library(e1071)
library(pbapply)
library(readODS)

set.seed(123)

Mode <- function (x) names(sort(-table(runs)))[1]

data <- read_ods("ElectoralCollege.ods")

stopifnot(!is.na(data$Ptrump) & data$Ptrump >= 0 & data$Ptrump <= 1)
stopifnot(sum(data$Votes) == 538)

one_run <- function () {
    ind <- runif(nrow(data)) <= data$Ptrump

    sum(ind * data$Votes)
}

runs <- pbreplicate(10000000, one_run())

cat(paste0("P(Trump) = ", mean(runs >= 269), "\n"))
cat(paste0("E(Trump EVs) = ", mean(runs), "\n"))
cat(paste0("Median(Trump EVs) = ", median(runs), "\n"))
cat(paste0("Mode(Trump EVs) = ", Mode(runs), "\n"))
cat(paste0("Skewness(Trump EVs) = ", skewness(runs), "\n"))

# plots

svg("plot1.svg", width=10, height=6)

hist(runs, breaks=50, probability=TRUE,
     main="Distribution of Trump Electoral College Votes",
     xlab="Electoral Votes", col="lightblue", border="white")

lines(density(runs), col="blue", lwd=2)

abline(v=269, col="red", lwd=2, lty=2)

dev.off()

svg("plot2.svg", width=10, height=6)

plot(ecdf(runs), do.points=F, verticals=T,
     main="Cumulative Distribution of Trump Electoral College Votes",
     xlab="Electoral Votes", ylab="Cumulative Probability")

abline(v=269, col="red", lwd=2, lty=2)

dev.off()