[Rd] pbirthday() for larger number of classes
Marius Hofert
m@r|u@@ho|ert @end|ng |rom uw@ter|oo@c@
Fri May 22 13:44:50 CEST 2020
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
More information about the R-devel
mailing list