Hi,
Thanks for the help. Your way works! I still carious about the meaning of the parameters. My question is how can I get the distribution of correlation?
model='''
data {
int D; //number of dimensions
int K; //number of gaussians
int N; //number of data
vector[D] y[N]; //data
}
parameters {
simplex[K] theta; //mixing proportions
ordered[D] mu[K]; //mixture component means
cholesky_factor_corr[D] L[K]; //cholesky factor of correlation matrix
vector<lower=0>[D] sigma[K]; // standard deviations
}
transformed parameters{
cholesky_factor_cov[D,D] cov[K];
for (k in 1:K) cov[k] = diag_pre_multiply(sigma[k],L[k]);
}
// real sx[K];
// real sy[K];
// real ro[K];
// cov=L;
// for (i in 1:K){
// sx[i]=sqrt(cov[i,1,1]);
// sy[i]=sqrt(cov[i,2,2]);
// ro[i]=cov[i,1,2]/sqrt(cov[i,1,1])/sqrt(cov[i,2,2]);
// }
//}
model {
real ps[K];
theta ~ dirichlet(rep_vector(2.0, D));
for(k in 1:K){
mu[1,k] ~ normal(3.67,1);//uniform(340/100,380/100);//
mu[2,k]~ normal(31.80,1);//uniform(3160/100,3190/100);//
sigma[k] ~ normal(1,2);
L[k] ~ lkj_corr_cholesky(2);
}
for (n in 1:N){
for (k in 1:K){
ps[k] = log(theta[k])+multi_normal_cholesky_lpdf(y[n] | mu[k], L[k]); //increment log probability of the gaussian
}
target += log_sum_exp(ps);
}
}
'''
Result:
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
theta[1] 0.85 1.6e-3 0.1 0.62 0.79 0.87 0.92 0.98 3731.0 1.0
theta[2] 0.15 1.6e-3 0.1 0.02 0.08 0.13 0.21 0.38 3731.0 1.0
mu[1,1] 1.14 9.2e-3 0.34 0.55 0.93 1.13 1.34 1.75 1345.0 1.0
mu[2,1] 31.22 0.01 0.85 29.51 30.67 31.23 31.81 32.83 3741.0 1.0
mu[1,2] 29.19 4.1e-3 0.3 28.61 28.99 29.2 29.4 29.78 5352.0 1.0
mu[2,2] 32.38 0.01 0.81 30.78 31.85 32.38 32.92 33.99 5504.0 1.0
L[1,1,1] 1.0 nan 0.0 1.0 1.0 1.0 1.0 1.0 nan nan
L[2,1,1] 1.0 nan 0.0 1.0 1.0 1.0 1.0 1.0 nan nan
L[1,2,1] 0.98 3.2e-3 0.06 0.93 0.99 0.99 1.0 1.0 358.0 1.0
L[2,2,1] -4.0e-3 5.7e-3 0.45 -0.81 -0.35 5.9e-4 0.33 0.82 6045.0 1.0
L[1,1,2] 0.0 nan 0.0 0.0 0.0 0.0 0.0 0.0 nan nan
L[2,1,2] 0.0 nan 0.0 0.0 0.0 0.0 0.0 0.0 nan nan
L[1,2,2] 0.13 4.5e-3 0.09 0.06 0.09 0.11 0.15 0.36 433.0 1.01
L[2,2,2] 0.88 3.7e-3 0.14 0.49 0.83 0.94 0.99 1.0 1415.0 1.0
sigma[1,1] 2.01 0.03 1.41 0.1 0.89 1.8 2.88 5.36 2978.0 1.0
sigma[2,1] 2.01 0.02 1.39 0.11 0.89 1.76 2.87 5.24 3984.0 1.0
sigma[1,2] 1.99 0.03 1.45 0.07 0.83 1.75 2.89 5.35 3098.0 1.0
sigma[2,2] 2.03 0.03 1.43 0.08 0.89 1.81 2.91 5.42 3054.0 1.0
cov[1,1,1] 2.01 0.03 1.41 0.1 0.89 1.8 2.88 5.36 2978.0 1.0
cov[2,1,1] 2.01 0.02 1.39 0.11 0.89 1.76 2.87 5.24 3984.0 1.0
cov[1,2,1] 1.96 0.03 1.43 0.06 0.8 1.72 2.84 5.27 2911.0 1.0
cov[2,2,1] -0.01 0.02 1.11 -2.58 -0.5 2.1e-4 0.47 2.43 4070.0 1.0
cov[1,1,2] 0.0 nan 0.0 0.0 0.0 0.0 0.0 0.0 nan nan
cov[2,1,2] 0.0 nan 0.0 0.0 0.0 0.0 0.0 0.0 nan nan
cov[1,2,2] 0.26 0.01 0.3 8.1e-3 0.09 0.2 0.35 0.91 903.0 1.0
cov[2,2,2] 1.79 0.02 1.31 0.08 0.75 1.56 2.57 4.92 2908.0 1.0
lp__ -371.9 0.09 2.82 -378.4 -373.5 -371.5 -369.8 -367.5 1069.0 1.0