Hi Ben,
Since I want to know the distribution of rho in that matrix 
. I write below code but it says initialization failure. So do you mean it is impossible to know the distribution of factor rho? Thanks.
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{
real ro[K];
cholesky_factor_cov[D,D] cov[K];
for (k in 1:K) {cov[k] = diag_pre_multiply(sigma[k],L[k]);
ro[k]=cov[1,2,k]/sqrt(cov[1,1,k])/sqrt(cov[2,2,k]);}
}
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,0.1);//uniform(3160/100,3190/100);//
sigma[k] ~ uniform(0.01,0.89);//normal(1,2);
ro[k] ~ uniform(-0.999,0.999);
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);
}
}
'''