最尤推定


1. はじめに

最尤推定法はデータから母集団のパラメータを推定する強力な手法です。
例えばコインを投げて表が出る確率を推定したいとします。
10回投げて7回表が出たとしたら、表が出る確率は7/10=0.7だと推測するでしょう。
これが最尤推定の基本的な考え方です。

しかし現実の問題はもっと複雑です。
多くのパラメータを同時に推定する必要があったり、データの分布が複雑だったりします。
そんな時に威力を発揮するのが以下の概念です。

2. 主要な概念の整理

2.1 尤度関数と最尤推定量

まず尤度関数(Likelihood function)について説明します。

尤度関数 L(θ) は、パラメータ θ が与えられたときに観測されたデータが得られる確率(密度)を表します。
最尤推定法は、この尤度関数を最大にするパラメータ θ を見つける方法です。

数学的には、最尤推定量 θ^ は以下のように定義されます。

L(θ^)L(θ),θΩ

ここで、Ω はパラメータ空間(パラメータ θ がとり得る全ての値の集合)を表します。

2.2 スコア関数

スコア関数(Score function)は、対数尤度関数の偏微分として定義されます。

U(θ)=l(θ)θ

ここで、l(θ)=log[L(θ)] は対数尤度関数です。

スコア関数は、パラメータの推定値がどの程度良いのかを示す指標として解釈できます。
スコア関数がゼロに近いほど、そのパラメータ値は尤度を最大化する値に近いと考えられます。

2.3 フィッシャー情報行列

フィッシャー情報行列(Fisher information matrix)は、スコア関数の分散共分散行列として定義されます。

I(θ)=E[U(θ)UT(θ)]

フィッシャー情報行列は、パラメータの推定精度に関する情報を表します。
直観的にはフィッシャー情報行列の値が大きいほど、そのパラメータの推定精度が高いと解釈できます。

スコア関数は確率変数ですが, Fisher 情報量はパラメータが決まれば定数になります.

3. 理論

3.1 最尤推定量の漸近正規性

最尤推定量の重要な性質の一つに漸近正規性があります。
これは、サンプルサイズが大きくなるにつれ最尤推定量の分布が正規分布に近づく性質です。

サンプルサイズ N が十分に大きいとき、最尤推定量 θ^ は以下の分布に収束します:
θ^dN(θ,I1)
ここで、d は分布収束を表し、I はフィッシャー情報行列です。

【証明】
確率変数Yi は独立な確率分布(同一でなくてもよい)に従うとします (i=1,,N)。
パラメータ θ を適当な推定方程式i=1NM(yi;θ)=0の解 θ^ として求めることを考えます。
A=E[M(Yi;θ)θ|θ0], B=E[M(Yi;θ0)M(Yi;θ0)T]とします。
パラメータの値 θ0E[M(Yi;θ0)]=0を満たす値として定義し、M¯(θ)=1Ni=1NM(yi;θ)とします。

テイラー展開により、M¯(θ^)θ0 の周りで展開すると、
M¯(θ^)=M¯(θ0)+M¯(θ0)(θ^θ0)+remainder
ここで、M¯(θ^)=0 (推定方程式の定義より)なので、
0=M¯(θ0)+M¯(θ0)(θ^θ0)+remainder

M¯(θ) が逆行列 M¯1(θ) を持つと仮定すると、N(θ^θ0)=M¯1(θ0)NM¯(θ0)+ remainder が得られます。

M¯(θ0)=1Ni=1N[θM(yi;θ)|θ0] は各項の平均です。
大数の法則により、サンプルサイズ N が大きくなるにつれて、この平均は期待値 A=E[M(Yi;θ)θ|θ0]に確率収束します。
M¯(θ0)=1Ni=1NM(yi;θ0) も各項の平均です。
中心極限定理により、NM¯(θ0) は正規分布に収束します。
その平均は 0E[M(Yi;θ0)]=0 より)、分散共分散行列は B=E[M(Yi;θ0)M(Yi;θ0)T] となります。
ゆえにN(θ^θ0) は確率収束する統計量と正規分布に収束する統計量の線型結合となります。
このような確率変数は, Slutskyの定理によって正規分布に分布収束します。

適当な正則条件の下で剰余項は 0 に収束します。
これらの結果を合わせることで、θ^dN[θ0,N1A1B(A1)T]という漸近正規性が成り立ちます。
もし θ0 がパラメータの真値であるとき、 M(yi;θ)不偏な推定関数(unbiased estimating equation)といいます。

最尤推定量はスコア方程式1Ni=1NU(yi;θ^)=0の解です。
スコア関数は不偏な推定関数ですので、上記の議論からθ^dN[θ,A1B(A1)T] となります。
最尤推定量の場合は特殊なケースで, 定義からB は Fisher 情報行列Iそのものです。
また, A は対象 i のスコア U(yi;θ) を微分してマイナスをとったものですから, 対象1つあたりの Fisher 情報行列 N1I に一致します。
つまり NA=B=I が成り立つので、 最尤推定量の漸近正規性θ^dN(θ,I1)が示されます。
Slutskyの定理
サンプルサイズ N のデータから計算された統計量 LNMN があったとします。
LNN が無限大のとき λ に確率収束するとします。
また, MNN(MNμ)dN(0,σ2)という分布収束を满たすとします。
このとき, LNMN の和と積は以下のように収束します.
N(LN+MNλμ)dN(0,σ2)N(LNMNλμ)dN(0,λ2σ2)

この定理は、最尤推定量が一致性(真の値に収束する性質)と漸近有効性(サンプルサイズが大きい時、漸近的に最小分散にむかう性質)を持つことを示しています。

3.2 スコア関数、フィッシャー情報行列、観測情報行列の関係

漸近正規性という便利な性質のおかげで、われわれは 正規分布を使って θ^ の信頼区間やp値を計算することができます。
しかし実際にはFisher 情報行列を求めなければなりません。
Fisher 情報行列の推定においてよく用いられるのは観測情報行列(observed information matrix)で代替する方法です.
観測情報行列 i(θ)i(θ)=2l(θ)θθTと定義されます.

対数尤度関数について, 微分と積分が交換可能と仮定します。
このとき, スコア関数 U(θ) の期待値はゼロ、つまりE[U(θ)]=0です。

また, スコア関数, Fisher 情報行列I(θ)=E[U(θ)U(θ)T],観測情報行列i(θ)=2l(θ)θθTには, 以下の関係が成り立ちます。
E[U(θ)θ]=I(θ)=E[i(θ)]

【証明】
簡単のため, パラメータが 1 個でデータが連続分布のときについてのみ考えます。

確率密度関数の性質からp(y;θ)dy=1です。
これを両辺について θ で微分します.
ddθp(y;θ)dy=d(1)dθ=0
ここで左辺は, 微分と積分が交換可能のときddθp(y;θ)dy=dp(y;θ)dθdyとなり, これはスコア関数の期待値なのでE[U(θ)]=0です。

次に, 上の式の左辺を微分すると, 微分と積分が交換可能であればdE[U(θ)]dθ=[ddθdl(θ)dθp(y;θ)]dy=[dl(θ)dθdp(y;θ)dθ+d2l(θ)dθ2p(y;θ)]dy=E[(dl(θ)dθ)2]+E[d2l(θ)dθ2]=I(θ)+E[dU(θ)dθ]
これはゼロなので、I(θ)=E[dU(θ)dθ]となります。

4. 具体例:二項分布

二項分布を例に、これらの概念を具体的に見てみます。

二項分布の尤度関数は

L(π)=(Ny)πy(1π)Nyとなります。

ここで、N は試行回数、y は成功回数、π は成功確率、(nk)は二項係数です。

パラメータが 1 個のとき、最尤推定量の漸近正規性からVar(θ^)1I(θ)ですから、標準誤差はSE(θ^)=[d2l(θ)dθ2|θ^]1/2となります。
これを 2 項尤度l(π)=Ylog(π)+(Ny)log(1π)に適用します。

l(π) を微分すると, 1 次の導関数は以下のようになります。

dl(π)dπ=yπNy1π

フィッシャー情報行列を求めるために2 次の導関数のマイナスをとると、d2l(π)dπ2=yπ2+Ny(1π)2が得られます。
これは π が含まれているので実際に計算できません。
最尤推定量でπを置き換えるとyπ^2+Ny(1π^)2=Nπ^(1π^)となります。
この逆数の平方根をとればSE=π^(1π^)Nが導かれます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom
import matplotlib as mpl

# 日本語フォントの設定
plt.rcParams['font.family'] = 'Yu Gothic'  

def plot_binomial_estimation(n, k):
    p_hat = k / n
    se = np.sqrt(p_hat * (1 - p_hat) / n)
    ci_lower = max(0, p_hat - 1.96 * se)
    ci_upper = min(1, p_hat + 1.96 * se)

    p_range = np.linspace(0, 1, 1000)
    likelihood = binom.pmf(k, n, p_range)
    
    plt.figure(figsize=(12, 6))
    plt.plot(p_range, likelihood, label='尤度関数')
    plt.axvline(p_hat, color='r', linestyle='--', label='最尤推定量')
    plt.axvline(ci_lower, color='g', linestyle=':', label='95%信頼区間')
    plt.axvline(ci_upper, color='g', linestyle=':')
    
    plt.xlabel('確率 p', fontsize=4)
    plt.ylabel('尤度', fontsize=4)
    plt.title(f'二項分布の最尤推定 (n={n}, k={k})', fontsize=4)
    plt.legend(fontsize=4)
    plt.grid(True, alpha=0.3)
    
    plt.text(0.5, plt.ylim()[1], f'推定値: {p_hat:.3f}, 95%CI: ({ci_lower:.3f}, {ci_upper:.3f})', 
             horizontalalignment='center', verticalalignment='bottom', fontsize=4)

    plt.show()

# 例:100回中60回表が出た場合
plot_binomial_estimation(100, 60)
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom
import matplotlib as mpl

# 日本語フォントの設定
plt.rcParams['font.family'] = 'Yu Gothic' 

def plot_binomial_estimation(n, k,fontsize):
    p_hat = k / n
    se = np.sqrt(p_hat * (1 - p_hat) / n)
    ci_lower = max(0, p_hat - 1.96 * se)
    ci_upper = min(1, p_hat + 1.96 * se)

    p_range = np.linspace(0, 1, 1000)
    likelihood = binom.pmf(k, n, p_range)
    
    plt.figure(figsize=(12, 6))
    plt.plot(p_range, likelihood, label='尤度関数')
    plt.axvline(p_hat, color='r', linestyle='--', label='最尤推定量')
    plt.axvline(ci_lower, color='g', linestyle=':', label='95%信頼区間')
    plt.axvline(ci_upper, color='g', linestyle=':')
    
    plt.xlabel('確率 p', fontsize=fontsize)
    plt.ylabel('尤度', fontsize=fontsize)
    plt.title(f'二項分布の最尤推定 (n={n}, k={k})', fontsize=fontsize)
    plt.legend(fontsize=fontsize)
    plt.grid(True, alpha=0.3)
    
    plt.text(0.5, plt.ylim()[1], f'推定値: {p_hat:.3f}, 95%CI: ({ci_lower:.3f}, {ci_upper:.3f})', 
             horizontalalignment='left', verticalalignment='center_baseline', fontsize=fontsize)
    plt.savefig('binomial_estimation.png', dpi=300, bbox_inches='tight')
    plt.close()

# 例:100回中60回表が出た場合
plot_binomial_estimation(100, 60,fontsize=12)

5. 限界

最尤推定法には以下のような限界や課題があります:

  1. モデルの誤特定:正しくないモデルを仮定すると、推定結果にバイアスが生じる
  2. 小標本での性能:サンプルサイズが小さい場合、漸近理論が適用できない
  3. 計算の複雑さ:複雑なモデルでは、尤度関数の最大化が計算上困難になることがある

参考文献

医学のための因果推論Ⅰ 一般化線形モデル

最終更新: