plotPvalues <- function() {
  l.norm <- rnorm(5000)
  shap.n <- function(n) {
    replicate(2000, unlist(shapiro.test(sample(l.norm, n, replace = TRUE))[c(1,2)]))
  }
  shap.l <- mapply(shap.n, seq(5, 50, 5), SIMPLIFY = FALSE)
  plot(seq(0,1,length.out = 10), seq(0.7, 1, length.out = 10), type = "n",
       xlab = "P",
       ylab = "W",
       main = "Simulated Empirical CDF for n = 5, 10, 15, ..., 50") 
  abline(v = seq(0, 1, length.out = 11), lty = 3, col = "grey")
  abline(h = seq(0.7, 1, length.out = 7), lty = 3, col = "grey")
  text(rep(0.21, 4), c(0.84, 0.88, 0.91, 0.98), labels = c("n = 5", "n = 10", "n = 15", "n = 50"), cex = 0.9)
  for (n in 1:10) {
    lines(sort(shap.l[[n]][2, ]), sort(shap.l[[n]][1, ]), col = rgb(red = (1 - 5*n/50), green = 0, blue = 5*n/50))
  }
}