Question

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?

Was it helpful?

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

OTHER TIPS

integrate:

Like all numerical integration routines, these evaluate the function on a finite set of points. If the function is approximately constant (in particular, zero) over nearly all its range it is possible that the result and error estimate may be seriously wrong.

So, the answer is no.

You really need to know something about the function to compute the integral correctly - for any automated algorithm which detects support there is a function for which it fails.

PS (7 years later). For any deterministic algorithm, and any error, there is a function, such that this algorithm will make this error on it.

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top