# Monte Carlo Simulations

Mimic process of picking at random we use the `sample()`

.

```
# probability red 2/5 or 40%
# and of blue 3/5 or 60%
(beads <- rep(c('red', 'blue'), times = c(2, 3)))
```

`## [1] "red" "red" "blue" "blue" "blue"`

```
# pick random sample
ifelse(sample(beads, 1)=='blue', 1, 0)
```

`## [1] 1`

```
# pick random sample
ifelse(sample(beads, 1)=='blue', 1, 0)
```

`## [1] 1`

```
# pick random sample
ifelse(sample(beads, 1)=='blue', 1, 0)
```

`## [1] 1`

`replicate()`

function can be used to repeat a task any number of times.

```
B <- 10000
# sample from beads 10000 times
events <- replicate(B, sample(beads, 1))
```

```
# use table to see the distribution
(tab <- table(events))
```

```
## events
## blue red
## 6021 3979
```

```
# prop.table() gives the proportion
prop.table(tab)
```

```
## events
## blue red
## 0.6021 0.3979
```

Monte Carlo simulatin gave a very good approximation.

## Combinations and Permutations

```
number <- "Three"
suit <- "Hearts"
paste(number, suit)
```

`## [1] "Three Hearts"`

`paste(letters[1:5], as.character(1:5))`

`## [1] "a 1" "b 2" "c 3" "d 4" "e 5"`

`expand.grid()`

gives all the combinations of two list

```
expand.grid(pants=c("blue", "black"),
shirt=c("white", "grid", "plaid"))
```

```
## pants shirt
## 1 blue white
## 2 black white
## 3 blue grid
## 4 black grid
## 5 blue plaid
## 6 black plaid
```

```
suits <- c("Diamonds", "Clubs", "Hearts", "Spades")
numbers <- c("Ace", "Deuce", "Three", "Four", "Five", "Six", "Seven",
"Eight", "Nine", "Ten", "Jack", "Queen", "King")
deck <- expand.grid(number=numbers, suit=suits)
( deck <- paste(deck$number, deck$suit) )
```

```
## [1] "Ace Diamonds" "Deuce Diamonds" "Three Diamonds" "Four Diamonds"
## [5] "Five Diamonds" "Six Diamonds" "Seven Diamonds" "Eight Diamonds"
## [9] "Nine Diamonds" "Ten Diamonds" "Jack Diamonds" "Queen Diamonds"
## [13] "King Diamonds" "Ace Clubs" "Deuce Clubs" "Three Clubs"
## [17] "Four Clubs" "Five Clubs" "Six Clubs" "Seven Clubs"
## [21] "Eight Clubs" "Nine Clubs" "Ten Clubs" "Jack Clubs"
## [25] "Queen Clubs" "King Clubs" "Ace Hearts" "Deuce Hearts"
## [29] "Three Hearts" "Four Hearts" "Five Hearts" "Six Hearts"
## [33] "Seven Hearts" "Eight Hearts" "Nine Hearts" "Ten Hearts"
## [37] "Jack Hearts" "Queen Hearts" "King Hearts" "Ace Spades"
## [41] "Deuce Spades" "Three Spades" "Four Spades" "Five Spades"
## [45] "Six Spades" "Seven Spades" "Eight Spades" "Nine Spades"
## [49] "Ten Spades" "Jack Spades" "Queen Spades" "King Spades"
```

```
## create a vector that contains the four ways we can get a king.
kings <- paste("King", suits)
## propoortion that are kings
mean(deck %in% kings)
```

`## [1] 0.07692`

We are going to be using the `combination()`

and `permutation()`

functions from the `gtools`

package.

The permutation functions computes for any list of size `n`

all the different ways we can select `R`

items.

`library(gtools)`

`permutations(5, 2)`

```
## [,1] [,2]
## [1,] 1 2
## [2,] 1 3
## [3,] 1 4
## [4,] 1 5
## [5,] 2 1
## [6,] 2 3
## [7,] 2 4
## [8,] 2 5
## [9,] 3 1
## [10,] 3 2
## [11,] 3 4
## [12,] 3 5
## [13,] 4 1
## [14,] 4 2
## [15,] 4 3
## [16,] 4 5
## [17,] 5 1
## [18,] 5 2
## [19,] 5 3
## [20,] 5 4
```

Optionally we can add a vector to this function that defines the range of digits. The code below samples 5 random 7 digit phone numbers out all possible phone numbers.

```
all_phone_numbers <- permutations(10, 7, v = 0:9)
n <- nrow(all_phone_numbers)
index <- sample(n ,5)
all_phone_numbers[index,]
```

```
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 6 0 3 7 2 4 1
## [2,] 7 2 1 8 4 3 9
## [3,] 3 6 9 0 1 4 7
## [4,] 2 8 3 9 6 7 4
## [5,] 8 9 7 4 3 0 6
```

To compute all possible ways that we can choose 2 cards when the **order** matters, we simply type the following piece of code.

```
hands <- permutations(52, 2, v=deck)
first_card <- hands[,1]
second_card <- hands[,2]
```

```
# check how many cases when first card is a king
sum(first_card %in% kings)
```

`## [1] 204`

```
# What fraction of these 204 also have a king as the second card
sum(first_card %in% kings & second_card %in% kings) /
sum(first_card %in% kings)
```

`## [1] 0.05882`

The above code is the `R`

version of the multiplication rule.

\[\Pr(A\mid B)=\frac{\Pr(A\ and\ B)}{\Pr(A)}\]

`permutations(3, 2)`

```
## [,1] [,2]
## [1,] 1 2
## [2,] 1 3
## [3,] 2 1
## [4,] 2 3
## [5,] 3 1
## [6,] 3 2
```

`combinations(3, 2)`

```
## [,1] [,2]
## [1,] 1 2
## [2,] 1 3
## [3,] 2 3
```

Now, what if the order does not matter? For example, in blackjack, if you get an ace and a face card or a 10, it’s called a natural 21, and you win automatically. If we want to compute the probability of this happening, we want to enumerate the combinations, not permutations, since the order doesn’t matter.

So if we get an **A** and a **king**, **king** and an **A**, it’s still a 21. We don’t want to count that twice. So notice the difference between the permutations functions, which lists all permutations, and the combination function, where order does not matter. This means that 2, 1 doesn’t appear because 1, 2 already appeared. Similarly, 3, 1 and 3, 2 don’t appear. So to compute the probability of a natural 21 in blackjack, we can do this.

```
aces <- paste("Ace", suits)
facecard <- c("King", "Queen", "Jack", "Ten")
facecard <- expand.grid(number=facecard, suit=suits)
facecard <- paste( facecard$number, facecard$suit)
hands <- combinations(52, 2, v = deck)
```

How often do you get aces and facecard

`mean(hands[,1] %in% aces & hands[,2] %in% facecard)`

`## [1] 0.04827`

`head(hands[,1])`

`## [1] "Ace Clubs" "Ace Clubs" "Ace Clubs" "Ace Clubs" "Ace Clubs" "Ace Clubs"`

We can also use a Monte Carlo to estimate this probably. In this case, we draw two cards over and over and keep track of how many 21’s we get. We can use the function sample to draw a card with a replacement like this.

```
# Here is one hand
(hand <- sample(deck, 2))
```

`## [1] "Jack Spades" "King Clubs"`

```
B <- 10000
results <- replicate(B, {
hand <- sample(deck, 2)
(hand[1] %in% aces & hand[2] %in% facecard) |
(hand[2] %in% aces & hand[1] %in% facecard)
})
mean(results)
```

`## [1] 0.0469`

## The Birthday Problem

Suppose you’re in a classroom with 50 people. If we assume this is a randomly selected group, what is the chance that at least two people have the same birthday? Although it is somewhat advanced, we can actually deduce this mathematically, and we do this later, but now, we’re going to use Monte Carlo simulations. For simplicity, we assumed that nobody was born on February 29th.

This actually doesn’t change the answer much. All right, first, note that birthdays can be represented as numbers between 1 and 365.

So a sample of 50 random birthdays can be obtained simply using the sample function, like this.

```
n <- 50
bdays <- sample(1:365, n, replace = TRUE)
```

To check if, in this particular set of 50 people we have at least two with the same birthday, we can use the function `duplicated`

, which returns true whenever an element of a vector has already appeared in that vector.

So here’s an example. If I type duplicated 1, 2, 3, 1, 4, 3, 5,

`duplicated(c(1,2,3,1,4,3,5))`

`## [1] FALSE FALSE FALSE TRUE FALSE TRUE FALSE`

And in the case that we just generated, we actually do have this happening. We did have two people, at least two people, with the same birthday. We get true. Now, to estimate the probability, we’re going to repeat this experiment. We’re going to run a Monte Carlo simulation over and over again. So what we do is we do it 10,000 times.

We use the replicate function like this.

```
B <- 10000
results <- replicate(B, {
bdays <- sample(1:365, n, replace = TRUE)
any(duplicated(bdays))
})
mean(results)
```

`## [1] 0.9706`

And when we’re done, we get that the probability of having two people, at least two people, with the same birthday in a group of 50 people is about 0%.

Were you expecting this to be so high? People tend to underestimate these probabilities, so it might be an opportunity for you to bet and make some money off of your friends. To get an intuition as to why this is so high, think of what happens as you get close to 365 people. At this stage, we run out of days, and the probability is 1. In the next video, we’re going to actually compute this probability for different group sizes and see how fast it gets close to 1.

## sapply

Say you want to use what you’ve just learned about the birthday problem to bet with friends about two people having the same birthday in a group of people. When are the chances larger than 50%? Larger than 75%? Let’s create a lookup table. We can quickly create a function to compute this for any group. We write the function like this. We’ll call it `compute_prob`

, and we’ll basically make the calculations for the probability of two people having the same birthday.

```
compute_prob <- function(n, B=10000){
same_day <- replicate(B, {
bdays <- sample(1:365, n, replace=TRUE)
any(duplicated(bdays))
})
mean(same_day)
}
```

Well, we use a small **Monte Carlo** simulation to do it. Now that we’ve done this, we want to compute this function, we want to apply this function to several values of n, let’s say from 1 to 60. Let’s define n as a sequence starting at 1 and ending at 60. Now, we can use a for loop to apply this function to each value in `n`

, but it turns out that for loops are rarely the preferred approach in R.

`n <- seq(1, 60)`

In general, we try to perform operations on entire vectors. Arithmetic operations, for example, operate on vectors in an element wise fashion. We saw this when we learned about `R`

. So if we type x equals 1 through 10, now X is the vector starting at 1 and ending at 10, and we compute the square root of x, it actually computes the square root for each element.

```
x <- 1:10
sqrt(x)
```

`## [1] 1.000 1.414 1.732 2.000 2.236 2.449 2.646 2.828 3.000 3.162`

Equally, if we define y to be 1 through 10, and then multiply x by y, it multiplies each element 1 by 1–

1 times 1, 2 times 2, et cetera.

```
y <- 1:10
x * y
```

`## [1] 1 4 9 16 25 36 49 64 81 100`

So there’s really no need for for loops. But not all functions work this way. You can’t just send a vector to any function in R. For example, the function we just wrote does not work element-wise since it’s expecting a scalar, it’s expecting an end.

This piece of code does not do what we want. If we type compute prob and send it the vector `n`

, we will not get what we want. We will get just one number.

`compute_prob(n)`

`## [1] 0`

What we can do instead is use the function `sapply`

. `sapply`

permits us to perform element-wise operations on any function. Here’s how it works. We’ll define a simple example for the vector 1 through 10. If we want to apply the square roots of each element of x, we can simply type sapply x comma square root, and it’ll apply square root to each element of x.

```
x <- 1:10
sapply(x, sqrt)
```

`## [1] 1.000 1.414 1.732 2.000 2.236 2.449 2.646 2.828 3.000 3.162`

We don’t need to do this because square root already does that, but we are using it as a simple example. So for our case, we can simply type prob equals sapply `n`

.

`n`

is our vector

`prob <- sapply(n, compute_prob)`

`plot(n, prob)`

and then the function we define `compute_prob`

. And this will assign to each element of prob the probability of two people having the same birthday for that end. And now we can very quickly make a plot. We plot the probability of two people having the same birthday against the size of the group.

Now, let’s compute the exact probabilities rather than use Monte Carlo simulations. The function we just define uses a Monte Carlo simulation, but we can use what we’ve learned about probability theory to compute the exact value. Not only do we get the exact answer using math, but the computations are much faster since we don’t have to generate experiments.

We simply get a number. To make the math simpler for this particular problem, instead of computing the probability of it happening, we’ll compute the probability of it not happening, and then we can use a multiplication rule. Let’s start with the first person. The probability that person 1 has a unique birthday is 1, of course.

Now let’s move on to the second one. The probability that the second person has a unique birthday given that person 1 already took one other day is 364 divided by 365. Then for a person 3, given that the first two people already have unique birthdays, that leaves 363. So now that probability is 363 divided by 365.

If we continue this way and find the chances of all, say, 50 people having unique birthdays, we would multiply 1 times 364 divided by 365, times 363 divided by 365, dot dot dot, all the way to the 50th element.

Here’s the equation: \[1\times\frac{364}{365}\times\frac{363}{365}\cdots\frac{365-n+1}{365}\]

Now, we can easily write a function that does this. This time we’ll call it `exact_prob`

. It takes `n`

as a argument, and it computes this probability using the simple code. Now we can compute each probably for each `n`

using `sapply`

again, like this.

```
exact_prob <- function(n){
prob_unique <- seq(365, 365 - n + 1) / 365
1 - prod(prob_unique)
}
eprob <- sapply(n, exact_prob)
```

And if we plot it, we can see that the Monte Carlo simulations were almost exactly right.

```
plot(n, prob)
lines(n, eprob, col="red")
```

They were almost exact approximations of the actual value. Now, notice had it not been possible to compute the exact probabilities, something that sometimes happens, we would have still been able to accurately estimate the probabilities using Monte Carlo.

## How many Monte Carlo experiments are enough?

In practice, we won’t know what the answer is, so we don’t know if our Monte Carlo estimate is accurate. We know that the larger the number of experiments, we’ve been using the letter B to represent that, the better the approximation. But how big do we need it to be? This is actually a challenging question, and answering it often requires advanced theoretical statistics training. One practical approach we will describe here is to check for the stability of the estimate.

Here’s an example with the birthday problem. We’re going to use `n <- 22`

. There’s 22 people. So we’re going to run a simulation where we compute or estimate the probability of two people having a certain birthday using different sizes of the Monte Carlo simulations. So the value of `B`

going to go from 10, to 20, to 40, to 100, et cetera. We compute the simulation, and now we look at the values that we get for each simulation.

```
B <- 10^seq(1, 5, len=100)
compute_prob <- function(B, n=22){
same_day <- replicate(B, {
bdays <- sample(1:365, n, replace = TRUE)
any(duplicated(bdays))
})
mean(same_day)
}
```

Remember, each simulation has a different B, a different number of experiments. When we see this graph, we can see that it’s wiggling up and down.

```
prob <- sapply(B, compute_prob)
plot(log10(B), prob, type = "l")
```

That’s because the estimate is not stable yet. It’s not such a great estimate. But as `B`

gets bigger and bigger, eventually it starts to stabilize. And that’s when we start getting a feeling for the fact that now perhaps we have a large enough number of experiments.

## The Addition Rule

The additional rule tells us that the probability of `A`

or `B`

, we’re going to have to make this calculation now, because you can get to 21 in two ways. You can get either a face card and then an ace. Or you can get an ace and then a face card. The additional rule tells us the probability of `A`

or `B`

is the probably of `A`

plus the probably `B`

minus the probability of `A`

and `B`

.

To understand why this makes sense, think of a Venn diagram. If we’re computing the probability of this whole thing happening, A or B, we can add the blue circle plus the pink circle, and then subtract the middle because we have added it twice by adding the blue plus the pink.

So now let’s apply it to the natural 21. In the case of natural 21, the intersection is empty. Since both hands can’t happen, you can’t have both an ace and then a face card, and then at the same time have a face card and then an ace. Those two things can’t happen at the same time. So this will be a very easy application of the addition rule. The probably of an ace followed by a face card we know is \(\frac{1}{13}\times\frac{16}{51}\). And the probability of a face card followed by an ace is \(\frac{16}{52}\times\frac{4}{51}\). These are actually the same, which makes sense to the symmetry. These two values are actually the same. In any case, we get the same result that we got before for the natural 21, which is about 0.05.

## The Monty Hall Problem

We’re going to look at one last example from the discrete probability. It’s the Monty Hall problem. In the 1970s, there was a game show called “Let’s Make a Deal”. Monty Hall was the host– this is where the name of the problem comes from. At some point in the game, contestants were asked to pick one of three doors. Behind one door, there was a prize. The other two had a goat behind them. And this basically meant that you lost.

If the contestant did not pick the prize door on his or her first try, Monty Hall would open one of the two remaining doors and show the contestant that there was no prize behind that door. So you’re left with two doors, the one you picked and one door that you do not know what’s behind it. So then Monty Hall would ask, do you want to switch doors? What would you do? We can use probability to show that if you stick to the original door, your chances of winning a car or whatever big prize is 1 in 3. But if you switch, your chances double to 2 in 3. This seems counterintuitive. Many people incorrectly think both chances are 1 in 2, since there’s only two doors left and there’s a prize behind one of them. You can find details, explanations– we’re going to provide links of the mathematics of how you can calculate that this is, in fact, wrong, that you have a higher chance if you switch. But here, we’re going to use a Monte Carlo simulation to show you that this is the case. And this will help you understand how it all works. Note that the code we’re going to show you is longer than it needs to be, but we’re using it as an illustration to help you understand the logic behind what actually happens here. So let’s start by creating a simulation that imitates the strategy of sticking to the same door.

Let’s go through the code. We’re going to do 10,000 experiments. So we’re going to do this over and over 10,000 times. First thing we do is we assign the prize to a door. So the prize is going to be in one of the three doors that we have created. Then that’s going to be in prize door. Prize door contains the number of the door with the prize behind it. Then we’re going to imitate your pick by taking a random sample of these three doors. Now, in the next step, we’re going to decide which door you’re shown. You’re going to be shown not your door and not the door with a prize. You’re going to be shown the other one. You stick to your door. That’s what you stick to. Nothing changed. All this that we did right above doesn’t matter, because you’re sticking to your door. So now we do this over and over again, and at the end, we ask is the door you chose, the one you stuck to, is that the prize door? How often does this happen? We know it’s going to be a 1/3, because none of this procedure changed anything. You started picking one out of three, nothing changed. And now, you are asked, is the one I picked the door? If we run the simulation, we actually see it happening. We ran it 10,000 times and the probability of winning was 0.3357, about 1 in 3.

```
# stick strategy code
B <- 10000
stick <- replicate(B, {
doors <- as.character(1:3)
prize <- sample(c('car', 'goat', 'goat'))
prize_door <- doors[prize == 'car']
my_pick <- sample(doors, 1)
show <- sample(doors[!doors %in% c(my_pick, prize_door)], 1)
stick <- my_pick
stick == prize_door
})
mean(stick)
```

`## [1] 0.3338`

```
# switch strategy code
switch <- replicate(B, {
doors <- as.character(1:3)
prize <- sample(c("car","goat","goat"))
prize_door <- doors[prize == "car"]
my_pick <- sample(doors, 1)
show <- sample(doors[!doors %in% c(my_pick, prize_door)], 1)
stick <- my_pick
switch <- doors[!doors%in%c(my_pick, show)]
switch == prize_door
})
mean(switch)
```

`## [1] 0.6675`

### Exercise 1. The Cavs and the Warriors

Two teams, say the Cavs and the Warriors, are playing a seven game championship series. The first to win four games wins the series. The teams are equally good, so they each have a 50-50 chance of winning each game.

If the Cavs lose the first game, what is the probability that they win the series?

- Assign the number of remaining games to the variable
`n`

. - Use the
`list`

function to create a list of game outcomes, where 0 indicates a loss for the Cavs and 1 indicates a win for the Cavs. Assign this value to the variable`l`

. - Use the
`expand.grid`

function to create a data frame containing all the possibilities for outcomes of the remaining games. - Use the
`rep`

function within the`expand.grid`

function to indicate the number of columns the results data frame should contain. - Use the
`rowSums`

function to identify which combinations of game outcomes result in the Cavs winning the number of games necessary to win the series. - Use the
`mean`

function to calculate the proportion of outcomes that result in the Cavs winning the series.

```
# Assign a variable 'n' as the number of remaining games.
n <- 6
# Assign a variable 'l' to a list of possible game outcomes, where 0 indicates a loss and 1 indicates a win for the Cavs.
l <- list(0:1)
# Create a data frame named 'possibilities' that contains all possible outcomes for the remaining games.
possibilities <- expand.grid(rep(l, n))
# Create a vector named 'results' that indicates whether each row in the data frame 'possibilities' contains enough wins for the Cavs to win the series.
results <- rowSums(possibilities)>=4
# Calculate the proportion of 'results' in which the Cavs win the series. Print the outcome to the console.
mean(results)
```

`## [1] 0.3438`

### Exercise 2. The Cavs and the Warriors - Monte Carlo

Confirm the results of the previous question with a Monte Carlo simulation to estimate the probability of the Cavs winning the series after losing the first game.

- Use the
`replicate`

function to replicate the sample code for`B <- 10000`

simulations. - Use the
`sampl`

e function to simulate a series of 6 games with random, independent outcomes of either a loss for the Cavs (0) or a win for the Cavs (1). - Use the
`sum`

function to indicate which of the simulated series contained at least 4 wins for the Cavs. - Use the
`mean`

function to find the proportion of simulations in which the Cavs win at least 4 of the remaining games.

```
# The variable `B` specifies the number of times we want the simulation to run. Let's run the Monte Carlo simulation 10,000 times.
B <- 10000
# Use the `set.seed` function to make sure your answer matches the expected result after random sampling.
set.seed(1)
# Create an object called `results` that replicates the sample code for `B` iterations and tallies the number of simulated series that contain at least four wins for the Cavs.
B <- 10000
results <- replicate(B, {
mean(sum(sample(0:1, 6, replace = TRUE)) >= 4)
})
# Calculate the frequency out of `B` iterations that the Cavs won at least four games in the remainder of the series. Print your answer to the console.
mean(results)
```

`## [1] 0.3453`

### Exercise 3. A and B play a series - part 1

Two teams, \(A\) and \(B\), are playing a seven series game series. Team \(A\) is better than team \(B\) and has a \(p > 0.5\) chance of winning each game.

- Use the function
`sapply`

to compute the probability, call it Pr of winning for`p <- seq(0.5, 0.95, 0.025)`

. - Then plot the result
`plot(p, prob)`

.

```
# Let's assign the variable 'p' as the vector of probabilities that team A will win.
p <- seq(0.5, 0.95, 0.025)
# Given a value 'p', the probability of winning the series for the underdog team B can be computed with the following function based on a Monte Carlo simulation:
prob_win <- function(p){
B <- 10000
result <- replicate(B, {
b_win <- sample(c(1,0), 7, replace = TRUE, prob = c(1-p, p))
sum(b_win)>=4
})
mean(result)
}
# Apply the 'prob_win' function across the vector of probabilities that team A will win to determine the probability that team B will win. Call this object 'Pr'.
Pr <- sapply(p, prob_win)
# Plot the probability 'p' on the x-axis and 'Pr' on the y-axis.
plot(p, Pr)
```

```
# Apply the 'prob_win' function across the vector of probabilities that team A will win to determine the probability that team B will win. Call this object 'Pr'.
Pr <- sapply(p, prob_win)
```

```
# Plot the probability 'p' on the x-axis and 'Pr' on the y-axis.
plot(p, Pr)
```

### Exercise 4. A and B play a series - part 2

Repeat the previous exercise, but now keep the probability that team \(A\) wins fixed at `p <- 0.75`

and compute the probability for different series lengths. For example, wins in best of 1 game, 3 games, 5 games, and so on through a series that lasts 25 games.

- Use the
`seq`

function to generate a list of odd numbers ranging from 1 to 25. - Use the function
`sapply`

to compute the probability, call it`Pr`

, of winning during series of different lengths. - Then plot the result
`plot(p, Pr)`

.

```
# Given a value 'p', the probability of winning the series for the underdog team B can be computed with the following function based on a Monte Carlo simulation:
prob_win <- function(N, p=0.75){
B <- 10000
result <- replicate(B, {
b_win <- sample(c(1,0), N, replace = TRUE, prob = c(1-p, p))
sum(b_win)>=(N+1)/2
})
mean(result)
}
# Assign the variable 'N' as the vector of series lengths. Use only odd numbers ranging from 1 to 25 games.
N <- seq(1, 25, 2)
# Apply the 'prob_win' function across the vector of series lengths to determine the probability that team B will win. Call this object `Pr`.
Pr <- sapply(N, prob_win)
# Plot the number of games in the series 'N' on the x-axis and 'Pr' on the y-axis.
plot(N, Pr)
```

`sessionInfo()`

```
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.1 LTS
##
## Matrix products: default
## BLAS: /home/michael/anaconda3/lib/R/lib/libRblas.so
## LAPACK: /home/michael/anaconda3/lib/R/lib/libRlapack.so
##
## locale:
## [1] en_CA.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] gtools_3.8.1 RevoUtils_11.0.1 RevoUtilsMath_11.0.0
##
## loaded via a namespace (and not attached):
## [1] bookdown_0.7 Rcpp_0.12.18 codetools_0.2-15 digest_0.6.15
## [5] rprojroot_1.3-2 backports_1.1.2 magrittr_1.5 evaluate_0.11
## [9] blogdown_0.9.8 stringi_1.2.4 rmarkdown_1.10 tools_3.5.1
## [13] stringr_1.3.1 xfun_0.4.11 yaml_2.2.0 compiler_3.5.1
## [17] htmltools_0.3.6 knitr_1.20
```

```
## Adding cites for R packages using knitr
knitr::write_bib(.packages(), "packages.bib")
```

R Core Team. 2018. *R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.