廿TT

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

指数分布に従う確率変数の和の分布の最尤推定

X をパラメータ 1/λ の指数分布に従う確率変数, Y をパラメータ 1/μ の指数分布に従う確率変数とすると, 和 Z=X+Y は密度関数

\displaystyle f_Z(z) = \frac{1}{\lambda-\mu}(e^{-\frac{1}{\lambda} z} - e^{-\frac{1}{\mu} z})

を持つ分布に従います(指数分布に従う確率変数の和の分布をたたみこみで求める - 廿TT).

この分布のパラメータ λμ最尤推定してみます.

和の分布の尤度方程式は解析的には解けないため, R の optim 関数を使って対数尤度を最大化するパラメータを数値的に求めることにします(手法 1).

一方, Z だけでなく X, Y の実現値も観測されている場合, λ, μ最尤推定量は X の標本平均, Y の標本平均です(手法 2).

1000回のシミュレーションの結果, やはり和 Z から求めた推定値はだいぶバラつきが大きくなってしまうようです.

f:id:abrahamcow:20150627044434p:plain

f:id:abrahamcow:20150627044443p:plain

箱ひげ図はそれぞれの手法による推定値、点線がシミュレーションで仮定した真値です。

なんとか改善の方法はないものでしょうか?

#R のコード
xbar <- numeric(1000)
ybar <- xbar
est1 <-  xbar
est2 <-  xbar

ll <- function(par){
  lambda <- par[1]
  mu <- par[2]
  out <- sum( log(1/(lambda-mu)) + log(exp(-z/lambda)-exp(-z/mu)) )
  ifelse(is.nan(out),-Inf,out)
}

for(i in 1:1000){
  x <- rexp(100, 1/8)
  y <- rexp(100, 1/2)
  z <- x+y
  xbar[i] <- mean(x)
  ybar[i] <- mean(y)
  res <-optim(c(xbar[i],ybar[i]),ll,control =list(fnscale=-1))
  est1[i] <- res$par[1]
  est2[i] <- res$par[2]
}

boxplot(est1,xbar,xlab="method",main=expression(hat(lambda)))
axis(1,1:2,1:2)
abline(h=8,lty=2,lwd=2,col="red2")
boxplot(est2,ybar,xlab="method",main=expression(hat(mu)))
axis(1,1:2,1:2)
abline(h=2,lty=2,lwd=2,col="red2")

関連エントリ:abrahamcow.hatenablog.com