线性回归误差项方差的估计

摘要

之前在文章线性回归系数的几个性质 中，我们证明了线性回归系数项的几个性质。在这篇短文中，我们介绍线性回归模型中误差项方差的估计。

线性回归误差项概念的回顾

我们先来回忆一下什么是线性回归的误差项。在文章线性回归系数的几个性质 中，我们指出，对于单变量的线性回归模型，$$y_i={\beta }_1{x}_i+{\beta }_0+{ϵ}_i,i=1,2,\cdots ,n$$，
其中 ${ϵ}_i$ 为误差项 (error)，每个 ${ϵ}_i$ 均是一个随机变量，独立且都服从 一个均值为0, 方差为 ${\sigma }^2$ 的概率分布。这篇短文介绍的就是对这个 ${\sigma }^2$ 的估计。

残差平方和 (residual sum of squares)

我们定义 $S{S}_{\mathrm{Res}}=\sum _{i=1}^n{e}_i^2=\sum _{i=1}^n(y_i-\widehat{y_i}{)}^2$。

首先，我们证明，
$$S{S}_{\mathrm{Res}}=\sum _{i=1}^n{y}_i^2-n{\bar{y}}^2-\widehat{{\beta }_1}{S}_{xy}$$。

我们把 $\widehat{y_i}=\widehat{{\beta }_0}+\widehat{{\beta }_1}\cdot x_i$ 代入 $S{S}_{\mathrm{Res}}$ 的表达式。

我们有，

$\begin{array}{rl}S{S}_{\mathrm{Res}}& =\sum _{i=1}^n(y_i-\widehat{y_i}{)}^2=\sum _{i=1}^n(y_i-(\widehat{{\beta }_0}+\widehat{{\beta }_1}\cdot x_i){)}^2\\ & =\sum _{i=1}^n{y}_i^2-2\sum _{i=1}^n{y}_i(\widehat{{\beta }_0}+\widehat{{\beta }_1}\cdot x_i)+\sum _{i=1}^n(\widehat{{\beta }_0}+\widehat{{\beta }_1}\cdot x_i{)}^2\\ & =\sum _{i=1}^n{y}_i^2-2\widehat{{\beta }_0}\cdot n\bar{y}-2\widehat{{\beta }_1}\sum _{i=1}^n{x}_i{y}_i+n{\widehat{{\beta }_0}}^2+2\widehat{{\beta }_0}\widehat{{\beta }_1}\sum _{i=1}^n{x}_i+{\widehat{{\beta }_1}}^2\sum _{i=1}^n{x}_i^2\\ & =\sum _{i=1}^n{y}_i^2-2(\bar{y}-\widehat{{\beta }_1}\bar{x})n\bar{y}-2\widehat{{\beta }_1}\sum _{i=1}^n{x}_i{y}_i+n(\bar{y}-\widehat{{\beta }_1}\bar{x}{)}^2+2(\bar{y}-\widehat{{\beta }_1}\bar{x})\widehat{{\beta }_1}\sum _{i=1}^n{x}_i+{\widehat{{\beta }_1}}^2\sum _{i=1}^n{x}_i^2\\ & =\sum _{i=1}^n{y}_i^2-2n{\bar{y}}^2+2n\widehat{{\beta }_1}\bar{x}\bar{y}-2\widehat{{\beta }_1}\sum _{i=1}^n{x}_i{y}_i+n{\bar{y}}^2-2n\widehat{{\beta }_1}\bar{x}\bar{y}+\\ & n{\widehat{{\beta }_1}}^2{\bar{x}}^2+2\bar{y}\widehat{{\beta }_1}\sum _{i=1}^n{x}_i-2{\widehat{{\beta }_1}}^2\bar{x}\sum _{i=1}^n{x}_i+{\widehat{{\beta }_1}}^2\sum _{i=1}^n{x}_i^2\\ & =\sum _{i=1}^n{y}_i^2-n{\bar{y}}^2-n{\widehat{{\beta }_1}}^2{\bar{x}}^2+{\widehat{{\beta }_1}}^2\sum _{i=1}^n{x}_i^2-2\widehat{{\beta }_1}\sum _{i=1}^n{x}_i{y}_i+2n\widehat{{\beta }_1}\bar{x}\bar{y}\\ & =\sum _{i=1}^n{y}_i^2-n{\bar{y}}^2+{\widehat{{\beta }_1}}^2(\sum _{i=1}^n{x}_i^2-n{\bar{x}}^2)-2\widehat{{\beta }_1}(\sum _{i=1}^n{x}_i{y}_i-n\bar{x}\bar{y})\\ & =\sum _{i=1}^n{y}_i^2-n{\bar{y}}^2+{\widehat{{\beta }_1}}^2{S}_{xx}-2\widehat{{\beta }_1}{S}_{xy}\\ & =\sum _{i=1}^n{y}_i^2-n{\bar{y}}^2+\widehat{{\beta }_1}\left(\frac{S_{xy}}{S_{xx}}\right)S_{xx}-2\widehat{{\beta }_1}{S}_{xy}\\ & =\sum _{i=1}^n{y}_i^2-n{\bar{y}}^2-\widehat{{\beta }_1}{S}_{xy}\end{array}$

于是，我们就证明了 $$S{S}_{\mathrm{Res}}=\sum _{i=1}^n{y}_i^2-n{\bar{y}}^2-\widehat{{\beta }_1}{S}_{xy}$$。

残差平方和的期望

下面我们来看残差平方和的期望，即 $\mathbb{E}[S{S}_{\mathrm{Res}}]$ 。

$\mathbb{E}[S{S}_{\mathrm{Res}}]=\mathbb{E}\left[\sum _{i=1}^n{y}_i^2\right]-n\mathbb{E}\left[{\bar{y}}^2\right]-\mathbb{E}\left[\widehat{{\beta }_1}{S}_{xy}\right]$。

我们分项来求。

我们先来求 $\mathbb{E}\left[\sum _{i=1}^n{y}_i^2\right]$。我们知道，
$\mathbb{E}\left[y_i\right]=\mathbb{E}[{\beta }_1{x}_i+{\beta }_0+{ϵ}_i]={\beta }_1{x}_i+{\beta }_0$。注意这里我们用到了 $\mathbb{E}\left[{ϵ}_i\right]=0$。
而 $\mathrm{Var}\left[{ϵ}_i\right]={\sigma }^2$。于是，我们知道 $\mathrm{Var}\left[y_i\right]={\sigma }^2$。从而，$\mathbb{E}\left[y_i^2\right]=\mathrm{Var}\left[y_i\right]+(\mathbb{E}[y_i]{)}^2=({\beta }_1{x}_i+{\beta }_0{)}^2+{\sigma }^2$。$\mathbb{E}\left[\sum _{i=1}^n{y}_i^2\right]=\sum _{i=1}^n(({\beta }_1{x}_i+{\beta }_0{)}^2+{\sigma }^2)=\sum _{i=1}^n({\beta }_1{x}_i+{\beta }_0{)}^2+n{\sigma }^2$。

对于 $\mathbb{E}\left[{\bar{y}}^2\right]$，我们采用一样的方法。因为 $\mathbb{E}\left[\bar{y}\right]=\frac{1}{n}\sum _{i=1}^n({\beta }_1{x}_i+{\beta }_0)={\beta }_0+{\beta }_1\bar{x}$。$\mathrm{Var}\left[\bar{y}\right]=\frac{1}{n^2}\sum _{i=1}^n\mathrm{Var}\left[y_i\right]=\frac{1}{n^2}\cdot n{\sigma }^2=\frac{{\sigma }^2}{n}$。
于是，$\mathbb{E}\left[{\bar{y}}^2\right]=\mathrm{Var}\left[\bar{y}\right]+\left(\mathbb{E}\right[\bar{y}]{)}^2=\frac{{\sigma }^2}{n}+({\beta }_0+{\beta }_1\bar{x}{)}^2$。

对于最后一项，$\mathbb{E}\left[\widehat{{\beta }_1}{S}_{xy}\right]$，因为 $\widehat{{\beta }_1}=\frac{S_{xy}}{S_{xx}}$，而 $S_{xx}$ 是常数，所以我们须要求 $\mathbb{E}\left[S_{xy}^2\right]$。
注意到，$S_{xy}=\sum _{i=1}^n(x_i-\bar{x})(y_i-\bar{y})=\sum _{i=1}^n(x_i-\bar{x})y_i$。

我们有 $\mathbb{E}[\sum _{i=1}^n(x_i-\bar{x})y_i]=\sum _{i=1}^n(x_i-\bar{x})\cdot \mathbb{E}\left[y_i\right]=\sum _{i=1}^n(x_i-\bar{x})({\beta }_0+{\beta }_1{x}_i)$。
另外，$\mathrm{Var}[\sum _{i=1}^n(x_i-\bar{x})y_i]=\sum _{i=1}^n(x_i-\bar{x}{)}^2\cdot \mathrm{Var}[y_i]=\sum _{i=1}^n(x_i-\bar{x}{)}^2{\sigma }^2$。这里我们用到了 $y_i,i=1,2,\cdots n$ 是非相关的 (uncorrelated)。
从而，
$\begin{array}{rl}\mathbb{E}\left[S_{xy}^2\right]& =\mathbb{E}\left[\right(\sum _{i=1}^n(x_i-\bar{x})y_i{)}^2]\\ & =\mathrm{Var}[\sum _{i=1}^n(x_i-\bar{x})y_i]+{\left(\mathbb{E}[\sum _{i=1}^n(x_i-\bar{x})y_i]\right)}^2\\ & =\sum _{i=1}^n(x_i-\bar{x}{)}^2{\sigma }^2+{\left(\sum _{i=1}^n(x_i-\bar{x})({\beta }_0+{\beta }_1{x}_i)\right)}^2\\ & =S_{xx}{\sigma }^2+{\left(\sum _{i=1}^n{\beta }_1{x}_i(x_i-\bar{x})\right)}^2\\ & =S_{xx}{\sigma }^2+{\beta }_1^2{S}_{xx}^2\end{array}$

从而，$\mathbb{E}\left[\widehat{{\beta }_1}{S}_{xy}\right]=\frac{S_{xx}{\sigma }^2+{\beta }_1^2{S}_{xx}^2}{S_{xx}}={\sigma }^2+{\beta }_1^2{S}_{xx}$。

我们把，$\mathbb{E}\left[\sum _{i=1}^n{y}_i^2\right],\mathbb{E}\left[{\bar{y}}^2\right],\left[S_{xy}^2\right]$ 这三项代入 $\mathbb{E}[S{S}_{\mathrm{Res}}]=\mathbb{E}\left[\sum _{i=1}^n{y}_i^2\right]-n\mathbb{E}\left[{\bar{y}}^2\right]-\mathbb{E}\left[\widehat{{\beta }_1}{S}_{xy}\right]$。
我们有
$\begin{array}{rl}\mathbb{E}[S{S}_{\mathrm{Res}}]& =\mathbb{E}\left[\sum _{i=1}^n{y}_i^2\right]-n\mathbb{E}\left[{\bar{y}}^2\right]-\mathbb{E}\left[\widehat{{\beta }_1}{S}_{xy}\right]\\ & =\sum _{i=1}^n({\beta }_1{x}_i+{\beta }_0{)}^2+n{\sigma }^2-n\left(\frac{{\sigma }^2}{n}+({\beta }_0+{\beta }_1\bar{x}{)}^2\right)-({\sigma }^2+{\beta }_1^2{S}_{xx})\\ & =(n-2){\sigma }^2\end{array}$

也就是说，$\frac{S{S}_{\mathrm{Res}}}{n-2}$ 是 ${\sigma }^2$ 的一个无偏估计。

实验验证

class sigmasqu_estimation:
    
    def __init__(self, arr_x: np.array, beta1: float, beta0: float, epsilon: float):
        #self.N = N
        self.X = arr_x
        self.beta1 = beta1
        self.beta0 = beta0
        self.epsilon = epsilon
        
        self.Sxx = ((self.X - self.X.mean()) ** 2).sum()
        self.X_bar = self.X.mean()
        
    def estimate_sigmasqu(self, N: int) -> tuple:
        res_sigmasqu_esti = []
        for i in range(N):
            #print(i)
            cur_error = np.random.normal(0, self.epsilon, arr_x.shape)
            cur_y = self.beta0 + self.beta1 * self.X + cur_error
            cur_y_bar = cur_y.mean()
            Sxy = ((self.X - self.X.mean()) * (cur_y - cur_y_bar)).sum()
            cur_beta1 = Sxy / self.Sxx
            cur_beta0 = cur_y_bar - cur_beta1 * self.X_bar
            cur_y_hat = cur_beta1 * self.X + cur_beta0 
            SS_res = ((cur_y - cur_y_hat) ** 2).sum()
            res_sigmasqu_esti.append(SS_res / (self.X.shape[0] - 2))
        return np.array(res_sigmasqu_esti)

arr_x = np.array(range(1, 11))
a = sigmasqu_estimation(arr_x, 2, 3, 1)
res = a.estimate_sigmasqu(10 ** 5)

res
np.mean(res)

1.0004110047596488

我们发现，$\frac{S{S}_{\mathrm{Res}}}{n-2}$ 的均值非常接近 ${\sigma }^2$

事实上，$(n-2)S{S}_{\mathrm{Res}}/{\sigma }^2$ 服从的是 ${\chi }_{n-2}^2$ 的分布 [1]。

plt.figure(figsize=(8, 6), dpi=100)
plt.hist(res, bins=50, density=True);
line_vert = [[1, c] for c in np.linspace(0, 1, 100)]
plt.plot([c[0] for c in line_vert], [c[1] for c in line_vert], '-', linewidth=4)
plt.xlabel("estimated $\sigma^2$ value", fontsize=20)
plt.ylabel("count", fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20);

参考文献

[1] Introduction to linear regression analysis, Douglas C. Montgomery, Elizabeth A. Peck, G. Geoffrey Vining, John Wiley & Sons (2021)