integrate a very peaked function

integrate, r

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

Problem

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?

Original source