Your function does NOT have a mode at zero. It is bimodal with 2 modes: -500 and 500. So, my approach still works. Zero is a stationary point where the gradient is zero, but it is not a mode (the second derivative at zero is not negative but it is zero).
Ravi. ------------------------------------------------------- Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edu -----Original Message----- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Hans W Borchers Sent: Friday, December 03, 2010 3:30 AM To: r-help@r-project.org Subject: Re: [R] Integral of PDF That does not remedy the situation in any case, take the following function fun <- function(x) dnorm(x, -500, 50) + dnorm(x, 500, 50) that has a 'mode' of 0 again. Interestingly, if I transform it by 1/x, to integrate I again have to reduce the error tolerance to at least 1e-10: f <- function(x) fun(1/x)/x^2 integrate(f, 0, 1) #-> 0, should be 1 while cubature::adaptIntegrate does not require to change the tolerance: require(cubature) adaptIntegrate(f, 0, 1) #-> 1 What I find *dangerous* about 'integrate' is not that it returns a wrong result, but that it --- w/o warning --- gives an error estimate that leads completely astray. Hans Werner But that only postpones the misery, Hans Werner! You can always make the mean large enough to get a wrong answer from `integrate'. integrate(function(x) dnorm(x, 700,50), -Inf, Inf, subdivisions=500, rel.tol=1e-11) integrate(function(x) dnorm(x, -700,50), -Inf, Inf, subdivisions=500, rel.tol=1e-11) The problem is that the quadrature algorithm is not smart enough to recognize that the center of mass is at `mean'. The QUADPACK algorithm (Gauss-Kronrod quadratutre) is adaptive, but it does not look to identify regions of high mass. Algorithms which can locate regions of highest density, such as those developed by Alan Genz, will do much better in problems like this. Genz and Kass (1997). J Computational Graphics and Statistics. There is a FORTRAN routine DCHURE that you might want to try for infinite regions. I don't think this has been ported to R (although I may be wrong). There may be other more recent ones. A simple solution is to locate the mode of the integrand, which should be quite easy to do, and then do a coordinate shift to that point and then integrate the mean-shifted integrand using `integrate'. Ravi. ------------------------------------------------------- Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edu That does not remedy the situation in any case, take the following function fun <- function(x) dnorm(x, -500, 50) + dnorm(x, 500, 50) that has a 'mode' of 0 again. Interestingly, if I transform it by 1/x, to integrate I again have to reduce the error tolerance to at least 1e-10: f <- function(x) fun(1/x)/x^2 integrate(f, 0, 1) #-> 0, should be 1 while cubature::adaptIntegrate does not require to change the tolerance: require(cubature) adaptIntegrate(f, 0, 1) #-> 1 What I find *dangerous* about 'integrate' is not that it returns a wrong result, but that it --- w/o warning --- gives an error estimate that leads completely astray. Hans Werner Ravi Varadhan wrote: > > But that only postpones the misery, Hans Werner! You can always make the > mean large enough to get a wrong answer from `integrate'. > > integrate(function(x) dnorm(x, 700,50), -Inf, Inf, > subdivisions=500, rel.tol=1e-11) > > integrate(function(x) dnorm(x, -700,50), -Inf, Inf, > subdivisions=500, rel.tol=1e-11) > > The problem is that the quadrature algorithm is not smart enough to > recognize that the center of mass is at `mean'. The QUADPACK algorithm > (Gauss-Kronrod quadratutre) is adaptive, but it does not look to identify > regions of high mass. Algorithms which can locate regions of highest > density, such as those developed by Alan Genz, will do much better in > problems like this. > > Genz and Kass (1997). J Computational Graphics and Statistics. > > There is a FORTRAN routine DCHURE that you might want to try for infinite > regions. I don't think this has been ported to R (although I may be > wrong). > There may be other more recent ones. > > A simple solution is to locate the mode of the integrand, which should be > quite easy to do, and then do a coordinate shift to that point and then > integrate the mean-shifted integrand using `integrate'. > > Ravi. > > ------------------------------------------------------- > Ravi Varadhan, Ph.D. > Assistant Professor, > Division of Geriatric Medicine and Gerontology School of Medicine Johns > Hopkins University > > Ph. (410) 502-2619 > email: rvarad...@jhmi.edu > -- View this message in context: http://r.789695.n4.nabble.com/Integral-of-PDF-tp3070243p3070768.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ 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. ______________________________________________ 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.