blackjack <- c(1, 19, 30, 52, 73, 113, 159, 209, 222, 227, 228, 282, 309, 318, 387)

xobs <- diff(c(0, blackjack)) 
xobs

xbarobs <- mean(xobs) 
phat <- 1/xbarobs
phat

n <- length(xobs) 
Ihat <- xbarobs^3/(xbarobs - 1) 
Ihat

Lxobs <- phat - 1/sqrt(n*Ihat) * qnorm(0.975) 
Uxobs <- phat - 1/sqrt(n*Ihat) * qnorm(0.025) 
round(c(L = Lxobs, U = Uxobs), 3)


ell <- function(p) {
stopifnot(all(0 < p & p < 1)) 
n * (log(p) + (xbarobs - 1) * log(1 - p)) 
}
pp <- seq(from = 0.001, to = 0.1, length = 201)
plot(pp, ell(pp), type = "l", main = "Log-likelihood function", xlab = "Values for p", ylab = "", bty = "n") 


init <- 0.1
opt <- optim(par = init, fn = ell, method = "L-BFGS-B", lower = 0.001, upper = 0.999, control = list(fnscale = -1), hessian = TRUE)
opt

phatopt <- opt$par
nIhatopt<- -opt$hessian[1, 1]
Lopt <- phatopt - sqrt(1/nIhatopt) * qnorm(0.975) 
Uopt <- phatopt - sqrt(1/nIhatopt) * qnorm(0.025) 
round(c(L = Lopt, U = Uopt), 3)

tauhat <- 1/phat
Ihattau <- Ihat/(-1/phat^2)^2
Ltau <- tauhat - 1/sqrt(n*Ihattau) * qnorm(0.975)
Utau <- tauhat - 1/sqrt(n*Ihattau) * qnorm(0.025)
round(c(L = Ltau, U = Utau), 3)

pp <- seq(from = 0.01, to = 0.08, length = 201) 
plot(pp, ell(pp), type = "l", main = "Log-likelihood function", xlab = "Values for p", ylab = "", bty = "n")
abline(h = ell(phat) - 2, col = "red", lty = 2) 


CIp <- range(pp[ell(pp)>ell(phat) - 2]) 
round(CIp, 3)

elltau <- function(tau) {
ell(1/tau)
}
tt <- seq(from = 10, to = 60, length = 201) 
plot(tt, elltau(tt), type = "l", main = "Log-likelihood function for tau = 1/p", xlab = "Values for tau", ylab = "", bty = "n") 
abline(h = elltau(tauhat) - 2, col = "red", lty = 2) 


CItau <- range(tt[elltau(tt) > elltau(tauhat) - 2])
round(CItau, 1) 

round(rev(1/CIp), 1