[R] Integrate to 1? (gauss.quad)

[Doran, Harold]

Thank you, Doug. I am still missing something here. Should this simply be sum(f(x_i) * w_i) where x_i is node i and w_i is the weight at node i? So, my function f(x) = (1/(s*sqrt(2*pi)))  * exp(-((qq$nodes-mu)^2/(2*s^2))) would then be multiplied only by the weights, qq$weights

Here is some reproducible code:

> library(statmod)
> Q <- 5
> mu <- 0; s <- 1
> qq <- gauss.quad(Q, kind = 'hermite')
> sum((1/(s*sqrt(2*pi)))  * exp(-((qq$nodes-mu)^2/(2*s^2))) * qq$weights)
[1] 0.5775682

> sum(qq$weights)
[1] 1.772454

> (1/(s*sqrt(2*pi)))  * exp(-((qq$nodes-mu)^2/(2*s^2)))
[1] 0.05184442 0.25198918 0.39894228 0.25198918 0.05184442

> dnorm(qq$nodes)
[1] 0.05184442 0.25198918 0.39894228 0.25198918 0.05184442

> -----Original Message-----
> From: dmbates at gmail.com [mailto:dmbates at gmail.com] On Behalf Of Douglas Bates
> Sent: Sunday, November 14, 2010 2:28 PM
> To: Doran, Harold
> Cc: r-help at r-project.org
> Subject: Re: [R] Integrate to 1? (gauss.quad)
> 
> I don't know about the statmod package and the gauss.quad function but
> generally the definition of Gauss-Hermite quadrature is with respect
> to the function that is multiplied by exp(-x^2) in the integrand.  So
> your example would reduce to summing the weights.
> 
> On Sun, Nov 14, 2010 at 11:18 AM, Doran, Harold <HDoran at air.org> wrote:
> > Does anyone see why my code does not integrate to 1?
> >
> > library(statmod)
> >
> > mu <- 0
> > s <- 1
> > Q <- 5
> > qq <- gauss.quad(Q, kind='hermite')
> > sum((1/(s*sqrt(2*pi)))  * exp(-((qq$nodes-mu)^2/(2*s^2))) * qq$weights)
> >
> > ### This does what's it is supposed to
> > myNorm <- function(theta) (1/(s*sqrt(2*pi)))  * exp(-((theta-mu)^2/(2*s^2)))
> > integrate(myNorm, -Inf, Inf)
> > ______________________________________________
> > R-help at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >