[R] Some Help Needed

[peter dalgaard]

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