CODEWITHSUNDEEP
rintegral
Kam
Kam

Reputation: 21

Evaluate bivariate normal distribution with Gauss-Hermite quadrature

By using the library(statmod), I can evaluate a univariate normal distribution using Gauss-Hermite quadrature. But how can I evaluate a bivariate normal distribution using Gauss-Hermite quadrature?

Any help will be appreciated. Thanks in advance.

I have given the codes that I used for evaluating univariate normal using Gauss-Hermite 5 points. But how I can I do for bivariate normal distribution?

library(statmod)
## generating Gauss-Hermite quadrature points and weights
q=gauss.quad(n=5,kind="hermite")

## defining univariate normal function
mu=0
sigma=2
norm=function(b){
    M=((2*pi*sigma)^(-1/2))*exp(-(1/2)*(b^2/sigma^2))
    return(M)
}

## approximating the integral of norm(b) using Gauss-Hermite method
sum(q$weights*norm(q$nodes)*exp(q$nodes^2))

Upvotes: 1

Views: 1172

Answers (1)

Maurits Evers
Maurits Evers

Reputation: 50678

Here is a minimal worked-through example, based on a bivariate standard normal (probability) density with a given covariance matrix sigma and the quadrature rules implemented in the mvQuad library.

Note that we use the Gauss-Legendre quadrature rule, which allows integration over an arbitrary bounded domain because in the Gauss-Hermite quadrature the domain is unbounded from (-∞, +∞). Since we're working with a bivariate standard normal probability density, the integral over the unbounded domain trivially equals 1. Generally, mvQuad::createNIGrid allows the implementation of various quadrature rules, including the Gauss-Hermite quadrature (see ?createNIGrid for details).

  1. Define a covariance matrix for the a bivariate standard normal probability density

    library(mvtnorm)
sigma <- matrix(c(1, 0.2, 0.2, 1), ncol = 2)
dens <- function(x) dmvnorm(x, sigma = sigma)

    We are interested in the integral of dens in the domain x ϵ [-1, 2] and y ϵ [-1, 2].

  2. We follow instructions from the mvQuad vignette to create a grid and rescale to the domain of interest

    library(mvQuad)
grid <- createNIGrid(dim = 2, type = "GLe", level = 6)
rescale(grid, domain = rbind(c(-1, 2), c(-1, 2)))

  3. Calculate the integral of the bivariate normal in the domain x ϵ [-1, 2] and y ϵ [-1, 2]

    quadrature(function(x) dmvnorm(x, sigma = sigma), grid = grid)
#[1] 0.6796583

  4. This value is in good agreement with the value from pmvnorm (which computes the distribution function of the multivariate normal for arbitrary limits and covariance matrices)

    pmvnorm(lower = c(-1, -1), upper = c(2, 2), sigma = sigma)
#[1] 0.6796584
#attr(,"error")
#[1] 1e-15
#attr(,"msg")
#[1] "Normal Completion"

Upvotes: 3

Related Questions