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

Enrico Schumann enricoschumann at yahoo.de
Tue Nov 16 09:29:30 CET 2010

```
gauss--hermite weights the function f(x) to be integrated by exp(-nodes^2),
ie, it uses the 'trick' that exp(nodes^2) * exp(-nodes^2) = 1. so your sum
should be exp(qq\$nodes^2) * exp(-qq\$nodes^2) * f(qq\$nodes^2) =
exp(qq\$nodes^2) * qq\$weights * f(qq\$nodes). hence

sum(exp(qq\$nodes^2) * (1/(s*sqrt(2*pi))) * exp(-((qq\$nodes-mu)^2/(2*s^2))) *
qq\$weights)

or

sum(exp(qq\$nodes^2) * qq\$weights * myNorm(qq\$nodes))

should work.

hth,
enrico

> -----Ursprüngliche Nachricht-----
> Von: r-help-bounces at r-project.org
> [mailto:r-help-bounces at r-project.org] Im Auftrag von Doran, Harold
> Gesendet: Montag, 15. November 2010 18:11
> An: Douglas Bates
> Cc: r-help at r-project.org
> Betreff: Re: [R] Integrate to 1? (gauss.quad)
>
> 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)
>  0.5775682
>
> > sum(qq\$weights)
>  1.772454
>
> > (1/(s*sqrt(2*pi)))  * exp(-((qq\$nodes-mu)^2/(2*s^2)))
>  0.05184442 0.25198918 0.39894228 0.25198918 0.05184442
>
> > dnorm(qq\$nodes)
>  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
> > > http://www.R-project.org/posting-guide.html
> > > and provide commented, minimal, self-contained, reproducible code.
> > >
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help