| # File src/library/stats/R/birthday.R |
| # Part of the R package, https://www.R-project.org |
| # |
| # Copyright (C) 1995-2012 The R Core Team |
| # |
| # This program is free software; you can redistribute it and/or modify |
| # it under the terms of the GNU General Public License as published by |
| # the Free Software Foundation; either version 2 of the License, or |
| # (at your option) any later version. |
| # |
| # This program is distributed in the hope that it will be useful, |
| # but WITHOUT ANY WARRANTY; without even the implied warranty of |
| # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| # GNU General Public License for more details. |
| # |
| # A copy of the GNU General Public License is available at |
| # https://www.R-project.org/Licenses/ |
| |
| |
| qbirthday <- function(prob = 0.5, classes = 365, coincident = 2) |
| { |
| k <- coincident |
| c <- classes |
| p <- prob |
| if (p <= 0) return(1) |
| if (p >= 1) return(c*(k-1)+1) |
| ## We need smallest n with pbirthday(n, c, k) >= prob |
| ## This is a crude inversion of Diaconis & Mosteller expression (7.5), |
| ## usually an underestimate. |
| N <- exp(((k-1)*log(c) + lgamma(k+1) + log(-log1p(-p)))/k) |
| N <- ceiling(N) |
| if(pbirthday(N, c, k) < prob) { |
| N <- N+1 |
| while(pbirthday(N, c, k) < prob) N <- N+1 |
| } else if (pbirthday(N-1, c, k) >= prob) { |
| N <- N-1 |
| while(pbirthday(N-1, c, k) >= prob) N <- N-1 |
| } |
| N |
| } |
| |
| pbirthday <- function(n, classes = 365, coincident = 2) |
| { |
| k <- coincident |
| c <- classes |
| if (k < 2) return(1) |
| if (k == 2) return( 1 - prod((c:(c-n+1))/rep(c, n)) ) |
| if (k > n) return(0) |
| if (n > c*(k-1)) return(1) |
| ## use Diaconis & Mosteller expression (7.5) on log scale |
| LHS <- n * exp(-n/(c*k))/(1 - n/(c*(k+1)))^(1/k) |
| lxx <- k*log(LHS) - (k-1)*log(c) - lgamma(k+1) |
| -expm1(-exp(lxx)) |
| } |
| |