Finding region of rejection with likelihood ratio test Let $X_1,\ldots,X_n$ be i.i.d. from a Gamma distribution with p.d.f. $f(x;\theta) = \theta^{-2} x e^{-x/\theta}$ for $x>0$ where $\theta$ is an unknown parameter.
I would like to test the hypothesis $H_0 : \theta = \theta_0$ against $H_a = \theta \neq \theta_0$ where $\theta_0$ is known.
In particular, I want to use the LRT test to determine the region of rejection at a level of $\alpha$.
I know I need to compute $L_0/L_1$ where $L_0$ is the max likelihood when parameters are restricted and $L_1$ is the max value when they aren't resreicted, but I'm really struggling to do so. I tried doing log-likelihood and so on. This is an exercise that I can't solve; I've seen related examples but can't figure this out for the specified Gamma distribution.
I would like to learn how to solve this by hand (without any software like R)
My attempt:
I got the likelihood function is $L(\theta) = \theta^{-2n} e^{-n\overline{x}/\theta} \cdot \prod_i x_i$ where $\overline{x} = \sum_i x_i$. Now taking a log gives $-2n\log(\theta) -n\overline{x}/\theta + n\overline{x}$. Differentiating and equating to $0$ gives $-2n/\theta + n\overline{x}/\theta^2 = 0 $ which means $\theta = \overline{x}/2$ where $\overline{x} = n^{-1}\sum_i x_i$.  Now I'm not sure how exactly to compute $L_0/L_1$. I get
$$\frac{L_0}{L_1} = \frac{\theta_0^{-2n}e^{-n\overline{x}/\theta_0}}{(\overline{x} /2)^{-2n} e^{-2n}} $$
Is this right? If so, how do I get my rule for rejecting $H_0$ in terms of this expression? I think it has something to do with a chi-squared distribution but not entirely sure.
 A: Likelihood function given the sample $\boldsymbol x=(x_1,\ldots,x_n)$ is
$$L(\theta\mid \boldsymbol x)=\theta^{-2n}\left(\prod_{i=1}^n x_i \right)\exp\left\{-\frac1{\theta}\sum_{i=1}^n x_i\right\}\mathbf1_{(0,\infty)^n}(x_1,\ldots,x_n)\quad,\,\theta>0$$
This gives the ML estimate $$\hat\theta=\frac1{2n}\sum_{i=1}^n x_i$$
As you have found, the LR test criterion takes the form
$$\Lambda(\boldsymbol x)=\frac{L(\theta_0\mid \boldsymbol x)}{L(\hat\theta\mid \boldsymbol x)}=\left(\frac1{2n\theta_0}\sum_{i=1}^n x_i\right)^{2n}\exp\left\{-\frac1{\theta_0}\sum_{i=1}^nx_i+2n\right\}=g(T)\,,$$
where $T(\boldsymbol x)=\sum\limits_{i=1}^n x_i$ and $g(x)=cx^{2n}e^{-x/\theta_0}$ for $x>0$ and some positive constant $c$.
If you roughly sketch the function $g$ (by checking the signs of $g'$ or $g''$ for example), you will find that it looks like a bell-shaped curve.
We reject $H_0$ for small values of $\Lambda$. So if $\Lambda<k$, i.e. if $g(T)<k$, keeping in mind the nature of the function $g$, the critical region has the form "$T<k_1$ or $T>k_2$" where $k_1<k_2$ are such that $$P_{\theta_0}(T(\boldsymbol X)<k_1)+P_{\theta_0}(T(\boldsymbol X)>k_2)=\alpha$$
and $$g(k_1)=g(k_2)$$
Now if $X$ has the pdf $f_X$ (say) in the question, then $Y=\frac{2}{\theta}X$ has pdf
$$f_Y(y)=f_X\left(\frac{\theta y}{2}\right)\left|\frac{\mathrm dx}{\mathrm dy}\right|=\frac14 ye^{-y/2}\mathbf1_{(0,\infty)}(y)$$
This is the density of a $\chi^2_4$ distribution, whence $$\frac2{\theta_0}\sum_{i=1}^n X_i=\frac2{\theta_0}T(\boldsymbol X)\stackrel{H_0}\sim \chi^2_{4n}$$
You can take $P_{\theta_0}(T(\boldsymbol X)<k_1)=P_{\theta_0}(T(\boldsymbol X)>k_2)=\frac{\alpha}2$ for convenience and write $k_1,k_2$ in terms of $\chi^2_{4n}$ quantiles but that is not necessarily the only possible solution.
A: This is the same as asking for a $(1-\alpha)\times 100\%$ confidence interval for the Likelihood Ratio Test.  Anything outside this set would lead to a rejection of the null.
These regions are often computed numerically.  Let $G^2 = F(\theta, \theta_0)$ be the LRT statistic.  You would need to solve for
$$ F(\theta, \theta_0) - \chi^2_{df, \alpha} = 0$$
Here, $\chi^2_{df,a}$ is the critical value of the test. Finding the roots of the equation above will give a confidence interval for $\theta$ and the compliment of this set is the rejection region.
Here is an example in R:



true_theta = 2.0
null_theta = 1.0
n_samples = 10

set.seed(0)
data = rgamma(n_samples, shape=1, rate = true_theta)

lrt_statistic = function(theta_hat) {
  log_lik_estimated = dgamma(data, shape=1, rate=theta_hat, log = T)
  log_lik_null = dgamma(data, shape=1, rate=null_theta, log = T)
  
  G2 = 2*sum(log_lik_estimated) - 2*sum(log_lik_null)
  
  G2
}

crit_stat = qchisq(0.05, n_samples-1)

obj_func = Vectorize(function(x) lrt_statistic(x) - crit_stat)


tt = seq(0, 6, 0.02)
plot(tt, obj_func(tt), type='l')
abline(h=0)

theta_est = 1/mean(data)
left_root = uniroot(obj_func, interval = c(0, theta_est ))$root
right_root = uniroot(obj_func, interval = c(theta_est, 6 ))$root

root = c(left_root, right_root)
[1] 1.319097 5.361789

```

