最小二乗法 | physics-asks

01合う直線をどう決めるか

点と直線の縦の差 $y_i - (a + b\,t_i)$ が、いちばん小さくなる直線を探したい。でも差は点の数だけあります。どう 1 つの数にまとめるか。

差を足す

+ と − が打ち消し合い、合っていない直線でも 0 になります。平均の点 $(\bar t, \bar y)$ を通る直線なら、どれでも 0 です。

点の平均の高さに引いた横の直線と、各点との差。上の差に +、下の差に − が付き、和は 0

絶対値を足す

打ち消し合いませんが、 $|x|$ は角があって微分できず、「微分が 0 の所が底」で探せません。答えが 1 本に決まらないこともあります。

2 乗を足す

負にならず、角もありません。差を 1 辺とする正方形の面積の和です。直線を動かしてみてください。

切片 a 傾き b

2 乗の和

S = \sum_{i=1}^{n} \bigl(y_i - f(t_i)\bigr)^2

これをいちばん小さくするのが最小二乗法です。

$t$ は正確で $y$ にだけ誤差があり、どの点の誤差も同じくらいで互いに関係がない、とします(誤差が点ごとに違うなら、差をその誤差で割ってから 2 乗)。誤差が正規分布なら、これがいちばんありそうな値(最尤推定)になります。

02直線を計算で求める

01 節の図のボタン「2 乗の和がいちばん小さい直線」が出した直線を、計算で求めます。 $S$ は $a, b$ の 2 次式で、 $a$ と $b$ を軸にとると、谷が 1 つのお椀の形です(異なる $t_i$ が 2 つ以上あるとき)。

01 節の 7 点に直線を合わせるときの S を、切片 a と傾き b の平面に等高線で描いた図。楕円が 1 つの点を囲む — 01 節の 7 点の S の等高線。谷の底(点)が、ボタンが出す直線

谷の底では、 $a$ で微分しても $b$ で微分しても 0 です(多変数の極値)。

\frac{\partial S}{\partial a} = -2\sum_i (y_i - a - b\,t_i) = 0,\qquad \frac{\partial S}{\partial b} = -2\sum_i t_i\,(y_i - a - b\,t_i) = 0

\begin{gathered}\frac{\partial S}{\partial a} = -2\sum_i (y_i - a - b\,t_i) = 0,\\[4pt]\frac{\partial S}{\partial b} = -2\sum_i t_i\,(y_i - a - b\,t_i) = 0\end{gathered}

これは $a, b$ の連立 1 次方程式で、解くと( $\bar t, \bar y$ は平均)

近似直線の係数

b = \frac{\sum_i (t_i - \bar t)(y_i - \bar y)}{\sum_i (t_i - \bar t)^2},\qquad a = \bar y - b\,\bar t

\begin{gathered}b = \frac{\sum_i (t_i - \bar t)(y_i - \bar y)}{\sum_i (t_i - \bar t)^2},\\[4pt]a = \bar y - b\,\bar t\end{gathered}

f(t) = a + b\,t

\bar t = \sum_i t_i/n

\bar y = \sum_i y_i/n

例題 14 点 (0, 0.9), (1, 3.1), (2, 4.9), (3, 7.1) に直線を合わせる。

\bar t = 1.5,\qquad \bar y = 4.0

\begin{gathered}\bar t = 1.5,\\[4pt]\bar y = 4.0\end{gathered}

平均との差 $(t_i - \bar t,\ y_i - \bar y)$ は $(-1.5, -3.1)$ 、 $(-0.5, -0.9)$ 、 $(0.5, 0.9)$ 、 $(1.5, 3.1)$ なので

\begin{aligned} &\textstyle\sum (t_i - \bar t)(y_i - \bar y) \\ &= 4.65 + 0.45 + 0.45 + 4.65 = 10.2 \end{aligned}

\textstyle\sum (t_i - \bar t)^2 = 2.25 + 0.25 + 0.25 + 2.25 = 5

b = \frac{10.2}{5} = 2.04

a = 4.0 - 2.04 \times 1.5 = 0.94

直線は $y = 0.94 + 2.04\,t$ 。残差は $-0.04,\ 0.12,\ -0.12,\ 0.04$ です。

$a + b\,t + c\,t^2$ のように $t$ について曲がっていても、 $a, b, c$ について 1 次なら同じように解けます( $S$ がやはり 2 次式になるため)。この場合は、異なる $t_i$ が 3 つ以上あれば答えはただ 1 つです。

演習 1

3 点 $(t, y) = (0, 1.0),\ (1, 2.9),\ (2, 5.1)$ に直線 $y = a + bt$ を合わせよう。

答えを見る

$\bar t = 1,\ \bar y = 3.0$ 。

平均との差は $(t_i - \bar t,\ y_i - \bar y) = (-1, -2.0),\ (0, -0.1),\ (1, 2.1)$ なので

$b = \frac{2.0 + 0 + 2.1}{1 + 0 + 1} = 2.05$

$a = 3.0 - 2.05 \times 1 = 0.95$

残差は $0.05,\ -0.10,\ 0.05$ で、足すと 0 になる( $\partial S/\partial a = 0$ の式そのもの)。

03直線に直せる式

合わせたい式が直線とは限りません。たとえば減っていく量 $y = a\,e^{-bt}$ は曲線です。でも対数をとると $\ln y = \ln a - b\,t$ で、 $t$ と $\ln y$ の直線です。02 節の式がそのまま使えます。

例題 23 点 (0, 5.0), (1, 2.9), (2, 1.9) に $y = a\,e^{-bt}$ を合わせる。

$\ln y$ は $1.6094,\ 1.0647,\ 0.6419$ (平均 $1.1053$ )。

$\ln y$ の平均との差は $0.5041,\ -0.0406,\ -0.4635$ 、 $t_i - \bar t$ は $-1, 0, 1$ なので

-b = \frac{-0.5041 + 0 - 0.4635}{1 + 0 + 1} = -0.4838

\ln a = 1.1053 + 0.4838 \times 1 = 1.5891

よって $b = 0.484,\ a = e^{1.5891} = 4.90$ 。 $\ln y$ の残差は $0.020,\ -0.041,\ 0.020$ です。

対数をとるので、 $y$ がすべて正のときだけ使えます。また、 $y$ の誤差 $\sigma$ が一定でも、 $\ln y$ の誤差は( $\sigma \ll y$ のとき)およそ $\sigma / y$ で、 $y$ が小さい点ほど大きくなります。きちんと合わせるなら、差を $\sigma / y_i$ で割ってから 2 乗します。

04直線に直せない式

対数でも直線にならない式、たとえば $y = \cos(\omega t)$ はどうするか。 $\omega$ について 1 次でないので一度では解けず、 $S$ の谷がいくつもできることもあります。

cos(ωt) を合わせるときの S を ω に対して描いた図。深い谷が 1 つと浅い谷がいくつもある — 深い谷(点)の外から始めると、近くの浅い谷で止まる

$a, b$ や $\omega$ のように、式の中で求めたい数をパラメータと呼び、まとめて $p$ と書きます。今の値の近くで式を 1 次で近似します。すると、ずらす量 $\Delta p$ は 02 節と同じように解けます。

f(t_i;\,p + \Delta p) \approx f(t_i;\,p) + \sum_j J_{ij}\,\Delta p_j,\qquad J_{ij} = \frac{\partial f(t_i;\,p)}{\partial p_j}

\begin{gathered}f(t_i;\,p + \Delta p) \approx f(t_i;\,p) + \sum_j J_{ij}\,\Delta p_j,\\[4pt]J_{ij} = \frac{\partial f(t_i;\,p)}{\partial p_j}\end{gathered}

J^{\mathsf T}J\,\Delta p = J^{\mathsf T} r,\qquad r_i = y_i - f(t_i;\,p)

\begin{gathered}J^{\mathsf T}J\,\Delta p = J^{\mathsf T} r,\\[4pt]r_i = y_i - f(t_i;\,p)\end{gathered}

$J$ は偏微分を並べた表(行列)、 $J^{\mathsf T}$ はその行と列を入れかえたもの、 $r$ は今の残差です。

これをくり返すのが Gauss–Newton 法です(scipy.optimize.curve_fit などはその改良版を使う)。反復の初期値はグラフから大まかに読んで、深い谷に置きます。たとえば振動なら、山の間隔から $\omega$ を読みます(減衰振動のデータ解析)。

05近似直線の誤差

点に誤差があるので、近似直線の傾きと切片にも誤差があります。「 $b = 2.04 \pm 0.06$ 」のように書くためのものです。どの点も同じ大きさ $\sigma$ の、互いに関係のない誤差をもつとき

\sigma_b = \frac{\sigma}{\sqrt{\sum_i (t_i - \bar t)^2}},\qquad \sigma_a = \sigma\sqrt{\frac{1}{n} + \frac{\bar t^{\,2}}{\sum_i (t_i - \bar t)^2}}

\begin{gathered}\sigma_b = \frac{\sigma}{\sqrt{\sum_i (t_i - \bar t)^2}},\\[4pt]\sigma_a = \sigma\sqrt{\frac{1}{n} + \frac{\bar t^{\,2}}{\sum_i (t_i - \bar t)^2}}\end{gathered}

$\sigma$ が分からないときは残差から見積もります( $m$ はパラメータの個数。直線なら 2)。

\sigma^2 \approx \frac{S_{\min}}{n - m}

例題 3例題 1 の直線の、傾きと切片の誤差を求める。

残差 $-0.04,\ 0.12,\ -0.12,\ 0.04$ から

S_{\min} = 0.032,\qquad \sigma^2 \approx \frac{0.032}{4 - 2} = 0.016,\qquad \sigma \approx 0.1265

\begin{gathered}S_{\min} = 0.032,\\[4pt]\sigma^2 \approx \frac{0.032}{4 - 2} = 0.016,\\[4pt]\sigma \approx 0.1265\end{gathered}

\sigma_b = \frac{0.1265}{\sqrt{5}} = 0.057,\qquad \sigma_a = 0.1265\sqrt{\frac{1}{4} + \frac{1.5^2}{5}} = 0.106

\begin{gathered}\sigma_b = \frac{0.1265}{\sqrt{5}} = 0.057,\\[4pt]\sigma_a = 0.1265\sqrt{\frac{1}{4} + \frac{1.5^2}{5}} = 0.106\end{gathered}

よって $b = 2.04 \pm 0.06,\ a = 0.94 \pm 0.11$ 。

一般の式では $\operatorname{Cov}(p) \approx \sigma^2\,(J^{\mathsf T} J)^{-1}$ の対角成分の平方根(直線のときは上の式と同じ)。式がパラメータについて 1 次で $\sigma$ が分かっていれば近似なし、それ以外は見積もりです。

06残差を見る

最小二乗法は、どんな式でも「いちばんまし」な値を出してしまいます。式の形が合っているかは、合わせたあとの残差 $y_i - f(t_i)$ (データと式の差)を描いて確かめます。

左: 残差が 0 のまわりにでたらめに散らばる。右: 残差が U 字の形に並ぶ — 左は式の形が合っているとき。右は曲がったデータに直線を合わせたときで、残差が U の字に並ぶ

形が見えたら、合わせる値ではなく式の形そのものを見直します。