Past exam of the mathematics course of the University of Cambridge / 2004 / iii / Paper 42 (Statistical Theory)

A kernel for density estimation is an integrable function K of unit integral, ordinarily nonnegative and centered. Square integrability is needed for integrated variance; a symmetric second-order kernel for density estimation also has $\int uK(u)du=0$ and finite nonzero second moment ${\mu }_2(K)$. Its scaled form is $K_h(u)=h^{-1}K(u/h)$, and the kernel density estimator isHigher-order signed kernels are possible, but need not give a nonnegative probability density function estimate.

The mean integrated squared error is $\mathbb{E}\int ({\hat{f}}_h(x)-f(x){)}^{2}dx$. Put $g_h=K_h\ast f$. Independence gives the pointwise mean $g_h$ and variance $n^{-1}[(K_h^2\ast f)-g_h^2]$. Integrating, use $\int f=1$ and Fubini to obtain $\int (K_h^2\ast f)=R(K_h)$. The bias-variance decomposition proves the exact formulaEquivalently this is $R(f)-2\int f(K_h\ast f)+(1-n^{-1})R(K_h\ast f)+R(K)/(nh)$. These identities assume the displayed integrals are finite; they are exact, rather than a small-smoothing bandwidth expansion.

For a symmetric second-order kernel for density estimation, Taylor expansion gives leading bias $h^2{\mu }_2(K)f^{\prime \prime }/2$, and the leading integrated variance is $R(K)/(nh)$. Thus the asymptotic mean integrated squared error balances inverse smoothing bandwidth against fourth-power bias, with optimal smoothing bandwidth $[R(K)/(n{\mu }_2(K{)}^{2}R(f^{\prime \prime })){]}^{1/5}$. Its appearance of the integrated squared density curvature motivates the scale comparison.

For $f_a(x)=af(ax)$, $a>0$, two differentiations give $f_a^{\prime \prime }(x)=a^3{f}^{\prime \prime }(ax)$. A change of variables therefore gives $R(f_a^{\prime \prime })=a^{5}R(f^{\prime \prime })$, which tends to zero as $a\downarrow 0$ if the original curvature integral is finite. The distribution with probability density function $f_a$ is that of X divided by a, so $\sigma (f_a)=\sigma (f)/a$. ConsequentlyThis scale-invariant density curvature measures shape without dependence on the units used to measure X.

For the final minimization, write $c=35/32$ and use the triweight kernel $f_0=c(1-x^2{)}^3$ on $(-1,1)$, zero outside. Integration gives total mass one, mean zero and variance $1/9$. Its value and its first two derivatives vanish at both endpoints, so its extension is twice continuously differentiable. On the interior,In particular $f_0^{\prime \prime \prime }(1)=-48c$ and $f_0^{\prime \prime \prime }(-1)=48c$.

Put $e=h-f_0$. The mass and moment constraints give $\int e=\int x^{2}e=0$, since both means vanish and their variances match. If $R(h^{\prime \prime })=\infty$ the claim is immediate; otherwise expand the squared norm and integrate its cross term twice by parts over $[-1,1]$. The first endpoint term vanishes because $f_0^{\prime \prime }(\pm 1)=0$, but the next one must be retained:For the second equality use the zero mass and second moment of e and $e=h\ge 0$ outside the support. All terms in the last expression are nonnegative, including the endpoint values of a probability density function. It follows thatThis proves that the triweight density minimizes integrated squared curvature under the specified constraints. Equality forces $e^{\prime \prime }=0$; a globally affine integrable difference of probability density functions is zero, so equality occurs only at $h=f_0$.

The empirical distribution function is ${\hat{F}}_n(x)=n^{-1}{\sum }_i{1}_{X_i\le x}$. For continuous F, the probability integral transform gives independent uniform variables $U_i=F(X_i)$. Strict increase of F is unnecessary. More explicitly, generate the observations by $X_i=F^{-1}(U_i)$ with the generalized quantile inverse. Then $1_{X_i\le x}=1_{U_i\le F(x)}$, so ${\hat{F}}_n(x)=G_n(F(x))$, where $G_n$ is the uniform empirical distribution.

Continuity and the limits of F ensure that its range contains $(0,1)$. The endpoint discrepancies for $G_n$ are zero. Hencewhich depends only on n. Equivalently it is ${\max }_{1\le i\le n}\{i/n-U_{(i)},U_{(i)}-(i-1)/n\}$, a function of uniform order statistics. This proves the distribution-free property of the Kolmogorov-Smirnov statistic even when F has flat portions.

Choose $c_{n,\alpha }$ from this common distribution so that $\mathbb{P}(D_n\le c_{n,\alpha })=1-\alpha$. The Kolmogorov-Smirnov confidence band isThe event that all these inequalities hold is exactly the event used to choose the critical value, giving simultaneous coverage $1-\alpha$. The Glivenko-Cantelli theorem states that for an independent identically distributed sample from any distribution function, continuous or not, ${\sup }_x|{\hat{F}}_n(x)-F(x)|\to 0$ almost surely.

For an integrable r-argument U-statistic kernel, replace it by $h_s(x_1,\dots ,x_r)=(r!{)}^{-1}{\sum }_{\pi }h(x_{\pi (1)},\dots ,x_{\pi (r)})$. Every permutation has the same expectation by exchangeability, so this symmetrization of a U-statistic kernel preserves the target functional. For a symmetric U-statistic kernel h and $n\ge r$, the U-statistic isEach summand has expectation theta, hence the statistic is unbiased.

For the population variance, use $h(x,y)=(x-y{)}^2/2$. Independence gives $\mathbb{E}(X_1-X_2{)}^2=2\mathrm{Var}(X_1)$ when the second moment is finite. The corresponding statistic simplifies using ${\sum }_{i<j}(X_i-X_j{)}^2=n{\sum }_i(X_i-\overline{X}{)}^2$:Thus variance as a U-statistic is precisely the usual unbiased sample variance, for $n\ge 2$.

Let the summands have mean mu, variance ${\sigma }^2>0$, and cumulants ${\kappa }_j$. For the standardized sum $Z_n=(S_n-n\mu )/(\sigma \sqrt{n})$, expansion of its log-characteristic function gives a normal leading term and successively smaller cumulant terms. Fourier inversion turns these terms into Probabilists' Hermite polynomials. Under sufficient moments and a smooth nonlattice condition, the Edgeworth expansion beginsfor bounded z. Higher orders include the fourth cumulant and products of lower cumulants. This is a central expansion in powers of $n^{-1/2}$; it need not remain positive, normalized, or accurate in far tails after truncation.

For a saddlepoint density approximation, require the summand cumulant-generating function K to be finite on an interval about zero. At a sum value $s=ny$, find the interior saddle $\hat{t}$ satisfying $K^{\prime }(\hat{t})=y$. Exponential tilting makes the tilted mean y, so the tilted sum is approximated locally by its normal probability density at its mean. Undoing the tilt givesEquivalently this follows by deforming the transform-inversion contour through its stationary point. Suitable analytic and nonlattice regularity gives relative error $O(n^{-1})$ at fixed interior y. This approximation incorporates a nonquadratic tail exponent and is positive, while an Edgeworth truncation uses a normal probability density times a polynomial. It uses more information from the distribution and stronger exponential-moment assumptions. Its integral may differ from one by a higher-order amount, so positivity does not imply exact normalization.

For the centered Laplace distribution, $K(t)=t^2+t^4/2+\cdots$, hence ${\kappa }_2=2$ and ${\kappa }_3=0$. Its skewness vanishes, so the order-$n^{-1/2}$ Edgeworth correction is zero. With $z=s/\sqrt{2n}$, the requested expansion, omitting terms of order $n^{-1}$ in the standardized correction, isFor bounded z the exact probability density function is this expression times $1+O(n^{-1})$. The prefactor is $n^{-1/2}$ because this is the probability density function of the unstandardized sum; that change of units should not be confused with the orders inside the Edgeworth expansion.

For the saddlepoint, $K^{\prime }(t)=2t/(1-t^2)$. Write $y=s/n$ and $d=\sqrt{1+y^2}$. The root in $(-1,1)$ isThese formulas include $s=0$ without a division by y. Substitution yieldsThis is the leading saddlepoint approximation, without its order-$n^{-1}$ relative correction. Expanding d near one recovers the same central normal approximation, while retaining d captures the Laplace-sum large-deviation exponent away from the center.

Assume a symmetric second-order probability density function kernel for density estimation and enough smoothness for the displayed derivative functionals. The benchmark smoothing bandwidth is $h_0=[R(K)/(n{\mu }_2(K{)}^{2}R(f^{\prime \prime })){]}^{1/5}$, balancing $R(K)/(nh)$ against $h^4{\mu }_2(K{)}^{2}R(f^{\prime \prime })/4$. Three widely used selection methods are as follows.

For least-squares density cross-validation, expand integrated squared error as $R({\hat{f}}_h)-2\int {\hat{f}}_{h}f+R(f)$. Estimate the middle integral by leaving out the observation at which the estimate is evaluated, and minimizeThe left-out observation is independent of its estimator, so the second term has expectation $-2\int (K_h\ast f)f$. Thus $\mathbb{E}\mathrm{CV}(h)=\mathrm{MISE}(h)-R(f)$ exactly. The subtraction of the self-contribution is essential. Computation uses $R({\hat{f}}_h)=n^{-2}{\sum }_{i,j}(K_h\ast {\tilde{K}}_h)(X_i-X_j)$, where ${\tilde{K}}_h(u)=K_h(-u)$. This method directly targets integrated error without estimating a curvature functional first. Under standard smoothness and a suitable search interval its selected smoothing bandwidth is ratio-consistent with the MISE optimum, and hence with $h_0$. Relative fluctuations are ordinarily $O_p(n^{-1/10})$, so the selected values can be quite variable in finite samples.

For biased density cross-validation, estimate curvature from the second-derivative probability density function estimator and subtract its diagonal variance term:Differentiating the kernel for density estimation estimator and applying the integrated-variance identity gives the exact expectation $\mathbb{E}R({\hat{f}}_h^{\prime \prime })=R(K^{\prime \prime })/(n{h}^5)+(1-n^{-1})R(K_h\ast f^{\prime \prime })$. Under the derivative regularity assumptions, $R(K_h\ast f^{\prime \prime })=R(f^{\prime \prime })+O(h^2)$, so $\mathbb{E}{\hat{A}}_h=R(f^{\prime \prime })+O(h^2)+O(n^{-1})$. Minimize BCV over a reasonable smoothing bandwidth range. It is a smoothed criterion directed at AMISE; the name reflects its finite-sample bias rather than an exact unbiased-risk identity. Standard conditions again give ratio consistency and relative fluctuation order $n^{-1/10}$. It often has a smoother objective than ordinary least-squares cross-validation, though numerical behavior and relative variance depend on the probability density function and kernel for density estimation.

For plug-in bandwidth selection, obtain an independent pilot scale g, estimate $A=R(f^{\prime \prime })$, then setThe pilot should be positive and stable; consistent positive estimates make this expression ratio-consistent by continuity. A simple normal-reference initialization uses $R(f^{\prime \prime })=3/(8\sqrt{\pi }{\sigma }^5)$ under a normal reference model, substituting a sample scale for sigma. For a Gaussian density kernel this gives $h\approx 1.06\hat{\sigma }{n}^{-1/5}$. A pure normal-reference rule is quick but is asymptotically correctly calibrated only when its assumed curvature matches the actual probability density function; it is useful as the first stage of a data-driven plug-in procedure.

A more adaptive fourth-derivative pilot estimate of density curvature uses $R(f^{\prime \prime })=\int f{f}^{(4)}$, obtained by two integrations by parts, and estimates it byThe diagonal expectation contributes $K^{(4)}(0)/(n{g}^5)$; the leading smoothing bias is $-g^2{\mu }_2(K)R(f^{\prime \prime \prime })/2$. Balancing these terms produces the stated pilot scale $g\asymp R(f^{\prime \prime \prime }{)}^{-1/7}{n}^{-1/7}$, with another pilot or reference model supplying $R(f^{\prime \prime \prime })$. Under sufficient derivatives and tail control, its error is bounded byat that pilot scale. Taylor expansion of the fifth-root formula gives ${\hat{h}}_{\mathrm{PI}}/h_0-1=-(\hat{A}-A)/(5A)+o_p(|\hat{A}-A|)$, so it inherits that bound. A calibrated cancellation of leading pilot bias can give a faster rate. This illustrates why an appropriately constructed plug-in selector can converge more quickly than cross-validation, while requiring more smoothness and pilot choices. For each ratio-consistent method, its attained AMISE divided by the minimal AMISE tends to one.

A marginal statistical likelihood uses a statistic whose marginal sampling distribution depends on the parameter of interest but does not involve the nuisance parameter. It integrates out other observed coordinates. A profile likelihood instead keeps the full data and, for every fixed interest parameter, fits the nuisance parameter by maximizing the full likelihood. Profiling supplies a likelihood function but usually not a normalized sampling probability density function. Here the marginalization is of data coordinates; Bayesian model evidence averages parameters against a prior and is a separate construction.

Use shape-rate notation for the gamma distribution. The independent sums have $X\sim \mathrm{Gamma}(m,\psi \lambda )$ and $Y\sim \mathrm{Gamma}(n,\lambda )$, so their joint probability density for $x,y>0$ isFor $x=tu$, $y=u$ the absolute Jacobian determinant is u. ThereforeIntegrate u using the gamma integral ${\int }_0^{\infty }{u}^{m+n-1}{e}^{-bu}du=\Gamma (m+n)b^{-(m+n)}$. The lambda factors cancel, leavingThus T gives a nuisance-free marginal likelihood. Discarding only terms independent of psi, its logarithm isThe factor t in the logarithm is required, as in the original PDF.

The full log-likelihood is $m\log \psi +(m+n)\log \lambda -\lambda (\psi x+y)+c(x,y)$. Its lambda derivative vanishes uniquely atThe second derivative in lambda is negative. Substituting this maximum gives ${\ell }_p(\psi )=m\log \psi -(m+n)\log (\psi x+y)+c^{\prime }(x,y)$. Since $\psi x+y=y(1+\psi t)$, it agrees with ${\ell }_m(\psi ;t)$ up to a psi-independent constant. This proves the Gamma-ratio marginal and profile likelihood identity. Their common maximizer is $\hat{\psi }=m/(nt)=my/(nx)$, if an interest-parameter estimate is also wanted.

In a regular model let ${\ell }_N(\theta )$ be the log-likelihood, $U_N(\theta )=\nabla {\ell }_N(\theta )$ its score, $I_N(\theta )$ its total Fisher information, and $\hat{\theta }$ the unrestricted maximum-likelihood estimator. For a point null $\theta ={\theta }_0$ in a p-dimensional parameter space, the standard statistics areThese define the Wald test, score test and likelihood-ratio test. With an interior true parameter, increasing information, suitable differentiability and asymptotic normality, all have asymptotic chi-squared distribution with p degrees of freedom under the null. Observed information can replace expected information when it is asymptotically equivalent.

For q smooth independent restrictions $g(\theta )=0$, replace the Wald statistic by $g(\hat{\theta }{)}^T[G{I}_N(\hat{\theta }{)}^{-1}{G}^T{]}^{-1}g(\hat{\theta })$, where G is the derivative of g at the estimate; compare with chi-squared q. The likelihood-ratio statistic uses the difference between unrestricted and null-restricted maxima. The score statistic is computed at the restricted maximum and uses efficient information for the tested coordinates: with nuisance coordinates lambda this is $I_{\psi \psi }-I_{\psi \lambda }{I}_{\lambda \lambda }^{-1}{I}_{\lambda \psi }$. Inverting the acceptance inequalities for point nulls produces approximate coverage-$1-\alpha$ confidence regions; for example the Wald region consists of theta satisfying $(\hat{\theta }-\theta {)}^T{I}_N(\hat{\theta })(\hat{\theta }-\theta )\le {\chi }_{p,1-\alpha }^2$. Analogous score and likelihood-ratio inequalities give their confidence regions.

For the Poisson family, put $Y_{-1}=1$ to include the initial observation, and defineMultiplying the initial and conditional Poisson probabilities gives, on the supported observation sequences,The parameter-free term is ${\sum }_{i=0}^n[Y_i\log Y_{i-1}-\log (Y_i!)]$, interpreting $0\log 0$ as zero. A positive count after a zero preceding count is outside the support.

For $S_n>0$, strict concavity gives unconstrained maximum $\tilde{\theta }=S_n/T_n$, and the upper endpoint restriction givesThere is an endpoint qualification in the printed question. With probability $e^{-\theta }$, $Y_0=0$ and all later counts are zero. Then $S_n=0$, $T_n=1$ and the likelihood is $e^{-\theta }$, which has no maximum on the stated space $0<\theta \le 1$. The formula produces zero only as an extended estimate on the closure $[0,1]$. This is a likelihood supremum at an excluded boundary, and zero cannot be called an actual maximum in the original space.

The observed negative second derivative is $S_n/{\theta }^2$. Conditional expectation gives $\mathbb{E}{Y}_0=\theta$ and $\mathbb{E}{Y}_i=\theta \mathbb{E}{Y}_{i-1}={\theta }^{i+1}$. Therefore for $0<\theta <1$,This is the bounded information in a subcritical Poisson branching process. Observing more generations of one family does not provide unbounded information. Indeed $\mathbb{P}(Y_n>0)\le \mathbb{E}{Y}_n={\theta }^{n+1}\to 0$, and zero is absorbing, so extinction occurs at a finite generation almost surely. The total observed progeny $S_{\infty }$ is finite; eventually $S_n=S_{\infty }$, $T_n=1+S_{\infty }$, and the extended estimate equals $S_{\infty }/(1+S_{\infty })$, a discrete nondegenerate random variable rather than a consistent increasingly precise estimate.

Here is a direct proof of the failure of a chi-squared Wald limit which does not rely solely on the information bound. Under any interior null ${\theta }_0$, the event $Y_0=1$, $Y_1=0$ has probability ${\theta }_0{e}^{-2{\theta }_0}>0$ and forces every later count to be zero. On it the genuine constrained estimate is $1/2$ for every $n\ge 1$, so the usual expected-information Wald statistic isThus it has a persistent atom tending to $4(1/2-{\theta }_0{)}^2$. The fact that atoms obstruct a continuous distributional limit rules out a chi-squared limit, including when ${\theta }_0$ equals one half and that atom is at zero. Any definition adopted for the separate all-zero sample event leaves this positive-count counterexample intact. Hence the Wald statistic does not have the usual asymptotic chi-squared null law for this increasing observation horizon.