{"word_count": 615, "character_count": 6825}
Cochran Theorem - from The distribution of quadratic forms in a normal system, with applications to the analysis of covariance published in 1934 - is probably the most import one in a regression course. It is an application of a nice result on quadratic forms of Gaussian vectors. More precisely, we can prove that if [latex]\boldsymbol{Y}\sim\mathcal{N}(\boldsymbol{0},\mathbb{I}_d)[/latex] is a random vector with [latex]d[/latex] [latex]\mathcal{N}(0,1)[/latex] variable then (i) if [latex]A[/latex] is a (squared) idempotent matrix [latex]\boldsymbol{Y}^\top A\boldsymbol{Y}\sim\chi^2_r[/latex] where [latex]r[/latex] is the rank of matrix [latex]A[/latex], and (ii) conversely, if [latex]\boldsymbol{Y}^\top A\boldsymbol{Y}\sim\chi^2_r[/latex] then [latex]A[/latex] is an idempotent matrix of rank [latex]r[/latex]. And just in case, [latex]A[/latex] is an idempotent matrix means that [latex]A^2=A[/latex], and a lot of results can be derived (for instance on the eigenvalues). The prof of that result (at least the (i) part) is nice: we diagonlize matrix [latex]A[/latex], so that [latex]A=P\Delta P^\top[/latex], with [latex]P[/latex] orthonormal. Since [latex]A[/latex] is an idempotent matrix observe that[latex display="true"]A^2=P\Delta P^\top=P\Delta P^\top=P\Delta^2 P^\top[/latex]where [latex]\Delta[/latex] is some diagonal matrix such that [latex]\Delta^2=\Delta[/latex], so terms on the diagonal of [latex]\Delta[/latex] are either [latex]0[/latex] or [latex]1[/latex]'s. And because the rank of [latex]A[/latex] (and [latex]\Delta[/latex]) is [latex]r[/latex] then there should be [latex]r[/latex] [latex]1[/latex]'s and [latex]d-r[/latex] [latex]1[/latex]'s. Now write[latex display="true"]\boldsymbol{Y}^\top A\boldsymbol{Y}=\boldsymbol{Y}^\top P\Delta P^\top\boldsymbol{Y}=\boldsymbol{Z}^\top \Delta\boldsymbol{Z}[/latex]where [latex]\boldsymbol{Z}=P^\top\boldsymbol{Y}[/latex] that satisfies[latex]\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0},PP^\top)[/latex] i.e. [latex]\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0},\mathbb{I}_d)[/latex]. Thus [latex display="true"]\boldsymbol{Z}^\top \Delta\boldsymbol{Z}=\sum_{i:\Delta_{i,i}-1}Z_i^2\sim\chi^2_r[/latex]Nice, isn't it. And there is more (that will be strongly connected actually to Cochran theorem). Let [latex]A=A_1+\dots+A_k[/latex], then the two following statements are equivalent (i) [latex]A[/latex] is idempotent and [latex]\text{rank}(A)=\text{rank}(A_1)+\dots+\text{rank}(A_k)[/latex] (ii) [latex]A_i[/latex]'s are idempotents, [latex]A_iA_j=0[/latex] for all [latex]i\neq j[/latex].Now, let us talk about projections. Let [latex]\boldsymbol{y}[/latex] be a vector in [latex]\mathbb{R}^n[/latex]. Its projection on the space [latex]\mathcal V(\boldsymbol{v}_1,\dots,\boldsymbol{v}_p)[/latex] (generated by those [latex]p[/latex] vectors) is the vector [latex]\hat{\boldsymbol{y}}=\boldsymbol{V} \hat{\boldsymbol{a}}[/latex] that minimizes [latex]\|\boldsymbol{y} -\boldsymbol{V} \boldsymbol{a}\|[/latex] (in [latex]\boldsymbol{a}[/latex]). The solution is[latex display="true"]\hat{\boldsymbol{a}}=( \boldsymbol{V}^\top \boldsymbol{V})^{-1} \boldsymbol{V}^\top \boldsymbol{y} \text{ and } \hat{\boldsymbol{y}} = \boldsymbol{V} \hat{\boldsymbol{a}}[/latex]Matrix [latex] P=\boldsymbol{V} ( \boldsymbol{V}^\top \boldsymbol{V})^{-1} \boldsymbol{V}^\top[/latex] is the orthogonal projection on [latex]\{\boldsymbol{v}_1,\dots,\boldsymbol{v}_p\}[/latex] and [latex]\hat{\boldsymbol{y}} = P\boldsymbol{y}[/latex].Now we can recall Cochran theorem. Let [latex]\boldsymbol{Y}\sim\mathcal{N}(\boldsymbol{\mu},\sigma^2\mathbb{I}_d)[/latex] for some [latex]\sigma&#62;0[/latex] and [latex]\boldsymbol{\mu}[/latex]. Consider sub-vector orthogonal spaces [latex]F_1,\dots,F_m[/latex], with dimension [latex]d_i[/latex]. Let [latex]P_{F_i}[/latex] be the orthogonal projection matrix on [latex]F_i[/latex], then (i) vectors [latex]P_{F_1}\boldsymbol{X},\dots,P_{F_m}\boldsymbol{X}[/latex] are independent, with respective distribution [latex]\mathcal{N}(P_{F_i}\boldsymbol{\mu},\sigma^2\mathbb{I}_{d_i})[/latex] and (ii) random variables [latex]\|P_{F_i}(\boldsymbol{X}-\boldsymbol{\mu})\|^2/\sigma^2[/latex] are independent and [latex]\chi^2_{d_i}[/latex] distributed.We can try to visualize those results. For instance, the orthogonal projection of a random vector has a Gaussian distribution. Consider a two-dimensional Gaussian vectorlibrary(mnormt)r = .7s1 = 1s2 = 1Sig = matrix(c(s1^2,r*s1*s2,r*s1*s2,s2^2),2,2)SigY = rmnorm(n = 1000,mean=c(0,0),varcov = Sig)plot(Y,cex=.6)vu = seq(-4,4,length=101)vz = outer(vu,vu,function (x,y) dmnorm(cbind(x,y),mean=c(0,0), varcov = Sig))contour(vu,vu,vz,add=TRUE,col='blue')abline(a=0,b=2,col="red")Consider now the projection of points [latex]\boldsymbol{y}=(y_1,y_2)[/latex] on the straight linear with directional vector [latex]\overrightarrow{\boldsymbol{u}}[/latex] with slope [latex]a[/latex] (say [latex]a=2[/latex]). To get the projected point [latex]\boldsymbol{x}=(x_1,x_2)[/latex] recall that [latex]x_2=ay_1[/latex] and [latex]\overrightarrow{\boldsymbol{x},\boldsymbol{y}}\perp\overrightarrow{\boldsymbol{u}}[/latex]. Hence, the following code will give us the orthogonal projectionsp = function(a){x0=(Y[,1]+a*Y[,2])/(1+a^2)y0=a*x0cbind(x0,y0)}withP = p(2)for(i in 1:20) segments(Y[i,1],Y[i,2],P[i,1],P[i,2],lwd=4,col="red")points(P[,1],P[,2],col="red",cex=.7)Now, if we look at the distribution of points on that line, we get... a Gaussian distribution, as expected,z = sqrt(P[,1]^2+P[,2]^2)*c(-1,+1)[(P[,1]&#62;0)*1+1]vu = seq(-6,6,length=601)vv = dnorm(vu,mean(z),sd(z))hist(z,probability = TRUE,breaks = seq(-4,4,by=.25))lines(vu,vv,col="red")Or course, we can use the matrix representation to get the projection on [latex]\overrightarrow{\boldsymbol{u}}[/latex], or a normalized version of that vector actuallya=2U = c(1,a)/sqrt(a^2+1)U[1] 0.4472136 0.8944272matP = U %*% solve(t(U) %*% U) %*% t(U)matP %*% Y[1,][,1][1,] -0.1120555[2,] -0.2241110P[1,]x0         y0-0.1120555 -0.2241110 (which is consistent with our manual computation). Now, in Cochran theorem, we start with independent random variables,Y = rmnorm(n = 1000,mean=c(0,0),varcov = diag(c(1,1)))Then we consider the projection on [latex]\overrightarrow{\boldsymbol{u}}[/latex] and [latex]\overrightarrow{\boldsymbol{v}}=\overrightarrow{\boldsymbol{u}}^\perp[/latex]U = c(1,a)/sqrt(a^2+1)matP1 = U %*% solve(t(U) %*% U) %*% t(U)P1 = Y %*% matP1z1 = sqrt(P1[,1]^2+P1[,2]^2)*c(-1,+1)[(P1[,1]&#62;0)*1+1]V = c(a,-1)/sqrt(a^2+1)matP2 = V %*% solve(t(V) %*% V) %*% t(V)P2 = Y %*% matP2z2 = sqrt(P2[,1]^2+P2[,2]^2)*c(-1,+1)[(P2[,1]&#62;0)*1+1]We can plot those two projectionsplot(z1,z2)and observe that the two are indeed, independent Gaussian variables. And (of course) there squared norms are [latex]\chi^2_{1}[/latex] distributed.