# Conditional likelihood for Poisson
x <- 10
y <- 20
poisson.test(c(y, x))
binom.test(20, 30)
binom.test(20, 30)$conf.int |> binomial()$linkfun() |> exp()

# Setup
source('https://myweb.uiowa.edu/pbreheny/7110/f25/notes/fun.R')

# Loss of efficiency
set.seed(3)
library(data.table)
opt <- list(
Likelihood = c("profile", "conditional"),
rep = 1:1000,
n = 20,
sd = 10
)
out <- expand.grid(opt) |> as.data.table()
out[, est := NA_real_]

# Loop
for (i in 1:nrow(out)) {
o <- out[i]
if (i == 1 || o$rep != out$rep[i-1]) {
  x <- runif(o$n) |> rep(each = 2)
  y <- x + rnorm(length(x), sd = o$sd)
}

# Analyze
if (o$Likelihood == "profile") {
  fit <- lm(y ~ x)
  out[i, est := mean(residuals(fit)^2)]
} else {
  ymat <- matrix(y, nrow = 2)
  d <- ymat[1,] - ymat[2,]
  out[i, est := mean(d^2) / 2]
}  
}

# Table
out[, err := sqrt(est) - sd]
out[, .(Bias = mean(err), Variance = var(err), MSE = mean(err^2)), Likelihood] |>
_[, Likelihood := stringr::str_to_sentence(Likelihood)] |>
knitr::kable(digits = 2)

# Neyman-Scott Illustration
col <- c("#FF4E37FF", "#00B500FF", "#008DFFFF")
ll <- function(v, Y, m) {
  n <- length(m)
  -n * log(v) - 1 / (2 * v) * sum((Y[, 1] - m)^2 + (Y[, 2] - m)^2)
}
llv <- function(s2, Y) {
  n <- nrow(Y)
  v <- apply(Y, 1, diff) / sqrt(2)
  -(n / 2) * log(s2) - (1 / (2 * s2)) * sum(v^2)
}
lik_plots <- function(n, seed=3) {
  set.seed(seed)
  m <- rnorm(n)
  Y <- matrix(rnorm(2*n, rep(m, 2)), n, 2)
  vv <- seq(0.01, 2, len = 299)
  oll <- ll(vv, Y, m)
  pll <- ll(vv, Y, apply(Y, 1, mean))
  mll <- llv(vv, Y)
  plotL(
    vv,
    exp(pll - max(pll)),
    col = col[1],
    xlab = expression(sigma^2),
    ylab = expression("L" * (sigma^2)))
  plotL(vv, exp(mll - max(mll)), col = col[2], add = TRUE)
  plotL(vv, exp(oll - max(oll)), col = col[3], add = TRUE)
  mtext(paste0('n = ', n))
  rightlegend(legend = c('Profile', 'Conditional', 'Oracle'), lwd = 3, col = col)
}
lik_plots(20)