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