# unbuffer the output

setwd( 'C:/DD/Teaching/AMS-131/Summer-2018' )

print( population.data.set <- c( rep( -1, 37 ), 35 ) )

#  [1] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
# [28] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 35

n <- 1000

M <- 1000000

random.number.generator.seed <- 4391557

set.seed( random.number.generator.seed )

simulated.sums <- rep( NA, M )

system.time( {

  for ( m in 1:M ) {

    simulated.sample.data.set <- sample( population.data.set, n, replace = T )

    simulated.sums[ m ] <- sum( simulated.sample.data.set )

  }

} 

)

#  user  system elapsed       n = 1000   M = 100000
#   2.7     0.0     2.7

#  user  system elapsed       n = 1000   M = 1000000 
# 26.31    0.03   26.36

estimated.probability.of.coming.out.ahead <- mean( simulated.sums > 0 )

paste( 'n =', n, 'M =', M, 'seed =', random.number.generator.seed,
  'estimated probability of coming out ahead =',
  estimated.probability.of.coming.out.ahead )

hist( simulated.sums, prob = T, main = '', xlab = 'Your Net Gain',
  breaks = 'FD' )

abline( v = 0, lwd = 2, col = 'red' )

table( simulated.sums ) / M

# simulated.sums

#     -100      -64      -28        8       44       80      116      152      188      224      260 
# 0.069067 0.187200 0.251781 0.221875 0.145672 0.075289 0.032287 0.011754 0.003770 0.000962 0.000265 

#      296      332      368 
# 0.000066 0.000009 0.000003