廿TT

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

周辺尤度を重点サンプリングで計算するときの提案分布に変分事後分布を使う

今日の川柳

といいんじゃないかなーと思った。

f:id:abrahamcow:20190606210210j:plain
(計算統計学の方法 p.199)

ここでは分布 g のことを「提案分布」と呼ぶことにしました。

以下、R のコードです。

混合二項分布のパラメータを平均場近似による変分ベイズ法で推定してます。

事前分布はベータ分布とディリクレ分布です。

 p_l \sim \mathrm{beta}(a,b)
 \pi_l \sim \mathrm{Dirichlet}(\alpha)
 z \sim \mathrm{categorical}(\pi)
 x \sim \mathrm{binomial}(10,p_z)

混合数は3が正解です。

library(gtools)
VBmixbinom <- function(x,m,K,a=1,b=1,alpha=rep(1,K),maxit=1000){
  p_ini <- runif(K)
  pi_ini <- rep(1/K,K)
  N<-length(x)
  m_x <- m-x
  s <- tmp <- matrix(NA,N,K)
  p_old <- p_ini
  pi_old <- pi_ini
  logp <- log(p_ini)
  log1_p <- log(1-p_ini)
  logpi <- log(pi_ini)
  for(k in 1:K){
    tmp[,k] <- exp(logpi[k]+x*logp[k]+m_x*log1_p[k])
  }
  den <-apply(tmp,1,sum)
  s<-sweep(tmp,1,den,"/")
  for(i in 1:maxit){
    ahat <- apply(s*x,2,sum) + a
    bhat <- apply(s*m_x,2,sum) + b
    alphahat <- apply(s,2,sum) + alpha
    logp <- digamma(ahat)-digamma(ahat+bhat)
    log1_p <- digamma(bhat)-digamma(ahat+bhat)
    logpi <- digamma(alphahat)-digamma(sum(alphahat))
    for(k in 1:K){
      tmp[,k] <- exp(logpi[k]+x*logp[k]+m_x*log1_p[k])
    }
    den <-apply(tmp,1,sum)
    s<-sweep(tmp,1,den,"/")
  }
  list(ahat=ahat,bhat=bhat,alphahat=alphahat,s=s,iter=i)
}
ISml <- function(out,np=10000){
  K <- length(out$ahat)
  pi_s <- t(rdirichlet(np,out$alphahat))
  p_s <- matrix(rbeta(np*K,out$ahat,out$bhat),K,np)
  log(mean(exp(log(colSums(pi_s*dbinom(x,10,p_s)))-
                 log(ddirichlet(t(pi_s),out$alphahat))+colSums(dbeta(p_s,out$ahat,out$bhat,log = TRUE)))
           ))
}
set.seed(1111)
x <-c(rbinom(100,10,0.9),rbinom(50,10,0.1),rbinom(100,10,0.5))
hist(x)
out2 <- VBmixbinom(x,10,2)
ISml(out2)
out3 <- VBmixbinom(x,10,3)
ISml(out3)
out4 <- VBmixbinom(x,10,4)
ISml(out4)

すべてのパラメータが1のディリクレ分布やベータ分布は、一様分布と一緒なので事前分布の密度は省略しました。

3のとき一番大きい値がでてる。

> ISml(out2)
[1] 2.433358
> out3 <- VBmixbinom(x,10,3)
> ISml(out3)
[1] 3.800278
> out4 <- VBmixbinom(x,10,4)
> ISml(out4)
[1] 0.6527613

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)