# mvQuad: Methods for Multivariate Quadrature

link to CRAN: https://cran.r-project.org/web/packages/mvQuad/

This package provides a collection of methods for (potentially) multivariate quadrature in R. It’s especially designed for use in statistical problems, characterized by integrals of the form $$\int_a^b g(x)p(x) \ dx$$, where $$p(x)$$ denotes a density-function. Furthermore the methods are also applicable to standard integration problems with finite, semi-finite and infinite intervals.

In general quadrature (also named: numerical integration, numerical quadrature) work this way: The integral of interests is approximated by a weighted sum of function values.

$I = \int_a^b h(x) \ dx \approx A = \sum_{i=1}^n w_i \cdot h(x_i)$

The so called nodes ($$x_i$$) and weights($$w_i$$) are defined by the chosen quadrature rule, which should be appropriate (better: optimal) for the integration problem in hand1.

This principle can be transferred directly to the multivariate case.

The methods provided in this package cover the following tasks:

• creating an uni-/multivariate grid (grid: collection of nodes and weights) for a chosen quadrature rule
• examining the created grid
• rescaling the grid for appropriate/efficient use
• computing the approximated integral

## Quick Start

This section shows a typical workflow for quadrature with the mvQuad-package. More details and additional features of this package are provided in the subsequent sections. In this illustration the following two-dimensional integral will be approximated: $I = \int_1^2 \int_1^2 x \cdot e^y \ dx dy$

{r } library(mvQuad)

# create grid nw <- createNIGrid(dim=2, type=“GLe”, level=6)

# rescale grid for desired domain rescale(nw, domain = matrix(c(1, 1, 2, 2), ncol=2))

# define the integrand myFun2d <- function(x){ x[,1]*exp(x[,2]) }

# compute the approximated value of the integral A <- quadrature(myFun2d, grid = nw) 

Explanation Step-by-Step

1. mvQuad-package is loaded
2. with the createNIGrid-command a two-dimensional (dim=2) grid, based on Gauss-Legendre quadrature rule (type="GLe") with a given accuracy level (level=6) is created and stored in the variable nw

The grid created above is designed for the domain $$[0, 1]^2$$ but we need one for the domain $$[1, 2]^2$$

1. the command rescale changes the domain-feature of the grid; the new domain is passed in a matrix (domain=...)

2. the integrand is defined

3. the quadrature-command computes the weighted sum of function values as mentioned in the introduction

## Supported Rules (build in)

The choice of quadrature rule is heavily related to the integration problem. Especially the domain/support ($$[a, b]$$ (finite), $$[a, \infty)$$ (semi-finite), $$(-\infty, \infty)$$ (infinite)) determines the choice.

The mvQuad-packages provides the following quadrature rules.

• cNC1, cNC2, ..., cNC6: closed Newton-Cotes Formulas of degree 1-6 (1=trapezoidal-rule; 2=Simpson’s-rule; …), finite domain

• oNC0, oNC1, ..., oNC3: open Newton-Cote Formula of degree 0-3 (0=midpoint-rule; …), finite domain

• GLe, GKr: Gauss-Legendre and Gauss-Kronrod rule, finite domain

• nLe: nested Gauss-Legendre rule (finite domain) [@Petras2003]

• Leja: Leja-Points (finite domain)

• GLa: Gauss-Laguerre rule (semi-finite domain)

• GHe: Gauss-Hermite rule (infinite domain)

• nHe: nested Gauss-Hermite rule (infinite domain) [@GenzKeister1996]

• GHN, nHN: (nested) Gauss-Hermite rule as before but weights are pre-multiplied by the standard normal density ($$\hat{w}_i = w_i * \phi(x_i)$$).2

For each rule grids can be created of different accuracy. The adjusting screw in the createNIGrid is the level-option. In general, the higher level the more precise the approximation. For the Newton-Cotes rules an arbitrary level can be chosen. The other rules uses lookup-tables for the nodes and weights and are therefore restricted to a maximum level (see help(QuadRules))

1. A rigorous introduction to numerical integration can be found in Davis/Rabinowitz [-@DavisRabinowitz1984]↩︎

2. Those rules are computationally more efficent for integrands of the form: $$\int_{-\infty}^{\infty} g(x)\phi(x)dx$$ because the approximation reduces to $$\sum \hat{w}_i \cdot g(x)$$.↩︎