- 追加された行はこの色です。
- 削除された行はこの色です。
- 最小二乗法 へ行く。
#freeze
*最少二乗法：The least-squares method [#nbc81473]
偏差(ズレ 隔たり)の二乗和を最小にするように、推定式の係数（パラメータ）を決める方法を最小二乗法と呼ぶ。
偏差(ズレ 隔たり)の尺度として、二乗和を採ることを提案したのは、ルジャンドルであった。1805年の論文「天体の軌道決定の新しい方法」において、「係数は、誤差の大きさを決めるが、誤差二乗和を最小となるように決定する必要がある。この方法が一般的に言って簡単な方法であることが分かった。」と述べている。
これには、深い根拠がある。
誤差二乗和を最小にすると定式化することで、最適化問題は、正規方程式と呼ばれる係数に関する連立一次方程式を解いて、その解を最適な係数とすることに帰着される。
偉大なガウスは、最小二乗法の確率論的基礎を導いている。最小二乗法は、測定誤差の正規分布と深く関係しており、その根拠を最尤法を用いてしめした。ガウスは、「被説明変数の推定誤差の二乗和は、期待値がゼロ、分散が最小になるように係数が決められねばならない」と考えた。そして、最小二乗法で推定したパラメータは、不偏性をもつこと、言い換えればパラメータの期待値が真のパラメータに一致することを明らかにしている。
"Least squares" means that the overall solution minimizes the sum of the squares of the errors made in solving every single equation.
The most important application is in data fitting. The best fit in the least-squares sense minimizes the sum of squared residuals, a residual being the difference between an observed value and the value provided by a model.
The least-squares method was first described by Carl Friedrich Gauss around 1794.
*歴史 [#sc1159d3]
-ラプラスが「最小二乗法」を発表（１７９９）
-ガウスが小惑星Ceresの軌道予測（１８０１）
--子午線上の測地データをもとに地球の楕円度を計算．（１７９９）
--9月にガウスが軌道を計算,予測,12月7日に予測通りに再発見。その後最小二乗法を用いて軌道を精密に計算。
-ルジャンドル,最小二乗法を発表（１８０５）
-ガウス,最小二乗法の原理を説明（１８０９）
--1795年に発見したと主張。「なぜ二乗か？」をガウス分布を用いて説明。ガウス,最小二乗法の論文を発表（１８２３）。ガウス分布に限らず一般の誤差分布に対する最適性を示す。
*問題の定式化 [#ce8a77b5]
A simple data set consists of n points (data pairs) (xi,yi), i = 1, ..., n, where xi is an independent variable and yi is a dependent variable whose value is found by observation. The model function has the form f(x,β), where the m adjustable parameters are held in the vector β. The goal is to find the parameter values for the model which "best" fits the data. The least squares method finds its optimum when the sum, S, of squared residuals
S = ∑ ri^2 i = 1, ..., n
is a minimum. A residual is defined as the difference between the value of the dependent variable and the model value
ri = yi − f(xi,β).
An example of a model is that of the straight line. Denoting the intercept as β0 and the slope as β1, the model function is given by f(x,β)=β0+β1・x .
*線形最少二乗法：Linear least squares [#y2d2ec72]
See [[Wikipedia>http://en.wikipedia.org/wiki/Linear_least_squares]].
*１次方程式の場合 [#u70deb22]
いま、観測してデータが(x1,y1),(x2,y2),・・・,(xn,Yn)を得た。
これを求めたい1次方程式(モデル)を
y=ax+b+e ,eは測定誤差の確率変数。
とする。
このとき誤差は
ri=yi-(axi+b) i = 1, ..., n
で表わされる。誤差の平方和を最少にするように、a,bの係数を求める。
S=∑ ri^2 --->Min(a,b)
sub. to ri=yi-(axi+b) i = 1, ..., n
極値なので、Sをそれぞれの係数aとbで、微分すれば0となる必要がある。
Σ2ri・xi=0
Σ2ri=0
ri=yi-(axi+b)
これを正規方程式とよぶ。この解が求めるパラメータ(a.b)である。
--標本誤差の総和はゼロであり、標本誤差と説明変数は直交していることが残差平方和最小化の条件となっていることが理解できる。(Σri・xi=0)
-aとbは次式で表わされる。
#ref(http://upload.wikimedia.org/math/3/0/9/309c7bbb62647d11dc9b285910d8b987.png)
#ref(http://upload.wikimedia.org/math/c/c/8/cc8b10aa58d755331990f5030a9b78a7.png)
これを、を使って表わすと
a=S(x,y)/S(x,x)
b=E(y)-(S(x,y)/S(x,x))・E(x)
但し、S(X,Y)=∑(xi-E(x))・(yi-E(y))/n ,E(x)=∑xi/n
*相関係数Rとは [#he503037]
上記の1次方程式の場合、xiに対するyの推定値y*i=aXi+bと観測されたyiについて
次式を定義する。
R^2=∑(y*i-E(y))^2 / ∑(yi-E(y))^2
これに、y*i=aXi+bを代入して整理すれば
R^2=S(x,y)^2/(s(x,x)・S(y,y))
すなわち
R=S(x,y)/(√s(x,x)・√S(y,y))
相関係数はxとyの共分散をxの標準偏差、yの標準偏差の積で割った値となる。
-R > 0 のとき x と y は 正の相関を持つ。傾きが正。
-R < 0 のとき x と y は 負の相関を持つ。傾きが負。
-R =1 のとき x と y は y=ax+b,a>0のモデル式で完全に説明できる。1次従属で推定誤差なし。
--正の強い相関がある場合の図
#ref(soukan.JPG)
*相関係数と1次回帰式 [#c037a2d4]
いま、観測してデータが(x1,y1),(x2,y2),・・・,(xn,Yn)を得た。
これを求めたい1次方程式(モデル)をy=ax+bとする。
この回帰係数は、先に示した。書き直すと
y=a・（x-E(x))＋E(y)
a=R・σy/σx
-性質１：傾き a
a=S(x,y)/S(x,x)= Σ(xi - E(x))(yi-E(y))/Σ(xi-E(x))2 = R・σy/σx
すなわち、傾き a は、相関係数とσy/σxの積で表わされる。
Rは、主として傾きの程度を規格化したものと考えられ、y方向とx方向の標準偏差の比に比例して傾きが決まる。
-性質２：X=E(x)のとき y=E(y) になる。点（E(x),E(y))を通る直線
*決定係数とは：決定係数＝相関係数の2乗 [#ib7965a3]
yのデータをy=ax+bで予測する時、次式が成り立つ。
Σ(yi-E(y))^2 = Σ(yi-y*i)^2 ＋ Σ(y*i-E(y))^2
これは
予測データの全変動＝予測誤差の全変動＋予測された全変動
と見ることができる。
--何故か？。
yi-E(y)=yi-y*i + y*i-E(y)
であるので、両辺の2乗をとると
(yi-E(y))^2=(yi-y*i)^2 + (y*i-E(y))^2+2・(yi-y*i)(y*i-E(y))
上式をi=1~nで合計したものが決定係数である。第3項の2乗和は零になる。
第3項の2乗和=2Σ[ri・(y*i-E(y))]^2=2Σ[ri・R・σy/σx（xi-E(x))]
=2R・σy/σx・Σ[rixi-E(x))]
--最初の正規方程式の直交条件Σrixi=0から、第3項は零になる。別の言い方では、回帰による偏差 (y - y*)は，残差(y* -E(y)) と無相関（独立)になるように求められるからである。
''決定係数とは、全変動の内、何割が予測されたかを表わす指標である。''
R^2=Σ(y*i-E(y))^2 / Σ(yi-E(y))^2
総変動のうち回帰式で説明できる変動の割合を表わす。
-証明
--予測式は y*i=(R・σy/σx)（xi-E(x))＋E(y)なので
決定係数 =[Σ(R・σy/σx)（xi-E(x))]^2/σy^2
=(R・σy/σx)^2・σx^2/σy^2
=R^2
*決定係数、相関係数、誤差二乗和の関係式 [#rf00b049]
決定係数=(相関係数)^2 = 1 - (yの誤差二乗和)/(yの分散：全変動)
--誤差二乗和を最小にすることは、相関係数が最も大きくなるように係数を決めることと等価である。
--全変動の内、何割が予測されたかを表わす指標である。1から標準化した予測誤差二乗和を引いた値である。
*推定値a,bの性質：期待値と分散 [#p938fd84]
パラメータ推定値をa*,b*とするとき、どちらもデータの関数であるので、観測データによって値が異なる。
本当のパラメータa,bに近い値が、データを十分に多くすれば得られるだろうか？。
そこで、残差項のみが確率変数である次のモデルを考える。
y=ax+b+e
（１）説明変数Xは確率変数ではない。
（２）誤差eの平均はゼロである。（すべてのiについて、E(ei)=0）
（３）誤差eの分散は均一である。（すべてのiについて、V(ei)=σ2（一定））
（４）誤差eに、自己相関がない。（すべてのi,jについて、Cov(ei,ej)=0）
このとき、パラメータ推定値の性質はどのようであろうか。
-不偏性
''推定値a*はａの不偏推定量である。E(a*)=a''
--証明:yi=axi+b+riを代入してみる。
E(a*)=E[(Σxiyi/n-E(x)E(y))/(Σ(xi-E(x))^2/n)
=E[(Σxi(axi+b+ri)/n-E(x)E(axi+b+ri))/(Σ(xi-E(x))^2/n)
={aE[Σxi^2/n-E(x)E(x)]+bE[Σxi/n-E(x)]+E[Σxiri/n-E(x)E(ri)}/(Σ(xi-E(x))^2/n)
=aE[Σxi^2/n-E(x)E(x)]/(Σ(xi-E(x))^2/n)
ゆえに E(a*)=a
同様に、''推定値b*はbの不偏推定量である。E(b*)=b 。''
-分散と共分散
--推定値a*の分散は
E[(a*-a)^2]=σ^2/(Σ(xi-E(x))^2/n)=σ^2/(E(x^2)-E(x)^2)
ただしσ^2はモデルの誤差の分散
--推定値b*の分散は
E[(b*-b)^2]=σ^2・E(x^2))/(E(x^2)-E(x)^2)
--共分散は、
Cov(a*,b*)=E((a*-E(a*))(b*-E(b*))
=E((a*-a)(b*-b))
=σ^2E(X)/(E(x^2)-E(x)^2)
となる。
''これらの分散の式を使って、パラメータの信頼性の検定ができる。''
*多変数最小二乗法 [#a6a194cd]
二つの標本XとYの間に従属（回帰）関係があって、
y=a0+a1xi+a2x2+・・・・+amxm +e 但し、eは誤差を表わす確率変数。
な多変数の関係が成り立つと考える。
標本について、
yi=a0+Σajxij+ei （総和は、j=1,･･･,m）
であることとなる。説明変数xは非確率変数で多重共線性はなく、誤差eに対するガウス・マルコフの仮定は満たされているとしよう。n個の標本のこの関係式を行列で表して、
Y=Xa+e
但し、Yは観測値のn次元列ベクトル、eは誤差の列ベクトル、Xは要素｛xij｝をもつm+1行n列のマトリックスでaはm次元のパラメータベクトル。ちなみにXのn行目は、Xn=（1,xn1,･･･,xnm）である。
すると誤差項の2乗の和（残差平方和）は、転地記号tを用いて
E^2=et・e
=(Y-Xa)t・(Y-Xa)
=Yt・Y-2at・Xt・Y+at・Xt・Xa
となるので、これの最小化条件は、aで偏微分し、ゼロとおく。
∂E2/∂a=-2Xt・Y+2Xt・Xa=0
故に、
(Xt・X)a=XtY
という正規方程式を得る。
もし、Xt・Xのm次元行列が正則ならば、逆行列が存在する。
未知パラメータaの推定値は、正規方程式の解である次式で表わされる。
a = [Xt・X]-1・XtY
1次式の場合に確認した残差平方和の最小化の条件式、Σeixi=0について確認する。
これは直交条件とも呼ばれていた。
Xte=Xt(Y-Xa)
=XtY-XtXa
=0
すなわち、正規方程式が直交条件そのものを表わしていることがわかる。
-推定パラメータの期待値と分散
続けて推定されたパラメータの期待値と分散を求めておこう。変数名が紛らわしくなるが、真の値のベクトルをａ、推定値をa*とすれば、推定値の期待値は、
E(a*)=E([Xt・X]-1・XtY)
=E([Xt・X]-1・Xt(Xa＋e))
=[Xt・X]-1(Xt・X)E(a)
=a
このように推定値の期待値が真の値になることを不偏推定量という。
推定値の分散は、スカラーcを伴う分散公式V(cα)=c2V(α)とV(c+α)=V(α)を使って、
V(a*)=V((XtX)-1XtY)
=(XtX)-1XtV(Y)((XtX)-1Xt)t
=(XtX)-1XtV(Xa+e)((XtX)-1Xt)t
=(XtX)-1Xtσ2I((XtX)-1Xt)t
=(XtX)-1σ2
となる。解となる分散の式は列ベクトルでなく行列となる。
うまくデータを取れば、分散を小さくできるであろうが、その情報はXtXにある。
そこで''XtXなるマリックスを情報行列と呼ぶ''。
*直交射影と最小二乗法 [#x86471b7]
--線形空間での射影
y観測値のベクトルは所与である。立体空間にある点Yから観測データXの列ベクトルが張る空間へ射影を考える。ベクタ y が X の列空間にある場合とは，等価的に X の列ベクタの線型和で書けるので，y=xaが厳密に成立します．そうでない場合は， y から 「Xの列ベクタで書けない成分」を抜き取ります．これが残差であります。空間上のデータをy*=xaで変換してyの空間上でyとy*の間にはy-y*なるズレが認められるとする。そこでe=y-y*なる残差のベクトルを考える。
この残差のベクトルをXの列ベクトルが張る空間に直交するように、言い換えればaを変化させて、eの長さ（ノルム）が最小になるようにyを射影することで、最適なaが求められます。
この射影では、ベクトルaは変数であるからいかほどでも伸び縮みすることができるので、Y-Xaという差分のベクトルのノルムが最小となるようにaを調整し定めるのが、最適化である。
最適化の結果、y-y*=eなる誤差ベクトルとデータXの張る空間は直交している。
最小二乗法では、Yの終点からXaへの最短距離が垂線であるから、YのXaへの射影を求めれば、Y-Xaというベクトルのノルムを最小とするベクトルとなる。
--ピタゴラスの定理は直角三角形の定理であるが
｜Y｜^2=｜Xa｜^2+｜Y-Xa｜^2
が、成立するためにはaは正規方程式を満たす必要があるとも言える。
#ref(projection.JPG)
*参考 [#u7e02430]
-[[一般化最小誤差法に基づく解法>http://www.ism.ac.jp/editsec/toukei/pdf/53-2-331.pdf]]
-[[離散システムのための逐次推定法>http://dickie.ees.kyushu-u.ac.jp/HP/wadalab/wadalab2004/undergd/miyakawa.pdf]]
-[[Exelおける回帰分析（最小二乗法）の手順と出力.pdf>http://keijisaito.info/pdf/excel_ols.pdf]]
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
