################################################################################ # Examples involving Gauss-Hermite quadrature # ################################################################################ ################################################################################ #Example on p. 270 of McCulloch and Searle for Gauss-Hermite library(pracma) save.GH <- gaussHermite(n = 3) save.GH$x save.GH$w x <- save.GH$x w <- save.GH$w sum((1 + x^2) * w) sum((1 + x^6) * w) save.GH <- gaussHermite(n = 4) x <- save.GH$x w <- save.GH$w sum((1 + x^6) * w) fx1 <- function(x) { (1 + x^2) * exp(-x^2) } integrate(f = fx1, lower = -Inf, upper = Inf) fx2 <- function(x) { (1 + x^6) * exp(-x^2) } integrate(f = fx2, lower = -Inf, upper = Inf) ################################################################################ # X ~ N(0, sigma^2) ###################################### # Show integrate to 1 #Original function fx1 <- function(x, sigma = 10) { 1/(sqrt(2*pi) * sigma) * exp(-x^2 /(2*sigma^2)) } fx1(x = 0) dnorm(x = 0, mean = 0, sd = 10) integrate(f = fx1, lower = -Inf, upper = Inf) #u substitution fu1 <- function(u) { 1/sqrt(pi) * exp(-u^2) } integrate(f = fu1, lower = -Inf, upper = Inf) #Gauss-Hermite save.GH <- gaussHermite(n = 2) x <- save.GH$x w <- save.GH$w sum(1/sqrt(pi) * w) ###################################### # E(X) fx2 <- function(x, sigma = 10) { x * 1/(sqrt(2*pi) * sigma) * exp(-x^2 /(2*sigma^2)) } integrate(f = fx2, lower = -Inf, upper = Inf) #u substitution fu2 <- function(u, sigma = 10) { (sqrt(2) * sigma * u)/sqrt(pi) * exp(-u^2) } integrate(f = fu2, lower = -Inf, upper = Inf) #Gauss-Hermite save.GH <- gaussHermite(n = 2) u <- save.GH$x w <- save.GH$w sigma <- 10 sum((sqrt(2) * sigma * u)/sqrt(pi) * w) ###################################### # E(X^2) fx3 <- function(x, sigma = 10) { x^2 * 1/(sqrt(2*pi) * sigma) * exp(-x^2 /(2*sigma^2)) } integrate(f = fx3, lower = -Inf, upper = Inf) #u substitution fu3 <- function(u, sigma = 10) { (sqrt(2) * sigma * u)^2 / sqrt(pi) * exp(-u^2) } integrate(f = fu3, lower = -Inf, upper = Inf) #Gauss-Hermite save.GH <- gaussHermite(n = 2) u <- save.GH$x w <- save.GH$w sigma <- 10 sum((sqrt(2) * sigma * u)^2 / sqrt(pi) * w) # #