(100/1000)*(.4 -.5)

1.96*sqrt(.25/2177)

coin <- c(0,1)

sum(sample(coin, 10, replace=T))/10
sum(sample(coin, 10, replace=T))/10
sum(sample(coin, 10, replace=T))/10
sum(sample(coin, 10, replace=T))/10
sum(sample(coin, 10, replace=T))/10


#The sample command lets us assign probabilities to each instance in coin
sample(coin, 10, prob = c(.52,.48), replace=T)
sample(coin, 10, prob = c(.52,.48), replace=T)
sample(coin, 10, prob = c(.52,.48), replace=T)
sample(coin, 10, prob = c(.52,.48), replace=T)

#The sample command lets us assign probabilities to each instance in coin
sample(coin, 10, prob = c(.43,.57), replace=T)
sample(coin, 10, prob = c(.43,.57), replace=T)
sample(coin, 10, prob = c(.43,.57), replace=T)
sample(coin, 10, prob = c(.43,.57), replace=T)

set.seed(19104)
samp <- rbinom(10000,1,.525)
head(samp, 50)

#The cumsum() function adds up every value that occurs previous in a vector

cumulative <- cumsum(samp)

#Divide by the position to get the average at each value

relative.error <- cumulative/1:length(samp)

head(samp)
head(cumulative)
head(relative.error)


plot(1:10000,relative.error, xlab="Sample Size",ylab="Fetterman Percent", type="l",
     ylim=c(0.3,0.7), xlim=c(0,1000))
abline(h=.525, lty=2, col="firebrick")

plot(1:10000,relative.error, xlab="Sample Size",ylab="Fetterman Percent", type="l",
     ylim=c(0.3,0.7))
abline(h=.525, lty=2, col="firebrick")


samp[1:10]



set.seed(19104)
binomial.sample <- matrix(NA, nrow=10000, ncol=1000)

for(i in 1:1000){
  binomial.sample[,i]<- cumsum(rbinom(10000,1,.525))/1:10000
}

plot(1:10000,binomial.sample[,1], xlab="Sample Size",ylab="Fetterman Percent", type="n",
     ylim=c(0.3,0.7))
for(i in 1:1000){
  points(1:10000,binomial.sample[,i], type="l", col="gray60")
}
abline(h=.525, lty=2, col="firebrick")
points(1:10000, binomial.sample[,1], col="magenta", type="l")
abline(v=c(100,500,1000,5000,10000), lty=2, col="darkblue")

means <- rep(NA,1000)

for(i in 1:1000){
means[i] <- mean(binomial.sample[,i])
}

plot(1:1000, means, pch=16, 
     ylim=c(0.3,0.7), xlab="Average of Averages", ylab="Sample Size")
abline(h=.525, col="firebrick", lty=2)


par(mfrow=c(2,3))
plot(density(binomial.sample[100,]), main="Sampling Distribution at n=100", xlim=c(0.3,0.7))
abline(v=.525, col="firebrick", lty=2)
plot(density(binomial.sample[500,]), main="Sampling Distribution at n=500", xlim=c(0.3,0.7))
abline(v=.525, col="firebrick", lty=2)
plot(density(binomial.sample[1000,]), main="Sampling Distribution at n=1000", xlim=c(0.3,0.7))
abline(v=.525, col="firebrick", lty=2)
plot(density(binomial.sample[5000,]), main="Sampling Distribution at n=5000", xlim=c(0.3,0.7))
abline(v=.525, col="firebrick", lty=2)
plot(density(binomial.sample[10000,]), main="Sampling Distribution at n=10000", xlim=c(0.3,0.7))
abline(v=.525, col="firebrick", lty=2)


var(binomial.sample[100,])
var(binomial.sample[500,])
var(binomial.sample[1000,])
var(binomial.sample[5000,])
var(binomial.sample[10000,])

#(Note that the dnorm function wants sd, not var, so we do sqrt(9)=3 in that spot)
plot(seq(55,73,.01), dnorm(seq(55,73,.01), mean = 64, sd = 3), type="l", 
     main="Female Height in the US", xlab="Height in Inches", ylab = "Probability of Occuring")

x <- seq(1,20,1)
mean(x*3)==3*mean(x)


x <- 1:20
z <- 21:40
mean(x+z)==mean(x) + mean(z)

x <- seq(1,20,1)
var(x*3)==3*var(x)
var(x*3)==(3^2)*var(x)

x <- NA
for(i in 1:1000){
x[i]  <- rnorm(1, mean=0, sd=1)
}
var(x)

par(mfrow=c(2,3))
plot(density(binomial.sample[100,]), main="Sampling Distribution at n=100", xlim=c(0.3,0.7))
abline(v=.525, col="firebrick", lty=2)
plot(density(binomial.sample[500,]), main="Sampling Distribution at n=500", xlim=c(0.3,0.7))
abline(v=.525, col="firebrick", lty=2)
plot(density(binomial.sample[1000,]), main="Sampling Distribution at n=1000", xlim=c(0.3,0.7))
abline(v=.525, col="firebrick", lty=2)
plot(density(binomial.sample[5000,]), main="Sampling Distribution at n=5000", xlim=c(0.3,0.7))
abline(v=.525, col="firebrick", lty=2)
plot(density(binomial.sample[10000,]), main="Sampling Distribution at n=10000", xlim=c(0.3,0.7))
abline(v=.525, col="firebrick", lty=2)


mean(binomial.sample[100,])
mean(binomial.sample[500,])
mean(binomial.sample[1000,])
mean(binomial.sample[5000,])
mean(binomial.sample[10000,])


var(binomial.sample[100,])
var(binomial.sample[500,])
var(binomial.sample[1000,])
var(binomial.sample[5000,])
var(binomial.sample[10000,])

#Calculated variances
.249/100
.249/500
.249/1000
.249/5000
.249/10000


par(mfrow=c(2,3))
x <- seq(0,1,.0001)

plot(density(binomial.sample[100,]), main="Sampling Distribution at n=100", xlim=c(0.3,0.7))
points(x, dnorm(x, mean=.525, sd=.499/sqrt(100)), col="darkblue", type="l")
abline(v=.525, col="firebrick", lty=2)
plot(density(binomial.sample[500,]), main="Sampling Distribution at n=500", xlim=c(0.3,0.7))
points(x, dnorm(x, mean=.525, sd=.499/sqrt(500)), col="darkblue", type="l")
abline(v=.525, col="firebrick", lty=2)
plot(density(binomial.sample[1000,]), main="Sampling Distribution at n=1000", xlim=c(0.3,0.7))
points(x, dnorm(x, mean=.525, sd=.499/sqrt(1000)), col="darkblue", type="l")
abline(v=.525, col="firebrick", lty=2)
plot(density(binomial.sample[5000,]), main="Sampling Distribution at n=5000", xlim=c(0.3,0.7))
points(x, dnorm(x, mean=.525, sd=.499/sqrt(5000)), col="darkblue", type="l")
abline(v=.525, col="firebrick", lty=2)
plot(density(binomial.sample[10000,]), main="Sampling Distribution at n=10000", xlim=c(0.3,0.7))
points(x, dnorm(x, mean=.525, sd=.499/sqrt(10000)), col="darkblue", type="l")
abline(v=.525, col="firebrick", lty=2)


ev <- .43
var <- (.43*.57)/1000

plot(seq(.3,.56,.001), dnorm(seq(.3,.56,.001), mean = ev, sd = sqrt(var)), type="l",
     main="Sampling Distribution for Trump Approval", 
     xlab="Proportion Approving Trump", ylab = "Probability of Getting a Sample with this Answer")
abline(v=.43, lty=2, col="firebrick")

ev <- .43
var <- (.43*.57)/500

plot(seq(.3,.56,.001), dnorm(seq(.3,.56,.001), mean = ev, sd = sqrt(var)), type="l",
     main="Sampling Distribution for Trump Approval", 
     xlab="Proportion Approving Trump", ylab = "Probability of Getting a Sample with this Answer")
abline(v=.43, lty=2, col="firebrick")

source("https://raw.githubusercontent.com/marctrussler/IIS-Data/main/NormalShader.R")

low <- qnorm(.025, mean = ev, sd = sqrt(var))
high <- qnorm(.025, mean = ev, sd = sqrt(var), lower.tail=F)

plot(seq(.3,.56,.001), dnorm(seq(.3,.56,.001), mean = ev, sd = sqrt(var)), type="l",
     main="Sampling Distribution for Trump Approval", 
     xlab="Proportion Approving Trump", ylab = "Probability of Getting a Sample with this Answer")
normal.shader(min=low, max=high, value.one = low, value.two=high, mu=ev, sd=sqrt(var), between=T, color="orange")
abline(v=.43, lty=2, col="firebrick")



qnorm(.025, mean = ev, sd = sqrt(var))
qnorm(.025, mean = ev, sd = sqrt(var), lower.tail=F)

qnorm(.025, mean = ev, sd = sqrt(var)) - ev
qnorm(.025, mean = ev, sd = sqrt(var), lower.tail=F) - ev

.0434/sqrt(var)

var <- (.43*.57)/1000

ev + 1.96*sqrt(var)
ev - 1.96*sqrt(var)


qnorm(.025, mean = ev, sd = sqrt(var))
#Yes

pi <- seq(0,1, .01)
one.minus.pi <- 1-pi
se <- sqrt((pi*one.minus.pi)/500)

plot(pi, se, type="l")
abline(v=.5, lty=2, col="firebrick")

ev <- .45
var <- (.5*.5)/1000

plot(seq(.3,.56,.001), dnorm(seq(.3,.56,.001), mean = ev, sd = sqrt(var)), type="l",
     main="Sampling Distribution for Trump Approval", 
     xlab="Proportion Approving Trump", ylab = "Probability of Getting a Sample with this Answer")
abline(v=.45, lty=2, col="firebrick")



1.96*sqrt(var)

moe <- sqrt((.5*.5)/1000)

result <- NA

for(i in 1:10000){
sample <- rbinom(1000, 1, .43)
prop <- mean(sample)
low <- prop - 1.96*moe
high <- prop + 1.96*moe

#Test to see whether truth is in interval
result[i] <- low<=.43 & high>=.43
}

mean(result)