coin <- c(0,1)

sample(coin,7,replace=T)

sample(coin,7,replace=T)
sample(coin,7,replace=T)
sample(coin,7,replace=T)

sum(sample(coin,7,replace=T))

for(i in 1:100000){
  sum(sample(coin,7,replace=T))
}

set.seed(19104)
num.heads <- rep(NA, 100000)
for(i in 1:100000){
 num.heads[i] <-  sum(sample(coin,7,replace=T))
}
head(num.heads)

table(num.heads)
prop.table(table(num.heads))

set.seed(19104)
dice <- c(1,2,3,4,5,6)
sample(dice,1,replace=T)

set.seed(19104)
dice <- c(1,2,3,4,5,6)
result <- rep(NA,10000)
  for(i in 1:10000){
    result[i] <- sample(dice,1,replace=T)
  }
table(result)

#Estimating the probability of guessing a locker combination

#Set Seed 
set.seed(19104)

#Draw prime combination
prime <- sample(seq(1,50),3, replace=F)

#How do we determine if two vectors are exactly equal?
c(2,34,29)==prime

#Use the all command to determine if all numbers match
all(c(28,21,37)==prime)

#Use a loop to repeatedly sample and check

results <- rep(NA, 1E6)

for(i in 1:length(results)){
  results[i] <- all(sample(seq(1,50),3, replace=F)==prime)
}

prop.table(table(results))

#Approximately equivalent to what we calculated via math. 

choose(70,5)

num <- factorial(365)/factorial(345)
num

(choose(365,20)*factorial(20))/365^20

set.seed(19104)
birthdays <- seq(1,365,1)
s <- sample(birthdays, 20, replace=T)
s
anyDuplicated(s)>0

set.seed(19104)
birthdays <- seq(1,365,1)
result <- rep(NA, 10000)
for(i in 1:10000){
  s <- sample(birthdays, 20, replace=T)
  result[i] <- anyDuplicated(s)>0
}
table(result)

set.seed(19104)
birthdays <- seq(1,365,1)
class.size <- seq(5,100,1)
prob.sim <- rep(NA, length(class.size))
prob.math <- rep(NA, length(class.size))

for(j in 1:length(class.size)){
result <- rep(NA, 10000)
for(i in 1:10000){
  s <- sample(birthdays, class.size[j], replace=T)
  result[i] <- anyDuplicated(s)>0
}
prob.sim[j] <- mean(result)
prob.math[j] <- 1-(choose(365,class.size[j])*factorial(class.size[j]))/365^class.size[j]
}

plot(class.size, prob.math, col="dodgerblue", type="b", xlab="Class Size",
     ylab="Probability of Shared Birthday", pch=16)
points(class.size, prob.sim, col="firebrick", type="b", pch=16)
legend("topleft", c("Math","Simulation"), pch=c(1,1), col=c("dodgerblue", "firebrick"))


anes <- read.csv("https://raw.githubusercontent.com/marctrussler/IIS-Data/main/ANES2020Clean.csv")

table(anes$race)
prop.table(table(anes$race))

table(anes$gender)

table(anes$race[anes$gender=="Female"])
prop.table(table(anes$race[anes$gender=="Female"]))

joint.prob <- prop.table(table(race=anes$race, gender=anes$gender))
joint.prob

sum(joint.prob)

sum(joint.prob[,1])
mean(anes$gender=="Female",na.rm=T)

#(There is some rounding error...)

sum(joint.prob[6,])
mean(anes$race=="White, non-Hispanic", na.rm=T)

joint.prob[6,1]/sum(joint.prob[,1])
mean(anes$gender=="Female" & anes$race=="White, non-Hispanic",na.rm=T)/mean(anes$gender=="Female",na.rm=T)

#First I will write one game without a loop:

#Pick a door, which unbeknownst to us contains a goat or a car:
pick <- sample(c("goat","goat","car"),1)

#If the player stays they get their initial selection
stay = pick

#If they switch, what happens depends on what their initial pick is

if(pick=="car"){
  #If they initially picked car, switching necessarily gets them a goat.
  switch = "goat"
} else {
  #If they initially picked a goat, switching necessarily gets them the car. 
  #This is the key to the whole thing. Monty never reveals the car, so if you pick
  #a goat he *must* reveal the other goat. 
  switch="car"
}
stay
switch

#Now put this whole thing in a loop, capturing what happens if you stay vs switch for 
#each initial selection:
stay <- NA
switch <- NA
for(i in 1:10000){
  pick <- sample(c("goat","goat","car"),1)
  
  stay[i] = pick
  
  if(pick=="car"){
    switch[i] = "goat"
  } else {
    switch[i]="car"
  }
}
mean(stay=="car")
mean(switch=="car")


#One Monty Hall game
set.seed(19104)
#Randomply place the Cars and goats
placement <- sample(c("Goat", "Goat","Car"),3, replace=F)
#In this case the placement is Goat, Goat, Car

#Person randomly selects a door:
door.selection <- sample(seq(1,3),1)
#In this case the person has selected door 1, that has a goat

#Now Monty must open a door
#He can't open the door we have opened, and can't open the door with the car
#Record which doors have goats:
goats <- which(placement=="Goat")
#Can't open the door you've selected
goats <- goats[goats!=door.selection]
#We are going to use an if statement here to represent the two situations
#In this situation there is only one door Monty can open, but if we had chosen the car first
#there would be two doors that monty can open and he has to choose one
if(length(goats)==1){
  open <- goats
} else {
  open <- sample(goats,size=1)
}
#
#The door that is switched to is the remaining door
switch.selection <- c(1,2,3)[c(-door.selection,-open)]

#Finally we see if the car is won either in the stay or switch condition
#Stay
placement[door.selection]=="Car"
#Switch
placement[switch.selection]=="Car"

#In this case the car is won when switching. 
#It should be clear from this that if you pick a goat initially and switch you *necessarily* win 
#the car. Because you have a 2/3 probability of picking a goat, you win the car 2/3 of the time if you
#switch. Is that what R says?

win.stay <- rep(NA, 1000)
win.switch <- rep(NA, 1000)

for(i in 1:1000){
#Placement
placement <- sample(c("Goat", "Goat","Car"),3, replace=F)
#Person randomly selects a door:
door.selection <- sample(seq(1,3),1)

#Now Monty must open a door
goats <- which(placement=="Goat")
goats <- goats[goats!=door.selection]
if(length(goats)==1){
  open <- goats
} else {
  open <- sample(goats,size=1)
}

switch.selection <- c(1,2,3)[c(-door.selection,-open)]



win.stay[i] <- placement[door.selection]=="Car"
win.switch[i] <- placement[switch.selection]=="Car"
}

table(win.stay)
table(win.switch)



set.seed(19104)
coin <- c(0,1)

all.heads <- rep(NA, 1000000)

for(i in 1:length(all.heads)){
  s <- sample(coin, 15, replace=T)
  all.heads[i] <- sum(s)==15
}

sum(all.heads)


set.seed(19104)
coin <- c(0,1)
s <- sample(coin, 40, replace=T)
s

rle(s)$lengths
rle(s)$values

rle(s)$lengths[rle(s)$values==1]

sort(rle(s)$lengths[rle(s)$values==1], decreasing = T)

sort(rle(s)$lengths[rle(s)$values==1], decreasing = T)[1]>=15

str.res <- NA
for(i in 1:6000){
  s <- sample(coin, 40, replace=T)
  str.res[i] <- sort(rle(s)$lengths[rle(s)$values==1], decreasing = T)[1]>=15
}
table(str.res)
any(str.res)

adv.streak <- function(n.adv, n.yrs, streak){
  str.res <- NA
  for(i in 1:n.adv){
    s <- sample(coin, n.yrs, replace=T)
    str.res[i] <- sort(rle(s)$lengths[rle(s)$values==1], decreasing = T)[1]>=streak
  }
  return(any(str.res))
}

adv.streak(n.adv=6000, n.yrs = 40, streak=15)

any.lucky <- NA
for(i in 1:100){
any.lucky[i] <- adv.streak(n.adv=6000, n.yrs = 40, streak=15)
}
table(any.lucky)

any.lucky <- NA
for(i in 1:1000){
any.lucky[i] <- adv.streak(n.adv=5, n.yrs = 10, streak=10)
}
table(any.lucky)

set.seed(19104)
dice <- c(1,2,3,4,5,6)
result <- rep(NA, 10000)
for(i in 1:10000){
  s <- sample(dice,4, replace=T)
  result[i] <- var(s)==0
}

mean(result)

set.seed(19104)
dice <- c(1,2,3,4,5,6)

result <- rep(NA, 10000)

for(i in 1:10000){
  roll <- sample(dice,2, replace=T)
  result[i] <- sum(roll)==10
}

mean(result)

set.seed(19104)
dice <- c(1,2,3,4,5,6)

result <- rep(NA, 10000)

for(i in 1:10000){
 roll <- sample(dice,2, replace=T)
result[i] <- sum(roll)
}

table(result)

barplot(prop.table(table(result)),
        xlab="Sum of Two Dice",
        ylab="P(Sum)")




set.seed(19104)
bag1 <- c(rep("R",3), rep("B",5))
bag2 <- c(rep("R",4), rep("B",4))

bags <- cbind(bag1,bag2)

result <- rep(NA,10000)

for(i in 1:length(result)){
  bag <- sample(1:2,1)  
  ball <- sample(bags[,bag], 1)  
  result[i] <- ball=="R"
}

mean(result)


philly.truth <- c(rep(1,6), rep(0,4))
result <- rep(NA,10000)
for(i in 1:length(result)){
  result[i] <- mean(sample(philly.truth, 7, replace=T))>.5
}
mean(result)

philly.truth <- c(rep(1,6), rep(0,4))
series.length <- seq(7,61,2)
percent.series.win <- rep(NA, length(series.length))

for(j in 1:length(percent.series.win)){
 result <- rep(NA,10000)
 for(i in 1:length(result)){
   result[i] <- mean(sample(philly.truth, series.length[j], replace=T))>.5
 }
percent.series.win[j] <-  mean(result)
}

plot(series.length, percent.series.win, type="b", pch=16)