Univariate Gauss-Hermite quadrature rule

Create a univariate Gauss-Hermite quadrature rule.

GHrule(ord, asMatrix = TRUE)

Arguments

ord

scalar integer between 1 and 100 - the order, or number of nodes and weights, in the rule. When the function being multiplied by the standard normal density is a polynomial of order \(2k-1\) the rule of order \(k\) integrates the product exactly.

asMatrix

logical scalar - should the result be returned as a matrix. If FALSE a data frame is returned. Defaults to TRUE.

Value

a matrix (or data frame, is asMatrix is false) with ord rows and three columns which are z the node positions, w the weights and ldnorm, the logarithm of the normal density evaluated at the nodes.

Details

This version of Gauss-Hermite quadrature provides the node positions and weights for a scalar integral of a function multiplied by the standard normal density.

Originally based on package SparseGrid's hidden GQN(), then on fastGHQuad's gaussHermiteData(.).

References

Qing Liu and Donald A. Pierce (1994). A Note on Gauss-Hermite Quadrature. Biometrika 81(3), 624–629. doi:10.2307/2337136

See also

a different interface is available via GQdk().

Examples

(r5  <- GHrule( 5, asMatrix=FALSE))
#>               z          w     ldnorm
#> 1 -2.856970e+00 0.01125741 -5.0000774
#> 2 -1.355626e+00 0.22207592 -1.8377997
#> 3  3.865099e-17 0.53333333 -0.9189385
#> 4  1.355626e+00 0.22207592 -1.8377997
#> 5  2.856970e+00 0.01125741 -5.0000774
(r12 <- GHrule(12, asMatrix=FALSE))
#>            z            w     ldnorm
#> 1  -5.500902 1.499927e-07 -16.048898
#> 2  -4.271826 4.837185e-05 -10.043187
#> 3  -3.223710 2.203381e-03  -6.115091
#> 4  -2.259464 2.911669e-02  -3.471528
#> 5  -1.340375 1.469670e-01  -1.817241
#> 6  -0.444403 3.216644e-01  -1.017686
#> 7   0.444403 3.216644e-01  -1.017686
#> 8   1.340375 1.469670e-01  -1.817241
#> 9   2.259464 2.911669e-02  -3.471528
#> 10  3.223710 2.203381e-03  -6.115091
#> 11  4.271826 4.837185e-05 -10.043187
#> 12  5.500902 1.499927e-07 -16.048898

## second, fourth, sixth, eighth and tenth central moments of the
## standard Gaussian N(0,1) density:
ps <- seq(2, 10, by = 2)
cbind(p = ps, "E[X^p]" = with(r5,  sapply(ps, function(p) sum(w * z^p)))) # p=10 is wrong for 5-rule
#>       p E[X^p]
#> [1,]  2      1
#> [2,]  4      3
#> [3,]  6     15
#> [4,]  8    105
#> [5,] 10    825
p <- 1:15
GQ12 <- with(r12, sapply(p, function(p) sum(w * z^p)))
cbind(p = p, "E[X^p]" = zapsmall(GQ12))
#>        p E[X^p]
#>  [1,]  1      0
#>  [2,]  2      1
#>  [3,]  3      0
#>  [4,]  4      3
#>  [5,]  5      0
#>  [6,]  6     15
#>  [7,]  7      0
#>  [8,]  8    105
#>  [9,]  9      0
#> [10,] 10    945
#> [11,] 11      0
#> [12,] 12  10395
#> [13,] 13      0
#> [14,] 14 135135
#> [15,] 15      0
## standard R numerical integration can do it too:
intL <- lapply(p, function(p) integrate(function(x) x^p * dnorm(x),
                                        -Inf, Inf, rel.tol=1e-11))
integR <- sapply(intL, `[[`, "value")
cbind(p, "E[X^p]" = integR)# no zapsmall() needed here
#>        p E[X^p]
#>  [1,]  1      0
#>  [2,]  2      1
#>  [3,]  3      0
#>  [4,]  4      3
#>  [5,]  5      0
#>  [6,]  6     15
#>  [7,]  7      0
#>  [8,]  8    105
#>  [9,]  9      0
#> [10,] 10    945
#> [11,] 11      0
#> [12,] 12  10395
#> [13,] 13      0
#> [14,] 14 135135
#> [15,] 15      0
all.equal(GQ12, integR, tol=0)# => shows small difference
#> [1] "Mean relative difference: 5.910896e-15"
stopifnot(all.equal(GQ12, integR, tol = 1e-10))
(xactMom <- cumprod(seq(1,13, by=2)))
#> [1]      1      3     15    105    945  10395 135135
stopifnot(all.equal(xactMom, GQ12[2*(1:7)], tol=1e-14))
## mean relative errors :
mean(abs(GQ12  [2*(1:7)] / xactMom - 1)) # 3.17e-16
#> [1] 3.172066e-16
mean(abs(integR[2*(1:7)] / xactMom - 1)) # 9.52e-17 {even better}
#> [1] 9.516197e-17