기사의 목적

포아송 혼합 모델의 깁스 샘플링을 R로 구현합니다.
참고 : 베이즈 추론에 의한 기계 학습 입문

목차

0. 라이브러리
1. 모델 설명
2. 추정할 분포
3. 파라미터 초기값 설정
4. 사후 분포
5. 확인

0. 라이브러리

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

1. 모델 설명

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)

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)