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

# Quadratic approx
n <- 20
x <- 6
theta <- seq(0, 1, len=499)
plotL(theta, dbinom(x, n, theta))
mtext('Binomial (n=20, x=6)')
plotl(theta, dbinom(x, n, theta, log=TRUE))
mtext('Binomial (n=20, x=6)')
theta <- seq(0.10, 0.55, len=499)
plotl(theta, dbinom(x, n, theta, log=TRUE))
mtext('Binomial (n=20, x=6)')
plotl(theta, dbinom(x, n, theta, log=TRUE))
mtext('Binomial (n=20, x=6)')
Info <- n/(6/20*14/20)
lines(theta, -Info*(theta-6/20)^2/2, lty=2)
m <- 30
x <- 14
k <- 45
theta <- 45:200
plotL(theta, dhyper(x, m, theta-m, k))
mtext('Hypergeometric (m=30, x=14, k=45)')
plotl(theta, dhyper(x, m, theta-m, k, log=TRUE))
mtext('Hypergeometric (m=30, x=14, k=45)')
theta <- 70:160
plotl(theta, dhyper(x, m, theta-m, k, log=TRUE))
mtext('Hypergeometric (m=30, x=14, k=45)')
l <- function(theta) dhyper(x, m, theta-m, k, log=TRUE)
plotl(theta, l(theta))
mtext('Hypergeometric (m=30, x=14, k=45)')
mle <- theta[which.max(l(theta))]
Hess <- ((l(mle+1) - l(mle)) - (l(mle) - l(mle-1)))
lines(theta, Hess*(theta-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)))

# 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
hessian(l, mle)

# Binomial n=10/100
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)')
Info <- n/(8/10*2/10)
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)
mean(Cov)
n <- 100
x <- 80
theta <- seq(0.7, 0.88, len=499)
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)
mean(Cov)
n <- 1000
x <- 800
theta <- seq(0.77, 0.83, len=499)
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)
mean(Cov)

# Normal means
plot(1, -3, bty='n', las=1, pch=16, xlim=c(-3,5), ylim=c(-7,1), xlab=expression(theta[1]), ylab=expression(theta[2]))
theta <- seq(0, 2*pi, len=101)
r <- sqrt(qchisq(0.95, 2))
lines(1+r*cos(theta), -3 + r*sin(theta), lty=2, lwd=3, col='gray')

# Persp
x <- seq(-3, 5, len=99)
y <- seq(-7, 1, len=99)
z <- outer(x, y, function(x,y) {-(x-1)^2 - (y+3)^2})
rgl::persp3d(x, y, z, col='slateblue')
rgl::planes3d(0, 0, 1, qchisq(0.95, 2), col='gray')

# AIC vs chisq
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))
toplegend(legend=c(expression(chi^2), '2p'), col=pal(2), lwd=3)