R语言计算多元积分_r语言偏积分-CSDN博客

本文链接：https://blog.csdn.net/weixin_46111814/article/details/129809117

0引言

在R语言中我们经常使用integrate函数计算一元积分，例如：《R语言【integrate】函数》¹. 但是在很多情形下我们需要计算多元积分（如归一化多元概率密度函数，求不规则物体的体积等），本文介绍R语言里的多元积分函数（或多重积分）。

一、包的载入与认识

包cubature主要围绕 Steven G. Johnson 的 cubature C 库的 R 包装器，用于超立方体上的自适应多元积分，以及 Thomas Hahn 的 Cuba C 库，用于确定性和蒙特卡罗积分。提供了用于容积和古巴例程的标量和矢量接口;强烈建议使用矢量接口，如封装小插图所示。

>  library(cubature) 
> ls("package:cubature")
[1] "adaptIntegrate" "cubintegrate"   "cuhre"          "default_args"  
[5] "divonne"        "hcubature"      "pcubature"      "suave"         
[9] "vegas"

上面给出了包cubature中的主要函数。R语言中可以使用adaptIntegrate函数求解多元积分的问题。该函数在超立方体上执行（可能的）向量值积分的自适应多维积分（立方体）。该函数包括一个向量接口，其中积分可以在单个调用中在数百个点上进行评估。下面看一下这个函数进行积分时的具体例子。，

二、使用实例

2.1 参数解析

先看一下函数主要参数:

> adaptIntegrate
function (f, lowerLimit, upperLimit, ..., tol = 1e-05, fDim = 1, 
    maxEval = 0, absError = 0, doChecking = FALSE, vectorInterface = FALSE, 
    norm = c("INDIVIDUAL", "PAIRED", "L2", 
        "L1", "LINF"))

每个参数的含义是：

f：要被积分的函数(被积函数)
lowerLimit：积分的下限，一个超立方体的向量
upperLimit：积分的上限，一个超立方体的向量
.．.所有其他参数传递给函数f
tol:最大公差，默认为1e-5。
fDim:被积函数的维数默认为1，与超立方体的维数无关
maxEval:所需函数求值的最大数量，默认0表示没有限制。注意，实际执行的函数求值的数量只能大致保证不超过这个数字。
absError:允许的最大绝对误差
doChecking:从版本2.0开始，这个标志将被忽略，并将在后续版本中删除
vectorInterface:指示是否使用vector接口的标志，默认为FALSE。详见下文
norm:对于向量值被积函数，范数指定用于测量误差和确定收敛性的范数。见下文。

2.2 运行实例1：多元正态分布的积分

> dmvnorm <- function (x, mean, sigma, log = FALSE) {
+     if (is.vector(x)) {
+         x <- matrix(x, ncol = length(x))
+     }
+     if (missing(mean)) {
+         mean <- rep(0, length = ncol(x))
+     }
+     if (missing(sigma)) {
+         sigma <- diag(ncol(x))
+     }
+     if (NCOL(x) != NCOL(sigma)) {
+         stop("x and sigma have non-conforming size")
+     }
+     if (!isSymmetric(sigma, tol = sqrt(.Machine$double.eps),
+         check.attributes = FALSE)) {
+         stop("sigma must be a symmetric matrix")
+     }
+     if (length(mean) != NROW(sigma)) {
+         stop("mean and sigma have non-conforming size")
+     }
+     distval <- mahalanobis(x, center = mean, cov = sigma)
+     logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
+     logretval <- -(ncol(x) * log(2 * pi) + logdet + distval)/2
+     if (log)
+         return(logretval)
+     exp(logretval)
+ }
> 
> m <- 3
> sigma <- diag(3)
> sigma[2,1] <- sigma[1, 2] <- 3/5 ; sigma[3,1] <- sigma[1, 3] <- 1/3
> sigma[3,2] <- sigma[2, 3] <- 11/15
> hcubature(dmvnorm, lower=rep(-0.5, m), upper=c(1,4,2),
+                         mean=rep(0, m), sigma=sigma, log=FALSE,
+                maxEval=10000)
$integral
[1] 0.3341125

$error
[1] 4.185435e-06

$functionEvaluations
[1] 10065

$returnCode
[1] 0

2.3 运行实例2：二元函数积分

计算如下函数的积分
$sin(4x_1+1)*cos(4x_2)*x_1*(x_1*(x_1*x_1)^2-x_2*(x_2*x_2-x_1)+2)$

> I.2d <- function(x) {
+   x1 = x[1]
+   x2 = x[2]
+   sin(4*x1+1) * cos(4*x2) * x1 * (x1*(x1*x1)^2 - x2*(x2*x2 - x1) +2)
+ }
> 
> hcubature(I.2d, rep(-1, 2), rep(1, 2), maxEval=10000)
$integral
[1] -0.01797993
$error
[1] 7.845607e-07
$functionEvaluations
[1] 10013
$returnCode
[1] 0
> pcubature(I.2d, rep(-1, 2), rep(1, 2), maxEval=10000)
$integral
[1] -0.01797993
$error
[1] 2.092601e-10
$functionEvaluations
[1] 1089
$returnCode
[1] 0

三、总结与展望

本文主要介绍求解多元定积分或者多重定积分的函数，并给出了两个实例。但是在实际使用中，当多元函数的维度比较多时，计算速度会变慢，但是可以通过调节收敛精度tol加快积分的计算和提升收敛速度。
希望本文的内容可以帮助到大家，您的批评是我前进的动力，欢迎评论区留言讨论本文的相关内容。注：如果想寻求进一步的合作，欢迎移步统计学小王子咨询中心处申请。