統計学 学習ノート

はじめに

「尤度が書けないときは、どうやって推定するの？」

前回の記事で、最尤推定の答案の書き方を学びました。とても強力な方法ですが、前提が1つあります。分布の形を、最後まで決めておかなければならないのです。

例えば、コンビニの売上を説明変数で予測する場面を考えましょう。誤差の平均が0で、ばらつきが一定だとは言えても、誤差の分布の形までは分からないことがよくありますよね。精密な地図がなくても、方角と距離さえ分かれば目的地に向かえるように、少ない仮定でできる推定があるはずです。

実は、推定の方法は、どこまで仮定するかで使い分けられます。平均と分散しか仮定しなくても、「この中では最良である」と証明できる推定量があるのです。

そこで登場するのが、今回の主役 線形推定（BLUE） です。BLUE は最良線形不偏推定量のことで、線形で不偏な推定量の中で分散が最小のものです。あわせて、モーメント法、最小二乗法、そして その他の手法 も扱います。

前回の記事では、最尤法を学びました。今回の記事では、モーメント推定量の漸近分散の求め方、最小二乗法の幾何学的な証明、BLUE の導出、そして最尤法以外の推定法について解説します。

モーメント法：推定量の分布まで求める

はじめに、モーメント法の性質の調べ方について解説します。

モーメント法 は、標本モーメント $m_k = \frac{1}{n}\sum_i X_i^k$ と母集団のモーメント $\mu_k(\theta)$ を等しいとおいて、母数について解く方法でした。準1級では、推定値を計算するところまで学びました。1級では「作って終わり」ではなく、その推定量がどれくらいよいのかまで述べます。

モーメント法は、目測の職人さんのようなものです。巻き尺で寸法を2、3か所測るだけで、すぐに見積もりを出してくれます。ただ、精密機械の測定（最尤法）と比べて、どのくらい精度が落ちるのかは知っておきたいですよね。

【モーメント推定量の性質を調べる手順】

1. 解いた形を、標本モーメントの関数として $\tilde{\theta} = h(m_1, \dots, m_p)$ と書く
2. 一致性： 大数の法則で $m_k \to \mu_k$ と確率収束し、$h$ が連続なら $\tilde{\theta} \to \theta$
3. 漸近正規性： 中心極限定理とデルタ法で、$\sqrt{n}(\tilde{\theta} - \theta)$ の極限分布を求める

1個のモーメントを使う場合、手順3は次の形です。デルタ法の証明は、漸近的性質の記事で行います。

$$\sqrt{n}\,\big(h(\bar{X}) - h(\mu)\big) \ \xrightarrow{d}\ N\big(0,\ h'(\mu)^2\, \sigma^2\big)$$

例題：モーメント推定量の漸近分散

密度 $f(x;\theta) = \theta x^{\theta-1}$（$0 < x < 1$、$\theta > 0$）からの標本について、$\theta$ のモーメント推定量 $\tilde{\theta}$ の漸近分散を求めよ。さらに、漸近分散が $\theta^2/n$ である最尤推定量に対する相対的な効率を求め、$\theta = 2$ の値を答えよ。

【答案】

1. $E[X] = \theta/(\theta + 1)$、$E[X^2] = \theta/(\theta + 2)$ より、$V[X] = \dfrac{\theta}{(\theta+1)^2(\theta+2)}$ である
2. $\bar{X} = \theta/(\theta+1)$ を解くと、$\tilde{\theta} = h(\bar{X})$、$h(m) = m/(1-m)$ である
3. $h'(m) = 1/(1-m)^2$ に $1 - \mu = 1/(\theta+1)$ を入れると、$h'(\mu) = (\theta+1)^2$ である
4. デルタ法より、$\sqrt{n}(\tilde{\theta} - \theta)$ の漸近分散は次のようになる

$$(\theta+1)^4 \times \frac{\theta}{(\theta+1)^2(\theta+2)} = \frac{\theta(\theta+1)^2}{\theta+2}$$

5. 最尤推定量との比は、$\dfrac{\theta^2}{\theta(\theta+1)^2/(\theta+2)} = \dfrac{\theta(\theta+2)}{(\theta+1)^2}$ である。分子は分母より1小さいので、この値は1未満である

答えは、漸近分散が $\theta(\theta+1)^2/\{n(\theta+2)\}$ で、効率は $\theta = 2$ のとき $8/9 \approx 0.889$ です。職人さんの目測は、精密機械の約9割の精度でした。

せっかくなので、別の例も見てみましょう。$N$ と $p$ がともに未知の二項分布 $B(N, p)$ では、$E[X] = Np$、$V[X] = Np(1-p)$ です。標本分散を $v$ とすると、$\tilde{p} = 1 - v/\bar{x}$、$\tilde{N} = \bar{x}/\tilde{p}$ となります。例えば、$\bar{x} = 6$、$v = 2.4$ なら、$\tilde{p} = 0.6$、$\tilde{N} = 10$ です。ただし、$v > \bar{x}$ のデータが出ると $\tilde{p} < 0$ になってしまいます。ありえない値が出ることは、モーメント法の弱点です。

もう1つの弱点は、十分統計量の関数になるとは限らないことです。$U(0, \theta)$ の $2\bar{X}$ は、十分統計量 $X_{(n)}$ の関数ではありません。この点は、次回の記事で改良します。

以上が、モーメント法についての解説です。

最小二乗法：直交性で最小を証明する

次に、最小二乗法が本当に最小を与えることの証明を説明します。

最小二乗法 は、線形模型 $\boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ で、$\|\boldsymbol{y} - X\boldsymbol{\beta}\|^2$ を最小にする方法でした。準1級では、微分して0とおき、正規方程式 $X^\top X\hat{\boldsymbol{\beta}} = X^\top\boldsymbol{y}$ を導きました。答案では、微分を使わずに最小を直接示せます。

川まで最短距離で歩きたいとき、川に垂直に歩きますよね。斜めに歩くと、ピタゴラスの定理の斜辺の分だけ遠くなります。最小二乗法も同じです。

$\hat{\boldsymbol{\beta}}$ を正規方程式の解とします。任意の $\boldsymbol{\beta}$ について、$\boldsymbol{y} - X\boldsymbol{\beta} = (\boldsymbol{y} - X\hat{\boldsymbol{\beta}}) + X(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})$ と分けて2乗します。交差項は $2(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\top X^\top(\boldsymbol{y} - X\hat{\boldsymbol{\beta}}) = 0$ なので、次の式が成り立ちます。

$$\|\boldsymbol{y} - X\boldsymbol{\beta}\|^2 = \|\boldsymbol{y} - X\hat{\boldsymbol{\beta}}\|^2 + \|X(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})\|^2 \ \ge\ \|\boldsymbol{y} - X\hat{\boldsymbol{\beta}}\|^2$$

等号は $X(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) = \boldsymbol{0}$ のときだけです。$X$ の列が線形独立なら、これは $\boldsymbol{\beta} = \hat{\boldsymbol{\beta}}$ を意味するので、最小点はただ1つです。

次に、誤差分散の不偏推定量を証明しましょう。$H = X(X^\top X)^{-1}X^\top$ を ハット行列 と呼びます。$\hat{\boldsymbol{y}} = H\boldsymbol{y}$ で、$H$ は $X$ の列空間への射影です。

【ハット行列の性質】

• 対称で、$H^2 = H$（影の影は同じ影）
• $HX = X$ なので、$(I - H)X = O$
• $\text{tr}\,H = \text{tr}\{(X^\top X)^{-1}X^\top X\} = p$

残差平方和は、$\text{RSS} = \boldsymbol{y}^\top(I - H)\boldsymbol{y} = \boldsymbol{\varepsilon}^\top(I - H)\boldsymbol{\varepsilon}$ です。スカラーをトレースとして扱うと、次のように計算できます。

$$E[\text{RSS}] = \text{tr}\{(I - H)\,E[\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^\top]\} = \sigma^2\,\text{tr}(I - H) = \sigma^2 (n - p)$$

したがって、$s^2 = \text{RSS}/(n - p)$ は $\sigma^2$ の不偏推定量です。正規分布の仮定は使っていません。

例えば、定数だけのモデルなら $p = 1$ で、$\text{RSS} = \sum_i (y_i - \bar{y})^2$ です。不偏分散の分母が $n - 1$ である理由は、$\text{tr}(I - H) = n - 1$ だったのです。

例題：制約のある最小二乗法

三角形の内角 $\theta_1, \theta_2, \theta_3$（和は180度）を別々に測り、測定値 $y_1, y_2, y_3$ を得た。測定誤差は互いに無相関で、平均0、分散 $\sigma^2$ とする。(1) 最小二乗推定量を導け。(2) $y = (61, 59, 63)$ のときの値を求めよ。(3) $V[\hat{\theta}_1]$ を求めよ。

【答案】

1. 制約 $\sum_j \theta_j = 180$ のもとで $\sum_j (y_j - \theta_j)^2$ を最小にする。ラグランジュ乗数 $\kappa$ を使い、$Q = \sum_j (y_j - \theta_j)^2 + 2\kappa(\sum_j \theta_j - 180)$ とおく
2. $\partial Q/\partial\theta_j = 0$ より $\theta_j = y_j - \kappa$ である。制約に代入すると、$\kappa = (\sum_j y_j - 180)/3$ である
3. 目的関数は狭義の凸関数で制約は線形なので、この停留点が最小点である。$\hat{\theta}_j = y_j - (\sum_l y_l - 180)/3$ となる
4. 値は、$\sum_j y_j = 183$ より $\kappa = 1$ で、$\hat{\theta} = (60, 58, 62)$ である
5. $\hat{\theta}_1 = (2y_1 - y_2 - y_3 + 180)/3$ なので、$V[\hat{\theta}_1] = (4 + 1 + 1)\sigma^2/9 = 2\sigma^2/3$ である

答えは (60, 58, 62) で、分散は $2\sigma^2/3$ です。合計のズレ3度を、3つの角に1度ずつ配る形になりました。「和は180度」という知識を使うだけで、分散は $\sigma^2$ から3分の2に減ります。

以上が、最小二乗法についての説明です。

線形推定（BLUE）：分散最小を導く

続いて、線形推定について解説します。

観測値の線形結合 $\sum_i a_i Y_i$ で表される推定量を、線形推定量 と呼びます。線形推定量の中で、不偏で分散が最小のものが 最良線形不偏推定量（BLUE） です。必要な仮定は、平均と分散の構造だけです。1つずつ見ていきましょう。

1. 精度の違う測定をまとめる

同じ量 $\mu$ を、精度の違う方法で測ったとします。友達の意見を聞くとき、信頼できる人の話ほど重く聞きますよね。その重みを計算で決めます。

$X_1, \dots, X_n$ は無相関で、$E[X_i] = \mu$、$V[X_i] = \sigma_i^2$（既知）とします。$T = \sum_i a_i X_i$ が不偏であるには、$\sum_i a_i = 1$ が必要です。分散 $\sum_i a_i^2\sigma_i^2$ の下限は、コーシー・シュワルツの不等式で出ます。

$$1 = \left(\sum_i a_i\right)^2 = \left(\sum_i a_i\sigma_i \cdot \frac{1}{\sigma_i}\right)^2 \le \left(\sum_i a_i^2\sigma_i^2\right)\left(\sum_i \frac{1}{\sigma_i^2}\right)$$

したがって $V[T] \ge 1/\sum_i (1/\sigma_i^2)$ で、等号は $a_i\sigma_i$ が $1/\sigma_i$ に比例するときに成り立ちます。

$$a_i = \frac{1/\sigma_i^2}{\sum_j 1/\sigma_j^2}, \qquad V[T] = \frac{1}{\sum_j 1/\sigma_j^2}$$

重みは分散の逆数に比例します。例えば、分散1の体重計と分散4の体重計で測って、50.0kgと52.5kgだったとします。重みは0.8と0.2で、BLUE の値は $0.8 \times 50 + 0.2 \times 52.5 = 50.5$kg、分散は0.8です。単純な平均の分散は $(1 + 4)/4 = 1.25$ で、BLUE の約1.56倍です。

2. ガウス・マルコフの定理を線形結合で証明する

線形模型 $\boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\varepsilon}$、$E[\boldsymbol{\varepsilon}] = \boldsymbol{0}$、$V[\boldsymbol{\varepsilon}] = \sigma^2 I$ で、$\boldsymbol{c}^\top\boldsymbol{\beta}$ を推定します。準1級では行列 $D$ を使って証明しました。ここでは、さきほどの直交分解を使います。

【証明】

1. $\boldsymbol{a}^\top\boldsymbol{y}$ が不偏であるとは、すべての $\boldsymbol{\beta}$ で $\boldsymbol{a}^\top X\boldsymbol{\beta} = \boldsymbol{c}^\top\boldsymbol{\beta}$ となること、すなわち $X^\top\boldsymbol{a} = \boldsymbol{c}$ である
2. $\boldsymbol{a}_0 = X(X^\top X)^{-1}\boldsymbol{c}$ とおくと、$X^\top\boldsymbol{a}_0 = \boldsymbol{c}$ で、$\boldsymbol{a}_0^\top\boldsymbol{y} = \boldsymbol{c}^\top\hat{\boldsymbol{\beta}}$ である
3. $\boldsymbol{d} = \boldsymbol{a} - \boldsymbol{a}_0$ は $X^\top\boldsymbol{d} = \boldsymbol{0}$ を満たす。$\boldsymbol{a}_0$ は $X$ の列の線形結合なので、$\boldsymbol{a}_0^\top\boldsymbol{d} = 0$ である
4. したがって $V[\boldsymbol{a}^\top\boldsymbol{y}] = \sigma^2\|\boldsymbol{a}\|^2 = \sigma^2(\|\boldsymbol{a}_0\|^2 + \|\boldsymbol{d}\|^2) \ge \sigma^2\|\boldsymbol{a}_0\|^2 = V[\boldsymbol{c}^\top\hat{\boldsymbol{\beta}}]$ である

ここでもピタゴラスの定理が活躍しました。余分な成分 $\boldsymbol{d}$ は、分散を増やすだけで何も得をしません。

3. 分散が等しくないとき

誤差の分散共分散行列が $\sigma^2 V$（$V$ は既知の正定値行列）のときはどうでしょうか？ $V = PP^\top$ と分解して、モデル全体に $P^{-1}$ を掛けます。すると誤差の分散が $\sigma^2 I$ になるので、ガウス・マルコフの定理が使えます。結果は次の 一般化最小二乗法 です。

$$\hat{\boldsymbol{\beta}}_G = (X^\top V^{-1} X)^{-1} X^\top V^{-1}\boldsymbol{y}, \qquad V[\hat{\boldsymbol{\beta}}_G] = \sigma^2 (X^\top V^{-1} X)^{-1}$$

具体的に見てみましょう。原点を通る回帰 $y_i = \beta x_i + \varepsilon_i$ で、$V[\varepsilon_i] = \sigma^2 x_i$（$x_i > 0$）とします。$V = \text{diag}(x_i)$ を入れると、$\hat{\beta}_G = \sum_i y_i/\sum_i x_i$ です。合計の比という、とても素朴な推定量が BLUE だったのです。

$x = 1, 2, 3$ なら、分散は $\sigma^2/6 \approx 0.167\sigma^2$ です。通常の最小二乗推定量 $\sum_i x_i y_i/\sum_i x_i^2$ の分散は $\sigma^2\sum_i x_i^3/(\sum_i x_i^2)^2 = 36\sigma^2/196 \approx 0.184\sigma^2$ で、確かに負けています。なお、1つ目の精度の違う測定は、$X = \boldsymbol{1}$ の場合にあたります。

試験のポイント：BLUE の「最良」は、線形かつ不偏な推定量の中での話です。誤差の分布を仮定せず、平均と分散の構造だけで証明できる点が強みです。

以上が、線形推定についての解説です。

その他の手法：道具箱の残りを開ける

最後に、その他の手法を見ていきましょう。

ここまでの方法が使えない場面でも、推定は必要です。工具箱に金づちしかないと、ネジは回せませんよね。代表的な道具を4つ紹介します。

【その他の手法】

1. 分位点法： 母集団と標本の分位点をそろえる
2. 最小カイ二乗法： 度数のズレを最小にする
3. M推定： 外れ値の影響を抑えた方程式を解く
4. ベイズ推定： 事前の知識とデータを合わせる

1つずつ見ていきましょう。

1. 分位点法

モーメントの代わりに、分位点をそろえる方法です。モーメントが存在しない分布でも使えます。

例えば、コーシー分布（位置 $\theta$、尺度 $\gamma$）は期待値が存在しないので、モーメント法は使えません。しかし、中央値は $\theta$ で、四分位数は $\theta \pm \gamma$ です。標本の第1四分位数が2、中央値が5、第3四分位数が8なら、$\tilde{\theta} = 5$、$\tilde{\gamma} = (8 - 2)/2 = 3$ です。

2. 最小カイ二乗法

度数データで、観測度数 $O_j$ と期待度数 $E_j(\theta)$ のズレを最小にする方法です。

$$\chi^2(\theta) = \sum_{j} \frac{\{O_j - E_j(\theta)\}^2}{E_j(\theta)}$$

前回の遺伝子型の例（度数 30, 50, 20）で最小点を数値で求めると約0.55で、最尤推定値とほぼ同じ値になります。一般に、最小カイ二乗推定量は最尤推定量と漸近的に同等です。

3. M推定

$\sum_i \psi(x_i - \theta) = 0$ の解を推定値とする方法です。最尤法を一般化したもので、$\psi$ の選び方で性格が変わります。

【$\psi$ の選び方】

• $\psi(u) = u$： 標本平均
• $\psi(u) = \text{sign}(u)$： 標本中央値
• $\psi(u) = \max\{-k, \min(u, k)\}$： ハンバーの推定量

ハンバーの推定量は、1人の発言力に上限 $k$ を設ける議長のようなものです。例えば、データ 1, 2, 3, 4, 20 に $k = 3$ を使いましょう。20だけが上限に当たると考えて、$(1 - \theta) + (2 - \theta) + (3 - \theta) + (4 - \theta) + 3 = 0$ を解くと、$\theta = 3.25$ です。このとき、他の4個の差はすべて3以内で、20との差は3を超えるので、仮定と合っています。平均の6と中央値の3の間で、外れ値に振り回されない値になりました。

漸近分散は $E[\psi^2]/(E[\psi'])^2$ を $n$ で割った値です。$\psi(u) = -\frac{d}{du}\log f(u)$ とすると、位置母数の最尤推定量に一致します。

4. ベイズ推定

母数に事前分布を置き、事後分布の平均などを推定値とする方法です。例えば、$p$ の事前分布を $\text{Beta}(a, b)$ とすると、$n$ 回中 $x$ 回の成功を見たあとの事後平均は $(a + x)/(a + b + n)$ です。$a = b = 1$、$n = 10$、$x = 3$ なら $4/12 = 1/3$ で、最尤推定値の0.3より、事前の平均0.5に少し引き寄せられます。詳しくはベイズ法の記事で扱います。

ここまでの方法を整理しましょう。

方法	必要な仮定	強み	弱み
最尤法	分布の形	漸近的に最良	分布を間違えると弱い
モーメント法	モーメントの式	手軽	効率が落ちる。ありえない値も出る
最小二乗法・BLUE	平均と分散の構造	分布の仮定が不要	線形不偏の中でのみ最良
分位点法	分位点の式	モーメントが不要	効率が落ちる
M推定	$\psi$ の選択	外れ値に強い	計算が要る

以上が、その他の手法についての解説です。

細かい手法を忘れてしまっても問題ありません。ここでは、「仮定が少ないほど使える場面は広く、仮定が多いほど精度は高い」という点だけ押さえておきましょう。

今回のまとめ

それでは、今回の内容を振り返ります。

今回は、最尤法以外の推定法について紹介しました。目測の職人さんのモーメント法、垂線を下ろす最小二乗法、分散の逆数で重みをつける BLUE と、それぞれ必要な仮定が違います。問題文がどこまで仮定しているかを読むことが、方法を選ぶ第一歩です。

今回のポイントを振り返りましょう。

振り返り

• モーメント推定量の性質は、大数の法則とデルタ法で調べる。最尤推定量より効率が落ちることが多い。
• 最小二乗法の最小性は直交分解で示せる。$\text{tr}(I - H) = n - p$ が $s^2$ の分母の理由である。
• 精度の違う測定の BLUE は、分散の逆数に比例する重みの平均である。
• ガウス・マルコフの定理も、余分な成分が直交することから示せる。分散が等しくなければ一般化最小二乗法を使う。
• 分位点法、最小カイ二乗法、M推定、ベイズ推定は、最尤法が使いにくい場面の道具である。

キーワードのおさらい

• モーメント法: 標本モーメントと母集団のモーメントを等しいとおいて母数を解く方法
• 最小二乗法: 残差の2乗和を最小にする方法
• 線形推定（BLUE）: 観測値の線形結合による推定。その中で不偏かつ分散最小のものが最良線形不偏推定量
• ハット行列: $\boldsymbol{y}$ を予測値へ写す射影行列 $X(X^\top X)^{-1}X^\top$
• 一般化最小二乗法: 誤差の共分散行列の逆行列で重みをつける最小二乗法
• その他の手法: 分位点法、最小カイ二乗法、M推定、ベイズ推定など
• M推定: $\sum_i \psi(x_i - \theta) = 0$ の解を推定値とする方法

さて、推定量の作り方がそろいました。では、できあがった推定量が「よい」ことは、どうやって証明すればよいのでしょうか？

次回は、不偏性、一致性、十分性、有効性を数学的に立証する方法について紹介します。

それでは、次の記事へ進んでいきましょう。