# Test on lambda, Poisson model ----

# H0: lambda = lambda0
set.seed(123)
y <- rpois(20, 1)

n <- length(y)
lambda_hat <- mean(y)      

lambda0 <- c(1, 1.5, 2)

# signed squared-root version of the test statistics
z_lr    <- sign(lambda_hat - lambda0) *
  sqrt(2 * n * (lambda_hat * log(lambda_hat / lambda0) - (lambda_hat - lambda0)))
z_wald  <- (lambda_hat - lambda0) / sqrt(lambda_hat / n)
z_score <- (lambda_hat - lambda0) / sqrt(lambda0 / n)


# 1) H1: lambda != lambda0 -----
tab_chi <- data.frame(
  lambda0 = lambda0,
  LRT     = 1-pchisq(z_lr^2,    df = 1),
  Wald    = pchisq(z_wald^2,  df = 1, lower.tail = FALSE),
  Score   = pchisq(z_score^2, df = 1, lower.tail = FALSE)
)
print(round(tab_chi, 4), row.names = FALSE)

# 2) H1: lambda != lambda0 ----
# from signed squared-root version 
tab_norm <- data.frame(
  lambda0 = lambda0,
  LRT     = 2 * pnorm(-abs(z_lr)),
  Wald    = 2 * pnorm(-abs(z_wald)),
  Score   = 2 * pnorm(-abs(z_score))
)
print(round(tab_norm, 4), row.names = FALSE)


# 3) on-sided alternatives ----
#    H1: lambda > lambda0  -> p = P(Z > z)
#    H1: lambda < lambda0  -> p = P(Z < z)

tab_right <- data.frame(
  lambda0 = lambda0,
  LRT     = pnorm(z_lr,    lower.tail = FALSE),
  Wald    = pnorm(z_wald,  lower.tail = FALSE),
  Score   = pnorm(z_score, lower.tail = FALSE)
)
tab_left <- data.frame(
  lambda0 = lambda0,
  LRT     = pnorm(z_lr),
  Wald    = pnorm(z_wald),
  Score   = pnorm(z_score)
)

print(round(tab_right, 4), row.names = FALSE)

print(round(tab_left, 4), row.names = FALSE)