岭估计的定义及性质¶
岭估计提出的想法很自然。当 $x_{1}, x_{2}, \cdots, x_{m}$ 间存在多重共线性关系时, 也就是 $\left|X^{T} X\right| \approx 0$, 我们设想给 $X^{T} X$ 加上一个正常数矩阵 $k I(k>0)$, 那么 $\left(X^{T} X+k I\right)$ 接近奇异的可能性要比 $X^{T} X$ 接近奇异的可能性小得多, 因此用 $$ \hat{\boldsymbol{\beta}}(\boldsymbol{k})=\left(X^{T} X+\boldsymbol{k} \boldsymbol{I}\right)^{-1} X^{T} \boldsymbol{Y} $$ 作为 $\beta$ 的估计应该比最小二乘估计稳定。为此, 我们给出岭估计的如下定义。
定义1 设 $0 \leq k<+\infty$, 满足 (12.16) 式的 $\hat{\beta}(k)$ 称为 $\beta$ 的岭估计。由 $\beta$ 的岭估计建立的回驲方程称为岭回归方程。其中 $k$ 称为岭参数。对于回归 系数 $\hat{\boldsymbol{\beta}}(k)=\left[b_{0}(k), b_{1}(k), \cdots, b_{m}(k)\right]^{T}$ 的分量 $b_{j}(k)(j \geq 1)$ 来说, 在直角坐标系 $\left(k, b_{j}(k)\right)$ 的图像是 $m$ 条曲线, 称为岭迹。 当 $k=0$ 时, $\hat{\beta}(0)$ 即为原来的最小二乘估计。
下面介绍岭估计的一些重要性质。
性质1 岭估计不再是无偏估计, 即 $E(\hat{\boldsymbol{\beta}}(k)) \neq \boldsymbol{\beta}_{\text {。 }}$
性质2 岭估计是压缩估计, 即 $\|\hat{\boldsymbol{\beta}}(\boldsymbol{k})\| \leq\|\hat{\boldsymbol{\beta}}\|$ 。
岭参数 $k$ 的选择¶
岭参数 $k$ 的选择有很多种方法, 我们这里只介绍两种容易理解的方法。
（1）岭迹法
观察岭迹曲线, 原则上应该选取使 $\hat{\boldsymbol{\beta}}(\boldsymbol{k})$ 稳定的最小 $k$ 值, 同时残差平 方和也不增加太多。
（2）均方误差法
岭估计的均方误差 $\operatorname{mse}(\hat{\beta}(k))=E\|\hat{\boldsymbol{\beta}}(k)-\boldsymbol{\beta}\|^{2}$ 是 $k$ 的函数, 可以证明它能 在某处取得最小值。计算并观察 $m \operatorname{se}(\hat{\boldsymbol{\beta}}(\boldsymbol{k}))$, 开始它将下降, 到达最小值后 开始上升。取它最小处的 $k$ 作为岭参数。
岭回归的应用¶
$$\begin{array}{ccccc|ccccc} \hline \text { 年份 } & x_{1} & x_{2} & x_{3} & y & \text { 年份 } & x_{1} & x_{2} & x_{3} & y \\ \hline 1949 & 149.3 & 4.2 & 108.1 & 15.9 & 1955 & 202.1 & 2.1 & 146.0 & 22.7 \\ 1950 & 171.5 & 4.1 & 114.8 & 16.4 & 1956 & 212.4 & 5.6 & 154.1 & 26.5 \\ 1951 & 175.5 & 3.1 & 123.2 & 19.0 & 1957 & 226.1 & 5.0 & 162.3 & 28.1 \\ 1952 & 180.8 & 3.1 & 126.9 & 19.1 & 1958 & 231.9 & 5.1 & 164.3 & 27.6 \\ 1953 & 190.7 & 1.1 & 132.1 & 18.8 & 1959 & 239.0 & 0.7 & 167.6 & 26.3 \\ 1954 & 202.1 & 2.2 & 137.7 & 20.4 & & & & & \\ \hline \end{array}$$下表 是 Malinvand 于 1966 年提出的研究法国经济问题的 一组数据。所考虑的因变量为进口总额 $y$, 三个解释变量分别为: 国内总产 值 $x_{1}$, 储存量 $x_{2}$, 总消费量 $x_{3}$ (单位均为 10 亿法郎)。
求岭回驲方程。
求解¶
画出的岭迹图如图12.1所示, 从岭迹图可以看出选 $\boldsymbol{k}=\mathbf{0 . 4}$ 较好。 对应的标准化岭回归方程为
$$
\hat{y}^{*}=0.2518 x_{1}^{*}+0.2262 x_{2}^{*}+0.6902 x_{3}^{*},
$$
将标准化回驲方程还原后得
$$
\hat{y}=-9.5320+0.0410 x_{1}+0.6231 x_{2}+0.1520 x_{3} \text {. }
$$
拟合优度 $R^{2}=0.9843$
代码
import numpy as np; import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge, RidgeCV
from scipy.stats import zscore
#plt.rc('text', usetex=True) #没装LaTeX宏包把该句注释
a=np.loadtxt("data/economic.txt")
n=a.shape[1]-1 #自变量的总个数
aa=zscore(a) #数据标准化
x=aa[:,:n]; y=aa[:,n] #提出自变量和因变量观测值矩阵
b=[] #用于存储回归系数的空列表
kk=np.logspace(-4,0,100) #循环迭代的不同k值
for k in kk:
md=Ridge(alpha=k).fit(x,y)
b.append(md.coef_)
st=['s-r','*-k','p-b'] #下面画图的控制字符串
for i in range(3):
plt.plot(kk,np.array(b)[:,i],st[i]);
plt.legend([r'$x_1$',r'$x_2$',r'$x_3$'],fontsize=15);plt.show()
mdcv=RidgeCV(alphas=np.logspace(-4,0,100)).fit(x,y);
print("最优alpha=",mdcv.alpha_)
#md0=Ridge(mdcv.alpha_).fit(x,y) #构建并拟合模型
md0=Ridge(0.4).fit(x,y) #构建并拟合模型
cs0=md0.coef_ #提出标准化数据的回归系数b1,b2,b3
print("标准化数据的所有回归系数为：",cs0)
mu=np.mean(a,axis=0); s=np.std(a,axis=0,ddof=1) #计算所有指标的均值和标准差
params=[mu[-1]-s[-1]*sum(cs0*mu[:-1]/s[:-1]),s[-1]*cs0/s[:-1]]
print("原数据的回归系数为：",params)
print("拟合优度：",md0.score(x,y))
在程序中我们使用了函数 RidgeCV 确定最佳 $k$ 值, 但拟合出 的 $x_{1}$ 系数是负的, 最终还是主观确定 $k=0.4$ 。
