最小二乗法 (.pdf)

最小二乗法
2008年3月24日(月)
最小二乗法(最小自乗法)

実験からn個のデータ (x1, y1), (x2, y2), …, (xn, yn)を取得

これらのデータに最もフィットする直線を f(x)=a0+a1x と
したとき a0, a1 の値は?
y
y
f(x)=a0+a1x
最小二乗法
(x4, y4) (x5, y5)
(x3, y3)
(x1, y1)
0
(x2, y2)
x
直線になるはずが、誤差によりバラつきが生じる
0
x
本来はどのような直線になるのか?
求める f(x)

各測定データと直線との誤差をdiとする

G=d12+d22+…+dn2 としたとき、Gを最小にするf(x)を
求める関数とする
di2
G=
N
∑
i= 1
=
N
∑
i= 1
( yi − f ( xi ))
f(x)=a0+a1x
2
yi
( yi − a0 − a1 xi )
2
f(xi)
a0
(xi, yi)
di
n=3 の例 (1/3)
G=
3
∑
i= 1
( yi − a0 − a1 xi ) 2
= ( y1 − a0 − a1 x1 ) 2 + ( y2 − a0 − a1 x2 ) 2 + ( y3 − a0 − a1 x3 ) 2
2
2
2
2
= ( y1 + a0 + a1 x1 − 2a0 y1 − 2a1 x1 y1 + 2a0 a1 x1 )
2
2
2
2
2
2
2
2
+ ( y2 + a0 + a1 x2 − 2a0 y2 − 2a1 x2 y2 + 2a0 a1 x2 )
+ ( y3 + a0 + a1 x3 − 2a0 y3 − 2a1 x3 y3 + 2a0 a1 x3 )
A
2
2
2
2
B
2
2
2
= ( y1 + y2 + y3 ) + a1 ( x1 + x2 + x3 ) + 3a0
2
− 2a0 ( y1 + y2 + y3 ) − 2a1 ( x1 y1 + x2 y2 + x3 y3 ) + 2a0 a1 ( x1 + x2 + x3 )
C
D
E
2
2
G = A + a1 B + 3a0 − 2a0C − 2a1 D + 2a0 a1 E
n=3 の例 (2/3)

xi, yi は実験により得られる値=A, B, C, D, E は定数
Gはa0, a1の関数
G
0
G
a0_min
a0
0
a1_min
a0
a1を固定すると、Gはa0の2次関数
a0を固定すると、Gはa1の2次関数
→Gが最小となるa0_minが存在
→Gが最小となるa1_minが存在
n=3 の例 (3/3)

最小値を探すため、微分して0とおく
 ∂G
 ∂ a = 6a0 − 2C + 2a1 E = 0
 0

 ∂ G = 2a B − 2 D + 2a E = 0
1
0
 ∂ a1

 3a0 − C + a1 E = 0

 a1 B − D + a0 E = 0
これを解くと、
C− D
3D − CE
a0 =
, a1 =
2
3− E
3E − E
一般的に (1/2)
G=
n
∑
i= 1
A n
=
∑
i= 1
( yi − a0 − a1 xi ) 2
B n
2
y i + a0
2
2
∑
i= 1
C n
D n
E n
xi + na1 − 2a1 ∑ yi − 2a0 ∑ xi yi + 2a0 a1 ∑ xi
2
2
i= 1
i= 1
i= 1
2
= A + a0 B + na1 − 2a1C − 2a0 D + 2a0 a1 E
微分して
 ∂G
 ∂ a = 2 Ba0 − 2 D + 2a1 E = 0
 0

 ∂ G = 2na − 2C + 2a E = 0
1
0
 ∂ a1
 Ba0 − D + a1 E = 0

 na1 − C + a0 E = 0
一般的に (2/2)
連立方程式を解いて
nD − CE
BC − DE
a=
, b =
2
nB − E
nB − E 2
元の値を代入して
n
a=
n∑ xi yi −
i= 1
n
n∑
i= 1
n
n
∑ x∑
i
i= 1
i= 1


xi −  ∑ xi 
 i= 1 
2
n
n
yi
2
, b =
∑
i= 1
xi
2
n
∑
i= 1
n
n∑
i= 1
yi −
n
∑
i= 1
xi yi ∑ xi
 n 
xi −  ∑ xi 
 i= 1 
2
n
i= 1
2
参考: n次式での近似(1/2)

2次式や3次式の多項式で近似したい場合も同様
…しかし次数が増えるにつれ計算が困難に
ガウスの消去法(連立方程式を解くプログラム)を用いる
f ( x) = a0 + a1 x + ... + an x
G=
N
∑
i= 1
n
( yi − f ( xi )) 2
n 次の多項式 f(x) を G に代入し、
∂G
∂G
∂G
= 0, = 0, ..., = 0 を考える
∂ a0
∂ a1
∂ an
参考: n次式での近似(2/2)
s をデータの個数とし、行列 A、行列 a、行列 b の各成分をそれぞれ
n
n
∑∑
i= 0 j= 0
Aij =
s
∑
k= 0
xk
n
∑
i+ j
i= 0
bi =
s
∑
k= 0
i
xk y k
n
∑
i= 0
ai = ai
とすると、n 個の連立1次方程式は
Aa = b
と表される。
これは未知数が n 個あり、n 個の方程式からなる連立1次方程式であるから、
ガウスの消去法で解ける
参考: ガウスの消去法(掃き出し法)

変数の数と方程式の本数が一致した連立一次方程式を
解くための解法

前進消去、後退代入の2ステップからなる
(例)
 2x + 4 y + 2z = 8

 4 x + 10 y + 3z = 17
 3x + 7 y + z = 11

を解くことを考える
 2 4 2 


 4 10 3  
 3 7 1 


x  8 
  
y  =  17 
z   11 
参考: ガウスの消去法 - 前進消去
 2 4 2 


4
10
3


 3 7 1 


x  8 
  
y  =  17 
z   11 
1行目を1/2倍
 1 2 1 


4
10
3


 3 7 1 


1行目を-4倍して2行目に加算
1行目を-3倍して3行目に加算
 1 2 1 


 0 2 − 1 
 0 1 − 2 


2行目を1/2倍
1 
1 2


 0 1 − 12  
 0 1 − 2 


x   4 

y =  1 
 2

z  − 1


2行目を-1倍して3行目に加算
3行目を-2/3倍
1 
1 2


 0 1 − 12  
0 0
1  

x   4 

y =  1 
 2
z   1 


x  4 
  
y  =  17 
z   11 
x  4 
  
y =  1 
z   − 1
参考: ガウスの消去法 - 後退代入
1 
1 2


1
0
1
−

2
0 0
1  

x   4 

y =  1 
 2 

z
1 
3行目を1/2倍して2行目に加算
3行目を-1倍して1行目に加算
2行目を-2倍して1行目に加算
 1 2 0   x   3

   
0
1
0

 y =  1
 0 0 1  z   1 

   
 1 0 0   x   1

   
0
1
0

  y  =  1
 0 0 1   z   1

   
∴ x = 1, y = 1, z = 1