r - 积分一个非常峰值的函数

Question

我正在使用 R 中的积分函数来积分一个非常峰值的函数。假设该函数是对数正态密度:

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

从图中，我知道峰值在 0.3-0.4 左右。

如果我在其支持上集成此密度函数(增加 abs.tol 和增加的分割)，integrate() 给我零，这不应该是真的。

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

但是，如果我将间隔限制为 0.3 - 0.4，它就会给出正确的答案。

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

有没有办法在不手动选择间隔的情况下积分这个密度？

Answer 1

不确定这是否有帮助 - 可能对于 dlnorm 来说过于具体，但您可以分区 [0, Inf[，特别是如果您很清楚峰值将在哪里结束:

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