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.