ywa136 - 8 months ago 48

R 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

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]))
}
P <- adaptIntegrate(NormalPDF, c(3,3), c(6,6))
> P
$integral
[1] 0.1577581
$error
[1] 1.046855e-06
$functionEvaluations
[1] 119
$returnCode
[1] 0
```

This integrates that density over the square with "lower left hand corner" at (3,3) and "upper right corner" at (6,6). The invocation in the question would have always returned 0, since the domain was a single point. Seems reasonable that the result (whaich would need to be extracted from the list with `P$integral`

if you were going to do anything with it) was less than 0.25 since we were only evaluating in the quarter-plane from the maximum at (3,3).

Source (Stackoverflow)