永井 忠一 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=0}^{N-1} x_i^{0 + 0} & \sum\limits_{i=0}^{N-1} x_i^{0 + 1} \\ \sum\limits_{i=0}^{N-1} x_i^{1 + 0} & \sum\limits_{i=0}^{N-1} x_i^{1 + 1} \end{pmatrix}\begin{pmatrix} w_0 \\ w_1 \end{pmatrix} &= \begin{pmatrix} \sum\limits_{i=0}^{N-1} x_i^0 y_i \\ \sum\limits_{i=0}^{N-1} x_i^1 y_i \end{pmatrix} \\ \begin{pmatrix} \sum\limits_{i=0}^{N-1} 1 & \sum\limits_{i=0}^{N-1} x_i \\ \sum\limits_{i=0}^{N-1} x_i & \sum\limits_{i=0}^{N-1} x_i^2 \end{pmatrix}\begin{pmatrix} w_0 \\ w_1 \end{pmatrix} &= \begin{pmatrix} \sum\limits_{i=0}^{N-1} y_i \\ \sum\limits_{i=0}^{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=0}^{N-1} 1 & \sum\limits_{i=0}^{N-1} x_i \\ \sum\limits_{i=0}^{N-1} x_i & \sum\limits_{i=0}^{N-1} x_i^2 \end{pmatrix}^{-1}\begin{pmatrix} \sum\limits_{i=0}^{N-1} y_i \\ \sum\limits_{i=0}^{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=0}^{N-1} 1 & \sum\limits_{i=0}^{N-1} x_i \\ \sum\limits_{i=0}^{N-1} x_i & \sum\limits_{i=0}^{N-1} x_i^2 \end{pmatrix}\begin{pmatrix} w_0 \\ w_1 \end{pmatrix} &= \begin{pmatrix} \sum\limits_{i=0}^{N-1} y_i \\ \sum\limits_{i=0}^{N-1} x_i y_i \end{pmatrix} \\ \begin{pmatrix} w_0 \\ w_1 \end{pmatrix} &= \begin{pmatrix} \sum\limits_{i=0}^{N-1} 1 & \sum\limits_{i=0}^{N-1} x_i \\ \sum\limits_{i=0}^{N-1} x_i & \sum\limits_{i=0}^{N-1} x_i^2 \end{pmatrix}^{-1}\begin{pmatrix} \sum\limits_{i=0}^{N-1} y_i \\ \sum\limits_{i=0}^{N-1} x_i y_i \end{pmatrix} \\ \mathbf w &= A^{-1} \mathbf{b} \end{align} \]
（上記の式は０オリジン。Octave、R言語は１オリジンなので注意）
スクリプトの実行結果
| Octave |
|---|
|
| R |
|
（同じテストデータを使用して、）同じ係数が求まっている
Maxima を使用して、記号計算で解く（連立方程式の形で解いているが、逆行列を求めていることと同じ）
| Maxima |
|---|
|
（Maxima で「.」演算子は行列の積）
Maxima の fortran() 関数で FORTRAN の式として書き出したものを、C言語の式に直して（エディタの置換コマンドで「,」を「][」に変換しただけ）、そのまま C/C++ 言語で実装
| C/C++言語 |
|---|
|
（同じ係数を計算できている）
| C/C++言語 |
|---|
|
手計算と同じように解いている
ピボット選択を pivoting() 関数で行っている。行の交換だけを行うものを部分選択、行と列も両方入れ替える方法を完全選択という（ここでは、partial pivoting を採用）
自前の row_reduction() 関数を、教科書を見ながら書いた関数「gauss_elimination()」に入れ替えて、最小二乗法を解く
（上記 row_reduction() 関数は、手計算のときと同じようなアニメーションを作成するために、本当は省略できる計算や処理もしているので、無駄がある）
| C/C++言語 |
|---|
|
（同じ係数を、計算できている）
（LU分解や、QR分解について）《積み残し》
（天下り）\[ \mathbf t = \Phi \mathbf w \]ここで、\( \mathbf t = \begin{pmatrix} y_0 \\ y_1 \\ \vdots \\ y_{n-1} \end{pmatrix} \)、\( \Phi \) は計画行列
（以前のページ → w の最尤推定解）《省略》
（以前のページ → gnuplot fit コマンド）
（以前のページ → R言語の lm() などの関数を使用）
多項式フィッティング
| Octave | Python NumPy |
|---|---|
|
|
（結果の係数は、降べきの順）
© 2024 Tadakazu Nagai