Can I trust the MCMC results?

Modeling
1

Hi there,

A quite short question

I used a very simple model, Y=\tau_n+\tau_s+\tau_c+N(0,5), and I wish the three parameters following the dependent prior: \tau_n<\tau_s<\tau_c, so all posterior samples of the three parameters should obey the rule.

In fitting, I used three different Stan program scripts

Choice 1

functions{
  real mylpdf_lpdf(vector X, vector D, vector sigma_square, real tau_n, real tau_s, real tau_c){
    real lprob;
    if(tau_n<tau_s && tau_s<tau_c){lprob=normal_lpdf(X| D , sqrt(sigma_square) );   }
      else{
         lprob=negative_infinity();
          }
    return lprob;
  }
}
data{
    int<lower=1> N;
    vector[N] Y;   //responses
}
parameters{
    real<lower=0,upper=20> tau_n;
    real<lower=0,upper=20> tau_s;
    real<lower=0,upper=20> tau_c;
}
model{
    vector[N] mu;
    vector[N] sigma_square;
    for ( i in 1:N ) {
        sigma_square[i] =5;
    }
    for ( i in 1:N ) {
        mu[i] = tau_n+tau_s+tau_c;
    }
    Y ~ mylpdf( mu , sigma_square,tau_n,tau_s,tau_c );

Choice 2

functions{
    real dependent_flat_lpdf(vector y){   
            if(y[1]<y[2] && y[2]<y[3]) 
                return 1;  
            else   
                return negative_infinity();
    }
}
data{
    int<lower=1> N;
    vector[N] Y;   //responses
}
parameters{
    real<lower=0,upper=20> tau_n;
    real<lower=0,upper=20> tau_s;
    real<lower=0,upper=20> tau_c;
}
transformed parameters{
    vector[3] taus;
    taus[1]=tau_n;
    taus[2]=tau_s;
    taus[3]=tau_c;
}
model{
    vector[N] mu;
    vector[N] sigma_square;
    taus ~ dependent_flat();
    for ( i in 1:N ) {
        sigma_square[i] =5;
    }
    for ( i in 1:N ) {
        mu[i] = tau_n+tau_s+tau_c;
    }
    Y ~ normal( mu , sqrt(sigma_square) );
}

Choice 3

data{
    int<lower=1> N;
    vector[N] Y;   //responses
}
parameters{
    real<lower=0,upper=20> tau_n;
    real<lower=0,upper=20> tau_s;
    real<lower=0,upper=20> tau_c;
}
model{
    vector[N] mu;
    vector[N] sigma_square;
    for ( i in 1:N ) {
        sigma_square[i] =5;
    }
    for ( i in 1:N ) {
        mu[i] = tau_n+tau_s+tau_c;
    }
    if(tau_n<tau_s && tau_s<tau_c){
      Y ~ normal( mu , sqrt(sigma_square));}
    else{   target+=negative_infinity();             }
}

To my knowledge, the three Stan scripts are the same and will generate the similar results, and indeed they are. The ESS and \hat{R} and even traceplots are also good. The 95% credibility intervals contain the real true values as well. However, I was told
There were ``` divergent transitions after warmup. Increasing adapt_delta above 0.85 may help.
and the pairs() shows

pairs
pairs1080×690 117 KB

Can I trust the MCMC results? If not, how can I improve the model without the change to the assumption \tau_n<\tau_s<\tau_c?

Thanks in advance for any tips and kind help.

2

instead of sampling the taus, sample tau1, the increment to tau2, and the increment to tau3

make these increments positive parameters… voila

1 Like
3

Use positive_ordered:

data {
    int <lower = 1> N;
    vector [N] Y;   //responses
}

transformed data {
    real <lower = 0> sigma_square = 5;
    real <lower = 0> sigma = sqrt(sigma_square);
}

parameters {
    positive_ordered [3] tau;
}

model {
    vector [N] mu = rep_vector(tau[1] + tau[2] + tau[3], N);
    Y ~ normal(mu, sigma);
}

I don’t think this model is identified though. Why are the parameters constrained to be less than 20?

1 Like
4

No. The link provided in the message

https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

has more details. In this case, you are going to need a different parameterization. Anytime there is an if (...) condition that depends on unknown variables, the sampler will struggle. In this case, you should be able to do just fine with

parameters {
  positive_ordered[3] tau;
}

and then indexing the three elements of tau as \tau_n, \tau_s, \tau_c.

2 Likes
5

Thanks for your help!

This is just an assumption, which is out of curiosity, but I think the assumption may make the coding messy and discourage me to use increment for tau.

By the way, if I use positive_ordered [3] tau in Stan, what is the corresponding (default) joint prior of them? Is it flat/uniform anywhere in the region tau[1]<tau[2]<tau[3]?

6

Yes (it is uniform over the y_k values), but this is not a very intuitive space to reason about.