Computing Probabilities

Learn to simulate many times with replicate()

Student writing code to simulate dice rolls
  • Review of set.seed(), sample()
  • replicate() for repeatedly running the same code many times, storing all the results in a vector when each run returns a single value
  • Don’t confuse replicate() with rep()! (rep()is for creating a vector of many repeated copies of something, as a vector)

Review: set.seed()

If you’re writing code that does random things (like sample), you may get a different answer every time you run your code. Run the cell below a few times, you’ll see that the results change.

If you want to make sure that you always get the SAME random results each time, use set.seed(any number) right before the code that uses randomness:

Run this cell many times - always the same result

Review: sample()

This is how to draw random elements from a vector, like drawing from a box of tickets. We’ll use these three arguments:

  1. x: the vector to be sampled from, this must be specified
  2. size: the number draws to take from that vector.
  3. replace: whether we replace a drawn item before we draw again. If we don’t set this argument, it defaults to FALSE (drawing without replacement).

Example: one sample of size 2 (WITHOUT replacement) from a box with tickets from 1 to 10

Read, then run the code:

Notice two things:

  1. The first argument doesn’t have a name. sample assumes the first argument is x, the vector to sample from. We could have said sample(x = box, size = 2) and it would have done the same thing. Edit the code above to try it if you like.
  2. The replace argument is not specified, so it falls back on its default value (FALSE in this case). This is what we want anyway, so we can just skip specifying it. We could have done sample(x = box, size = 2, replace = FALSE) and gotten the same behavior.

Example: one sample of size 2 (WITH replacement) from a box with tickets from 1 to 10

Read, then run the code:

The usefulness of replicate()


BUT FIRST: replicate() ≠ rep()!

rep() is for creating a vector of many repeated copies of something. E.g. rep(1, times = 5) gives you a vector of five 1’s: c(1, 1, 1, 1, 1).

On to replicate()

Let’s say we want two people to play a game. Each person rolls a die 5 times, then adds the results up as their “score.” We want to store their results in a vector. For this two-player case, it’s easy to do with familiar code:

Read, then run the code:

…but this quickly gets unwieldy if we want to have 100 people do this. Our code would have to look something like:

die = 1:6
score1 = sum(sample(die, size = 5, replace = TRUE))
score2 = sum(sample(die, size = 5, replace = TRUE))
score3 = sum(sample(die, size = 5, replace = TRUE))
# ... 100 lines of these!
scores = c(score1, score2, score3, ..., score100)

This is a headache. Let alone trying to do it for 1,000,000 people! Instead, we can use replicate() to do all of this in one line of code.

replicate() takes two required arguments:

  1. n: the number of times to do something
  2. expr: the thing we want to do (code to run)

The result of replicate() is a vector with n values, representing the results of each run of the code we specified. Let’s try it out.

Now we can rewrite our code in dramatically cleaner fashion. Again, we can skip naming these arguments, as replicate() assumes the first argument is n and the second is expr.

Read, then run the code:

Or even in one line by moving the die vector directly into the sample() function!

Read, then run the code:

You try: Simulate a game where 50 people play. Each person draws two tickets, without replacement, from a box with tickets numbered 1 through 20. Store all the sums of everyone’s draws in a vector called sums.

Check your answer

Another example of where replicate() is useful

Let’s say we want to not only run the experiment above, but also plot the distribution of scores. It would be nice to see the distribution to get a sense of what a good/bad/average score would look like.

BUT before we use ggplot(), we need a dataframe, not a vector! This is easy enough, as you can just call data.frame(sums) to create a one-column dataframe, with one column named sums with all of our data. Then we can pipe it to ggplot().

First, observe the simple dataframe creation. (Read, then run the code.)

Now let’s put it all together in one pipe to get plotting! (Read, then run the code.)

(Advanced): Simulating de Mere’s dice game

This code is harder to parse, but it’s absolutely within your power. See if you can follow.

Game 1: “At least one six in four rolls of a die”

Read, then run the code.

Game 2: “At least one double-six in 24 rolls of two dice”

Read, then run the code.

This requires TWO uses of replicate! Again, much harder code to follow, but give it a try.