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

[Enrico Schumann]

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)
> [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.
> > >
> 
> ______________________________________________
> 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.
> No virus found in this incoming message.
> Checked by AVG - www.avg.com
> Version: 9.0.869 / Virus Database: 271.1.1/3256 - Release 
> Date: 11/14/10 08:34:00
>