## —-echo=FALSE———————————————————-
set.seed(2312)
## ————————————————————————
pbinom(3.5, 10, 0.5)
pbinom(3, 10, 0.5) # same answer as above! (why?)
## ————————————————————————
1 – pbinom(3.5, 10, 0.25)
## ————————————————————————
T <- rbinom(200, 10, 0.5) # simulate under H0
alpha <- sum(T < 3.5) / length(T) # estimate
alpha
## ------------------------------------------------------------------------
T1 <- rbinom(200, 10, 0.25) # simulate under H1
beta <- sum(T1 > 3.5) / length(T1) # estimate
beta
## ————————————————————————
alpha + c(-1 ,1) * 1.96 * sqrt(alpha * (1 – alpha) / 200) # CI for alpha
beta + c(-1, 1) * 1.96 * sqrt(beta * (1 – beta) / 200) # CI for beta
## —-fig.width=6, fig.height=4——————————————-
K1 <- function(p)
1 - pbinom(12, 25, p)
p <- seq(0, 1, 0.01)
K <- K1(p)
plot(p, K, type = "l", ylab = "Power, K(p)")
## ------------------------------------------------------------------------
K1(0.4)
## ------------------------------------------------------------------------
1 - K1(0.6)
## ----fig.width=6, fig.height=4-------------------------------------------
K2 <- function(p)
1 - pbinom(12, 30, p)
curve(K1, from = 0, to = 1, xlab = "p", ylab = expression(1 - beta))
curve(K2, from = 0, to = 1, add = TRUE, col = 2, lty = 2)
## ------------------------------------------------------------------------
data <- read.table("candies.txt", header = TRUE) # load the data
## ------------------------------------------------------------------------
y <- data[, 1]
n <- data[, 2]
p <- y / n
z <- abs(p - 0.2) / sqrt(0.2 * 0.8 / n)
sum(z > 1.96) / length(z) # proportion
which(z > 1.96) # this shows *which* students rejected the null
## ————————————————————————
b1 <- p - 1.96 * sqrt(p * (1 - p) / n)
b2 <- p + 1.96 * sqrt(p * (1 - p) / n)
sum((b1 <= 0.2) & (0.2 <= b2)) / length(p)
## ------------------------------------------------------------------------
x <- sum(y)
N <- sum(n)
prop.test(x, N, p = 0.2, alternative = "two.sided")
## ------------------------------------------------------------------------
x <- c( 25, 50) # successes
n <- c(351, 389) # sample sizes
prop.test(x, n, alternative = "two.sided")