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]]

______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Reply via email to