求解与多元正态密度相关的二重积分

Question

我正在尝试用已知的均值向量和协方差矩阵求解与多元正态密度相关的二重积分：

library(cubature)

mu1 <- matrix(c(3,3), nrow=2)
sigma1 <- rbind(c(4,-1), c(-1,6))

quadratic <- function(a,b) {
  X <- matrix(c(a,b),nrow=2)
  Q <- (-1/2)*t(X-mu1)%*%solve(sigma1)%*%(X-mu1)
}

NormalPDF <- function(x1,x2) {
  f <- (1/(2*pi))*(1/sqrt(det(sigma1)))*exp(quadratic(x1,x2))
}

# Solving for P(1 < X1 < 3, 1 < X2 < 3)
P <- adaptIntegrate(NormalPDF(x1,x2), c(1,3), c(1,3))

但是，它一直给我错误：

Error in matrix(c(a, b), nrow = 2) : object 'x1' not found

我的代码有没有明显的错误？

Answer 1

HubertL 指出第一个参数应该是函数，而不是带参数的函数调用。假定该函数将接受一个 "x" 参数，一个长度为 2 的向量，因此 NormalPDF 函数需要在其参数和对辅助函数的调用中进行修改。另一个错误是如何设置限制。

考虑一下：

library(cubature)

mu1 <- matrix(c(3,3), nrow=2)
sigma1 <- rbind(c(4,-1), c(-1,6))

quadratic <- function(a,b) {
  X <- matrix(c(a,b),nrow=2)
  Q <- (-1/2)*t(X-mu1)%*%solve(sigma1)%*%(X-mu1)
}

NormalPDF <- function(x) {
  f <- (1/(2*pi))*(1/sqrt(det(sigma1)))*exp(quadratic(x[1],x[2]))
}
# Solving for P(1 < X1 < 3, 1 < X2 < 3)
P <- adaptIntegrate( NormalPDF, lowerLimit= c(1,1),  upperLimit=c(3,3))
P
#==============
$integral
[1] 0.09737084

$error
[1] 1.131395e-08

$functionEvaluations
[1] 17

$returnCode
[1] 0

这将在 (1,1) 处 "lower left hand corner" 和 (3,3) 处 "upper right corner" 的正方形上对密度进行积分。问题中的调用总是返回 0，因为域是一个点。如果您要对其进行任何操作，则需要使用 P$integral 从列表中提取 "numerical"。结果小于 0.25 似乎是合理的，因为我们仅从 (3,3) 处的最大值在四分之一平面中进行评估。

求解与多元正态密度相关的二重积分

Solving a double integral associated with Multivariate Normal desity

r

numerical-integration