winning.probability.single.number <- function( n ) { return( 1 - pbinom( n / 36, n, 1 / 38 ) ) } cbind( 1:100, winning.probability.single.number( 1:100 ) ) [,1] [,2] [1,] 1 0.02631579 [2,] 2 0.05193906 [3,] 3 0.07688803 [4,] 4 0.10118045 [5,] 5 0.12483360 [6,] 6 0.14786429 [7,] 7 0.17028892 [8,] 8 0.19212342 [9,] 9 0.21338333 [10,] 10 0.23408377 [11,] 11 0.25423946 [12,] 12 0.27386473 [13,] 13 0.29297356 [14,] 14 0.31157952 [15,] 15 0.32969584 [16,] 16 0.34733543 [17,] 17 0.36451081 [18,] 18 0.38123421 [19,] 19 0.39751752 [20,] 20 0.41337232 [21,] 21 0.42880989 [22,] 22 0.44384121 [23,] 23 0.45847697 [24,] 24 0.47272758 [25,] 25 0.48660317 [26,] 26 0.50011361 [27,] 27 0.51326851 [28,] 28 0.52607724 [29,] 29 0.53854889 [30,] 30 0.55069234 [31,] 31 0.56251622 [32,] 32 0.57402896 [33,] 33 0.58523872 [34,] 34 0.59615349 [35,] 35 0.60678103 [36,] 36 0.24460566 [37,] 37 0.25440891 [38,] 38 0.26421932 [39,] 39 0.27402973 [40,] 40 0.28383334 [41,] 41 0.29362373 [42,] 42 0.30339479 [43,] 43 0.31314076 [44,] 44 0.32285621 [45,] 45 0.33253597 [46,] 46 0.34217522 [47,] 47 0.35176936 [48,] 48 0.36131411 [49,] 49 0.37080542 [50,] 50 0.38023949 [51,] 51 0.38961276 [52,] 52 0.39892189 [53,] 53 0.40816378 [54,] 54 0.41733551 [55,] 55 0.42643438 [56,] 56 0.43545786 [57,] 57 0.44440364 [58,] 58 0.45326954 [59,] 59 0.46205357 [60,] 60 0.47075391 [61,] 61 0.47936888 [62,] 62 0.48789694 [63,] 63 0.49633670 [64,] 64 0.50468691 [65,] 65 0.51294643 [66,] 66 0.52111425 [67,] 67 0.52918948 [68,] 68 0.53717134 [69,] 69 0.54505915 [70,] 70 0.55285233 [71,] 71 0.56055039 [72,] 72 0.29446066 [73,] 73 0.30166309 [74,] 74 0.30887352 [75,] 75 0.31608923 [76,] 76 0.32330754 [77,] 77 0.33052585 [78,] 78 0.33774163 [79,] 79 0.34495241 [80,] 80 0.35215579 [81,] 81 0.35934946 [82,] 82 0.36653114 [83,] 83 0.37369865 [84,] 84 0.38084986 [85,] 85 0.38798271 [86,] 86 0.39509521 [87,] 87 0.40218542 [88,] 88 0.40925149 [89,] 89 0.41629161 [90,] 90 0.42330405 [91,] 91 0.43028712 [92,] 92 0.43723923 [93,] 93 0.44415881 [94,] 94 0.45104438 [95,] 95 0.45789449 [96,] 96 0.46470778 [97,] 97 0.47148291 [98,] 98 0.47821864 [99,] 99 0.48491374 [100,] 100 0.49156707 plot.winning.probability.single.number.1 <- function( n ) { plot( 1:n, winning.probability.single.number( 1:n ), type = 'l', lwd = 2, col = 'red', xlab = 'Number of Plays', ylab = 'P( coming out ahead)', main = 'Single Number' ) abline( v = 35, lwd = 1, col = 'blue', lty = 2 ) abline( h = 0.60678103, lwd = 1, col = 'blue', lty = 2 ) abline( h = 0.5, lwd = 1, col = 'black', lty = 1 ) } ############################################################# plot.winning.probability.single.number.1( 100 ) plot.winning.probability.single.number.1( 500 ) plot.winning.probability.single.number.1( 1000 ) ############################################################# set.seed( 1 ) x.star.n.35 <- rbinom( 1000000, 35, 1 / 38 ) table( x.star.n.35 ) x.star.n.35 0 1 2 3 4 5 6 7 8 393069 372494 170902 50502 10975 1788 236 33 1 hist( x.star.n.35, breaks = ( 0:7 ) - 0.5, prob = T, main = '', xlab = 'x.star' ) 35 / 36 # [1] 0.9722222 1 - pnorm( 35 / 36, mean( x.star.n.35 ), sd( x.star.n.35 ) ) # [1] 0.4781048 1 - pnorm( 0.5, mean( x.star.n.35 ), sd( x.star.n.35 ) ) # [1] 0.6716528 s.star.n.35 <- 36 * x.star.n.35 - 35 table( s.star.n.35 ) s.star.n.35 -35 1 37 73 109 145 181 217 253 393069 372494 170902 50502 10975 1788 236 33 1 mean( s.star.n.35 ) # [1] -1.869236 population.mean <- mean( c( rep( -1, 37 ), 35 ) ) 35 * population.mean # [1] -1.842105 naive.normal.approximation.single.number <- function( n ) { mu <- n / 38 sigma <- sqrt( n * ( 1 / 38 ) * ( 1 - 1 / 38 ) ) return( 1 - pnorm( ( ( n / 36 ) - mu ) / sigma ) ) } plot.winning.probability.single.number.2 <- function( n ) { plot( 1:n, winning.probability.single.number( 1:n ), type = 'l', lwd = 2, col = 'red', xlab = 'Number of Plays', ylab = 'P( coming out ahead)', main = 'Single Number' ) abline( v = 35, lwd = 1, col = 'blue', lty = 2 ) abline( h = 0.60678103, lwd = 1, col = 'blue', lty = 2 ) abline( h = 0.5, lwd = 1, col = 'black', lty = 1 ) lines( 1:n, naive.normal.approximation.single.number( 1:n ), lwd = 2, col = 'cyan4' ) } ############################################################# plot.winning.probability.single.number.2( 100 ) plot.winning.probability.single.number.2( 500 ) plot.winning.probability.single.number.2( 1000 ) plot.winning.probability.single.number.2( 5000 ) ############################################################# winning.probability.roulette <- function( n, p, w ) { return( 1 - pbinom( n / ( w + 1 ), n, p ) ) } cbind( 1:100, winning.probability.roulette( 1:100, 1 / 38, 35 ) ) ############################################################# naive.normal.approximation.roulette <- function( n, p, w ) { mu <- n * p sigma <- sqrt( n * p * ( 1 - p ) ) return( 1 - pnorm( ( ( n / ( w + 1 ) ) - mu ) / sigma ) ) } plot.winning.probability.roulette <- function( n, p, w, title ) { plot( 1:n, winning.probability.roulette( 1:n, p, w ), type = 'l', lwd = 2, col = 'red', xlab = 'Number of Plays', ylab = 'P( coming out ahead)', main = title ) maximum.winning.probability <- max( winning.probability.roulette( 1:n, p, w ) ) n.maximum.winning.probability <- ( 1:n )[ winning.probability.roulette( 1:n, p, w ) == maximum.winning.probability ] abline( v = n.maximum.winning.probability, lwd = 1, col = 'blue', lty = 2 ) abline( h = maximum.winning.probability, lwd = 1, col = 'blue', lty = 2 ) abline( h = 0.5, lwd = 1, col = 'black', lty = 1 ) lines( 1:n, naive.normal.approximation.roulette( 1:n, p, w ), lwd = 2, col = 'cyan4' ) } plot.winning.probability.roulette( 500, 1 / 38, 35, 'Single Number' ) plot.winning.probability.roulette( 500, 2 / 38, 17, 'Split' ) plot.winning.probability.roulette( 500, 18 / 38, 1, 'Red And Black' ) # need a better normal approximation #############################################################