Past exam of the mathematics course of the University of Cambridge / 2010 / iii / Paper 32 (Statistical Theory)

Put $G=X^{T}X$ and $H=X{G}^{-1}{X}^T$. The rank of a matrix assumption makes $G$ invertible, and $H$ is the orthogonal projection matrix onto the column space of the design matrix. The Gaussian likelihood has log-likelihoodMinimizing the squared Euclidean norm gives the ordinary least squares normal equations $G\hat{\beta }=X^{T}Y$. Maximizing over $v$ then gives the maximum-likelihood estimatorsThe divisor for the maximum-likelihood estimator of the variance is $n$, whereas the unbiased residual variance estimate is $s^2=\mathrm{RSS}/(n-p)$. The exceptional event $\mathrm{RSS}=0$ has probability zero under the Gaussian linear model with ${\sigma }^2>0$; on that event the likelihood is unbounded as $v\downarrow 0$.

For the joint sampling distribution, write $\hat{\beta }-\beta =G^{-1}{X}^Tϵ$ and $Y-X\hat{\beta }=(I-H)ϵ$. Their cross-covariance matrix vanishes because $X^T(I-H)=0$. These two vectors are jointly multivariate normal, so they are independent random variables. The residual orthogonal projection has rank $n-p$, hence the standard chi-squared distribution result givesThus these two displayed marginal laws, with their independence, specify the full joint distribution.

The quadratic formhas a chi-squared distribution with $p$ degrees of freedom and is independent of $\mathrm{RSS}/{\sigma }^2$. Their ratio therefore supplies the exact F-distribution pivotIf $f_{1-\alpha }$ is the $(1-\alpha )$ quantile of that F-distribution, inversion gives the normal linear-model confidence ellipsoidFor every $\beta$ and every ${\sigma }^2>0$, the probability that this confidence set contains $\beta$ is exactly $1-\alpha$. To test the null hypothesis, reject precisely when ${\beta }_0$ lies outside this confidence ellipsoid, equivalently when $T({\beta }_0)>f_{1-\alpha }$. Under the null hypothesis the rejection probability is exactly $\alpha$, independent of the unknown variance.

For the numerical calculation, the second column of the design matrix sums to zero and its sum of squares is $60$. Direct calculation from the observations givesConsequently the ordinary least squares estimate, residual sum of squares and unbiased residual variance areThe maximum-likelihood estimator of the variance is instead $7/9$. At the proposed null value the F-distribution statistic isTherefore reject the null hypothesis at the 5% level. The relevant critical value has $2$ numerator and $7$ residual degrees of freedom.

One sufficient maximum likelihood estimation statistical consistency theorem is the following. Let the observations be independent random variables with a common density $f_{\theta }$ relative to a fixed measure, let $\Theta$ be a compact metric parameter set containing the true value ${\theta }_0$, and suppose that the expected log-likelihoodis finite, continuous, and uniquely maximized at ${\theta }_0$. Assume also the uniform law of large numbersand assume that a measurable global maximum-likelihood estimator exists. Then every such maximum-likelihood estimator converges in probability to the true parameter. A useful sufficient condition for this uniform law of large numbers is almost-sure continuity of $\theta \mapsto \log f_{\theta }(Y_1)$ on the compact parameter set and an integrable envelope ${\sup }_{\theta }|\log f_{\theta }(Y_1)|\le M(Y_1)$, ${\mathbb{E}}_{{\theta }_0}M(Y_1)<\infty$. Under these conditions, identifiability supplies uniqueness through the Kullback-Leibler divergence identitywith equality only at $\theta ={\theta }_0$.

To see the statistical consistency conclusion, fix $\epsilon >0$. Compactness and the unique maximum give a strictly positive separation gapwhenever the set outside the neighbourhood is nonempty. If ${\Delta }_n$ denotes the supremum discrepancy in the uniform law of large numbers, maximization of the empirical log-likelihood yieldsHence the probability of $d({\hat{\theta }}_n,{\theta }_0)\ge \epsilon$ is at most $\mathbb{P}(2{\Delta }_n\ge {\delta }_{\epsilon })$, which tends to zero. This is the argmin consistency under uniform convergence in probability argument applied to the negative log-likelihood. These are sufficient conditions, not necessary ones; parameter-dependent support can require a direct argument, and independent observations need not have a common distribution in other models.

The likelihood function for the shifted exponential distribution isIt increases strictly up to the smallest observation and is zero beyond it. Thus shifted exponential maximum likelihood givesFor $t\ge 0$, independence and the survival function of the exponential distribution give the exact order statistic lawIn particular, ${\hat{\theta }}_n\ge \theta$ almost surely and ${\mathbb{P}}_{\theta }(|{\hat{\theta }}_n-\theta |>\epsilon )=e^{-n\epsilon }\to 0$. This proves statistical consistency. The exact endpoint limit for a shifted exponential distribution is stronger than an asymptotic approximation:The limit is one-sided and nonnormal; the moving support endpoint explains why a regular square-root-$n$ asymptotic normality of a maximum likelihood estimator theorem does not apply here.

Set $S_n={\sum }_{i=1}^n{Y}_i$ and $a_n={\sum }_{i=1}^n{2}^{-i}=1-2^{-n}$. For the Poisson distribution observations, factors independent of $\theta$ can be removed from the likelihood function, leavingWhen $S_n>0$, the log-likelihood derivative is $S_n/\theta -a_n$ and its second derivative is negative. The unique positive maximum-likelihood estimator is thereforeThere is a genuine boundary qualification: if $S_n=0$, the likelihood decreases strictly on $\theta >0$, so no maximum is attained in the stated open parameter space. On the closure $\theta \ge 0$, its maximizer is $0$. Defining the usual extended maximum-likelihood estimator by the same displayed formula for every $S_n$ allows its statistical consistency to be examined; the all-zero event cannot be ignored because its probability does not tend to zero.

Indeed, sums of independent Poisson distribution variables give $S_n\sim \mathrm{Poi}(\theta a_n)$, soThus the extended maximum-likelihood estimator is inconsistent. More precisely, $S_n$ increases to $S_{\infty }={\sum }_{i\ge 1}{Y}_i$, and $\mathbb{E}{S}_{\infty }=\theta <\infty$ shows that this total count is finite almost surely. Its distribution is Poisson distribution with mean $\theta$, by taking limits of the finite-sum laws. HenceThis finite-exposure Poisson inconsistency leaves a nondegenerate random limit even with infinitely many observations. Correspondingly, the total Fisher information is $a_n/\theta \to 1/\theta$, rather than diverging. The common-distribution hypothesis in the introductory maximum likelihood estimation theorem is absent here.

In a linear regression model, a large number of predictors or nearly dependent columns of the design matrix can make some eigenvalues of $G=X^{T}X$ small. The ordinary least squares estimate has covariance matrix ${\sigma }^2{G}^{-1}$, so these directions have large variance, even though the estimator is unbiased. Ridge regression reduces this variance by shrinking the estimate towards zero, at the cost of introducing bias. This bias-variance tradeoff can improve the mean squared error when the signal is not too large in the shrunken directions.

For the usual preprocessing with a regression intercept, separate the constant column from the predictors and leave the intercept unpenalized. Subtract the sample mean from the response and from each nonconstant predictor, then divide each centered predictor column by its positive scale, for example its Euclidean norm. Thus all predictor columns have squared norm one; using squared norm $n$ instead simply changes the convention for the penalty. Centering is performed on the slope columns, not on the constant column. A full-rank design containing an intercept retains full rank on its centered slope columns. A model without an intercept should retain its specified mean structure unless an intercept is deliberately introduced.

Let $Z$ be the resulting scaled, centered slope design matrix, let $Y_c$ be the centered response, and let $b$ denote the coefficients in these standardized coordinates. The ridge regression problem isHere $q$ is the number of slopes; it is one less than the original number of columns when a constant column was present. The gradient of this objective is $2(Z^{T}Z+\lambda I)b-2Z^T{Y}_c$. Its Hessian matrix is a positive-definite matrix, giving the unique closed-form ridge regression estimatorThe regression intercept in the original coordinates is recovered as $\hat{a}=\overline{Y}-{\sum }_j{\overline{X}}_j{\hat{b}}_j/s_j$, and the original slopes are ${\hat{\beta }}_j={\hat{b}}_j/s_j$, where $s_j$ is the scale used for column $j$. This is the unpenalized intercept in ridge regression convention. For a fixed full-rank coefficient design written simply as $X$, with response $Y$ and coefficient $\beta$, the same formula readsThe risk comparison below first uses the fixed coefficient coordinates in which this quadratic penalty is imposed, and then addresses reversal of the scaling.

Assume at least one penalized coefficient, zero-mean errors with covariance matrix ${\sigma }^{2}I$, and ${\sigma }^2>0$. No normal distribution assumption is needed for this mean squared error calculation. Put $G=X^{T}X$ and $A_{\lambda }=(G+\lambda I{)}^{-1}$. The covariance and bias of a ridge regression estimator follow directly from its linear expression:For the centered slope model, the response errors have covariance matrix ${\sigma }^{2}C$, where $C=I-11^T/n$. Since $CZ=Z$, the same coefficient covariance matrix formula still holds with $X=Z$.

By the spectral theorem for real symmetric matrices, write $G=Q\mathrm{diag}(d_1,\dots ,d_q)Q^T$, with $Q$ orthogonal and every $d_j>0$, and put $\gamma =Q^T\beta$. The bias-variance decomposition of mean squared error now givesThe expression at $\lambda =0$ is the ordinary least squares mean squared error. Its right derivative isTherefore, for every fixed coefficient vector, all sufficiently small positive penalties strictly improve mean squared error.

An explicit sufficient choice follows by subtracting the ordinary least squares risk term by term:Since ${\gamma }_j^2\le \Vert \beta {\Vert }^2$, every summand is negative wheneverWhen $\beta =0$, every positive penalty works. This is a pointwise result: the permitted penalty may depend on the true coefficient vector. A fixed positive penalty cannot improve risk for all arbitrarily large signals, because its squared bias grows quadratically with the signal. The same distinction appears in uniform directional risk improvement by ridge regression and unbounded directional risk of fixed ridge shrinkage. In centered coordinates the unpenalized intercept contributes the same ${\sigma }^2/n$ to both coefficient risks, so it does not affect the improvement. Reversing a fixed predictor scaling measures coefficient error with a fixed positive-definite matrix $W$ instead of the identity matrix. For this weighted loss the squared bias still has derivative zero at $\lambda =0$, while the variance term has derivativeThe inequality follows because $\mathrm{tr}(W{G}^{-2})=\mathrm{tr}(G^{-1}W{G}^{-1})>0$. Thus sufficiently small penalties also improve the original-coordinate mean squared error, though the explicit sufficient bound above was for the unweighted coordinates. Recovering an uncentered intercept adds a fixed nonnegative quadratic form to the slope loss, and its sample-mean error is uncorrelated with the centered slope errors, so the same weighted argument applies. A model with only an unpenalized intercept gives identical estimators. If the noise variance were zero, the strict improvement assertion would fail; positive noise variance is essential.

Let $R$ be the total number of rejected null hypotheses, and let $V$ count rejected null hypotheses that are in fact true. The false discovery proportion and false discovery rate areThus the false discovery proportion is zero when there are no rejections, and the false discovery rate averages the realized proportion of erroneous rejections.

For the Benjamini-Hochberg procedure, order the p-values as $P_{(1)}\le \cdots \le P_{(m)}$ and defineReject every hypothesis whose p-value satisfies $P_i\le \alpha R/m$ when $R>0$, and reject none otherwise. There are exactly $R$ such hypotheses: if a further p-value also met the threshold, its larger rank would satisfy its own, still larger, threshold, contradicting maximality. Under the stated independence and uniform true-null p-values, this procedure has false discovery rate $\alpha m_0/m\le \alpha$.

For the first modified Benjamini-Hochberg procedure, remove $P_1$, order the remaining $m-1$ p-values as $Q_{(1)},\dots ,Q_{(m-1)}$, and use the shifted critical values $\alpha (j+1)/m$. Define its rejection count byEquivalently, insert a zero in place of $P_1$ and run the ordinary Benjamini-Hochberg procedure on $m$ values; its total rejection count is $R^{(1)}=R_{-1}+1$. For $r=1,\dots ,m$, the Benjamini-Hochberg leave-one-out identity in event form isTo prove it, suppose first that the left event holds. Then $P_1$ is among the rejected values. Lowering it to zero leaves all ordered p-values at ranks greater than $r$ unchanged, so no new rank can meet its threshold. Rank $r$ still meets its threshold, giving $R^{(1)}=r$. Conversely, if $R^{(1)}=r$ and $P_1\le \alpha r/m$, the $r-1$ rejected remaining values and $P_1$ all meet the original rank-$r$ threshold. Hence the original count is at least $r$. Lowering a p-value cannot decrease the rejection count, so the original count is also at most $R^{(1)}=r$. This proves equality of the events, including possible ties.

For the second modified Benjamini-Hochberg procedure, remove both $P_1,P_2$, order the remaining p-values as $T_{(j)}$, and use shifted critical values $\alpha (j+2)/m$. Its count isReplacing $P_1,P_2$ by zeros in the full procedure gives a total of $R^{(12)}=R_{-12}+2$ rejections. The same argument gives the Benjamini-Hochberg leave-two-out identity, for $r=2,\dots ,m$:In the forward direction both replaced values already lie among the first $r$, so ordered ranks greater than $r$ remain unchanged. In the reverse direction reinserting the two values below the rank-$r$ threshold still leaves at least $r$ qualifying values, while monotonicity under lowering bounds the original count above by $r$. For $r=1$, the event of two p-values meeting the displayed threshold while only one is rejected is empty. For $m=1$ only the first modification is needed.

To calculate the second moment, write $I_i=1\{i\text{ rejected}\}$ for a true null. Expanding the square givesFor a true-null index $i$, let $R_{-i}$ denote the first modified count after removing $P_i$. It is a function only of the remaining p-values, and so is independent of the uniform variable $P_i$. The first event identity therefore yieldsFor distinct true-null indices $i,j$, the two independent uniform p-values are independent of the second modified count $R_{-ij}$. The second event identity givesThe last sum is one because $R_{-ij}$ always takes a value in $\{0,\dots ,m-2\}$. The ordered pairs in the expansion of $V^2$ number $m_0(m_0-1)$, with no factor of one half.

The true-null p-values are identically distributed and independent, so permutations among their indices leave the law of the full collection unchanged. Thus every $R_{-i}$ with $i\le m_0$ has the same distribution as $R_{-1}$, even though the false-null p-values need not be identically distributed. Combining the diagonal and off-diagonal contributions proves the second moment of the Benjamini-Hochberg false discovery proportion:If the first modified count is instead defined to include the inserted zero, the same answer is $A=1/R^{(1)}$. When $m_0=0$, both sides are zero; when $m_0=1$, there are no off-diagonal terms. Finally, repeating the first calculation with denominator $R\vee 1$ rather than its square gives $\mathbb{E}[I_i/(R\vee 1)]=\alpha /m$, which verifies the claimed exact false discovery rate.