Lecture 12 Sampling Distributions and Standard Errors

The model of the day is $X_{1},\dots,X_{n}$ i.i.d. from the exponential distribution with rate $\theta$, where the maximum likelihood estimator is $\hat{\theta}=1/\bar{X}$ by example 5.2. The command rexp can be used to simulate such a sample, and as in lecture 6 we package one run as a function. Before you run the code below, predict the two histograms: for $n=10$, will the estimates form a symmetric bell around the true value $\theta=1$? Will the peak sit above, below or at $1$? How will the picture change for $n=100$?

As in lecture 9, the histograms are drawn with freq = FALSE, so that a density can be drawn on top of them. The command curve can be used to plot any function of x on top of the existing plot, if the argument add = TRUE is used; we use it to overlay the density. The density we overlay is a normal approximation which we can derive from appendix A: the observations have mean $1/\theta$ and variance $1/\theta^{2}$, and thus the central limit theorem (theorem A.15) gives $\sqrt{n}\,(\bar{X}-1/\theta)\xrightarrow{\ \mathrm{d}\ }N(0,1/\theta^{2})$. Applying the delta method (proposition A.18) with $g(m)=1/m$ and $g^{\prime}(1/\theta)=-\theta^{2}$ we find

so that for large $n$ the estimator $\hat{\theta}=1/\bar{X}$ is approximately normal with mean $\theta$ and standard deviation $\theta/\sqrt{n}$. The command qqnorm can be used to plot the sorted estimates against the quantiles of a standard normal distribution, and qqline can be used to add the straight line which the points would follow if the estimates were exactly normally distributed. Points above the line at the right-hand end indicate a right tail which is heavier than that of the normal distribution, and points above the line at the left-hand end indicate a left tail which is lighter.

Figure 12.1 shows that for $n=10$ the estimates are not normally distributed: the histogram has a long right tail, the peak is to the left of $1$ and the quantile plot shows quantiles diverging from the line. This is no contradiction: the central limit theorem and the delta method only promise the normal shape for large $n$. For $n=100$ the distortion is gone.

The standard deviation $\theta/\sqrt{n}$ of the normal curves in figure 12.1 can also be read off the Fisher information from lecture 11: by example 11.8 and proposition 11.4 the sample contains information $\mathcal{I}_{n}(\theta)=n/\theta^{2}$, and thus $1/\sqrt{\mathcal{I}_{n}(\theta)}=\theta/\sqrt{n}$. This is not a coincidence. Lecture 17 will show that every maximum likelihood estimator is, under conditions, approximately normal with spread $1/\sqrt{\mathcal{I}_{n}(\theta)}$ for large $n$, and that $\theta$ may be replaced by the estimate inside the information. We use these two facts on trust in this session; the delta method above has verified the first of them for the exponential rate. Since the value of $\theta$ is normally unknown, the standard error of the estimate is given by

This is a measure of the likely size of the estimation error, computed from the data. Since the estimator is approximately normally distributed with this spread, the Wald interval $\hat{\theta}\pm 1.96\,\mathop{\mathrm{se}}\nolimits(\hat{\theta})$ covers the true value with probability approximately $95\%$. This is a confidence interval of the kind we have seen in MATH2701, and is used informally here and more formally in lecture 25. Before you try the code, make a prediction: will the interval cover the true value $1$?

The interval $[0.407,0.861]$ does not contain the true value $1$. This is not a programming error: this seed happens to give one of the misses, and in the long run a $95\%$ interval misses the truth in about one sample out of twenty. From a single sample it is not possible to know whether this is one of the misses; how often the recipe is wrong in the long run can be answered by simulation.

Recall from lecture 9 that a probability is the expected value of an indicator, and thus by the law of large numbers (proposition A.14) it can be estimated by the mean of a logical vector over many independent repetitions. The same principle applies to all probabilities we estimate by simulation: for the coverage probability of an interval, the indicator is that the interval contains the truth, and for the rejection rate of a test (see lecture 24) it is that the test rejects.

We apply this to the coverage of the Wald interval: the box below contains a complete script, with its lines shuffled, except that the lines of the function are kept together. Before you continue, arrange the lines on paper into a working script which estimates, from $10000$ samples of size $n=10$, the probability that the Wald interval contains $\theta$. Several orders satisfy the rules below.

Names must be defined before they are used and the seed set before the first random draw (lecture 6): the function, theta and the seed come first, in any order, then replicate, then the summary. One detail is new: the function returns the two endpoints, which replicate collects as the columns of a two-row matrix, so that ci[1, ] and ci[2, ] hold the lower and upper endpoints, and & combines two logical vectors elementwise (see appendix B). Below, the reconstructed script is extended by a loop over $n=10$, $30$ and $100$. Before running it, make a prediction about the misses at $n=10$: the standard error $\hat{\theta}/\sqrt{n}$ is proportional to the estimate, so will an interval which misses the truth more often lie entirely below or entirely above $\theta$?

The coverage is close to the nominal $0.95$ for all three sample sizes, even for $n=10$ where the sampling distribution was far from normal. Since $10000$ repetitions estimate a probability near $0.95$ with a standard error of about $0.002$, the differences between the three values are within simulation error. The total hides an imbalance: small estimates get short intervals, and at $n=10$ most misses are intervals whose upper endpoint lies below the truth. Check this yourself, using the logical vectors ci[2, ] < theta and ci[1, ] > theta in place of the coverage indicator.

The following task is in the style of a short practical report. Choose a rate $\theta$, let $n=50$, set the random number generator seed to your student ID, and simulate one sample from an exponential distribution. Determine the estimate $\hat{\theta}=1/\bar{x}$, the standard error, and the Wald interval for this sample. Write down your results in one sentence each. Then, using $10000$ samples, estimate the coverage of the interval at $n=50$. Also determine how often the upper endpoint of the interval lies below $\theta$ and how often the lower endpoint lies above $\theta$, and report the two proportions separately. In two sentences, discuss whether you would trust the interval for this sample size.

Before the next session, type the commands used in this session into R, and check that your answers match the text.