rintegrate

Integrate function returning roundoff error after working previously


When using integrate to integrate a lognormal density function from 2000 -> Inf, I am returned with an error. I had used a very similar equation previously with no problems.

I have tried disabling stop on error, and setting rel.tol lower. I am fairly new and unfamilar with r so I apologize if neither of those are expected to have done anything.

> integrand = function(x) {(x-2000)*(1/x)*(1/(.99066*((2*pi)^.5)))*exp(-((log(x)-7.641)^2)/((2*(.99066)^2)))}
> integrate(integrand,lower=2000,upper=Inf)
1854.002 with absolute error < 0.018
#returns value fine

> integrand = function(x) {(x-2000)*(1/x)*(1/(1.6247*((2*pi)^.5)))*exp(-((log(x)-9.0167)^2)/((2*(1.6247)^2)))}
> integrate(integrand,lower=2000,upper=Inf)
Error in integrate(integrand, lower = 2000, upper = Inf) : 
  roundoff error is detected in the extrapolation table
#small change in the mu and sigma in the lognormal density function results in roundoff error

> integrate(integrand,lower=1293,upper=Inf)
29005.08 with absolute error < 2
#integrating on lower bound works fine, but having lower=1294 returns a roundoff error again
> integrate(integrand,lower=1294,upper=Inf)
Error in integrate(integrand, lower = 1294, upper = Inf) : 
  roundoff error is detected in the extrapolation table

I should be getting returned a value, no? I struggle to see how very slightly altering the values would cause the function to no longer integrate.


Solution

  • First of all, I believe you are complicating when you define the integrand by writing down the entire expression, it seems better to use the built-in dlnorm function.

    g <- function(x, deduce, meanlog, sdlog){
      (x - deduce) * dlnorm(x, meanlog = meanlog, sdlog = sdlog)
    }
    
    curve(g(x, deduce = 2000, meanlog = 9.0167, sdlog = 1.6247), 
          from = 1294, to = 1e4)
    

    enter image description here

    As for the integration problem, package cubature generally does a better job when integrate fails. All of the following produce the results, with no errors.

    library(cubature)
    
    cubintegrate(g, lower = 1293, upper = Inf, method = "pcubature", 
                 deduce = 2000, meanlog = 9.0167, sdlog = 1.6247)
    
    cubintegrate(g, lower = 1294, upper = Inf, method = "pcubature", 
                 deduce = 2000, meanlog = 9.0167, sdlog = 1.6247)
    
    cubintegrate(g, lower = 2000, upper = Inf, method = "pcubature", 
                 deduce = 2000, meanlog = 9.0167, sdlog = 1.6247)