最小二乗法

HOME>>メモ>>アルゴリズム>>最小二乗法

最小二乗法

最小二乗法(最小自乗法とも書く)は、実験データのような誤差を含む値から、 それにフィットするような関数を求める手法の一つです。

例えば等速直線運動運動の観測実験を行なった結果、 次のような実験データが得られたとします。

実験結果
時間(s)位置(m)
0.01.0
1.02.3
2.03.4
3.04.4
4.05.9
5.07.0

この実験結果を時間をx軸、位置をy軸にとってグラフにプロットすると次のようになります。

実験結果
この運動は等速直線運動なので、すべての点を通る直線が一本引けるはずです。 このような直線を引ければ、速度や初期位置がわかるので、いろいろと便利です。

しかし現実にはどうしても測定誤差が発生してしまうため、 どんなに頑張って直線を引いても、図のようにデータポイントと直線は少しずれてしまいます。 仕方がないので、できるだけデータポイントから外れないような直線を引くことにしましょう。 最小二乗法はこのような直線を引く代表的な手法の一つです。

線形回帰

「できるだけデータポイントから外れない」とはどういうことでしょう？ もう少し数学的に表してみることにします。

i番目のデータポイントを$(x_i, y_i)$、求めたい直線を$y=ax+b$と 定義します。データポイントの数はnとします。

誤差
すると、実際の測定値と予想される理論値との誤差を$y_i-(ax_i+b)$と表すことができます。 この誤差を小さくすればデータポイントと直線が近くなりそうです。

一つのデータポイントにだけ近づけても仕方がないのでこの誤差の総和を取ります。 しかし、ただ単純に足すだけだとプラスとマイナスが打ち消しあってしまいます。 絶対値を取ってから足すのもいいのですが、絶対値というのは数学的に少し扱いが難しいので、 二乗して総和をとることにしましょう。

\[J=\sum_{i=1}^n (y_i-(ax_i+b))^2\]

このJを最小化することが最小二乗法の目標です。

Jを最小にするには「極値となる点では、微分した結果が0になる」ということを 利用します。Jは下に凸な二次関数なのでこれで最小値が求まるはずです。 今具体的に求めたいのはaとbの値なので、それぞれについて偏微分してみましょう。

\begin{eqnarray*} \frac{\partial J}{\partial a} &=&-2\sum_{i=1}^nx_i(y_i-(ax_i+b))\\ \frac{\partial J}{\partial b} &=&-2\sum_{i=1}^n(y_i-(ax_i+b)) \end{eqnarray*}

それぞれをイコール0と置いて整理すると、aとbに関する連立方程式となります。

\begin{eqnarray*} a\sum_{i=1}^nx^2 + b\sum_{i=1}^nx_i &=&\sum_{i=1}^nx_iy_i\\ a\sum_{i=1}^nx^2 + bn &=&\sum_{i=1}^ny_i\\ \end{eqnarray*}

あとはこれを消去法なり代入法なりクラメルの公式なりを使って解くだけです。 結論だけ書くと次のようになります。

\begin{eqnarray*} a&=&\frac{\displaystyle n\sum_{i=1}^nx_iy_i-\sum_{i=1}^nx_i\sum_{i=1}^ny_i} {\displaystyle n\sum_{i=1}^nx_i^2-\left(\sum_{i=1}^nx_i\right)^2}\\ b&=&\frac{\displaystyle\sum_{i=1}^nx_i^2\sum_{i=1}^nx_iy_i-\sum_{i=1}^nx_i\sum_{i=1}^nx_iy_i} {\displaystyle n\sum_{i=1}^nx_i^2-\left(\sum_{i=1}^nx_i\right)^2}\\ \end{eqnarray*}

もっと複雑な関数で近似

もっと複雑な関数で近似してみることを考えてみましょう。 求めたい関数$f(x)$がm個の既知の関数$g_k(x)$の 線形結合で表されているとします。

\[ f(x) = \sum_{k=1}^m a_kg_k(x) \]

ここで、$a_k$は関数を決定するためのパラメータです。 この$a_k$を求めることが最終目標です。

$f(x)$ですが、行列を使って次のように書くことができます。

$f(x)=(g_1(x),g_2(x),\cdots,g_m(x))(a_1,a_2,\cdots,a_m)^T$

これを縦にn個だけ並べれば、$x=x_i$における関数$f(x)$の値を 簡単に表すことが出来ます。

\begin{eqnarray*} \left(\begin{array}{c} f(x_1)\\f(x_2)\\\vdots\\f(x_n) \end{array}\right) &=& \left(\begin{array}{cccc} g_1(x_1) & g_2(x_1) & \cdots & g_m(x_1) \\ g_1(x_2) & g_2(x_2) & \cdots & g_m(x_2) \\ \vdots & \vdots & \ddots & \vdots \\ g_1(x_n) & g_2(x_n) & \cdots & g_m(x_n) \\ \end{array}\right) \left(\begin{array}{c} a_1\\a_2\\\vdots\\a_m \end{array}\right) \\ &=&\bm{G}\bm{a} \end{eqnarray*}

ここで、 $\bm{a}=(a_1,a_2,\cdots,a_m)^T$ \[\bm{G}=\left(\begin{array}{cccc} g_1(x_1) & g_2(x_1) & \cdots & g_m(x_1) \\ g_1(x_2) & g_2(x_2) & \cdots & g_m(x_2) \\ \vdots & \vdots & \ddots & \vdots \\ g_1(x_n) & g_2(x_n) & \cdots & g_m(x_n) \\ \end{array}\right)\] と置きました。

実際に測定した値を$\bm{y}=(y_1,y_2,\cdots,y_n)^T$とすると 誤差の総和Jは次のように表せます。

\begin{eqnarray*} J&=&(\bm{G}\bm{a}-\bm{y})^T(\bm{G}\bm{a}-\bm{y})\\ &=&\bm{a}^T\bm{G}^T\bm{G}\bm{a}-2\bm{a}^T\bm{G}^T\bm{y}+\bm{y}^T\bm{y} \end{eqnarray*}

これを最小化するためにパラメータaについて偏微分し、0とと置きます。

\[\frac{\partial J}{\partial\bm{a}}=2\bm{G}^T\bm{G}\bm{a}-2\bm{G}^T\bm{y}=\bm{0}\]

整理すると、解かなければならない連立方程式は \[ \bm{G}^T\bm{G}\bm{a}=\bm{G}^T\bm{y} \] となります。

実際にやってみる

SVGを使ってグラフを描くJavascriptを作ってみました。 黒い線が元の関数(直線, 2次関数, 三角関数から選択可)、 赤い点が元関数にランダムノイズを加えたもの、 青い線が最小二乗法によって求めた近似曲線です。

SVGの描画にはRaphaëlを 行列計算にSylvesterを使用しています。

元関数:
誤差: 測定点の個数:
近似関数: