Statistical Analysis

Start with the Statistical Model

Least squares becomes a statistical estimator only after we state $how the observations were generated$ . Let

y=X\beta+\varepsilon

In this model, $X\in\mathbb{R}^{m\times p}$ is a fixed design matrix, $\beta\in\mathbb{R}^p$ is the unknown parameter vector, and $y\in\mathbb{R}^m$ contains the observations. For ordinary Gaussian least squares, assume

\varepsilon\sim\mathcal N(0,\sigma^2I_m),\qquad \sigma^2>0

The second assumption is that $X$ has full column rank. Under these two assumptions, the observation vector $y$ follows $\mathcal N(X\beta,\sigma^2I_m)$ . Its mean is $X\beta$ , and its covariance matrix is $\sigma^2I_m$ .

Fisher Information and Estimator Covariance

Up to a constant that does not depend on $\beta$ , the log-likelihood is

\ell(\beta)=-\frac1{2\sigma^2}\lVert y-X\beta\rVert_2^2

Its score, Hessian, and Fisher information are

\nabla\ell(\beta)=\frac1{\sigma^2}X^T(y-X\beta)

$X^TX$ is the Fisher information only in the standardized case $\sigma^2=1$ . Solving the score equation gives

\hat\beta=(X^TX)^{-1}X^Ty

Because a linear transformation of a Gaussian vector is Gaussian,

\hat\beta\sim\mathcal N\left(\beta,\sigma^2(X^TX)^{-1}\right)

The covariance is $exactly the inverse Fisher information$ under this model:

\operatorname{Cov}(\hat\beta)=\sigma^2(X^TX)^{-1}=\mathcal I(\beta)^{-1}

The diagonal entry $c_{jj}$ is the variance of $\hat\beta_j$ . An off-diagonal entry $c_{jk}$ measures how two coefficient estimates vary together, and their correlation is

\rho_{jk}=\frac{c_{jk}}{\sqrt{c_{jj}c_{kk}}}

Confidence Intervals When the Variance Is Unknown

If $\sigma^2$ is unknown and $m>p$ , estimate it from the residuals:

s^2=\frac{\lVert y-X\hat\beta\rVert_2^2}{m-p}

Under the Gaussian model, an exact two-sided $100(1-\alpha)\%$ confidence interval for coefficient $\beta_j$ is

\hat\beta_j\pm t_{m-p,\,1-\alpha/2}\,s\sqrt{[(X^TX)^{-1}]_{jj}}

The interval is a $repeated-sampling statement$ about the procedure. It is not the probability that a fixed parameter lies in one already computed interval.

Stable Computation from an Orthogonal Decomposition

Let $X=QR$ be the reduced $\mathrm{QR}$ decomposition, where $Q$ has orthonormal columns and $R$ is nonsingular upper triangular. Then

\hat\beta=R^{-1}Q^Ty

In numerical code, $solve the triangular system$ $R\hat\beta=Q^Ty$ directly. This avoids forming $R^{-1}$ . $\mathrm{QR}$ also avoids squaring the condition number through $X^TX$ .

Weighted and Generalized Least Squares

Suppose instead that the errors have known positive-definite covariance $\Omega$ :

\varepsilon\sim\mathcal N(0,\Omega)

Generalized least squares minimizes

(y-X\beta)^T\Omega^{-1}(y-X\beta)

The generalized least-squares solution, Fisher information, and estimator covariance are

\hat\beta_{\mathrm{GLS}}=(X^T\Omega^{-1}X)^{-1}X^T\Omega^{-1}y

For independent errors with variances $\sigma_i^2$ , use $\Omega=\operatorname{diag}(\sigma_1^2,\ldots,\sigma_m^2)$ . The objective becomes

\sum_{i=1}^m\frac{(y_i-x_i^T\beta)^2}{\sigma_i^2}

Each observation $i$ receives $weight$ $1/\sigma_i^2$ . Whitening with $W=\Omega^{-1/2}$ converts the problem to ordinary least squares for $\widetilde X=WX$ and $\widetilde y=Wy$ .

Fitting a Weighted Linear Model

Consider three observations for an intercept and slope:

X=\begin{pmatrix}1&0\\1&1\\1&2\end{pmatrix},\qquad y=\begin{pmatrix}1\\2\\5\end{pmatrix}

The third observation has standard deviation $2$ , so it receives only $one quarter of the weight$ of either of the first two. The weighted normal equations are

X^T\Omega^{-1}X=\begin{pmatrix}\frac94&\frac32\\\frac32&2\end{pmatrix}

Writing the two equations and eliminating the intercept makes the parameter calculation explicit:

\begin{aligned}9\hat\beta_0+6\hat\beta_1&=17\\3\hat\beta_0+4\hat\beta_1&=9\\6\hat\beta_1&=27-17=10\\\hat\beta_1&=\frac53,\qquad\hat\beta_0=\frac{9-4(5/3)}3=\frac79\end{aligned}

Observation	Standard deviation	Inverse-variance weight
1	1	1
2	1	1
3	2	$1/4$

Solving the weighted normal equations gives

\hat\beta_{\mathrm{GLS}}=\begin{pmatrix}\frac79\\\frac53\end{pmatrix}

The weighting has a probabilistic meaning only if the variance model is defensible. Incorrect weights can produce misleading standard errors, and exact Gaussian confidence intervals also require the linear mean, covariance, and normal-error assumptions to be adequate.