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

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) 处的“左下角”和 (3,3) 处的“右上角”进行积分。问题中的调用始终返回 0，因为域是单个点。如果您要使用它执行任何“数字”操作，则需要使用 P$integral 从列表中提取。结果小于 0.25 似乎是合理的，因为我们仅在四分之一平面中从 (3,3) 处的最大值进行评估。