Likelihood Ratio Test for Common Variance from Two Normal Distribution Samples

Question

$\newcommand{\szdp}[1]{\!\left(#1\right)} \newcommand{\szdb}[1]{\!\left[#1\right]}$ Problem Statement: Let $S_1^2$ and $S_2^2$ denote, respectively, the variances of independent random samples of sizes $n$ and $m$ selected from normal distributions with means $\mu_1$ and $\mu_2$ and common variance $\sigma^2.$ If $\mu_1$ and $\mu_2$ are unknown, construct a likelihood ratio test of $H_0: \sigma^2=\sigma_0^2$ against $H_a:\sigma^2=\sigma_a^2,$ assuming that $\sigma_a^2>\sigma_0^2.$

Note 1: This is Problem 10.89 in Mathematical Statistics with Applications, 5th Ed., by Wackerly, Mendenhall, and Sheaffer.

Note 2: This is cross-posted here.

My Work So Far: Let $X_1, X_2,\dots,X_n$ be the sample from the normal distribution with mean $\mu_1,$ and let $Y_1, Y_2,\dots,Y_m$ be the sample from the normal distribution with mean $\mu_2.$ The likelihood is \begin{align*} L(\mu_1,\mu_2,\sigma^2) =\szdp{\frac{1}{\sqrt{2\pi}}}^{\!\!(m+n)} \szdp{\frac{1}{\sigma^2}}^{\!\!(m+n)/2} \exp\szdb{-\frac{1}{2\sigma^2}\szdp{\sum_{i=1}^n(x_i-\mu_1)^2 +\sum_{i=1}^m(y_i-\mu_2)^2}}. \end{align*} We obtain $L\big(\hat{\Omega}_0\big)$ by replacing $\sigma^2$ with $\sigma_0^2$ and $\mu_1$ with $\overline{x}$ and $\mu_2$ with $\overline{y}:$ \begin{align*} L\big(\hat{\Omega}_0\big) =\szdp{\frac{1}{\sqrt{2\pi}}}^{\!\!(m+n)} \szdp{\frac{1}{\sigma_0^2}}^{\!\!(m+n)/2} \exp\szdb{-\frac{1}{2\sigma_0^2}\szdp{\sum_{i=1}^n(x_i-\overline{x})^2 +\sum_{i=1}^m(y_i-\overline{y})^2}}. \end{align*} The MLE for the common variance in exactly this scenario is: $$\hat\sigma^2=\frac{1}{m+n}\szdb{\sum_{i=1}^n(x_i-\overline{x})^2 +\sum_{i=1}^m(y_i-\overline{y})^2}.$$ So this estimator plugged into the likelihood yields \begin{align*} L\big(\hat{\Omega}\big) &=\szdp{\frac{1}{\sqrt{2\pi}}}^{\!\!(m+n)} \szdp{\frac{1}{\hat\sigma^2}}^{\!\!(m+n)/2} \exp\szdb{-\frac{1}{2\hat\sigma^2}\szdp{\sum_{i=1}^n(x_i-\overline{x})^2 +\sum_{i=1}^m(y_i-\overline{y})^2}}. \end{align*} It follows that the ratio is \begin{align*} \lambda &=\frac{L\big(\hat{\Omega}_0\big)}{L\big(\hat{\Omega}\big)}\\ &=\szdp{\frac{\hat\sigma^2}{\sigma_0^2}}^{\!\!(m+n)/2} \exp\szdb{\frac{(\sigma_0^2-\hat\sigma^2)(m+n)}{2\sigma_0^2}}.\\ -2\ln(\lambda) &=(m+n)\szdb{\frac{\hat\sigma^2}{\sigma_0^2} -\ln\szdp{\frac{\hat\sigma^2}{\sigma_0^2}}-1}. \end{align*} Now the function $f(x)=x-\ln(x)-1$ first decreases, then increases. It has a global minimum of $0$ at $x=1.$ Note also that the original inequality becomes: \begin{align*} \lambda&<k\\ 2\ln(\lambda)&<2\ln(k)\\ -2\ln(\lambda)&>k'. \end{align*} As the test is for $\sigma_a^2>\sigma_0^2,$ we will expect the estimator $\hat\sigma^2>\sigma_0^2.$ We can, evidently, use Theorem 10.2 and claim that $-2\ln(\lambda)$ is $\chi^2$ distributed with d.o.f. $1-0.$ So we reject $H_0$ when $$(m+n)\szdb{\frac{\sum_{i=1}^n(x_i-\overline{x})^2 +\sum_{i=1}^m(y_i-\overline{y})^2}{(m+n)\sigma_0^2} -\ln\szdp{\frac{\sum_{i=1}^n(x_i-\overline{x})^2 +\sum_{i=1}^m(y_i-\overline{y})^2}{(m+n)\sigma_0^2}}-1} >\chi^2_{\alpha},$$ or $$\frac{\sum_{i=1}^n(x_i-\overline{x})^2 +\sum_{i=1}^m(y_i-\overline{y})^2}{\sigma_0^2} -\ln\szdp{\frac{\sum_{i=1}^n(x_i-\overline{x})^2 +\sum_{i=1}^m(y_i-\overline{y})^2}{\sigma_0^2}}-(m+n) >\chi^2_{\alpha}.$$

My Questions:

1. Is my answer correct?
2. My answer is not the book's answer. The book's answer is simply that $$\chi^2=\frac{(n-1)S_1^2+(m-1)S_2^2}{\sigma_0^2}$$ has a $\chi_{(n+m-2)}^2$ distribution under $H_0,$ and that we reject $H_0$ if $\chi^2>\chi_a^2.$ How is this a likelihood ratio test? It's not evident that they went through any of the steps of forming the likelihood ratio with all the necessary optimizations. Their estimator is not the MLE for $\sigma^2,$ is it?

Answer 1

Let $x_1,...,x_n\sim N(\mu_1,\sigma^2)$ and $y_1,...,y_m\sim N(\mu_2,\sigma^2)$. Under $H_0$, the likelihood function is

$$L_0=\left(\frac{1}{\sqrt{2\pi\sigma_0^2}}\right)^nexp\left( \frac{-1}{2\sigma_0^2}\left( \sum_{i=1}^{n}(x_i-\bar{x})^2 + \sum_{i=1}^{m}(y_i-\bar{y})^2 \right) \right)=\left(\frac{1}{\sqrt{2\pi\sigma_0^2}}\right)^nexp\left( \frac{-1}{2\sigma_0^2}\left( (n-1)S_1^2 + (m-1)S_2^2 \right) \right).$$

Under $H_a$ we get:

$$L_a=\left(\frac{1}{\sqrt{2\pi\sigma_a^2}}\right)^nexp\left( \frac{-1}{2\sigma_a^2}\left( (n-1)S_1^2 + (m-1)S_2^2 \right) \right).$$

Inside the latter exponent, we have the same sum of $S_1^2$ and $S_2^2$ but multiplied by $\frac{1}{2\sigma_a^2}$, which can be presented as $\frac{\sigma_0^2}{\sigma_a^2}\frac{1}{2\sigma_0^2}$. The LRT is then

$$\lambda=\frac{L_0}{L_a}=\left( \frac{\sqrt{\sigma_a}}{\sqrt{\sigma_0}} \right)^nexp\left( \frac{-1}{2\sigma_0^2}( (n-1)S_1^2 + (m-1)S_2^2) + \frac{\sigma_0^2}{\sigma_a^2}\frac{1}{2\sigma_0^2}( (n-1)S_1^2 + (m-1)S_2^2) \right)=\left( \frac{\sqrt{\sigma_a}}{\sqrt{\sigma_0}} \right)^nexp\left(-\frac{1}{2} \left( 1- \frac{\sigma_0^2}{\sigma_a^2}\right)\frac{1}{\sigma_0^2}( (n-1)S_1^2 + (m-1)S_2^2) \right).$$

Given that $\sigma_a>\sigma_0$, the expression $\left( \frac{\sqrt{\sigma_a}}{\sqrt{\sigma_0}} \right)^n$ is larger that 1 and $\left( 1- \frac{\sigma_0^2}{\sigma_a^2}\right)$ is positive, for all $n$. Thus, our test depends on $\frac{1}{2\sigma_0^2}( (n-1)S_1^2 + (m-1)S_2^2)$ (note I ignored the minus 0.5 there, as it is a constant), which is a $\chi^2$ variable with $n+m-2$ degrees of freedom (why? simple to find).

Ultimately, our test is in the form of

$$\left\{ \frac{1}{\sigma_0^2}\left( (n-1)S_1^2 + (m-1)S_2^2\right) > \chi^2_{\alpha,~n+m-2} \right\}$$

Answer 2

Ok, part of this is that you are using a theorem about asymptotic distribution of log likelihood ratios and they are using the exact distribution of the residual sum of squares. I agree that the book answer is a little underexplained.

Personally, I would have taken the the logarithm first, then subtracted, but I think your calculations are right. You find that the log likelihood is a function of $\hat \sigma^2/\sigma^2_0$, and use the asymptotic results to get a $\chi^2_1$ reference distribution. I think you're wrong in that the asymptotic theory is without the restriction $\sigma^2_a>\sigma^2_0$, so that with the restriction you should have $\frac{1}{2}\chi^2_1$ as the asymptotic null sampling distribution. The probability is halved because you only reject when $\hat\sigma^2>\sigma^2_0$, even though you can just as easily get a large likelihood ratio with $\hat\sigma^2<\sigma^2_0$.

Once you find that the loglikelihood ratio depends on the data only through $\hat \sigma^2/\sigma^2_0$, you can look for other functions of this statistic that have more tractable null sampling distributions. Any rejection region you define based on the log likelihod ratio plus the restriction $\hat\sigma^2>\sigma^2_0$ can also be defined based on $\hat \sigma^2/\sigma^2_0$; it's the same set of rejection regions, so it's the same test. In more complicated cases you aren't going to be able to find an exact sampling distribution so there's no point and you should just stick with the loglikelihod ratio as the test statistic. This case is simple enough that there might be a tractable exact distribution -- especially as it's a textbook example.

Since the data are Normal with variance $\sigma^2_0$ under the null, the residual sum of squares after estimating two means has a null $\sigma^2_0\chi^2_{n+m-2}$ distribution. Thus $\mathrm{RSS}/\sigma^2_0$ has a $\chi^2_{n+m-2}$ null distribution. You will reject $H_0$ if the ratio is large, so the test is as described in the book.

Their test will have closer to the nominal size in small (and moderate) sample sizes, eg

m<-6
n<-10
sigma2_0<-17
mu1<-420
mu2<-69


sim_once<-function(){
y<-rnorm(m+n, m=rep(c(mu1,mu2), c(m,n)),s=sqrt(sigma2_0))

xbar1<-mean(y[1:m])
xbar2<-mean(y[-(1:m)])
xbars<-rep(c(xbar1,xbar2),c(m,n))

ell_0 <- dnorm(y,xbars, s=sqrt(sigma2_0),log=TRUE)

RSS<- sum((y-xbars)^2)
hatsigma2<-max(sigma2_0, RSS/(m+n))

ell_1 <- dnorm(y,xbars, s=sqrt(hatsigma2), log=TRUE)

lrt<-sum(ell_1-ell_0)
c(lrt, RSS/sigma2_0)
}

results<-replicate(10000,sim_once())

gives

> mean(results[1,]*2>qchisq(0.95,1))
[1] 0.0088

for the asymptotic distribution of $2\ell$ that you use

> mean(results[1,]*2>  0.5*qchisq(0.95,1)+0.5*0)
[1] 0.0334

for the one-tailed version of the asymptotic distribution of $2\ell$ and

> mean(results[2,]>qchisq(0.95,m+n-2))
[1] 0.0501

for the exact LRT from the book.