How to reparameterize a Gaussian mixture model?

Modeling
12

Maybe a non-centered parameterization of the mixture components works better?

// Guassian mixture model
data {
  int<lower=1> N;        // number of data points
  real y[N];             // observed values
  real uncertainties[N]; // uncertainties
}
parameters {
  real<lower=0,upper=1> p;   // mixing proportion
  ordered[2] mu;             
  real<lower=0> sigma;       // spread of mixture components
  vector[N] z1;
  vector[N] z2;
}
transformed parameters {
  vector[N] y_mix1 = mu[1] + sigma*z1;
  vector[N] y_mix2 = mu[2] + sigma*z2; // I corrected this in an edit!!
}
model {
  sigma ~ normal(0, 0.1);
  mu[1] ~ normal(-1.4, 1);
  mu[2] ~ normal(-1.0, 1);
  p ~ beta(2, 2);
  
  z1 ~ std_normal();
  z2 ~ std_normal();

  for (n in 1:N) {
    target += log_mix(p,
                      normal_lpdf(y[n] | y_mix1[n], uncertainties), // corrected (see edit2 below)
                      normal_lpdf(y[n] | y_mix2[n], uncertainties)); // corrected (see edit2 below)
  }
}
generated quantities{
  vector[N] y_true;
  
  for (n in 1:N)
    y_true[n] = bernoulli_rng(p) == 0 ? y_mix1[n] : y_mix2[n];

}

Can’t make anything of the results though… Do they make sense?

Edit: I corrected something in the code above (see comment in code): the auxillary variable for y_mix2 should be be z2…

Edit2: As pointed out by @Evgenii, y_mix1 and y_mix2 need to be indexed as well in the log_mix call \rightarrow changed to y_mix1[n] and y_mix2[n].

4 Likes