Lecture 9 Consistency by Simulation

Consistency, definition 2.8, is a statement about the sequence of estimators $\hat{\theta}_{1},\hat{\theta}_{2},\dots$ as the sample size increases. The most direct way to visualise this sequence is to simulate a long sample and to re-compute the estimate after each observation. 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, and we will see that this estimator is consistent for $\theta$, by theorems 2.9 and A.17.

The command cumsum can be used to compute the running totals $X_{1}$, $X_{1}+X_{2}$ and so on in one step, and dividing the count $n$ by these totals gives all $2000$ running estimates at once. Before you look at the output, predict the five printed numbers. How far from $2$ do you expect the estimate after one observation to be? How far after $2000$?

The first estimate is $1/x_{1}$, which could be any number; in this example it is more than twice the truth. By $n=1000$ the estimate is within $0.1$ of the truth. To get a better understanding of the path taken by the estimate, we plot the estimate against $n$. The argument log = "x" for the plot command tells R to plot $n$ on a logarithmic axis, so that we can still see the early part of the path. Before we look at the plot, try to sketch the path you expect the estimate to take. Where is the curve most volatile? Once the curve enters the band $\theta\pm 0.1$, does it stay there?

Figure 9.1 shows the expected behaviour: the path is wild for small $n$ and ends inside the band. It does not converge in the sense of a sequence of numbers where, after some point, no term is outside the band; the path leaves and re-enters the band at random times. What decreases is the probability of being outside the band at a given $n$, as given in definition 2.8.

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$. To see how each decreases with $n$ we repeat the study for several sample sizes, and this is the everyday shape of a simulation study in the practical report: 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 collects the $10000$ repetitions into a matrix, a rectangular table of numbers with two rows, one per estimator, 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. Before you read the script, you should try 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, i.e. theta and ns must come first, the seed of the random number generator must be set before the loop starts, and the names and the display (using the signif function to round the numbers) require the matrix to be complete. Before you try the reconstructed script, try to predict the last row of the output. What will the two mean squared errors be for $n=1000$? For $n=10$ the second mean squared error was approximately half of the first one. Will the ratio between the two mean squared errors stay close to two, or will it change?

The row for $n=10$ coincides with the two numbers from lecture 6 up to simulation error (the samples are not the same, since the loop has already sampled the $n=5$ values first), and the values for each column decrease as $n$ increases, as expected from the sufficient condition for consistency given in corollary 2.11. The ratio between the mean squared errors is not close to two anymore: for $n=1000$ the mean squared error for the maximum is smaller by a factor of approximately $170$. The exact values from exercises 2.2 and 4.4 show that this is expected, since $\theta^{2}/(3n)$ decreases like $1/n$ while $2\theta^{2}/\bigl{(}(n+1)(n+2)\bigr{)}$ decreases like $1/n^{2}$. If we plot the data on log-log axes, i.e. if we use the argument log = "xy" for the plot, then a quantity proportional to $n^{-k}$ can be plotted as a straight line with slope $-k$. What will the two sets of points look like on a log-log plot? 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 usual rate of decrease for a mean squared error is like $1/n$. All the estimators we will consider in the R strand of this module will have slope $-1$. The faster rate for the maximum is a consequence of the support of the uniform distribution moving with $\theta$. We will consider this example again in lecture 13.

We can also consider the full sampling distribution of the maximum at $n=50$, since in lecture 8 we know the exact density: by proposition 8.9 and example 8.10, the maximum of $n$ observations from $\text{Uniform}(0,1)$ has density $f(y)=ny^{n-1}$ for $0\leq y\leq 1$. To plot this density on top of a histogram, we need to plot the histogram using the density scale. This can be achieved by using the argument freq = FALSE of hist. Then the command curve can be used to plot any function of x on the horizontal axis, on top of the existing plot, if add = TRUE is used. Sketch the density $50y^{49}$ before the plot appears. Where is the density largest? What is the density at $y=1$?

The curve closely follows the tops of the bars, and reaches the value $50$ at the boundary of the interval. The histogram is not a bell-shaped distribution as in figure 6.1: a sampling distribution can be concentrated around the truth, without being normal. In figure B.1 in appendix B you can see the same construction for a normal sample.

To conclude these notes, we consider an example for use in lecture 8: for Poisson samples, the total number of counts $\sum_{i}X_{i}$ is sufficient by example 8.4. The factorisation theorem 8.2 states that for data sets with the same total count, the likelihoods are proportional, differing by a constant shift on the log scale. We consider three data sets of size $n=5$, where the first two have the same total count. We compute the log-likelihood for these data sets using the dpois function in R, with the log = TRUE option to get the log-likelihoods, and without any manual simplification of the results. The sapply command can be used to evaluate a function on every element of a vector; here we use it to evaluate the log-likelihoods on a grid of $\theta$-values. Will the difference between the first two curves be the same at every grid point? Will the same be true for the first and third curves? Do the three maxima coincide?

The difference between the first two curves is constant at every grid point, taking the value $-2.197$. This is equal to $\log(32/288)$, the logarithm of the ratio of the two factors $1/\prod_{i}x_{i}!$, which do not depend on $\theta$. The difference to the third curve varies between $-9.1$ and $3.3$. The first two maxima are both at $\bar{x}=2$ while the third maximum is at $0.8$. The likelihood does not allow us to infer anything about $\theta$ beyond what is already known from the total count: for x1 and x2, the information about $\theta$ is the same, once the total count is known.

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.2, 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}$.

Just like it is important to understand the syntax of a programming language, it is also important to be able to read the output of R commands. Before the next session, type in the commands from this session, and check that the output matches what is shown in the text.