Link Copied
Doc:
There are a lot of humans watching this forum, in 98 minutes 38 people have viewed this last missive, which is one every 2.5 minutes. I mean that is a high hit rate for a small forum of limited interest to a very narrow group of elderly Fortran programmers.
Continuing the :: am I having fun:
Do not steal other people's code unless you are going to check it.
Case in point, I needed to do a quick linear regression analysis on linear and then log mapped data. I grabbed this code,
SUBROUTINE MAP_UTIL_ANALY(J8,J9,CONST,SLOPE,CORR,X,Y)
REAL X(J8),Y(J8),AX,AY,AX2,AY2,AXY
REAL CONST,SLOPE,CORR
INTEGER I1,J8,J9,I
I1=J8
J8=J8-J9
AX=0.0D+00
AY=0.0D+00
AX2=0.0D+00
AY2=0.0D+00
AXY=0.0D+00
DO 1 I=1,J8
AX=AX+X(I)
AY=AY+Y(I)
AXY=AXY+X(I)*Y(I)
AX2=AX2+X(I)*X(I)
AY2=AY2+Y(I)*Y(I)
1 CONTINUE
CONST=(AY*AX2-AX*AXY)/(J8*AX2-AX*AX)
SLOPE=((J8*AXY-AX*AY)/(J8*AX2-AX*AX))
CORR=(J8*AXY-AX*AY)/(SQRT((J8*AX2-AX*AX)*(J8*AY2-AY*AY)))
J8=I1
write(*,100)CONST,SLOPE,CORR,CORR*CORR
100 Format(" Constant m :: ",F20.15,/,&
" Slope m/m :: ",F20.15,/,&
" Correlation R :: ",F10.5,/,&
" Correlation RSq :: ",F10.5,/)
RETURN
END
I added the write statement.
I am at the stage of checking the answer, so I use the tried and true EXCEL for simple data set.
The code gives me R as a return value. One is so used to the traditional Rsquared being returned and this is a single independent variable model that is all you want. I would never expect the R instead of the Rsquared.
So I grabbed this code to check the map code.
!****************************************************
! 相関係数計算
!
! date name version
! 2018.12.18 mk-mode.com 1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
module const
! SP: 単精度(4), DP: 倍精度(8)
integer, parameter :: SP = kind(1.0)
integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
end module const
module comp
use const
implicit none
private
public :: calc_r
contains
! 相関係数計算
!
! :param(in) real(8) x(:): X 値配列
! :param(in) real(8) y(:): Y 値配列
! :return real(8) r: 相関係数
function calc_r(x, y) result(r)
implicit none
real(DP), intent(in) :: x(:), y(:)
real(DP) :: r
integer(SP) :: size_x, size_y, i
real(DP) :: mean_x, mean_y, cov, var_x, var_y
size_x = size(x)
size_y = size(y)
if (size_x == 0 .or. size_y == 0) then
print *, "[ERROR] array size == 0"
stop
end if
if (size_x /= size_y) then
print *, "[ERROR] size(X) != size(Y)"
stop
end if
r = 0.0_DP
mean_x = sum(x) / size_x
mean_y = sum(y) / size_y
cov = sum((x(1:size_x) - mean_x) * (y(1:size_x) - mean_y))
var_x = sum((x(1:size_x) - mean_x) * (x(1:size_x) - mean_x))
var_y = sum((y(1:size_x) - mean_y) * (y(1:size_x) - mean_y))
r = (cov / sqrt(var_x)) / sqrt(var_y)
end function calc_r
end module comp
program correlation_coefficient
use const
use comp
implicit none
real(DP) :: x(10), y(10), r
integer(SP) :: i
x = (/(i, i=1,10)/)
y = (/(i, i=1,10)/)
r = calc_r(x, y)
print '("X = (", 10F6.2, ")")', x
print '("Y = (", 10F6.2, ")")', y
print '("r = ", F11.8)', r
print *
y = (/2, 3, 3, 4, 6, 7, 8, 9, 10, 11/)
r = calc_r(x, y)
print '("X = (", 10F6.2, ")")', x
print '("Y = (", 10F6.2, ")")', y
print '("r = ", F11.8)', r
print *
y = (/15, 13, 12, 12, 10, 10, 8, 7, 4, 3/)
r = calc_r(x, y)
print '("X = (", 10F6.2, ")")', x
print '("Y = (", 10F6.2, ")")', y
print '("r = ", F11.8)', r
end program correlation_coefficient
The Chinese comment says Pearson coefficient, so this is also R as one observes by looking at the equations.
As some one observed to me yesterday, how do you know the code is correct? I looked and said, you check against known data and methods and you look at the graphs and say - does it look ok.
Fun - yes a pleasant interlude, like reading a murder mystery.
of course I should not have been lazy and use the NR in Fortran.
