[R] Gauss-Laguerre using statmod

[Ravi Varadhan]

Harold,

I don't think there is any error in your code.  The problem is with using
Gauss-Laguerre quadrature for this "integrand".  

I changed your function `ff' slightly so that it admits a closed-form
integral (I took the means and sigmas to be the same):


ff2 <- function(x) pnorm(x, mu, sigma) * dnorm(x, mu, sigma)

# exact anti-derivative of this is 0.5 * (pnorm(x, mu, sigma))^2

options(digits=6)

xt <- seq(200, 400, by=20)
results <- matrix(NA, length(xt), 4)
colnames(results) <- c("target", "exact", "integrate", "laguerre")
i <- 0
for (target in xt) {
i <- i + 1
num1 <- integrate(ff2, lower= -Inf, upper=target)$val
exact <- 0.5 * pnorm(target, mu, sigma)^2
num2 <- falsePos(target = target, m = 300, mu = 300, sigma = 50, sigma_i =
50)
results[i, ] <- c(target, exact, num1, num2)
}

 results 

xt <- seq(200, 400, by=20)
results2 <- matrix(NA, length(xt), 4)
colnames(results2) <- c("target", "exact", "integrate", "laguerre")
i <- 0
for (target in xt) {
i <- i + 1
num1 <- integrate(ff2, lower= target, upper=Inf)$val
exact <- 0.5 - 0.5 * pnorm(target, mu, sigma)^2
num2 <- truePos(target = target, m = 300, mu = 300, sigma = 50, sigma_i =
50)
results2[i, ] <- c(target, exact, num1, num2)
}

 results2

These 2 experiments clearly show that you cannot accurately compute the
integral that you want using "Gauss-Laguerre" quadrature.  This problem
perists even after you increase the number of quadrature points from Q=30 to
Q=100.

It is also clear that `integrate' does a good job.

Ravi.

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Doran, Harold
Sent: Friday, August 07, 2009 9:34 AM
To: r-help at r-project.org
Subject: [R] Gauss-Laguerre using statmod

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

______________________________________________
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.