source("http://myweb.uiowa.edu/pbreheny/7210/f15/notes/fun.R")

# The data
Data <- read.delim("http://myweb.uiowa.edu/pbreheny/data/Pike1966.txt")
d <- sum(Data$Death)
v <- sum(Data$Time/365.25)

# Definitely doesn't follow an exponential distribution
require(survival)
fit <- survfit(Surv(Time/365.25, Death) ~ 1, Data)
plot(fit, mark.time=FALSE, col=pal(2)[1], lwd=3, bty='n', las=1,
     xlab='Time on study (Years)', ylab='Survival', conf.int=FALSE)
tt <- seq(0, 1, len=199)
lines(tt, pexp(tt, d/v, low=FALSE), lwd=3, col=pal(2)[2])

# Setup
lam <- seq(0.5, 3, len=199)
l <- function(ll) {d*log(ll) - ll*v - d*log(d/v) + d}
baseplot <- function() {plot(lam, l(lam), col=pal(2)[1], lwd=3, type='l', bty='n',
     xlab=expression(lambda), ylab=expression(loglik(lambda)))}
baseplot()

# Score
baseplot()
lines(c(0, 2), l(1) + c(-1,1) * (d-v), lwd=3, col='gray')
mtext("Slope: 11")
d - v
sqrt(d)
(d-v)/sqrt(d) # z

# Wald
baseplot()
arrows(1, 0, d/v, 0, lwd=3, col='gray', angle=90, code=3, xpd=TRUE)
mtext(expression(hat(theta)-1: " "*0.44))
arrows(1, 0, d/v, 0, lwd=3, col='gray', angle=90, code=3, xpd=TRUE)
mtext(expression(hat(theta)-1: " "*0.44))
lines(lam, -(lam-d/v)^2/(2*d/v^2), col=pal(2)[2], lwd=3)
d/v-1
sqrt(d)/v
(d/v-1)/(sqrt(d)/v) # z

# LRT
baseplot()
arrows(d/v, 0, d/v, l(1), lwd=3, col='gray', angle=90, code=3, xpd=TRUE)
mtext(expression(loglik(hat(theta))-loglik(1): " "*2.14))
arrows(d/v, 0, d/v, l(1), lwd=3, col='gray', angle=90, code=3, xpd=TRUE)
mtext(expression(loglik(hat(theta))-loglik(1): " "*2.14))
abline(h=-qchisq(.95,1)/2, col=pal(2)[2], lwd=3)
l(1)
pchisq(-2*l(1), 1, low=FALSE)

# "Exact"
vv <- seq(0, 60, len=199)
plot(vv, dgamma(vv, shape=d, rate=1), type='l', lwd=3, bty='n', col=pal(2)[1],
     xlab='v', ylab='Density', las=1)
abline(v=v, lwd=3, col='gray')
pgamma(v, shape=d, rate=1)

# Bayes
require(caTools)
lam <- sort(c(seq(0, 3, len=199), 1))
nc <- trapz(lam, exp(l(lam)))
posterior <- exp(l(lam))/nc
prior <- rep(1/2.5, length(lam))
plot(lam, prior, col='gray', lwd=3, type='l', bty='n', ylim=c(0, 1.75),
     xlab=expression(lambda), ylab=expression(p(lambda)))
lines(lam, posterior, col=pal(2)[2], lwd=3)
ind <- which(lam <= 1)
trapz(lam[ind], posterior[ind])

# Gamma prior
prior <- dgamma(lam, 2, 2)
plot(lam, prior, col='gray', lwd=3, type='l', bty='n', ylim=c(0, 1.75),
     xlab=expression(lambda), ylab=expression(p(lambda)))
nc <- trapz(lam, exp(l(lam))*prior)
posterior <- exp(l(lam))*prior/nc
lines(lam, posterior, col=pal(2)[2], lwd=3)
ind <- which(lam <= 1)
trapz(lam[ind], posterior[ind])