Gaussian Quadrature with Probability Distributions
Calculate nodes and weights for Gaussian quadrature in terms of probability distributions.
gauss.quad.prob(n, dist = "uniform", l = 0, u = 1,
mu = 0, sigma = 1, alpha = 1, beta = 1)n |
number of nodes and weights |
dist |
distribution that Gaussian quadrature is based on, one of |
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 probability interpretation is explained by Smyth (1998).
For details on the underlying quadrature rules, see gauss.quad.
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, using Netlib Fortran code http://www.netlib.org/go/gaussq.f, and including corrections suggested by Spencer Graves
Smyth, G. K. (1998). Polynomial approximation. In: Encyclopedia of Biostatistics, P. Armitage and T. Colton (eds.), Wiley, London, pages 3425-3429. http://www.statsci.org/smyth/pubs/PolyApprox-Preprint.pdf
# the 4th moment of the standard normal is 3 out <- gauss.quad.prob(10,"normal") sum(out$weights * out$nodes^4) # the expected value of log(X) where X is gamma is digamma(alpha) out <- gauss.quad.prob(32,"gamma",alpha=5) sum(out$weights * log(out$nodes))
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.