ネムネコの部屋: 最小２乗法による近似式の決定

最小２乗法による近似式の決定

§１　最小２乗法による近似多項式の求め方

n個のデータ

があるとする。

このデータがm次の多項式（m≦n+1）

に沿っているとし、

が最小になるようにm次多項式（１）の係数

を定める方法を最小２乗法という。

最小２乗法では、

とおいたとき、（２）式が最小になるように

を定めるので、

にならなければならない。

したがって、

でなければならない。

ところで、

となるので、（４）式は

ここで、

とおくと、（５）式は

となる。

したがって、この連立方程式（７）を解くことによって、多項式（１）の係数を定めることができる。

特に、m=1のとき、すなわち、

のとき、（７）は

ここで、

となるので、（８）は

これを解くと、

ここで、

したがって、統計学の回帰直線と一致する。

§２　プログラム

n個のデータ

をもとに、最小２乗法を用いてm次の近似多項式を求める連立方程式は（６）、（７）式から与えられる。

（６）、（７）から与えられる連立方程式をガウス消去法を用いて解くことにすると、例えば、次のようなプログラムを作ることができる。

real x(20),y(20) !xとyのデータは２０個まで
real c(0:10) ! 近似多項式の係数（１０次まで）
real y_c(20)
data x/20*0/ ! x(i)の初期化
data y/20*0/ ! y(i)の初期化

write(6,*) "データの個数nを入力"
read(5,*) n

write(6,*) "データの入力"
read(5,*) (x(i), i=1,n) ! xのデータをn個入力
read(5,*) (y(i), i=1,n) ! yのデータをn個入力

write(6,*) "近似多項式の次数mを入力"
read(5,*) m

call lsm(x,y,n,m,c)

write(6,*)
! 近似式の係数の表示
write(6,*) "近似多項式の係数"
do i=0, m
    write(6,100) 'c(',i, ')=', c(i)
end do

write(6,*)
! 計算結果
write(6,*) " x(i)       y(i)        計算値"
do i= 1, n
    s = 0
    do j=0, m
        s = s + c(j)*x(i)**j
    end do
    y_c(i)=s
    write(6,110) x(i), y(i), s
end do

open(1,file='lsm.dat', status='replace')
do i=1,n
    write(1,120) x(i),y(i),y_c(i)
end do

100 format(a,i2,a3,f10.7)
110 format(f10.7,2x,f10.7,2x,f10.7)
120 format(f10.7,x,f10.7,x,f10.7)
end

! ここから下は変えないほうがよい(^^)
! 最小二乗法
subroutine lsm(x,y,n,m,c)
real a(0:50,0:51)
real x(20), y(20), c(0:10)
real s(0:20), t(0:10)

! 係数の初期化
do i=0, 2*m
    s(i) = 0
end do
do i=0, m
    t(i)=0;
end do

! 係数を計算
do i=1, n
    do j=0, 2*m
        s(j) = s(j)+x(i)**j
    end do
    do j=0, m
        t(j)= t(j)+x(i)**j * y(i)
    end do
end do

! 係数行列に値をセット
do i= 0, m
    do j=0, m
        a(i,j) = s(i+j)
    end do
    a(i,m+1) = t(i)
end do

! ガウス消去法で連立方程式を解く
call Gauss(a,m)

do i=0, m
    c(i) = a(i,m+1)
end do

end

subroutine Gauss(a,n) ! ガウス消去法
real a(0:50,0:51)

! foward marching process
do k=0, n-1
! ピボット選択
    i_max = k
    do j=k+1, n
        if (abs(a(k,i_max)).lt.abs(a(k,j))) i_max = j
    end do
    if (i_max.ne.k) then
        do j=k,n+1
            w=a(i_max, j)
            a(i_max, j) = a(k,j)
            a(k,j) = w
        end do
        det=-det
    end if
    if (abs(a(k,k)).lt.eps) then
        det =0.
        return
    end if
! ピボット選択終わり
    do i=k+1, n
        t=a(i,k)/a(k,k)
        do j=k+1, n+1
            a(i,j)=a(i,j)-t*a(k,j)
        end do
    end do
end do

! backward marching process
do i=n,0, -1
    d = a(i,n+1)
    do j= n, i+1, -1
        d=d- a(i,j)*a(j,n+1)
    end do
    a(i,n+1)=d/a(i,i)
end do
end

たとえば、次の３個のデータ

(1,1) , (2,5), (3,8)

をもとに２次の近似多項式を求めると、次のようになる。