Integrate_1d errors in stan

Algorithms
,
1

Hi all,

I am trying to run following stan code with CmdStanR to estimate parameters.
These are the several equations in my model

D_{1,t_{0}}=\beta(t_{0})\int_{0}^{t_{0}}w(x)Q(t_{0}-x)dx
D_{k,t_{0}}=\beta(t_{k-1})\Bigg\{\sum_{i=0}^{k-2}[1-\beta(t_{i})]...[1-\beta(t_{k-2})]\int_{t_{i-1}}^{t_{i}}w(x)Q(t_{k-1}-x)dx +\int_{t_{k-2}}^{t_{k-1}}w(x)Q(t_{k-1}-x)dx \Bigg\}
I_{k,t_{0}}(t)=\sum_{i=0}^{k-1}[1-\beta(t_{i})]...[1-\beta(t_{k-1})]\int_{t_{i-1}}^{t_{i}}w(x)q(t-x)dz +\int_{t_{K-1}}^{t}w(x)q(t-x)dz, \hspace{0.5cm} \forall t \in (t_{k-1},t_{k})

This is a sample data table used in my coding,
Data Table.csv (217 Bytes)

and stan file is

model.stan (5.0 KB)

I am getting errors because of the integrations in my model.

Chain 1 Exception: integrate: error estimate of integral 3.53534 exceeds the given relative tolerance times norm of integral

Could you please suggest a way to solve this issue, I have already simplified a bit according to previous suggestions by @bbbales2 . I want to use following uniform distributions for my initial values .

init_func <- function(){
list(b0=runif(1,0,4),mu=runif(1,4,4.5),sigma2=runif(1,0.01,0.05),lambda=runif(1,0.01,0.5),alpha=runif(1,1.5,4))
}

Here is the code

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

MT<- as.matrix(DataTable)
head(MT)
n=nrow(MT)
m=ncol(MT)
tstar = mean(MT[,1])
K=2

model_data<-list(MT=MT,n=n,m=m,K=K, tstar=tstar)
mod<-cmdstan_model("C:/Users/nanayaer/Documents/.cmdstanr/cmdstan-2.24.1/model.stan")

fit <- mod$sample(
  data = model_data,
  init = init_func,
  chains = 1,
  refresh = 500,
  iter_warmup = 200,
  max_treedepth=15 , 
  adapt_delta = 0.99
)

THANK YOU !

2

These integral problems can be tricky, but the thing you need to do is:

  1. Figure out the input to integrate_1d that causes the error

  2. Verify that the function you are integrating actually is integrable by Stan (it needs to be continuous)

  3. Make sure that you don’t have any bugs in your code

The easiest way to figure out #1 is to add print debugging statements before you integrate_1d function call.

So like with this:

integrate_1d(WQ_integral,t[1],t[2],{mu, sigma2,lambda, alpha, t[2]}, {0.0}, x_i, 1e-8))

Maybe add print(t[1], t[2], mu, sigma2, lambda, alpha); before it, record the values of these variables right before the error throws, and then make plots of the function.

It would be easiest if you compressed your model to just one integral – that way you know exactly the integral that is causing the problems and there isn’t any confusion.

Once you have identified what the arguments are that are causing the problem, you’ll want to compute a reference integral somehow – if Stan is going to compute this it should be computable in R or something.

You also want to make plots of your function and make sure that there isn’t anything crazy going on (that it looks like something the numerical quadrature would work on).

One way of doing this is using the expose_stan_functions function in rstan to expose the function that you have written to R (https://mc-stan.org/rstan/reference/expose_stan_functions.html). This way you can directly make plots using the Stan code to see if there are things that might cause problems.

Once you can compute your integral in R and it looks like your integrand is formulated in Stan correctly, you’ll want to try to integrate your function in Stan. Again you can do this by writing the function in Stan and exposing it with expose_stan_functions.

At some point in this process you should find your problem! Hopefully it’s easy to fix :D.

2 Likes
3

Hi @bbbales2,

Thank you again for the instructions provided here.
Yes, this is a very tricky integration. So, I moved to use numerical integrations in my stan code.

One question, do you how to use these codes in CmdStanR.

MCMCsummary(fit)
MCMCtrace(fit)
MCMCplot(fit)
MCMCchains(fit)

THANK YOU!

4

Summary for cmdstanr is just:

fit$summary()

Or just fit in an R terminal should give you a summary printout.

You can do traceplots with with bayesplot (https://mc-stan.org/bayesplot/) using mcmc_trace(fit$draws()) (https://cran.r-project.org/web/packages/bayesplot/vignettes/plotting-mcmc-draws.html).

I’m not sure autocorrelations (I’d just go by the effective sample size estimates in the summary output instead of worrying about those), but you can do density plots in Bayesplot too.

1 Like
5

Hi Ben, @bbbales2

Great, Thank you very much.
I did it.

Best. :)

6

Hi Ben @bbbales2,

When the $sample() method of a CmdStanModel object runs, how we can change the number of iterations ? Is default 2000?

So then, please explain how to use iter_warmup and iter_sampling with custom number of iterations. Suppose I want to run 10000 iterations.

THANK YOU.

8

If you want 10000 posterior samples and are running four chains, run with iter_sampling = 2500.

iter_sampling * num_chains is how many posterior samples you get.

9

Hi Ben, @bbbales2

Wow… thank you for the ideas.
I used Trapezoidal rule for above integrations and was able to get results similar to what I wanted.

Not with Integrate_1d. I feel that it is so difficult and tricky task. But, later hope to come back to it.

THANK YOU.

😊

2 Likes
10

If you end up having to do something like this, you can try to monitor your error by doing the same integration in generated quantities with a higher resolution, maybe like 5x smaller trapezoids? Then compare your solution computed in the model itself. Generated quantities runs much less than the model itself so it’s possible to do this.

If this error gets large, you know you have a problem with the integration.

2 Likes