最小二乗法の理論

9. 最小 2 乗法
田中雅博
最適化プログラミング、教科書 pp.105–113
1 次関数の当てはめ
1
1.1
変数が 1 つの場合
いま x とそれに関係する y のデータが N 組ある。すなわち、{(x1 , y1 ), . . . , (xN , yN )} とする。
これらが直線の上に乗っていれば、y = ax + b と書ける。しかし、一般にはばらつく。そこで、y
を x の 1 次式で近似するという意味で、誤差の 2 乗を最小にすることを考える。すなわち、
1∑
2
(yi − (axi + b)) → min
2 i=1
N
J=
(1)
これを、a 及び b を調整して最小化するということからは、
∂J
= 0,
∂a
∂J
=0
∂b
とすればよい。
(1) 式を直接そのまま a で微分して 0 とおくことにより、
N
∑
(yi − (axi + b)) (−xi ) = 0
i=1
これより
−
N
∑
yi xi +
i=1
N
∑
axi xi +
i=1
N
∑
bxi = 0
i=1
すなわち
a
N
∑
x2i + b
i=1
N
∑
xi =
i=1
N
∑
yi xi
(2)
i=1
を得る。同様に、(1) 式を直接そのまま a で微分して 0 とおくことにより、
N
∑
(yi − (axi + b)) = 0
i=1
よって
a
N
∑
xi + bN =
i=1
N
∑
i=1
1
yi
(3)
を得る。(2) 式と (3) 式を連立させることにより
) ( ∑
)
( ∑
)(
∑N
N
N
2
a
i=1 xi yi
i=1 xi
i=1 xi
=
∑N
∑N
b
N
i=1 yi
i=1 xi
(4)
となり、これ（正規方程式と呼ばれている）を解くことにより、a と b が求まる。
( ∑
)
∑N
N
2
i=1 xi
i=1 xi
∑N
N
i=1 xi
が逆行列をもつとき
(
a
b
)
( ∑
N
x2i
∑i=1
N
i=1 xi
=
∑N
i=1 xi
N
)−1 ( ∑
N
xi yi
∑i=1
N
i=1 yi
)
(5)
となる。
教科書はこの形を示しているが、行列とベクトルを使った表現の方が簡潔でわかりやすい。先に
述べたように、y は厳密に ax + b とは書けないので、各データ (xi , yi ) ごとに誤差 ei を持つものと
しよう。そうすると、
yi = axi + b + ei ,
i = 1, . . . , N
と書ける。データを縦に並べたベクトルや、それに関連する行列を




y1
x1 1
(
)




 y2 
 x2 1 
a




y =  . , M =  .
,
..  , θ =
b
 .. 
 ..
. 
yN
xN
1



e=


e1
e2
..
.






eN
とおくと、
y = Mθ + e
と書ける。誤差の 2 乗和は
N
∑
(6)
e2i = eT e
i=1
であるから、
J=
1 T
e e = (y − M θ)T (y − M θ)
2
となる。内積表現を用いると
J=
1
1
(e, e) = ((y − M θ), (y − M θ))
2
2
1
1
1
1
(y, y) − (y, M θ) − (M θ, y) + (M θ, M θ)
2
2
2
2
1
1
= (y, y) − (M θ, y) + (M θ, M θ)
2
2
1
1
= (y, y) − (θ, M T y) + (θ, M T M θ)
2
2
となる。ここで、パラメータに関して最小化を行うために、勾配ベクトル ∇J(θ) = 0 とおく。教
科書定理 1.5, 1.6 より
=
∇J = −M T y + M T M θ = 0
2
となるので、M T M が逆行列を持てば
θ = (M T M )−1 M T y
となる。この式は (5) 式と同等である。
’ 例題 4.2’
X=[4 15 30 100 200]’;
M=[X ones(5,1)];
Y=[ -17 -4 -7 50 70]’;
th=inv(M’*M)*M’*Y
plot(X’,Y’,’o’)
hold on
X=(0:250)’;
M=[X ones(251,1)];
Yest=M*th;
figure(1)
plot(X’,Yest’,’r’)
xlabel(’x’)
ylabel(’y’)
hold off
120
100
80
y
60
40
20
0
-20
0
50
100
150
200
250
x
図 1: 変数が 1 つの場合の最小二乗法による推定値
1.2
変数が 2 つ以上の場合
いま，変数 y の値を，l 種類の説明変数 x1 , . . . , xl を用いて推定したいとする．x1 , . . . , xl の値は
容易に取得できるものとする．このとき，y はどのようにして推定できるのであろうか．もっとも
3
簡単なモデルは，線形回帰モデル
y = a0 + a1 x1 + · · · + al xl
(7)
とするものである．係数 a0 , a1 , . . . , al は各変数がそれぞれ y に及ぼす影響の強さを示す．この場
合，y は誤差を含まないが，一般的に，
「真」の y は x1 , . . . , xl で正確に決まることはほとんどなく，
誤差を含むことから，(7) 式を改めて
y = a0 + a1 x1 + · · · + al xl + e
(8)
と書く．ここに，e は誤差項である．これを、重回帰モデルともいう。
重回帰モデルを推定するには，まず，個々のデータに対する (8) 式を縦に並べて改めてベクトル
を太字，行列を大文字で表現した式
y = Mθ + e
(9)
を定義しよう．ただし，上付き y (j) は，j 番目のデータに対する y を示す．他の変数の上付きも同
様である．また

y (1)
 (2)
 y
y=
 ..
 .
y (N )


1



 1
, M =  .

 .

 .
(1)
x2
(2)
x2
..
.
(N )
x2
x1
x1
..
.
1 x1
とおいた．ここで，誤差の 2 乗和
∑N
(1)
(1)
···
xl
(2)
···
..
.
xl
..
.
(N )
· · · xl
(j) 2
j=1 (e )
(2)





, θ = 



(N )


e(1)
a0
 (2)
 e
.. 
, e =  .
. 
 .
 .
al
e(N )






= (y − M θ)⊤ (y − M θ) を最小にすることを考えよ
う．2 乗誤差の極値を与えるパラメータを、1 変数の場合と同じ手順で求めると、
θ̂ = (M ⊤ M )−1 M ⊤ y
(10)
で与えられる．
この場合に，注意が必要なのは，データの個数と独立性である．逆行列を構成する行列 M のサ
イズは N × (l + 1) なので，N < l + 1 の場合は，M の階数が高々N であり，M ⊤ M は決して逆行
列をもたない．つまり，データの個数が，データの次元に満たないときは，回帰モデルは最小２乗
法では係数を一意的に決定できない．
では，データ個数 N が N ≥ l + 1 のときは，いつも θ が求まるのであろうか．実はここで，多
重共線性といわれる問題がある．つまり，説明変数 1, x1 , . . . , xl のうちで一次従属な変数がある
と，いくらデータをたくさん集めても（すなわち，N を大きくしても），M の階数は l 以下にな
り，M ⊤ M は逆行列をもたない．また，たとえ逆行列をもつ場合でも，行列式が 0 に近ければ，得
られた結果は信頼がおけない．
2
多項式の当てはめ
N 個のデータ {(x1 , y1 ), . . . , (xN , yN )} に 2 次式を当てはめよ。
当てはめる式を y = ax2 + bx + c とすると、今までの議論同様、
)2
1 ∑(
yi − (ax2i + bxi + c) → min
J=
2 i=1
N
4
これを最小化するために、(1.2) 同様
y = ax2i + bxi + c + e
とし、全部のデータを縦に並べると


y1
 . 
. 
Y =
 . ,
yN

x21
 .
.
M =
 .
x2N

1
.. 

. 
1
x1
..
.
xN
θ = (M T M )−1 M T Y
X=[-1 0 0 1]’;
Y=[0 -2 -1 0]’;
M=[X.^2 X ones(4,1)];
th=inv(M’*M)*M’*Y
figure(3)
plot(X’,Y’,’o’)
X=(-1.5:0.1:1.5)’;
M=[X.^2 X ones(31,1)];
hold on
Yest=M*th;
plot(X’,Yest’,’r’)
xlabel(’x’)
ylabel(’y’)
hold off
2
1.5
1
y
0.5
0
-0.5
-1
-1.5
-2
-1.5
-1
-0.5
0
0.5
x
図 2: 二次曲線近似
5
1
1.5
表 1: 課題用データ
x
y
z
課題
0.4387
0.2513
1.3825
0.3816
-0.2449
1.1936
0.7655
0.0060
1.5376
0.7952
0.1991
1.6164
0.1869
0.3909
1.2481
0.4898
0.4593
1.4806
0.4456
0.0472
1.3261
0.6463
-0.3614
1.3440
0.7094
-0.3507
1.3913
0.7547
-0.2425
1.4555
0.2760
0.3407
1.2954
0.6797
-0.2457
1.4021
0.6551
0.3143
1.5529
0.1626
-0.2565
1.0369
0.1190
0.4293
1.2121
0.4984
-0.1500
1.3038
0.9597
-0.3034
1.5808
0.3404
-0.2489
1.1636
0.5853
0.1160
1.4445
0.2238
-0.0267
1.1487
表１のデータ (http://carnation.is.konan-u.ac.jp/opt9kadai.txt として保存してある）にお
いて、最小 2 乗法により z を x と y の 1 次式による回帰モデルとして表現せよ。
[ファイルの読み方] まず、データをブラウザで読み、MATLAB のカレントディレクトリ（プロ
グラムが置いてあるところ）に保存する。次に、
data=load(’opt9kadai.txt);
とすれば、データが配列 data に入る。あとは、その 1 列目を x、2 列目を y、3 列目を z として使
えばよい。
[提出物] Word に、以下の内容を示せ。
• 学籍番号、氏名、課題出題日
• 作成した MATLAB プログラムリスト
• 得られた回帰式
• データ z、その推定値、誤差の 2 乗平均（並べて表示）
• 感想
6

Download Report