最小二乗法のソースを表示

[[ファイル:Linear least squares2.png|right|thumb|データセットを4次関数で最小二乗近似した例]]
{{回帰分析}}
'''最小二乗法'''（さいしょうにじょうほう、さいしょうじじょうほう；'''最小自乗法'''とも書く、{{Lang-en-short|least squares method}}）は、誤差を伴う測定値の処理において、その誤差の二乗の和を最小にするようにし、最も確からしい関係式を求める方法である。測定で得られた数値の組を、適当なモデルから想定される[[1次関数]]、[[対数]]曲線など特定の[[関数 (数学)|関数]]を用いて[[近似]]するときに、想定する関数が測定値に対してよい近似となるように、[[残差平方和]]を最小とするような係数を決定する方法<ref name="ut"/><ref name="lh"/><ref name="bj"/>、あるいはそのような方法によって近似を行うことである<ref name="ut"/><ref name="lh">{{cite book |author=Lawson, Charles L.; Hanson, Richard J. |year=1995 |url=https://epubs.siam.org/doi/10.1137/1.9781611971217 |title=Solving Least Squares Problems |publisher=Society for Industrial and Applied Mathematics |doi=10.1137/1.9781611971217}}{{要購読}}</ref><ref name="bj">{{cite book |author=Björck, Åke |year=1996 |url=https://epubs.siam.org/doi/abs/10.1137/1.9781611971484 |title=Numerical Methods for Least Squares Problems |publisher=Society for Industrial and Applied Mathematics |doi=10.1137/1.9781611971484}}</ref>。

== 歴史 ==
[[1805年]]に[[アドリアン＝マリ・ルジャンドル]]が出版したのが初出である。しかし、[[1809年]]に[[カール・フリードリヒ・ガウス]]が出版した際に[[1795年]]には最小二乗法を考案済みだったと主張したことで、最小二乗法の発明者が誰であるかについては不明になっている。

== 計算の概要 ==
=== 前提条件 ===
最小二乗法では測定データ <math>y</math> はモデル関数 <math>f(x)</math> と誤差 <math>\varepsilon</math> の和で
:<math>y=f(x)+\varepsilon</math>
と表せるとする。物理現象の測定データには、[[誤差]]が含まれ、それは[[系統誤差]]と[[誤差#偶然誤差|偶然誤差]]を含んでいる。この内、偶然誤差は、測定における信号経路の[[微視的]]現象に由来するならば、[[正規分布]]であると期待されることが多い。また、社会調査などの誤差理由の特定が困難な場合でも誤差が正規分布になると期待する[[正規分布#正規分布の適用|考え方]]もある。

誤差が正規分布に従わない場合、最小二乗法によって得られたモデル関数はもっともらしくないことに注意する必要がある。偶然誤差が正規分布していない場合、系統誤差が無視できない位大きくそれをモデル関数に含めていない場合、測定データに正規分布から大きく外れた[[外れ値]]を含む場合などが該当する。

上記を含め、最小二乗法の理論的基盤には次のような前提が設けられている<ref name="ut">{{cite|和書 |author=中川徹|author2=小柳義夫 |title=最小二乗法による実験データ解析 |publisher=東京大学出版会 |year=1982 |isbn=4-13-064067-4 |page=30}}</ref><ref name="lh"/><ref name="bj"/>。
* 測定値の誤差には偏りがない。すなわち誤差の[[平均値]]は 0 である。
* 測定値の誤差の[[分散 (確率論)|分散]]は既知である。ただし測定データごとに異なる値でも良い。
* 各測定は互いに[[独立 (確率論)|独立]]であり、誤差の[[共分散]]は 0 である。
* 誤差は正規分布する。
* ''<math>m</math>'' 個のパラメータ<ref group="注釈">''<math>m</math>'' は、測定データの数よりも小さいとする。</ref>（フィッティングパラメータ）を含むモデル関数 <math>f</math> が知られていて、測定量の真の値を近似誤差なく再現することのできるパラメータが存在する。

=== 基礎的な考え方 ===
話を簡単にするため、測定値は <math>x, y</math> の二次元の平面に分布するものとし、想定される分布（モデル関数）が <math>y=f(x)</math> の形である場合を述べる。想定している関数 <math>f</math> は、既知の関数 <math>g(x)</math> の[[線型結合]]で表されていると仮定する。すなわち、

{{Indent|<math>f(x)=\sum_{k=1}^{m} a_k g_k(x)</math>}}

例えば、 <math>g_k(x)=x^{k-1}</math> は、多項式近似であり、特に ''<math>m = 2</math>'' の時は <math>f(x)= a_1 + a_2 x</math> という直線による近似（[[線形回帰]]）になる。<!--図は多項式近似で ''m=5'' の例。-->

今、測定で得られた、次のような数値の組の集合があるとする。

{{Indent|<math>\{(x_1, y_1),\ (x_2, y_2),\ \ldots ,\ (x_n, y_n)\}</math>}}

これら <math>(x, y)</math> の分布が、<math>y = f(x)</math> というモデル関数に従うと仮定した時、想定される理論値は

{{Indent|<math>(x_1, f(x_1)), (x_2, f(x_2)), ..., (x_n, f(x_n))</math>}}

ということになり、実際の測定値との残差は、各 <math>i</math> につき <math>|y_i - f(x_i)|</math> ということになる。

この残差の大きさは、<math>xy</math> 平面上での <math>(x_i, y_i)</math> と <math>(x_i, f(x_i))</math> との距離でもある。

ここで、理論値からの誤差の[[分散 (確率論)|分散]]の推定値は残差の平方和

{{Indent|<math>J = \sum_{i=1}^n (y_i - f(x_i))^2</math>}}

で与えられるから、''<math>J</math>'' が最小になるように想定分布 <math>f</math> を（すなわち <math>a_k</math> を）、定めればよいということになる。

それには、上式は ''<math>a_k</math>'' を変数とする関数と見なすことができるので、''<math>J</math>'' を ''<math>a_k</math>'' について偏微分したものを 0 と置く。こうして得られた ''<math>m</math>'' 個の連立方程式（正規方程式）を解き、''<math>a_k</math>'' を決定すればよい。

=== 一次方程式の場合 ===
さらに簡単な例として、モデル関数を1次関数とし、
:<math>f(x)=ax+b\,</math>
とおくと、''<math>a, b</math>'' は次式で求められる。
:<math>a=\frac{n\textstyle\sum\limits_{k=1}^n x_k y_k -\sum\limits_{k=1}^n x_k \sum\limits_{k=1}^n y_k}{n\textstyle\sum\limits_{k=1}^n {x_k}^2 -\left( \sum\limits_{k=1}^n x_k \right)^2}</math>
:<math>b=\frac{\textstyle\sum\limits_{k=1}^n {x_k}^2 \sum\limits_{k=1}^n y_k -\sum\limits_{k=1}^n x_k y_k \sum\limits_{k=1}^n x_k}{n\textstyle\sum\limits_{k=1}^n {x_k}^2 -\left( \sum\limits_{k=1}^n x_k \right)^2}</math>
:

== 解法例 ==
当てはめたい関数 ''<math>f</math>'' は、
{{Indent|<math>f(x)= (g_1(x), g_2(x), \ldots, g_m(x))(a_1, a_2, \ldots, a_m)^\textrm{T}</math>}}
と表すことができる。上付き添字 T は[[転置行列]]を表す。最小にすべき関数 ''<math>J</math>'' は
:<math>
\begin{align}
J (\boldsymbol{a}) &= (G \boldsymbol{a}-\boldsymbol{y})^\textrm{T} \, (G \boldsymbol{a}-\boldsymbol{y}) \\
&= \left(\begin{bmatrix} G & \boldsymbol{y} \end{bmatrix} \begin{bmatrix}  \boldsymbol{a}\\ -1\end{bmatrix} \right)^\textrm{T} \, \left(\begin{bmatrix} G & \boldsymbol{y} \end{bmatrix} \begin{bmatrix}  \boldsymbol{a}\\ -1\end{bmatrix} \right)
\end{align}
</math>
と表される。ここに ''<math>G</math>'' は、<math>G_{ij} = g_j(x_i)</math> なる成分を持つ ''<math>n \times m</math>'' 行列、<math> \boldsymbol{y}=(y_1, y_2, \ldots, y_n)^\textrm{T}</math>、係数<math>\boldsymbol{a}=(a_1, a_2, \ldots, a_m)^\textrm{T}</math> である。

これの最小解 <math>\boldsymbol{a}</math> は、<math>\begin{bmatrix} G & \boldsymbol{y} \end{bmatrix} ^\textrm{T} \begin{bmatrix} G & \boldsymbol{y} \end{bmatrix}  = \tilde{R}^\textrm{T} \tilde{R}</math> を満たす上三角行列<math>\tilde{R}</math>の計算を経て、解 <math>\boldsymbol{a}</math> を得ることができ、全体の計算量に無駄が少ない。下記の表式を用いると <math>\tilde{R} = \begin{bmatrix} R  & Q^\textrm{T} \boldsymbol{y} \\ \boldsymbol{0}^\textrm{T} & \alpha \end{bmatrix} </math> が得られ、<math>R \boldsymbol{a}= Q^\textrm{T} \boldsymbol{y}</math> から係数解 <math>\boldsymbol{a}</math> を求める<ref group="注釈"><math>\alpha^2 = \min J (\boldsymbol{a})</math>。<math>G^\textrm{T} G</math> は正則行列と仮定。</ref>。<!--この<math>\min J (\boldsymbol{a})</math>の解は、<math>\begin{bmatrix} G & \boldsymbol{y} \end{bmatrix} ^\textrm{T} \begin{bmatrix} G & \boldsymbol{y} \end{bmatrix}  = \tilde{R}^\textrm{T} \tilde{R}</math>を満たす正方行列<math>\tilde{R}</math>を得れば、<math>J (\boldsymbol{a}) = \begin{bmatrix}  \boldsymbol{a}\\ -1\end{bmatrix}^\textrm{T}  \tilde{R}^\textrm{T} \tilde{R} \begin{bmatrix}  \boldsymbol{a}\\ -1\end{bmatrix} </math>を最小化して得られる。<math>\tilde{R}</math>を上[[三角行列]]とすれば計算量に無駄がなく<ref><math>\begin{bmatrix} G & \boldsymbol{y} \end{bmatrix} ^\textrm{T} \begin{bmatrix} G & \boldsymbol{y} \end{bmatrix}  = \tilde{R}_1^\textrm{T} \tilde{R}_2</math>を満たすLU分解で上[[三角行列]]を得ても良く、その<math>\tilde{R}_1</math>を使っても<math>\tilde{R}_2</math>を使っても、係数解<math>\boldsymbol{a}</math>を計算できる。</ref>、-->

また前節で述べたように ''J'' を<math>\boldsymbol{a}</math>のそれぞれの成分で偏微分して 0 と置いた ''<math>m</math>'' 個の式（正規方程式）は行列を用いて、

{{Indent|<math>G^\textrm{T} G \boldsymbol{a}= G^\textrm{T} \boldsymbol{y}</math>}}

と表される。これを'''正規方程式'''（{{Lang-en-short|normal equation}}）と呼ぶ。この正規方程式を解けば係数解 <math>\boldsymbol{a}</math> が求まる。

係数解 <math>\boldsymbol{a}</math> の解法には以下のようないくつかの方法がある。

* 逆行列で正規方程式を解く
*: 行列 <math>G^TG</math> が[[正則行列]]（つまりフルランク）である場合は、解 <math>\boldsymbol{a} = (G^\textrm{T} G)^{-1} G^\textrm{T} \boldsymbol{y}</math> は一意に求まる。ただし <math>G^TG</math> の逆行列を明示的に求めることは通常は良い方法ではない。
:: 計算量が小さい方法として[[コレスキー分解]]<ref name="Yamamoto1"/>（<math>G^\textrm{T} G = R^\textrm{T} R</math>、<math>R</math> は ''<math>m \times m</math>'' 上[[三角行列]]）による三角行列分解を経て、最終的に <math>R \boldsymbol{a}= R^{-\textrm{T}} G^\textrm{T} \boldsymbol{y}</math> を解けばよい。
:: [[数値的安定性]]確保のためには、積 <math>G^TG</math> を経ない三角行列分解が良い。すなわち以下と同じく[[QR分解]]（直交分解）による <math>G = Q R</math> から、<math>R \boldsymbol{a}= Q^\textrm{T} \boldsymbol{y}</math> を解く。
<!--:: ''R''およびこの右辺は、まとめて三角行列分解を行い<math>\begin{bmatrix} G & \boldsymbol{y} \end{bmatrix} = Q \, \begin{bmatrix} R & Q^\textrm{T} \boldsymbol{y} \end{bmatrix} </math>と得ることもできる。-->

* 直交分解で正規方程式を解く
*: コレスキー分解の方法よりも計算量が大きいが、数値的に安定かつ汎用な方法として、[[QR分解]]や[[特異値分解]]（SVD）を用いる方法がある<ref name="Yamamoto1">{{Cite book |和書 |author=山本哲朗 |title=数値解析入門 |edition=増訂版 |date=2003-06 |publisher=[[サイエンス社]] |series=サイエンスライブラリ 現代数学への入門 14 |ISBN=4-7819-1038-6}}</ref>。これらの方法では計算の過程で積 <math>G^TG</math> を必要としないため[[数値的安定性]]が高い。また <math>G^TG</math> が正則行列でない（ランク落ちしている）場合は正規方程式の解が不定となるが、その場合でも、これらの手法では解 ''a'' のうち[[ノルム]]が最も小さいものを求めることができる。特異値分解を用いる場合は、特異値のうち極めて小さい値を 0 とみなして計算することで数値計算上の大きな誤差の発生を防ぐことができる({{Lang-en-short|truncated SVD}})<ref>{{Cite journal |author=Hansen, Per Christian |year=1987 |url=https://doi.org/10.1007/BF01937276 |title=The truncated SVD as a method for regularization |journal=BIT Numerical Mathematics |volume=27 |pages=534-553 |publisher=Springer |doi=10.1007/BF01937276}} {{要購読}}</ref>。
* [[擬似逆行列]]を使う方法もあるが<ref>安川章. (2017). 科学実験/画像変換の近似計算に便利な 「疑似逆行列」 入門 できる人が使っている最小二乗法の一発フィット. インターフェース= Interface, 43(8), 142-146.</ref>、計算効率が悪いため、特殊な場合（解析的な数式が必要な場合など）を除いてあまり用いられない。

== 拡張 ==
=== 多次元 ===
想定される分布が[[媒介変数]] ''t'' を用いて <math>(x, y) = (f(t), g(t))</math> の形（あるいは <math>f, g</math> は複数の媒介変数によって決まるとしても同様）であっても考察される。

すなわち、測定値 <math>(x_i, y_i)</math> がパラメータ <math>t_i</math> に対する <math>(f(t_i), g(t_i))</math> を理論値として近似されているものと考えるのである。

この場合、各点の理論値 <math>(f(t_i), g(t_i))</math> と測定値 <math>(x_i, y_i)</math> の間に生じる残差は
{{Indent|<math>\sqrt{(x_i - f(t_i))^2 + (y_i - g(t_i))^2}</math>}}
である。故に、[[残差平方和]]は

{{Indent|<math>\sum_{i=1}^{n}\left\{(x_i - f(t_i))^2 + (y_i -g(t_i))^2\right\}</math>}}

となるから、この値が最小であるように、<math>f, g</math> を決定するのである。

このように、''<math>n</math>'' 組の <math>(x, y)</math> の測定値 <math>(x_i, y_i)(i = 1, 2, ..., n)</math> を''<math>n</math>'' 組の <math>(x_1, x_2, ..., x_m)</math> の測定値 <math>(x_{1i}, x_{2i}, ..., x_{mi})(i = 1, 2, ..., n)</math> に拡張したものも考察することができる。

=== 測定の誤差が既知の場合 ===
''<math>n</math>'' 回の測定における[[誤差]]があらかじめ分かっている場合を考える。異なる測定方法で測定した複数のデータ列を結合する場合などでは、測定ごとに誤差が異なることはしばしばある。誤差が正規分布していると考え、その[[標準偏差]] <math> \sigma_i (i=1,2,\ldots,n )</math> で、誤差の大きさを表す。すると、誤差が大きい測定より、誤差が小さい測定の結果により重みをつけて近似関数を与えるべきであるから、

{{Indent|<math>J'= \sum_{i=1}^n \frac{(y_i - f(x_i))^2}{\sigma_i^2}</math>}}

を、最小にするように ''f'' を定める方がより正確な近似を与える。

毎回の測定が[[独立 (確率論)|独立]]ならば、測定値の[[尤度]]は ''<math>exp(-J')</math>'' に比例する。そこで、上記の ''<math>J'</math>'' を最小にする ''<math>f</math>'' は、[[最尤法|最尤推定値]]であるとも解釈できる。また、''<math>J'</math>'' は[[自由度]] ''<math>n-m</math>'' の[[カイ二乗分布]]<ref>Weisstein, Eric W. "Chi-Squared Distribution." From MathWorld--A Wolfram Web Resource. {{URL|https://mathworld.wolfram.com/Chi-SquaredDistribution.html}}</ref>に従うので、それを用いてモデル ''<math>f</math>'' の妥当性を検定することもできる。

毎回の測定誤差が同じ場合、''<math>J'</math>'' を最小にするのは ''<math>J</math>'' を最小にするのと同じ意味になる。

=== 非線形最小二乗法 ===
{{main|非線形最小二乗法}}
もし、''<math>f</math>'' が、''<math>a_k</math>'' の線型結合で表されないときは、正規方程式を用いた解法は使えず、[[反復解法]]を用いて数値的に ''<math>a_k</math>'' の近似値を求める必要がある。例えば、[[ガウス・ニュートン法]]<ref>{{cite book |author=MAGREÑÁN, Alberto; Argyros, Ioannis |year=2018 |title=A contemporary study of iterative methods: convergence, dynamics and applications |publisher=Academic Press |url=https://books.google.co.jp/books?hl=ja&lr=lang_ja|lang_en&id=62U6DwAAQBAJ&oi=fnd&pg=PP1&dq=Magrenan,+A.+A.,+%26+Argyros,+I.+(2018).+A+contemporary+study+of+iterative+methods:+convergence,+dynamics+and+applications.+Academic+Press.&ots=gozfhHu53F&sig=oEUmo3nfc5tWdGWQguzrEVkDQM4#v=onepage&q=Magrenan%2C%20A.%20A.%2C%20%26%20Argyros%2C%20I.%20(2018).%20A%20contemporary%20study%20of%20iterative%20methods%3A%20convergence%2C%20dynamics%20and%20applications.%20Academic%20Press.&f=false}}</ref>や[[レーベンバーグ・マーカート法]]<ref>
Weisstein, Eric W. "Levenberg-Marquardt Method." From MathWorld--A Wolfram Web Resource. {{URL|https://mathworld.wolfram.com/Levenberg-MarquardtMethod.html}}</ref><ref name="more">{{Cite journal |author=Moré Jorge J. |year=1978 |url=https://link.springer.com/chapter/10.1007/BFb0067700 |title=The Levenberg-Marquardt algorithm : Implementation and theory |journal=Lecture Notes in Mathematics |volume=630 |pages=105-116 |doi=10.1007/BFb0067700}} {{要購読}}</ref><ref name="yu">Yu, H., & Wilamowski, B. M. (2011). Levenberg-marquardt training. Industrial electronics handbook, 5(12), 1.</ref><ref name="rang">{{Cite journal |author=Ranganathan, Ananth |year=2004 |url=https://sites.cs.ucsb.edu/~yfwang/courses/cs290i_mvg/pdf/LMA.pdf |format=PDF |title=The levenberg-marquardt algorithm |journal=Tutoral on LM algorithm |volume=11 |issue=1 |pages=101-110}}</ref><ref>{{Cite journal|和書|author=山下信雄, 福島雅夫 |date=2000-10 |url=https://hdl.handle.net/2433/64462 |title=Levenberg-Marquardt法の局所収束性について (最適化の数理科学) |journal=数理解析研究所講究録 |ISSN=1880-2818 |publisher=京都大学数理解析研究所 |volume=1174 |pages=161-168 |hdl=2433/64462 |CRID=1050001201691367552}}</ref>が用いられる。とくにレーベンバーグ・マーカート法（{{Lang-en-short|Levenberg-Marquardt Method}}）は多くの多次元非線形関数でパラメータを発散させずに効率よく収束させる（探索する）方法として知られている<ref name="more"/><ref name="yu"/><ref name="rang"/>。
<!--以下の説明はパラメータaについて非線形になっていないのでコメントアウトします。
あてはめの曲線の関数値''y''<sub>''n''</sub>を求める。

{{Indent|<math>\tilde{y}_n = a_0 + a_1 x_n + a_2 x_n^2 + \cdots a_N x_n^N, \qquad n = 1,~2,~\cdots N, \qquad N \geqq n + 1</math>}}

これは、ベクトルを利用して以下のように表せる。

{{Indent|<math>\tilde{y}_n = \boldsymbol{x}_n \cdot \boldsymbol{a}</math>}}

ただし、

{{Indent|<math>\boldsymbol{x}_n = (1,~x_n,~x_n^2,\cdots x_n^N)</math>}}
{{Indent|<math>\boldsymbol{a} = (a_0,~a_1,~a_2,\cdots a_N)^\textrm{T}</math>}}

二乗誤差の総和をεとすると、

{{Indent|
<math>
\epsilon = \sum_{n=1}^{N} (y_n - \tilde{y}_n)^2 = \sum_n (y_n - \boldsymbol{x}_n \cdot \boldsymbol{a}) = (\boldsymbol{y} - \hat{X} \cdot \boldsymbol{a})
</math>
}}

ただし、

{{Indent|<math>\hat{X} = (\boldsymbol{x}_1,~\boldsymbol{x}_2,\cdots \boldsymbol{x}_N)^\textrm{T}</math>}}
{{Indent|<math>\boldsymbol{y} = (y_1,~y_2,\cdots y_N)^\textrm{T}</math>}}

ここで、εの最小値を計算するために、微分して0とおく。

{{Indent|<math>
\frac{\partial \epsilon}{\partial \boldsymbol{a}} = -2 \hat{X} \cdot (\boldsymbol{y} - \hat{X} \cdot \boldsymbol{a}) = 0
</math>}}
{{Indent|<math>
\boldsymbol{a} = (\hat{X}^\textrm{T} \hat{X})^{-1} \hat{X} \boldsymbol{y}
</math>}}

これを計算することで、当てはめ曲線多項式の定数を定めることができる。
-->

=== 異常値の除去 ===
前提条件の節で述べたように、測定データを最小二乗法によって近似する場合、外れ値または異常値が含まれていると極端に近似の尤もらしさが低下することがある。また、様々な要因によって異常値を含む測定はしばしば得られるものである。

誤差が正規分布から極端に外れた異常値を取り除くための方法として修正トンプソン-τ法が用いられる<ref>Michele Rienzner (2020). Find Outliers with Thompson Tau ({{URL|https://www.mathworks.com/matlabcentral/fileexchange/27553-find-outliers-with-thompson-tau}}), [[MATLAB]] Central File Exchange. Retrieved May 17, 2020.</ref>。

== 関連項目 ==
{{Commonscat|Least squares}}
{{Wikibooks|最小二乗法}}
* {{仮リンク|総最小二乗法|en|Total least squares}} (<ref>Van Huffel, S., & Vandewalle, J. (1991). The total least squares problem: computational aspects and analysis (Vol. 9). SIAM.</ref><ref>{{Cite journal |author=Golub, Gene H; Van Loan, Charles F |year=1980 |url=https://doi.org/10.1137/0717073 |title=An analysis of the total least squares problem |journal=SIAM journal on numerical analysis |volume=17 |issue=6 |pages=883-893 |publisher=SIAM |doi=10.1137/0717073}}</ref>)<!--Perpendicular offsets?-->
* [[ガウス＝マルコフの定理]]<ref>Drygas, H. (2012). The coordinate-free approach to Gauss-Markov estimation (Vol. 40). Springer Science & Business Media.</ref>
* [[曲線あてはめ]]<ref>Motulsky, H., & Christopoulos, A. (2004). Fitting models to biological data using linear and nonlinear regression: a practical guide to curve fitting. Oxford University Press.</ref>

== 脚注 ==
{{脚注ヘルプ}}

=== 注釈 ===
<references group="注釈" responsive="" />

=== 出典 ===
<references responsive="" />

{{authority control}}

{{DEFAULTSORT:さいしようししようほう}}
[[Category:曲線あてはめ]]
[[Category:最適化]]
[[Category:回帰分析]]
[[Category:数学に関する記事]]
[[Category:最小二乗法|*]]