Integral approximation using gauss-hermit quadrature method: mvQuad package

user_6403044 0
Viewed 81

I want to approximate the integral corresponds to $E(XY)$ where X and Y are independent and X~N(0.5,1) and Y~N(0.5,1), using gauss-hermit quadrature using mvQuad package.

Since the two random variables are normal distributed, the approximated value should be close to 0.25.

But I am getting a different results.

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

Can anybody helps me to figure out what I am doing wrong?

Thank you

integral numerical-integration r
1 Answers

You should use the product of PDFs for both X and Y in function myFun2d 

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

and you will see

> 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
user_121587570
Related