Past exam of the mathematics course of the University of Cambridge / 2023 / iii / Paper 218 (Statistical Learning in Practice)

For individual $i$, let $Y_i$ be the reported count, $g_i\in \{0,1\}$ the gender indicator, and $m_i\in \{0,1\}$ the minority indicator. The fitted Poisson regression isIts log-likelihood isThe maximum-likelihood estimator isHolding minority status fixed, changing the gender indicator from zero to one multiplies the fitted conditional expected value by $e^{-0.1916}=0.826$. Thus the fitted mean count for men is about $17.4\%$ lower than that for women with the same minority status.

The researcher computed the Pearson chi-squared statisticand compared it with a ${\chi }_{1305}^2$ distribution, using $1308-3=1305$ residual degrees of freedom. Under an adequate large-sample Poisson regression, $X_P^2$ should be roughly the residual degrees of freedom. The reported tail probability rounds numerically to zero and gives strong evidence of overdispersion.

Possible causes include unobserved heterogeneity or omitted covariates, dependence among respondents, excess zeros, or an incorrect mean function. The conclusion that the Poisson variance assumption fails is well supported, although this test alone does not identify the cause or establish that the Quasi-Poisson regression variance $\mathrm{Var}(Y_i)=\varphi {\mu }_i$ is correct.

The reported Pearson dispersion estimator is $\hat{\varphi }=X_P^2/1305=1.719346$, soup to rounding.

With $X$ having rows $x_i^T=(1,g_i,m_i)$ and ${\mu }_i(\beta )=e^{x_i^T\beta }$, the quasi-score equation isIt is the same coefficient equation as for Poisson maximum likelihood, explaining why the two models have identical coefficient estimates.

Let $W=\mathrm{diag}({\hat{\mu }}_1,\dots ,{\hat{\mu }}_n)$ and let coordinate $2$ denote gender. The model-based standard error is

The Quasi-Poisson standard errors are trustworthy only if observations are independent, the log-linear mean is correct, and the variance is proportional to the mean with one common dispersion. Dependence, zero inflation, or covariate-dependent dispersion can invalidate this covariance formula.

A parametric bootstrap under model 1 proceeds as follows. Fit the Poisson model once and retain $X$ and the fitted means ${\hat{\mu }}_i$. For bootstrap repetition $b$, independently drawrefit the same Poisson regression to $(X,Y^{\ast (b)})$, and save its gender estimate ${\hat{\beta }}_1^{\ast (b)}$. The sample standard deviation of these estimates over many repetitions estimates the model-1 standard error. This bootstrap deliberately measures uncertainty under the fitted Poisson model; it does not repair real overdispersion unless the resampling model is enlarged to represent its cause.

The generalized linear mixed model tries to explain overdispersion by replacing the fixed minority coefficient with a Gaussian random intercept. Conditional on the group effect $b_m$,

This is a poor use of a random effect because minority has only two levels. Two realized intercepts contain almost no information about a random-effects distribution or its variance, and the two levels are substantively fixed categories rather than a sample from a population of groups. The model also has worse Akaike information criterion than model 1, $1132.8>1122.3$, and its fit does not establish that the original overdispersion has disappeared.

Because models 1 and 3 are full likelihood models for the same response, one can compare their Akaike information criterion values; that favors model 1. One can also compare held-out count prediction by K-fold cross-validation, using a common loss such as Poisson deviance or negative log predictive density.

The AIC comparison cannot include model 2 because a Quasi-Poisson fit specifies only mean and variance and has no full likelihood. Cross-validation can compare model 2 with model 3 if all predictions are scored by the same proper out-of-sample loss.

A process $(X_t{)}_{t\in \mathbb{Z}}$ is weakly stationary when it has finite second moments, a time-independent mean $\mathbb{E}{X}_t=\mu$, and an autocovariance functionthat depends only on the lag.

At lag $h$, the plot shows the sample autocorrelation functionUnder a white noise process, each fixed nonzero-lag sample autocorrelation is approximately $N(0,1/n)$, so the dashed pointwise $95\%$ reference lines are approximately $\pm 1.96/\sqrt{n}$.

The first nonzero-lag bar is well above the upper line, which contradicts the zero autocorrelation expected from white noise. Since the plot then largely cuts off, an moving-average process of order one is a plausible model; with the sampling interval as the time unit this is an $\mathrm{ARMA}(0,1)$ model.

The selected zero-mean autoregressive moving-average process isThe reported maximum-likelihood estimates areUsing the displayed asymptotic standard error gives the Wald confidence intervalThis normal interval is unreliable and likely too narrow because the series has only about twenty observations, the moving-average estimate is near the noninvertibility boundary $\theta =1$, and the same data were used to select the model. The finite-sample likelihood is consequently skewed and model-selection uncertainty is omitted.

The autoregressive process of order one is causal exactly whenbecause then $X_t={\sum }_{j\ge 0}{\varphi }^j{\epsilon }_{t-j}$ converges in mean square. Its autocovariance is

Under the intended assumption that the two white-noise sequences are mutually uncorrelated at every pair of times, $X$ and $W$ are uncorrelated. Their sum is therefore weakly stationary with

Strictly, the printed condition $\mathbb{E}[{\epsilon }_t{W}_t]=0$ only at equal times is insufficient. For example, $W_t=(-1{)}^t{\epsilon }_{t-1}$ is itself white noise and is contemporaneously uncorrelated with ${\epsilon }_t$, but the cross-covariance contribution can depend on $t$. The displayed answer therefore uses the standard intended cross-series white-noise assumption $\mathbb{E}[{\epsilon }_t{W}_s]=0$ for all $s,t$.

Using $X_t-\varphi X_{t-1}={\epsilon }_t$ givesUnder the cross-series uncorrelatedness used in part ii, terms in $U_t$ and $U_{t-h}$ involve disjoint white-noise times whenever $|h|>1$. HenceFor completeness,

Part iii and the stated characterization imply that $U$ has a moving-average process of order one representation $U_t={\eta }_t+\theta {\eta }_{t-1}$. Sincewhere $B$ is the backshift operator,Thus $Y$ is a causal autoregressive moving-average process of order $(1,1)$.

Let $T_i>0$ be record time and let $x_i=(1,c_i,d_i{)}^T$ contain standardized climb and distance. Model 1 is the normal linear modelModel 2 applies the same model after a logarithmic transformation:so $T_i$ is conditionally log-normal. Model 3 is a Gamma regression with logarithmic link:

Model 1's residuals have a systematic curved pattern and a spread that grows strongly with fitted time. This indicates an incorrect linear mean on the original scale and heteroscedasticity; several observations are also influential or outlying. The constant-variance assumption is therefore implausible.

Logging time greatly stabilizes the spread and removes most of the mean pattern, so model 2 is much more compatible with constant conditional variance and linearity. A few conspicuous residuals remain. A residual-versus-fitted plot alone does not check independence or fully establish normality.

For a likelihood $L(\theta )$ with $k$ estimated parameters, the Akaike information criterion iswhere $\hat{\theta }$ is the maximum-likelihood estimator. Among likelihoods for the same observed response and reference measure, smaller AIC estimates smaller expected out-of-sample Kullback-Leibler divergence up to a model-independent constant.

The printed values appear to favor model 2 because $-12.529<336.381$. They are not directly comparable: model 2 reports the Gaussian likelihood of $Z_i=\log T_i$, whereas model 3 reports a density for $T_i$. The change-of-variables formula for a probability density givesso the transformed model's AIC on the original response scale isBecause the standardized predictors have zero sample means and the ordinary-least-squares residuals sum to zero,Thereforewhich is slightly worse than model 3's $336.381$.

The ordinary least squares equations for model 2 areFor the Gamma generalized linear model, $V(\mu )={\mu }^2$ and $d\mu /d\eta =\mu$, so its score equation under the logarithmic link isWhen $T_i$ is close to ${\mu }_i$,by the first-order Taylor expansion of the exponential function. The Gamma score equations then become the model-2 normal equations, so their coefficient estimates are close.

Model 3 directly specifies the conditional mean and variance of the positive response on its observed scale. Its coefficients give multiplicative effects on mean record time, its prediction intervals concern time itself, and its likelihood can be compared directly with other original-scale models.

Model 2 instead models the mean of $\log T$. Back-transformation does not give the mean time without a retransformation bias correction; for a log-normal model, $\mathbb{E}(T\mid x)=\exp (x^T\alpha +{\tau }^2/2)$. The Gamma model also represents the increasing original-scale variance seen in the diagnostic plot, so it is preferable here.

With $f=1_{\{\mathrm{font}=\mathrm{serif}\}}$ and $d=1_{\{\mathrm{display}=\mathrm{popup}\}}$, the model-matrix input isThe feedforward neural network has three inputs, a fully connected layer of two ReLU units, and a fully connected two-class softmax output. Algebraically,where $W\in {\mathbb{R}}^{2\times 3}$, $b\in {\mathbb{R}}^2$, $V\in {\mathbb{R}}^{2\times 2}$, and $c\in {\mathbb{R}}^2$. The coding is $0$ for no click and $1$ for a click. There areparameters.

The fitted classifier isand training minimizes the empirical categorical cross-entropy loss

The softmax non-identifiability already proves that the coefficient vector is not unique. For any $a\in {\mathbb{R}}^2$ and $r\in \mathbb{R}$, replace both rows of $V$ by $V_{k\cdot }+a^T$ and both output biases by $c_k+r$. Every logit then gains the same value $a^{T}h+r$, so all softmax probabilities, classifications, and losses remain unchanged. Permuting the two hidden units supplies another non-uniqueness.

For the one-layer softmax fit, write its class logits as $z_k=a_k+b_k^{T}x$. ThenThus an equivalent logistic regression classifier predicts a click exactly whenUsing one reference class removes the common-logit-shift non-identifiability. Subject to the usual full-rank and no-separation conditions, the logistic parameter is identifiable. It uses four effective parameters rather than an eight-parameter redundant softmax representation, has a convex loss, and gives directly interpretable log-odds coefficients, so it is preferable for this binary linear classifier.

The first observation has $x=(1,0,0{)}^T$ and one-hot label $y=(0,1{)}^T$ in the output order $(\mathrm{no},\mathrm{yes})$. If every kernel weight and bias initially equals one, each hidden preactivation is two, so $h=(2,2{)}^T$. Both logits equal five and $p=(1/2,1/2{)}^T$.

For stochastic gradient descent on one cross-entropy observation,Hence the output-weight gradient has first row $(1,1)$ and second row $(-1,-1)$, while the output-bias gradient is $(1/2,-1/2)$. With learning rate one,Using the old output weights for backpropagation givesso every entry of $W$ and $b$ remains equal to one after this batch.

For fixed ${\lambda }_2$, letting ${\lambda }_1\downarrow 0$ gives ridge regression, including ordinary least squares when ${\lambda }_2=0$; letting ${\lambda }_1\to \infty$ forces every coefficient to zero. For fixed ${\lambda }_1$, letting ${\lambda }_2\downarrow 0$ gives the Lasso. Letting both penalties vanish gives an ordinary least squares solution, unique when $X$ has full column rank and otherwise potentially nonunique or path-dependent.

When ${\lambda }_2>0$, the term ${\lambda }_2\Vert \beta {\Vert }_2^2$ is strictly convex. Its sum with the convex squared loss and ${\ell }^1$ penalty is strictly convex and coercive, so the elastic net solution exists and is unique.

For coordinate $j$, hold all other coefficients fixed and form the partial residualThe one-coordinate coordinate descent problem isIts exact update isCyclically update coordinates and their residuals until the objective or coefficients converge. Convexity makes every limit point a global minimizer, and ${\lambda }_2>0$ makes it the unique minimizer.

When $X^{T}X=I_p$, expanding the objective separates it by coordinates:The soft thresholding solution is therefore

Model m1 is ridge regression, while m2 is the Lasso. Ridge shrinks but normally retains every coefficient; the Lasso's ${\ell }^1$ penalty sets many coefficients exactly to zero. Model m3 combines sparsity with the grouping effect of the elastic net: correlated predictors tend to enter together and receive more similar coefficients.

Accordingly, m3 keeps the weak variables age, lcp, and gleason at zero as m2 does, but retains lweight, lbph, svi, and pgg45 as the ridge fit does. The correlation $0.54$ between svi and lcavol explains why m3 keeps both with substantial coefficients, whereas m2 selects lcavol and discards svi. This lies between the dense ridge behavior and the more aggressively sparse Lasso behavior.

Choose $({\lambda }_1,{\lambda }_2)$ on a grid by K-fold cross-validation, comparing the same held-out prediction loss and optionally applying the one-standard-error rule for a simpler model. A genuinely untouched test set can then estimate final prediction error.

Ordinary model-based intervals after selecting nonzero coefficients ignore selection and are generally invalid. Valid approaches include Debiased Lasso or a selective-inference procedure under its assumptions, sample splitting followed by an unpenalized refit and inference on the independent half, or a bootstrap that repeats both tuning and fitting and is interpreted with care near the nonsmooth zero threshold.

The binary regression functions are the conditional class probabilitiesUnder zero-one loss, the Bayes classifier iswith arbitrary tie breaking. Its Bayes decision boundary is

Let $N_k(x)$ index the $k$ closest training covariates to $x$. The K-nearest neighbors algorithm estimatesand predicts one when this average is at least $1/2$.

Small $k$ gives low smoothing bias but high sampling variance and a jagged decision boundary. Larger $k$ averages more labels, reducing variance and producing a smoother boundary, but it mixes increasingly distant covariates and raises bias. The optimal balance depends on sample size, dimension, and smoothness of ${\eta }_1$, and is commonly selected by cross-validation.

For $\varphi (x)=x$ and centered data, the objective is the Rayleigh quotientThis is the first-direction optimization in principal component analysis. Therefore $\hat{u}$ is any unit eigenvector of the sample covariance matrix $S$ corresponding to its largest eigenvalue.

Write $\Phi \alpha ={\sum }_i{\alpha }_i\varphi (X_i)$. By the definition of the kernel matrix,andThe irrelevant positive factor $1/n$ gives exactly the stated constrained optimization.

Let $K{v}_1={\lambda }_1{v}_1$, where ${\lambda }_1>0$ is the largest eigenvalue and $\Vert v_1{\Vert }_2=1$. The generalized Rayleigh quotient is maximized byup to sign and addition of a vector in $\ker K$, which does not change $\hat{u}=\Phi \hat{\alpha }$. Repeated leading eigenvectors give further principal directions.

Kernel principal component analysis diagonalizes the centered $n$ by $n$ kernel matrix instead of an explicit covariance operator in a possibly infinite-dimensional feature space. A new point has coordinateso the leading coordinates require only evaluations of the positive-semidefinite kernel.

Running the K-nearest neighbors algorithm in a truncated collection of these coordinates can remove low-variance noise, reduce effective dimension, and allow a nonlinear boundary in the original covariates. This is the kernel trick: every feature-space inner product needed for fitting and projection is replaced by $k(x,z)=⟨\varphi (x),\varphi (z){⟩}_{\mathcal{H}}$ without constructing $\varphi (x)$ explicitly.