On Apr 24, 2012, at 05:52 , li li wrote:
> Dear all,
> I need to do some calculation where the code used are below. I get
> error message when I choose k to be large, say greater than 25.
> The error message is
> "Error in integrate(temp, lower = 0, upper = 1, k, x, rho, m) :
> the integral is probably divergent".
> Can anyone give some help on resolving this. Thanks.
For the case at hand, adding stop.on.error=FALSE to the integrate() call
probably works, but you need to be somewhat suspicious of the results and check
carefully.
As a general matter, numerical integrators are confused by near-discontinuous
behaviour. If you plot your integrand at one of the values where the integration
fails:
curve(temp(y, 30, x=.06, rho, m), xname="y")
you'll see that it is 1 for y=0 then drops rapidly to zero, but not so
quickly that it is completely flat for y>0.
These things are tricky to work around, but it usually involves some sort of
hinting to tell the integrator where "things are happening", which in
turn may require mathematical analysis (pen-and-paper style) of the integrand,
or crude guesswork. E.g., you could eyeball it and split into two integrals:
> integrate(temp, lower = 0, upper=.1, k=30, x=.06, rho,m)
0.000809987 with absolute error < 1.6e-05> integrate(temp, lower = 0.1, upper=1, k=30, x=.06, rho,m)
1.444079e-09 with absolute error < 1.1e-11
(Notice that the result is well below your target of 0.05, which is why I said
that you could probably get away with using stop.on.error -- all you really need
is that the integral is too small for x=.06)
You could also utilze the fact that the integrand is a rather nicer function of
log(y)
> curve(temp(y, 30, x=.06, rho, m), xname="y", log="x",
from=1e-6 )
and consider integration by substitution:
> temp2 <- function(u,...) {y <- exp(u); y*temp(y,...)}
> integrate(temp2, lower = -Inf, upper=0, k=30, x=.06, rho, m)
0.0008099758 with absolute error < 3.7e-05
> Hannah
>
> m <- 100
>
> alpha <- 0.05
>
> rho <- 0.1
>
>
>
> F0 <- function(y,x,rho){
>
> pnorm((qnorm(x, lower.tail = F)-sqrt(rho)*qnorm(y, lower.tail = F))/sqrt(1-
> rho), lower.tail = F)
>
> }
>
> temp <- function(y,k,x,rho,m) {
>
> t <- F0(y=y, x=x, rho=rho)
>
> pbinom(q=k, size=m, prob=t, lower.tail = F, log.p = FALSE)
>
> }
>
>
>
> est <- function(k,x,rho,m) {
>
> integrate(temp, lower = 0, upper=1, k,x, rho,m)$value-alpha}
>
>
>
> estf <- function(x){est(x, k=30, rho=rho, m=m)}
>
> thre <- uniroot(estf, c(0,1), tol=1/10^12)$root
>
> [[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.
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com