r - 极限内变量函数的多重积分

Question

我需要对以下内容进行数值积分:

我尝试使用 cubature 和 pracma 但它们似乎不支持功能集成限制。我发现尝试通过以下方式使用cubature:

library(cubature)

integrand <- function(arg) {
   x <- arg[1]
   y <- arg[2]
   z <- arg[3]
   w <- arg[4]
   v<- arg[5]
   ff <- dnorm(x, 10,2)*dnorm(y, 10,2)*dnorm(z, 10,2)*dnorm(w, 10,2)* dnorm(v, 10,2)* (x+y+z+w+v<=52)
   return(ff)
}

R <- cuhre(f = integrand, 
           lowerLimit=c(0,0,0,0,0), 
           upperLimit=c(20,20,20,20,20), 
           relTol = 1e-5, absTol= 1e-5)

但作者不保证这样做是正确的。

有没有办法在 R 中对具有函数极限的多个积分进行数值积分？

Answer 1

积分域是按因子 42 缩放的规范单纯形。要计算单纯形上的积分，请使用 SimplicialCubature 包:

integrand <- function(arg) {
  x <- arg[1]
  y <- arg[2]
  z <- arg[3]
  w <- arg[4]
  v <- arg[5]
  dnorm(x, 10, 2) * 
    dnorm(y, 10, 2) * 
    dnorm(z, 10, 2) * 
    dnorm(w, 10, 2) * 
    dnorm(v, 10, 2)
}

library(SimplicialCubature)
Simplex <- 42 * CanonicalSimplex(5)

这是要运行的命令:

adaptIntegrateSimplex(integrand, S = Simplex)
# $integral
# [1] 0.03252553
# 
# $estAbsError
# [1] 0.3248119
# 
# $functionEvaluations
# [1] 9792
# 
# $returnCode
# [1] 1
# 
# $message
# [1] "error: maxEvals exceeded - too many function evaluations"

算法已达到最大函数求值次数，估计绝对误差为0.3248119，而积分估计值为0.03252553。这是一个很大的错误。

我们可以增加允许的函数评估的最大数量。以1e6为例，计算有点慢，我们得到:

adaptIntegrateSimplex(integrand, S = Simplex, maxEvals = 1e6)
# $integral
# [1] 0.03682535
# 
# $estAbsError
# [1] 0.001004083
# 
# $functionEvaluations
# [1] 999811
# 
# $returnCode
# [1] 1
# 
# $message
# [1] "error: maxEvals exceeded - too many function evaluations"

估计误差已降至 0.001004083，相当好。

请注意，我们可以通过模拟来近似该积分，因为该积分是多元正态分布下单纯形的测度:

library(mvtnorm)
Sigma <- 2^2 * diag(5)
Mean <- rep(10, 5)
set.seed(666)
sims <- rmvnorm(1e6, mean = Mean, sigma = Sigma)

f <- function(X){ # test whether 0 < x < 42, 0 < x + y < 42, 0 < x + y + z < 42, ...
  all(X > 0 & cumsum(X) < 42)
}

mean(apply(sims, 1, f))
# 0.037083