Past exam of the mathematics course of the University of Cambridge / 2009 / iii / Paper 34 (Statistical Theory)

For every nonzero $v\in {\mathbb{R}}^{n-p}$, full column rank of $L$ gives $Lv\ne 0$, and henceThus the Gram matrix $L^{T}L$ is a positive-definite matrix and is invertible. Using $BA=0$ and $A=L{L}^T$,so

A linear image of a multivariate normal vector is again multivariate normal. Therefore the stacked vector has mean $(B\mu ,L^T\mu {)}^T$ and covarianceIn particular,The two diagonal blocks are positive-definite by the rank assumptions. The two blocks are jointly Gaussian and have zero cross covariance. The fact that uncorrelated jointly Gaussian variables are independent makes $BY$ and $L^{T}Y$ independent. Sinceis a measurable function of the second block, independence is preserved under this transformation. This proves Gaussian independence of linear and quadratic statistics:

The log-likelihood in the normal linear model, up to an additive constant, isFull column rank makes $X^{T}X$ positive-definite, so minimizing the residual sum of squares gives the normal equations $X^{T}X\hat{\beta }=X^{T}Y$. Maximizing over ${\sigma }^2$ afterward gives the maximum-likelihood estimatorsThe residual sum of squares is positive almost surely because $n-p>0$ and ${\sigma }^2>0$, so this maximization has an interior variance solution almost surely. The first estimator is a linear Gaussian statistic with mean $\beta$ and covariance ${\sigma }^2(X^{T}X{)}^{-1}$, giving

Let $P_X=X(X^{T}X{)}^{-1}{X}^T$ be the orthogonal projection onto the design column space, and let $A=I-P_X$. This is a symmetric idempotent matrix of rank $n-p$. Choose an orthonormal basis for its range and collect its columns in $L$, so $A=L{L}^T$, $L^{T}L=I_{n-p}$ and $L^{T}X=0$. ThenThe definition of the chi-squared distribution therefore givesThe maximum-likelihood denominator is $n$; the unbiased residual-variance estimator uses $n-p$.

Finally put $B=(X^{T}X{)}^{-1}{X}^T$. It has rank $p$ and $BA=0$, because $X^T(I-P_X)=0$. Part (a) implies that $BY=\hat{\beta }$ is independent of $Y^{T}AY=n{\hat{\sigma }}^2$. Multiplication by a fixed constant preserves independence, so the normal linear model maximum-likelihood sampling distributions satisfy

Put $Q_n(\theta )=n^{-1}{\sum }_{i=1}^n\log f(Y_i;\theta )$ and $Q(\theta )={\mathbb{E}}_{{\theta }_0}\log f(Y_1;\theta )$. The integrable absolute-log envelope makes these expectations finite. Compactness and continuity of the likelihood ensure that a maximum exists, and maximizing it is equivalent almost surely to maximizing $Q_n$.

First identify the unique maximizer of $Q$. Let $r(y)=f(y;\theta )/f(y;{\theta }_0)$ on the support of the true density. Since $\log r\le r-1$,Equality requires $r=1$ almost surely under the true distribution and no mass outside its support, so the two probability distributions coincide. By identifiability, this happens only for $\theta ={\theta }_0$. Equivalently $Q({\theta }_0)-Q(\theta )$ is the Kullback-Leibler divergence from the true distribution to the candidate.

Fix $\epsilon >0$ and let $K_{\epsilon }=\{\theta \in \Theta :\Vert \theta -{\theta }_0\Vert \ge \epsilon \}$. If this set is empty there is nothing to prove. Otherwise it is compact, so continuity and the strict maximum implyOn the event ${\sup }_{\Theta }|Q_n-Q|<{\delta }_{\epsilon }/3$, every point of $K_{\epsilon }$ has $Q_n(\theta )<Q({\theta }_0)-2{\delta }_{\epsilon }/3$, whereas $Q_n({\theta }_0)>Q({\theta }_0)-{\delta }_{\epsilon }/3$. Hence no maximizer can belong to $K_{\epsilon }$. The supplied uniform law of large numbers therefore proves compact-parameter consistency of maximum likelihood:

Under the usual regularity conditions, with nonsingular one-observation Fisher information matrix $I({\theta }_0)$, the asymptotic normality of a maximum likelihood estimator iswhere $I({\theta }_0)={\mathbb{E}}_{{\theta }_0}[S_{{\theta }_0}(Y_1)S_{{\theta }_0}(Y_1{)}^T]$ and $S_{\theta }={\nabla }_{\theta }\log f(\cdot ;\theta )$ is the score function.

For the uniform model, let $M_n={\max }_i{Y}_i$. For a sample in its support, the likelihood is ${\theta }^{-n}{1}_{\{\theta \ge M_n\}}$, decreasing on its nonzero region. Thus the uniform endpoint maximum-likelihood estimator isThe all-zero sample is a probability-zero exception and may be assigned an arbitrary positive estimate. Independence gives ${\mathbb{P}}_{\theta }(M_n\le y)=(y/\theta {)}^n$ on $[0,\theta ]$. For $0<\epsilon <\theta$,for $\epsilon \ge \theta$ this probability is zero. This proves consistency directly.

Let $T_n=n(\theta -M_n)/\theta$. Its exact cumulative distribution function isFor fixed $t\ge 0$, this converges to $1-e^{-t}$, so $T_n\stackrel{d}{\to }\mathrm{Exp}(1)$. In particular the errors are of order $1/n$. Explicitly, for every fixed $\eta >0$ and sufficiently large $n$,ConsequentlyOne violated regularity condition is parameter-independent support: the density is supported on $[0,\theta ]$. Differentiation under the density integral is not valid in the usual regular-model form. Indeed the interior score equals $-1/\theta$, whose expectation is not zero; the changing endpoint contributes the omitted boundary term. Thus the uniform model is handled by its exact maximum distribution, rather than by the regular asymptotic normality theorem.

Put $P=I-p^{-1}{1}_p{1}_p^T$, $Z=PX=X-\bar{X}{1}_p$, $Q=\Vert Z{\Vert }^2$ and $a=p-3$. Componentwise,Thus the James–Stein shrinkage toward the sample mean leaves the average unchanged and multiplies every contrast about it by the same data-dependent factor $1-a/Q$. If $Q>a$ the contrasts shrink without changing sign; if $Q<a$ their signs reverse, and if $Q<a/2$ their magnitudes increase. The event $Q=0$ has probability zero, and the estimator can be assigned $\bar{X}{1}_p$ there.

The sample mean is a linear normal statistic with mean $\bar{\theta }$ and variance $1/p$, soThe error splits into orthogonal mean and contrast components:Since the second term is perpendicular to $1_p$, its cross product with the first is zero for every sample. Taking squared norms and expectations, $p\mathbb{E}(\bar{X}-\bar{\theta }{)}^2=1$, which yields the risk function

To compute the contrast risk, choose an orthonormal basis matrix $U$ for the range of $P$, with $U^{T}U=I_m$, $U{U}^T=P$ and $m=p-1\ge 3$. Then $W=U^{T}X\sim N_m(\eta ,I)$ with $\eta =U^T\theta$, $Q=\Vert W{\Vert }^2$, and the contrast norm is the norm of $(1-a/Q)W-\eta$. Expand it:For $g_j(w)=w_j/\Vert w{\Vert }^2$, the Gaussian Stein identity givesThis application at a singularity is valid: near zero the reciprocal squared radius is integrable in dimension $m>2$, since its radial integral is ${\int }_0^1{r}^{m-3}dr$. Apply Gaussian integration by parts outside a ball of radius $\epsilon$; the inner boundary contribution is bounded by a constant times ${\epsilon }^{m-2}$ and vanishes. The Gaussian tail eliminates the boundary at infinity. Therefore all required expectations are finite and the displayed identity holds.

Since $a=m-2=p-3$, the Stein risk identity now gives contrast risk $m-a^2\mathbb{E}(1/Q)$. Adding the unit mean risk proves

Finally $Q$ has the noncentral chi-squared distribution ${\chi }_m^2(\lambda )$ with $\lambda =\Vert P\theta {\Vert }^2$. To establish where its reciprocal expectation is largest, complete the square in the normal representation $Q=\Vert W{\Vert }^2$ to obtain its Laplace transform:Since $1/q={\int }_0^{\infty }{e}^{-tq}dt$, the Tonelli theorem gives the reciprocal moment of a noncentral chi-squared variableIts integrand is positive and strictly decreasing in $\lambda$ for every $t>0$. The expectation is therefore strictly largest at $\lambda =0$, whereBecause this expectation is subtracted in the risk, the minimizing set is exactly $P\theta =0$, that is, all coordinates of $\theta$ equal. On this line the risk is $p-(p-3)=3$:

Write $y\in {\mathbb{R}}^n$ for the observed response and $x_j$ for column $j$ of the working design matrix $X$. The usual LARS convention gives every covariate unit Euclidean norm. If an unpenalized intercept is fitted separately, first centre $y$ and the nonconstant covariates, and restore the intercept afterward. For an intercept-free model no centring is required. The algebra below also applies to an unscaled full-rank design; only the literal equal-angle interpretation uses unit-length columns. The penalty is defined on the regression coefficients of whichever working design is used, since rescaling covariates changes the units of an L1 penalty.

Initialize ${\hat{\beta }}^0=0$, ${\hat{\mu }}^0=X{\hat{\beta }}^0=0$ and residual $r^0=y$. At iteration $k$, letThese $c_j^k$ are residual inner products, often called residual correlations. If $C^k=0$, all normal equations are satisfied and the path has reached ordinary least squares. Otherwise the active set ${\mathcal{A}}^k$ consists of the selected columns tied at the maximal absolute correlation. Initially select the maximizing column, or all such columns in a tie. Ordinary least angle regression retains previously selected columns and adds each new one when it ties the active correlations. Put $s_j^k=\mathrm{sign}(c_j^k)$ for active $j$.

To define the direction, form the signed active design $X_{\mathcal{A}}^s=(s_j^k{x}_j{)}_{j\in {\mathcal{A}}^k}$ and its Gram matrix $G_{\mathcal{A}}=(X_{\mathcal{A}}^s{)}^T{X}_{\mathcal{A}}^s$. Full column rank makes this matrix positive-definite. Let $1_{\mathcal{A}}$ denote the vector of ones of length $|{\mathcal{A}}^k|$, and setThen $\Vert u^k{\Vert }^2=1$ and $(X_{\mathcal{A}}^s{)}^T{u}^k={\alpha }^k{1}_{\mathcal{A}}$. This is the equiangular direction in least angle regression. Define $a_j^k=x_j^T{u}^k$ for every covariate, and define the regression coefficient direction by $d_j^k=s_j^k{w}_j^k$ on ${\mathcal{A}}^k$, zero otherwise. It satisfies $X{d}^k=u^k$.

Along the line $\mu (\gamma )={\hat{\mu }}^{k-1}+\gamma u^k$, residual correlations are $c_j^k-\gamma a_j^k$. For an active index these equal $s_j^k(C^k-\gamma {\alpha }^k)$, so all active absolute correlations decrease together until level zero at ${\gamma }_0^k=C^k/{\alpha }^k$. An inactive covariate first joins when its correlation equals either sign of the common level. Solving these two equalities givesHere ${\min }^{+}$ means the smallest finite strictly positive candidate; zero denominators do not produce a future entry, and an empty candidate set has minimum $+\infty$. Candidates beyond the zero-correlation endpoint are capped by that endpoint. TakeAt an entry add the newly tied covariate, or all simultaneous ties, and recompute the direction. Once all covariates are active there is no inactive candidate: the final move is to correlation level zero and the full least-squares fit. These definitions specify every term in the stated step formula and the entry knot in least angle regression.

For definiteness use the following penalty normalization for the Lasso estimator:Full column rank makes the criterion strictly convex, so the minimizer is unique. With the alternative squared-error normalization $1/(2n)$, the penalty parameter is divided by $n$. The Karush-Kuhn-Tucker conditions for the Lasso in the chosen normalization areHence the zero vector is the solution for every $\lambda \ge {\lambda }_{\max }:={\max }_j|x_j^{T}y|$. Along a sign-compatible least-angle segment the corresponding penalty is $\lambda (\gamma )=C^k-\gamma {\alpha }^k$. Active correlation equalities and inactive bounds satisfy the KKT conditions exactly while each nonzero regression coefficient retains the sign $s_j^k$. Ordinary LARS can let a regression coefficient cross zero while its correlation sign stays unchanged; continuing past that crossing would violate the KKT conditions.

The required modification is therefore to monitor regression coefficient zero crossings. Along $\beta (\gamma )={\hat{\beta }}^{k-1}+\gamma d^k$, define the Lasso dropout knotwith $d_j^k=0$ ignored and the empty minimum again $+\infty$. Stop atIf a nonzero active regression coefficient reaches zero before the next entry, remove it from the active direction and recompute; the fitted value is continuous at this event. If entry comes first, add that covariate as before. A removed covariate is allowed to re-enter later, so the active set need not increase monotonically. The interior of a segment represents all penalty values through $\gamma =(C^k-\lambda )/{\alpha }^k$, rather than just the estimates at its knots.

For complete treatment of simultaneous events, full column rank alone does not imply that entries and dropouts occur singly. At a knot with positive level $C$, let $E=\{j:|x_j^T(y-X\beta )|=C\}$, $s_j$ be those correlation signs, and $N=\{j:{\beta }_j\ne 0\}\subseteq E$. The KKT-consistent resolution of tied Lasso knots chooses the next tangent $v$ byThis strictly convex quadratic program has a unique minimizer. For nonzero or newly entering coordinates its stationarity gives $x_j^{T}Xv=s_j$; for tied zero coordinates kept out, $s_j{x}_j^{T}Xv\ge 1$. Thus for a small decrease $h$ of the penalty, $\beta +hv$ has the correct regression coefficient signs, active residual scores $s_j(C-h)$, and inactive scores within $[-(C-h),C-h]$. Strictly inactive coordinates remain feasible by continuity. Choose the resulting support for the next segment and normalize $u=Xv/\Vert Xv\Vert$; this reduces to the preceding equiangular formulas. It handles ties, including a newly tied zero regression coefficient whose unconstrained direction would point against its correlation sign, without perturbing the data.

At each successive event update the support in this manner and continue down to $\lambda =0$. Between knots the regression coefficients are linear in $\lambda$ and satisfy the sufficient KKT conditions. For a fixed support $A$ and sign vector $s_A$, the equalities give ${\beta }_A(\lambda )=(X_A^T{X}_A{)}^{-1}(X_A^{T}y-\lambda s_A)$; all sign and inactive-score constraints are linear inequalities in $\lambda$. Each support/sign pattern therefore has an interval of validity, and there are only finitely many such patterns, so the complete path has finitely many linear segments. Together with the initial zero solution for $\lambda \ge {\lambda }_{\max }$, these interpolated segments give