与极限变量函数的多重集成

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