plot(0:14, dbinom(0:14, 14, .1), pch=16, xlab="Number of Left Handed Presidents", ylab="P(X=x)") plot(seq(-100,100,.001),dnorm(seq(-100,100,.001), 7.5,18.5), type="l", xlab="Percent Return", ylab=("P(X=x)")) plot(seq(-100,100,.001),dnorm(seq(-100,100,.001), 7.5,18.5), type="l", xlab="Percent Return", ylab=("P(X=x)")) abline(v=-20, col="firebrick", lty=2) library(rio) sm.pa <- import("https://raw.githubusercontent.com/marctrussler/IIS-Data/main/PAFinalWeeks.csv") head(sm.pa) nrow(sm.pa[sm.pa$senate.topline %in% c("Democrat","Republican"),]) nrow(sm.pa[sm.pa$senate.topline %in% c("Democrat"),]) 2432/4369 set.seed(19140) heads <- c(0,1,2,3) true.prob <- c(.125,.375,.375,.125) sim.heads<- sample(heads,size=15, prob=true.prob,replace=T) sim.prob <- as.numeric(prop.table(table(sim.heads))) plot(heads, true.prob, pch=16, ylim=c(0,.5), xlim=c(-1,4), axes=F, xlab="",ylab="Prob") points(heads, sim.prob, pch=16, col="firebrick") axis(side=1, at=c(0,1,2,3)) axis(side=2, las=2) set.seed(19140) heads <- c(0,1,2,3) true.prob <- c(.125,.375,.375,.125) sim.heads<- sample(heads,size=1000, prob=true.prob,replace=T) sim.prob <- as.numeric(prop.table(table(sim.heads))) plot(heads, true.prob, pch=16, ylim=c(0,.5), xlim=c(-1,4), axes=F, xlab="",ylab="Prob") points(heads, sim.prob, pch=16, col="firebrick") axis(side=1, at=c(0,1,2,3)) axis(side=2, las=2) 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="Sample Index", ylab="Average of Averages") 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) y <-c(.475, .525) plot(c(0,1), y, ylim=c(0,1), pch=16) segments(c(0,1), c(0,0), c(0,1), y) var(binomial.sample[100,]) var(binomial.sample[500,]) var(binomial.sample[1000,]) var(binomial.sample[5000,]) var(binomial.sample[10000,]) 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) (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) se.table <- data.frame( n = c(100, 500, 1000, 5000, 10000), simulation = c(sd(binomial.sample[100,]), sd(binomial.sample[500,]), sd(binomial.sample[1000,]), sd(binomial.sample[5000,]), sd(binomial.sample[10000,])), math = NA ) se.table se.table$math <- sqrt(.249 / se.table$n) se.table 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) qnorm(.005) interval.high <- .525 + 2.57*(.499/sqrt(1:10000)) interval.low <- .525 - 2.57*(.499/sqrt(1:10000)) plot(1:10000,binomial.sample[,1], xlab="Sample Size",ylab="Fetterman Percent", type="n", ylim=c(0.3,0.7)) points(1:10000, interval.high, col="darkorange", type="l", lwd=2) points(1:10000, interval.low, col="darkorange", type="l", lwd=2) 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") points(1:10000, interval.high, col="darkorange", type="l", lwd=2) points(1:10000, interval.low, col="darkorange", type="l", lwd=2) evaluation.values <- seq(.4,.6,.0001) pdf <- dnorm(evaluation.values, mean=.5,sd=.007) plot(evaluation.values,pdf, main="Sampling Distribution for Tied Race", xlim=c(.4,.6), type="l") abline(v=.556, col="firebrick", lwd=2, lty=2) set.seed(19104) samp <- rbinom(1,5000,.5) samp.mean <- samp/5000 samp.mean samp <- rbinom(5000,1,.5) head(samp,100) set.seed(19104) simulation.means <- rbinom(10000,5000,.5)/5000 plot(evaluation.values,pdf, main="Sampling Distribution for Tied Race", xlim=c(.4,.6), type="l") abline(v=.556, col="firebrick", lwd=2, lty=2) points(density(simulation.means), type="l", col="magenta") #install.packages("fGarch") #You can install that package but you don't have to know anything about fgarch or #The dsnorm command. I just wanted to generate a skewed distribution to show you. library(fGarch) inc <- seq(0,100000) plot(inc, dsnorm(inc, mean = 30000, sd = 20000, xi = 2.2), type="l", main="Population", xlab="Income", ylab="Density") legend("topright", c("Mean=30k", "SD=20K")) set.seed(19104) x <- rsnorm(100, mean = 30000, sd = 20000, xi = 2.2) head(x) plot(density(x), main="Sample vs Population", xlab="Income", ylim=c(0,4e-5), xlim=c(0,100000)) points(inc, dsnorm(inc, mean=30000, sd=20000, xi=2.2), col="firebrick", type="l") par(mfrow=c(3,3)) for(i in 1:9){ samp <- rsnorm(100, mean=30000, sd=20000, xi=2.2) plot(density(samp), main="n=100", xlab="Income", ylim=c(0,4e-5), xlim=c(0,100000)) points(inc, dsnorm(inc, mean=30000, sd=20000, xi=2.2), col="firebrick", type="l") } par(mfrow=c(3,3)) for(i in 1:9){ samp <- rsnorm(1000, mean=30000, sd=20000, xi=2.2) plot(density(samp), main="n=1000", xlab="Income", ylim=c(0,4e-5), xlim=c(0,100000)) points(inc, dsnorm(inc, mean=30000, sd=20000, xi=2.2), col="firebrick", type="l") } evaluation.values <- seq(20000,40000,1) pdf <- dnorm(evaluation.values,30000,20000/sqrt(100)) plot(evaluation.values, pdf, type="l", main="Sampling Distribution of Sample Mean, n=100") income.sample.means <- rep(NA, 10000) for(i in 1:10000){ income.sample.means[i]<- mean(rsnorm(100,mean = 30000, sd = 20000, xi = 2.2)) } plot(evaluation.values, pdf, type="l", main="Sampling Distribution of Sample Mean, n=100") points(density(income.sample.means), col="darkorange", type="l") interval.high <- 30000 + 2.57*(20000/sqrt(1:10000)) interval.low <- 30000 - 2.57*(20000/sqrt(1:10000)) plot(c(0,10000),c(0,60000), type="n", main="Sampling Distribution of Mean of Income", xlab="Sample Size", ylab="Sample Mean") points(1:10000, interval.high, col="darkorange", lwd=2, type="l") points(1:10000, interval.low, col="darkorange", lwd=2, type="l") abline(h=30000, lty=2, col="firebrick") income.sample <- import("https://github.com/marctrussler/IIS-Data/raw/main/SkewedSamples.rds") interval.high <- 30000 + 2.57*(20000/sqrt(1:10000)) interval.low <- 30000 - 2.57*(20000/sqrt(1:10000)) plot(c(0,10000),c(0,60000), type="n", main="Sampling Distribution of Mean of Income", xlab="Sample Size", ylab="Sample Mean") for(i in 1:ncol(income.sample)){ points(1:10000, income.sample[,i], type="l", col="gray80") } points(1:10000, interval.high, col="darkorange", lwd=2, type="l") points(1:10000, interval.low, col="darkorange", lwd=2, type="l") points(1:10000, income.sample[,1], col="magenta", type="l") abline(h=30000, lty=2, col="firebrick") #Randomly sample and take a mean means <- rep(NA, 10000) for(i in 1:length(means)){ means[i] <- mean(rnorm(100, mean=35, sd=20)) } #Is estimator unbiased? plot(density(means-35), main="Bias of the Sample Mean") abline(v=mean(means-35), lty=2, col="firebrick") #Yes first.obs <- rep(NA, 10000) for(i in 1:length(first.obs)){ first.obs[i] <- mean(rnorm(1, mean=35, sd=20)) } #Is estimator unbiased? plot(density(first.obs-35), main="Bias of the Sample Mean") abline(v=mean(first.obs-35), lty=2, col="firebrick") plot(density(first.obs-35), main="Comparison of the two estimators") abline(v=mean(first.obs-35), lty=2, col="firebrick") points(density(means-35), type="l") rmse <- rep(NA, 10000) for(i in 1:10000){ rmse[i] <- sqrt(var(binomial.sample[i,]) + mean( binomial.sample[i,]-.525 )^2) } plot(1:10000, rmse, main="RMSE of Sample Mean", xlab="Number of observation", type="l") plot(100:10000, rmse[100:10000], main="RMSE of Sample Mean (100-10000", xlab="Number of observation", type="l")