在 R 中使用 integrate 函数进行集成

Question

library(pbivnorm)

rho <- 0.5
f1 <- function(x, y) {
  pbivnorm(log(x)-10, log(y)-10, rho)*(exp(-(log(x)-10)^2/2)/(sqrt(2*pi)*x))*(exp(-(log(y)-10)^2/2)/(sqrt(2*pi)*y))
}
integration1 <- round(integrate(function(y) {
  sapply(y, function(y) {
    integrate(function(x) f1(x,y), 0, Inf, rel.tol = 1e-12)$value
  })
}, 0, Inf, rel.tol = 1e-12)$value, 10)

这个积分应该在0.3左右，但是R给出了0。谁能指出问题所在？ R中积分的最佳功能是什么？非常感谢。

Answer 1

包cubature可以解决问题，给出预期的结果。该函数必须重写为单参数函数，并在函数体中设置 x 和 y 的值。

library(cubature)

f2 <- function(X) {
  x <- X[1]
  y <- X[2]
  pbivnorm(log(x)-10, log(y)-10, rho)*(exp(-(log(x)-10)^2/2)/(sqrt(2*pi)*x))*(exp(-(log(y)-10)^2/2)/(sqrt(2*pi)*y))
}

hcubature(f2, c(0, 0), c(Inf, Inf))
#$integral
#[1] 0.2902153
#
#$error
#[1] 2.863613e-06
#
#$functionEvaluations
#[1] 7599
#
#$returnCode
#[1] 0

编辑。

在之后，这里是用 hcubature

计算的积分

f3 <- function(x) {
  pnorm(log(x)-10.2)*(exp(-(log(x)-10)^2/2)/(sqrt(2*pi)*x)) 
} 

hcubature(f3, lowerLimit = 0, upperLimit = Inf, tol = 1e-12)$integral 
#[1] 0.4437685

在 R 中使用 integrate 函数进行集成

Integration in R with integrate function

r

integral