custom {custom.gauss.quad} | R Documentation |
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)
Arguments
which.f |
a list specifying the nonnegative integrable
weight function |
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.
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
.
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:
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