Example: Failure Times

An engineer is interested in understanding the failure time of a particular engine component. He undertakes an experiment where 10 independent copies of a component are tested separately and their lifetime (in hours) measured, leading to the following data:

      90 255  40 143  30 239 484  28  39  15

What conclusions can be made about the mean failure time?

Our first decision is what probability model to use for our data. This is an important part of the modelling process; formal model choice procedures will be described later in the module.

For now, our informal procedure, recalled from MATH235 is: plot the data first, and check the model fit after finding the MLE (e.g. using a QQ plot).

We first take a look at the data; for example we could use the R code

fail<-c(90,255,40,143,30,239,484,28,39,15)
stripchart(fail,pch=4,xlab="failure times (in hours)")

This leads to the strip plot in Figure 2.1.

For lifetime data of this sort, for reasons that were developed in the MATH230 course, a common model for such data (not least because the failure times are bound to be non-negative) is the exponential distribution, with rate parameter $\theta$ and mean ${\theta }^{-1}$.

Thus, we will assume that our data $x_1,x_2,\mathrm{\dots },x_n$ (where $n=10$, but we’ll work more generally for the moment) are independent realisations of independent random variables, $X_1,X_2,\mathrm{\dots },X_n$, each having the $\text{Exponential}(\theta )$ distribution.

Assuming the exponential distribution, for each of the $X_i$,

$$f(x_i|\theta )=\theta e^{-\theta x_i}\mathit{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.25em}}(x_i>0),\text{with }\theta \in \mathrm{\Omega }=(0,\mathrm{\infty }).$$

Hence:

	$L(\theta )$	$=$	${\displaystyle \prod _{i=1}^n}\theta e^{-\theta x_i}\mathit{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}}(\theta >0)$
		$=$	${\theta }^n{e}^{-\theta {\sum }_{i=1}^n{x}_i}$

and

$$\mathrm{\ell }(\theta )=n\log \theta -\theta \sum _{i=1}^n{x}_i=n\log \theta -\theta n\overline{x}.$$

So

$$S(\theta )={\mathrm{\ell }}^{\prime }(\theta )=\frac{n}{\theta }-n\overline{x}$$

and for a stationary point ${\mathrm{\ell }}^{\prime }(\widehat{\theta })=0$. Hence,

$$\widehat{\theta }=\frac{1}{\overline{x}}$$

which is verified as a maximum since

Thus, we obtain $\widehat{\theta }=1/\overline{x}$ as the maximum likelihood estimator of $\theta$: for our data this leads to $\widehat{\theta }=0.00734$.

The 95% confidence interval for $\theta$ based on the asymptotic normality of $\widehat{\theta }$ is given by

$$\widehat{\theta }\pm z_{\alpha /2}{[I_E(\widehat{\theta })]}^{-1/2}$$

where $z_{\alpha /2}=1.96$ and

	$I_E(\theta )$	$=$	$E\left[-{\displaystyle \frac{d^2\mathrm{\ell }(\theta ;\overrightarrow{X})}{d{\theta }^2}}\right]$
		$=$	$E\left[{\displaystyle \frac{n}{{\theta }^2}}\right]$
		$=$	${\displaystyle \frac{n}{{\theta }^2}}.$

So, replacing $I_E(\theta )$ by $I_E(\widehat{\theta })$ the interval becomes

$$0.00734\pm 1.96/\sqrt{10/{(0.00734)}^2}=(0.00279,0.0119).$$

For comparison, we can now calculate an approximate 95% confidence interval based on the deviance.

The deviance is given by

	$D(\theta )$	$=$	$2\left\{n\log {\displaystyle \frac{1}{\overline{x}}}-{\displaystyle \frac{n\overline{x}}{\overline{x}}}-n\log \theta +\theta n\overline{x}\right\}$
		$=$	$2\left\{10\log {\displaystyle \frac{1}{\overline{x}}}-10-10\log \theta +10\theta \overline{x}\right\}.$

From ${\chi }_1^2$ tables, for a 95% interval, we need to solve $D(\theta )=3.84$. This is a non-linear equation in $\theta$; to solve it we plot $D(\theta )$ and look for the intersections with the $D(\theta )=3.84$ line.

From Figure Figure 2.2 (Link), this leads to the interval (0.0037,0.013).

At this point, recall that inferences are required not about $\theta$, but about the mean of the population, $\varphi =1/\theta$, i.e. $g(\theta )=1/\theta$. We could go right back to the start of the problem, replace $\theta$ with $1/\varphi$, and do the whole thing again working with $\varphi$ instead of $\theta$.

However, because of the invariance property of the likelihood that’s not necessary: we get exactly the same results by simple transformation $\varphi =1/\theta$ of both the maximum likelihood estimate and the deviance-based confidence interval.

In other words,

$$\widehat{\varphi }=1/\widehat{\theta }=\overline{x}=136.3.$$

and the confidence interval becomes

$$({\tilde{\varphi }}_l,{\tilde{\varphi }}_u)=(1/{\tilde{\theta }}_u,1/{\tilde{\theta }}_l)=(76.9,270.3).$$

For comparison, the confidence interval for the mean parameter $\varphi$ based on the asymptotic normality of the MLE $\widehat{\varphi }$ is $(51.8,221)$.

Exercise 5: Verify that the confidence interval for $\varphi$ based on asymptotic normality of the MLE $\widehat{\varphi }$ is indeed $(51.8,221)$.

xbar<-mean(fail)

# function to calculate D(theta)
dev<-function(theta){
     2*(10*log(1/xbar) - 10 - 10*log(theta) + theta*10*xbar)
}

# a sequence of x-values at which to evaluate D(theta)
theta<-seq(from=0.001,to=0.02,length=1000)

plot(theta,dev(theta),type="l",ylab="deviance",xlim=c(0,0.02))


# add horizontal line
abline(h=3.84)

# find exact solutions via...
devroot<-function(theta){
     dev(theta)-3.84
}

# lower interval
uniroot(devroot,c(0.002,0.005))$root

# upper interval
uniroot(devroot,c(0.010,0.015))$root