1.Logistic回归

Logistic回归是实践中解决分类问题的最重要方法，方法简单、容易实现、效果良好、易于解释，应用范围不只是分类，亦可用于推荐系统。

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

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

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

——《机器学习实战》 $P_{74}$

Logistic回归的一般过程：

收集数据：采用任意方法收集数据。

准备数据：由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳。

分析数据：采用任意方法对数据进行分析。

训练算法：大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数。

测试算法：一旦训练步驟完成，分类将会很快。

使用算法：首先，我们需要输入一些数据，并将其转换成对应的结构化数值；接着，基于训练好的回归系数就可以对这些数值进行简单的回归计算，判定它们属于哪个类别，在这之后，我们就可以夺输出的类别上做一些其他分析工作。

——《机器学习实战》 $P_{73}$

对于二分类问题，其输出为0或1，最理想的模型是”单位阶跃函数”（unit-step function），如下图右侧红字所示。而对数几率函数（Logistic function）正是可替代”单位阶跃函数”并连续的函数，如下式；该函数是一种“Sigmoid函数”（即形似S的函数）。

y = \frac{1}{1 + e^{- z}}

$y=\frac{1}{1+e^{-z}}$
这里写图片描述

为了实现Logistic回归分类器，对每个特征都乘以各自的回归系数，结果相加，再带入到对数几率函数的z值中，如下所示，也可写成向量形式 $z = ω^TX$ 。为了算出最佳系数，需要用到最优化理论的知识。

z = ω_{0} x_{0} + ω_{1} x_{1} + ω_{2} x_{2} + \dots \dots + ω_{0} n x_{n}

$z=ω_0x_0+ω_1x_1+ω_2x_2+……+ω_0nx_n$

2.训练常用算法实现

2.1梯度上升

梯度上升法的思想是，沿着函数的梯度方向探寻，能找到函数的最大值。梯度上升算法的迭代公式如下，α为步长：

ω := ω + α ▽_{ω} f (ω)

$\omega:=\omega+α▽_ωf(ω)$
而梯度下降算法与梯度上升算法相同，只是梯度下降算法用来求函数的最小值，而梯度上升算法用来求最大值，迭代公式如下：

ω := ω - α ▽_{ω} f (ω)

$\omega:=\omega-α▽_ωf(ω)$
代码实现如下：

def sigmoid(inX):
    '''
    sigmoid函数
    :param inX: 输入值
    :return: sigmoid结果
    '''
    return 1.0/(1+exp(-inX))

def gradAscent(dataMatIn, classLabels):
    '''
    梯度上升算法
    :param dataMatIn: 数据集
    :param classLabels: 数据标签
    :return: 回归系数
    '''
    dataMatrix = mat(dataMatIn)             # 转换为NumPy matrix
    labelMat = mat(classLabels).transpose() # 转换为NumPy matrix，转置为列向量
    m,n = shape(dataMatrix)                 # 返回dataMatrix行数和列数
    alpha = 0.001                           # 移动步长
    maxCycles = 500                         # 迭代次数
    weights = ones((n,1))                   # 回归系数单位矩阵初始化
    for k in range(maxCycles):              # 迭代求解
        h = sigmoid(dataMatrix*weights)     # 求sigmoid函数值
        error = (labelMat - h)              # 向量相减，计算误差
        weights = weights + alpha * dataMatrix.transpose()* error # 梯度上升算法公式
    return weights.getA()                   # 将矩阵转换为数组，返回权重数组

2.2随机梯度上升

由2.1可知梯度上升算法在每次更新回归系数时都需要遍历整个数据集，如果有较多的样本或特征，则计算复杂度会很高，改进的方法为每次仅用一个样本点来更新回归系数。由于是单个样本更新，则可用于增量式的学习，有新样本就可以对分类器进行更新，即一个在线学习算法。

代码实现如下：

def stocGradAscent0(dataMatrix, classLabels):
    '''
    随机梯度上升法
    :param dataMatIn: 数据集
    :param classLabels: 数据标签
    :return: 回归系数
    '''
    m,n = shape(dataMatrix)                         # 返回dataMatrix行数和列数
    alpha = 0.01                                    # 移动步长
    weights = ones(n)                               # 回归系数单位矩阵初始化
    for i in range(m):                              # 每次使用一个样本
        h = sigmoid(sum(dataMatrix[i]*weights))     # 求sigmoid函数值
        error = classLabels[i] - h                  # 向量相减，计算误差
        weights = weights + alpha * error * dataMatrix[i]  # 随机梯度上升算法公式
    return weights

2.3改进的随机梯度上升

在随机梯度上升基础上，加入迭代求解，并使alpha随着迭代不断减少，但不会为0；以保证随着迭代的进行，不会产生大幅度的波动，从而收敛到某个值，加快收敛速度。

代码实现如下：

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    '''
    改进的随机梯度上升算法
    :param dataMatIn: 数据集
    :param classLabels: 数据标签
    :param numIter: 迭代次数
    :return: 回归系数
    '''
    m,n = shape(dataMatrix)                         # 返回dataMatrix行数和列数
    weights = ones(n)                               # 回归系数单位矩阵初始化
    for j in range(numIter):                        # 迭代求解
        dataIndex = list(range(m))                  # 数据行数索引
        for i in range(m):                          # 样本迭代
            alpha = 4/(1.0+j+i)+0.0001              # 随着迭代alpha不断减小，但不会为0
            randIndex = int(random.uniform(0,len(dataIndex)))   # 随机选取更新
            h = sigmoid(sum(dataMatrix[randIndex]*weights))     # 求sigmoid函数值
            error = classLabels[randIndex] - h      # 向量相减，计算误差
            weights = weights + alpha * error * dataMatrix[randIndex] # 随机梯度上升算法公式
            del(dataIndex[randIndex])               # 删除已经使用过的样本
    return weights

3.实战：鸢尾花分类

3.1使用sklearn中的LogisticRegression()

实现代码如下：

#!/usr/bin/python
# -*- coding:utf-8 -*-

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.pipeline import Pipeline

if __name__ == "__main__":
    path = 'iris.data'  # 数据文件路径
    data = pd.read_csv(path, header=None)
    # 使用pandas的Categorical对标签进行编码
    data[4] = pd.Categorical(data[4]).codes
    # 将数据分割为特征集和标签集
    x, y = np.split(data.values, (4,), axis=1)
    lr = Pipeline([('sc', StandardScaler()),                    # 去均值和方差归一化
                   ('poly', PolynomialFeatures(degree=2)),      # 对数据进行2阶多项式转换
                   ('clf', LogisticRegression()) ])             # 逻辑回归分类器
    lr.fit(x, y.ravel())                                        # 训练模型
    y_hat = lr.predict(x)                                       # 得到预测集
    np.set_printoptions(suppress=True)
    print('准确度：%.2f%%' % (100*np.mean(y_hat == y.ravel())))

输出结果如下：

准确度：98.00%

尝试直接用LogisticRegression()而不进行归一化和多项式转换，准确度降为96%。

若仅使用前两列特征（花萼宽度，花萼长度），进行计算，准确度：81.33%；画出分类效果图，如下所示：
这里写图片描述

作图代码如下：

    # 画图
    N, M = 500, 500                                 # 横纵各采样多少个值
    x1_min, x1_max = x[:, 0].min(), x[:, 0].max()   # 第0列的范围
    x2_min, x2_max = x[:, 1].min(), x[:, 1].max()   # 第1列的范围
    t1 = np.linspace(x1_min, x1_max, N)
    t2 = np.linspace(x2_min, x2_max, M)
    x1, x2 = np.meshgrid(t1, t2)                    # 生成网格采样点
    x_test = np.stack((x1.flat, x2.flat), axis=1)   # 测试点

    mpl.rcParams['font.sans-serif'] = ['simHei']
    mpl.rcParams['axes.unicode_minus'] = False
    cm_light = mpl.colors.ListedColormap(['#77E0A0', '#FF8080', '#A0A0FF'])
    cm_dark = mpl.colors.ListedColormap(['g', 'r', 'b'])
    y_hat = lr.predict(x_test)                      # 预测值
    y_hat = y_hat.reshape(x1.shape)                 # 使之与输入的形状相同
    plt.figure(facecolor='w')
    plt.pcolormesh(x1, x2, y_hat, cmap=cm_light)    # 预测值的显示
    plt.scatter(x[:, 0], x[:, 1], c=np.squeeze(y), edgecolors='k', s=50, cmap=cm_dark)    # 样本的显示
    plt.xlabel(u'花萼长度', fontsize=14)
    plt.ylabel(u'花萼宽度', fontsize=14)
    plt.xlim(x1_min, x1_max)
    plt.ylim(x2_min, x2_max)
    plt.grid()
    patchs = [mpatches.Patch(color='#77E0A0', label='Iris-setosa'),
              mpatches.Patch(color='#FF8080', label='Iris-versicolor'),
              mpatches.Patch(color='#A0A0FF', label='Iris-virginica')]
    plt.legend(handles=patchs, fancybox=True, framealpha=0.8)
    plt.title(u'鸢尾花Logistic回归分类效果 - 标准化', fontsize=17)
    plt.show()

3.2使用sklearn中的LogisticRegressionCV()并降维处理

实现代码如下：

# -*- coding:utf-8 -*-

import pandas as pd
import numpy as np
from sklearn.decomposition import PCA
from sklearn.feature_selection import SelectKBest, SelectPercentile, chi2
from sklearn.linear_model import LogisticRegressionCV
from sklearn import metrics
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches


def extend(a, b):
    return 1.05*a-0.05*b, 1.05*b-0.05*a


if __name__ == '__main__':
    stype = 'pca'
    pd.set_option('display.width', 200)
    data = pd.read_csv('.\\iris.data', header=None)
    columns = np.array(['花萼长度', '花萼宽度', '花瓣长度', '花瓣宽度', '类型'])
    data.rename(columns=dict(list(zip(np.arange(5), columns))), inplace=True)
    data['类型'] = pd.Categorical(data['类型']).codes
    print(data.head(5))
    x = data[columns[:-1]]
    y = data[columns[-1]]

    if stype == 'pca':
        pca = PCA(n_components=2, whiten=True, random_state=0)      # 利用PCA降为2维
        x = pca.fit_transform(x)
        print('各方向方差：', pca.explained_variance_)
        print('方差所占比例：', pca.explained_variance_ratio_)
        x1_label, x2_label = '组分1', '组分2'
        title = '鸢尾花数据PCA降维'
    else:
        fs = SelectKBest(chi2, k=2)                                 # 选取两个特征
        # fs = SelectPercentile(chi2, percentile=60)
        fs.fit(x, y)
        idx = fs.get_support(indices=True)
        print('fs.get_support() = ', idx)
        x = x[columns[idx]]
        x = x.values                    # 为下面使用方便，DataFrame转换成ndarray
        x1_label, x2_label = columns[idx]
        title = '鸢尾花数据特征选择'
    print(x)

    # 拆分数据为训练集（70%）和测试集（30%）
    x, x_test, y, y_test = train_test_split(x, y, test_size=0.3)
    model = Pipeline([
        ('poly', PolynomialFeatures(degree=2, include_bias=True)),                          # 对数据进行2阶多项式转换
        ('lr', LogisticRegressionCV(Cs=np.logspace(-3, 4, 8), cv=5, fit_intercept=False))   # 逻辑回归分类器（交叉验证）
    ])
    # 训练模型
    model.fit(x, y)
    print('最优参数：', model.get_params('lr')['lr'].C_)
    # 模型预测
    y_hat = model.predict(x)
    print('训练集精确度：', metrics.accuracy_score(y, y_hat))
    y_test_hat = model.predict(x_test)
    print('测试集精确度：', metrics.accuracy_score(y_test, y_test_hat))

    # 画图
    N, M = 500, 500     # 横纵各采样多少个值
    x1_min, x1_max = extend(x[:, 0].min(), x[:, 0].max())   # 第0列的范围
    x2_min, x2_max = extend(x[:, 1].min(), x[:, 1].max())   # 第1列的范围
    t1 = np.linspace(x1_min, x1_max, N)
    t2 = np.linspace(x2_min, x2_max, M)
    x1, x2 = np.meshgrid(t1, t2)                    # 生成网格采样点
    x_show = np.stack((x1.flat, x2.flat), axis=1)   # 测试点
    y_hat = model.predict(x_show)  # 预测值
    y_hat = y_hat.reshape(x1.shape)  # 使之与输入的形状相同
    plt.figure(facecolor='w')
    plt.pcolormesh(x1, x2, y_hat, cmap=cm_light)  # 预测值的显示
    plt.scatter(x[:, 0], x[:, 1], s=30, c=y, edgecolors='k', cmap=cm_dark)  # 样本的显示
    plt.xlabel(x1_label, fontsize=12)
    plt.ylabel(x2_label, fontsize=12)
    plt.xlim(x1_min, x1_max)
    plt.ylim(x2_min, x2_max)
    plt.grid(b=True, ls=':', color='k')
    patchs = [mpatches.Patch(color='#77E0A0', label='Iris-setosa'),
              mpatches.Patch(color='#FF8080', label='Iris-versicolor'),
              mpatches.Patch(color='#A0A0FF', label='Iris-virginica')]
    plt.legend(handles=patchs, fancybox=True, framealpha=0.8, loc='lower right')
    plt.title('鸢尾花Logistic回归分类效果', fontsize=15)
    plt.show()

A.使用PCA降维，输出结果如下：

各方向方差： [ 4.22484077  0.24224357]
方差所占比例： [ 0.92461621  0.05301557]
最优参数： [ 10.  10.   1.]
训练集精确度： 0.961904761905
测试集精确度： 0.955555555556

这里写图片描述

B.使用SelectKBest特征选择，输出结果如下：

fs.get_support() =  [2 3]
最优参数： [  0.1  10.   10. ]
训练集精确度： 0.942857142857
测试集精确度： 1.0