How to code a joint posterior $p(\theta | y,c) \propto p(c|y) p(y,c | \theta)$

Modeling
1

My posterior is a joint distribution of location(y) and concentration. c is concentration. The concentration at each point y obey exponential distribution. The probability of the location obey a mixture gaussian model. The model described below:

p(\theta | y,c) \propto p(c|y) p(y,c | \theta).
p(c|y)\sim exponential(1/f(y,c))
p(y,c|\theta)\sim mixture\space gaussian(\theta)
Likelihood(y,c|\theta) = \prod \space exp(1/f(y,c))
Now I write the code for

p(\theta | y) \propto p(\theta|y) p(y| \theta).
The likelihood function is

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);

How could I add the concentration information to the whole model and likelihood function. Now I give each point a transform parameter lamda. Lamda is caculated by \theta. But it gives me errors.

model='''
data {
 int D; //number of dimensions
 int K; //number of gaussians
 int N; //number of data
 vector[D] y[N]; //data
 real con[N]; //co ncentration

}

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
}

generated quantities {
real ro[K];
  matrix[D,D] Omega[K];
  matrix[D,D] Sigma[K];
  for (k in 1:K){
  Omega[k] = multiply_lower_tri_self_transpose(L[k]);
  Sigma[k] = quad_form_diag(Omega[k], sigma[k]); 
  ro[k]=Omega[k,2,1]; 
  }
}
transformed parameters{    
	cholesky_factor_cov[D,D] cov[K];  
  vector[N] lamba;
    for (k in 1:K) {cov[k] = diag_pre_multiply(sigma[k],L[k]);
    }
    for (i in 1:N){lamba[i]=1/(theta[1]*(1./2./Sigma[1,1,1]/Sigma[1,2,2]/sqrt(1-ro[1]*ro[1])*exp(-1./2./(1-ro[1]*ro[1])*((y[i,1]-mu[1,1])*(y[i,1]-mu[1,1])/Sigma[1,1,1]/Sigma[1,1,1])+(y[i,2]-mu[1,2])*(y[i,2]-mu[1,2])/Sigma[1,2,2]/Sigma[1,2,2]-2.*ro[1]*(y[i,1]-mu[1,1])*(y[i,2]-mu[1,2])/Sigma[1,1,1]/Sigma[1,2,2]))+theta[2]*(1./2./Sigma[2,1,1]/Sigma[2,2,2]/sqrt(1-ro[2]*ro[2])*exp(-1./2./(1-ro[2]*ro[2])*((y[i,1]-mu[2,1])*(y[i,1]-mu[2,1])/Sigma[2,1,1]/Sigma[2,1,1])+(y[i,2]-mu[2,2])*(y[i,2]-mu[2,2])/Sigma[2,2,2]/Sigma[2,2,2]-2.*ro[2]*(y[i,1]-mu[2,1])*(y[i,2]-mu[2,2])/Sigma[2,1,1]/Sigma[2,2,2])));
    }
}

model {
 real ps[K];
 theta ~ dirichlet(rep_vector(2.0, D));
 for(k in 1:K){
 mu[1,k] ~ normal(3.67,0.1);// uniform(340/100,380/100);//
 mu[2,k]~  normal(31.80,0.1);//uniform(3160/100,3190/100);//
 sigma[k] ~ exponential(0.5);//beta(2,5);//uniform(0.1,0.89); //////normal(1,2);
 L[k] ~ lkj_corr_cholesky(2);// contain rho 
 con ~ exponential(lamba);
 }
 
 for (n in 1:N){
 for (k in 1:K){
 ps[k] = log(theta[k])+multi_normal_cholesky_lpdf(y[n] | mu[k], cov[k]); //increment log probability of the gaussian
 }
 target += log_sum_exp(ps);
 }
   for(i in 1:N){
   target +=   - lamba[i]*con[i];
  }
 }

The error is:

ValueError: Failed to parse Stan model 'anon_model_e06a3361734d2035ca027d4b7fd124de'. Error message:
PARSER EXPECTED: whitespace to end of file.
FOUND AT line 28:
2

The generated quantities block should go at the end of your stan file and there is one too many curly brackets in the generated quantities.

3

Hi,
Thanks for replying. But when I put the generated quantities block at the end, the lamda line shows undefined since I am going to use Sigma[1,1,1] etc parameters in generated quantities. So I don’t know how to deal with it. I thinkcurly brackets is not a problem.

for (i in 1:N){lamba[i]=1/(theta[1]*(1./2./Sigma[1,1,1]/Sigma[1,2,2]/sqrt(1-ro[1]*ro[1])*exp(-1./2./(1-ro[1]*ro[1])*((y[i,1]-mu[1,1])*(y[i,1]-mu[1,1])/Sigma[1,1,1]/Sigma[1,1,1])+(y[i,2]-mu[1,2])*(y[i,2]-mu[1,2])/Sigma[1,2,2]/Sigma[1,2,2]-2.*ro[1]*(y[i,1]-mu[1,1])*(y[i,2]-mu[1,2])/Sigma[1,1,1]/Sigma[1,2,2]))+theta[2]*(1./2./Sigma[2,1,1]/Sigma[2,2,2]/sqrt(1-ro[2]*ro[2])*exp(-1./2./(1-ro[2]*ro[2])*((y[i,1]-mu[2,1])*(y[i,1]-mu[2,1])/Sigma[2,1,1]/Sigma[2,1,1])+(y[i,2]-mu[2,2])*(y[i,2]-mu[2,2])/Sigma[2,2,2]/Sigma[2,2,2]-2.*ro[2]*(y[i,1]-mu[2,1])*(y[i,2]-mu[2,2])/Sigma[2,1,1]/Sigma[2,2,2])));
    }
4

If you want to use the parameters Sigma, Omega and ro in the model or transformed parameter block, they should not go in the generated quantities but in the model or transformed parameter block.

5

Hi,

After I put those parameters to transformed parameter block, it works.
One thing I want to make sure is how could I let the model know the likelihood function is a exponential distribution? Also, how could I let the model know the posterior distribution is a exponential multiply mixture gaussian?

Is my present way right?

 for (n in 1:N){
 for (k in 1:K){
 ps[k] = log(theta[k])+multi_normal_cholesky_lpdf(y[n] | mu[k], cov[k]); //increment log probability of the gaussian
 }
 target += log_sum_exp(ps);
 }
   for(i in 1:N){
   target +=   - lamba[i]*con[i];
  }