[SOLVED] integrate a very peaked function

integrate a very peaked function

I am using integrate function in R to integrate a very peaked function. Say that function is a log-normal density:

 xs <- seq(0,3,0.00001)
 fun <- function(xs) dlnorm(xs, meanlog=-1.057822,sdlog=0.001861871)
 plot(xs,fun(xs),type="l")

From the plot, I know that the peak is at around 0.3-0.4.

If I integrate this density function over its support (with increased abs.tol and increased subdivisions) the integrate() gives me zero, which should not be true.

integrate(fun,lower=0,upper=Inf,subdivisions=10000000,abs.tol=1e-100) 
0 with absolute error < 0

However, if I restrict the interval to 0.3 - 0.4, it gives me the correct answer.

integrate(fun,lower=0.3,upper=0.4,subdivisions=10000000,abs.tol=1e-100) 
1 with absolute error < 1.7e-05

Is there a way to integrate this density without manually choosing the interval?

Solution

Not sure whether this is helpful -- might be too specific to dlnorm, but you can partition [0, Inf[, especially if you have a good idea of where the peak will end up:

integrate.dlnorm <- function(mu=0, sd=1, width=2) {
    integral.l <- integrate(f=dlnorm, lower=0, upper=exp(mu - width * sd), meanlog=mu, sdlog=sd)$value
    integral.m <- integrate(f=dlnorm, lower=exp(mu - width * sd), upper=exp(mu + width * sd), meanlog=mu, sdlog=sd)$value
    integral.u <- integrate(f=dlnorm, lower=exp(mu + width * sd), upper=Inf, meanlog=mu, sdlog=sd)$value
    return(integral.l + integral.m + integral.u)
}

integrate.dlnorm()  # 1
integrate.dlnorm(-1.05, 10^-3)  # .97
integrate.dlnorm(-1.05, 10^-3, 3)  # .998