title: “Empirical Probabilities, Discrete Distributions, Visualizations” author: “Evan Robinson 118933641” output: html_document: number_sections: yes toc: yes pdf_document: toc: yes
#Functions for Coin Toss and Die Rolls coin_toss <- function(p=0.5){ coin <- c(0,1) toss <- sample(coin, size = 1, replace = TRUE, prob = c(1-p, p)) return(toss) } coin_toss()
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) }
die_roll() ## Experiments involving Coin tosses and die rolls # 2.1 Die Rolls two_rolls <- function(){ rolls <- replicate(2, die_roll()) return(rolls) } sum_two_rolls <- function(){ rolls <- two_rolls() return(sum(rolls)) } sum_die <- replicate(10000,sum_two_rolls())
round(table(sum_die)/10000, 2)
unfair_die_roll <- function(){ die <- 1:6 probs <- c(1/10,1/10,1/10,1/10,1/10,1/2) roll <- sample(die, size = 5,replace = TRUE, prob = probs) return(roll) } E_happens <- function(){ rolls <- unfair_die_roll() sum <- sum(rolls) if (sum == 13 ){ return(1) } else { return(0) } } G_happens <- function(){ rolls <- unfair_die_roll() four_max <- all(rolls <= 4) if (four_max){ return(1) } else{ return(0) } } G_happens() #2.3 Coin Tosses heads_in_two_tosses <- function(){ head_count <- 0 toss_count <- 0 while(toss_count < 2){ toss <- coin_toss(p=0.8) if(toss == 0){ toss_count = toss_count + 1 } else{ head_count = head_count + 1 toss_count = toss_count + 1 } } return(head_count) }
heads_in_three_tosses <- function(){ head_count <- 0 toss_count <- 0 while(toss_count < 3){ toss <- coin_toss(p=0.8) if(toss == 0){ toss_count = toss_count + 1 } else{ head_count = head_count + 1 toss_count = toss_count + 1 } } return(head_count) }
heads_in_four_tosses <- function(){ head_count <- 0 toss_count <- 0 while(toss_count < 4){ toss <- coin_toss(p=0.8) if(toss == 0){ toss_count = toss_count + 1 } else{ head_count = head_count + 1 toss_count = toss_count + 1 } } return(head_count) } #2.4 Coin Tosses with Infinite Sample Space tails_until_first_head <- function(){ tail_count <- 0 head_count <- 0 while(head_count < 1){ toss <- coin_toss(p=0.8) if (toss ==0){ tail_count = tail_count + 1 } else{ head_count = head_count + 1 } } return(tail_count) }
tails_until_second_head <- function(){ tail_count <- 0 head_count <- 0 while(head_count < 2){ toss <- coin_toss(p=0.8) if (toss ==0){ tail_count = tail_count + 1 } else{ head_count = head_count + 1 } } return(tail_count) }
tails_until_third_head <- function(){ tail_count <- 0 head_count <- 0 while(head_count < 3){ toss <- coin_toss(p=0.8) if (toss ==0){ tail_count = tail_count + 1 } else{ head_count = head_count + 1 } } return(tail_count) }
plot_binom <- function(n,p){ n <- 40 p_values <-
c(0.05,0.1,0.4,0.6,0.9,0.95) par(mfrow = c(3,2)) for (p in
p_values){
x <- 0:n binom <- dbinom(x, size = n, prob = p) plot(x,binom,type
=‘h’) } } # 3.2 Poisson Distribution n <- 100 plot_pois <-
function(mu){ x <- 0:n pois <- dpois(x, lambda = mu) mu_X <- mu
plot(x, pois, type = “h”, lwd = 5, main = paste(“Poisson density: mu =”,
mu), xlab = “x”, ylab = “P(X=x)”) abline(v = mu_X, col = ‘red’, lwd = 4)
} plot_pois(0.5) plot_pois(1) plot_pois(5) plot_pois(8) plot_pois(20)
plot_pois(50)
cat(“As mu increases, the Poisson density plots go from right skewed to more and more left skewed”)
#4.1.2 mu <- 5 sigmas <- c(0.2,0.4, 0.8,1,1.3,1.8,2)
x <- seq(mu-4sigmas[1], mu+4sigmas[1], length.out = 1000) nor <- dnorm(x, mean = mu, sd = sigmas[1]) plot(x,nor,type =‘l’) for (i in 2:length(sigmas)){ nor <- dnorm(x, mean = mu, sd = sigmas[i]) lines(x, nor, col = i) } cat(“as the sigma values increase, the distribution gets wider”) #4.1.2.1 mu <- 5 sigmas <- c(0.2,0.4, 0.8,1,1.3,1.8,2) a <- mu - 4max(sigmas) b <- mu+4max(sigmas) x <- seq(a, b, length.out=1000) nor <- pnorm(x, mean =mu, sd=sigmas[1]) plot(x,nor, type =‘l’, xlim = c(a,b), ylim = c(0,1)) for (i in 2:length(sigmas)){ nor <- pnorm(x, mean =mu, sd=sigmas[i]) lines(x, nor, col = i) } }
#4.1.3 mus <- c(-1,-0.5,0,0.4,0.9,2.5,4) sigma <- 0.4 a <- max(mus) - 4sigma b <- max(mus)+4sigma x <- seq(mu[1]-4sigma, mu[1] +4sigma, length.out = 1000) nor <- dnorm(x, mean = mus[1], sd = sigma) plot(x,nor,type =‘l’,xlim =c(a,b)) for (i in 2:length(mus)){ nor <- dnorm(x, mean = mus[i], sd = sigma) lines(x, nor, col = i) } cat(“as mu values increase, plot shifts to the right”) #4.1.3.2 mus <- c(-1,-0.5,0,0.4,0.9,2.5,4) sigma <- 0.4 a <- max(mus) - 4sigma b <- max(mus)+4sigma x <- seq(mu[1]-4sigma, mu[1] +4sigma, length.out = 1000) nor <- pnorm(x, mean = mus[1], sd = sigma) plot(x,nor,type =‘l’,xlim =c(a,b), ylim =c(0,1)) for (i in 2:length(mus)){ nor <- pnorm(x, mean = mus[i], sd = sigma) lines(x, nor, col = i) }
alpha <- 1 beta_values <- c(0.2, 0.6, 1, 1.5, 2) x <- seq(0,10, length.out = 1000) gamm <- dgamma(x, shape = alpha, scale = beta_values[1]) plot(x, gamm, type = ‘l’, xlim = c(mu-4max(sigmas),mu +4max(sigmas))) for (i in 2:length(beta_values)){ gamm <- dgamma(x, shape = alpha, scale = beta_values[i]) lines(x, gamm, col = i) }
#4.2.2 alpha <- 0.6 beta_values <- c(0.2, 0.6, 1, 1.5, 2) x <- seq(0,10, length.out = 1000) gamm <- dgamma(x, shape = alpha, scale = beta_values[1]) plot(x, gamm, type = ‘l’) for (i in 2:length(beta_values)){ gamm <- dgamma(x, shape = alpha, scale = beta_values[i]) lines(x, gamm, col = i) } #4.2.3 alpha <- 2 beta_values <- c(0.2, 0.6, 1, 1.5, 2) x <- seq(0,10, length.out = 1000) gamm <- dgamma(x, shape = alpha, scale = beta_values[1]) plot(x, gamm, type = ‘l’) for (i in 2:length(beta_values)){ gamm <- dgamma(x, shape = alpha, scale = beta_values[i]) lines(x, gamm, col = i) } #4.2.4 alpha <- 5 beta_values <- c(0.2, 0.6, 1, 1.5, 2) x <- seq(0,10, length.out = 1000) gamm <- dgamma(x, shape = alpha, scale = beta_values[1]) plot(x, gamm, type = ‘l’) for (i in 2:length(beta_values)){ gamm <- dgamma(x, shape = alpha, scale = beta_values[i]) lines(x, gamm, col = i) } #4.3 Beta Distribution alpha_values <- c(1, 5, 0.3, 2, 1) beta_values <- c(3, 1, 0.6, 2, 1) x <- seq(0,1, length.out = 1000) bet <- dbeta(x, shape1 = alpha_values[1], shape2 = beta_values[1]) plot(x, bet, type = ‘l’) for (i in 2:length(alpha_values)){ bet <- dbeta(x, shape1 = alpha_values[i], shape2 = beta_values[i]) lines(x, bet, col = i) } legend(‘top’, legend = paste(‘alpha =’, alpha_values, ‘beta =’, beta_values), col = 1:length(alpha_values)) ## 5 Simulating Monty-Hall Problem #5.1 set.seed(123) monty_3doors_noswitch <- function(){ prize_door <- sample(1:3, 1) first_choice <- sample(1:3, 1) if (first_choice == prize_door){ return(1) } else { return(0) } } monty_3doors_switch <- function(){ prize_door <- sample(1:3, 1) first_choice <- sample(1:3, 1) remaining_doors <- setdiff(1:3, c(first_choice, prize_door)) revealed <- sample(remaining_doors, 1) door_after_switch <- setdiff(1:3, c(first_choice, revealed)) if (door_after_switch[[1]] == prize_door){ return(1) } else { return(0) } } noswitch <- mean(replicate(10000, monty_3doors_noswitch())) switch <- mean(replicate(10000, monty_3doors_switch())) cat(“Empirical win probability with switch:”,switch ) cat(“Empirical win probability without switching:”,noswitch)
#5.2 set.seed(123) monty_10doors_noswitch <- function(){ prize_door <- sample(1:10, 1) first_choice <- sample(1:10, 1) if (first_choice == prize_door){ return(1) } else { return(0) } } monty_10doors_switch <- function(){ prize_door <- sample(1:10, 1) first_choice <- sample(1:10, 1) remaining_doors <- setdiff(1:10, c(first_choice, prize_door)) revealed <- sample(remaining_doors, 8) door_after_switch <- setdiff(1:10, c(first_choice, revealed)) if (door_after_switch[[1]] == prize_door){ return(1) } else { return(0) } } noswitch_ten <- mean(replicate(10000, monty_10doors_noswitch())) switch_ten <- mean(replicate(10000, monty_10doors_switch())) cat(“Empirical win probability with switch:”,switch_ten ) cat(“Empirical win probability without switching:”,noswitch_ten)
set.seed(123)
monty_10doors_mod_noswitch <- function(){ prize_door <- sample(1:10, 1) first_choice <- sample(1:10, 1) if (first_choice == prize_door){ return(1) } else { return(0) } } monty_10doors_mod_switch <- function(){ prize_door <- sample(1:10, 1) first_choice <- sample(1:10, 1) remaining_doors <- setdiff(1:10, c(first_choice, prize_door)) revealed <- sample(remaining_doors, 7) switch_choices <- setdiff(1:10, c(first_choice, revealed)) door_after_switch <- sample(switch_choices,1) if (door_after_switch[[1]] == prize_door){ return(1) } else { return(0) } }
noswitch_ten_mod <- mean(replicate(10000, monty_10doors_mod_noswitch())) switch_ten_mod <- mean(replicate(10000, monty_10doors_mod_switch())) cat(“Empirical win probability with switch:”,switch_ten_mod ) cat(“Empirical win probability without switching:”,noswitch_ten_mod)