Issues in my stan code

Hi guys,

I am new to stan, there are a lot of problems when I ran the code below, I have no idea what is wrong.
Could you please help with this? Thanks a lot

modelString1 = "
parameters {
  real<lower=0> theta; 
  real<lower=0> sigma;
}
model {
  theta ~ normal (0.5,0.1) ;     
  sigma ~ normal (1, 0.3) ;  
}
generated quantities {
  real<lower=1,upper=100> n ;
  real<lower=0,upper=1> poweraj ;
  real<lower=0,upper=1> pcrit;
  n= 2 * pow((.84 +1.96) * sigma / theta, 2);
  poweraj= Phi( sqrt(31.5)* theta/sigma  -1.96  );
  pcrit= 1- poweraj>=0.7;
}
"
#initsList =  list( theta = 0.5, sigma = 1*0.9+0.1)
writeLines( modelString1 , con="TEMPmodel.stan" )
stanModel <- stan(
  file = "TEMPmodel.stan",  
  chains = 4,            
  warmup = 5000,         
  iter = 25000,
  control = list(adapt_delta = 0.99)
)
[[1]]
Stan model 'TEMPmodel' does not contain samples.

[[2]]
Stan model 'TEMPmodel' does not contain samples.

[[3]]
Stan model 'TEMPmodel' does not contain samples.

Warning messages:
In .local(object, ...) :
  some chains had errors; consider specifying chains = 1 to debug

Hi and welcome. I cleaned up your post so it was easier to read. By using ``` to enclose blocks of model or code it will make it easier to read.

Thanks. But do you have any clue about how to slove the issues. Many thanks.

If I run your model I get the errors:

n is 111.684, but must be less than or equal to 100

This is because your constraints in the generated quantities block do not match the parameters being passed in. Using constraints in generated quantities does not influence the sampling of parameters, it just serves as a means of input-checking. If you want n to be constrained to a given range, then you need to appropriately constrain either theta or sigma.

Your current constraint for n is:

1 < 2*(\frac{2.8*\sigma}{\theta})^2 < 100

Double-check that’s what you intended. If it is, then we can rearrange/simplify to give the appropriate constraint for sigma.

Divide all terms by 2:

0.5 < (\frac{2.8*\sigma}{\theta})^2 < 50

Take the square root of all terms:

\sqrt{0.5} < \frac{2.8*\sigma}{\theta} < \sqrt{50}

Multiply all terms by theta:

\sqrt{0.5}*\theta < 2.8*\sigma < \sqrt{50}*\theta

Divide all terms by 2.8:

\frac{\sqrt{0.5}*\theta}{2.8} < \sigma < \frac{\sqrt{50}*\theta}{2.8}

If I then update your model with this constraint:

modelString1 = "
parameters {
  real<lower=0> theta; 
  real<lower=((sqrt(0.5)*theta)/2.8),
       upper=((sqrt(50)*theta)/2.8)> sigma;
}
model {
  theta ~ normal (0.5,0.1) ;     
  sigma ~ normal (1, 0.3) ;  
}
generated quantities {
  real<lower=1,upper=100> n ;
  real<lower=0,upper=1> poweraj ;
  real<lower=0,upper=1> pcrit;
  n= 2 * pow((.84 +1.96) * sigma / theta, 2);
  poweraj= Phi( sqrt(31.5)* theta/sigma  -1.96  );
  pcrit= 1- poweraj>=0.7;
}
"

Everything samples without issue

2 Likes

Thanks very much! It is clear. But I have the other question regarding this,
actually I have jags model like this: theta ~ dnorm(0.5,1/(0.1^2))I(0,)
sigma ~ dnorm(1,1/(0.3^2))I(0,) . It looks like the results still off from the Jags.

I am think if I can do this below. but looks like a lot of issues raises again.
modelString1 = "
parameters {
real theta;
real sigma;
real<lower=0> L;
}
model {
theta ~ normal (0.5,0.1) T[L,] ;
sigma ~ normal (1, 0.3) T[L,];
}
generated quantities {
real<lower=1> n ;
real<lower=0,upper=1> poweraj ;
real<lower=0,upper=1> pcrit;
n= 2 * pow((.84 +1.96) * sigma / theta, 2);
poweraj= Phi( sqrt(31.5)* theta/sigma -1.96 );
pcrit= 1- poweraj>=0.7;
}
"

In your JAGS model you’re truncating at 0, but your Stan model is truncating at a random variable in the the range [0,\infty]. For the models to be the same you need to truncate at 0:

theta ~ normal (0.5,0.1) T[0,] ;
sigma ~ normal (1, 0.3) T[0,];
1 Like

Really appreciate your help. I got the other two warnings below, do you know how to slove those? Sorry, this is simple code with many errors…But it can help me understand how stan work out…

options(mc.cores = parallel::detectCores())

modelString1 = "
parameters {
real theta;
real sigma;
}
model {
theta ~ normal (0.5,0.1) T[0,] ;
sigma ~ normal (1, 0.3) T[0,];
}
generated quantities {
real n ;
real poweraj ;
real pcrit;
n= 2 * pow((.84 +1.96) * sigma / theta, 2);
poweraj= Phi( sqrt(31.5)* theta/sigma -1.96 );
pcrit= poweraj<0.7;
}
"
writeLines( modelString1 , con=“TEMPmodel.stan” )
stanModel ← stan(
file = “TEMPmodel.stan”,
chains = 4,
warmup = 5000,
iter = 25000,
control = list(adapt_delta = 0.99)
)

2: There were 178 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
3: Examine the pairs() plot to diagnose sampling problems

When using truncation you also need to have the respective bound specified for the parameter, otherwise you can run into sampling issues (more information in this section of the User’s Guide).

Try the following model:

parameters {
  real<lower=0> theta;
  real<lower=0> sigma;
}
model {
  theta ~ normal (0.5,0.1) T[0,];
  sigma ~ normal (1, 0.3) T[0,];
}
generated quantities {
  real n;
  real poweraj;
  real pcrit;
  n= 2 * pow((.84 +1.96) * sigma / theta, 2);
  poweraj= Phi( sqrt(31.5)* theta/sigma -1.96 );
  pcrit= poweraj<0.7;
}
3 Likes

Thank you very much!!!

1 Like