2 Estimation

2.1 Data and questions

Data set 2.1 Wages in the USA

These data are taken from Dyer, D. (1981), Canadian Journal of Statistics, 9. They comprise annual wages (in multiples of US$100) of 30 production line workers in a large American firm.

Table 2.1 US wages

112 123 128 116 119
154 103 132 140 111
119 115 107 108 101
108 107 151 105 157
112 125 103 158 112
156 119 104 104 115

□

Data set 2.2 Severe ideopathic respiratory distress syndrome

The data were published by van Vliet, P.K. and Gupta, J.M. (1973), Sodium bicarbonate in ideopathic respiratory distress syndrome, Arch. Diseases in Childhood, 48, 249-255. The condition, known by its acronym SIRDS, is very serious and can result in the death. The data comprise birth weights (kg) of 50 infants who displayed SIRDS

. Table 2.2 Birth weights (kg) of infants with severe ideopathic respiratory distress syndrome 1.050* 2.500* 1.890* 1.760 2.830 1.175* 1.030* 1.940* 1.930 1.410 1.230* 1.100* 2.200* 2.015 1.715 1.310* 1.185* 2.270* 2.090 1.720 1.500* 1.225* 2.440* 2.600 2.040 1.600* 1.262* 2.560* 2.700 2.200 1.720* 1.295* 2.730* 2.950 2.400 1.750* 1.300* 1.130 3.160 2.550 1.770* 1.550* 1.575 3.400 2.570 2.275* 1.820* 1.680 3.640 3.005 * child died

There are important questions to be asked about these data. Can we identify children at risk quickly and accurately? Is it possible to relate the risk of death to birthweight? □ Data set 2.3 Silver content of Byzantine coins

A number of coins from the reign of King Manuel I, Comnemus (1143 - 80) were discovered in Cyprus. They arise from four different coinages at intervals throughout his reign. The question of interest is whether there is any significant difference in their silver content with the passage of time; there is a suspicion that it was deliberately and steadily reduced. The data give the silver content (%Ag) of the coins.

Table 2.3 Silver content of coins

First Second Third Fourth
5.9  6.9  4.9  5.3
6.8  9.0  5.5  5.6
6.4  6.6  4.6  5.5
7.0  8.1  4.5  5.1
6.6  9.3  6.2
7.7  9.2  5.8
7.2  8.6  5.8

6.9 6.2

On the face of it the suspicion could be correct in that the fourth coinage would seem to have lower silver content than, say, the first coinage, but there is a need for firm statistical evidence if it is to be confirmed. □

Data set 2.4 Radiocarbon dating

Radiocarbon datings have given historians a powerful tool for determining more precisely the periods during which ancient civilisations flourished. The following set of data came from the Lake Lamoka site and were contributed by S.R. Wilson from the Australian National University to Andrews, D.F. and Herzberg, A.M. (1985), Data, Springer-Verlag: New York.

Table 2.4 Radiocarbon dating

Sample Radiocarbon age
number determination
C-288  2419
M-26  2485
C-367  3433
M-195  2575
M-911  2521
M-912  2451
Y-1279  2550
Y-1280  2540

□

2.2 Desirable properties of estimators

Suppose we have X = X_{1}, X_{2}, …, X_{n} drawn from a distribution with some parameter θ. Often X_{1}, X_{2}, …, X_{n} form a sample.

Definition 2.0 Estimators

An estimator θ̂_{n} of θ is just a function of the observed data which (we hope) forms a useful approximation to the parameter:

θ̂_{n} = g(X_{1}, X_{2}, …, X_{n}).

Note that θ̂_{n} can depend only on the observed data, and not on any unknown parameters. You have already met a number of estimators in Mods, such as. the sample mean and sample variance. The estimator is a function of random variables, so is itself a random variable, with a distribution, mean, and variance, etc.

Here is a familiar property which you have already met.

Definition 2.1 Unbiasedness θ̂_{n} is said to be unbiased for θ if

E \binom{}{θ_{n}} = θ, ∀θ ∈ Θ.

□

Next a new property, but one which appeals to common-sense. The idea is that, the larger the amount of data, the closer the estimate should be to the parameter to be estimated. This is expressed in terms of the estimator converging in probability to the parameter value.

Convergence in probability (borrowed from probability course) A sequence of random variables Z_{1}, Z_{2}, …, Z_{n}, ... is said to converge in probability to a constant z if for any ϵ > 0, as n → ∞

P (|Z_{n} − z| > ϵ) → 0.

Definition 2.2 Consistency θ̂_{n} is said to be consistent for θ if

θ̂_{n} →^{P} θ.

□

Definition 2.3 Efficiency θ̂_{A} is said to be more efficient than θ̂_{B} if

V \binom{}{θ_{A}} < V \binom{}{θ_{B}} , ∀θ ∈ Θ.

□

Again, this appeals to common-sense.

2.3 Revision and extension of maximum likelihood estimation

The basic idea starts with the joint distribution of X = X_{1}, X_{2}, …, X_{n} depending upon a parameter θ,

f (x; θ) = f (x_{1}, x_{2}, …, x_{n}; θ) .

For fixed θ, probability statements can be made about X. If we have observations, x, but θ is unknown, we regard information about θ as being contained in the likelihood

L(θ; x) = f (x; θ) ,

where L is regarded as a function of θ with x fixed.

Example 2.1

Suppose X = X_{1}, X_{2}, …, X_{n} are independent Bernoulli random variables with parameter θ ∈ [0, 1].

i.e. P (X_{i} = 1) = θ, P (X_{i} = 0) = 1 − θ.

Observations are x = (1, 0, 0, 1, 0, 1, 1) and

L(θ; x) = ∏_{i=1}^{7}f (x_{i}; θ)

= ∏_{i=1}^{7}θ^{x_{i}} (1 − θ)^{1−x_{i}}

= θ^{4} (1 − θ)^{3} .

In general, for a sample size n,

L(θ; x) = ∏_{i=1}^{n} θ^{x_{i}} (1 − θ)^{1−x_{i}} = θ^{∑ x_{i}} (1 − θ)^{n−∑ x_{i}} .

■

Maximum likelihood estimation

The value of θ which maximises L(θ) is called the maximum likelihood estimate of θ.

Example 2.1 (continued)

In the previous example, we can find this by differentiating L(θ), or equivalently by differentiating l(θ) = log L(θ).

\begin{aligned}l(θ) & = & ∑_{i} x_{i} log θ + (n − ∑_{i} x_{i}) log (1 − θ) \\ \frac{∂l(θ)}{∂θ} & = & ∑_{i} x_{i}/ θ − (n − ∑_{i} x_{i})/ (1 − θ) \end{aligned}

Putting \frac{∂l(θ)}{∂θ} = 0 we obtain

θ = ∑_{i} x_{i}/ n.

This leads us to θ̂_{n}, the maximum likelihood estimator based on a sample of size n

θ̂_{n} = ∑_{i} X_{i}/ n.

■ Example 2.2 Suppose X = X_{1}, X_{2}, …, X_{n} are independent normal random variables, N (μ, σ^{2}).

\begin{aligned}L (μ, σ^{2}) & = & ∏_{i=1}^{n} \frac{1}{\sqrt{2πσ^{2}}} exp [−\frac{(x_{i} − μ)^{2}}{2σ^{2}} ] \\ & = & (2πσ^{2})^{− n/2} exp [−\frac{1}{2σ^{2}} ∑_{i=1}^{n} (x_{i} − μ)^{2}] \end{aligned}

and

\begin{aligned}l (μ, σ^{2}) & = & log L (μ, σ^{2}) \\ & = & −\frac{n}{2} log (2πσ^{2}) − \frac{1}{2σ^{2}} ∑_{i=1}^{n} (x_{i} − μ)^{2} \end{aligned}

Differentiating,

\begin{aligned}\frac{∂l (μ, σ^{2})}{∂μ} & = & \frac{1}{σ^{2}} ∑_{i=1}^{n} (x_{i} − μ) , \\ \frac{∂l (μ, σ^{2})}{∂σ^{2}} & = & −\frac{n}{2σ^{2}} + \frac{1}{2σ^{4}} ∑_{i=1}^{n} (x_{i} − μ)^{2} \end{aligned}

Equating these derivatives to zero results in

μ̂ = \overline{X}, σ̂^{2} = \frac{1}{n} ∑_{i=1}^{n} (X_{i} − μ̂)^{2} = \frac{1}{n} ∑_{i=1}^{n} (X_{i} − \overline{X})^{2} .

■

2.4 Order statistics

To expand our armoury of useful and descriptive statistics we shall need to know about order statistics.

Definition 2.4 Order statistics

Suppose that X_{1} , X_{2}, . . , X_{n} is a random sample then the full set of order statistics is the re-ordering X_{(1)} , X_{(2)}, . . , X_{(n)} such that X_{(1)} ≤ X_{(2)} ≤ . . ≤ X_{(n)}. In particular the minimum order statistic is X_{(1)} = min X_{i} and the maximum order statistic is X_{(n)} = max X_{i}. If the random sample is derived from a continuous distribution then all the inequalities will be strict. □

Example 2.3 US wages

The data in Table 2.1 are shown below, with the data having been ordered.

Table 2.1: US wages

101 107 112 119 140
103 107 112 119 151
103 108 115 123 154
104 108 115 125 156
104 111 116 128 157
105 112 119 132 158

You can see that the realisation of the minimum order statistic is 101 and the realisation of the maximum order statistic is 158. ■ The distributions of X_{(1)} = min {X_{i}} and X_{(n)} = max {X_{i}} are straightforward to derive. The c.d.f. of X_{(1)} is

F_{(1)}(x) = 1 − P (X_{(1)} > x)) = 1 − P (X_{1} > x, …, X_{n} > x)
= 1 − [1 − F (x)]^{n}

so that, differentiating,

f_{(1)}(x) = F′_{(1)}(x) = n [1 − F (x)]^{n−1} f(x).

Similarly

F_{(n)}(x) = P (X_{(n)} ≤ x) = F (x)^{n}

giving

f_{(n)}(x) = nF (x)^{n−1}f(x).

Example 2.4 Suppose X_{i} ∼ U(0, θ). Then

F (x) = \{ \begin{matrix}0, & x < 0, \\ \frac{x}{θ}, & 0 ≤ x < θ, \\ 1, & x ≥ θ, \end{matrix} and f(x) = \{ \begin{matrix}\frac{1}{θ}, & x ∈ (0, θ), \\ 0, & otherwise. \end{matrix}

The p.d.f.’s of X_{(1)} = min {X_{i}} and X_{(n)} = max {X_{i}} are given by

f_{(1)}(x) = \frac{n (θ − x)^{n−1}}{θ^{n}} , x ∈ (0, θ)

and

f_{(n)}(x) = \frac{nx^{n−1}}{θ^{n}} , x ∈ (0, θ).

■ For a continuous random sample, each with density function f(x) and cumulative distribution function F (x), the joint density function of the order statistics is

f(x_{(1)}, x_{(2)}, .., x_{(n)}) = n! ∏_{i=1}^{n} f(x_{(i)}), x_{(1)} < x_{(2)} < ⋯ < x_{(n)},

and the marginal distribution for X_{(r)} is

f_{(r)}(x) = \frac{n!}{(r − 1)!(n − r)!}F (x)^{r−1}(1 − F (x))^{n−r}f(x).

Lemma 2.1 The p.d.f. of X_{(r)} is given by

f_{(r)}(x) = \frac{n!}{(r − 1)!(n − r)!}F (x)^{r−1} [1 − F (x)]^{n−r} f(x).

□ Proof Putting r = 1 and r = n in the expression gives

f_{(1)}(x) = n [1 − F (x)]^{n−1} f(x), f_{(n)}(x) = nF (x)^{n−1}f(x).

The c.d.f. F_{(r)} of X_{(r)} is

F_{(r)}(x) = P (X_{(r)} ≤ x) = ∑_{j=r}^{n} \binom{n}{j} [F (x)]^{j} [1 − F (x)]^{n−j}

i.e. the probability that at least r of the X_{i}’s are less than or equal to x. Therefore

F_{(r)}(x) − F_{(r+1)}(x) = \binom{n}{r} [F (x)]^{r} [1 − F (x)]^{n−r}

and, by differentiation,

f_{(r+1)}(x) = f_{(r)}(x) − \binom{n}{r} [F (x)]^{r−1} [1 − F (x)]^{n−r−1} [r − nF (x)] f(x)

= \frac{n!}{r!(n − r)!}F (x)^{r−1} [1 − F (x)]^{n−r−1} (n − r)F (x)f(x)

= \frac{n!}{r! (n − (r + 1))!}F (x)^{(r+1)−1} [1 − F (x)]^{n−(r+1)} f(x),

so the proposition is proved inductively. ■

There are many uses of the order statistics in terms of analysing and describing data, three in particular being the median, the lower quartile and the upper quartile.

Definition 2.5 Median

The median of a random sample is defined by

X_{m} = X_{(_{\frac{1}{2}} (n+1)).}

This is the central order statistic if the size of the sample is odd and is the average of the two central order statistics if it is even, so that

if n = 2k + 1, then X_{m} = X_{(k+1)}, and if n = 2k, then X_{m} = \frac{1}{2} (X_{(k)} + X_{(k+1)}) .

□

Example 2.3 (revisited) US wages For the wages data, the realisation is x_{m} = \frac{1}{2} (115 + 115) . Thus 50% of the data have a smaller value than the median and 50% of the data have a higher value. ■

Definition 2.6 Quartiles

The lower quartile is defined by

q_{L} = X(_{\frac{1}{4}}_{(n+1)}),

and the upper quartile is defined by

q_{U} = X(_{\frac{3}{4}}_{(n+1)}).

□ Example 2.3 (revisited) US wages The wages data comprises 30 values and 31/4 = 7\frac{3}{4}, so the realisation of the lower quartile is q_{L} = x_{(7 \frac{3}{4})} = 107 + \frac{3}{4} (107 − 107) = 107. The realisation of the upper quartile is q_{U} = x_{(23 \frac{1}{4})} = 128 + \frac{1}{4} (132 − 128) = 129. ■

You will see one of the main uses of order statistics in the next sub-section, where we look at how to produce useful plots of the data.

2.4.1 Plotting the data

For univariate data there are many useful plots each of which attempts to give a graphical representation of the data. These include histograms, boxplots and plots of empirical cumulative distribution functions (see Rice, Chapter 10 for examples and descriptions of these plots). For paired data, scatter plots are a very quick way of visually checking for association between the paired observations.

Boxplots

Boxplots are particularly useful in giving an intuitive feel for the way the data are distributed. They are very easy to construct, but their main advantage is that all computer packages can produce them on-screen instantly. Let us look at how they are constructed.

Figure 2.1 Boxplot construction

With one exception the plot is largely self-explanatory; that exception is the idea of adjacent values. These are calculated by first calculating the inter-quartile range q_{U} −q_{L}; then the upper adjacent value is the largest value not exceeding that distance from the upper quartile, the lower adjacent value is the smallest value not less than that distance from the lower quartile. Any point which lies outside these boundaries is depicted as a point and considered as a possible outlier (i.e. a value which is non-typical for some reason or other). Thus we have a picture which, in some sense, represents the way the data are distributed.

Example 2.3 (revisited) US wages A boxplot of the wage data, produced using the R package, is shown in Figure 2.2.

Figure 2.2 Boxplot of US wages data

It shows some skewness in the data and some possible outliers, but overall you can see that most of the workers at that time earned between $10,000 and $14,000 with a median wage around $11,500 and 50% of workers earning between just over $10,500 and just under $13,000. ■

But boxplots are really useful when it comes to comparison.

Example 2.5 Infants with SIRDS

The data set in Table 2.2 is reproduced below.

. Table 2.2 Birth weights (kg) of infants with severe ideopathic respiratory distress syndrome 1.050* 2.500* 1.890* 1.760 2.830 1.175* 1.030* 1.940* 1.930 1.410 1.230* 1.100* 2.200* 2.015 1.715 1.310* 1.185* 2.270* 2.090 1.720 1.500* 1.225* 2.440* 2.600 2.040 1.600* 1.262* 2.560* 2.700 2.200 1.720* 1.295* 2.730* 2.950 2.400 1.750* 1.300* 1.130 3.160 2.550 1.770* 1.550* 1.575 3.400 2.570 2.275* 1.820* 1.680 3.640 3.005 * child died

One way of comparing the birthweights of infants who died with the birthweights of those who survived is to use a graph with the boxplots compared against the same scale.

Figure 2.3 Boxplots of SIRDS data

The first thing you notice about this plot is that the median value of birtheweight for children who died is less than the lower quartile for those who survived. The plot makes it immediately clear that the two groups have typically different birthweights. ■

Probability Plots

These are sometimes called quantile-quantile plots. They can be used to examine the plausibility of a density function from which a set of data could have been drawn. If such a known distribution is found then analysis of the data is performed more readily, particularly if this distribution should turn out to be normal. We need some theory before we develop the graphical method.

Definition 2.7

Let X be a continuous random variable with cumulative distribution function F_{X}(x), then the probability integral transform is

Y = F_{X}(X).

□

The transformed random variable Y has a particularly useful property.

Lemma 2.2 For F_{X}(x) strictly monotonic increasing on (a, b) , X taking values only

in (a, b) , then if Y = F_{X}(X),

Y ∼ U(0, 1).

Y so defined is uniformly distributed on (0, 1). □

Proof For 0 ≤ y ≤ 1,

P (Y ≤ y) = P (F_{X}(X) ≤ y) = P (X ≤ F_{X}^{−1}(y)) = F_{X} (F_{X}^{−1}(y)) = y,

giving the cumulative distribution function of the uniform distribution on (0, 1). ■

We shall also need a result concerning the expected values of order statistics from a uniform distribution.

Lemma 2.3 If Y_{(1)}, Y_{(2)}, …, Y_{(n)} are the order statistics of a random sample drawn from a uniform distribution U(0, 1), then

E (Y_{(k)}) = \frac{k}{n + 1} = P (Y_{i} ≤ \frac{k}{n + 1}) .

[Proof left as an exercise] □

We conclude that, for the uniform distribution, E (Y_{(1)}) , E (Y_{(2)}) , …, E(Y_{(n)}) splits the

interval (0, 1) into equal bites of probability, each of length \frac{1}{n + 1}.

Using the results of the two lemmas, it is now possible to give the method of construction for a probability plot. The question we wish to address is: ‘Is the data set plausibly drawn from a distribution with cdf F (x)?’ To construct a probability plot, you take the following steps:

(i) Arrange the data as x_{(1)} < x_{(2)} < ⋯ < x_{(n)}. From Lemma 2.2, if the supposition is correct, F (x_{(1)}), F (x_{(2)}), …, F (x_{(n)}) are realisations of order statistics from U(0, 1).

(ii) Solve

F (z_{(k)}) = \frac{k}{n + 1}

for z_{(k)}: remember that E (F (X_{(k)}) = \frac{k}{n + 1}.

(iii) Plot the ordered pairs (x_{(k)}, z_{(k)}).

If the supposition is reasonable we should get an approximate straight line, z = x.

If the distribution contains some unknown parameters we need to adapt the method, if possible.

Example 2.6 Suppose we wish to check whether a distribution is exponential with mean μ. Then we note F (x) = 1 − e^{−x/μ}, so that

\frac{k}{n + 1} ≈ 1 − e^{−x_{(k)}/μ}.

Re-arranging

x_{(k)} ≈ −μ log (1 − \frac{k}{n + 1}) .

Hence, if the distribution is plausibly exponential, plot (x_{(k)}, log (1 − \frac{k}{n+1})) and expect to see an approximate straight line. The slope will provide an estimate of μ. ■

Example 2.3 (revisited) US wages

Many financial modellers assume that wages have a Pareto distribution with c.d.f.

F (x) = 1 − (\frac{α}{x})^{θ} , x > α,

where α represents the minimum wage and both α and θ are unknown parameters. Writing

\frac{k}{n + 1} ≈ 1 − (\frac{α}{x_{(k)}} )^{θ}

and re-arranging,

x_{(k)} ≈ α (1 − \frac{k}{n + 1})^{−1/θ} .

Clearly we cannot produce a probability plot because we do not know the value of θ, but taking logs of both sides gives

log x_{(k)} ≈ log α − \frac{1}{θ} log (1 − \frac{k}{n + 1}) ,

so plotting ordered pairs

(log x_{(k)}, log (1 − \frac{k}{n + 1}))

should produce a straight line if the distribution is plausible as a model for the US wages data.

Figure 2.4 Pareto probability plot for US wages data

As you can see, a Pareto probability model doesn’t look very convincing. ■

In practice the distribution we wish to consider is frequently the normal distribution, which has two unknown parameters. We do a little more adjustment.

Definition 2.8 Normal scores

The numbers

z_{(k)} = Φ^{−1} (\frac{k}{n + 1}) ,

where Φ is the c.d.f. of Z ∼ N(0, 1), are called the normal scores. They are the \frac{k}{n + 1} quantiles of the standard normal distribution. □

We use the normal scores to construct a normal probability plot. Suppose that we wish to check whether or not the data could have been drawn from a normal distribution N (μ, σ^{2}) . As above take the ordered data x_{(1)} < x_{(2)} < ⋯ < x_{(n)}. Note that, if X ∼ N (μ, σ^{2}), then

\frac{X − μ}{σ} ∼ N(0, 1).

Plot the pairs(x_{(k)}, z_{(k)}) . If the data are plausibly normal then we should have an approximate straight line with slope 1/σ and intercept μ/σ.

In practice you will find that most statistical computer packages will generate normal scores and produce probability plots: in fact, most packages will do this for a variety of distributions. The plots in the examples which follow were generated using the R package.

Example 2.7 Silver content of Byzantine coins

The data are reproduced below.

Table 2.3 Silver content of coins

First Second Third Fourth
5.9  6.9  4.9  5.3
6.8  9.0  5.5  5.6
6.4  6.6  4.6  5.5
7.0  8.1  4.5  5.1
6.6  9.3  6.2
7.7  9.2  5.8
7.2  8.6  5.8

6.9 6.2

The graph shows a normal probability plot of the silver content of the first coinage.

Figure 2.5 Silver content of coins

It looks as if a normal distribution is appropriate. ■ Example 2.8 Radio-carbon dating

The data from Table 2.4 are reproduced below.

Table 2.4 Radiocarbon dating

Sample Radiocarbon age
number determination
C-288  2419
M-26  2485
C-367  3433
M-195  2575
M-911  2521
M-912  2451
Y-1279  2550
Y-1280  2540

Radioactive-carbon dating was undertaken on 8 samples from a single early site. There is one rather obvious outlier which is easily seen from the plot.

Figure 2.6 Radiocarbon dating

If possible we would check out the outlier with the laboratory and the collector but, since this is not possible, we will drop it and consider whether the remaining data points could have been drawn from a normal distribution.

Figure 2.7 Radiocarbon dating

As can be seen from Figure 2.7, the 7 remaining data points do give an approximate straight line on a normal probability plot. ■

2.5 Comparing estimators

Remember that we began with three desirable properties of estimators. These were

Unbiasedness E \binom{}{θ_{n}} = θ, ∀θ ∈ Θ.
Consistency θ̂_{n} →^{P} θ.
Efficiency  θ̂_{A} is more efficient than θ̂_{B} if V \binom{}{θ_{A}} < V \binom{}{θ_{B}} , ∀θ ∈ Θ.

Even if we assume unbiasedness and consistency to be desirable, it is possible to have more than one such estimator.

Example 2.9

Consider the linear estimator

θ̂_{n} = ∑_{i=1}^{n} a_{i}X_{i}

where E (X_{i}) = θ, V (X_{i}) = σ^{2} for 1 ≤ i ≤ n.

E \binom{}{θ_{n}} = ∑_{i=1}^{n} a_{i}E (X_{i}) = θ ∑_{i=1}^{n} a_{i},

so the estimator is unbiased provided

∑_{i=1}^{n} a_{i} = 1.

For i.i.d. random variables,

V \binom{}{θ_{n}} = ∑_{i=1}^{n} a_{i}^{2}σ^{2}.

Now

∑_{i=1}^{n} (a_{i} − \frac{1}{n})^{2} ≥ 0
⇒ ∑_{i=1}^{n} a_{i}^{2} − \frac{2}{n} ∑_{i=1}^{n} a_{i} + \frac{1}{n} ≥ 0
and, if ∑_{i=1}^{n} a_{i} = 1,
⇒  ∑_{i=1}^{n} a_{i}^{2} ≥ \frac{1}{n}.

Equality occurs iff a_{i} = \frac{1}{n}, 1 ≤ i ≤ n.

The conclusion then is that, if θ̂_{n} is a linear unbiased estimator of the form ∑_{i=1}^{n} a_{i}X_{i}

and if \overline{X} = \frac{1}{n} ∑_{i=1}^{n} X_{i}, then

V (\overline{X}) ≤ V \binom{}{θ_{n}} .

Of all linear unbiased estimators, \overline{X} is most efficient. ■

2.5.1 Food for thought

1. Is the sample mean the “best” estimator of the distribution mean?

2. If \overline{X} is unbiased for θ, is g (\overline{X}) unbiased for g (θ)?

3. Is there any reason to doubt the principle of seeking the unbiased estimator with the minimum variance?

4. Can a biased estimator be “better” (whatever that means!) than an unbiased estimator?

Is the sample mean the “best” estimator of the distribution mean?

Let X_{1}, X_{2}, …, X_{n} be a random sample from a uniform distribution U(θ − \frac{1}{2}, θ + \frac{1}{2}).

E (\overline{X}) = θ, V (\overline{X}) = \frac{1}{12n}.

Now consider θ̂_{n} = \frac{1}{2} (X_{(1)} + X_{(n)}). Obviously θ̂_{n} is symmetrically distributed about θ so

E \binom{}{θ_{n}} = θ.

Now

\begin{aligned}f_{X_{(1)}X_{(n)}}(u, v) & = & n(n − 1) [F (v) − F (u)]^{n−2} f(u)f(v), \\ & & θ − \frac{1}{2} < u < v < θ + \frac{1}{2}. \end{aligned}

where

F (x) = \{ \begin{matrix}0, & x < θ − \frac{1}{2}, \\ x + \frac{1}{2} − θ, & θ − \frac{1}{2} ≤ x < θ + \frac{1}{2}, \\ 1, & x ≥ θ + \frac{1}{2}, \end{matrix}

and

f(x) = F′(x),

so that

\begin{aligned}f_{X_{(1)}X_{(n)}}(u, v) & = & n(n − 1) [v − u]^{n−2} , \\ & & θ − \frac{1}{2} < u < v < θ + \frac{1}{2}. \end{aligned}

You have my word as a gentleman that, after a lot of boring slog,

V \binom{}{θ_{n}} = \frac{1}{2(n + 1)(n + 2)}

and

V \binom{}{θ_{n}} < V (\overline{X}) for n > 2.

Question: Does this mean we can sometimes do better by ignoring some of the data values?

If \overline{X} is unbiased for θ, is g (\overline{X}) unbiased for g (θ)? Let

f(x) = \frac{1}{θ}e^{− x/θ}, x ≥ 0.

Then

E (\overline{X}) = \frac{1}{n} ∑_{i=1}^{n} E (X_{i}) = θ, V (\overline{X}) = \frac{1}{n^{2}} ∑_{i=1}^{n} V (X_{i}) = \frac{θ^{2}}{n} .

Now consider g(x) = x^{2}.

E [g (\overline{X})] = E (\overline{X}^{2})

= V (\overline{X}) + E (\overline{X})^{2}

= \frac{θ^{2}}{n} + θ^{2} = θ^{2} (1 + \frac{1}{n}) .

Is there any reason to doubt the principle of seeking the unbiased estimator with the minimum variance? Finals question 1979 − IV − 15, part (iv).

X_{1}, X_{2}, …, X_{n} is a random sample from a Poisson distribution with mean λ. Hence

Y = ∑_{i=1}^{n} X_{i} ∼ P oisson (nλ) .

Part (iv) asks for an unbiased minimum variance estimator of e^{−2nλ}.

Let this be g(Y). Then

E [g(Y)] = ∑_{y=0}^{∞}g(y)\frac{e^{−nλ}(nλ)^{y}}{y!} = e^{−2nλ}

or

∑_{y=0}^{∞}g(y)\frac{(nλ)^{y}}{y!} = e^{−nλ}.

But

e^{−x} = ∑_{k=0}^{∞} (−1)^{k} \frac{x^{k}}{k!}

so that

g(y) = (−1)^{y} , y = 0, 1, …

i.e.

g(y) = \{ \begin{matrix}−1, & y odd, \\ 1, & y even. \end{matrix}

THIS IS COMPLETE AND UTTER RUBBISH!

Can a biased estimator be “better” than an unbiased estimator? Look at the distributions of the following two estimators of θ.

Figure 2.8

Which do you prefer?

Even though θ̂_{1} is biased,

E [(θ̂_{1} − θ)^{2}] < V \binom{}{θ_{2}} .

θ̂_{1} has smaller mean squared error.

Are m.l.e.’s the answer? The maximum likelihood principle is intuitively appealing and does not involve worries about bias.

(i) Maximum likelihood estimators are asymptotically unbiased.

(ii) If θ̂ is the m.l.e. of θ, then g \binom{}{θ}is the m.l.e. of g (θ).

This is the invariance property of maximum likelihood estimators.

Are there any problems with m.l.e.’s? Consider X_{1}, X_{2}, …, X_{n} where X_{i} ∼ U(θ− \frac{1}{2}, θ + \frac{1}{2}), 1 ≤ i ≤ n, and note that θ − \frac{1}{2} ≤ x_{(1)} < x_{(n)} ≤ θ + \frac{1}{2}.

The likelihood function is

L(θ; x) = \{ \begin{matrix}1, & x_{(n)} − \frac{1}{2} ≤ θ ≤ x_{(1)} + \frac{1}{2}, \\ 0, & otherwise. \end{matrix}

Clearly

max_{θ} \{L(θ; x)\} = 1 ∀θ ∈ [x_{(n)} − \frac{1}{2}, x_{(1)} + \frac{1}{2}]

i.e. the likelihood function is maximised for every statistic θ̂(X_{1}, X_{2}, …, X_{n}) such that θ̂(X) ∈ [X_{(n)} − \frac{1}{2}, X_{(1)} + \frac{1}{2}] . Conclusion: An m.l.e. may not be unique.

2.6 More about likelihood

2.6.1 Invariance property of m.l.e.’s Lemma 2.4 If θ̂ is an m.l.e. of θ and if g is a function, then g \binom{}{θ} is an m.l.e. of g(θ). □

Proof If g is one-to-one, then

L(θ) = L (g^{−1} (g(θ)))

are both maximised by θ̂, so

θ̂ = g^{−1} \binom{}{g(θ)}

or

g \binom{}{θ} = ĝ(θ).

If g is many-to-one, then θ̂ which maximises L(θ) still corresponds to g \binom{}{θ}, so g \binom{}{θ} still corresponds to the maximum of L(θ) ■ Example 2.10 Suppose X_{1}, X_{2}, …, X_{n} is a random sample from a Bernoulli distribution B(1, θ). Consider m.l.e.’s of the mean, θ, and variance, θ(1 − θ). Note, by the way, that θ(1 − θ) is not a 1-1 function of θ. The log-likelihood is

l(θ) = ∑_{i} x_{i} log θ + (n − ∑_{i} x_{i}) log(1 − θ)

and

\frac{dl(θ)}{dθ} = ∑_{i} x_{i}/ θ − (n − ∑_{i} x_{i})/ (1 − θ)

so it is easily shown that the m.l.e. of θ is θ̂ = \overline{X.} Putting ν = θ(1 − θ),

\frac{dl(ν)}{dν} = \frac{dl (ν(θ))}{dθ} · \frac{dθ}{dν}

so it is easily seen that, since \frac{dθ}{dν} is not, in general, equal to zero,

ν̂ = ν \binom{}{θ} = \overline{X} (1 − \overline{X}) .

■

2.6.2 Relative likelihood

If sup_{θ} L(θ) < ∞, the relative likelihood is

RL(θ) = \frac{L(θ)}{sup_{θ} L(θ)}; 0 ≤ RL(θ) ≤ 1.

Relative likelihood is invariant to known 1-1 transformations of x, for if y is a 1-1 function of x,

f_{Y} (y; θ) = f_{X} (x(y); θ) |\frac{dx}{dy}| .

|\frac{dx}{dy}| is independent of θ, so

RL_{X}(θ) = RL_{Y} (θ).

2.6.3 Likelihood summaries

Realistic statistical problems often have many parameters. These cause problems because it can be hard to visualise L(θ), and it becomes necessary to use summaries.

Key idea In large samples, log-likelihoods are often approximately quadratic near the maximum.

Example 2.11

Suppose X_{1}, X_{2}, …, X_{n} is a random sample from an exponential distribution with parameter λ. i.e.

f_{X}(x) = λe^{−λx}, x ≥ 0.

Then

\begin{matrix}l(λ) = & n log λ − λ ∑_{i} x_{i}, & \frac{dl(λ)}{dλ} = & n/ λ − ∑_{i} x_{i}, \\ \frac{d^{2}l(λ)}{dλ^{2}} = & − n/ λ^{2}, & \frac{d^{3}l(λ)}{dλ^{3}} = & 2n/ λ^{3}. \end{matrix}

The log-likelihood has a maximum at λ̂ = n /∑_{i} x_{i}, so

RL(λ) = (\frac{λ}{λ̂})^{n} e^{n−λ ∑_{i} x_{i}}
= (\frac{λ}{λ̂}e^{1− λ/λ̂})^{n} , λ > 0.
→ 1 as λ → λ̂.

Now, what does the likelihood look like in the neighbourhood of λ̂, as n → ∞?

log RL(λ) = l(λ) − l(λ̂)

= l(λ̂) + l′(λ̂) (λ − λ̂) + \frac{1}{2}l′′(λ̂) (λ − λ̂)^{2} − l(λ̂)

+ O(λ − λ̂)^{3}

using Taylor series. Now l′(λ̂) = 0 and l′′(λ̂) = −n /λ̂^{2} , so

log RL(λ) ≃ −\frac{n (λ − λ̂)^{2}}{2λ̂^{2}} → −∞ as n → ∞

unless λ = λ̂. Thus, as n → ∞,

RL(λ) → \{ \begin{matrix}1, & λ = λ̂, \\ 0, & otherwise. \end{matrix}

Conclusion Likelihood becomes more concentrated about the maximum as n → ∞, and values far from the maximum become less and less plausible. ■

In general We call the value θ̂ which maximises L(θ) or, equivalently, l(θ) = log L(θ) the maximum likelihood estimate, and

J(θ) = −\frac{∂^{2}l(θ)}{∂θ^{2}}

is called the observed information. Usually J(θ) > 0 and J(θ̂) measures the concentration of l(θ) at θ̂. Close to θ̂, we summarise

l(θ) ≃ l(θ̂) − \frac{1}{2} (θ − θ̂)^{2} J(θ̂).

2.6.4 Information In a model with log-likelihood l(θ), the observed information is

J(θ) = −\frac{∂^{2}l(θ)}{∂θ^{2}}.

When observations are independent, L(θ) is a product of densities so

l(θ) = ∑_{i} log f (x_{i}; θ)

and

J(θ) = − ∑_{i} \frac{∂^{2}}{∂θ^{2}} log f (x_{i}; θ) .

Since

l(θ) ≃ l(θ̂) − \frac{1}{2} (θ − θ̂)^{2} J(θ̂),

for θ near to θ̂, we see that large J(θ̂) implies that l(θ) is more concentrated about θ̂. This means that the data are less ambiguous about possible values of θ, i.e. we have more information about θ.

2.6.5 Expected information Univariate distributions Before an experiment is conducted, we have no data so that we cannot evaluate J(θ). But we can find its expected value

I(θ) = E (−\frac{∂^{2}l(θ)}{∂θ^{2}}) .

This is called the expected information or Fisher’s information. If the observations are a random sample, then the whole sample expected information is

I(θ) = ni(θ)

where

i(θ) = E (−\frac{∂^{2}}{∂θ^{2}} log f (X_{i}; θ)) ,

the single observation Fisher information.

Sir Ronald Aylmer Fisher (1890 - 1962)

Example 2.12

Suppose X_{1}, X_{2}, …, X_{n} is a random sample from a Poisson distribution with parameter θ.

L(θ) = ∏_{i=1}^{n} \frac{θ^{x_{i}}e^{−θ}}{x_{i}!} ,

giving

l(θ) = log L(θ) = ∑_{i} x_{i} log θ − nθ − ∑_{i} log x_{i}!

Thus

J(θ) = −\frac{∂^{2}l(θ)}{∂θ^{2}} = ∑_{i} x_{i}/ θ^{2}.

To find I(θ), we need E (X_{i}) = θ and

I(θ) = \frac{1}{θ^{2}} ∑_{i} E (X_{i}) = \frac{n}{θ}.

■

Multivariate distributions If θ is a (p × 1) vector of parameters, then I(θ) and J(θ) are (p × p) matrices.

\{J(θ)\}_{rs} = −\frac{∂^{2}l(θ)}{∂θ_{r}∂θ_{s}} and \{I(θ)\}_{rs} = E (−\frac{∂^{2}l(θ)}{∂θ_{r}∂θ_{s}} ) .

These matrices are obviously symmetric.

We can also write the above as

J(θ) = −\frac{∂^{2}l(θ)}{∂θ∂θ^{T}} and I(θ) = E (−\frac{∂^{2}l(θ)}{∂θ∂θ^{T}} ) .

Example 2.13 X_{1}, X_{2}, …, X_{n} is a random sample from a normal distribution with parameters μ and σ^{2}. We have already seen that

L (μ, σ^{2}) = (2πσ^{2})^{− n/2} exp [−\frac{1}{2σ^{2}} ∑_{i} (x_{i} − μ)^{2}] ,

so

l (μ, σ^{2}) = −\frac{n}{2} log 2π − \frac{n}{2} log σ^{2} − \frac{1}{2σ^{2}} ∑_{i} (x_{i} − μ)^{2} .

and

\begin{aligned}\frac{∂l}{∂μ} & = & \frac{1}{σ^{2}} ∑_{i} (x_{i} − μ) , \\ \frac{∂l}{∂σ^{2}} & = & −\frac{n}{2σ^{2}} + \frac{1}{2σ^{4}} ∑_{i} (x_{i} − μ)^{2} , \\ \frac{∂^{2}l}{∂μ^{2}} & = & −\frac{n}{σ^{2}} , \\ \frac{∂^{2}l}{∂μ∂σ^{2}} & = & −\frac{1}{σ^{4}} ∑_{i} (x_{i} − μ) . \\ \frac{∂^{2}l}{∂ (σ^{2})^{2}} & = & \frac{n}{2σ^{4}} − \frac{1}{σ^{6}} ∑_{i} (x_{i} − μ)^{2} . \end{aligned}

J(μ, σ^{2}) = ( \begin{matrix}\frac{n}{σ^{2}} & \frac{1}{σ^{4}} ∑_{i} (x_{i} − μ) \\ \frac{1}{σ^{4}} ∑_{i} (x_{i} − μ) & \frac{1}{σ^{6}} ∑_{i} (x_{i} − μ)^{2} − \frac{n}{2σ^{4}} \end{matrix} ) .

To find I(μ, σ^{2}), use

E (X_{i}) = μ,
V (X_{i}) = E [(X_{i} − μ)^{2}] = σ^{2},

so that

I(μ, σ^{2}) = E (J(μ, σ^{2})) = ( \begin{matrix}\frac{n}{σ^{2}} & 0 \\ 0 & \frac{n}{2σ^{4}} \end{matrix} ) .

Example 2.14 Censored exponential data Lifetimes of n components, safety devices, etc. are observed for a time c, when r have failed and (n − r) are still OK. We have two kinds of observation:

(i) Exact failure times x_{i} observed if x_{i} ≤ c, so that f(x; λ) = λe^{−λx}, x ≥ 0;

(ii) x_{i} unobserved if x_{i} > c,

P (X > c) = e^{−λc}.

Data are therefore x_{1}, …, x_{r}, \underbrace{c, …, c}_{n−r times}

The (n−r) components, safety devices, etc. which have not failed are said to be censored. The likelihood is

L(λ) = ∏_{i=1}^{r} λe^{−λx_{i}} ∏_{i=r+1}^{n} e^{−λc}

= λ^{r} exp [−λ (∑_{i=1}^{r} x_{i} + (n − r)c)] .

\begin{aligned}l(λ) & = & r log λ − λ (∑_{i=1}^{r} x_{i} + (n − r)c) \\ l′(λ) & = & r/ λ − (∑_{i=1}^{r} x_{i} + (n − r)c) \\ l′′(λ) & = & −r /λ^{2} . \end{aligned}

Thus J(λ) = r /λ^{2} > 0 if r > 0 so we must observe at least one exact failure time.

I(λ) = E (r /λ^{2} ) = \frac{1}{λ^{2}} E (#X_{i} observed exactly.)

Now P (X_{i} observed exactly) = P (X_{i} ≤ c) = 1 − e^{−λc}, so

I_{c}(λ) = \frac{n (1 − e^{−λc})}{λ^{2}}.

No censoring if c → ∞, giving

I_{∞}(λ) = \frac{n}{λ^{2}} > I_{c}(λ)

as one might expect.

The asymptotic efficiency when there is censoring at c relative to no censoring is

I_{c}(λ) /I_{∞}(λ) = 1 − e^{−λc}.

■

Revision Example 2.15 Events in a Poisson process

Events are observed for period (0, T).

n events occur at times 0 < t_{1} < t_{2} < … < t_{n} < T

Two observers A and B. A records exact times, B uses an automatic counter and goes to the pub (i.e. B merely records how many events there are).

A knows exact times, and times between events are independent and exponentially distributed, so

L_{A}(λ) = λe^{−λt_{1}} × λe^{−λ(t_{2}−t_{1})} × ⋯ × λe^{−λ(t_{n}−t_{n−1})} × e^{−λ(T −t_{n})}
= λ^{n}e^{−λT} .

B merely observes the event [N = n], where N ∼ P oi(λT), so

L_{B}(λ) = \frac{(λT)^{n} e^{−λT}}{n!} .

Log-likelihoods are

l_{A}(λ) = n log λ − λT,
l_{B}(λ) = n log λ + n log T − λT − log n!

and

J_{A}(λ) = J_{B}(λ) = n /λ^{2} .

E(N) = λT, so I_{A}(λ) = I_{B}(λ) = T /λ, and both observers get the same information. As usual, the one who went to the pub did the right thing. ■

2.6.6 Maximum likelihood estimates

The maximum likelihood estimate θ̂ of θ maximises L(θ) and often (but not always) satisfies the likelihood equation

\frac{∂l}{∂θ} \binom{}{θ} = 0,

with

J(θ̂) = −\frac{∂^{2}l}{∂θ^{2}} \binom{}{θ} > 0

for a maximum.

In the vector case, θ̂ solves simultaneously

\frac{∂l}{∂θ_{r}} \binom{}{θ} = 0, r = 1, …, p,

with J(θ̂) positive definite This implies

det J(θ̂) > 0

If the likelihood equation has many solutions, we find them all and check L(θ) for each. Usually, the equation has to be solved numerically. One way is by Newton-Raphson.

Suppose we have a starting value θ_{0}. Then

0 = \frac{∂l}{∂θ} \binom{}{θ} ≃ \frac{∂l}{∂θ} (θ_{0}) + \frac{∂^{2}l}{∂θ^{2}} (θ_{0}) (θ̂− θ_{0})

which may be re-arranged to

θ̂ = θ_{0} + \frac{U(θ_{0})}{J(θ_{0})},

where

U(θ) = \frac{∂l}{∂θ} is the score function,
J(θ) = −\frac{∂^{2}l}{∂θ^{2}} is the observed information.

Now we iterate using θ_{0} as a starting value and

θ_{n+1} = θ_{n} + \frac{U(θ_{n})}{J(θ_{n})}.

Example 2.16 Extreme value (Gumbel) distribution

This distribution is used to model such things as annual maximum temperature. Data due to Bliss on numbers of beetles killed by exposure to carbon disulphide are fitted by this model. The c.d.f. is

F (x) = exp (−e^{−(x−η)}) , x ∈ ℝ, η ∈ ℝ,

and the density is

f(x) = exp [−(x − η) − e^{−(x−η)}] , x ∈ ℝ, η ∈ ℝ.

The sample log-likelihood is

l(η) = − ∑_{i}(x_{i} − η) − ∑_{i} e^{−(x_{i}−η)},

so that

U(η) = n − ∑_{i} e^{−(x_{i}−η)},
J(η) = ∑_{i} e^{−(x_{i}−η)}.

Starting at η_{0} = \overline{x}, iterate using

η_{n+1} = η_{n} + \frac{n − ∑_{i} e^{−(x_{i}−η_{n})}}{∑_{i} e^{−(x_{i}−η_{n})}.}

■ Fisher scoring This simply involves replacing J(θ) with I(θ).

Example 2.17 Extreme value distribution We need

I(η) = E [J(η)] = ∑_{i} E [e^{−(X_{i}−η)}]
= n ∫_{−∞}^{∞} e^{−(x−η)} exp [−(x − η) − e^{−(x−η)}] dx.

Put u = e^{−(x−η)} and the integral becomes

I(η) = n ∫_{0}^{∞} ue^{−u}du = n,

so Fisher scoring gives the iteration

η_{n+1} = η_{n} + 1 − \frac{1}{n} ∑_{i} e^{−(x_{i}−η_{n})}.

■ 2.6.7 Sufficient statistics You have already seen a likelihood which cannot be summarised by a quadratic. Example 2.18

f (x_{i}; θ) = θ^{−1}, 0 < x_{i} < θ,

so

L(θ) = θ^{−n}, 0 < max {x_{i}} < θ.

Figure 2.9

Clearly a quadratic approximation is useless here. Definition 2.9 Sufficient statistic If S = s (X) is such that the conditional density f_{X|S} (x|s; θ) is independent of θ, then S is a sufficient statistic for θ. The important question is:

Does s(x) reduce the dimensionality of the problem?□

As we will see, the definition is equivalent to saying that the likelihood L(θ) only depends upon data x through s (x). So maximum likelihood inference also only depends on s (x) .

Example 2.19 Suppose X_{1}, X_{2} ∼ B(n, θ) and consider

P (X_{1} = x|X_{1} + X_{2} = r)

= \frac{P (X_{1} = x, X_{1} + X_{2} = r)}{P (X_{1} + X_{2} = r)}

= \frac{P (X_{1} = x, X_{2} = r − x)}{P (X_{1} + X_{2} = r)}

= \frac{\binom{n}{x}θ^{x} (1 − θ)^{n−x} \binom{n}{r − x}θ^{r−x} (1 − θ)^{n−r+x}}{\binom{2n}{r}θ^{r} (1 − θ)^{2n−r}}

= \begin{matrix}\binom{n}{x} & \binom{n}{r − x} \\ ( & {2n \atop r} & ) \end{matrix}.

This does not contain θ, so that X_{1} + X_{2} is a sufficient statistic for θ. ■

Theorem 2.1 Factorization Theorem s(X) is a sufficient statistic for θ if and only if there exist functions g and h such that

f(x; θ) = g (s(x); θ) h(x)

for all x ∈ ℝ^{n}, θ ∈ Θ. □ Proof for discrete random variables

(i) Let s(x) = a and suppose the factorization condition to be satisfied, so that f(x; θ) = g (s(x); θ) h(x). Then

P (s(X) = a) = ∑_{y∈s^{−1}(a)} p(y) = g(a; θ) ∑_{y∈s^{−1}(a)} h(y).

Hence

P (X = x | s(X) = a) = \frac{h(x)}{∑_{y∈s^{−1}(a)} h(y)}

and this does not depend upon θ.

(ii) Let s(X) be a sufficient statistic for θ. Then

P (X = x) = P (X = x | s(X) = a) P (s(X) = a) .

But sufficiency ⇒ P (X = x | s(X) = a) does not depend upon θ, so writing

P (s(X) = a) = g(a; θ) and P (X = x | s(X) = a) = h(x) gives the result.

The proof in the continuous case requires measure theory and is beyond the scope of this course. ■

Example 2.20

X_{1}, X_{2}, …, X_{n} ∼ U(0, θ), so that

L(θ) = θ^{−n}, 0 < x_{1}, …, x_{n} < θ.
L(θ) = θ^{−n}, θ > x_{(n)}.

This factorizes with s(x) = x_{(n)} and h (x) = 1, so that X_{(n)} is a sufficient statistic for θ.■

Example 2.21

Suppose X_{1}, X_{2}, …, X_{n} is a random sample from a Bernoulli distribution. Then

p(x; θ) = θ^{∑_{i} x_{i}}(1 − θ)^{n−∑_{i} x_{i}}.

Trivially this factorizes with s(x) = ∑_{i} x_{i} and h (x) = 1 ■ Example 2.22

Suppose X_{1}, X_{2}, …, X_{n} is a random sample from a N(μ, σ^{2}) distribution, where (μ, σ^{2})^{T} is a vector of unknown parameters. Then

f(x; μ, σ^{2}) = (2πσ^{2})^{− n/2} exp [−\frac{1}{2σ^{2}} ∑_{i}(x_{i} − μ)^{2}]

= (2πσ^{2})^{− n/2} exp [−\frac{1}{2σ^{2}} ∑_{i}(x_{i} − \overline{x})^{2} + n(\overline{x} − μ)^{2}] .

Again this factorizes where s(x) = (\overline{x}, ∑_{i}(x_{i} − \overline{x})^{2})^{T}, a vector valued function. ■ 2.6.8 Large sample distribution of θ̂ From the data summary point of view, the m.l.e. θ̂ and J \binom{}{θ}have been thought of in terms of a particular set of data. We now wish to think of θ̂ in terms of repeated sampling (i.e. as a random variable).

Main results In many situations and subject to regularity conditions

θ̂ →^{D} N(θ, I(θ)^{−1}),

We can apply the approach to obtaining confidence intervals developed in Mods, for the CLT. Then an approximate 95% confidence interval for θ is given by

θ̂± 1.96I(θ̂)^{−1/2} .

[or θ̂± 1.96J(θ̂)^{−1/2} , regarded by many as better, but not in the books].

In the multivariate case,

θ̂ →^{D} N(θ, I(θ)^{−1}).

Example 2.23 Exponential distribution For an exponential distribution with mean θ,

L(θ) = θ^{−n}e^{− ∑_{i} x_{i}/θ} , θ > 0,
l(θ) = −n log θ − ∑_{i} x_{i} /θ ,

so that

U(θ) = −\frac{n}{θ} + \frac{∑_{i} x_{i}}{θ^{2}} , J(θ) = −\frac{n}{θ^{2}} + \frac{2 ∑_{i} x_{i}}{θ^{3}} .

Thus

θ̂ = \overline{x}, J(θ̂) = \frac{n}{\overline{x}^{2}}

and an approximate 95% confidence interval is

\overline{x} ± 1.96\overline{x} /\sqrt{n}

■ Example 2.24 Normal distribution For a normal random sample,

μ̂ = \overline{x}, σ̂^{2} = n^{−1} ∑_{i} (x_{i} − \overline{x})^{2} ,

J (μ, σ^{2}) = ( \begin{matrix}n /σ^{2} & σ^{−4} ∑_{i} (x_{i} − μ) \\ σ^{−4} ∑_{i} (x_{i} − μ) & σ^{−6} ∑_{i} (x_{i} − μ)^{2} − n /2σ^{4} \end{matrix} ) .

Therefore

I (μ̂, σ̂^{2}) = ( \begin{matrix}n /σ̂^{2} & 0 \\ 0 & n /2σ̂^{4} \end{matrix} ) .

An approximate 95% confidence interval for μ is

\overline{x} ± 1.96σ̂ /\sqrt{n} ,

and for σ^{2} is

σ̂^{2} ± 1.96σ̂^{2}\sqrt{\frac{2}{n}}.

Note that the estimators μ̂ and σ̂^{2} are asymptotically uncorrelated. ■

Lemma 2.5 If θ̂ is the maximum likelihood estimator of a parameter θ based on a random sample, under suitable regularity conditions

θ̂ →^{D} N(θ, I(θ)^{−1}),

where I(θ) is Fisher’s information for the sample. □

Sketch Proof Suppose X_{1}, X_{2}, …, X_{n} is a random sample from a distribution with p.d.f. f(x; θ). Then the log-likelihood, score and observed information are

l(θ) = ∑_{i} log f(x_{i}; θ),
U(θ) = ∑_{i} \frac{∂}{∂θ} log f(x_{i}; θ),

J(θ) = − ∑_{i} \frac{∂^{2}}{∂θ^{2}} log f(x_{i}; θ).

Let U_{i}(θ) be the random variable U_{i}(θ) = \frac{∂}{∂θ} log f(X_{i}; θ), and, provided that conditions are such that integration and differentiation are interchangeable,

E [U_{i}(θ)] = ∫ f(x; θ)\frac{∂}{∂θ} log f(x; θ)dx

= ∫ \frac{∂}{∂θ}f(x; θ)dx

= \frac{∂}{∂θ} ∫ f(x; θ)dx = \frac{∂}{∂θ}1 = 0

and

0 = \frac{∂}{∂θ} ∫ f(x; θ)\frac{∂}{∂θ} log f(x; θ)dx

= ∫ f(x; θ)\frac{∂^{2}}{∂θ^{2}} log f(x; θ)dx + ∫ \frac{∂}{∂θ}f(x; θ)\frac{∂}{∂θ} log f(x; θ)dx

= E [\frac{∂^{2}}{∂θ^{2}} log f(X; θ)] + ∫ f(x; θ) (\frac{∂}{∂θ} log f(x; θ))^{2} dx.

So

0 = −i(θ) + E [U_{i}(θ)^{2}]

and, therefore, V [U_{i}(θ)] = i(θ). It follows that E [U(θ)] = 0, V [U(θ)] = ni(θ) = I(θ), and the Central Limit Theorem shows that

U(θ) →^{D} N (0, I(θ)) .

Now the m.l.e. is a solution of U(θ̂) = 0, so that, Taylor expanding about θ,

U(θ) + U′(θ)(θ̂− θ) ≃ 0

or

U(θ) − J(θ)(θ̂− θ) ≃ 0.

Re-arranging,

\sqrt{I(θ)}(θ̂− θ) ≃ U(θ)\frac{\sqrt{I(θ)}}{J(θ)} = \frac{U(θ)}{\sqrt{I(θ)}} /\frac{J(θ)}{I(θ)} .

From the CLT,

\frac{U(θ)}{\sqrt{I(θ)}} →^{D} N(0, 1)

and (from WLLN which you will meet in probability)

\frac{J(θ)}{I(θ)} →^{P} 1.

Slutsky’s Theorem states that, if S_{n} →^{D} S and U_{n} →^{P} k, where k is a constant, then

S_{n} + U_{n} →^{D} S + k, S_{n}U_{n} →^{D} Sk.

It therefore results in

\sqrt{I(θ)}(θ̂− θ) →^{D} N(0, 1)

or

θ̂ →^{D} N(θ, I(θ)^{−1}).

■ Requirements of this proof (i) The true value of θ is interior to the parameter space. (ii) Differentiation under the integral is valid, so that E [U(θ)] = 0 and V [U(θ)] =

ni(θ). This allows a central limit theorem to apply to U(θ). (iii) Taylor expansions are valid for the derivatives of the log-likelihood, so that higher

order terms may be neglected. (iv) A weak law of large numbers applies to J(θ).

2.7 The δ-method

We often want the asymptotic distribution of some function of a random variable. The basic method of doing this is the δ-method. Let X_{1}, X_{2}, …, X_{n}, be a random sample from a distribution with mean μ, variance σ^{2}. Then E (\overline{X} − μ) = 0, V (\overline{X}) = σ^{2}/n. By Taylor’s Theorem

g (\overline{X}) = g (μ + (\overline{X} − μ))
= g(μ) + (\overline{X} − μ)g′(μ) + \frac{1}{2}(\overline{X} − μ)^{2}g′′(δ_{n}),
0 < |δ_{n} − μ| < |\overline{X} − μ| ,

and as n becomes large

\begin{aligned}E [g (\overline{X})] & = & g(μ) + O (n^{−1}) , \\ V [g (\overline{X})] & = & n^{−1}σ^{2}g′(μ)^{2} + O (n^{− 3/2}) . \end{aligned}

Example 2.25 Exponential distribution If

f_{X}(x) = λe^{−λx}, x ≥ 0,

then

E(X) = λ^{−1}, V (X) = λ^{−2},

so that

E(\overline{X}) = λ^{−1}, V (\overline{X}) = 1/ nλ^{2}.

Suppose we want E(Y) and V (Y), where Y = log \overline{X}. Then

g(x) = log x, g′(x) = 1/ x

and

E(Y) ≃ − log λ, V (Y) ≃ \frac{1}{nλ^{2}} (\frac{1}{x})^{2}|_{x= 1/λ} = \frac{1}{n}.

■

Example 2.26 Variance stabilising transformations

The aim of a variance stabilising transformation is to find a transformation such that g (\overline{X}) has a variance which is approximately constant, i.e.

σ^{2}g′(μ)^{2} ≃ c, a constant.

Suppose we have a random sample X_{1}, X_{2}, …, X_{n} with mean μ, variance V (μ). Then

\begin{aligned}& g′(μ) & = & cV (μ)^{− 1/2} \\ ⇒ & g(μ) & = & ∫ ^{μ} V (u)^{− 1/2}du. \end{aligned}

In practice, the transformation is usually applied to the sample data directly to produce transformed data g(X_{1}), g(X_{2}), …, g(X_{n}). For a Poisson distribution, E(X) = μ, V (X) = μ so that g(μ) = μ^{1/2}. Thus, if X ∼ P oisson(μ), then

Y = \sqrt{X}

has a variance which is approximately constant.

For an exponential distribution, V (μ) = μ^{2}, so

g(μ) = log μ.

■