# Integral of PDF

11 messages
Open this post in threaded view
|

## Integral of PDF

 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]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: 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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Integral of PDF

 In reply to this post by Doran, Harold To really understaand it you will have to look at the fortran code underlying integrate.  I tracked it back through a couple of layers (dqagi, dqagie, ...  just use google, these are old netlib subroutines) then decided I ought to get back to grading papers :-) It looks like the integral is split into the sum of two integrals, (-Inf,0] and [0,Inf).   > integrate(function(x) dnorm(x, 350,50), 0, Inf) 1 with absolute error < 1.5e-05 > integrate(function(x) dnorm(x, 400,50), 0, Inf) 1.068444e-05 with absolute error < 2.1e-05 > integrate(function(x) dnorm(x, 500,50), 0, Inf) 8.410947e-11 with absolute error < 1.6e-10 > integrate(function(x) dnorm(x, 500,50), 0, 1000) 1 with absolute error < 7.4e-05 The integral from 0 to Inf is the lim_{x -> Inf} of the integral over [0,x].  I suspect that the algorithm is picking an interval [0,x], evaluating the integral over that interval, and then performing some test to decide whether to expand the interval.  When the initial interval doesn't contain much, expanding a little won't add enough to trigger another iteration.   albyn On Thu, Dec 02, 2010 at 04:21:34PM -0500, 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]] > > ______________________________________________ > [hidden email] 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. > -- Albyn Jones Reed College [hidden email] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Integral of PDF

 In reply to this post by Doran, Harold You can use trace() to see what is happening > trace(dnorm, quote(print(round(x)))) > integrate(function(x) dnorm(x, 500,50), -Inf, Inf) Tracing dnorm(x, 500, 50) on entry  [1]   1 233   0  38   0  14   0   7   0   4   0   2   0   2   1 Tracing dnorm(x, 500, 50) on entry  [1]   -1 -233    0  -38    0  -14    0   -7    0   -4    0   -2    0 -2   -1 Tracing dnorm(x, 500, 50) on entry  [1]   3 467   1  78   1  29   1  14   1   9   2   6   2   4   2 Tracing dnorm(x, 500, 50) on entry  [1]   -3 -467   -1  -78   -1  -29   -1  -14   -1   -9   -2   -6   -2 -4   -2 Tracing dnorm(x, 500, 50) on entry  [1] 0 1 0 1 0 1 0 1 0 1 0 1 0 0 0 Tracing dnorm(x, 500, 50) on entry  [1]  0 -1  0 -1  0 -1  0 -1  0 -1  0 -1  0  0  0 Tracing dnorm(x, 500, 50) on entry  [1]   7 935   3 156   3  58   3  30   4  18   4  12   5   9   6 Tracing dnorm(x, 500, 50) on entry  [1]   -7 -935   -3 -156   -3  -58   -3  -30   -4  -18   -4  -12   -5 -9   -6 Tracing dnorm(x, 500, 50) on entry  [1] 2 3 1 3 1 3 1 3 1 2 1 2 1 2 1 Tracing dnorm(x, 500, 50) on entry  [1] -2 -3 -1 -3 -1 -3 -1 -3 -1 -2 -1 -2 -1 -2 -1 8.410947e-11 with absolute error < 1.6e-10 You are asking integrate to find the relatively tiny portion of the the floating point real line (from c. -10^300 to 10^300) on which dnorm(x) is bigger than c. 10^-300.  It found a few points where it was that big, but apparently not enough to home on the answer.  You need to give it better hints abouts dnorm's support region.  E.g., > integrate(function(x) dnorm(x, 500,50), -100, 900) Tracing dnorm(x, 500, 50) on entry  [1] 400 -87 887 -33 833  60 740 183 617 326 474 -98 898 -65 865  10 790 119 681 [20] 253 547 Tracing dnorm(x, 500, 50) on entry  [1] 150 -93 393 -66 366 -20 320  42 258 113 187 -99 399 -83 383 -45 345 9 291 [20]  76 224 Tracing dnorm(x, 500, 50) on entry  [1] 650 407 893 434 866 480 820 542 758 613 687 401 899 417 883 455 845 509 791 [20] 576 724 Tracing dnorm(x, 500, 50) on entry  [1] 525 403 647 417 633 440 610 471 579 506 544 401 649 409 641 427 623 455 595 [20] 488 562 Tracing dnorm(x, 500, 50) on entry  [1] 775 653 897 667 883 690 860 721 829 756 794 651 899 659 891 677 873 705 845 [20] 738 812 1 with absolute error < 4.0e-07 Making the limits +-10^4 works ok also, but +-Inf is apparently asking for too much. Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com   > -----Original Message----- > From: [hidden email] > [mailto:[hidden email]] On Behalf Of Doran, Harold > Sent: Thursday, December 02, 2010 1:22 PM > To: [hidden email] > Subject: [R] Integral of PDF > > 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]] > > ______________________________________________ > [hidden email] 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. > ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Integral of PDF

Open this post in threaded view
|

## Re: Integral of PDF

Open this post in threaded view
|

## Re: Integral of PDF

 In reply to this post by Ravi Varadhan On Thu, Dec 02, 2010 at 06:23:45PM -0500, Ravi Varadhan wrote: > 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. Translation:  think before you compute! albyn ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Integral of PDF

Open this post in threaded view
|

## Re: Integral of PDF

 In reply to this post by Doran, Harold Hi, You need to be careful with -Inf and Inf in R, I suppose they are some large numbers arbitrarily chosen by R, but not infinite quantities. So the function integrate can return errors if the function to integrate doesn't have negligible values beyond those large numbers. I ran the following to convince myself: > integrate(function(x){return(ifelse(x<500,0,1/(x^2)))},-Inf,Inf) 0 with absolute error < 0 > integrate(function(x){return(ifelse(x<500,0,1/(x^2)))},500,Inf) 0.001999999 with absolute error < 0.000011 The second result shouldn't be greater that the first one but... And to solve your problem: > integrate(function(x) dnorm(x, 500,50), 500-1000, 500+1000) 1 with absolute error < 0.000074 I centred the interval of integration on 500 instead of on 0 if you integrate between -Inf and Inf. HTH, Samuel -----Original Message----- From: [hidden email] [mailto:[hidden email]] On Behalf Of Doran, Harold Sent: 02 December 2010 21:22 To: [hidden email] Subject: [R] Integral of PDF 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]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code. __________ Information from ESET NOD32 Antivirus, version of virus signature database 5668 (20101202) __________ The message was checked by ESET NOD32 Antivirus. http://www.eset.com__________ Information from ESET NOD32 Antivirus, version of virus signature database 5669 (20101203) __________ The message was checked by ESET NOD32 Antivirus. http://www.eset.com______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.