Hi, pbirthday(, coincident = 2) starts to issue warnings (see (*) below) for larger number of classes (R 4.0.0, R-devel ./src/library/stats/R/birthday.R:47).
The default coincident = 2 is computed as 1 - prod((c:(c - n + 1))/rep(c, n)) where c = classes. Using exp(log(...)), one can derive the return value if(n > 0) 1 - exp(sum(log1p(-(0:(n-1))/c))) else 0. Simplifying this a bit further one obtains if(n >= 2) 1 - exp(sum(log1p(-(1:(n-1))/c))) else 0. For large c, sum(log1p(-(1:(n-1))/c)) is close to 0, so a more robust version would be to return if(n >= 2) -expm1(sum(log1p(-(1:(n-1))/c))) else 0 in the default case 'coincident = 2' (internally: if (k == 2) ...). ## Auxiliary function *just* considering 'coincident = 2' pbirthday2 <- function(n, classes = 365) { c <- classes # as pbirthday() if(n >= 2) -expm1(sum(log1p(-(1:(n-1))/c))) else 0 # return value suggested for the case 'if (k == 2) ...' } pbirthday (3, classes = 2) == pbirthday2(3, classes = 2) # identical pbirthday (3, classes = 2^53) == pbirthday2(3, classes = 2^53) # not identical anymore... stopifnot(all.equal(pbirthday (3, classes = 2^53), pbirthday2(3, classes = 2^53))) # ... but numerically equal pbirthday (3, classes = 2^54) # warnings start to appear (*) pbirthday2(3, classes = 2^54) # fine pbirthday (3, classes = 2^56) == 0 # numerically indistinguishable from 0 pbirthday2(3, classes = 2^56) # fine Cheers, M ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel