我在 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/