Mixture of Two Normals

Question

Suppose we have a data which consists of two normals,

x = rnorm(50,mean=1,sd=2)
y = rnorm(50,mean=2,sd=3)
z = sample( c(x,y) , size = 100, replace=FALSE )

The goal is to estimate the two normal parameters and then classify which of the two normal families each of the data points belong to.

Here is an approach to this problem that I was thinking of. Let $p = (\mu_1,\sigma_1,\mu_2,\sigma_2)$ denote the 4-dimensional vector of real numbers, and $I$ a binary string of length $100$, i.e. for example, $I = (1,1,...,1,0,0,..0)$, here $I$ is categorical and indicates which normal family each data point belongs to. We can define the following likelihood function,

L = function(p,I){
likelihood = 0
         for(i in 1:100){
         likelihood = likelihood + dnorm(z[i],p[1],p[2])*I[i] + dnorm(z[i],p[3],p[4])*(1 - I[i]) 
} }

Now that this function has been defined it remains to sample from it. It is possible to write a custom Metropolis-algorithm here, however, I am not sure how efficient it would be, and wondered if there is a way to do this using Stan? According to the Stan manual, the software does not allow discrete parameters, so I would be unable to define a binary parameter $I$ as a vector of $100$ components. Is there a workaround to this that will allow Stan to sample from this high-dimensional distribution?

Answer 1

If I remember correctly, in cases where Stan does not allow you to use discrete parameters, the most straightforward solution is to marginalize over the latent discrete variable. I am quoting the pymc3 documentation (from: https://pymc3-testing.readthedocs.io/en/rtd-docs/notebooks/marginalized_gaussian_mixture_model.html)

---

A natural parameterization of the Gaussian mixture model is as the latent variable model:

$$ \begin{split}\begin{align*} \mu_1, \ldots, \mu_K & \sim N(0, \sigma^2) \\ \tau_1, \ldots, \tau_K & \sim \textrm{Gamma}(a, b) \\ \boldsymbol{w} & \sim \textrm{Dir}(\boldsymbol{\alpha}) \\ z\ |\ \boldsymbol{w} & \sim \textrm{Cat}(\boldsymbol{w}) \\ x\ |\ z & \sim N(\mu_z, \tau^{-1}_z). \end{align*}\end{split} $$

An alternative, equivalent parameterization that addresses these problems is to marginalize over z. The marginalized model is:

$$ \begin{split}\begin{align*} \mu_1, \ldots, \mu_K & \sim N(0, \sigma^2) \\ \tau_1, \ldots, \tau_K & \sim \textrm{Gamma}(a, b) \\ \boldsymbol{w} & \sim \textrm{Dir}(\boldsymbol{\alpha}) \\ f(x\ |\ \boldsymbol{w}) & = \sum_{i = 1}^K w_i\ N(x\ |\ \mu_i, \tau^{-1}_z), \end{align*}\end{split} $$

Marginalizing $z$ out of the model generally leads to faster mixing and better exploration of the tails of the posterior distribution. Marginalization over discrete parameters is a common trick in the Stan community, since Stan does not support sampling from discrete distributions.

---

I think this post might be useful: https://jeremy9959.net/Blog/StanMixture/