Lecture 32 The Exponential Model: From Start to Finish

Assume that the lifetime of a component of a machine is exponentially distributed with rate $\theta>0$. Then the lifetimes $X_{1},\dots,X_{n}$ of $n$ components are i.i.d. with density $f(x;\theta)=\theta e^{-\theta x}$ for $x\geq 0$. The mean lifetime is $\mathbb{E}_{\theta}(X_{i})=1/\theta$, the variance is $1/\theta^{2}$ (table A.2). In example 7.4 we have seen that this is a one-parameter exponential family with natural statistic $T(x)=x$ and natural parameter $\eta(\theta)=-\theta$. Thus, by lemma 7.5, the joint density of the sample is

for $x_{1},\dots,x_{n}\geq 0$. Since this density depends on the data only through $T=\sum_{i=1}^{n}X_{i}$, we can use the factorisation theorem (proposition 8.3) to conclude that $T$ is sufficient for $\theta$. Since the family has full rank in the sense of definition 7.8 (the natural parameter space is the open half-line), we can use proposition 10.6 to conclude that $T$ is also complete. Thus, from lectures 7 to 10 we know that all our results can be based on $T$ alone, and that any unbiased function of $T$ will be the best unbiased estimator for the corresponding quantity.

As data we have observed $n=10$ lifetimes (in years) of the components, and the sum of these values is $\sum_{i}x_{i}=25.0$, i.e. we have $\bar{x}=2.5$. The manufacturer claims that the mean lifetime is two years, i.e. $\theta_{0}=0.5$, and we want to test whether the components of the machine last longer than the manufacturer claims.

The log-likelihood is $\ell(\theta)=n\log\theta-\theta T$, with derivative $\ell^{\prime}(\theta)=n/\theta-T$ and second derivative $-n/\theta^{2}<0$, and thus the likelihood equation has the unique root $\hat{\theta}=n/T=1/\bar{X}$, which is the maximum. Thus $\hat{\theta}$ is the MLE for $\theta$ (example 5.2). For our data we find $\hat{\theta}=10/25=0.4$. By the invariance of the MLE, proposition 5.6, the MLE for the mean lifetime $1/\theta$ is $\bar{X}=2.5$ years.

The MLE for the rate is biased: exercise 2.3 shows $\mathbb{E}_{\theta}(\hat{\theta})=n\theta/(n-1)$, so that $\hat{\theta}$ overestimates $\theta$ by a factor of $n/(n-1)$. The corrected estimator $\tilde{\theta}=(n-1)/T$ is unbiased. Since $\tilde{\theta}$ is an unbiased function of the complete sufficient statistic $T$, the Lehmann–Scheffé theorem 10.5 shows that this is the UMVUE for $\theta$. The variance of the UMVUE can be found as $\theta^{2}/(n-2)$ (exercise 5.3). For the data considered here, we find $\tilde{\theta}=9/25=0.36$.

What is the optimal performance of an estimator? From example 11.8 we know that the Fisher information of one observation is $\mathcal{I}(\theta)=1/\theta^{2}$ and thus we have $\mathcal{I}_{n}(\theta)=n/\theta^{2}$. Since the conditions of section 11.1 hold, the Cramér–Rao inequality (theorem 13.2) gives the bound $\theta^{2}/n$ for the variance of any unbiased estimator satisfying (13.1). The variance of the UMVUE is $\theta^{2}/(n-2)$, which is larger than the bound and by the discussion after corollary 13.6 no unbiased estimator can attain the bound, since the bound is only attained for an affine function of the natural statistic $T$ and no such function is unbiased for the rate. The efficiency $(n-2)/n$ of $\tilde{\theta}$ tends to $1$ as $n\to\infty$. The situation for the mean $1/\theta$ is different, since $\bar{X}$ is unbiased with variance $1/(n\theta^{2})$. This variance equals the bound for the mean (exercise 11.6) and thus $\bar{X}$ is efficient.

For large samples, theorems 16.3 and 17.1 apply, because the regularity conditions are satisfied: $\hat{\theta}_{n}$ is consistent and $\sqrt{n}(\hat{\theta}_{n}-\theta)\xrightarrow{\ \mathrm{d}\ }N(0,\theta^{2})$. From corollary 17.3 we know that the standard error is

For the data at hand we get $\mathop{\mathrm{se}}\nolimits(\hat{\theta})=0.4/\sqrt{10}=0.126$. This is the number we would expect to see next to the estimate.

Proposition 26.8 turns the standard error into the approximate $95\%$ confidence interval $\hat{\theta}\pm 1.96\,\mathop{\mathrm{se}}\nolimits(\hat{\theta})$. For the data considered here, this confidence interval is $0.4\pm 0.248$, i.e. $[0.152,0.648]$. The interval is symmetric about the estimate, and the probability of the interval covering the correct value is only approximately $95\%$, since the interval is based on the normal limit at $n=10$.

For the exponential distribution we can do better, because an exact pivot is available. In exercise 13.1 we have seen that $2\theta T\sim\chi^{2}_{2n}$ for every value of $\theta$, so $Q=2\theta T$ is a pivot in the sense of lecture 26. Solving $\chi^{2}_{2n,\alpha/2}\leq 2\theta T\leq\chi^{2}_{2n,1-\alpha/2}$ for $\theta$ gives the exact interval

with coverage exactly $1-\alpha$ for every $n$. With $\chi^{2}_{20,0.025}=9.591$ and $\chi^{2}_{20,0.975}=34.170$ we find the interval $[9.591/50,34.170/50]=[0.192,0.683]$, as in exercise 26.1. The exact interval is shifted to the right relative to the Wald interval and is not symmetric about $0.4$: it makes use of the fact that the sampling distribution of $\hat{\theta}$ is skewed, which the normal approximation does not. Both intervals are $95\%$ confidence intervals in the sense of definition 26.1: the statement is about the procedure, and for the interval in front of us the parameter is either inside the interval or it is not.

We now test the manufacturer’s claim. We take $H_{0}\colon\theta=\theta_{0}=0.5$ and start, as in lecture 23, with a simple alternative $H_{1}\colon\theta=\theta_{1}$ where $\theta_{1}<\theta_{0}$ (a longer mean lifetime). The likelihood ratio is

which is an increasing function of $T$, since $\theta_{0}-\theta_{1}>0$. By the Neyman–Pearson lemma, theorem 23.2, the most powerful test at level $\alpha$ rejects $H_{0}$ if $T>c$, where $c$ is chosen such that $\mathbb{P}_{\theta_{0}}(T>c)=\alpha$. Under $H_{0}$ the pivot gives $2\theta_{0}T\sim\chi^{2}_{2n}$, and thus $c=\chi^{2}_{2n,1-\alpha}/(2\theta_{0})$. For $\alpha=0.05$ we find $c=31.410/1=31.41$. The critical region $\{T>31.41\}$ does not depend on $\theta_{1}$ and thus, by proposition 23.4 with the inequalities reversed, the same test is most powerful against all alternatives $\theta<\theta_{0}$. Since the observed value $T=25.0$ is less than the critical value, we do not reject $H_{0}$ at the $5\%$ level: the data do not show that the components have a longer lifetime than claimed. The $p$-value from definition 22.8 is $\mathbb{P}_{\theta_{0}}(T\geq 25.0)=\mathbb{P}(\chi^{2}_{20}\geq 25.0)=0.20$.

The power function of definition 22.4 can be found by evaluating the same pivot at the alternative: $\beta(\theta)=\mathbb{P}_{\theta}(T>c)=\mathbb{P}(\chi^{2}_{2n}>2\theta c)$. For example $\beta(0.25)=\mathbb{P}(\chi^{2}_{20}>15.71)=0.73$. If the true mean lifetime of the components is four years, then a sample of ten components has probability $0.73$ of detecting this.

For the two-sided alternative $\theta\neq\theta_{0}$ there is no uniformly most powerful test, by the argument of proposition 23.5: the two one-sided tests reject on opposite tails of $T$. We use the generalised likelihood-ratio test of definition 25.1 instead. The statistic is

and by Wilks’ theorem 25.3 it is approximately $\chi^{2}_{1}$-distributed under $H_{0}$, since one parameter is fixed. For our data, $W=20\bigl{(}\log 0.8-1+1.25\bigr{)}=0.54$, well below $\chi^{2}_{1,0.95}=3.84$, and thus again $H_{0}$ is not rejected. Exercise 32.1 verifies the formula for $W$.

Since we can assume $\theta$ to be random with prior density $\pi(\theta)$, and since the likelihood function is $\theta^{n}e^{-\theta T}$ as a function of $\theta$, a prior of the same shape as the likelihood, namely the Gamma$(a,b)$ density $\pi(\theta)\propto\theta^{a-1}e^{-b\theta}$, is conjugate. This is a special case of proposition 28.4: the posterior distribution is given by

which is again the Gamma$(a+n,b+T)$ distribution. The update increases $a$ by $n$ and $b$ by $T$. For the natural parameter $\eta=-\theta$ the Jacobian is $1$, and thus the Gamma$(a,b)$ prior is the prior $\pi_{b,a-1}$ of proposition 28.4. In the example of the manufacturer’s claim about the lifetime of the components, we centre the prior at the claimed rate by choosing $a=2$ and $b=4$, so that the prior mean is $a/b=0.5=\theta_{0}$. This prior is weak, since the update adds the sample size $n=10$ to the small shape $a=2$, and thus the data outweigh the prior. The posterior distribution is then Gamma$(12,29)$ and, using lemma 28.7, the Bayes estimator for $\theta$, using squared error loss, is given by

This is a weighted average of the prior mean and the MLE, where the weight of the data increases as $T$ increases.

The credible interval of definition 29.1 comes from the same chi-squared fact as before: for integer $2a$, if $\theta\sim\text{Gamma}(a+n,b+T)$ then $2(b+T)\theta\sim\chi^{2}_{2(a+n)}$, and thus the equal-tailed $95\%$ credible interval is $[\chi^{2}_{24,0.025},\chi^{2}_{24,0.975}]/(2\cdot 29)=[0.214,0.679]$. For the manufacturer’s claim, the posterior probability of long lifetimes is $\mathbb{P}(\theta<0.5\mskip 1.0mu|\mskip 1.0mux)=\mathbb{P}(\chi^{2}_{24}<29)=0.78$, up from the prior probability $0.59$. Unlike the $p$-value, this is a probability about the parameter, computed after the data have been seen.

The Jeffreys prior from definition 28.8 is $\pi_{J}(\theta)\propto\sqrt{\mathcal{I}(\theta)}=1/\theta$. This is the improper limit $a,b\to 0$ of the Gamma prior. Using this prior, the posterior distribution is Gamma$(n,T)$. The credible interval is $[\chi^{2}_{2n,\alpha/2},\chi^{2}_{2n,1-\alpha/2}]/(2T)$. This is the same as the exact confidence interval introduced in section 32.3. Also, $\mathbb{P}(\theta\geq\theta_{0}\mskip 1.0mu|\mskip 1.0mux)=\mathbb{P}(\chi^{2}_{2n}\geq 2\theta_{0}T)$. This is the $p$-value of the one-sided test. While the numbers coincide, the statements are different: the confidence interval covers the fixed $\theta$ in $95\%$ of repetitions, whereas the credible interval covers the random $\theta$ with probability $0.95$ given the data; and the $p$-value is a probability about the data under $H_{0}$, whereas $\mathbb{P}(\theta\geq\theta_{0}\mskip 1.0mu|\mskip 1.0mux)$ is a probability about the hypothesis. This illustrates the comparison from section 29.3, for a concrete model, and it is special that the numbers coincide for the Jeffreys prior and one-sided hypotheses.

Every question about the exponential model we have been able to answer so far has used a small number of different tools: the exponential-family structure of the model allowed us to find a complete sufficient statistic, the likelihood function allowed us to find the estimator, the Fisher information allowed us to find the standard error and the bound, a pivot or the normal limit allowed us to find the interval, the Neyman–Pearson lemma and Wilks’ theorem allowed us to find tests, and Bayes’ theorem allowed us to find the posterior distribution. The same list of tools applies to most regular models, and it is this list which we can use for a model not considered in these notes: write down the likelihood function, find the sufficient statistic, solve the likelihood equation and check that the root is a maximum, compute the information, and use the results to compute the standard error, the interval, and the tests. The exact pivot is available because of the special form of the exponential distribution. In general we will only have the large-sample results from lectures 16, 17 and 25, and some care is needed when these results are used.

This is also the material we will study after this module. The standard errors computed by all regression packages are the inverse Fisher information from lecture 14, the tests which compare fitted models from two different analyses are the likelihood-ratio tests from lecture 25, mixture models and hidden Markov models are fitted using the EM algorithm from lecture 19, and Bayesian computation, from the simple conjugate updating to the simulation methods of later modules, is based on the posterior distribution from lecture 28. The models there will be much larger, but the questions we ask of these models will be similar to the ones we considered in this module.