最小二乗法: なぜ分散の逆数で重み付けすると良いのか？

はじめに

投稿論文の査読コメントに返信するのに、久しぶりに5次元の連立微分方程式を解く機会があったのですが、計算能力がかなり落ちていることに気づきました。
思い返すとD論を書き始めた辺りから計算が面倒で論文に載っている式をそのまま使うようになった気がします。学部時代は学費と生活費を稼ぐためにバイトで忙しかったため、歩きながら頭の中で微積分を復習するとかやっていました。今は簡単な計算すらやり方を忘れて鉛筆を握ったまま固まり、数秒後ネット検索をすることが増えました。
実用上は、調べて多少時間かけても解ければ研究を行う上では問題はないです。
が、論文/G検定がひと段落したため、計算練習のためにこれまで証明をせずに結果だけ知っていた事柄の証明を行なっていこうと思います。ただし、厳密な数学の説明は避けます。
初回は、最小二乗法を行う特になぜ分散(観測誤差の２乗)の逆数で重み付けするのか？です。

最小二乗法

最小二乗法は得られたデータに最もよく当てはまる関数を求める方法です。"最もよく当てはまる"の定義はたくさんありますが、一般的にユークリッド距離で定義された目的関数と近似関数との残差の２乗が最小になる時が最も良いと判断します。具体的に１次関数で近似すると下記の$${\sigma^2}$$を最小にすることを言います。

$${\sigma^2 = \sum_{i=1}^n (y_i - ax_i - b)^2}$$

$${\sigma^2}$$を最小にすると言うことは、a, bで偏微分した値が０になるa, bを求めれば良いため(ラグランジュの未定乗数法)

$${\frac{\partial \sigma}{\partial a} = -\sum_{i=1}^n x_i(y_i - ax_i - b)=0}$$
$${\frac{\partial \sigma}{\partial b} = -\sum_{i=1}^n (y_i - ax_i - b)=0}$$

であり、上記の連立方程式は下記の行列式で表現できます。

$$
\begin{pmatrix}\sum_{i=1}^n x_i^2 & \sum_{i=1}^n x_i \\ \sum_{i=1}^nx_i&n \\\end{pmatrix}\begin{pmatrix}a \\b \\\end{pmatrix} = \begin{pmatrix} \sum_{i=1}^n x_iy_i \\ \sum_{i=1}^n y_i \\\end{pmatrix}
$$

左辺の逆行列を両編にかけると

$$
\begin{pmatrix}a \\b \\\end{pmatrix} = \frac{1}{n\sum_{i=1}^n x_i^2 - (\sum_{i=1}^nx_i)^2}  \begin{pmatrix}n & -\sum_{i=1}^n x_i \\ -\sum_{i=1}^nx_i& \sum_{i=1}^n x_i^2 \\\end{pmatrix} \begin{pmatrix} \sum_{i=1}^n x_iy_i \\ \sum_{i=1}^n y_i \\\end{pmatrix}
$$

最尤推定法

同様に最尤推定法を用いても上記と同じ結果が得られることが知られています。誤差が平均0、分散$${\alpha^2}$$の確率密度関数$${f_(x; 0, \alpha)}$$に従うとすると

$$
{f_(x; 0, \alpha) = \frac{1}{\sqrt{2\pi}{\alpha}}} exp{-\frac{\sigma^2}{2\alpha^2}}
$$

と$${x_i}$$を生成する確率密度関数は表現できる。ここからn回$${x_i}$$を生成した時の尤度関数$${L(\theta)}$$は

$$
{{L(\theta)} = \prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}{\alpha}}} exp{-\frac{\sigma_i^2}{2\alpha^2}}
$$

この尤度関数を最大化するa, bの組み合わせを考えれば良いため、

$${\log(L(\theta))=-\frac{n}{2}\log(2\pi{\alpha}^2)- \sum_{i=1}^{n}{-\frac{(y_i - ax_i - b)^2}{2\alpha^2}}}$$

をa,bで微分すれば良い訳です。右辺第一項はa,bで微分すると消えるため、第二項のみを考えればよく、結局ユークリッド距離で残差を定義した時の結果と合流するのがわかると思います。

それぞれの誤差の分散が異なる場合の最尤推定法

上記より、ユークリッド距離で残差を評価した時の最小二乗法と最尤推定法の結果が一致することがわかりました。必然的に一致するのか私は理解していませんが、ユークリッド距離で残差を定義すると中心極限定理より正規分布になるからではないかと思います(確証はないため注意してください)。
上記では、それぞれの誤差の分散が全て同一であることを仮定して、最尤推定法を行いました。しかし、実際のデータを扱うとデータによって誤差が異なることはままあります。
例えば、光を使用して大気気温を測定する観測機(ライダー)のデータを私は解析していました。光が減衰するため遠い場所(高い高度)の観測値の誤差は増え、近い場所は逆に小さくなる特性があります。高度1 km : 10±0.5℃, 高度2 km: 0±1℃のようになります。
直感的に、高度2 kmより高度1 kmのデータに重みを大きくした方が良い推定ができそうな気がすると思います。つまり、データ毎の誤差の分散を考慮して重み付けして最小二乗法を解くのが良い訳です。では、どのような重みをつけるのが良いか？と言う答えがタイトルです。

分散の変化を考慮した時の尤度関数は、下記のように表すことができます。

$$
{{L(\theta)} = \prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}{\alpha_i}}} exp{-\frac{\sigma_i^2}{2\alpha_i^2}}
$$

分散$${\alpha}$$がデータ毎に変化するため、上記の最尤関数に添字iを加えただけです。対数を取ると

$${\log(L(\theta))=-\frac{n}{2}\log(2\pi{\alpha_i}^2)- \sum_{i=1}^{n}{-\frac{(y_i - ax_i - b)^2}{2\alpha_i^2}}}$$

ここまで来ると結果が見えてくると思いますが、念の為a, bで微分してみましょう。

$${\frac{\partial\log(L(\theta))}{\partial a}= \sum_{i=1}^{n}{-a\frac{(y_i - ax_i - b)}{\alpha_i^2}} = 0}$$
$${\frac{\partial\log(L(\theta))}{\partial b}= \sum_{i=1}^{n}{-\frac{(y_i - ax_i - b)}{\alpha_i^2}} = 0}$$

上記の式が「最小二乗法」の節で導出した式と非常によく似た結果になっています。が、$$frac{1}{\alpha_i^2}$$が残っています。
これは、それぞれのデータの残差に分散(観測誤差の２乗)の逆数を掛けた値の合計値が最小になる時の結果と一致します。

結論

誤差がそれぞれのデータ点で異なる場合は、その分散で重みをつけて最小二乗法を解いた方が尤度が高くなると言うことです。このため、誤差の２乗の逆数で重みをつける訳です。