Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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?

like image 517
FairyOnIce Avatar asked Apr 22 '26 06:04

FairyOnIce


2 Answers

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
like image 77
Adrian Avatar answered Apr 24 '26 20:04

Adrian


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.

like image 20
sds Avatar answered Apr 24 '26 19:04

sds



Donate For Us

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!