Rstan is running too slowly, seeking advice on ways to improve speed.

[huiying]

Dear people,

I'm rather new to R and very new to Rstan.

I'm attempting to reproduce someone else's work. However having troubles with rstan that runs very very slowly

This is one example of the time taken to complete one chain.

#  Elapsed Time: 541.943 seconds (Warm-up)

#                552.663 seconds (Sampling)

#                1094.61 seconds (Total)

The attached syntax is meant for one age group.

I have to repeat the same procedures for another 4 age groups, each with 5-chains with 60000 iterations. That means it will take more than 8 hours to run the full syntax.

The syntax was originally done by someone else, and I don't think they had the same problems. The syntax i'm trying to run is almost identical to theirs, except changing the following command line 

for (i in 1:W)

y[i] ~ binomial(n[i], pi[i]) 

to y ~ binomial(n, pi)

Please kindly advise if there's any error I've overlooked or just any comments are welcomed.

With thanks,

Huiying

[Krzysztof Sakrejda]

On Tuesday, April 19, 2016 at 6:14:37 AM UTC-4, huiying wrote:

Dear people,

I'm rather new to R and very new to Rstan.

I'm attempting to reproduce someone else's work. However having troubles with rstan that runs very very slowly

This is one example of the time taken to complete one chain.

#  Elapsed Time: 541.943 seconds (Warm-up)

#                552.663 seconds (Sampling)

#                1094.61 seconds (Total)

The attached syntax is meant for one age group.

I have to repeat the same procedures for another 4 age groups, each with 5-chains with 60000 iterations.

60000 iterations is likely very excessive.  Have you tried seeing if using 1000 iterations would be enough?  

Krzysztof

[Bob Carpenter]

Also, it's impossible to diagnose a Stan program
with just a single line.

You can run chains in parallel. We usually run 4
chains.

And rather than running multiple models by age group, have you
thought of making some kind of hierarchical model and putting
them all together?

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[huiying]

Hey ! Thanks for taking time to reply.

When running in parallel, the time taken is around the same as that for the first chain.

Anyway, I've tried to run 1000 iterations with 4 chains. It runs in seconds but it gave me the following error. 

Warning messages:

1: There were 1964 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. 

2: Examine the pairs() plot to diagnose sampling problems

Increasing max_treedepth to 15 took a good 20 minutes or so for 4 chains and isn't good enough. I got similar warning message to increase max_treedepth above 15. I would think that increasing max_treedepth further will take even longer time to run.

Is it correct to interpret this as non-convergence and can I increase number of iterations/chains instead of increasing max_treedepth?

How can I decide how many iterations and chains is sufficient?

Yes, thought of doing a hierarchical model, but haven't figured out exactly how to do that and I was thinking to start off with one age group before writing out a larger syntax, and not knowing the source of error when any error occurs.

I'm also concerned that a hierarchical model be more time-consuming than running it separately for each age group, but I guess I can worry about this later.

Huiying

[James Camac]

Exceeding the maximum tree depth doesn't mean the model hasn't converged. It just means the sampler becomes less efficient as when max tree depth is reached it reverts to a random walk.

To check convergence you can look at the traceplots with rstan.

Another option (ideally you do both) is to look at the Rhats for each parameter. They need to be less than 1.1 to indicate convergence.

The other thing to look at is the effective sample size you are getting for each parameter. I tend to ensure that each parameter has at least 100 samples.

You might also want to consider decreasing the initial stepsize from the default of 1 to something smaller e.g. 0.5 (control(list(stepsize= 0.5)), especially if you are getting a lot of numerical issues during warmup.

I'm guessing you're not getting divergent iterations based on the error messages you're receiving.

Sorry haven't looked at your model.

[huiying]

Hi James,

Thanks for the helpful explanation and guidelines.

With 10000 iterations and 4 chains, the output for Rhats was > 1.1 and effective sample size for some samples were < 100. Traceplots clearly showed evidence of divergence.

With 1000 iterations, it is even more poorly converged.

Do I increase the number of iterations until the output shows desired Rhats and ESS? Is that the proper way of conduct? 

With 20000 iterations, 4 chains, Rhats were less than 1.1 and ess were above 100 (at borderline for some parameters), however I still get warning message concerning max_treedepth. Would it be appropriate to ignore the warning message? (I doubt so?)

However, with further increased # of iterations (tried with 30,000 and 40,000), some parameters had very low ESS, while some had way above 100. 

What are the possible reasons for increased number of iterations required to converge and also for exceeding max treedepth?

Can rstan run on multiple thread (in case I really have to run 60,000 iterations!)?

I've also tried decreasing initial stepsize to 0.5, but that did not seem to improve performance. What does this do actually?

:(

Huiying

[Krzysztof Sakrejda]

Could you explain the purpose of the replication here? Are you currently working with the same model and the same data as it was originally applied to?  Simulated data? Another data set?  I think you already have ample evidence that this model will not work well for this data (comments on this below) and rather than fiddling with sampler parameters the density needs to be reparameterized.

On Wednesday, April 20, 2016 at 2:57:44 AM UTC-4, huiying wrote:

Hi James,

Thanks for the helpful explanation and guidelines.

With 10000 iterations and 4 chains, the output for Rhats was > 1.1 and effective sample size for some samples were < 100. Traceplots clearly showed evidence of divergence. 

With 1000 iterations, it is even more poorly converged.

If 10000 iterations doesn't fix a problem visible at 1000 iterations you will likely never run enough iterations to get good results out of this model.

 

Do I increase the number of iterations until the output shows desired Rhats and ESS? Is that the proper way of conduct? 

Fix the model.  Look for strongly correlated parameters,  poorly identified sets of parameters, or parameters with correlations changing as a function of other parameters, look for parameters on boundaries, these are all things you can fix but they need to be found first.  

 

With 20000 iterations, 4 chains, Rhats were less than 1.1 and ess were above 100 (at borderline for some parameters), however I still get warning message concerning max_treedepth. Would it be appropriate to ignore the warning message? (I doubt so?)

Treedepth limits how many steps the integrator running HMC will take before giving up (whether it finds a u-turn or not).  I'm sure there are models where it's ok and you can live with the reduced efficiency but I haven't run into one yet.  Typically this means that in at least one direction in one portion of parameter space the scale the sampler needs to move on is much larger than whatever part of the space it was tuned in (during adaptation).  If it's not a practically relevant part of the parameter space then you can make it go away by setting tighter priors, or transforming parameters or something similar.

 

However, with further increased # of iterations (tried with 30,000 and 40,000), some parameters had very low ESS, while some had way above 100. 

The 20k iteration chains missed pathologies in the model---not surprising since they had less time to find them.

 

What are the possible reasons for increased number of iterations required to converge and also for exceeding max treedepth?

Can rstan run on multiple thread (in case I really have to run 60,000 iterations!)?

You could run many chains in parallel but each one would still have to run until convergence so this is not likely to be helpful.

[Krzysztof Sakrejda]

So in the model, the -log(W) terms don't matter (they're constants), but taking a parameter that's the min of other parameters is confusing (and likely to cause problems).  What's this supposed to be a model fror?

On Tuesday, April 19, 2016 at 6:14:37 AM UTC-4, huiying wrote:

[James Camac]

I had a quick look at the model. I have to admit I'm a little lost at what you are trying to do.

I agree with Kryzysztof, you'll need to consider a reparameterisation of the model and that just upping the number of iterations is not going to help here.

rstan's pairs plot function might help with diagnosing where the issue is.

A few extra notes:

- You don't need to set a parameters lower bound as something like 0.0000015625 just use zero.

- You specify theta0 as a bounded parameter between 0 and 1, and the prior you use is beta(0.5, 1129.5). This prior suggests you expect the mean posterior of this parameter to be very very close to zero. I could be wrong... but I think this could lead to issues as the sampler will struggle to sample so close to a boundary.

[Bob Carpenter]

> On Apr 20, 2016, at 6:02 PM, James Camac <james...@gmail.com> wrote:
>
> ...

> - You specify theta0 as a bounded parameter between 0 and 1, and the prior you use is beta(0.5, 1129.5). This prior suggests you expect the mean posterior of this parameter to be very very close to zero. I could be wrong... but I think this could lead to issues as the sampler will struggle to sample so close to a boundary.

Yes, this can be a problem.

- Bob

[huiying]

Dear all,

Thanks very much for taking time to address all my questions.

I am reproducing the results, with an intention to further improve the model. Before I could do that, I obviously have lots to learn. Yes, I am working on the same model and same data as it was originally applied to.

So my understanding from the discussion is that the definition for the parameters is the main issue, and I will work on it next. 

I have to admit that I'm not very certain with the command at the last part of the model

  Pr_S <- binomial_coefficient_log(est_in, S)

  + multiply_log(S, est_is[1])

  + multiply_log((est_in - S), est_is[2]); 

  increment_log_prob((Pr_S));

  increment_log_prob(-log(W));

  Pr_D <- binomial_coefficient_log(est_in, D)

  + multiply_log(D, est_if[1])

  + multiply_log((est_in - D), est_if[2]); 

  

  increment_log_prob((Pr_D));

  increment_log_prob(-log(W));

I tried to google for binomial_coefficient_log but had limited search results on forum/mailing list. I am going to ask the person who did this previously, but here is what I think it is which may not be correct,

Pr_D is the probability of death, and will take the value est_in which is estimated from posterior distribution of pi from the same model.

binomial_coefficient_log is used in view of unknown est_in and also non-informative/weak prior of est_if (risk of death), defined as having a beta distribution (0.5, 0.5) which is rather random?

I'm sorry if this doesn't make sense.

Could I have used the number of death, D ~ binomial(est_in, est_if)?

theta gets a lower bound 0.0000015625 based on the knowledge that 20 cases were observed in the population of size 12.8e6.

theta0 takes a prior beta(0.5, 1129.5) based on the knowledge that the 0 out of 1129 samples was found positive in a previous study and alpha for beta distribution has to be >0. It is also set to bound between 0 and 1 because it is a parameter for baseline cumulative incidence, which is a proportion.

if prior for pi is known to be beta(0.5, 193.5), can I incorporate this into the model when at the same time, pi is defined as theta0 + theta*x1?

Thank you for your patience.

Huiying

[Bob Carpenter]

> On Apr 20, 2016, at 10:42 PM, huiying <hyc.t...@gmail.com> wrote:
>
> Dear all,
>
> Thanks very much for taking time to address all my questions.
>

> I am reproducing the results, with hope to further improve the model. Before I could do that, I obviously have lots to learn. Yes, I am working on the same model and same data as it was originally applied to.

>
> So my understanding from the discussion is that the definition for the parameters is the main issue, and I will work on it next.
>
> I have to admit that I'm not very certain with the command at the last part of the model

Where did it come from?

> Pr_S <- binomial_coefficient_log(est_in, S)
> + multiply_log(S, est_is[1])
> + multiply_log((est_in - S), est_is[2]);
>
> increment_log_prob((Pr_S));
> increment_log_prob(-log(W));
>
> Pr_D <- binomial_coefficient_log(est_in, D)
> + multiply_log(D, est_if[1])
> + multiply_log((est_in - D), est_if[2]);
>
> increment_log_prob((Pr_D));
> increment_log_prob(-log(W));
>
> I tried to google for binomial_coefficient_log

The index in the manual includes all functions and their
entries have definitions.

> but had limited search results on forum/mailing list. I am going to ask the person who did this previously, but here is what I think it is which may not be correct,
>
> Pr_D is the probability of death, and will take the value est_in which is estimated from posterior distribution of pi from the same model.
>
> binomial_coefficient_log is used in view of unknown est_in and also non-informative/weak prior of est_if (risk of death), defined as having a beta distribution (0.5, 0.5) which is rather random?

It's probably someone trying to write out all the of
terms in a density like the binomial.

Is "rather random" a meta-comment?

> I'm sorry if this doesn't make sense.
>
> Could I have used the number of death, D ~ binomial(est_in, est_if)?

binomial requires a count as first argument and probaiblity in (0, 1)
as the second one. It's just the usual binomial.

>
> theta gets a lower bound 0.0000015625 based on the knowledge that 20 cases were observed in the population of size 12.8e6.

You should just remove that lower bound.

> theta0 takes a prior beta(0.5, 1129.5) based on the knowledge that the 0 out of 1129 samples was found positive

This isn't the same as the 20 out of 1e7?

But if you started from a uniform beta(1, 1), the posterior
after 0 out of 1129 samples would be beta(1, 1130).

If that study's the same as this one, you can either combine
the data, or start with that posterior if nothing else is going on.

> in a previous study and alpha for beta distribution has to be >0.

beta(1, 1) is uniform; beta(0.5, 0.5) is very U-shaped, with density
going to infinity as variate approaches 0 or 1.

> It is also set to bound between 0 and 1 because it is a parameter for baseline cumulative incidence, which is a proportion.

Yes, you need that constraint.

> if prior for pi is known to be beta(0.5, 193.5), can I incorporate this into the model when at the same time, pi is defined as theta0 + theta*x1?

You should think about priors on theta0 and theta rather
than a function of them.

- Bob

[huiying]

Hi Bob,

It comes from someone who was working on this project previously. 

I said "rather random" as I do not see from where the beta(0.5, 0.5) for priors est_if and est_is come from.

I think I understand D ~ binomial(est_in, est_if). I have the information on D, but can we use it when est_in is going to be estimated from the model and we have very weak prior est_if?

prior beta(0.5, 1129.5) for theta0 is different from 20 out of 12.8e6.

theta0 is defined as the baseline cumulative incidence, which is the proportion of positive before an outbreak. and it was observed that 0 out of 1129 samples was tested positive before the outbreak. This is a finding from a separate study that took place at a different city.

My current prior for theta0 is beta(0.5, 1129.5). So do you mean it is better to set it at beta(1, 1130) instead, as it takes the same shape as beta(1,1)?

theta is defined as the actual cumulative incidence that we aim to find out, and it was observed that 20 out of 12.8e6 samples were tested positive. Thus, theta is set to have the value of at least 20/12.8e6.

My current prior for theta is beta(0.5, 0.5). 

Both theta and theta0 have their own prior.

The only one that has no prior incorporated yet is the parameter pi. We believe that pi takes the value theta0+theta*x1, and prior for pi was also calculated from 0/1129, using weighted age distribution to obtain age-specific prior. that gives beta(0.5, 193.5).

Do you mean I should keep beta(0.5, 193.5) for pi instead and not the function for pi?

Huiying

[Bob Carpenter]

I went back and looked at that model and there are all kinds
of problems with it, I'm afraid.

This:

for (i in 1:W)

pi_0[i] <- theta0 + theta * x1[i];

for (i in 1:W)

pi[i] <- fmin(pi_0[i], 1);

It doesn't make sense to truncate like this. Instead, you
want to do a logistic regression and do this:

pi[i] <- inv_logit(pi_0[i]);

But rather than doing that, you can just use binomial_logit.
But little steps first.

And you don't want to put a prior on pi_0; you want to
express the priors on theta0 and theta. If x is centered,
you can put the informative prior on theta0 (on the logit
scale).

I don't like the variables names "est_is" for estimated severity
risk; the variable is just the severity risk and you estimate
it conditional on data.

Then I'm not sure why the code has the binomial spelled out
the way it is.

I'd suggest trying to write the model down mathematically first
as a GLM with appropirate priors, then simulate data and fit it,
then apply it to real data.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <modelcode.stan><Rstan_test.R>

[huiying]

Hi Bob,

Thanks for taking time to go through with me my command line by line. 

I was trying to digest them and agree that I should rewrite them.

This is what I've got so far.

data {

  int W; // Number of weeks samples were collected

  

  int n[W]; // Number of serum samples collected at week i (n will be a vector of length W)

  

  int y[W]; // Number of positive serum samples at week i

  int S; // Total number of severe infections

  int D; // Total number of deaths

  real x1[W]; // Number of expected scaled seroprevalence

}

parameters {

  real<lower=0, upper=1> theta; // Cumulative incidence of infection in each age group

  real<lower=0, upper=1> theta0; // Baseline cumulative incidence in all age group

  real<lower=0, upper=1> est_is; // Estimated infection-severity risk 

  real<lower=0, upper=1> est_if; // Estimated infection-fatality risk

}

transformed parameters {

  real pi_0[W];

  

  real<lower=0, upper=1> pi[W]; // Probability of being infected

  real est_in; // Estimated number of infections

  

  est_in <- 12800000*theta; 

  

  for (i in 1:W)

  pi_0[i] <- theta0 + theta * x1[i]; 

  

  

  for (i in 1:W)

  pi[i] <- logit(pi_0[i]); //logit to transform pi_0(real numbers) to a value in the interval[0, 1]

}

model {    

  theta ~ beta(0.5, 0.5); // Prior information on theta

  theta0 ~ beta(0.5, 193.5); // Prior information on age-specific theta0

  y ~ binomial(n, pi); // Binomial likelihood of a positive serum sample; 

  S ~ binomial(est_in, est_is); // Binomial likelihood of a severe infection 

  D ~ binomial(est_in, est_if); // Binomial likelihood of a fatal infection

}

It gives me errors pointing at the last few lines that I'm not sure how to fix.

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

No matches for: 

  int ~ binomial(real, real)

Available argument signatures for binomial:

  int ~ binomial(int, real)

  int ~ binomial(int, real[])

  int ~ binomial(int, vector)

  int ~ binomial(int, row vector)

  int ~ binomial(int[], real)

  int ~ binomial(int[], real[])

  int ~ binomial(int[], vector)

  int ~ binomial(int[], row vector)

  int[] ~ binomial(int, real)

  int[] ~ binomial(int, real[])

  int[] ~ binomial(int, vector)

  int[] ~ binomial(int, row vector)

  int[] ~ binomial(int[], real)

  int[] ~ binomial(int[], real[])

  int[] ~ binomial(int[], vector)

  int[] ~ binomial(int[], row vector)

ERROR at line 50

 48:      y ~ binomial(n, pi); // Binomial likelihood of a positive serum sample; 

 49:    

 50:      S ~ binomial(est_in, est_is); // Binomial likelihood of a severe infection 

                                       ^

 51:

The issue is est_in is a real number. I tried est_in <- round(12800000*theta) but it is still being recognised as real number.

I could potentially use floor but I read on the manual that says smth like the a less efficient sampler will be resulted etc. If floor can work, I believe I could also use round which is also stated in the manual.

I wrote S ~ binomial(est_in, est_is) as the binomial likelihood for S is L(est_is|est_in, S) = est_is^S x (1 - est_is)^(est_in - S), so does D.

However, the first argument for binomial (int ~ binomial(real, real)) has to be an integer.

Any advice how I can go about this?

Huiying

 

[Bob Carpenter]

The binomial has to be an integer. You have to know it
ahead of time. Stan doesn't allow integer unknowns (and there's
no good way to convert a real to a rounded integer --- it breaks
all the derivatives in the systems).

- Bob

[huiying]

Theta is estimated by the model. 

So, instead of using the mean value of theta (cumulative incidence) to estimate total number of infection (est_in) using the following equation:

est_in =  total population size * theta

Can I treat est_in as a stochastic node to express its uncertainties? Something like est_in ~ binomial(12800000, theta), which I could not do because of the above reason mentioned e.g unknown est_in.

Is there other way I can consider or is it good enough if I have confidence interval expressed along est_in calculated by total population size * theta?

Huiying

[huiying]

I've tried to run two different models. The first is to estimate age-specific cumulative incidence, theta.

data {

  int N; // Number of weeks samples were collected

  int n[N]; // Number of serum samples collected at week i (n will be a vector of length W)

  int y[N]; // Number of positive serum samples at week i

  int S; // Total number of severe infections

  int D; // Total number of deaths

  real x1[N]; // Number of expected scaled seroprevalence

  int M; // population size for each age group

  real ub; //Upper bound for age-specific prior theta0

}

parameters {

  real<lower=0, upper=1> theta; // Cumulative incidence of infection in each age group

  real<lower=0, upper=1> theta0; // Baseline cumulative incidence in all age group

}

transformed parameters {

  real<lower=0, upper=1> p[N]; // Probability of being infected

  

  real est_inf;// Estimated number of infections

  

  for (i in 1:N)

  p[i] <- inv_logit(theta0 + theta * x1[i]); // binomial_logit to transform p to values in the interval [0,1]

  

  est_inf <- theta*M; // M is to be replaced by general population size of each age group

}

  

model {

  theta ~ beta(0.5, 0.5); // Noninformative prior on theta

  theta0 ~ beta(0.5, ub); // Prior information on age-specific theta_0

  for(i in 1:N)

  y[i] ~ binomial(n[i], p[i]); //Binomial likelihood of a positive serum sample

}

The output from the first model will be used to calculate population cumulative incidence and the estimated total number of cases, est_inf.

Est_inf will then be used as an input for the 2nd model, as below, to estimate infection-severity risk (est_isr) and infection-fatality risk (est_ifr). 

data {

  int est_inf; // Estimated number of infected cases

  int S; // Number of severe cases

  int D; // Number of deaths

}

parameters {

  real<lower=0, upper=1> est_isr; // infection-severity risk

  real<lower=0, upper=1> est_ifr; // infection-fatality risk

}

model {

  est_isr ~ beta(0.5, 0.5); // prior for est_isr

  est_ifr ~ beta(0.5, 0.5); // prior for est_ifr

  S ~ binomial(est_inf, est_isr); // binomial likelihood of severe cases

  D ~ binomial(est_inf, est_ifr); // binomial likelihood of fatal cases

}

Doing two models do not seem to be very ideal. However, with my limited knowledge, I guess this is one way I can get around with the limitations I obtained when attempting to everything in one model. 

Please let me know what you think of this approach. I'm sorry if this seems stupid.

Thanks !

Huiying

[Bob Carpenter]

You almost always want to put everything in a single joint
model. There's no good way to chain models together without
doing that.

- Bob

[huiying]

Hi Bob !

I decided to take your and try to put this whole block in the generated quantities section, as these parameters to be estimated (est_isr, est_ifr and est_inf) will take the value of other parameters returned by rstan.

generated quantities {

  real<lower=0, upper=1> est_isr; // Estimated infection-severity risk

  real<lower=0, upper=1> est_ifr; // Estimated infection-fatality risk

  real est_inf;// Estimated number of infections

  int int_inf;

  int S; // Total number of severe infections

  int D; // Total number of deaths

 

  est_isr <- beta_rng(0.5, 0.5); // prior for estimated isr

  

  est_ifr <- beta_rng(0.5, 0.5); // prior for estimated ifr

  

  est_inf <- theta*M; 

  

  int_inf <- 1;

  while ((int_inf+1) < floor(est_inf))

  int_inf <- int_inf+1; // force est_inf from real to integer

  

  S <- binomial_rng(int_inf, est_isr); // Binomial likelihood of a severe case

  

  D <- binomial_rng(int_inf, est_ifr); // Binomial likelihood of a death case

}

However, I'm not sure why in the final output, S and D are also returned and seem to be estimated by the model as well. I have data for S and D specified in my data list that is supposed to be taken by. 

When I put the whole block in model section, it runs much slower with warning message concerning max_treedepth and the results also differ (which could be related to the max_treedepth?)???

Would appreciate if you or anyone could enlighten me.

Huiying

[Bob Carpenter]

It's very hard to imagine what the rest of your model is from
the generated quantities block.

You can't declare S and D in multiple blocks, so your statement
about the model can't be right (or maybe I misunderstood it).

S and D are declared and defined in the generated quantities block, so they
show up in the output.

- Bob

> On Apr 26, 2016, at 6:32 AM, huiying <hyc.t...@gmail.com> wrote:
>
> Hi Bob !
>
> I decided to take your and try to put this whole block in the generated quantities section, as these parameters to be estimated (est_isr, est_ifr and est_inf) will take the value of other parameters returned by rstan.
>
> generated quantities {
> real<lower=0, upper=1> est_isr; // Estimated infection-severity risk
> real<lower=0, upper=1> est_ifr; // Estimated infection-fatality risk
> real est_inf;// Estimated number of infections
> int int_inf;
> int S; // Total number of severe infections
> int D; // Total number of deaths
>
> est_isr <- beta_rng(0.5, 0.5); // prior for estimated isr
>
> est_ifr <- beta_rng(0.5, 0.5); // prior for estimated ifr
>
> est_inf <- theta*M;
>
> int_inf <- 1;
> while ((int_inf+1) < floor(est_inf))
> int_inf <- int_inf+1; // force est_inf from real to integer
>
> S <- binomial_rng(int_inf, est_isr); // Binomial likelihood of a severe case
>
> D <- binomial_rng(int_inf, est_ifr); // Binomial likelihood of a death case
> }
>
> However, I'm not sure why in the final output, S and D are also returned and seem to be estimated by the model as well. I have data for S and D specified in my data list that is supposed to be taken by.
>

> When I put the whole block in model section, it runs much slower.

>
> Would appreciate if you or anyone could enlighten me.
>
> Huiying
>

[huiying]

Apologise I haven't been very clear.

I've attached the two different model codes, one with S and D in the generated quantities block, the other in model block. Hope that clarifies.

The results given by the two are quite different, in terms of theta, est_isr and est_ifr.

I'm not sure which is the right way to do, I have the feeling its test1?

Also, I have a separate question. I realised rstan does not allow variable to replace value for alpha in the prior theta0 (model block). 

What I mean is I can write theta0 ~ beta(1.5, beta) with beta defined in datalist. but with theta0 ~ beta(alpha, beta), a warning message says alpha variable does not exist, even though I have it defined in datalist?

Huiying

[Bob Carpenter]

I can't follow what you're trying to do, but it's unlikely test1.stan
is right.

All this does in generated quantities:

real<lower=0, upper=1> est_ifr;

est_ifr <- beta_rng(0.5, 0.5);

is to assign est_ifr a random beta(0.5, 0.5) draw. So the posterior
mean should always be 0.5.

Almost any conversion you do from real to integer in Stan outside
a piecewise approximation of a function with smooth knots is
going to be wrong.

- Bob

> On Apr 26, 2016, at 9:24 PM, huiying <hyc.t...@gmail.com> wrote:
>
> Apologise I haven't been very clear.
>
> I've attached the two different model codes, one with S and D in the generated quantities block, the other in model block. Hope that clarifies.
>
> The results given by the two are quite different, in terms of theta, est_isr and est_ifr.
>
> I'm not sure which is the right way to do, I have the feeling its test1?
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <test.R><test1.stan><test2.stan>

[huiying]

Hi Bob,

Thanks for explaining!

What I'm trying to do can be divided into 3 parts:

Part 1: to estimate posterior distribution of theta (cumulative incidence). 

----Part 1 is pretty straight forward, and is expressed from data and parameters block to model block, as in test1.stan.

Part 2: Then, using posterior theta to estimate est_inf (number of infection), using this equation: est_inf = theta x population size.

Part 3: Then, with est_inf, estimate est_isr (infection-severity risk) and est_ifr (infection-fatality risk), which are expressed as binomial likelihoods in the following equations:

S ~ binomial(est_inf, est_isr) and D ~ binomial(est_inf, est_ifr), where S and D are known.

Part 2 and 3 are written in generated quantities block because they take posterior estimate of theta (from part 1) as input (we say they are conditional on theta?).

It is mainly the part 3 that I'm struggling with, because,

1. est_inf has to be integer. From what I gathered from the mailing list discussion, a while loop can convert it from real to integer, if I really need to? I do understand the limitations associated with converting real to integer.

2. I can't write the usual format D ~ binomial(est_inf, est_ifr), as ~ is only allowed in the model block. 

3. it seems like est_ifr <- beta_rng(0.5, 0.5) is not the right way to specify a prior beta distribution for est_ifr? If so, how should I write?

Is my understanding correct? If you have any solutions to the challenges I'm facing, please kindly advise.

Huiying