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: rvaradhan at 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: rvaradhan at 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.