机器学习基础 第二章 预测算法

1 一元线性回归

1.1 为什么用回归

 

 

 

图1.1.1 Google的票房与搜索量的关系

图1.1显示的是Google发布的电影的搜索量与票房的关系。如何用历史的信息预测票房就是（线性）回归问题。

1.2 一元线性回归模型

1 数学描述

图1.1.1中的横、纵轴分别用用{xi,yi}图1.1.1中的横、纵轴分别用用{xi,yi}表示，表示，1≤i≤N。假设图1.1中使用的一元线性模型的形式为：1≤i≤N。假设图1.1中使用的一元线性模型的形式为：

 

 

t=ω0+ω1x−−−−−（1.2−1）t=ω0+ω1x−−−−−（1.2−1）

 

显然只要求出ω0,ω1线性模型就可以确定了。为了求解系数ω0,ω1需要构造一个目标函数（损失函数），如下显然只要求出ω0,ω1线性模型就可以确定了。为了求解系数ω0,ω1需要构造一个目标函数（损失函数），如下

 

 

E(ω0,ω1)=12∑i=1N(ti−yi)2=12∑i=1N(ω0+ω1xi−yi)2−−−−−（1.2−2）E(ω0,ω1)=12∑i=1N(ti−yi)2=12∑i=1N(ω0+ω1xi−yi)2−−−−−（1.2−2）

 

只要最小化式（1.2-2），就可以求出系数a,ba,b。这种做法非常直观，就是要使预测的结果和真值之间的差最小。

 

 
图1.2.1E(ω0,ω1)函数的几何解释图1.2.1E(ω0,ω1)函数的几何解释

 

E(ω0,ω1)E(ω0,ω1)是一个非负值，最小值为0，它的几何解释如图1.2.1，就是要使yn,tnyn,tn的距离平方和最小，回归函数要穿过真实数据ynyn。

2 矩阵表示

对于N各数据点，式（1.2-1）有N个等式，并用线性代数表示为

 

 

⎡⎣⎢⎢t1⋮tN⎤⎦⎥⎥=⎡⎣⎢⎢ω0+ω1x1⋮ω0+ω1xN⎤⎦⎥⎥=Xω−−−−−（1.2−3）[t1⋮tN]=[ω0+ω1x1⋮ω0+ω1xN]=Xω−−−−−（1.2−3）

 

其中X=[XT1 ⋮ XTN ]X=[X1T ⋮ XNT ]，，Xi=[1,xi]TXi=[1,xi]T，，ω=[ω0,ω1]Tω=[ω0,ω1]T

此时

 

 

⎡⎣⎢⎢t1−y1⋮tN−yN⎤⎦⎥⎥=⎡⎣⎢⎢ω0+ω1x1−y1⋮ω0+ω1xN−y1⎤⎦⎥⎥=Xω−Y−−−−−−（1.2−4）[t1−y1⋮tN−yN]=[ω0+ω1x1−y1⋮ω0+ω1xN−y1]=Xω−Y−−−−−−（1.2−4）

 

其中Y=⎡⎣⎢⎢y1⋮yN⎤⎦⎥⎥Y=[y1⋮yN]

所以式（1.2.2）又可以表示为

 

 

E(ω0,ω1)=12 (Xω−Y)T(Xω−Y)=ωTXTXω−ωTXTY−YTXω+YTY−−−−−（1.2−5）E(ω0,ω1)=12 (Xω−Y)T(Xω−Y)=ωTXTXω−ωTXTY−YTXω+YTY−−−−−（1.2−5）

 

【因为，对于zϵRn，zTz=∑iz2i】【因为，对于zϵRn，zTz=∑izi2】

3 目标函数最小化

在高等数学中，使函数一阶导数为0，且二阶导数要大于0的点为函数的最小值点。式（1.2-5）所表示的是二次函数且开口向上，只要求一阶导数为0即可

 

 

∂E(ω0,ω1)∂ω=XTXω−XTY=0−−−−−−（1.2−6）∂E(ω0,ω1)∂ω=XTXω−XTY=0−−−−−−（1.2−6）

 

【矩阵求导： ∂(XTAX)∂X=AX+ATX ∂(XTAX)∂X=AX+ATX；∂(XTA)∂X=A∂(XTA)∂X=A；∂(AX)∂X=AT∂(AX)∂X=AT，∂E(ω0,ω1)∂ω=∂(ωTXTXω)∂ω−∂(ωTXTY)∂ω−∂(YTXω)∂ω∂E(ω0,ω1)∂ω=∂(ωTXTXω)∂ω−∂(ωTXTY)∂ω−∂(YTXω)∂ω，设第一项中XTX=AXTX=A，第二项中XTY=BXTY=B，第三项中YTX=CYTX=C，所以∂E(ω0,ω1)∂ω=[XTXω+(XTX)Tω]−XTY−(YTX)T∂E(ω0,ω1)∂ω=[XTXω+(XTX)Tω]−XTY−(YTX)T，所以，∂E(ω0,ω1)∂ω=2XTXω−2XTY∂E(ω0,ω1)∂ω=2XTXω−2XTY

1 详细推导

 

∂(XTAX)∂X∂(XTAX)∂X

 

 

 

XTAX=∑i=1N∑j=1NAijXiXjXTAX=∑i=1N∑j=1NAijXiXj

 

 

∂(XTAX)∂Xk=∂∂Xk∑i=1N∑j=1NAijXiXj=∑i=1NAikXi+∑j=1NAkjXj=ATX+AX∂(XTAX)∂Xk=∂∂Xk∑i=1N∑j=1NAijXiXj=∑i=1NAikXi+∑j=1NAkjXj=ATX+AX

 

2 详细推导 ∂(AX)∂X∂(AX)∂X

AXAX是n×1n×1维的，第ii个元素为

 

 

[AX]i=∑j=1NAijXj[AX]i=∑j=1NAijXj

 

 

 

AX=[∑j=1NA1jXj,∑j=1NA2jXj,⋯,∑j=1NAnjXj]AX=[∑j=1NA1jXj,∑j=1NA2jXj,⋯,∑j=1NAnjXj]

 

 

 

∂(AX)∂X=[∂∑Nj=1A1jXj∂X,∂∑Nj=1A2jXj∂X,⋯,∂∑Nj=1AnjXj∂X]∂(AX)∂X=[∂∑j=1NA1jXj∂X,∂∑j=1NA2jXj∂X,⋯,∂∑j=1NAnjXj∂X]

 

 

 

∂∑Nj=1AkjXj∂X=⎡⎣⎢⎢⎢⎢Ak1Ak2⋮AkN⎤⎦⎥⎥⎥⎥∂∑j=1NAkjXj∂X=[Ak1Ak2⋮AkN]

】

 

所以线性模型所对应的最优参数ωω表示如下，也称为正则解或者闭形式解

 

 

ω=(XTX)−1XTY−−−−−（1.2−7）ω=(XTX)−1XTY−−−−−（1.2−7）

 

一个简单的例子，代码见文件夹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]

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8

第二步：用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

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11

第三步：将第二步中计算的ωω和第一步中生成的数据，传递给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
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16

 

 
图1.2.2 线性回归的例子

 

在图1.2.2中散点代表的是训练数据，训练数据是由程序随机生成，没有实际意义，直线是回归直线，并标出了直线方程，在运行程序时直线结果可能与图中的结果稍有不同，因为训练数据是随机生成的缘故。

2 最优化方法-梯度下降法

在第一节的第3部分介绍了，将损失函数表示成矩阵形式，然后求导方法，求出最优的ww，这种方法对线性问题可以求出最优解，称为闭形式解或者解析解。本节介绍的梯度下降法是数值最优化方法，普适性更强，对于非线性问题依然可以求解。

梯度下降法是最常用、也是最容易理解的最优化方法。学会了梯度下降法，其它基于梯度的改进方法：共轭梯度法、牛顿法、拟牛顿法等，就比较容易理解。

1 盲人是如何下山的

第一步：左踩一脚，右踩一脚，如果发现这两脚在在高度上没有差别，此时他所面对的应该是山顶或者山脚，反之盲人面对的应该是山脊。（计算偏导数）

第二步：上踩一脚，下踩一脚，脚低的那个方向就对着山脚。（计算偏导数）

第三步：四个脚中，高度最低的那个方向就是山脚，从当前位置向下夸一小步，向着山脚进发。（确定步长，学习率）

重复第一、二、三步，直到山脚。

2 梯度下降法

梯度法就和盲人下山类似，就两个步骤：首先确定下山方，然后再确定的方向上按照一定的步长下山。

下面介绍最优化问题。

单目标、无约束、多维最优化问题的数学描述：

 

 

minxf(x)minxf(x)

 

其中，x∈RNx∈RN。

梯度下降法算法流程如下：

1）给定初值x(0)x(0)，精度ε>ε>0，并令k=1k=1。

2）计算梯度下降方向（搜索方向）v(k)=−∇f(x(k))，∇f(x(k))表示f(x)在x(k)处的梯度计算梯度下降方向（搜索方向）v(k)=−∇f(x(k))，∇f(x(k))表示f(x)在x(k)处的梯度。

【∇f(x(k))∇f(x(k))所表示的是数值梯度，求法如下：

 

 

∇f(x(k))=⎡⎣⎢⎢g1⋮gN⎤⎦⎥⎥∇f(x(k))=[g1⋮gN]

 

 

 

gi=f(x(k)1,⋯,x(k)i+,⋯,x(k)N)−f(x(k)1,⋯,x(k)i,⋯,x(k)N)(x(k)i+)−x(k)igi=f(x1(k),⋯,xi(k)+,⋯,xN(k))−f(x1(k),⋯,xi(k),⋯,xN(k))(xi(k)+)−xi(k)

 

其中，=0.000001=0.000001】

3）若∣∣v(k)∣∣≤ε|v(k)|≤ε，则停止计算，否则从x(k)x(k)出发，沿着v(k)v(k)一维搜索，即求λkλk，使的f(x(k)+λkv(k))=minλ>0f(x(k)+λkv)f(x(k)+λkv(k))=minλ>0f(x(k)+λkv)，此处的一维搜索可以用黄金分割法或者二次差值法等。

4）令x(k+1)=x(k)+λkv(k)x(k+1)=x(k)+λkv(k)，k=k+1k=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_

• 1
• 2
• 3
• 4
• 5

第三步：完成数值梯度的定义，按照梯度下降法中的第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

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8

第四步：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

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29

第五步：将闭形式的ωω和梯度下降法的ωω，以及数据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()

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26

 

 
图2.1 解析解与梯度下降法解的对比图

 

3 基函数

3.1 多项式回归

如果有如图2.1的数据，依然采用式（1.2-1）的模型，则回归模型如图2.2。从图2.2中可以看出，用式（1.2-1）所表示的模型无法拟合这种带多个峰的数据。一个很直观的想法是增加式（1.2-1）中的项数，用多项式拟合这种多个峰的数据

 

 

t=ω0+ω1x+ω2x2+ω3x3+⋯−−−−−（3.1−1）t=ω0+ω1x+ω2x2+ω3x3+⋯−−−−−（3.1−1）

式（3.1）写成矩阵形式为

 

 

 

t=[ω0,ω1,⋯,ωK]⎡⎣⎢⎢⎢⎢1x⋮xK⎤⎦⎥⎥⎥⎥−−−−−（3.1−2）t=[ω0,ω1,⋯,ωK][1x⋮xK]−−−−−（3.1−2）

 

按照1.2节中的方法，也可以得到

 

ω=[ω0,ω1,⋯,ωK]ω=[ω0,ω1,⋯,ωK]

的解，这里就不做详细推导，直接给出结论：

 

 

 

ω=(X¯TX¯)−1X¯TY−−−−−（3.1−3）ω=(X¯TX¯)−1X¯TY−−−−−（3.1−3）

 

为了与式（1.2-7）加以区别，用X¯X¯代替了式（1.2-7）中的XX，这里的X¯X¯表示为如下形式

 

 

X¯=⎡⎣⎢⎢⎢⎢⎢11⋮1x1x2⋮xNx21⋯xK1x22⋯xK2⋮⋯⋮x2N⋯xKN⎤⎦⎥⎥⎥⎥⎥−−−−−（3.1−4）X¯=[1x1x12⋯x1K1x2x22⋯x2K⋮1⋮xN⋮⋯⋮xN2⋯xNK]−−−−−（3.1−4）

 

其中，KK为多项式中自变量的最高次数。

 

 
图3.1.1 x+0.3∗sin(2∗pi∗x)加上随机噪声的数据x+0.3∗sin(2∗pi∗x)加上随机噪声的数据

 

 

 
图3.1.2 按照式（3.1-3），分别用3、5、10、12阶多项式拟合数据，结果如图2.3。图（d）中的拟合曲线的末端上翘与数据不吻合了。这是过拟合导致的。过拟合问题将会在下一节中介绍。

 

实现图3.1.3的代码如下：

第一步：生成x+0.3∗sin(2∗pi∗x)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]

• 1
• 2
• 3
• 4
• 5

第二步：按照式（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

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14

第三步：用第三步中的参数ωω画出拟合的

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

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17

 

 
图3.1.3 多项式拟合结果

 

3.2 回归模型中的基函数

式（3.1-1）更一般化的表示为 

 

t=ω0+∑j=1K−1ωjϕj(x)−−−−−（3.2−1）t=ω0+∑j=1K−1ωjϕj(x)−−−−−（3.2−1）

 

其中的 ϕj(x)ϕj(x) 称为基函数，引入基函数是为了对数据进行非线性变换，以解决非线性问题。

多项式回归中的 x,x2,⋯,xKx,x2,⋯,xK 都可看成是基函数。

其它形式的常用基函数有：高斯基函数，逻辑蒂斯基函数，它们的表达式分别如下

高斯基函数： 

 

ϕj(x)=exp(−(x−μj)22s2)−−−−（3.2−2）ϕj(x)=exp⁡(−(x−μj)22s2)−−−−（3.2−2）

 

其中，μjμj 控制着基函数的位置，ss控制着基函数的形状。

 

 

ϕj(x)=σ(x−μjs)−−−−−（3.2−2）ϕj(x)=σ(x−μjs)−−−−−（3.2−2）

 

其中 σ(a)=11+exp(−a)σ(a)=11+exp(−a)

4 欠拟合与过拟合

4.1 欠拟合

忽略严格的数学定义，从一般的直观理解欠拟合，概念如下。

欠拟合：模型过于简单，无法捕获数据中所存在的规律，图3.1.2所示的情况就是欠拟合，因为采用的模型为t=ω0+ω1xt=ω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(ω)E(ω) 上加上罚项，对高阶项进行处罚，达到降低高阶项前面的系数 ωiωi，ωiωi 变小，说明模型中 ωiωi 所对应的那项对模型的影响程度就会较低，达到简化模型的目的，常用的正则化后的损失函数为

 

 

E(ω)=12∑n=1N{tn−ωϕ(xn)}2+α2ωTω−−−−−（4.2−1）E(ω)=12∑n=1N{tn−ωϕ(xn)}2+α2ωTω−−−−−（4.2−1）

按照前文介绍的方法，依然可以得到，正则化后的模型参数如下 

 

ω=(αI+X¯TX¯)−1X¯TY−−−−−（4.2−2）ω=(αI+X¯TX¯)−1X¯TY−−−−−（4.2−2）

 

其中，αα 为罚参数，II 为单位阵。αα 越大惩罚越强，如果 αα 非常大，则 ωω 会趋向于0。

按照式（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

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17

图4.2.2是不同正则化参数下的回归曲线，从中可以看出 λλ 越大，惩罚越强，α=0.1α=0.1 时，多项式回归模型已经欠拟合了，对应于多项式中的高次项的系数接近于0了。随着 λλ 变小，惩罚程度减弱，高次项的系数基本上没有变小，见（d）。

 

 
图4.2.2 不同罚参数αα下的拟合曲线

 

5 多元线性回归

以上介绍的是一元回归，如果自变量的个数不止一个，就会要求使用多元回归，多元线性回归最简单的形式如下

 

 

t=ω0+∑j=1Dωjxj−−−−−（5−1）t=ω0+∑j=1Dωjxj−−−−−（5−1）

 

其中DD为自变量的个数。但是这种形式的多元回归模型有很多限制。

所以和一元回归类似，也可以引入基函数的概念，引入基函数后的表达如下

 

 

t=ω0+∑j=1K−1ωjϕj(x)−−−−−（5−2）t=ω0+∑j=1K−1ωjϕj(x)−−−−−（5−2）

 

其中，x=(x1,⋯,xM)Tx=(x1,⋯,xM)T，KK为基函数的个数

同理可以得到模型的参数ωω为

 

 

ω=(ΦTΦ)−1ΦTY−−−−−（5−3）ω=(ΦTΦ)−1ΦTY−−−−−（5−3）

 

其中

 

 

Φ=⎡⎣⎢⎢⎢⎢⎢ϕ0(x1)ϕ0(x2)⋮ϕ0(xN)ϕ1(x1)ϕ1(x2)⋮ϕ1(xN)⋯ϕK−1(x1)⋯ϕK−1(x2)⋱⋯⋮ϕK−1(xN)⎤⎦⎥⎥⎥⎥⎥Φ=[ϕ0(x1)ϕ1(x1)⋯ϕK−1(x1)ϕ0(x2)ϕ1(x2)⋯ϕK−1(x2)⋮ϕ0(xN)⋮ϕ1(xN)⋱⋯⋮ϕK−1(xN)]

 

 

 

ϕ0(x)=1ϕ0(x)=1

 

下面看一个二元回归的例子。

有如图5.1所示的曲面，建立一个回归模型拟合这个曲面，由于这个曲面是个二次曲面，所以在选择基函数时可以选择到二次或者更高次，比如选择ϕ1=xϕ1=x，ϕ2=x2ϕ2=x2，此时

 

 

Φ=⎡⎣⎢⎢⎢⎢⎢11⋮1x1x2⋮xNx21x22⋮x2N⎤⎦⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢11⋮1x11x12x21x22⋮xN1xN2x211x212x221x222⋮x2N1x2N2⎤⎦⎥⎥⎥⎥⎥−−−−−（5.4）Φ=[1x1x121x2x22⋮1⋮xN⋮xN2]=[1x11x12x112x1221x21x22x212x222⋮1⋮xN1xN2⋮xN12xN22]−−−−−（5.4）

 

按照式（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函数产生的曲面