理学療法士がまとめたノート

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

ポアソン分布の最尤推定

忘れないうちに書いときます
ポアソン分布とは…
二項分布において、n(→∞)が大きく、p(→0)が小さい場合の確率分布.

f:id:yoshida931:20170524162954p:plain



平均1,3,5,7のポアソン分布
y<-0:15
prob<-dpois(y,lambda = 1)  #ラムダ=1
plot(y,prob,type = "b",lty=2,xlim=c(0,15),ylim=c(0,0.4),xlab = "",ylab = "")
prob<-dpois(y,lambda = 3)
lines(y,prob,type = "b",lty=2,pch=2)
prob<-dpois(y,lambda = 5)
lines(y,prob,type = "b",lty=2,pch=4)
prob<-dpois(y,lambda = 7)
lines(y,prob,type = "b",lty=2,pch=5)
legend(12,0.35,legend=c("λ=1","λ=3","λ=5","λ=7"),pch = c(1,2,4,5))

f:id:yoshida931:20170524153211p:plain

 

ポアソン分布のパラメータλの最尤推定の可視化
準備として平均2.8になるポアソン分布の乱数を50個用意します
x<-rpois(50,2.8)

尤度対数
logL<-function(m)sum(dpois(x,m,log=TRUE))  
# 対数尤度関数 sum(dpois(x,m,log=TRUE))
# ポアソン分布の密度関数(dpois)を対数変換(log=TRUE)したもの

lambda<-seq(1,5,0.1)
#λ(平均) 1~5、0.1間隔

plot(lambda,sapply(lambda,logL),type="l",xlab = "λ",ylab = "log likelihood") 
#sapply() logLにlambda<-seq(2,5,0.1) を挿入した値を返す

f:id:yoshida931:20170524153325p:plain

平均2.8になるポアソン分布の乱数なので、もちろん最終推定値は2.8付近になります.
計算では対数尤度関数をλで偏微分して、=0となる値を求めます.

 

参考資料
久保拓哉:確率分布と一般化線形モデル(GLM)
http://hosho.ees.hokudai.ac.jp/~kubo/stat/2012/kobe/k2/kubostat2012k2.pdf