Coin Game Lab

R notebook quick tutorial

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

Simulation Experiments

Simulation provides a straightforward way of approximating probabilities. One simulates a particular random experiment a large number of times, and a probability of an outcome is approximated by the relative frequency of the outcome in the repeated experiments.

The use of simulation experiments to better understand probability patterns is called the Monte Carlo method. We focus on two particular R functions that simplify the process of programming simulation experiments. The sample() function will take samples with or without replacement from a set, and the replicate() function is helpful in repeating a particular simulation experiment.

NOTE: remember to use the set.seed(2020) function in each code chunk in which you are generating random numbers to obtain the desired results.

Simulating a Game of Chance of coin tossing

Peter and Paul play a simple game involving repeated tosses of a fair coin. In a given toss, if heads is observed, Peter wins 1€ from Paul; otherwise if tails is tossed, Peter gives 1€ to Paul. If Peter starts with zero €, we are interested in his fortune as the game is played for 50 tosses.

We can simulate this game using the R sample() function. Peter’s winning on a particular toss will be 1€ or -1€ with equal probability. His winnings on 50 repeated tosses can be considered to be a sample of size 50 selected with replacement from the set {-1, 1}.

set.seed(2020)
FlipCoin <- function(n) sample(c(-1, 1), n, replace = TRUE)
flips = FlipCoin(50)
print(flips)

 [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  1  1  1  1  1  1
[34] -1  1 -1  1  1  1 -1  1  1  1 -1  1  1  1 -1  1 -1

Desired Result

  [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  1  1  1  1  1  1 -1  1 -1  1  1  1 -1  1  1  1 -1  1  1  1 -1  1 -1

In the outcome represented above, Peter lost the first two tosses, won the third, lost the fourth, and so on….

Exploring Cumulative Winnings

Suppose Peter is interested in his cumulative winnings as he plays this game. We store his individual toss winnings in the variable win. The function cumsum() will compute the cumulative winnings of the individual values and the cumulative values are stored in cum_win.

just adds the individual wins/losses up:

cum_win=cumsum(flips)
cum_win

 [1]  1  2  1  2  3  2  1  2  3  4  3  2  3  4  5  6  7  8  9  8  9 10  9 10  9 10 11 12 13 14 15 16 17
[34] 16 17 16 17 18 19 18 19 20 21 20 21 22 23 22 23 22

Desired Result

 [1]  1  2  1  2  3  2  1  2  3  4  3  2  3  4  5  6  7  8  9  8  9 10  9 10  9 10 11 12 13 14 15 16 17 16 17 16 17 18 19 18 19 20 21 20 21 22 23 22 23 22

We extend this and plot the sequence of cumulative winnings for 4 games. In order to plot the cumulative winnings we will use a line plot, this kind of plots can be printed using the plot() function, including the cumulative winnings as the vector we want to represent and type = "l" as a parameter to obtain a line plot.

# Extend this and plot the sequence of cumulative winnings for 4 games

set.seed(2020)
par(mfrow=c(2,2))
matrix4games <- data.frame(replicate(4,FlipCoin(50)))
colnames(matrix4games) <- c('First', 'Second', 'Third','Fourth')
for (col in names(matrix4games)){
  plot(apply(matrix4games[col],2,cumsum),type = "l",main = col,xlab="Index",
        ylab="cumsum(win)",)
  abline(h=0, col="black")
}

Desired Result Below you can see the four plots you have to obtain:

Is evidently much variability. How do we address:

• 1. What is prob of Peter breaking even after 50 games?
• 2. How often Peter will be on the lead?
• 3. What will be value of Peter’s fortune?

First question can be calculated exactly, can only approximate 2) and 3) using Monte Carlo methods

R Function to Implement a Monte Carlo Experiment

One can obtain approximate answers to these questions by a Monte Carlo methods. In this type of experiment, we will simulate the random process and computes the statistic or statistics of interest. By repeating the random process many times, one obtains a collection of the statistics. We will then use the collection to approximate probabilities or expectations that answer the questions.

Let’s first consider Peter’s fortune F at the end of the game (question 3). We write a function peter_paul that simulates the fortunes for the 50 tosses and computes F that is equal to the sum of the individual fortunes. To make this function more general, we define n to be the number of tosses and let the default value of n be 50.

# Write a peter_paul function that simulates the fortunes for 50 tosses and calculate F, 
# that is the sum of the individual fortunes
set.seed(2020)
peter_paul=function(n=50){
  flip=FlipCoin(n)
  sum(flip)
}
F = peter_paul()
F

[1] 22

Desired Result

 [1] 22

In this game, Peter finished with a fortune of 22€. To repeat this for 1000 games, we use the replicate function with arguments 1000, the number of games, and the name of the function peter_paul() to repeat. The output of replicate is the vector of fortunes for the 1000 games that we assign to the variable F. Print the header of variable F to see the result.

# Use the replicate() function to repeat the experiment 1000 times
set.seed(2020)
F=replicate(1000, peter_paul())
F[1:6]

[1] 22 -2  2 -8 -2 -4

Desired Result

[1] 22 -2  2 -8 -2 -4

Summarizing the Monte Carlo Results

Since Peter’s fortune is integer-valued, a convenient way of summarizing the collection of values of F using table().

# Use the function table() to sumarize the values of F
table(F)

F
-22 -20 -18 -16 -14 -12 -10  -8  -6  -4  -2   0   2   4   6   8  10  12  14  16  18  20  22 
  4   3   5  15   8  31  31  69  72 111 100 120  95 109  79  61  36  20  14   5   5   6   1 

Desired Result

F
-22 -20 -18 -16 -14 -12 -10  -8  -6  -4  -2   0   2   4   6   8  10  12  14  16  18  20  22 
  4   3   5  15   8  31  31  69  72 111 100 120  95 109  79  61  36  20  14   5   5   6   1 

We can display the frequencies of F using plot() applied to the table output, and see what will be the value of Peter’s best fortune (question 3)

# Plot
plot(table(F))

Desired Result Below you can see the plot you have to obtain:

What is Peter’s chance of breaking even? (question 1) It is the ratio (approximated) of Peter finishing with 0 out of the 1000

But that answer can also be determined exactly. He breaks even if there are exactly n/2 heads in a binomial experiment of n trials with a probability of success equal to 0.50. So if n = 50 we can apply the dbinom function to obtain the exact value:

# Compute P(X=25) for X Binomial(n=50,p=0.5)
?dbinom
dbinom(25, 50, 0.5)

[1] 0.1122752

print(120/1000)

[1] 0.12

Note how close this exact number is to our approximated answer from the Monte Carlo simulation, that should be 120/1000 = 0.12

Modifying the Experiment to Learn About

New Statistics

We can add additional lines of code to our function peter_paul to compute several statistics of interest in our experiment. To answer our questions, we focus on the final fortune F, the number of times Peter is in the lead L, and the maximum cumulative winning M.

In the function, we define the vector of cumulative winnings cum_win. Here the output of the function is a vector consisting of the values of F, L, and M.

By naming the components (using, for example, F=sum(win)), we get more attractive output.

# Create a function to obtain the final fortune (F), the number of times Peter is in the lead (L) and 
# the maximum cumulative winning (W)
set.seed(2020)
peter_paul=function(n=50){
  flip=FlipCoin(n)
  cum_win=cumsum(flip)
  c(F =sum(flip), L=sum(cum_win>0),  M=max(cum_win))    
}

We simulate the game once

peter_paul()

 F  L  M 
22 50 23 

Desired Result

 F  L  M 
22 50 23 

In this game (with seed = 2020), Peter’s final fortune was 22€, he was in the lead for 50 plays, and his maximum total winning during the game was 23€. To repeat this game 1000 times, we again use replicate and store output in variable S.

Since the output of peter_paul is a vector, S will be a matrix of 3 rows and 1000 columns, where the rows correspond to the simulated draws of F, L, and M. We can verify the dimension of the matrix of S using the dim() function.

# Simulate 1000 times
set.seed(2020)
peter_paul=function(n=50){
  flip=FlipCoin(n)
  cum_win=cumsum(flip)
  c(F =sum(flip), L=sum(cum_win>0),  M=max(cum_win))    
}
S = replicate(1000,peter_paul())
dim(S)

[1]    3 1000

Desired Result

[1]    3 1000

How many times Peter in the lead? The likely answer to this question (we do not know for sure) is in the L row. Create the times_in_lead variable with all the values and print the header

# How many times is Peter in lead?
times_in_lead = S["L",]
times_in_lead[1:6]

[1] 50  4 47 17 17  9

Desired Result

[1] 50  4 47 17 17  9

We tabulate the simulated values using the table function, and the prop.table function will find the corresponding relative frequencies. We plot the result in a line graph:

# plot

plot(prop.table(table(times_in_lead)))

Desired Result Below you can see the plot you have to obtain:

The pattern of this graph for L isn’t what most people expect. Since it is a fair game, one might think it would be most likely for Peter to be ahead in 25 out of 50 tosses.

Instead, the most likely values are the extreme values L = 0,1,50, and the remaining values of L appear equally likely. So actually it is relatively common for Peter to always be losing or always winning in this game.

Last, let’s consider the distribution of M, Peter’s maximum winning during the 50 plays. We store the 1000 simulated values of M in the variable maximum_lead, tablulate the values using table(), and then plot them

# Plot the Maximum winings in all the experiments
maximum_lead = S["M",]
#maximum_lead[1:6]
plot(prop.table(table(maximum_lead)))

Desired Result Below you can see the plot you have to obtain:

From the plot, we can see that it is most likely for Peter to have maximum winnings of 1 or 2, but values between 3 and 6 are all relatively likely.

To compute the approximate probability that Peter will have a maximum winning of 10 or more, we find the number of values of maximum_lead that are 10 or greater and divide this sum by the number of simulation iterations 1000:

# What is the number of times Peter have a maximum lead greater than 10?
sum(maximum_lead >= 10)/1000

[1] 0.158

Since this probability is only about 16%, it is relatively rare for Peter to have a maximum winning of 10 or more.

Apply what you have learned

If you have managed to get this far you have secured 80% of the grade.

In this last section of the lab you will have to be creative and generate new questions about the coin game that you can later solve using R-simulation experiments.

Generate at least one new question and solve it.

Question:

What are the probabilities of Peter to win consecutively for 50 tosses ?

To demonstrate the idea, we will toss the coin only for 10 times in only 1 game

set.seed(2020)
flip=FlipCoin(10)
flip

 [1]  1  1 -1  1  1 -1 -1  1  1  1

We can count the consecutive by using function rle. Here lengths represent number of sequence of both -1 and 1. For example:

• 1. It shows here that only -1 has sequence of 1
• 2. -1 has 1 frequency of having 1 consecutive (lost)
• 3. -1 has 1 frequency of having 2 consecutive (lost)
• 4. 1 has 2 frequencies of having 2 consecutive(win)
• 5. 1 has 1 frequency of having 3 consecutive(win)

flip_st<-rle(flip)
table(flip_st)

       values
lengths -1 1
      1  1 0
      2  1 2
      3  0 1

Now we will calculate first the consecutive winning and sum the rest and assign zero value to show events without consecutive winning

win2con = flip_st$lengths[which(flip_st$values==1)]
other_win2con = flip_st$lengths[which(flip_st$values== - 1)]
other_win2con = sum(other_win2con)
other_win2con

[1] 3

win2con

[1] 2 2 3

el_win1orlost = 10 - sum(win2con)
el_win1orlost

[1] 3

We will need to create entries for event oof no consecutive win or losses and assign this variable win1orlost. We use function rep to add zero for those event. We combine both win2more and win1orlost and assign new table final_tab

win1orlost<- rep(0,el_win1orlost)
final_tab = append(win2con,win1orlost)
final_tab

[1] 2 2 3 0 0 0

Now let us calculate the frequencies of the consecutive winning and have the probability table

table (final_tab)

final_tab
0 2 3 
3 2 1 

prop.table(table (final_tab))

final_tab
        0         2         3 
0.5000000 0.3333333 0.1666667 

By tossing the coin by 10 times, we know that we going to have 3 consecutive win only once, 2 consecutive win twice and others which consist of only winning once or lost by thrice.

Let us put this in a **function and run 50 tossing coins and replicate this 1000 times to get the distribution probability of consecutive winning.

# Write your code here
set.seed(2020)
peter_paul_streak <- function(n=50){
  
  flip=FlipCoin(n)
  flip_st<-rle(flip)
  win2con = flip_st$lengths[which(flip_st$values==1)]
  other_win2con = flip_st$lengths[which(flip_st$values== - 1)]
  other_win2con = sum(other_win2con)
  el_win1orlost = n - sum(win2con)
  ##assigned zero for other condition (no consecutive win or lost)
  win1orlost<- rep(0,el_win1orlost)
  final_tab = append(win2con,win1orlost)
  
  
}
# We use Monte Carlo approach by running the function 1000 times for 50 
Winning2More = replicate(1000,peter_paul_streak(50))
# we do unlist to handle different size of list from the consecutive result
table(unlist(Winning2More))

    0     1     2     3     4     5     6     7     8     9    10    11    12 
25102  6462  3188  1564   767   350   191    97    47    16    16     3     1 

Table above shows the frequency of the winning in sequence from 1 to 12, meanwhile, zero represent no consecutive win or lost. We can see that the number of maximum 12 consecutive win which happen only once. And we can see that the number of consecutive winning is decreasing

round(prop.table(table(unlist(Winning2More))),3)

    0     1     2     3     4     5     6     7     8     9    10    11    12 
0.664 0.171 0.084 0.041 0.020 0.009 0.005 0.003 0.001 0.000 0.000 0.000 0.000 

Table here shows the probability of the consecutive winning. We can see that the probability of winning 6 or more consecutively is very low (below 0.00)

plot(prop.table(table(unlist(Winning2More))), main = paste("Probability of Winning in Sequence"),xlab="Frequencies",
ylab="probability")

Here is the chart