Week 3 — Random samples, order statistics, and exponential-family structure

Where this week starts

Week 2 gave you three routes to the distribution of something computed from a random variable: the distribution-function route, the change-of-variables route with its Jacobian, and the moment-generating-function route. Each acted on one variable at a time. This week the object being transformed is the whole sample. You have nn independent draws X1,,XnX_1, \dots, X_n from one density f(xθ)f(x \mid \theta) — a random sample — and the question is what that collection knows that no single draw does.

Two answers pull in opposite directions. Sorting is itself a transformation of the whole sample, and the sorted values, the order statistics, have a joint distribution you can write down exactly. The smallest and largest observation have different distributions, and they report on edges of the population the sample mean throws away. The second answer runs the other way: for a large and familiar class of models the joint density of all nn observations depends on the parameter only through a handful of sums, and that handful does not grow with nn. Those are the exponential families.

The hinge is what Week 2 left unfinished: change of variables required a one-to-one map, and sorting is not one-to-one, since all n!n! orderings of the same numbers give the same sorted vector.

By Thursday you should be able to write the joint density of the order statistics, derive the marginal density of the kk-th smallest two independent ways, carry the uniform case through, and read moments off a cumulant function. Three things are deferred: what “all the information about θ\theta” means is Week 10, the limit on any estimator is Week 12, and mean squared error is Week 6.

Why this matters downstream

Picture a reliability group estimating the ceiling of a timing mechanism: the wait before a device reports is modelled as uniform on (0,θ)(0, \theta). Someone takes the largest observed wait, quotes a standard error, and reports a symmetric interval around it. Every part of that is wrong in a way this week makes visible. The maximum is below θ\theta with probability one, so it is biased downward by a fixed fraction rather than by luck, and its sampling distribution, rescaled correctly, converges to a one-sided exponential law. The symmetric interval puts roughly half its mass on values the data have already excluded.

The exponential-family half carries a different stake, and it is worth separating two things that get bundled together under the word “regularity”. One is structural. For a full-rank exponential family the natural statistic iT(Xi)\sum_i T(X_i) is automatically minimal sufficient and complete, so the reduction theorems of Weeks 10 and 11 arrive free the moment you recognise the form. That is a sufficient condition and not a necessary one: Uniform(0,θ)(0, \theta) sits outside the class and still has a complete sufficient statistic, namely X(n)X_{(n)}, which is exactly what will certify the bias-corrected maximum as the best unbiased estimator of θ\theta. The other thing is analytic — differentiating under the integral sign, which the score identity of Week 7 and the information bound of Week 12 both rest on, and which the uniform genuinely breaks because its region of integration moves with θ\theta. Confusing the two costs you either a theorem you were entitled to use or a bound you were never entitled to quote.

What you will be able to do

  • Write the joint density of the order statistics of a continuous sample, and say where n!n! comes from.
  • Derive the marginal density of X(k)X_{(k)} by counting, then check it against the distribution-function route at k=1k = 1 and k=nk = n.
  • Compute the exact distribution, mean, and variance of a uniform sample’s largest observation, and state the size and sign of its bias.
  • Recognise a kk-parameter exponential family, exhibit hh, TT, η\eta and AA, and describe the natural parameter space.
  • Differentiate a cumulant function for the mean and covariance of the natural statistic, verifying each against a moment you know.
  • Explain why Uniform(0,θ)(0, \theta) is outside the class, name a later theorem that therefore does not apply to it, and name one that still does.

Words worth owning

Term What it means in this course
Random sample X1,,XnX_1, \dots, X_n independent, each with density f(xθ)f(x \mid \theta).
Order statistics The same values sorted upward, X(1)X(n)X_{(1)} \le \cdots \le X_{(n)}.
Sample extremes The minimum and maximum, governed by the tails of FF, not its centre.
Support Where the density is strictly positive. Whether it moves with θ\theta is this week’s decisive question.
Exponential family A family factoring as h(x)h(x) times an exponential linear in θ\theta-free statistics, on a support fixed in θ\theta.
Natural parameter The coefficient ηj\eta_j multiplying Tj(x)T_j(x), taken as the parameter itself. Also called canonical.
Cumulant function The normalising term A(η)A(\eta). Its first derivative is a mean, its second a variance.
Natural statistic The vector T(x)T(x) in the exponent; for a sample, iT(xi)\sum_i T(x_i), whose dimension does not grow with nn.

What sorting does to a random sample

Throughout, X1,,XnX_1, \dots, X_n are independent with distribution function FF and continuous density ff. Continuity buys one convenience: ties have probability zero, so the sorted values are strictly increasing almost surely. Everything below needs care for discrete data.

From the sample to its sorted values

Sorting sends the observations to the sorted vector, X(1)X_{(1)} smallest through X(n)X_{(n)} largest. The picture shows one sample of eight going through it.

Eight labelled observations in a top row joined by crossing arrows to the same eight numbers in increasing order below; the largest value, 9.2, is highlighted as it moves to the final position.

Sorting a sample of eight, each recorded value joined by an arrow to its position in the sorted row.

Two features do the work. No number is created or destroyed, so the sorted vector determines the multiset exactly. But the arrows are lost: nothing afterwards records that 9.29.2 came from the fourth run. That is the course’s first data reduction, and it rests on an assumption — if the observations really are i.i.d. the labels carried nothing about θ\theta, and if the fourth run followed a recalibration, sorting destroyed the evidence.

For the joint density, split n\mathbb{R}^n into the n!n! open regions on which the coordinates are in a fixed strict order. On each, sorting is a permutation of coordinates: one-to-one, with Jacobian determinant of absolute value one. Week 2’s formula applies region by region, each contributing the same product of densities, so

fX(1),,X(n)(y1,,yn)=n!i=1nf(yi),y1<y2<<yn,f_{X_{(1)}, \dots, X_{(n)}}(y_1, \dots, y_n) = n! \prod_{i=1}^{n} f(y_i), \qquad y_1 \lt y_2 \lt \cdots \lt y_n,

and zero elsewhere. Audit it: by symmetry among the orderings, the region {y1<<yn}\{y_1 \lt \cdots \lt y_n\} carries probability 1/n!1/n! under the original product density, so the factor n!n! restores total mass one. Notice what the support says, too. This density lives on that ordered region rather than on a product of intervals, and a density supported on a non-product set cannot describe independent coordinates. The order statistics are dependent by construction.

The marginal law of one order statistic

The fastest correct route to a single X(k)X_{(k)} is counting. For X(k)X_{(k)} to fall in [y,y+dy)[y, y + dy), exactly one observation lands there, exactly k1k - 1 fall below, and nkn - k fall above. The number of ways to assign labelled observations to those roles is n!/{(k1)!1!(nk)!}n!/\{(k-1)!\,1!\,(n-k)!\}, and each assignment has probability F(y)k1f(y)dy{1F(y)}nkF(y)^{k-1} f(y)\,dy\,\{1 - F(y)\}^{n-k}, so

f(k)(y)=n!(k1)!(nk)!F(y)k1{1F(y)}nkf(y).f_{(k)}(y) = \frac{n!}{(k-1)!\,(n-k)!}\, F(y)^{k-1} \{1 - F(y)\}^{n-k} f(y).

The second route never mentions infinitesimals. The event {X(k)y}\{X_{(k)} \le y\} says at least kk of the observations are at most yy, and that count is Binomial(n,F(y))(n, F(y)), so

F(k)(y)=j=kn(nj)F(y)j{1F(y)}nj.F_{(k)}(y) = \sum_{j=k}^{n} \binom{n}{j} F(y)^{j} \{1 - F(y)\}^{n-j}.

Differentiating in yy makes almost every term cancel against its neighbour, and what survives is the density above. Two independent derivations agreeing is the corroboration this course asks for: a slip in one is unlikely to be repeated by the other.

The extremes are worth rebuilding rather than memorising. All observations are at most yy exactly when the largest is, so F(n)(y)=F(y)nF_{(n)}(y) = F(y)^n; all exceed yy exactly when the smallest does, so F(1)(y)=1{1F(y)}nF_{(1)}(y) = 1 - \{1 - F(y)\}^{n}. Differentiating gives nF(y)n1f(y)n F(y)^{n-1} f(y) and n{1F(y)}n1f(y)n \{1 - F(y)\}^{n-1} f(y), matching the general formula at k=nk = n and k=1k = 1.

The uniform sample, done completely

Take X1,,XnX_1, \dots, X_n uniform on (0,θ)(0, \theta), so F(y)=y/θF(y) = y/\theta and f(y)=1/θf(y) = 1/\theta on [0,θ][0, \theta]. Substituting and cancelling powers of θ\theta,

f(k)(y)=n!(k1)!(nk)!yk1(θy)nkθn,0yθ.f_{(k)}(y) = \frac{n!}{(k-1)!\,(n-k)!} \cdot \frac{y^{k-1} (\theta - y)^{n-k}}{\theta^{n}}, \qquad 0 \le y \le \theta.

With θ=1\theta = 1 this is the Beta(k,nk+1)(k, n-k+1) density, so every uniform order statistic is a Beta variable and the Beta moments apply:

𝔼[X(k)]=kn+1θ,Var(X(k))=k(nk+1)(n+1)2(n+2)θ2.\mathbb{E}[X_{(k)}] = \frac{k}{n+1}\,\theta, \qquad \operatorname{Var}(X_{(k)}) = \frac{k\,(n-k+1)}{(n+1)^2 (n+2)}\,\theta^{2}.

The mean formula says something concrete: the nn order statistics cut the interval into n+1n+1 pieces of equal expected length.

Three density curves on the unit interval for a sample of ten: the smallest draw piles up against zero, the fifth is a broad hump near 0.45, and the largest piles up against one; triangles mark the three means.

Marginal densities of the smallest, the fifth, and the largest of ten uniform draws, with their means marked.

The figure makes three points formulas tend not to. The curves have different shapes, so the order statistics are not identically distributed. The extremes are far more concentrated than the middle one — the variance formula gives 10/(12112)=0.006910/(121 \cdot 12) = 0.0069 for X(10)X_{(10)} against 30/(12112)=0.020730/(121 \cdot 12) = 0.0207 for X(5)X_{(5)}, a visible factor of three. And the extreme densities are sharply skewed against a hard wall, which is already why a symmetric error bar around a sample maximum is the wrong shape.

Dependence can be quantified. For iji \le j the covariance of two uniform order statistics on (0,1)(0,1) is i(nj+1)/{(n+1)2(n+2)}i(n-j+1)/\{(n+1)^2 (n+2)\}, strictly positive for every pair. Check the smallest case: with n=2n = 2 the formula gives 1/361/36, while directly 𝔼[X(1)X(2)]=𝔼[X1X2]=1/4\mathbb{E}[X_{(1)} X_{(2)}] = \mathbb{E}[X_1 X_2] = 1/4, since the product of the sorted values is the product of the original ones, and 1/4(1/3)(2/3)=1/361/4 - (1/3)(2/3) = 1/36. The sign makes sense: a sample landing high pushes every order statistic up together.

One guard against over-generalising, stated with the hypothesis that actually makes it true. Central order statistics do behave conventionally, but the index has to approach its target level fast enough, and that condition is the part students drop. Let ξp\xi_p be the population quantile, so that F(ξp)=pF(\xi_p) = p with 0<p<10 \lt p \lt 1, and suppose ff is positive and continuous at ξp\xi_p. If the indices satisfy the rate condition n(kn/np)0\sqrt{n}\,(k_n/n - p) \to 0 — of which kn=npk_n = \lceil np \rceil is the ordinary case — then

n(X(kn)ξp)dN(0,p(1p)f(ξp)2).\sqrt{n}\,\big(X_{(k_n)} - \xi_p\big) \xrightarrow{d} N\Big(0, \; \frac{p(1-p)}{f(\xi_p)^{2}}\Big).

The weaker requirement kn/npk_n/n \to p will not do, and what fails is the centring rather than the spread. Take kn=np+n3/4k_n = \lceil np + n^{3/4} \rceil. Then kn/npk_n/n \to p, so the weak hypothesis holds, but X(kn)X_{(k_n)} concentrates near the quantile at level kn/nk_n/n, which sits roughly n1/4/f(ξp)n^{-1/4}/f(\xi_p) above ξp\xi_p; multiplying by n\sqrt{n} leaves a centre growing like n1/4/f(ξp)n^{1/4}/f(\xi_p), and the sequence converges to nothing at all. The general statement keeps the variance and tracks the drift: if n(kn/np)c\sqrt{n}\,(k_n/n - p) \to c for a finite cc, the limit is N(c/f(ξp),p(1p)/f(ξp)2)N\big(c/f(\xi_p),\, p(1-p)/f(\xi_p)^{2}\big), and c=0c = 0 recovers the clean version above. The sample median is the case to keep in mind — kn=n/2k_n = \lceil n/2 \rceil and p=1/2p = 1/2 satisfy the rate condition, and the limiting variance is 1/{4f(ξ1/2)2}1/\{4 f(\xi_{1/2})^{2}\}: a density that is flat where it is cut gives an imprecise median, because neighbouring values of the quantile are then nearly as plausible. Only the extremes, where the count of observations beyond the point of interest stays bounded instead of growing with nn, need a theory of their own.

Exponential families and what they buy you

The second half of the week asks a different question: what shape makes a whole sample easy to handle at once?

The k-parameter form, written carefully

A family {f(xθ):θΘ}\{f(x \mid \theta) : \theta \in \Theta\} is a kk-parameter exponential family when

f(xθ)=h(x)exp{j=1kηj(θ)Tj(x)A(θ)},x𝒳,f(x \mid \theta) = h(x) \exp\Big\{ \sum_{j=1}^{k} \eta_j(\theta)\, T_j(x) - A(\theta) \Big\}, \qquad x \in \mathcal{X},

where h0h \ge 0 and each TjT_j depend on the data alone, each ηj\eta_j on the parameter alone, and — the condition that does the damage later — the set 𝒳\mathcal{X} where the density is positive is the same for every θ\theta. The term A(θ)A(\theta) makes the density integrate to one.

Reparameterise so the coefficients are the parameters. With η=(η1,,ηk)\eta = (\eta_1, \dots, \eta_k),

f(xη)=h(x)exp{ηT(x)A(η)},A(η)=log𝒳h(x)eηT(x)dx,f(x \mid \eta) = h(x) \exp\{ \eta^{\top} T(x) - A(\eta)\}, \qquad A(\eta) = \log \int_{\mathcal{X}} h(x)\, e^{\eta^{\top} T(x)}\, dx,

with a sum replacing the integral in the discrete case. The natural parameter space \mathcal{H} is the set of η\eta making that integral finite; it is convex, which matters because the differentiation below needs an interior point.

For nn independent draws the joint density is a product, and a product of exponentials adds exponents:

fn(x1,,xnη)={i=1nh(xi)}exp{j=1kηji=1nTj(xi)nA(η)}.f_n(x_1, \dots, x_n \mid \eta) = \Big\{ \prod_{i=1}^{n} h(x_i) \Big\} \exp\Big\{ \sum_{j=1}^{k} \eta_j \sum_{i=1}^{n} T_j(x_i) - n A(\eta) \Big\}.

The sample lies in an exponential family of the same order kk, with natural statistic iT(xi)\sum_i T(x_i) and cumulant function nA(η)nA(\eta). The parameter touches the data through kk numbers however large nn becomes: a thousand Poisson observations reach the likelihood through their sum and nothing else. A fixed-dimensional summary of arbitrarily much data is exactly the kind of structure a sample can have and a single observation cannot exhibit.

Differentiating the cumulant function

Take k=1k = 1. By construction eA(η)=h(x)eηT(x)dxe^{A(\eta)} = \int h(x)\, e^{\eta T(x)}\, dx. Differentiate both sides, moving the derivative inside the integral:

A(η)eA(η)=T(x)h(x)eηT(x)dx.A'(\eta)\, e^{A(\eta)} = \int T(x)\, h(x)\, e^{\eta T(x)}\, dx.

Divide by eA(η)e^{A(\eta)} and recognise the density on the right:

A(η)=T(x)h(x)eηT(x)A(η)dx=𝔼η[T(X)].A'(\eta) = \int T(x)\, h(x)\, e^{\eta T(x) - A(\eta)}\, dx = \mathbb{E}_\eta[T(X)].

Differentiate once more. The derivative of eηT(x)A(η)e^{\eta T(x) - A(\eta)} in η\eta is {T(x)A(η)}\{T(x) - A'(\eta)\} times itself, so

A(η)={T(x)A(η)}T(x)h(x)eηT(x)A(η)dx=𝔼η[T2]{𝔼η[T]}2=Varη(T(X)).A''(\eta) = \int \{T(x) - A'(\eta)\}\, T(x)\, h(x) e^{\eta T(x) - A(\eta)}\, dx = \mathbb{E}_\eta[T^2] - \{\mathbb{E}_\eta[T]\}^2 = \operatorname{Var}_\eta(T(X)).

In kk dimensions the same computation makes the gradient of AA the mean vector of TT and the Hessian its covariance matrix. Since a covariance matrix is positive semidefinite, AA is convex.

NoteWhere the interchange is licensed, and where it is not

Moving /η\partial/\partial\eta inside the integral is the step to be suspicious of. It is legal here because at an interior point of \mathcal{H} the integrand is smooth in η\eta and locally dominated by an integrable function, and — crucially — the region of integration does not move when η\eta moves. Week 7’s score identity 𝔼θ[/θ]=0\mathbb{E}_\theta[\partial \ell / \partial \theta] = 0, in which \ell is the log-likelihood, is the same interchange in other clothes, and the uniform breaks it precisely because its region of integration does move.

It is worth seeing A=𝔼[T]A' = \mathbb{E}[T] rather than only reading it, on the Bernoulli cumulant function A(η)=log(1+eη)A(\eta) = \log(1 + e^{\eta}).

An increasing, gently curving function of the natural parameter eta, with dashed tangent lines at eta equal to minus 1.5 and eta equal to 0.8; the tangent slopes are labelled 0.18 and 0.69, the matching success probabilities.

The Bernoulli cumulant function with tangent lines at two natural parameter values, slopes labelled as means.

At η=0.8\eta = 0.8 the tangent slope is e0.8/(1+e0.8)=0.690e^{0.8}/(1 + e^{0.8}) = 0.690, the success probability there, and the curve’s bend is 0.690×0.310=0.2140.690 \times 0.310 = 0.214, the variance. At η=1.5\eta = -1.5 the slope is 0.1820.182 and the bend 0.1490.149. Two structural facts are visible without algebra: the curve never bends downward, since a variance is never negative, and far out in either direction it straightens as pp saturates.

Say what that straightening means carefully, because the loose version of the sentence is false. The bend is A(η)=p(1p)A''(\eta) = p(1-p), which is the Fisher information about the natural parameter η\eta, and it does collapse toward zero as pp approaches 00 or 11: one observation then moves your view of η\eta hardly at all. Notice, though, that there is no edge being approached. Saturation means η\eta \to -\infty or η+\eta \to +\infty, and the natural parameter space here is the whole line, so “near the boundary of the parameter space” is not even a description of what is happening in the picture. Information about pp runs the other way. A one-line computation from logf(xp)=xlogp+(1x)log(1p)\log f(x \mid p) = x \log p + (1-x)\log(1-p) gives I(p)=1/{p(1p)}I(p) = 1/\{p(1-p)\}, which diverges as pp saturates, and correspondingly the Cramér-Rao bound p(1p)/np(1-p)/n for unbiased estimators of pp shrinks there rather than growing. Both statements are correct because information belongs to a model and a parameterisation together, the two being tied by the chain rule I(p)=I(η)(dη/dp)2I(p) = I(\eta)\,(d\eta/dp)^{2}; here dη/dp=1/{p(1p)}d\eta/dp = 1/\{p(1-p)\}, and the square of that factor converts p(1p)p(1-p) into 1/{p(1p)}1/\{p(1-p)\} exactly. So the honest reading of the flattening curve is that a saturated Bernoulli is uninformative about its log-odds, not that it is uninformative full stop. Week 12 makes the transformation rule quantitative and shows which questions the bound answers.

Family h(x)h(x) T(x)T(x) η\eta A(η)A(\eta) A(η)A'(\eta) A(η)A''(\eta)
Bernoulli(p)(p) 11 on {0,1}\{0,1\} xx log{p/(1p)}\log\{p/(1-p)\} log(1+eη)\log(1 + e^{\eta}) pp p(1p)p(1-p)
Poisson(λ)(\lambda) 1/x!1/x! xx logλ\log \lambda eηe^{\eta} λ\lambda λ\lambda
Exponential(rate λ\lambda) 11 on (0,)(0, \infty) xx λ-\lambda log(η)-\log(-\eta) 1/λ1/\lambda 1/λ21/\lambda^2

Check at least one row. For the exponential, A(η)=log(η)A(\eta) = -\log(-\eta) gives A(η)=1/η=1/λA'(\eta) = -1/\eta = 1/\lambda, the mean, and A(η)=1/η2=1/λ2A''(\eta) = 1/\eta^2 = 1/\lambda^2, the variance. The natural parameter is λ-\lambda, so the natural parameter space is the negative half-line and a sign slip flips both derivatives.

Why the uniform sits outside the family

Suppose Uniform(0,θ)(0, \theta) could be written in exponential form. The exponential factor is strictly positive wherever defined, so for each θ\theta the set where the density is positive equals the set where hh is positive. The second set does not involve θ\theta; the first is the interval from 00 to θ\theta, which differs for θ=1\theta = 1 and θ=2\theta = 2. Contradiction — and not one a cleverer parameterisation could dodge, since common support is part of the definition.

A box holds eight named families including Bernoulli, Poisson and normal; two sit outside it: the uniform on zero to theta, whose support moves with the parameter, and the Cauchy, which admits no such factorization.

Familiar families inside the exponential-family form, with the uniform and the Cauchy outside it for different reasons.

Almost everything familiar is inside, which is why the class feels universal until you meet a member of the complement, and the two exclusions drawn fail for different reasons. The uniform’s support moves, and the argument above disposes of it without looking at the density. The Cauchy location family has the whole real line as support for every parameter value, so the support test passes, and it is still not an exponential family: no finite factorisation of the required shape exists for f(xμ)=1/[π{1+(xμ)2}]f(x \mid \mu) = 1/[\pi\{1 + (x - \mu)^2\}]. Test the support first, then hunt for the factorisation.

The consequences are why the course keeps returning to the uniform. The interchange of derivative and integral fails, so the score has no reason to have mean zero (Week 7); the information bound of Week 12 has no force, and the maximum beats the rate it would suggest; the likelihood is not differentiable at its maximum, so the score equation is useless.

Now do not overcorrect, because the list of casualties is shorter than students expect. Everything lost above is analytic, a consequence of the moving region of integration. The reduction theory survives intact. The joint density of the sample is θn\theta^{-n} when 0<x(1)0 \lt x_{(1)} and x(n)θx_{(n)} \le \theta, and zero otherwise, which factors through x(n)x_{(n)} alone — an indicator doing the work an exponent usually does — so Week 10’s criterion makes X(n)X_{(n)} sufficient, and in fact minimal sufficient. It is also complete, and the argument is short enough to run here. Suppose 𝔼θ[g(X(n))]=0\mathbb{E}_\theta[g(X_{(n)})] = 0 for every θ>0\theta \gt 0. The density of X(n)X_{(n)} is ntn1/θnn t^{n-1}/\theta^{n} on [0,θ][0, \theta], the k=nk = n case of the marginal formula above, so the assumption says 0θg(t)tn1dt=0\int_0^{\theta} g(t)\, t^{n-1}\, dt = 0 for every θ>0\theta \gt 0; differentiating in θ\theta gives g(θ)θn1=0g(\theta)\,\theta^{n-1} = 0 for almost every θ\theta, hence g=0g = 0 almost everywhere, which is completeness. Sufficiency plus completeness is precisely the input Week 11’s Lehmann-Scheffé theorem wants, and it will hand back n+1nX(n)\tfrac{n+1}{n}X_{(n)} as the unique unbiased estimator of θ\theta with smallest variance — the model with no information bound and no usable score nonetheless has a clean optimality result. Neither sufficiency nor completeness is the private property of exponential families; the smooth theory built on top of them is.

Worked example — the largest observation in a uniform sample

Setting. A workshop instruments a batch-testing rig that reports after a wait modelled as uniform on (0,θ)(0, \theta) seconds, where θ\theta is an unknown polling ceiling. Eight independent runs give waits of 2.62.6, 7.17.1, 4.84.8, 9.29.2, 1.31.3, 8.48.4, 5.55.5 and 3.93.9. The estimand is θ\theta, so n=8n = 8 and the observed maximum is 9.29.2.

Step 1 — the distribution function. The largest observation is at most tt exactly when every observation is, and they are independent, so for 0tθ0 \le t \le \theta,

P(X(8)t)=i=18P(Xit)=(tθ)8.P(X_{(8)} \le t) = \prod_{i=1}^{8} P(X_i \le t) = \left(\frac{t}{\theta}\right)^{8}.

Step 2 — the density. Differentiating in tt gives f(8)(t)=8t7/θ8f_{(8)}(t) = 8 t^{7}/\theta^{8} on [0,θ][0, \theta]. Cross-check with the general marginal density at k=n=8k = n = 8: the coefficient is 8!/(7!0!)=88!/(7!\,0!) = 8, F(t)7F(t)^{7} is t7/θ7t^7/\theta^7, {1F(t)}0\{1 - F(t)\}^{0} is one, and f(t)=1/θf(t) = 1/\theta. The product is the same.

Step 3 — the mean, and the bias.

𝔼[X(8)]=0θt8t7θ8dt=8θ8θ99=89θ.\mathbb{E}[X_{(8)}] = \int_0^{\theta} t \cdot \frac{8 t^{7}}{\theta^{8}}\, dt = \frac{8}{\theta^{8}} \cdot \frac{\theta^{9}}{9} = \frac{8}{9}\,\theta.

In general 𝔼[X(n)]=nθ/(n+1)\mathbb{E}[X_{(n)}] = n\theta/(n+1), so the bias is θ/(n+1)-\theta/(n+1), here θ/9-\theta/9. The sign is no accident of the algebra: X(n)θX_{(n)} \le \theta for every possible sample, so the estimator cannot average out to θ\theta. Multiplying by (n+1)/n(n+1)/n removes the bias exactly, and 98×9.2=10.35\tfrac{9}{8} \times 9.2 = 10.35 is the corrected estimate. That corrected estimator is not merely one unbiased option among many: since X(n)X_{(n)} is complete and sufficient here, the Lehmann-Scheffé argument of Week 11 makes n+1nX(n)\tfrac{n+1}{n}X_{(n)} the unbiased estimator of θ\theta with the smallest variance there is. Whether to insist on being unbiased at all is a different and harder question, a mean-squared-error one, and that is Week 6’s business.

Step 4 — the variance. The second moment is

𝔼[X(8)2]=0θt28t7θ8dt=8θ8θ1010=45θ2,\mathbb{E}[X_{(8)}^{2}] = \int_0^{\theta} t^{2} \cdot \frac{8 t^{7}}{\theta^{8}}\, dt = \frac{8}{\theta^{8}} \cdot \frac{\theta^{10}}{10} = \frac{4}{5}\,\theta^{2},

so the variance is 45θ26481θ2=8810θ20.00988θ2\tfrac{4}{5}\theta^2 - \tfrac{64}{81}\theta^2 = \tfrac{8}{810}\theta^2 \approx 0.00988\,\theta^2 and the standard deviation about 0.0994θ0.0994\,\theta. The general formula nθ2/{(n+2)(n+1)2}n\theta^2/\{(n+2)(n+1)^2\} returns 8θ2/8108\theta^2/810, the same number. If θ\theta is near ten seconds, the maximum of eight runs has a standard deviation of about one second, and all of that spread lies below θ\theta.

Step 5 — what the numbers say. With X(8)=9.2X_{(8)} = 9.2 we know with certainty that θ9.2\theta \ge 9.2, a guarantee no other estimator here supplies, and we know the estimate runs low by about eleven percent on average. As a calibration, P(X(8)0.9θ)=0.98=0.430P(X_{(8)} \le 0.9\,\theta) = 0.9^{8} = 0.430: in more than two samples out of five the maximum lands below ninety percent of the ceiling.

What this assumed. Independence, a genuinely uniform wait, and a ceiling that does not drift. Uniformity does more work than it looks: (t/θ)n(t/\theta)^n is a statement about the far right tail of FF, which is where a modelling assumption is least constrained by data near the centre. If the true density thins near θ\theta, every number here changes while the sample mean barely notices.

The same reasoning, transferred

Let X1,,XnX_1, \dots, X_n be independent Exponential draws with rate λ\lambda, so P(Xi>t)=eλtP(X_i \gt t) = e^{-\lambda t}, and look at the smallest observation. The smallest exceeds tt exactly when all of them do, so

P(X(1)>t)=i=1nP(Xi>t)=enλt,P(X_{(1)} \gt t) = \prod_{i=1}^{n} P(X_i \gt t) = e^{-n \lambda t},

the survival function of an Exponential with rate nλn\lambda. Hence 𝔼[X(1)]=1/(nλ)\mathbb{E}[X_{(1)}] = 1/(n\lambda) and Var(X(1))=1/(nλ)2\operatorname{Var}(X_{(1)}) = 1/(n\lambda)^2.

The engine stayed the same: an event about one extreme became an event about all nn observations at once, independence turned it into a product of identical factors, and the distribution followed. What changed is which tail is probed. The uniform maximum left its family and became a Beta variable; the exponential minimum stays exponential, with the rate multiplied by nn.

The consequence warns against a lazy generalisation. Since 𝔼[nX(1)]=1/λ\mathbb{E}[n X_{(1)}] = 1/\lambda, the rescaled minimum is exactly unbiased for the mean. But its variance is n2/(nλ)2=1/λ2n^2/(n\lambda)^2 = 1/\lambda^2 for every nn: it never shrinks, so the estimator is unbiased and not consistent. The uniform maximum, by contrast, has standard deviation shrinking like 1/n1/n, faster than the usual 1/n1/\sqrt{n}. Same class of statistic, opposite behaviour: the shape of FF near the boundary decides, not the word “extreme”.

Second worked example — the two-parameter normal in natural form

Setting. Let X1,,XnX_1, \dots, X_n be i.i.d. Normal(μ,σ2)(\mu, \sigma^2), both parameters unknown. This exercises a different facet of the week: not a sampling distribution, but canonical form and the moments read off a cumulant function.

Step 1 — expand the exponent.

(xμ)22σ2=x22σ2+μxσ2μ22σ2.-\frac{(x - \mu)^2}{2\sigma^2} = -\frac{x^2}{2\sigma^2} + \frac{\mu x}{\sigma^2} - \frac{\mu^2}{2\sigma^2}.

Step 2 — read off the ingredients. Absorb (2π)1/2(2\pi)^{-1/2} into hh and the θ\theta-dependent leftovers into AA:

f(xμ,σ2)=12πexp{μσ2x12σ2x2μ22σ212logσ2}.f(x \mid \mu, \sigma^2) = \frac{1}{\sqrt{2\pi}} \exp\Big\{ \frac{\mu}{\sigma^2} x - \frac{1}{2\sigma^2} x^{2} - \frac{\mu^2}{2\sigma^2} - \frac{1}{2}\log \sigma^{2} \Big\}.

So h(x)=(2π)1/2h(x) = (2\pi)^{-1/2}, T(x)=(x,x2)T(x) = (x, x^2), η1=μ/σ2\eta_1 = \mu/\sigma^2, η2=1/(2σ2)\eta_2 = -1/(2\sigma^2), and A=μ2/(2σ2)+12logσ2A = \mu^2/(2\sigma^2) + \tfrac{1}{2}\log\sigma^2. The support is the whole line for every parameter value, so the structural condition holds.

Step 3 — invert to natural coordinates. From η2=1/(2σ2)\eta_2 = -1/(2\sigma^2) we get σ2=1/(2η2)\sigma^2 = -1/(2\eta_2), forcing η2<0\eta_2 \lt 0, and then μ=η1σ2=η1/(2η2)\mu = \eta_1\sigma^2 = -\eta_1/(2\eta_2). So =×(,0)\mathcal{H} = \mathbb{R} \times (-\infty, 0) and

A(η)=η124η212log(2η2).A(\eta) = -\frac{\eta_1^{2}}{4 \eta_2} - \frac{1}{2}\log(-2\eta_2).

Step 4 — the gradient gives the means.

Aη1=η12η2=μ,Aη2=η124η2212η2=μ2+σ2.\frac{\partial A}{\partial \eta_1} = -\frac{\eta_1}{2\eta_2} = \mu, \qquad \frac{\partial A}{\partial \eta_2} = \frac{\eta_1^{2}}{4\eta_2^{2}} - \frac{1}{2\eta_2} = \mu^{2} + \sigma^{2}.

Those are 𝔼[X]\mathbb{E}[X] and 𝔼[X2]\mathbb{E}[X^2], exactly as the general theorem promised.

Step 5 — the Hessian gives the covariances.

2Aη12=12η2=σ2,2Aη1η2=η12η22=2μσ2,2Aη22=η122η23+12η22=4μ2σ2+2σ4.\frac{\partial^2 A}{\partial \eta_1^{2}} = -\frac{1}{2\eta_2} = \sigma^{2}, \qquad \frac{\partial^2 A}{\partial \eta_1 \partial \eta_2} = \frac{\eta_1}{2\eta_2^{2}} = 2\mu\sigma^{2}, \qquad \frac{\partial^2 A}{\partial \eta_2^{2}} = -\frac{\eta_1^{2}}{2\eta_2^{3}} + \frac{1}{2\eta_2^{2}} = 4\mu^{2}\sigma^{2} + 2\sigma^{4}.

Audit all three against moments obtained another way. The first is Var(X)=σ2\operatorname{Var}(X) = \sigma^2. For the second, the normal third moment is μ3+3μσ2\mu^3 + 3\mu\sigma^2, so Cov(X,X2)=μ3+3μσ2μ(μ2+σ2)=2μσ2\operatorname{Cov}(X, X^2) = \mu^3 + 3\mu\sigma^2 - \mu(\mu^2 + \sigma^2) = 2\mu\sigma^2. For the third, the fourth moment is μ4+6μ2σ2+3σ4\mu^4 + 6\mu^2\sigma^2 + 3\sigma^4, and subtracting (μ2+σ2)2(\mu^2 + \sigma^2)^2 leaves 4μ2σ2+2σ44\mu^2\sigma^2 + 2\sigma^4. All three match.

Step 6 — put numbers through it. Fix μ=3\mu = 3 and σ2=4\sigma^2 = 4, so η1=0.75\eta_1 = 0.75 and η2=0.125\eta_2 = -0.125.

Quantity From the natural form Value From the ordinary form Value
AA η12/(4η2)12log(2η2)-\eta_1^2/(4\eta_2) - \tfrac12\log(-2\eta_2) 1.8181.818 μ2/(2σ2)+12logσ2\mu^2/(2\sigma^2) + \tfrac12\log\sigma^2 1.8181.818
A/η1\partial A/\partial\eta_1 η1/(2η2)-\eta_1/(2\eta_2) 33 𝔼[X]=μ\mathbb{E}[X] = \mu 33
A/η2\partial A/\partial\eta_2 η12/(4η22)1/(2η2)\eta_1^2/(4\eta_2^2) - 1/(2\eta_2) 1313 𝔼[X2]=μ2+σ2\mathbb{E}[X^2] = \mu^2 + \sigma^2 1313
2A/η12\partial^2 A/\partial\eta_1^2 1/(2η2)-1/(2\eta_2) 44 Var(X)=σ2\operatorname{Var}(X) = \sigma^2 44
2A/η1η2\partial^2 A/\partial\eta_1\partial\eta_2 η1/(2η22)\eta_1/(2\eta_2^2) 2424 Cov(X,X2)=2μσ2\operatorname{Cov}(X, X^2) = 2\mu\sigma^2 2424
2A/η22\partial^2 A/\partial\eta_2^2 η12/(2η23)+1/(2η22)-\eta_1^2/(2\eta_2^3) + 1/(2\eta_2^2) 176176 Var(X2)=4μ2σ2+2σ4\operatorname{Var}(X^2) = 4\mu^2\sigma^2 + 2\sigma^4 176176

Step 7 — go to the sample. For nn observations the natural statistic becomes (iXi,iXi2)\big(\sum_i X_i, \sum_i X_i^2\big) and the cumulant function nA(η)nA(\eta), so 𝔼[iXi]=nμ\mathbb{E}\big[\sum_i X_i\big] = n\mu, 𝔼[iXi2]=n(μ2+σ2)\mathbb{E}\big[\sum_i X_i^2\big] = n(\mu^2 + \sigma^2), and the covariance matrix is nn times the Hessian above. Week 4 builds the exact distributions of X\bar{X} and S2S^2 from that pair, and Week 10 proves it minimal sufficient.

Reading the same two derivatives for the Bernoulli

The Bernoulli is the one-parameter rehearsal. Writing f(xp)=px(1p)1xf(x \mid p) = p^{x}(1-p)^{1-x} on {0,1}\{0, 1\} and taking logs inside an exponential,

f(xp)=exp{xlogp1p+log(1p)},f(x \mid p) = \exp\Big\{ x \log \frac{p}{1-p} + \log(1-p) \Big\},

so h(x)=1h(x) = 1 on {0,1}\{0,1\}, T(x)=xT(x) = x, the natural parameter is the log-odds η=log{p/(1p)}\eta = \log\{p/(1-p)\}, and A=log(1p)A = -\log(1-p). Inverting, p=eη/(1+eη)p = e^{\eta}/(1 + e^{\eta}), so A(η)=log(1+eη)A(\eta) = \log(1 + e^{\eta}) — the curve in the figure above. One derivative returns pp and two return p(1p)p(1-p), and for a sample the statistic iXi\sum_i X_i has mean npnp and variance np(1p)np(1-p) read straight off nAnA' and nAnA''.

A numerical check worth running

Neither result needs simulation to be true, but running one is how you catch a dropped factor. This block checks the uniform maximum against Steps 3, 4 and 5.

set.seed(75063)
theta <- 10; n <- 8; reps <- 50000
mx <- replicate(reps, max(runif(n, 0, theta)))

mean(mx)                  # near 8.889  = n * theta / (n + 1)
sd(mx)                    # near 0.994  = theta * sqrt(n) / ((n + 1) * sqrt(n + 2))
mean((n + 1) / n * mx)    # near 10     = theta, the bias-corrected estimator
mean(mx <= 0.9 * theta)   # near 0.430  = 0.9^8

The second differences AA rather than differentiating it, so an algebra slip shows up.

cumulant <- function(e) log(1 + exp(e))
eta <- 0.8; h <- 1e-4

(cumulant(eta + h) - cumulant(eta - h)) / (2 * h)                    # near 0.6900
(cumulant(eta + h) - 2 * cumulant(eta) + cumulant(eta - h)) / h^2    # near 0.2139

exp(eta) / (1 + exp(eta))                                            # 0.6900

A close match at one parameter value with fifty thousand replicates rules out gross algebra errors and sign flips, which is most of what goes wrong. It is not a proof, though: a wrong derivation can agree with a simulation by coincidence.

The misreading to avoid

Here is the sentence as graduate students write it: “the largest observation is the natural estimator of the upper endpoint, so it should be unbiased for θ\theta, and with a reasonable sample size I can attach a symmetric two-standard-error interval to it.” Each clause is wrong, instructively so.

The estimator is not unbiased and cannot be. Every possible sample satisfies X(n)θX_{(n)} \le \theta, so the sampling distribution lives below θ\theta and its mean is strictly below it, short by exactly θ/(n+1)\theta/(n+1). That follows from never being able to overshoot; a bigger study does not fix it.

The normal approximation is not merely inaccurate; it is the wrong shape at every nn. For 0tn0 \le t \le n,

P(n(θX(n))θt)=P(X(n)θ(1tn))=(1tn)net.P\Big( \frac{n(\theta - X_{(n)})}{\theta} \ge t \Big) = P\Big( X_{(n)} \le \theta\big(1 - \tfrac{t}{n}\big)\Big) = \Big(1 - \frac{t}{n}\Big)^{n} \longrightarrow e^{-t}.

The correct normalising factor is nn, not n\sqrt{n}, and the limit is Exponential with rate one: one-sided, right-skewed, with no left tail. The exact distribution offers a simple one-sided interval instead, since X(n)/θX_{(n)}/\theta has distribution function unu^{n} on (0,1)(0,1) and so does not involve θ\theta. Taking u=α1/nu = \alpha^{1/n} gives P(θX(n)α1/n)=1αP\big(\theta \le X_{(n)}\,\alpha^{-1/n}\big) = 1 - \alpha; for the eight observations above with α=0.05\alpha = 0.05, the factor is 0.051/8=1.4540.05^{-1/8} = 1.454 and the interval runs from 9.29.2 to 13.3813.38. It is asymmetric because the situation is. Week 4 names the object that made this work: a pivot.

Behind the specific error sits a general one: treating the order statistics as though they were themselves an i.i.d. sample. They are not identically distributed, as the three Beta densities plotted earlier show, and not independent, as the covariance i(nj+1)/{(n+1)2(n+2)}i(n-j+1)/\{(n+1)^2(n+2)\} shows. Any argument that quietly adds variances of order statistics, or applies the central limit theorem to one without asking whether it is central or extreme, is unsound. The earlier qualification stands in both directions: the median, whose index n/2\lceil n/2 \rceil meets the rate condition, really is asymptotically normal at rate n\sqrt{n}, and the failure is specific to extremes and to indices that drift away from their level too slowly.

A companion misreading sits on the exponential-family side: “the families I know are all exponential families, so the regularity conditions are safe to assume.” Membership can fail with a moving support or a fixed one, and the check takes a single line.

Practice on your own

Work these with a pencil first; they are for your own checking.

  1. Differentiate F(k)(y)=j=kn(nj)Fj(1F)njF_{(k)}(y) = \sum_{j=k}^{n} \binom{n}{j} F^{j} (1-F)^{n-j} in yy. Show the terms cancel in pairs, that what remains is the marginal density above, and that it collapses correctly at k=1k = 1 and k=nk = n.
  2. For a uniform sample on (0,1)(0,1) with n=3n = 3, compute the three variances from the Beta formula, add them, and compare with 3Var(X1)=1/43\operatorname{Var}(X_1) = 1/4. Account for the difference using the covariance formula, and say why it has the sign it does.
  3. A classmate argues: “sorting only relabels, so kX(k)=iXi\sum_k X_{(k)} = \sum_i X_i; therefore Var(X(1))++Var(X(n))\operatorname{Var}(X_{(1)}) + \cdots + \operatorname{Var}(X_{(n)}) equals nVar(X1)n \operatorname{Var}(X_1).” Locate the step that fails and use item 2 to say by how much.
  4. Membership hunt, one line of justification each: Uniform(θ,θ+1)(\theta, \theta + 1); Normal(μ,μ2)(\mu, \mu^{2}) with μ>0\mu \gt 0; Binomial(m,p)(m, p) with mm known; Binomial(m,p)(m, p) with both unknown. For any that fail, name the requirement that fails.
  5. Describe a simulation, in enough detail for someone else to code, separating the two candidate limit laws for the uniform maximum: histogram n(θX(n))/θn(\theta - X_{(n)})/\theta against an Exponential density with rate one, then histogram n(X(n)θ)\sqrt{n}(X_{(n)} - \theta) as nn grows. Say in advance what each picture should show if this unit is right.

Where to read more

Where this goes next

Week 4 takes the normal model apart exactly, starting from the pair (iXi,iXi2)\big(\sum_i X_i, \sum_i X_i^2\big) this unit produced: XN(μ,σ2/n)\bar{X} \sim N(\mu, \sigma^2/n), the scaled sample variance as a chi-square variable, and the independence of the two, which is a peculiarly normal phenomenon. The order-statistics half gets its sequel in Week 6’s mean-squared-error comparison, Week 7’s regularity counterexample, and Week 12’s bound that does not apply. Read week-04.qmd next; week-02.qmd set up the change-of-variables argument, and the notes index lists the sequence.

One forward connection is worth planting now. If the model is an exponential family with natural parameter η\eta and you take a prior π(η)exp{ητνA(η)}\pi(\eta) \propto \exp\{\eta^{\top}\tau - \nu A(\eta)\} for fixed constants τ\tau and ν\nu, then multiplying by the likelihood of a sample adds iT(xi)\sum_i T(x_i) to τ\tau and nn to ν\nu, and the posterior belongs to the same family. Conjugacy is not a coincidence and not a philosophy; it is this form closing under multiplication. That is the week’s Connect step, developed in Week 9: the structure you learned to spot here is what makes Bayesian updating tractable.