【西瓜书笔记】——第三章：线性模型

1、基本形式

定义：最简单的线性函数为：
$$f(x)=wx+b$$
将其扩展为矩阵形式，其中$$x_i=(x_{i1},x_{i2},x_{i3},\dots ,x_{ij}{)}^T$$表示对于单个样本xi，它拥有j个不同的特征，而针对每个特征的重要度不同，自然有权重矩阵：
$$w=(w_1,w_2,w_3,\dots ,w_j)$$
因此 wx 结果为一个数，再加上b得到f(xi)的评判结果。

---

2、线性回归

目的是通过线性模型将离散属性通过连续化转化为连续值，在线性模型中确定w和b的值是关键，因此需要衡量f(x)的预估结果于y真实值之间的差别。均方误差法通常被用来进行性能度量，试图让均方误差最小化来求解最合适的w和b值。
$$(w^{\ast },b^{\ast })=argmin(w,b)\sum _{i=1}^m(f(x_i)-y_i{)}^2=argmi{n}_{(w,b)}\sum _{i=1}^m(y_i-w{x}_i-b{)}^2$$
实为欧式距离计算公式或“最小二乘法”。几何意义为找一条直线使所有样本到直线上的欧式距离之和最小。
令
$$E_{(w,b)}={\sum }_{i=1}^n(y_i-w{x}_i-b{)}^2$$
因此，求解w和b使得E最小化的过程，称为线性回归模型的最小二乘参数估计。可以利用求导的方式解得w和b（结果见书p54页）。

对于多元线性回归，将原式改为矩阵形式可得：
$$f(x_i)=w^T{x}_i+b即：{\hat{w}}^{\ast }=argmi{n}_{\hat{w}}(y-X\hat{w}{)}^T(y-X\hat{w})同样的，令E为0，对w求导可知：当X^{T}X为满秩矩阵或正定矩阵时，令E为0可得{\hat{w}}^{\ast }=(X^{T}X{)}^{-1}{X}^{T}y$$
但是现实任务种X^TX往往不是满秩矩阵，此时可解出多个w，它们都能使均方误差最小化，因此选择哪一个解作为输出将由学习算法的归纳偏好决定，常见的作法是引用正则化。

---

2.1 正则化

向模型加入某些规则，加入先验知识缩小解空间，减小求出错误解的可能性。将所用哟䣌数学知识数字化，告诉模型，对代价函数来说就是加入对模型“好坏”的评判标准。

数学解释
1、通俗定义
就是给平面不可约代数曲线以某种形式的全纯参数表示。
即
$$对于P{C}^2中的不可约代数曲线C，寻找一个紧Riemann面C^{\ast }和一个全纯映射\sigma :C^{\ast }\to P{C}^2,使得\sigma (C^{\ast })=C$$

2、严格定义
设C是不可约平面代数曲线，S是C的奇点的集合。
$$如果存在紧Riemann面C^{\ast }及全纯映射\sigma :C^{\ast }\to P{C}^2$$
使得
$$\sigma (C^{\ast })=C$$
$${\sigma }^{-1}(S)是有限点集$$
$$\sigma :\frac{C^{\ast }}{{\sigma }^{-1}}$$
(S)→C\S是一对一的映射
则称(C*,σ)为C的正则化。不至于混淆的时候，也可以称C*为C的正则化。正则化的做法，实际上是在不可约平面代数曲线的奇点处，把具有不同切线的曲线分支分开，从而消除这种奇异性。

---

2.2 从狭义模型到广义

线性模型虽然简单，但是却又很多的变换形式。常用的指数尺度变化模型为：
$$y=w^{T}x+b$$
也可以利用对数的形式将线性回归模型的预测值与真是标记联系起来：
$$lny=w^{T}x+b，其本质上也是线性回归，只是将y转化为了y=e^{w^{T}x+b}的形式$$
为了更一般的考虑，也可以利用单调可微函数g(·)：
$$y=g^{-1}(w^{T}x+b)$$
这个被称为广义线性模型，其中g(·)被称为“联系函数”，很显然，ln(·)是g(·)的一个特例。

---

3、对数几率回归

在处理分类任务的情况时，只需找一个单调可微函数将分类任务的真实标记y与线性回归模型的预测值联系起来。

3.1 单位阶跃函数

考虑二分类任务，其输出标记y∈{0, 1}，而线性回归模型产生的预测值z=wx+b是实值，于时，我们需要将实值z转换为0/1值。最理想的方式就是单位阶跃函数。
$$y=\{\begin{array}{l}0,z<0;\\ 0.5,z=0;\\ 1,z>0\end{array}$$
若预测值z大于0则判断为正例，小于零则为反例，预测值为临界值则可以任意判别，这也被称为Heaviside函数。

但是阶跃函数不连续，我们希望找到一个一定程度上近似单位阶跃函数的替代函数：
$$y=\frac{1}{1+e^{-z}}$$
并希望它可以微调，这便是更多情况下选用的对数几率函数。

3.2对数几率函数

对数几率函数就是我们常说的Sigmoid函数种的一种，它将z值转化为一个接近0或1的y值，并且其输出值在z=0附近变化很显著。

---

Logistic回归

逻辑回归的主要作用是分类，其思想是：根据现有数据对分类边界线建立回归共识，以此进行分类。旨在找到最佳拟合参数集，所运用的方法便是常见的一些最优化算法。

逻辑回归的一般过程：

1、收集数据：采用任意方法收集数据；
2、准备数据：因算法涉及欧氏距离计算，故要求数据类型为数值型（结构化数据格式最好）；
3、分析数据：采用任意方法分析数据；
4、训练数据：目的是找到最佳拟合参数集（回归系数）；
5、测试算法：训练完成后对算法的泛化性能做评估；
6、使用算法：对任意新的待分类数据集进行实际应用。注意对数据进行结构化处理。

此部分将对逻辑回归的基本原理做出解释说明，并介绍相关的最优化算法如梯度上升和随机梯度上升。最后将附上测试代码。

---

1、基于逻辑回归和Sigmoid函数的分类

逻辑回归的特点：
【优点】：计算代价不高，容易理解；
【缺点】：容易出现欠拟合，分类精度不高；
【数据类型】：数值型和标称型。

常用的用于进行二分类的函数有海维赛德阶跃函数和单位阶跃函数，都是用于将输入结果映射为0或1的输出结果。而上述两个函数有各自的缺点，例如瞬间跳跃无法处理、对整个输入数据集的数据不具备普遍操作性等。因此，为了数学上便与处理，Sigmoid函数的出现起到了很重要的作用。

【*注：虽然在神经网络中Sigmoid函数存在梯度弥散的问题。但是此处不对双曲正切Tanh和ReLU做过多说明】
$$\sigma (z)=\frac{1}{1+e^{-z}}$$
Sigmoid函数的特点在于其在x取0时，S(x)=0.5。因此，以0.5作为阈值出现了x<0 → 0而x>0 → 1的情况。

为了实现Logistic回归，我们可以在每个特征上乘上一个回归系数，然后把所有的结果值相加，将总和代入Sigmoid函数中，从而得到一个范围在0~1间的数值。这也是逻辑回归常被视为一种概率估计的原因。

显然，在函数模型确定的情况下，回归系数是多少才是逻辑回归的关键。

---

2、基于最优化方法的最佳回归系数确定

首先需要明白一点：逻辑回归的输入并不是原始的样例特征值输入。
$$input={\theta }^{T}X$$
$${\hat{Y}}_{predict}=sigmoid(input)=\frac{1}{1+e^{-input}}$$
其中，θ就是回归系数矩阵，而X就是原始的样例输入集合。为了得到最佳拟合系数矩阵θ，我们常用梯度上升法这一最优化方法，利用python可以绘制决策边界图，用以直观地秒回梯度与分类精度之间的变化关系。

---

2.1 梯度上升法

【概念】：要找到某函数的最大值，最好的方法就是沿着该函数的梯度方向探寻。即对函数求导。（沿不同方向便是分别对不同方向的参数求偏导）梯度算子在图像中（梯度方向上）总是指向数值增长最快的方向，即与该点处的切线始终垂直，而每一次梯度移动的长度被称为步长α。

由此概念我们可以得到梯度迭代公式：
$$w:=w+\alpha {\Delta }_{w}f(w)$$
对系数函数求导*步长+原系数=新系数值。该公式将一直被迭代直到达到终止条件，例如某一个指定数值或某一个指定误差范围。

【*注】梯度上升的公式中是“＋”，用于求函数最大值；而梯度下降公式中是“－”，用于求最小值。

---

2.2 实战：使用梯度上升找到最佳参数

样本空间为100个样例，有两个特征A1和A2。我们将对该数据集进行测试，利用梯度上升法找到最佳回归系数，也就是拟合出最佳的Logistic回归模型。
伪代码如下：

每个回归系数初始化为1
重复R次：
	计算整个数据集的梯度
	使用“步长*梯度”来更新回归系数
返回回归系数

程序清单
logRegres.py

import numpy as np

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 = alpha*dataMatrix.transpose() * error + weights
    return weights

这里设置了步长α=0.001，设置了迭代次数maxCycles=500。而梯度上升公式也变成了：
$$w_{new}:=\alpha X^{T}Err+w_{original}$$
这里的公式推导如下：

---

2.3 决策边界的可视化

程序清单
logRegres.py

def plotBestFit(weights):
    import matplotlib.pyplot as plt
    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, c='red', s=30, marker = 's')
    ax.scatter(xcord2, ycord2, c='green', s=30)
    x = np.arange(-3.0, 3.0, 0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]
    x = x.reshape(np.shape(y))
    ax.plot(x, y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

这里预设了sigmoid的阈值为0。本代码经过修改为python3版本，经过可视化操作发现效果不错，但是计算成本仍然很高。

---

2.4 梯度随机上升法

在样本空间规模不够庞大时，使用梯度上升算法是很有效的，但是当有数十亿样本和成千上万特征时，该方法的计算复杂度就显得很高。因为梯度上升算法每一次更新回归系数时都需要遍历整个数据集。

随机梯度上升法：一次仅用一个样本点来更新回归系数。同时也是一个在线学习算法（在新样本到来时对分类器进行增量式更新），伪代码如下：

所有回归系数初始化为1
对数据集中每个样本
		计算该样本的梯度
		使用α*Δf(x)的方式来更新回归系数值
返回回归系数值

程序清单
logRegres.py

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

同梯度上升算法相比，随机梯度上升法的h、error都是由向量形式而非数值形式表现；且在计算权重的同时没有对dataMatrix做转置处理。经过可视化处理后发现，该方法在此样本空间中的分类精度似乎不及梯度上升法。但事实上并非如此，在随机梯度上升法中，当迭代次数增多时，回归系数的迭代也会趋于稳定，这时再将回归系数带入待定函数中进行求解效果会更好。

当然，随机梯度下降法也可以再做一定的改进：
程序清单：
logRegres.py

def stocGradAscent1(dataMatrix, classLabels, numIter = 150):
	m, n = np.shape(dataMatrix)
	weights = np.ones(n)
	for j in range(numIter):
		dataIndex = list(range(m))
		for i in range(m):
			alpha = 4/(1.0 + j + i) + 0.01
			randIndex = int(np.random.uniform(0, len(dataIndex)))
			h = sigmoid(sum(dataMatrix[randIndex] * weights))
			error = classLabels[randIndex] - h
			weights = weights + alpha * error * dataMatrix[randIndex]
			del(dataIndex[randIndex])
	return weights

显然，上述代码添加了alpha = 4/(1.0 + j + i) + 0.01的部分来设置步长的自动迭代，同时也避免了步长随着迭代次数的增多而减少到0（因为末尾有常数项0.01）。这样的做法是为了避免数据陷入梯度中而无法走出的情况，即使梯度达到了最小，函数也能从中”走“出来。关于α，j是迭代次数而i是样本下标。

同时，randIndex的设置可以避免周期性的波动（因为一般的随机梯度上升方法在α稳定前后会经历或大或小的波动甚至是周期性波动详见书82-83页图5-6进而5-7）。这样做可以节省计算时间和步骤。

---

3、实例：疝气疾病预测病马死亡率

【X = 368, Attribute = 28】
回顾一下，首要任务是找到合适的回归系数，然后来使用算法。

---

3.1 处理数据中的缺省值

数据的价值往往超乎想象，有时候面对缺省的数据集，用之不妥，弃之可惜，这时候往往就需要对缺省值进行处理，而我们处理缺省值的目的就是充分利用数据，尽可能弥补数据的不足，常用的处理数据缺省的方法有：

1、用缺省值所在特征的均值填补；
2、用特殊值填补，例如-1；
3、忽略该样本；
4、用相似样本的值来填补此处的值；
5、换一种机器学习算法

再用numpy处理时，必须保证填补的缺省值为实数值，这里的logistic回归公式为
$$weights=weights+alpha\ast error\ast dataMatrix[randIndex]$$
本实例中采用0来填补缺省值，这是因为sigmoid函数本身的优点。另外，如果在数据中发现缺省的是标签值，那么该样例就一定要舍弃，因为无法对利用监督学习模型分析非监督学习的样例。

---

3.2 用logistic回归进行分类

程序清单
logRegres.py

def classifyVector(inX, weights):
     prob = sigmoid(sum(inX * weights))
     if prob > 0.5: 
         return 1.0
     else: 
         return 0.0

def colicTest():
    frTrain = open('horseColicTraining.txt')
    frTest = open('horseColicTest.txt')
    trainingSet = []; trainingLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    trainWeights = stocGradAscent1(np.array(trainingSet), trainingLabels, 500)
    errorCount = 0; numTestVec = 0.0
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(np.array(lineArr), trainWeights))!=int(currLine[21]):
            errorCount += 1
    errorRate = (float(errorCount)/numTestVec)
    print("the error rate of this test is: %f" % errorRate)
    return errorRate

def multiTest():
    numTests = 10; errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print("after %d iterations the average error rate is: %f" % (numTests, errorSum/float(numTests)))

---

小结

Logistic回归的目的是寻找一个非线性函数Sigmoid的最佳拟合参数，求解过程由最优化方法（梯度）来完成。而适用于本处的梯度上升法又可以改进为随机梯度上升法。

而关于数据的预处理阶段，此处仅涉及到了对缺省值的处理，有不同的处理方法。同样在将来的数据预处理实验中还会涉及到对数据的清洗，包括对噪声数据的处理，这将在接下来的博客中有所涉及。

【*注：从概率角度解释梯度下降】
满足正态分布：$N(\mu ,\sigma )$，其中$f(x)=\frac{1}{2\sigma }{e}^{\frac{(x-\mu )}{2{\sigma }^2}}$。
约束条件：$argma{x}_{w}N(y|w^{T}x,{\beta }^{-1})N(w|0,{\lambda }^{-1})$
分布满足：$p(w|x,y,\beta )\propto e^{-\frac{\beta }{2}{\sum }_{n=1}^N(y_n-w^T{x}_n{)}^2}\cdot e^{-\frac{\lambda }{2}{w}^{T}w}$
再对w求导。

---

满足拉普拉斯算子：
$$p(y|w^{T}x,{\beta }^{-1})=e^{-\beta |y-w^{T}w|}$$