My model is a mixture model Y\sim pNorm(\mu_1,\sigma_1)+(1-p)Norm(\mu_2,\sigma_2), and I set
mu1=10;mu2=100;sd1=4;sd2=10;n=1000;p=0.3
with my stan code
... ...
parameters {
ordered[2] mu;
vector<lower=0>[2] sigma;
real<lower=0, upper=1> p;
}
model {
mu[1]~normal(120,1000);
mu[2]~normal(50,1000);
sigma ~ inv_gamma(0.1, 0.1);
p ~ beta(1, 1);
for (n in 1:N)
target += log_mix(p,
normal_lpdf(y[n] | mu[1], sigma[1]),
normal_lpdf(y[n] | mu[2], sigma[2]));
}
But my result is
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu[1] 10.14 0.01 0.31 9.50 9.93 10.15 10.35 10.75 561 1
mu[2] 10.26 0.01 0.29 9.68 10.06 10.26 10.45 10.80 714 1
sigma[1] 86.96 0.07 2.17 82.78 85.52 87.01 88.29 91.12 908 1
sigma[2] 3.36 0.01 0.21 2.98 3.21 3.35 3.50 3.78 816 1
p 0.76 0.00 0.02 0.73 0.75 0.76 0.77 0.79 743 1
lp__ -5526.89 0.09 1.62 -5530.72 -5527.75 -5526.59 -5525.70 -5524.68 331 1
If I re-run the same model, i will get another result
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu[1] 10.03 0.01 0.22 9.62 9.89 10.03 10.19 10.45 986 1
mu[2] 100.50 0.01 0.40 99.71 100.23 100.50 100.75 101.29 1575 1
sigma[1] 3.73 0.01 0.16 3.43 3.63 3.73 3.83 4.08 956 1
sigma[2] 10.42 0.01 0.27 9.89 10.22 10.42 10.61 10.96 1150 1
p 0.28 0.00 0.01 0.25 0.27 0.28 0.29 0.31 1315 1
lp__ -4065.13 0.09 1.62 -4069.11 -4065.92 -4064.79 -4063.97 -4063.05 346 1
DO SOMEONE KNOW THE REASON? THANKS!