共役分布 2 カテゴリカル分布の場合

カテゴリカル分布　多項分布　の世界

　カテゴリカルは世界をK個に分ける。K個のサブ世界の比率は　λ1　λ2　λ3　...　λｋ

で、λ1+λ2+λ3+　...　+λｋ＝1　であり　それぞれは　0または正　かつ　1または1以下である。1000個のデータを（λ1, λ2, λ3, λ4)＝(0.1, 0.2, 0.4, 0.3) の場合

のもとで　rmultinom()　を用いて抽出できる。

---

set.seed(126)
n<-1000
lambda<-c(0.1,0.2,0.4,0.3)
k<-length(lambda)
catk<-rmultinom(n,1 ,prob=lambda)
row.names(catk)<-as.character(1:k)
barplot(apply(catk,1,sum)/n,col=3,main="Categolical distribution n=1000")

---

このときデータが（λ1, λ2, λ3, λ4)＝(0.1, 0.2, 0.4, 0.3)の世界から生み出されたであろう確率を計算できる。ただし定数項

gamma(n)/(gamma(n1)…gamma(nk))*Π λ k^nK

を無視している。

---

PrCat<-function(x,lambda){
  k<-nrow(x)
  n<-ncol(x)
  p<-1
    for(ki in 1:k){ p<- p*(lambda[ki])^(sum(x[ki,])/n)}
  p}

PrCat(x=catk,lambda=lambda)

---

[1] 0.2594558

---

ここでlambdaをいろいろ変えて　PrCatを計算すれば　データXにふさわしいλがみつかるはずである。

lambdaを過不足なくそろえるのが難しい。

---

seek4lambda<-function(bigin=c(0,0,0,0),end=c(1,1,1,1),step=c(0.1,0.1,0.1,0.1)){
  la<-NULL
  for(l4 in seq(bigin[4],end[4], step[4])){
    for(l3 in seq(bigin[3] ,min(1-l4,end[3]),step[3])){
      for(l2 in seq(bigin[2], min(1-l4-l3,end[2]),step[2])){
        l1=1-l4-l3-l2
        if (l1>=bigin[1] && l1<=end[1]) la<-rbind(la,c(l1,l2,l3,l4))
  }}}
  la}
#####
la<-seek4lambda(step=c(0.025,0.025,0.025,0.025))
head(la)

---

##       [,1]  　[,2] 　[,3] [,4]
## [1,] 1.000 0.000    0    0
## [2,] 0.975 0.025    0    0
## [3,] 0.950 0.050    0    0
## [4,] 0.925 0.075    0    0
## [5,] 0.900 0.100    0    0
## [6,] 0.875 0.125    0    0

---

このlambdaを順にPrCat(x, lambda)で計算できる。

---

Likelifood.Cat<-function(x,lam){
  Lk<-NULL
  for (i in 1:nrow(lam)){
    p<-PrCat(x=x,lambda=lam[i,])
    if( p==Inf) p<-0  #
    Lk<-rbind(Lk,c(lam[i,],p))
  }
  Lk<-data.frame(Lk)
  sp<-sum(Lk[,5]) # 1/sp = CONSTANT
  Lk[,5]<-Lk[,5]/sp
  names(Lk)<-c("L1","L2","L3","L4","Likelihood")
  Lk
}
#####
library(scatterplot3d)
col.n<-100
ch.col<-topo.colors(col.n+1,alpha=0.2)
#####
#n=1000
LkCat<-Likelifood.Cat(x=catk,lam=la)
ch.max<-max(LkCat$Likelihood)
col<-ch.col[(1-LkCat$Likelihood/ch.max)*col.n+1]
scatterplot3d(LkCat[,1:3],color=col,pch=20,angle=140,main="Likelihood distribution n=1000")

---

ここでLikelihoodが最大になる点を求めることにより、確かに世界に隠されたlamdaをデータXから求めることができる。ML（マキシマムライクリフッド）である。

---

LkCat[LkCat$Likelihood==max(LkCat$Likelihood),]

##       L1  L2  L3  L4   Likelihood

## 8199 0.1 0.2 0.4 0.3 0.0001633475

さてカテゴリカル分布の共役分布はDirichlet 分布である。

Drch<-function(a,lambda){

  PrD<-NULL  for(n in 1:nrow(lambda)){

    lam<-lambda[n,]

    p<-1

    for(ki in 1:length(a)){

        p<-p*lam[ki]^(a[ki])}

    PrD<-rbind(PrD,c(lam,p))}

  sp<-sum(PrD[,5])

  PrD[,5]<-PrD[,5]/sp

  PrD<-data.frame(PrD)

  names(PrD)<-c("L1","L2","L3","L4","prior")

  PrD}   

#

a<-c(0.1,0.2,0.4,0.3)*5

prior<-Drch(a=a,lambda=la)

ch.max<-max(prior$prior)

col<-ch.col[(1-prior$prior/ch.max)*col.n+1]

scatterplot3d(prior[,1:3],color=col,angle=140,pch=20,main="Dirichlet distribution")

Bayes’ ruleによれば

　　Posterior = Liklihood * Prior / Evidence

ここで

　　Liklihood = PrCat(x,λ); Pr(x|λ)
 　  Prior = Drch(a,λ)       ; Pr(λ)
      Evidence = ∫ PrCat(x,λ)dλ

---

Posterior.Cat<-function(x,a,lam){
  k<-length(a)
  n<-ncol(x)
  p<-Likelifood.Cat(x,lam)$Likelihood*Drch(a=a,lambda=lam)$prior
  p<-p/sum(p)
  Post<-data.frame(lam,p)
  names(Post)<-c("L1","L2","L3","L4","posterior")
  Post}
#
Post<-Posterior.Cat(x=catk,a=a,lam=la)
ch.max<-max(Post$posterior)
col<-ch.col[(1-Post$posterior/ch.max)*col.n+1]
scatterplot3d(Post[,1:3],color=col,pch=20,angle=140,main="Posterior distribution n=1000")

---

ここでposteriorが最大になる点を求めれば、これがMAPと言われているものである。

Post[Post$posterior==max(Post$posterior),]

##       L1  L2  L3  L4    posterior

## 8199 0.1 0.2 0.4 0.3 0.0005428246

参照　http://rpubs.com/HN-JP/194782