初识分类（鸢尾花卉数据集）

注：本文用到的模块如Numpy，Scipy，matplotlib等都可以在 pythonlibs 里下载并使用pip安装，如果没有安装过这里有安装教程。

本文提供的代码基于windows的python2.x，数据和代码都可以在github上打包下载。

机器是否能够识别出图像中的花朵种类？从机器学习的角度来说，我们可以通过以下方式解决这个问题：先让机器学习一下每种花朵的样本数据，然后让它根据这些信息，对未标识出花朵种类的图像进行分类。这个过程就叫做分类（或者叫监督学习），这是一个已经研究了几十年的经典问题。

Iris数据集（Iris dataset，也称鸢尾花卉数据集）是源自20世纪30年代的经典数据集。它是最早应用统计分类的现代示例之一。

数据中包含有不同种类的Iris花朵的数据，这些种类可以通过它们的形态来识别。时至今日，我们已经可以通过基因签名（genomic signature）来识别这些分类。但在20世纪30年代，人们还不确定DNA是不是基因信息的载体。

测量的是每个花朵的以下四个属性：

花萼长度

花萼宽度

花瓣长度

花瓣宽度

一般来说，我们把数据中所有的测量结果都叫做特征。

此外，每个花朵所属的种类都已经标出。现在的问题是：如果看到这种植物的一个新花朵，我们能否通过它的四个特征来预测它的种类？

这就是监督学习或分类问题; 对于给定的带有类别标签的样本，我们要设计出一种规则，实现对其他样本的预测。这跟垃圾邮件分类问题是一样的：根据用户标出的垃圾邮件和非垃圾邮件样本，我们能否判定一个新来的信息是否是垃圾邮件？

在这个时候，Iris数据集可以很好的满足我们的需求。它的规模很小（只包含150个样本，每个样本有4个特征，共有3种分类），很容易可视化地显示出来并进行处理。

我们通过以下代码，对这四个特征进行两两组合（共6种组合）来展示两个特征之间的相关性。

figure1.py

from matplotlib import pyplot as plt
from sklearn.datasets import load_iris
import numpy as np

data = load_iris()
features = data['data']
feature_names = data['feature_names']
target = data['target']

pairs = [(0,1), (0,2), (0,3), (1,2), (1,3), (2,3)]
for i,(p0,p1) in enumerate(pairs):
    plt.subplot(2,3,i+1)
    for t, marker, c in zip(xrange(3), '>ox', 'rgb'):
        plt.scatter(features[target == t, p0],
                    features[target == t, p1],
                    marker=marker,
                    c=c)
    plt.xlabel(feature_names[p0])
    plt.ylabel(feature_names[p1])
    plt.xticks([])
    plt.yticks([])

plt.savefig("01.png")

下图中的每一个子图都画出了所有数据点在其中两个维度上的映射。外面的那一组点（三角形）代表的是山鸢尾花（Iris Setosa），中间的（圆圈）代表的是变色鸢尾花（Iris Versicolor），用x标记的是维吉尼亚鸢尾花（Iris Virginica）。

构建第一个分类模型

如果我们的目标是区分这三种花朵的类型，根据上图不难发现，根据花瓣（petal）长度似乎就可以将山鸢尾花（三角形表示）跟其他两类花朵区分开。我们写一点代码来寻找切分点在哪里，如下所示：

simple_threshold.py

import numpy as np
from sklearn.datasets import load_iris

data=load_iris()
features = data['data']
labels = data['target_names'][data['target']]

plength = features[:,2]
is_setosa = (labels == 'setosa')
print("Maxinum of setosa: {0}." .format(plength[is_setosa].max()))
print("Minimum of others: {0}." .format(plength[~is_setosa].min()))

输出结果如下：

Maxinum of setosa: 1.9.
Minimum of others: 3.0.

现在我们可以通过花瓣长度将山鸢尾花和其他两类花区分开，而且效果非常好，不会出现任何错误。然而，我们却无法立即找到区分维吉尼亚鸢尾花（Iris Virginica）和变色鸢尾花（Iris Vesicolor）的最佳阈值。我们甚至会发现永远都不可能找到完美的划分。不过，我们还是可以尝试一下。下面通过pred（预测结果）和virginica（真实标签）的比较次数的平均值，可以得到正确结果所占的比例，即正确率，我们寻找后两种花分类的最佳阈值。

stump.py

from matplotlib import pyplot as plt
from sklearn.datasets import load_iris

data = load_iris()
features = data['data']
labels = data['target_names'][data['target']]

setosa = (labels == 'setosa')
features = features[~setosa]
labels = labels[~setosa]
virginica = (labels == 'virginica')

best_acc = -1.0
for fi in range(features.shape[1]):
    thresh = features[:,fi].copy()
    for t in thresh:
        pred = (features[:,fi] > t)
        acc = (pred == virginica).mean()
        if acc > best_acc:
            best_acc = acc
            best_fi = fi
            best_t = t
print("Best cut is {0} on feature {1}, which achieves accuracy of {2:.1%}."
      .format(best_t, best_fi,best_acc))

输出结果如下：

Best cut is 1.6 on feature 3, which achieves accuracy of 94.0%.

我们可以通过matplotlib很直观的展示根据feature3（即petal width特征）简单阈值的分类情况。

figure2.py

from matplotlib import pyplot as plt
from sklearn.datasets import  load_iris

data = load_iris()
features = data['data']
feature_names = data['feature_names']
species = data['target_names'][data['target']]

setosa = (species == 'setosa')
features = features[~setosa]
species = species[~setosa]
virginica = (species == 'virginica')

t = 1.75
p0, p1 = 3, 2

area1c = (1.,1, 1)
area2c = (.7, .7, .7)

x0, x1 =[features[:,p0].min()*.9, features[:,p0].max()*1.1]
y0, y1 = [features[:p1].min()*.9, features[:,p1].max()*1.1]
plt.fill_between([t,x1],[y0,y0],[y1,y1],color=area2c)
plt.fill_between([x0,t],[y0,y0],[y1,y1],color=area1c)
plt.plot([t,t],[y0,y1],'k--',lw=2)
plt.plot([t-.1,t-.1],[y0,y1],'k:',lw=2)
plt.scatter(features[virginica,p0],features[virginica,p1],c='b',marker='o')
plt.scatter(features[~virginica,p0],features[~virginica,p1],c='r',marker='x')
plt.ylim(y0,y1)
plt.xlim(x0,x1)
plt.xlabel(feature_names[p0])
plt.ylabel(feature_names[p1])
plt.show()
plt.savefig('02.png')

评估：留存数据和交叉验证

上节中我们讨论了一个简单模型：它在训练集上达到了94%的正确率。然而，这个评估也许过于乐观了。因为我们用这些数据去确定阈值，同时又用它评估这个模型。这在逻辑上犯了循环论证的错误。

我们真正想做的事情是衡量模型对新样本的泛化能力。所以，这里应该用训练中未出现过的数据来评估模型的性能。对此，我们把数据分成两部分，一部分用于训练模型，一部分（从训练集拿出来的数据）用于测试模型效果。为此，我们抽象出三个方法，learn_model, apply_model, accuracy 优化代码。如下所示：

threshold.py

import numpy as np
def learn_model(features, labels):
    best_acc = -1.0
    for fi in range(features.shape[1]):
        thresh = features[:,fi].copy()
        thresh.sort()
        for t in thresh:
            pred = (features[:,fi] > t)
            acc = (pred == labels).mean()
            if acc > best_acc:
                best_acc = acc
                best_fi = fi
                best_t = t
    return best_t, best_fi

def apply_model(features, model):
    t, fi = model
    return features[:,fi] > t

def accuracy(features, labels, model):
    preds = apply_model(features, model)
    return np.mean(preds == labels)

heldout.py

from matplotlib import  pyplot as plt
import numpy as np
from sklearn.datasets import load_iris
from threshold import learn_model, apply_model, accuracy

data = load_iris()
features = data['data']
labels = data['target_names'][data['target']]

setosa = (labels == 'setosa')
features = features[~setosa]
labels = labels[~setosa]
virginica = (labels == 'virginica')

testing = np.tile([True,False], 50)
training = ~testing
model = learn_model(features[training], virginica[training])
train_error = accuracy(features[training],virginica[training],model)
test_error = accuracy(features[testing],virginica[testing],model)

print('''\
Training accuracy was {0:.1%}.
Testing accuracy was {1:.1%} (N = {2}).
'''.format(train_error,test_error,testing.sum()))

输出结果：

Training accuracy was 96.0%.

Testing accuracy was 90.0% (N = 50).

训练数据上的误差叫做训练误差，它对算法效果的估计常常过于乐观。我们总是应该测量和报告测试误差，也就是在未用于训练的样本集合上的误差。

随着模型越来越复杂，这些概念也变得越来越重要。在这个实例中，这两种误差的差距并不是很大。但当使用一个复杂的模型时，训练集上的正确率很可能达到100%，但在测试集上的效果却跟随机猜测差不多。

我们之前采用了从训练集中留存数据的方式，这里有一个潜在的问题，即在训练中只使用了部分数据（在这个例子中，我们使用了一半的数据）。另一方面，如果我们在测试中使用了过少的数据，误差估计将只在很少的一部分样本上进行。在理想情况下，我们希望在训练和测试中都能使用所有的数据。

我们可以通过交叉验证（cross-validation）达到类似的效果。交叉验证的一个极端（但有时很有用）形式叫做去一法（leave-one-out）.从训练集中拿出一个样本，并在缺少这个样本的数据上训练一个模型，然后看模型是否能对这个样本正确分类：

error = 0.0
for ei in range(len(features)):
    training = np.ones(len(features),bool)
    training[ei] = False
    testing = ~training
    model = learn_model(features[training],virginica[training])
    predictions = apply_model(features[testing],virginica[testing],model)
    error += np.sum(predictions != virginica[testing])

交叉验证去一法最主要的问题是工作量大，我们可以通过x折交叉验证。5折和10折是比较常见的。

在生成数据折的时候，你需要谨慎地保持数据分布的平衡。这里不再深入介绍如何来做这个事情，因为机器学习工具包可处理好此事。

构建更复杂的分类器

一个分类模型是由什么组成的？我们可以把它分成三部分。

1.模型结构 在这里我们采用一个阈值在一个特征上进行划分

2.搜索过程 在这里我们尽可能多的尝试所有特征和阈值的组合

3.损失函数 我们通过损失函数来确定哪些可能性不会太差。我们可以用训练误差或者其他方式定义这一点，比如我们想要的正确率。一般来说，人们希望损失函数最小化。

作为选择，我们可以有不同的损失函数。一种可能的情况是，一种错误比另一种的代价更大。在医疗领域，假阴性和假阳性的代价并不相同。假阴性（当检查见过是阴性，但实际上是假的）可能会导致患有严重疾病的患者无法及时得到治疗。假阳性（当检查结果是阳性，即使患者并没有真正患有这种病）可能会导致患者进行更多不必要的诊断和治疗（这仍然是有代价的，例如治疗的副作用）。因此，根据具体情况，选择不同的折中是可以理解的。在一个极端情况下，如果疾病是致命的，而治疗费用很低，副作用也很少，那么我们会希望尽可能使假阴性最少。在垃圾邮件过滤问题中我们面临同样的问题；误删一个正常邮件对用户来说是十分危险的，而让一个垃圾邮件通过，只会带来一个小麻烦。

更复杂的数据集和更复杂的分类器

数据集seeds.tsv

小麦种子的测量数据集。
包括七个特征：面积（A），周长（P），紧密度（C=4πA/（P*P）），谷粒的长度，宽度，偏度系数，谷粒槽长度。
分为三个类别：Canadian，Rosa，Kama

数据集的特征都是从数字图像上自动计算生成的。

图像模式识别是这样实现的：你得到一些数字格式的图像，从中计算出一些相关特征，然后使用一般的分类系统进行分类。后面会有一篇计算机视觉相关的文，从图像中提取特征。

特征和特征工程

这些特征的一个有趣的地方是，紧密度特征实际上并不是一个新的测量值，而是之前两个特征，面积和周长所组成的函数。这个通用领域叫做特征工程（featur engineering）。人们有时认为它没有算法那样富有魅力，但它对系统性能也许有很大影响（一个简单算法在精心选择的特征上的效果比一个漂亮算法在较差的特征上的效果还要好）。

在这里，原作者计算了“紧密度”这个特征，这是一个典型的形状特征（也叫做“圆度”）。如果有两个谷粒，一个的大小是另一个的两倍，但形状一样，那么这个特征将会具有相同的值。然而，对于非常圆的谷粒（特征值接近1）和扁长的谷粒（特征值接近0），它的值会很不相同。

好特征的目标是在重要的地方取不同值，而在不重要的地方不变。例如，紧密度不会随大小而改变，但会随形状变化。在实践中，要同时完美的达到这两个目标可能很困难，但我们希望能逼近这种理想情况。

你需要借助背景知识通过直觉来判断哪些是好特征。幸运的是，在很多问题领域，已经有很多文献介绍了可能用到的特征和特征类型。对于图像，之前提到的所有特征都是很典型的，计算机视觉库可以帮你计算出来。基于文本的问题也是如此，有很多标准的解决方案可以混合搭配。尽管，你通常可以利用特定问题中的领域知识来设计一个特定的特征。

甚至在获取数据之前，你必须决定哪些数据值得收集。然后，你需要将所有的特征交给机器去评估和计算最优分类器。

一个很自然就会想到的问题是，我们能否自动地把好特征选取出来。这个问题叫做特征选择（feature selection）。人们已经提出了很多方法来解决这个问题，但在实践中，极简单的想法可能已经做的很好。在这些小数据集上使用特征选择没有什么意义，但是如果你有几千个特征，那扔掉其中大多数特征将会大大加快后续的流程。

最邻近分类

对于这个数据集，我们首先用简单的阈值对Canadian和非Canadian进行分类：

load.py

import numpy as np
def load_dataset(dataset_name):
    data = []
    labels = []
    with open(r'data\{0}.tsv'.format(dataset_name)) as ifile:
        for line in ifile:
            tokens = line.strip().split('\t')
            data.append([float(tk) for tk in tokens[:-1]])
            labels.append(tokens[-1])
    data = np.array(data)
    labels = np.array(labels)
    return data, labels

seeds_threshold.py

from load import load_dataset
import numpy as np
from threshold import learn_model, apply_model, accuracy

features, labels =load_dataset('seeds')
labels = labels == 'Canadian'

error = 0.0
for fold in range(10):
    training = np.ones(len(features), bool)
    training[fold::10] = 0
    testing = ~training
    model = learn_model(features[training],labels[training])
    test_error =accuracy(features[testing],labels[testing],model)
    error += test_error
error /= 10.0
print('Ten fold cross-validated error was {0:.1%}'.format(error))

输出结果：

Ten fold cross-validated error was 73.8%

可以看到，即使采用之前的方法只能把两个类区分出来，但不能得到很好的结果。因此，让我们介绍一个新的分类器：最邻近分类器。

考虑到每个样本是由它的特征所表示的（用数学语言来讲，它是N维空间中的点），我们可以计算样本之间的距离，而且可以选择不同方法来计算这个距离，例如：

def distance(p0, p1):
    return np.sum((p0-p1)**2)

在分类的时候，我们采用一个简单的规则：对于一个新样本，我们在数据集中寻找最接近它的点（它最近的近邻），并查看它的标签：

def nn_classify(training_set, training_labels, new_example):
    dists = np.array([distance(t, new_example) for t in training_set])
    nearest = dists.argmin()
    return training_labels[nearest]

在这种情况下，我们的模型不使用任何训练数据及标签，就能在分类阶段计算出所有的结果。一个更好的实现方法是，在学习阶段对这些数据做索引，从而加速分类的计算，但这个实现是一个很复杂的算法。

现在，我们注意到这个模型在训练数据上表现的很完美。在每个数据点上，它的最近近邻就是它自己，所以它的标签一定完美匹配（除非两个样本具有相同的特征但标签不同，而这是可能发生的）。因此，采用交叉验证来进行测试是很必要的。

我们在这个数据集上应用这个算法，并进行10折交叉验证，可以得到89.5%的正确率，正如之前几节中讨论的，交叉验证正确率会低于训练正确率，但这是对模型性能更可靠的估计。

我们现在来看一下判别边界。对此，我们必须简化一下，只考虑两个维度（这样我们就可以把它画在纸上了）。

在前图中，Canadian样本是用菱形表示的，Kama种子是用圆圈表示的，Rosa种子是用三角形表示的。它们的区域分别用白色、黑色和灰色表示。你可能会对为何这些区域都是水平方向感到疑惑，觉得这非常古怪。这里的问题在于，x轴（面积）的值在10到22之间，而y轴（紧密度）在0.75到1之间。这意味着，x值的一小点改变实际上比y值的一小点变化大得多。所以，在用之前的函数计算距离的时候，我们多半只把x轴考虑进去了。

如果你了解一些物理学背景知识，你可能已经注意到我们已经把长度、面积和无量纲的量加了起来，把各种单位混合在一起了（我们从不会在物理系统中这样做）。我们需要把所有特征都归一化到一个公共尺度上。这个问题有很多解决方法，一个简单的方法是把它们归一到Z值（Z-score）。Z值表示的是特征值离它的平均值有多远，它用标准方差的数量来计算。它可以归结到以下这一对简单的操作上：

features -= features.mean(axis = 0)
features /= features.atd(axis = 0)

转换为Z值之后，0就是平均值，正数是高于平均值的值，负数是低于平均值的值，它独立于原始值。

现在每个特征都采用了同样的单位（严格来说，每个特征现在都是无量纲的，它没有单位），因此我们可以更放心地混合不同维度的特征。事实上，如果现在运行最邻近分类器，我们可以得到94.3%的正确率！

knn.py

import numpy as np
def learn_model(k, features, labels):
    return k, features.copy(),labels.copy()

def plurality(xs):
    from collections import defaultdict
    counts = defaultdict(int)
    for x in xs:
        counts[x] += 1
    maxv = max(counts.values())
    for k,v in counts.items():
        if v == maxv:
            return k

def apply_model(features, model):
    k, train_feats, labels = model
    results = []
    for f in features:
        label_dist = []
        for t,ell in zip(train_feats, labels):
            label_dist.append( (np.linalg.norm(f-t), ell) )
        label_dist.sort(key=lambda d_ell: d_ell[0])
        label_dist = label_dist[:k]
        results.append(plurality([ell for _,ell in label_dist]))
    return np.array(results)

def accuracy(features, labels, model):
    preds = apply_model(features, model)
    return np.mean(preds == labels)

seeds_knn.py

from load import load_dataset
import numpy as np
from knn import learn_model, apply_model, accuracy

features,labels = load_dataset('seeds')

def cross_validate(features, labels):
    error = 0.0
    for fold in range(10):
        training = np.ones(len(features), bool)
        training[fold::10] = 0
        testing = ~training
        model = learn_model(1, features[training], labels[training])
        test_error = accuracy(features[testing], labels[testing], model)
        error += test_error

    return error/ 10.0

error = cross_validate(features, labels)
print('Ten fold cross-validated error was {0:.1%}.'.format(error))

features -= features.mean(0)
features /= features.std(0)
error = cross_validate(features, labels)
print('Ten fold cross-validated error after z-scoring was {0:.1%}.'.format(error))

输出结果：

Ten fold cross-validated error was 89.5%.
Ten fold cross-validated error after z-scoring was 94.3%.

再次看看这两个维度的决策空间，它如下图所示：

分界面变得非常复杂，而且在两个维度之间有一个相互交叉。如果使用数据全集，那么所有这些都发生在七维空间中，它非常难以可视化地显示出来。但是道理是一样的：之前是少数维度占据统治地位，而现在所有维度都具有相同的重要性。

最邻近分类器虽然简单，但有时它的效果已经足够好。我们可以把它泛化成一个k邻近分类器，不仅考虑最邻近的点，还要考虑前k个最临近点。所有这k个邻近点通过投票方式来选择标签。k一般是一个小数字，比如5，但它也可以更大，特别是当数据规模非常大的时候。

附灰度图和彩色图的代码如下：

figure4_5.py

COLOUR_FIGURE = True

from matplotlib import pyplot as plt
from matplotlib.colors import ListedColormap
from load import load_dataset
import numpy as np
from knn import learn_model,apply_model,accuracy

feature_names = [
    'area',
    'perimeter',
    'compactness',
    'length of kernel',
    'width of kernel',
    'asymmetry coefficien',
    'length of kernel groove'
]

def train_plot(features, labels):
    y0, y1 = features[:,2].min()*.9, features[:,2].max()*1.1
    x0, x1 = features[:,0].min()*.9, features[:,0].max()*1.1
    X = np.linspace(x0,x1,100)
    Y = np.linspace(y0,y1,100)
    X,Y = np.meshgrid(X,Y)

    model = learn_model(1, features[:,(0,2)], labels)
    C = apply_model(np.array(np.vstack([X.ravel(), Y.ravel()]).T), model).reshape(X.shape)
    if COLOUR_FIGURE:
        cmap = ListedColormap([(1,.6,.6),(.6,1.,.6),(.6,.6,1.)])
    else:
        cmap = ListedColormap([(1.,1.,1.),(.2,.2,.2),(.6,.6,.6)])
    plt.xlim(x0,x1)
    plt.ylim(y0,y1)
    plt.xlabel(feature_names[0])
    plt.ylabel(feature_names[2])
    plt.pcolormesh(X,Y,C, cmap=cmap)
    if COLOUR_FIGURE:
        cmap = ListedColormap([(1.,.0,.0),(.0,1.,.0),(.0,.0,1.)])
        plt.scatter(features[:,0],features[:,2], c =labels,cmap=cmap)
    else:
        for lab,ma in zip(range(3),"Do^"):
            plt.plot(features[labels == lab, 0],features[lab==lab,2],ma,c=(1.,1.,1.))

features, labels = load_dataset('seeds')
names = sorted(set(labels))
labels = np.array([names.index(ell) for ell in labels])

train_plot(features, labels)
plt.show()
plt.savefig('04.png')

features -= features.mean(0)
features /= features.std(0)
train_plot(features,labels)
plt.show()
plt.savefig('05.png')

二分类和多分类

我们看到的第一个分类器，阈值分类器，是一个简单的二类分类器（由于数据点不是高于阈值就是低于阈值，所以分类结果不是第一个类，就是第二个类）。我们用的第二个分类器，最邻近分类器，天然就是一个多类分类器（它的输出可以是多个类别中的一个）。

构建一个二分类方法通常要比构建一个解决多分类问题的方法更加简单。然而，我们可以将多分类问题细化成一系列二分决策。这就是之前我们在Iris数据集上顺带做出的；我们观察到，将原始类别中的一个类别分离出来很容易，我们需要专注于另外两个类别的区分，而这些可以退化成几个二分类决策。

总结

我们用简单示例介绍了很多一般性概念。然而这只是一个规模很小的问题。它的优点在于能让我们把它画出来，看到我们具体在做什么。当我们换一个维度高、样本多的问题时，这一优点就不见了。

分类意味着对样本进行归纳，从而构建出一个模型（这是一个能够自动对新的、未分类的数据进行分类的规则）。这是机器学习的一个基础工具。

我们还学习到，对于模型效果，训练误差是一个有误导性的、过于乐观的估计。相反，我们必须使用未用于训练的测试数据来评估效果。为了在测试中不浪费过多的样本，交叉验证计划可以帮我们兼得两者的优势（以更多的计算作为代价）。

我们还探究了一下特征工程问题。特征并不是天生就为你预备的，但选择和设计特征却是设计机器学习流程的一个组成部分。事实上，这通常是一个能够获得最大正确率提升的地方，这是因为更好的特征数据往往可以击败更漂亮的方法。

接下来，我们来看下当数据中没有预设的类别信息时应当如何处理。