I have an issue computing the integral of a function in rstan.
With “integrate_1d”, the integral cannot be computed over an interval, whereby it throws up “error estimate of integral above zero exceeds the given relative tolerance times norm of integral above zero” .
However, if I do the integration with “integrate” in R, then I have no such issue over that interval.
I am looking for
- why the two functions behave differently
- how I can get “integrate_1d” to behave more like “integrate”, OR
- how I can use the function “integrate” in my rstan code.
I paste some code below to reproduce the issue. The density function correspond to a dynamic model of choice (Linear Ballistic Accumulator model), and the integral computes the probability that one value is more than the other at a given time.
I also attach a file with the code for convenience.forum question code.R (3.7 KB)
stancode <- 'functions{
real lbaX_pdf(real X, real t, real A, real v, real s){
//PDF of the LBA model
real b_A_tv_ts;
real b_tv_ts;
real term_1;
real term_2;
real pdf;
b_A_tv_ts = (X - A - t*v)/(t*s);
b_tv_ts = (X - t*v)/(t*s);
term_1 = Phi(b_A_tv_ts);
term_2 = Phi(b_tv_ts);
pdf = (1/A)*(-term_1 + term_2);
return pdf;
}
real lbaX_cdf(real X, real t, real A, real v, real s){
//CDF of the LBA model
real b_A_tv;
real b_tv;
real ts;
real term_1;
real term_2;
real term_3;
real term_4;
real cdf;
b_A_tv = X - A - t*v;
b_tv = X - t*v;
ts = t*s;
term_1 = b_A_tv * Phi(b_A_tv/ts);
term_2 = b_tv * Phi(b_tv/ts);
term_3 = ts * exp(normal_lpdf(b_A_tv/ts|0,1));
term_4 = ts * exp(normal_lpdf(b_tv/ts|0,1));
cdf = (1/A)*(- term_1 + term_2 - term_3 + term_4);
return cdf;
}
// proba that 1 is ranked before 2 in first race
real rank_density(real x, // Function argument
real xc, // Complement of function argument
// on the domain (defined later)
real[] theta, // parameters
real[] x_r, // data (real)
int[] x_i) { // data (integer)
real t = theta[1];
real A = theta[2];
real v1 = theta[3];
real v2 = theta[4];
real s = theta[5];
real v;
v=lbaX_pdf(x, t, A, v1, s)*lbaX_cdf(x, t, A, v2, s);
return v;
}
// cumulative proba of 1 before 2 in first race, knowing neither reached b
real order(real down, real up, real[] theta, data real[] x_r) {
int x_i[0];
real v;
v=integrate_1d(rank_density, down, up, theta, x_r, x_i,1e-8);
// }
return v;
}
}
'
library(rstan)
expose_stan_functions(stanc(model_code = stancode))
t<-0.5
b<-1
A<-0.5
v1<-1
v2<-1
s<-1
# This compute the probability that 1 and 2 are less than b but 1 is more than 2 at time t
order(-10,b,theta = c(t,A,v1,v2,s), x_r = double())
order(-10,0.66,theta = c(t,A,v1,v2,s), x_r = double())
order(-10,0.67,theta = c(t,A,v1,v2,s), x_r = double())
order(-10,0.70,theta = c(t,A,v1,v2,s), x_r = double())
order(-10,0.76,theta = c(t,A,v1,v2,s), x_r = double())
order(-10,0.77,theta = c(t,A,v1,v2,s), x_r = double())
# gives out following error, which occurs for x between 0.67 and 0.76
# Exception: integrate: error estimate of integral above zero 9.59331e-10
# exceeds the given relative tolerance times norm of integral above zero
# (in 'unknown file name' at line 74)
# However, if I now use the function integrate in R, I do not have any problems in that domain:
f_1_2<-function(x,t,A,v1,v2,s){
p<-lbaX_pdf(x, t, A, v1, s)*lbaX_cdf(x, t, A, v2, s)
return(p)}
integrate(Vectorize(f_1_2),-10,0.66,t=t,A=A,v1=v1,v2=v2,s=s)$value
integrate(Vectorize(f_1_2),-10,0.67,t=t,A=A,v1=v1,v2=v2,s=s)$value
integrate(Vectorize(f_1_2),-10,0.70,t=t,A=A,v1=v1,v2=v2,s=s)$value
integrate(Vectorize(f_1_2),-10,0.76,t=t,A=A,v1=v1,v2=v2,s=s)$value
integrate(Vectorize(f_1_2),-10,0.77,t=t,A=A,v1=v1,v2=v2,s=s)$value
# here is the whole graph of the function "order" between -1 and 3:
N<-400
x<-(1:N)/100-1
ft<-rep(NA,N)
for (i in (1:N)) {ft[i]<-integrate(Vectorize(f_1_2),-10,x[i],t=t,A=A,v1=v1,v2=v2,s=s)$value}
plot(x,ft)