第5章
CHAPTER 5


线 性 回 归




线性模型是一类统计模型的总称,它包括了线性回归模型、方差分析模型、协方差分析模型和线性混合效应模型(或称方差分量模型)等。许多生物、医学、经济、管理、地质、气象、农业、工业、工程技术等领域的现象都可以用线性模型来近似描述。因此线性模型成为现代统计学中应用最为广泛的模型之一。
5.1概述
给定样本x,用列向量表示该样本x=(x(1),x(2),…,x(n))T。样本有n种特征,我们用x(i)表示样本x的第i个特征。线性模型(Linear model)的形式为: 


f(x)=w·x+b

其中w=(w(1),w(2),…,w(n))T为每个特征对应的权重生成的权重向量,称为权重向量,权重向量直观地表达了各个特征在预测中的重要性。
线性模型中的“线性”其实就是一系列一次特征的线性组合,在二维空间中是一条直线,在三维空间中是一个平面,然后推广到n维空间,这样可以理解为广义的线性模型。
线性模型非常简单,易于建模,应用广泛,主要包括岭回归、lasso回归、Elastic Net、逻辑回归、线性判别分析等。
5.2普通线性回归
线性回归是一种回归分析技术。回归分析本质上就是一个函数估计的问题(函数估计包括参数估计和非参数估计两类),就是找出因变量和自变量之间的因果关系。回归分析的因变量应该是连续变量,如果因变量为离散变量,则问题转化为分类问题,回归分析是一个有监督学习的问题。
5.2.1基本概述
给定数据集T={(x1,y1),(x2,y2),…,(xN,yN)},xi∈XRn,yi∈YR,i=1,2,…,N,其中xi=(x(1)i,x(2)i,…,x(n)i)T。需要学习的模型为: 


f(x)=w·x+b
也即根据已知的数据集T来计算参数w和b。




对于给定的样本x,其预测值为y^i=f(xi)=w·xi+b。采用平方损失函数,则在训练集T上,模型的损失函数为: 


L(f)=∑Ni=1(y^i-yi)2=∑Ni=1(w·xi+b-yi)2
我们的目标是损失函数的最小化,即: 


(w*,b*)= argminw,b∑Ni=1(w·xi+b-yi)2

可以采用梯度下降法来求解上述最优化问题的数值解。在使用梯度下降法时,要注意特征归一化(Feature Scaling)。
特征归一化有两个好处: 
(1) 提升模型的收敛速度,比如两个特征x1和x2,x1的取值为0~2000,而x2的取值为1~5,假如只有这两个特征,对其进行优化时,会得到一个窄长的椭圆形,导致在梯度下降时,梯度的方向为垂直等高线的方向而走“之”字形路线,这样会使迭代很慢。相比之下,归一化之后,是一个圆形,梯度的方向为直接指向圆心,迭代就会很快。可见,归一化可以大大减少寻找最优解的时间。
(2) 提升模型精度,归一化的另一个好处是提高精度,这在涉及一些距离计算的算法时效果显著,比如算法要计算欧几里得距离,上面x2的取值范围比较小,涉及距离计算时其对结果的影响远比x1带来的小,所以这就会造成精度的损失。
因此归一化很有必要,它可以让各个特征对结果做出的贡献相同。在求解线性回归的模型时,还有一个问题要注意,那就是特征组合问题,比如房子的长度和宽度作为两个特征参与模型的构造,不如将二者相乘得到的面积作为一个特征来进行求解,这样就在特征选择上做了减少维度的工作。
上述最优化问题实际上是有解析解的,可以用最小二乘法求解解析解,该问题称为多元线性回归(multivariate linear regression)。令: 


w~=(w(1),w(2),…,w(n),b)T=(wT,b)T

x~=(x(1),x(2),…,x(n),1)T=(xT,1)T

y=(y1,y2,…,yN)T
则: 


∑Ni=1(w·xi+b-yi)2=(y-(x~1,x~2,…,x~N)Tw~)T(y-(x~1,x~2,…,x~N)Tw~)
令: 


x=(x~1,x~2,…,x~N)T=x~T1

x~T2

︙

x~TN=
x(1)1x(2)1…x(n)11

x(1)2x(2)2…x(n)21

︙︙︙︙

x(1)Nx(2)N…x(n)N1
则: 


w~*= argminw~(y-xw~)T(y-xw~)
令Ew~=(y-xw~)T(y-xw~),求它的极小值。对w~求导令导数为零,得到解析解: 


Ew~w~=2xT(xw~-y)=0xTxw~-xTy

 当xTx为满秩矩阵或者正定矩阵时,可得: 


w~*= (xTx)-1xTy
其中(xTx)-1为xTx的逆矩阵。于是得到的多元线性回归模型为: 


f(x~i)=x~Tiw~*

 xTx不是满秩矩阵。比如N<n(样本数量小于特征种类的数量),根据x的秩小于或等于(N,n)中的最小值,即小于或等于N(矩阵的秩一定小于或等于矩阵的行数和列数); 而矩阵xTx是n×n大小的,它的秩一定小于或等于N,因此不是满秩矩阵。此时存在多个解析解。常见的做法是引入正则化项,如L1正则化或者L2正则化。以L2正则化为例: 


w~*= argminw~(y-xw~)T(y-xw~)+λ‖w~‖22
其中,λ>0调整正则化项与均方误差的比例; ‖…‖2为L2范数。
根据上述原理,得到多元线性回归算法。
 输入: 数据集T={(x1,y1),(x2,y2),…,(xN,yN)},xi∈XRn,yi∈YR,i=1,2,…,N,正则化项系数λ>0。
 输出: 


f(x)=w·x+b

 算法步骤。

① 令: 


w~=(w(1),w(2),…,w(n),b)T=(wT,b)T

x~=(x(1),x(2),…,x(n),1)T=(xT,1)T

y=(y1,y2,…,yN)T
计算: 


x=(x~1,x~2,…,x~N)T=x~T1

x~T2

︙

x~TN=
x(1)1x(2)1…x(n)11

x(1)2x(2)2…x(n)21

︙︙︙1

x(1)Nx(2)N…x(n)N1
② 求解: 


w~*= argminw~(y-xw~)T(y-xw~)+λ‖w~‖22

③ 最终的模型: 


f(x~i)=x~Tiw~*
5.2.2Python实现
前面对普通线性回归的求解过程进行了介绍,下面直接通过Python来实现线性回归的求解。
【例51】Python实现线性回归模型。

#导入必要的编程库

import numpy as np

from pylab import *

"""训练得到w和b的向量"""

def train_wb(X, y):

"""

:param X:N*D的数据

:param y:X对应的y值

:return: 返回(w,b)的向量

"""

if np.linalg.det(X.T * X) != 0:

wb = ((X.T.dot(X).I).dot(X.T)).dot(y)

return wb

def test(x, wb):

return x.T.dot(wb)

"""获得数据的函数"""

def getdata():

x = []; y = []

file = open("ex0.txt", 'r')

for line in file.readlines():

temp = line.strip().split("\t")

x.append([float(temp[0]),float(temp[1])])

y.append(float(temp[2]))

return (np.mat(x), np.mat(y).T)

"""画图函数,分别把训练用的数据的散点图还有回归直线画出来了"""

def draw(x, y, wb):

#画回归直线y = wx+b

a = np.linspace(0, np.max(x))#横坐标的取值范围

b = wb[0] + a * wb[1]

plot(x, y, '.')

plot(a, b)

show()

X, y = getdata()

wb = train_wb(X, y)

draw(X[:, 1], y, wb.tolist())

运行程序,效果如图51所示。


图51线性回归


【例52】以表51为例,针对运输里程、运输次数与运输总时间的关系,建立多元线性回归模型。


表51运输里程、运输次数与运输总时间的关系



运 输 里 程运 输 次 数运输总时间

10049.3
5034.8
10048.9
10026.5
5024.2
8026.2
7537.4
6546.0
9037.6
9026.1

Python的实现代码为: 

import numpy as np

from sklearn import datasets,linear_model

#定义训练数据

x = np.array([[100,4,9.3],[50,3,4.8],[100,4,8.9],

[100,2,6.5],[50,2,4.2],[80,2,6.2],

[75,3,7.4],[65,4,6],[90,3,7.6],[90,2,6.1]])

print(x)

X = x[:,:-1]

Y = x[:,-1]

print(X,Y)

#训练数据

regr = linear_model.LinearRegression()

regr.fit(X,Y)

print('coefficients(b1,b2...):',regr.coef_)

print('intercept(b0):',regr.intercept_)

#预测

x_test = np.array([[102,6],[100,4]])

y_test = regr.predict(x_test)

print(y_test)

运行程序,输出如下: 

[[100.4.9.3]

[ 50.3.4.8]

[100.4.8.9]

[100.2.6.5]

[ 50.2.4.2]

[ 80.2.6.2]

[ 75.3.7.4]

[ 65.4.6. ]

[ 90.3.7.6]

[ 90.2.6.1]]

[[100. 4.]

[ 50. 3.]

[100. 4.]

[100. 2.]

[ 50. 2.]

[ 80. 2.]

[ 75. 3.]

[ 65. 4.]

[ 90. 3.]

[ 90. 2.]] [9.3 4.8 8.9 6.5 4.2 6.2 7.4 6.7.6 6.1]

coefficients(b1,b2...): [0.06113460.92342537]

intercept(b0): -0.868701466781709

[10.907579818.93845988]


如果特征向量中存在分类型变量,例如车型,则需要进行特殊处理,数据如表52所示。


表52车型、运输里程、运输次数与运输总时间的关系



运 输 里 程输 出 次 数车型隐 式 转 换运输总时间

100410109.3
50301004.8
100410108.9
100220016.5
50220014.2续表



运 输 里 程输 出 次 数车型隐 式 转 换运输总时间

80210106.2
75310107.4
65401006.0
90301007.6
100410109.3
50301004.8
100410108.9
100220016.5

Python的实现代码为: 

import numpy as np

from sklearn.feature_extraction import DictVectorizer

from sklearn import linear_model

#定义数据集

x = np.array([[100,4,1,9.3],[50,3,0,4.8],[100,4,1,8.9],

[100,2,2,6.5],[50,2,2,4.2],[80,2,1,6.2],

[75,3,1,7.4],[65,4,0,6],[90,3,0,7.6],

[100,4,1,9.3],[50,3,0,4.8],[100,4,1,8.9],[100,2,2,6.5]])

x_trans = []

for i in range(len(x)):

x_trans.append({'x1':str(x[i][2])})

vec = DictVectorizer()

dummyX = vec.fit_transform(x_trans).toarray()

x = np.concatenate((x[:,:-2],dummyX[:,:],x[:,-1].reshape(len(x),1)),axis=1)

x = x.astype(float)

X = x[:,:-1]

Y = x[:,-1]

print(x,X,Y)

#训练数据

regr = linear_model.LinearRegression()

regr.fit(X,Y)

print('coefficients(b1,b2...):',regr.coef_)

print('intercept(b0):',regr.intercept_)

运行程序,输出如下: 

[[100.4.0.1.0.9.3]

[ 50.3.1.0.0.4.8]

[100.4.0.1.0.8.9]

[100.2.0.0.1.6.5]

[ 50.2.0.0.1.4.2]

[ 80.2.0.1.0.6.2]

[ 75.3.0.1.0.7.4]

[ 65.4.1.0.0.6. ]

[ 90.3.1.0.0.7.6]

[100.4.0.1.0.9.3]

[ 50.3.1.0.0.4.8]

[100.4.0.1.0.8.9]

[100.2.0.0.1.6.5]] [[100. 4. 0. 1. 0.]

[ 50. 3. 1. 0. 0.]

[100. 4. 0. 1. 0.]

[100. 2. 0. 0. 1.]

[ 50. 2. 0. 0. 1.]

[ 80. 2. 0. 1. 0.]

[ 75. 3. 0. 1. 0.]

[ 65. 4. 1. 0. 0.]

[ 90. 3. 1. 0. 0.]

[100. 4. 0. 1. 0.]

[ 50. 3. 1. 0. 0.]

[100. 4. 0. 1. 0.]

[100. 2. 0. 0. 1.]] 

[9.3 4.8 8.9 6.5 4.2 6.2 7.4 6.7.6 9.3 4.8 8.9 6.5]

coefficients(b1,b2...): [ 0.054525070.70930079 -0.180196420.60821607 -0.42801964]

intercept(b0): 0.19899589563177766

5.3广义线性模型
考虑单调可导函数h(·),令h(y)=wTx+b,这样得到的模型称为广义线性模型(generalized linear model)。
广义线性模型的一个典型的例子就是对数线性回归。当h(·)=ln(·)时的广义线性模型就是对数线性回归,即


lny=wTx+b

它是通过exp(wTx+b)拟合y的。它虽然称为广义线性回归,但实质上是非线性的。
【例53】本例中使用一个2次函数加上随机的扰动来生成500个点,然后尝试用1、2、100次方的多项式对该数据进行拟合。
拟合的目的是使得根据训练数据能够拟合出一个多项式函数,这个函数能够很好地拟合现有数据,并且能对未知的数据进行预测。

import matplotlib.pyplot as plt

import numpy as np

import scipy as sp

from scipy.stats import norm

from sklearn.pipeline import Pipeline

from sklearn.linear_model import LinearRegression

from sklearn.preprocessing import PolynomialFeatures

from sklearn import linear_model

''''' 数据生成 '''

x = np.arange(0, 1, 0.002)

y = norm.rvs(0, size=500, scale=0.1)

y = y + x**2

''''' 均方误差根 '''

def rmse(y_test, y):

return sp.sqrt(sp.mean((y_test - y) ** 2))

''''' 与均值相比的优秀程度,介于0~1。0表示不如均值,1表示完美预测。这个版本的实现是参考scikit-learn官网文档'''

def R2(y_test, y_true):

return 1 - ((y_test - y_true)**2).sum() / ((y_true - y_true.mean())**2).sum()

def R22(y_test, y_true):

y_mean = np.array(y_true)

y_mean[:] = y_mean.mean()

return 1 - rmse(y_test, y_true) / rmse(y_mean, y_true)

plt.scatter(x, y, s=5)

degree = [1,2,100]

y_test = []

y_test = np.array(y_test)

for d in degree:

clf = Pipeline([('poly', PolynomialFeatures(degree=d)),

('linear', LinearRegression(fit_intercept=False))])

clf.fit(x[:, np.newaxis], y)

y_test = clf.predict(x[:, np.newaxis])

print(clf.named_steps['linear'].coef_)

print('rmse=%.2f, R2=%.2f, R22=%.2f, clf.score=%.2f' %

(rmse(y_test, y),

R2(y_test, y),

R22(y_test, y),

clf.score(x[:, np.newaxis], y))) 

plt.plot(x, y_test, linewidth=2)

plt.grid()

plt.legend(['1','2','100'], loc='upper left')

plt.show()

运行程序,输出如下,效果如图52所示。

[-0.176469091.01838941]

rmse=0.12, R2=0.86, R22=0.63, clf.score=0.86

[-0.028726930.1283764 0.8917966 ]

rmse=0.10, R2=0.91, R22=0.69, clf.score=0.91

[-1.11632925e-014.29980995e+01 -4.76106227e+032.49103305e+05

-7.46452218e+061.41084996e+08 -1.78775541e+091.58133296e+10

-1.00114941e+114.58895566e+11 -1.51772245e+123.53223913e+12

…

4.46276417e+11 -5.96923498e+11 -1.51906488e+12 -1.90374020e+12

-1.63247529e+12 -7.23567011e+115.91292805e+111.81499049e+12

2.17224445e+121.00237307e+12 -1.37732005e+12 -3.02020489e+12

1.81820159e+12]

rmse=0.09, R2=0.91, R22=0.70, clf.score=0.91



图52广义线性回归效果1


这里要注意几点: 
(1) 误差分析。做回归分析,常用的误差主要有均方误差根(RMSE)和R平方(R2)。RMSE是预测值与真实值的误差平方根的均值。这种度量方法很流行(Netflix机器学习比赛的评价方法),是一种定量的权衡方法。
R2方法是将预测值跟使用均值的情况下相比,看能好多少。R2通常在(0,1)区间。0表示还不如什么都不预测,直接取均值的情况,而1表示所有预测跟真实结果完美匹配的情况。我们看到多项式次数为1的时候,虽然拟合得不太好,但R2也能达到0.82。二次多项式提高到了0.88。而次数提高到100次时,R2也只提高到了0.89。
(2) 过拟合。使用100次方多项式做拟合,效果确实是好了一些,但该模型的据测能力很差。而且注意看多项式系数,出现了大量的大数值,甚至达到10的12次方。这里修改代码,将500个样本中的最后2个从训练集中移除,但在测试中仍然会测试所有500个样本。


clf.fit(x[:498, np.newaxis], y[:498])


这样修改后的多项式拟合结果如下,效果如图53所示。

[-0.164906190.97634264]

rmse=0.13, R2=0.83, R22=0.59, clf.score=0.83

[ 0.0046522-0.049212131.03174524]

rmse=0.10, R2=0.90, R22=0.68, clf.score=0.90

……

rmse=0.35, R2=-0.31, R22=-0.15, clf.score=-0.31



彩色图片
图53




图53广义线性回归效果2


仅仅只是缺少了最后2个训练样本,绿线(100次方多项式拟合结果)的预测发生了剧烈的偏差,R2也急剧下降到0.57。而反观一次多项式和二次多项式的拟合结果,R2反而略微上升了。
这说明高次多项式过度拟合了训练数据,包括其中大量的噪声,导致其完全丧失了对数据趋势的预测能力。前面也看到,100次多项式拟合出的系数数值无比巨大。人们自然想到通过在拟合过程中限制这些系数数值的大小来避免生成这种畸形的拟合函数。
其基本原理是将拟合多项式的所有系数绝对值之和(L1正则化)或者平方和(L2正则化)加入惩罚模型中,并指定一个惩罚力度因子w,来避免产生这种畸形系数。
这样的思想应用在了岭(Ridge)回归(使用L2正则化)、Lasso法(使用L1正则化)、弹性网(Elastic net,使用L1+L2正则化)等方法中,都能有效避免过拟合。
下面以岭回归为例看看100次多项式的拟合是否有效。将代码修改如下: 


clf = Pipeline([('poly', PolynomialFeatures(degree=d)),

('linear', linear_model.Ridge ())])

clf.fit(x[:400, np.newaxis], y[:400])


修改程序后,输出如下,效果如图54所示。

[0. 0.76953671]

rmse=0.15, R2=0.78, R22=0.53, clf.score=0.78

[0. 0.29611609 0.62106811]

rmse=0.11, R2=0.88, R22=0.65, clf.score=0.88

…

rmse=0.11, R2=0.88, R22=0.65, clf.score=0.88



图54岭回归效果


由图54可以看到,100次多项式的系数参数变得很小,大部分都接近于0。
另外值得注意的是,使用岭回归之类的惩罚模型后,一次多项式或二次多项式回归的R2值可能会稍微低于基本线性回归。
这样的模型即使使用100次多项式,在训练400个样本、预测500个样本的情况下不仅有更小的R2误差,还具备优秀的预测能力。
5.4逻辑回归
线性模型也可以用于分类。考虑二分类问题,给定数据集T={(x1,y1),(x2,y2),…,(xN,yN)},xi∈XRn,yi∈YR,i=1,2,…,N,其中xi= (x(1)i,x(2)i,…,x(n)i)T。我们需要知道P(y/x),这里用条件概率的原因是: 预测的时候都是已知x,然后需要判断此时对应的y值。
考虑到w·x+b取值是连续的,因此它不能拟合离散变量。可以考虑用它来拟合条件概率P(y=1/x),因为概率的取值也是连续的,但是对于w≠0(如果等于零向量,则没有求解的价值),w·x+b取值是从-∞~+∞,不符合概率取值为0~1,因此考虑采用广义线性模型,最理想的是单位阶跃函数: 


P(y=1/x)=0,z<0

0.5,z=0

1,z>0,z=w·x+b

但是阶跃函数不满足单调可导的性质。退而求其次,我们寻找一个可导的、与阶跃函数相似的函数。对数概率函数(logistic function)就是这样的一个替代函数: 


P(y=1/x)=11+e-z,z=w·x+b

由于P(y=0/x)=1-P(y=1/x),所以有lnP(y=1/x)P(y=0/x)=z=w·x+b。lnP(y=1/x)P(y=0/x)表示样本为正例的可能性与负例的可能性之比,称为概率(odds),反映了样本作为正例的相对可能性。概率的对数称为对数概率(log odds,也称为logit)。
下面给出逻辑回归模型参数估计: 给定训练数据集T={(x1,y1),(x2,y2),…,(xN,yN)},xi∈XRn,yi∈{0,1}。模型估计的原理是: 用极大似然估计法估计模型参数。
为了便于讨论,将参数b吸收进w中,令: 


w~=(w(1),w(2),…,w(n),b)T∈Rn+1,x~= (x(1),x(2),…,x(n),1)T∈Rn+1

令P(y=1/x)=π(x~)=exp(w~·x~)1+exp(w~·x~),P(y=0/x)=1-π(x~),则似然函数为: 


∏Ni=1[π(x~)]yi[1-π(x~)]1-yi
对数似然函数为: 


L(w~)=∑Ni=1[yilogπ(x~)+(1-yi)log(1-π(x~))]

=∑Ni=1yilogπ(x~)1-π(x~)+log(1-π(x~))

又由于π(x~)=exp(w~·x~)1+exp(w~·x~),因此: 


L(w~)=∑Ni=1[yi(w~·x~)-log(1+exp((w~·x~))]

对L(w~)求极大值,得到w~的估计值。设估计值为w~*,则逻辑回归模型为: 


P(Y=1/X= x~)=exp(w~*·x~)1+exp(w~*·x~)

P(Y=0/X= x~)=11+exp(w~*·x~)

提示: 通常用梯度下降法或者拟牛顿法来求解该最大值问题。
以上讨论的都是二分类的逻辑回归模型,可以推广到多分类逻辑回归模型。设离散型随机变量Y的取值集合为{1,2,…,K},则多分类逻辑回归模型为: 


P(Y=k/x~)=exp(w~k·x~)1+∑K-1k=1exp(w~k·x~),k=1,2,…,K-1

P(Y=K/x~)=11+∑K-1k=1exp(w~k·x~),x~∈Rn+1,w~k∈Rn+1

其参数估计方法与二分类逻辑回归模型类似。
【例54】对给定的数据进行逻辑回归分析。

from numpy import *

filename=('testSet.txt')#文件目录

def loadDataSet():#读取数据(这里只有两个特征)

dataMat = []

labelMat = []

fr = open(filename)

for line in fr.readlines():

lineArr = line.strip().split()

dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])#代码中的1.0,表示方程的常量比如两个特征X1,X2,共需要三个参数,W1+W2*X1+W3*X2

labelMat.append(int(lineArr[2]))

return dataMat,labelMat

def sigmoid(inX):#sigmoid函数

return 1.0/(1+exp(-inX))

def gradAscent(dataMat, labelMat):#梯度上升求最优参数

dataMatrix=mat(dataMat)#将读取的数据转换为矩阵

classLabels=mat(labelMat).transpose()#将读取的数据转换为矩阵

m,n = shape(dataMatrix)

alpha = 0.001#设置梯度的阈值,该值越大梯度上升幅度越大

maxCycles = 500#设置迭代的次数,一般看实际数据进行设定,

#有些可能200次就够了

weights = ones((n,1))#设置初始的参数,并都赋默认值为1。注意这里

#权重以矩阵形式表示三个参数

for k in range(maxCycles):

h = sigmoid(dataMatrix*weights)

error = (classLabels - h)#求导后差值


weights = weights + alpha * dataMatrix.transpose()* error#迭代更新权重

return weights

"""随机梯度上升,当数据量比较大时,每次迭代都选择全量数据进行计算,计算量会非常大。所以采用每次迭代中一次只选择其中的一行数据进行更新权重。"""

def stocGradAscent0(dataMat, labelMat):

dataMatrix=mat(dataMat)

classLabels=labelMat

m,n=shape(dataMatrix)

alpha=0.01

maxCycles = 500

weights=ones((n,1))

for k in range(maxCycles):

for i in range(m):#遍历计算每一行

h = sigmoid(sum(dataMatrix[i] * weights))

error = classLabels[i] - h

weights = weights + alpha * error * dataMatrix[i].transpose()

return weights

"""改进版随机梯度上升,在每次迭代中随机选择样本来更新权重,并且随迭代次数增加,权重变化越小。"""

def stocGradAscent1(dataMat, labelMat):

dataMatrix=mat(dataMat)

classLabels=labelMat

m,n=shape(dataMatrix)

weights=ones((n,1))

maxCycles=500

for j in range(maxCycles):#迭代

dataIndex=[i for i in range(m)]

for i in range(m):#随机遍历每一行

alpha=4/(1+j+i)+0.0001#随迭代次数增加,权重变化越小

randIndex=int(random.uniform(0,len(dataIndex)))#随机采样

h=sigmoid(sum(dataMatrix[randIndex]*weights))

error=classLabels[randIndex]-h

weights=weights+alpha*error*dataMatrix[randIndex].transpose()

del(dataIndex[randIndex])#去除已经抽取的样本

return weights

def plotBestFit(weights):#画出最终分类的图

import matplotlib.pyplot as plt

dataMat,labelMat=loadDataSet()

dataArr = array(dataMat)

n = shape(dataArr)[0]

xcord1 = []; ycord1 = []

xcord2 = []; ycord2 = []

for i in range(n):

if int(labelMat[i])== 1:

xcord1.append(dataArr[i,1])

ycord1.append(dataArr[i,2])

else:

xcord2.append(dataArr[i,1])

ycord2.append(dataArr[i,2])

fig = plt.figure()

ax = fig.add_subplot(111)

ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')

ax.scatter(xcord2, ycord2, s=30, c='green')

x = arange(-3.0, 3.0, 0.1)

y = (-weights[0]-weights[1]*x)/weights[2]

ax.plot(x, y)

plt.xlabel('X1')

plt.ylabel('X2')

plt.show()

def main():

dataMat, labelMat = loadDataSet()

weights=gradAscent(dataMat, labelMat).getA()

plotBestFit(weights)

if __name__=='__main__':

main()

运行程序,效果如图55所示。



图55逻辑回归分析


5.5岭回归
岭回归是在最小二乘法的基础之上的,所以要了解岭回归,首先要了解最小二乘的回归原理。设有多重线性回归模型y=xβ+ε,参数β的最小二乘估计为: 


β^=(xTx)-1xTy

当自变量间存在多重共线性,即|xTx|≈0时,设想给|xTx|加上一个正常数矩阵(k>0),那么|xTx|+kI接近奇异的程度就会比接近奇异的程度小得多。考虑到变量的量纲问题,要先对数据标准化,标准化捕捉设计矩阵仍用x表示,定义称为岭回归估计,其中,k称为岭参数。由于假设x已经标准化,所以就是自变量样本相关阵。y可以标准化也可以未标准化,如果y也经过标准化,那么计算的实际是标准化岭回归估计。k作为β的估计应比最小二乘估计稳定,k=0时的岭回归估计就是普通的最小二乘估计。因为岭参数k不是唯一确定的,所以得到的岭回归实际是回归参数的一个估计族。
岭回归的参数估计为: 


β^(k)=(xTx+kI)-1xTy

【例55】随机产生100组数据集,每组数据集包含25个点,每个点满足: y=sin(2πx)+ε,这里x∈{0.041×i,i=1,2,…,24},e是添加的高斯噪声(0,0.32)。在每组数据集上用具有不同λ的7阶多项式进行岭回归拟合。

import numpy as np 

import matplotlib.pyplot as plt

from tkinter import _flatten 

x_arange = 0.041 * np.arange(0, 25, 1)#每组数据的25个点

y_True = np.sin(2 * np.pi * x_arange)#每个数据点对应的值(没有添加噪声)

y_Noise = np.zeros(y_True.shape)#添加噪声的值

x_Prec = np.linspace(0, 24*0.041, 100)#画图范围

mu = 0#噪声的mu值

sigma = 0.3#噪声的sigma值

Num = 100#100组数据集

n = 8#7阶多项式

lamda = [np.exp(1), np.exp(0), np.exp(-5), np.exp(-10)]#不同的lambda值

phi = np.mat(np.zeros((x_arange.size, n)))#phi矩阵

x = np.mat(x_arange).T#输入数据矩阵

#phi矩阵运算

for i_n in range(n):

for y_n in range(x_arange.size):

phi[y_n, i_n] = x[y_n, 0] ** i_n 

plt.figure(figsize=(15, 10))

index = 221

for i_lamda in lamda:

plt.subplot(index)

index += 1

plt.title("lambda = %f" % i_lamda)

plt.plot(x_Prec, np.sin(2 * np.pi * x_Prec), color='g')

for k in range(Num):

for i in range(x_arange.size):

y_Noise[i] = y_True[i] + np.random.normal(mu, sigma)

y = np.mat(y_Noise).T

#求解w参数

W = (phi.T * phi + i_lamda*np.eye(n)).I * phi.T * y

ploy = list(_flatten(W.T.tolist()))

ploy.reverse()

p = np.poly1d(ploy)

if k%5==0:#只画20条曲线

plt.plot(x_Prec, p(x_Prec), color='r')

plt.show()

运行程序,效果如图56所示。


图56不同的lambda值的岭回归效果


5.6Lasso回归
Lasso回归与岭回归非常相似,它们的差别在于使用了不同的正则化项。最终都实现了约束参数从而防止过拟合的效果。Lasso之所以重要,还有另一个原因,即Lasso能够将一些作用比较小的特征的参数训练为0,从而获得稀疏解。也就是说,用这种方法,在训练模型的过程中实现了降维(特征筛选)的目的。
Lasso回归的代价函数为: 


J(θ)=12m∑mi=1(y(i)-(wx(i)+b))2+λ‖w‖1=12MSE(θ)+λ∑mi=1|θi|(51)
式中的w是长度为n的向量,不包括截距项的系数θ0,θ是长度为n+1的向量,包括截距项的系数θ0,m为样本数,n为特征数。‖w‖1表示参数w的L1范数,也是一种表示距离的函数。加入w表示三维空间中的一个点(x,y,z),那么‖w‖1=|x|+|y|+|z|,即各个方向上的绝对值(长度)之和。
式(51)的梯度为: 


θMSE(θ)+λsigm(θ1)

sigm(θ2)

︙

sigm(θn)

其中,sigm(θi)由θi的符号决定: θi>0,sigm(θi)=1; θi=0,sigm(θi)=0; θi<0,sigm(θi)=-1。
【例56】利用Lasso回归对自带的数据集进行回归分析。

import matplotlib.pyplot as plt

import numpy as np

from sklearn import datasets, linear_model, discriminant_analysis, cross_validation

def load_data():

diabetes = datasets.load_diabetes()

return cross_validation.train_test_split(diabetes.data, diabetes.target, test_size=0.25, random_state=0)

def test_lasso(*data):

X_train, X_test, y_train, y_test = data

lassoRegression = linear_model.Lasso()

lassoRegression.fit(X_train, y_train)

print("权重向量:%s,b的值为:%.2f"%(lassoRegression.coef_, lassoRegression.intercept_))

print("损失函数的值:%.2f" % np.mean((lassoRegression.predict(X_test) - y_test) ** 2))

print("预测性能得分: %.2f" % lassoRegression.score(X_test, y_test))

#测试不同的α值对预测性能的影响

def test_lasso_alpha(*data):

X_train, X_test, y_train, y_test = data

alphas = [0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000]

scores = []

for i, alpha in enumerate(alphas):

lassoRegression = linear_model.Lasso(alpha=alpha)

lassoRegression.fit(X_train, y_train)

scores.append(lassoRegression.score(X_test, y_test))

return alphas, scores

def show_plot(alphas, scores):

figure = plt.figure()

ax = figure.add_subplot(1, 1, 1)

ax.plot(alphas, scores)

ax.set_xlabel(r"$\alpha$")

ax.set_ylabel(r"score")

ax.set_xscale("log")

ax.set_title("lasso")

plt.show()

if __name__=='__main__':

X_train, X_test, y_train, y_test = load_data()

#使用默认的alpha

#test_lasso(X_train, X_test, y_train, y_test)

#使用自己设置的alpha

alphas, scores = test_lasso_alpha(X_train, X_test, y_train, y_test)

show_plot(alphas, scores)

运行程序,效果如图57所示。


图57lasso回归


5.7弹性网络
弹性网络(Elastic Net)结合了岭回归和Lasso回归,由两者加权平均所得。据介绍,这种方法在特征数大于训练集样本数或有些特征之间高度相关时比Lasso更加稳定。其代价函数为: 


J(θ)=12MSE(θ)+rλ∑mi=1|θi|+1-r2λ∑ni=1θ2i
其中,r表示L1所占的比例。
【例57】用Python实现弹性网络算法(多变量)。使用鸢尾花数据集,后3个特征作为特征,用来预测第一个特征。

#导入必要的编程库,创建计算图,加载数据集

import matplotlib.pyplot as plt 

import tensorflow as tf 

import numpy as np 

from sklearn import datasets 

from tensorflow.python.framework import ops 

ops.get_default_graph() 

sess = tf.Session() 

iris = datasets.load_iris() 

x_vals = np.array([[x[1], x[2], x[3]] for x in iris.data]) 

y_vals = np.array([y[0] for y in iris.data]) 

#声明学习率,批量大小,占位符和模型变量,模型输出

learning_rate = 0.001

batch_size = 50

x_data = tf.placeholder(shape=[None, 3], dtype=tf.float32)#占位符大小为3 

y_target = tf.placeholder(shape=[None, 1], dtype=tf.float32) 

A = tf.Variable(tf.random_normal(shape=[3,1])) 

b = tf.Variable(tf.random_normal(shape=[1,1])) 

model_output = tf.add(tf.matmul(x_data, A), b) 

#对于弹性网络回归算法,损失函数包括L1正则和L2正则

elastic_param1 = tf.constant(1.) 

elastic_param2 = tf.constant(1.) 


l1_a_loss = tf.reduce_mean(abs(A)) 

l2_a_loss = tf.reduce_mean(tf.square(A)) 

e1_term = tf.multiply(elastic_param1, l1_a_loss) 

e2_term = tf.multiply(elastic_param2, l2_a_loss) 

loss = tf.expand_dims(tf.add(tf.add(tf.reduce_mean(tf.square(y_target - model_output)), e1_term), e2_term), 0) 

#初始化变量,声明优化器,然后遍历迭代运行,训练拟合得到参数

init = tf.global_variables_initializer() 

sess.run(init) 

my_opt = tf.train.GradientDescentOptimizer(learning_rate) 

train_step = my_opt.minimize(loss) 

loss_vec = [] 

for i in range(1000): 

rand_index = np.random.choice(len(x_vals), size=batch_size) 

rand_x = x_vals[rand_index] 

rand_y = np.transpose([y_vals[rand_index]]) 


sess.run(train_step, feed_dict={x_data:rand_x, y_target:rand_y}) 

temp_loss = sess.run(loss, feed_dict={x_data:rand_x, y_target:rand_y}) 

loss_vec.append(temp_loss) 

if (i+1)%250 == 0: 

print('Step#' + str(i+1) +'A = ' + str(sess.run(A)) + 'b=' + str(sess.run(b))) 

print('Loss= ' +str(temp_loss))

#现在能观察到,随着训练迭代后损失函数已收敛

plt.plot(loss_vec, 'k--') 

plt.title('Loss per Generation') 

plt.xlabel('Generation') 

plt.ylabel('Loss') 

plt.show()

运行程序,输出如下,效果如图58所示。

Step#250A = [[1.3357686 ]

[0.7352283 ]

[0.12850094]]b=[[-1.4265894]]

Loss= [1.9901819]

Step#500A = [[1.3833795 ]

[0.6989296 ]

[0.03573647]]b=[[-1.2994616]]

Loss= [1.9257231]

Step#750A = [[ 1.3752384e+00]

[ 6.9431901e-01]

[-3.4378259e-05]]b=[[-1.1875116]]

Loss= [1.7963648]

Step#1000A = [[1.3545375e+00]

[6.8431729e-01]

[3.3413176e-05]]b=[[-1.0804418]]

Loss= [1.7084534]



图58弹性网络算法


5.8线性判别分析
线性判别分析(Linear Discriminant Analysis,LDA)的思想是: 
 训练时,设法将训练样本投影到一条直线上,使得同类样本的投影点尽可能地接近,异类样本的投影点尽可能地远离。要学习的就是这样的一条直线。
 预测时,将待预测样本投影到学到的直线上,根据它的投影点的位置来判定它的类别。
5.8.1线性判别二分类情况
回顾之前的逻辑回归方法,给定m个n维特征的训练样例,每个x(i),i=1,2,…,m对应一个类标签y(i)。我们就是要学习出参数θ,使得y(i)=g(θ)Tx(i)。
现在只考虑二分类情况,也就是y=1或y=0。为了方便表示,先换符号重新定义问题,给定特征为d维的N个样例,x(i)=(x(i)1,x(i)2,…,x(i)1),其中有N1个样例属于类别w1,另外N2个样例属于类别w2。
现在我们觉得原始特征数太多,想将d维特征降到只有一维,又要保证类别能够“清晰”地反映在低维数据上,也就是说,这一维就能决定每个样例的类别。我们将这个最佳的向量称为w(d维),那么样例x(d维)到w上的投影可以用下式来计算: 


y=wTx

这里得到的y值不是0/1值,而是x投影到直线上的点到原点的距离。若x是二维的,则要找一条直线(方向为w)来做投影,然后寻找最能使样本点分离的直线,如图59所示。



图59直线投影效果


从直观上看,右图比较好,可以很好地将不同类别的样本点分离。接着从定量的角度来找到这个最佳的w。
首先寻找每类样例的均值(中心点),这里i只有两个: 


μi=1Ni∑x∈wix

由于x到w投影后的样本点均值为: 


μ~i=1Ni∑y∈wiy=1Ni∑x∈wiwTx=wTμi



图510样本点分布图

由此可知,投影后的均值也就是样本中心点的投影。那什么是最佳直线(w)呢?首先发现,能够使投影后的两类样本中心点尽量分离的直线是好的直线,定量表示就是: 



J(w)=|μ~1-μ~2|=|wT(μ1-μ2)|

J(w)越大越好,但只考虑J(w)是不行的,观察图510。

样本点均匀分布在椭圆里,投影到横轴x1上时能够获得更大的中心点间距J(w),但是由于有重叠,x1不能分离样本点。投影到纵轴x2上,虽然J(w)较小,但是能够分离样本点。因此还需要考虑样本点之间的方差,方差越大,样本点越难以分离。
我们使用另外一个度量值,称作散列值(scatter),对投影后的类求散列值,如下: 


s~2i=∑y∈wi(y-μ~i)2

从上式可以看出,只是少除以样本数量的方差值,散列值的几何意义是样本点的密集程度,值越大,越分散; 反之,越集中。而我们想要的投影后的样本点是: 不同类别的样本点越分开越好,同类的越聚集越好,也就是均值差越大越好,散列值越小越好。因此,可以使用J(w)和s来度量,最终的度量公式为: 


J(w)=|μ~1-μ~2|2s~21+s~22
因此,只需寻找使J(w)最大的w即可。
先把散列值公式展开,


s~2i=∑y∈wi(y-μ~i)2=∑x∈wi(wTx-wTμi)2=∑x∈wiwT(x-μi)(x-μi)Tw(52)

定义上式的中间部分: 


si=∑x∈wi(x-μi)(x-μi)T(53)

这个公式就是样例数的协方差矩阵,称为散列矩阵(scatter matrix)。继续定义: 


sw=s1+s2

称sw为类内散度矩阵。回到式(52)中,使用si替换中间部分,得


s~2i=wTsiw

s~21+s~22=wTsww
然后,展开分子


(μ~1-μ~2)2=(wTμ1-wTμ2)2=wT(μ1-μ2)(μ1-μ2)Tw=wTsBw
其中sB=(μ1-μ2)(μ1-μ2)T称为类间散度矩阵,是两个向量的外积,为一个矩阵,秩为1。那么J(w)最终可以表示为


J(w)=wTsBwwTsww
在求导前,需要对分母进行归一化,因为不做归一的话,无论w扩大多少倍,上式都成立,也就无法确定w。所以令‖wTsww‖=1,加入拉格朗日乘子后,求导: 


c(w)=wTsBw-λ(wTsww-1)

dcdw=2sBw-2λsww=0

sBw=λsww

用到了矩阵微积分,求导时可以简单地把wTsww当作sww2看待。如果sw可逆,那么将求导后的结果两边都乘以s-1w,得


s-1wsBw=λw(54)

即w就是矩阵s-1wsB的特征向量了。这个公式称为费希尔判别准则。由sB=(μ1-μ2)(μ1-μ2)T,有sBw=(μ1-μ2)(μ1-μ2)Tw=(μ1-μ2)λw,将其代入公式(54)中得


s-1wsBw=s-1w(μ1-μ2)λw=λw

由于对w扩大缩小任何倍不影响结果,因此可以约去两边的未知常数λ和λw,得到: 


w=s-1w(μ1-μ2)

至此,我们只需要求出原始样本的均值和方差就可以求出最佳的方向w,这就是Fisher于1936年提出的线性判别分析。二维样本的投影结果如图511所示。



图511二维样本的投影结果


5.8.2线性判别多类情况
前面是针对只有两个类的情况,假设类别变成多个了,那么要怎么改变,才能保证投影后类别能够分离呢?之前讨论的是如何将d维降到一维,现在类别多了,一维可能已经不能满足要求。假设有C个类别,需要K维向量(或者叫作基向量)来做投影。

图512类间散列度和类内散

列度几何图


将这K维向量表示为W=(w1,w2,wK)。将样本点在这K维向量投影后结果表示为(y1,y2,…,yK),有以下公式成立: 


yi=wTix

y=wTx

仍然从类间散列度和类内散列度来考虑。当样本是二维时,从几何意义上考虑,如图512所示。


其中,μi和sw与5.8.1节的含义一样,sw1是类别1里的样本点相对于该类中心点μ1的散列度。sB1变成类别1中心点相对于样本点心点μ的协方差矩阵,即类1相对于μ的散列程度。sw为: 


sw=∑ci=1swi
swi的计算公式不变,仍然类似于类内部样本点的协方差矩阵: 


swi=∑x∈wi(x-μi)(x-μi)T

sB需要改变,原来度量的是两个均值点的散列情况,现在度量的是每类的均值点相对于样本中心的散列情况。类似于将μi看作样本点,μ是均值的协方差矩阵,如果某类里面的样本点较多,那么其权重稍大,权重用NiN表示,但由于J(w)对倍数不敏感,因此使用Ni。


sB=∑ci=1Ni(μi-μ)(μi-μ)T
其中,μ=1N∑xx=1N∑x∈wiNiμi是所有样本的均值。
上面讨论的都是在投影前的公式变化,但真正的J(w)的分子分母都是在投影后计算的。下面看一下样本点投影后的公式改变: 
这两个是第i类样本点在某基向量上投影后的均值计算公式: 


μ~i=1Ni∑y∈wiy

μ~=1N∑yy
下面两个是在某基向量上投影后的sw和sB: 


s~w=∑ci=1∑y∈wi(y-μ~i)(y-μ~i)T

s~B=∑ci=1Ni(μ~i-μ~)(μ~i-μ~)T
其实就是将μ换成μ~。
综合各个投影向量(w)上的s~w和s~B,更新这两个参数,得到: 


s~w=wTsww

s~B=wTsBw
w是基向量矩阵,s~w是投影后的各个类内部的散列矩阵之和,s~B是投影后各个类中心相对于全样本中心投影的散列矩阵之和。
在5.8.1节的公式J(w),分子是两类中心距,分母是每个类自己的散列度。现在投影方向是多维了,分子需要做一些改变,我们不是求两样本中心距之和,而是求每类中心相对于全样本中心的散列度之和。然而,最后的J(w)形式为: 


J(w)=|s~B||s~w|=|wTsBw|wTsww
由于得到的分子分母都是散列矩阵,要将矩阵变成实数,需要取行列式。又因为行列式的值实际上是矩阵特征值的积,一个特征值可以表示在该特征向量上的发散程度。因此使用行列式来计算。整个问题又回归求J(w)的最大值了,固定分母为1,然后求导,得出最后结果: 


sBwi=λswwi

与5.8.1节得出的结论一样: s-1wsBwi=λwi。因此还是归结到求矩阵的特征值,首先求出s-1wsB的特征值,然后取前K个特征向量组成w矩阵即可。
5.8.3线性判别分析实现
5.8.1节和5.8.2节分别介绍了线性判别分析对二分类、多分类过程进行了介绍,下面直接通过例子来演示线性判别分析的实现。
【例58】演示线性判别分析对随机数据的判别。

'''

LDA算法实现

'''

import os

import sys

import numpy as np

from numpy import *

import operator

import matplotlib

import matplotlib.pyplot as plt

def createDataSet():

group1=mat(random.random((2,4))+10)

group2=mat(random.random((2,4))+2)

return group1, group2

#计算样本均值



#参数samples为n×m维矩阵,其中n表示维数,m表示样本个数

def compute_mean(samples):

mean_mat=mean(samples, axis=1)

return mean_mat

#计算样本类内离散度



"""参数samples表示样本向量矩阵,大小为nxm,其中n表示维数,m表示样本个数

参数mean表示均值向量,大小为1xd,d表示维数,大小与样本维数相同,即d=m"""

def compute_withinclass_scatter(samples, mean):

#获采样本维数,样本个数

dimens,nums=samples.shape[:2]

#将所有样本向量减去均值向量

samples_mean=samples-mean

#初始化类内离散度矩阵

s_in=0

for i in range(nums):

x=samples_mean[:,i]

s_in+=dot(x,x.T)

return s_in

if __name__=='__main__':

group1,group2=createDataSet()

print ("group1 :\n",group1)

print ("group2 :\n",group2)

mean1=compute_mean(group1)

print ("mean1 :\n",mean1)

mean2=compute_mean(group2)

print ("mean2 :\n",mean2)

s_in1=compute_withinclass_scatter(group1, mean1)

print ("s_in1 :\n",s_in1)

s_in2=compute_withinclass_scatter(group2, mean2)

print ("s_in2 :\n",s_in2)

#求总类内离散度矩阵

s=s_in1+s_in2

print( "s :\n",s)

#求s的逆矩阵

s_t=s.I

print ("s_t :\n",s_t)

#求解权向量

w=dot(s_t, mean1-mean2)

print ("w :\n",w)

#判断(2,3)是在哪一类

test1=mat([1,1])

g=dot(w.T, test1.T-0.5*(mean1-mean2))

print("g(x) :",g)

#判断(4,5)是在哪一类

test2=mat([10,10])

g=dot(w.T, test2.T-0.5*(mean1-mean2))

print("g(x) :",g)

运行程序,输出如下: 

[2.81742415 2.64585055 2.16746183 2.80708945]]

mean1 :

[[10.54752722]

[10.58633559]]

mean2 :

[[2.70746802]

[2.6094565 ]]

s_in1 :

[[ 0.70944616 -0.24457606]

[-0.244576060.10698882]]

s_in2 :

[[0.11233544 0.03508707]

[0.03508707 0.27899314]]

s :

[[ 0.8217816-0.20948899]

[-0.209488990.38598196]]

s_t :

[[1.41226396 0.76649631]

[0.76649631 3.00680512]]

w :

[[17.18648137]

[29.99429734]]

g(x) : [[-139.82117865]]

g(x) : [[284.80582973]]

5.9习题
1. 线性模型是一类统计模型的总称,它包括了、、和等。
2. 特征归一化有两个好处,是什么?
3. 拟合时,需要注意几点?
4.  Lasso重要的原因有哪些?
5. 创建一组随机三维数据,绘制方程z=ax+by+c的空间平面。