Efficiently sampling from mixture distribution posterior

Question

I have the following model: $$ \begin{align} \pi_1\sim & \text{Unif}(0,1)\\ \lambda_1,\lambda_2\sim & \text{Ga}(1,1)\\ z_i\sim & \pi_1^{1(z_i=1)}\pi_2^{1(z_i=2)}\\ p(y_i|\lambda_1,\lambda_2,z_i)=&[\lambda_1\exp(-\lambda_1 y_i)]^{1(z_i=1)}[\lambda_2\exp(-\lambda_2 y_i)]^{1(z_i=2)} \end{align} $$ where $\pi_2=1-\pi_1$.

The full conditionals that I used for the Gibbs samplers are as follows: $$ \begin{align} p(\lambda_j|y,z,\pi)= & \text{Ga}\left(1+n_j,1+\sum_{i:z_i=j} y_i\right)\\ p(z_i=j|\lambda,\pi,y)= & \frac{1}{C} \lambda_j\exp(-\lambda_j y_i)\pi_j\\ p(\pi) = & \text{Dir}(1+n_1,1+n_2), \end{align} $$ where $n_j$ is the number of observations currently assigned to cluster $j$, $C$ is the normalising constant of a categorical distribution such that $p(z_i=1|\lambda,\pi,y)+p(z_i=2|\lambda,\pi,y)=1$ and $\text{Dir}$ is a Dirichlet distribution.

For data generated from N=250, $\lambda=[2, 6]$ and $\pi=[0.4, 0.6]$ and using a Gibbs sampler, the joint posterior that I get on $\lambda_1,\lambda_2$ are as follows:

My question is considering that after I marginalise out all $z_i$, $$p(y)=\int \sum_{j=1}^2 \pi_j\lambda_j\exp(-\lambda_j y_i) d\lambda_1 d\lambda_2,$$ I would expect the posterior on $\lambda$ to be symmetric about $\lambda_1=\lambda_2$. If I were to take these $\lambda$ samples and swap their dimensions would it give me the correct posterior?

I'm assuming that considering the relatively high number of samples, the Gibbs sampler gets stuck in local modes.

Answer 1

You are correct that there is an unidentifiability in mixture models due to the problem of label-switching, i.e. you get the same likelihood if the mixture components are labeled 1,2 or 2,1. For convergence in a Gibbs sampler, you need the labels to swap back and forth. But if the modes are well separated, then it can take a long time for the labels to swap.

In univariate problems like yours, the easiest approach is to enforce identifiability in the prior by ordering the mean parameters. So, you could change your prior on $\lambda$ to be $$ p(\lambda) \propto Ga(\lambda_1;1,1) Ga(\lambda_2;1,1) 1(\lambda_1<\lambda_2) $$ which modifies the full conditional distributions to be $$ \lambda_j|\lambda_{-j},y,z \sim Ga\left(1+n_j,1+\sum_{i:z_i=j} y_i\right)1(\lambda_1<\lambda_2). $$ Note that the $\lambda_j$ are no longer conditionally independent.