# Confidence Intervals

# initiate constants
n.samples <- 100
n <- 100
p1 <- .28

# create samples and histogram of phats
distrib1 <- matrix(rbinom(n.samples*n, size=1, 
	prob=p1), ncol=n)
sample1 <- rowMeans(distrib1[,1:n])
hist(sample1)

# compute std. dev. of histogram (2 std errors)
summary(sample1)
sd(sample1)
diff <- 1.96 * sqrt(p1*(1-p1)/n)
lower <- p1 - diff
upper <- p1 + diff
print("95% CI is")
lower ; upper
hist(sample1)
abline(v=lower,lwd=4)
abline(v=upper,lwd=4)


# count proportion of CI's not capturing p
diff1 <- 1.96 * sqrt(sample1 * (1 - sample1)/n)
lower1 <- sample1 - diff1
upper1 <- sample1 + diff1
test1 <- which(lower1>p1 | upper1<p1)
length(test1)
cat("proportion of ci's not capturing", 
	p1, "are", length(test1)/n.samples,"\n")
cat("sample interval numbers not capturing p are",
	which(lower1>p1 | upper1< p1), "\n")


# plot ci's 

x1 <- 1:n.samples
y1 <- lower1
y2 <- upper1
x <- c(x1,x1)
y <- c(y1,y2)
plot(x,y)
abline(p1,0)
lines(x1,y1)
lines(x1,y2)