読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

指数分布の標本平均と不偏推定量

R 確率分布

命題

指数分布(密度関数は f_X(x)=\lambda e^{-\lambda x})のパラメータ λ の m 個の標本平均  \bar{X}_j =\sum^n_{i=1}X_i/n をもとに, パラメータ λ の不偏推定量を求める.

\displaystyle \hat{\lambda} = \frac{1}{m}\sum^m_{j=1}\frac{n-1}{n} \frac{1}{\bar X_j}
これは不偏.

証明

まず  E[1/\bar X] を求める.

UMVUE(一様最小分散不偏推定量)について - BIGLOBEなんでも相談室 の解法は教育的だと思う.

平均 \theta=1/\lambda の指数分布(密度関数は f_X(x)= \lambda e^{-\lambda x})に従う確率変数の和は, 形状パラメータ n, 尺度パラメータ 1/λ のガンマ分布に従う.

すなわち,
 Y = \sum _ {i=1}^n X_i確率密度関数  f_Y(y) は、
 \displaystyle f_Y(y) = \frac{\lambda^{n}}{\Gamma (n)} y^{n-1} e^{-\lambda y}.

これより k 次モーメントを求める.
\displaystyle E[Y^k] = \int_0^{\infty} y^k f_Y(y) \,dy \\
\displaystyle = \int_0^{\infty} \frac{\lambda^{-k}}{\Gamma (n)} (\lambda y)^{n+k-1} e^{-\lambda y}\lambda \,dy

\lambda y = t と置き,

\displaystyle E[Y^k] = \int_0^{\infty} {\lambda^{-k}} \frac{1}{\Gamma(n)} t^{n+k-1} e^{-t} \,dt\\
\displaystyle= \lambda^{-k} \frac{1}{\Gamma(n)} \int_0^{\infty} t^{n+k-1} e^{-t} \,dy\\
\displaystyle= \lambda^{-k} \frac{1}{\Gamma(n)} \Gamma(n+k)

上式で k=-1 とすると, E[Y^{-1}]=\lambda\frac{1}{(n-1)} より,

\displaystyle E[1/\bar X]=E[nY^{-1}] =\frac{n}{n-1}\lambda

結果,
\displaystyle E[ \hat{\lambda} ]= E\left[\frac{1}{m}\sum^m_{j=1}\frac{n-1}{n} \frac{1}{\bar X_j}\right] \\
\displaystyle= \frac{1}{m}\sum^m_{j=1}\frac{n-1}{n} \frac{n}{n-1}\lambda \\
\displaystyle=\lambda.

シミュレーション

平均の逆数 \displaystyle \frac{1}{\bar X} を取った場合の平均 result1 と, \displaystyle \frac{n-1}{n} \frac{1}{\bar X_j} = \frac{n-1}{\sum_{i}^{n}X_i} の平均 result2 を比較する. 試行回数は10000回.

f:id:abrahamcow:20140720033747p:plain
赤い点線はシミュレーションで仮定した真値. result1 はバイアスが生じていることがわかる.

R のコード:

res1 = numeric(m)
res2 = numeric(m)
for(j in 1:10000){
n=100
m=10
ans1 = numeric(m)
ans2 = numeric(m)
for(i in 1:m){
  ans1[i]=1/mean(rexp(n,3))
  ans2[i]=(n-1)/sum(rexp(n,3))
}
res1[j]=mean(ans1)
res2[j]=mean(ans2)
}

#グラフ
library(ggplot2)
library(reshape2)
result <- melt(data.frame(result1 =res1,result2 =res2))
ggplot(result, aes(x=variable,y=value)) + 
  geom_boxplot() +
  geom_abline(intercept = 3, slope = 0, colour=2,  linetype = 2)