廿TT

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

EM アルゴリズムによるゼロ切断ポアソン分布のパラメータ推定

最尤推定と条件付き最尤推定

ゼロ切断ポアソン分布のパラメータの最尤推定 - 廿TT の続きです。

上記のエントリでは、0 でない観測の数 n を所与として、n 個の観測が得られるために必要な試行回数 N は観測されないとしていました。

しかし現実的には n が固定されておらず N が固定されている例のほうが多いかもしれません。

この場合、ゼロ切断ポアソン分布のパラメータの最尤推定 - 廿TT で与えた尤度関数、

\displaystyle \log L_1 (\lambda|n) = \log(\lambda)\sum_{i=1}^{n} k_i - n\lambda - n\log(1-e^{-\lambda})+\mbox{定数}

は、 n が得られたという条件の下での条件付き対数尤度関数です。

完全な対数尤度関数は、n が得られる確率も考慮して、

\displaystyle\log L_2(\lambda)=\log(\lambda) \sum_{i=1}^{n}k_i -n \lambda -n \log(1-e^{-\lambda})\\
\displaystyle~~~+n \log(1-e^\lambda)+(N-n)\log e^{-\lambda}+\mbox{定数}\\
\displaystyle = \log(\lambda) \sum_{i=1}^{n} k_i -N \lambda+\mbox{定数}

となります。(上式の後半部は二項分布の確率関数です。)

一回微分して λ について解くと、最尤推定量は、

\displaystyle \hat \lambda =\frac{\sum_{i=1}^{n}}{N}

です。

最尤推定と条件付き最尤推定をシミュレーションで比較してみます。

library(tidyr)
library(cowplot)
N=100
lambda=1
iter=10000
#条件付き尤度関数
logL0 <- function(ysum,n){
  function(lambda) {log(lambda)*ysum - n*lambda - n*log(1-exp(-lambda))}
}
lambdahat1 <- lambdahat2<- numeric(iter)
for(i in 1:iter){
  y <- rpois(N,lambda)
  y2 <- y[y != 0]
  ysum <- sum(y2)
  n=length(y2)
  logL1 <- logL0(ysum,n)
  opt1 <- optimise(logL1,c(0,100),maximum = TRUE)
  lambdahat1[i] <-opt1$maximum
  lambdahat2[i] <-ysum/N
}
est_df <- data.frame(lambdahat1,lambdahat2)

ggplot(gather(est_df),aes(x=value))+
  geom_histogram(aes(fill=key,colour=key),alpha=0.3,position="identity")

f:id:abrahamcow:20160806010425p:plain

赤いヒストグラムが条件付き最尤推定量、青いヒストグラム最尤推定量の分布です。

λ=1 と設定してシミュレートしました。

最尤推定量のほうが分散が小さく抑えられています。

EM アルゴリズム

N がわかっていれば最尤推定量を使えばいいのですが、0 が観測されない分布で N がわかっている状況はまずありません。

でもなんとか N の情報も入れて推定してやりたい。どうしたらいいか、と考えて、EM アルゴリズムを使ったらどうだろうと思いつきました。

0 でない値が n 回出るまでの試行回数の分布は負の二項分布になるので、N の期待値は

 \displaystyle \frac{n}{1-\exp(-\lambda)}

です(E ステップ)。

N も含めた完全な観測が得られたという条件下での尤度関数と、その最大化は上述の通りです(M ステップ)。

EM アルゴリズムによる推定量と、条件付き最尤推定量をシミュレーションで比較します。

myEM <- function(lambda0,ysum,n){
  for(i in 1:1000){
    EN <- n/(1-exp(-lambda0)) #E step
    lambda2 <-ysum/EN #M step
    if(abs(lambda0-lambda2)<1e-6){break}
    lambda0 <- lambda2
  }
  return(c(par=lambda2,iter=i))
}

lambdahat3 <- lambdahat4 <- numeric(iter)
for(i in 1:iter){
  y <- rpois(N,lambda)
  y2 <- y[y != 0]
  ysum <- sum(y2)
  n=length(y2)
  logL1 <- logL0(ysum,n)
  opt1 <- optimise(logL1,c(0,100),maximum = TRUE)
  opt2 <- myEM(1,ysum,n)
  lambdahat3[i] <-opt1$maximum
  lambdahat4[i] <-opt2[1]
}
est_df <- data.frame(lambdahat3,lambdahat4)

ggplot(est_df,aes(x=lambdahat3,y=lambdahat4))+
  geom_point(alpha=0.3)

結果、両者はかなりぴったり一致しました。意味ねー。

f:id:abrahamcow:20160806011621p:plain

ゼロ切断ポアソン分布によるセッション数のモデル化とデータ拡大による潜在利用者数の推定 - 廿TT に続く。