# The t-test as a likelihood ratio test (Normal model, unknown sigma)
# Y_i ~ N(mu, sigma^2),   H0: mu = mu0
# Notation: r_P = signed root of the LRT based on the profile likelihood

# Data and null hypothesis
y   <- rnorm(10, 4, 1)
mu0 <- 5
n   <- length(y)

mu_hat     <- mean(y)                 # MLE of mu
sigma2     <- mean((y - mu_hat)^2)    # MLE of sigma^2 (/n)
sigma2_mu0 <- mean((y - mu0)^2)       # MLE of sigma^2 under H0

# 1) Profile log-likelihood: sigma^2 replaced by its MLE for each mu
lp <- function(s2) -n/2 * log(s2) - n/2

# Difference of profile log-likelihoods:
# lP(mu_hat) - lP(mu0) = n/2 log(sigma2_mu0) - n/2 log(sigma2)
diff_lp <- lp(sigma2) - lp(sigma2_mu0)
diff_lp
n/2 * log(sigma2_mu0) - n/2 * log(sigma2)

# 2) Using sigma2_mu = sigma2 + (mu_hat - mu)^2, we get
#    r_P(mu0) = sign(mu_hat - mu0) * sqrt( n * log(1 + (mu_hat - mu0)^2 / sigma2) )
rP <- sign(mu_hat - mu0) * sqrt(n * log(1 + (mu_hat - mu0)^2 / sigma2))
rP

# 3) Classical t statistic: T = (ybar - mu0) / sqrt(s^2 / n), with s^2 using divisor n-1
t_obs <- (mu_hat - mu0) / sqrt(var(y) / n)
t_obs

# 4) r_P is a monotonic increasing function of T:
#    (mu_hat - mu0)^2 / sigma2 = t^2 / (n-1)
t_grid  <- seq(-4, 4, length.out = 200)
rP_grid <- sign(t_grid) * sqrt(n * log(1 + t_grid^2 / (n - 1)))
plot(t_grid, rP_grid, type = "l", lwd = 2,
     xlab = "T", ylab = expression(r[P]),
     main = "r_P is a monotonic increasing function of T")

# 5) By monotonicity, Pr{R >= r_obs} = Pr{T >= t_obs} under H0:
#    the LRT p-value is obtained from Student's t with n-1 d.f.
#    (alternative H1: mu > mu0)
pt(t_obs, df = n - 1, lower.tail = FALSE)

# 6) normal approximation on r_P (only valid for large n)
pnorm(rP, lower.tail = FALSE)