两种数值积分方法的两种不同结果？

Question

我计算了高斯密度和某个函数的乘积的积分。

首先，我使用函数 int2()（rmutil 包）完成了它。然后，我用 Gauss-Hermite 点来做。我得到的两个结果是不同的。我是否应该认为 Gauss-Hermite 方法是好的，数值积分是一种近似？

我在下面提供一个例子：

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 函数中计算点的方式。

我将 Sigma 矩阵的分解更改为乘以 2 的平方根的 Cholesky 分解。并且两种方法的结果变得非常相似。下面是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))
}

两种数值积分方法的两种不同结果？

Two different results with two methods of numerical integration?

r

gaussian

numerical-methods