使用 gauss-hermit 求积法的积分逼近:mvQuad 包

Integral approximation using gauss-hermit quadrature method: mvQuad package

我想近似对应于 $E(XY)$ 的积分,其中 X 和 Y 是独立的,X~N(0.5,1) 和 Y~N(0.5,1),使用 gauss-hermit 正交使用mvQuad 包。

由于两个随机变量呈正态分布,近似值应该接近0.25。

但我得到了不同的结果。

> require(mvQuad)
> nw <- createNIGrid(dim=2, type="GHe", level=c(10,10))
> m=c(0.5,0.5)
> c=matrix(c(1,0,0,1),nrow = 2,byrow = F)
> rescale(nw, m = m, C = c, dec.type = 0)
> 
> myFun2d <- function(x){
+     (0.5+sqrt(2)*x[,1])*(0.5+sqrt(2)*x[,2])/pi
+ }
> 
> quadrature(myFun2d, grid = nw)
[1] 58.71479
> 

谁能帮我弄清楚我做错了什么?

谢谢

您应该在函数 myFun2d

中为 XY 使用 PDF 的乘积
myFun2d <- function(x) {
  x[, 1] * exp(-0.5 * (x[, 1]-0.5)^2) / sqrt(2 * pi) * x[, 2] * exp(-0.5 * (x[, 2]-0.5)^2) / sqrt(2 * pi)
}

你会看到

> library(mvQuad)

> # create grid
> nw <- createNIGrid(dim = 2, type = "GHe", level = c(10, 10))

> m <- c(0.5, 0.5)

> c <- matrix(c(1, 0, 0, 1), nrow = 2, byrow = F)

> rescale(nw, m = m, C = c, dec.type = 0)

> # define the integrand
> myFun2d <- function(x) {
+   x[, 1] * exp(-0.5 * (x[, 1] - 0.5)^2) / sqrt(2 * pi) * x[, 2] * exp(-0.5 * (x[, 2] - 0.5)^2) / .... [TRUNCATED]

> # compute the approximated value of the integral
> (A <- quadrature(myFun2d, grid = nw))
[1] 0.25