統計モデルのパラメーターを決定する際の一手段として最尤推定法を使う事が多々あるが、この手法だとモデルを複雑化する（≒パラメーターを増やす）だけでも、与えられたデータに対する最大対数尤度を改善することは可能である。しかし、これはたまたま与えられたデータへのあてはまり度合を高めているだけであって、モデル自身の予測性能を向上させているわけではない。
統計モデルとしての予測の良さ、すなわち「観測データで推定・構築されたモデルが、また新たに取得してきたデータに対してどの程度正確にあてはまるのか」を評価する量は平均対数尤度\[ E\left[ \log(L)\right] \]である。
ここで\( bias \)という量を平均対数尤度と最大対数尤度を使って
\[ bias := \log(L^{*}) - E\left[ \log(L)\right] \]
と定義すると、\( bias \)の推定量はモデルパラメーターの個数(\( k \)と書く)に等しくなることが数理統計学により示されている。従って、平均対数尤度を最大化するようにモデルを構築したい場合には平均対数尤度の代わりに\[ \log(L^{*}) - k \]を最大化すればよく、実はこれがAIC（赤池情報量規準）でモデルの良さを評価した統計モデル構築法となっている。
前置きが随分と長くなってしまったが、ここではより単純な話で\( bias \)の推定量はきちんとモデルパラメーター数になっているのか？という点を確かめてみたいと思う。
手順は以下の通り。ここではポアソン分布の強度をモデルパラメーターとしている。従って、モデルパラメーターは１となる。
まずはbiasを計算する関数を書く。こんな感じ。（１の内容）
bias <- function(lambda.true, sample.size){
#１サンプルセット（sample.size個）のデータを生成し、ポアソン分布の強度推定
sample.rpois <- rpois(sample.size, lambda.true)
fit <- glm(sample.rpois~1, family=poisson)
#glm推定結果からモデルパラメーター（ポアソン分布の強度（推定））を計算
#(モデル：log(lambda) = beta)
lambda.estimated <- exp(coef(fit))
#また別に本物のパラメーター(lambda.true)から
#サンプルセット（sample.size個）を200セットサンプリング
#平均対数尤度をポアソン分布の強度（推定）から計算
likelihood.mean <- mean(sapply(1:200, function(i){sum(log(dpois(rpois(N, lambda.true), lambda.estimated)))}))
#bias(最大尤度-平均尤度）を返却
logLik(fit) - likelihood.mean
}
次にこれを使ってbaisを計算、平均値も算出。（２＆３の内容）
ちゃんと\( bias \)の推定値がモデルパラメーター数（１）に近い値をとっていることがわかる。
# 各種パラメーター設定 １サンプルセットに含まれるサンプル数
N <- 50
# 真のポアソン分布の強度
lambda.true <- 8
# 1000回biasの推定を繰り返す
bias.sampled <- sapply(1:1000, function(i) bias(lambda.true, N))
mean(bias.sampled)
## [1] 1.053
分布も合わせてみてみる。（ggplot2だとこんなに簡単に美しくPLOTできる！）
library(ggplot2)
qplot(bias.sampled, geom = "blank") + geom_histogram(aes(y = ..density..),
colour = "black", fill = "white") + geom_density(alpha = 0.2, fill = "#6666FF") +
geom_vline(aes(xintercept = mean(bias.sampled)), color = "red", linetype = "dashed",
size = 2)