Past exam of the mathematics course of the University of Cambridge / 2007 / iii / Paper 49 (Nonparametric Statistical Theory)

A cubic spline belongs to $C^2[a,b]$ and is a polynomial of degree at most three on each interval separated by the spline knots. A natural cubic spline is additionally an affine function on the exterior intervals $[a,x_1]$ and $[x_n,b]$. On $[x_1,x_n]$ this is equivalent to the endpoint conditions $s^{\prime \prime }(x_1)=s^{\prime \prime }(x_n)=0$; these conditions remain part of the definition when an outer knot coincides with $a$ or $b$.

Let $s$ be the natural cubic spline interpolant of the prescribed values, and let $v$ be any other interpolating $C^2$ function. Set $u=v-s$, so $u(x_i)=0$. On the interval $[x_i,x_{i+1}]$, integration by parts twice givesbecause the fourth derivative of each cubic piece is zero. The second boundary term is zero because $u$ vanishes at the spline knots. After summing the first terms, the interior contributions cancel by continuity of $s^{\prime \prime }$ and $u^{\prime }$, and the two remaining terms vanish by the natural endpoint conditions. On the exterior intervals $s^{\prime \prime }=0$. Thus ${\int }_a^b{s}^{\prime \prime }{u}^{\prime \prime }=0$, and the second derivative roughness penalty satisfiesIf equality holds, the continuous nonnegative function $(u^{\prime \prime }{)}^2$ has zero integral, so $u^{\prime \prime }=0$ everywhere. Hence $u$ is an affine function, and its zeros at two distinct spline knots force $u=0$. This proves the minimum roughness property of the natural cubic spline interpolant and its equality case: the unique minimum-roughness interpolant is $s$.

For the fixed-design nonparametric regression fit, use the penalized least squares criterionLet $s_z$ be the natural cubic spline interpolant with value vector $z=(z_1,\dots ,z_n{)}^T$. Replacing $v$ by $s_z$ for $z_i=v(x_i)$ leaves the residual sum unchanged and reduces the second derivative roughness penalty, strictly unless $v=s_z$. It therefore remains to minimize the finite-dimensional quadratic formThe supplied spline roughness penalty matrix $K$ is positive semidefinite. For every nonzero vector $w$,Consequently $I+\lambda K$ is a positive-definite matrix, and $Q_{\lambda }$ is strictly convex. Its normal equations for linear least squares have the unique solutionThe preceding strict reduction also proves uniqueness over the whole class $C^2[a,b]$, not just over the natural cubic splines. This is the cubic smoothing spline. If the residual criterion is divided by $n$ instead, the corresponding formula is $(I+n\lambda K{)}^{-1}Y$, reflecting a change in the numerical convention for the smoothing parameter.

The independence and specified error mean and variance are not needed for existence or uniqueness of this minimizer. They do imply, with smoothing matrix $A_{\lambda }=(I+\lambda K{)}^{-1}$ and $g_X=(g(x_i){)}_i$, that $\mathbb{E}\hat{z}=A_{\lambda }{g}_X$ and $\mathrm{Cov}(\hat{z})={\sigma }^2{A}_{\lambda }^2$. The null space of $K$ consists of sampled affine functions: these have zero roughness, and conversely zero roughness forces their natural cubic spline interpolant to be an affine function. Thus these functions are reproduced without bias of an estimator, while other components are smoothed.

At a target $t$, local polynomial regression fits a polynomial of degree at most $p$ by weighted least squares:Define the response vector $Y=(Y_1,\dots ,Y_n{)}^T$, the design matrix $X$ by $X_{i,j+1}=(x_i-t{)}^j$ for $j=0,\dots ,p$, and the diagonal matrix $W$ by $W_{ii}=K((x_i-t)/h)$. The normal equations for linear least squares giveThis formula assumes the local polynomial Gram matrix $X^{T}WX$ is invertible. It is invertible whenever at least $p+1$ distinct sites have positive weights: a nonzero degree-$p$ polynomial cannot vanish at all of them. With fixed $p$, a nonzero continuous regression kernel, and $nh\to \infty$, this condition holds eventually at the target points considered here. Using weights $h^{-1}K((x_i-t)/h)$ instead gives exactly the same fit.

Use the usual unit-integral normalization ${\mu }_0(K)=1$. A regression kernel need not intrinsically have unit integral, but multiplying all weights by a positive constant leaves local polynomial regression unchanged. With unnormalized raw moments the interior coefficient below is ${\mu }_2(K)/{\mu }_0(K)$; the version written using ${\mu }_2(K)$ alone requires this normalization. A zero regression kernel cannot define an estimator.

Write the local moment sums asThe local constant estimator is ${\hat{m}}_h(t;0)={\sum }_{i}K((x_i-t)/h)Y_i/{\sum }_{i}K((x_i-t)/h)$. Since the errors have zero mean, its bias of an estimator is the corresponding weighted average of $m(x_i)-m(t)$. By the second-order Taylor expansion, uniformly for contributing sites with $|x_i-t|\le h$,where ${\eta }_n\to 0$. This uniform remainder follows from uniform continuity of $m^{\prime \prime }$ on $[0,1]$. Nonnegativity of the regression kernel therefore giveswhenever $A_0$ tends to a positive limit and $A_2=O(h^2)$.

For fixed interior $x$, the moment approximations give $A_0(x)=1+o(1)$, $A_1(x)=O(n^{-1})$, and $A_2(x)=h^2{\mu }_2(K)+o(h^2)$. The hypothesis $n{h}^2\to \infty$ makes $n^{-1}=o(h^2)$, so the odd-moment contribution is negligible at second order. ThusThe symmetry of the regression kernel is what removes the first-order interior bias of an estimator.

At $z_n=\alpha h$, the rescaled observation window is truncated to $[-\alpha ,1]$. Put ${\mu }_{r,\alpha }={\int }_{-\alpha }^1{u}^{r}K(u)du$ and ${\beta }_r={\mu }_{r,\alpha }/{\mu }_{0,\alpha }$. Symmetry, nonnegativity, and unit mass give ${\mu }_{0,\alpha }\ge 1/2$. The supplied moment approximations already imply the leading boundary bias of local constant regression:To retain the second-order term rigorously, we need slightly more precise control of the first moment. The rescaled grid $u_i=i/(nh)-\alpha$ has spacing $1/(nh)$. Each truncated function $u^{r}K(u)1_{[-\alpha ,1]}$ has bounded variation, because $K$ is continuously differentiable on its support. The Riemann sum error is bounded by this variation times the grid spacing: summing the oscillation bounds on the grid cells provesEndpoint cells, including the omitted site at zero, contribute the same error order. ConsequentlyBoth discretization errors are $o(h^2)$ under $n{h}^2\to \infty$. The preceding Taylor expansion therefore proves the second-order boundary expansion for local constant regression:Equivalently, expanding the derivatives about zero givesIn general the truncated first moment is nonzero, so the boundary bias of an estimator has order $h$, although it can be smaller when $m^{\prime }(0)=0$ or the omitted part of the regression kernel carries no mass.

For local linear regression, the polynomial reproduction property of local polynomial regression removes the constant and linear terms exactly. Let $\Delta =A_0{A}_2-A_1^2$. Its effective kernel weights areso ${\sum }_i{l}_i=1$, ${\sum }_i{l}_i(x_i-t)=0$, andAt either an interior or the specified boundary target, $\Delta$ is asymptotic to a positive constant times $h^2$. At the boundary that constant is ${\mu }_{0,\alpha }{\mu }_{2,\alpha }-{\mu }_{1,\alpha }^2>0$, the mass squared times the variance of the continuous normalized truncated regression kernel. Moreover ${\sum }_i|l_i|=O(1)$, as follows by bounding $|A_2-A_1(x_i-t)|=O(h^2)$ in the window. Thus the Taylor remainder still contributes $o(h^2)$, even though some effective kernel weights may be negative. In particular,At an interior point this coefficient reduces to ${\mu }_2(K)$ by symmetry and normalization. Hence local linear regression has bias $O(h^2)$ at both interior and boundary points; its linear reproduction removes the generic order-$h$ boundary term.

A nondegenerate distribution function $G$ is a max-stable distribution if, for every integer $r\ge 1$, there are $a_r>0$ and $b_r\in \mathbb{R}$ such thatEquivalently, the suitably centered and positively scaled sample maximum of $r$ independent and identically distributed random variables with law $G$ again has law $G$. The extremal types theorem states that every nondegenerate limiting law of normalized sample maxima is a max-stable distribution. Thus this property identifies all possible limiting types, although a particular starting distribution function need not admit any such limiting normalization.

The three standard types are the Gumbel distributionthe Fréchet distribution, with $\alpha >0$,and the negative Weibull distribution, with $\alpha >0$,They are max-stable distributions, as follows respectively by taking $(a_r,b_r)=(1,\log r)$, $(r^{1/\alpha },0)$, and $(r^{-1/\alpha },0)$. The necessary and sufficient classification is: $G$ is max-stable exactly when $G(x)=H((x-b)/a)$ for $a>0$, $b\in \mathbb{R}$, and one of these three types $H$. A positive affine map sends a max-stable distribution to another max-stable distribution.

For the requested threshold equivalence, set $p_n=1-F(u_n)$. By independence,If $n{p}_n\to \tau <\infty$, then $p_n\to 0$ and the expansion $\log (1-p)=-p+O(p^2)$ gives $n\log (1-p_n)\to -\tau$, since $n{p}_n^2=(n{p}_n)p_n\to 0$. Exponentiation proves the claimed limit. Conversely, if $(1-p_n{)}^n\to e^{-\tau }>0$, then $p_n\to 0$: any subsequence with $p_n\ge \delta >0$ would have its powers tending to zero. Taking logarithms now gives $-n\log (1-p_n)\to \tau$. Since $-\log (1-p)/p\to 1$, including the continuous extension at $p=0$, it follows that $n{p}_n\to \tau$. Hence

The exceedance count $S_n$ has the binomial distribution with parameters $n,p_n$. For $\tau >0$ and each fixed nonnegative integer $s$, its probability mass function satisfieswhere $(n{)}_s=n(n-1)\cdots (n-s+1)$ and $(n{)}_0=1$. Indeed $(n{)}_s/n^s\to 1$, $(n{p}_n{)}^s\to {\tau }^s$, and $(1-p_n{)}^{n-s}\to e^{-\tau }$. If $\tau =0$, the bound $\mathbb{P}(S_n\ge 1)\le \mathbb{E}{S}_n=n{p}_n\to 0$ gives the degenerate Poisson distribution at zero. Summing the finitely many terms proves, in both cases,This derives the needed instance of the Poisson limit theorem directly.

Finally fix $x$ with $G(x)>0$ and put $u_n=a_{n}x+b_n$. The assumed limit for the normalized sample maximum and the proved equivalence give $n{p}_n\to -\log G(x)$. The second-largest order statistic is at most $u_n$ exactly when at most one observation is strictly greater than $u_n$. This event identity holds even if observations have ties. Apply the preceding Poisson limit theorem with $k=1$ to obtainThe endpoint $G(x)=1$ corresponds to $\tau =0$ and is already included. This is the case $r=2$ of the exceedance count limit for extreme order statistics.

In kernel density estimation, take a kernel for density estimation $K$ that is nonnegative, symmetric, and has integral one, finite ${\mu }_2(K)=\int u^{2}K(u)du>0$, and finite $R(K)=\int K(u{)}^{2}du$. Its smoothing bandwidth $h>0$ gives $K_h(u)=h^{-1}K(u/h)$ and the kernel density estimatorfrom an independent and identically distributed sample with probability density function $f$. Nonnegativity and unit mass make the estimate itself a probability density function. We compare kernel for density estimation shapes by their optimized mean integrated squared error, rather than by their performance at an arbitrarily common numerical smoothing bandwidth.

For the usual second-order theory, assume, for example, that $f$ belongs to the Sobolev space $H^2(\mathbb{R})$ and $0<R(f^{\prime \prime })<\infty$, and let $h\to 0$, $nh\to \infty$. The expectation is the convolution $K_h\ast f$. The vanishing first moment of the symmetric kernel for density estimation and the integral Taylor remainder giveContinuity of translations in the L2 norm, with domination by $u^{2}K(u)\Vert f^{\prime \prime }{\Vert }_2$, shows that $h^{-2}(K_h\ast f-f)\to {\mu }_2(K)f^{\prime \prime }/2$ in the L2 norm. Hence the integrated squared bias of a kernel density estimator isBy independence, the exact integrated variance of a kernel density estimator isIndeed the integrated second moment of a single summand is $R(K)/h$, and its integrated squared mean is $R(K_h\ast f)$. The latter quantity is bounded by $R(f)$, by Young's convolution inequality, so its contribution is $O(n^{-1})$. The bias-variance decomposition of mean squared error therefore yields the leading criterionThis asymptotic mean integrated squared error displays the bias-variance tradeoff: a small smoothing bandwidth increases integrated variance, while a large one increases squared bias of an estimator.

Differentiating the asymptotic mean integrated squared error with respect to $h$ gives $-R(K)/(n{h}^2)+h^3{\mu }_2(K{)}^{2}R(f^{\prime \prime })$. It changes sign once, from negative to positive, so its unique minimum isAt this minimum the integrated variance term is four times the squared bias of an estimator term. Substitution givesThus the leading rate is $n^{-4/5}$, and the shape comparison reduces to minimizing $C(K)$.

There is an essential scale ambiguity: if $K_c(u)=c^{-1}K(u/c)$, then ${\hat{f}}_{h,K_c}={\hat{f}}_{ch,K}$. Changing the scale of the kernel for density estimation merely changes the numerical convention for the smoothing bandwidth. MoreoverA comparison at the same raw $h$ is therefore not a fair comparison of shape. One can normalize the second moment to one, but that does not make the optimal smoothing bandwidth independent of shape.

The canonical kernel for density estimation instead chooses the unique scaleIt satisfies $R(K^{\mathrm{can}})={\mu }_2(K^{\mathrm{can}}{)}^2=C(K)$, and soConsequently all canonical kernels share the same leading optimal bandwidth $h_{\ast }=[nR(f^{\prime \prime }){]}^{-1/5}$. Only the multiplicative constant $C(K)$ depends on shape. This is precisely the separation of shape choice and scale choice that makes the canonical kernel normalization useful. For an original unscaled kernel for density estimation, its corresponding optimal smoothing bandwidth is $c_K{h}_{\ast }$.

To prove optimal shape, scale a candidate kernel for density estimation to have second moment one. Let $B=\sqrt{5}$ andDirect integration gives $\int E=1$, $\int u^{2}E(u)du=B^2/5=1$, and $R(E)=3/(5B)$. Thus $E$ is a rescaled Epanechnikov kernel. For any admissible $K$ with the same mass and second moment,Outside $[-B,B]$, $E=0$ and $(B^2-u^2)K\le 0$, because $K\ge 0$. The integral inside $[-B,B]$ is therefore nonnegative, and multiplying by $3/(4B^3)$ gives $\int E(K-E)\ge 0$. Expanding the square now provesEquality forces $K=E$ almost everywhere. At fixed second moment, minimizing $R(K)$ is equivalent to minimizing $C(K)$. Since $C(K)$ is scale invariant, this proves the optimality of the Epanechnikov kernel over all nonnegative symmetric second-order kernels for density estimation, including the equality case up to rescaling.

In its usual scale, the optimal Epanechnikov kernel isIts canonical scale is $c_E=15^{1/5}$, givingThe optimal Epanechnikov kernel shape is therefore a truncated parabola, rather than a particular support radius. The canonical kernel and unit-second-moment version represent the same shape but different scale conventions.

The numerical gains over other reasonable kernels for density estimation are modest. Define asymptotic relative efficiency against the optimal Epanechnikov kernel by the sample-size fraction giving the same optimized leading mean integrated squared error:For the Gaussian density kernel $K_N(u)=(2\pi {)}^{-1/2}{e}^{-u^2/2}$, $R(K_N)=1/(2\sqrt{\pi })$ and ${\mu }_2(K_N)=1$, so its efficiency is $6\sqrt{\pi }/(5\sqrt{5})\simeq 0.9512$. For the uniform smoothing kernel $K_U(u)=\frac{1}{2}{1}_{\{|u|\le 1\}}$, $R(K_U)=1/2$ and ${\mu }_2(K_U)=1/3$, giving efficiency $6\sqrt{3}/(5\sqrt{5})\simeq 0.9295$. The Gaussian density kernel is infinitely smooth, whereas the optimal Epanechnikov kernel has compact support and allows each observation to affect only a finite window. Such computational and smoothness considerations may outweigh a small leading-efficiency difference.

In practice $R(f^{\prime \prime })$ is unknown. A plug-in estimator estimates this roughness functional from a pilot estimate of the second derivative and substitutes it into the optimal smoothing bandwidth formula. Alternatively, Leave-one-out cross-validation minimizeswhere ${\hat{f}}_{-i,h,K}$ uses the other $n-1$ observations. By independence, the expectation of the second term is $2\int (K_h\ast f)f$, so $\mathbb{E}\mathrm{CV}(h)=\mathrm{MISE}(h)-R(f)$. Thus it targets integrated error without knowing $f$. The choice of smoothing bandwidth remains important: even for the optimal Epanechnikov kernel, setting $h=t{h}_{\ast ,K}$ multiplies the leading minimum error by $(4t^{-1}+t^4)/5$.

The preceding optimum concerns nonnegative second-order kernels for density estimation and twice-smooth densities on the real line. Signed kernels of order ell can cancel higher moments for smoother densities and yield different rates, but may produce negative density estimates. A support boundary or a density jump also requires separate bias of a kernel density estimator analysis. These settings do not contradict the proved second-order optimality of the Epanechnikov kernel.