theta =.7
x = sample(c(1,0),25, prob=c(theta, 1-theta),replace=T)
#
#likelihood
#
p=seq(0,1,.01)
n=25
t = 18
lik=p^t*(1-p)^(n-t)
plot(p,lik)
m = which.max(lik)
m
p[m]

#binomial dist

prob = dbinom(0:25,25,theta)
plot(0:25/25,prob)

#CI from normal approx
#that = sum(x)/n
that = 18/25
shat = sqrt(theta*(1-theta)/n)
shat
that - 1.96*shat
that + 1.96*shat

#bootstrap

that=.72
tstar = rep(0,1000)
for (i in 1:1000)
tstar[i] = mean(sample(c(1,0),25, prob=c(that, 1-that),replace=T))
hist(tstar)
mean(tstar)
sd(tstar)

#bootstrap CI

hist(tstar - that)
q = quantile(tstar-that,probs = c(.025, .975))
q
lower=that-q[2]
upper = that-q[1]
lower
upper