python实现 logistic 回归 二分类算法 （通俗讲解逻辑回归本质与由来）

logistic回归

• 将数据样本看作是欧式空间的点，尝试找到一个超平面，将空间分成两部分，如果样本点在”正面“，则它被分为0类；如果样本点在”负面“，则它被分为1类。
• 怎么判断样本点在超平面的哪一面？将样本点坐标$x$代入超平面方程$a^{T}x=0$的等式左边$a^{T}x$，如果大于0，则在”阳面“；小于0，则在”阴面“；等于0，则在超平面上。
• 根据上面大于0 or 小于0 已经能判断属于哪一类了，再代入符号函数(机器学习里叫阶跃函数，数学人觉得就是个符号函数)，直接输出类别0 or 1，就更好了。这样就有了一个完整模型，即$\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(a^{T}x)$。
• 为什么非要套一个符号函数呢？我直接用大于0，小于0判断不行嘛？不行，因为后面要进行优化选参数$a$，有完整模型才能建立评价指标，比如风险函数，然后这个风险函数作为优化目标，取其极值点不是吗？
• $\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(a^{T}x)$这个模型并不是很好，因为后续要进行优化，绝大多数优化算法依靠的是导数信息(虽然现在有很多非光滑优化算法)，参数$a$作为变量，函数$\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(a^{T}x)$最好要是光滑的，但是它不是，因为原点是符号函数的跳跃间断点，它甚至都不连续啊…性质太差了！
• 怎么解决上面的问题？光滑化呗，取一个光滑函数去近似符号函数，下图是符号函数图像：
  
  我想让它变成下面这样光滑去近似符号函数：
  
• 正好，sigmoid函数完美符合要求，连续光滑，长得像符号函数。
• 把上面模型的符号函数替换成sigmoid函数就是逻辑回归模型了。

sigmoid函数：
$$\sigma (x)=\frac{1}{1+e^{-x}}$$
顺便记录一下它求导过程：
$${\sigma }^{{}^{\prime }}(x)=\frac{e^{-x}}{(1+e^{-x}{)}^2}=\frac{1}{1+e^{-x}}-\frac{1}{(1+e^{-x}{)}^2}={\sigma }^{{}^{\prime }}(x)-({\sigma }^{{}^{\prime }}(x){)}^2={\sigma }^{{}^{\prime }}(x)(1-{\sigma }^{{}^{\prime }}(x))$$

• 注： 以上很多名词只是为了自己方便记忆随便乱取的。
• 由于sigmoid函数取值不是离散的0，1，而是一个连续区间(0,1)，所以，可以看作是一种概率估计，概率大于0.5分入1类，小于0.5分入0类。

梯度下降法

太熟了，不写。
机器学习里的优化好像并没有给出明确的停机准则，都是靠最大迭代次数来停止的。

训练/优化：使用梯度下降找到模型的最佳参数

有关评价准则，风险函数，损失函数等有时间再补上，下面先只给出训练程序。

• 下面的梯度grad是对风险函数求导得到，有机会再补充，这本书也没有详细说明。
• 注意：dataMtrix第一列是1，这是对应于线性函数的截距项。

import numpy as np


## 逻辑回归
# 自定义sigmoid函数
## 逻辑回归
# 自定义sigmoid函数
def sigmoid(inX):
    res = np.zeros(inX.shape)
    for i in range(inX.shape[0]):
        if inX[i] >= 0:
            res[i] = 1 / (1 + np.exp(-inX[i]))
        else:
            res[i] = np.exp(inX[i]) / (1 + np.exp(inX[i]))
    return res


# 获取数据
def loadDataSet():
    """这里数据集使用python list存储"""
    dataMat = []
    labelMat = []
    # 一行一行读取文本，获取数据，读到空串的时候就是EOF
    f_obj = open('testSet.txt')
    lines = f_obj.readlines()
    for line in lines:
        line = line.strip()  # 去掉末尾换行符
        line_splited_list = line.split()
        labelMat.append(int(line_splited_list[-1]))
        dataMat.append([1.0, float(line_splited_list[0]), float(line_splited_list[1])])  # 注意不是用extend，dataMat每个元素就是样本列表
    return dataMat, labelMat


# 使用梯度下降法 训练数据 得到最佳参数
def gradDscent(dataMatIn, labelMatIn):
    # 将数据转换成numpy ndarray
    dataMatrix = np.mat(dataMatIn)
    labelMat = np.mat(labelMatIn).transpose()  # 转置
    m, n = dataMatrix.shape
    alpha = 0.001  # 迭代步长/学习率
    maxItNum = 500  # 最大迭代次数
    weights = np.ones((n, 1))  # 初始化参数为全1，是列向量
    for i in range(maxItNum):
        grad = -((labelMat - sigmoid(dataMatrix * weights)).transpose() * dataMatrix).transpose()
        weights = weights - alpha * grad
    return weights


w = gradDscent(*loadDataSet())
print(w)

结果

[[ 4.12414349]
 [ 0.48007329]
 [-0.6168482 ]]

画出决策边界

def plotBestFit(weights, dataSet, labels):
    # 先划分dataSet中属于不同类的点
    weights = weights.A
    dataMat = np.mat(dataSet)
    lablesArr = np.array(labels)
    aux_1 = lablesArr == 1
    aux_2 = ~aux_1
    class_1_cord1 = (dataMat[aux_1, :][:, 1]).A;
    class_1_cord2 = dataMat[aux_1, :][:, 2].A;
    class_2_cord1 = dataMat[aux_2, :][:, 1].A;
    class_2_cord2 = dataMat[aux_2, :][:, 2].A;
    fig = plt.figure(1)  # 画布对象
    ax = fig.add_subplot(111)  # 坐标对象
    # 画散点图
    ax.scatter(class_1_cord1, class_1_cord2, s=20, c='b', marker='.')
    ax.scatter(class_2_cord1, class_2_cord2, s=20, c='r', marker='v')
    # 画直线
    x = np.arange(-3, 3, 0.1)
    y = -(weights[1]*x + weights[0])/weights[2]
    plt.plot(x,y,'-g')
    plt.show()


plotBestFit(w, *loadDataSet())

报错记录

1. sigmoid 函数出现：python overflow encountered in exp。

这是由于如果x的值是负数，且绝对值非常大，那么exp(-x)会非常大，就会出现溢出的现象。
解决办法： 当x<0时，对sigmoid函数分子分母同乘以exp(-x)，函数值不变。
注意： 上面sigmoid函数是可以处理向量和矩阵的，就是对向量的每一个分量求函数值，所以判断的时候必须一个值一个值判断。

2. plot只接受一维的x,y。所以传入列向量(是matrix)会报错。可以用mat.A 或者mat.getA()将mat装换成ndarray。
3. numpy中 numpy.ndarray和numpy.matrix是不一样的！
4. suqueeze函数或者方法，可以去掉ndarray中维数是1的那个多余的维度。对matrix不管用！！
5. matrix must be 2-dimensional numpy中matrix类型数据只能是2维的。
6. 一个数组的元素全是逻辑值，对所有元素取否定不能前面加 -，而是加 ~。