廿TT

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

ポアソン混合モデルの変分ベイズによる推定をRで(ベイズ推論による機械学習入門)

今日の川柳

モデル

パラメータ \lambdaポアソン分布の確率関数を {\rm Poi}(x|\lambda) と書くことにする.

x_n の確率関数を,
\displaystyle \prod_{k=1}^{K} {\rm Poi}(x_n|\lambda_k)^{s_{n,k}}
とする.

ここで  \boldsymbol{s}_n = \{s_{n,1},\ldots,s_{n,k}\} はカテゴリカル分布にしたがう変数とする.  \boldsymbol{s}_n は観測されない潜在変数である.

\lambda_k の事前分布にパラメータ a, b のガンマ分布を仮定する.

カテゴリカル分布のパラメータ \boldsymbol{\pi} の事前分布には, パラメータ \boldsymbol{\alpha} のディリクレ分布を仮定する.

R による計算例

 \boldsymbol{s}_n, \lambda_k, \boldsymbol{\pi} の期待値を推定するために変分ベイズを使う.

以下のようなコードを書いた.

#x: 観測されたデータ
#a,b: ガンマ分布のパラメータ
#alpha: ディリクレ分布のパラメータ
#lambda_ini: ポアソン分布のパラメータの初期値
#pi_ini: カテゴリカル分布のパラメータの初期値
VBPoissonMix <- function(x,a,b,alpha,K,lambda_ini,pi_ini,tol=1e-6,maxit=1000){
  N<-length(x)
  s <- tmp <- matrix(NA,N,K)
  lambda_old <- lambda_ini
  pi_old <- pi_ini
  loglambda <- log(lambda_ini)
  logpi_old <- log(pi_ini)
  for(k in 1:K){
    tmp[,k] <- exp(x*loglambda[k]-lambda_old[k]+logpi_old[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,2,sum) + b
    alphahat <- apply(s,2,sum) + alpha
    lambda_new <- ahat/bhat
    loglambda <- digamma(ahat)-log(bhat)
    logpi_new <- digamma(alphahat)-digamma(sum(alphahat))
    for(k in 1:K){
      tmp[,k] <- exp(x*loglambda[k]-lambda_new[k]+logpi_new[k])
    }
    den <-apply(tmp,1,sum)
    s<-sweep(tmp,1,den,"/")
    if(all(abs(lambda_old-lambda_new)<tol) & all(abs(logpi_old-logpi_new)<tol)){
      break
    }
    lambda_old <- lambda_new
    logpi_old <- logpi_new
  }
  list(lambda=lambda_new,pi=alphahat/sum(alphahat),s=s,iter=i)
}

なぜこういうアルゴリズムになるのかは本を見てください. 式の展開とかがめっちゃていねいです。

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

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

適当に以下のような乱数を生成してそこからパラメータを推定してみる.

f:id:abrahamcow:20171224231040p:plain

set.seed(1)
N <- 1000
true_lambda <- c(5,15)
true_s <- rbinom(N,1,0.5)
x <-rpois(N,true_lambda[true_s+1])

hist(x,breaks = "scott",main = "",xlab = "")
est1 <-VBPoissonMix(x = x,a = 0.001,b = 0.001,alpha = c(1,1),K = 2,
                   lambda_ini = c(3,14),pi_ini = c(0.5,0.5))

だいたい元のパラメータに近い値が求まる.

> est1$iter
[1] 27
> est1$lambda
[1]  5.016922 14.953575
> est1$pi
[1] 0.5331708 0.4668292

 \boldsymbol{s}_n の期待値の分布を見てみる. シミュレーションで真の  \boldsymbol{s}_n がわかっているので, それで層別してヒストグラムを書いた.

library(tidyverse)
library(cowplot)
checkdf <-data.frame(true_s,probability=est1$s[,1])

p1 <-ggplot(checkdf,aes(x=probability))+
  geom_histogram(colour="black", fill="white", bins=10)+
  facet_grid(true_s~.)

print(p1)

f:id:abrahamcow:20180318123928p:plain

けっこうはっきり傾向が分かれている.

 \boldsymbol{s}_n の期待値の平均をヒストグラムのビンごとに平均したもので色分けして, x_nヒストグラムを書いてみる.

h1<-hist(x,plot = FALSE,breaks = "scott")

df1 <-data.frame(x=x,s=est1$s[,1]) %>% 
  mutate(cut=cut(x,h1$breaks,include.lowest=TRUE,right=TRUE)) %>% 
  group_by(cut) %>% 
  summarise(count=n(),s=mean(s)) %>% 
  mutate(mids=h1$mids)

p2 <-ggplot(df1,aes(x=mids,y=count,fill=s))+
  geom_col(width = h1$breaks[2]-h1$breaks[1],colour="black")+
  scale_fill_gradient2(low = "cornflowerblue",mid = "white",high = "orange",midpoint = 0.5)+
  xlab("")

print(p2)

f:id:abrahamcow:20180318123515p:plain

こんな感じだ.