[R] integrate: Don't do this?

Göran Broström gb at stat.umu.se
Wed Feb 6 13:36:23 CET 2013


To answer my own question: I thought for a while that centering is the 
solution to the problem, but it is not:

 > integrate(dnorm, 0, Inf, mean = 25)
3.187474e-05 with absolute error < 5.9e-05
 > integrate(dnorm, -25, Inf, mean = 0)
3.187474e-05 with absolute error < 5.9e-05

Here "Don't do this" is the wrong advise. Better:

 > integrate(dnorm, -25, 10, mean = 0)
1 with absolute error < 1.5e-06

So, in my slightly more complicated example below,
the integration limits should be

max(0, nv - 10), nv + 10 (or so)

i.e., a finite interval of suitable length. Maybe the advice in the 
documentation about "integrating over infinite intervals" needs some 
modification?

Göran

On 02/05/2013 07:42 PM, Göran Broström wrote:
> When I run the following function
>
> HQ2 <- function(n) {
>        nv <- 6 * sqrt(n)
>        fcn <- function(z) {
>            pchisq(z^2 / 36, n - 1) * dnorm(nv - z)
>        }
>        ## I want the integral from 0 to infinity:
>        f.Inf <- integrate(fcn, 0, Inf)
>        ## Doc: "Don't do this":
>        f.100 <- integrate(fcn, 0, 100)
>        cbind(f.Inf, f.100)
> }
>
> I get, for n = 9 and n = 10:
>
>    > HQ2(9)
>                 f.Inf        f.100
> value        0.6531951    0.6531951
> abs.error    5.509777e-05 2.931259e-07
> subdivisions 10           5
> message      "OK"         "OK"
> call         Expression   Expression
>
> which indicates that the "Don't do this" works. With n = 10,
>
>    > HQ2(10)
>                 f.Inf        f.100
> value        3.292159e-05 0.6451915
> abs.error    6.121266e-05 4.176584e-08
> subdivisions 2            5
> message      "OK"         "OK"
> call         Expression   Expression
>
> the "recommended" way fails badly, but "Don't do this" works! Any
> recommendations? (Or explanations?)
>
> Göran
>
>    > sessionInfo()
> R version 2.15.2 Patched (2012-12-07 r61250)
> Platform: x86_64-unknown-linux-gnu (64-bit)
>
> locale:
>     [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>     [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>     [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
>     [7] LC_PAPER=C                 LC_NAME=C
>     [9] LC_ADDRESS=C               LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> loaded via a namespace (and not attached):
> [1] tools_2.15.2
>    >
>



More information about the R-help mailing list