共役分布 3 ポアソン分布の場合

ポアソンな世界の有様

　ポアソンな世界は単位時間または空間で仕切られ、その単位の中で事象は平均λで生起する。ターム方式のストラテジーゲームである。

　k回起こる確率を式であらわせば　p(k) = λ^k * e^ λ  / k!

　Rでは　rpois() で生成（random generation）できる。

---

set.seed(125)
n<-1000
lambda<-5
poik<-rpois(n, lambda)
barplot(table(poik)/n,col=3,main="Poisson distribution n=1000")

poi20<-rpois(20, lambda)
barplot(table(poi20)/10,col=4,main="Categolical distribution n=20")

---

データからもっともらしい世界像を思い描く

　ポアソン世界から得た　1000個のデータ　と　20個のデータ　から　今のところ知られていないλを　推定しなければならないとしよう。

　ポアソン世界からのデータであることだけはわかっているとすれば、そのようなデータが得られるはずのλを推定できることになる。（λ値を変えながら今あるデータが得られる確率を計算できれば、どのようなλ値の時今あるデータを得る可能性が最も高いかわかるであろう。）

---

PrPoi<-function(x,lambda){
  p<-1
  for(ki in 0:max(x)){
    p<-p*exp(-lambda)*lambda^ki/factorial(ki)  }
  p　}
#
Likelifood.Poi<-function(x){
  Lk<-NULL
  for( lam in 0:max(x)){
     p<-PrPoi(x=x,lambda=lam)
     Lk<-rbind(Lk,c(lam,p))  }
  Lk[,2]<-Lk[,2]/sum(Lk[,2])
  Lk<-data.frame(Lk)
  names(Lk)<-c("lambda","Likelihood")
  Lk　}
#####
#n=1000
LkPoi<-Likelifood.Poi(x=poik)
plot(LkPoi,col=3,type="l",main="Likelihood")
#n=20
LkPoi<-Likelifood.Poi(x=poi20)
lines(LkPoi,col=4）

---

共役分布を用いてベイズの企みに参加する　

想像可能な世界と実現されている世界を繋ぐ

ところでガンマ分布は

 dgamma(y, shape, rate) = rate^shape/gamma(shape) *y^(shape-1)*exp(-rate* y)
 dgamma(y, shape, rate＝１)  = y^(shape-1)*exp(-y)/gamma(shape) 

とポアソン分布と形が似ていることから共役分布と呼ばれる。

---

m<-max(c(poi20,poik))
prior<-dgamma(x=0:m,shape=5,rate=1)
plot(x=0:m,y=prior,col=2,type="l", main="Gamma distribution",xlab="lambda")

---

ここでピークは4だが平均は5であることが面白い。

---

#Max
which(prior==max(prior))-1

[1] 4

#Mean
sum(0:m*prior)

[1] 4.926971

---

ここでベイス　の定理を思い起こせば

　Posterior = Liklihood * Prior / Evidence

　　　　　  =PrPoi(x,lambda) * dgamma(y, shape, rate=1)  / Evidence
                 = k(x,shape) * dgamma(y, x+shape, rate=2)   
                    where k(x,shape)＝1/{x!* (shape+1)!} *Evidence

---

Posterior.Poi<-function(x,a){
  Post<-NULL
  for (l in 0:max(x)){
    p<-1
    for(k in 0:max(x)){
      p<-p*dgamma(l,shape=a+k,rate=2)
    }
    Post<-rbind(Post,c(l,p))   }
  Post<-data.frame(Post)
  names(Post)<-c("Lambda","posterior")
  sp<-sum(Post$posterior)
  Post$posterior<-Post$posterior/sp
  Post }
#####
Post<-Posterior.Poi(x=poik,a=5)
plot(Post,col=3,type="l", main="Posterior distribution")
#
Post<-Posterior.Poi(x=poi20,a=5)
lines(Post,col=4)

---

Likelihoodよりそれらしい範囲に収まってはいるが、あくまでも事前分布がもっともらしければということである。