Lecture 9 Consistency by Simulation

In this section we simulate one long sample and watch the estimate change as the observations arrive one by one.

We will consider the exponential distribution with rate $\theta=2$, where the maximum likelihood estimator for the rate is $\hat{\theta}_{n}=1/\bar{X}_{n}$, as found in example 5.2. This estimator is consistent for $\theta$ by example 2.12 and proposition A.17. The most direct way to see the estimator at work is to simulate one long sample and to re-compute the estimate after each observation. The command cumsum can be used to compute the running totals $X_{1}$, $X_{1}+X_{2}$, $X_{1}+X_{2}+X_{3},\ldots$ in one step, and dividing the count $n$ by these totals gives all $2000$ running estimates at once. Before you look at the plot, sketch the path you expect. Where is it most volatile?

Figure 9.1 shows the expected behaviour: the path fluctuates strongly for small $n$, crosses the edge of the band $\theta\pm 0.1$ several times and ends inside it. A single path, however, cannot tell us how likely it is to be outside the band at a given $n$, and this is what consistency is about.

In this section we compare the sampling distribution at three sample sizes and which fraction of estimates misses the truth by more than a fixed amount.

By definition 2.9, consistency means that for every band around the truth, here $\theta\pm 0.1$, the probability $\mathbb{P}_{\theta}\bigl{(}|\hat{\theta}_{n}-\theta|\geq 0.1\bigr{)}$ of missing the band goes to zero as $n$ increases. A probability is a property of many repetitions, and thus, as in lecture 6, we use replicate to compute $10000$ independent estimates, at each of the sample sizes $n=10$, $n=100$ and $n=1000$. The fraction of estimates outside the band estimates the probability; as in appendix B, the mean of a condition is the fraction of cases where it holds. Before you run the code, write down three guesses for the three fractions. Will the fraction for $n=1000$ be close to zero?

The probability of missing the band decreases from $0.88$ at $n=10$ to $0.12$ at $n=1000$, as consistency requires, but at $n=1000$ it is still far from zero: the band is narrow, and consistency only says that the probability goes to zero eventually. To see the whole sampling distribution we draw the three histograms with a common horizontal axis, as in lecture 6. New is the argument freq = FALSE, which draws the histogram on the density scale, so that the total area of the bars is one and a narrower histogram must be a taller one. Values outside the horizontal range are not shown; for $n=10$ these are about $8\%$ of the estimates, almost all of them in the long right tail. The dotted lines mark the band.

Figure 9.2 is the picture promised in lectures 3 and 6: the histogram gets tighter around the truth as $n$ grows. For $n=10$ it is also visibly skewed to the right, since $1/\bar{X}_{n}$ is large whenever the sample mean is small, while for $n=1000$ it is close to symmetric.

In this section we estimate the mean squared error of two estimators for a range of sample sizes, collect the results in a table and plot them.

For a report a single number per sample size is often more useful than a histogram. In lecture 6 we estimated the mean squared error of the two uniform estimators $2\bar{X}$ and $\max_{i}X_{i}$ for $n=10$. We now repeat this study for several sample sizes, which is the usual shape of a simulation study: an outer loop over a setting, an inner Monte Carlo run and a table of results. We first put both estimates into one function, so that both are computed from the same sample and returned as a vector of length two. The command replicate then collects the $10000$ repetitions into a matrix, a rectangular table of numbers with two rows, one per estimator, and $10000$ columns, one per repetition, and the command rowMeans can be used to average each row.

The outer loop runs over a vector ns of sample sizes, using seq_along(ns) to count through its positions. The command matrix(0, nrow = length(ns), ncol = 2) can be used to set aside a matrix for the results, mse[j, ] refers to the $j$-th row of this matrix, and row and column names, set using rownames and colnames, can be used to make the matrix a readable table. As in lecture 6, the box below shows the complete script, with the lines of the loop kept together for clarity, but with the lines in a shuffled order.

Exercise. Before you read the script, you need to work out the correct order of lines to make the script work.

The rule for writing R scripts is still in force: all names must be defined before they are used, so ns must come before the matrix is created and theta before the loop, the seed must be set before the loop starts, and the display (using signif to round the numbers) must come last. As in lecture 6, several orders satisfy these rules.

For $n=10$ the mean squared error of the maximum was approximately half of the one for $2\bar{X}$. Before you run the reconstructed script, make a prediction. Will this ratio stay close to two for $n=1000$?

The values in each column decrease towards zero as $n$ increases, and thus both estimators are consistent by lemma 2.11. The ratio between the two columns is not close to two anymore: for $n=1000$ the mean squared error for the maximum is smaller by a factor of approximately $170$. This is expected, since the exact mean squared error $\theta^{2}/(3n)$ of exercise 2.2 decreases like $1/n$, while the exact value $2\theta^{2}/\bigl{(}(n+1)(n+2)\bigr{)}$ of exercise 4.4 decreases like $1/n^{2}$. On log-log axes, which the argument log = "xy" selects, a quantity proportional to $n^{-k}$ is a straight line with slope $-k$. The command points can be used to add a second set of points to an existing plot, with the plotting symbol chosen by the argument pch, and the command legend can be used to add a key. Will the two sets of points be parallel?

Both sets of points lie on straight lines: the slope is $-1$ for $2\bar{X}$ and the slope is $-2$ for the maximum. The rate $1/n$ is the usual rate of decrease for a mean squared error; it is the rate of the sample mean in example 2.8. The faster rate for the maximum is a consequence of the support of the uniform distribution moving with $\theta$, the same feature which, by the remark after example 7.4 in lecture 7, keeps the uniform distribution out of the exponential families. We will consider this example again in lecture 13.

In this section we check by simulation that the conditional distribution of a Poisson sample given its total does not depend on $\theta$.

To conclude this lecture, we consider an example from lecture 8: for Poisson samples, the total $\sum_{i}X_{i}$ is sufficient by example 8.6. By definition of sufficiency, given the total, the conditional distribution of the sample does not depend on $\theta$. We can see this by simulation, without any algebra. We generate many Poisson samples of size $n=5$, keep only those whose total equals $10$, and look at the first observation in the kept samples; we do this for two different values of $\theta$.

The command replicate returns a matrix with $5$ rows and one column per sample, colSums computes the column totals, and samples[1, totals == 10] selects from the first row the columns where the condition holds. The command table can be used to count how often each value occurs, and dividing by the number of kept samples turns the counts into proportions. Before you run the code, make a prediction. Will the two tables of proportions be different, since the two values of $\theta$ are different? What will the mean of the first observation be in each case?

The total equals $10$ for about $9\%$ of the samples when $\theta=1.5$ and for about $5\%$ when $\theta=3$, so the two values of $\theta$ clearly matter for how often the total is $10$. Once the total is fixed, however, the two tables agree up to simulation error, and both means are close to $10/5=2$: given the total, the first observation has the same distribution whether $\theta=1.5$ or $\theta=3$. This is definition 8.1 in action. All that the data can tell us about $\theta$ is carried by the total, and the distribution of the individual counts around it is the same for every $\theta$. (The proportions in the table are themselves estimates; with about $9000$ and $5000$ kept samples the simulation error is between $0.005$ and $0.007$ for the larger entries.)

The following task is a review of the topics of this lecture, in the style of a practical report. Following the structure of section 9.3, for the exponential rate with $\theta=2$ and sample sizes given in ns, estimate the bias and mean squared error of the maximum likelihood estimator $1/\bar{X}$ using $10000$ repetitions each, in a table with named rows and columns. Create a plot of the mean squared error against $n$, using logarithmic axes. From the plot, determine the slope. Compare the empirical bias to the exact value from exercise 2.3. Write three sentences about what the table and plot show about the consistency of the estimator $1/\bar{X}$.

Before the next session, type in the commands from this session and check that the output matches what is shown in the text.