--- title: "mvQuad - Package for multivariate Quadrature" author: "Constantin Weiser" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{mvQuad_intro} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} references: - id: DavisRabinowitz1984 title: Methods of numerical integration author: - family: Davis given: P.J. - family: Rabinowitz given: P. publisher: Academic Press type: book issued: year: 1984 - id: Gerstner1998 title: Numerical integration using sparse grids author: - family: Gerstner given: T. - family: Griebel given: M. container-title: Numerical Algorithms volume: 18 page: 209-232 type: article-journal issued: year: 1998 - id: Heiss2008 title: Likelihood approximation by numerical integration on sparse grids author: - family: Heiss given: F. - family: Winschel given: V. container-title: Journal of Econometrics volume: 144 page: 62-80 type: article-journal issued: year: 2008 - id: Jackel2005 title: A note on multivariate Gauss-Hermite quadrature author: - family: Jäckel given: P. type: article-journal container-title: mimeo issued: year: 2005 - id: GenzKeister1996 title: Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight author: - family: Genz given: A. - family: Keister given: B.D. container-title: Journal of Computational and Applied Mathematics volume: 71 page: 299-309 type: article-journal issued: year: 1996 - id: Petras2003 title: Smolyak cubature of given polynomial degree with few nodes for increasing dimension author: - family: Petras given: K. container-title: Numerische Mathematik volume: 93 page: 729-753 type: article-journal issued: year: 2003 --- ## Introduction 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 hand^[A rigorous introduction to numerical integration can be found in Davis/Rabinowitz [-@DavisRabinowitz1984]]. 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** 0. `mvQuad`-package is loaded 1. 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$ 2. the command `rescale` changes the domain-feature of the grid; the new domain is passed in a matrix (`domain=...`) 3. the integrand is defined 4. 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)$).^[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)$.] 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)`) ## User-defined Rules Through the very general functionality of quadrature, `mvQuad` allows to use user-defined quadrature rules. This can be of interest for specific integration problems. There are two ways to hand over user-defined rules to `createNIGrid`. (1) Via a function, which takes a value for the level and returns a grid. The returned grid have to be a list with three objects (n: nodes, w: weights, features: initial domain; see the example below). (2) Via a text-file, which contains the nodes and weights for potentially several levels. **Variant 1** ```{r } # via an user-defined function myRule.fun <- function(l){ n <- seq(1, 2*l-1, by=2)/ (l*2) w <- rep(1/(l), l) initial.domain <- matrix(c(0,1), ncol=2) return(list(n=as.matrix(n), w=as.matrix(w), features=list(initial.domain=initial.domain))) } nw.fun <- createNIGrid(d=1, type = "myRule.fun", level = 10) print(data.frame(nodes=getNodes(nw.fun), weights=getWeights(nw.fun))) ``` **Variant 2** ```{r } # via a text-file myRule.txt <- readRule(file=system.file("extdata", "oNC0_rule.txt", package = "mvQuad")) nw.txt <- createNIGrid(d=1, type = myRule.txt, level = 10) print(data.frame(nodes=getNodes(nw.txt), weights=getWeights(nw.txt))) ``` The text-file for variant 2 have to be formatted like the following example. The first line have to declare the domain `initial.domain a b`, where `a` and `b` denotes the lower and upper-bound for the integration domain. This can be either a number or '-Inf'/'Inf' (for example `initial.domain 0 1` or `initial.domain 0 Inf`) Every following line contains one single node and weight belonging to one level of the rule (format: 'level' 'node' 'weight'). This example shows the use for the "midpoint-rule" (levels: 1 - 3). ```{r echo=FALSE} txt <- readLines(system.file("extdata", "oNC0_rule.txt", package = "mvQuad")) for (i in 1:9) { cat(i,"\t", txt[i], "\n") if (i==9) { cat(" ... \n") } } ``` ## Construction of multivariate grids The quadrature rules described above are inherent 1-dimensional. There exists several construction principles for multivariate grids, based on 1D-quadrature rules. The traditional one is the so called "product-rule", which leads to an even designed grid (see figure "Product-Rule"). The main drawback of this principle is the fast growing number of grid points for an increasing number of dimensions. The total number of grid points $N = n^d$, where $n$ denotes the number of points in the underlying 1d-rule. This property make this principle _unfeasible_ for higher dimensions (d>4). `createNIGrid` can be set to use this principle with the `ndConstruction="product"` argument ```{r fig.show="hold", fig.width=4.5} nw <- createNIGrid(dim=2, type="cNC1", level=5, ndConstruction = "product") plot(nw, main="Example: Product-Rule") ``` Another principle is the "combination technique", which leads to an more sophisticated compilation of grid points. Therefore the construction is a bit more complex [@Gerstner1998, @Heiss2008], but the number of grid points grow much slower with increasing dimensionality, with almost the same accuracy. Therefore _sparse grids_ are feasible for more then 20 dimensions `creatNIGrid` can be forced to use this principle with the `ndConstruction="sparse" argument. ```{r fig.show="hold", fig.width=4.5} nw <- createNIGrid(dim=2, type="cNC1", level=5, ndConstruction = "sparse") plot(nw, main="Example: Combination Technique") ``` ## Rescaling `createNIGrid` generate an un-adjusted grid. Un-adjusted in the sense that it maybe doesn't fit the domain needed for the integration problem in hand or it is not efficient. Typically build in grids for finite intervals are adjusted for $[0, 1]^d$ and grids for infinite intervals for a density with zero mean and variance equals one. A first example for rescaling was shown in the Quick-Start section. Another one is shown here. A bivariate normal density with mean-vector $m=(2,1)'$ and covariance matrix $C=\begin{pmatrix} 1.0 & 0.6 \\ 0.6 & 2.0 \end{pmatrix}$ should be integrated for $[-\infty, \infty]^2$. The following four figures shows the adjusted grid added by the contour plot of the integrand. ```{r, fig.height=6, fig.width=6, echo=FALSE, message=FALSE} nw <- createNIGrid(dim=2, type="GHe", level=3) C = matrix(c(1,0.6,0.6, 2),2) m = c(2, 1) dmvnorm <- function(x,m,C){ C.inv <- solve(C) tmp <- apply(x, MARGIN = 1, FUN = function(x){ 1/sqrt((2*pi)^2 * det(C)) * exp(-0.5 * t(x-m)%*%C.inv%*%(x-m)) }) return(unlist(tmp)) } xx <- expand.grid(seq(-5,5,length.out=100),seq(-5,5,length.out=100)) yy <- dmvnorm(xx, m, C) yy <- matrix(yy,100) par(mfrow=c(2,2), mar=c(2,2,3,1), oma=c(0,0,0,0)) plot(nw, main="initial grid", xlim=c(-6,6), ylim=c(-6,6), pch=8) contour(seq(-5,5,length.out=100),seq(-5,5,length.out=100),yy,asp=1, labels="", add = T, col="gray") rescale(nw, m = m, C = C, dec.type = 0) plot(nw, main="no dec.", xlim=c(-6,6), ylim=c(-6,6), pch=8) contour(seq(-5,5,length.out=100),seq(-5,5,length.out=100),yy,asp=1, labels="", add = T, col="gray") rescale(nw, m = m, C = C, dec.type = 1) plot(nw, main="spectral dec.", xlim=c(-6,6), ylim=c(-6,6), pch=8) contour(seq(-5,5,length.out=100),seq(-5,5,length.out=100),yy,asp=1, labels="", add = T, col="gray") rescale(nw, m = m, C = C, dec.type = 2) plot(nw, main="Cholesky dec.", xlim=c(-6,6), ylim=c(-6,6), pch=8) contour(seq(-5,5,length.out=100),seq(-5,5,length.out=100),yy,asp=1, labels="", add = T, col="gray") ``` The upper left figure shows the initial grid (optimal for zero mean an variance equals 1). The upper right figure shows a rescaled grid, where for the mean and variances (diagonal elements of $C$) is accounted for. In both lower pictures it's also accounted for the covariance, left via the spectral-decomposition, right via the Cholesky decomposition [@Jackel2005]. Here is the code for the rescaling: ```{r, fig.height=6, fig.width=6, eval=FALSE} nw <- createNIGrid(dim=2, type="GHe", level=3) # no rescaling plot(nw, main="initial grid", xlim=c(-6,6), ylim=c(-6,6), pch=8) rescale(nw, m = m, C = C, dec.type = 0) plot(nw, main="no dec.", xlim=c(-6,6), ylim=c(-6,6), pch=8) rescale(nw, m = m, C = C, dec.type = 1) plot(nw, main="spectral dec.", xlim=c(-6,6), ylim=c(-6,6), pch=8) rescale(nw, m = m, C = C, dec.type = 2) plot(nw, main="Cholesky dec.", xlim=c(-6,6), ylim=c(-6,6), pch=8) ``` ## Examinig the Grid For technical or didactic reasons there is sometimes the need to examine a created grid. The `mvQuad`-package provides the methods illustrated in the following example-session ```{r } nw <- createNIGrid(dim=2, type="GHe", level=6, ndConstruction = "sparse") plot(nw) print(nw) dim(nw) size(nw) ``` ## Miscellaneous 1. What `mvQuad` (V. 1.0.1) can't do: + `mvQuad` can't select an appropriate quadrature rule for your problem. + The grid create by 'createNIGrid' is static. There is no refinement nor adaption. + There is no error estimate, which gives you a feeling of the achieved accuracy. 2. Save created grids. The creation of high-dimensional grids is 'time-consuming'. Ones you are satisfied by an grid for a concrete task, save the grid (i.e. `save(nw, file=....)`) for replications. 3. Ex ante it's not trivial to find a balance between 'standard grids' and 'highly adapted grids'. Both approaches 'brute force' and 'high specialization' can cost a lot of CPU- and/or live-time. 4. All comments regarding the package or the documentation are highly appreciated. ## References