机器学习入门学习笔记：（2.2）线性回归python程序实现

上一篇博客中，推导了线性回归的公式，这次试着编程来实现它。（机器学习入门学习笔记：（2.1）线性回归理论推导）
我们求解线性回归的思路有两个：一个是直接套用上一篇博客最后推导出来的公式；另一个是使用梯度下降法来求解极值。如果数据量很大不建议采用第一个，采用后者能更有效地减小计算量。这篇博客后面的程序也采用的是后者。
事先声明，这篇博客是我自己的学习笔记，代码参考了《机器学习实战》一书，小部分代码有所修改。代码以及数据集我打包上传好了，感兴趣的可以下载下来看看：
线性逻辑回归代码
（今天才发现csdn下载居然没有了0分资源，最少都要1分，有点坑啊）

代码不难，简要介绍下，程序中有注释。
按照顺序将代码加入logRegres.py。

导入数据集

# 导入数据集
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

创建了列表dataMat和labelMat，分别用于存储特征变量对应的值和标签对应的值；
数据集存放在testSet.txt，从中读取数据；
逐行读取，遍历所有行；读取的每一行做一下分割，每行的最后一个数据是标签，前面的都是特征向量的数据。将其分开后，分别加入列表dataMat和labelMat中。

激励函数

# sigmoid函数
def sigmoid(inX):return 1.0 / (1 + exp(-inX))

我在以前的博客中介绍了激励函数的概念。（机器学习入门学习笔记：（1）BP神经网络原理推导及程序实现）
如果只考虑两个类的情况下，我们的结果可以用0或1来表示，这种函数叫做阶跃函数。但是由于阶跃函数存在一个0到1的跳变的过程，这会导致不可导、不光滑，有时在数学上难以处理。所以我们采用sigmoid函数来代替激励函数，以便能在数学上更易处理。
直接给出sigmoid的公式，不做赘述了：

梯度上升法

通常我们更多地接触到的是梯度下降法，这两种方法本质都是一样的。由于梯度是变化最快的方向，梯度下降法就是向最快的方向减小，而梯度上升法就是增大，最后两者分别收敛于局部最小值和局部最大值。
什么时候用梯度下降法或者是梯度上升法呢？具体问题具体分析，如果我们用的是一个误差函数，我们希望这个误差尽可能小，那么就选择梯度下降法了，很简单的逻辑。

# 梯度上升法
def gradAscent(dataMatIn, classLabels):# 转换为numpy矩阵数据类型dataMatrix = mat(dataMatIn)labelMat = mat(classLabels).transpose()# 预设初始参数m, n = shape(dataMatrix)# 求出矩阵大小alpha = 0.001   # 步长maxCycles = 500 # 迭代500次# 使用梯度上升找到最优值weights = ones((n, 1))  # 权值for k in range(maxCycles):h = sigmoid(dataMatrix*weights)error = labelMat - hweights = weights + alpha * dataMatrix.transpose() * errorreturn weights

首先，输入的两个参数dataMatIn和classLabels默认是列表类型，后面要使用numpy进行矩阵运算，所以调用numpy的mat()函数将其转换为numpy下的矩阵数据类型；
然后预设一些相关的参数，步长不宜太大，迭代次数设为500次；
接着就是迭代了，不断更新权值，使用如下公式：

上面这个是梯度上升的公式，如果是梯度下降的公式把加号改成减号就行。
关于算法推导这里不做赘述，我以前在其他博客中也有有关梯度下降法的介绍。梯度下降法和梯度上升法本质上是一样的。（四旋翼姿态解算——梯度下降法理论推导）
算法不难，也可以自行查阅相关资料。
只是看公式的话还是有点懵，梯度上升法的伪代码如下：

初始化权值，使所有回归系数为1
重复n次（n表示迭代次数）： 计算整个数据集的梯度使用alpha（步长）*gradient（梯度）套上面公式更新回归系数的向量
返回回归系数

跟给出的代码中对一对，就很容易搞懂了。在命令行测试一下看看：

另外还有一点，我们看看代码不难发现，可以很轻易地将这段程序改成梯度下降法，也就是将程序中的梯度更新代码的几个符号改一下，结果还是一样的。

画出边界

只是得到了参数还不够直观，我们接着使用matplotlib来是数据可视化。

# 画出数据集和logistic回归的最佳你和直线的函数
def plotBestFit(weights):import matplotlib.pyplot as pltdataMat, labelMat = loadDataSet()dataArr = array(dataMat)n = 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 = 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()

调用库，没什么好说的。无非就是把点全部画出来，最后又画了一条边界。

权值返回来是一个numpy的矩阵类型，我们要先调用getA()方法将其转换为数组类型。

分类结果挺不错的，只有很少的点被分错了。

随机梯度上升

前面使用的梯度上升法也叫作批量梯度上升法(batch gradient ascent)，它每一次更新回归系数时都要遍历整个数据集，如果数据集很大，计算量也会非常大。
那么一个改进方法就是一次仅仅使用一个样本点来更新回归系数，这个方法叫随机梯度上升法。
随机梯度上升法的伪代码：

所有回归系数初始化为1
遍历数据中的每个样本：计算样本的梯度使用alpha（步长）*gradient（梯度）套上面公式更新回归系数的向量
返回回归系数值

编写代码：

# 随机梯度上升算法
def stocGradAscent(dataMatIn, classLabels):# 转换为numpy矩阵数据类型dataMatrix = mat(dataMatIn)labelMat = mat(classLabels).transpose()m,n = shape(dataMatrix)alpha = 0.01weights = ones((n, 1))for i in range(m):h = sigmoid(sum(dataMatrix[i] * weights))error = labelMat[i] - hweights = weights + alpha *dataMatrix[i].transpose() * error  return weights

这部分代码大体上跟前面梯度上升法是一样的，不做赘述。唯一的区别就是，没有了迭代次数，这里是遍历了整个数据集，使用每个样本来对参数进行更新而不是每个样本都计算一遍。

测试一下：

结果不太好，我们还需要修改程序。

改进算法

# 改进的随机梯度上升算法
def stocGradAscent_optimized(dataMatIn, classLabels, numIter=150):# 转换为numpy矩阵数据类型dataMatrix = mat(dataMatIn)labelMat = mat(classLabels).transpose()m, n = shape(dataMatrix)weights = ones((n, 1))for j in range(numIter):dataIndex = range(m)for i in range(n):alpha = 4 / (1.0 + j + i) + 0.01randIndex = int(random.uniform(0, len(dataIndex)))h = sigmoid(sum(dataMatrix[randIndex] * weights))error = classLabels[randIndex] - hweights = weights + alpha * dataMatrix[randIndex].transpose() * errordel(dataIndex[randIndex])return weights

首先，步长alpha不再是一个常数，随着迭代次数的增大不断减小。步长alpha的值后面还加了一个0.01，因为不能让它趋近于0，如果趋近于0了新数据的更新也几乎没有作用，乘积都是0。
然后，这里样本选取不再是按照默认顺序选取了。随机选取数据集进行更新，反复迭代多次。

试试看结果：

效果跟批量梯度下降差不多了，但是随机梯度下降的计算量要少很多。这里使用的默认迭代次数只有150次。

测试算法

从疝气病症预测病马死亡的概率：

分析特征变量

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

以回归系数和输入的特征变量计算对应的sigmoid函数的值，如果大于0.5，就预测为1；如果小于0.5，就预测为0。

测试结果

def colicTest():# 打开测试集和数据集frTrain = open('horseColicTraining.txt')frTest = open('horseColicTest.txt')# 导入训练集的特征变量和标签trainingSet = []trainingLabels = []for line in frTrain.readlines():currentLine = line.strip().split('\t')lineArr = []for i in range(21):lineArr.append(float(currentLine[i]))trainingSet.append(lineArr)trainingLabels.append(float(currentLine[21]))# 使用优化过后的随机梯度上升法计算权重trainWeights = stocGradAscent_optimized(trainingSet, trainingLabels, 1200)# 导入测试集，并检验效果errorCount = 0numTestVec = 0for line in frTest.readlines():numTestVec += 1currentLine = line.strip().split('\t')lineArr = []for i in range(21):lineArr.append(float(currentLine[i]))if int(classifyVector(array(lineArr), trainWeights)) != int(currentLine[21]):errorCount += 1errorRate = float(errorCount) / numTestVecprint('The error rate of this test is : %f'%errorRate)return errorRate

直接调用前面写好的函数了，不多说了。

调用colicTest()10次，比较结果：

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