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!
 A: 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
