[R] Gauss-Laguerre using statmod

[Doran, Harold]

I believe this may be more related to analysis than it is to R, per se.
Suppose I have the following function that I wish to integrate:

ff <- function(x) pnorm((x - m)/sigma) * dnorm(x, observed, sigma)

Then, given the parameters:

mu <- 300
sigma <- 50
m <- 250
target <- 200
sigma_i <- 50

I can use the function integrate as:

> integrate(ff, lower= -Inf, upper=target)
0.002169851 with absolute error < 4.4e-05

I would like to also use Gauss-Laguerre methods to also integrate this
function. In doing so, I believe the only change of variable needed when
integrating from -Inf to target is x = target - y_i where y_i is node i.
As such, I can implement the following:

library(statmod)

falsePos <- function(target, m, mu, sigma, sigma_i, Q = 30){
   gq <- gauss.quad(Q, kind="laguerre")
   nodes <- gq$nodes
   whts  <- gq$weights
   y <- pnorm((target - nodes - m)/sigma_i) * dnorm(target - nodes, mu,
sigma)
   sum(y * exp(nodes)* whts)
}

> falsePos(target = 200, m = 250, mu = 300, sigma = 50, sigma_i = 50)
[1] 0.002169317

Which yields the same value as the internal R function. Now suppose I
want to integrate in the opposite direction going from target to Inf
using the same parameters as previously used. Again, the internal
integrate function yields:

> integrate(ff, lower=target, upper=Inf)
0.7580801 with absolute error < 7.2e-05

Now, my understanding of the change of variable needed for
Gauss-Laguerre in this instance is simple, x = y_i + target. As such,
the integration should be

truePos <- function(target, m, mu, sigma, sigma_i, Q = 30){
   gq <- gauss.quad(Q, kind="laguerre")
   nodes <- gq$nodes
   whts  <- gq$weights
   y <- pnorm((nodes + target - m)/sigma_i) * dnorm(nodes + target, mu,
sigma)
   sum(y * exp(nodes)* whts)
}

> truePos(target = 200, m = 250, mu = 300, sigma = 50, sigma_i = 50)
[1] 0.2533494

Clearly, there is not a match in this instance. Looking at the density
we can see that R's internal function is correct:

plot(ff, 0, 500)

I've used Gauss-Laguerre in the past with this same change of variable
and obtained correct results for the integral. However, I seem to have
an error in this situation that I can't seem to identify. Can anyone
point out whether I have an error in the way I am approaching the
problem mathematically or is there a programming error I seem to be
missing?

Many thanks for your help.

Harold

> sessionInfo()
R version 2.9.0 (2009-04-17) 
i386-pc-mingw32 

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base


other attached packages:
[1] MiscPsycho_1.4     statmod_1.3.8      xtable_1.5-5
lme4_0.999375-28   Matrix_0.999375-25
[6] VAM_0.8-5          lattice_0.17-22   

loaded via a namespace (and not attached):
[1] grid_2.9.0  tools_2.9.0