How to derive the likelihood ratio test for the Poisson distribution?

Question

I am stumped as to how to find the one-sided LRT for the Poisson distribution. The majority of the examples I am finding are for two-sided tests, which does not give all that much help. If you have links to any websites, those would be very useful.

Answer 1

If anything, the one-sided case is simpler than the two-sided case, since the evidence in favour of the alternative hypothesis is only in one direction with respect to the sufficient statistic. Let's work through the problem to have a look a the test in this case.

Since $X_1,...,X_n \sim \text{IID Pois}(\theta)$ you can easily establish that the log-likelihood function is:

$$\ell_\mathbf{x}(\theta) = T(\mathbf{x}) \log(\theta) -n \theta \quad \quad\ \quad \text{for all } \theta>0,$$

where $T(\mathbf{x}) \equiv \sum x_i$ is the sample total, which is a sufficient statistic for the parameter $\theta$. Moreover, using a well-known result for the Poisson distribution, it can be shown that $T(\mathbf{X}) \sim \text{Pois}(n\theta)$.

---

Now, suppose you have hypotheses $H_0: \theta \leqslant \theta_0$ and $H_1: \theta > \theta_0$. Then from the Neyman-Pearson lemma we know that the UMP test is the likelihood-ratio test, which in this case is equivalent to a monotone test on the sample mean. Since the alternative hypothesis is for a larger rate parameter, a higher sample mean constitutes evidence in favour of the alternative hypothesis. Hence, the p-value function for the test is:

$$\begin{equation} \begin{aligned} p(\mathbf{x}) \equiv \sup_{H_0} \mathbb{P}(T(\mathbf{X}) \geqslant T(\mathbf{x})|\theta) &= \sup_{\theta \leqslant \theta_0} \mathbb{P}(T \geqslant T(\mathbf{x})|T \sim \text{Pois}(n\theta)) \\[6pt] &= \mathbb{P}(T \geqslant T(\mathbf{x})|T \sim \text{Pois}(n\theta_0)) \\[12pt] &= 1-\mathbb{P}(T < T(\mathbf{x})|T \sim \text{Pois}(n\theta_0)) \\[6pt] &= 1-\sum_{t = 0}^{T(\mathbf{x})-1} \frac{\theta_0^x e^{-\theta_0}}{x!}. \\[6pt] \end{aligned} \end{equation}$$

The p-value function fully determines the test, and so this constitutes a full specification of the test. As can be seen, the p-value is just the upper-tail of the Poisson distribution with parameter $n \theta_0$.

---

Implementing the test: It is relatively simple to program this test in R using the standard method for programming a hypothesis test. Here we will program a more general version of the test that allows you to stipulate the alternative to be greater or less than the null. (For brevity, we will not include the two-sided test.)

Poisson.rate.test <- function(X, theta0, alternative = "greater") {

#Check validity of inputs
if(!is.numeric(X))      { stop("Error: Data should be numeric"); }
if(length(X) == 0)      { stop("Error: You require at least one observation"); }
if(all(X != as.integer(X)))  { stop("Error: Data should be non-negative integers"); }
if(min(X) < 0)          { stop("Error: Data should be non-negative integers"); }
if(!is.numeric(theta0)) { stop("Error: Parameter theta0 should be numeric"); }
if(length(theta0) != 1) { stop("Error: Parameter theta0 should be a scalar"); }
if(theta0 <= 0)         { stop("Error: Parameter theta0 should be positive"); }
if(!(alternative %in% c("greater", "less"))) 
                        { stop("Error: Alternative must be 'greater' or 'less'"); }

#Set description of test and data
method      <- "Poisson rate test (one-sided)";
data.name   <- paste0(deparse(substitute(X)));

#Set null hypothesis value
null.value  <- theta0;
attr(null.value, "names") <- "rate parameter";

#Calculate test statistics
n           <- length(X);
estimate    <- mean(X);
attr(estimate, "names") <- "estimated rate";
statistic   <- n*estimate;
attr(statistic, "names") <- "sample total";

#Calculate p-value
if (alternative == "less") { 
  p.value   <- ppois(statistic, lambda = n*null.value, 
                     lower.tail = TRUE,  log.p = FALSE); }
if (alternative == "greater") { 
  p.value   <- ppois(statistic, lambda = n*null.value, 
                     lower.tail = FALSE, log.p = FALSE); }
attr(p.value, "names") <- NULL;

#Create htest object
TEST        <- list(method = method, data.name = data.name,
                    null.value = null.value, alternative = alternative,
                    estimate = estimate, statistic = statistic, p.value = p.value);
class(TEST) <- "htest";
TEST; }

We can implement this with some data as an example. I will create a data vector with $n=40$ values that I generated from a Poisson distribution (with a parameter value that I will keep secret). We will test to see if the parameter I used is greater than $\theta_0=15$.

#Create a data vector
DATA <- c(21L, 25L, 13L, 21L,  8L,  9L, 14L, 12L, 14L, 14L, 
          11L, 16L, 19L, 20L, 16L,  9L, 13L, 16L, 17L, 23L, 
          14L, 16L, 16L, 14L, 17L, 14L, 12L,  9L,  8L, 16L, 
          20L, 11L, 10L,  4L, 18L, 15L, 15L, 14L, 16L, 20L);

#Test the alternative hypothesis that theta > 15
TEST <- Poisson.rate.test(DATA, theta0 = 15, alternative = 'greater');
TEST;

        Poisson rate test (one-sided)

data:  DATA
sample total = 590, p-value = 0.6488
alternative hypothesis: true rate parameter is greater than 15
sample estimates:
estimated rate 
         14.75

We can see in this case that there is no evidence that the rate parameter for this data was greater than the hypothesised value of $\theta_0=15$.