gauss.quad.prob {statmod}R Documentation

Gaussian Quadrature with Probability Distributions


Calculate nodes and weights for Gaussian quadrature in terms of probability distributions.




n number of nodes and weights
dist distribution that Gaussian quadrature is based on, one of "uniform", "normal", "beta" or "gamma"
l lower limit of uniform distribution
u upper limit of uniform distribution
mu mean of normal distribution
sigma standard deviation of normal distribution
alpha positive shape parameter for gamma distribution or first shape parameter for beta distribution
beta positive scale parameter for gamma distribution or second shape parameter for beta distribution


This is a rewriting and simplification of gauss.quad in terms of probability distributions.

The expected value of f(X) is approximated by sum(w*f(x)) where x is the vector of nodes and w is the vector of weights. The approximation is exact if f(x) is a polynomial of order no more than 2n-1. The possible choices for the distribution of X are as follows:

Uniform on (l,u).

Normal with mean mu and standard deviation sigma.

Beta with density x^(alpha-1)*(1-x)^(beta-1)/B(alpha,beta) on (0,1).

Gamma with density x^(alpha-1)*exp(-x/beta)/beta^alpha/gamma(alpha).


A list containing the components

nodes vector of values at which to evaluate the function
weights vector of weights to give the function values


Gordon Smyth


Golub, G. H., and Welsch, J. H. (1969). Calculation of Gaussian quadrature rules. Mathematics of Computation 23, 221-230.

Golub, G. H. (1973). Some modified matrix eigenvalue problems. Siam Review 15, 318-334.

Smyth, G. K. (1998). Polynomial approximation. In: Encyclopedia of Biostatistics, P. Armitage and T. Colton (eds.), Wiley, London, pp. 3425-3429.

Stroud and Secrest (1966). Gaussian Quadrature Formulas. Prentice- Hall, Englewood Cliffs, N.J.

See Also

gauss.quad, integrate


out <- gauss.quad.prob(10,"normal")
sum(out$weights * out$nodes^4)
#  the 4th moment of the standard normal is 3

out <- gauss.quad.prob(32,"gamma",alpha=5)
sum(out$weights * log(out$nodes))
#  the expected value of log(X) where X is gamma is digamma(alpha)

[Package statmod version 1.2.4 Index]