机器学习算法与Python实践(2) - 逻辑回归

首先，先说一下机器学习算法的一般步骤：

对于一个问题，首先用数学语言进行描述，再建立一个数学模型。例如，用回归/分类模型来描述该问题

通过最大似然、最大后验概率或者最小化分类误差等等建立模型的代价函数，也就是一个最优化的问题。找到最优化问题的解，也就是能拟合我们的数据的最好的模型参数

求解代价函数，找到最优解。求解分很多情况：

a. 优化函数存在解析解。例如：求最值。首先对代价函数求导，找到导数为0的点，就可以找到最(大/小)值。如果代价函数能简单求导，并且求导后为0的式子存在解析解，就可以得出最有的参数了。

b. 1)优化函数很难求导。例如：函数之间存在隐含的变量/变量间存在耦合。2)求导后的式子很难求解。例如：未知参数的个数大雨已知方程组的个数。 这种情况下，需要借助迭代一步步接近最优解。

如果代价函数是凸函数，那么就存在全局最优解。如果是非凸函数，就存在多个局部最优解。

逻辑回归的一般步骤：

拿到一个分类或者回归问题 —–> 建立代价函数 —–> 求导或者迭代求出最优的模型参数 —–> 验证模型的好坏。找到最合适的模型即可。

逻辑回归的优缺点

优点：计算代价不高，易于理解和实现。

缺点：容易欠拟合，分类精度可能不高。

使用数据类型：数值型和标称型数据。

逻辑回归（LogisticRegression）

逻辑回归是当前机器学习领域比较常用的方法，常用于估计某种事物的可能性。例如：某用户购买某件商品的可能性，某病人得某种疾病的可能性等。

逻辑回归可以用来回归，也可以用来分类（主要用于二分类）。类似于SVM分类器，就是一个二分类的例子（0/1）。逻辑回归提供的就是，这个实例属于正类的可能性有多大。

假设样本为{x,y}，y为 0/1，表示正类或者负类，x为我们的m维的特征向量。那么样本x属于正类（y=1）的“概率”可以通过下面的逻辑函数来表示：



θ为模型参数(回归系数)，σ是sigmoid函数。

实际上这个函数是由下面的对数几率(也就是x属于正类的可能性和负类的可能性的比值的对数)变换得到的：



其中y为分类结果，x1,x2,x3,x4 … xm为特征向量，θ1,θ2,θ3,θ4 … θm 为回归系数(特征向量的系数/权值)

所以说从上述的logistic回归可以看出其就是一个线性分类模型，而它与线性回归的不同点在于：逻辑回归将很大的范围的数压缩到0～1之间，更直观，简单的进行表达。对于二分类来说，可以设为：y >= m 为正类，y < m 为负类（m为0～1之间任意数值）实际上，SVM的类概率就是样本到边界的距离，就是通过LogisticsRegression实现的。



所以，LogisticsRegression 就是一个比logistic方程归一化后的线性回归。

模型选择之后，模型的参数θ还是未知的，所以我们需要用我们收集到的数据来进行训练，得到θ。那么下一步就是建立我们的代价函数了。

LogisticsRegression最基本的学习算法是最大似然。//最大似然

假设我们有n个独立的训练样本

。那么每一个观察到的样本 (xi，yi) 出现的概率是：



从上面的式子中可以看出。当y = 1 时，后面的一项为1，只剩下x = 1 的概率；同样，当y = 0 时，前面的一项为1，只剩下x = 0 的概率。所以不管y是0还是1，上面得到的数，都是(x, y)出现的概率。那么整个样本集，也就是n个独立的样本出现的似然函数为（因为每个样本都是独立的，所以n个样本出现的概率就是他们各自出现的概率相乘）：



那最大似然法就是求模型中使得似然函数最大的系数取值θ*。这个最大似然就是我们的代价函数(cost function)了。

代价函数有了之后，下一步要做的就是优化求解了。我们先尝试对上面的代价函数求导，令导数为0，看可不可以解出来(有没有解析解)，如果没有就通过迭代了。

首先变换L(θ)：取自然对数，然后化简。注：有xi的时候，表示它是第i个样本，下面没有做区分了，得到：



(其中第三步到第四部用了以下变换：)



用L(θ)对θ求导，得到：



然后令该导数为0，发现其无法解析求解。

最后问题变成了，求解参数使方程L最大化，求解参数的方法梯度上升法（原理这里不解释了，看详细的代码的计算方式应该更容易理解些）。

根据这个转换公式



我们代入参数和特征，求P，也就是发生1的概率。



上面这个也就是常提及的sigmoid函数，俗称激活函数，最后用于分类（若P(y=1|x;Θ )大于0.5，则判定为1）。

下面是详细的逻辑回归代码，代码比较简单，主要是要理解上面的算法思想。个人建议，可以结合代码看一步一步怎么算的，然后对比上面推导公式，可以让人更加容易理解，并加深印象。

from numpy import *

filename = '../txtfile/test.txt'  # 文件目录

def loadDataSet():  # 读取数据（这里只有两个特征）
dataMat = []
labelMat = []
fr = open(filename)
for line in fr.readlines():
lineArr = line.strip().split()
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])  # 前面的1，表示方程的常量。比如两个特征X1,X2，共需要三个参数，W1+W2*X1+W3*X2
labelMat.append(int(lineArr[2]))
return dataMat, labelMat

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

def gradAscent(dataMat, labelMat):  # 梯度上升求最优参数
dataMatrix = mat(dataMat)  # 将读取的数据转换为矩阵
classLabels = mat(labelMat).transpose()  # 将读取的数据转换为矩阵
m, n = shape(dataMatrix)
alpha = 0.001  # 设置梯度的阀值，该值越大梯度上升幅度越大
maxCycles = 500  # 设置迭代的次数，一般看实际数据进行设定，有些可能200次就够了
weights = ones((n, 1))  # 设置初始的参数，并都赋默认值为1。注意这里权重以矩阵形式表示三个参数。
for k in range(maxCycles):
h = sigmoid(dataMatrix * weights)
error = (classLabels - h)  # 求导后差值
weights = weights + alpha * dataMatrix.transpose() * error  # 迭代更新权重
return weights

def stocGradAscent0(dataMat, labelMat):  # 随机梯度上升，当数据量比较大时，每次迭代都选择全量数据进行计算，计算量会非常大。所以采用每次迭代中一次只选择其中的一行数据进行更新权重。
dataMatrix = mat(dataMat)
classLabels = labelMat
m, n = shape(dataMatrix)
alpha = 0.01
maxCycles = 500
weights = ones((n, 1))
for k in range(maxCycles):
for i in range(m):  # 遍历计算每一行
h = sigmoid(sum(dataMatrix[i] * weights))
error = classLabels[i] - h
weights = weights + alpha * error * dataMatrix[i].transpose()
return weights

def stocGradAscent1(dataMat, labelMat):  # 改进版随机梯度上升，在每次迭代中随机选择样本来更新权重，并且随迭代次数增加，权重变化越小。
dataMatrix = mat(dataMat)
classLabels = labelMat
m, n = shape(dataMatrix)
weights = ones((n, 1))
maxCycles = 500
for j in range(maxCycles):  # 迭代
dataIndex = [i for i in range(m)]
for i in range(m):  # 随机遍历每一行
alpha = 4 / (1 + j + i) + 0.0001  # 随迭代次数增加，权重变化越小。
randIndex = int(random.uniform(0, len(dataIndex)))  # 随机抽样
h = sigmoid(sum(dataMatrix[randIndex] * weights))
error = classLabels[randIndex] - h
weights = weights + alpha * error * dataMatrix[randIndex].transpose()
del (dataIndex[randIndex])  # 去除已经抽取的样本
return weights

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()

def main():
dataMat, labelMat = loadDataSet()
weights = gradAscent(dataMat, labelMat).getA()
plotBestFit(weights)

if __name__ == '__main__':
main()

结果展示：



当然，还可以换随机梯度上升和改进的随机梯度上升算法试试，效果可能会更好一点。

test.txt数据下载