In order to run simulations of experiments with random outcomes we make use of R's ability to generate random numbers. More precisely it generates a 'pseudo-random' number, which behaves in many ways like a truly random number.
# sample(x,k) generates a random permutation of k objects from the vector x. That is, all k choices are different # start with a vector x with 5 elements. > x = 1:5 > print(x) [1] 1 2 3 4 5 # randomly sample 3 of the elements of x. > y = sample(x, 3) > print(y) [1] 2 5 1 # Your answer may be different. > y = sample(x, 3) > print(y) [1] 3 1 4 # Your answer may be different. > sample(x, 3) > print(y) [1] 1 2 5 # Your answer may be different. # Note every element is chosen at most once. # randomly sampling 5 elements of x is a permutation of all 5 elements. > y = sample(x, 5) > print(y) [1] 1 3 2 5 4 # Your answer may be different. > y = sample(x, 5) > print(y) [1] 2 3 1 4 5 # Your answer may be different. # You can't pick more than 5 different elements from a set of 5 things, so R gives an error. > sample(x, 6) Error in sample.int(length(x), size, replace, prob) : cannot take a sample larger than the population when 'replace = FALSE' Allowing repeated elements # Sometimes you want to allow repeats. For example when we roll a die repeatedly we expect to see numbers repeating. We can think of this as picking a random element from a set and then putting it back, i.e. replacing it, so it can be drawn again. # In R we do this by setting the optional argument replace=TRUE # Note that we get can get repeated values > y = sample(x, 3, replace=TRUE) > print(y) [1] 3 1 4 # Your answer may be different. > y = sample(x, 3, replace=TRUE) > print(y) [1] 4 5 4 # Your answer may be different. > y = sample(x, 3, replace=TRUE) > print(y) [1] 2 1 5 # Your answer may be different. > y = sample(x, 5, replace=TRUE) > print(y) [1] 5 1 5 1 4 # Your answer may be different. > y = sample(x, 5, replace=TRUE) > print(y) [1] 3 5 5 2 1 # Your answer may be different. > y = sample(x, 5, replace=TRUE) > print(y) [1] 4 2 3 2 1 # Your answer may be different. # Now there is no problem asking for more than 5 things from a set of 5 elements. > y = sample(x, 12, replace=TRUE) > print(y) [1] 1 1 2 1 4 2 3 4 1 2 4 5 # Your answer may be different. > y = sample(x, 12, replace=TRUE) > print(y) [1] 3 2 5 1 4 2 4 1 4 5 3 1 # Your answer may be different. # To generate an m x n array of random values we can use the function sample function followed by the matrix function. # Let's simulate rolling a die. # We use 1:6 to make vector (1,2,3,4,5,6). # Generate a 3 x 4 array of random dice rolls. > y = sample(1:6, 12, replace=TRUE) > z = matrix(y, nrow=3, ncol=4) > print(z) [,1] [,2] [,3] [,4] [1,] 1 3 6 1 [2,] 2 2 3 6 [3,] 2 1 4 3 # Or we could make it a 2 x 6 array. > z = matrix(y, nrow=2, ncol=6) > print(z) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 2 2 6 4 6 [2,] 2 3 1 3 1 3 # Your answer may be different.
# First let's simulate rolling 1 die 3 times and checking if one of the rolls was a 6. # First we use sample to generate 3 samples from 1:6 > x = sample(1:6, 3, replace=TRUE) > print(x) [1] 6 4 1 # Your answer may be different. # To check if any of the 3 rolls are 6 we use the following command. It returns a vector of TRUE or FALSE depending on whether that entry of x is 6 or not. Note the use of the double equal sign. We can't use a single equal sign because that would mean 'set the value of x to 6'. Compare the result with the value of x above. > print(x == 6) [1] TRUE FALSE FALSE # Your answer may be different. # We can also see which elements of x are less than 6 > print(x < 6) [1] FALSE TRUE TRUE # Your answer may be different. # Now let's roll the die 1000 times and see what fraction of the rolls give 6. We expect that about 1/6 of them will be. # Simulate 1000 rolls > x = sample(1:6, 1000, replace=TRUE) # Your answer may be different. # x == 6 gives a vector of TRUE or FALSE. R is clever, when we sum the vector: each TRUE counts as 1 and each FALSE counts as 0. So the sum is the number of TRUE's. In this case that means the number of 6's, which happens to be 168. > s = sum(x == 6) > print(s) [1] 168 # Your answer may be different. # Divide by 1000 to get the fraction of 6's. > a = sum(x == 6)/1000 > print(a) [1] 0.168 # Your answer may be different. # Compare that with the theoretical value of 1/6. > print(1/6) [1] 0.1666667 # Not bad! # Now let's estimate the probability of getting at least one 6 in 4 rolls. # Goal: estimate the probability of getting at least one 6 in 4 rolls. # Experiment: roll 1 die 4 times and check for a 6. # Repeat the experiment 10 times and see what fraction of times this happens. # So you can see all the commands at once, we'll show them all and then explain them later. For commenting, we'll put a command number before each '>' 1. > x = matrix(sample(1:6, 4*10, replace=TRUE), nrow=4, ncol=10) 2. > print(x) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 5 6 1 1 3 6 1 4 6 1 [2,] 3 6 6 4 3 2 2 1 6 5 [3,] 2 3 5 2 5 3 4 6 2 5 [4,] 2 6 6 5 5 3 1 4 5 3 # Your answer may be different. 3. > y = (x==6) 4. > print(y) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE [2,] FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE [3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE [4,] FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE # Your answer may be different. 5. > z = colSums(y) 6. > print(z) [1] 0 3 2 0 0 1 0 1 2 0 # Your answer may be different. 7. > print(z > 0) [1] FALSE TRUE TRUE FALSE FALSE TRUE FALSE TRUE TRUE FALSE 8. > s = sum(z > 0) 9. > print(s) [1] 5 # Your answer may be different. 10. > m = mean(z > 0) 11. > print(m) [1] 0.5 # Your answer may be different. # Command 1: Generate a 4 by 10 random array. Each column represents one experimental trial of 4 rolls. The 10 columns represent the 10 trials. # Command 2: Display x # Command 3: See which of the rolls are 6's. Note the double equals to check which rolls were 6. The result y is an array of TRUE or FALSE. You can check that it picks out the 6's. # Command 4: Display y # Command 5: Sum each of the columns. The result is then number of 6's in each trial (column). # Command 6: Display z # Command 7: If the column sum is greater than 0 there was at least one TRUE in the column, that is at least one 6. This command prints a vector of TRUE or FALSE representing whether or not the experiment yielded a 6. # command 8: sum the vector in command 7 to get the number of trials that yielded a 6. We see that 5 out of 10 trials did. (This is random. Your answer may be different.) # Command 9: Display s # command 10: The mean function is just the sum divided by the number of trials. This is just 5/10 = .5. Half the trials yielded a 10. # command 11: Disply the mean # Let's repeat this but with 1000 trials instead of 10. Without all the comments it's pretty short. > x = matrix(sample(1:6, 4*1000, replace=TRUE), nrow=4, ncol=1000) > y = (x==6) > z = colSums(y) > print(sum(z > 0)) [1] 525 > print(mean(z>0)) [1] 0.525 Our estimate of the probability of at least one 6 in 4 rolls is .525. This is pretty close to the theoretical value of .518. The dim() function # We can always check that x has 4 rows and 1000 columns using the dim() function. > print(dim(x)) [1] 4 1000
# Goal: estimate the probability of getting a sum of 7 when rolling two dice. 1. > ntrials = 10000 2. > x = matrix(sample(1:6, 2*ntrials, replace=TRUE), nrow=2, ncol=ntrials) 3. > y = colSums(x) 4. > print(mean(y == 7)) [1] 0.1658 # Command 1: We assign the number of trials to the variable ntrials. Writing code like this makes it easier to modify later. If we want to change the number of trials to 7 we just have to change this one line of code. # Command 2: we create 10000 columns with 2 rows. That is, we run 10000 experiments of rolling 2 dice. # Command 3: we sum each of the columns, i.e., we sum the two dice. # Command 4: we find the fraction of sums that are 7. # Note, this is essentially the exact answer of 1/6. Exercise: Try to estimate the probability of two sixes when rolling two dice.
18.05 Introduction to Probability and Statistics
Spring 2022
Author: Jeremy Orloff
For information about citing these materials or our Terms of Use, visit: https://ocw.mit.edu/terms