トレンド除去

x(i) のデータがある時、１次近似で x = a*i+b と a, b を 最小自乗法で求め (calc_ab) 、それを使って、線形のトレンドからの 偏差を x に戻す ( detrend ) サブルーチンの例。 テストのための主プログラムをとりあえずつけてある。

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

implicit none

integer, parameter :: num=100
real :: x(num)
integer :: i, j, n

  i=1
  do while( 1==1 )
    read(*,*, END=900) x(i)
    i=i+1
  end do
900 n=i-1

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

  call detrend( x, n)

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

end


subroutine detrend( x, n )
 implicit none
 integer :: n
 real    :: x(n)
 real    :: a, b
 integer :: i

 call calc_ab( x, n, a, b )

 write(*,*) 'a=', a, 'b=', b

 do i=1, n
  x(i) = x(i) - (a*i+b)
 end do
end


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

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

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

 d = sx2 * n - sx * sx
 a = (n*sxy-sx*sy)/d   ! (2008.1.31 にはここにタイプミスがありました)
 b = (-sx*sxy+sx2*sy)/d
end

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

戻る

玉川 2008.1.20