Integral of PDF
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: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Doran, Harold Sent: Thursday, December 02, 2010 1:22 PM To: r-help at r-project.org 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]]
______________________________________________ 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.