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?

  • $\begingroup$ I think you can estimate a mixture model with two states each constituting two seperate gamma distributions. Where $P$ is the probability of belonging to the first gamma distribution, you can test the null that $P=1$ against the alternative that $P\neq 1$. The parameters of the first distriburion can be fixed inputs from a seperate model or estimated simultaniously in the mixture model. $\endgroup$ Jan 23, 2016 at 2:56

2 Answers 2


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.

  • $\begingroup$ (+1) Note that BIC has many alternatives which sometimes have better properties or seem to perform better, including its predecessor AIC or newer ones like WAIC. Here is a good, though lengthy, overview of some alternatives. $\endgroup$
    – Danica
    Jan 28, 2016 at 19:25
  • $\begingroup$ @Dougal Choice of information criteria for model selection can go either way. I chose BIC here because of its relationship to the marginal likelihood which allows you to compute a posterior probability. This cannot be done with AIC and it’s variants. Andrew Gelman (one of the authors in the paper you reference) argues against the use of posterior probabilities (and BIC) and has written multiple papers on the issue. However, many other statisticians encourage it’s use including Raftery, Hoeting, West, and others. (Similar debate exists on the validity of all frequentist hypothesis testing) $\endgroup$ Jan 29, 2016 at 1:57
  • $\begingroup$ Doesn't the BIC comparison just boil down to a comparison of the likelihood ratio, which is (as you explained) problematic when one of the models is at a boundary of the parameter range. $\endgroup$ Apr 22 at 7:03

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$.1 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 gammas2, 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)

  • 3
    $\begingroup$ I think this is a good answer. The one concern I would have is that because $p=1$ occurs on the boundry of the parameter space, the asymptotic chi - squared no longer holds. I have seen this problem come up in econometrics literature. Whether or not it is a big concern here, I do not know. $\endgroup$ Jan 23, 2016 at 9:32
  • $\begingroup$ @ZacharyBlumenfeld I didn't realize that would be a concern for the LR test. A source to demonstrate the problem would be great, if you have one in mind $\endgroup$ Jan 23, 2016 at 14:58
  • $\begingroup$ This is really helpful, but I'm switching to @ZacharyBlumenfeld's answer because of the points he linked to about the boundary issue in LR tests $\endgroup$ Jan 28, 2016 at 19:22
  • $\begingroup$ Doesn't the log likelihood ratio just follow approximately a mixture of the chi-squared distribution and a pointmass at 0? E.g. if we estimate a normal distribution $X_i \sim N(\mu,1)$ and $\mu \in [0,\infty)$ then $\hat{\mu}_{H_0}=0$ and $\hat{\mu}_{H_a}=max(\bar{x},0)$. Half the time the likelihoods will be equal the other half of the time they follow a chi-squared distribution. So the solution to the problem is to divide the p-value by two. Additionally this could be estimated with a simulation approach. This might be anyways better since the $\chi^2$ distribution is only an approximation. $\endgroup$ Apr 22 at 6:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.