[Rd] proposed pbirthday fix

[Martin Maechler]

>>>>> "ken" == ken knoblauch <knoblauch at lyon.inserm.fr>
>>>>>     on Mon, 23 Jan 2006 09:43:28 +0100 writes:

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

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

Well, I liked your first fix better -- thank you for it! --
since it's always good practice to formulate such as to avoid
overflow when possible. 
All things considered, I think I'd go for

   upper <- min( exp(k * log(n) - (k-1) * log(c)), 1, na.rm = TRUE)

Martin 

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

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

    Ken> I was able to get it to work by using logs, however, as
    Ken> 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
    >> }