혼합 모델(Mixture model)

다음의 데이터의 히스토그램을 확인하고, 이를 적합할 분포에 대하여 생각해 보자.

 

In:

data = read.csv('../input/mixture.csv', header = F)

y = data$V1
n = length(y)

hist(y, breaks = 20)

 

Out:

 

▷ 일반적인 분포와는 다르게 -2와 1 근처에 두 개의 봉우리를 형성하고 있는 것을 확인할 수 있다. 우리가 알고 있는 정규분포, 지수분포, 감마분포는 단봉형태를 띄고 있다. 따라서 이를 위의 데이터에 적합할 경우, 두 개의 봉우리에 대해 적합한 분포를 얻을 수 없다.

 

▶ 혼합 모델(Mixture model)은 두 개이상의 분포를 합쳐서 만든 모델로 기존 분포에 비하여 자유도 높은 적합이 가능하다는 장점을 가지고 있다. 또한 복잡한 분포의 적합뿐만 아니라, 군집화에 사용되기도 한다.

 

위의 데이터에 대하여 혼합 모델을 이용하여 분포를 적합시켜보자.

 

In:

library('rjags')

mod_string = " model {
  for (i in 1:length(y)) {
    y[i] ~ dnorm(mu[z[i]], prec[z[i]])
    z[i] ~ dcat(omega)
  }
  
  mu[1] ~ dnorm(-1.0, 1.0/100.0)
  mu[2] ~ dnorm(1.0, 1.0/100.0) T(mu[1],)
  
  prec[1] ~ dgamma(1.0/2.0, 1.0*1.0/2.0)
  prec[2] ~ dgamma(1.0/2.0, 1.0*1.0/2.0)
  
  sig[1] = sqrt(1.0/prec[1])
  sig[2] = sqrt(1.0/prec[2])
  
  omega ~ ddirich(c(1.0, 1.0))
} "

data_jags = list(y = y)

params = c('mu', 'sig', 'omega')

mod = jags.model(textConnection(mod_string), 
                 data = data_jags, 
                 n.chains = 3)

update(mod, 1e3)

mod_sim = coda.samples(model = mod,
                       variable.names = params,
                       n.iter = 5e3)

mod_comb_sim = as.mcmc(do.call(rbind, mod_sim))

colMeans(mod_comb_sim)

 

Out:

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 200
   Unobserved stochastic nodes: 205
   Total graph size: 817

Initializing model

  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
  |**************************************************| 100%
  |**************************************************| 100%
     mu[1]      mu[2]   omega[1]   omega[2]     sig[1]     sig[2] 
-2.0624518  1.5038066  0.3987204  0.6012796  1.1804035  1.1344625 

 

▷ 혼합 모델의 구성요소로써 두 개의 정규분포를 이용하였다. 위의 결과로부터 얻은 혼합 모델의 식은 다음과 같다.

 

위의 식을 이용하여 혼합 모델의 결과를 확인하여 보자.

 

In:

dmixnorm = function(x) {
  return(0.40 * dnorm(x, -2.06, 1.18) + 0.60 * dnorm(x, 1.50, 1.13))
}

x = seq(-6, 6, by = 0.01)

plot(x, dmixnorm(x), type = 'l', col = 'blue')

hist(y, breaks = 20, freq = F, add = TRUE)

 

Out:

 

▷ 혼합 모델이 두 개의 봉우리를 가진 데이터의 분포에 잘 적합한 것을 확인할 수 있다.

 

---

Reference:

"Bayesian Statistics: From Concept to Data AnalysisTechniques and Models," Coursera, https://www.coursera.org/learn/bayesian-statistics/.