集成非常峰值的函数

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

Answer 2

integrate

：

像所有数值积分例程一样，这些例程评估函数 在有限的点集上。如果函数近似恒定 （特别是零）在几乎所有范围内，有可能 结果和误差估计可能是严重错误的。

所以，答案是不。

您确实需要了解有关正确计算积分的函数的一些信息 - 对于任何检测支持的自动算法，都有一个函数会失败。

PS（7年后）。对于任何确定性算法，以及任何错误，都有一个函数，使得该算法会在其上犯这个错误。