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

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

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)

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

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

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

 

	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

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

	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

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

Out[8]:

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

#封装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

rSquare(aba, ridgeRegres)

Out[15]:

0.5274036413294183

rSquare(aba, standRegres)

Out[14]:

0.5276299399919837

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

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)

ridge.coef_#返回各行系数

Out[17]:

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

ridge.intercept_#返回截距

Out[18]:

3.0265413865908872

绘制岭迹图

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

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

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算法

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)

ridge.coef_#返回各行系数

Out[23]:

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

ridge.intercept_#返回截距

Out[24]:

3.0265413865908872