X1 <- seq(-3,3, length = 300)

Y1 <- dnorm(X1)

plot(X1, Y1, type = "n", xlab=" ", ylab = " ")

x0 <- min(which(X1 >= -2.58))

x1 <- min(which(X1 >= -1.96))

x2 <- min(which(X1 >= -1))

x3 <- max(which(X1 <= 1))

x4 <- max(which(X1 <= 1.96))

x5 <- max(which(X1 <= 2.58))

polygon(x = c(X1[c(1, 1:x0, x0)]), y = c(0, Y1[1:x0], 0), col = "yellow", border = NA)

polygon(x = c(X1[c(x0, x0:x1, x1)]), y = c(0, Y1[x0:x1], 0), col = "yellow", border = NA)

polygon(x = c(X1[c(x1, x1:x2, x2)]), y = c(0, Y1[x1:x2], 0), col = "orange", border = NA)

polygon(x = c(X1[c(x2, x2:x3, x3)]), y = c(0, Y1[x2:x3], 0), col = "orange", border = NA)

polygon(x = c(X1[c(x3, x3:x4, x4)]), y = c(0, Y1[x3:x4], 0), col = "orange", border = NA)

polygon(x = c(X1[c(x4, x4:x5, x5)]), y = c(0, Y1[x4:x5], 0), col = "yellow", border = NA)

polygon(x = c(X1[c(x5, x5:300, 300)]), y = c(0, Y1[x5:300], 0), col = "yellow", border = NA)

abline(v = 0, lwd = 0.5, lty = 1)

points(X1, Y1, type = "l")

abline(v = 1.96, lwd = 1, lty = 2)

abline(v = -1.96, lwd = 1, lty = 2)

text(0.5, 0.3, "1-alpha = 0.95 so alpha = 0.05", cex = 1, col = "black")