More on Likelihood Ratio Test, the following problem is originally from Casella and Berger (2001), exercise 8.12.
Problem
For samples of size n=1,4,16,64,100 from a normal population with mean \mu and known variance \sigma^2, plot the power function of the following LRTs (Likelihood Ratio Tests). Take \alpha = .05.- H_0:\mu\leq 0 versus H_1:\mu>0
- H_0:\mu=0 versus H_1:\mu\neq 0
Solution
-
The LRT statistic is given by
\lambda(\mathbf{x})=\frac{\displaystyle\sup_{\mu\leq 0}\mathcal{L}(\mu|\mathbf{x})}{\displaystyle\sup_{-\infty<\mu<\infty}\mathcal{L}(\mu|\mathbf{x})}, \;\text{since }\sigma^2\text{ is known}.
The denominator can be expanded as follows: \begin{aligned} \sup_{-\infty<\mu<\infty}\mathcal{L}(\mu|\mathbf{x})&=\sup_{-\infty<\mu<\infty}\prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{(x_i-\mu)^2}{2\sigma^2}\right]\\ &=\sup_{-\infty<\mu<\infty}\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\displaystyle\sum_{i=1}^{n}\frac{(x_i-\mu)^2}{2\sigma^2}\right]\\ &=\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\displaystyle\sum_{i=1}^{n}\frac{(x_i-\bar{x})^2}{2\sigma^2}\right],\\ &\quad\text{since }\bar{x}\text{ is the MLE of }\mu.\\ &=\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\frac{n-1}{n-1}\displaystyle\sum_{i=1}^{n}\frac{(x_i-\bar{x})^2}{2\sigma^2}\right]\\ &=\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\frac{(n-1)s^2}{2\sigma^2}\right],\\ \end{aligned}while the numerator is evaluated as follows: \begin{aligned} \sup_{\mu\leq 0}\mathcal{L}(\mu|\mathbf{x})&=\sup_{\mu\leq 0}\prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{(x_i-\mu)^2}{2\sigma^2}\right]\\ &=\sup_{\mu\leq 0}\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\displaystyle\sum_{i=1}^{n}\frac{(x_i-\mu)^2}{2\sigma^2}\right]. \end{aligned}Above expression will attain its maximum if the value inside the exponential function is small. And for negative values of \mu\in(-\infty,0) the quantity (x_i-\mu)^2 would be large, implies that the exponential term would become small. Therefore, the only value that will give us the supremum likelihood is \mu=\mu_0=0. Hence, \begin{aligned} \sup_{\mu\leq 0}\mathcal{L}(\mu|\mathbf{x})&=\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\displaystyle\sum_{i=1}^{n}\frac{(x_i-\mu_0)^2}{2\sigma^2}\right]\\ =\frac{1}{(2\pi\sigma^2)^{1/n}}&\exp\left[-\displaystyle\sum_{i=1}^{n}\frac{(x_i-\bar{x}+\bar{x}-\mu_0)^2}{2\sigma^2}\right]\\ =\frac{1}{(2\pi\sigma^2)^{1/n}}&\exp\left\{-\displaystyle\sum_{i=1}^{n}\left[\frac{(x_i-\bar{x})^2+2(x_i-\bar{x})(\bar{x}-\mu_0)+(\bar{x}-\mu_0)^2}{2\sigma^2}\right]\right\}\\ =\frac{1}{(2\pi\sigma^2)^{1/n}}&\exp\left[-\frac{(n-1)s^2+n(\bar{x}-\mu_0)^2}{2\sigma^2}\right], \\ &\text{since the middle term is 0.}\\ =\frac{1}{(2\pi\sigma^2)^{1/n}}&\exp\left[-\frac{(n-1)s^2+n\bar{x}^2}{2\sigma^2}\right], \text{since }\mu_0=0.\\ \end{aligned}So that \begin{equation} \label{eq:lrtre} \begin{aligned} \lambda(\mathbf{x})&=\frac{\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\frac{(n-1)s^2+n\bar{x}^2}{2\sigma^2}\right]}{\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\frac{(n-1)s^2}{2\sigma^2}\right]}\\ &=\exp\left[-\frac{n\bar{x}^2}{2\sigma^2}\right].\\ \end{aligned} \end{equation}And we reject the null hypothesis if \lambda(\mathbf{x})\leq c, that is \begin{aligned} \exp\left[-\frac{n\bar{x}^2}{2\sigma^2}\right]&\leq c\\ -\frac{n\bar{x}^2}{2\sigma^2}&\leq \log c\\ \frac{\lvert\bar{x}\rvert}{\sigma/\sqrt{n}}&\geq\sqrt{-2\log c}=c'. \end{aligned}
Figure 1: Plot of Likelihood Ratio Test Statistic for n = 4,\sigma = 1.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characterslibrary(magrittr) library(lattice) x3 <- seq(-2, .5, length.out = 100) lrt_fun1 <- function (x, n = 4, s = 1) { exp(- (n * x ** 2) / (2 * s ** 2)) } { lrt_fun1(x4) ~ x4 } %>% xyplot( xlab = expression(bar(x)), ylab = expression( paste( 'exp', bgroup("[", frac(- n * bar(x) ** 2, 2 * sigma ** 2),"]") ) ), xlim = c(-2.2, 2.2), ylim = c(-0.05, 1.05), col = 'black', lwd = 2, type = 'l', panel = function (x, y, ...) { panel.grid(h = -1, v = -1) panel.xyplot(x, y, ...) from.z <- .5 to.z <- 2 S.x <- c(from.z, seq(from.z, to.z, 0.01), to.z) S.y <- c(0, exp(- (n * (seq(from.z, to.z, 0.01)) ** 2) / (2 * s ** 2)), 0) panel.polygon(S.x, S.y, col = "gray", border = "white") panel.segments(.5, 0, 2, 0) panel.segments(.5, 0, .5, .61) panel.abline(v = 0, lty = 2) panel.arrows(0.5, exp(- (n * .5 ** 2) / (2 * s ** 2)), -1.25, exp(- (n * .5 ** 2) / (2 * s ** 2)), angle = 20) panel.text(-1.5, exp(- (n * .5 ** 2) / (2 * s ** 2)), labels = "c = .61") panel.arrows(.7, .2, 1.2, .4, angle = 20) panel.text(1.5, .46, labels = expression(paste(frac(bar(x), sigma/sqrt(n)) >= c, "'", sep = ""))) panel.arrows(1, exp(- (n * (-1) ** 2) / (2 * s ** 2)) + .02, 1.4, .19, angle = 20) panel.text(1.7, .22, labels = expression(lambda(bold(x)) <= c)) panel.curve(exp(- (n * x ** 2) / (2 * s ** 2)), .5, 2, lwd = 2, lty = 2, add = TRUE) from.z <- -2 to.z <- -.5 S.x <- c(from.z, seq(from.z, to.z, 0.01), to.z) S.y <- c(0, exp(- (n * (seq(from.z, to.z, 0.01)) ** 2) / (2 * s ** 2)), 0) panel.polygon(S.x, S.y, col = "gray", border = "white") panel.segments(-.5, 0, -2, 0) panel.segments(-.5, 0, -.5, .61) panel.arrows(-.7, .2, -1.2, .4, angle = 20) panel.text(-1.25, .46, labels = expression(paste(frac(bar(x), sigma/sqrt(n)) <= -c, "'", sep = ""))) panel.arrows(-1, exp(- (n * (-1) ** 2) / (2 * s ** 2)) + .02, -1.4, .19, angle = 20) panel.text(-1.7, .22, labels = expression(lambda(bold(x)) <= c)) panel.curve(exp(- (n * x ** 2) / (2 * s ** 2)), -2, -.5, lwd = 2, lty = 2, add = TRUE) } )
So that the power function, that is the probability of rejecting the null hypothesis given that it is true (the probability of Type I error) is, \begin{aligned} \beta(\mu)&=\mathrm{P}\left[\frac{\bar{x}-\mu_0}{\sigma/\sqrt{n}}\geq c'\right],\quad\mu_0=0\\ &=1-\mathrm{P}\left[\frac{\bar{x}+\mu-\mu-\mu_0}{\sigma/\sqrt{n}}< c'\right]\\ &=1-\mathrm{P}\left[\frac{\bar{x}-\mu}{\sigma/\sqrt{n}} + \frac{\mu-\mu_0}{\sigma/\sqrt{n}}< c'\right]\\ &=1-\mathrm{P}\left[\frac{\bar{x}-\mu}{\sigma/\sqrt{n}}< c'+ \frac{\mu_0-\mu}{\sigma/\sqrt{n}}\right]\\ &=1-\Phi\left[c'+ \frac{\mu_0-\mu}{\sigma/\sqrt{n}}\right]. \end{aligned}Values taken by \Phi are negative and so it decreases, but since we subtracted it to 1, then \beta(\mu) is an increasing function. So that for \alpha=.05, \begin{aligned} \alpha&=\sup_{\mu\leq \mu_0}\beta(\mu)\\ .05&=\beta(\mu_0)\Rightarrow\beta(\mu_0)=1-\Phi(c')\\ .95&=\Phi(c')\Rightarrow c'=1.645. \end{aligned}Since, \begin{aligned} \Phi(1.645)=\int_{-\infty}^{1.645}\frac{1}{\sqrt{2\pi}}\exp\left[-\frac{x^2}{2}\right]\operatorname{d}x=.9500151. \end{aligned}Therefore for c'=1.645,\mu_0=0,\sigma=1, the plot of the power function as a function of \mu for different sample size, n, is shown in Figure 2. For example, for n=1 we compute for the function \begin{equation} \label{eq:powcomp} \begin{aligned} \beta(\mu)&=1-\Phi\left[c'+ \frac{\mu_0-\mu}{\sigma/\sqrt{n}}\right]\\ &=1-\Phi\left[1.645+ \frac{0-\mu}{1/\sqrt{1}}\right]\\ &=1-\int_{-\infty}^{\left(1.645+ \frac{0-\mu}{1/\sqrt{1}}\right)}\frac{1}{\sqrt{2\pi}}\exp\left[-\frac{x^2}{2}\right]\operatorname{d}x. \end{aligned} \end{equation}The obtained values would be the y. For n = 64 \begin{aligned} \beta(\mu)&=1-\Phi\left[c'+ \frac{\mu_0-\mu}{\sigma/\sqrt{n}}\right]\\ &=1-\Phi\left[1.645+ \frac{0-\mu}{1/\sqrt{64}}\right]\\ &=1-\int_{-\infty}^{\left(1.645+ \frac{0-\mu}{1/\sqrt{64}}\right)}\frac{1}{\sqrt{2\pi}}\exp\left[-\frac{x^2}{2}\right]\operatorname{d}x, \end{aligned}and so on.Figure 2: Power Function for Different Values of n.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characterslibrary(magrittr) library(lattice) phi_fun1 <- function (x, mu0 = 0, s = 1, cr = 1.645, n = c(1, 4, 16, 64, 100), i = 1) { 1 - pnorm(cr + (mu0 - x) / (s / sqrt(n[i]))) } x1 <- seq(-2, 5, length.out = 100) { phi_fun1(x1) ~ x1 } %>% xyplot( xlab = expression(mu), ylab = expression( paste( "1 - ", Phi, bgroup("[", c + frac(mu[0]-mu, sigma/sqrt(n)), "]"), sep = "")), col = 'black', lwd = 2, type = 'l', panel = function (x, y, ...) { panel.grid(h = -1, v = -1) for (i in 2:5) { panel.curve(phi_fun1(x, i = i), -2, 5, add = TRUE, lty = i, lwd = 2) } panel.xyplot(x, y, ...) }, key = list(corner = c(.9, .1), text = list(c('n = 1', 'n = 4', 'n = 16', 'n = 64', 'n = 100')), lines = list(lty = c(1, 2, 3, 4, 5), lwd = 2) ) ) -
The LRT statistic is given by
\lambda(\mathbf{x})=\frac{\displaystyle\sup_{\mu= 0}\mathcal{L}(\mu|\mathbf{x})}{\displaystyle\sup_{-\infty<\mu<\infty}\mathcal{L}(\mu|\mathbf{x})}, \;\text{since }\sigma^2\text{ is known}.
The denominator can be expanded as follows: \begin{aligned} \sup_{-\infty<\mu<\infty}\mathcal{L}(\mu|\mathbf{x})&=\sup_{-\infty<\mu<\infty}\prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{(x_i-\mu)^2}{2\sigma^2}\right]\\ &=\sup_{-\infty<\mu<\infty}\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\displaystyle\sum_{i=1}^{n}\frac{(x_i-\mu)^2}{2\sigma^2}\right]\\ &=\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\displaystyle\sum_{i=1}^{n}\frac{(x_i-\bar{x})^2}{2\sigma^2}\right],\\ &\quad\;\text{since }\bar{x}\text{ is the MLE of }\mu.\\ &=\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\frac{n-1}{n-1}\displaystyle\sum_{i=1}^{n}\frac{(x_i-\bar{x})^2}{2\sigma^2}\right]\\ &=\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\frac{(n-1)s^2}{2\sigma^2}\right],\\ \end{aligned}and the numerator is evaluated as follows: \begin{aligned} \sup_{\mu=0}\mathcal{L}(\mu|\mathbf{x})&=\sup_{\mu=0}\prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{(x_i-\mu)^2}{2\sigma^2}\right]\\ &=\sup_{\mu=0}\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\displaystyle\sum_{i=1}^{n}\frac{(x_i-\mu)^2}{2\sigma^2}\right]\\ &=\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\displaystyle\sum_{i=1}^{n}\frac{(x_i-0)^2}{2\sigma^2}\right]\\ &=\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\frac{(n-1)s^2+n\bar{x}^2}{2\sigma^2}\right], \end{aligned}we skip some lines in the above simplification since we've done this already in part (a). And by Equation (1), \lambda(\mathbf{x})=\exp\left[-\frac{n\bar{x}^2}{2\sigma^2}\right]. So that \lambda(\mathbf{x})\leq c would be \begin{aligned} \exp\left[-\frac{n\bar{x}^2}{2\sigma^2}\right]&\leq c\\ -\frac{n\bar{x}^2}{2\sigma^2}&\leq \log c\\ \frac{\lvert\bar{x}-\mu_0\rvert}{\sigma/\sqrt{n}}&\geq\sqrt{-2\log c}=c',\quad \mu_0=0. \end{aligned}So rejecting the null hypothesis if \lambda(\mathbf{x})\leq c' is equivalent to rejecting H_0 if \frac{\lvert\bar{x}\rvert}{\sigma/\sqrt{n}}\geq c'. And since H_1 is two-sided, then we reject H_0 if \frac{\bar{x}}{\sigma/\sqrt{n}}\geq c' or \frac{\bar{x}}{\sigma/\sqrt{n}}\leq -c'. To illustrate this, consider Figure 3 where the two shaded regions are the lower and upper rejection regions.
Figure 3: Plot of Likelihood Ratio Test Statistic for n = 4,\sigma = 1.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characterslibrary(magrittr) library(lattice) x4 <- seq(-.5, .5, length.out = 100) { lrt_fun1(x4) ~ x4 } %>% xyplot( xlab = expression(bar(x)), ylab = expression( paste( 'exp', bgroup("[", frac(- n * bar(x) ** 2, 2 * sigma ** 2),"]") ) ), xlim = c(-2.2, 2.2), ylim = c(-0.05, 1.05), col = 'black', lwd = 2, type = 'l', panel = function (x, y, ...) { panel.grid(h = -1, v = -1) panel.xyplot(x, y, ...) from.z <- .5 to.z <- 2 S.x <- c(from.z, seq(from.z, to.z, 0.01), to.z) S.y <- c(0, exp(- (n * (seq(from.z, to.z, 0.01)) ** 2) / (2 * s ** 2)), 0) panel.polygon(S.x, S.y, col = "gray", border = "white") panel.segments(.5, 0, 2, 0) panel.segments(.5, 0, .5, .61) panel.abline(v = 0, lty = 2) panel.arrows(0.5, exp(- (n * .5 ** 2) / (2 * s ** 2)), -1.25, exp(- (n * .5 ** 2) / (2 * s ** 2)), angle = 20) panel.text(-1.5, exp(- (n * .5 ** 2) / (2 * s ** 2)), labels = "c = .61") panel.arrows(.7, .2, 1.2, .4, angle = 20) panel.text(1.5, .46, labels = expression(paste(frac(bar(x), sigma/sqrt(n)) >= c, "'", sep = ""))) panel.arrows(1, exp(- (n * (-1) ** 2) / (2 * s ** 2)) + .02, 1.4, .19, angle = 20) panel.text(1.7, .22, labels = expression(lambda(bold(x)) <= c)) panel.curve(exp(- (n * x ** 2) / (2 * s ** 2)), .5, 2, lwd = 2, lty = 2, add = TRUE) from.z <- -2 to.z <- -.5 S.x <- c(from.z, seq(from.z, to.z, 0.01), to.z) S.y <- c(0, exp(- (n * (seq(from.z, to.z, 0.01)) ** 2) / (2 * s ** 2)), 0) panel.polygon(S.x, S.y, col = "gray", border = "white") panel.segments(-.5, 0, -2, 0) panel.segments(-.5, 0, -.5, .61) panel.arrows(-.7, .2, -1.2, .4, angle = 20) panel.text(-1.25, .46, labels = expression(paste(frac(bar(x), sigma/sqrt(n)) <= -c, "'", sep = ""))) panel.arrows(-1, exp(- (n * (-1) ** 2) / (2 * s ** 2)) + .02, -1.4, .19, angle = 20) panel.text(-1.7, .22, labels = expression(lambda(bold(x)) <= c)) panel.curve(exp(- (n * x ** 2) / (2 * s ** 2)), -2, -.5, lwd = 2, lty = 2, add = TRUE) } ) Notice that \Phi_1 is an increasing function, while \Phi_2 is decreasing as a function of \mu. We expect this since the alternative hypothesis is a two-sided one, so does the power. To see this, consider Figure 4 for different values of n.Figure 4: Two-Sided Power Function for Different n.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characterslibrary(magrittr) library(lattice) phi_fun2 <- function (x, mu0 = 0, s = 1, cr = 1.645, n = c(1, 4, 16, 64, 100), i = 1) { 1 - pnorm(cr - x / (s / sqrt(n[i]))) + pnorm(-cr - x / (s / sqrt(n[i]))) } x2 <- seq(-5, 5, length.out = 100) { phi_fun2(x2) ~ x2 } %>% xyplot( xlab = expression(mu), ylab = expression( paste( "1 - ", Phi, bgroup('[', c - frac(mu, sigma/sqrt(n)), ']'), ' + ' , Phi, bgroup('[', -c - frac(mu, sigma/sqrt(n)), ']') ) ), col = 'black', lwd = 2, type = 'l', panel = function (x, y, ...) { panel.grid(h = -1, v = -1) for (i in 2:5) { panel.curve(phi_fun2(x, i = i), -5, 5, add = TRUE, lty = i, lwd = 2) } panel.xyplot(x, y, ...) }, key = list(corner = c(.9, .1), text = list(c('n = 1', 'n = 4', 'n = 16', 'n = 64', 'n = 100')), lines = list(lty = c(1, 2, 3, 4, 5), lwd = 2) ) )