Решение двойного интеграла, связанного с многовариантным нормальным числом

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

Хьюбертл указал, что первый аргумент должен быть функцией, а не вызовом функции с аргументами. Предполагается, что функция будет принимать аргумент "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).