Bias in the MLE of variance component in a multivariate Gaussian?

Given an $n$-vector $y$ (responses) and a design matrix $X$, I wish to fit them with a simple linear regression model $$y=X\beta+e,$$ where $e\sim\mathcal{N}(0, \sigma^2I)$. Then, we have $$y\sim\mathcal{N}(X\beta, \sigma^2I).$$ Then the maximum likelihood estimations (MLE) of $\beta$ and $\sigma^2$ are just $$\hat\beta=(X^TX)^{-1}X^Ty$$ and $$\hat{\sigma^2}=\frac{(y-X\hat\beta)^T(y-X\hat\beta)}{n}.$$

I understand that $\sigma^2$'s estimator is biased in that $E\{\hat{\sigma^2}\}\neq\sigma^2$.

On Page 4 of this lecture notes, the author claims $E\{\hat{\sigma^2}\}=\frac{n-r}{n}\sigma^2$ without specifying what $r$ is. I guess $r$ is the DoF loss while estimating $\beta$, i.e., the dimension of $\beta$. But how do we derive it?

I tried to prove it myself but got confused at "over which distribution are we taking the expectation?"

Yes, $r$ is indeed the length/dimension of $\beta$. Define the 'hat' matrix $H=X(X^TX)^{-1}X^T$ so that $Hy=X\hat\beta=\hat y$. Then $E((y-X\hat\beta)^T(y-X\hat\beta))=E((y-Hy)^T(y-Hy))=E(y^T(I-H)^T(I-H)y)$. Now $I-H$ is an orthogonal projection so $(I-H)^T(I-H)=I-H$ and hence the expectation above reduces to $E(y^Ty)-E((Hy)^T(Hy))$. Given that $y\sim\mathcal{N}(X\beta, \sigma^2I)$, it follows that $E(y^Ty)=\sigma^2trace(I)+(X\beta)^T(X\beta)=n\sigma^2+(X\beta)^T(X\beta)$. It also follows that the transformed variable $Hy\sim\mathcal{N}(HX\beta,\sigma^2HH^T)=\mathcal{N}(X\beta,\sigma^2H)$. Hence $E((Hy)^T(Hy))=\sigma^2trace(H)+(X\beta)^T(X\beta)$. Because $H$ is an orthogonal projection, its trace is equal to its rank, which equals the rank of $X$, which is just $dim(\beta)$. So, combining the obtained equations for $E(y^Ty)$ and $E((Hy)^T(Hy))$ leads to $E((y-X\hat\beta)^T(y-X\hat\beta))=(n-dim(\beta))\sigma^2$.

• (+1) It took me so long to realize how awesome this answer is. Now only one thing confuses me: I obtained $E(y^Ty)$ from "sum of $N$ independent standard Gaussians is a chi-square distribution of DoF $N$". Here $N$ happens to be $tr(I)$. However, how can we directly know we can compute $E((Hy)^T(Hy))$ by simply replacing $tr(I)$ with $tr(H)$? Because to me, $\sigma^2H$ may not even be diagonal, implying the Gaussians may not be independent! In this case, we no longer have chi-square distribution! So I presume $\sigma^2H$ is diagonal? If so, why? If not, how did you compute that? Thanks a lot Sep 25 '15 at 11:03
• @SibbsGambling Thanks! The matrix $H$ is not necessarily diagonal. Rather than trying to characterise the distribution of $(Hy)^T(Hy)$, you can compute $E((Hy)^T(Hy))$ directly, using the mean and covariance matrix of $Hy$. Writing $z=Hy$, $E((Hy)^T(Hy))=E(z^T z)=\sum_i E[z_i^2]$. Now use the fact that $E[z_i^2]=var(z_i)+E(z_i)^2=\sigma^2 H_{ii}+(X\beta)_i^2$. Sep 26 '15 at 16:08
• This is beautiful! Thanks! I also found another way of computing it: Theorem 6.4 in nyu.edu/econ/dept/courses/peracchi/gstat6.pdf. But definitely your solution is simpler! Sep 26 '15 at 16:18

Just make our life easier, suppose $X$ is just one random variable and your model has no intercept then the MLE of $\sigma^2$ is $\hat{\sigma^2}=\frac{\sum_{i=1}^n (x_i-\hat{\mu})^2}{n}$ (This is just a simpler case for $\hat{\sigma^2}=\frac{(y-X\hat\beta)^T(y-X\hat\beta)}{n}$ and $\hat\mu=X\hat\beta=\bar{X}$)

We know $\hat{\sigma^2}=\frac{\sum_{i=1}^n (x_i-\hat{\mu})^2}{n}$ is a biased estimator of $\sigma^2$, the unbiased one should be $\frac{\sum_{i=1}^n (x_i-\hat{\mu})^2}{n-1}$

Therefore, $E(\frac{n}{n-1} \hat{\sigma^2})=E(\frac{n}{n-1}*\frac{\sum_{i=1}^n (x_i-\hat{\mu})^2}{n})=\sigma^2$

$\Rightarrow E(\hat{\sigma^2)}=\frac{n-1}{n}\sigma^2$

you can see $r=1$ here.

For multivariate situation I guess it should be the same, hope someone else can help.