[R] Integral of PDF

[Ravi Varadhan]

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


-----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: Thursday, December 02, 2010 5:16 PM
To: r-help at r-project.org
Subject: Re: [R] Integral of PDF


You can dive into the thread "puzzle with integrate over infinite range"
from September this year. The short answer appears to be: Increase the
error tolerance.

    integrate(function(x) dnorm(x, 500,50), -Inf, Inf,
              subdivisions=500, rel.tol=1e-11)
    # 1 with absolute error < 1.1e-12

Hans Werner



Doran, Harold wrote:
> 
> The integral of any probability density from -Inf to Inf should equal 1,
> correct? I don't understand last result below.
> 
>> integrate(function(x) dnorm(x, 0,1), -Inf, Inf)
> 1 with absolute error < 9.4e-05
> 
>> integrate(function(x) dnorm(x, 100,10), -Inf, Inf)
> 1 with absolute error < 0.00012
> 
>> integrate(function(x) dnorm(x, 500,50), -Inf, Inf)
> 8.410947e-11 with absolute error < 1.6e-10
> 
>> all.equal(integrate(function(x) dnorm(x, 500,50), -Inf, Inf)$value, 0)
> [1] TRUE
> 
>> sessionInfo()
> R version 2.10.1 (2009-12-14)
> i386-pc-mingw32
> 
> locale:
> [1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United
> States.1252
> [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
> [5] LC_TIME=English_United States.1252
> 
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
> 
> other attached packages:
> [1] statmod_1.4.6      mlmRev_0.99875-1   lme4_0.999375-35  
> Matrix_0.999375-33 lattice_0.17-26
> 
> loaded via a namespace (and not attached):
> [1] grid_2.10.1   nlme_3.1-96   stats4_2.10.1 tools_2.10.1
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> 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.
> 
> 
-- 
View this message in context:
http://r.789695.n4.nabble.com/Integral-of-PDF-tp3070243p3070315.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.