naive.implementation <- function( n, k ) { return( 1 - factorial( n ) / ( factorial( n - k ) * n^k ) ) } naive.implementation( 365, 2 ) [1] NaN Warning messages: 1: In factorial(n) : value out of range in 'gammafn' 2: In factorial(n - k) : value out of range in 'gammafn' factorial( 100 ) [1] 9.332622e+157 factorial( 200 ) [1] Inf Warning message: In factorial(200) : value out of range in 'gammafn' 1e308 [1] 1e+308 1.8 * 1e308 [1] Inf .Machine $double.eps [1] 2.220446e-16 $double.neg.eps [1] 1.110223e-16 $double.xmin [1] 2.225074e-308 $double.xmax [1] 1.797693e+308 $double.base [1] 2 $double.digits [1] 53 $double.rounding [1] 5 $double.guard [1] 0 $double.ulp.digits [1] -52 $double.neg.ulp.digits [1] -53 $double.exponent [1] 11 $double.min.exp [1] -1022 $double.max.exp [1] 1024 $integer.max [1] 2147483647 $sizeof.long [1] 4 $sizeof.longlong [1] 8 $sizeof.longdouble [1] 16 $sizeof.pointer [1] 8 method.1 <- function( n, k ) { product <- 1 for ( i in 0:( k - 1 ) ) { product <- product * ( n - i ) / n } return( 1 - product ) } method.1( 365, 93 ) [1] 0.9999974 system.time( print( method.1( 365, 93 ) ) ) [1] 0.9999974 user system elapsed 0 0 0 n <- 365 k.max <- 150 match.probability <- rep( NA, k.max ) for ( k in 1:k.max ) { match.probability[ k ] <- method.1( n, k ) } cbind( 1:k.max, match.probability ) cbind( 1:k.max, match.probability ) match.probability [1,] 1 0.000000000 [2,] 2 0.002739726 [3,] 3 0.008204166 [4,] 4 0.016355912 [5,] 5 0.027135574 [6,] 6 0.040462484 [7,] 7 0.056235703 [8,] 8 0.074335292 [9,] 9 0.094623834 [10,] 10 0.116948178 [11,] 11 0.141141378 [12,] 12 0.167024789 [13,] 13 0.194410275 [14,] 14 0.223102512 [15,] 15 0.252901320 [16,] 16 0.283604005 [17,] 17 0.315007665 [18,] 18 0.346911418 [19,] 19 0.379118526 [20,] 20 0.411438384 [21,] 21 0.443688335 [22,] 22 0.475695308 [23,] 23 0.507297234 [24,] 24 0.538344258 [25,] 25 0.568699704 [26,] 26 0.598240820 [27,] 27 0.626859282 [28,] 28 0.654461472 [29,] 29 0.680968537 [30,] 30 0.706316243 [31,] 31 0.730454634 [32,] 32 0.753347528 [33,] 33 0.774971854 [34,] 34 0.795316865 [35,] 35 0.814383239 [36,] 36 0.832182106 [37,] 37 0.848734008 [38,] 38 0.864067821 [39,] 39 0.878219664 [40,] 40 0.891231810 [41,] 41 0.903151611 [42,] 42 0.914030472 [43,] 43 0.923922856 [44,] 44 0.932885369 [45,] 45 0.940975899 [46,] 46 0.948252843 [47,] 47 0.954774403 [48,] 48 0.960597973 [49,] 49 0.965779609 [50,] 50 0.970373580 [51,] 51 0.974431993 [52,] 52 0.978004509 [53,] 53 0.981138113 [54,] 54 0.983876963 [55,] 55 0.986262289 [56,] 56 0.988332355 [57,] 57 0.990122459 [58,] 58 0.991664979 [59,] 59 0.992989448 [60,] 60 0.994122661 [61,] 61 0.995088799 [62,] 62 0.995909575 [63,] 63 0.996604387 [64,] 64 0.997190479 [65,] 65 0.997683107 [66,] 66 0.998095705 [67,] 67 0.998440043 [68,] 68 0.998726391 [69,] 69 0.998963666 [70,] 70 0.999159576 [71,] 71 0.999320753 [72,] 72 0.999452881 [73,] 73 0.999560806 [74,] 74 0.999648644 [75,] 75 0.999719878 [76,] 76 0.999777437 [77,] 77 0.999823779 [78,] 78 0.999860955 [79,] 79 0.999890668 [80,] 80 0.999914332 [81,] 81 0.999933109 [82,] 82 0.999947953 [83,] 83 0.999959646 [84,] 84 0.999968822 [85,] 85 0.999975997 [86,] 86 0.999981587 [87,] 87 0.999985925 [88,] 88 0.999989280 [89,] 89 0.999991865 [90,] 90 0.999993848 [91,] 91 0.999995365 [92,] 92 0.999996521 [93,] 93 0.999997398 [94,] 94 0.999998061 [95,] 95 0.999998560 [96,] 96 0.999998935 [97,] 97 0.999999215 [98,] 98 0.999999424 [99,] 99 0.999999578 [100,] 100 0.999999693 [101,] 101 0.999999777 [102,] 102 0.999999839 [103,] 103 0.999999884 [104,] 104 0.999999917 [105,] 105 0.999999940 [106,] 106 0.999999957 [107,] 107 0.999999970 [108,] 108 0.999999979 [109,] 109 0.999999985 [110,] 110 0.999999989 [111,] 111 0.999999993 [112,] 112 0.999999995 [113,] 113 0.999999996 [114,] 114 0.999999998 [115,] 115 0.999999998 [116,] 116 0.999999999 [117,] 117 0.999999999 [118,] 118 0.999999999 [119,] 119 1.000000000 [120,] 120 1.000000000 [121,] 121 1.000000000 [122,] 122 1.000000000 [123,] 123 1.000000000 [124,] 124 1.000000000 [125,] 125 1.000000000 [126,] 126 1.000000000 [127,] 127 1.000000000 [128,] 128 1.000000000 [129,] 129 1.000000000 [130,] 130 1.000000000 [131,] 131 1.000000000 [132,] 132 1.000000000 [133,] 133 1.000000000 [134,] 134 1.000000000 [135,] 135 1.000000000 [136,] 136 1.000000000 [137,] 137 1.000000000 [138,] 138 1.000000000 [139,] 139 1.000000000 [140,] 140 1.000000000 [141,] 141 1.000000000 [142,] 142 1.000000000 [143,] 143 1.000000000 [144,] 144 1.000000000 [145,] 145 1.000000000 [146,] 146 1.000000000 [147,] 147 1.000000000 [148,] 148 1.000000000 [149,] 149 1.000000000 [150,] 150 1.000000000 plot( 1:k.max, match.probability, type = 'l', xlab = 'Number of People in AMS 131 Class', ylab = 'P( at least one birthday match )', lwd = 2, col = 'blue' ) abline( h = 0.5, lwd = 2, col = 'red' ) abline( v = 22.5, lwd = 2, col = 'red' ) text( 60, 0.2, '>= 23 for probability > 0.5' ) method.2 <- function( n, k ) { return( 1 - exp( ( n - k + 0.5 ) * ( log( n ) - log( n - k ) ) - k ) ) } method.2( 365, 93 ) method.3 <- function( n, k ) { p <- 1 - exp( lgamma( n + 1 ) - lgamma( n - k + 1 ) - k * log( n ) ) return( p ) } method.3( 365, 93 ) ####################################################################### sum.1 <- function( x, y ) { z <- x + y return( z ) } sum.2 <- function( x, y ) { z <- x + y return( z ) } sum.3 <- function( x, y ) { z <- x + y return( z ) } #######################################################################