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

# Gamma example: Setup
set.seed(1)
n <- 50
x <- rgamma(n, shape=2, rate = 1)

# MLE
f <- function(a, mx, mlx) log(a) - log(mx) - digamma(a) + mlx
a.mle <- a <- uniroot(f, c(1e-5, 10), mx=mean(x), mlx=mean(log(x)))$root
b.mle <- b <- a/mean(x)

# Profile likelihood illustration
reml <- function(b) {
  uniroot(function(a, b) {log(b) - digamma(a) + mean(log(x))}, c(0.01, 10), b=b)$root
}
l <- function(a, b) {
  n <- length(x)
  n*(a*log(b) - log(gamma(a))) + (a-1)*sum(log(x)) - b*sum(x)
}
f1 <- function(b) {
  a <- reml(b)
  2*(l(a, b) - l(a.mle, b.mle))
}
f2 <- function(b) {
  l(a.mle, b) - l(a.mle, b.mle)
}
bb <- seq(0.75, 3, len=99)

plot(bb, exp(sapply(bb, f1)), lwd=2, type='l', col=pal(2)[2], bty='n', las=1, xlab=expression(beta),
     ylab=expression(L(beta)))
plot(bb, sapply(bb, f1), lwd=2, type='l', col=pal(2)[2], bty='n', las=1, xlab=expression(beta),
     ylab=expression("\u2113"*(beta)))
plot(bb, exp(sapply(bb, f1)), lwd=2, type='l', col=pal(2)[2], bty='n', las=1, xlab=expression(beta),
     ylab=expression(L(beta)))
lines(bb, exp(sapply(bb, f2)), lwd=2, col=pal(2)[1])
plot(bb, sapply(bb, f1), lwd=2, type='l', col=pal(2)[2], bty='n', las=1, xlab=expression(beta),
     ylab=expression("\u2113"*(beta)))
lines(bb, sapply(bb, f2), lwd=2, col=pal(2)[1])
toplegend(legend=c("Plug-in MLE", "Profile"), col=pal(2), lwd=2)
plot(bb, exp(sapply(bb, f1)), lwd=2, type='l', col=pal(2)[2], bty='n', las=1, xlab=expression(beta),
     ylab=expression(L(beta)))
lines(bb, exp(sapply(bb, f2)), lwd=2, col=pal(2)[1])
abline(h=0.15, col='gray', lwd=2, lty=2)
plot(bb, sapply(bb, f1), lwd=2, type='l', col=pal(2)[2], bty='n', las=1, xlab=expression(beta),
     ylab=expression("\u2113"*(beta)))
lines(bb, sapply(bb, f2), lwd=2, col=pal(2)[1])
toplegend(legend=c("Plug-in MLE", "Profile"), col=pal(2), lwd=2)
abline(h=-qchisq(0.95, 2)/2, col='gray', lwd=2, lty=2)

# Neyman-Scott
ll <- function(v, Y, m) {
  n <- length(m)
  -n*log(v) - 1/(2*v)*sum((Y[,1] -m)^2 + (Y[,2] -m)^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))
  plotL(vv, exp(pll-max(pll)), col=pal(2)[1], xlab=expression(sigma^2), ylab=expression("L"*(sigma^2)))
  plotL(vv, exp(oll-max(oll)), col=pal(2)[2], add=TRUE)
  mtext(paste0('n = ', n))
}
lik_plots(10)
lik_plots(20)
lik_plots(40)
lik_plots(1000)