Integral of PDF
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.