RはExcelからRStudio - 詳細表示 - OneSegの長寿と健康を願うひとりがき

ポアソンな世界の有様

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

　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)