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
>
______________________________________________
R-help@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.
