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

来源：互联网发布：linux jdk设置jvm内存编辑：程序博客网时间：2024/05/24 01:06

1 一元线性回归

1.1 为什么用回归

这里写图片描述
图1.1.1 Google的票房与搜索量的关系

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

1.2 一元线性回归模型

1 数学描述

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

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

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

E (ω 0, ω 1) = 1 2 \sum i = 1 N (t i - y i) 2 = 1 2 \sum i = 1 N (ω 0 + ω 1 x i - y i) 2 - - - - - （ 1.2 - 2 ）

只要最小化式（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个等式，并用线性代数表示为

⎡ ⎣ ⎢ ⎢ t 1 ⋮ t N ⎤ ⎦ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ω 0 + ω 1 x 1 ⋮ ω 0 + ω 1 x N ⎤ ⎦ ⎥ ⎥ = X ω - - - - - （ 1.2 - 3 ）

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

此时

⎡ ⎣ ⎢ ⎢ t 1 - y 1 ⋮ t N - y N ⎤ ⎦ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ω 0 + ω 1 x 1 - y 1 ⋮ ω 0 + ω 1 x N - y 1 ⎤ ⎦ ⎥ ⎥ = X ω - Y - - - - - - （ 1.2 - 4 ）

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

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

E (ω 0, ω 1) = 1 2 (X ω - Y) T (X ω - Y) = ω T X T X ω - ω T X T Y - Y T X ω + Y T Y - - - - - （ 1.2 - 5 ）

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

3 目标函数最小化

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

\partial E ( ω 0 , ω 1 ) \partial ω = X T X ω - X T Y = 0 - - - - - - （ 1.2 - 6 ）

【矩阵求导： ∂(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 详细推导

\partial ( X T A X ) \partial X

X T A X = \sum i = 1 N \sum j = 1 N A ij X i X j

\partial ( X T AX ) \partial X k = \partial \partial X k \sum i = 1 N \sum j = 1 N A ij X i X j = \sum i = 1 N A i k X i + \sum j = 1 N A k j X j = A T X + A X

2 详细推导 ∂(AX)∂X

AX是n×1维的，第i个元素为

[AX] i = \sum j = 1 N A ij X j

A X = ⎡ ⎣ \sum j = 1 N A 1 j X j, \sum j = 1 N A 2 j X j, \dots, \sum j = 1 N A nj X j ⎤ ⎦

\partial ( AX ) \partial X = ⎡ ⎣ \partial \sum N j = 1 A 1 j X j \partial X, \partial \sum N j = 1 A 2 j X j \partial X, \dots, \partial \sum N j = 1 A nj X j \partial X ⎤ ⎦

\partial \sum N j = 1 A kj X j \partial X = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ A k 1 A k 2 ⋮ A kN ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

】

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

ω = (X T X) - 1 X T Y - - - - - （ 1.2 - 7 ）

一个简单的例子，代码见文件夹1_regression。
第一步：用synthic_data.py中的linearSamples方法生成数据

import numpy as npimport randomdef 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 梯度下降法

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

下面介绍最优化问题。

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

m i n x f (x)

其中，x∈RN。

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

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

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

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

\nabla f (x (k)) = ⎡ ⎣ ⎢ ⎢ g 1 ⋮ g N ⎤ ⎦ ⎥ ⎥

g i = f ( x ( k ) 1 , \dots , x ( k ) i + , \dots , x ( k ) N ) - f ( x ( k ) 1 , \dots , x ( k ) i , \dots , x ( k ) N ) ( x ( k ) i + ) - x ( k ) i

其中，=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）中的项数，用多项式拟合这种多个峰的数据

t = ω 0 + ω 1 x + ω 2 x 2 + ω 3 x 3 + \dots - - - - - （ 3.1 - 1 ）

式（3.1）写成矩阵形式为

t = [ω 0, ω 1, \dots, ω K] ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ 1 x ⋮ x K ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ - - - - - （ 3.1 - 2 ）

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

ω = [ω 0, ω 1, \dots, ω K]

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

ω = (X ¯ T X ¯) - 1 X ¯ T Y - - - - - （ 3.1 - 3 ）

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

X ¯ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ 11 ⋮ 1 x 1 x 2 ⋮ x N x 21 \dots x K 1 x 22 \dots x K 2 ⋮ \dots ⋮ x 2 N \dots x K N ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ - - - - - （ 3.1 - 4 ）

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

图3.1.1

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)样本的代码

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）更一般化的表示为

t = ω 0 + \sum j = 1 K - 1 ω j ϕ j (x) - - - - - （ 3.2 - 1 ）

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

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

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

高斯基函数：

ϕ j (x) = exp ⎛ ⎝ ⎜ ⎜ - ( x - μ j ) 2 2 s 2 ⎞ ⎠ ⎟ ⎟ - - - - （ 3.2 - 2 ）

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

ϕ j (x) = σ (x - μ j s) - - - - - （ 3.2 - 2 ）

其中 σ(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 所对应的那项对模型的影响程度就会较低，达到简化模型的目的，常用的正则化后的损失函数为

E (ω) = 1 2 \sum n = 1 N {t n - ω ϕ (x n)} 2 + α 2 ω T ω - - - - - （ 4.2 - 1 ）

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

ω = (α I + X ¯ T X ¯) - 1 X ¯ T Y - - - - - （ 4.2 - 2 ）

其中，α 为罚参数，I 为单位阵。α 越大惩罚越强，如果 α 非常大，则 ω 会趋向于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

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

图4.2.2 不同罚参数

α下的拟合曲线

5 多元线性回归

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

t = ω 0 + \sum j = 1 D ω j x j - - - - - （ 5 - 1 ）

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

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

t = ω 0 + \sum j = 1 K - 1 ω j ϕ j (x) - - - - - （ 5 - 2 ）

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

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

ω = (Φ T Φ) - 1 Φ T Y - - - - - （ 5 - 3 ）

其中

Φ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ϕ 0 (x 1) ϕ 0 (x 2) ⋮ ϕ 0 (x N) ϕ 1 (x 1) ϕ 1 (x 2) ⋮ ϕ 1 (x N) \dots ϕ K - 1 (x 1) \dots ϕ K - 1 (x 2) ⋱ \dots ⋮ ϕ K - 1 (x N) ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥

ϕ 0 (x) = 1

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

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

Φ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ 11 ⋮ 1 x 1 x 2 ⋮ x N x 21 x 22 ⋮ x 2 N ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ 11 ⋮ 1 x 11 x 12 x 21 x 22 ⋮ x N 1 x N 2 x 211 x 212 x 221 x 222 ⋮ x 2 N 1 x 2 N 2 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ - - - - - （ 5.4 ）

按照式（5.4）代入式（5.3）可以求得模型参数ω。

图5.1 二维曲面

下面看看用以上思路能否拟合出图5.1的曲面。

第一步：生成图5.1所示的数据，在x，y方向分别等距离采集数据点40个，Python代码如下：
第二步：用生成的数据点，按照式（5-3）求出参数ω
第三步：用第二步中计算的参数ω，并且用第一步中的方法生成测试数据，在x，y方向上分别采样20个点，用这20个点作为测试数据，输入到模型中，看模型预测的值和真实的值之间的误差。图5.2中，将真实值和预测点按照对应的一行一行的展开，分别做成两个一维向量，这样便于比对。从图中看出，预测的结果还是精确的。