기사의 목적

목차

0. 라이브러리

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

1. 모델 설명

2. 추정할 분포

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. 확인

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)