写在最前面：Logistic回归通过Sigmoid函数接受输入然后进行预测
首先，介绍一下什么是Sigmoid函数。大家一定听过海维赛德阶跃函数（Heaviside step function），什么？没听过，好吧，换个名字，单位阶跃函数，这个认识吧!
这个函数的问题在于该函数在跳跃点上从0瞬间跳跃到1，这样处理起来不是很方便
然鹅我们还有另一个具有类似性质的函数Sigmoid函数，计算公式如下：
f(z) = 1/1+e^-z
当x轴的尺度足够大
import matplotlib.pyplot as plt
import numpy as np
def sigmoid(x):
# 直接返回sigmoid函数
return 1. / (1. + np.exp(-x))
def plot_sigmoid():
# param:起点，终点，间距
x = np.arange(-30, 30, 0.2)
y = sigmoid(x)
plt.plot(x, y)
plt.show()
if __name__ == '__main__':
plot_sigmoid()
Sigmoid函数的输入记为z，给出输入z的公式：
z = w0x0 + w1x1 + w2x2 + w3x3 +…+ wnxn
用向量去理解的话，就是w^T * x，w是一个列向量，大家都知道行向量和列向量相乘得到一个数。其中向量x是分类器的输入数据，向量w是我们要找的最佳参数。
我们通过梯度上升法寻找最佳回归系数。梯度上升法的基本思想是：想要找到某函数的最大值，最好的方法是沿着该函数的梯度放向探寻。
记梯度为▽，则函数f(x,y)的梯度▽f(x,y) =
这个梯度意味着要沿着x的方向移动
然后沿着y的方向移动
这里引入移动量的大小α，用向量表示梯度上升算法的迭代公式如下：
w := w + α▽wf(w) 大家肯定会疑惑这里的 := 什么意思，其实就是编程语言里的赋值，编程里一般用 = ，但是这里用 = 会有歧义。
该公式将会一直迭代，直到达到某个停止条件为止。
好了，基础介绍到这里
下面举个栗子
def loadDataSet(filename):
dataMat = []; labelMat = []
fr = open(filename)
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
该函数的作用是打开文本
-0.017612 14.053064 0
-1.395634 4.662541 1
-0.752157 6.538620 0
-1.322371 7.152853 0
0.423363 11.054677 0
0.406704 7.067335 1
0.667394 12.741452 0
-2.460150 6.866805 1
0.569411 9.548755 0
-0.026632 10.427743 0
0.850433 6.920334 1
1.347183 13.175500 0
1.176813 3.167020 1
.
.
.
每行前两个值分别是x1，x2，第三个值是数据对应的类别标签。
第二个函数是梯度上升算法，该函数有两个参数，第一个参数是dataMatin，是一个二维Numpy数组，每列分别代表每个不同的特征，每行则代表每个训练样本。dataMatin存放的事100x3的矩阵，alpha是目标移动的步长，maxCycless是迭代次数。
由于txt文件中包含x1,x2，我们还将x0设为1.0，那么将三个特征值作为二维矩阵的列，txt中的100个数据作为行，这就是一个100x3的矩阵，第二个参数是一个1x100的矩阵，列是gradAscent的第二个参数classLabels，为了便于计算，我们使用transpose函数进行转置，转成100x1的列向量
labelMat = mat(classLabels).transpose()
然后shape函数计算这个100x3的矩阵的维度，m=100,n=3
m,n = shape(dataMatrix)
步长默认为0.001
alpha = 0.001
迭代次数为500次
maxCycles = 500
ones函数构造了一个全1矩阵，n行1列的矩阵
weights = ones((n,1))
相当于将回归系数全部初始化为1
循环500次，将100x3的矩阵和nx1的矩阵相乘，得到100x1的矩阵，再把这个矩阵作为参数输入Sigmoid函数，error是标签减去函数的输出，这里计算的是真实类别和预测类别的差值，毫无疑问，h是我们的预测类别。按照该差值的方向调整回归系数。
for k in range(maxCycles): #heavy on matrix operations
h = sigmoid(dataMatrix*weights) #matrix mult
error = (labelMat - h) #vector subtraction
weights = weights + alpha * dataMatrix.transpose()* error #matrix mult
下面是这个函数的完整代码
def gradAscent(dataMatIn, classLabels):
dataMatrix = mat(dataMatIn) #convert to NumPy matrix
labelMat = mat(classLabels).transpose() #convert to NumPy matrix
m,n = shape(dataMatrix)
alpha = 0.001
maxCycles = 500
weights = ones((n,1))
for k in range(maxCycles): #heavy on matrix operations
h = sigmoid(dataMatrix*weights) #matrix mult
error = (labelMat - h) #vector subtraction
weights = weights + alpha * dataMatrix.transpose()* error #matrix mult
print(weights)
return weights
计算出回归系数：
[[ 4.12414349]
[ 0.48007329]
[-0.6168482 ]]
画出数据集和Logistic最佳拟合直线的函数
def plotBestFit(weights):
import matplotlib.pyplot as plt
dataMat,labelMat=loadDataSet('testSet.txt')
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()
这个例子用到了100x3的数据集，数据量很小，但是却需要迭代300次，那么当样本达到千万甚至亿的级别的时候，遍历数据集的成本和代价会非常大。
下面介绍一种改进方法，随机梯度上升算法
该算法每次仅用一个样本点更新回归系数
def stocGradAscent0(dataMatrix, classLabels):
m,n = shape(dataMatrix)
alpha = 0.01
weights = 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
可以看到，随机梯度上升算法和梯度上升算法差别并不大
将每一行和回归系数相乘， weights是一个1x3的Numpy数组，初始化为1
[1. 1. 1.]
然后和每一行的特征值相乘，得出一个numpy型的float
<class 'numpy.float64'>
最后绘制图像
今天先到这儿吧，后续看情况更新，大家周末愉快~