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

廿TT

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

ゼロ切断ポアソン分布のパラメータの最尤推定

R 確率分布

確率質量関数

fポアソン分布の確率質量関数とすると、ゼロ切断ポアソン分布(zero-truncated Poisson)の確率質量関数は、

\displaystyle g(k;\lambda )=P(X=k\mid X>0) \\ 
\displaystyle~~~={\frac {f(k;\lambda )}{1-f(0;\lambda )}} \\
\displaystyle~~~={\frac {\lambda ^{k}e^{-\lambda }}{k!\left(1-e^{-\lambda }\right)}}\\
\displaystyle~~~={\frac {\lambda ^{k}}{(e^{\lambda }-1)k!}}

期待値は、

\displaystyle E[X]=\frac {\lambda }{1-e^{-\lambda }}=\frac {\lambda e^{\lambda }}{e^{\lambda }-1}

Zero-truncated Poisson distribution - Wikipedia, the free encyclopedia

最尤推定

サンプルサイズを n とすると対数尤度関数は、

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

これを最大化するパラメータ λ は解析的に求まりそうで求まらない。

対数尤度関数を一回微分すると、

\displaystyle \frac{\partial}{\partial \lambda} \log L (\lambda) = \frac{1}{\lambda}\sum_{i=1}^{n} k_i - n - n\frac{e^{-\lambda}}{1-e^{-\lambda}}

二回微分すると、

\displaystyle \frac{\partial^2}{\partial \lambda^2} \log L (\lambda) =  - \frac{1}{\lambda^2}\sum_{i=1}^{n} k_i  + n\frac{e^{-\lambda}}{1-e^{-\lambda}}

期待値を取ると、フィッシャー情報量が求まり、

\displaystyle I= E \left( \frac{\partial^2}{\partial \lambda^2} \log L (\lambda)\right) =  -\frac{n}{1-e^{-\lambda}}\left( \frac{1}{\lambda} - \frac{e^{-\lambda}}{1-e^{-\lambda}}\right)

漸近的にはこれの負の逆数が最尤推定量の分散になる。

シミュレーション

R です。

パラメータ λ=1 のゼロ切断ポアソン乱数を 100 個生成して、最尤推定を行うことを 10000 回繰り返す。

推定量は平均 1、分散 -1/I正規分布に従っていることがわかる。

library(tidyr)
library(cowplot)
lambda <- 1   # Poisson mean
n <- 100      # Sample size
# Define log of the likelihood function
logL0 <- function(ysum,n){
  function(lambda) {log(lambda)*ysum - n*lambda - n*log(1-exp(-lambda))}
}
lambdahat <- N <- numeric(10000)
for(i in 1:10000){
  y <- rpois(n*30,lambda) # 多めに出しておく
  N[i] <-which.max(cumsum(y!=0)==100)
  y2 <- y[y != 0][1:n]
  ysum <- sum(y2)
  logL <- logL0(ysum,n)
  opt1 <- optimise(logL,c(0,100),maximum = TRUE)
  lambdahat[i] <-opt1$maximum
}
hist(lambdahat,freq = FALSE)
NFI <- (n/(1-exp(-lambda)))*(1/lambda - exp(-lambda)/(1-exp(-lambda)))
curve(dnorm(x,lambda,sqrt(1/NFI)),add = TRUE)

ヒストグラム正規分布の密度関数を重ねてプロットした。

f:id:abrahamcow:20160805025758p:plain

この事実から正規近似による信頼区間を求めることができる。

ZTPCI <- function(lambda,n,level=0.95){
  NFI <- (n/(1-exp(-lambda)))*(1/lambda - exp(-lambda)/(1-exp(-lambda)))
  alpha = 1-level
  lower=qnorm(alpha/2,lambda,sqrt(1/NFI))
  upper=qnorm(1-alpha/2,lambda,sqrt(1/NFI))
  cbind(lower,upper)
}
CIs <- ZTPCI(lambdahat,n)
CP <-mean(CIs[,"lower"]<= lambda & lambda<=CIs[,"upper"])
cat("coverage probability is",round(CP,2))

coverage probability(信頼区間が真の値を含む割合)は約 0.95 で名目上の信頼水準とほぼ一致することがわかった。

隠れた試行回数の推定

λ の推定値より、0 でない値が n 個出るまでに要した試行回数が何回であったかを推定することができるだろうと考えた。

ポアソン分布で 0 の出る確率は  \exp(-\lambda)、0 以外の値の出る確率は  1-\exp(-\lambda)、0 でない値が 100 回出るまでの試行回数の分布は負の二項分布になるので、期待値は

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

分散は、

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

である。

h<-hist(N,freq = FALSE)
points(h$breaks,dnbinom(h$breaks-100,100,1-exp(-lambda)),type = "b")

ヒストグラムと負の二項分布の確率質量関数を重ねてプロットした。

f:id:abrahamcow:20160805032243p:plain

期待値の式に λ の推定値をプラグインしてやると隠れた試行回数(N としよう)が求まる。

Nhat <-n/(1-exp(-lambdahat))
Ndf <-data.frame(N,Nhat)
ggplot(Ndf,aes(N,Nhat))+
  geom_point(alpha=0.2)+
  xlim(range(Ndf))+ylim(range(Ndf))+
  geom_abline()

横軸に N、縦軸に N の推定値をとってプロットした。

f:id:abrahamcow:20160805032709p:plain

もっとぴったり一致してくれるかとおもったけどそうでもない。

NN の推定値のヒストグラム

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

f:id:abrahamcow:20160805032753p:plain

分布をみるとピークの位置はだいたい一致してるけど、N の推定値のほうが右すそが重い分布になっていて、バイアスもある。

> (EN <- n/(1-exp(-lambda)))
[1] 158.1977
> mean(Nhat)
[1] 160.1209
> mean(N)
[1] 158.1816