################################################################# library( moments ) # unbuffer the output ################################################################# population <- c( 35, rep( -1, 37 ) ) print( mu <- mean( population ) ) # [1] -0.05263158 print( sigma <- sqrt( mean( population^2 ) - mu^2 ) ) # [1] 5.762617 print( skewness <- mean( ( ( population - mu ) / sigma )^3 ) ) # [1] 5.918364 skewness( population ) # [1] 5.918364 print( kurtosis <- mean( ( ( population - mu ) / sigma )^4 ) - 3 ) # [1] 33.02703 kurtosis( population ) # [1] 36.02703 ################################################################# seed <- 1 set.seed( seed ) M <- 100000 S.star.1 <- rep( NA, M ) n <- 35 system.time( S.star <- apply( matrix( sample( population, n * M, replace = T ), M, n ), 1, sum ) ) plot( density( S.star ), main = 'PDF of S', xlab = 'Sum (S)', lwd = 2 ) E.s <- n * mu SD.s <- sigma * sqrt( n ) s.grid <- seq( E.s - 3 * SD.s, E.s + 3 * SD.s, length = 500 ) lines( s.grid, dnorm( s.grid, E.s, SD.s ), lwd = 2, lty = 2, col = 'red' ) text( E.s, 0.8 * max( dnorm( s.grid, E.s, SD.s ) ), 'CLT', col = 'red' ) ################################################################# c( skewness( S.star ), kurtosis( S.star ) - 3 ) # [1] 0.5856234 0.3145210 mean( S.star > 0 ) # [1] 0.492046 1 - pnorm( - E.s / SD.s ) # [1] 0.4636141 mean( S.star ) # [1] -1.790144 ################################################################# 1 - ( 1 - 1 / 38 )^34 # [1] 0.5961535 ################################################################# win.probability.analytic <- function( n ) { return( 1 - sum( dbinom( 0:floor( n / 36 ), n, 1 / 38 ) ) ) } n.max <- 1000 win.p.analytic <- rep( NA, n.max ) normal.approximation <- rep( NA, n.max ) for ( i in 1:n.max ) { win.p.analytic[ i ] <- win.probability.analytic( i ) normal.approximation[ i ] <- 1 - pnorm( ( 19 + 2 * i ) / ( 36 * sqrt( 37 * i ) ) ) } par( mfrow = c( 1, 1 ) ) plot( 1:n.max, win.p.analytic, type = 'l', xlab = 'Number of $1 bets ( n )', ylab = 'P (Come Out Ahead)', ylim = c( 0, 0.61 ) ) lines( 1:n.max, normal.approximation, lty = 2, lwd = 2, col = 'red' ) abline( h = 0.5, lwd = 2, col = 'blue' ) ############################################################# cbind( 1:100, win.p.analytic[ 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.61712890 <- maximized at n = 36 [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.56815296 [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 ############################################################# par( mfrow = c( 1, 1 ) ) plot( 1:n.max, - ( 1:n.max ) / 19, type = 'l', lwd = 2, xlab = 'Number of $1 bets', ylab = 'Your Net Gain ($)', ylim = c( -500, 500 ) ) lines( 1:n.max, - ( 1:n.max ) / 19 - qnorm( 0.975 ) * ( 36 * sqrt( 37 * 1:n.max ) / 38 ), lty = 2, lwd = 2, col = 'red' ) lines( 1:n.max, - ( 1:n.max ) / 19 + qnorm( 0.975 ) * ( 36 * sqrt( 37 * 1:n.max ) / 38 ), lty = 2, lwd = 2, col = 'red' ) abline( h = 0, lwd = 2, col = 'blue' ) ############################################################# population.data.set <- c( rep( -1, 37 ), 35 ) win.probability.simulation <- function( n, M ) { p <- 0 for ( m in 1:M ) { y <- sample( population.data.set, n, replace = T ) p <- p + ( sum( y ) > 0 ) } return( p / M ) } seed <- 1 set.seed( seed ) M <- 100000 win.probability.simulation( 1, M ) win.p.simulation <- rep( NA, 100 ) for ( i in 1:100 ) { win.p.simulation[ i ] <- win.probability.simulation( i, M ) print( i ) } plot( 1:100, win.p.simulation, type = 'l' ) n <- 36 y <- matrix( NA, M, n ) for ( i in 1:M ) { y[ i, ] <- sample( population.data.set, n, replace = T ) } S <- apply( y, 1, sum ) hist( S, breaks = 1000, prob = T ) table( S ) # -36 0 36 72 108 144 180 216 # 38289 37509 17438 5350 1139 245 25 5 ############################################################# w or lw or llw or lllw or ... or ll ... llw (34 losses in a row) p <- 1 / 38 for ( i in 1:34 ) { p <- p + ( 37 / 38 )^i / 38 } print( p ) [1] 0.606781 ############################################################# population.data.set <- c( rep( -1, 37 ), 35 ) M <- 5000 n <- 1000 S <- NULL S.grid <- seq( -600, 500, length = 500 ) set.seed( 3 ) for ( i in 1:M ) { S <- c( S, sum( y <- sample( population.data.set, n, replace = T ) ) ) hist( S, breaks = 18, prob = T, xlim = c( -600, 500 ), ylim = c( 0, 0.003 ), xlab = 'Your Net Gain ($)', main = '' ) lines( S.grid, dnorm( S.grid, - n / 19, 36 * sqrt( 37 * n ) / 38 ), lty = 1, lwd = 2, col = 'red' ) } #############################################################