我在 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
有没有办法在不手动选择间隔的情况下整合这个密度?

最佳答案

不确定这是否有帮助——可能对 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

关于r - 积分一个非常尖的函数,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/20665689/

10-12 13:58
查看更多