%%R
y <- c(35,34,38,35,37)
y

%%R
n <- length(y)
n

%%R
estimate <- mean(y)
estimate

%%R
se <- sd(y)/sqrt(n)
se

%%R
int.50 <- estimate + qt(c(.25,.75),n-1)*se
int.50

%%R
int.95 <- estimate + qt(c(.025,.975),n-1)*se
int.95

"""
Explanation: 2.3: Classical confidence intervals
CI for continuous data, Pg 18
End of explanation
"""


%%R
y <- 700
y

%%R
n <- 1000
n

%%R
estimate <- y/n
estimate

%%R
se <- sqrt(estimate*(1-estimate)/n)
se

%%R
int.95 <- estimate + qnorm(c(.025,.975))*se
int.95

"""
Explanation: CI for proportions, Pg 18
End of explanation
"""


%%R
y <- rep(c(0,1,2,3,4), c(600,300,50,30,20))
y

%%R
n <- length(y)
n

%%R
estimate <- mean(y)
estimate

%%R
se <- sd(y)/sqrt(n)
se

%%R
int.50 <- estimate + qt(c(.25,.75),n-1)*se
int.50

%%R
int.95 <- estimate + qt(c(.025,.975),n-1)*se
int.95

"""
Explanation: CI for discrete data, Pg 18
End of explanation
"""


%%R
# Data is available in death.polls directory of ARM_Data
# Note: data seems to come from Gallup
# See: http://www.gallup.com/poll/1606/Death-Penalty.aspx
#
# Columns are
# 1- year
# 2- month
# 3- percentage support
# 4- percentage against
# 5- percentage no opinion
polls <- matrix(scan("../../ARM_Data/death.polls/polls.dat"),
                ncol=5, byrow=TRUE)
polls

%%R
# Note: this is the percentage support **among those that have an opinion** (pg 18)
#
# -- "No opinion" percentages are ignored
# -- !s a result, our plot below looks different from the Gallup plot
support <- polls[,3]/(polls[,3]+polls[,4])
support

%%R
year <-  polls[,1] + (polls[,2]-6)/12
year

%%R
# confidence intervals in plot assume sample size N=1000
# itervals are 68% confidence
par (mar=c(5,5,4,2)+.1)

plot(year, support*100, xlab="Year", ylim=c(min(100*support)-1, max(100*support)+1),
     ylab="Percentage support for the death penalty", cex=1.1, cex.main=1.2,
     cex.axis=1.1, cex.lab=1.1, pch=20)

for (i in 1:nrow(polls))
  lines (rep(year[i],2), 100*(support[i]+c(-1,1)*sqrt(support[i]*(1-support[i])/1000)))

"""
Explanation: Plot Figure 2.3, Pg 19
End of explanation
"""


%%R
N <- c(66030000, 81083600, 60788845)
p <- c(0.55, 0.61, 0.38)
se <- c(0.02, 0.03, 0.03)

%%R
w.avg <- sum(N*p)/sum(N)
w.avg

%%R
se.w.avg <- sqrt(sum((N*se/sum(N))^2))
se.w.avg

%%R
# this uses +/- 2 std devs
int.95 <- w.avg + c(-2,2)*se.w.avg
int.95

"""
Explanation: Weighted averages, Pg 19
The example R-code for this part is incomplete, so I'll make up N, p and
se loosely related to the text on page 19.
End of explanation
"""


%%R
n.men <- 500
n.men

%%R
p.hat.men <- 0.75
p.hat.men

%%R
se.men <- sqrt (p.hat.men*(1-p.hat.men)/n.men)
se.men

%%R
n.women <- 500
n.women

%%R
p.hat.women <- 0.65
p.hat.women

%%R
se.women <- sqrt(p.hat.women*(1-p.hat.women)/n.women)
se.women

%%R
n.sims <- 10000
n.sims

%%R
p.men <- rnorm(n.sims, p.hat.men, se.men)
p.men[1:10] # show first ten

%%R
p.women <- rnorm(n.sims, p.hat.women, se.women)
p.women[1:10] # show first ten

%%R
ratio <- p.men/p.women
ratio[1:10] # show first ten

%%R
int.95 <- quantile(ratio, c(.025,.975))
int.95

"""
Explanation: CI using simulations, Pg 20
End of explanation
"""


