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

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) 处的最大值在四分之一平面中进行评估。