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

# binomial
n_b <- 20
x_b <- 6
theta_b <- seq(0, 1, len=499)
L_b <- plotL(theta_b, dbinom(x_b, n_b, theta_b))

# hypergeometric
m <- 30
x_h <- 14
k <- 45
theta_h <- 45:200
L_h <- plotL(theta_h, dhyper(x_h, m, theta_h-m, k))

# point and interval
mle_b <- theta_b[which.max(L_b)]
mle_h <- theta_h[which.max(L_h)]
ci_b <- range(theta_b[L_b > 0.15])
ci_h <- range(theta_h[L_h > 0.15])

# likelihood intervals
par(mfrow=c(1,2))
plotL(theta_b, dbinom(x_b, n_b, theta_b))
abline(h=0.15, lty=2, col='gray', lwd=2)
plotL(theta_h, L_h)
abline(h=0.15, lty=2, col='gray', lwd=2)

# invariance
par(mfrow=c(1,2))
L <- plotL(theta_b, dbinom(x_b, n_b, theta_b))
plotL(theta_b, dbinom(x_b, n_b, theta_b), xaxt='n', xlab=expression(phi1), ylab=expression(L(phi1)))
tt <- seq(0, 1, 0.25)
axis(1, at=tt, label=round(binomial()$linkfun(tt), 2))

# hypothesis test discrepancy
x <- 6
n <- 20
d <- dbinom(0:n, n, 0.5)
p1 <- sum(d[d <= dbinom(x, n, 0.5)])
d <- dnbinom(0:100, x, 0.5)
p2 <- sum(d[d <= dnbinom(n-x, x, 0.5)])