R devel - Jan 2006 - proposed pbirthday fix

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

Actually, since NaN's are also detected in na.action operations, a 
simpler fix
might just be to use the na.rm = TRUE option of min

upper <- min(n^k/(c^(k - 1)), 1, na.rm = TRUE)

> 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
> }
	[[alternative text/enriched version deleted]]