polyCub: polyCub.exact.Gauss – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

polyCub.exact.Gauss

Quasi-Exact Cubature of the Bivariate Normal Density

Description

The bivariate Gaussian density can be integrated based on a triangulation of the (transformed) polygonal domain, using formulae from the Abramowitz and Stegun (1972) handbook (Section 26.9, Example 9, pp. 956f.). This method is quite cumbersome because the A&S formula is only for triangles where one vertex is the origin (0,0). For each triangle of the tristrip we have to check in which of the 6 outer regions of the triangle the origin (0,0) lies and adapt the signs in the formula appropriately: (AOB+BOC-AOC) or (AOB-AOC-BOC) or (AOB+AOC-BOC) or (AOC+BOC-AOB) or .... However, the most time consuming step is the evaluation of pmvnorm.

Usage

polyCub.exact.Gauss(polyregion, mean = c(0, 0), Sigma = diag(2), plot = FALSE)

Arguments

`polyregion`	a `"gpc.poly"` polygon or something that can be coerced to this class, e.g., an `"owin"` polygon (via `owin2gpc`), an `"sfg"` polygon (via `sfg2gpc`), or – given rgeos is available – a `"SpatialPolygons"` object.
`mean, Sigma`	mean and covariance matrix of the bivariate normal density to be integrated.
`plot`	logical indicating if an illustrative plot of the numerical integration should be produced. Note that the `polyregion` will be transformed (shifted and scaled).

Value

The integral of the bivariate normal density over polyregion. Two attributes are appended to the integral value:

`nEval`	number of triangles over which the standard bivariate normal density had to be integrated, i.e. number of calls to `pmvnorm` and `pnorm`, the former of which being the most time-consuming operation.
`error`	Approximate absolute integration error stemming from the error introduced by the `nEval` `pmvnorm` evaluations. For this reason, the cubature method is in fact only quasi-exact (as is the `pmvnorm` function).

Note

The package gpclib is required to produce the tristrip, since this is not implemented in rgeos (as of version 0.3-25). The restricted license of gpclib (commercial use prohibited) has to be accepted explicitly via gpclibPermit() prior to using polyCub.exact.Gauss.

References

Abramowitz, M. and Stegun, I. A. (1972). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. New York: Dover Publications.

Examples

## a function to integrate (here: isotropic zero-mean Gaussian density)
f <- function (s, sigma = 5)
    exp(-rowSums(s^2)/2/sigma^2) / (2*pi*sigma^2)

## a simple polygon as integration domain
hexagon <- list(
    list(x = c(7.33, 7.33, 3, -1.33, -1.33, 3),
         y = c(-0.5, 4.5, 7, 4.5, -0.5, -3))
)

## quasi-exact integration based on gpclib::tristrip() and mvtnorm::pmvnorm()
if (requireNamespace("mvtnorm") && gpclibPermit()) {
    hexagon.gpc <- new("gpc.poly", pts = lapply(hexagon, c, list(hole = FALSE)))
    plotpolyf(hexagon.gpc, f, xlim = c(-8,8), ylim = c(-8,8))
    print(polyCub.exact.Gauss(hexagon.gpc, mean = c(0,0), Sigma = 5^2*diag(2),
                              plot = TRUE), digits = 16)
}

polyCub

Cubature over Polygonal Domains

v0.8.0

GPL-2

Authors

Sebastian Meyer [aut, cre, trl] (<https://orcid.org/0000-0002-1791-9449>), Leonhard Held [ths], Michael Hoehle [ths]

Initial release

2021-01-26