参考1
参考2

0.前言

Logistic regression （逻辑回归）是当前业界比较常用的机器学习方法，用于估计某种事物的可能性。之前在经典之作《数学之美》中也看到了它用于广告预测，也就是根据某广告被用户点击的可能性，把最可能被用户点击的广告摆在用户能看到的地方，然后叫他“你点我啊！”用户点了，你就有钱收了。这就是为什么我们的电脑现在广告泛滥的原因了。
还有类似的某用户购买某商品的可能性，某病人患有某种疾病的可能性啊等等。这个世界是随机的（当然了，人为的确定性系统除外，但也有可能有噪声或产生错误的结果，只是这个错误发生的可能性太小了，小到千万年不遇，小到忽略不计而已），所以万物的发生都可以用可能性或者几率（Odds）来表达。“几率”指的是某事物发生的可能性与不发生的可能性的比值。
Logistic regression可以用来回归，也可以用来分类，主要是二分类。如向量机SVM，它就是个二分类的，例如，它可以将两个不同类别的样本给分开，思想是找到最能区分它们的那个分类超平面。但当你给一个新的样本给它，它能够给你的只有一个答案，你这个样本是正类还是负类。例如你问SVM，某个女生是否喜欢你，它只会回答你喜欢或者不喜欢。这对我们来说，显得太粗鲁了，要不希望，要不绝望，这都不利于身心健康。那如果它可以告诉我，她很喜欢、有一点喜欢、不怎么喜欢或者一点都不喜欢，你想都不用想了等等，告诉你她有49%的几率喜欢你，总比直接说她不喜欢你，来得温柔。而且还提供了额外的信息，她来到你的身边你有多少希望，你得再努力多少倍，知己知彼百战百胜，哈哈。Logistic regression就是这么温柔的，它给我们提供的就是你的这个样本属于正类的可能性是多少。
假设我们的样本是{x, y}，y是0或者1，表示正类或者负类，x是我们的m维的样本特征向量。那么这个样本x属于正类，也就是y=1的“概率”可以通过下面的逻辑函数来表示：

1.what is the logistic regression?

许多人对线性回归都比较熟悉，但知道逻辑回归的人可能就要少的多。从大的类别上来说，逻辑回归是一种有监督的统计学习方法，主要用于对样本进行分类。

在线性回归模型中，输出一般是连续的，例如

逻辑回归与线性回归的关系

逻辑回归也被称为广义线性回归模型，它与线性回归模型的形式基本上相同，都具有 ax+b，其中a和b是待求参数，其区别在于他们的因变量不同，多重线性回归直接将ax+b作为因变量，即y = ax+b，而logistic回归则通过函数S将ax+b对应到一个隐状态p，p = S(ax+b)，然后根据p与1-p的大小决定因变量的值。这里的函数S就是Sigmoid函数。

sigmoid函数的图像

通过函数S的作用，我们可以将输出的值限制在区间[0， 1]上，p(x)则可以用来表示概率p(y=1|x)，即当一个x发生时，y被分到1那一组的概率。可是，等等，我们上面说y只有两种取值，但是这里却出现了一个区间[0, 1]，这是什么鬼？？其实在真实情况下，我们最终得到的y的值是在[0, 1]这个区间上的一个数，然后我们可以选择一个阈值，通常是0.5，当y>0.5时，就将这个x归到1这一类，如果y<0.5就将x归到0这一类。但是阈值是可以调整的，比如说一个比较保守的人，可能将阈值设为0.9，也就是说有超过90%的把握，才相信这个x属于1这一类。了解一个算法，最好的办法就是自己从头实现一次。下面是逻辑回归的具体实现。

逻辑回归模型的代价函数

逻辑回归一般使用交叉熵作为代价函数。关于代价函数的具体细节，请参考代价函数
这里只给出交叉熵公式

hθ(x)：用参数θ和x预测出来的y值；

y：原训练样本中的y值，也就是标准答案

上角标(i)：第i个样本

优化求解(以0.前言中最前面的最大似然函数举例)

梯度下降(gradient descent)

Gradient descent 又叫 steepest descent，是利用一阶的梯度信息找到函数局部最优解的一种方法，也是机器学习里面最简单最常用的一种优化方法。它的思想很简单，要找最小值，我只需要每一步都往下走（也就是每一步都可以让代价函数小一点），然后不断的走，那肯定能走到最小值的地方，例如下图所示

但，我同时也需要更快的到达最小值啊，怎么办呢？我们需要每一步都找下坡最快的地方，也就是每一步我走某个方向，都比走其他方法，要离最小值更近。而这个下坡最快的方向，就是梯度的负方向了。
对logistic Regression来说，梯度下降算法新鲜出炉，如下：

其中，参数α叫学习率，就是每一步走多远，这个参数非常关键的。如果设置的太多，那么很容易就在最优值附加徘徊，因为你步伐太大了。例如要从广州到上海，但是你的一步的距离就是广州到北京那么远，没有半步的说法，自己能迈那么大步，是幸运呢？还是不幸呢？事物总有两面性嘛，它带来的好处是能很快的从远离最优值的地方回到最优值附近，只是在最优值附近的时候，它有心无力了。但如果设置的太小，那收敛速度就太慢了，向蜗牛一样，虽然会落在最优的点，但是这速度如果是猴年马月，我们也没这耐心啊。所以有的改进就是在这个学习率这个地方下刀子的。我开始迭代是，学习率大，慢慢的接近最优值的时候，我的学习率变小就可以了。所谓采两者之精华啊！
梯度下降算法的伪代码如下：

初始化回归系数为1
重复下面步骤直到收敛{
        计算整个数据集的梯度
        使用alpha x gradient来更新回归系数
}
返回回归系数值

  注：因为本文中是求解的Logit回归的代价函数是似然函数，需要最大化似然函数。所以我们要用的是梯度上升算法。但因为其和梯度下降的原理是一样的，只是一个是找最大值，一个是找最小值。找最大值的方向就是梯度的方向，最小值的方向就是梯度的负方向。另外，最大似然可以通过取负对数，转化为求最小值。

随机梯度下降SGD (stochastic gradient descent)

梯度下降算法在每次更新回归系数的时候都需要遍历整个数据集（计算整个数据集的回归误差），该方法对小数据集尚可。但当遇到有数十亿样本和成千上万的特征时，就有点力不从心了，它的计算复杂度太高。改进的方法是一次仅用一个样本点（的回归误差）来更新回归系数。这个方法叫随机梯度下降算法。由于可以在新的样本到来的时候对分类器进行增量的更新（假设我们已经在数据库A上训练好一个分类器h了，那新来一个样本x。对非增量学习算法来说，我们需要把x和数据库A混在一起，组成新的数据库B，再重新训练新的分类器。但对增量学习算法，我们只需要用新样本x来更新已有分类器h的参数即可），所以它属于在线学习算法。与在线学习相对应，一次处理整个数据集的叫“批处理”。
随机梯度下降算法的伪代码如下：

初始化回归系数为1

重复下面步骤直到收敛{

        对数据集中每个样本

               计算该样本的梯度

                使用alpha xgradient来更新回归系数

 }

返回回归系数值

改进的随机梯度下降

评价一个优化算法的优劣主要是看它是否收敛，也就是说参数是否达到稳定值，是否还会不断的变化？收敛速度是否快？

上图展示了随机梯度下降算法在200次迭代中（请先看第三和第四节再回来看这里。我们的数据库有100个二维样本，每个样本都对系数调整一次，所以共有200*100=20000次调整）三个回归系数的变化过程。其中系数X2经过50次迭代就达到了稳定值。但系数X1和X0到100次迭代后稳定。而且可恨的是系数X1和X2还在很调皮的周期波动，迭代次数很大了，心还停不下来。产生这个现象的原因是存在一些无法正确分类的样本点，也就是我们的数据集并非线性可分，但我们的logistic regression是线性分类模型，对非线性可分情况无能为力。然而我们的优化程序并没能意识到这些不正常的样本点，还一视同仁的对待，调整系数去减少对这些样本的分类误差，从而导致了在每次迭代时引发系数的剧烈改变。对我们来说，我们期待算法能避免来回波动，从而快速稳定和收敛到某个值。
对随机梯度下降算法，我们做两处改进来避免上述的波动问题：
1）在每次迭代时，调整更新步长alpha的值。随着迭代的进行，alpha越来越小，这会缓解系数的高频波动（也就是每次迭代系数改变得太大，跳的跨度太大）。当然了，为了避免alpha随着迭代不断减小到接近于0（这时候，系数几乎没有调整，那么迭代也没有意义了），我们约束alpha一定大于一个稍微大点的常数项。
2）每次迭代，改变样本的优化顺序。也就是随机选择样本来更新回归系数。这样做可以减少周期性的波动，因为样本顺序的改变，使得每次迭代不再形成周期性。
改进的随机梯度下降算法的伪代码如下：

初始化回归系数为1
重复下面步骤直到收敛{

       对随机遍历的数据集中的每个样本

              随着迭代的逐渐进行，减小alpha的值

              计算该样本的梯度

              使用alpha x gradient来更新回归系数

    }

返回回归系数值

比较原始的随机梯度下降和改进后的梯度下降，可以看到两点不同：

1）系数不再出现周期性波动。2）系数可以很快的稳定下来，也就是快速收敛。这里只迭代了20次就收敛了。而上面的随机梯度下降需要迭代200次才能稳定。

2.数据实验演示

-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

上面的数据一共是3列10行，其中前两列为x1和x2的值，第3列表示y的值；10行表示取了10个样本点。我们可以将这些数据当做训练模型参数的训练样本。

见到训练样本就可以比较直观的理解算法的输入，以及我们如何利用这些数据来训练逻辑回归分类器，进而用训练好的模型来预测新的样本（检测样本）。

从逻辑回归的参数形式，式子（1）我们可以看到逻辑回归模型中有两个待定参数a（x的系数）和b（常数项），我们现在给出来的数据有两个特征x1, x2，因此整个模型就增加了一项：ax1 + cx2 + b。为了形式上的统一，我们使用带下标的a表示不同的参数（a0表示常数项b并作x0的参数令 <x0=1>，a1、a2分别表示x1和x2的参数），就可以得到：

注：数据是二维的，也就是说这组观察样本中有两个自变量，即两个特征（feature）。

3.训练分类器

就像上面说的，训练分类器的过程，就是根据已经知道的数据（训练样本）确定一个使得代价函数的值最小的a（参数向量/回归系数）的过程。逻辑回归模型属于有监督的学习方法，上面示例数据中的第3列其实是训练样本提供的"标准答案"。也就是说，这些数据是已经分好类的（两类，0或者1）。在训练阶段，我们要做的就是利用训练样本和（2）式中的模型，估计一个比较合适的参数a，使得仅通过前面两列数据（观察值/测量值）就可以估计一个值h(a)，这个值越接近标准答案y，说明我们的模型预测的越准确。

import numpy as np
import os
import matplotlib.pyplot as plt
path = 'D:/Python basic/logistic regression/data'
training_sample = 'Logistic_Regression-trainingSample.txt'
testing_sample = 'Logistic_Regression-testingSample.txt'

# 从文件中读入训练样本的数据
def loadDataSet(p, file_n):
    dataMat = []
    labelMat = []
    fr = open(os.path.join(p, file_n))
    for line in fr.readlines():# readlines() 方法用于读取所有行并返回列表list 该列表可以由 Python 的 for... in ... 结构进行处理
        lineArr = line.strip().split() #s.strip(rm) 删除s字符串中开头、结尾处
        # 三个特征x0, x1, x2, x0=1
        # print(lineArr[0])
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))  # 样本标签y
    return dataMat, labelMat

def sigmoid(X):# sigmoid函数
    return 1.0/(1+np.exp(-X))

# 梯度下降法求回归系数a
def gradAscent(dataMatIn, classLabels):
    dataMatrix = np.mat(dataMatIn)             # 转换成numpy中的矩阵, X, 90 x 3
    labelMat = np.mat(classLabels).transpose()  # 转换成numpy中的矩阵, y, 90 x 1
    m, n = np.shape(dataMatrix)  # m=90, n=3
    alpha = 0.001  # 学习率
    maxCycles = 1000
    weights = np.ones((n, 1))  # 初始参数, 3 x 1
    for k in range(maxCycles):              # heavy on matrix operations
        h = sigmoid(np.dot(dataMatrix, weights))     # 模型预测值, 90 x 1, 矩阵乘法
        error = h - labelMat              # 真实值与预测值之间的误差, 90 x 1
        temp = np.dot(dataMatrix.transpose(), error)  # 所有参数的偏导数, 3 x 1, np.dot矩阵乘法
        weights = weights - alpha * temp  # 更新权重
    return weights

# 分类效果展示，参数weights就是回归系数
def plotBestFit(weights):
    dataMat,labelMat=loadDataSet(path, training_sample)
    dataArr = np.array(dataMat)
    n = np.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 = np.arange(-3.0, 3.0, 0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]  # x2 = f(x1)
    ax.plot(x.reshape(1, -1), y.reshape(1, -1))
    plt.xlabel('X1'); plt.ylabel('X2');
    plt.show()
# 测试函数
def test_logistic_regression():
    dataArr, labelMat = loadDataSet(path, training_sample)  # 读入训练样本中的原始数据
    A = gradAscent(dataArr, labelMat)  # 回归系数a的值
    h = sigmoid(np.dot(dataArr, A)) # 预测结果h(a)的值, 矩阵乘法
    # print(dataArr, labelMat)
    print("回归系数：",A)
    print("预测值",h)
    plotBestFit(A)

test_logistic_regression()