# setup
source('https://myweb.uiowa.edu/pbreheny/7110/f23/notes/fun.R')
library(numDeriv)

# Quadratic approximation, Setup
par(mfrow=c(1,2), mar=c(4.3, 4.3, 1, 0.5))
nb <- 20
xb <- 6
tb <- seq(0, 1, len=499)
m <- 30
k <- 45
xh <- 14
th <- 45:200

# Plot likelihood
plotL(tb, dbinom(xb, nb, tb))
mtext('Binomial (n=20, x=6)')
plotL(th, dhyper(xh, m, th-m, k))
mtext('Hypergeometric (m=30, x=14, k=45)')

# Plot log-likelihood
plotl(tb, dbinom(xb, nb, tb, log=TRUE))
mtext('Binomial (n=20, x=6)')
plotl(th, dhyper(xh, m, th-m, k, log=TRUE))
mtext('Hypergeometric (m=30, x=14, k=45)')

# Zoom in
tb <- seq(0.1, 0.55, len=499)
th <- 70:160
plotl(tb, dbinom(xb, nb, tb, log=TRUE))
mtext('Binomial (n=20, x=6)')
plotl(th, dhyper(xh, m, th-m, k, log=TRUE))
mtext('Hypergeometric (m=30, x=14, k=45)')

# Zoom in
plotl(tb, dbinom(xb, nb, tb, log=TRUE))
mtext('Binomial (n=20, x=6)')
Info <- nb/(6/20*14/20)
lines(tb, -Info*(tb-6/20)^2/2, lty=2)
lh <- function(theta) dhyper(xh, m, theta-m, k, log=TRUE)
plotl(th, lh(th))
mtext('Hypergeometric (m=30, x=14, k=45)')
mle <- th[which.max(lh(th))]
Hess <- ((lh(mle+1) - lh(mle)) - (lh(mle) - lh(mle-1)))
lines(th, Hess*(th-mle)^2/2, lty=2)

# Random score functions
set.seed(5)
N <- 20
col <- "#4D4D4D66"
par(mfrow=c(2,2), mar=c(4, 4, 1, 0))
tt <- seq(2, 8, len=99)
S <- sapply(rnorm(N, 4, 1/sqrt(10)), function(x) 10*(x-tt)/1)
matplot(tt, S, type='l', lty=1, las=1, bty='n', col=col, xlab=expression(theta), ylab='Score')
abline(h=0, lty=2, col='slateblue'); abline(v=4, lty=2, col='slateblue')
mtext('Normal')
tt <- seq(2, 8, len=99)
S <- matrix(rpois(N*10, 4), 10, N) |> apply(2, mean) |> sapply(function(x) 10*(x-tt)/tt)
matplot(tt, S, type='l', lty=1, las=1, bty='n', col=col, xlab=expression(theta), ylab='Score')
abline(h=0, lty=2, col='slateblue'); abline(v=4, lty=2, col='slateblue')
mtext('Poisson')
tt <- seq(0.1, 0.8, len=99)
S <- sapply(rbinom(N, 10, 0.4)/10, function(x) {10*(x-tt)/(tt*(1-tt))})
matplot(tt, S, type='l', lty=1, las=1, bty='n', col=col, xlab=expression(theta), ylab='Score')
abline(h=0, lty=2, col='slateblue'); abline(v=0.4, lty=2, col='slateblue')
mtext('Binomial')
tt <- seq(-14, 24, len=99)
s <- function(x, tt) {
  out <- numeric(length(tt))
  for (j in 1:length(tt)) {
    out[j] <- sum((2*(x-tt[j]))/(1 + (x-tt[j])^2))
  }
  out
}
S <- matrix(rcauchy(N*10, 4), 10, N) |> apply(2, s, tt=tt)
matplot(tt, S, type='l', lty=1, las=1, bty='n', col=col, xlab=expression(theta), ylab='Score')
abline(h=0, lty=2, col='slateblue'); abline(v=4, lty=2, col='slateblue')
mtext('Cauchy')

# Cauchy likelihood
tt <- seq(-14, 24, len=99)
l <- rcauchy(10, 4) |> sapply(function(x) dcauchy(x, location=tt, log=TRUE)) |> apply(1, sum)
plotl(tt, l)
plotL(tt, exp(l))

# out.width='90%'
# Information
set.seed(1)
ylim <- c(-10, 0)
par(mfrow=c(1,2), mar=c(4.5, 4.5, 1, 0))
tt <- seq(2, 8, len=99)
x <- rpois(10, 4)
l1 <- sapply(x, function(x) dpois(x, tt, log=TRUE)) |> apply(1, sum)
plotl(tt, l1, ylim=ylim)
Info <- 10 / mean(x)
mtext(paste0('n=10, Info=', round(Info, 1)))
tt <- seq(2, 8, len=99)
l2 <- rpois(40, 4) |> sapply(function(x) dpois(x, tt, log=TRUE)) |> apply(1, sum)
plotl(tt, l2, ylim=ylim)
Info <- 40 / mean(x)
mtext(paste0('n=40, Info=', round(Info, 1)))

# Info of max
l <- function(theta) {4*pnorm(3.5-theta, log=TRUE) + dnorm(3.5-theta, log=TRUE)}
mle <- optimize(l, c(0, 10), maximum=TRUE)$maximum
grad(l, mle)  # Of course
(h <- hessian(l, mle) |> drop())

# Binomial: Sim 1
n <- 10
x <- 8
theta <- seq(0.45, 0.97, len=499)
par(mar=c(4.3, 4.3, 1, 0.5))
plotl(theta, dbinom(x, n, theta, log=TRUE))
mtext('Binomial (n=10, x=8)')
p <- x/n
Info <- n/(p*(1-p))
lines(theta, -Info*(theta-x/n)^2/2, lwd=2, lty=2)
abline(h=-qchisq(0.95, 1)/2, col='gray', lwd=3)
N <- 10000
x <- rbinom(N, n, 0.8)
p <- x/n
se <- sqrt(p*(1-p)/n)
Cov <- (p + qnorm(.025)*se <= 0.8) & (p + qnorm(.975)*se >= 0.8)

# Binomial: Sim 2
n <- 100
x <- 80
theta <- seq(0.7, 0.88, len=499)
par(mar=c(4.3, 4.3, 1, 0.5))
plotl(theta, dbinom(x, n, theta, log=TRUE))
mtext('Binomial (n=100, x=80)')
p <- x/n
Info <- n/(p*(1-p))
lines(theta, -Info*(theta-x/n)^2/2, lwd=2, lty=2)
abline(h=-qchisq(0.95, 1)/2, col='gray', lwd=3)
N <- 10000
x <- rbinom(N, n, 0.8)
p <- x/n
se <- sqrt(p*(1-p)/n)
Cov <- (p + qnorm(.025)*se <= 0.8) & (p + qnorm(.975)*se >= 0.8)

# Binomial: Sim 3
n <- 1000
x <- 800
theta <- seq(0.77, 0.83, len=499)
par(mar=c(4.3, 4.3, 1, 0.5))
plotl(theta, dbinom(x, n, theta, log=TRUE))
mtext('Binomial (n=1000, x=800)')
p <- x/n
Info <- n/(p*(1-p))
lines(theta, -Info*(theta-x/n)^2/2, lwd=2, lty=2)
abline(h=-qchisq(0.95, 1)/2, col='gray', lwd=3)
N <- 10000
x <- rbinom(N, n, 0.8)
p <- x/n
se <- sqrt(p*(1-p)/n)
Cov <- (p + qnorm(.025)*se <= 0.8) & (p + qnorm(.975)*se >= 0.8)

# AIC vs chisq
par(mar=c(4.3, 4.3, 1, 0.5))
d <- 1:20
Cutoff <- cbind(qchisq(0.95, d), 2*d)
matplot(d, Cutoff, lty=1, type='l', bty='n', las=1, lwd=3, col=pal(2), ylim=c(0, 40), xlim=c(0, 20), ylab="Cutoff / Penalty")
toplegend(legend=c(expression(chi^2), '2p'), col=pal(2), lwd=3)