custom.gauss.quad (version 1.0.0)

custom: Custom-Made Gauss Quadrature Nodes and Weights

Description

For the nonnegative weight function specified by which.f and given number n of nodes, the function custom computes the Gauss quadrature nodes asNumeric.nodes and corresponding weights asNumeric.weights which are double precision vectors.

Usage

custom(which.f, n)

Value

A list with the following elements: which.f, n, nbits.vec, list.Gauss.nodes, list.Gauss.weights, mat.timings, max.abs.diffs.nodes, sum.abs.diffs.weights, L.nodes, L.weights, asNumeric.nodes, asNumeric.weights.

Arguments

which.f

a list specifying the nonnegative integrable weight function \(f\), with the following three components: (i) name (in the form of a character string), (ii) support specified by a 2-vector of the endpoints of the interval, (iii) parameter vector when \(f\) belongs to a family of weight functions and is specified by the value of this parameter vector (if \(f\) is already fully specified then the parameter vector is set to NULL)

n

number of Gauss quadrature nodes

Details

Suppose that we wish to evaluate $$\int_{-\infty}^{\infty} g(x) f(x) dx,$$ where \(f\) is a specified nonnegative integrable weight function. The Gauss quadrature approximation to this integral has the form $$\sum_{i=1}^n \lambda_i \, g(\tau_i),$$ where \(\tau_1, \dots, \tau_n\) are called the nodes and \(\lambda_1, \dots, \lambda_n\) are called the corresponding weights. This approximation is exact whenever \(g\) is a polynomial of degree less than or equal to \(2n - 1\).

If \(f\) takes a form that leads to Gauss quadrature rules with nodes that are the roots of classical orthogonal polynomials of a continuous variable then these rules (such as Gauss Legendre, Gauss Hermite and Gauss Laguerre) are readily accessible to statisticians via R packages such as statmod. If, however, \(f\) does not take one of these particular forms then the Gauss quadrature nodes and weights need to be custom-made.

custom computes the Gauss quadrature nodes and weights, for given n, using a user-supplied formula for the \(r\)'th moment $$\int_{-\infty}^{\infty} x^r f(x) dx$$ for all nonnegative integers \(r\). This formula must be inserted by the user into the code for the function moments and must be able to be computed to an arbitrary number of bits (nbits) of precision using the R package Rmpfr.

To obtain some assurance that the Gauss quadrature nodes and weights are computed to sufficient precision for subsequent double precision computations in R, these nodes and weights are computed for the increasing numbers of bits of precision given in the 5-vector nbits.vec and the results compared. This comparison results in the criteria max.abs.diffs.nodes, sum.abs.diffs.weights, L.nodes and L.weights described in detail by Kabaila (2022). The execution times for various parts of the code are stored in mat.timings whose components are described by Kabaila (2022).

list.Gauss.nodes[[i]] and list.Gauss.weights[[i]] are the n-vectors of Gauss quadrature nodes and weights, respectively, computed using nbits.vec[i] bits of precision (i=1,...,5).

The most accurate approximations to the Gauss quadrature nodes and weights are list.Gauss.nodes[[5]] and list.Gauss.weights[[5]]. These are converted to double precision by applying the asNumeric function from the R package Rmpfr, resulting in asNumeric.nodes and asNumeric.weights, respectively.

References

Kabaila, P. (2022) Custom-made Gauss quadrature for statisticians. arXiv:2211.04729

See Also

moments

Examples

# Suppose that the weight function f is the probability density
# function of a random variable with the same probability
# distribution as R divided by the square root of m, where R has a
# chi distribution with m degrees of freedom.
# Also suppose that we wish to compute the Gauss quadrature nodes
# and weights, for number of nodes n = 5, when the parameter m = 160.
# The r th moment can be computed to an arbitrary number of bits of
# precision using the R package Rmpfr. We describe the weight function
# f using the following R commands:

m <- 160
which.f <- list(name="scaled.chi.pdf", support=c(0, Inf),
parameters=m)

# Here, "scaled.chi.pdf" is the name (a character string) that we
# have given to the weight function f. The R function moments includes
# the code needed to compute the r th moment to an arbitrary number
# of bits of precision using the R package Rmpfr.
# We compute the Gauss quadrature node and weight, for the toy example
# with number of nodes n=1, using the following R commands:

n <- 1
gauss.list <- custom(which.f, n)
old <- options(digits = 17)
gauss.list$asNumeric.nodes
gauss.list$asNumeric.weights
options(old)

# These commands take less than 1 second to run. The resulting
# of node and corresponding weight in double precision are:
# > gauss.list$asNumeric.nodes
# [1] 0.99843873022375829
# > gauss.list$asNumeric.weights
# [1] 1

# The computation times for number of nodes n=5, 17 and 33 are roughly
# 160 seconds, 31 minutes and 5 hours,respectively.
#
# We compute the Gauss quadrature nodes and weights, for number of
# nodes n=5, using the following R commands:
# \donttest{
n <- 5
gauss.list <- custom(which.f, n)
old <- options(digits = 17)
gauss.list$asNumeric.nodes
gauss.list$asNumeric.weights
options(old)
# }
# These commands take roughly 3 minutes to run. The resulting vectors
# of nodes and corresponding weights in double precision are:
# > gauss.list$asNumeric.nodes
# [1] 0.84746499810651410 0.92785998378868118 1.00262691212158761
# [4] 1.07930375924992528 1.16628363226782716
# > gauss.list$asNumeric.weights
# [1] 0.0144433732487188448 0.2483585328946608384 0.5305446123744097520
# [4] 0.1977278905956056654 0.0089255908866048821