Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

Rでポアソン混合モデル実装 (ギブスサンプリング)

記事の目的

ポアソン混合モデルのギブスサンプリングをRで実装します。
参考:ベイズ推論による機械学習入門

目次

0. ライブラリ
1. モデルの説明
2. 推定する分布
3. パラメータ初期値設定
4. 事後分布
5. 確認

0. ライブラリ

library(dplyr)
library(mvtnorm)
library(MCMCpack)
library(gganimate)
set.seed(100)

1. モデルの説明

IMG_0178.jpeg

2. 推定する分布

50が平均のポアソン分布と、80が平均のポアソン分布が1:3の割合で混合しているモデルです。

N <- 1000
s.true  <- rmultinom(N, 1, c(1/4, 3/4)) %>% apply(2, which.max)
X <- append(rpois(N, 50)[s.true==1], rpois(N, 80)[s.true==2])
NULL %>% ggplot(aes(x=X)) + geom_histogram(position = "identity", alpha=0.7)
N <- length(X)

image.png

3. パラメータの初期値設定

#ハイパーパラメータ設定
K <- 2
alpha0 <- rep(1/K, K)
a0 <- rep(1, K)
b0 <- rep(1, K)
#パラメータ初期値設定
pi <- rep(1/K, K)
s <- rep(1:K, N)[1:N]
lambda <- rep(100, K)
alpha <- {}
a <- {}
b <- {}
s.Data1 <- data.frame(s, iter=rep(0, N))

4. 事後分布

max.iter <- 30
eata <- matrix(NA, nrow = N, ncol=K)
for(t in 1:max.iter){
  #sのサンプリング
  for(i in 1:N){
    eata[i, ] <- (X[i]*log(lambda) - lambda + log(pi)) %>% exp()
    s[i] <- rmultinom(1, 3, eata[i, ]) %>% apply(2, which.max)
  }
  s.Data2 <- data.frame(s, iter=rep(t, N))
  s.Data1 <- rbind(s.Data1, s.Data2)
  #lambdaのサンプリング
  for(k in 1:K){
    a[k] <- sum(X[s==k]) + a0[k]
    b[k] <- (X[s==k] %>% length()) + b0[k]
    lambda[k] <- rgamma(1, a[k], b[k])
  }
  #piのサンプリング
  for(k in 1:k){
    alpha[k] <- (X[s==k] %>% length()) + alpha0[k]
  }
  pi <- rdirichlet(1, alpha)
}

5. 確認

iter15でだいたい分類できていることが分かります。

s.Data <- s.Data1 %>% filter(iter==0|iter==5|iter==15|iter==20|iter==25|iter==30)
s.Data <- data.frame(X=rep(X, 6), s.Data)
s.Data$s <- s.Data$s %>% as.factor()
s.Data %>% ggplot(aes(x=X, fill=s)) + geom_histogram(position = "identity", alpha=0.7) +
  transition_states(iter, transition_length = 2, state_length = 1)

a5.gif

Tatsuki-Oike
I am a third year university student in Japan. I started qiita as a memorandum.Nice to meet you.
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away