Probability Foundations

Combinatorics and Simulations in R

Introduction

In this tutorial, we will see two uses of R to help us compute probabilities

• R has built in factorial, exponential, and n choose k functions so we can get numbers for expressions like

\[ \frac{\binom{10}{4}\,2^{6}}{3^{10}} \approx 0.228\]

• R has built in random number generators. It can simulate sampling a vector with or without replacement. We can make R repeat an experiment many times and get frequencies of an event to compare with our computed probabilities.

R as a Calculator: Factorial, N choose K, and Exponential

In R, the code for “5 to the power of 4” is an up arrow

In R, the code for 6 factorial is factorial(6)

In R, the code for “N Choose K” is choose(N,K)

Let’s put this all together to get a value for the penguin bowtie expression from before,

\[\frac{\binom{10}{4}\,2^{6}}{3^{10}}\]

Sampling With and Without Replacement using sample()

Our first step toward simulating experiments is introducing randomness in R. The following three functions are a good start.

sample()

This function takes in a vector of possible values (the state space), and randomly picks a value from the list with or without replacement. We can also pick from the list more than one time. The inputs are:

• Arguments
  • x: the vector to be sampled from, this must be specified
  • size: the number of items to be sampled, the default value is the length of x
  • replace: whether we replace a drawn item before we draw again, the default value is FALSE, indicating that we would draw without replacement.

Example: Simulate two dice rolls, sampling with replacement

What would happen if we don’t specify values for size and replace?

By default, size was set to 6 and replace was FALSE, so this is exactly a permutation of our vector, one of the \(6!\) possible orderings where each entry appears only once. (Go hit run a few more times, you will see different outputs!)

Note the input list does not have to be numbers. Recall our friend the penguin who picks one of 3 bowties a day for 10 days, let’s simulate a single 10 day output for him by picking a bowtie 10 times

Let’s check how many red bowties the penguin wore in this sample, and see if it equals 3

So here we have computed how many red bowties did the penguin wear, and a logical output that tells us did the penguin wear exactly 3 red bowties in this sample.

set.seed()

The random number generator in R is called a “Pseudo Random Number Generator”, because the process can be controlled by a “seed number”. These numbers are not actually “perfectly random”, but are emulating randomness. We can make a “random” R code chunck produce the same number every time by setting what is called the seed.

The seed is the starting value for an algorithmic random number generator, which means that if you provide the same seed (a starting number), R will generate the same sequence of random numbers. This makes it easier to debug your code, and reproduce your results if needed.

• Arguments
  • n: the seed number to use. You can use any number you like, for example 1, or 31415 etc You might have noticed that each time you run sample in the code chunk above, it gives you a different sample. Sometimes we want it to give the same sample so that we can check how the code is working without the sample changing each time. We will use the set.seed function for this, which ensures that we will get the same random sample each time we run the code.

Example: Roll a Die Twice

Example: Roll Die with the same seed

Notice that we get the same sample. You can try to run sample(die) without using set.seed() and see what happens.

Though we used set.seed() twice here to demonstrate its purpose, generally, you will only need to run set.seed() one time per document. This is a line of code that fits perfectly at the beginning of your work, when you are also loading libraries and packages.

seq()

Above, we created the vector die using die <- c(1, 2, 3, 4, 5, 6), which is fine, but this method would be tedious if we wanted to simulate a 20-sided die, for instance. The function seq() allows us to create any sequence we like by specifying the starting number, how we want to increment the numbers, and either the ending number or the length of the sequence we want.

• Arguments
  • from: where to start
  • by: size of jump from number to number (the increment)

You can end a sequence in one of two ways: - to: at what number should the sequence end - length: how long should the sequence be

Example: sequence with the to argument

Example: sequence with the length argument

The For Loop: Running a Line of Code Many, Many Times

In our simulations, we want to

• Repeat an experiment many times
• Check for each individual trial if the event occurred or not
• Compute the proportion of times the event occurred to estimate probability

To run the same line of code in R, we use a for loop

The for loop has three parts to it:

• A block of code we want to run many times
• The number of times we want to run that code, called N
• An indicator variable, usually called i, that will increase by 1 in each run of the code block. So i is 1, then 2, then 3, up to the final run when it is N.

(Note: you don’t have to call them i and N, you can use whatever variable names you want, but I will use them here for consistency.)

The code structure of a for loop is: (note this chunck will not do anything since I put no code in the loop)

So we say for, then in smooth brackets () we say how many times we want to code to run, and the name of the indicator variable. Then we enclose the code we want to run in {} curly brackets.

Let’s see an example of adding up the numbers 1+2+3+4+5+6

The first time this code runs, \(i\) is 1, total_count is 0, and 0+1 is 1. The second time it runs, i is now 2, total_count is 1, and 1+2 is 3. And so on.

The reason we want to use this is to simulate many repeated experiments and track the total number of times an event occurs. For example, let’s roll a die 100 times and check how many times it is even.

This is the general code structure to simulate many possible outcomes and see if an event occurs. What Null does is it makes a vector but doesn’t put any information in it yet, so it can be filled in later. Here we will see code that computes the probabilities from the 4 practice exampeles in the notes:

1.

Say I roll 2 dice and add up the sum of the dice. So rolling a 2 and a 4 adds up to 6. Simulate 10,000 trials and track the number that add up to 7. (From notes: should be close to 0.167)

2.

Say I am baking a cake.

• The flavor can be vanilla, chocolate, or strawberry (3 options)
• The icing can be buttercream, cream cheese, fondant, or whipped cream (4 options)
• The cake type can be an ice cream cake or a baked cake (2 options)

I make a selection of each option, each equally likely. Simulate 10,000 trials and track the number of times I get buttercream icing or strawberry flavor (From notes: should be close to 0.5)

3.

Say I have an urn with 4 red balls and 5 blue balls. I reach into the urn 3 times and pick a ball, and I do not replace the ball each time.Simulate 10,000 trials and count What is the probability I draw 3 red balls in a row? (From notes: should be close to 0.048)

4.

Every morning for 10 days, a penguin reaches into its wardrobe and picks one of 3 bowties — red, blue, or green — at random, with each equally likely and each day independent of the others. Simulate 10,000 trials and compute What is the probability the penguin wears the red bowtie on exactly 4 of the 10 days? (From notes: should be close to 0.228)

Summary

So, we can see 4 examples of general sampling code here using a for loop where you can plug in your own state space, sampling line, and event check. The structure should be mostly the same for the kinds of experiments you would like to simulate.

The sample() function is the main base random number generator, it can pick numbers from a vector either with or without replacement as many times as we like.