推定量と分散について

投稿日: 作成者:

統計的な推定問題における推定量の分布について考えてみたいと思います。

例として \(\{X_1,\; \cdots,\; X_n\}\) という独立な指数分布からのサンプルを考えます。
指数分布とは \(\lambda>0\) というパラメータを持つ確率分布で密度関数は

\begin{align*}
f(x|\lambda)=\lambda\exp(-\lambda x)
\end{align*}
のようなものです。したがって尤度関数 \(L(\lambda)\) は

\begin{align*}
L(\lambda)&=\prod_{i=1}^n\lambda\exp(-\lambda X_i)\\
&=\lambda^n\exp\bigg(-\lambda\sum_{i=1}^n X_i\bigg)
\end{align*}
となりますが、n乗根 \(L(\lambda)^{1/n}\) を最大化しても同じことなので，新しい尤度関数 \(\hat{L}(\lambda):=L(\lambda)^{1/n}\) は

\begin{align*}
\hat{L}(\lambda)&=\lambda\exp\bigg(-\lambda\cdot\frac{1}{n}\sum_{i=1}^n X_i\bigg)\\
&=\lambda\exp(-\lambda\bar X)
\end{align*}
となります。このことは平均値さえ知っていれば最尤推定ができることを意味しています。つまり、
\begin{align*}\hat\lambda=1/\bar X\end{align*}という推定量が得られます。

ほかの推定方法として、q-分位点 T の統計量（たとえば中央値）より \(\lambda\) を使うことでも推定可能です：
\begin{align*}
&q=1-F(\tilde T)=1-\exp(-\lambda\tilde T)\\
\Leftrightarrow \quad&
\hat\lambda=-\frac1{\tilde T}\log(1-q)
\end{align*}これも先ほどの平均値から得られる推定量とは別のものですが、わりと自然な考え方で推定は可能な気がします。
この推定量をここでは 1/q-分位推定量とよぶことにしましょう。

ここで

  • 最尤推定量
  • 二分位推定量
  • 四分位推定量
  • 八分位推定量

の性能を比較するために簡単な実験をしてみましょう。


それぞれの推定量に対して、

  • 100個の独立なλ=2の指数分布の乱数を発生させる
  • 推定量を計算
  • という操作を10000回繰り返して推定量の分布を生成する

という手順を通して、推定量を比較してみます。

結果は次の図のようになりました。
EstiDist.jpg
すべての推定量は真のパラメータ 2 の周囲に分布しますが、そのばらつきがまったく異なりますね。
推定量の分散はフィッシャー情報量の逆行列よりも小さくなることはない、というクラメル・ラオの定理がありますが、分位点による推定量は明らかに最尤推定量よりも大きな分散を持っているため、効率が悪いと考えることができるでしょう。

それでは最尤推定量の分散はフィッシャー情報量の逆数よりも大きいのか？それとも理論的な分散の下限であるフィッシャー情報量の逆数を達成しているのか？ということが知りたくなりますが、このブログでは理論的にではなく数値的にチェックしてみる、というアプローチをとってみます。

（次回へつづく）

library(MASS)
plot.new()
par(mfrow=c(4,1))
xlabel="estimator"
ylabel="density"
# estimator by mean
truehist(sapply(1:10000, function(x) 1/mean(rexp(100,2))),xlim=c(1,4),ylim=c(0,2),col="gray",xlab=xlabel,ylab=ylabel)
title("Distribution of Maximum likelihood estimator")
# estimator by median
truehist(sapply(1:10000, function(x) -1/median(rexp(100,2))*log(0.5)),xlim=c(1,4),ylim=c(0,2),col="gray",xlab=xlabel,ylab=ylabel)
title("Distribution of median estimator")
# estimator by 1/4-quantile
truehist(sapply(1:10000, function(x) -1/quantile(rexp(100,2),0.25)*log(0.75)),xlim=c(1,4),ylim=c(0,2),col="gray",xlab=xlabel,ylab=ylabel)
title("Distribution of 1/4-quantile estimator")
# estimator by 1/8-quantile
truehist(sapply(1:10000, function(x) -1/quantile(rexp(100,2),0.125)*log(0.875)),xlim=c(1,4),ylim=c(0,2),col="gray",xlab=xlabel,ylab=ylabel)
title("Distribution of 1/8-quantile estimator")