########## R code ############################## n <- 365 k <- 93 product <- 1 for ( i in 0:( k - 1 ) ) { product <- product * ( n - i ) / n } 1 - product # 0.9999974 ################################################ 1 - exp( ( n - k + 1 / 2 ) * ( log( n ) - log( n - k ) ) - k ) # 0.9999974 1 - exp( lgamma( n + 1 ) - lgamma( n - k + 1 ) - k * log( n ) ) # 0.9999974 birthday <- function( n, k ) { return( 1 - exp( lgamma( n + 1 ) - lgamma( n - k + 1 ) - k * log( n ) ) ) } birthday( 365, 93 ) # 0.9999974 birthday( 365, 10 ) # 0.1169482 # 0.4114384 plot( 1:40, birthday( 365, 1:40 ), type = 'l', xlab = 'Number of People', ylab = 'Probability of a Match' ) abline( h = 0.5, lwd = 2, col = 'red' ) birthday( 365, 23 ) # 0.5072972 abline( v = 23, lwd = 2, col = 'red' )