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$.

  • $\begingroup$ (+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 $\endgroup$ Sep 25 '15 at 11:03
  • $\begingroup$ @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$. $\endgroup$ Sep 26 '15 at 16:08
  • $\begingroup$ 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! $\endgroup$ 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.