R integrate: returns wrong solution (is using wrong quadrature points?)

Question

I have a function in R which I am trying to integrate, but for some (extreme) values of the function parameters, integrate returns the incorrect solution. I believe the issue may be that integrate selects improper quadrature points for some of these extreme values, but first I will provide demonstrate the issue.

The function I wish to integrate is the following.

integrandFunc_F <- function(x, func_u, func_u_lowerBar, 
  func_u_upperBar, func_mean_v, func_sigma_v, func_sigma_epsilon, 
  func_sigma_y, func_gamma, func_rho) {
#print(x);
p <- 1 - pnorm(func_u_upperBar,x,func_sigma_y);
q <- pnorm(func_u_lowerBar,x,func_sigma_y);
p <- p*(1-func_rho); q <- q*(1-func_rho);
alpha <- ifelse(func_gamma*(p+q) == 0, 0, pmax((func_gamma*p-q)/(func_gamma*(p+q)), 0));
g <- ifelse(x > func_u, dnorm(x,func_mean_v,sqrt(func_sigma_v^2 + func_sigma_epsilon^2))/(1-pnorm(func_u,func_mean_v,sqrt(func_sigma_v^2 + func_sigma_epsilon^2))), 0);
output <- alpha*g;
output
}

When I try to calculate the following, I get the correct solution of 1:

integrate(integrandFunc_F, lower=-Inf, upper=Inf, func_u= 8, func_u_lowerBar= 8, 
  func_u_upperBar= 8, func_mean_v= 30, func_sigma_v= .1, func_sigma_epsilon= 2, 
  func_sigma_y= 1, func_gamma= 1/1.1, func_rho= .05)

However, when I try to calculate the following, I get the incorrect solution of 0:

integrate(integrandFunc_F, lower=-Inf, upper=Inf, func_u= 8, func_u_lowerBar= 8, 
  func_u_upperBar= 8, func_mean_v= 50, func_sigma_v= .1, func_sigma_epsilon= 2, 
  func_sigma_y= 1, func_gamma= 1/1.1, func_rho= .05)

Above I indicated I believe the issue may have to do with the selection of quadrature points. If you uncomment #print(x) in the function above, you can see in the func_mean_v = 30 case, integrate settles on quadrature points that are relatively large/near 30. However, in the func_mean_v=50 case, after a few iterations integrate selects quadrature points that are near 0. Quadrature points near 0 are inappropriate to evaluate this function, which incorporates a normal distribution with mean at func_mean_v.

Any ideas on how to address this issue? Why would integrate iterate to quadrature points near 0 in some cases? Note, the choices of func_mean_v = 30 and func_mean_v = 50 are admittedly extreme parameters for this function, however I need to be able to calculate such cases correctly.

Answer 1

you could shift the integration variable to centre the peak,

wrapper <- function(x, func_mean_v, ...)
   integrandFunc_F(x+func_mean_v, func_mean_v=func_mean_v, ...)


integrate(wrapper, rel.tol = 1e-8, lower=-Inf, upper=Inf, func_u= 8, func_u_lowerBar= 8, 
          func_u_upperBar= 8, func_mean_v= 50, func_sigma_v= .1, func_sigma_epsilon= 2, 
          func_sigma_y= 1, func_gamma= 1/1.1, func_rho= .05)
# 1 with absolute error < 1.3e-09