# Obtaining posteriors for multivariate Normal mixture models

So I want to fit a mixture model

$$f(y) = \pi_1 f_1 (y) + \pi_2 f_2 (y)$$

where $$\pi_k = P(S = k)$$ and $$S_i$$ is a latent unobserved variable.

I assume that, conditional on $$S=k$$, we have the model

$$y_i = \mu_k + \beta_k x_i + \epsilon_i$$

where $$\epsilon_i \sim N_p (0, \Sigma_k)$$.

My goal is to estimate parameters $$\mu_k, \beta_k, \Sigma_k, \pi_k$$ for each $$k$$.

Firstly, As I understand it, under this model we have likelihood function

$$L(\theta |y) = \prod_i \left( \pi_1 \phi_{1,p} \left(y_i-\mu_1 -\beta_1 x_i\right) + (1-\pi_1) \phi_{2,p} \left(y_i-\mu_2 -\beta_2 x_i\right) \right)$$

where $$\phi_{k,p} \sim N_p(0, \Sigma_k)$$.

Secondly, in order to obtain point estimates of the parameters in a Bayesian framework, I need to obtain a posterior sample from the posterior distribution

$$p(\mu, \beta, \Sigma, \pi | y)$$

But as this is difficult, an alternative would be to use Gibbs sampling with posteriors

$$p(\mu | \beta, \Sigma, \pi , y)$$

$$p(\beta | \mu, \Sigma, \pi , y)$$

$$p(\Sigma | \mu, \beta, \pi , y)$$

$$p(\pi | \mu, \Sigma, \beta , y)$$

What I don't understand, is how to get from the given model equation and likelihood function to those posteriors. Once I have those posteriors, I am confident I can impliment Gibbs sampling.

Since my likelihood is not multivariate normal, it is a mixture of multivariate normals, I am unsure how to proceed, as there are no conjugates listed on the Wikipedia page for mixture distributions. Any help/advice is much appreciated!

You are quite correct there is no conjugate for the mixture model.

One way you can proceed. You have $$N$$ data points of $$K$$ dimensional multivariate data $$X$$ ($$N \times K$$) (I'm assuming you already think that there are 2 components? I'll call it $$M$$ here to be a bit more general.). You define a latent parameter $$Z$$ ($$N \times M$$) that is an assignment of each data point into one of the $$M$$ components. Initialize your $$\mu$$ ($$M \times K$$), $$\Sigma$$ ($$M \times K \times K$$)and $$\pi$$ ($$M \times K$$). I'm sorry, I didn't quite understand what purpose $$\beta$$ serves.

The basic idea is you assign each data point to one of the components, and conditional on these assignments, you can sample parameters for each of the $$M$$ components using the data assigned to it using conditional posteriors for a single multivariate normal.

Each Gibbs iteration consists of 2 steps:

1) Re-assign (re-sample) component membership, conditional on current parameter values using the probability that the $$j$$th data is assigned into the $$m$$th component:

$$p(z_{jm} = 1) = \frac{\pi_m N(x_j | \mu_m, \Sigma_m)}{\sum_l \pi_l N(x_j | \mu_l, \Sigma_l)}$$,

where $$N(..)$$ stands for the normal density.

2) Conditional on component membership, sample parameters for each component, using the conditional posteriors as you wrote above, except with it also being conditional on ONLY the data that are currently assigned to $$m$$th component. Your conditional posterior had it on ALL the data.

$$\mu_m \sim p(\mu_m | \Sigma_m, \pi_m, \{\forall j, x_j | z_{jm} = 1\})$$

Here, each component parameter sampling is independent of parameters of all other components.

Repeat.

========

There are ways to also infer # of components in a Bayesian manner, but this is a fairly simple algorithm that you can try first. For the priors, you can use standard Gaussian prior for each of the $$M$$ components in step 2.

This is roughly the procedure developed here:

https://www.jstor.org/stable/2345907?seq=1#page_scan_tab_contents

• The thing that confuses me is the model I'm supposed to fit specifically conditional densities $f(y|S=k) = x_i \beta_k + \epsilon_i$ where $\epsilon_i \sim N_p (\mu_k , \Sigma_k)$ and covariate $x_i$. Without this beta, I understand your answer well. But with it, I am completely lost as to how I'm supposed to obtain estimates for $\beta$. I think (?) this is called a Mixture Structural Equation Model, but tbh our lecture notes do not go into detail in the slightest and have no info on fitting it... I could use EM but then I'm maximising a function of 20ish parameters which seems unworldly. Oct 13 '18 at 16:44
• If I understand correctly, estimation of all parameters except $\mu$ and $\beta$ are the same regardless of whether $\beta$ is included or not. But, with the inclusion of $\beta$, neither $\mu$ nor $\beta$ are the mean of the normal densities, and so I'm not sure how to obtain congujate priors for each of them in this situation. Oct 13 '18 at 16:51
• What does the data look like when scatter plotted? x_i are the raw data points?
– ken
Oct 13 '18 at 16:51
• The data is simulated and follows that mixture model with the specified conditional densities. $x_i$ are generated covariates used to construct the response variables $y_i$ from known multivariate normals so that the model holds. It's just estimation that has me completely confused Oct 13 '18 at 16:52
• can you post a pic of your simulated data? I'm having a bit of hard time picturing what's going on.
– ken
Oct 13 '18 at 16:54