ps = seq(0.01, 0.99, length.out=100)
ys = dbinom(12, 18, ps)

plot(ps, ys, type='l', xlab='p', ylab='dbinom(X, N, p)')
abline(v=2/3, col='red')
grid()

X = 12
N = 18
phat = X/N
ci = phat + c(-1.96*sqrt(phat*(1-phat)/N), 1.96*sqrt(phat*(1-phat)/N))
print(ci)

res = prop.test(12, 18)
print(res)