data { int N; vector[N] y; } parameters { ordered[2] mu; real sigma[2]; real theta; } model { sigma ~ normal(0, 2); mu ~ normal(0, 2); theta ~ beta(5, 5); for (n in 1:N) target += log_mix(theta, normal_lpdf(y[n] | mu[1], sigma[1]), normal_lpdf(y[n] | mu[2], sigma[2])); }