最小二乗法

データセットを4次関数で最小二乗近似した例

最小二乗法(さいしょうにじょうほう,さいしょうじじょうほう；最小自乗法とも書く、英: least squares method)は、測定で得られた数値の組を、適当なモデルから想定される1次関数、対数曲線など特定の関数を用いて近似するときに、想定する関数が測定値に対してよい近似となるように、残差の二乗和を最小とするような係数を決定する方法、あるいはそのような方法によって近似を行うことである。

前提条件

測定データはモデル関数と誤差の和で表せる。モデル関数が正しければ誤差は偶然誤差だから、測定データとモデル関数の差の分散を最小にすることによって、モデル関数を決定することができる。物理現象の測定データには、誤差が含まれ、それは系統誤差と偶然誤差を含んでいる。この内、偶然誤差は、測定における信号経路の微視的現象に由来するならば、正規分布であると期待されることが多い。また、社会調査などの誤差理由の特定が困難な場合でも誤差が正規分布になると期待する考え方もある。

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

基礎的な考え方

話を簡単にするため、測定値は x, y の二次元の平面に分布するものとし、想定される分布が y = f(x) の形である場合を述べる。想定している関数 f は、既知の関数 g(x) の線型結合で表されていると仮定する。すなわち、

$f(x)=\sum_{k=1}^{m} a_k g_k(x)$

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

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

$(x, y) = (x_1, y_1), (x_2, y_2), \ldots , (x_n, y_n)$

これら (x, y) の分布が、y = f(x) という関数に従うと仮定した時、想定される理論値は (x₁, f(x₁)), (x₂, f(x₂)), ..., (x_n, f(x_n)) ということになり、実際の測定値との残差は、各 i につき |y_i - f(x_i)| ということになる。 この残差の大きさは、xy-平面上での (x_i, y_i) と (x_i, f(x_i)) との距離でもある。

ここで、理論値からの誤差の分散の推定値は残差の平方和

$J = \sum_{i=1}^n (y_i - f(x_i))^2$

で与えられるから、J が最小になるように想定分布 f を(すなわち a_kを)、定めればよいということになる。

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

一次方程式への近似

今、

$(x, y) = (x_1, y_1), (x_2, y_2), \ldots , (x_n, y_n)$

という測定結果が得られたとする。求めたい1次方程式の式を

$y=ax+b\,$

とおくと、aとbは次式で求められる。

$a=\frac{\displaystyle n\sum_{k=1}^n x_ky_k-\sum_{k=1}^n x_k\sum_{k=1}^n y_k}{\displaystyle n\sum_{k=1}^n x^2_k-\left( \sum_{k=1}^n x_k \right)^2}$

$b=\frac{\displaystyle \sum_{k=1}^n x^2_k\sum_{k=1}^n y_k-\sum_{k=1}^n x_ky_k\sum_{k=1}^n x_k}{\displaystyle n\sum_{k=1}^n x^2_k-\left( \sum_{k=1}^n x_k \right)^2}$

解法

当てはめたい関数 f は、

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

と、行列で表すことができる。ここに上付き添字 T は転置行列を表す。すると、最小にする関数 J は

$J=(Ga-y)^T (Ga-y)$

と、表される。ここに、G は、$G_{ij} = g_j(x_i)$ なる成分を持つ行列、$a=(a_1, a_2, \ldots, a_m)^T$ で、$y=(y_1, y_2, \ldots, y_m)^T$。

前節で述べたように J を a のそれぞれの成分で偏微分してゼロと置いた m 個の式（正規方程式）は行列を用いて、

$G^T Ga= G^T y$

と表される。これを正規方程式 (normal equation) と呼ぶ。この正規方程式を解けば係数 a が求まる。

逆行列で正規方程式を解く

行列 G^T G が正則行列（つまりフルランク）である場合は、解 a は一意に求まるため、コレスキー分解を使って解くことができる。ただし、以下の方法のほうが一般に数値的に安定である。

直交分解で正規方程式を解く

上記の方法よりも時間がかかるが、数値的に安定かつ汎用な方法として、QR分解や特異値分解 (SVD) を用いる方法がある。これらの方法では計算の過程で積 G^T G を必要としないため数値的安定性が高い。また G^T G が正則行列でない（ランク落ちしている）場合は正規方程式の解が不定となるが、その場合でも、これらの手法では解 a のうちノルムが最も小さいものを求めることができる。特異値分解を用いる場合は、特異値のうち極めて小さい値を0とみなして計算することで数値計算上の大きな誤差の発生を防ぐことができる (truncated SVD)。

疑似逆行列を使う方法もあるが、効率が悪いため特殊な場合（解析的な数式が必要な場合など）を除いてあまり用いられない。

拡張

多次元

想定される分布が媒介変数 t を用いて (x, y) = (f(t), g(t)) の形(あるいは f, g は複数の媒介変数によって決まるとしても同様)であっても考察される。

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

この場合、各点の理論値 (f(t_i), g(t_i)) と測定値 (x_i, y_i) の間に生じる残差は

$\sqrt{(x_i - f(t_i))^2 + (y_i - g(t_i))^2}$

である。故に、残差平方和は

$\sum_{i=1}^{n}\left\{(x_i - f(t_i))^2 + (y_i -g(t_i))^2\right\}$

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

このように、n 組の (x, y) の測定値 (x_i, y_i) (i = 1, 2, ..., n) を n 組の (x₁, x₂, $\ldots$ x_m) の測定値 (x_1i, x_2i, ..., x_mi) (i = 1, 2, ..., n) に拡張したものも考察することができる。

測定の誤差が既知の場合

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

$J'= \sum_{i=1}^n \frac{(y_i - f(x_i))^2}{\sigma_i^2}$

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

毎回の測定が独立ならば、測定値の尤度は exp(-J') に比例する。そこで、上記の J' を最小にする f は、最尤推定値であるとも解釈できる。また、J' は自由度 n-m のカイ二乗分布に従うので、それを用いてモデル f の妥当性を検定することもできる。

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

非線型最小二乗法

もし、f が、a_k の線型結合で表されないときは、正規方程式を用いた解法は使えず、数値的に a_k の近似値を求める必要がある。例えば、ガウス・ニュートン法やLevenberg-Marquardt法が用いられる。とくにLevenberg-Marquardt法は多くの多次元非線形関数でパラメータを発散させずに効率よく収束させる(探索する)方法として知られている。

異常値の除去

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

誤差が正規分布から極端に外れた異常値を取り除くための方法として修正トンプソン－τ法が用いられる。

関連項目

• 曲線あてはめ