永井 忠一 2024.11.3
二乗和誤差を最小化する\[ J(a, b) = \sum_{i = 0}^{N-1} \left( y_i - (a x_i + b) \right)^2 \rightarrow \min \]係数 a, b を求める
評価関数（誤差関数） J を、係数ベクトル \( \begin{pmatrix} a \\ b \end{pmatrix} \) で微分（偏微分）して 0 になるとして、（連立方程式を）解く\[ \begin{gather} \frac{\partial J(a, b)}{\partial a} = 0, \quad \frac{\partial J(a, b)}{\partial b} = 0 \end{gather} \]
《M 次多項式（一般の場合）として、連立一次方程式を解いた例 → 以前のページ》
（Maxima で記号計算）
データ点数 N = 10 の場合
| Maxima |
|---|
|
解を確認
| Python | |
|---|
解
（最小二乗法の、M 次式の解より）ベクトルと行列で表現\[ \begin{align} A\mathbf w &= \mathbf{b} \\ \begin{pmatrix} \sum\limits_i^{N-1} x_i^{0 + 0} & \sum\limits_i^{N-1} x_i^{0 + 1} \\ \sum\limits_i^{N-1} x_i^{1 + 0} & \sum\limits_i^{N-1} x_i^{1 + 1} \end{pmatrix}\begin{pmatrix} w_0 \\ w_1 \end{pmatrix} &= \begin{pmatrix} \sum\limits_i^{N-1} x_i^0 y_i \\ \sum\limits_i^{N-1} x_i^1 y_i \end{pmatrix} \\ \begin{pmatrix} \sum\limits_i^{N-1} 1 & \sum\limits_i^{N-1} x_i \\ \sum\limits_i^{N-1} x_i & \sum\limits_i^{N-1} x_i^2 \end{pmatrix}\begin{pmatrix} w_0 \\ w_1 \end{pmatrix} &= \begin{pmatrix} \sum\limits_i^{N-1} y_i \\ \sum\limits_i^{N-1} x_i y_i \end{pmatrix} \\ \mathbf w &= A^{-1} \mathbf{b} \\ \begin{pmatrix} b \\ a \end{pmatrix} = \begin{pmatrix} w_0 \\ w_1 \end{pmatrix} &= \begin{pmatrix} \sum\limits_i^{N-1} 1 & \sum\limits_i^{N-1} x_i \\ \sum\limits_i^{N-1} x_i & \sum\limits_i^{N-1} x_i^2 \end{pmatrix}^{-1}\begin{pmatrix} \sum\limits_i^{N-1} y_i \\ \sum\limits_i^{N-1} x_i y_i \end{pmatrix} \end{align} \]
（Maxima による、記号計算）
| Maxima |
|---|
|
（同じテストデータを使用して、）解を表示して確認
| Python | |
|---|
解
\( \Sigma \) が展開されないようにして、解を求める
| Maxima |
|---|
|
（Maxima では、シングルクォート演算子「'」は、シンボルの評価を抑制する）
（同じテストデータを使用して、）検算
| IPython |
|---|
|
解が等しいことを確認（求めた係数 \( a = 0.4719094236134641,\ b = -0.31682418530985723 \)）
（Python のリスト内包表記と zip() 関数を利用して、sum() を計算）
（Octave、R 言語を使用）
| Octave | R 言語 |
|---|---|
|
|
上記の式を、式のとおりに解く\[ \begin{align} A\mathbf w &= \mathbf{b} \\ \begin{pmatrix} \sum\limits_i^{N-1} 1 & \sum\limits_i^{N-1} x_i \\ \sum\limits_i^{N-1} x_i & \sum\limits_i^{N-1} x_i^2 \end{pmatrix}\begin{pmatrix} w_0 \\ w_1 \end{pmatrix} &= \begin{pmatrix} \sum\limits_i^{N-1} y_i \\ \sum\limits_i^{N-1} x_i y_i \end{pmatrix} \\ \begin{pmatrix} w_0 \\ w_1 \end{pmatrix} &= \begin{pmatrix} \sum\limits_i^{N-1} 1 & \sum\limits_i^{N-1} x_i \\ \sum\limits_i^{N-1} x_i & \sum\limits_i^{N-1} x_i^2 \end{pmatrix}^{-1}\begin{pmatrix} \sum\limits_i^{N-1} y_i \\ \sum\limits_i^{N-1} x_i y_i \end{pmatrix} \\ \mathbf w &= A^{-1} \mathbf{b} \end{align} \]
スクリプトの実行結果
| Octave |
|---|
|
| R |
|
（同じテストデータを使用して、）同じ係数が求まっている
Maxima を使用して、記号計算で解く（連立方程式の形で解いているが、逆行列を求めていることと同じ）
| Maxima |
|---|
|
（Maxima で「.」演算子は非可換の乗算、行列の積）
Maxima の fortran() 関数で FORTRAN の式として書き出したものを、C 言語の式に直して（エディタの置換コマンドで「,」を「][」に変換しただけ）、そのまま C/C++ 言語で実装
| C/C++ 言語 |
|---|
|
（同じ係数を計算できている）
... ToDo ...
... ToDo ...
... ToDo ...
掃き出し法、ピボット選択
... ToDo ...
... ToDo ...
... ToDo ...
（天下り）\[ \mathbf t = \Phi \mathbf w \]ここで、\( \mathbf t = \begin{pmatrix} y_0 \\ y_1 \\ \vdots \\ y_{n-1} \end{pmatrix} \)、\( \Phi \) は計画行列
... W.I.P. ...
... W.I.P. ...
... W.I.P. ...
... W.I.P. ...
... W.I.P. ...
... W.I.P. ...
... W.I.P. ...
... W.I.P. ...
... W.I.P. ...
© 2024 Tadakazu Nagai