surveillance (version 1.5-4)

polyCub: Two-Dimensional Numerical Integration over a Polygonal Domain

Description

The functions documented here numerically integrate a two-dimensional function over a polygonal domain. The first method uses the classical midpoint rule converting the polygons to binary pixel images using the as.im.function method from package spatstat. The second one uses the more sophisticated product Gauss cubature proposed by Sommariva & Vianello (2007). The third one is a quasi-exact method specific to the integration of the bivariate Gaussian density over polygonal domains and is based on formulae from Chapter 26 of the famous Abramowitz & Stegun (1970) handbook, i.e. triangulation of the polygonal domain and appropriate evaluations of pmvnorm from package mvtnorm. See Section 3.2 of Meyer (2010) for a more detailed description and benchmark experiment of those three (among others).

Usage

polyCub(polyregion, f, method = c("SV", "midpoint", "exact.Gauss"), ...,
        plot = FALSE)
polyCub.midpoint(polyregion, f, ...,
                 eps = NULL, dimyx = NULL, plot = FALSE)
polyCub.SV(polyregion, f, ...,
           nGQ = 20, alpha = NULL, rotation = FALSE, plot = FALSE)
polyCub.exact.Gauss(polyregion, mean = c(0,0), Sigma = diag(2),
                    plot = FALSE)

Arguments

polyregion

a polygon representing the integration domain. For polyCub.midpoint, polyregion can be anything coercible to the spatstat class "owin" (by

f

two-dimensional function to be integrated. As its first argument the function must take a coordinate matrix, i.e. a matrix with two columns. For the polyCub.exact.Gauss method, f is ignored since it is specific to the biv

method

choose one of the implemented cubature methods (partial argument matching is applied). Defaults to using the product Gauss cubature implemented in polyCub.SV.

...

further arguments passed to f.

eps

numeric scalar specifying the (approximate) width/height of the pixels for the generation of the image.

dimyx

vector of length two specifying the number of subdivisions in each dimension for the generation of the image. See as.im.function and as.mask.

nGQ

degree of the one-dimensional Gauss-Legendre quadrature rule (default: 20). See gauss.quad in package statmod, on which this function depends.

alpha

base-line of the (rotated) polygon at $x = alpha$ (see Sommariva & Vianello (2007) for an explication). If NULL (default), the midpoint of the x-range of the polygon is chosen if no rotation is performed, and otherwis

rotation

logical (default: FALSE) or a list of points "P" and "Q" describing the preferred direction. If TRUE, the polygon is rotated according to the vertices "P" and "Q", w

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.

Value

• The approximated value of the integral of f over polyregion. For the polyCub.exact.Gauss method, two attributes are appended to the integral value: [object Object],[object Object]

encoding

latin1

References

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

Meyer, S. (2010): Spatio-Temporal Infectious Disease Epidemiology based on Point Processes. Master's Thesis, Ludwig-Maximilians-Universit�{ae}t M�{ue}nchen. Available as http://epub.ub.uni-muenchen.de/11703/ Sommariva, A. and Vianello, M. (2007): Product Gauss cubature over polygons based on Green's integration formula. Bit Numerical Mathematics, 47 (2), 441-453.

Examples

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

# simple polygonal integration domain
disc.owin <- discpoly(c(3,2), 5, npoly=8, class="owin")

# plot image of the function and integration domain
gr <- seq(-8,8,by=0.2)
vals <- matrix(f(expand.grid(gr,gr)), nrow=length(gr))
image(x=gr, y=gr, z=vals)
library("spatstat")    
plot.owin(disc.owin, hatch=TRUE, add=TRUE)

# Two-dimensional midpoint rule
testmidpoint <- function (eps, main=paste("Two-dimensional midpoint rule with eps =",eps))
{
    image(x=gr, y=gr, z=vals, main=main)
    plot.owin(disc.owin, hatch=TRUE, add=TRUE)
    # add binary mask to plot
    #plot(as.mask(disc.owin, eps=eps), add=TRUE, box=FALSE)
    # add evaluation points to plot
    with(as.mask(disc.owin, eps=eps), points(expand.grid(xcol, yrow), col=m, pch=20))
    polyCub.midpoint(disc.owin, f, eps=eps)
}
testmidpoint(5)
testmidpoint(3)
testmidpoint(0.5)
testmidpoint(0.2)

# Product Gauss cubature by Sommariva & Vianello (2007)
if (require("statmod")) {
    for (nGQ in c(1:5,10,20,60)) {
        cat("nGQ =", sprintf("%2i",nGQ), ": ",
            format(polyCub.SV(disc.owin, f, nGQ=nGQ), digits=16),
            "")
    }
}

# quasi-exact cubature by triangulation (gpclib) and using mvtnorm::pmvnorm()
if (require("mvtnorm") && require("gpclib")) {
    oopt <- surveillance.options(gpclib=TRUE)
    print(polyCub.exact.Gauss(disc.owin, mean=c(0,0), Sigma=5^2*diag(2),
                              plot=TRUE), digits=16)
    surveillance.options(oopt)
}