Exact Gaussian regression inference and prediction

Let $W=\mathrm{col}X$. The columns of $X$ are a basis of $W$, so $\dim W=p$. By [F1], ${\mathbb{R}}^n=W\oplus W^{\perp }$. Since $W^{\perp }$ is a subspace of the finite-dimensional space ${\mathbb{R}}^n$, [F4] makes it finite-dimensional; [F5] then gives $n=\dim W+\dim W^{\perp }$, so $\dim W^{\perp }=n-p$. By [F2], choose an orthonormal basis of $W^{\perp }$ and index its $n-p$ vectors as $z_{p+1},\dots ,z_n$. The columns of $A=[X{z}_{p+1}\cdots z_n]$ form a basis of ${\mathbb{R}}^n$: their $X$-part spans $W$, their remaining part spans $W^{\perp }$, and the direct sum in [F1] makes the union linearly independent. Thus $A$ is invertible.

Apply [F3] to write $A=UT$, where $U$ is orthogonal and $T$ is upper triangular with positive diagonal. Let $Q$ be the first $p$ columns of $U$, let $Q_0$ be its remaining columns, and let $R$ be the leading $p\times p$ block of $T$. Because $T_{ij}=0$ when $i>j$, every entry below row $p$ in its first $p$ columns is zero. The first $p$ columns of $A$ are therefore $X=QR$. The matrix $R$ is upper triangular with positive diagonal, while $U=[Q{Q}_0]$ gives $Q^{T}Q=I_p$ and the required orthonormal completion.

Since $R$ is invertible, $X^{T}X=R^{T}R$ and $(X^{T}X{)}^{-1}=R^{-1}{R}^{-T}$. Hence $P_X=QR(R^{T}R{)}^{-1}{R}^T{Q}^T=Q{Q}^T$, because $R(R^{T}R{)}^{-1}{R}^T=I_p$.

The Gaussian experiment has the fixed-design first moment in [F14], and [F15] defines its least-squares solution set. Apply [F1] and write $X=QR$, with $[Q{Q}_0]$ orthogonal. Since $R$ is invertible, $X^{T}X=R^{T}R$ and $(X^{T}X{)}^{-1}=R^{-1}{R}^{-T}$. Therefore $\widehat{\beta }=R^{-1}{Q}^{T}Y$ and $P_X=QR(R^{T}R{)}^{-1}{R}^T{Q}^T=Q{Q}^T$. For any $b\in {\mathbb{R}}^p$, orthogonality of the columns of $Q$ and $Q_0$ gives $\mathrm{\parallel }Y-Xb{\mathrm{\parallel }}^2=\mathrm{\parallel }{Q}^{T}Y-Rb{\mathrm{\parallel }}^2+\mathrm{\parallel }{Q}_0^{T}Y{\mathrm{\parallel }}^2$. The unique minimizer is $b=R^{-1}{Q}^{T}Y$, consistent with [F16]; the residual is $Q_0{Q}_0^{T}Y$.

Set $W_U=\mathrm{col}X$, $A=W_U\cap W_R^{\perp }$, and $B=W_U^{\perp }$, and let $P_A,P_B$ be the orthogonal projections onto $A,B$. By [F21], ${\mathbb{R}}^n=W_U\oplus B$. For any $u\in W_U$, write $u=P_{W_R}u+(u-P_{W_R}u)$. The first term belongs to $W_R\subseteq W_U$, while the second lies in $W_R^{\perp }$ by [F11] and in $W_U$ because it is a difference of vectors in $W_U$. Thus $u\in W_R+A$. Conversely $W_R,A\subseteq W_U$, and $W_R\cap A=\{0\}$ because $A\subseteq W_R^{\perp }$ and [F21] gives $W_R\cap W_R^{\perp }=\{0\}$, so $W_U=W_R\oplus A$. Combining the two decompositions gives ${\mathbb{R}}^n=W_R\oplus A\oplus B$, with pairwise orthogonal summands. By [F23], these subspaces are finite-dimensional. The columns of $X$ are a basis of $W_U$, so $\dim W_U=p$. If $q=p$, then $W_R=\{0\}$ and $\dim W_R=0=p-q$; if $q<p$, [F25] gives $\dim W_R=\mathrm{crank}(X_R)=\mathrm{rrank}(X_R)=\mathrm{rank}(X_R)=p-q$. Applying [F24] gives $\dim A=q$ and $\dim B=n-p=\nu$. By [F22], choose orthonormal bases of $W_R,A,B$, using the empty basis when $q=p$. Let $O$ be the matrix whose columns are these bases in that order; then $O$ is orthogonal. For every $y=r+a+b$ in the unique decomposition, [F11] gives $P_{W_R}y=r$, $P_{A}y=a$, $P_{B}y=b$, and $P_{W_U}y=r+a$. Hence $y-P_{W_R}y=P_{A}y+P_{B}y$ and $y-P_{W_U}y=P_{B}y$. Standardize the errors to $\eta =ϵ\mathrm{/}\sigma$. By [F9], each ${ϵ}_i\mathrm{/}\sigma$ is standard normal, and [F19] preserves independence under coordinatewise scaling, so $\eta$ has independent standard-normal coordinates. Write $O^T\eta =(Z_R,Z_A,Z_B)$. By [F2], these transformed coordinates are independent standard normals, so [F18] makes $Z_A$ and $Z_B$ independent. Under the null, $\mu =X\beta \in W_R$, so $P_A\mu =P_B\mu =0$. The two residual vectors above are therefore $P_Aϵ+P_Bϵ$ and $P_Bϵ$; since $P_X=P_{W_U}$ and the chosen columns are orthonormal, $RS{S}_R-RS{S}_F=\mathrm{\parallel }{P}_Aϵ{\mathrm{\parallel }}^2={\sigma }^2\mathrm{\parallel }{Z}_A{\mathrm{\parallel }}^2$ and $RS{S}_F=\mathrm{\parallel }{P}_Bϵ{\mathrm{\parallel }}^2={\sigma }^2\mathrm{\parallel }{Z}_B{\mathrm{\parallel }}^2$. Thus [F3] gives independent ${\chi }_q^2$ and ${\chi }_{\nu }^2$ variables; both degrees are positive, [F4] makes the denominator positive almost surely, and the ratio definition [F5] gives the stated $F_{q,\nu }$ law.

Put $\eta =ϵ\mathrm{/}\sigma$ and use $U=[Q{Q}_0]$ from [F1]. The density in [F9] becomes standard normal under this scaling, and [F19] preserves independence under coordinatewise scaling; hence $\eta$ has independent standard-normal coordinates. By [F2], $U^T\eta =(G,H)$ has independent standard-normal coordinates, and [F18] groups the first $p$ and last $\nu$ coordinates into independent vectors. Step 1.1 gives $\widehat{\beta }-\beta =\sigma R^{-1}G$ and $(I-P_X)Y=\sigma Q_{0}H$. By [F13], these linear and squared-norm maps are Borel, so [F19] preserves the independence of their values. From [F20], $\widehat{\beta }\sim N_p(\beta ,{\sigma }^2{R}^{-1}{R}^{-T})=N_p(\beta ,{\sigma }^2(X^{T}X{)}^{-1})$; [F3] gives $RS{S}_F\mathrm{/}{\sigma }^2=\mathrm{\parallel }H{\mathrm{\parallel }}^2\sim {\chi }_{\nu }^2$, and $s>0$ almost surely by [F4]. For $c\mathrm{\ne }0$, positive definiteness gives $v_c=c^T(X^{T}X{)}^{-1}c>0$. The linear image $Z_c=(R^{-T}c{)}^{T}G\mathrm{/}\sqrt{v_c}$ is $N(0,1)$ by [F20], since $\mathrm{\parallel }{R}^{-T}c{\mathrm{\parallel }}^2=v_c$. It is independent of $\mathrm{\parallel }H{\mathrm{\parallel }}^2$ by [F18] and [F19]. Therefore, almost surely, $T_c=Z_c\mathrm{/}\sqrt{\mathrm{\parallel }H{\mathrm{\parallel }}^2\mathrm{/}\nu }$ has law $t_{\nu }$ by [F3] and [F5]; the defined value on $s=0$ does not change its law. Finally, $\mathrm{row}(X)={\mathbb{R}}^p$, so [F17] makes every coefficient contrast estimable; explicitly, $X^T[X(X^{T}X{)}^{-1}c]=c$.

For fixed $v={\sigma }^2>0$, independence and the normal density [F9] make the joint likelihood as a function of $b$ proportional to $\exp (-\mathrm{\parallel }Y-Xb{\mathrm{\parallel }}^2\mathrm{/}(2v))$. Step 1.1 shows its unique maximizer is $\widehat{\beta }$.

Define $T_{c,{\theta }_0}$ by replacing $c^T\beta$ in the numerator of $T_c$ with ${\theta }_0$. This statistic is measurable by [F10] and [F13], since its piecewise formula uses Borel linear and norm maps. Under $H_0:c^T\beta ={\theta }_0$, step 2.1 gives the symmetric continuous $t_{\nu }$ law in [F6]; the two tails beyond $\pm t_{\nu ,1-\alpha \mathrm{/}2}$ each have probability $\alpha \mathrm{/}2$. The contrast null and alternative are both nonempty: $\beta ={\theta }_{0}c\mathrm{/}(c^{T}c)$ is in the null, and $\beta =({\theta }_0+1)c\mathrm{/}(c^{T}c)$ is in its alternative. The partial F statistic is Borel because $RS{S}_R$ and $RS{S}_F$ are continuous functions of $Y$, and its ratio on $\{RS{S}_F>0\}$ with the constant branch on $\{RS{S}_F=0\}$ is Borel by [F13], so the two threshold rejection events are Borel. Step 1.2 and the unique $F_{q,\nu }$ quantile in [F6] give upper-tail rejection probability $\alpha$. The F-test null contains $\beta =0$; step 1.2 gives a $q$-dimensional subspace $A\subseteq \mathrm{col}X$ with $q\ge 1$, so some $a\in A\setminus \{0\}$ equals $X\beta$ for a coefficient vector $\beta$. Since $A\cap W_R=\{0\}$, this mean is outside $W_R$, so that parameter is in the alternative. Each Borel rejection event $E$ defines a measurable deterministic kernel by $K(y,\{1\})={\mathbf{1}}_E(y)$ and $K(y,\{0\})=1-{\mathbf{1}}_E(y)$, so [F7] applies. Both tests reject with probability $\alpha$ at every null parameter, so each has size $\alpha$ by [F8].

Write $r=RS{S}_F$. At $b=\widehat{\beta }$, the joint log-likelihood for $v>0$ is $-\frac{n}{2}\log (2\pi v)-\frac{r}{2v}$. If $r>0$, its derivative $(r-nv)\mathrm{/}(2v^2)$ is positive for $v<r\mathrm{/}n$ and negative for $v>r\mathrm{/}n$, so the unique maximum is $v=r\mathrm{/}n$. Step 2.1 shows $r>0$ almost surely. If an observed vector has $r=0$, the likelihood tends to infinity as $v\downarrow 0$ and has no maximizer on $(0,\mathrm{\infty })$. Therefore the joint MLE is $(\widehat{\beta },RS{S}_F\mathrm{/}n)$ almost surely, while $s^2=RS{S}_F\mathrm{/}(n-p)$ differs because $n>p$ and $RS{S}_F>0$ almost surely. Together with steps 2.1 and 3.1 this proves all asserted distributional and test-calibration claims, with the likelihood exception stated for every data vector.

Fix $c\mathrm{\ne }0$. Then $v_c>0$ because $X^{T}X$ is positive definite. For each full parameter $\theta =(\beta ,{\sigma }^2)$, the statistic $q_{\theta }(Y)=(c^T\widehat{\beta }-c^T\beta )\mathrm{/}(s\sqrt{v_c})$, assigned value $0$ when $s=0$, is a pivot by [F1] and [F2], with common law $t_{\nu }$. By [F5], $B=[-t^{\ast },t^{\ast }]$ has probability $1-\alpha$. Define $D_c(Y)=\{({\beta }^{\mathrm{\prime }},{\sigma }^{\mathrm{\prime }2})\in \mathrm{\Theta }:c^T{\beta }^{\mathrm{\prime }}\in C_c(Y)\}$. For each fixed $({\beta }^{\mathrm{\prime }},{\sigma }^{\mathrm{\prime }2})\in \mathrm{\Theta }$, its membership event is $\mathrm{\mid }{c}^T\widehat{\beta }-c^T{\beta }^{\mathrm{\prime }}\mathrm{\mid }\le t^{\ast }s\sqrt{v_c}$, a closed, hence Borel, set because $\widehat{\beta }$ and $s$ are continuous functions of $Y$. Thus [F3] applies to $D_c$. Pivot inversion [F4] gives exact coverage for $D_c^{\mathrm{\prime }}(Y)=\{{\theta }^{\mathrm{\prime }}\in \mathrm{\Theta }:\mathrm{\mid }{q}_{{\theta }^{\mathrm{\prime }}}(Y)\mathrm{\mid }\le t^{\ast }\}$. On $s>0$, $D_c(Y)=D_c^{\mathrm{\prime }}(Y)$; the exceptional event $s=0$ has probability zero for every parameter by [F1]. Therefore $D_c$ has exact coverage $1-\alpha$, which is precisely ${\mathbb{P}}_{\beta ,{\sigma }^2}\{c^T\beta \in C_c(Y)\}=1-\alpha$.

If $x_0\mathrm{\ne }0$, the same argument with $c=x_0$ gives the displayed mean-response interval because $h_0=v_{x_0}$. If $x_0=0$, then $x_0^T\beta =x_0^T\widehat{\beta }=h_0=0$ for every parameter and data vector, so the interval is $\{0\}$ and contains the true mean with probability $1$.

Write $Z=(\widehat{\beta }-\beta )\mathrm{/}\sigma$ and $G_{\ast }={ϵ}_0^{\ast }\mathrm{/}\sigma$. By [F1] and [F8], $Z$ and $G_{\ast }$ are independent centered normal vectors/scalars with covariance $(X^{T}X{)}^{-1}$ and variance $1$. Their joint law is the product law by [F7], so [F6] gives $\frac{Y_0^{\ast }-x_0^T\widehat{\beta }}{\sigma \sqrt{1+h_0}}=\frac{G_{\ast }-x_0^{T}Z}{\sqrt{1+h_0}}\sim N(0,1)$, since the numerator has variance $1+x_0^T(X^{T}X{)}^{-1}{x}_0=1+h_0$. Because $G_{\ast }$ is independent of the training data, [F7] gives $\mathcal{L}(Z,RS{S}_F,G_{\ast })=\mathcal{L}(Z,RS{S}_F)\otimes \mathcal{L}(G_{\ast })$. Since $Z$ is a measurable affine function of $\widehat{\beta }$, [F1] and [F7] imply that $Z$ is independent of $RS{S}_F$; the joint law therefore factors as $\mathcal{L}(Z)\otimes \mathcal{L}(RS{S}_F)\otimes \mathcal{L}(G_{\ast })$, so $(Z,G_{\ast })$ is independent of $RS{S}_F$. The displayed standard normal is a measurable function of $(Z,G_{\ast })$, so it is independent of $RS{S}_F\mathrm{/}{\sigma }^2\sim {\chi }_{\nu }^2$ by [F7]. By [F9] and [F10], the ratio $(Y_0^{\ast }-x_0^T\widehat{\beta })\mathrm{/}(s\sqrt{1+h_0})$ has law $t_{\nu }$ with positive denominator almost surely. Its central acceptance interval has probability $1-\alpha$ by [F5], which is exactly the event $Y_0^{\ast }\in C_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}(Y;x_0)$.