最小二乗法の計算法

POINT

  • 最小二乗法の計算を解説.

最小二乗法の計算について紹介します.微分法による極値問題の一例としても良い題材です.

【注】この記事は(まだ?)細かい部分を詰めることはせず,ざっくりとした計算の流れを整理することを主眼としています.

Excel関数を用いる(実用的な)方法は以下を参照してください:

1次関数(直線)

$N$組のデータ点$(x_1,y_1),...,(x_N,y_N)$を最小二乗法で直線フィッティングすることを考えます.つまり,
\begin{align}
S(a,b)=\sum_{i=1}^N\bigl[y_i - (\underset{\text{直線の式}}{\underline{ax_i+b}} )\bigr]^2
\end{align}を最小にするような$a,b$の値を求める問題を考えます.



【方針】
求める$a,b$は,$\displaystyle \frac{\partial S}{\partial a}(a,b)=0$, $\displaystyle \frac{\partial S}{\partial b}(a,b)=0$を満たします.したがって,この2式を連立して解くことで,$a,b$を計算できます.


【解】
偏微分は,合成関数の微分法を使うことで
\begin{align}
\frac{\partial S}{\partial a}(a,b)
&=\sum_{i=1}^N 2\left[y_i - (ax_i+b)\right]\cdot (-x_i) \\
\frac{\partial S}{\partial b}(a,b)
&=\sum_{i=1}^N 2\left[y_i - (ax_i+b)\right]
\end{align}と計算できます(補足:*1).したがって,$\displaystyle \frac{\partial S}{\partial a}(a,b)=0$, $\displaystyle \frac{\partial S}{\partial b}(a,b)=0$はそれぞれ
\begin{align}
\frac{\partial S}{\partial a}(a,b)&=0:\\
&\sum_{i=1}^N \left(x_i \right)^2 \cdot a +\sum_{i=1}^N x_i \cdot b = \sum_{i=1}^N x_i y_i\\
\frac{\partial S}{\partial b}(a,b)&=0: \\
&\sum_{i=1}^N x_i\cdot a+Nb=\sum_{i=1}^N y_i
\end{align}となります.


ここで,それぞれ$\{x_i\}_i$, $\{y_i\}_i$を当確率で値にとる確率変数$X$, $Y$を導入しましょう(*2).すると

  • $X$の平均値:$\displaystyle\langle X\rangle=\frac{1}{N}\sum_{i=1}^N x_i$
  • $Y$の平均値:$\displaystyle\langle Y\rangle=\frac{1}{N}\sum_{i=1}^N y_i$
  • $X$の分散:
    \begin{align}
    \sigma_X^2
    &=\biggl\langle \bigl(X-\langle X\rangle\bigr)^2
    \biggr\rangle \\
    &=\langle X^2\rangle
    -\langle X\rangle^2
    \end{align}
  • $X,Y$の共分散:
    \begin{align}
    \mathrm{Cov}\left(\mathrm{X,Y}\right)
    &=\biggl\langle \bigl(X-\langle X\rangle\bigr)
    \bigl(Y-\langle Y\rangle\bigr)
    \biggr\rangle \\
    &=\langle XY\rangle
    -\langle X\rangle\langle Y\rangle
    \end{align}
という量を考えられます.これらを使うと,得られた方程式は
\begin{align}
\frac{\partial S}{\partial a}(a,b)&=0:\\
& \left(\sigma_X^2 +\langle X\rangle^2 \right) a +\langle X\rangle b \\
&\qquad =\mathrm{Cov}\left(\mathrm{X,Y}\right) + \langle X\rangle\langle Y\rangle \\
\frac{\partial S}{\partial b}(a,b)&=0: \\
& \langle X\rangle a+b=\langle Y\rangle
\end{align}と書き直すことができます.


2式を連立させて解くことで
\begin{align}
a&=\frac{\mathrm{Cov}\left(\mathrm{X,Y}\right)}{\sigma_X^2} \\
b&=-a\langle X\rangle+\langle Y\rangle \\
&=\frac{\langle X^2\rangle \langle Y\rangle - \langle XY\rangle \langle X\rangle}{\sigma_X^2}
\end{align}が得られます.//

【注意】
求めた直線は
\begin{align}
y = a(x-\langle X\rangle) + \langle Y\rangle
\end{align}という形をしています.したがって,この直線は平均値からなる点$(\langle X\rangle,\langle Y\rangle)$を通ります.

多変数の線形回帰

(上の場合を含むような)一般化した問題として,複数のパラメータについて線形な関数によるフィッティングを考えます.

フィッティング関数として
\begin{align}
f(x^1,…,x^m)
=\sum_{i=1}^m \theta_m x^m
\end{align}を考えます(ここで,$x^i$は冪乗ではなく「単なる添字」を意味します).この式は,ベクトル$\boldsymbol{\theta}={}^t(\theta_1,...,\theta_m)$, $\boldsymbol{x}={}^t(x^1,...,x^m)$を用いて
\begin{align}
f(\boldsymbol{x})
=\boldsymbol{x} \cdot \boldsymbol{\theta}
\end{align}と表すことができます.

それでは$N$組のデータ点$(\boldsymbol{x}_1,y_1),…,(\boldsymbol{x}_N,y_N)$のフィッティングを考えましょう.$N$組の式
\begin{align}
y_i - \boldsymbol{x}_i \cdot \boldsymbol{\theta}
\end{align}は行列$A=(a_{ij})=x^j_{\:i}$($x^j_{\:i}$は$\boldsymbol{x}_i$の$j$成分)を用いて
\begin{align}
\boldsymbol{y} - A \boldsymbol{\theta}
\end{align}と書き直せます.

したがって,二乗和は$\left\| \boldsymbol{y} - A \boldsymbol{\theta}\right\|^2$$={}^t(\boldsymbol{y} - A \boldsymbol{\theta})(\boldsymbol{y} - A \boldsymbol{\theta})$となります.ここでは,(データに相関がある場合にも扱えるよう)重み付きの二乗和
\begin{align}
S(\boldsymbol{\theta})=
{}^t(\boldsymbol{y} - A \boldsymbol{\theta}) W (\boldsymbol{y} - A \boldsymbol{\theta})
\end{align}を最小化することを考えます(*3).

1変数の場合と同じように$\displaystyle\frac{\partial S}{\partial \theta_i}(\boldsymbol{\theta})=0$$\,(i=1,...,m)$をみたす$\boldsymbol{\theta}$が求めるものです.これは簡単に計算することができ(*4
\begin{align}
\frac{\partial S}{\partial \theta_i}(\boldsymbol{\theta})&=0: \\
&-2\:{}^tA W (\boldsymbol{y} - A \boldsymbol{\theta})=0
\end{align}となります.したがって
\begin{align}
\boldsymbol{\theta}=\left({}^tA WA \right)^{-1} \left( {}^tA W\boldsymbol{y} \right)
\end{align}が求める値です.

参考文献

データ・誤差解析の基礎

データ・誤差解析の基礎

*1:$g(x)=x^2$, $h_i(a,b)=y_i - (ax_i+b)$により,$S(a,b)=\sum_i g\left(h_i(a,b)\right)$と表されるので,合成関数の微分法より \begin{align} &\frac{\partial S}{\partial a}(a,b) \\ &=\sum_{i=1}^N \frac{\mathrm{d} g}{\mathrm{d} x}\left(h_i(a,b)\right) \frac{\partial h_i}{\partial a}(a,b) \\ &=\sum_{i=1}^N 2\left[y_i - (ax_i+b)\right]\cdot (-x_i) \end{align}と計算できます.$\partial S/\partial b$についても同様です.

*2:あるいは, \begin{align} & \langle X\rangle =\frac{1}{N}\sum_{i=1}^N x_i \\ & \Delta x_i =x_i - \langle X\rangle \\ & \langle f(\Delta X,\Delta Y) \rangle =\frac{1}{N}\sum_{i=1}^N f(\Delta x_i, \Delta y_i) \end{align}という量を「定義」したと考えても良いです($f$は関数).この場合,最終的な結果は \begin{align} \sigma_X^2 &\rightarrow \langle (\Delta X)^2\rangle \\ \mathrm{Cov}\left(\mathrm{X,Y}\right) &\rightarrow \langle \Delta X \Delta Y\rangle \\ \end{align}と置き換えたものになります.

*3:$W=(w_{ij})$は重みを表す対称行列です.ベクトル$\boldsymbol{a}$に対して \begin{align} {}^t\boldsymbol{a} W\boldsymbol{a} =\sum_{i,j} a_i w_{ij} a_j \end{align}となります.

*4:関数$f(\boldsymbol{\theta})=[{}^t\boldsymbol{a}(\boldsymbol{\theta})] W [\boldsymbol{b}(\boldsymbol{\theta})]$を考えましょう.ベクトル・行列の成分で表せば \begin{align} f(\boldsymbol{\theta})=\sum_{jk}a_j w_{jk}b_k \end{align}なので,積の微分法より \begin{align} \frac{\partial f}{\partial \theta_i} &=\sum_{jk}w_{jk}\left( \frac{\partial a_j}{\partial \theta_i}\ b_k +a_j \frac{\partial b_k}{\partial \theta_i} \right) \\ &={}^t \left( \partial \boldsymbol{a} / \partial \theta_i \right) W \boldsymbol{b} +{}^t \boldsymbol{a} W \left( \partial \boldsymbol{b} / \partial \theta_i \right) \end{align}となります.特に,$\boldsymbol{a}=\boldsymbol{b}$の場合には \begin{align} \frac{\partial f}{\partial \theta_i} &={}^t \left( \partial \boldsymbol{a} / \partial \theta_i \right) (W + {}^t W)\boldsymbol{a} \end{align}と変形することができます(スカラー量は転置をとっても値を変えないので,${}^t\boldsymbol{a} W\boldsymbol{b}$$={}^t\left({}^t\boldsymbol{a} W\boldsymbol{b}\right)$$={}^t\boldsymbol{b} {}^tW\boldsymbol{a}$.).