[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