[Rd] pbinom gotcha (PR#1569)

[p.dalgaard@biostat.ku.dk]

This came up due to a question from Anders Hald:

Bernoulli calculated an approximation to the smallest n so that

P(0.58 <= x/n <= 0.62) >= 1000/1001

What is the exact value?

Now try

n <- 6350:6500
Pr <- function(n)pbinom(0.62*n,n,0.6) -  
  pbinom(0.58*n,n,0.6) + dbinom(0.58*n,n,0.6)
plot(n,Pr(n),type="b")
abline(h=1000/1001)
min(n[Pr(n)>1000/1001])

Next, try

Pr <- function(n)pbinom(0.62*n,n,0.6) - 
        pbinom(0.58*n + 1e-7,n,0.6) + dbinom(0.58*n,n,0.6)
plot(n,Pr(n),type="b")
abline(h=1000/1001)
min(n[Pr(n)>1000/1001])

At least on Linux, you get quite different results.
The reason is of course that 

> 0.58*6400-3712
[1] -4.547474e-13

and dbinom has a fuzz factor, but pbinom unceremoniously takes
floor(x). I suppose we ought to build the same fuzz into pbinom? 

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._