推定量と分散について 2 (クラメル・ラオの不等式)

  • ????????????????????

先日、指数分布のパラメータλのさまざまな推定方法に対して数値実験をしました(前回のエントリー:推定量と分散について)。
100個の独立な確率変数から推定を10000回行った結果、それぞれの推定量に対して次のような分布が得られました(推定量は確率変数の関数なので、やはり確率変数となり、広がりのある分布を持ちます)。
EstiDist.jpg
それぞれの推定量の意味は前回のエントリーを見てください。

グラフをみてすぐにわかることは、明らかに最尤推定量(一番上のグラフ)のバラつきが一番小さくなっていることがわかると思います(ちなみに真のパラメータはλ=2です)。

一方、統計学の有名な定理で クラメル・ラオの不等式 というものがあります。

【定理】クラメル・ラオ

正則統計モデルの n 個のサンプルを使った、不偏推定量\(\theta_n\)の分散に関して以下の不等式が成立する。
\begin{align*}V[\theta_n]\ge\{nI(\theta)\}^{-1}\end{align*}
ただし、\(I(\theta)\) はフィッシャー情報量で以下のように定義される。
\begin{align*}I(\theta)&=E\bigg[\frac{\partial}{\partial \theta}\log p(X|\theta)\bigg]^2
=-E\bigg[\frac{\partial^2}{\partial \theta^2}\log p(X|\theta)\bigg]\end{align*}

二番目の等号は微積の交換性などを仮定すれば導出できます。また、パラメータ空間が二次元以上でもこの式は成立します。その場合、不等号は半正定値の意味になります。

【定義】有効推定量

不偏推定量 \(\theta_n\) の分散がクラメル・ラオ不等式の下限を達成する(等号が成立する)とき有効推定量という。
ということで、先ほどのグラフの4つの推定量のなかで、有効推定量になっているものがあるとするならば、一番上の最尤推定量ということになるでしょう(本当はすべての推定量について不偏性があるかチェックしなければいけないのですが…)。


さて、指数分布に対してフィッシャー情報量を計算してみましょう。フィッシャー情報量の二回微分のほうを使うと、簡単に計算できて
\begin{align*}
I(\theta)&=\displaystyle\int_0^\infty \frac1{\lambda^2}p(x|\lambda)dx=\frac1{\lambda^2}
\end{align*}
となります。つまり、推定量の分散の下限は
\begin{align*}\frac{\lambda^2}{n}\end{align*}となります。

ここで、この理論的な値と、4つの推定量の分散を比較してみましょう。
分散はサンプル数に依存するので、n=100, 200, 300, 400, 500 とします。
各推定量に対して10000回ずつ、推定を行って、その標本分散を計算したものが次のグラフです。
CramerRao2.png

最尤推定量が理論的な分散の下限とほぼ一致していることがわかります。その他の推定量は下限を達成していません。注意としては、一般に最尤推定量はいつでも有効推定量になるわけではないということです。このあたりについては稲垣宣生「数理統計学」に証明も含めて書いてあります。

数理統計学 (数学シリーズ) 数理統計学 (数学シリーズ)
稲垣 宣生

裳華房
売り上げランキング : 366872

Amazonで詳しく見る

最後に上のグラフを計算するための R で書かれたスクリプトを示しておきます。

# 指数分布の推定量の分散とクラメル・ラオの下限を比較する R スクリプト
plot.new()
par(mfrow=c(1,2))
lambda=2
K=c(100,200,300,400,500)
thr.var = lambda^2/K
mle.var = sapply( K, function(k){
var(sapply(1:10000, function(x) 1/mean(rexp(k,lambda)))) } )
med.var = sapply( K, function(k){ var(sapply(1:10000,
function(x) -1/median(rexp(k,lambda))*log(0.5))) } )
q4.var  = sapply( K, function(k){ var(sapply(1:10000,
function(x) -1/quantile(rexp(k,lambda),0.25)*log(0.75))) } )
q8.var  = sapply( K, function(k){ var(sapply(1:10000,
function(x) -1/quantile(rexp(k,lambda),0.125)*log(0.875))) } )
plot.new()
par(mfrow=c(1,2))
plot(K,thr.var,ylim=c(0.,0.4),type="o",
xlab="# of samples used",ylab="Variance of estimator")
points(K,mle.var,pch=2,type="o",col=2)
points(K,med.var,pch=3,type="o",col=3)
points(K,q4.var,pch=4,type="o",col=4)
points(K,q8.var,pch=5,type="o",col=5)
legend("topright",legend=c("theoretical lower bound","maximum likelihood",
"median","1/4-quantile","1/8-quantile"),col=1:5,pch=1:5,lt=1,cex=0.8)
title("Relationship between # of samples \nand variance of estimators")
はてなブックマーク - 推定量と分散について 2 (クラメル・ラオの不等式)
Pocket