How to test H0: “this sample is drawn from a gamma distribution” against HA: “this sample is drawn from two different gamma distributions”

Question

This is a targeted follow-up to " all of these data points come from the same distribution." How to test?

I have a sample of strictly positive data. I fit a gamma distribution to it and obtain sample estimates for the shape and rate.

How can I test the hypothesis that all of the data points are drawn from this gamma distribution, against the alternative that some of the data points are drawn from a gamma distribution with different parameters?

Answer 1

Sorry it took me a while to post this answer.

To re-phrase the question again; we are given the sample $X_1,...,X_n$ where $X_i \sim p \Gamma(\alpha_1, \beta_1) + (1-p) \Gamma(\alpha_2, \beta_2)$. The hypothesis test of interest is $$ H_0:\; p=1\;\;\;\;\;H_1:\; p \in (0,1) $$

A possible issue with using a likelihood ratio test in this setting is that the null hypothesis, $p=1$, lies on the boundary of the parameter space $p \in (0,1]$*, causing the asymptotic Chi-squared distribution to no longer hold. This is a well researched topic and there are many papers discussing this issue and proposing new, more robust, hypothesis tests. See here, here,here, and here for a few examples (google "likelihood ratio test boundary of parameter space" and many more results will pop-up). The literature on this topic is very abstruse, but it still may be worth your time if the problem is important to you.

*In fact, the hypotheses are technically not even nested since the true parameter space for a two state mixture model is $p \in (0,1)$ (if $p$ was allowed to be exactly 1 then $\alpha_2$ and $\beta_2$ could change value without effecting the likelihood causing an identification problem) .

The purpose of my answer is to propose 2 alternative testing techniques which, although maybe not exactly what you had in mind, are easier to understand and circumvent the boundary issue.

1. Bayesian Model Comparison

This is not hypothesis testing in the way most people think about it. However, it can be used to accomplish the same objectives and many practitioners actually prefer this technique over traditional hypothesis testing, especially for problems related to the one you're asking about.

Instead of obtaining a p-value which tells you the probability of obtaining data that is "as extreme or more extreme" then that observed conditional on the null hypothesis; Bayesian model comparison uses the posterior probability of the hypothesis/model, $P(H_i|X)$.

The nice thing about Bayesian model comparison is that $P(H_i|X)$ is much easier to interpret than a p-value, it is simply the probability that the hypothesis/model $H_i$ is true given the data. The bad thing about this approach is that estimating the exact value of $P(H_i|X)$ can be extremely computationally burdensome and in some cases virtually impossible.

Fortunately you can approximate $P(H_i|X)$ with the Bayesian Information Criterion (BIC). In the case of just two models, $H_0$ and $H_1$, The posterior probabilities can be approximated as $$ P(H_0|X) \approx \bigg[1+\exp\bigg(\frac{B_0 -B_1}{2}\bigg)\bigg]^{-1}\;\;\;\;\;P(H_1|X)=1-P(H_0|X) $$ where $B_0$ is the BIC of the model estimated under $H_0$ (a single gamma) and $B_1$ is the BIC of the model estimated under $H_1$ (the 2 gamma mixture).

The BIC approximation is an asymptotic one meaning that it gets better as the amount of data increases and may be very unreliable when the number of observations is close to the number of parameters. See This question and answer and the references cited therein for more information on this technique and it's assumptions.

2. Simulate The Null Distribution

Where $\Theta$ is the parameter space, $\Theta_0 \subset \Theta$ is the parameter space under the null hypothesis, and $\hat \theta$ are the parameter estimates; the distribution under the null hypothesis is $f(\hat \theta | \theta \in \Theta_0)$.

The basic idea here is that if we can simulate draws from $f(\hat \theta | \theta \in \Theta_0)$ we can empirically estimate a p-value.

For this to work here, you have to make a more restrictive null hypothesis then the one given above, specifically

$$ H_0:\; p=1,\;\; \alpha_1=\alpha_0,\;\; \beta_1=\beta_0 $$

You would proceed with the simulation as follows

1. draw a sample of size $n$ from $\Gamma(\alpha_0,\beta_0)$
2. fit the sample with the mixture model $p \Gamma(\alpha_1, \beta_1) + (1-p) \Gamma(\alpha_2, \beta_2)$ and collect the point estimate $\hat p$
3. Repeat the above 2 steps until you have collected a large sample of estimates $\hat p_1,...,\hat p_G$

Once you fit the mixture model to your original data $X_1,..,X_n$ and collect the resulting point estimate $\hat p_{*}$, the p-value can be approximated as; $$ Pr(\hat p_{*}|p=1)=\int_0^{\hat p_{*}} f(t|p=1)dt \approx \frac{1}{G}\sum_{g=1}^G I(\hat p_{*}>\hat p_g) $$ The term on the right most side is called the Empirical Distribution Function. $I()$ is an indicator function that is equal to 1 when the statment inside it is true and zero otherwise.

Notice that if you were to replace $\hat p_{*}$ with $\hat p_{g}$ there is a $(1-\alpha)$% chance of falsely rejecting the null hypothesis for $\hat p_{g}$, reflecting a proper type one error rate.

The obvious downsides of this technique are it's computational burden and restriction of the null hypothesis. This method is crude when compared to the techniques proposed in the references above.

Answer 2

We can phrase this problem as follows, following Zachary's suggestion:

• We are given a set of sampled points $X = \{x_1, \dots, x_n\}$.
• We assume that $X \sim p \Gamma(\alpha_1, \beta_1) + (1-p) \Gamma(\alpha_2, \beta_2)$, where $p \in (0, 1]$ is the mixture weight of the first component. Denote the parameters of the mixture by $\theta = (p, \alpha_1, \beta_1, \alpha_2, \beta_2)$, which is in the space of possible parameters $\Theta$.
• Our null hypothesis is $H_0: \theta \in \Theta_0$, where $\Theta_0$ is the set of parameters with $p = 1$. The alternative is $H_1: \theta \not\in \Theta_0$.

This fits the setup for a composite-hypothesis likelihood ratio test. That is, we want to find the likelihood ratio $$ \Lambda(X) = \frac{\sup_{\theta \in \Theta_0} L(\theta; X)}{\sup_{\theta \in \Theta} L(\theta; X)} .$$ Let's call the numerator $\lambda_0$ and the denominator $\lambda$ for convenience. Note that $\lambda \ge \lambda_0$: any value that $L(\theta; X)$ can take on with $\theta \in \Theta_0$ is also a valid result for $\theta \in \Theta$. So the mixture model will necessarily fit the data better; the question is, how much better?

Well, it turns out that as $n \to \infty$, $-2\log\Lambda$ becomes $\chi^2$-distributed under $H_0$, with degrees of freedom $\dim \Theta - \dim \Theta_0 = 5 - 2 = 3$.¹ So we can compute $-2 \log \Lambda$, compare to the (say) 95th percentile of the $\chi^2(3)$ distribution (let's call that threshold $c$), and then reject if $-2 \log \Lambda > c$ (i.e. the mixture model is much more likely than the simple model).

Now, how to compute $\Lambda(X)$?

First, the numerator $\lambda_0$. $\theta \in \Theta_0$ corresponds to a single gamma distribution, so $\sup_{\theta \in \Theta_0} p(X; \theta)$ is the likelihood of the maximum likelihood estimator for the gamma distribution. This is easy enough to find numerically to arbitrary precision.

Next, the denominator $\lambda$. Same thing, but for mixtures of two gammas. This is a harder problem; you should be able to get pretty good results with say EM, or there's some variants specifically for mixtures of gammas², but you can't be sure that you got the best one. Instead, you get some $\hat{\lambda} \le \lambda$, and who knows how big that gap is.

Luckily, though, we don't actually need to find the best $\theta \in \Theta$ to be confident we should reject $H_0$. If we have a $\hat\lambda$ such that $-2 \log \frac{\lambda_0}{\hat \lambda} > c$, then since $\hat\lambda \le \lambda$ we know that $-2 \log \frac{\lambda_0}{\lambda} > c$. Not finding the best mixture reduces the power of our test but doesn't increase the probability of a false rejection.

---

^{1. As I've written it, $\dim \Theta_0$ is actually 4, but clearly the model under $H_0$ should be two-dimensional. We could alternatively have defined $\Theta$ that if $p = 1$, necessarily $\alpha_2 = \beta_2 = 0$, or something along those lines.}

^{2. Schwander and Nielsen. Fast learning of Gamma mixture models with k-MLE. SIMBAD 2013. (doi) (author's site)}