机器学习入门学习笔记：（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 - h
weights = weights + alpha * dataMatrix.transpose() * error
return weights

  首先，输入的两个参数dataMatIn和classLabels默认是列表类型，后面要使用numpy进行矩阵运算，所以调用numpy的mat()函数将其转换为numpy下的矩阵数据类型；

  然后预设一些相关的参数，步长不宜太大，迭代次数设为500次；

  接着就是迭代了，不断更新权值，使用如下公式：



  上面这个是梯度上升的公式，如果是梯度下降的公式把加号改成减号就行。

  关于算法推导这里不做赘述，我以前在其他博客中也有有关梯度下降法的介绍。梯度下降法和梯度上升法本质上是一样的。（四旋翼姿态解算——梯度下降法理论推导 ）

  算法不难，也可以自行查阅相关资料。

  只是看公式的话还是有点懵，梯度上升法的伪代码如下：

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

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



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

画出边界

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

# 画出数据集和logistic回归的最佳你和直线的函数
def plotBestFit(weights):
import matplotlib.pyplot as plt
dataMat, 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.01
weights = ones((n, 1))
for i in range(m):
h = sigmoid(sum(dataMatrix[i] * weights))
error = labelMat[i] - h
weights = 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.01
randIndex = int(random.uniform(0, len(dataIndex)))
h = sigmoid(sum(dataMatrix[randIndex] * weights))
error = classLabels[randIndex] - h
weights = weights + alpha * dataMatrix[randIndex].transpose() * error
del(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 1
else:
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 = 0
numTestVec = 0
for line in frTest.readlines():
numTestVec += 1
currentLine = 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 += 1
errorRate = float(errorCount) / numTestVec
print('The error rate of this test is : %f'%errorRate)
return errorRate

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

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

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