廿TT

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

変分ベイズによる近似事後分布は相当おおざっぱかもしれない

以下のポアソン回帰モデルを考えます。

 y_i \sim \mathrm{Poisson}(\prod_{j=1}^{K}\lambda_j^{x_{ij}})
 \lambda_j \sim \mathrm{Gamma}(a,b)

x_{ij} (i=1,\ldots,n; j= 1,\ldots,K) はダミー変数で 0 か 1 の値を取るとします。

平均場近似による \lambda_j の近似事後分布はガンマ分布です。

 q(\lambda_j)=\mathrm{Gamma}(\hat a_j,\hat b_j)
\hat a_j = \sum_{i=1}^{n} y_i x_{ij} + a
\hat b_j = \sum_{i=1}^{n} x_{ij}\prod_{j' \neq j} \lambda_{j'}^{x_{ij'}}

です。導出の詳細は省略しますが、気が向いたら追記するかもしれません。

ギブスサンプリングのための条件付き分布も同様です。

ギブスサンプリングによる事後分布と平均場近似による近似事後分布を比較します。

R のコードです。

doGibbs <- function(y,X,a=1,b=1,iter=2000){
  K <- ncol(X)
  sxy = t(X)%*%y
  lambda <- matrix(0,iter,K)
  lmd <- rgamma(K,1,1)
  for (i in 1:iter) {
    for(j in 1:K){
      lmd[j] = rgamma(1,sxy[j]+a, X[,j]%*%exp(X[,-j,drop=FALSE]%*%log(lmd[-j]))+b)
    }
    lambda[i,] <- lmd
  }
  return(lambda)  
}
doVB <- function(y,X,a=1,b=1,iter=200){
  K <- ncol(X)
  sxy = t(X)%*%y
  lambda <- rgamma(K,1,1)
  ahat <- rep(a,K)
  bhat <- rep(b,K)
  for (i in 1:iter) {
    for(j in 1:K){
      ahat[j] = sxy[j]+a
      bhat[j] = X[,j]%*%exp(X[,-j,drop=FALSE]%*%log(lambda[-j]))+b
      lambda[j] <- ahat[j]/bhat[j]
    }
  }
  return(list(lambda=lambda,ahat=ahat,bhat=bhat))
}
set.seed(1)
x <- rbinom(100,1,0.5)
X <- model.matrix(~x)
y <- rpois(100,exp(X%*%c(1,1)))
outGibbs <- doGibbs(y,X)
outVB <- doVB(y,X)
burnin <- 1:1000
hist(outGibbs[-burnin,1],freq = FALSE,breaks = "FD",ylim = c(0,4))
curve(dgamma(x,outVB$ahat[1],outVB$bhat[1]),add=TRUE,col="red")

hist(outGibbs[-burnin,2],freq = FALSE,breaks = "FD",ylim = c(0,4))
curve(dgamma(x,outVB$ahat[2],outVB$bhat[2]),add=TRUE,col="red")

f:id:abrahamcow:20190314223334p:plain

f:id:abrahamcow:20190314223348p:plain

かんたんなモデルなのでサンプルサイズが100もあれば両者はだいたい一致するんじゃないかと思っていましたが、ぜんぜん違いました。
やってみてよかったと思っています。

ギブスサンプリングや変分ベイズの更新式の導出方法について知りたい方には『ベイズ推論による機械学習入門』がおすすめです。

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

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