最小二乗法による線形モデルのパラメータ決定中野武雄

2005 年 6 月 4 日
概 要
線形最小二乗法の定義と、そのモデルパラメータの決定法について述べる。得られたモデルパラメータの誤差の求め方についても記し、最後にいくつかのケースに関して具体的な表式を示す。

1

線形最小二乗法

1.1

測定と解析モデルの定義

ある実験において n 回の測定が行われた場合を考える。一つ一つの測定に際しては、l 個の測定条件が (数値によって) 指定されており、1 つの測定値が得られるものとする。i 番目の測定
(1)

(2)

(l)

(1 < i ≤ n) における条件パラメータの組を xi , xi , . . . , xiのように表し、得られる測定値をyi とする。測定値 yi は、おのおの分散 σi を持って、真の値 yi 0 の周りを正規分布しているもの(j)とする。また xi には誤差は含まれないものとする。
(j)
このとき、yi 0 と xi が、あらかじめ導入したモデル(1)(2)(l)yi 0 = f (xi , xi , . . . , xi ; θ1 , θ2 , . . . , θm )によって誤差なく関係付けられるものと仮定したとき、このモデルのパラメータ θkを決定することを考える。パラメータ推定値はベクトルの形式に
noθ̂ ≡ θˆ1 , . . . , θˆm(1)(k = 1 ∼ m)(2)と書き、この推定値からモデルを用いて計算される値を

yˆi=
≡(1)
(l)
f (xi , . . . , xi ; θˆ1 , . . . , θˆm ) (i = 1 ∼ n)fi (θˆ1 , . . . , θˆm ) ≡ fi (θ̂)(3)のように表わすものとする。
最小二乗法とは、このパラメータ推定法のひとつであり、具体的には以下の最小二乗条件を満たすように θ̂ を定める方法である。

S(θ̂) ≡nX[yi − fi (θ̂)]2σi 2i=1= min(4)ここで、測定値と計算値との差

vi ≡ yi − yˆi = yi − fi (θ̂) (i = 1 ∼ n)1(5)を残差と呼ぶ。また次の式で各測定値 yi の重み wi を定義する。

wi ≡ σ0 2 /σi 2(6)ここで σ0 は任意に選べる比例定数であり、測定 i における分散 σi 2 が各測定についてそれぞれ既知である場合には、単純に 1 と置いて構わない。もし σi 2 が未知のときは、代表的と思われる値を設定しておき、式 6 に従って各測定値における分散の逆数の相対値を重みとすれば良い。
式 5 と式 6 を用いると、式 4 の最小二乗条件は、

S 0 (θ) ≡ σ0 2 S(θ) =nX2vi (θ) wi = min(7)i=1と書き直すことができる。
最小二乗法がパラメータ θ̂ の最尤推定法であるためには、以下の前提が必要となる。

1. 測定値の誤差 ²i = yi − yi 0 は不偏である。すなわち、h²i i = 0 である。
2. 測定値の誤差の分散は既知である。
3. 各測定は互いに独立であり、共分散は 0 である。
4. 誤差の分布系は正規分布である。
実際には条件 2 が満たされることは困難であるが、最小二乗法では測定値の誤差の相対値のみでモデルパラメータを決定することができる。ただしその場合、モデルの有意性やパラメータの決定誤差を求めることはできない。

1.2線形モデル

線形モデルとは、式 1 において、理論式 f がモデルのパラメータ θk に対する一次結合(1)(l)f (xi , . . . , xi ; θ1 , . . . , θm )≡ fi (θ1 , . . . , θm )= Xi1 θ1 + · · · + Xim θmmX=
Xij θj (i = 1 ∼ n)(8)j=1(1)(l)で表されるものをいう。ここで各係数 Xij は、条件パラメータ xi , . . . , xiを使って表される既知の定数であり、モデルパラメータ θi には依存しないものとする。つまり『線形』モデルとは、あくまでモデルのパラメータに対しての意味であることに注意する。1
この線形モデルを最小二乗条件である式 4 に代入すると、S(θ) =nXi=1yi −mX
j=1

2

Xij θj  /σi 21 例えばフーリエ展開で展開係数を計算するのも線形最小二乗問題の一種であると言える。

2(9)が得られる。式 9 が各モデルパラメータ θj に関して最小になるためには、各パラメータに対する偏微分が 0 になれば良い。

mnXX∂S(θ)
yi −Xij 0 xj 0  · Xij /σi 20== −2∂θj0i=1(j = 1 ∼ m)

j =1

これを未知のモデルパラメータ θj 0 に対して整理すると、
Ã n!
mn
XXX
¢
¡
2
Xij 0 Xij /σi θj 0 =Xij /σi 2 yi
j 0 =1i=1i=1(j = 1 ∼ m)

(10)(11)となる。これは θj 0 に関する以下のような連立一次方程式になる。

A11 θ1
 A21 θ1
ただしここで、Am1 θ1+A12 θ2+..
.

A22 θ2+Am2 θ2

Ajj 0 ≡nX

aj ≡+ ··· ++ ··· +i=1i=1=a1A2m θm=a2+ · · · + Amm θm

Xij 0 Xij /σi 2nX¡

A1m θm(12)= am(j, j 0 = 1 ∼ m)

(13)

(j = 1 ∼ m)

(14)

¢Xij /σi 2 yiである。これらの量は、既知の値 (測定条件パラメータと測定値) で表されていることに注意。式 10を正規方程式 (normal equation) と呼び、この方程式の解が最小二乗法で決定されるパラメータ推定値 θ̂ となる。

1.3

行列表式

以上の表式を簡潔に表わすためには、行列形式を用いるのが便利である。測定値や測定パラメー
タなど、これまで出てきた組になっている量を、今後列ベクトル (縦ベクトル) で表わす。ベクトルは小文字の太字、行列は大文字の太字で表し、ベクトルや行列の転置には ˜ を用いる。
測定値のベクトルを y とし、その真値を y0 とおいて、誤差ベクトル ² を次のように定義する。
これらのベクトルの次元は測定回数である n となる。

² ≡ y − y0(15)さらに、この n 個の測定に関する誤差行列 Σ をΣ ≡ h²²̃iと定義する。すなわち Σii0 ≡ h²i ²i0 i対角成分を共分散と呼ぶ。

(16)(i, i0 = 1 ∼ n) である。この誤差行列の対角成分を分散、非最小二乗法における、各測定値の誤差に関する要請から、誤差ベクトルと誤差行列は次の性質を持っていることがわかる。

h²i = 03(17)


Σ=


σ1 20σ2 2..

.
σn 20



(18)

さらに、式 6 で定義されていた、重み wi を対角成分とする重み行列 W を、誤差行列 Σ の逆行列を用いてW ≡ σ0 2 Σ−1(19)と定義する。測定値間に相関があるとき (共分散があり、Σ が非対角行列となるとき) も、式 19 はそのまま成立し、この後の表式を用いれば良い。
パラメータの列ベクトル (m 元) を θ 、モデル理論式の列ベクトル (n 次元) を f (x; θ) ≡ f (θ)と書くことにする。すると最小二乗問題は、誤差 σ (より一般的には誤差行列 Σ) のもとで、y∼= f (θ)(20)に対する最適パラメータ θ̂ を求めることであると表現できる。

θ̂ の解は、式 4 を一般化した以下の最小二乗条件
S(θ) ≡ ṽ(θ)Σ−1 v(θ) = σ0 −2 ṽ(θ)W v(θ) = min(21)を満たすものとして定義されることになる。
ここから線形モデルの話になる。式 8 で表される線形モデルは、行列の表現ではf (θ) ≡ Xθ(22)と表される。ここで X は式 8 の係数 Xij を要素とする n × m 次の長方行列で、ヤコビアン行列と呼ばれる。一般に m 個のパラメータを決定するにはそれ以上の回数の測定が必要であるから、n > m であることに注意する。
より一般的には、ヤコビアン行列の要素 Xij は、モデル f i (θ) のパラメータ θi に対する偏微分係数であると定義されており、

Xij ≡∂fi (θ)∂θj(i = 1 ∼ n, j = 1 ∼ m)(23)である。これは線形モデルの場合には、測定の条件パラメータととモデルとから既知となる定数である。
最小二乗条件 (式 9) を微分して得られる正規方程式は、式 10 を行列形式で表し、³´X̃W X θ = X̃W y(24)と書ける。あるいは、式 12 の形式で書けば、Ax = a(25)A ≡ X̃W X(26)a ≡ X̃W y(27)となる。ただしここで、4である。
したがって最小二乗解 θ̂ は、

θ̂ = A−1 a = (X̃W X)−1 X̃W y(28)である。ここで、m × n の行列 C を³´−1X̃W
C ≡ A−1 X̃W = X̃W X(29)θ̂ = Cy(30)のように定義すれば、結局式 28 はとなり、結局線形最小二乗法の解 θ̂ は、測定値 y の線形結合で表されることがわかる。

2

モデル有意性の検定
線形最小二乗解を求める際に利用したモデルが妥当であるかどうかについては、各測定値の分散がわかっている場合は χ2 検定によって検定することができる。(To be written: 参考文献の §3.5を参照)

3

決定パラメータの誤差前提として、式 15 で定義される誤差ベクトル ² は既知とする。また測定値の真値 y 0 と、パラメータ推定値の真値 θ 0 は y 0 = Xθ 0 という関係にあるとする。このとき

Cy 0 = CXθ 0 = (X̃W X)−1 X̃W Xθ 0 = θ 0(31)だから、式 30 に y 0 を用いれば θ 0 が得られることがわかる。
このときパラメータ推定値 θ̂ の真値に対するずれ δθ は、

δθ ≡ θ̂ − θ 0 = Cy − Cy 0 = C²(32)となる。ちなみに最小二乗法の適用条件 (式 17) より hδ θ̂i = 0 となるので、パラメータ推定値の期待値も 0 となり、これが不偏推定量になっていることも理解できる。

いまこのパラメータ推定値の分散・共分散からなる誤差行列 Σ ˆ を次のように定義する。
θ(Σ ˆ )jjθ≡ σ ˆ 2 = hhθ̂j − θj 0 i2 iθj(Σ ˆ )jj 0θ≡ h(θ̂j − xj 0 )(θ̂j 0 − xj 0 0 )i (j 6= j 0 , j, j 0 = 1 . . . m)

(j = 1 . . . m)(33)(34)

これらをまとめて表し、32 を代入すると、

Σˆ
θ

≡

˜h(θ̂ − θ 0 )(θ̂ − θ 0 )i

= hC²²̃C̃i = Ch²²̃iC̃ = CΣC̃が得られる。この関係は、測定値からパラメータへの誤差伝播そのものといえる。

5(35)式 35 に 19 と 29 を代入して整理すると (W̃ = W に注意)、

Σˆ
θ

´−1³
´−1
¡
¢
X̃W X
X̃W σ0 2 W −1 W X X̃W X
³
´−1 ³
´−1
= σ0 2 X̃W X= X̃Σ−1 X≡³(36)と書ける。すなわち、パラメータ推定値の誤差は、測定点の取り方とモデル (X) 及び測定精度 (Σ)
から決まり、測定値そのもの (y) には依存しない。

測定誤差の分散がすべて等しい場合

3.1

最小二乗法が適用できるためには 1.1 で与えた 4 つの条件が必要であったが、さらに特殊な場合として、5. 各測定値の分散がすべて等しい (σi = const ≡ σ)。
ことを期待できる場合も多い。この場合、測定数がある程度大きければ、最小二乗解 θ̂ を使って
2 +*mX
σ 2 = yi −Xij θ̂j (37)j=1とみなせる。このとき誤差行列は Σ = σ 2 I (I は単位行列) となり、これを式 35 に代入すればパラメータの推定誤差が得られる。

4ケーススタディ

4.1y = θx

このときモデルパラメータはひとつだけだから n = 1 で θ = θ。いま m 回の測定を行って(xi , yi )(j = 1 . . . m) というデータセットが得られたとすると、ヤコビアン行列 X は
x1 x2 
X= .  .. (38)xmとなる。ここで式 19 のような重み行列を用い、式 30 に代入すれば θ が得られる。
式 18 のように、誤差行列が対角成分しか持たない (つまり測定間の共分散が 0) ならば、式 30 は´³Pmj=1

θ= Pmj=1xj y jσj 2³xj 2σj 2´(39)´(40)となる。またこのパラメータ推定値の誤差は、式 32 を用いて

1³σθ = Pmj=1となる。

6xj 2σj 2

4.2y = θ1 + θ2 x

同様に考え、このときのヤコビアン行列は

1
 1X= ..
 .

x1x2..
.

1 xm

(41)である。
このケースを GNU octave (http://www.octave.org) によって実装した例を示す。以下のスクリプトでは、データセットを標準入力から与えると、θ ベクトルとその誤差行列とが出力される。

Tx を変更すれば、他の場合に適用するのも容易だろう。なお、誤差行列の成分は分散であり、標準偏差ではない点に注意すること。
なお octave には gls, ols という自前の最小二乗フィット関数もあるので、そちらを利用するのも良いだろう。ただしこれらの関数ではパラメータの誤差は評価されない。

#!/usr/bin/octave -qf
line=fgets(stdin);while ( !feof(stdin) )[xi, yi, dyi] = sscanf(line, "%f %f %f", "C");
x(i)
y(i)

= xi;
= yi;

dy(i) = dyi;i++;line=fgets(stdin);endwhile[rx, cx] = size(x);Tx = [ones(rx,1), x];
EM = diag(dy.^2);
WM = inv(EM);C = inv(Tx’ * WM * Tx) * Tx’ * WM;
Theta = C * yvMat = inv(Tx’ * WM * Tx)

4.3

片対数・両対数プロットの場合

測定値 y˜i に測定誤差 δ y˜i がある場合、yi ≡ log10 y˜i の誤差 σi は、誤差伝播を用いて

σi==

∂yi
δ y˜i
∂ y˜i
1 δ y˜iloge 10 y˜i7(42)となる。この yi と σi を用いて、ここまで述べてきた計算を行えばよろしい。絶対誤差 (δ y˜i ) が一定ならば、小さな測定値 y˜i では yi が小さくなるために係数決定に寄与する重みは小さくなる。
一方相対誤差 (δ y˜i /yi ) が一定ならば、それぞれの測定点は同じ重みを持つことがわかるであろう。

参考文献[1] 中川 徹, 小柳義夫「最小二乗法による実験データ解析」東京大学出版会 UP 応用数学選書 7(1982)

8