当上限为 Inf 时,R 中的积分失败

问题描述 投票:0回答:1

我需要对两个对数正态随机变量之和的对数方差进行数值近似。我想在 R 中执行此操作,这是一个示例:

################
## Setup data ##
################

sdX1 = 0.33
sdX2 = 0.70
muX1 = log(32765) - 0.5 * sdX1^2
muX2 = log(52650) - 0.5 * sdX2^2

####################################
## PDF for sum of 2 lognormal RVs ##
####################################

d2lnorm = function(z, muX1, muX2, sdX1, sdX2){

    #PDFs
    L1 = distr::Lnorm(meanlog = muX1, sdlog = sdX1)
    L2 = distr::Lnorm(meanlog = muX2, sdlog = sdX2)

    #Convlution integral
    L1plusL2 = distr::convpow(L1 + L2, 1)

    #Density function
    f.Z = distr::d(L1plusL2)

    #Evaluate
    return(f.Z(z))

}

############################################
## Expectation for sum of 2 lognormal RVs ##
############################################    

ex2lnorm = function(muX1, muX2, sdX1, sdX2){    

    #E(g(x)) = integral of g(x) * f(x) w.r.t x
    integrate(function(z) log(z) * d2lnorm(z, muX1 = muX1, muX2 = muX2, sdX1 = sdX1, sdX2 = sdX2), lower = 0, upper = +Inf)$value

}

##############
## Run code ##
##############

ex2lnorm(muX1, muX2, sdX1, sdX2)

但是,积分计算结果为 0。如果我将

integrate
的上限更改为较大的数字,它就会起作用。但是,我正在运行大量模拟,并且我无法摆弄上限以使其每次都能正常工作。还有其他方法可以让它持续工作吗?

r numerical-integration
1个回答
1
投票

如果找到不为 0 的最大上限值,您可以做什么:

ex2lnorm = function(muX1, muX2, sdX1, sdX2){    

  f <- function(z) log(z) * d2lnorm(z, muX1 = muX1, muX2 = muX2, sdX1 = sdX1, sdX2 = sdX2)

  upper <- 2^10
  cond.lower <- FALSE
  repeat {
    new <- integrate(f, lower = 0, upper = upper)$value
    if (new > 0) cond.lower <- TRUE
    if (cond.lower && new == 0) break
    upper <- upper * 2
    prev <- new
  }

  prev
}
© www.soinside.com 2019 - 2024. All rights reserved.