获取数据

# from sklearn.datasets import fetch_mldata
# from sklearn import datasets

# mnist = fetch_mldata('MNIST   original') 
# mnist

好像下载不到它的数据,直接从网上找到它的数据,放到当面目录下的dataset目录下。

from sklearn.datasets import fetch_mldata
from sklearn import datasets
import numpy as np

mnist = fetch_mldata('mnist-original', data_home = './datasets/') 
mnist
{'DESCR': 'mldata.org dataset: mnist-original',
 'COL_NAMES': ['label', 'data'],
 'target': array([0., 0., 0., ..., 9., 9., 9.]),
 'data': array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint8)}

```

网上很多的说法是错误的,只有我这个才是正解。

X, y = mnist['data'], mnist['target']
print(X.shape)
print(y.shape)
(70000, 784)
(70000,)

从上面看出来,X是一个\(7000\times784\)的一个矩阵,一般来说,7000行表示有7000个样本,784列,表示样本有784这么多个属性。

%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt

some_digit = X[36000]
some_digit_image = some_digit.reshape(28,28)

plt.imshow(some_digit_image, cmap=matplotlib.cm.binary, interpolation="nearest")
plt.axis('off')
plt.show()

分类-MNIST-LMLPHP

说个数看起来像是5,我觉得更像是6,我们可查看一下它的标签。

y[36000]
5.0
# EXTRA
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 = matplotlib.cm.binary, **options)
    plt.axis("off")
plt.figure(figsize=(9,9))
example_images = np.r_[X[:12000:600], X[13000:30600:600], X[30600:60000:590]]
plot_digits(example_images, images_per_row=10)
# save_fig("more_digits_plot")
plt.show()

分类-MNIST-LMLPHP

可能这个标签写错了都不一定,我们得新写一下这个标签,说不定可以提高模型的准确率呢。这只是我个人在这里开玩笑说的,不用当真哈。

在做数据的训练前,应该找出测试集,这里MNIST已经帮我们把测试集做好了。

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

MNIST的数据是按数字大小顺序排列的,所我们先要打乱它的顺序,这样可以保证我们的交叉验证是每一次都是相似的。

import numpy as np

shuffle_index = np.random.permutation(60000)
shuffle_index
array([52603, 56601, 42625, ..., 17778, 24267, 29358])
X_train, y_train = X_train[shuffle_index],y_train[shuffle_index]

训练一个二分类器

先不做一个多类器,我们先不去识别里面的手写数字是0~10中的某一个数。目前做一个最简单的,判断它是否是5,即将数据分成两个类别:“5”和“非5”

# 这是一个逻辑数组,5:True, 非5:False
y_train_5 = (y_train == 5)
y_test_5 = (y_test == 5)

现在开始用一个分类器去训练它。用随机梯度下降分类器SGD。用Scikit-Learn的SGDClassifier类。这个分类器有一个好处是能够高效地处理非常大的数据集。部分原因是它每次只处理一条数据。

from sklearn.linear_model import SGDClassifier

sgd_clf = SGDClassifier(random_state = 32)
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=None,
       n_iter=None, n_iter_no_change=5, n_jobs=None, penalty='l2',
       power_t=0.5, random_state=32, shuffle=True, tol=None,
       validation_fraction=0.1, verbose=0, warm_start=False)
sgd_clf.predict([some_digit])
array([ True])

这个模型的准确度好像受随机种子的影响比较大,如果我将模型的随机种改为42,我们再来看一下它预测的结果是不是正确的

sgd_clf = SGDClassifier(random_state = 42)
sgd_clf.fit(X_train, y_train_5)
sgd_clf.predict([some_digit])
array([ True])

对性能的评估

可以到这个时候它又预测错了。下面来整体评估一下这个分类的性能。

使用交叉验证测量准确性

在交叉验证过程中,有时候我们会需要更多的控制权,相较于函数cross_val_score()或者其他相似函数所提供的功能。下面代码做了和cross_val_score()相同的事情

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

skfolds = StratifiedKFold(n_splits = 3, random_state = 42)
clone_clf = clone(sgd_clf)
for train_index, test_index in skfolds.split(X_train, y_train_5):
    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.9612

0.9531

0.9688

下面直接使用sklearn中的库进行交叉评估。使用cross_val_score函数来评估SGDClassifier模型。

from sklearn.model_selection import cross_val_score

cross_val_score(sgd_clf, X_train, y_train_5, cv = 3, scoring = "accuracy")
array([0.9612, 0.9531, 0.9688])

这精度看起来还不错,有大于95%的精度,有点让人兴奋,感觉做个分类还是挺容易的,不难。
我们再来看下一个非常简单的分类器去分类,看看它在“非5”这个类上的表现。

from sklearn.base import BaseEstimator
# 这个模型的预测的策略就是将所有的数据都认为是'非5'
class Never5Classifier(BaseEstimator):
    def fit(self,X,y=None):
        pass
    def predict(self,X):
        return np.zeros((len(X),1), dtype=bool)
np.zeros((2,1), dtype=bool)
array([[False],
       [False]])
never_5_clf = Never5Classifier()
cross_val_score(never_5_clf, X_train, y_train_5, cv = 3, scoring = "accuracy")
array([0.90815, 0.9124 , 0.9084 ])

这么一个简单的分类器也有90%的精度,这是因为只有10%的样本是5,其它都是非5,所以只我们一直猜这个图像不是5,当然有90%的精度,这叫数据不平衡。就像我们如果在日本,站到大街上,见到人就猜他是一个日本人,我们几乎肯定是正确的。

所以精度并不是一个好的性能度量指标,特别是在我们数据不平衡的时候。

混淆矩阵

对一般分类器来说,一人好得多的性能评估指标是混淆矩阵。大体思路是:输出类别A被分成类别B的次数。

为了计算混淆矩阵,首先你需要有一系列的预测值,这样才能将预测值与真实值做比较。你或许想在测试集上做预测。

from sklearn.model_selection import cross_val_predict

y_train_pred = cross_val_predict(sgd_clf, X_train, y_train_5, cv = 3)
from sklearn.metrics import confusion_matrix

confusion_matrix(y_train_5, y_train_pred)
array([[54306,   273],
       [ 2065,  3356]], dtype=int64)

混淆矩阵中的每一行表示一个实际的类,而每一列表一个预测的类。该矩阵的第一行认为"非5"中的53993张被正确地归类为非5(这被称为真反例,true negatives),而其余586被错误归类为5(这被称为假正例,false positive),其余3905正确分类为"5"类(真正例,true positive)。一个完美的分类器将只有真反例和真正例,所混淆矩阵的非零值仅在其主对角线(左上至右下)。

# confusion_matrix(y_train_5, y_train_perfect_predictions)

混淆矩阵可以提供很多信息。有时候你会想要更加简明的指标。一个有趣的指标是正例预测的精度,也叫做分类器的准确率(precision)

\[precision = \frac{TP}{TP + FP}\tag{3-1}\]

其中\(TP\)真正例的数目,\(FP\)假正例的数目。

以准确率一般会伴随另一个指标一起使用,这个指标叫做召回率(recall),也叫做敏感度(sensitivity)或者真正例率(true positive rate, TPR)。这是正例被分类器正确探测出的比率。

\[recall = \frac{TP}{TP+FN}\tag{3-2}\]

\(FN\)是假反例的数目。

from sklearn.metrics import precision_score, recall_score

print(precision_score(y_train_5, y_train_pred))
print(recall_score(y_train_5, y_train_pred))
0.924772664645908
0.6190739715919572

这样看起,这个分类器的准确率并不高,只有56.8%左右,而且只是分成两类的一个分类器,这跟我们猜差不多。

通常结合准确率和召回率会更加方便,这个指标叫做F1值,特别是当你需要一个简单的方法去比较两个分类器的优劣的时时候。F1值是准确率和召回率的调和平均

\[F1 = \frac{2}{\frac{1}{precision}+\frac{1}{recall}} = 2 \times \frac{precision \times recall}{precision + recall} = \frac{TP}{TP + \frac{FN+FP}{2}}\tag{3-3}\]

计算F1值,简单调用f1_score()即可。

from sklearn.metrics import f1_score

f1_score(y_train_5, y_train_pred)
0.7416574585635358

F1支持那些有着相近准确率和召回率的分类(意思是只有当准确率和召回率一样大的时个,F1值才会大)。但并不是所的时候,我们都关心F1值,有时候我们只关心准确率(precision),或者有时候我们只关心召回率(recall)。

在这里,我以自己的理解,举两个例子,比如公司想找个人当总经理,有一群人来应聘它。我们这时候的目标是,找到的这个人肯定是能够当总经理的,就算有的人看起来像是能当总经理,但是为了确保万无一失,我们要找一个看起来非常非常像能够当总经理的人。这个时候我们当然有着很高的准确率,因为我们找的人几乎肯定是能够当总经理的,但是此时,我们会犯另一个错误,就是有些人确实有能力当总经理,只是我们没有看出来(人不可貌像),所以我们拒绝他,因此我们有低的召回率,这在统计学上被称为犯了第一类错误,即弃真。这样做是合理的,因为即使弃真,但我们保真了。

另一种情况是,比如警察在一群人中想找出几个犯罪的人,这个时候我们就不能要超高的准确率了,因为有可能把真正的犯人放走。找犯人的原则一般是,只要他看起来像个犯人,都应该审查一下,即使最后真像大白后,他真的不是一个犯人。我们平时听到的宁可错杀一千,不可放走一个说的就是这个道理,因此这有着比较低的准确率,但是有高的召回率,这在统计学上被称为犯了第二类错误,即取伪

准备率/召回率之间的折中

y_scores = sgd_clf.decision_function([some_digit])
y_scores
array([15905.22111141])
threshold = 0
y_some_digit_pred = (y_scores > threshold)
y_some_digit_pred
array([ True])
y_scores = cross_val_predict(sgd_clf, X_train,y_train_5,cv=3,
                             method = "decision_function")
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")
    plt.plot(thresholds, recalls[:-1], "g-", label = "Recall")
    plt.xlabel("Threshold")
    plt.legend(loc="upper left")
    plt.ylim([0,1.1])
    
plot_precision_recall_vs_threshold(precisions,recalls,thresholds)
plt.grid()
plt

分类-MNIST-LMLPHP

ROC曲线

受试者工作特征(ROC)曲线是另一个二分类器常用的工具。它非常类似与准确率/召回率曲线,但不是画出准确率对召回率的曲线,,ROC曲线是真正例率(true positive rate,另一个名字叫做召回率)对假正例率(false positive rate, FPR)的曲线。FPR是反例被错误分成正例的比率。它等于1减去真反例率(true negative rate,TNR)。TNR是反例被正确分类的比率。TNR也叫做特异性。

为了画出ROC曲线,你首先需要计算各种不同阈值下的TPR、FPR,使用roc_curve()函数:

from sklearn.metrics import roc_curve
fpr, tpr, thresholds = roc_curve(y_train_5, y_scores)

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')
    plt.ylabel('True Positive Rate')

plot_roc_curve(fpr,tpr)
plt

分类-MNIST-LMLPHP

一个比较分类器之间优劣的方法是:测量ROC曲线下的面积(AUC)。一个完美的分类器的 ROC AUC 等于1,而一个纯随机分类器的ROC AUC等于0.5。Scikit-Learn提供了一个函数来计算ROC AUC:

from sklearn.metrics import roc_auc_score
roc_auc_score(y_train_5,y_scores)
0.9623990527630832
from sklearn.ensemble import RandomForestClassifier

forest_clf = RandomForestClassifier(random_state = 42)
y_probas_forest = cross_val_predict(forest_clf, X_train, y_train_5, cv=3, method = "predict_proba")
y_scores_forest = y_probas_forest[:,1]
fpr_forest, tpr_forest, thresholds_forest=roc_curve(y_train_5,y_scores_forest)
plt.plot(fpr,tpr,"b:",label="SGD") 
plot_roc_curve(fpr_forest,tpr_forest,"Random Forest") 
plt.legend(loc="bottom  right") 
plt

分类-MNIST-LMLPHP

# 将概率大于0.5的,置为true, 否则为false
print(precision_score(y_train_5, y_scores_forest > 0.5))
print(recall_score(y_train_5, y_scores_forest > 0.5))
0.9844298245614035
0.8280760007378712

可以看出来,它的准确率还可,挺高的。

下面我们将分类出更多的数字,而不仅仅是5。

多类分类

二分类器只能分出两个类,而多分类器能分出多于两个类别的类。

一些算法(比如随机森林分类器或者朴素贝叶斯分类器)可以直接处理多类分类问题。其他一些算法(比如SVM分类器或者线性分类器)则是严格的二分类器,然后有许多策略可以让你用二分类器去执行多类分类。

Scikit-Learn可以探测出你想使用一个二分类器去完成多分类的任务,它会自动地执行OvA(除了SVM分类器,它使用OvO)。让我们试一下SGDClassifier

sgd_clf.fit(X_train, y_train)
sgd_clf.predict([some_digit])
array([5.])

你可以调用decision_function()方法。不是返回每个样例的一个数值,而是返回10个数值,一个数值对应于一个类。

some_digit_scores = sgd_clf.decision_function([some_digit])
some_digit_scores
array([[-253639.46707377, -425198.63904333, -354213.80127786,
        -229676.13263264, -376404.48500382,   15905.22111141,
        -564592.12430579, -194289.65607053, -748913.30208666,
        -597652.52038338]])

最高的数值对应类别5

np.argmax(some_digit_scores)
5
sgd_clf.classes_
array([0., 1., 2., 3., 4., 5., 6., 7., 8., 9.])

如果你想强制Scikit-Learn使用OvO策略或者OvA策略,你可以使用OneVsOneClassifier类或者OneVsRestClassifier类。创建一个样例,传递一个二分类器给它的构造函数。举例子,下面的代码会创建一个多类分类器,使用OvO策略,基于SGDClassifier

from sklearn.multiclass import OneVsOneClassifier

ovo_clf = OneVsOneClassifier(SGDClassifier(random_state=42))
ovo_clf.fit(X_train, y_train)
ovo_clf.predict([some_digit])
array([5.])

训练一个RandomForestClassifier同样简单:

forest_clf.fit(X_train,y_train)
forest_clf.predict([some_digit])
array([5.])

这次Scikit-Learn没有必要去运行OvO或者OvA, 因为随机森林分类器能够直接将一个样例分到多个类别。你可调用predict_proba(),得到样例对应的类别的概率值的列表:

forest_clf.predict_proba([some_digit])
array([[0. , 0. , 0. , 0. , 0. , 0.9, 0. , 0. , 0.1, 0. ]])

接下来,我们当然想评估一下这些分类器。像以前一样,想便用交叉验证。让我们用cross_val_score来评估SGDClassifier的精度。

cross_val_score(sgd_clf, X_train, y_train,cv = 3, scoring = "accuracy")
array([0.86002799, 0.8760438 , 0.88093214])

我们可以看到这个分类器有86.3%的精度,这个精度还不错,比我们随便乱猜的精度要高出不少(如果我们随机猜,那么精度只有10%)。看起来也并不差,这里可以使输入正则化,得到更高的精度,可以将其精度提高到90%以上。

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_scaled, y_train, cv = 3, scoring="accuracy")
array([0.9080184 , 0.91049552, 0.91043657])

误差分析

分析模型产生的误差,首先,我们可以检查混淆矩阵。需要使用cross_val_predict()做出预测,然后调用confusion_matrix()函数,像以前做的那样

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([[5739,    3,   22,    8,    9,   50,   43,    7,   38,    4],
       [   2, 6451,   50,   23,    6,   46,    5,   14,  133,   12],
       [  58,   38, 5348,   87,   76,   26,   83,   56,  169,   17],
       [  50,   40,  134, 5300,    2,  267,   37,   64,  140,   97],
       [  25,   26,   36,    7, 5356,    9,   54,   32,   83,  214],
       [  68,   37,   34,  179,   74, 4617,  106,   30,  171,  105],
       [  35,   21,   42,    2,   39,   98, 5630,    6,   44,    1],
       [  27,   18,   66,   27,   52,   10,    7, 5793,   17,  248],
       [  58,  150,   68,  140,   16,  156,   51,   29, 5050,  133],
       [  43,   29,   24,   84,  158,   36,    3,  194,   83, 5295]],
      dtype=int64)

这里是一堆数字,使用Matplotlib的matshow()函数,将混淆矩阵以图像的方式呈现,将会更加方便。

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

分类-MNIST-LMLPHP

可以看到,几乎所有的图片都在对角线上,这意味着分类几乎全部正确。现我们只看看其误差的图像

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)
plt.show()

分类-MNIST-LMLPHP

现在可以清楚看出分类器的各类误差,其中行代表实际类别,列代表预测的类别。第8、9列很亮,这说明很多图片被误分成数字8或者数字9。

分析混淆矩阵通常可以提供深刻的见解去改善分类器。回顾这幅图,看样子应该努力改善分类器在数字8和数字9上的表现,和纠正3/5的混淆。举例子,你可以尝试去收集更多的数据,或者你可以构造新的、有助于分类器的特征(新的分类器的特征,我们可以在数据里面加一个新的列———这相当添加了一个新的属性,比如字数8有两个环,数字6有一个,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()

分类-MNIST-LMLPHP

多标签分类

到目前为止,所有的样例都总是被分配到仅一个类(比如我们前面训练的分类,要么输出是1,要么是2,3,...,9,一次只能输出一个类别)。有些情况下,你也许想让你的分类器给一个样例输出多个类别。比如有时候我们想识别某个人脸,想判断它的性别,还有是否为中国人,这就有两个类别了([gender, isChinese])。。这种输出多个二值标签的分类系统被叫做多标签分类系统。

目前不打算深入脸部识别。我们可以先看一个简单点的例子。

from sklearn.neighbors import KNeighborsClassifier
y_train_large = (y_train >=7)
y_train_odd = (y_train % 2 == 1)
y_multilabel = np.c_[y_train_large, y_train_odd]

knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train, y_multilabel)
KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski',
           metric_params=None, n_jobs=None, n_neighbors=5, p=2,
           weights='uniform')

这段代码创造了一个y_multilabel数组,里面包含两个目标标签。第一个标签指出这个数字是否为大数(即是否为7,8,9),第二个标签指示这个数字是否为奇数

knn_clf.predict([some_digit])
array([[False,  True]])

这个预测器预测对,我们输入的数据代表5,5不是一个大数,但是是一个奇数。

# y_train_knn_pred = cross_val_predict(knn_clf, X_train, y_train, cv = 3)
# f1_score(y_train, y_train_knn_pred, average="macro")

多输出分类

我们即将讨论最后一种分类任务,被叫做"多输出-多分类"(或者简称多输出分类)。在这里每一个标签可以是多类别的(比如我们前面所举的例子)

为了说明这点,我们建立一个系统,它可以去除图片当中的噪音。它将一张混有噪音的图片作为输入,期待它输出一张干净的数字图片,用一个像素强度的数组表示,就像 MNIST图片那样。注意到这个分类器的输出是多标签的(一个像素一个标签)和每个标签可以有多个值 (像素强度取值范围从0到255)。所以它是一个多输出分类系统的例子。

我们从MNIST的图版创建训练集和测试集开始,然后给图片的像素强度添加噪声,这里是用NumPy的randint()函数。目标图像是原始图像。

noise = np.random.randint(0, 100, (len(X_train), 784))
X_train_mod = X_train + noise
noise = np.random.randint(0, 100, (len(X_test), 784))
X_test_mod = X_test + noise
y_train_mod = X_train
y_test_mod = X_test
def plot_digit(data):
    image = data.reshape(28, 28)
    plt.imshow(image, cmap = matplotlib.cm.binary,
               interpolation="nearest")
    plt.axis("off")
    
some_index = 5500

plt.subplot(121); plot_digit(X_test_mod[some_index])
plt.subplot(122); plot_digit(y_test_mod[some_index])
# save_fig("noisy_digit_example_plot")
plt.show()

分类-MNIST-LMLPHP

knn_clf.fit(X_train_mod, y_train_mod)
clean_digit = knn_clf.predict([X_test_mod[some_index]])
plot_digit(clean_digit)
# save_fig("cleaned_digit_example_plot")

分类-MNIST-LMLPHP

上面的图片看起来还行,比较接近原图片,去噪的效果还可以。

到这里,分类的知识学得差不多了。

11-04 11:54