[R] Gauss-Laguerre using statmod
Doran, Harold
HDoran at air.org
Fri Aug 7 15:33:42 CEST 2009
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
More information about the R-help
mailing list