From knoblauch at lyon.inserm.fr  Mon Jan 23 01:18:54 2006
From: knoblauch at lyon.inserm.fr (Ken Knoblauch)
Date: Sun, 22 Jan 2006 23:18:54 -0100 (CET)
Subject: [Rd] proposed pbirthday fix
Message-ID: <51940.82.231.93.240.1137968334.squirrel@webmail.lyon.inserm.fr>

Recent news articles concerning an article from The Lancet with fabricated
data indicate
that in the sample containing some 900 or so patients, more than 200 had
the same
birthday.  I was curious and tried out the p and q birthday functions but
pbirthday
could not handle 250 coincidences with n = 1000.  The calculation of upper
prior
to using uniroot produces NaN,

upper<-min(n^k/(c^(k-1)),1)

I was able to get it to work by using logs, however, as in the following
version

function(n, classes = 365, coincident = 2){
    k <- coincident
    c <- classes
    if (coincident < 2) return(1)
    if (coincident > n) return(0)
    if (n > classes * (coincident - 1)) return(1)
    eps <- 1e-14
    if (qbirthday(1 - eps, classes, coincident) <= n)
        return(1 - eps)
    f <- function(p) qbirthday(p, c, k) - n
    lower <- 0
    upper <- min( exp( k * log(n) - (k-1) * log(c) ), 1 )
    nmin <- uniroot(f, c(lower, upper), tol = eps)
    nmin$root
}

Ken



[Rd] proposed pbirthday fix
Ken Knoblauch