mvQuad: Methods for Multivariate Quadrature

link to CRAN:

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:

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.

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)\).↩︎