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

______________________________________________

[hidden email] mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel