1.MNIST

import sklearn
import numpy as np

# 随机种子
np.random.seed(42)

%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rc('axes', labelsize=14)
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)

# 下载手写数字集
from sklearn.datasets import fetch_openml
mnist = fetch_openml('mnist_784', version=1)
mnist.keys()
type(mnist)

【输出】：

dict_keys(['data', 'target', 'feature_names', 'DESCR', 'details', 'categories', 'url'])
sklearn.utils.Bunch

Bunch本质上的数据类型是dict，属性有：

1. DESCR：数据描述。
2. target_names：标签名。可自定义，默认为文件夹名。
3. filenames：文件名。
4. target：文件分类。
5. data：数据数组。

2.查看数据基本情况

2.1 查看数据集特征

# 查看数据特征
X= mnist["data"]
print(X.shape)

# 将数据特征的784数组形式，转换成28*28的图片形式
some_digit = X[0]  
some_digit_image = some_digit.reshape(28, 28)
plt.imshow(some_digit_image, cmap=mpl.cm.binary)
plt.axis("off")
plt.show()

【输出】：

(70000, 784)

2.2 查看数据集标签

y= mnist["target"]
print(y.shape)
print(y[0])       # 查看第1张照片所对应的标签
print（type(y[0])）

# 提前将，标签的数据类型转换成8位整型
y = y.astype(np.uint8)
print(type(y[0])）

【输出】：

(70000,)
'5'
str

2.3 定义相关函数

# 定义显示一张图片的函数
def plot_digit(data):
    image = data.reshape(28, 28)
    plt.imshow(image, cmap = mpl.cm.binary,
               interpolation="nearest")
    plt.axis("off")
	
# 定义显示多张图片的函数
def plot_digits(instances, images_per_row=10, **options):
    size = 28
    images_per_row = min(len(instances), images_per_row)
    images = [instance.reshape(size,size) for instance in instances]
    n_rows = (len(instances) - 1) // images_per_row + 1
    row_images = []
    n_empty = n_rows * images_per_row - len(instances)
    images.append(np.zeros((size, size * n_empty)))
    for row in range(n_rows):
        rimages = images[row * images_per_row : (row + 1) * images_per_row]
        row_images.append(np.concatenate(rimages, axis=1))
    image = np.concatenate(row_images, axis=0)
    plt.imshow(image, cmap = mpl.cm.binary, **options)
    plt.axis("off")

# 显示多张图片
plt.figure(figsize=(9,9))
example_images = X[:100]
plot_digits(example_images, images_per_row=10)
plt.show()

3.创建训练集与测试集

数据集本身就已经是打乱的：

X_train, X_test, y_train, y_test = X[:60000], X[60000:], y[:60000], y[60000:]

4.简化问题：二分类

y_train_5 = (y_train == 5) # 返回布尔类型
y_test_5 = (y_test == 5)
y_train_5[:10]             # 查看前10个

array([ True, False, False, False, False, False, False, False, False,
       False])

5.SGD训练的线性分类器

5.1 模型创建及训练

该估计器通过随机梯度下降(SGD)学习实现正则化线性模型：一次估计每个样本的损失梯度，并随着学习率的降低而更新模型。

from sklearn.linear_model import SGDClassifier

sgd_clf = SGDClassifier(max_iter=1000, tol=1e-3, random_state=42) # max_iter：最大的迭代次数；random_state=42：让随机数起点一样
sgd_clf.fit(X_train, y_train_5)

【输出】：

SGDClassifier(alpha=0.0001, average=False, class_weight=None,
              early_stopping=False, epsilon=0.1, eta0=0.0, fit_intercept=True,
              l1_ratio=0.15, learning_rate='optimal', loss='hinge',
              max_iter=1000, n_iter_no_change=5, n_jobs=None, penalty='l2',
              power_t=0.5, random_state=42, shuffle=True, tol=0.001,
              validation_fraction=0.1, verbose=0, warm_start=False)

5.2 模型预测

# 取出第一条数据
some_digit = X[0]

# 将其格式转换，显示照片
some_digit_image = some_digit.reshape(28, 28)
plt.imshow(some_digit_image, cmap=mpl.cm.binary)
plt.axis("off")
plt.show()

# 模型预测结果
sgd_clf.predict([some_digit])

【输出】：

array([ True])

some_digit = X[1]
some_digit_image = some_digit.reshape(28, 28)
plt.imshow(some_digit_image, cmap=mpl.cm.binary)
plt.axis("off")
plt.show()

sgd_clf.predict([some_digit])

【输出】：

array([False])

5.3 交叉验证

# 利用sklearn中的cross_val_score进行交叉验证
from sklearn.model_selection import cross_val_score

cross_val_score(sgd_clf, X_train, y_train_5, cv=3, scoring="accuracy")  # 3折

【输出】：

array([0.95035, 0.96035, 0.9604 ])

【说明】：3次正确率都超过了95% 【ps：当然也可以自己写一个交叉验证】

from sklearn.model_selection import StratifiedKFold
from sklearn.base import clone

# 分层抽样
skfolds = StratifiedKFold(n_splits=3, random_state=42)

for train_index, test_index in skfolds.split(X_train, y_train_5):
    clone_clf = clone(sgd_clf)
    X_train_folds = X_train[train_index]
    y_train_folds = y_train_5[train_index]
    X_test_fold = X_train[test_index]
    y_test_fold = y_train_5[test_index]

    clone_clf.fit(X_train_folds, y_train_folds)
    y_pred = clone_clf.predict(X_test_fold)
    n_correct = sum(y_pred == y_test_fold)
    print(n_correct / len(y_pred))

【输出】：

0.95035
0.96035
0.9604

6.以准确率评估分类器的缺点

根据上面的交叉验证发现准确率都超过了95%，我们的模型看似很完美，但是：

from sklearn.base import BaseEstimator

# 定义傻瓜分类器
class Never5Classifier(BaseEstimator):
    def fit(self, X, y=None):
        pass
    def predict(self, X):
        return np.zeros((len(X), 1), dtype=bool)

# 创建傻瓜分类器模型
never_5_clf = Never5Classifier()

# 交叉验证
cross_val_score(never_5_clf, X_train, y_train_5, cv=3, scoring="accuracy")

array([0.91125, 0.90855, 0.90915])

【说明】：发现傻瓜分类器的准确率竟然也高达90%！这是因为只有大约10%的图像是5，所以如果你一直猜不是5，90%的可能性你是对的！所以，准确率无法成为分类器的首要性能指标。

7.混淆矩阵、ROC曲线、

# 利用sklearn中的cross_val_predict进行交叉验证
# 但是与cross_val_socre不同，它返回的是每个折叠的预测
from sklearn.model_selection import cross_val_predict

y_train_pred = cross_val_predict(sgd_clf, X_train, y_train_5, cv=3)
y_train_pred

【输出】：

array([ True, False, False, ...,  True, False, False])

7.1 混淆矩阵

from sklearn.metrics import confusion_matrix

# 获取混淆矩阵
confusion_matrix(y_train_5, y_train_pred)

【输出】：

array([[53892,   687],
       [ 1891,  3530]], dtype=int64)

# 下面是最理想的混淆矩阵
y_train_perfect_predictions = y_train_5
confusion_matrix(y_train_5, y_train_perfect_predictions)

array([[54579,     0],
       [    0,  5421]], dtype=int64)

7.2 查准率与查全率

precision：查准率、精度

recall：查全率、召回率

# 先利用sklearn中现成的算出，再用算式验证一下：
from sklearn.metrics import precision_score, recall_score

print('查准率：', precision_score(y_train_5, y_train_pred))
print('验证一下对不对：', 3530/(687+3530))

print('查全率：', recall_score(y_train_5, y_train_pred))

查准率： 0.8370879772350012
验证一下对不对： 0.8370879772350012

查全率： 0.6511713705958311

7.3 F1分数

F1分数是精度和召回率的谐波平均值，只有当两者都很高的时候，分类器才能得到较高的F1分数。但是，并不是所有情况都一直符合你的要求，可能你的更关心精度、可能你更关心召回率，鱼和熊掌不可兼得。

from sklearn.metrics import f1_score

f1_score(y_train_5, y_train_pred)

0.7325171197343846

7.4 精度/召回率的权衡关系

我们创建的二分类模型是SGDClassifier，它在每次预测时会基于决策函数计算出一个分值，如果该分值大于阈值，则将该实例判为正类，否则判为负类。如下图，我们从左到右依次提升阈值，不难分析出模型的精度会随着阈值的升高而升高，但是召回率会随着阈值的升高而下降，这就是两者的权衡关系。

7.5 阈值该如何选择

# 取出第一条数据
some_digit = X[0]

# 利用模型计算出第一条数据的得分
y_scores = sgd_clf.decision_function([some_digit])
print('得分：', y_scores)

# 间接设置阈值
threshold = 0
y_some_digit_pred = (y_scores > threshold)
print('阈值为0时的判断结果：', y_some_digit_pred)

# 间接设置阈值
threshold = 2500
y_some_digit_pred = (y_scores > threshold)
print('阈值为2500时的判断结果：', y_some_digit_pred)

【输出】：

得分： [2164.22030239]
阈值为0时的判断结果： [ True]
阈值为2500时的判断结果： [False]

# 利用cross_val_predict获取训练集中所有实例的分数
# 但这次返回的是预测分数，并不是预测结果！
y_scores = cross_val_predict(sgd_clf, X_train, y_train_5, cv=3,
                             method="decision_function")
print(y_scores)
print(y_scores.shape)

【输出】：

array([  1200.93051237, -26883.79202424, -33072.03475406, ...,
        13272.12718981,  -7258.47203373, -16877.50840447])
(60000,)

画出决策阈值与精度、召回率的图像：

# 首先，先通过刚刚获取的决策分数，利用precision_recall_curve算出所有可能的阈值的精度和召回率
from sklearn.metrics import precision_recall_curve

precisions, recalls, thresholds = precision_recall_curve(y_train_5, y_scores)

# 然后，根据手中的阈值、精度、召回率描绘图像
def plot_precision_recall_vs_threshold(precisions, recalls, thresholds):
    plt.plot(thresholds, precisions[:-1], "b--", label="Precision", linewidth=2)
    plt.plot(thresholds, recalls[:-1], "g-", label="Recall", linewidth=2)
    plt.legend(loc="center right", fontsize=16)
    plt.xlabel("Threshold", fontsize=16)
    plt.grid(True)
    plt.axis([-50000, 50000, 0, 1])

plt.figure(figsize=(8, 4))
plot_precision_recall_vs_threshold(precisions, recalls, thresholds)
plt.show()

print('精度的个数：', len(precisions))
print('精度大于0.90，返回对应的布尔类型数组：', precisions >= 0.90)
print('返回精度大于0.90的最大的最大索引', np.argmax(precisions >= 0.90))
print('返回精度大于0.90的最大的最大索引位置对应的召回率：', recalls[np.argmax(precisions >= 0.90)])
print('返回精度大于0.90的最大的最大索引位置对应的召回率：', recall_score(y_train_5, y_train_pred_90))   # 两个表达一样
print('返回精度大于0.90的最大的最大索引位置对应的决策阈值：', thresholds[np.argmax(precisions >= 0.90)])
print('返回精度大于0.90的最大的最大索引位置对应的查准率：', precision_score(y_train_5, y_train_pred_90))

精度的个数： 59967
精度大于0.90，返回对应的布尔类型数组： [False False False ...  True  True  True]
返回精度大于0.90的最大的最大索引 57075
返回精度大于0.90的最大的最大索引位置对应的召回率： 0.4799852425751706
返回精度大于0.90的最大的最大索引位置对应的召回率： 0.4799852425751706
返回精度大于0.90的最大的最大索引位置对应的决策阈值： 3370.0194991439557
返回精度大于0.90的最大的最大索引位置对应的查准率： 0.9000345901072293

7.6 ROC曲线及AUC面积

ROC曲线也经常与二元分类器一起使用，叫做受试者工作特征曲线，面积越大模型效果越好。

# 首先，利用roc_curve计算多种阈值的真正类率TPR（召回率）、假正类率FPR
from sklearn.metrics import roc_curve

fpr, tpr, thresholds = roc_curve(y_train_5, y_scores)

# 然后，利用得到的FPR、TPR描绘图像
def plot_roc_curve(fpr, tpr, label=None):
    plt.plot(fpr, tpr, linewidth=2, label=label)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.axis([0, 1, 0, 1])
    plt.xlabel('False Positive Rate (Fall-Out)', fontsize=16)
    plt.ylabel('True Positive Rate (Recall)', fontsize=16)
    plt.grid(True)

plt.figure(figsize=(8, 6))
plot_roc_curve(fpr, tpr)
plt.show()

# 利用roc_auc_score计算ROC曲线的AUC面积
from sklearn.metrics import roc_auc_score

roc_auc_score(y_train_5, y_scores)

0.9604938554008616

【说明】：AUC面积的理想值是1，纯随机分类器ROC曲线的AUC是0.5，用于参照。

7.7 创建一个随机森林分类器比较ROC曲线

from sklearn.ensemble import RandomForestClassifier

# 随机森林分类器
forest_clf = RandomForestClassifier(n_estimators=100, random_state=42)
y_probas_forest = cross_val_predict(forest_clf, X_train, y_train_5, cv=3,
                                    method="predict_proba")  # predict_proba：概率
y_scores_forest = y_probas_forest[:, 1]

# 计算随机森林分类器的多种阈值的真正类率TPR（召回率）、假正类率FPR
fpr_forest, tpr_forest, thresholds_forest = roc_curve(y_train_5,y_scores_forest)

# 描绘图像
plt.figure(figsize=(8, 6))
plt.plot(fpr, tpr, "b:", linewidth=2, label="SGD")
plot_roc_curve(fpr_forest, tpr_forest, "Random Forest")
plt.grid(True)
plt.legend(loc="lower right", fontsize=16)
plt.show()

8.多分类问题的解决思路

一些算法支支持严格的二分类，如：支持向量机分类器、线性分类器等，这种情况也可以使用二元分类器实现多分类问题：

1. OvA一对多策略：例如：0与其他、1与其他、2与其他......共10个分类器。
2. OvO一对一策略：为每一对数字训练一个二分类器，例如：区分0和1、区分0和2、区分1和2、....共45个二分类器。

还有一些算法本身就支持多分类，如：随机森林分类器、朴素贝叶斯分类器等。

【ps：sklearn可以检测到你尝试使用二元分类器进行多分类任务，自动运行OvA(SVM分类器除外)】

9.支持向量机、随机梯度下降SGD分类器

9.1 支持向量机

自动实现一对多

# 创建支持向量机分类器
from sklearn.svm import SVC

svm_clf = SVC(gamma="auto", random_state=42)
svm_clf.fit(X_train[:1000], y_train[:1000])  # 只选了1000个训练

# 取出第一条数据进行预测
some_digit = X[0]
print('第一条数据预测为：', svm_clf.predict([some_digit]))
some_digit_scores = svm_clf.decision_function([some_digit])
print('第一条数据的预测分数为：', some_digit_scores)
print('第一条数据的预测分数中最大值所对应的索引：', np.argmax(some_digit_scores))
print('类别：', svm_clf.classes_)

# 查看分类器的数量
from sklearn.multiclass import OneVsRestClassifier

ovr_clf = OneVsRestClassifier(SVC(gamma="auto", random_state=42))
ovr_clf.fit(X_train[:1000], y_train[:1000])

【输出】：

第一条数据预测为： [5]
第一条数据的预测分数为： [[ 2.81585438  7.09167958  3.82972099  0.79365551  5.8885703   9.29718395
   1.79862509  8.10392157 -0.228207    4.83753243]]
第一条数据的预测分数中最大值所对应的索引： 5
类别： [0 1 2 3 4 5 6 7 8 9]

强制实现一对一

# 强制一对一
from sklearn.multiclass import OneVsRestClassifier

ovr_clf = OneVsRestClassifier(SVC(gamma="auto", random_state=42))
ovr_clf.fit(X_train[:1000], y_train[:1000])

# 预测
print('预测第一条数据的结果：', ovr_clf.predict([some_digit]))

# 分类器数量
print('分类器数量：', len(ovr_clf.estimators_))

【输出】：

预测第一条数据的结果： [5]
分类器数量： 10

9.2 随机梯度下降SGD分类器

from sklearn.linear_model import SGDClassifier

# 随机梯度下降SGD分类器
sgd_clf = SGDClassifier(max_iter=1000, tol=1e-3, random_state=42)
sgd_clf.fit(X_train, y_train)

# 模型预测
some_digit = X_train[0]
print('预测的类别是：', sgd_clf.predict([some_digit]))
print('每个类别的评分是：', sgd_clf.decision_function([some_digit]))

# 交叉验证
cross_val_score(sgd_clf, X_train, y_train, cv=3, scoring='accuracy')

【输出】：

预测的类别是： [3]
每个类别的评分是： [[-31893.03095419 -34419.69069632  -9530.63950739   1823.73154031
  -22320.14822878  -1385.80478895 -26188.91070951 -16147.51323997
   -4604.35491274 -12050.767298  ]]
array([0.87082583, 0.87089354, 0.88628294])

【说明】：纯随机分类器的准确率大概是10%，所以目前这个结果还算可以，如果想再提高一些，可以标准化特征。

# 将特征简单缩放，准确率还会上升一点
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train.astype(np.float64))
cross_val_score(sgd_clf, X_train, y_train, cv=3, scoring='accuracy')

【输出】：

array([0.89707059, 0.8960948 , 0.90693604])

10.错误分析

获取多分类的混淆矩阵

# 获取多分类的混淆矩阵
y_train_pred = cross_val_predict(sgd_clf, X_train_scaled, y_train, cv=3)
conf_mx = confusion_matrix(y_train, y_train_pred)
conf_mx

array([[5578,    0,   22,    7,    8,   45,   35,    5,  222,    1],
       [   0, 6410,   35,   26,    4,   44,    4,    8,  198,   13],
       [  28,   27, 5232,  100,   74,   27,   68,   37,  354,   11],
       [  23,   18,  115, 5254,    2,  209,   26,   38,  373,   73],
       [  11,   14,   45,   12, 5219,   11,   33,   26,  299,  172],
       [  26,   16,   31,  173,   54, 4484,   76,   14,  482,   65],
       [  31,   17,   45,    2,   42,   98, 5556,    3,  123,    1],
       [  20,   10,   53,   27,   50,   13,    3, 5696,  173,  220],
       [  17,   64,   47,   91,    3,  125,   24,   11, 5421,   48],
       [  24,   18,   29,   67,  116,   39,    1,  174,  329, 5152]])

【说明】：由于多分类的混淆矩阵不太明显，我们可以利用可视化将他描述出来。

def plot_confusion_matrix(matrix):
    fig = plt.figure(figsize=(8,8))
    ax = fig.add_subplot(111)
    cax = ax.matshow(matrix)
    fig.colorbar(cax)

plt.matshow(conf_mx, cmap=plt.cm.gray)
plt.show()

【说明】：由于我们模型判断正确的很多，混淆矩阵对角线上格外明亮，而对角线以外的部分数值较小接近黑色。因此，我们需要做一下处理：

1. 首先，将矩阵中的每个值除以相对应类别的图片数目(这样比较的就是错误率而不是错误的绝对值，相对公平)
2. 然后，用0填充对角线，重新绘制结果

# 第一步
row_sums = conf_mx.sum(axis=1, keepdims=True)
norm_conf_mx = conf_mx / row_sums

# 第二布
np.fill_diagonal(norm_conf_mx, 0)
plt.matshow(norm_conf_mx, cmap=plt.cm.gray)
save_fig("confusion_matrix_errors_plot", tight_layout=False)
plt.show()

【说明】：每行代表实际类别，每列代表预测类别。第8列和第9列整体比较明亮，很多照片被错误的分成了8、9

# 查看数字3、5的例子
cl_a, cl_b = 3, 5
X_aa = X_train[(y_train == cl_a) & (y_train_pred == cl_a)]
X_ab = X_train[(y_train == cl_a) & (y_train_pred == cl_b)]
X_ba = X_train[(y_train == cl_b) & (y_train_pred == cl_a)]
X_bb = X_train[(y_train == cl_b) & (y_train_pred == cl_b)]

# 把图片显示出来
plt.figure(figsize=(8,8))
plt.subplot(221); plot_digits(X_aa[:25], images_per_row=5)
plt.subplot(222); plot_digits(X_ab[:25], images_per_row=5)
plt.subplot(223); plot_digits(X_ba[:25], images_per_row=5)
plt.subplot(224); plot_digits(X_bb[:25], images_per_row=5)
save_fig("error_analysis_digits_plot")
plt.show()

【说明】：左侧两个55的矩阵显示了被分类成3的图片，右侧两个55的矩阵显示了被分成5的图片。出错的原因在于，我们使用简单的SGDClassifier模型是一个线性模型，它所做的就是为每个像素分配一个各个类别的权重，当它看见新的图像时，将加权后的像素强度汇总，从而得到一个分数进行分类。但是图中的3、5只是在一部分像素上有区别，所以分类器很容易将它们弄混。

数字3和数字5之间的主要区别在于连接顶线和下方弧线的中间那段小线条的位置。如果你写的数字3将连接点略往左移，分类器可能就将其分类为数字5，反之亦然。总之，这个分类器对图像的移位和旋转非常的敏感。因此，减少3和5的混淆方法之一，就是对图片进行预处理，保证它们位于中心位置并且没有旋转。