相関係数

x(n), y(n) のデータからシンプルに、それぞれの平均と標準偏差、 相関係数を求めるサブルーチン。
ついでに、線形回帰も
テストのための主プログラムとして、1000行までの x y の値を書いた テキストファイルを標準入力から読み込むものをとりあえずつけてある。

プログラム (テスト用メインプログラム付き)

implicit none

integer, parameter :: num=1000
real :: x(num), y(num)
integer :: i, j, n
real :: xm, ym, sigx, sigy, cor
real :: a, b

  i=1
  do while( 1==1 )
    read(*,*, END=900) x(i), y(i)
    i=i+1
    if(i>num) then
      write(0,*) 'too many data. Change num'
      stop
    end if
  end do
900 n=i-1

  do i=1, n
   write(0,*) i, x(i), y(i)
  end do

 call correlation(x, y, n, xm, ym, sigx, sigy, cor)

  write(*,*) 'Mean x = ', xm, 'Sig x =', sigx
  write(*,*) 'Mean y = ', ym, 'Sig y =', sigy
  write(*,*) 'Correlation = ', cor

 call calc_ab( x, y, n, a, b)
  write(*,*) 'y=ax+b', 'a=', a, 'b=', b

end


subroutine correlation(x, y, n, xm, ym, sigx, sigy, cor)

 implicit none
 integer :: n, i
 real    :: x(n), y(n)
 real :: xm, ym, sigx, sigy, cor

 real(kind=8) :: sx, sy, sxy, sx2, sy2
  
 sx=0; sy=0; sxy=0; sx2=0; sy2=0


if( n > 2 ) then
 do i=1, n
  sx = sx +x(i)
  sy = sy +y(i)
 end do
  xm = sx/n; ym = sy/n
 do i=1, n
  sx2 = sx2 + (x(i)-xm)*(x(i)-xm)
  sy2 = sy2 + (y(i)-ym)*(y(i)-ym)
  sxy = sxy + (x(i)-xm)*(y(i)-ym)
 end do

  sigx = sqrt(sx2/n); sigy = sqrt(sy2/n)
  if( sigx > 0 .AND. sigy > 0 ) then
  cor = (sxy/n)/(sigx*sigy)
  end if
end if
end


subroutine calc_ab( x, y, n, a, b)
!   y = ax+b 
 implicit none
 integer :: n, i
 real :: x(n), y(n), a, b
 double precision :: sx2, sx, sxy, sy, d

 sx2 = 0 ;  sx = 0 ;  sxy = 0;  sy = 0

 do i=1,n
   sx2 = sx2 + x(i)*x(i) ! sx2 = sx2 + i*i
   sx = sx + x(i)        ! sx  = sx + i
   sxy = sxy + x(i)*y(i) !sxy = sxy + i*x(i)
   sy = sy + y(i)        !sy  = sy  + x(i)
 end do

 d = sx2 * n - sx * sx
 a = (n*sxy-sx*sy)/d
 b = (-sx*sxy+sx2*sy)/d
end

プログラム例をダウンロードするのはここ

戻る

玉川 2011.12.1