1. 多元线性回归要解决的问题

如果有下面的训练数据,$x_1$和$x_2$是特征数据,$y$是数据标签,观察到$y$的取值要么为1,要么为0.

那么问题来了, 已知有上面的观察数据,如果有一组新的输入$x_1 = 1.5,x_2=2.4$,要如何预测$y$的值呢?

2. sigmoid函数

要解决上面的预测问题,我们需要找到一个函数,输入参数$x_i,x_2$之后即可获得预测的$y$值.

根据$y$值的特殊性($y \in {0,1}$),有下面的函数非常类似我们要找的函数,这个函数为:
$$g(z) = \frac{1}{1 + e^{-z}} \tag{1}$$
这个函数的图像长成下面这个样子

由函数图像可见,这个函数在$x\lt-5$或者$x \gt 5$就已经无限接近0和1了.

有了这个函数,就与我们的目标接近一大步了.然而,这个函数只接受一个参数$z$,而我们的有两个参数$x_1,x_2$(实际问题中可能有更多的参数).
解决这个问题,只需要再构造一个函数
$$h_\theta(x)= g(\mathbf{\vec{\theta}}^\mathrm{T} \cdot \vec{x}) = \frac{1}{1 + e^{\mathbf{\vec{\theta}}^{\mathrm{T}} \cdot \vec{x}}} \tag{2}$$

这里$\vec{{\theta}}^\mathrm{T}$表示向量$\vec{\theta}$的转置$={\theta_1,\theta_2}$,于是
$$\mathbf{\vec{\theta}}^\mathrm{T} \cdot \vec{x} = \begin{bmatrix} \theta_1 & \theta_2 \end{bmatrix} \cdot \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \theta_1 \cdot x_1 + \theta_2 \cdot x_2 \\ \tag{3}$$

将公式(3)带入(2)可得$h_\theta(x)$是关于$x_1,x_2$的函数.

实际计算时,对于系数向量$\vec{\theta}$我们还会为线性方程增加一个常量项$\theta_0$,因此,对于向量$\vec{x}$也需要增加一列$x_0$参与计算,因为
对应的$\theta_0$是常数项,所有的$x_0$取值=1即可. 因此开头的表格在实际参与计算的时候变成下面的样子.

因此公式(3)变成:
$$z_{\theta} = \mathbf{\vec{\theta}}^\mathrm{T} \cdot \vec{x} = \begin{bmatrix} \theta_0 & \theta_1 & \theta_2 \end{bmatrix} \cdot \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} = \theta_0 \cdot x_0 + \theta_1 \cdot x_1 + \theta_2 \cdot x_2 \tag{4}$$

将(4)带入(2)得:
$$h_\theta(x) = \frac{1}{1 + e^{- \left (\begin{bmatrix} \theta_0 & \theta_1 & \theta_2 \end{bmatrix} \cdot \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} \right )}} \tag{5}$$

$$h_\theta(x) = \frac{1}{1 + e^{-(\theta_0 \cdot x_0 + \theta_1 \cdot x_1 + \theta_2 \cdot x_2)}} \tag{6}$$

然而,有了公式(5)还是不能愉快的求出要预测的$y$,因为(5)还有参数$\theta_0,\theta_1,\theta_2$是未知的.
接下来看如何求解参数$\theta_0,\theta_1,\theta_2$

3. 极大似然估计

函数$h_\theta(x)$的值有特殊的含义,如果我们把上文表格中的每一行数据看成一个事件,则$h_\theta(x)$的值就是条件概率$P(1|\vec{x};\theta)$,
于是有:
$$P(1|\vec{x};\theta) = h_\theta(x) \tag{7}$$
因为y取值非0即1,所以又有
$$P(0|\vec{x};\theta) = 1 – h_\theta(x) \tag{8}$$

综合起来就有:
$$P(y|\vec{x};\theta) = [h_\theta(\vec{x})]^{y} \cdot [1 – h_\theta(\vec{x})]^{1-y} \tag{9}$$

根据公式(9)就可以计算上文表格中每一行数据出现的概率.假设每一行训练数据之间是独立的,那么三行数据同时出现的概率为
$$L(\vec{\theta}) = p = \prod_{i=1}^3 [h_\theta(\vec{x}^{(i)})]^{y(i)} \cdot [1 – h_\theta(\vec{x}^{(i)})]^{1-y(i)} \tag{10}$$

因为表中的数据已经实实在在的发生了,所以应用最大似然估计原理,需要使得使用公式(10)计算出的概率应该尽可能的大.
因此可求的式(10)中的参数$\vec{\theta}$的最优解,即这一组参数能使得(10)式取得最大值.
为简化计算,(10)式两边取自然底数的对数得到:

$$l(\vec{\theta}) = \ln{L(\vec{\theta})} = \sum_{i=1}^3 y^{(i)}\ln{h_\theta(\vec{x}^{(i)})} + (1-y^{(i)})\ln{[1-h_\theta(\vec{x}^{(i)})]} \tag{11}$$

公式(11)是$\vec{\theta}$的函数,现在需要求出一组$\theta_0,\theta_1,\theta_2$使得函数能取得最大值即可.注意,我们并不关心这个最大值式多少,也不会去计算这个最大值.

4. 梯度上升法求系数$\vec{\theta}$的最优解

先举一个简单的例子说明梯度上升算法的思想.

假如有一个函数$f( \theta ) =\theta_0 + a \theta_1^3 + b \theta_2^2$,其中$a,b$是常数,求函数取得最大值时$\theta_0,\theta_1,\theta_2$的值,然后可以将它们这个函数求出最大值.

1. 首先任取一组初始的$\theta_0,\theta_1,\theta_2$值,不妨去$(1,1,1)$,这里函数要在这个的初始点可微.
2. 求出这一点的梯度,分别对$\theta_0,\theta_1,\theta_2$求偏导数,可得梯度为$(1,3a\theta_1^2,2b\theta_2)$,梯度是一个方向,当$\theta_0,\theta_1,\theta_2$
   往这个方向移动时,获得一组新的$\theta_0,\theta_1,\theta_2$,将这组值带入函数后,求的的函数的值是增大的.假如每次迭代时,朝梯度方向移动$\alpha$(很小的步长)个单位后,
   得到的新的一组
   $$ \begin{aligned} \theta_0 &:=\theta_0 + 1 &\cdot \quad \alpha \\ \theta_1 &:= \theta_1 + 3a\theta_1^2 & \cdot \quad \alpha \\ \theta_2 &:= \theta_2 + 2b\theta_2 &\cdot \quad \alpha \end{aligned} \tag{12}$$
3. 重复上面的步骤,直到达到迭代次数或者符合要求的精度

同理,回到我们现在的问题,要使得上述(11)式取得最大值,需要依次对(11)式求$\;{\theta_0,\theta_1,\theta_2}\;$的偏导数,然后得到梯度.

根据复合函数的求导法则(同济六版高等数学上册92页)有:

如果$u=g(x)$在点 $x$ 可导,而$y=f(u)$在点 $u=g(x)$ 可导, 则复合函数$y=f[g(x)]$在点 $x$ 可导,且其导数为
$$\frac{\mathrm{d}y}{\mathrm{d}x} = f^\prime(u) \cdot g^\prime(x) \;或\;\frac{\mathrm{d}y}{\mathrm{d}x} = \frac{\mathrm{d}y}{\mathrm{d}u} \cdot \frac{\mathrm{d}u}{\mathrm{d}x}$$

这也是传说中的链式求导法则.同时求偏导$\frac{\partial l_{(\theta)}}{\partial \theta_1}$时将$\theta_0,\theta_2$看成常数,再根据一元复合函数求导规则求导即可
(同济六版高等数学下册64页).

$$\frac{\partial l({\vec{\theta}})}{\partial \theta_{j}} = \frac{\partial l({\vec{\theta}})}{\partial g(\vec{\theta}^\mathrm{T} \cdot \vec{x})} \quad {\cdot} \quad \frac{\partial g(\vec{\theta}^\mathrm{T} \cdot \vec{x})}{\partial \vec{\theta}^\mathrm{T} \cdot \vec{x}} \quad \cdot \quad \frac{\partial \left ( \vec{\theta}^\mathrm{T} \cdot \vec{x}\right )}{\partial \theta_j} \tag{13}$$

因为(11)式是求和,而求和的每一项i表示一组训练数据,为简化计算我们暂时不写求和符号$\sum$和$i$

下面分别(13)式等号右边三部分:

• 第1部分

$$\frac{\partial l({\vec{\theta}})}{\partial g(\vec{\theta}^\mathrm{T} \cdot \vec{x})} = y \cdot \frac{1}{g(\mathbf{\vec{\theta}}^\mathrm{T} \cdot \vec{x})} + (1- y) \cdot \frac{1}{1 – g(\mathbf{\vec{\theta}}^\mathrm{T} \cdot \vec{x})} \tag{14.1}$$

• 第2部分,将$(\vec{\theta}^\mathrm{T} \cdot \vec{x})$看成一个整体由$z$表示,这里转化成求(1)式的导数$g(z) = \frac{1}{1 + e^{-z}}$,容易求得,
  $g^\prime(z) = g(z) \cdot (1 – g(z))$,可见$g$的导数可以由它自己表示,可简化后面的计算,于是有:

$$\frac{\partial g(\vec{\theta}^\mathrm{T} \cdot \vec{x})}{\partial \vec{\theta}^\mathrm{T} \cdot \vec{x}} = g(\vec{\theta}^\mathrm{T} \cdot \vec{x}) \cdot \left (1 – g(\vec{\theta}^\mathrm{T} \cdot \vec{x}) \right ) \tag{14.2}$$

• 第3部分

$$\frac{\partial \left ( \vec{\theta}^\mathrm{T} \cdot \vec{x}\right )}{\partial \theta_j} = \frac{\partial (\theta_0 \cdot x_0 + \theta_1 \cdot x_1 + \theta_2 \cdot x_2)}{\partial \theta_j} = x_j \tag{14.3}$$

这里j取值为(0,1,2).

接下来,将(14.1),(14.2),(14.3)带入(13)式,并加上求和符号$\sum$和$i$(第$i$行训练数据),化简后求的偏导数为:

$$\frac{\partial l({\vec{\theta}})}{\partial \theta_{j}} =\sum_i^3 \left (y^{(i)} – h\theta (x^{(i)}) \right ) \cdot x_j \tag{15}$$

终于,我们得到了梯度上升迭代公式:

$$\theta_j := \theta_j + \alpha \cdot \sum_i^3 \left (y^{(i)} – h\theta (x^{(i)}) \right ) \cdot x_j \tag{16}$$

NOTE: 这里每一个带$i$标的表示从训练数据第$i$行拿到对应的数据参与计算,而$h\theta (x^{(i)})$时根据上一次计算出来的${\theta_0,\theta_1,\theta_2}$
带入(5)式之后计算的值,$\alpha$是根据实际需要自定的一个迭代步长,所以(16)式中每一个值都算出来,接下来就可以进行迭代计算了.

至此,终于理解了78页中的一行代码了. 感觉下面这行代码是这一章(第5章)中最核心的东西,而然作者对这一行代码的解释竟然只字未提,呵呵...

weights = weights + alpha * dataMatrix.transpose() * error

5. 实战代码

#!/usr/bin/env python3

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import sklearn.datasets as skDs
import sklearn as skl

def main():
 print(__file__)
 dataArr, labelMat = loadDataSet()
 # print(np.array(dataArr))
 # print(np.array(labelMat))
 # inX = np.mat([1, 2, 3])
 # print(sigmoid(inX))
 dataMatIn, classLabels = loadDataSet()
 weights = gradAscent(dataMatIn, classLabels)
 print(weights)
 plotBestFit(weights.getA())
 weights2 = stocGradAscent1(dataMatIn, classLabels)
 plotBestFit(weights2)

def loadDataSet():
 dataMat = []
 labelMat = []
 fr = open('testSet.txt')
 for line in fr.readlines():
 lineArr = line.strip().split()
 dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
 labelMat.append(int(lineArr[2]))
 return dataMat, labelMat

def sigmoid(inX):
 return 1.0 / (1 + np.exp(-inX))

def gradAscent(dataMatIn, classLabels):
 dataMatrix = np.mat(dataMatIn)
 labelMat = np.mat(classLabels).transpose()
 m, n = np.shape(dataMatrix)
 alpha = 0.001
 maxCycles = 500
 weights = np.ones((n, 1))
 for k in range(maxCycles):
 h = sigmoid(dataMatrix * weights)
 error = (labelMat - h)
 weights = weights + alpha * dataMatrix.transpose() * error
 return weights

def stocGradAscent0(dataMatrix, classLabels):
 m, n = np.shape(dataMatrix)
 alpha = 0.01
 weights = np.ones(n) # initialize to all ones
 for i in range(m):
 h = sigmoid(sum(dataMatrix[i] * weights))
 error = classLabels[i] - h
 weights = weights + alpha * error * dataMatrix[i]
 return weights

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
 m, n = np.shape(dataMatrix)
 weights = np.ones(n) # initialize to all ones
 for j in range(numIter): # 迭代预定的次数
 dataIndex = list(range(m)) # 获取一个训练数据集行标的数组
 for i in range(m):
 alpha = 4 / (1.0 + j + i) + 0.0001
 # apha decreases with iteration, does not go to 0 because of the constant
 randIndex = int(np.random.uniform(0, len(dataIndex)))
 h = sigmoid(sum(dataMatrix[randIndex] * weights))
 error = classLabels[randIndex] - h
 weights = weights + alpha * error * np.array(dataMatrix[randIndex])
 # 这里跟书上不一样, 改成 np.array(dataMatrix[randIndex]
 # 不然报错 can't multiply sequence by non-int of type 'numpy.float64'
 # 应该是Python和numpy的版本跟书中用的不一样了
 del (dataIndex[randIndex])
 return weights

def plotBestFit(weights):
 dataMat, labelMat = loadDataSet()
 dataArr = np.array(dataMat)
 n = np.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 = np.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()

if __name__ == '__main__':
 main()

6. 代码下载

机器学习实战Ch05 多元线性回归 2.55 KB

Download

本文固定链接: https://www.jack-yin.com/coding/ml/3143.html | 边城网事