# Setup
source("https://pbreheny.org/7110/f26/notes/fun.R")
library(numDeriv)

# Quadratic approximation: Setup
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"
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))

# Information
set.seed(1)
ylim <- c(-10, 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)))

# Information of maximum
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 simulation #1
n <- 10
x <- 8
theta <- seq(0.45, 0.97, len = 499)
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(0.025) * se <= 0.8) & (p + qnorm(0.975) * se >= 0.8)

# Binomial simulation #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(0.025) * se <= 0.8) & (p + qnorm(0.975) * se >= 0.8)

# Binomial simulation #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(0.025) * se <= 0.8) & (p + qnorm(0.975) * se >= 0.8)