simulations/analysis.R

In chrisroberts2112/msProject: Functions for evaluating Rivest Estimator

#
#
# ## create power curves ##
#
# # individual plots
#
# rejections_phi0 <- read.csv("simulations/rejections_phi0d.csv")
# plot_sim(dat = rejections_phi0, legend = "bottom")
#
# rejections_phi1 <- read.csv("simulations/rejections_phi1d.csv")
# plot_sim(dat = rejections_phi1)
#
# rejections_phi2 <- read.csv("simulations/rejections_phi2d.csv")
# plot_sim(dat = rejections_phi2)
#
# rejections_phi3 <- read.csv("simulations/rejections_phi3d.csv")
# plot_sim(dat = rejections_phi3)
#
# rejections_phi4 <- read.csv("simulations/rejections_phi4d.csv")
# plot_sim(dat = rejections_phi4)
#
# rejections_phi5 <- read.csv("simulations/rejections_phi5d.csv")
# plot_sim(dat = rejections_phi5)
#
#
# # plot assumption violated plots togather
#
# plot1 <- plot_sim(dat = rejections_phi1) + labs(title = expression(paste(phi, " = 0.1")))
# plot2 <- plot_sim(dat = rejections_phi2) + labs(title = expression(paste(phi, " = 0.2")))
# plot3 <- plot_sim(dat = rejections_phi3) + labs(title = expression(paste(phi, " = 0.3")))
# plot4 <- plot_sim(dat = rejections_phi4) + labs(title = expression(paste(phi, " = 0.4")))
# plot5 <- plot_sim(dat = rejections_phi5) + labs(title = expression(paste(phi, " = 0.5")))
#
# library(cowplot)
#
# plots <- plot_grid(plot1, plot2, plot3, plot4, plot5, ncol = 1)
# legend <- get_legend(plot_sim(dat = rejections_phi1, legend = "bottom"))
#
# plot_grid(plots, legend, rel_heights = c(1, .2), ncol = 1)
#
# round(mean(rejections_phi0$V3), 4)*100
#
#
# # Ni for 50% and 80% power for FET
#
# Ni <- seq(10, 1000, by = 10)
#
# plot_sim(dat = rejections_phi5, legend = "bottom")
#
# tmp <- as.data.frame(sapply(rejections_phi5, FUN = f))
# power1 <- Ni[c(nearest(tmp$V1, .5), nearest(tmp$V1, .8),
#                nearest(tmp$V2, .5), nearest(tmp$V2, .8),
#                nearest(tmp$V3, .5), nearest(tmp$V3, .8))]
#
#
# vec <- tmp$V1
# vec[1:length(vec)]
#
# f <- function(vec){
#   out <- zoo::rollmean(vec, k = 7)
#   return(out)
# }
#
# MAs <- rbind(NA, NA, NA,
#              as.data.frame(sapply(dat, FUN = function(x) zoo::rollmean(x, k = 7))),
#              NA, NA, NA)
#
#
#
# tmp <- rbind(NA, NA, NA, as.data.frame(sapply(rejections_phi0, FUN = f)), NA, NA, NA)
# plot_sim(dat = rejections_phi0, legend = "bottom")
#
# library(ggplot2)
# ggplot(data = as.data.frame(sapply(rejections_phi0, FUN = f)), aes(x = X)) +
#   geom_line(aes(y = V1)) +
#   geom_line(aes(y = V2)) +
#   geom_line(aes(y = V3))
#
#
#
# movingAvg <- function(vec){
#   vec[1:length(vec)]
# }
#
# power1 <- Ni[c(nearest(zoo::rollmean(rejections_phi1$V2, k = 9), .5), nearest(rejections_phi1$V3, .8))]
# power2 <- Ni[c(nearest(zoo::rollmean(rejections_phi2$V2, k = 9), .5), nearest(rejections_phi2$V3, .8))]
# power3 <- Ni[c(nearest(zoo::rollmean(rejections_phi3$V2, k = 9), .5), nearest(rejections_phi3$V3, .8))]
# power4 <- Ni[c(nearest(zoo::rollmean(rejections_phi4$V2, k = 9), .5), nearest(rejections_phi4$V3, .8))]
# power5 <- Ni[c(nearest(zoo::rollmean(rejections_phi5$V2, k = 9), .5), nearest(rejections_phi5$V3, .8))]
#
# rbind(power1, power2, power3, power4, power5)
#
#
#
# ## estimate accuracy ##
#
# # phi = 0
#
# library(ggplot2)
# theme_set(theme_bw())
#
# phi0 <- read.csv("simulations/results_phi0c.csv")
#
# mean(abs(phi0$pctDiff))
#
#
# plot_pctDiff(phi0, phi = 0.0, prop = .3)
#
# ggplot(data = phi0) +
#   geom_line(aes(x = X.1, y = pctDiff)) +
#   scale_y_continuous(limits = c(-100, 100))
#
# phi0_means <- data.frame(Ni = seq(10, 1000, by = 10),
#                          Ttrue = tapply(phi0[,"Ttrue"], INDEX = factor(phi0$i), FUN = mean),
#                          That = tapply(phi0[,"That"], INDEX = factor(phi0$i), FUN = mean))
#
# sum((phi0_means$Ttrue - phi0_means$That)^2)
#
# ggplot(data = phi0_means, aes(x = Ni)) +
#   geom_line(aes(y = Ttrue), color = "blue") +
#   geom_line(aes(y = That), color = "red") +
#   scale_y_continuous(limits = c(1000000, 1500000))
#
# plot_meanTs(phi0, legend = c(.8, .8))
#
#
# tmp <- na.omit(phi0$pctDiff[abs(phi0$pctDiff) < 50])
# hist(tmp, breaks = 50)
# sum(tmp^2)
#
#
# # phi = 1
#
#
#
# phi1 <- na.omit(read.csv("simulations/results_phi1c.csv"))
# phi1_means <- data.frame(Ni = seq(10, 1000, by = 10),
#                          Ttrue = tapply(phi1[,"Ttrue"], INDEX = factor(phi1$i), FUN = mean),
#                          That = tapply(phi1[,"That"], INDEX = factor(phi1$i), FUN = mean))
#
#
#
# plot_pctDiff(phi1, phi = 0.1, prop = .3)
#
# ggplot(data = phi1_means, aes(x = Ni)) +
#   geom_line(aes(y = Ttrue), color = "blue") +
#   geom_line(aes(y = That), color = "red") +
#   scale_y_continuous(limits = c(1000000, 1500000))
#
#
# tmp <- na.omit(phi1$pctDiff[abs(phi1$pctDiff) < 100])
# hist(tmp, breaks = 50)
# sum(tmp^2)
#
#
# # phi = 2
#
# phi2 <- na.omit(read.csv("simulations/results_phi2c.csv"))
# phi2_means <- data.frame(Ni = seq(10, 1000, by = 10),
#                          Ttrue = tapply(phi2[,"Ttrue"], INDEX = factor(phi2$i), FUN = mean),
#                          That = tapply(phi2[,"That"], INDEX = factor(phi2$i), FUN = mean))
#
# ggplot(data = phi2) +
#   geom_line(aes(x = X.1, y = pctDiff)) +
#   scale_y_continuous(limits = c(-100, 100))
#
# ggplot(data = phi2_means, aes(x = Ni)) +
#   geom_line(aes(y = Ttrue), color = "blue") +
#   geom_line(aes(y = That), color = "red") +
#   scale_y_continuous(limits = c(1000000, 1500000))
#
#
# tmp <- na.omit(phi2$pctDiff[abs(phi2$pctDiff) < 50])
# hist(tmp, breaks = 50)
# sum(tmp^2)
#
#
#
# # phi = 3
#
# phi3 <- na.omit(read.csv("simulations/results_phi3c.csv"))
# phi3_means <- data.frame(Ni = seq(10, 1000, by = 10),
#                          Ttrue = tapply(phi3[,"Ttrue"], INDEX = factor(phi3$i), FUN = median),
#                          That = tapply(phi3[,"That"], INDEX = factor(phi3$i), FUN = median))
# sum((phi3_means$Ttrue - phi3_means$That)^2)
#
# ggplot(data = phi3) +
#   geom_line(aes(x = X.1, y = pctDiff)) +
#   scale_y_continuous(limits = c(-100, 100))
#
# ggplot(data = phi3_means, aes(x = Ni)) +
#   geom_line(aes(y = Ttrue), color = "blue") +
#   geom_line(aes(y = That), color = "red") +
#   scale_y_continuous(limits = c(1000000, 1500000))
#
# mean(abs(phi3$pctDiff))
#
# tmp <- na.omit(phi3$pctDiff[abs(phi3$pctDiff) < 50])
# hist(tmp, breaks = 50)
# sum(tmp^2)
#
# sum((phi3_means$Ttrue - phi3_means$That)^2)
#
# # phi = 4
#
# phi4 <- na.omit(read.csv("simulations/results_phi4c.csv"))
# phi4_means <- data.frame(Ni = seq(10, 1000, by = 10),
#                          Ttrue = tapply(phi4[,"Ttrue"], INDEX = factor(phi4$i), FUN = mean),
#                          That = tapply(phi4[,"That"], INDEX = factor(phi4$i), FUN = mean))
#
# ggplot(data = phi4) +
#   geom_line(aes(x = X.1, y = pctDiff)) +
#   scale_y_continuous(limits = c(-100, 100))
#
# ggplot(data = phi4_means, aes(x = Ni)) +
#   geom_line(aes(y = Ttrue), color = "blue") +
#   geom_line(aes(y = That), color = "red") +
#   scale_y_continuous(limits = c(1000000, 1500000))
#
# tmp <- na.omit(phi4$pctDiff[abs(phi4$pctDiff) < 50])
# hist(tmp, breaks = 50)
# sum(tmp^2)
#
#
#
# # phi = 5
#
# phi5 <- na.omit(read.csv("simulations/results_phi5c.csv"))
# phi5_means <- data.frame(Ni = seq(10, 1000, by = 10),
#                          Ttrue = tapply(phi5[,"Ttrue"], INDEX = factor(phi5$i), FUN = median),
#                          That = tapply(phi5[,"That"], INDEX = factor(phi5$i), FUN = median))
# sum((phi5_means$Ttrue - phi5_means$That)^2)
#
#
# ggplot(data = phi5) +
#   geom_line(aes(x = X.1, y = pctDiff)) +
#   scale_y_continuous(limits = c(-100, 100))
#
# ggplot(data = phi5_means, aes(x = Ni)) +
#   geom_line(aes(y = Ttrue), color = "blue") +
#   geom_line(aes(y = That), color = "red") +
#   scale_y_continuous(limits = c(1000000, 1500000))
#
# mean(abs(phi5$pctDiff))
#
# tmp <- na.omit(phi5$pctDiff)
# hist(tmp, breaks = 50)
# sum(tmp^2)
#
# sum((phi5_means$Ttrue - phi5_means$That)^2)
#
# pd0 <- na.omit(phi0$pctDiff)
# pd1 <- na.omit(phi1$pctDiff)
# pd2 <- na.omit(phi2$pctDiff)
# pd3 <- na.omit(phi3$pctDiff)
# pd4 <- na.omit(phi4$pctDiff)
# pd5 <- na.omit(phi5$pctDiff)
#
# plot(0:5, lapply(list(pd0, pd1, pd2, pd3, pd4, pd5), FUN = function(x) sum(x^2)))
#
# boxplot(pd0, pd1, pd2, pd3, pd4, pd5)
#
#
#
# wilcox.test(c(phi0_means$That - phi0_means$Ttrue),
#        c(phi5_means$That - phi5_means$Ttrue), paired = TRUE)
#
# diffs <- data.frame(diffs = c(phi0_means$That - phi0_means$Ttrue,
#                               phi1_means$That - phi1_means$Ttrue,
#                               phi2_means$That - phi2_means$Ttrue,
#                               phi3_means$That - phi3_means$Ttrue,
#                               phi4_means$That - phi4_means$Ttrue,
#                               phi5_means$That - phi5_means$Ttrue),
#                     phi = factor(rep(0:5, times = rep(100, 6))))
#
# kruskal.test(diffs ~ phi, data = diffs)
#
#
#
#
# ## plot_meanTs
#
# phi0 <- read.csv("simulations/results_phi0d.csv")
# phi1 <- read.csv("simulations/results_phi1d.csv")
# phi2 <- read.csv("simulations/results_phi2d.csv")
# phi3 <- read.csv("simulations/results_phi3d.csv")
# phi4 <- read.csv("simulations/results_phi4d.csv")
# phi5 <- read.csv("simulations/results_phi5d.csv")
#
# plot_meanTs(phi1, legend = c(.7, .2))
# plot_pctDiff(phi0, phi = "0", prop = .2)
#
# round(mean(phi0$pctDiff), 4)
#
# plot_meanTs(phi1, legend = c(.7, .2))
# plot_pctDiff(phi1, phi = "0.1", prop = .2)
#
# -round(mean(phi1$pctDiff), 4)
#
# plot_meanTs(phi2, legend = c(.7, .2))
# plot_pctDiff(phi2, phi = "2", prop = .2)
#
# round(mean(phi2$pctDiff), 4)
#
# plot_meanTs(phi3, legend = c(.7, .2))
# plot_pctDiff(phi3, phi = "3", prop = .2)
#
# round(mean(phi3$pctDiff), 4)
#
# plot_meanTs(phi4, legend = c(.7, .2))
# plot_pctDiff(phi4, phi = "4", prop = .2)
#
# round(mean(phi4$pctDiff), 4)
#
# plot_meanTs(phi5, legend = c(.7, .2))
# plot_pctDiff(phi5, phi = "5", prop = .2)
#
# round(mean(phi5$pctDiff), 4)
#
#
#
# # table 3
# tab <- data.frame(phi = 0:5, SSE = c(
#   sum((na.omit(phi0$Ttrue - phi0$That))^2),
#   sum((na.omit(phi1$Ttrue - phi1$That))^2),
#   sum((na.omit(phi2$Ttrue - phi2$That))^2),
#   sum((na.omit(phi3$Ttrue - phi3$That))^2),
#   sum((na.omit(phi4$Ttrue - phi4$That))^2),
#   sum((na.omit(phi5$Ttrue - phi5$That))^2))
# )
#
#
# ggplot(data = na.omit(phi0)) +
#   geom_line(aes(x = ntrue, y = Ttrue - That))
#
# mean(na.omit(phi5$pctDiff[phi5$ntrue < 200]))
#