proposed pbirthday fix

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

proposed pbirthday fix

knoblauch
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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel