最小二乗法

目的

最小二乗法[Ordinary Least Squares;OLS]の目的は複雑な構造を持つデータを少数の式の重ね合わせ[線形結合,linear combination]で表現することにあります.このように大量のデータを少数のデータで表現することを音声や画像に応用する場合はデータ圧縮[data compression]と呼ばれます.最小二乗法はデータ圧縮の一つともいえます.

最良線形予測

$y$ と $x$ という2つのデータがあって,いま,$x$ の値が分かっている場合に,$y$ の値を推測(予測)するということを考えてみます.$x$ と $y$ というデータは次のように表されます.\[\left(\begin{array}{cc}x_{1}\ y_{1}\\x_{2}\ y_{2}\\\vdots \ \vdots\\x_{N}\ y_{N}\end{array}\right)\]$x$ が与えられたとき,$y$ の条件付き期待値 $E(y|x)$ は,$y$ の平均2乗予測誤差を最小にする $x$ の関数となります.

一般的にいうと,$E(y|x)$ は $x$ の非線形関数になります.しかし,ここでは簡単に次のような線形関数を考えることとします.\[f(x)=\beta_{0}+\beta_{1}x\]そうすると,平均2乗予測誤差は,\[E\{(y-b_{0}-b_{1}x)^{2}\}\]と表されます.

この平均最小2乗予測誤差を最小化します.

ここで新たに,\[y=\beta_{1}+\beta_{2}x+u\]となる $u$ という変数を定義します.

この式を,平均2乗予測誤差の\[E\{(y-\beta_{0}-\beta_{1}x)^{2}\}\]に代入します.

\[E(u)=0,E(xu)=0\]という関係を利用すると,\[\begin{align}E\{(y-\beta_{0}-\beta_{1}x)^{2}\}&=E\{u+(\beta_{1}-b_{1})+(\beta_{2}-b_{2})x\}^{2}]\\&=E(u^{2})+[\beta_{1}-b_{1}\ \beta_{2}-b_{2}]\left[\begin{array}{cc}1\ E(x)\\E(x)\ E(x^{2})\end{array}\right]\left[\begin{array}{c}\beta_{1}-b_{1}\\\beta_{2}-b_{2}\end{array}\right]\end{align}\]となります.

さらに,\[x=[1 \ x]'\]とおくと,\[\left[\begin{array}{cc}1\ E(x)\\E(x)\ E(x^{2})\end{array}\right]=E[xx']\]となります.

つまり,\[\begin{align}E\{(y-\beta_{0}-\beta_{1}x)^{2}\}&=E\{u+(\beta_{1}-b_{1})+(\beta_{2}-b_{2})x\}^{2}]\\&=E(u^{2})+[\beta_{1}-b_{1}\ \beta_{2}-b_{2}]\left[\begin{array}{cc}1\ E(x)\\E(x)\ E(x^{2})\end{array}\right]\left[\begin{array}{c}\beta_{1}-b_{1}\\\beta_{2}-b_{2}\end{array}\right]\\&=E(u^{2})+[\beta_{1}-b_{1}\ \beta_{2}-b_{2}] E(xx') \left[\begin{array}{c}\beta_{1}-b_{1}\\\beta_{2}-b_{2}\end{array}\right]\end{align}\]となります.

ここで,任意の2次元ベクトル $a$ に対しては,\[a'xx'a=(a'x) \geq 0\]となり,\[a'E(xx')a=E\{(a'x)^{2}\}\]となるので,\[E(xx')\]は非負定符号行列[nonnegative definite matrix]であり,固有値は全て非負になります.

また,$E(xx')$ の行列式,\[|E(xx')|=E(x^{2})-\{E(x)\}^{2}=var(x)\]となるため,$var(x) > 0$ ならば,$E(xx')$ は正則ということになります.

以上の結果から,\[E\{(y-\beta_{0}-\beta_{1}x)^{2}\}=E(u^{2})+[\beta_{1}-b_{1}\ \beta_{2}-b_{2}] E(xx') \left[\begin{array}{c}\beta_{1}-b_{1}\\\beta_{2}-b_{2}\end{array}\right]\]の右辺第2項は,\[b_{1}=\beta_{1},b_{2}=\beta_{2}\]のときを除いて常に正となります.

よって,$\beta_{1},\beta_{2}$ が最小解であることが分かりました.

この解は最良線形予測[best linear prediction]と呼ばれますが,このように,条件付き期待値とは一般には一致しません.

最小2乗推定量(OLS estimation)

\[S(\beta_{1},\beta_{2})=E\{(y-\beta_{1}-\beta_{2}x)^{2}\}\]を $\beta_{1},\beta_{2}$ で偏微分すると,$S(\beta_{1},\beta_{2})$ が最小化されるためには,\[\frac{\partial S(\beta_{1},\beta_{2})}{\partial \beta_{1}}=0\]\[\frac{\partial S(\beta_{1},\beta_{2})}{\partial \beta_{2}}=0\]が成り立っている必要があります.これは,最小化の1階条件[first-order condition]と言われます.

最小化の1階条件の式を展開すると,合成関数の微分から,\[\frac{\partial S(\beta_{1},\beta_{2})}{\partial \beta_{1}}=2(y-\beta_{1}-\beta_{2}x)=0\]\[\frac{\partial S(\beta_{1},\beta_{2})}{\partial \beta_{1}}=2(y-\beta_{1}-\beta_{2}x)x=0\]となります.

両辺を $2$ で割ると,\[\begin{array}y-\beta_{1}-\beta_{2}x&=0\\(y-\beta_{1}-\beta_{2}x)x&=0\end{array}\]

ここで,先ほどと同じように,\[X=[1 \ x]'\]とおくと,\[\begin{array}(y-\beta_{1}-\beta_{2}x)&=0\\(y-\beta_{1}-\beta_{2}x)x&=0\end{array}\]という2つの式は1つにまとめることが出来て,\[X'(y-X\beta)=0\]と表すことが出来ます.この関係式を正規方程式(normal equation)と言われます.

正規方程式[normal equation]を変形すると,\[X'y-X'X\beta=0\]となり,さらに変形すると,\[X'X\beta=X'y\]となり,$X$ の階数が列数に等しいとすると,$X'X$ は正則であって逆行列が存在します.

従って,\[\beta=(X'X)^{-1}X'y\]となります.

この $\beta$ を最小2乗推定量(OLS estimation)といいます.

ここで,$x$ の平均は $E(x)$ なので,\[(y-\beta_{1}-\beta_{2}x)\]という式は,\[(y-\beta_{1}-\beta_{2}x+\beta_{2}E(x)-\beta_{2}E(x))\]と変形出来ます.

この式を整理すると,\[(y-\beta_{1}-\beta_{2}(x-E(x))-\beta_{2}E(x))\]となります.

この式の期待値をとると,\[E(y-\beta_{1}-\beta_{2}(x-E(x))-\beta_{2}E(x))\]となりますので,さらに変形して,\[E(y)-E(\beta_{1})-E(\beta_{2}x)+E(\beta_{2}E(x))-E(\beta_{2}E(x))\]となり,\[E(y)-\beta_{1}-\beta_{2}E(x)+\beta_{2}E(x)-\beta_{2}E(x)\]つまり,\[E(y)-\beta_{1}-\beta_{2}E(x)\]となるので,\[E(y)-\beta_{1}-\beta_{2}E(x)=0\]から,\[\beta_{1}=E(y)-\beta_{2}E(x)\]となります.

続いて,\[y-\beta_{1}-\beta_{2}x=0\]であることから,\[(x-E(x))(y-\beta_{1}-\beta_{2}x)=0\]となるので,ここに,\[\beta_{1}=E(y)-\beta_{2}E(x)\]を代入すると,\[(x-E(x))(y-(E(y)-\beta_{2}E(x))-\beta_{2}x)=0\]つまり,\[(x-E(x))(y-E(y)+\beta_{2}E(x)-\beta_{2}x)=0\]となるので,さらに,\[(x-E(x))(y-E(y)+\beta_{2}(E(x)-x))=0\]そして,\[(x-E(x))(y-E(y))+(x-E(x))(\beta_{2}(E(x)-x))=0\]となるので,\[(x-E(x))(y-E(y))-\beta_{2}(E(x)-x)^{2}=0\]すなわち,\[\beta_{2}=\frac{(x-E(x))(y-E(y))}{(E(x)-x)^{2}}\]となるので,結局,$\beta_{1},\beta_{2}$ は,\[\begin{align}\beta_{1}&=E(y)-\beta_{2}E(x)\\\beta_{2}&=\frac{E[\{x-E(x)\}\{y-E(y)\}]}{E[\{x-E(x)\}^{2}]}\end{align}\]となります.

Pythonでの例


# coding: utf-8
from __future__ import division
#Python 2.7でimportをすればPython 3風に,Python 2.7では割り算「/」の結果は切り捨てでしたが
#Python 3以降ではちゃんと小数点以下まで保持されるようになります

from matplotlib import pyplot as plt

import numpy as np


def main():
# データ点数は 100 にする
N = 100
# 0 ~ 30 (-15 ~ +15) の範囲で,ランダムに値を散らばらせる
R = 30

#N個のデータを生成して x に代入
x = np.arange(N)

# y = 50 + x の直線に (-5 ~ +5) でランダム性をもたせる
#0x@Ρの一様乱数をN個生成する
y = np.arange(N) + np.random.rand(N) * R - R // 2 + 50

# 平均値と標準偏差
xmean, xsigma2 = x.mean(), x.var()
# 平均値
ymean = y.mean()

# 共分散
sxy = np.sum(x * y) / N - xmean * ymean

# 回帰係数ベータ
beta = sxy / xsigma2
# 回帰係数アルファ
alpha = ymean - beta * xmean

# 得られた回帰直線の一次式
print('y = {0} + {1}x'.format(alpha, beta))

# グラフに回帰直線をひく
ols = np.array([alpha + beta * xi for xi in x])
plt.plot(x, ols, color='r')

# 元データもプロットする
plt.scatter(x, y)

# グラフを表示する
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()


if __name__ == '__main__':
main()

Vita brevis, ars longa. Omnia vincit Amor.





















確率測度の拡張 情報量 世界最古の都市 Tipologia Edilizia(都市類型学) 層としての都市 都市構造の普遍性