机器学习-逻辑回归_python 负二项回归

机器学习-逻辑回归

本节简单介绍逻辑回归（Logistic Regression）的原理以及Python实现，包含以下内容，

1. 概念
2. Sigmoid函数
3. 最大似然估计
4. 梯度上升算法
5. 梯度上升算法求解J(θ)的最大θ值
6. Python实现
7. 随机梯度上升算法
8. 预测

部分内容引用自《Machine Learning in Action》

概念

逻辑回归（Logistic Regression）是一种用于解决二分类问题的机器学习算法，用于估计某种事物的可能性，即存在/不存在，是/不是等问题。比如某病人患有某种疾病的可能性，某个员工是否会离职，明天是否下雨等等。 注意，这里用的是“可能性”不是数学上的“概率”，逻辑回归的结果并非数学定义中的概率值，不可以直接当做概率值来用。

逻辑回归假设数据服从伯努利分布（0-1分布），其通过最大似然估计的方法，运用梯度上升算法来求解参数，最终达到数据二分类的目的。

逻辑回归的两个假设：

1）假设数据服从伯努利分布
2）假设模型的输出值是样本为正的概率。

逻辑回归与多重线性回归实际上有很多相同之处，最大的区别就在于它们的因变量不同，其他的基本都差不多。正是因为如此，这两种回归可以归于同一个家族，即广义线性模型（generalizedlinear model）。

这一家族中的模型形式基本上都差不多，不同的就是因变量不同。

• 如果是连续的，就是多重线性回归
• 如果是二项分布，就是Logistic回归
• 如果是Poisson分布，就是Poisson回归
• 如果是负二项分布，就是负二项回归

逻辑回归的主要用途：

• 寻找危险因素：寻找某一疾病的危险因素等
• 预测：根据模型，预测在不同的自变量情况下，发生某病或某种情况的可能性
• 判别：实际上跟预测有些类似，也是根据模型，判断某人属于某病或属于某种情况的概率有多大，也就是看一下这个人有多大的可能性是属于某病

逻辑回归主要在流行病学中应用较多，比较常用的情形是探索某疾病的危险因素，根据危险因素预测某疾病发生的概率，等等。例如，想探讨胃癌发生的危险因素，可以选择两组人群，一组是胃癌组，一组是非胃癌组，两组人群肯定有不同的体征和生活方式等。这里的因变量就是是否胃癌，即“是”或“否”，自变量就可以包括很多了，例如年龄、性别、饮食习惯、幽门螺杆菌感染等。自变量既可以是连续的，也可以是分类的。

Sigmoid函数

Sigmoid函数又叫激活函数，其函数表达式如下：

函数图像：

可知其定义域为全体实数，值域为（0, 1），Sigmoid函数的作用是将任何实数映射到（0, 1）。

逻辑回归的假设函数如下：

即，

其中，$\large x$是输入，即特征值，$\large \theta$是待求的回归系数。（使用此激活函数来拟合样本空间中的样本点，即，样本空间中的样本点都有一个真实的类别，通过寻找一个激活函数，该激活函数能够很好的计算（拟合）样本空间中的样本值。）

因为Sigmoid函数的值域是(0,1]，所以其在二分类问题上很适用。

逻辑回归要求我们根据一系列的学习样本（即上面公式中的$\large x$，多维矩阵）得到一个$\large \theta$向量，然后来预测新的元素分类。

最大似然估计

逻辑回归采用最大似然估计来作为其损失函数（这部分参考自：Logistic回归基础篇之梯度上升算法）：

注：如果采用常见的误差平方和来当做损失函数，即如下函数，

这是一个非凸函数，这就意味着损失函数有着许多的局部最小值，不利于求解（参考：逻辑回归与最大似然估计推导_逻辑回归极大似然估计推导-CSDN博客）。所以，逻辑回归中使用的是最大最大似然估计来作为其损失函数。

根据sigmoid函数的特性，我们可以做出如下的假设：

上式即为在已知样本x和参数θ的情况下，样本x属性正样本(y=1)和负样本(y=0)的条件概率。理想状态下，根据上述公式，求出各个点的概率均为1，也就是完全分类都正确。但是考虑到实际情况，样本点的概率越接近于1，其分类效果越好。比如一个样本属于正样本的概率为0.51，那么我们就可以说明这个样本属于正样本。另一个样本属于正样本的概率为0.99，那么我们也可以说明这个样本属于正样本。但是显然，第二个样本概率更高，更具说服力。我们可以把上述两个概率公式合二为一：

合并出来的Cost，我们称之为代价函数(Cost Function)。当y等于1时，(1-y)项(第二项)为0；当y等于0时，y项(第一项)为0。为了简化问题，我们对整个表达式求对数，(将指数问题对数化是处理数学问题常见的方法)：

这个代价函数，是对于一个样本而言的。给定一个样本，我们就可以通过这个代价函数求出，样本所属类别的概率，而这个概率越大越好，所以也就是求解这个代价函数的最大值。既然概率出来了，那么最大似然估计也该出场了。假定样本与样本之间相互独立，那么整个样本集生成的概率即为所有样本生成概率的乘积，再将公式对数化，便可得到如下代价函数公式：

其中，m为样本的总数，y(i)表示第i个样本的类别，x(i)表示第i个样本，需要注意的是θ是多维向量，x(i)也是多维向量。

综上所述，满足J(θ)的最大的θ值即是我们需要求解的模型。（求解θ，使得该θ可以让样本空间中的样本点出现其真实类别的概率最大，后验概率。）

怎么求解使J(θ)最大的θ值呢？因为是求最大值，所以我们需要使用梯度上升算法。如果面对的问题是求解使J(θ)最小的θ值，那么我们就需要使用梯度下降算法。面对我们这个问题，如果使J(θ) := -J(θ)，那么问题就从求极大值转换成求极小值了，使用的算法就从梯度上升算法变成了梯度下降算法，它们的思想都是相同的，学会其一，就也会了另一个。本文使用梯度上升算法进行求解。

梯度上升算法

梯度上升算法用于求解函数的最大值，梯度下降算法用于求解函数的最小值，其原理都一样。因为逻辑回归中求解的是某个元素属于某个分类的最大可能性，求解的是最大值，因此需要使用梯度上升算法。公式：

​

其中，$\large \alpha$是步长，$\large \frac{\partial f(X_{i})) }{X_{i}}$是梯度，也就是偏导数组成的向量，表示$\large X$应该沿此方向运动方能取到最大值。下面举一个实例，

求以下函数的最大值：

$\large f(x) = -x^{2} + 2x + 2$

函数图像：

​

可知，当$\large x=1$时取最大值为$\large y=3$。

下面通过Python实现梯度上升算法，

import math
 
 
def f(x):
    return -x * x + 2 * x + 2
 
 
def ff(x):
    return -2 * x + 2
 
 
def gradient_ascent(x_0=1, steps=0.1):
    while True:
        x_1 = x_0 + ff(x_0) * steps
        if math.fabs(x_0 - x_1) <= 0.000001:
            return x_1
        x_0 = x_1
 
 
x = gradient_ascent(-1)
print("Max value x: %r" % x)
print("Max value y: %r" % f(x))

 运行结果：

D:\work\python_workspace\machine_learning\venv\Scripts\python.exe D:/work/python_workspace/machine_learning/logistic_regression/gradient_ascent.py
Max value x: 0.9999961687611478
Max value y: 2.9999999999853215
 
Process finished with exit code 0

使用梯度上升算法求解J(θ)的最大θ值

由上小节可知J(θ)为：

sigmoid函数为：

那么，现在我只要求出J(θ)的偏导，就可以利用梯度上升算法，求解J(θ)的极大值了。

那么现在开始求解J(θ)对θ的偏导，求解如下(数学推导)：

其中，是步长系数，可以是0.1,0.01,0.001等；是原始样本的类别值，可以是0,1等，通过原始数据集直接得到；表示通过Sigmoid函数在当前场景下计算出来的估计值；$x_{j}^{(i)}$表示原始样本的特征值，通过原始数据集直接得到。

知道了，梯度上升迭代公式，我们就可以自己编写代码，计算最佳拟合参数了。

Python实现

数据集文件：test_set.txt，

-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
-1.781871	9.097953	0
-0.566606	5.749003	1
0.931635	1.589505	1
-0.024205	6.151823	1
-0.036453	2.690988	1
-0.196949	0.444165	1
1.014459	5.754399	1
1.985298	3.230619	1
-1.693453	-0.557540	1
-0.576525	11.778922	0
-0.346811	-1.678730	1
-2.124484	2.672471	1
1.217916	9.597015	0
-0.733928	9.098687	0
-3.642001	-1.618087	1
0.315985	3.523953	1
1.416614	9.619232	0
-0.386323	3.989286	1
0.556921	8.294984	1
1.224863	11.587360	0
-1.347803	-2.406051	1
1.196604	4.951851	1
0.275221	9.543647	0
0.470575	9.332488	0
-1.889567	9.542662	0
-1.527893	12.150579	0
-1.185247	11.309318	0
-0.445678	3.297303	1
1.042222	6.105155	1
-0.618787	10.320986	0
1.152083	0.548467	1
0.828534	2.676045	1
-1.237728	10.549033	0
-0.683565	-2.166125	1
0.229456	5.921938	1
-0.959885	11.555336	0
0.492911	10.993324	0
0.184992	8.721488	0
-0.355715	10.325976	0
-0.397822	8.058397	0
0.824839	13.730343	0
1.507278	5.027866	1
0.099671	6.835839	1
-0.344008	10.717485	0
1.785928	7.718645	1
-0.918801	11.560217	0
-0.364009	4.747300	1
-0.841722	4.119083	1
0.490426	1.960539	1
-0.007194	9.075792	0
0.356107	12.447863	0
0.342578	12.281162	0
-0.810823	-1.466018	1
2.530777	6.476801	1
1.296683	11.607559	0
0.475487	12.040035	0
-0.783277	11.009725	0
0.074798	11.023650	0
-1.337472	0.468339	1
-0.102781	13.763651	0
-0.147324	2.874846	1
0.518389	9.887035	0
1.015399	7.571882	0
-1.658086	-0.027255	1
1.319944	2.171228	1
2.056216	5.019981	1
-0.851633	4.375691	1
-1.510047	6.061992	0
-1.076637	-3.181888	1
1.821096	10.283990	0
3.010150	8.401766	1
-1.099458	1.688274	1
-0.834872	-1.733869	1
-0.846637	3.849075	1
1.400102	12.628781	0
1.752842	5.468166	1
0.078557	0.059736	1
0.089392	-0.715300	1
1.825662	12.693808	0
0.197445	9.744638	0
0.126117	0.922311	1
-0.679797	1.220530	1
0.677983	2.556666	1
0.761349	10.693862	0
-2.168791	0.143632	1
1.388610	9.341997	0
0.317029	14.739025	0

创建模块 log_regres.py，并输入以下代码：

import numpy as np
 
 
def load_data():
    data_mat = []
    label_mat = []
    with open('test_set.txt') as f:
        for line in f.readlines():
            line_array = line.strip().split()
            data_mat.append([1.0, float(line_array[0]), float(line_array[1])])
            label_mat.append(int(line_array[2]))
    return data_mat, label_mat
 
 
def sigmoid(in_X):
    return 1.0 / (1 + np.exp(-in_X))
 
 
def grad_ascent(data_mat_in, class_labels):
    data_matrix = np.mat(data_mat_in)
    label_mat = np.mat(class_labels).transpose()
    m, n = np.shape(data_matrix)
    alpha = 0.001
    max_cycles = 500
    weights = np.ones((n, 1))
    for k in range(max_cycles):
        h = sigmoid(data_matrix * weights)
        error = label_mat - h
        weights = weights + alpha * data_matrix.transpose() * error
    return weights
 
 
if __name__ == '__main__':
    data, label = load_data()
    weight = grad_ascent(data, label)
    print(weight)

注意下面几行代码，就是通过梯度上升算法得到的公式的具体实现，

h = sigmoid(data_matrix * weights)
error = label_mat - h
weights = weights + alpha * data_matrix.transpose() * error

输出结果：

D:\work\python_workspace\machine_learning\venv\Scripts\python.exe D:/work/python_workspace/machine_learning/logistic_regression/log_regres.py
[[ 4.12414349]
 [ 0.48007329]
 [-0.6168482 ]]
 
Process finished with exit code 0

因此，该逻辑回归的决策函数为：

$\large z = 4.12414349 + 0.48007329 * x_{1} -0.6168482 * x_{2}$

下面我们根据该决策函数画出决策边界，创建模块 plot_best_fit.py，并输入以下代码：

import matplotlib.pyplot as plt
import numpy as np
import logistic_regression.log_regres as lr
 
def plot_best_fit(weights):
    data, label = lr.load_data()
    data_array = np.array(data)
    n = np.shape(data_array)[0]
    xcord1 = []
    ycord1 = []
    xcord2 = []
    ycord2 = []
    for i in range(n):
        if int(label[i]) == 1:
            xcord1.append(data_array[i, 1])
            ycord1.append(data_array[i, 2])
        else:
            xcord2.append(data_array[i, 1])
            ycord2.append(data_array[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 = np.arange(-3.0, 3.0, 0.1)
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x, np.array(y)[0])
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()
 
 
if __name__ == '__main__':
    data, label = lr.load_data()
    weights = lr.grad_ascent(data, label)
    plot_best_fit(weights)

运行结果：

​

注意，假设以0.5为决策阈值，通过Sigmoid可知$z$的值为0，可得到$\large x_{1}$和$\large x_{2}$的以下关系式，

$\large x_{2} = \frac{4.12414349 + 0.48007329 * x_{1}}{0.6168482}$

这就是该示例的决策边界函数。

随机梯度上升算法

上面的梯度算法需要进行大量的计算，在样本点和特征值很多的情况下会有一定的性能问题，我们可以通过随机梯度算法进行优化。随机梯度算法和一般梯度算法类似，只是每次计算时是随机的从样本中选取一个样本点，而不是像一般梯度算法中每次迭代时把所有样本点都计算一次。且随机梯度算法每次迭代时都要逐渐减小步长。

创建模块 stoc_grad_ascent1.py，并输入以下代码：

import numpy as np
import random
import matplotlib.pyplot as plt
import logistic_regression.log_regres as lr
import logistic_regression.plot_best_fit as pf
 
 
def stoc_grad_ascent1(data_matrix, class_label, num_iter=150):
    m, n = np.shape(data_matrix)
    weights = np.ones(n)
    for j in range(num_iter):
        data_index = list(range(m))
        for i in range(m):
            alpha = 4 / (1.0 + j + i) + 0.01
            rand_index = int(random.uniform(0, len(data_index)))
            h = lr.sigmoid(np.sum(data_matrix[rand_index] * weights))
            error = class_label[rand_index] - h
            weights = weights + alpha * error * data_matrix[rand_index]
            del (data_index[rand_index])
    return weights
 
 
def plot_best_fit(weights):
    data, label = lr.load_data()
    data_array = np.array(data)
    n = np.shape(data_array)[0]
    xcord1 = []
    ycord1 = []
    xcord2 = []
    ycord2 = []
    for i in range(n):
        if int(label[i]) == 1:
            xcord1.append(data_array[i, 1])
            ycord1.append(data_array[i, 2])
        else:
            xcord2.append(data_array[i, 1])
            ycord2.append(data_array[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 = np.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()
 
 
if __name__ == '__main__':
    data_arr, label_mat = lr.load_data()
    weights = stoc_grad_ascent1(np.array(data_arr), label_mat)
    print(weights)
    plot_best_fit(weights)

运行结果：

D:\work\python_workspace\machine_learning\venv\Scripts\python.exe D:/work/python_workspace/machine_learning/logistic_regression/stoc_grad_ascent1.py
[13.66512034  0.9332941  -2.05701081]
 
Process finished with exit code 0

图像：

​

预测

一旦得到θ向量，我们就可以得到逻辑回归的假设函数的准确表达式，

将目标变量的特征值带入到此表达式，得到$h_{\theta }(x)$的值，再通过设定阈值（比如0.5）的判断，即可得到目标变量所属的类别。