# Ridge Regression and Bayesian Regression 学习笔记

Ridge Regression

概念

• 数学形式及解读

给定数据集D=(x1,y1),(x2,y2),...,(xm,ym),其中x∈Rd,y∈R.我们考虑最简单的线性回归模型，以平方误差为损失函数，则优化目标为

minw∑i=1m(yi−wTxi)2

当样本特征很多，而样本数相对较少时，上式很容易陷入过拟合，为缓解过拟合问题，对上式引入正则化项，若使用L2范数正则化，则有

minw∑i=1m(yi−wTxi)2+λ||w||22

其中正则化参数λ>0。上式称为“岭回归”，通过引入正则化，确实能降低过拟合的风险。

• 算法适用范围

岭回归最先用来处理特征数多于样本数的情况，现在也用于在估计中加入偏差，从而得到更好的估计。这里通过引入λ来限制了所有w之和，通过引入该惩罚项，能够减少不重要的参数，这个技术在统计学中也叫缩减。

• 岭回归中的岭是什么？
  岭回归使用了单位矩阵乘以常量λ，我们观察其中的单位矩阵I，可以看到值I贯穿整个对角线，其余元素全是0，形象地，在0构成的平面上有一条1组成的“岭”，这就是岭回归中的“岭”的由来。

缩减方法可以去掉不重要的参数，因此能更好地理解数据。此外，与简单的线性回归相比，缩减法能取得更好的预测效果。

2. 算法代码

简单来说，岭回归就是在矩阵xTx上加一个λx从而使得矩阵非奇异，进而能对xTx+λx求逆。其中矩阵x是一个m*m的单位矩阵，对角线上元素全为1，其他元素全为0.而λ是一个用户定义的数值，后面会做介绍。在这种情况下，回归系数的计算公式将变成:

wˆ=(XTX+λI)−1XTy

公式推导

#岭回归def ridgeRegress(xMat,yMat,lam = 0.2) #lambda = 0.2    xTx = xMat.T*xMat    denom = xTx + eye(shape(xMat)[1]*lam    if linalg.set(denom) == 0.0:        print "This matrix id singular, cannot do inverse"        return    ws = denom.I*(xMat.T*yMat)    return wsdef ridgeTest(xArr,yArr):    xMat = mat(xArr)    yMat = mat(yArr).T    yMean = mean(yMat,0)    yMat = yMat - yMean    XMeans = mean(xMat,0)    xVar = var(xMat,0)    xMat = (xMat - xMeans)/xVar    #以上步骤为标准化    numTestPts = 30 #30个lambda    wMat = zeros((numTestPts,shape(xMat)[1]))    for i in range(numTestPts):        ws = ridgeRegress(xMat,yMat,lam = exp(i-10))        wMat[i:] = ws.T    return wMat    #这样可以看出lambda在取非常小的值时和取非常大的值时分别对结果造成的影响

在sklearn中的文档代码：

import numpy as npimport matplotlib.pyplot as pltfrom sklearn import linear_model# X is the 10x10 Hilbert matrix#生成数据X = 1. / (np.arange(1, 11) + np.arange(0, 10)[:, np.newaxis])y = np.ones(10)# print(X)# print(y)n_alphas = 200 #lambdaalphas = np.logspace(-10, -2, n_alphas) #从-10到-2两百个数clf = linear_model.Ridge(fit_intercept=False)#初定义，不需要截距coefs = []for a in alphas:    clf.set_params(alpha=a)    clf.fit(X, y)    coefs.append(clf.coef_)#10个因子plt.figure(figsize=(10,8))ax = plt.gca()#gca() returns the current axesax.plot(alphas, coefs)ax.set_xscale('log')#注意这一步，alpha是对数化了的ax.set_xlim(ax.get_xlim()[::-1])  # reverse axisplt.xlabel('alpha')plt.ylabel('weights')plt.title('Ridge coefficients as a function of the regularization')plt.axis('tight')#图边界的设置plt.show()

得到结果

Bayesian Regression

1.概念
贝叶斯回归
贝叶斯线性回归的引入主要是在最大似然估计很难决定模型的复杂程度，ridge回归加入的惩罚参数其实也是解决这个问题的，同时可以采用的方法还有对数据的进行正规化处理，另一个可以解决此问题的方法就是采用贝叶斯方法。

贝叶斯回归模型假设先验概率，似然函数和后验概率都是正态分布。先验概率是假设模型输出Y是符合均值为Xw的正态分布，正则化参数α被看作是一个需要从数据中估计得到的随机变量。

p(y|X,w,α)=N(y|Xw,α)

• 贝叶斯方法中，控制模型参数分布的参数，被称为超参数（可以理解为参数的参数，因为参数本身符合一个分布，而这个分布中的参数就是超参数）

贝叶斯回归的优点

• 它适应于手上拥有的数据
• 在估计的过程中可以加上正则化参数

贝叶斯回归的缺点

• 模型的参数估计耗费时间较长

  1. Bayesian Ridge Regression
     回归系数w的先验分布规律为球形正态分布，超参数为λ。我们需要通过最大化边际似然函数来估计超参数α和λ，以及回归系数w。
     
     p(w|λ)=N(w|0,λ−1Ip)
     
     α和λ的先验概率被选择为gamma分布.这些先验常常是无信息的，我们需要通过最大化边际似然函数来估计。

3.算法适用
如果我们的数据有很多缺失或者矛盾的病态数据，可以考虑BayesianRidge类，它对病态数据鲁棒性很高，也不用交叉验证选择超参数。但是极大化似然函数的推断过程比较耗时，一般情况不推荐使用。

• 计算机科学中，健壮性（英语：Robustness）是指一个计算机系统在执行过程中处理错误，以及算法在遭遇输入、运算等异常时继续正常运行的能力。

代码实现

import numpy as npimport matplotlib.pyplot as pltfrom scipy import statsfrom sklearn.linear_model import BayesianRidge, LinearRegression

Generating simulated data with Gaussian weights

np.random.seed(0)n_samples, n_features = 100, 100 #一共100个特征X = np.random.randn(n_samples, n_features)  # Create Gaussian data# Create weights with a precision lambda_ of 4.lambda_ = 4.w = np.zeros(n_features)# Only keep 10 weights of interestrelevant_features = np.random.randint(0, n_features, 10) #只选取10个特征有值for i in relevant_features:    w[i] = stats.norm.rvs(loc=0, scale=1. / np.sqrt(lambda_))# Create noise with a precision alpha of 50.alpha_ = 50.noise = stats.norm.rvs(loc=0, scale=1. / np.sqrt(alpha_), size=n_samples)# Create the targety = np.dot(X, w) + noise

Fit the Bayesian Ridge Regression and an ordinary least squares for comparison

clf = BayesianRidge(compute_score=True) #贝叶斯回归估计clf.fit(X, y)ols = LinearRegression() #线性估计ols.fit(X, y)