A-岭回归的python'实现

虽然在线性回归的求解过程中，通过借助最⼩⼆乘法能够迅速找到全域最优解，但最⼩⼆乘法本身的使⽤条件较为苛刻，必须要求当 为满秩矩阵或正定阵时才可进⾏逆矩阵或⼴义逆矩阵的求解， 在实际应⽤中经常会遇⻅矩阵不存在逆矩阵或⼴义逆矩阵的情况，并且当X的各列存在线性相关关系 （即多重共线性）的时候，最⼩⼆乘法的求解结果不为⼀，这⾥需要注意的是，在进⾏数据采集过程 中，数据集各列是对同⼀个客观事物进⾏描述，很难避免多重共线性的存在，因此存在共线性是很多 数据集的⼀般情况。当然更为极端的情况则是数据集的列⽐⾏多，此时最⼩⼆乘法⽆法对其进⾏求 解。 总的来说，解决共线性的问题的⽅法主要有以下几种： 其⼀是在建模之前对各特征进⾏相关性检验，若存在多重共线性，则可考虑进⼀步对数据集进⾏ SVD分解或PCA主成分分析，在SVD或PCA执⾏的过程中会对数据集进⾏正交变换，最终所得数 据集各列将不存在任何相关性。当然此举会对数据集的结构进⾏改变，且各列特征变得不可解 释。 其⼆则是采⽤逐步回归的⽅法，以此选取对因变量解释⼒度最强的⾃变量，同时对于存在相关性 的⾃变量加上⼀个惩罚因⼦，削弱其对因变量的解释⼒度，当然该⽅法不能完全避免多重共线性 的存在，但能够绕过最⼩⼆乘法对公线性较为敏感的缺陷，构建线性回归模型； 其三则是在原有的算法基础上进⾏修改，放弃对线性⽅程参数⽆偏估计的苛刻条件，使其能够容 忍特征列存在多重共线性的情况，并且能够顺利建模，且尽可能的保证RSS取得最⼩值。 通常来说，能够利⽤⼀个算法解决的问题尽量不⽤多个算的组合来解决，因此此处我们主要考虑后两 个解决⽅案，特别第三种方案我们用到岭回归算法和Lasso算法。 本文来自CDA数据分析学院（吴昊天）老师

 

所谓岭回归就是在原有⽅程系数计算公式中添加了⼀个扰动项aI,原先⽆法求⼴义逆的情况变成可以求出其⼴义逆，使得问题稳定并得以求解，

其中 a是⾃定义参数， I则是单位矩阵

In [ ]:

w=(x.T * x +aI).I * X.T * y 

In [3]:

#对于对单位矩阵的⽣成，需要借助NumPy中的eye函数，该函数需要输⼊对⻆矩阵规模参数
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
np.eye(3)

Out[3]:

array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

In [4]:

def ridgeRegres(dataSet, lam=0.2):
    xMat = np.mat(dataSet.iloc[:, :-1].values)
    yMat = np.mat(dataSet.iloc[:, -1].values).T
    xTx = xMat.T * xMat
    denom = xTx + np.eye(xMat.shape[1])*lam
    ws = denom.I * (xMat.T * yMat)
    return ws

In [6]:

aba = pd.read_table('abalone.txt', header = None)#该数据集源于UCI，记录了鲍⻥的⽣物属性，⽬标字段是该⽣物的年龄
aba.head()

 

Out[6]:

	0	1	2	3	4	5	6	7	8
0	1	0.455	0.365	0.095	0.5140	0.2245	0.1010	0.150	15
1	1	0.350	0.265	0.090	0.2255	0.0995	0.0485	0.070	7
2	-1	0.530	0.420	0.135	0.6770	0.2565	0.1415	0.210	9
3	1	0.440	0.365	0.125	0.5160	0.2155	0.1140	0.155	10
4	0	0.330	0.255	0.080	0.2050	0.0895	0.0395	0.055	7

In [7]:

#其中最后⼀列为标签列，同时，由于数据集第⼀列是分类变量，且没有⽤于承载截距系数的列，因此
#此处将第⼀列的值全部修改为1，然后再进⾏建模
aba.iloc[:, 0] = 1
aba.head()

Out[7]:

	0	1	2	3	4	5	6	7	8
0	1	0.455	0.365	0.095	0.5140	0.2245	0.1010	0.150	15
1	1	0.350	0.265	0.090	0.2255	0.0995	0.0485	0.070	7
2	1	0.530	0.420	0.135	0.6770	0.2565	0.1415	0.210	9
3	1	0.440	0.365	0.125	0.5160	0.2155	0.1140	0.155	10
4	1	0.330	0.255	0.080	0.2050	0.0895	0.0395	0.055	7

In [8]:

rws = ridgeRegres(aba)
rws
#返回各列的系数同样的第一行为截距

Out[8]:

matrix([[  3.0178564 ],
        [ -0.02075972],
        [ 11.4007625 ],
        [ 11.02315972],
        [  8.71725704],
        [-19.64613377],
        [ -8.96257395],
        [  9.18180982]])

In [9]:

#封装R**2
def rSquare(dataSet, regres):#设置参数为 数据集 与 回归方法
    sse = sseCal(dataSet, regres) 
    y = dataSet.iloc[:, -1].values
    sst = np.power(y - y.mean(), 2).sum()
    return 1 - sse / sst

In [10]:

#封装最小二乘
def standRegres(dataSet):
    xMat = np.mat(dataSet.iloc[:, :-1].values)#将DataFrame转换成array 再将array转换成matrix(矩阵) 因为只有矩阵可以进行计算
    yMat = np.mat(dataSet.iloc[:, -1].values).T
    xTx = xMat.T*xMat
    if np.linalg.det(xTx) == 0:#判断xTx是否是满秩矩阵，若不满秩，则⽆法对其进⾏求逆矩阵的操作
        print("This matrix is singular, cannot do inverse")
        return
    ws = xTx.I * (xMat.T*yMat)
    return ws
#这⾥需要注意的是，当使⽤矩阵分解来求解多元线性回归时，必须添加⼀列全为1的列，⽤于表征线性⽅程截距b。

In [12]:

#将SSE做一个封装
def sseCal(dataSet, regres):#设置参数为 数据集 与 回归方法
    n = dataSet.shape[0] 
    y = dataSet.iloc[:, -1].values
    ws = regres(dataSet)
    yhat = dataSet.iloc[:, :-1].values * ws
    yhat = yhat.reshape([n,])
    rss = np.power(yhat - y, 2).sum()
    return rss

In [15]:

rSquare(aba, ridgeRegres)

Out[15]:

0.5274036413294183

In [14]:

rSquare(aba, standRegres)

Out[14]:

0.5276299399919837

 

调⽤Scikit—Learn中岭回归算法验证⼿动建模有效性

In [16]:

from sklearn import linear_model
ridge = linear_model.Ridge(alpha=.2)
ridge.fit(aba.iloc[:, :-1], aba.iloc[:, -1])

Out[16]:

Ridge(alpha=0.2, copy_X=True, fit_intercept=True, max_iter=None,
   normalize=False, random_state=None, solver='auto', tol=0.001)

In [17]:

ridge.coef_#返回各行系数

Out[17]:

array([  0.        ,  -0.04157241,  11.39836699,  11.01582006,
         8.71877234, -19.64361088,  -8.9576113 ,   9.1872849 ])

In [18]:

ridge.intercept_#返回截距

Out[18]:

3.0265413865908872

 

绘制岭迹图

 

绘制岭迹图基本函数基本思路为从⼩到⼤依次取 a，然后查看各特征列系数的衰减速度，从中查看是否存在公线性，以及 的最佳取值

In [19]:

def ridgeTest(dataSet):
    xMat = np.mat(dataSet.iloc[:, :-1].values)
    yMat = np.mat(dataSet.iloc[:, -1].values).T
    yMean = np.mean(yMat, axis = 0)
    yMat = yMat - yMean
    xMeans = np.mean(xMat, axis = 0) 
    xVar = np.var(xMat,axis = 0) 
    xMat = (xMat - xMeans)/xVar
    numTestPts = 30
    wMat = np.zeros((numTestPts,xMat.shape[1]))
    for i in range(numTestPts):
        ws = ridgeRegres(dataSet,np.exp(i-10))
        wMat[i,:]=ws.T
    return wMat

In [20]:

aba = pd.read_table('abalone.txt', header = None)
ridgeWeights = ridgeTest(aba)
plt.plot(ridgeWeights)
plt.xlabel('log(lambda)')
plt.ylabel('weights')

Out[20]:

Text(0,0.5,'weights')

 

 

把所有回归系数的岭迹都绘制在一张图上，如果这些曲线比较稳定，如上图所示，利用最小二乘估计会有一定的把握。

 

在Scikit-Learn中执⾏Lasso算法

In [22]:

from sklearn import linear_model
las = linear_model.Lasso(alpha = 0.01)
las.fit(aba.iloc[:, :-1], aba.iloc[:, -1])
 

Out[22]:

Lasso(alpha=0.01, copy_X=True, fit_intercept=True, max_iter=1000,
   normalize=False, positive=False, precompute=False, random_state=None,
   selection='cyclic', tol=0.0001, warm_start=False)

In [23]:

ridge.coef_#返回各行系数

Out[23]:

array([  0.        ,  -0.04157241,  11.39836699,  11.01582006,
         8.71877234, -19.64361088,  -8.9576113 ,   9.1872849 ])

In [24]:

ridge.intercept_#返回截距

Out[24]:

3.0265413865908872