Practical — Estimation Theory (STS-205, Section A)

Topics Covered

Section A — The Five Practicals

Practical 1 — Jackknife Estimates

1. Problem

For the sample \(3, 5, 7, 11, 14\), the variance is estimated by the divisor-\(n\) (biased) variance \(\hat\theta = \frac{1}{n}\sum_i (x_i - \bar x)^{2}\). Compute the jackknife estimate of its bias, the bias-corrected jackknife estimate, and the jackknife standard error.

2. Aim

To reduce the bias of an estimator and estimate its standard error by the jackknife, which leaves out one observation at a time.

3. Formula

\(\hat\theta_{(-i)}\) is the same statistic computed with \(x_i\) omitted. The pseudo-values \(\tilde\theta_i\) average to \(\hat\theta_{\text{jack}}\), and their sample variance gives the standard error. If the bias has the form \(a_1/n + a_2/n^{2} + \cdots\), the construction cancels the \(a_1/n\) term exactly; nothing about the distribution is assumed (Unit 2, section 4).

Applying it:

1. Compute \(\hat\theta\) on the full sample.
2. Compute the \(n\) leave-one-out values \(\hat\theta_{(-i)}\) and their mean.
3. Estimate the bias and correct for it.
4. Form the pseudo-values; the standard error comes from them, not from the leave-one-out values directly.

4. Calculation

\(\bar x = 40/5 = 8\), and \(\sum_i (x_i - 8)^{2} = 25 + 9 + 1 + 9 + 36 = 80\), so \(\hat\theta = 80/5 = 16\).

The unbiased variance is \(s^{2} = 80/(5-1) = 20\). The pseudo-values have mean \(20\), as they must, and their sample variance divided by \(n\) gives \(\widehat{SE}_{\text{jack}} = 7.925434\).

x <- c(3, 5, 7, 11, 14); n <- length(x)
biased_var <- function(v) mean((v - mean(v))^2)
th <- biased_var(x)
loo <- sapply(seq_len(n), function(i) biased_var(x[-i]))
ps  <- n * th - (n - 1) * loo            # pseudo-values
loo; ps
c(estimate = mean(ps), se = sd(ps) / sqrt(n))
var(x)                                    # 20 -- the jackknife recovered it

[1] 12.1875 17.1875 19.6875 17.1875  8.7500
[1] 31.25 11.25  1.25 11.25 45.00
 estimate        se 
20.000000  7.925434 
[1] 20

5. Result

\(\widehat{\text{bias}} = -4\), \(\hat\theta_{\text{jack}} = 20\) and \(\widehat{SE}_{\text{jack}} = 7.93\). The jackknife reproduced the unbiased variance \(s^{2} = 20\) exactly, having been told nothing about the correct divisor. This is not a coincidence — the bias of the divisor-\(n\) variance is exactly \(-\sigma^{2}/n\), of the form the jackknife is built to remove. The same four lines apply to a statistic whose bias nobody has worked out, and give a standard error where no formula exists.

Practical 2 — Bootstrap Estimates

1. Problem

For the sample \(2, 5, 11\) and the statistic \(\hat\theta = \bar x\), compute the bootstrap estimates of bias and standard error by enumerating every resample.

2. Aim

To estimate the bias and standard error of a statistic from its bootstrap distribution, computed here exactly by complete enumeration.

3. Formula

Resampling from the data, with replacement, is resampling from the empirical distribution function \(F_n\), justified by the Glivenko–Cantelli lemma (Unit 2, section 5). With \(n = 3\) there are exactly \(3^{3} = 27\) resamples, so the bootstrap distribution can be written down in full and no randomness enters; the answer is exact and reproducible in an examination.

Applying it:

1. Compute \(\hat\theta\) on the original sample.
2. List all \(n^{n}\) resamples and compute \(\hat\theta^{*}\) on each.
3. The bias is the mean of the \(\hat\theta^{*}\) less \(\hat\theta\); the standard error is the square root of their variance (divisor: the number of resamples).

4. Calculation

\(\bar x = (2 + 5 + 11)/3 = 6\). Over all 27 resamples, the mean of the bootstrap means is \(\bar{\hat\theta}^{*} = 6.000000\), exactly \(\bar x\), and

The biased sample variance is \(s_n^{2} = \{(2-6)^{2} + (5-6)^{2} + (11-6)^{2}\}/3 = (16 + 1 + 25)/3 = 14\), and \(s_n^{2}/n = 14/3 = 4.666667\) — exactly the bootstrap variance. The classical \(\sqrt{s^{2}/n} = \sqrt{21/3} = 2.645751\).

s <- c(2, 5, 11)
g <- expand.grid(s, s, s)                # all 27 resamples
bm <- rowMeans(g)
c(n_resamples = nrow(g), mean = mean(bm), bias = mean(bm) - mean(s),
  var = mean((bm - mean(bm))^2), se = sqrt(mean((bm - mean(bm))^2)))

n_resamples        mean        bias         var          se 
  27.000000    6.000000    0.000000    4.666667    2.160247 

5. Result

The bootstrap bias is \(0\) — correctly, since the sample mean is unbiased — and \(\widehat{SE}_{\text{boot}} = 2.160\), against the classical 2.646. The two differ by exactly \(\sqrt{(n-1)/n} = \sqrt{2/3} = 0.8165\), negligible for realistic \(n\): for the mean the bootstrap reproduces a known formula exactly, and for a statistic with no known formula the same procedure runs unchanged.

A practical warning. \(n^{n}\) resamples is \(27\) at \(n = 3\) but \(10^{10}\) at \(n = 10\). Enumeration is an examination device for small \(n\); for real data the resamples are drawn at random, and \(B = 1000\) or more is used with the seed recorded.

Practical 3 — MLE by the Method of Scoring, Cauchy Location

1. Problem

A sample of \(n = 7\) from Cauchy\((\theta, 1)\) is

Find the maximum likelihood estimate of \(\theta\) by the method of scoring, its standard error, and an approximate 95% confidence interval.

2. Aim

To solve a likelihood equation that has no closed form, by Fisher's method of scoring.

3. Formula

For the Cauchy there is no closed-form MLE and no useful moment estimator — Distribution Theory, Unit 1 proves the mean does not exist, so the sample mean estimates nothing. The information is \(\tfrac12\) per observation, a constant — which is what makes scoring simpler than Newton–Raphson here: the divisor never has to be recomputed.

Applying it:

1. Start from the sample median, \(0.72\), the fourth of seven ordered values. Never start from the sample mean \(0.840000\): for a Cauchy it is not consistent for anything, and a poor start can send the iteration to a local maximum, of which the Cauchy likelihood can have several.
2. Compute the score at the current value and add \(U/3.5\).
3. Repeat until the score is effectively zero.
4. Give the standard error from the information, and the interval \(\hat\theta \pm 1.96\,\widehat{SE}\).

4. Calculation

Each step is about \(0.26\) times the previous one, so the error falls by roughly a quarter per iteration — linear convergence, which is what scoring gives (Newton–Raphson would be quadratic but needs the observed information recomputed each time). After five iterations \(\hat\theta = 0.666520\) and the score is \(-1.66 \times 10^{-4}\), effectively zero.

y <- c(-1.94, -0.44, 0.10, 0.72, 1.35, 2.18, 3.91)
score <- function(t) sum(2 * (y - t) / (1 + (y - t)^2))
th <- median(y)                           # NOT mean(y)
for (k in 1:5) th <- th + score(th) / (length(y) / 2)
c(mle = th, score = score(th), se = sqrt(2 / length(y)))
th + c(-1, 1) * 1.96 * sqrt(2 / length(y))

          mle         score            se 
 0.6665195849 -0.0001647852  0.5345224838 
[1] -0.3811445  1.7141837

5. Result

\(\hat\theta = 0.6665\), with standard error 0.5345 and approximate 95% interval \((-0.381,\ 1.714)\).

Interpretation. The MLE \(0.6665\) sits below the sample mean \(0.8400\), which the single large observation \(3.91\) has dragged upwards; the likelihood downweights it automatically, because a Cauchy expects occasional extreme values. The interval is very wide for seven observations — \(\sqrt{2/n}\) falls only as \(n^{-1/2}\), and the Cauchy is an uninformative distribution to sample from.

Practical 4 — Confidence Limits for a Normal Population

1. Problem

A sample of \(n = 16\) from \(N(\mu, \sigma^{2})\), with both parameters unknown, has \(\bar x = 24.5\) and \(s = 3.2\). Find 95% confidence limits for \(\mu\), \(\sigma^{2}\) and \(\sigma\).

2. Aim

To find exact confidence limits for the mean and the variance of a normal population from the t and chi-square pivots.

3. Formula

Applying it:

1. Compute the standard error \(s/\sqrt n\).
2. Read \(t_{15,\,0.025}\) and form the interval for \(\mu\).
3. Read the two 2.5% points of \(\chi^{2}_{15}\) and form the interval for \(\sigma^{2}\).
4. Take square roots for \(\sigma\), which is legitimate because the transformation is monotone.

4. Calculation

\(t_{15,\,0.025} = 2.131449\) (check: \(P(|t_{15}| > 2.131449) = 0.050000\)), so

The 2.5% points of \(\chi^{2}_{15}\) are \(6.262\) and \(27.488\), and \((n-1)s^{2} = 15 \times 10.24 = 153.6\), so

nn <- 16; xb <- 24.5; s4 <- 3.2
round(xb + c(-1, 1) * qt(0.975, nn - 1) * s4 / sqrt(nn), 4)       # mu
round((nn - 1) * s4^2 / qchisq(c(0.975, 0.025), nn - 1), 4)       # sigma^2
round(sqrt((nn - 1) * s4^2 / qchisq(c(0.975, 0.025), nn - 1)), 4) # sigma

[1] 22.7948 26.2052
[1]  5.5878 24.5284
[1] 2.3639 4.9526

5. Result

With 95% confidence, \(22.79 < \mu < 26.21\), \(5.59 < \sigma^{2} < 24.53\) and \(2.36 < \sigma < 4.95\).

Interpretation. The interval for \(\mu\) is symmetric and narrow; the one for \(\sigma\) is neither, running from \(0.74\) to \(1.55\) times the estimate \(s = 3.2\). Sixteen observations locate a mean well and a standard deviation poorly, and any report that quotes \(s\) without an interval is hiding that.

Practical 5 — Large-Sample Confidence Limits

1. Problem

Find approximate 95% confidence limits for: (a) a binomial \(p\), with \(x = 88\) successes in \(n = 400\) trials; (b) a Poisson \(\lambda\), with \(132\) events observed over \(50\) units of exposure; (c) an exponential rate \(\theta\), with \(n = 100\) and \(\bar x = 4.8\).

2. Aim

To find large-sample confidence limits from the asymptotic normality of the maximum likelihood estimator.

3. Formula

All three cases use the asymptotic normality of the MLE from Unit 2, with the information evaluated at the estimate because \(\theta\) is unknown.

Applying it:

1. Write down \(I(\theta)\) and the MLE.
2. Substitute \(\hat\theta\) and take the reciprocal square root of \(nI\).
3. Form \(\hat\theta \pm 1.96\,\widehat{SE}\).

4. Calculation

(a) Binomial. \(\hat p = 88/400 = 0.22\); \(\widehat{SE} = \sqrt{\hat p(1 - \hat p)/n} = \sqrt{0.22 \times 0.78/400}\). The numerator is \(0.1716\), so \(\widehat{SE} = \sqrt{0.000429} = 0.020712\) and

(b) Poisson.

(c) Exponential.

A check worth doing, and what it reveals. The mean \(1/\theta\) can be handled either way. Directly, \(\widehat{SE}(\bar X) = \bar x/\sqrt n = 0.48\), giving \(4.8 \pm 1.96 \times 0.48 = (3.859200,\ 5.740800)\). Inverting that interval and reversing the endpoints gives \((0.174192,\ 0.259121)\) for \(\theta\) — not the same as the interval in (c), which was \((0.167500,\ 0.249167)\).

ph <- 88/400;  ph + c(-1, 1) * 1.96 * sqrt(ph * (1 - ph) / 400)   # (a)
lh <- 132/50;  lh + c(-1, 1) * 1.96 * sqrt(lh / 50)             # (b)
tb <- 1/4.8;   tb + c(-1, 1) * 1.96 * tb / sqrt(100)            # (c)
rev(1 / (4.8 + c(-1, 1) * 1.96 * 0.48))                         # (c) on the 1/theta scale

[1] 0.1794039 0.2605961
[1] 2.189626 3.090374
[1] 0.1675000 0.2491667
[1] 0.1741918 0.2591211

5. Result

95% limits: (a) \(0.179 < p < 0.261\); (b) \(2.19 < \lambda < 3.09\); (c) \(0.1675 < \theta < 0.2492\).

Interpretation. The discrepancy found in (c) is real and instructive. A confidence interval is not invariant under a non-linear transformation: building it on the \(\theta\) scale and building it on the \(1/\theta\) scale give different answers, because the normal approximation is being made at different places. Both are valid to the same asymptotic order, and they agree as \(n\) grows — here they overlap over most of their length. The point estimate is invariant, by property (M1) of the MLE in Unit 2; the interval is not. Which scale to work on is decided by where the normal approximation is better, and for a positive parameter that is usually the log scale.

What the Practical Record Should Contain