本記事は、データ解析のための統計モデリング入門（通称緑本）の2章を参考にしています。
（もし内容に不備等あればご一報ください）
R,Pythonコード、使用しているデータはこちら
最尤推定法とは
「あてはあまりの良さ」を表す統計量（尤度）を最大にするようなパラメーターの値を探すパラメータ推定法
（尤度とは、あるパラメータ $\theta$ をパラメータとする確率分布から観測データ $y_i$ が発生する確率 $p(y_i|\theta)$ の積で定義される）
と書いてもよく分からないので、例題にそって最尤推定法を実施してみる。
例題：種子の統計モデリング
目的
以下のデータを解析する。
- 架空の植物50個体からなる集団を調査する。
- 各個体の種子数を数えたものがデータだとする。
植物A：種子数2
植物B：種子数3
植物C：種子数6
⋮
（合計50個体）
植物50個体それぞれの種子の数 $\{y_i\}$
{2, 2, 4, 6, 4, 5, 2, 3, 1, 2, 0, 4, 3, 3, 3, 3, 4, 2, 7, 2, 4, 3, 3, 3, 4,
3, 7, 5, 3, 1, 7, 6, 4, 6, 5, 2, 4, 7, 2, 2, 6, 2, 4, 5, 4, 5, 1, 3, 2, 3}
この種子数データを統計モデルとして表現するために、ふさわしい確率分布を考える
データの観察
データの分布（種子数のとグラム）を図示してみると以下の図のようになる。
Rコード
# データから読み込み
load("data.RData")
data <- c(2, 4, 6, 4, 5, 2, 3, 1, 2, 0, 4, 3, 3, 3, 3, 4, 2, 7, 2, 4,
3, 3, 3, 4, 3, 7, 5, 3, 1, 7, 6, 4, 6, 5, 2, 4, 7, 2, 2, 6,
2, 4, 5, 4, 5, 1, 3, 2, 3)
#ヒストグラム
hist(data, breaks = seq(-0.5, 9.5, 1))pythonコード
# データの読み込み
data = pd.read_csv(r'./data.csv', header=None)
# ヒストグラム
plt.hist(data, bins=np.arange(-0.5, 8.5, 1))ヒストグラムやデータの背景から以下の情報が読み取れる。
- データの値は種子数 $y_i$ は1個,2個,･･･と数えられる非負の整数である（カウントデータである）
- $y_i$ に下限（ゼロ）はあるみたいだが、上限があるか不明
- この観測データでは種子数の標本平均と標本分散がだいたい等しい（3.56と2.98）
- 分布の形状（一山のグラフになっている）
以上の条件より、種子数データを統計モデルとして表現するための確率分布としてポアソン分布が適切だと予想できる。
ポアソン分布について（補足）
ちなみに、ポアソン分布は以下の式で定義される。
$$ p(y|\lambda) = \frac{\lambda \exp(-\lambda)}{y!}$$
ここの $p(y|\lambda)$ は、平均が $\lambda$ である時に、ポアソン分布に従う確率変数が ある値 $y$ になる確率を示す。
例えば、平均( $\lambda$ ) 3.0 （分散も3.0）のポアソン分布を描くと以下のようになり、
この場合、種子数が4($y=4$)になる確率は、下記式で 0.168と求められる。
$$ p(4|3) = \frac{3.0^4 \exp(-3.0)}{4!} \approx 0.168 $$
最尤推定法でパラメータを推定する（本題）
最尤推定法とは
「あてはあまりの良さ」を表す統計量（尤度）を最大にするようなパラメーターの値を探すパラメータ推定法
（今回の場合｜尤度が最大となるような $\lambda$ の値を推定する）
尤度とは
あるパラメータ $\theta$ をパラメータとする確率分布から観測データ $y_i$ が発生する確率 $p(y_i|\theta)$ の積のこと
（「どういうパラメータを設定すると、得られたデータをうまく再現できるか」というイメージ）
$$ L(\theta|\boldsymbol{Y}) = \prod_{i}p(y_i|\theta) $$
通常は計算上の都合などから、尤度を対数変換した対数尤度がよく用いられる。
$$ \log L(\theta|\boldsymbol{Y}) = \sum_{i} p(y_i|\theta)$$
この対数尤度が最大になるような $\hat{\theta}$ を探すのが最尤推定
例題）
下記のような 3個体分の種子数データが得られた場合の尤度は $\lambda$=3、 $\lambda$=6の場合でそれぞれ次のような値となる。
$$ {y_1, y_2, y_3} = {2,2,4} $$
尤度 = ($y_1=2$ となる確率) × ($y_2=2$ となる確率) × ($y_3=4$ となる確率)なので
$\lambda = 3$ の場合
$$
\begin{align}
{L(\lambda=3|\boldsymbol{Y})}&=
L(\lambda=3|y_1=2) \times L(\lambda=3|y_2=2) \times L(\lambda=3|y_3=4)\\
&=
0.224 \times 0.224 \times 0.168 \approx 0.0084
\end{align}
$$
$\lambda = 6$ の場合
$$
\begin{align}
{L(\lambda=6|\boldsymbol{Y})}&=
L(\lambda=6|y_1=2) \times L(\lambda=6|y_2=2) \times L(\lambda=6|y_3=4)\\
&= 0.045 \times 0.045 \times 0.134 \approx 0.0003
\end{align}
$$
$\lambda = 3$ の場合と$\lambda = 6$ の場合を比べると前者の方が尤度が大きいので、よりデータをうまく再現できていると言える。
Rコード
# 平均3,6のポアソン分布を図示してみる。
y <- 0:9
prob <- dpois(y, lambda = 3)
plot(y, prob, type = "b", lty = 3)
prob <- dpois(y, lambda = 6)
plot(y, prob, type = "b", lty = 3)Pythonコード
# 平均3.56のポアソン分布を図示してみる。
y = range(10)
prob = poisson.pmf(y, mu=3.56)
plt.plot(y, prob, "o--")また、$\lambda$ と$\log{L}$ の関係を図示すると以下のようになり、対数尤度が最大になるのはおおよそ3.5付近だと予想することができる。
Rコード
# 対数尤度LogL(lambda)とlambdaの関係を調べる
logL <- function(m) sum(dpois(data, m, log=TRUE))
lambda <- seq(2, 5, 0.1)
plot(lambda, sapply(lambda, logL), type= "l")Python
# 対数尤度 logL(lambda)とlambdaの関係を調べる
logL = lambda lambda_: sum(poisson.logpmf(data,mu=lambda_))
lambda_list = np.arange(2, 5, 0.1)
logL_list = [logL(lambda_) for lambda_ in lambda_list]
plt.plot(x, logL_list)解析的に最尤推定量を求める
対数尤度関数が最大となるのは対数尤度関数の傾きが0になる点なので、関数を $\lambda$ で微分することで、最尤推定量が求められる。
$$ \frac{\partial \log{L}(\lambda)}{\partial\lambda} = \sum_i{\frac{y!}{\lambda}} = \frac{1}{\lambda}\sum_iyi - 50$$
$$\hat{\lambda}=\frac{1}{50}\sum_i{y_i}=\frac{\sum_iy_i}{50} = (データの標本平均) = 3.56$$
ゆえに、最尤推定量値 $\hat{\lambda}$ は3.56と求められる
補足（ポアソン分布の尤度関数）
$$ L(\lambda|y_i) = \prod_{i}\frac{\lambda^{y_i}\exp(-\lambda)}{y_{i}!} $$
$$ \log L(\lambda|y_i) = \sum_{i}(y_i\log{\lambda}-\lambda-\sum_{k}^{y_i}\log{k})$$
実際のデータ解析で使うモデルはより複雑なので、今回のように解析的に（式の変形で）最尤推定値を求めるのは難しい。
そのため、通常は数値的な試行錯誤により、対数尤度ができるだけ大きくなるパラメータを探す。
以上
参考
書籍
