最尤推定について
本記事は、データ解析のための統計モデリング入門(通称緑本)の2章を参考にしています。
(もし内容に不備等あればご一報ください)
R,Pythonコード、使用しているデータはこちら
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
最尤推定法とは
「あてはあまりの良さ」を表す統計量(尤度)を最大にするようなパラメーターの値を探すパラメータ推定法 (尤度とは、あるパラメータ
と書いてもよく分からないので、例題にそって最尤推定法を実施してみる。
例題:種子の統計モデリング
目的
以下のデータを解析する。
- 架空の植物50個体からなる集団を調査する。
- 各個体の種子数を数えたものがデータだとする。
植物A:種子数2
植物B:種子数3
植物C:種子数6
⋮
(合計50個体)
植物50個体それぞれの種子の数
{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} \{y_i\}
この種子数データを統計モデルとして表現するために、ふさわしい確率分布を考える
データの観察
データの分布(種子数のとグラム)を図示してみると以下の図のようになる。
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))
ヒストグラムやデータの背景から以下の情報が読み取れる。
- データの値は種子数
は1個,2個,・・・と数えられる非負の整数である(カウントデータである)y_i -
に下限(ゼロ)はあるみたいだが、上限があるか不明y_i - この観測データでは種子数の標本平均と標本分散がだいたい等しい(3.56と2.98)
- 分布の形状(一山のグラフになっている)
以上の条件より、種子数データを統計モデルとして表現するための確率分布としてポアソン分布が適切だと予想できる。
ポアソン分布について(補足)
ちなみに、ポアソン分布は以下の式で定義される。
ここの
例えば、平均(
最尤推定法でパラメータを推定する(本題)
最尤推定法とは
「あてはあまりの良さ」を表す統計量(尤度)を最大にするようなパラメーターの値を探すパラメータ推定法 (今回の場合|尤度が最大となるような
尤度とは
あるパラメータ
(「どういうパラメータを設定すると、得られたデータをうまく再現できるか」というイメージ)
通常は計算上の都合などから、尤度を対数変換した対数尤度がよく用いられる。
この対数尤度が最大になるような
例題)
下記のような 3個体分の種子数データが得られた場合の尤度は
尤度 = (
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--")
また、
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になる点なので、関数を
ゆえに、最尤推定量値
補足(ポアソン分布の尤度関数)
実際のデータ解析で使うモデルはより複雑なので、今回のように解析的に(式の変形で)最尤推定値を求めるのは難しい。 そのため、通常は数値的な試行錯誤により、対数尤度ができるだけ大きくなるパラメータを探す。
以上
参考
書籍
- データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- StanとRでベイズ統計モデリング (Wonderful R)
- Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド
以上
Discussion