1 一元线性回归
1.1 为什么用回归
图1.1.1 Google的票房与搜索量的关系
图1.1显示的是Google发布的电影的搜索量与票房的关系。如何用历史的信息预测票房就是(线性)回归问题。
1.2 一元线性回归模型
1 数学描述
图1.1.1中的横、纵轴分别用用{xi,yi} 表示, 1≤i≤N。假设图1.1中使用的一元线性模型的形式为:
显然只要求出ω0,ω1线性模型就可以确定了。为了求解系数ω0,ω1需要构造一个目标函数(损失函数),如下
只要最小化式(1.2-2),就可以求出系数 a,b 。这种做法非常直观,就是要使预测的结果和真值之间的差最小。
图1.2.1E(ω0,ω1)函数的几何解释
E(ω0,ω1) 是一个非负值,最小值为0,它的几何解释如图1.2.1,就是要使 yn,tn 的距离平方和最小,回归函数要穿过真实数据 yn 。
2 矩阵表示
对于N各数据点,式(1.2-1)有N个等式,并用线性代数表示为
其中
X=[XT1 ⋮ XTN ]
,
此时
其中 Y=⎡⎣⎢⎢y1⋮yN⎤⎦⎥⎥
所以式(1.2.2)又可以表示为
【因为,对于zϵRn,zTz=∑iz2i】
3 目标函数最小化
在高等数学中,使函数一阶导数为0,且二阶导数要大于0的点为函数的最小值点。式(1.2-5)所表示的是二次函数且开口向上,只要求一阶导数为0即可
【矩阵求导: ∂(XTAX)∂X=AX+ATX ; ∂(XTA)∂X=A ; ∂(AX)∂X=AT , ∂E(ω0,ω1)∂ω=∂(ωTXTXω)∂ω−∂(ωTXTY)∂ω−∂(YTXω)∂ω ,设第一项中 XTX=A ,第二项中 XTY=B ,第三项中 YTX=C ,所以 ∂E(ω0,ω1)∂ω=[XTXω+(XTX)Tω]−XTY−(YTX)T ,所以, ∂E(ω0,ω1)∂ω=2XTXω−2XTY
1 详细推导
2 详细推导 ∂(AX)∂X
AX 是 n×1 维的,第 i 个元素为
所以线性模型所对应的最优参数 ω 表示如下,也称为正则解或者闭形式解
一个简单的例子,代码见文件夹1_regression。
第一步:用synthic_data.py中的linearSamples方法生成数据
import numpy as np
import random
def linearSamples(n = 20):
a = 0.5
b = 1.0
r = [i + 2.0*random.random() for i in xrange(n)]
return [range(0, len(r)), r]
第二步:用linear_regression.py中的lR方法,完成式(1.2-7),最终的结果为 ω 的值
def lR(x, y):
x = np.matrix(x)
if x.shape[0] == 1:
x = x.transpose()
y = np.matrix(y)
if y.shape[0] == 1:
y = y.transpose()
one = np.ones((x.shape[0], 1))
x = np.hstack([one, x])
w = inv((x.transpose()).dot(x)).dot(np.transpose(x)).dot(y)
return w
第三步:将第二步中计算的 ω 和第一步中生成的数据,传递给plotLM方法,画出的数据点和回归直线如图1.2.2
def plotLM(w, x,y):
xx = [i for i in np.arange(0.0,20.0,0.5)]
yy = [w[0,0] + w[1,0] * i for i in xx]
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x,y, '.')
ax.plot(xx,yy)
s = 'y = %s + %s * x' %(str(w[0,0])[0:7], str(w[1, 0])[0:7])
ax.annotate(s, xy=(12.5, 13.3), xycoords='data',
xytext=(-180, 30), textcoords='offset points',
bbox=dict(boxstyle="round", fc="0.8"),
arrowprops=dict(arrowstyle="->", connectionstyle="angle,angleA=0,angleB=90,rad=10"))
plt.xlabel('x')
plt.ylabel('y')
plt.legend(('training sampes','regression line'))
plt.show()
图1.2.2 线性回归的例子
在图1.2.2中散点代表的是训练数据,训练数据是由程序随机生成,没有实际意义,直线是回归直线,并标出了直线方程,在运行程序时直线结果可能与图中的结果稍有不同,因为训练数据是随机生成的缘故。
2 最优化方法-梯度下降法
在第一节的第3部分介绍了,将损失函数表示成矩阵形式,然后求导方法,求出最优的 w ,这种方法对线性问题可以求出最优解,称为闭形式解或者解析解。本节介绍的梯度下降法是数值最优化方法,普适性更强,对于非线性问题依然可以求解。
梯度下降法是最常用、也是最容易理解的最优化方法。学会了梯度下降法,其它基于梯度的改进方法:共轭梯度法、牛顿法、拟牛顿法等,就比较容易理解。
1 盲人是如何下山的
第一步:左踩一脚,右踩一脚,如果发现这两脚在在高度上没有差别,此时他所面对的应该是山顶或者山脚,反之盲人面对的应该是山脊。(计算偏导数)
第二步:上踩一脚,下踩一脚,脚低的那个方向就对着山脚。(计算偏导数)
第三步:四个脚中,高度最低的那个方向就是山脚,从当前位置向下夸一小步,向着山脚进发。(确定步长,学习率)
重复第一、二、三步,直到山脚。
2 梯度下降法
梯度法就和盲人下山类似,就两个步骤:首先确定下山方,然后再确定的方向上按照一定的步长下山。
下面介绍最优化问题。
单目标、无约束、多维最优化问题的数学描述:
其中, x∈RN 。
梯度下降法算法流程如下:
1)给定初值 x(0) ,精度 ε> 0,并令 k=1 。
2) 计算梯度下降方向(搜索方向)v(k)=−∇f(x(k)),∇f(x(k))表示f(x)在x(k)处的梯度 。
【 ∇f(x(k)) 所表示的是数值梯度,求法如下:
其中, =0.000001 】
3)若 ∣∣v(k)∣∣≤ε ,则停止计算,否则从 x(k) 出发,沿着 v(k) 一维搜索,即求 λk ,使的 f(x(k)+λkv(k))=minλ>0f(x(k)+λkv) ,此处的一维搜索可以用黄金分割法或者二次差值法等。
4)令 x(k+1)=x(k)+λkv(k) , k=k+1 ,转2)。
3 基于梯度下降法的线性回归
下面用用梯度下降法,优化目标函数:式(1.2.2),并给出相应的代码解释
第一步:依然使用linearSamples生成数据,代码见前文
第二步:完成目标函数的定义,即式(1.2-2)
def obj(x, y, w):
t = x.dot(w) - y
t = np.multiply(t, t)
sum_ = 0.5 * np.sum(t)
return sum_
第三步:完成数值梯度的定义,按照梯度下降法中的第2)点介绍,完成代码编写
def gradient(fun, x, y, w, delta = 1e-6, *args):
l = len(w)
g = []
for i in range(0, l):
delta_w = deepcopy(w)
delta_w[i] = delta_w[i] + delta
g.append(-(obj(x, y, delta_w) - obj(x, y, w))/delta)
return g
第四步:gdLR方法将实现,梯度下降法中介绍的流程1),2),4),忽略了第三步,其中的学习率由手动调整。计算结束后返回最优的 ω 。
def gdLR(fun, x, y, step = 0.0007,tol = 1e-6):
#preprocess the data
x = np.matrix(x)
if x.shape[0] == 1:
x = x.transpose()
y = np.matrix(y)
if y.shape[0] == 1:
y = y.transpose()
one = np.ones((x.shape[0], 1))
x = np.hstack([one, x])
w = [0.0, 0.0]
w = np.matrix(w)
if w.shape[0] == 1:
w = w.transpose()
l = len(w)
k = 1
while(True):
step1 = step / k
#1)compute negative gradient
g = gradient(fun, x, y, w)
err = linalg.norm(g)
print err
if err < tol or k > 200:
break
#2)updata the parameters
w = [w[i,0] + step * g[i] for i in range(0, l)]
w = np.matrix(w).transpose()
k = k + 1
return w
第五步:将闭形式的 ω 和梯度下降法的 ω ,以及数据x,y传递给方法plotGdLM画出对比图,见图2.1。
def plotGdLM(cf_w,gd_w, x,y):
xx = [i for i in np.arange(0.0,20.0,0.5)]
cf_yy = [cf_w[0,0] + cf_w[1,0] * i for i in xx]
gd_yy = [gd_w[0,0] + gd_w[1,0] * i for i in xx]
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x,y, '.')
ax.plot(xx,cf_yy,color = 'g', linewidth=3)
s = 'y = %s + %s * x' %(str(cf_w[0,0])[0:7], str(cf_w[1, 0])[0:7])
ax.annotate(s, xy=(12.5, 13.3), xycoords='data',
xytext=(-180, 30), textcoords='offset points',
bbox=dict(boxstyle="round", fc='g', ec='g'),
arrowprops=dict(arrowstyle="->",fc='g', ec='g',
connectionstyle="angle,angleA=0,angleB=90,rad=10"))
ax.plot(xx,gd_yy, color = 'r', linewidth=3)
s = 'y = %s + %s * x' %(str(gd_w[0,0])[0:7], str(gd_w[1, 0])[0:7])
ax.annotate(s, xy=(8.5, 9.3), xycoords='data',
xytext=(-180, 30), textcoords='offset points',
bbox=dict(boxstyle="round", fc='r', ec='r'),
arrowprops=dict(arrowstyle="->", fc='r', ec='r',
connectionstyle="angle,angleA=0,angleB=90,rad=10"))
plt.xlabel('x')
plt.ylabel('y')
plt.legend(('training sampes', 'closed-form regression','gradient descent regression'),loc='upper center')
plt.show()
图2.1 解析解与梯度下降法解的对比图
3 基函数
3.1 多项式回归
如果有如图2.1的数据,依然采用式(1.2-1)的模型,则回归模型如图2.2。从图2.2中可以看出,用式(1.2-1)所表示的模型无法拟合这种带多个峰的数据。一个很直观的想法是增加式(1.2-1)中的项数,用多项式拟合这种多个峰的数据
式(3.1)写成矩阵形式为
按照1.2节中的方法,也可以得到
为了与式(1.2-7)加以区别,用
X¯
代替了式(1.2-7)中的
X
,这里的
其中, K 为多项式中自变量的最高次数。
图3.1.1
图3.1.2 按照式(3.1-3),分别用3、5、10、12阶多项式拟合数据,结果如图2.3。图(d)中的拟合曲线的末端上翘与数据不吻合了。这是过拟合导致的。过拟合问题将会在下一节中介绍。
实现图3.1.3的代码如下:
第一步:生成 x+0.3∗sin(2∗pi∗x) 样本的代码
def nlSamples(n = 100):
t = np.arange(0, 1.0, 1.0 / n)
y = [ti + 0.3 * math.sin(2 * math.pi * ti)+random.random()*0.01 for ti in t]
t = list(t)
return [t, y]
第二步:按照式(3.1-3)计算模型的参数 ω
def bFLR(x, y, rank = 2):
x = np.matrix(x)
if x.shape[0] == 1:
x = x.transpose()
y = np.matrix(y)
if y.shape[0] == 1:
y = y.transpose()
one = np.ones((x.shape[0], 1))
tmp = np.zeros((x.shape[0], rank))
for i in xrange(rank):
tmp[:,i] = np.power(x.A, i + 1).transpose()
xx = np.hstack([one, tmp])
w = inv((xx.transpose()).dot(xx)).dot(np.transpose(xx)).dot(y)
return w
第三步:用第三步中的参数 ω 画出拟合的
def plotBFLR(w, x, y, rank = 2):
xx = [i for i in np.arange(0.0,1.0,1.0/20)]
w = w.A.transpose()
yy = [w.dot(xlist(i, rank))[0,0] for i in xx]
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x,y, 'ro')
ax.plot(xx,yy)
plt.xlabel('x')
plt.ylabel('y')
plt.title(str(rank) + ' order regression')
plt.legend(('training sampes','regression line'))
plt.show()
def xlist(i, rank):
l = [np.power(i ,ii) for ii in xrange(rank+1)]
l = np.array([l]).transpose()
return l
图3.1.3 多项式拟合结果
3.2 回归模型中的基函数
式(3.1-1)更一般化的表示为
其中的 ϕj(x) 称为基函数,引入基函数是为了对数据进行非线性变换,以解决非线性问题。
多项式回归中的 x,x2,⋯,xK 都可看成是基函数。
其它形式的常用基函数有:高斯基函数,逻辑蒂斯基函数,它们的表达式分别如下
高斯基函数:
其中, μj 控制着基函数的位置, s 控制着基函数的形状。
其中 σ(a)=11+exp(−a)
4 欠拟合与过拟合
4.1 欠拟合
忽略严格的数学定义,从一般的直观理解欠拟合,概念如下。
欠拟合:模型过于简单,无法捕获数据中所存在的规律,图3.1.2所示的情况就是欠拟合,因为采用的模型为 t=ω0+ω1x 形式,这种形式的模型只能拟合x和y成线性关系的数据,对于非线性的数据,应该采用更高阶的回归模型。
欠拟合的解决办法:增加模型的复杂度,如将一次多项式模型,增加到3阶或者更高阶。对比图3.1.2和图3.1.3即可发现其变化过程。
4.2 过拟合
同样也可以给出过拟合的概念如下。
过拟合:和欠拟合相对,指模型过于复杂,模型在训练数据上的训练误差很小,而在测试数据的测试误差很大,即泛化能力很差。图3.1.3中的(d)就是过拟合现象,12阶的多项式模型对于(d)中的数据复杂度太高了,其实用3阶多项式模型就能取得不错的效果。
过拟合的解决办法:
1)不改变模型,增加数据
当过拟合时,不改变模型,增加数据可以改善过拟合问题,图3.1.3中的(d)只有20个数据点,现在将数据点增加到2000个,依然用12阶的多现实拟合,结果下图
图4.2.1 增加数据点后的12阶回归模型
2)改变模型:正则化
依然对图3.1.3中的(d)的问题,如果没有足够的数据点,则可以减少模型中的特征,即将12阶模型降低为更低阶的模型,有一种方法称之为正则化,正则化方法通过在损失函数 E(ω) 上加上罚项,对高阶项进行处罚,达到降低高阶项前面的系数 ωi , ωi 变小,说明模型中 ωi 所对应的那项对模型的影响程度就会较低,达到简化模型的目的,常用的正则化后的损失函数为
按照前文介绍的方法,依然可以得到,正则化后的模型参数如下
其中,
α
为罚参数,
I
为单位阵。
按照式(4.2-2),实现的代码如下
def rTLR(x, y, lamda = 0.5,rank = 2):
x = np.matrix(x)
if x.shape[0] == 1:
x = x.transpose()
y = np.matrix(y)
if y.shape[0] == 1:
y = y.transpose()
one = np.ones((x.shape[0], 1))
tmp = np.zeros((x.shape[0], rank))
for i in xrange(rank):
tmp[:,i] = np.power(x.A, i + 1).transpose()
xx = np.hstack([one, tmp])
dim = xx.shape[1]
I = lamda * np.diag(np.ones(dim))
w = inv(I + (xx.transpose()).dot(xx)).dot(np.transpose(xx)).dot(y)
return w
图4.2.2是不同正则化参数下的回归曲线,从中可以看出 λ 越大,惩罚越强, α=0.1 时,多项式回归模型已经欠拟合了,对应于多项式中的高次项的系数接近于0了。随着 λ 变小,惩罚程度减弱,高次项的系数基本上没有变小,见(d)。
图4.2.2 不同罚参数 α 下的拟合曲线
5 多元线性回归
以上介绍的是一元回归,如果自变量的个数不止一个,就会要求使用多元回归,多元线性回归最简单的形式如下
其中 D 为自变量的个数。但是这种形式的多元回归模型有很多限制。
所以和一元回归类似,也可以引入基函数的概念,引入基函数后的表达如下
其中, x=(x1,⋯,xM)T , K 为基函数的个数
同理可以得到模型的参数
其中
下面看一个二元回归的例子。
有如图5.1所示的曲面,建立一个回归模型拟合这个曲面,由于这个曲面是个二次曲面,所以在选择基函数时可以选择到二次或者更高次,比如选择 ϕ1=x , ϕ2=x2 ,此时
按照式(5.4)代入式(5.3)可以求得模型参数 ω 。
图5.1 二维曲面
下面看看用以上思路能否拟合出图5.1的曲面。
第一步:生成图5.1所示的数据,在x,y方向分别等距离采集数据点40个,Python代码如下:
第二步:用生成的数据点,按照式(5-3)求出参数
ω
第三步:用第二步中计算的参数
ω
,并且用第一步中的方法生成测试数据,在x,y方向上分别采样20个点,用这20个点作为测试数据,输入到模型中,看模型预测的值和真实的值之间的误差。图5.2中,将真实值和预测点按照对应的一行一行的展开,分别做成两个一维向量,这样便于比对。从图中看出,预测的结果还是精确的。
图5.2 真实值和预测值得对比图
实际上线性回归模型的非线性拟合能力较差,对于图5.3的多峰值函数就无能为力,本来这里给出的例子想用图5.3,结果,解释采用20阶以上的多项式也拟合不出图5.3的样子,无奈采用图5.1,相对简单一些,不过神经网络可以对高度非线性数据进行拟合,详细的Python代码和参考文献可以见
http://blog.csdn.net/zc02051126/article/details/9337319
图5.3 Matlab中的peaks函数产生的曲面
6 应用实例-Google票房预测模型
“相关推荐”对你有帮助么?
-
非常没帮助
-
没帮助
-
一般
-
有帮助
-
非常有帮助