[R] Integral of PDF

[Hans W Borchers]

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.