From p.dalgaard@biostat.ku.dk  Fri May 17 17:30:11 2002
From: p.dalgaard@biostat.ku.dk (p.dalgaard@biostat.ku.dk)
Date: Fri, 17 May 2002 18:30:11 +0200 (MET DST)
Subject: [Rd] pbinom gotcha (PR#1569)
Message-ID: <200205171630.SAA21501@pubhealth.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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._


[Rd] pbinom gotcha (PR#1569)
p.dalgaard@biostat.ku.dk