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
polyCub.exact.Gauss(polyregion, mean = c(0, 0), Sigma = diag(2), plot = FALSE)
mean and covariance matrix of the bivariate normal density to be integrated.
logical indicating if an illustrative plot of the numerical
integration should be produced. Note that the
The integral of the bivariate normal density over
Two attributes are appended to the integral value:
number of triangles over which the standard bivariate normal density had to
be integrated, i.e. number of calls to
Approximate absolute integration error stemming from the error introduced by
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
Abramowitz, M. and Stegun, I. A. (1972). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. New York: Dover Publications.
circleCub.Gauss for quasi-exact cubature of the
isotropic Gaussian density over a circular domain.
## 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() ## Not run: ## (this example requires gpclib and acceptance of its license) 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) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.