How to write generated quantities syntax from a joint density function

Modeling
1

Hi,

I have a model syntax that successfully compiles and fits simulated data with no issues. I am trying to add the generated quantities block to compute model fit indices with the loo package.

Does anyone have any suggestions? Thank you.

The model syntax is below. Please ignore any redundancy in the code pointed out by @Bob_Carpenter in another post. I will deal with it later.

data{
  int  I;                 //  number of items
  int  P;                 //  number of individuals
  int  n_obs;             //  number of observed responses
  int  item[n_obs];       //  item id
  int  id[n_obs];         //  person id
  real Y[n_obs];          //  vector outcome
}

parameters {
  real beta[I];         // vector of beta parameters for
  real ln_alpha[I];     // vector of alpha parameters for I items     
  real omega[I];        // vector of omega parameters for I items     
  real gamma0[I];       // vector of gamma0 parameters for I items     
  real gamma1[I];       // vector of gamma1 parameters for I items     
  real theta[P];        // vector of theta parameters for P individuals      
}

transformed parameters{
  real<lower=0> alpha[I];
  for(i in 1:I){
       alpha[i] = exp(ln_alpha[i]);
  } 
}
   
model{
    
  beta      ~ normal(0,10);
  ln_alpha  ~ normal(0,10);
  omega     ~ normal(0,10);
  gamma0    ~ normal(0,10);
  theta     ~ normal(0,1);      
  
  for(i in 1:I){
    gamma1[i] ~ normal(0,10)T[gamma0[i],];
  }
  
  for(i in 1:n_obs) {
    
    real   k0 = inv_logit(gamma0[item[i]] - alpha[item[i]]*theta[id[i]]);
    real   k1 = inv_logit(gamma1[item[i]] - alpha[item[i]]*theta[id[i]]);
  
    if (Y[i] == 0) {
      target += log(k0);
    } else if (Y[i] == 1) {
      target += log1m(k1);                    
    } else {
      
      real  sh1 = exp((alpha[item[i]]*theta[id[i]] + beta[item[i]] + omega[item[i]])/2);
      real  sh2 = exp((-(alpha[item[i]]*theta[id[i]] + beta[item[i]]) + omega[item[i]])/2);

      target += log(k1 - k0) + beta_lpdf(Y[i] | sh1,sh2);   
    }	

  }
}
1 Like
2

Can you be more specific on what the issue is? I am not sure what the loo package requires (or what is the package, for that matter), but you can compute any quantity in the generated quantities block, and since it use samples from the posterior it will be from a joint density (as your question states). So at taken at face value your answer is simply to put whatever function you need to retrieve (that is not already in the transformed parameters) in that block and use it for whatever you need outside of the Stan program.