polyCub.iso: Cubature of Isotropic Functions over Polygonal Domains

In bastistician/polyCub: Cubature over Polygonal Domains

Cubature of Isotropic Functions over Polygonal Domains

Description

polyCub.iso numerically integrates a radially symmetric function f(x,y) = f_r(||(x,y)-\boldsymbol{\mu}||), with \mu being the center of isotropy, over a polygonal domain. It internally approximates a line integral along the polygon boundary using integrate. The integrand requires the antiderivative of r f_r(r)), which should be supplied as argument intrfr (f itself is only required if check.intrfr=TRUE). The two-dimensional integration problem thereby reduces to an efficient adaptive quadrature in one dimension. If intrfr is not available analytically, polyCub.iso can use a numerical approximation (meaning integrate within integrate), but the general-purpose cubature method polyCub.SV might be more efficient in this case. See Meyer and Held (2014, Supplement B, Section 2.4) for mathematical details.

.polyCub.iso is a “bare-bone” version of polyCub.iso.

Usage

polyCub.iso(polyregion, f, intrfr, ..., center, control = list(),
  check.intrfr = FALSE, plot = FALSE)

.polyCub.iso(polys, intrfr, ..., center, control = list(),
  .witherror = FALSE)

Arguments

polyregion	a polygonal domain. The following classes are supported: "owin" from package spatstat.geom, "gpc.poly" from gpclib, "SpatialPolygons", "Polygons", and "Polygon" from package sp, as well as "(MULTI)POLYGON" from package sf. (For these classes, polyCub knows how to get an xylist.)
f	a two-dimensional real-valued function. As its first argument it must take a coordinate matrix, i.e., a numeric matrix with two columns, and it must return a numeric vector of length the number of coordinates.
intrfr	a function(R, ...), which implements the (analytical) antiderivative of r f_r(r) from 0 to R. The first argument must be vectorized but not necessarily named R. If intrfr is missing, it will be approximated numerically via integrate(function(r, ...) r * f(cbind(x0 + r, y0), ...), 0, R, ..., control=control), where c(x0, y0) is the center of isotropy. Note that f will not be checked for isotropy.
...	further arguments for f or intrfr.
center	numeric vector of length 2, the center of isotropy.
control	list of arguments passed to integrate, the quadrature rule used for the line integral along the polygon boundary.
check.intrfr	logical (or numeric vector) indicating if (for which r's) the supplied intrfr function should be checked against a numeric approximation. This check requires f to be specified. If TRUE, the set of test r's defaults to a seq of length 20 from 1 to the maximum absolute x or y coordinate of any edge of the polyregion.
plot	logical indicating if an image of the function should be plotted together with the polygonal domain, i.e., plotpolyf(polyregion, f, ...).
polys	something like owin$bdry, but see xylist.
.witherror	logical indicating if an upper bound for the absolute integration error should be attached as an attribute to the result?

Value

The approximate integral of the isotropic function f over polyregion.
If the intrfr function is provided (which is assumed to be exact), an upper bound for the absolute integration error is appended to the result as attribute "abs.error". It equals the sum of the absolute errors reported by all integrate calls (there is one for each edge of polyregion).

Author(s)

Sebastian Meyer

The basic mathematical formulation of this efficient integration for radially symmetric functions was ascertained with great support by Emil Hedevang (2013), Dept. of Mathematics, Aarhus University, Denmark.

References

Hedevang, E. (2013). Personal communication at the Summer School on Topics in Space-Time Modeling and Inference (May 2013, Aalborg, Denmark).

Meyer, S. and Held, L. (2014). Power-law models for infectious disease spread. The Annals of Applied Statistics, 8 (3), 1612-1639. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1214/14-AOAS743")}

See Also

system.file("include", "polyCubAPI.h", package = "polyCub") for a full C-implementation of this cubature method (for a single polygon). The corresponding C-routine polyCub_iso can be used by other R packages, notably surveillance, via ‘⁠LinkingTo: polyCub⁠’ (in the ‘DESCRIPTION’) and ‘⁠#include <polyCubAPI.h>⁠’ (in suitable ‘/src’ files). Note that the intrfr function must then also be supplied as a C-routine. An example can be found in the package tests.

Other polyCub-methods: polyCub.SV(), polyCub.exact.Gauss(), polyCub.midpoint(), polyCub()

Examples

## we use the example polygon and f (exponential decay) from
example(plotpolyf)

## numerical approximation of 'intrfr' (not recommended)
(intISOnum <- polyCub.iso(letterR, f, center = fcenter))

## analytical 'intrfr'
## intrfr(R) = int_0^R r*f(r) dr, for f(r) = dexp(r), gives
intrfr <- function (R, rate = 1) pgamma(R, 2, rate) / rate
(intISOana <- polyCub.iso(letterR, f, intrfr = intrfr, center = fcenter,
                          check.intrfr = TRUE))
## f is only used to check 'intrfr' against a numerical approximation

stopifnot(all.equal(intISOana, intISOnum, check.attributes = FALSE))


### polygon area: f(r) = 1, f(x,y) = 1, center does not really matter

## intrfr(R) = int_0^R r*f(r) dr = int_0^R r dr = R^2/2
intrfr.const <- function (R) R^2/2
(area.ISO <- polyCub.iso(letterR, intrfr = intrfr.const, center = c(0,0)))

if (require("spatstat.geom")) { # check against area.owin()
    stopifnot(all.equal(area.owin(owin(poly = letterR)),
                        area.ISO, check.attributes = FALSE))
}

bastistician/polyCub documentation built on Feb. 3, 2024, 5:41 p.m.