[R] Integral of PDF

[William Dunlap]

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.
>