Package {quadriceps}

quadriceps: positive-weight cubature rules for the Gaussian weight and the cube

Description

Positive-weight cubature rules in d >= 1 dimensions for two weight functions:

Details

• GH (ghpos()): the Gaussian weight, by default the standard normal density N(0, I_d);

• Le (lepos()): the uniform weight, by default the uniform density on ⁠[0,1]^d⁠.

A rule of degree p integrates every polynomial of total degree ⁠<= p⁠ exactly. All weights are strictly positive. The rules shipped with the package are the smallest ones known to its author; available_rules() lists them, and LICENSE.note and the files RULES.md and NOTICE.md of the repository say where each one comes from.

The second argument q of ghpos() and lepos() is the number of nodes of the one-dimensional Gauss rule whose exactness is wanted, as in statmod::gauss.quad(q, kind = "hermite"); the rule returned has degree ⁠p = 2q - 1⁠. The degree can be given instead, as the argument p.

This package is the R twin of the Julia package Quadriceps.jl and of the Python package quadriceps; the three share their data and their conventions.

Author(s)

Maintainer: Joris Pinkse pinkse@gmail.com [copyright holder]

Authors:

• Joris Pinkse pinkse@gmail.com [copyright holder]

Other contributors:

• Moustapha Diallo (some of the stored rules for the cube descend from his published rules (MIT license, see LICENSE.note)) [copyright holder]

See Also

Useful links:

• https://github.com/NittanyLion/quadriceps-r

• Report bugs at https://github.com/NittanyLion/quadriceps-r/issues

Examples

r <- ghpos(3, 4)        # as exact as the 4x4x4 Gauss-Hermite grid (degree 7), with 27 nodes
dim(r$nodes)
r <- ghpos(3, p = 7)    # the same rule, requested by degree
r <- lepos(2, 5)        # degree 9 on the unit square: 17 nodes instead of 25

List the stored rules

Description

All stored rules of a family, one row per rule, sorted by dimension and degree. Tensor products, which pragmatic = TRUE builds on demand, are not listed.

Usage

available_rules(family)

Arguments

Value

A data frame with columns

• family, d (dimension), q, p (degree of exactness, ⁠p = 2q - 1⁠), n (number of nodes);

• moller: Möller's lower bound on n, or -1 where it is not tabulated; a rule with n == moller is proven minimal;

• relerr: largest relative monomial error of the stored double-precision rule, measured when the data were built (see exactness_error());

• minweight: smallest weight (always ⁠> 0⁠);

• interior: TRUE when every node lies inside the integration domain;

• origin: who the rule belongs to: "own", or "derived: ...", "same-rule: ...", "transcribed: ..." followed by the published source;

• source_id: the origin as the small integer stored in rules.bin (0: own; 3: derived from Diallo and Worku; 10 and up: a published rule; see FORMAT.md in the repository).

Examples

head(available_rules("gh")[, c("d", "q", "p", "n")])

Largest relative monomial error of a rule

Description

The maximum, over all monomials x^a of total degree ⁠<= p⁠, of ⁠|sum_i w_i x_i^a - E x^a| / max(sum_i |w_i| |x_i^a|, 1)⁠, for the normalized weight of family ("gh": N(0, I_d); "le": uniform on ⁠[0,1]^d⁠). The denominator is the scale of the sum being computed, so a value near machine epsilon means the rule is exact to rounding. Rules returned with normalize = FALSE must be checked in the normalized frame.

Usage

exactness_error(nodes, weights = NULL, p, family)

Arguments

Value

The error, a single number.

Examples

r <- ghpos(4, 5)                       # degree 2 * 5 - 1 = 9
exactness_error(r, p = 9, family = "gh")

Positive-weight cubature rule for the Gaussian weight

Description

The smallest positive-weight rule the package has for the Gaussian weight in d dimensions.

Usage

ghpos(d, q = NULL, p = NULL, normalize = TRUE, pragmatic = FALSE)

Arguments

Details

ghpos(d, q) follows one-dimensional Gauss-Hermite functions such as statmod::gauss.quad(q, kind = "hermite"): there, q is the number of nodes of the one-dimensional Gauss rule, which is exact to degree ⁠2q - 1⁠. Here, ghpos(d, q) returns a rule with the same exactness in d dimensions, degree ⁠p = 2q - 1⁠, that replaces the q^d-node product grid; for d = 1 it is the q-node Gauss-Hermite rule itself. Alternatively give the degree: ghpos(d, q) is ghpos(d, p = 2 * q - 1).

With rule <- ghpos(d, q), sum(rule$weights * f(rule$nodes)) (for a vectorized f of the rows) equals the integral of f against the weight for every polynomial f of total degree ⁠<= p⁠, up to rounding; ruleinfo() gives the measured error of each stored rule.

Rules are stored for ⁠2 <= d <= 5⁠ at odd degrees, up to a ceiling that depends on d (see available_rules()). A request is served by the smallest stored rule of degree ⁠>= p⁠, so an even p gets the rule for p + 1. For d = 1 the Gauss rule is computed on the fly (Golub-Welsch) and is never an error.

Value

A list with components nodes, an ⁠n x d⁠ matrix with one node per row (also for d = 1), and weights, a vector of n strictly positive weights.

See Also

lepos(), nnodes(), ruleinfo(), available_rules()

Examples

r <- ghpos(3, 4)                                   # degree 7: 27 nodes instead of 4^3 = 64
sum(r$weights * r$nodes[, 1]^2 * r$nodes[, 2]^4)   # E[Z1^2 Z2^4] = 3
identical(ghpos(3, p = 7), r)

try(ghpos(3, 21))                                  # nothing stored at that degree
r <- ghpos(3, 21, pragmatic = TRUE)                # tensor product of lower-dimensional rules
r <- ghpos(7, 5, pragmatic = TRUE)                 # d = 7 as (d = 2) x (d = 5)

Positive-weight cubature rule for the uniform weight on the cube

Description

The smallest positive-weight rule the package has for the uniform weight on a d-dimensional cube.

Usage

lepos(d, q = NULL, p = NULL, normalize = TRUE, pragmatic = FALSE)

Arguments

Details

lepos(d, q) follows one-dimensional Gauss-Legendre functions such as statmod::gauss.quad(q, kind = "legendre"): q is the number of nodes of the one-dimensional Gauss rule, and lepos(d, q) returns a rule with the same exactness in d dimensions, degree ⁠p = 2q - 1⁠, that replaces the q^d-node product grid; for d = 1 it is the q-node Gauss-Legendre rule itself. Alternatively give the degree p.

The details given for ghpos() apply here too. All nodes of every stored Le rule lie inside the cube.

Value

A list with components nodes (⁠n x d⁠ matrix, one node per row) and weights (n strictly positive weights).

See Also

ghpos(), nnodes(), ruleinfo(), available_rules()

Examples

r <- lepos(2, 5)                                   # degree 9: 17 nodes on [0,1]^2 instead of 25
sum(r$weights * r$nodes[, 1]^3 * r$nodes[, 2]^2)   # 1/4 * 1/3

r <- lepos(2, 5, normalize = FALSE)                # the same rule on [-1,1]^2, weights sum to 4
r <- lepos(6, 6, pragmatic = TRUE)                 # d = 6 from a product of stored rules

Number of nodes of a rule, without building it

Description

The number of nodes of the rule that ghpos() (family = "gh") or lepos() (family = "le") returns for the same arguments. Signals the same error when there is no rule.

Usage

nnodes(family, d, q = NULL, p = NULL, pragmatic = FALSE)

Arguments

Value

The node count, as a double (tensor products in high dimensions exceed the integer range).

Examples

nnodes("gh", 5, 7)
nnodes("gh", 7, 5, pragmatic = TRUE)
5^7                                    # the product grid

Describe a rule: its factors and their origins

Description

Describes the rule that ghpos() or lepos() returns for the same arguments.

Usage

ruleinfo(family, d, q = NULL, p = NULL, pragmatic = FALSE)

Arguments

Details

Use it to find out whom to cite: origin names the published source of every rule that is not the package author's own.

Value

A data frame with one row per tensor factor, in the order of the columns of nodes, and the columns of available_rules(). A one-dimensional Gauss factor has origin = "Gauss" and NA in the columns that do not apply. A request answered by a single stored rule gives one row.

Examples

ruleinfo("le", 3, p = 41)$origin
ruleinfo("gh", 7, 5, pragmatic = TRUE)[, c("d", "p", "n", "origin")]