# MCMC of a mixture and the label switching problem

I generated some data according to a mixture of two lognormals: $f(x) = p \cdot \mathcal LN(\mu_1, \sigma) + (1-p) \cdot \mathcal LN (\mu_2, \sigma)$.

given $p$ and $\sigma$, this is my code to find $\mu_1, \mu_2$ by MCMC with the mcmcPack R Package:

the mixture density:

d.lognorm.mix <- function(x,p=0.7,mu1,mu2,sigma=0.2) {
dens <- p*dlnorm(x,meanlog=mu1,sdlog=sigma) + (1-p)*dlnorm(x,meanlog=mu2,sdlog=sigma)
return(log(dens))


}

the loglikeli function, the data ist stored in the dat variable (with $\mu_1=2, \mu_2=8$)

loglikeli <- function(theta) {  #theta = mu1, mu2
mu1 <- theta[1]
mu2 <- theta[2]
return(sum(d.lognorm.mix(dat, p=0.7, mu1=mu1, mu2=mu2, sigma=0.2)))
}


find the two parameters

library(MCMCpack)
post.samp <- MCMCmetrop1R(loglikeli, theta.init=c(1,1),
thin=1, mcmc=20000, burnin=500,
#     tune=c(1.5, 1.5),
verbose=FALSE, logfun=TRUE)


but the result is way off:

summary(post.samp)

Iterations = 501:20500
Thinning interval = 1
Number of chains = 1
Sample size per chain = 20000
# Empirical mean and standard deviation for each variable, plus standard error of the mean:

Mean       SD  Naive SE Time-series SE
[1,] 2.0703 0.062859 4.445e-04      1.146e-02
[2,] 0.6896 0.003494 2.471e-05      7.299e-05

# Quantiles for each variable:

2.5%    25%    50%    75%  97.5%
var1 1.968 2.0745 2.0797 2.0844 2.0934
var2 0.683 0.6876 0.6896 0.6917 0.6959


Is this an example of the label switching problem? Is there an R package for MCMC of mixtures?

I'd also appreciate an explanation of what exactly the label switching problem is.

-

Caveat: I'm no expert in mixture modeling, this is just what I've absorbed by osmosis.

I'd also appreciate an explanation of what exactly the label switching problem is.

Imagine the simplest case of a mixture model in the context of your problem, where you have two components. Component $A$ we know has mean $\mu_1$. Component $B$ we know has mean $\mu_2$. So we can model this as you do: $$p(x)=\pi\times\mathcal{LN}(\mu_1,\sigma)+(1-\pi)\times\mathcal{LN}(\mu_2,\sigma)=\pi A+(1-\pi)B$$

And you want to do inference on $\pi$ (or also other parameters). But, without changing anything about your model, you could also have written $$p(x)=(1-\pi)A+\pi B$$ Without further constraints on $\pi$, these labeling schemes occur with equal posterior probability.

That is, in the two component problem, there are two modes which can occur with equal posterior probability. Part-way through sampling, the sampler may wander to the second mode and swap between exploring the first expression and the second. That is, it swaps $B$ and $A$. This pathology is called "label switching."

One solution is to constrain your parameters such that $\pi<1-\pi$. This can create its own problems, though; if the opposite is the "truth," and your $\pi$ samples will tend to "pile up" on the boundary. What it means is that your prior is in direct conflict with your data, so you should reverse the inequality. One way to do this is to just enforce the prior $\Pr(\pi<0.5)=1$ (or the opposite).

Obviously, as you increase the number of components, the number of modes can proliferate. I'm not sure what one does in that case, since the probability of guessing the correct ordering and correspondence of $\left\{ \pi_1, \pi_2, ..., \pi_i \right\}$ to $\left\{A,B,...N\right\}$ diminishes. But I'm sure some smarter people have looked into this question.

If label switching is the problem, I'm not sure that either thinning or burn-in will fix it, since at any step, the sampler may decide to wander to the other mode, since will be equal probability of the alternative mode.

-

The problem here is that you are not thinning and burning your posterior sample. The first $N$ iterations in your MCMC simulation have to be removed (because the MCMC works asymptotically) and then, you need to obtain a subsample (every $k$ iterations) in order to reduce correlation/dependence (MCMC samples are correlated). The MCMC itself doesn't seem to be problematic in this specific case. See the following code, based on your own code:

library(MCMCpack)

# Simulated data
dat <- vector()

sigma=0.2

for(i in 1:1000){
if(runif(1)<0.7) dat[i] = rlnorm(1,meanlog=2,sdlog=sigma)
else dat[i] = rlnorm(1,meanlog=8,sdlog=sigma)
}

hist(dat)

d.lognorm.mix <- function(x,p=0.7,mu1,mu2,sigma=0.2) {
dens <- p*dlnorm(x,meanlog=mu1,sdlog=sigma) + (1-p)*dlnorm(x,meanlog=mu2,sdlog=sigma)
return(log(dens))
}

loglikeli <- function(theta) {  #theta = mu1, mu2
mu1 <- theta[1]
mu2 <- theta[2]
return(sum(d.lognorm.mix(dat, p=0.7, mu1=mu1, mu2=mu2, sigma=0.2)))
}

library(MCMCpack)
post.samp <- MCMCmetrop1R(loglikeli, theta.init=c(1,1),
thin=25, mcmc=60000, burnin=10000,
#     tune=c(1.5, 1.5),
verbose=FALSE, logfun=TRUE)

summary(post.samp)

# Posterior samples burned and thinned N=10000, k=25

hist(post.samp[,1])
hist(post.samp[,2])


Here is a link to a paper on the label switching problem for your joy:

A. Jasra, C. C. Holmes and D. A. Stephens. Markov Chain Monte Carlo Methods and the Label Switching Problem in Bayesian Mixture Modeling. Statist. Sci. Volume 20, Number 1 (2005), 50-67.

-
thanks for your input. I think I can just as well use "thin = 20" and "burnin = 5000" in the MCMCmetrop1R function – spore234 Aug 31 '14 at 18:37
@spore234 Indeed, I have updated my answer in order to account for this. – Sony Aug 31 '14 at 20:50