Два разных результата с двумя методами численного интегрирования?

Question

Я рассчитал интеграл от произведения гауссовой плотности и некоторой функции.

Сначала я сделал это с помощью функции int2() (rmutil package).А потом я сделал это с очками Гаусса-Эрмита.Два результата, которые я получил, разные.Должен ли я считать, что метод Гаусса-Эрмита является хорошим, а численное интегрирование является приближенным?

Я приведу ниже пример:

1.rmutil :: int2 ()

library(rmutil)

Sig <- matrix (c(0.2^2, 0, 0, 0.8^2), ncol=2)
Mu<- c(2, 0) 


to.integrate <- function(B0, B1) {
  first.int= 1/0.8 * (1.2 * exp(B0 + B1 * 0.5))^(-1/0.8) * gamma(1/0.8) 
  B=matrix(c(B0, B1), ncol=1)
  multi.norm=1 / (2 * pi * det(Sig)^(1/2)) * 
    exp (- 0.5 * t( B - Mu ) %*% solve(Sig) %*%( B - Mu ) )
  return (first.int %*% multi.norm)
}

result_int2 <- int2(to.integrate, a=c(-Inf, -Inf), b=c(Inf, Inf), 
                     eps=1.0e-6, max=16, d=5)

2.Вычислить многомерные гауссовы квадратурные точки:

library(statmod)
mgauss.hermite <- function(n, mu, sigma) {
  dm  <- length(mu)
  gh  <- gauss.quad(n, 'hermite')
  gh  <- cbind(gh$nodes, gh$weights)
  idx <- as.matrix(expand.grid(rep(list(1:n), dm)))
  pts <- matrix(gh[idx, 1], nrow(idx), dm)
  wts <- apply(matrix(gh[idx, 2], nrow(idx), dm), 1, prod)
  eig <- eigen(sigma) 
  rot <- eig$vectors %*% diag(sqrt(eig$values))
  pts <- t(rot %*% t(pts) + mu)
  return(list(points=pts, weights=wts))
}

nod_wei <- mgauss.hermite(10, mu=Mu, sigma=Sig)
gfun <- function(B0, B1) {
  first.int <- 1/0.8 *(1.2 * exp(B0 + B1 * 0.5))^(-1/0.8)* gamma(1/0.8) 
  return(first.int)
}

result_GH <- sum(gfun(nod_wei$points[, 1], nod_wei$points[, 2]) * nod_wei$weights)/pi

result_int2
result_GH

Answer 1

Ошибка возникла из-за способа вычисления точек в функции mgauss.hermite.

Я изменил декомпозицию матрицы Сигмы для декомпозиции Холецкого с умножением на квадратный корень из 2. И результатыиз двух методов стали очень похожи.Ниже приведена коррекция функции mgauss.hermite

mgauss.hermite <- function(n, mu, sigma) {
  dm  <- length(mu)
  gh  <- gauss.quad(n, 'hermite')
  gh  <- cbind(gh$nodes, gh$weights)
  idx <- as.matrix(expand.grid(rep(list(1:n),dm)))
  pts <- matrix(gh[idx,1],nrow(idx),dm)
  wts <- apply(matrix(gh[idx,2],nrow(idx),dm), 1, prod)
  rot <- 2.0**0.5*t(chol(sigma))
  pts <- t(rot %*% t(pts) + mu)
  return(list(points=pts, weights=wts))
}