最小2乗法

図のようにV(電圧)とI(電流)の複数のデータから一定の関係を見つける問題を考える。

いいかえると次式の最確値a,bを決定する問題となる。

V = aI+b (a,b=定数)

このような問題の解法として最小2乗法がある。

まず1組の測定値が確実に得られるためにはどのようにするかを考えてみよう。

n回の測定値を

$x_1,x_2,\cdots, x_n$ とし、その最確値を
x_m とすると、 x_i の測定値をうる確率 $P_{xi}$

は、確率誤差関数に比例するから、

$P_{xi} \propto \frac{1}{\sigma\sqrt{2\pi}}\exp({-\frac{(x_i-x_m)^2}{2\sigma^2}})$

よって、 $x_1,x_2,\cdots,x_n$ なる測定値が得られる確率は、

$P_{x1,x2,\cdots,xn} \propto (\frac{1}{\sigma\sqrt{2\pi}})^n \exp({-\frac{1}{2\sigma^2}\bigl((x_1-x_m)^2 + (x_2 -x_m)^2 + \cdots + (x_n - x_m)^2\bigr)})$ (1-1)

一方、上の一組の測定値が現実に得られることは、得られる確率(1-1)式
において、 $P_{x1,x2,\cdots,xn}$ が最大になっているはずである。
そのためには分子が最小になればいいから、

$S=(x_1 - x_m)^2 + (x_2 - x_m)^2 + \cdots + (x_n - x_m)^2$
=最小 (1-2)

が成り立てばいい。

ここで例えば平均値 x_m を確実に得たいと考えた場合は、 x_m の偏微分を取ればいい。

$\frac{\partial~S}{\partial~{x_m}} = 0$ (1-3)

(1-3)式を計算すれば、
$x_m = \frac{1}{n}(x_1 + x_2 + \cdots + x_n)$ (1-4)

これは、測定値を全部足してそれを総数で割った平均と同じになる。

aとbの最確値を求める

次に実際にaとbの最確値を求めていこう。
測定量xとyとの間に、簡単のために次の関係があるとし、
xとyの測定値からaとbの最確値を決定する。

y = ax + b (2-1)

(2-1)式の直線を決めるためにn組の測定値 (x_i,y_i) が得られたとき、
平均 x_m から偏差を求めるように、 y_i からの偏差を考えると

$y_i - (ax_i + b) \neq 0$ (2-2)

となり、その確率誤差関数P(x)を求めると,

$P_{i} \propto \frac{1}{\sigma\sqrt{2\pi}}\exp({-\frac{(y_i-(ax_i+b))^2}{2\sigma^2}})$ (2-3)

のようになる。

ここから同じようにSが最小になるようにa,bを求めれば最確値a,bが決定する。

$S = \sum^{n}_{i=1}\{ y_i - (ax_i + b)\}^2$ (2-4)

a,bを最小化する条件は、

$\frac{\partial~S}{\partial~{a}} = 0$ 及び
$\frac{\partial~S}{\partial~{b}} = 0$ であり、

$\sum^{n}_{i=1}y_i = a\sum^{n}_{i=1}x_i + nb$ 　①
$\sum^{n}_{i=1}x_{i}y_i = a\sum^{n}_{i=1}x_i^2 + b\sum^{n}_{i=1}x_i$ ②　(2-5)

(メモ　①x $\sum{x_i}$ - ②xnでaを計算。)

そして上式を解いて、a,bは次のように決定される。

$a = \frac{ n\sum{x_i y_i} - (\sum{x_i})(\sum{y_i})}{ n\sum{x^2_i}-(\sum x_i)^2}$

$b = \frac{ (\sum{y_i})(\sum{x^2_i}) - (\sum{x_i y_i})(\sum{x_i})}{ n\sum{x^2_i}-(\sum x_i)^2}$

aが決定していればbは、①式より

$b=\frac{ \sum^{n}_{i=1}y_i - a\sum^{n}_{i=1}x_i }{n}$

とした方が簡単である。

(2-6)式

また標準偏差は、

$\sigma_y = \{ \frac{ \sum{(y_i - ax_i - b)^2}}{n} \}^{1/2}$ (2-7)

となる。

(例) xとyを測定してA及びBの最確値を求めよう。

y=Ax+B (A,Bは定数)においてxとyを測定して表1.を得たときのAとBの最確値を求める。

*x	1	2	3	4
*y	4.5	5.7	7.3	8.5

S=(4.5-A-B)^2 + (5.7-2A-B)^2 + (7.3-3A-B)^2 + (8.5-4A-B)^2

SをA,Bで編微分して0とおくと、次のようになる。

$\frac{dS}{dA} = -2(4.5-A-B)-4(5.7-2A-B)-6(7.3-3A-B)-8(8.5-4A-B)=0$
$\frac{dS}{dB} = -2(4.5-A-B)-2(5.7-2A-B)-2(7.3-3A-B)-2(8.5-4A-B)=0$

60A+20B-143.6=0.20

A+8B-52.0 = 0

このようにして得られた両式からA=1.36,B=3.10となる。

最小2乗法　Fortran90でのコード

 ! least squares method
 PROGRAM leastSquaresMethod
 implicit none

 integer::open_status ! determination open flag
 integer:: i,n
 real*8::x(100),y(100),sumx,sumy,sumxx,sumxy,a,b
 integer::buf

 write(*,*) "calculate using  least squares method"

 !input file
 open(9,file='input.txt',iostat=open_status) 
 if(open_status>0) stop "Input error"

 !input data

 write(*,*) "input data from input.txt"
 n=4
 read(9,*) (x(i),y(i),i=1,n)
 close(9)

 !!!!!!!!!!!!!!!!!!!!!!!!!!
 ! least squaares method
 !!!!!!!!!!!!!!!!!!!!!!!!!

 sumx=0
 sumy=0
 sumxx=0
 sumxy=0

 do i=1,n
    sumx = sumx + x(i)
    sumy = sumy + y(i)
    sumxx = sumxx + x(i)**2
    sumxy = sumxy + x(i)*y(i)
 end do

 a=(n*sumxy-sumx*sumy) / (n*sumxx-sumx**2)
 b=(sumy-a*sumx)/n

 write(*,777) 'a=',a,'b=',b
 777 format(' ',2(a5,f8.3))

 close(9)

 stop
 end PROGRAM leastSquaresMethod

結果

./test
calculate using  least squares method
input data from input.txt
a=   1.360   b=   3.100

最小2乗法

最小2乗法

aとbの最確値を求める

(例) xとyを測定してA及びBの最確値を求めよう。

最小2乗法 Fortran90でのコード

結果

最小2乗法　Fortran90でのコード