simulations/sim.R

In chrisroberts2112/msProject: Functions for evaluating Rivest Estimator

#
#
# library(caribouSim)
#
#
# ########################################################################################
#
# #################################### SIMULATIONS #######################################
#
# ########################################################################################
#
#
# T <- 1000000
# n <- seq(.00001*T, .001*T, by = 10)
# iters <- 100
#
#
# ###################### phi = 0 ######################
#
#
# results <- vector(mode = "list", length = length(n))
# rejections <- matrix(nrow = length(n), ncol = 3)
# for(i in seq(n)){
#   tmp <- iterate(iters = iters, T = T, n = n[i], phi = 0, r = .5, modelSim = "I")
#   print(tmp)
#   results[[i]] <- cbind(i, tmp$collars, tmp$abundance, tmp$p_values, tmp$rhat)
#   rejections[i,] <- tmp$reject_rate
# }
#
# # write.csv(do.call(rbind, results), "simulations/results_phi0d.csv")
# # write.csv(rejections, "simulations/rejections_phi0d.csv")
#
# rejections_phi0 <- read.csv("simulations/rejections_phi0d.csv")
# plot_sim(dat = rejections_phi0)
#
# ###################### phi = 0.1 ######################
#
# results <- vector(mode = "list", length = length(n))
# rejections <- matrix(nrow = length(n), ncol = 3)
# for(i in seq(n)){
#   tmp <- iterate(iters = iters, T = T, n = n[i], phi = 0.1, r = .5, modelSim = "I")
#   print(tmp)
#   results[[i]] <- cbind(i, tmp$collars, tmp$abundance, tmp$p_values, tmp$rhat)
#   rejections[i,] <- tmp$reject_rate
# }
#
# # write.csv(do.call(rbind, results), "simulations/results_phi1d.csv")
# # write.csv(rejections, "simulations/rejections_phi1d.csv")
#
# rejections_phi1 <- read.csv("simulations/rejections_phi1d.csv")
# plot_sim(dat = rejections_phi1)
#
# ###################### phi = 0.2 ######################
#
#   results <- vector(mode = "list", length = length(n))
#   rejections <- matrix(nrow = length(n), ncol = 3)
#   for(i in seq(n)){
#     tmp <- iterate(iters = iters, T = T, n = n[i], phi = 0.2, r = .5, modelSim = "I")
#     print(tmp)
#     results[[i]] <- cbind(i, tmp$collars, tmp$abundance, tmp$p_values, tmp$rhat)
#     rejections[i,] <- tmp$reject_rate
#   }
#
#   write.csv(do.call(rbind, results), "simulations/results_phi2d.csv")
#   write.csv(rejections, "simulations/rejections_phi2d.csv")
#
#   rejections_phi2 <- read.csv("simulations/rejections_phi2d.csv")
#   plot_sim(dat = rejections_phi2) + ylim(0,1)
#
#
#   ###################### phi = 0.3 ######################
#
#   results <- vector(mode = "list", length = length(n))
#   rejections <- matrix(nrow = length(n), ncol = 3)
#   for(i in seq(n)){
#     tmp <- iterate(iters = iters, T = T, n = n[i], phi = 0.3, r = .5, modelSim = "I")
#     print(tmp)
#     results[[i]] <- cbind(i, tmp$collars, tmp$abundance, tmp$p_values, tmp$rhat)
#     rejections[i,] <- tmp$reject_rate
#   }
#
#   write.csv(do.call(rbind, results), "simulations/results_phi3d.csv")
#   write.csv(rejections, "simulations/rejections_phi3d.csv")
#
#   rejections_phi3 <- read.csv("simulations/rejections_phi3d.csv")
#   plot_sim(dat = rejections_phi3) + ylim(0,1)
#
#
#   ###################### phi = 0.4 ######################
#
#   results <- vector(mode = "list", length = length(n))
#   rejections <- matrix(nrow = length(n), ncol = 3)
#   for(i in seq(n)){
#     tmp <- iterate(iters = iters, T = T, n = n[i], phi = 0.4, r = .5, modelSim = "I")
#     print(tmp)
#     results[[i]] <- cbind(i, tmp$collars, tmp$abundance, tmp$p_values, tmp$rhat)
#     rejections[i,] <- tmp$reject_rate
#   }
#
#   write.csv(do.call(rbind, results), "simulations/results_phi4d.csv")
#   write.csv(rejections, "simulations/rejections_phi4d.csv")
#
#   rejections_phi4 <- read.csv("simulations/rejections_phi4d.csv")
#   plot_sim(dat = rejections_phi4) + ylim(0,1)
#
# ###################### phi = 0.5 ######################
#
# results <- vector(mode = "list", length = length(n))
# rejections <- matrix(nrow = length(n), ncol = 3)
# for(i in seq(n)){
#   tmp <- iterate(iters = iters, T = T, n = n[i], phi = 0.5, r = .5, modelSim = "I")
#   print(tmp)
#   results[[i]] <- cbind(i, tmp$collars, tmp$abundance, tmp$p_values, tmp$rhat)
#   rejections[i,] <- tmp$reject_rate
# }
#
# write.csv(do.call(rbind, results), "simulations/results_phi5d.csv")
# write.csv(rejections, "simulations/rejections_phi5d.csv")
#
# rejections_phi5 <- read.csv("simulations/rejections_phi5d.csv")
# plot_sim(dat = rejections_phi5) + ylim(0,1)
#
#
#
#
#
###################### phi = 0.7 ######################

Theta <- 1000000
nu <- seq(10, 250, by = 10)
iters <- 100

results <- vector(mode = "list", length = length(nu))
rejections <- matrix(nrow = length(nu), ncol = 3)
for(i in seq(nu)){
  tmp <- iterate(iters = iters, Theta = Theta, nu = nu[i], phi = .7, r = .5, modelSim = "I")
  print(tmp)
  results[[i]] <- cbind(i, tmp$collars, tmp$abundance, tmp$p_values, tmp$rhat)
  rejections[i,] <- tmp$reject_rate
}

# write.csv(do.call(rbind, results), "simulations/results_phi7.csv")
# write.csv(rejections, "simulations/rejections_phi7.csv")

# write.csv(phi7, "simulations/results_phi7.csv")
# read.csv("simulations/results_phi7.csv")


nu <- seq(10, 250, by = 10)
means <- data.frame(nu = nu,
                    Ttrue = tapply(phi7[,"Ttrue"], INDEX = factor(phi7$i), FUN = mean),
                    That = tapply(phi7[,"That"], INDEX = factor(phi7$i), FUN = mean))

cols <- c("Estimated" = "#f04546", "True" = "#3591d1")
library(ggplot2)
ggplot(data = means, aes(x = nu)) +
  geom_line(aes(y = That/1000000, color = names(cols)[1])) +
  geom_line(aes(y = Ttrue/1000000, color = names(cols)[2])) +
  ylab("Mean Population Size (Millions)") + xlab(bquote(paste(nu, " (nu)"))) +
  scale_y_continuous(limits = c(1, 1.3)) +
  scale_colour_manual(name = "", values = cols) +
  theme_bw() +
  theme(legend.position = "bottom")


###################### phi = 1 ######################

Theta <- 1000000
nu <- seq(160, 250, by = 10)
iters <- 100

results <- vector(mode = "list", length = length(nu))
rejections <- matrix(nrow = length(nu), ncol = 3)
for(i in seq(nu)){
  tmp <- iterate(iters = iters, Theta = Theta, nu = nu[i], phi = 1, r = .5, modelSim = "I")
  print(tmp)
  results[[i]] <- cbind(i, tmp$collars, tmp$abundance, tmp$p_values, tmp$rhat)
  rejections[i,] <- tmp$reject_rate
}


# phi10b <- do.call(rbind, results)

write.csv(phi10b, "simulations/results_phi10b.csv")
phi10b <- read.csv("simulations/results_phi10b.csv")

phi10b$i <- dplyr::recode(as.factor(phi10b$i), "1" = "16", "2" = "17", "3" = "18", "4" = "19",
       "5" = "20", "6" = "21", "7" = "22", "8" = "23", "9" = "24", "10" = "25")

phi10d <- rbind(phi10, phi10b)

# write.csv(phi10d, "simulations/results_phi10d.csv")
# write.csv(rejections, "simulations/rejections_phi10.csv")

# rejections_phi10 <- read.csv("simulations/rejections_phi10.csv")
phi10 <- na.omit(read.csv("simulations/results_phi10d.csv"))

nu <- seq(10, 250, by = 10)
means <- data.frame(nu = nu,
                    Ttrue = tapply(phi10[,"Ttrue"], INDEX = factor(phi10$i), FUN = mean),
                    That = tapply(phi10[,"That"], INDEX = factor(phi10$i), FUN = mean))

cols <- c("Estimated" = "#f04546", "True" = "#3591d1")
library(ggplot2)
ggplot(data = means, aes(x = nu)) +
  geom_line(aes(y = That/1000000, color = names(cols)[1])) +
  geom_line(aes(y = Ttrue/1000000, color = names(cols)[2])) +
  ylab("Mean Population Size (Millions)") + xlab(bquote(paste(nu, " (nu)"))) +
  scale_y_continuous(limits = c(0.5, 1.3)) +
  scale_colour_manual(name = "", values = cols) +
  geom_text(aes(x = 200, y = 1.25,
                label = paste("Avg difference =",
                              round(mean((phi10$That - phi10$Ttrue)/ phi10$Ttrue), 4)*100,
                              "%")), cex = 5) +
  theme_bw() +
  theme(legend.position = "bottom")