Rで行うEMアルゴリズム~ゼロ過剰ポワソン分布を例にして~

今回は、EMアルゴリズムを用いて、ゼロ過剰ポワソン分布のパラメータの推定を行ってみようと思います。

過去記事
最尤推定
EMアルゴリズム
- アルゴリズム
- 実例

過去記事

ゼロ過剰ポワソン分布に関する詳しい説明はこちらをご覧ください。
ゼロ過剰ポアソン分布とは何か - バナナでもわかる話

前にゼロ過剰ポワソン分布に関する尤度の説明は記事にしましたが、今回は潜在変数を利用するので、また新しく説明しようと思います！
ゼロ過剰ポアソン分布の尤度関数 - バナナでもわかる話

最尤推定

ゼロ過剰ポワソン分布の尤度関数

ゼロ過剰ポワソン分布から得られたデータがn個あるとします。
ここで、普通のポワソン分布から出ている(観測されない未知の)サンプルサイズをMと置くことにします。
つまり、ゼロのみからなる標本空間から得られたサンプルがn-M個あると考えても問題ありません。

更に、今回は確率 $\omega$ でゼロのみからなる標本空間からサンプリングされると仮定します。
そうするとこの時の尤度関数 $l(\lambda,\omega)$ は次のようになりますね。ただ、Mは一旦固定しておきます。

$l(\lambda,\omega)= \omega^{n-M} \prod_{i=1}^M \{\frac{(1-\omega) \lambda^{y_i}exp(-\lambda)}{y_i !} \}$

最尤推定量

対数尤度 $ln(\lambda,\omega)$ を求めます。
$ln(\lambda,\omega)=(n-M)log(\omega) +M log(1-\omega)+\sum_{i=1}^M \{y_i log(\lambda)-\lambda -log(y_i !)\}$

Mはとりあえず放っておいて、最尤推定量を求めてみます。
$\frac{\partial ln(\lambda,\omega)}{\partial \omega} =\frac{n-M}{\omega}-\frac{M}{1-\omega}=0$
$\hat{\omega} = 1-\frac{M}{n}$

$\frac{\partial ln(\lambda,\omega)}{\partial \lambda} =-M+\frac{\sum_{i=1}^M y_i}{\lambda}=0$
$\hat{\lambda} = \frac{\sum_{i=1}^M y_i}{M}$

Mの理論値

とりあえず、1以上の値は全て普通のポワソン分布から出ていることがわかっています。そこで1以上の個数 $m=\#(Y_i≧1)$ と置くことにします。
この時点でM≧mであることが分かっています。
さて、何個の0が普通のポワソン分布に属しているかという話ですが、理論的にはポワソン分布から0が出る確率は $exp(-\lambda)$ なので、
サンプルサイズMのもとでの、0の個数に関する期待値は $M exp(-\lambda)$ ですね。そこで、Mの理論値は次のように書けます。

$M = m +M exp(-\lambda)$

Mの理論値にMを使ってることに少し違和感がありますが、これをMについて解いてやります。
$\hat{M}=\frac{m}{1-exp(-\lambda)}$

EMアルゴリズム

$M$ の推定には $\lambda$ を使うし、 $\lambda$ や $\omega$ の推定には $M$ を使います。よって、どちらかが分からない限り推定値を得ることができません。
そこで、EMアルゴリズムという手法を利用します。

EMアルゴリズムでは、 $M$ 又は $\lambda$ に初期値を与えてもう片方を推定し、その推定結果を利用して、またもう片方を推定する...という手順を繰り返し、逐次的に解を得ます。

アルゴリズム

今回は、まず初期値として $\lambda_{init} =\frac{\sum_{i=1}^n y_i}{n}$ を与えることにします。
そして、この $\lambda_{init}$ を利用して $M$ の初期値 $M_1$ を次のように得ます。
$M_1 = \frac{m}{1-exp(-\lambda_{init})}$
更に $M_1$ を利用して $\lambda_1 =\frac{\sum_{i=1}^M y_i}{M_1}$ を得られます。
※どうせポワソン分布以外から出てくるのは0なので、総和に変化はない。

また、 $\lambda_1$ を利用して $M_2$ を....といったようにそれぞれの数列
$M_1,M_2,\cdots$
$\lambda_1,\lambda_2,\cdots$
を得ます。

そして、 $\lambda$ の列が収束した段階で探索を止めます。

実例

シミュレーションデータ

prob=0.1
data=rpois(100000,lambda=1)*rbinom(100000,1,1-prob)

$\omega = 0.1,\lambda =1 , n=100000$

EMアルゴリズム

#lambda
init_lam=mean(data)
est_lam=init
alt_lam=100
#Mはポワソンから出ているサンプルの個数(潜在変数)
init_M=0
est_M=init_M
m=sum(data>0)
n=length(data)
S=sum(data)
while(abs(est_lam-alt_lam)>=0.001){
	alt_lam=est_lam
	est_M=m/(1-exp(-est_lam))
	est_lam=S/est_M
	if(abs(est_lam-alt_lam)<0.001){
		omega=1-est_M/n
	}
}

推定結果

> omega
[1] 0.09526407
> est_lam
[1] 0.9954507

かなり良い値が出てきましたね。