Experiments, Conditional Probability, and Bayes’ Theorem
Madeline Moldrem 118800318
- 1 Basic R manipulation.
- 2 Simulating Experiments
- 3 Experiments involving Coin Tosses and Die Rolls
- 4 Bayes’s Theorem and Witnesses
1 Basic R manipulation.
- In the code segment below, assign the value 2 to variable x, and the value 6 to variable y. Then compute the sum ‘x+y’:
x <- 2
y <- 6
print(x + y)## [1] 8- Construct a list of all even numbers between 1 and 100, assign this
to the variable ‘even’. Construct a list of 100 numbers, starting with
1, such that the difference between consecutive numbers is three, assign
this to the variable ‘threes’. Construct a list of 100 equally spaced
numbers between 0 and 2. Assign this to the variable
even = seq(from = 2, to = 100, by = 2)
threes = seq(from = 1, by = 3, length.out = 100)
x = seq(from = 0, to = 2, length.out = 100)- Define a function that computes the square of the input value, call this function square:
square <- function(a){
a * a
} - Plot trigonometric functions sine and cosine for values of x between 0 and 10. Also plot the ‘square’ function for input values between 0 and 2.
square <- function(a){
a * a
}
square(5)## [1] 25x = seq(0, 10, length.out = 100)
sin = sin(x)
cos = cos(x)
plot(x, sin, type = "l", main = "sin(x) function", xlab = "x",ylab = "sin(x)")
plot(x, cos, type = "l", main = "cos(x) function", xlab = "x",ylab = "cos(x)")
x = seq(0, 2, length.out = 100)
squarefunct = square(x)
plot(x, squarefunct, type = "l", main = "square function", xlab = "x",ylab = "x^2")
2 Simulating Experiments
In this problem we will simulate two important experiments: coin tosses and die rolls.
We will explore empirical probability of events using random samples. The two important functions to understand are “sample” and “replicate”.
We use the sample function to set up our experiment, and replicate function to repeat it as many times as we want.
2.1 Functions for Coin Toss and Die Rolls
Set up an function, called ‘coin_toss’, that simulates a single toss of a coin with \(P(H) = p\in (0,1)\). Note, the function will take input \(p\), and will output the outcome of a single toss of the coin. It might help to have the function return 1 if coin lands H, and return 0 otherwise.
coin <- c("heads", "tails")
coin_toss <- function(p){
sample(coin, size = 1, replace = TRUE, prob = c(p, 1-p))
}Set up an function, called ‘die_roll’, that simulates a single roll of a six-sided die, with a given probability distribution. Note, the function will take input a vector \(p\) with the six probabilities associated the respective faces of the die, and will output the outcome of a single roll of the die.
die_roll <- function(probs){
die <- 1:6
sample(die, size = 1, prob = probs)
}3 Experiments involving Coin Tosses and Die Rolls
3.1 Die Rolls
Setup the experiment where we roll two fair six 6-sided die independently, and calculate the sum of the numbers that appear for each of them.
probabilities <- c(1/6, 1/6, 1/6, 1/6, 1/6, 1/6)
roll1 <- die_roll(probabilities)
roll2 <- die_roll(probabilities)
sum_rolls <- roll1 + roll2Use the replicate function to repeat this experiment 10000 times, store the output of the experiment in a variable “sum_die”.
sum_die <- replicate(10000, die_roll(probabilities) + die_roll(probabilities))Use “sum_die” to calculate the empirical probability of getting the sums 2, 3, 4, …, 12. List these empirical probabilities into a table, using the ‘table’ function.
counts = table(sum_die)
empirical <- counts / sum(counts)Compare your answers with the theoretical values. P(sum = 2) Theoretical probability: 1/36, or 0.0278 Experiment result: 0.0293
P(sum = 3) Theoretical probability: 2/36, or 0.0556 #Experiment result: 0.0560
#P(sum = 4) #Theoretical probability: 3/36, or 0.0833 #Experiment result: 0.0814
#P(sum = 5) #Theoretical probability: 4/36, or 0.1111 #Experiment result: 0.1063
#P(sum = 6) #Theoretical probability: 5/36, or 0.1389 #Experiment result: 0.1329
#P(sum = 7) #Theoretical probability: 6/36, or 0.1667 #Experiment result: 0.1707
#P(sum = 8) #Theoretical probability: 5/36, or 0.1389 #Experiment result: 0.1426
#P(sum = 9) #Theoretical probability: 4/36, or 0.1111 #Experiment result: 0.1159
#P(sum = 10) #Theoretical probability: 3/36, or 0.0833 #Experiment result: 0.0827
#P(sum = 11) #Theoretical probability: 2/36, or 0.0556 #Experiment result: 0.0546
#P(sum = 11) #Theoretical probability: 1/36, or 0.0278 #Experiment result: 0.0276
Plot the following: relative frequency of getting sum equal to 8 as a function of number of trials, for this experiment.
trials <- 1:10000
eights <- cumsum(sum_die == 8) / trials
plot(trials, eights, type = "l", xlab = "trials", ylab = "frequency of 8", main = "Getting a sum of 8")
What do you notice about this graph as the number of trials increase? YOUR ANSWER: The theoretical probability of getting a sum of 8 is 0.1389, and as the number of trials increases, the empirical probability gets closer to the theoretical probability.
3.2 Calculating a Die Roll Probability
Suppose an experiment involves rolling a 6-sided die 5 times, the die has the following distribution: \[P(1) = P(2) = P(3) = P(4)= P(5)=0.1, \quad P(6) = 0.5\]
Suppose \[E\] is the event that the sum of all the die rolls is 13. And \[G\] is the event that the max of the 5 rolls is atmost 4.
Set up a function ‘E_happens’, that rolls the die 5 times and returns 1 if E happens, and returns zero otherwise
E_happens <- function(){
results <- sample(1:6, size = 5, replace = TRUE, prob = c(0.1, 0.1, 0.1, 0.1, 0.1, 0.5))
if(sum(results) == 13){
return (1)
} else {
return (0)
}
}Set up a function ‘G_happens’, that rolls the die 5 times and returns 1 if G happens, and returns 0 otherwise
G_happens <- function(){
results <- sample(1:6, size = 5, replace = TRUE, prob = c(0.1, 0.1, 0.1, 0.1, 0.1, 0.5))
if(max(results) < 5){
return (1)
} else {
return (0)
}
}Use the two functions you constructed, ‘E_happens’ and ‘G_happens’, and the ‘replicate’ function, to calculate the empirical approximations to \(P(E)\) and \(P(G)\).
E_results <- replicate(10000, E_happens())
counts <- table(E_results)
P_E <- counts / sum(counts)
G_results <- replicate(10000, G_happens())
counts <- table(G_results)
P_G <- counts / sum(counts)What are the approximate probabilities? YOUR ANSWER: P(E) = 0.0077 P(G) = 0.0111
3.3 Coin Tosses
Suppose we toss a coin with \(P(H) = 0.8\), \(n\) times independently. Let \(X\) be the the random variable that counts the number of heads in \(n\) tosses.
3.3.1 n = 2
Setup a function that simulates the random variable \(X\) when we toss the coin two times.
coin <- c(0, 1)
P_H = 0.8
toss_2 <- function(p){
X <- 0
for(i in 1:2){
toss <- sample(coin, size = 1, replace = TRUE, prob = c(p, 1 - p))
if(toss == 0){
X = X + 1
}
}
return (X)
}Use the above function to calculate the emprical probability distribution table for \(X\). Compare this with the theoretical probability distribution of \(X\).
exp <- replicate(1000, toss_2(P_H))
empirical <- table(exp) / 1000
#**P(X = 0)**
#**Theoretical probability: (0.2)(0.2) = 0.04**
#**Experiment result: 0.036**
#**P(X = 1)**
#**Theoretical probability: (0.2)(0.8) + (0.2)(0.8) = 0.32**
#**Experiment result: 0.322**
#**P(X = 2)**
#**Theoretical probability: (0.8)(0.8) = 0.64**
#**Experiment result: 0.642**Calculate empirical estimates for the expected value and the variance of \(X\). Compare these with the actual theoretical calculations.
E_X_empirical <- 0*0.036 + 1*0.322 + 2*0.642
E_X_the <- 2*(0.8)
#**The expected theoretical E(X) is 1.6, and the empirical E(X) is 1.606, closely modeling the theoretical value.**
V_X_empirical <- ((0*0.036 + 1*0.322 + 4*0.642) - (1.606 * 1.606))
V_X_the <- 2*(0.8)*(0.2)
#**The expected theoretical V(X) is 0.32, and the empirical V(X) is 0.310764, closely modeling the theoretical value.**3.3.2 n = 3
Setup a function that simulates the random variable \(X\) when we toss the coin three times.
coin <- c(0, 1)
P_H = 0.8
toss_3 <- function(p){
X <- 0
for(i in 1:3){
toss <- sample(coin, size = 1, replace = TRUE, prob = c(p, 1 - p))
if(toss == 0){
X = X + 1 #increments every time we get a 0, or "heads"
}
}
return (X)
}Use the above function to calculate the emprical probability distribution table for \(X\). Compare this with the theoretical probability distribution of \(X\).
exp <- replicate(1000, toss_3(P_H))
empirical <- table(exp) / 1000
#**P(X = 0)**
#**Theoretical probability: (0.2)(0.2)(0.2) = 0.008**
#**Experiment result: 0.009**
#**P(X = 1)**
#**Theoretical probability: 3(0.2)(0.2)(0.8) = 0.096**
#**Experiment result: 0.105**
#**P(X = 2)**
#**Theoretical probability: 3(0.8)(0.8)(0.2) = 0.384**
#**Experiment result: 0.382**
#**P(X = 3)**
#**Theoretical probability: (0.8)(0.8)(0.8) = 0.512**
#**Experiment result: 0.504**Calculate emprical estimates for the expected value and the variance of \(X\). Compre these with the actual theoretical calculations.
E_X_empirical <- 0*0.009 + 1*0.105 + 2*0.382 + 3*0.512
E_X_the <- 3*(0.8)
#**The expected theoretical E(X) is 2.4, and the empirical E(X) is 2.405, closely modeling the theoretical value.**
V_X_empirical <- ((0*0.009 + 1*0.105 + 4*0.382 + 9*0.512) - (2.405 * 2.405))
V_X_the <- 3*(0.8)*(0.2)
#**The expected theoretical V(X) is 0.48, and the empirical V(X) is 0.456975, closely modeling the theoretical value.**3.3.3 n = 4
Setup a function that simulates the random variable \(X\) when we toss the coin four times.
coin <- c(0, 1)
P_H = 0.8
toss_4 <- function(p){
X <- 0
for(i in 1:4){
toss <- sample(coin, size = 1, replace = TRUE, prob = c(p, 1 - p))
if(toss == 0){
X = X + 1 #increments when we get a 0
}
}
return (X)
}Use the above function to calculate the empirical probability distribution table for \(X\). Compare this with the theoretical probability distribution of \(X\).
exp <- replicate(1000, toss_4(P_H))
empirical <- table(exp) / 1000
#**P(X = 0)**
#**Theoretical probability: (0.2)(0.2)(0.2)(0.2) = 0.0016**
#**Experiment result: 0.000**
#**P(X = 1)**
#**Theoretical probability: 4(0.2)(0.2)(0.2)(0.8) = 0.0256**
#**Experiment result: 0.022**
#**P(X = 2)**
#**Theoretical probability: 6(0.2)(0.2)(0.8)(0.8) = 0.1536**
#**Experiment result: 0.151**
#**P(X = 3)**
#**Theoretical probability: 4(0.2)(0.8)(0.8)(0.8) = 0.4096**
#**Experiment result: 0.420**
#**P(X = 4)**
#**Theoretical probability: (0.8)(0.8)(0.8)(0.8) = 0.4096**
#**Experiment result: 0.407**Calculate emprical estimates for the expected value and the variance of \(X\). Compre these with the actual theoretical calculations.
E_X_empirical <- 0*0.000 + 1*0.022 + 2*0.151 + 3*0.420 + 4*0.407
E_X_the <- 4*(0.8)
#**The expected theoretical E(X) is 3.2, and the empirical E(X) is 3.212, closely modeling the theoretical value.**
V_X_empirical <- ((0*0.000 + 1*0.022 + 4*0.151 + 9*0.420 + 16*0.407) - (3.212 * 3.212))
V_X_the <- 4*(0.8)*(0.2)
#**The expected theoretical V(X) is 0.64, and the empirical V(X) is 0.601056, closely modeling the theoretical value.**3.4 Coin Tosses for with Infinite Sample Space.
Suppose we toss a coin with \(P(H) = 0.8\) until we get \(r\) heads.
3.4.1 r=1
Setup an experiment that counts the number of tails until the first head. You will need to use the ‘while’ loop.
Run 10000 simulations of this experiment, and use these to find estimate for the expected number of tails until the first head. Match this value with the theoretical/calculated value.
P_H = 0.8
coin <- c(0, 1)
toss_until_H <- function() {
a <- 2
num_tails <- 0
while (a == 2){
toss <- sample(coin, size = 1, replace = TRUE, prob = c(P_H, 1 - P_H))
if(toss == 0){ #first head!
return (num_tails)
} else {
num_tails <- num_tails + 1 #increment number of tails if toss is 1
}
}
}
toss_until_H_results <- replicate(10000, toss_until_H())
counts = table(toss_until_H_results)
empirical <- counts / sum(counts)
E_X_empirical = 0*0.7998 + 1*0.1581 + 2*0.0327 + 3*0.0068 + 4*0.0023 + 5*0.0003
E_X_the = (0.2)/(0.8)
#**The expected theoretical E(X) is 0.25, and the empirical E(X) is 0.2546, closely modeling the theoretical value.**3.4.2 r = 2
Setup an experiment that counts the number of tails until the second head. Run 10000 simulations of this experiment, and use these to find an estimate for the expected number of tails until the second head. Match this value with the theoretical/calculated value.
toss_until_2H <- function() {
num_heads <- 0
num_tails <- 0
while (num_heads < 2){ #repeat until second head
toss <- sample(coin, size = 1, replace = TRUE, prob = c(P_H, 1 - P_H))
if(toss == 0){
num_heads = num_heads + 1
} else {
num_tails <- num_tails + 1
}
}
return (num_tails)
}
toss_until_2H_results <- replicate(10000, toss_until_2H())
counts = table(toss_until_2H_results)
empirical <- counts / sum(counts)
E_X_empirical = 0*0.6424 + 1*0.2515 + 2*0.0786 + 3*0.0210 + 4*0.0048 + 5*0.0014 + 6*0.0003
E_X_the = (2*0.2)/(0.8)
#**The expected theoretical E(X) is 0.5, and the empirical E(X) is 0.4997, closely modeling the theoretical value.**3.4.3 r = 3
Setup an experiment that counts the number of tails until the third head. Run 10000 simulations of this experiment, and use these to find an estimate for the expected number of tails until the third head. Match this value with the theoretical/calculated value.
toss_until_3H <- function() {
num_heads <- 0
num_tails <- 0
while (num_heads < 3){ #repeat until third head
toss <- sample(coin, size = 1, replace = TRUE, prob = c(P_H, 1 - P_H))
if(toss == 0 ){
num_heads = num_heads + 1
} else {
num_tails = num_tails + 1
}
}
return (num_tails)
}
toss_until_3H_results <- replicate(10000, toss_until_3H())
counts = table(toss_until_3H_results)
empirical <- counts / sum(counts)
E_X_empirical = 0*0.8055 + 1*0.1538 + 2*0.0337 + 3*0.0058 + 4*0.0011 + 5*0 + 6*0.0001
E_X_the = (3*0.2)/(0.8)
#**The expected theoretical E(X) is 0.75, and the empirical E(X) is 0.2436**4 Bayes’s Theorem and Witnesses
Recall the hit-and-run example from the asynchronous lecture and HW4. We noted that the Bayes’ probabilities are calculated as the following formula \[ P(B|W) = \frac{pq}{pq + (1-p)(1-q)},\] where \(P(B) = q\) and \(P(W|B) = p\) (check notes). We also noted in HW4 that if \(P(W|B)=p > \frac{1}{2}\) we will have \[ P(B|W) > P(B).\]
In this problem, we will explore how changing the level of Witness-Reliability, that is \(p\), will affect the chance of getting a Blue Cab driver arrested.
Recall that \(P(B)\) is the probability that the Blue cab is at fault with no added information. We will now run multiple iteration of Bayes’ rule, and at every iteration, we will update the value of \(P(B)\) with the value of \(P(B|W)\) from earlier iteration. From our experience from the homework, if \(p > \frac{1}{2}\), we must have \(P(B) \longrightarrow 1\) as the number of iterations increase (intuitively, the probability that the Blue cab is at fault will go to 1 as the number of “reliable witness” who say that “the Blue Cab is at fault” increases).
Write a function called “bayes” that takes input \(q, p\) and outputs \(P(B|W)\) (use the formula above)
bayes <- function(q, p){
return(p*q/(p*q + (1-p)*(1-q)))
}
bayes(0.01, 0.8)## [1] 0.038834954.1 Reliable Witnesses with fixed probability
In this part, we will assume that we are using testimony of reliable witnesses with fixed reliability \(p\).
4.1.1 Initial Plots
Initialize \(p=0.6\), \(q=0.01\), and q_vals = c(), the empty vector. Run a for loop (1 in 1:40) inside it, update q at every iteration to \(q = \text{bayes}(p,q)\), and concatenate the vector q_vals with this new value of q.
p <- 0.6
q <- 0.01
q_vals <- c()
for (i in 1:40) {
q <- bayes(q, p)
q_vals <- c(q_vals, q)
}Construct a plot with x-axis 1:40 and y-axis q_vals.
plot(1:40, q_vals, type = "l", main = "Bayes after 40 updates", xlab = "x",ylab = "q_vals")
Repeat the above experiment with \(p=0.4\), \(q=0.99\).
p <- 0.4
q <- 0.99
q_vals <- c()
for (i in 1:40) {
q <- bayes(q, p)
q_vals <- c(q_vals, q)
}
plot(1:40, q_vals, type = "l", main = "Bayes after 40 updates", xlab = "x",ylab = "q_vals")
4.1.2 How many reliable Witnesses needed?
Suppose the police decide to arrest the Blue cab driver if \(P(B|W) \ge 0.9\).
Suppose we intialize \(p=0.6\) (witness reliability), \(q=0.01\) (proportion of Blue cabs in the city), how many witnesses with reliability level \(p\) would we need for the police to consider arresting the Blue cab?
You will need a counter variable initialized to zero and a while loop inside which you update \(q\) using the bayes function, and update the counter by 1 at every iteration of the while loop.
a <- 0
q <- 0.1
p <- 0.6
while(q < 0.9){
q <- bayes(q, p)
a = a + 1
}
#It takes 11 witness4.2 Working with two types of witnesses
In this problem we will work with two types of witnesses (1) reliable witnesses with fixed reliability level \(p_1 > \frac{1}{2}\) and (2) unreliable witness with fixed reliability level \(p_2 < \frac{1}{2}\).
4.2.1 Initial Plots
Initialize \(p_1=0.6\), \(p_2 = 0.45\), \(q=0.01\), and q_vals = c(), the empty vector. Run a for loop (1 in 1:90) inside it, update q at every even iteration to \(q = \text{bayes}(p_2,q)\), and at every odd iteration to \(q = \text{bayes}(p_1,q)\) (You will need to use a “if” command), and finally concatanating the vector q_vals with the new value of q at every iteration.
p1 <- 0.6
p2 <- 0.45
q <- 0.01
q_vals <- c()
for(i in 1:90){
if(i %% 2 == 0){ #even number
q <- bayes(q, p2)
} else {
q <- bayes(q, p1)
}
q_vals <- c(q_vals, q)
}Construct a plot with x-axis 1:90 and y-axis q_vals.
plot(1:90, q_vals, type = "l", main = "Bayes after 90 updates w 2 witnesses", xlab = "x",ylab = "q_vals")
What do you notice? (Compared to case with single witness) #The graph is more jagged, and although it follows a similar trend upward as the single witness example, it oscillates up and down.
4.2.2 How many witnesses needed?
Suppose the police decide to arrest the Blue cab driver if \(P(B|W) \ge 0.9\).
Initialize \(p_1=0.6\), \(p_2 = 0.45\), \(q=0.01\). How many witnesses will we have to have if we are going to sequentially work with one reliable witness followed by an unreliable witness?
p1 <- 0.6
p2 <- 0.45
q <- 0.01
a <- 0
b <- 1
while(q < .9){
if(b %% 2 == 0){ #even number
q <- bayes(q, p2)
} else {
q <- bayes(q, p1)
}
b = b + 1
a = a + 1
}
print(a)## [1] 65#It took 65 witnessesThis problem can be generalized to choosing \(p\) from any given probability distribution on \((0,1)\), we will be working with the Beta distribution later in the semester, which can be an interesting candidate to sample our values for \(p\) (that is witness reliability).
4.3 Working with four types of witnesses
In this problem we will work with two types of witnesses (1) two reliable witnesses with fixed reliability level \(p_1, p_2 > \frac{1}{2}\) and (2) two unreliable witnesses with fixed reliability level \(p_3, p_4 < \frac{1}{2}\).
4.3.1 Initial Plots
Initialize \(q=0.01\) and we will sample \(p\) from the vector c(0.6, 0.8, 0.45, 0.3) and suppose choosing a witness with either of these witness-reliabilities is equally likely.
Now let q_vals = c(), the empty vector and run a for loop (1 in 1:100) inside it, update q by sampling p from the vector c(0.6, 0.8, 0.45, 0.3), and concatanating the vector q_vals with the new value of q at every iteration.
q <- 0.01
p_vals <- c(0.6, 0.8, 0.45, 0.3)
q_vals <- c()
for(i in 1:100){
p <- sample(p_vals, 1)
q <- bayes(q, p)
q_vals <- c(q_vals, q)
}Construct a plot with x-axis 1:100 and y-axis q_vals.
plot(1:100, q_vals, type = "l", main = "Bayes after 100 updates w 4 witnesses", xlab = "x",ylab = "q_vals")
What do you notice? How does this compare with the earlier cases? #It follows a similar upward trend as as the earlier cases, but the curve is much less consistent. It’s not smooth #like the one-type witness example, or consistently jagged like the two-type example that alternated neatly. #Because the witness reliability is somewhat randomized, it likely not be as smooth as the previous examples.
4.3.2 How many witnesses needed?
Suppose the police decide to arrest the Blue cab driver if \(P(B|W) \ge 0.9\).
Initialize \(q=0.01\) and we will sample \(p\) from the vector c(0.6, 0.8, 0.45, 0.3) and suppose choosing a witness with either of these witness-reliabilities is equally likely. How many witnesses will we need to have if we are going to uniformly sample from each of these four types of of witnesses at every step?
q <- 0.01
a <- 0
b <- 1
while(q < 0.9){
if(b %% 4 == 1){ #0.6 reliability every four witnesses
q <- bayes(q, 0.6)
} else if (b %% 4 == 2){ #0.8 reliability every four witnesses
q <- bayes(q, 0.8)
} else if (b %% 4 == 3){#0.45 reliability every four witnesses
q <- bayes(q, 0.45)
} else {#0.3 reliability every four witnesses
q <- bayes(q, 0.3)
}
a = a + 1
b = b + 1
}
#We needed 30 witnesses4.4 BONUS: Working with witnesses coming from a Beta distribution
In this problem we will assume that the witness reliability of the people in the city follows the Beta distribution with parameters \(\alpha\), \(\beta\). Choose \(\alpha, \beta\) in such a way that the expected value of this Beta distribution is greater than .65.
4.4.1 Initial Plots
Initialize \(q=0.01\) and we will choose \(p\) by sampling the Beta distribution with parameters \(\alpha, \beta\) chosen above.
Now let q_vals = c(), the empty vector and run a for loop (1 in 1:100) inside it, update q by sampling p \(\text{Beta}(\alpha, \beta)\), and concatanating the vector q_vals with the new value of q at every iteration.
Construct a plot with x-axis 1:100 and y-axis q_vals.
What do you notice? How does this compare with the earlier cases?
4.4.2 How many witnesses needed?
Suppose the police decide to arrest the Blue cab driver if \(P(B|W) \ge 0.9\).
Initialize \(q=0.01\) and we will choose \(p\) by sampling the Beta distribution with parameters \(\alpha, \beta\) chosen above. How many witnesses will we need to have in this setting?