1 Basic R manipulation.

  1. 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
  1. 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 x.
#Here, we do the evens, threes, and x all given by instruction
even <- seq(2, 100, by = 2)
threes <- seq(1, 100, by = 3)
x <- seq(0, 2, length.out = 100)
  1. Define a function that computes the square of the input value, call this function square:
square <- function(num){
  return (num^2)
}
  1. 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.
#We make the original plot with sine and add cosine to the graph
x <- seq(0, 10, length.out = 1000)
plot(x, sin(x), type = "l", col = "red", xlim = c(0, 10))
lines(x, cos(x), col = "red")
abline( h = 0, v = 0)

xs <- seq(0, 2, length.out = 1000)
plot(xs, square(xs))

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 ‘cion_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_toss <- function(p) {
  coin <- c(0,1)
  toss <- sample(coin, size = 1, replace = TRUE, prob = c(1-p, p))
  return(toss)
}

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.

#This is a simple die roll, and the probabilities are all 1/6
die_roll <- function() {
  die <- 1:6
  probs <- c(1/6,1/6,1/6,1/6,1/6,1/6)
  roll <- sample(die, size = 1, replace = TRUE, prob = probs)
  return(roll)
}

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.

#This shows the experiment where we roll twice, using the die_roll function
two_rolls <- function() {
  roll1 <- replicate(1, die_roll())
  roll2 <- replicate(1, die_roll())
  return(roll1 + roll2)
}
replicate(1, two_rolls())
## [1] 8

Use the replicate function to repeat this experiment 10000 times, store the output of the experiment in a variable “sum_die”.

sum_die <- replicate(10000, two_rolls())

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.

#Used a table here to show the empirical probability where I divided the sum_dies by the length of sum_die
emp_prob <- round(table(sum_die) / length(sum_die), 2)
emp_prob
## sum_die
##    2    3    4    5    6    7    8    9   10   11   12 
## 0.03 0.05 0.08 0.11 0.14 0.17 0.14 0.11 0.09 0.05 0.03

Compare your answers with the theoretical values.

Plot the following: relative frequency of getting sum equal to 8 as a function of number of trials, for this experiment.

#To plot the relative frequency we have to replicate the sum_die equal to 8, and get that result done 1000 times.
B <- 1000
reps <- replicate(B, sum_die == 8)
#We must go through from 1 to 1000 and get the mean of the reps, so we can then plot them
x <- 1:B
y <- c()
for (i in 1:B) {
  y <- c(y, mean(reps[1:i]))
}

plot(x,y, type = 'l', main = "Relative frequencies for the Experiment", ylim = c(0,1), xlab = "Number of Trials", ylab = "Relative Frequencies")
abline(h = 0.15, col = 'red')
abline(h = 0, v = 0)

What do you notice about this graph as the number of trials increase? YOUR ANSWER: It seems to stabilize more compared to when there are fewer trials

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

#this time, our probabilites are different and unfair. so we have to get a sample of the new probabilites and see if the sum is 13.
E_happens <- function() {
rolls <- sample(1:6, 5, replace = TRUE, prob = c(0.1,0.1,0.1,0.1,0.1,0.5))
  if (sum(rolls) == 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

#This time, we do the same thing and see if the max of the 5 rolls is less than or equal to 4.
G_happens <- function() {
  rolls <- sample(1:6, 5, replace = TRUE, prob = c(0.1,0.1,0.1,0.1,0.1,0.5))
if (max(rolls) <= 4) {
  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)\).

#to get the empirical approx. we must replicate the experiments 10000 times and find the mean to get the average value.
emp_prob_E <- mean(replicate(10000, E_happens()))

emp_prob_G <- mean(replicate(10000, G_happens()))

What are the approximate probabilities? YOUR ANSWER: 0.0111 for E_happens and 0.0068 for G_happens

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.

#we have to sample the experiment if the probability for heads was 0.8, then we return the number of heads that were found in 2 tosses.
simulate_coin_toss <- function() {
  
  coin_tosses <- sample(c("H", "T"), 2, replace = TRUE, prob = c(0.8, 0.2))

  num_heads <- sum(coin_tosses == "H")
  
  return(num_heads)
}

print(simulate_coin_toss())
## [1] 2

Use the above function to calculate the emprical probability distribution table for \(X\). Compare this with the theoretical probability distribution of \(X\).

#replicate the experiment 10000 times, and then find the empirical probability through dividing the num heads by the amount of times we replicated it, 
#then find the theoretical probability by using a vector to show the calculations required. 
num_heads <- replicate(10000, simulate_coin_toss())

emp_x <- table(num_heads) / 10000
theo_prob <- c((1-0.8)^2, 2*0.8*0.2, 0.8^2)

theo_prob
## [1] 0.04 0.32 0.64
emp_x
## num_heads
##      0      1      2 
## 0.0402 0.3217 0.6381

Calculate empirical estimates for the expected value and the variance of \(X\). Compare these with the actual theoretical calculations.

#First, we replicate this experiment and to find the empirical expected value we do the mean of this experiment. For the variance, we simply use the var function.
num_heads <- replicate(10000, simulate_coin_toss())

empirical_exp_val <- mean(num_heads)
emp_variance <- var(num_heads)

#For the theoretical probabilities, we use the formula E(x) = np and V(x) = np(1-p)

n<-2
p<-0.8
theoretical_exp_val <- n * p
theoretical_variance <- n * p * (1 - p)

#Then we print to show they are the same if not very close.

theoretical_exp_val
## [1] 1.6
theoretical_variance
## [1] 0.32
empirical_exp_val
## [1] 1.6038
emp_variance
## [1] 0.3236579

3.3.2 n = 3

Setup a function that simulates the random variable \(X\) when we toss the coin three times.

#we have to sample the experiment if the probability for heads was 0.8, then we return the number of heads that were found in 3 tosses.
simulate_coin_toss <- function() {
  
  coin_tosses <- sample(c("H", "T"), 3, replace = TRUE, prob = c(0.8, 0.2))

  num_heads <- sum(coin_tosses == "H")
  
  return(num_heads)
}

Use the above function to calculate the emprical probability distribution table for \(X\). Compare this with the theoretical probability distribution of \(X\).

#replicate the experiment 10000 times, and then find the empirical probability through dividing the num heads by the amount of times we replicated it, 
#then find the theoretical probability by using a vector to show the calculations required. 
num_heads <- replicate(10000, simulate_coin_toss())
p <- 0.8
n <- 3

theoretical_probs <- c((1 - p)^n, choose(n, 1) * p * (1 - p)^(n - 1), choose(n, 2) * p^2 * (1 - p)^(n - 2), p^n)

theoretical_probs
## [1] 0.008 0.096 0.384 0.512
emp_x <- table(num_heads) / 10000

emp_x
## num_heads
##      0      1      2      3 
## 0.0072 0.1026 0.3765 0.5137

Calculate emprical estimates for the expected value and the variance of \(X\). Compre these with the actual theoretical calculations.

#First, we replicate this experiment and to find the empirical expected value we do the mean of this experiment. For the variance, we simply use the var function.
num_heads <- replicate(10000, simulate_coin_toss())

empirical_exp_val <- mean(num_heads)
emp_variance <- var(num_heads)

#For the theoretical probabilities, we use the formula E(x) = np and V(x) = np(1-p)
n <- 3
p <- 0.8
theoretical_exp_val <- n * p
theoretical_variance <- n * p * (1 - p)
#Then we print to show they are the same if not very close.
empirical_exp_val
## [1] 2.3914
emp_variance
## [1] 0.489255
theoretical_exp_val
## [1] 2.4
theoretical_variance
## [1] 0.48

3.3.3 n = 4

Setup a function that simulates the random variable \(X\) when we toss the coin four times.

#we have to sample the experiment if the probability for heads was 0.8, then we return the number of heads that were found in 4 tosses.
simulate_coin_toss <- function() {
  
  coin_tosses <- sample(c("H", "T"), 4, replace = TRUE, prob = c(0.8, 0.2))

  num_heads <- sum(coin_tosses == "H")
  
  return(num_heads)
}

Use the above function to calculate the empirical probability distribution table for \(X\). Compare this with the theoretical probability distribution of \(X\).

#replicate the experiment 10000 times, and then find the empirical probability through dividing the num heads by the amount of times we replicated it, 
#then find the theoretical probability by using a vector to show the calculations required. 
num_heads <- replicate(10000, simulate_coin_toss())

p <- 0.8
n <- 4

theoretical_probs <- c((1 - p)^n, choose(n, 1) * p * (1 - p)^(n - 1), choose(n, 2) * p^2 * (1 - p)^(n - 2), choose(n, 3) * p^3 * (1 - p)^(n - 3), p^n)

theoretical_probs
## [1] 0.0016 0.0256 0.1536 0.4096 0.4096
emp_x <- table(num_heads) / 10000
emp_x
## num_heads
##      0      1      2      3      4 
## 0.0015 0.0242 0.1514 0.4096 0.4133

Calculate emprical estimates for the expected value and the variance of \(X\). Compre these with the actual theoretical calculations.

#First, we replicate this experiment and to find the empirical expected value we do the mean of this experiment. For the variance, we simply use the var function.
num_heads <- replicate(10000, simulate_coin_toss())

empirical_exp_val <- mean(num_heads)
emp_variance <- var(num_heads)

#For the theoretical probabilities, we use the formula E(x) = np and V(x) = np(1-p)
n <- 4
p <- 0.8
theoretical_exp_val <- n * p
theoretical_variance <- n * p * (1 - p)
#Then we print to show they are the same if not very close.

empirical_exp_val
## [1] 3.2076
emp_variance
## [1] 0.6439666
theoretical_exp_val
## [1] 3.2
theoretical_variance
## [1] 0.64

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.

#This function goes and counts the heads and tails, using the coin_toss function. We only continue this function while we have counted less than one head,
#since we are only going until the first head.
tails_till_first_head <- function(p) {
  head_count <- 0
  tail_count <- 0
  while(head_count < 1) {
    toss <- coin_toss(p)
    if (toss == 0) {
      tail_count = tail_count + 1
    } else {
      head_count = head_count + 1
    }
  }
  return(tail_count)
}
#to get the expected number of tails we replicate this experiment and find the mean.
results <- replicate(10000, tails_till_first_head(0.8))
expected <- mean(results)
#Then we find the theo prob and compare the two values.
theo_val <- (1-0.8)/0.8
theo_val
## [1] 0.25
expected
## [1] 0.2462

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.

#This function goes and counts the heads and tails, using the coin_toss function. We only continue this function while we have counted less than 2 heads,
#since we are only going until the second head.
tails_till_first_head <- function(p) {
  head_count <- 0
  tail_count <- 0
  while(head_count <= 1) {
    toss <- coin_toss(p)
    if (toss == 0) {
      tail_count = tail_count + 1
    } else {
      head_count = head_count + 1
    }
  }
  return(tail_count)
}
#to get the expected number of tails we replicate this experiment and find the mean.
results <- replicate(10000, tails_till_first_head(0.8))
expected <- mean(results)
#Then we find the theo prob and compare the two values.
theo_val <- ((2-0.8) / 0.8) - 1
theo_val
## [1] 0.5
expected
## [1] 0.4988

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.

#This function goes and counts the heads and tails, using the coin_toss function. We only continue this function while we have counted less than 3 heads,
#since we are only going until the third head.
tails_till_first_head <- function(p) {
  head_count <- 0
  tail_count <- 0
  while(head_count <= 2) {
    toss <- coin_toss(p)
    if (toss == 0) {
      tail_count = tail_count + 1
    } else {
      head_count = head_count + 1
    }
  }
  return(tail_count)
}
#to get the expected number of tails we replicate this experiment and find the mean.
results <- replicate(10000, tails_till_first_head(0.8))
expected <- mean(results)
#Then we find the theo prob and compare the two values.
theo_val <- ((3-0.8) / 0.8) - 2
theo_val
## [1] 0.75
expected
## [1] 0.7593

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)

#This simply uses the formula above and calculates bayes.
bayes <- function(q, p){
  return(p*q/(p*q + (1-p)*(1-q)))
}

bayes(0.01, 0.8)
## [1] 0.03883495

4.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(q)

for(i in 1:40) {
  new_q <- bayes(q, p)
  q_vals <- c(q_vals, new_q)
  q <- new_q
}

length(q_vals)
## [1] 41

Construct a plot with x-axis 1:40 and y-axis q_vals.

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
q_data <- data.frame(x = 1:41, y = q_vals)

ggplot(data = q_data, aes(x=x, y=y)) + geom_point()+geom_abline(slope = 0, intercept = 0.9, col = 'red')

Repeat the above experiment with \(p=0.4\), \(q=0.99\).

#We must do the exact same task as above, but with different p and q values.
p <- 0.4
q <- 0.99
q_vals <- c(q)

for(i in 1:40) {
  new_q <- bayes(q, p)
  q_vals <- c(q_vals, new_q)
  q <- new_q
}

library(tidyverse)
q_data <- data.frame(x = 1:41, y = q_vals)

ggplot(data = q_data, aes(x=x, y=y)) + geom_point()+geom_abline(slope = 0, intercept = 0.1, col = 'red')

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.

#With the new p and q values, we increment a counter and count how many witnesses we need for the police to consider arresting the blue cab.
q <- 0.01
p <- 0.6
counter <- 0
while (q <= 0.9) {
  counter <- counter + 1
  q <- bayes(q,p)
}
counter
## [1] 17

4.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.

#Now, with two types of witnesses we have two p and one q. We do the same task except distinguish which task to perform using %, whether it is an
#odd or even iteration.
p1 <- 0.6
p2 <- 0.45
q <- 0.01
q_vals <- c(q)

for (i in 1:90) {
  if (i%%2 == 0) {
    new_q <- bayes(q, p2)
    q <- new_q
    q_vals <- c(q_vals, new_q)
  } else {
    new_q <- bayes(q, p1)
    q <- new_q
    q_vals <- c(q_vals, new_q)
  }
}

Construct a plot with x-axis 1:90 and y-axis q_vals.

q_data <- data.frame(x = 1:91, y = q_vals)

ggplot(data = q_data, aes(x=x, y=y)) + geom_line() + geom_abline(slope = 0, intercept = 0.9, col = 'red')

What do you notice? (Compared to case with single witness) It takes much longer to hit the threshold, which depends on the values of p1, p2, and p (from the case with the single witness).

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?

#Now we must perform the same task as before, counting witnesses, which may take longer since we have one reliable witness and one unreliable witness.
q <- 0.01
p_1 <- 0.6
p_2 <- 0.45
counter <- 0
while(q <= 0.9) {
  ps <- c(p_1, p_2)
  p <- sample(ps, size = 1, replace = TRUE)
  q <- bayes(q,p)
  counter <- counter + 1
}
counter
## [1] 102

This 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.

#We must create a vector of all the witness reliabilites, and loop through and perform bayes. Similar to the previous examples, but we just
#have 4 types of witnesses to go through instead of 2.
p <- c(0.6, 0.8, 0.45, 0.3)
q <- 0.01
q_vals <- c()
for (i in 1:100) {
  sampled_p <- sample(p, 1)
  
  new_q <- bayes(q, sampled_p)
  
  q_vals <- c(q_vals, new_q)
  q <- new_q
}

Construct a plot with x-axis 1:100 and y-axis q_vals.

q_data <- data.frame(x = 1:100, y = q_vals)
library(ggplot2)
ggplot(data = q_data, aes(x=x, y=y)) + geom_line() + geom_abline(slope = 0, intercept = 0.9, col = 'red')

What do you notice? How does this compare with the earlier cases? This time, the graph seems to be more unpredictable and rise faster. This is because there are 4 witnesses with all different reliability. This compares to the earlier cases by showing how it reaches the threshold much quicker, has steeper ups and downs, and flattens out after reaching the threshold quicker.

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?

#As before, we are performing the same task of counting every time our threshold isnt met, except this time we have a vector of 4 types of witnesses instead of #two.
q <- 0.01
p <- c(0.6, 0.8, 0.45, 0.3)
counter <- 0
while(q <= 0.9) {
  sample_p <- sample(p, 1)
  q <- bayes(q, sample_p)
  counter <- counter + 1
}
counter
## [1] 65

4.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?