各種推定法 — 最尤法とモーメント法

準1級-10　応用統計の手法

各種推定法 — 最尤法とモーメント法

最尤法
モーメント法
最小二乗法
線形模型

はじめに

「母数の値って、結局どうやって決めるの？」

2級では、母平均は標本平均で、母比率は標本比率で推定しました。どちらも当たり前に思えますよね。では、ガンマ分布の形状母数はどうでしょうか？標本の何を計算すれば、形状母数の推定値になるのか、すぐには思いつきません。

例えば、ガチャを10回引いて3回当たったとしましょう。当たる確率は30%だろうと、誰もが思います。でも、なぜ30%なのでしょうか？実は、この「なんとなく」の裏には、どんな分布にも使える一般的な手順が隠れています。

そこで登場するのが、今回の主役 最尤法 です。あわせて、もう2つの基本の作り方である モーメント法 と 最小二乗法 も扱います。

前回の記事では、情報を失わない要約である十分統計量を学びました。今回の記事では、推定量を作る3つの方法と、回帰分析の土台である線形模型について解説します。

モーメント法：標本と母集団のモーメントをそろえる

はじめに、最も素朴な作り方であるモーメント法について解説します。

モーメント法 は、標本のモーメントと母集団のモーメントを等しいとおいて、母数について解く方法です。標本の $k$ 次モーメント $m_k$ と、母集団の $k$ 次モーメントを、次のように書きます。

m_k = \frac{1}{n}\sum_{i=1}^{n} X_i^k, \qquad \mu_k(\theta) = E[X^k]

母数が $p$ 個あるなら、 $m_k = \mu_k(\theta)$ （ $k = 1, \dots, p$ ）という $p$ 本の連立方程式を解きます。大数の法則により $m_k$ は $\mu_k$ に近づくので、理屈も素直です。

洋服の仕立てに例えると、これは「お客さんの肩幅と着丈を測って、型紙の寸法をそろえる」作業です。形の型紙はもう決まっていて、寸法だけを合わせるイメージです。

具体的に見てみましょう。正規分布 $N(\mu, \sigma^2)$ では、 $E[X] = \mu$ 、 $E[X^2] = \sigma^2 + \mu^2$ です。方程式を解くと、次のようになります。

\hat{\mu} = \bar{X}, \qquad \hat{\sigma}^2 = m_2 - m_1^2 = \frac{1}{n}\sum_{i=1}^{n}(X_i - \bar{X})^2

分母が $n-1$ ではなく $n$ である点に注意しましょう。

せっかくなので、別の例も見てみましょう。形状 $\alpha$ 、尺度 $\beta$ のガンマ分布では、期待値は $\alpha\beta$ 、分散は $\alpha\beta^2$ です。分散を期待値で割れば $\beta$ が出るので、 $\hat{v} = m_2 - m_1^2$ として次の式が得られます。

\hat{\beta} = \frac{\hat{v}}{\bar{X}}, \qquad \hat{\alpha} = \frac{\bar{X}^2}{\hat{v}}

例えば、 $\bar{x} = 6$ 、 $\hat{v} = 12$ なら、 $\hat{\beta} = 12 \div 6 = 2$ 、 $\hat{\alpha} = 36 \div 12 = 3$ です。冒頭の疑問だった形状母数も、これで計算できました。

モーメント法は手軽ですが、弱点もあります。 $U(0, \theta)$ を考えましょう。 $E[X] = \theta/2$ なので、 $\hat{\theta} = 2\bar{X}$ です。データが 1, 2, 2, 3, 9 だったら、 $\bar{x} = 3.4$ で $\hat{\theta} = 6.8$ となります。でも、データに9があるので、 $\theta$ は少なくとも9のはずです。推定値が、ありえない値になってしまいました。

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

最尤法：データが最も起きやすい母数を選ぶ

次に、推定量の主役である最尤法について説明します。

最尤法（maximum likelihood method）は、手元のデータが起きる確率を最も大きくする母数を、推定値に選ぶ方法です。尤は「もっともらしい」という意味です。

探偵に例えると、分かりやすいです。現場に残った証拠を見て、「この人が犯人なら、この証拠が残る確率は高いか」を容疑者ごとに比べます。そして、証拠を最もよく説明する人を第一の容疑者にするのです。ここで証拠がデータ、容疑者が $\theta$ の候補です。

式で書きましょう。データ $x_1, \dots, x_n$ を固定して、同時確率（密度）を $\theta$ の関数と見たものを 尤度関数 と呼びます。その対数が 対数尤度 です。

L(\theta) = \prod_{i=1}^{n} f(x_i;\theta), \qquad \ell(\theta) = \log L(\theta) = \sum_{i=1}^{n} \log f(x_i;\theta)

$L(\theta)$ を最大にする $\hat{\theta}$ が、最尤推定量 です。対数は単調増加なので、最大になる場所は $L$ でも $\ell$ でも同じです。積が和に変わるので、ふつうは $\ell$ を使います。

【最尤推定量の求め方】

尤度関数を書き、対数をとる
$\theta$ で微分して0とおく（これを尤度方程式と呼ぶ）
解いて、最大であることを確かめる

冒頭のガチャを解いてみましょう。10回中3回当たる並びの確率は $L(p) = p^3(1-p)^7$ です。候補の値を入れて比べます。

$p$	0.2	0.3	0.5
$L(p)$	0.00168	0.00222	0.00098

0.3が最も大きいですね。微分でも確かめましょう。 $\ell(p) = 3\log p + 7\log(1-p)$ なので、次の式を解きます。

\ell'(p) = \frac{3}{p} - \frac{7}{1-p} = 0 \quad \Rightarrow \quad \hat{p} = 0.3

「なんとなく30%」の正体は、最尤推定値だったのです。

次に、母数が2つある例です。 $N(\mu, \sigma^2)$ の対数尤度は次の形です。

\ell(\mu, \sigma^2) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i - \mu)^2

$\mu$ と $\sigma^2$ でそれぞれ偏微分して0とおくと、 $\hat{\mu} = \bar{x}$ 、 $\hat{\sigma}^2 = \frac{1}{n}\sum_i (x_i - \bar{x})^2$ が得られます。ここでも分母は $n$ です。

微分が使えない例もあります。 $U(0, \theta)$ の尤度は、 $\theta \ge \max_i x_i$ のとき $\theta^{-n}$ 、それ以外では0です。 $\theta^{-n}$ は $\theta$ が小さいほど大きいので、許される範囲の左端を選びます。つまり、 $\hat{\theta} = \max_i X_i$ です。さきほどのデータなら $\hat{\theta} = 9$ で、モーメント法の6.8のような矛盾は起きません。前回の分解定理から分かるように、最尤推定量が1つに定まるなら、それは十分統計量の関数になります。

最尤法には、もう1つ便利な性質があります。 $\hat{\theta}$ が $\theta$ の最尤推定量なら、 $g(\theta)$ の最尤推定量は $g(\hat{\theta})$ です。これを 不変性 と呼びます。

例題：最尤推定とモーメント推定を比べる

$f(x;\theta) = \theta x^{\theta - 1}$ （ $0 < x < 1$ ）からの大きさ4の標本で、 $\sum_i \log x_i = -2$ 、 $\bar{x} = 0.65$ でした。 $\theta$ の最尤推定値と、モーメント法の推定値を求めましょう。

【解き方】

対数尤度： $\ell(\theta) = n\log\theta + (\theta - 1)\sum_i \log x_i$
微分して0とおく： $n/\theta + \sum_i \log x_i = 0$ なので、 $\hat{\theta} = -n/\sum_i \log x_i = 4 \div 2 = 2$
モーメント法： $E[X] = \int_0^1 \theta x^{\theta}\,dx = \theta/(\theta + 1)$ なので、 $\bar{x} = \theta/(\theta+1)$ を解いて $\tilde{\theta} = \bar{x}/(1 - \bar{x})$
代入する： $\tilde{\theta} = 0.65 \div 0.35 \approx 1.86$

答えは最尤推定値2、モーメント推定値約1.86です。同じデータでも、方法が違えば値も変わります。

例題：不変性を使う

電球の寿命が指数分布に従うとします。5個の寿命が 120, 80, 200, 60, 40 時間のとき、200時間より長くもつ確率の最尤推定値を求めましょう。

【解き方】

平均 $\theta$ の指数分布の最尤推定量は $\hat{\theta} = \bar{x}$ なので、 $\hat{\theta} = 500 \div 5 = 100$
求めたいのは $g(\theta) = P(X > 200) = e^{-200/\theta}$
不変性より、 $g(\hat{\theta}) = e^{-2} \approx 0.135$

答えは約0.135です。

補足（ここは読み飛ばしてOK）：ガンマ分布の形状母数のように、尤度方程式が式変形では解けないこともあります。その場合は、 $\theta_{\text{new}} = \theta - \ell'(\theta)/\ell''(\theta)$ を繰り返すニュートン法などの数値計算で求めます。

試験のポイント：範囲が母数で決まる分布では、尤度方程式は使えません。尤度の形を直接見て、最大の場所を探します。

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

最小二乗法：ズレの2乗を最小にする

続いて、最小二乗法について解説します。

2級の回帰分析で学んだように、最小二乗法 は、観測値とモデルの予測値のズレの2乗和を最小にする方法です。

S(\theta) = \sum_{i=1}^{n} \{y_i - g(x_i;\theta)\}^2

例えば、各データ点と直線をばねでつないだ場面を想像してください。ばねのエネルギーは伸びの2乗に比例するので、手を離すと、直線はエネルギーの合計が最小になる位置に落ち着きます。それが最小二乗法の答えです。

最も単純な $g = a$ （定数）の場合は、 $S(a) = \sum_i (y_i - a)^2$ を $a$ で微分すると $\hat{a} = \bar{y}$ が出ます。標本平均も最小二乗推定量なのです。

$g = a + bx$ の場合は、 $a$ と $b$ で偏微分して0とおきます。

\sum_{i}(y_i - a - bx_i) = 0, \qquad \sum_{i} x_i (y_i - a - bx_i) = 0

この2本を 正規方程式 と呼びます。解くと、2級で使った $\hat{b} = S_{xy}/S_{xx}$ 、 $\hat{a} = \bar{y} - \hat{b}\bar{x}$ になります。

ここがポイントです。最小二乗法は、誤差の分布を何も仮定していません。最尤法は分布を決めないと使えませんが、最小二乗法はズレの大きさだけで計算できます。

では、2つの方法は関係がないのでしょうか？実は、誤差が $N(0, \sigma^2)$ に従うと仮定すると、対数尤度は次の形になります。

\ell = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{S(\theta)}{2\sigma^2}

$\ell$ を最大にすることは、 $S(\theta)$ を最小にすることと同じです。つまり、誤差が正規分布なら、最小二乗推定量は最尤推定量に一致します。

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

線形模型：行列で一度に解く

最後に、線形模型を見ていきましょう。

線形模型（線形モデル）は、応答 $y$ を説明変数の線形結合と誤差の和で表すモデルです。 $n$ 個の観測をまとめて、行列で次のように書きます。

\boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad E[\boldsymbol{\varepsilon}] = \boldsymbol{0}, \quad V[\boldsymbol{\varepsilon}] = \sigma^2 I

$\boldsymbol{y}$ は $n$ 次元の応答、 $X$ は $n \times p$ の 計画行列、 $\boldsymbol{\beta}$ は $p$ 次元の未知の係数です。「線形」とは、係数 $\boldsymbol{\beta}$ について1次式という意味です。したがって、 $y = \beta_0 + \beta_1 x + \beta_2 x^2$ のような2次式も線形模型です。単回帰、重回帰、分散分析は、 $X$ の作り方が違うだけで、どれもこの1つの式に入ります。万能の型紙のようなものです。

ズレの2乗和は $S(\boldsymbol{\beta}) = (\boldsymbol{y} - X\boldsymbol{\beta})^\top(\boldsymbol{y} - X\boldsymbol{\beta})$ です。 $\boldsymbol{\beta}$ で微分して0とおくと、 $X^\top(\boldsymbol{y} - X\boldsymbol{\beta}) = \boldsymbol{0}$ となり、行列の形の正規方程式が得られます。 $X^\top X$ が逆行列をもつとき、解は次のとおりです。

X^\top X \hat{\boldsymbol{\beta}} = X^\top \boldsymbol{y}, \qquad \hat{\boldsymbol{\beta}} = (X^\top X)^{-1} X^\top \boldsymbol{y}

記号の意味を確認しましょう。 $X^\top$ は $X$ の転置です。 $X^\top(\boldsymbol{y} - X\hat{\boldsymbol{\beta}}) = \boldsymbol{0}$ は、残差が $X$ のすべての列と直交することを表しています。つまり、 $X\hat{\boldsymbol{\beta}}$ は、 $\boldsymbol{y}$ を $X$ の列が張る平面へ落とした影です。真上から光を当てたとき、棒の先の影が棒の先に最も近い平面上の点になるのと同じです。

推定量の性質は、次の3つです。

【線形模型の推定量の性質】

不偏： $E[\hat{\boldsymbol{\beta}}] = \boldsymbol{\beta}$
分散共分散行列： $V[\hat{\boldsymbol{\beta}}] = \sigma^2 (X^\top X)^{-1}$
誤差分散の不偏推定量： $s^2 = \text{RSS}/(n - p)$ （ $\text{RSS}$ は残差平方和）

なお、誤差が正規分布のとき、 $\sigma^2$ の最尤推定量は $\text{RSS}/n$ です。分母が違う点は、正規分布の $\hat{\sigma}^2$ と同じ話です。

例題：行列で回帰係数を求める

4組のデータ $(x, y)$ が (1, 2), (2, 3), (3, 5), (4, 6) のとき、 $y = \beta_0 + \beta_1 x + \varepsilon$ の $\hat{\boldsymbol{\beta}}$ と $s^2$ を求めましょう。

【解き方】

$X$ は1列目がすべて1、2列目が $x$ 。 $X^\top X$ と $X^\top \boldsymbol{y}$ を計算する

X^\top X = \begin{pmatrix} 4 & 10 \\ 10 & 30 \end{pmatrix}, \qquad X^\top \boldsymbol{y} = \begin{pmatrix} 16 \\ 47 \end{pmatrix}

行列式は $4 \times 30 - 10 \times 10 = 20$ なので、逆行列を掛ける

\hat{\boldsymbol{\beta}} = \frac{1}{20}\begin{pmatrix} 30 & -10 \\ -10 & 4 \end{pmatrix}\begin{pmatrix} 16 \\ 47 \end{pmatrix} = \begin{pmatrix} 0.5 \\ 1.4 \end{pmatrix}

予測値は 1.9, 3.3, 4.7, 6.1 で、残差は 0.1, −0.3, 0.3, −0.1
$\text{RSS} = 0.01 + 0.09 + 0.09 + 0.01 = 0.2$ なので、 $s^2 = 0.2 \div (4 - 2) = 0.1$

答えは $\hat{\beta}_0 = 0.5$ 、 $\hat{\beta}_1 = 1.4$ 、 $s^2 = 0.1$ です。傾きの分散の推定値は、 $0.1 \times 4 \div 20 = 0.02$ となります。

ここまでの3つの方法を比べてみましょう。

方法	考え方	必要な仮定	特徴
モーメント法	モーメントをそろえる	モーメントの式	計算が手軽。ありえない値が出ることも
最尤法	データが最も起きやすい母数	分布の形	汎用で性質がよい。数値計算が要ることも
最小二乗法	ズレの2乗和を最小に	回帰式の形	分布の仮定が不要。正規誤差なら最尤法に一致

試験のポイント： $\hat{\boldsymbol{\beta}} = (X^\top X)^{-1}X^\top\boldsymbol{y}$ と $V[\hat{\boldsymbol{\beta}}] = \sigma^2(X^\top X)^{-1}$ は、2次の逆行列を手で計算できるようにしておきましょう。

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

今回のまとめ

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

今回は、推定量を作る3つの方法と、線形模型について紹介しました。寸法をそろえるモーメント法、証拠から犯人を探す最尤法、ばねのエネルギーを最小にする最小二乗法と、発想はそれぞれ違います。その中で最尤法は、どんな分布でも同じ手順で使える、最も汎用な方法です。

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

振り返り

モーメント法は標本と母集団のモーメントを等しいとおく。手軽だが、ありえない値が出ることもある。
最尤法は尤度を最大にする。対数をとって微分するのが基本で、範囲が母数で決まる分布では直接考える。
最尤推定量には不変性があり、 $g(\theta)$ の推定量は $g(\hat{\theta})$ でよい。
最小二乗法は分布の仮定なしに使え、誤差が正規分布なら最尤法と一致する。
線形模型の解は $(X^\top X)^{-1}X^\top\boldsymbol{y}$ で、残差は $X$ の列と直交する。

キーワードのおさらい

モーメント法: 標本モーメントと母集団モーメントを等しいとおいて母数を解く方法
最尤法: 尤度関数を最大にする値を推定値とする方法
最小二乗法: 観測値と予測値のズレの2乗和を最小にする方法
線形模型: 応答を、係数について1次式の部分と誤差の和で表すモデル

さて、同じデータでも、方法によって推定値が変わることが分かりました。では、どの推定量が「よい」推定量なのでしょうか？分散の分母は $n$ と $n-1$ のどちらがよいのでしょうか？

次回は、不偏性や有効性など、推定量のよさを測る基準である「点推定の性質」について紹介します。

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

練習問題5問開いて挑戦する

問1. 形状

\alpha

、尺度

\beta

のガンマ分布（期待値

\alpha\beta

、分散

\alpha\beta^2

）からの標本で、

\bar{x} = 8

、

\hat{v} = m_2 - m_1^2 = 16

だった。モーメント法の推定値の組として正しいものはどれか。

問2. 幾何分布

P(X = x) = p(1-p)^{x-1}

（

x = 1, 2, \dots

）からの標本で、

\bar{x} = 4

だった。

p

の最尤推定値として正しいものはどれか。

問3.

U(0, \theta)

からの標本が 2, 5, 3, 7, 4 のとき、

\theta

の最尤推定値として正しいものはどれか。

問4. ポアソン分布からの標本で

\bar{x} = 1.2

だった。

P(X = 0)

の最尤推定値に最も近いものはどれか。

問5. 線形模型で

X^\top X = \begin{pmatrix} 5 & 10 \\ 10 & 30 \end{pmatrix}

、

X^\top \boldsymbol{y} = \begin{pmatrix} 20 \\ 50 \end{pmatrix}

のとき、

\hat{\boldsymbol{\beta}}

として正しいものはどれか。

一覧へ戻る