最尤推定

1. はじめに

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

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

2. 主要な概念の整理

2.1 尤度関数と最尤推定量

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

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

数学的には、最尤推定量 $\widehat{\boldsymbol{\theta}}$ は以下のように定義されます。

$$ L(\widehat{\boldsymbol{\theta}}) \geq L(\boldsymbol{\theta}), \quad \boldsymbol{\theta} \in \Omega $$

ここで、$\Omega$ はパラメータ空間（パラメータ $\boldsymbol{\theta}$ がとり得る全ての値の集合）を表します。

2.2 スコア関数

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

$$ \boldsymbol{U}(\boldsymbol{\theta})=\frac{\partial l(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} $$

ここで、$l(\boldsymbol{\theta}) = \log[L(\boldsymbol{\theta})]$ は対数尤度関数です。

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

2.3 フィッシャー情報行列

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

$$ \boldsymbol{I}(\boldsymbol{\theta})=\mathrm{E}\left[\boldsymbol{U}(\boldsymbol{\theta}) \boldsymbol{U}^T(\boldsymbol{\theta})\right] $$

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

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

3. 理論

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

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

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

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

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

4. 具体例：二項分布

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

二項分布の尤度関数は

$$ L(\pi) = \binom{N}{y} \pi^y (1-\pi)^{N-y} $$となります。

ここで、$N$ は試行回数、$y$ は成功回数、$\pi$ は成功確率、$\binom{n}{k}$は二項係数です。

パラメータが 1 個のとき、最尤推定量の漸近正規性から$\begin{aligned}\operatorname{Var}(\widehat{\theta}) \approx \frac{1}{I(\theta)}\end{aligned}$ですから、標準誤差は$\begin{aligned}\mathrm{SE}(\widehat{\theta})=\left[-\left.\frac{d^2 l(\theta)}{d \theta^2}\right|_{\widehat{\theta}}\right]^{-1 / 2}\end{aligned}$となります。
これを 2 項尤度$l(\pi)=Y \log (\pi)+(N-y) \log (1-\pi)$に適用します。

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

$$ \frac{d l(\pi)}{d \pi}=\frac{y}{\pi}-\frac{N-y}{1-\pi} $$

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

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. 計算の複雑さ：複雑なモデルでは、尤度関数の最大化が計算上困難になることがある

参考文献

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