[R] Integral of PDF

[Ravi Varadhan]

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


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Hans W Borchers
Sent: Friday, December 03, 2010 3:30 AM
To: r-help at 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: 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.

______________________________________________
R-help at 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.