【Rで統計入門】(4) 最尤法でパラメタ推定 - ケィオスの時系列解析メモランダム

　このブログでは、「データの背後にある法則を見つける」という問いを強調してきました。そして前回の記事では、データが従う法則の一つとして確率分布があり、その分布のパラメタを推定する問題について考えました。この問題に対する有力な解決策の一つが、最尤推定法（Maximum Likelihood Estimate (MLE) Method） です。最尤推定法は「観測されたデータが最も起こりやすくなるように確率分布のパラメタを設定する」という考え方に基づいています。

　この手法を確立したのが、ロナルド・フィッシャー（Ronald A. Fisher） です。以前の記事で、フィッシャーがモーメント法を開発したピアソンと対立関係にあったことについて触れました。フィッシャーは20世紀初頭に統計学の発展に多大な貢献を果たし、最尤法をはじめとする数多くの画期的なアイデアを生み出しました。彼は間違いなく天才と称される存在ですが、個人的には、能力の差を見せつけられてもなお、威張り散らすピアソンの方が、人間臭くて親近感がわきます。

　ということで、今回は最尤法の考え方を説明してみます。前回の記事は以下です。

【Rで統計入門】(3) モーメント法を使った確率分布のパラメタ推定 - ケィオスの時系列解析メモランダム

最尤法とは？
尤度のイメージ
尤度を表す式
尤度関数の対数をとる
最尤推定の例
最後に一言
付録：指数関数の尤度を最大化
- 尤度関数 (Likelihood Function)
- 対数尤度関数の最大化するパラメタの計算

最尤法とは？

　最尤推定法 (さいゆうすいていほう)とは、確率分布のパラメータを推定する方法のひとつです。ここからは、ちょっと短く「最尤法」と呼びます。最尤法では、"Likelihood"という考え方を使い、観測されたデータは最も起こりやすいものが起こったという仮定で、パラメタの値を求めます。

　英単語のLikelihoodには、「ありそうなこと、見込み、可能性」という意味があります。統計学では、Likelihoodの訳は「尤度（ゆうど）」になります。尤度って統計でしか使わないような堅い単語ですが、英語でのLikelihoodは、

There is a high likelihood of rain tomorrow. （明日は雨が降る可能性が高い。）
The likelihood of winning the lottery is very low. （宝くじに当たる可能性は非常に低い。）

のように使うことができます。日本訳で、

明日は雨が降る尤度が高い。
宝くじに当たる尤度は非常に低い。

とか言いませんので、尤度ってなんだか難しそうです。

　上の例文では、「尤度」を、「確率」と言い換えても意味が通じます。

明日は雨が降る確率が高い。
宝くじに当たる確率は非常に低い。

　しかし、統計では、「確率は未来の出来事について考えるもの」であり、「尤度は過去に得られたデータをもとに仮説のもっともらしさを評価するもの」として、意味が区別されています。とはいえ、尤度も、確率みたいなものです。

尤度のイメージ

　統計における「尤度（Likelihood）」は、

「あるパラメータを仮定したとき、観測データがどれくらい説明できるか？」

を評価するものです。

　尤度の大小の解釈は以下のようになります。

「尤度が高い」 = 観測データが、そのパラメータのもとで起こる可能性が高い
「尤度が低い」 = 観測データが、そのパラメータのもとで起こるのは不自然（起こる可能性が低い）

　たとえば、コイン投げで、「50回中40回表が出た」データに対して、

「このコインが公平である（表が出る確率 $p=0.5$ と仮定すると、このデータの尤度はどれくらいか？」
「このコインは表裏のバランスが崩れていて、表が出る確率 $p=0.8$ だったとすると、このデータの尤度はどれくらいか？」

と考えます。

　50回中40回も表が出たならば、このコインの表と裏が出る確率が同じとは思えないと思います。この「思えない」程度を尤度で定量的に評価します。

尤度を表す式

　確率分布は、パラメータと呼ばれる数値によって形が決まります。たとえば、正規分布（ガウス分布）は平均 $\displaystyle{ \mu }$ と標準偏差 $\displaystyle{ \sigma }$ というパラメタによって形が決まります。

　データが $n$ 個、

$\displaystyle{ x_1, x_2, \dots, x_n }$

のように観測されたとします。このとき、確率分布のパラメタ $\displaystyle{ \theta }$ がある値だとして、それらのデータの組み合わせが発生する「確率のようなもの」を考えます (「確率のようなもの」=「尤度」です)。

　確率密度関数 (Probability Density Function: PDF)を $\displaystyle{ f(x | \theta) }$ とすると、 $x _ 1$ が出る確率は、 $\displaystyle{ f(x _ 1| \theta) }$ に比例し、 $x _ 2$ が出る確率は、 $\displaystyle{ f(x _ 2| \theta) }$ に比例します。ですので、すべてのデータが得られる確率は、それらの確率密度の値の積に比例します。これを 尤度関数（Likelihood Function） といい、次の式で表されます。

$\displaystyle{ L(\theta) = \prod_{i=1}^{n} f(x_i | \theta) }$

ここで、 $\displaystyle{ x _ i }$ が観測された値で、 $\displaystyle{ \theta }$ がパラメタでです。

　 $\displaystyle{ f(x _ i | \theta) }$ は、もともと、パラメタが既知の $\displaystyle{ \theta }$ の下で、 $\displaystyle{ x _ i }$ が観測される確率密度を表す式です。

　しかし、尤度 $\displaystyle{ L(\theta) }$ の計算では、データ $\displaystyle{ x _ 1, x _ 2, \dots, x _ n }$ が既に観測されていますので、 $\displaystyle{ x _ i }$ が既知で、 $\displaystyle{ \theta }$ が未知です。この意味を踏まえれば、 $\displaystyle{ f(\theta|x _ i) }$ とするべきです。

　この尤度関数 $\displaystyle{ L(\theta) }$ が最も大きくなるような $\displaystyle{ \theta }$ を求めるのが最尤法です。

尤度関数の対数をとる

　多くの場合、実際の計算では、尤度関数の対数をとった対数尤度関数（Log-Likelihood Function）を使います。対数をとることで、積の計算が和に変わり、計算がしやすくなります。

$\displaystyle{ \ell(\theta) = \log L(\theta) = \sum_{i=1}^{n} \log f(x_i | \theta) }$

この対数尤度関数を最大化する $\displaystyle{ \theta }$ を求めることで、最尤推定ができます。

最尤推定の例

　最尤推定を指数分布と正規分布のパラメタ推定に使ってみます。

指数分布の場合

指数分布とは？

　指数分布（Exponential Distribution）は、ある事象が発生するまでの時間をモデル化する確率分布です。たとえば、故障が起こるまでの時間や、電話の着信間隔の時間などを表すのに使われます。

　指数分布の確率密度関数は次のようになります（パラメータは $\displaystyle{ \lambda }$ ）：

$\displaystyle{ f(x | \lambda) = \lambda e^{-\lambda x}, \quad (x > 0, \lambda > 0) }$

ここで、

$\displaystyle{ x }$ は観測データ（例：到着時間間隔）
$\displaystyle{ \lambda }$ は分布のパラメータ（大きいほど短い間隔で発生）

尤度関数を計算してみる

データ $\displaystyle{ x _ 1, x _ 2, \dots, x _ n }$ が得られたとき、尤度関数は

$\displaystyle{ L(\lambda) = \prod_{i=1}^{n} \lambda e^{-\lambda x_i} }$

となります。この対数をとると、

$\displaystyle{ \ell(\lambda) = \sum_{i=1}^{n} \log \lambda - \lambda \sum_{i=1}^{n} x_i }$

です。これを $\lambda$ について最大化すると、最尤推定の式 (最尤推定量) $\displaystyle{ \hat{\lambda} }$ がえられます。具体的な計算は、下の付録を見てください。

　結果は、これです。

$\displaystyle{ \hat{\lambda} = \frac{n}{\sum_{i=1}^{n} x_i} }$

　つまり、データの平均値の逆数が最尤推定値になります。

　(心の声) なんだよ、ピアソンのモーメント法の結果と同じじゃねーの。ピアソン頑張れ！

　モーメント法の結果は、以下の記事を参照してください。

【Rで統計入門】(3) モーメント法を使った確率分布のパラメタ推定 - ケィオスの時系列解析メモランダム

正規分布の場合

正規分布とは？

　正規分布（Normal Distribution）は、データが平均値を中心に左右対称に分布するモデルです。確率密度関数は次のようになります（パラメータは $\displaystyle{ \mu }$ と $\displaystyle{ \sigma }$ ：

$\displaystyle{ f(x | \mu, \sigma) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right) }$

ここで、
- $\displaystyle{ \mu }$ は平均（データの中心） - $\displaystyle{ \sigma }$ は標準偏差（データの広がり）

尤度関数を計算してみる

尤度関数は

$\displaystyle{ L(\mu, \sigma) = \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{(x_i - \mu)^2}{2\sigma^2}\right) }$

です。対数をとると、

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

です。

　これを最大化するために偏微分して計算すると、最尤推定値がえらます。

　結果はこれです。

$\displaystyle{ \begin{align} \hat{\mu} &= \frac{1}{n} \sum_{i=1}^{n} x_i \\ \hat{\sigma}^2 &= \frac{1}{n} \sum_{i=1}^{n} (x_i - \hat{\mu})^2 \end{align} }$

つまり、最尤推定では「データの標本平均が $\displaystyle{ \mu }$ の推定値」、「標本分散が $\displaystyle{ \sigma^ 2 }$ 」という当たり前に感じられる結果になります。

　(心の声) おーい、また、ピアソンのモーメント法の結果と同じじゃねーか。ピアソンが怒るのも当然だ。

と、結果だけをみて、最尤法を評価してはいけません。最尤法には深～い意味があり、結果だけを見て、なんだか面倒という印象を持たないでください。最尤法の本質を理解するには良い推定とは何かを整理して、理解する必要があります。最尤推定量は、下の表にある良い性質を持つことが証明されています。これについては、今回は詳しく説明しません。

性質	意味	実用的なメリット
一致性	データが増えると真の値に収束	サンプルが多ければ信頼性が上がる
漸近正規性	大きなサンプルでは正規分布に近づく	誤差の推定・信頼区間の計算が容易
有効性	最も分散が小さい推定量	精度の高い推定ができる

Rで最尤推定しません

　Rスクリプトで、指数分布と正規分布の最尤推定を数値実験しようと思いましたが、つまらないので止めます。

最後に一言

　若いころに学んだ基礎は、意外と覚えているものです。でも、計算のひらめきや暗算のスピードは確実に衰えていきます。

　昔は、講義中に計算の過程を予習しなくても、過去に一度解いた計算は覚えていたので、その場でスラスラ解けました。でも今では、たとえ事前に予習していても、その予習中に、これどうやって計算したっけ？と悩みます。講義で、計算過程を説明するときは、とても緊張します。

　なぜかって？それは、うっかりミスが増えたから。

数字の書き間違い、符号の見落とし…… どれも些細なことですが、講義中のミスで余計な時間を使い、学生に申し訳ないなと感じます。昔の自分なら、こんな間違いはしなかったのに —— そう思うことも増えました。

　それでも、学び続けることに意味があると信じています。間違いながらでも、考え続けること。それこそが、知識を深め、成長し続けるための秘訣なのかもしれません。

　ピアソンには、権威の座にこだわらず、フィッシャーの成果を正当に評価できるような、大きな器の人物になってほしかったです。

付録：指数関数の尤度を最大化

　ここでは、指数分布のパラメータ $\displaystyle{ \lambda }$ の最尤推定量の導出を、詳しく説明します。

尤度関数 (Likelihood Function)

　指数分布の尤度関数は以下です。

$\displaystyle{ \begin{aligned} L(\lambda) &= \prod_{i=1}^n f(x_i \mid \lambda) \\ &= \prod_{i=1}^n \lambda e^{-\lambda x_i} \\ &= \lambda^n \exp\Bigl(-\lambda \sum_{i=1}^n x_i\Bigr). \end{aligned} }$

したがって、対数尤度関数 (Log-Likelihood) は、

$\displaystyle{ \ell(\lambda) = \log L(\lambda) = \log \Bigl[\lambda^n \exp\Bigl(-\lambda \sum_{i=1}^n x_i\Bigr)\Bigr]. }$

です。この式を展開すると、

$\displaystyle{ \ell(\lambda) = \log \Bigl(\lambda^n\Bigr) + \log \Bigl(\exp\bigl(-\lambda \sum_{i=1}^n x_i\bigr)\Bigr) = n \log \lambda - \lambda \sum_{i=1}^n x_i. }$

対数尤度関数の最大化するパラメタの計算

　最尤推定量を求めるには、対数尤度 $\displaystyle{ \ell(\lambda) }$ を最大にする $\displaystyle{ \lambda }$ を探します。まずは、極値を探すために微分し、それが0の点を探します。

$\displaystyle{ \begin{align} \frac{d\,\ell(\lambda)}{d \lambda} &= \frac{d}{d \lambda} \Bigl(n \log \lambda - \lambda \sum_{i=1}^n x_i\Bigr)\\ &= \frac{n}{\lambda} - \sum_{i=1}^n x_i \end{align} }$

これを 0 にすると、

$\displaystyle{ \frac{n}{\lambda} - \sum_{i=1}^n x_i = 0 \quad \Longrightarrow \quad \frac{n}{\lambda} = \sum_{i=1}^n x_i \quad \Longrightarrow \quad \lambda = \frac{n}{\sum_{i=1}^n x_i} }$

　微分係数を 0 にするだけでは、その点が最大か最小か分かりません。そこで、二次微分係数を計算して、それが最大値であることを確認します。

　二次微分係数を計算すると、

$\displaystyle{ \frac{d^2\,\ell(\lambda)}{d \lambda^2} = \frac{d}{d \lambda} \Bigl(\frac{n}{\lambda} - \sum_{i=1}^n x_i\Bigr) = -\frac{n}{\lambda^2} < 0 }$

ですので、第二次微分が負ということは、対数尤度関数がその点で上に凸になっていることを意味します。

　つまり、先ほど求めた点は極大値になっていることが分かります。

　尤度の端点と無限大での挙動についても調べておく必要があります。

　 $\displaystyle{ \lambda \to 0^ + }$ のとき、 $\displaystyle{ \lambda^ n \to 0 }$ 、かつ、 $\displaystyle{ \exp \left(-\lambda \sum x _ i \right) \to 1 }$ となり、