source(file = "beta.R")

n <- 100
std <- 155
x_bar <- 262
mu <- 286

ub <- 292.38
lb <- 231.62 # rejection region is (-Inf, 231.62)

x <- (x_bar - n) : (x_bar + n)
y <- dnorm(x = x, mean = x_bar, sd = std / sqrt(n))
plot(x, y, type = 'n', xlab = "Calories", ylab = "Density", main = "Calorie Density Function")
lines(x, y, lwd = 2)
abline(v = mu, col = 'red', lty = 2, lwd = 3)
abline(v = lb, col = 'blue', lty = 2, lwd = 3)

x.shade <- min(x) : round(lb)
y.shade <- dnorm(x = x.shade, mean = x_bar, sd = std / sqrt(n))
polygon(c(0, x.shade), c(y.shade, -1), density = 20, angle = 50,  col = 'darkred')
legend("topright", legend = c("Population Mean", "Critical Z", "Rejection Region"), col = c("red", "blue", "darkred"), pch = c(16, 16, 16))


beta.x <- x
beta.y <- beta(beta.x, x_bar, std, n)
plot(beta.x, beta.y, type = 'n', xlab = "Alternative Mean", ylab = "Power", main = "Power Plot for Varying Alternative Means")
lines(beta.x, beta.y, lwd = 5, col = 'darkred')