Issues in my stan code

Modeling
,
1

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
2

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.

3

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

4

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
5

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;
}
"

6

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
7

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

8

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
9

Thank you very much!!!

1 Like