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-boun...@r-project.org > [mailto:r-help-boun...@r-project.org] Im Auftrag von Doran, Harold > Gesendet: Montag, 15. November 2010 18:11 > An: Douglas Bates > Cc: r-help@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: dmba...@gmail.com [mailto:dmba...@gmail.com] On Behalf Of > > Douglas Bates > > Sent: Sunday, November 14, 2010 2:28 PM > > To: Doran, Harold > > Cc: r-help@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 > <hdo...@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@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@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 > ______________________________________________ R-help@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.