Lecture 31 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 8.2 to conclude that $T$ is sufficient for $\theta$. Since the family has full rank in the sense of definition 7.10 (the natural parameter space is the open half-line), we can use proposition 10.5 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$, so the likelihood equation has the unique root $\hat{\theta}=n/T=1/\bar{X}$ and this root is the maximum: $\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, theorem 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–Scheffe theorem 10.4 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 model is regular (exercise 16.1), the Cramer–Rao inequality (theorem 13.2) gives the bound $\theta^{2}/n$ for the variance of any unbiased estimator. The variance of the UMVUE is $\theta^{2}/(n-2)$, which is larger than the bound and by proposition 13.5 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: $\bar{X}$ is unbiased with variance $1/(n\theta^{2})$. This variance equals the bound $1/\mathcal{I}_{n}(\mu)$ for the mean (exercise 11.5) 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 25.7 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, written there in terms of the mean $1/\theta$, 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 25. 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]$. 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 25.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 22, 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 22.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 22.4, 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 20.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 20.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)$, so 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 22.5: the two one-sided tests reject on opposite tails of $T$. We use the generalised likelihood-ratio test of definition 23.1 instead. The statistic is

and by Wilks’ theorem 23.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$, so again $H_{0}$ is not rejected. The exercises at the end of this lecture verify 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 theorem 26.2: the posterior distribution is given by

which is again the Gamma$(a+n,b+T)$ distribution. The situation in the example fits the result of proposition 26.4: the update process increases $a$ by $n$ and $b$ by $T$, so that the effect of the prior is as if we had observed $a$ additional values with total lifetime $b$. In the example of the manufacturer’s claim about the lifetime of the components, we can assume that the claim is equivalent to having observed two more prior observations with mean lifetime two years, i.e. $a=2$ and $b=4$. The prior mean is then $a/b=0.5$. The posterior distribution is then Gamma$(12,29)$ and, using theorem 26.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 28.1 comes from the same chi-squared fact as before: if $\theta\sim\text{Gamma}(a+n,b+T)$ then $2(b+T)\theta\sim\chi^{2}_{2(a+n)}$, so 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 26.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 31.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 28.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 for any model, and it is this list which you should use when you are given 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 a luxury of the exponential distribution; in general you will only have the large-sample results from lectures 16, 17 and 23, 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 23, 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 26. The models under consideration here will be much larger, but the questions we ask of these models will be similar to the ones we considered in this module.