Solving a double integral associated with Multivariate Normal desity

Question

I am trying to solve a double integral associated with the Multivariate Normal density with known mean vector and covariance matrix:

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))

However, it keeps giving me the error:

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

Is there any obvious error with my code?

Answer 1

HubertL has pointed out that the first argument should be a function, not a function-call with arguments. It is assumed that the function will accept an "x" argument, a single length 2 vector, so the NormalPDF function needs to be modified in its arguments and in its call to the helper function. Another error was is how the limits are set up.

Consider this:

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

This integrates that density over the square with "lower left hand corner" at (1,1) and "upper right corner" at (3,3). The invocation in the question would have always returned 0, since the domain was a single point. Would need to be extracted from the list with P$integral if you were going to do anything "numerical" with it. Seems reasonable that the result was less than 0.25 since we were only evaluating in the quarter-plane from the maximum at (3,3).