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

Я пытаюсь решить двойной интеграл, связанный с многомерной нормальной плотностью с известным средним вектором и ковариационной матрицей:

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

Есть ли очевидная ошибка в моем коде?


r numerical-integration
person ywa136    schedule 08.11.2016    source источник
comment
может быть adaptIntegrate(NormalPDF, c(1,3), c(1,3))   -  person HubertL    schedule 08.11.2016

Ответы (1)


arrow_upward
2
arrow_downward

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

person IRTFM    schedule 08.11.2016
comment
Тег SO для [интеграции] явно не связан с математическим использованием этого слова. Я не вижу тега для интеграла или интеграции, поэтому не вставлял его. - person IRTFM; 08.11.2016