Univariate Gauss-Hermite quadrature rule

Create a univariate Gauss-Hermite quadrature rule.

Usage

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

Liu Q, Pierce DA (1994). “A note on Gauss—Hermite quadrature.” Biometrika, 81(3), 624–629. doi:10.2307/2337136 .

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