原文：
annas-archive.org/md5/7039549ff2b32d189f96a3420dc66360译者：飞龙
第三章：降维
在本章中，我们将涵盖以下食谱：
-
使用 PCA 进行降维
-
使用因子分析进行分解
-
使用核 PCA 进行非线性降维
-
使用截断 SVD 进行降维
-
使用分解进行分类与 DictionaryLearning
-
使用流形进行降维——t-SNE
-
测试通过管道减少维度的方法
介绍
在本章中，我们将减少输入到机器学习模型中的特征或输入数量。这是一个非常重要的操作，因为有时数据集有很多输入列，减少列数可以创建更简单的模型，减少计算能力的需求以进行预测。
本节使用的主要模型是主成分分析（PCA）。由于 PCA 的解释方差，您无需知道可以将数据集减少到多少特征。一个性能相似的模型是截断奇异值分解（truncated SVD）。最好首先选择一个线性模型，允许您知道可以将数据集减少到多少列，例如 PCA 或截断 SVD。
在本章后面，查看现代方法t 分布随机邻居嵌入（t-SNE），它使特征在低维度中更容易可视化。在最后一个食谱中，您可以检查一个复杂的管道和网格搜索，找到由降维与多个支持向量机组成的最佳复合估计器。
使用 PCA 进行降维
现在是时候将数学提升到一个新层次了！PCA 是本书中讨论的第一个相对高级的技术。到目前为止，所有的内容都只是简单的统计学，而 PCA 将统计学与线性代数结合起来，产生一个预处理步骤，有助于减少维度，这是简化模型的敌人。
准备工作
PCA 是 scikit-learn 中分解模块的一个成员。还有几种其他的分解方法，稍后将在本食谱中介绍。我们将使用鸢尾花数据集，但如果你使用自己的数据会更好：
from sklearn import datasets
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
%matplotlib inline
iris = datasets.load_iris()
iris_X = iris.data
y = iris.target
如何实现…
- 导入
decomposition模块：
from sklearn import decomposition
- 实例化一个默认的 PCA 对象：
pca = decomposition.PCA()
pca
PCA(copy=True, iterated_power='auto', n_components=None, random_state=None,
svd_solver='auto', tol=0.0, whiten=False)
- 与 scikit-learn 中的其他对象相比，PCA 对象所需的参数相对较少。现在 PCA 对象（一个 PCA 实例）已经创建，只需通过调用
fit_transform方法来转换数据，iris_X作为参数：
iris_pca = pca.fit_transform(iris_X)
iris_pca[:5]
array([[ -2.68420713e+00, 3.26607315e-01, -2.15118370e-02,
1.00615724e-03],
[ -2.71539062e+00, -1.69556848e-01, -2.03521425e-01,
9.96024240e-02],
[ -2.88981954e+00, -1.37345610e-01, 2.47092410e-02,
1.93045428e-02],
[ -2.74643720e+00, -3.11124316e-01, 3.76719753e-02,
-7.59552741e-02],
[ -2.72859298e+00, 3.33924564e-01, 9.62296998e-02,
-6.31287327e-02]])
- 现在 PCA 对象已经拟合完成，我们可以看到它在解释方差方面的效果如何（将在接下来的*工作原理…*部分中进行说明）：
pca.explained_variance_ratio_
array([ 0.92461621, 0.05301557, 0.01718514, 0.00518309])
它是如何工作的…
PCA 有一个通用的数学定义，并在数据分析中有特定的应用案例。PCA 找到一组正交方向，这些方向表示原始数据矩阵。
通常，PCA 通过将原始数据集映射到一个新空间来工作，其中矩阵的每个新列向量都是正交的。从数据分析的角度来看，PCA 将数据的协方差矩阵转换为可以解释方差某些百分比的列向量。例如，使用鸢尾花数据集，92.5%的整体方差可以通过第一个分量来解释。
这非常有用，因为维度问题在数据分析中很常见。许多应用于高维数据集的算法会在初始训练时出现过拟合，从而失去对测试集的泛化能力。如果数据的绝大部分结构可以通过更少的维度忠实地表示，那么这通常被认为是值得的权衡：
pca = decomposition.PCA(n_components=2)
iris_X_prime = pca.fit_transform(iris_X)
iris_X_prime.shape
(150L, 2L)
现在我们的数据矩阵是 150 x 2，而不是 150 x 4。即便在减少维度为二之后，类别的可分性依然保持。我们可以看到这两维所表示的方差有多少：
pca.explained_variance_ratio_.sum()
0.97763177502480336
为了可视化 PCA 的效果，让我们绘制鸢尾花数据集的前两维，并展示 PCA 变换前后的对比图：
fig = plt.figure(figsize=(20,7))
ax = fig.add_subplot(121)
ax.scatter(iris_X[:,0],iris_X[:,1],c=y,s=40)
ax.set_title('Before PCA')
ax2 = fig.add_subplot(122)
ax2.scatter(iris_X_prime[:,0],iris_X_prime[:,1],c=y,s=40)
ax2.set_title('After PCA')
PCA对象也可以在一开始就考虑解释方差的数量。例如，如果我们希望至少解释 98%的方差，那么PCA对象将按如下方式创建：
pca = decomposition.PCA(n_components=.98)
iris_X_prime = pca.fit(iris_X).transform(iris_X)
pca.explained_variance_ratio_.sum()
0.99481691454981014
由于我们希望解释的方差稍微多于两个分量示例，因此包含了第三个分量。
即使最终数据的维度是二维或三维，这两三列也包含了所有四个原始列的信息。
还有更多…
建议在使用 PCA 之前进行缩放。操作步骤如下：
from sklearn import preprocessing
iris_X_scaled = preprocessing.scale(iris_X)
pca = decomposition.PCA(n_components=2)
iris_X_scaled = pca.fit_transform(iris_X_scaled)
这导致了如下图：
fig = plt.figure(figsize=(20,7))
ax = fig.add_subplot(121)
ax.scatter(iris_X_prime[:,0],iris_X_prime[:,1],c=y,s=40)
ax.set_title('Regular PCA')
ax2 = fig.add_subplot(122)
ax2.scatter(iris_X_scaled[:,0],iris_X_scaled[:,1],c=y,s=40)
ax2.set_title('Scaling followed by PCA')
看起来有点差。无论如何，如果你考虑使用 PCA，始终应该考虑使用缩放后的 PCA。最好能通过管道按如下方式进行缩放：
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
pipe = Pipeline([('scaler', StandardScaler()), ('pca',decomposition.PCA(n_components=2))])
iris_X_scaled = pipe.fit_transform(iris_X)
使用管道可以防止错误，并减少复杂代码的调试工作量。
使用因子分析进行分解
因子分析是我们可以用来减少维度的另一种技术。然而，因子分析有前提假设，而 PCA 没有。基本假设是存在一些隐含特征，它们决定了数据集的特征。
这个方法将提取样本中的显式特征，以期理解独立变量和因变量。
准备就绪
为了比较 PCA 和因子分析，我们再次使用鸢尾花数据集，但我们首先需要加载FactorAnalysis类：
from sklearn import datasets
iris = datasets.load_iris()
iris_X = iris.data
from sklearn.decomposition import FactorAnalysis
如何做到这一点…
从编程角度来看，因子分析与 PCA 没有太大区别：
fa = FactorAnalysis(n_components=2)
iris_two_dim = fa.fit_transform(iris.data)
iris_two_dim[:5]
array([[-1.33125848, -0.55846779],
[-1.33914102, 0.00509715],
[-1.40258715, 0.307983 ],
[-1.29839497, 0.71854288],
[-1.33587575, -0.36533259]])
将以下图与上一节中的图进行比较：
由于因子分析是一个概率变换，我们可以检查不同的方面，例如模型下观测值的对数似然性，甚至更好的是，比较不同模型的对数似然性。
因子分析并非没有缺点。原因在于你不是在拟合一个模型来预测结果，而是将模型作为准备步骤来拟合。这并不是什么坏事，但当你在训练实际模型时，错误会被累积。
工作原理…
因子分析与之前讲解的 PCA 相似，然而它们之间有一个重要的区别。PCA 是数据的线性变换，转到一个不同的空间，在这个空间里，第一个主成分解释了数据的方差，而每个后续主成分与第一个主成分正交。
例如，你可以将 PCA 想象成将一个N维的数据集降到某个M维的空间，其中M < N。
另一方面，因子分析假设只有M个重要特征，这些特征的线性组合（加噪声）创建了N维度的数据集。换句话说，你不是在对结果变量做回归，而是在对特征做回归，以确定数据集的潜在因子。
此外，一个大缺点是你不知道可以将数据降到多少列。PCA 会提供解释方差的指标，以指导你完成这一过程。
使用核 PCA 进行非线性降维
统计学中的大多数技术本质上是线性的，因此为了捕捉非线性，我们可能需要应用一些变换。PCA 当然是线性变换。在本步骤中，我们将看一下应用非线性变换，然后应用 PCA 进行降维。
准备工作
如果数据总是线性可分，生活会变得非常简单，但不幸的是，数据并非总是如此。核主成分分析（Kernel PCA）可以帮助解决这个问题。数据首先通过核函数进行处理，将数据投影到一个不同的空间；然后，执行 PCA。
为了熟悉核函数，一个很好的练习是思考如何生成能够被核 PCA 中可用的核函数分离的数据。在这里，我们将使用余弦核来完成。这个步骤的理论内容会比之前的更多。
在开始之前，加载鸢尾花数据集：
from sklearn import datasets, decomposition
iris = datasets.load_iris()
iris_X = iris.data
如何操作…
余弦核通过比较在特征空间中表示的两个样本之间的角度来工作。当向量的大小扰动了用于比较样本的典型距离度量时，它就显得很有用。提醒一下，两个向量之间的余弦通过以下公式给出：
这意味着A和B之间的余弦是这两个向量的点积，通过各自的范数的乘积进行归一化。向量A和B的大小对这个计算没有影响。
所以，让我们回到鸢尾花数据集，使用它进行视觉对比：
kpca = decomposition.KernelPCA(kernel='cosine', n_components=2)
iris_X_prime = kpca.fit_transform(iris_X)
然后，展示结果：
结果看起来稍微好一些，尽管我们需要进行测量才能确认。
它是如何工作的…
除了余弦核，还有几种不同的核可供选择。你甚至可以编写自己的核函数。可用的核如下：
-
多项式（Poly）
-
RBF（径向基函数）
-
Sigmoid
-
余弦
-
预计算
还有一些选项依赖于核的选择。例如，degree 参数将为多项式核、RBF 核和 Sigmoid 核指定度数；此外，gamma 将影响 RBF 或多项式核。
SVM 中的食谱将更详细地介绍 RBF 核函数。
核方法非常适合创建可分性，但如果使用不当，也可能导致过拟合。确保适当进行训练和测试。
幸运的是，现有的核是平滑的、连续的且可微的函数。它们不会像回归树那样产生锯齿状的边缘。
使用截断 SVD 来减少维度
截断 SVD 是一种矩阵分解技术，将矩阵M分解为三个矩阵U、Σ和V。这与 PCA 非常相似，不同之处在于，SVD 的分解是在数据矩阵上进行的，而 PCA 的分解则是在协方差矩阵上进行的。通常，SVD 在幕后用于找到矩阵的主成分。
准备工作
截断 SVD 不同于常规 SVD，它产生的分解结果列数等于指定的截断数。例如，对于一个n x n的矩阵，SVD 将产生n列的矩阵，而截断 SVD 将产生指定列数的矩阵。通过这种方式，维度得以减少。这里我们将再次使用鸢尾花数据集，供你与 PCA 结果进行比较：
from sklearn.datasets import load_iris
iris = load_iris()
iris_X = iris.data
y = iris.target
如何操作…
这个对象与我们使用的其他对象形式相同。
首先，我们将导入所需的对象，然后拟合模型并检查结果：
from sklearn.decomposition import TruncatedSVD
svd = TruncatedSVD(2)
iris_transformed = svd.fit_transform(iris_X)
然后，展示结果：
结果看起来相当不错。与 PCA 一样，有explained_variance_ratio_的解释方差：
svd.explained_variance_ratio_
array([ 0.53028106, 0.44685765])
它是如何工作的…
现在我们已经了解了在 scikit-learn 中如何执行，让我们看看如何只使用 SciPy，并在这个过程中学到一些东西。
首先，我们需要使用 SciPy 的linalg执行 SVD：
from scipy.linalg import svd
import numpy as np
D = np.array([[1, 2], [1, 3], [1, 4]])
D
array([[1, 2],
[1, 3],
[1, 4]])
U, S, V = svd(D, full_matrices=False)
U.shape, S.shape, V.shape
((3L, 2L), (2L,), (2L, 2L))
我们可以重构原始矩阵D，以确认U、S和V作为分解：
np.dot(U.dot(np.diag(S)), V)
array([[1, 2],
[1, 3],
[1, 4]])
实际上，由截断 SVD 返回的矩阵是U和S矩阵的点积。如果我们想模拟截断，我们将丢弃最小的奇异值及其对应的U列向量。因此，如果我们想要一个单一的组件，我们将执行以下操作：
new_S = S[0]
new_U = U[:, 0]
new_U.dot(new_S)
array([-2.20719466, -3.16170819, -4.11622173])
一般来说，如果我们想要截断到某个维度，例如t，我们会丢弃N - t个奇异值。
还有更多内容…
截断 SVD 有一些杂项内容值得注意，特别是在方法方面。
符号翻转
截断 SVD 有个陷阱。根据随机数生成器的状态，连续应用截断 SVD 可能会翻转输出的符号。为了避免这种情况，建议只进行一次截断 SVD 拟合，然后从那时起使用变换。这是管道方法的另一个好理由！
为了实现这一点，请执行以下操作：
tsvd = TruncatedSVD(2)
tsvd.fit(iris_X)
tsvd.transform(iris_X)
稀疏矩阵
截断 SVD 相对于 PCA 的一个优势是，截断 SVD 可以作用于稀疏矩阵，而 PCA 不能。这是因为 PCA 必须计算协方差矩阵，而这需要在整个矩阵上进行操作。
使用分解方法通过 DictionaryLearning 进行分类
在这个示例中，我们将展示如何使用分解方法进行分类。DictionaryLearning试图将数据集转化为稀疏表示。
准备工作
使用DictionaryLearning，其思路是特征是结果数据集的基础。加载 iris 数据集：
from sklearn.datasets import load_iris
iris = load_iris()
iris_X = iris.data
y = iris.target
此外，通过取iris_X和y的每隔一个元素来创建训练集。剩下的元素用来进行测试：
X_train = iris_X[::2]
X_test = iris_X[1::2]
y_train = y[::2]
y_test = y[1::2]
如何实现…
- 导入
DictionaryLearning：
from sklearn.decomposition import DictionaryLearning
- 使用三个组件来表示三种鸢尾花：
dl = DictionaryLearning(3)
- 变换每隔一个数据点，这样我们就可以在学习器训练后，在结果数据点上测试分类器：
transformed = dl.fit_transform(X_train)
transformed[:5]
array([[ 0\. , 6.34476574, 0\. ],
[ 0\. , 5.83576461, 0\. ],
[ 0\. , 6.32038375, 0\. ],
[ 0\. , 5.89318572, 0\. ],
[ 0\. , 5.45222715, 0\. ]])
- 现在通过简单输入以下命令来测试变换：
test_transform = dl.transform(X_test)
我们可以可视化输出。注意每个值是如何在x、y或z轴上定位的，并且与其他值和零一起显示；这被称为稀疏性：
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(14,7))
ax = fig.add_subplot(121, projection='3d')
ax.scatter(transformed[:,0],transformed[:,1],transformed[:,2],c=y_train,marker = '^')
ax.set_title("Training Set")
ax2 = fig.add_subplot(122, projection='3d')
ax2.scatter(test_transform[:,0],test_transform[:,1],test_transform[:,2],c=y_test,marker = '^')
ax2.set_title("Testing Set")
如果你仔细观察，你会发现有一个训练错误。某个类别被误分类了。不过，错误只发生一次并不算大问题。分类中也有错误。如果你记得其他一些可视化，红色和绿色类别经常出现在彼此接近的位置。
它是如何工作的…
DictionaryLearning的背景包括信号处理和神经学。其思路是，在任何给定时刻，只有少数特征是活跃的。因此，DictionaryLearning试图在大多数特征应该为零的约束下，找到数据的合适表示。
使用流形进行维度降维 – t-SNE
准备工作
这是一个简短而实用的示例。
如果你阅读了本章的其余部分，你会发现我们已经在使用 iris 数据集进行很多维度降维。我们继续这种模式进行额外的简便比较。加载 iris 数据集：
from sklearn.datasets import load_iris
iris = load_iris()
iris_X = iris.data
y = iris.target
加载PCA以及manifold模块中的一些类：
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE, MDS, Isomap
#Load visualization library
import matplotlib.pyplot as plt
%matplotlib inline
如何做…
- 对
iris_X应用所有变换。其中一个变换是 t-SNE：
iris_pca = PCA(n_components = 2).fit_transform(iris_X)
iris_tsne = TSNE(learning_rate=200).fit_transform(iris_X)
iris_MDS = MDS(n_components = 2).fit_transform(iris_X)
iris_ISO = Isomap(n_components = 2).fit_transform(iris_X)
- 绘制结果：
plt.figure(figsize=(20, 10))
plt.subplot(221)
plt.title('PCA')
plt.scatter(iris_pca [:, 0], iris_pca [:, 1], c=y)
plt.subplot(222)
plt.scatter(iris_tsne[:, 0], iris_tsne[:, 1], c=y)
plt.title('TSNE')
plt.subplot(223)
plt.scatter(iris_MDS[:, 0], iris_MDS[:, 1], c=y)
plt.title('MDS')
plt.subplot(224)
plt.scatter(iris_ISO[:, 0], iris_ISO[:, 1], c=y)
plt.title('ISO')
t-SNE 算法最近很受欢迎，但它需要大量的计算时间和算力。ISO 生成了一个有趣的图形。
此外，在数据的维度非常高（超过 50 列）的情况下，scikit-learn 文档建议在 t-SNE 之前进行 PCA 或截断 SVD。鸢尾花数据集较小，但我们可以编写语法，在 PCA 之后执行 t-SNE：
iris_pca_then_tsne = TSNE(learning_rate=200).fit_transform(iris_pca)
plt.figure(figsize=(10, 7))
plt.scatter(iris_pca_then_tsne[:, 0], iris_pca_then_tsne[:, 1], c=y)
plt.title("PCA followed by TSNE")
它是如何工作的…
在数学中，流形是一个在每个点局部欧几里得的空间，但它嵌入在更高维的空间中。例如，球体的外表面在三维空间中是一个二维流形。当我们在地球表面行走时，我们倾向于感知地面的二维平面，而不是整个三维空间。我们使用二维地图导航，而不是更高维的地图。
scikit-learn 中的manifold模块对于理解高维空间中的二维或三维空间非常有用。该模块中的算法收集关于某个点周围局部结构的信息，并试图保持这一结构。什么是一个点的邻居？一个点的邻居有多远？
例如，Isomap 算法试图在算法中保持所有点之间的测地距离，从最近邻搜索开始，接着是图搜索，再到部分特征值分解。该算法的目的是保持距离和流形的局部几何结构。多维尺度法（MDS）算法同样尊重距离。
t-SNE 将数据集中点对之间的欧几里得距离转化为概率。在每个点周围都有一个以该点为中心的高斯分布，概率分布表示任何其他点成为邻居的概率。相距很远的点，成为邻居的概率很低。在这里，我们将点的位置转化为距离，再转化为概率。t-SNE 通过利用两点成为邻居的概率，能够很好地保持局部结构。
从非常一般的角度来看，流形方法通过检查每个点的邻居来开始，这些邻居表示流形的局部结构，并试图以不同的方式保持该局部结构。这类似于你在邻里或街区上走动，构建你周围局部结构的二维地图，并专注于二维而不是三维。
使用管道进行降维的测试方法
在这里，我们将看到由降维和支持向量机组成的不同估算器的表现。
准备工作
加载鸢尾花数据集和一些降维库。对于这个特定的步骤来说，这是一个重要的步骤：
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.model_selection import GridSearchCV
from sklearn.pipeline import Pipeline
from sklearn.svm import SVC, LinearSVC
from sklearn.decomposition import PCA, NMF, TruncatedSVD
from sklearn.manifold import Isomap
%matplotlib inline
如何实现…
-
实例化一个包含两大部分的管道对象：
-
一个用于降维的对象
-
具有 predict 方法的估算器
-
pipe = Pipeline([
('reduce_dim', PCA()),
('classify', SVC())
])
- 请注意以下代码中，Isomap 来自
manifold模块，并且非负矩阵分解（NMF）算法利用 SVD 将矩阵分解为非负因子，其主要目的是与其他算法进行性能比较，但在自然语言处理（NLP）中非常有用，因为矩阵分解的结果不能为负。现在输入以下参数网格：
param_grid = [
{
'reduce_dim': [PCA(), NMF(),Isomap(),TruncatedSVD()],
'reduce_dim__n_components': [2, 3],
'classify' : [SVC(), LinearSVC()],
'classify__C': [1, 10, 100, 1000]
},
]
这个参数网格将允许 scikit-learn 通过一些降维技术与两种 SVM 类型结合：线性 SVC 和用于分类的 SVC。
- 现在运行一次网格搜索：
grid = GridSearchCV(pipe, cv=3, n_jobs=-1, param_grid=param_grid)
iris = load_iris()
grid.fit(iris.data, iris.target)
- 现在查看最佳参数，以确定最佳模型。使用 PCA 和 SVC 是最佳模型：
grid.best_params_
{'classify': SVC(C=10, cache_size=200, class_weight=None, coef0=0.0,
decision_function_shape=None, degree=3, gamma='auto', kernel='rbf',
max_iter=-1, probability=False, random_state=None, shrinking=True,
tol=0.001, verbose=False),
'classify__C': 10,
'reduce_dim': PCA(copy=True, iterated_power='auto', n_components=3, random_state=None,
svd_solver='auto', tol=0.0, whiten=False),
'reduce_dim__n_components': 3}
grid.best_score_
0.97999999999999998
- 如果你想创建一个结果的数据框，使用以下命令：
import pandas as pd
results_df = pd.DataFrame(grid.cv_results_)
- 最后，你可以通过
grid.predict(X_test)方法对未见实例进行预测，X_test是测试集。我们将在后续章节中进行几次网格搜索。
它是如何工作的……
网格搜索进行交叉验证，以确定最佳得分。在这种情况下，所有数据都用于三折交叉验证。对于本书的其余部分，我们将保留一些数据用于测试，以确保模型不会出现异常行为。
关于你刚才看到的管道的最后一点：sklearn.decomposition方法将用于管道中的第一个降维步骤，但并非所有流形方法都设计为适用于管道。
第四章：使用 scikit-learn 进行线性模型
本章包含以下几个步骤：
-
通过数据拟合一条直线
-
使用机器学习通过数据拟合一条直线
-
评估线性回归模型
-
使用岭回归克服线性回归的不足
-
优化岭回归参数
-
使用稀疏性正则化模型
-
采用更基础的方法使用 LARS 进行正则化
介绍
我推测我们天生能很好地感知线性函数。它们非常容易可视化、解释和说明。线性回归非常古老，可能是第一个统计模型。
在本章中，我们将采用机器学习方法进行线性回归。
请注意，这一章节与降维和 PCA 章节类似，涉及使用线性模型选择最佳特征。即使您决定不使用线性模型进行预测回归，也可以选择最有效的特征。
还要注意，线性模型提供了许多机器学习算法使用背后的直觉。例如，RBF 核 SVM 具有平滑边界，从近距离看，它们看起来像一条直线。因此，如果你记住你的线性模型直觉，解释 SVM 就会变得容易。
通过数据拟合一条直线
现在我们将从线性回归的基础建模开始。传统线性回归是第一个，因此可能是最基本的模型——数据的一条直线。
直观地说，对于大多数人来说很熟悉：一个输入变量的变化会按比例改变输出变量。许多人在学校、报纸的数据图表、工作中的演示中都见过它，因此你可以很容易向同事和投资者解释它。
准备工作
波士顿数据集非常适合用于回归分析。波士顿数据集包括波士顿多个地区的房屋中位数价格。它还包括可能影响房价的其他因素，例如犯罪率。首先，导入datasets模块，然后我们可以加载数据集：
from sklearn import datasets boston = datasets.load_boston()
如何做…
实际上，在 scikit-learn 中使用线性回归非常简单。线性回归的 API 基本上与你从前一章已经熟悉的 API 相同。
- 首先，导入
LinearRegression对象并创建一个对象：
from sklearn.linear_model import LinearRegression lr = LinearRegression()
- 现在，只需将独立变量和因变量传递给
LinearRegression的fit方法即可开始：
lr.fit(boston.data, boston.target)
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
- 现在，要获得预测结果，请执行以下操作：
predictions = lr.predict(boston.data)
- 你已经获得了线性回归生成的预测结果。现在，进一步探索
LinearRegression类。查看残差，即实际目标集和预测目标集之间的差异：
import numpy as np import pandas as pd import matplotlib.pyplot as plt #within an Ipython notebook:
%matplotlib inline
pd.Series(boston.target - predictions).hist(bins=50)
表达特征系数及其名称的常见模式是zip(boston.feature_names, lr.coef_)。
- 通过键入
lr.coef_找到线性回归的系数：
lr.coef_
array([ -1.07170557e-01, 4.63952195e-02, 2.08602395e-02,
2.68856140e+00, -1.77957587e+01, 3.80475246e+00,
7.51061703e-04, -1.47575880e+00, 3.05655038e-01,
-1.23293463e-02, -9.53463555e-01, 9.39251272e-03,
-5.25466633e-01])
所以，回到数据上，我们可以看到哪些因素与结果呈负相关，哪些因素呈正相关。例如，正如预期的那样，城镇的人均犯罪率的增加与波士顿的房价呈负相关。人均犯罪率是回归分析中的第一个系数。
- 你还可以查看截距，即当所有输入变量为零时目标的预测值：
lr.intercept_
36.491103280361955
- 如果你忘记了系数或截距属性的名称，可以输入
dir(lr)：
[... #partial output due to length
'coef_',
'copy_X',
'decision_function',
'fit',
'fit_intercept',
'get_params',
'intercept_',
'n_jobs',
'normalize',
'predict',
'rank_',
'residues_',
'score',
'set_params',
'singular_']
对于许多 scikit-learn 预测器，参数名以单词加一个下划线结尾，如 coef_ 或 intercept_，这些参数特别值得关注。使用 dir 命令是检查 scikit 预测器实现中可用项的好方法。
它是如何工作的…
线性回归的基本思想是找到满足 y = Xv 的 v 系数集合，其中 X 是数据矩阵。对于给定的 X 值，通常不太可能找到一个完全满足方程的系数集合；如果存在不精确的规格或测量误差，误差项将被添加进来。
因此，方程变为 y = Xv + e，其中 e 被假定为正态分布，并且与 X 的值独立。从几何角度来看，这意味着误差项与 X 垂直。虽然这超出了本书的范围，但你可能想亲自证明 E(Xv) = 0。
还有更多内容…
LinearRegression 对象可以自动对输入进行标准化（或缩放）：
lr2 = LinearRegression(normalize=True)
lr2.fit(boston.data, boston.target)
LinearRegression(copy_X=True, fit_intercept=True, normalize=True)
predictions2 = lr2.predict(boston.data)
使用机器学习拟合数据线
使用机器学习的线性回归涉及在未见过的数据上测试线性回归算法。在这里，我们将进行 10 折交叉验证：
-
将数据集分成 10 个部分
-
在 9 个部分上进行训练，剩下的部分用来测试
-
重复这个过程 10 次，以便每一部分都能作为测试集一次
准备工作
如前一节所示，加载你想要应用线性回归的数据集，这里是波士顿住房数据集：
from sklearn import datasets
boston = datasets.load_boston()
如何操作…
执行线性回归的步骤如下：
- 导入
LinearRegression对象并创建一个实例：
from sklearn.linear_model import LinearRegression
lr = LinearRegression()
- 将自变量和因变量传递给
LinearRegression的fit方法：
lr.fit(boston.data, boston.target)
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
- 现在，为了获得 10 折交叉验证的预测结果，执行以下操作：
from sklearn.model_selection import cross_val_predict
predictions_cv = cross_val_predict(lr, boston.data, boston.target, cv=10)
- 观察残差，即真实数据与预测值之间的差异，它们比前一节中没有交叉验证的线性回归的残差更接近正态分布：
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
#within Ipython
%matplotlib inline
pd.Series(boston.target - predictions_cv).hist(bins=50)
- 以下是通过交叉验证获得的新残差。与没有交叉验证时的情况相比，正态分布变得更加对称：
评估线性回归模型
在这个步骤中，我们将查看我们的回归如何拟合基础数据。我们在上一节中拟合了一个回归模型，但没有太关注我们实际的拟合效果。拟合模型后，首先要问的问题显然是：模型拟合得如何？在这个步骤中，我们将深入探讨这个问题。
准备工作
让我们使用 lr 对象和波士顿数据集——回到你在 数据拟合直线 章节中的代码。现在，lr 对象会有很多有用的方法，因为模型已经拟合完毕。
如何操作…
- 从 IPython 开始，导入多个库，包括
numpy、pandas和matplotlib用于可视化：
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
- 值得查看 Q-Q 图。我们这里使用
scipy，因为它有内置的概率图：
from scipy.stats import probplot
f = plt.figure(figsize=(7, 5))
ax = f.add_subplot(111)
tuple_out = probplot(boston.target - predictions_cv, plot=ax)
以下屏幕截图展示了概率图：
- 输入
tuple_out[1]，你将得到以下结果：
(4.4568597454452306, -2.9208080837569337e-15, 0.94762914118318298)
这是一个形式为 (slope, intercept, r) 的元组，其中 slope 和 intercept 来自最小二乘拟合，而 r 是决定系数的平方根。
- 在这里，我们之前看到的偏斜值变得更清晰了。我们还可以查看一些其他的拟合指标；均方误差 (MSE) 和 均值绝对偏差 (MAD) 是两种常见的指标。让我们在 Python 中定义每个指标并使用它们。
def MSE(target, predictions):
squared_deviation = np.power(target - predictions, 2)
return np.mean(squared_deviation)
MSE(boston.target, predictions)
21.897779217687503
def MAD(target, predictions):
absolute_deviation = np.abs(target - predictions)
return np.mean(absolute_deviation)
MAD(boston.target, predictions)
3.2729446379969205
- 现在你已经看到了使用 NumPy 计算误差的公式，你还可以使用
sklearn.metrics模块快速获取误差：
from sklearn.metrics import mean_absolute_error, mean_squared_error
print 'MAE: ', mean_absolute_error(boston.target, predictions)
print 'MSE: ', mean_squared_error(boston.target, predictions)
'MAE: ', 3.2729446379969205
'MSE: ', 21.897779217687503
它是如何工作的…
MSE 的公式非常简单：
它将每个预测值与实际值之间的偏差平方后再取平均，这实际上是我们优化的目标，用于找到线性回归的最佳系数。高斯-马尔科夫定理实际上保证了线性回归的解是最优的，因为系数具有最小的期望平方误差且是无偏的。在下一节中，我们将探讨当我们允许系数存在偏差时会发生什么。MAD 是绝对误差的期望误差：
在拟合线性回归时，MAD 并没有被使用，但它值得一看。为什么？想一想每个指标的作用，以及在每种情况下哪些错误更重要。例如，在 MSE 中，较大的错误会比其他项受到更多的惩罚，因为平方项的存在。离群值有可能显著地扭曲结果。
还有更多…
有一点被略过了，那就是系数本身是随机变量，因此它们具有分布。让我们使用自助法（bootstrapping）来查看犯罪率系数的分布。自助法是一种常见的技术，用来了解估计的不确定性：
n_bootstraps = 1000
len_boston = len(boston.target)
subsample_size = np.int(0.5*len_boston)
subsample = lambda: np.random.choice(np.arange(0, len_boston),size=subsample_size)
coefs = np.ones(n_bootstraps) #pre-allocate the space for the coefs
for i in range(n_bootstraps):
subsample_idx = subsample()
subsample_X = boston.data[subsample_idx]
subsample_y = boston.target[subsample_idx]
lr.fit(subsample_X, subsample_y)
coefs[i] = lr.coef_[0]
现在，我们可以查看系数的分布：
import matplotlib.pyplot as plt
f = plt.figure(figsize=(7, 5))
ax = f.add_subplot(111)
ax.hist(coefs, bins=50)
ax.set_title("Histogram of the lr.coef_[0].")
以下是生成的直方图：
我们还可能需要查看自助法（bootstrapping）生成的置信区间：
np.percentile(coefs, [2.5, 97.5])
array([-0.18497204, 0.03231267])
这很有趣；事实上，有理由相信犯罪率可能对房价没有影响。注意零值位于置信区间（CI）-0.18 到 0.03 之间，这意味着它可能不起作用。还值得指出的是，自助法可能会带来更好的系数估计，因为自助法的均值收敛速度比使用常规估计方法找出系数要快。
使用岭回归克服线性回归的不足
在本节中，我们将学习岭回归。它与普通的线性回归不同；它引入了一个正则化参数来收缩系数。当数据集包含共线性因素时，这非常有用。
岭回归在共线性存在的情况下非常强大，甚至可以建模多项式特征：向量 x、x²、x³，……这些特征高度共线且相关。
准备开始
让我们加载一个有效秩较低的数据集，并通过系数来比较岭回归和线性回归。如果你不熟悉秩，它是线性独立列和线性独立行中的较小者。线性回归的一个假设是数据矩阵是满秩的。
如何操作…
- 首先，使用
make_regression创建一个包含三个预测变量的简单数据集，但其effective_rank为2。有效秩意味着，尽管矩阵在技术上是满秩的，但许多列之间存在高度的共线性：
from sklearn.datasets import make_regression
reg_data, reg_target = make_regression(n_samples=2000,n_features=3, effective_rank=2, noise=10)
- 首先，让我们回顾一下上一章使用自助法的常规线性回归：
import numpy as np
n_bootstraps = 1000
len_data = len(reg_data)
subsample_size = np.int(0.5*len_data)
subsample = lambda: np.random.choice(np.arange(0, len_data),size=subsample_size)
coefs = np.ones((n_bootstraps, 3))
for i in range(n_bootstraps):
subsample_idx = subsample()
subsample_X = reg_data[subsample_idx]
subsample_y = reg_target[subsample_idx]
lr.fit(subsample_X, subsample_y)
coefs[i][0] = lr.coef_[0]
coefs[i][1] = lr.coef_[1]
coefs[i][2] = lr.coef_[2]
- 可视化系数：
- 使用岭回归执行相同的步骤：
from sklearn.linear_model import Ridge
r = Ridge()
n_bootstraps = 1000
len_data = len(reg_data)
subsample_size = np.int(0.5*len_data)
subsample = lambda: np.random.choice(np.arange(0, len_data),size=subsample_size)
coefs_r = np.ones((n_bootstraps, 3))
for i in range(n_bootstraps):
subsample_idx = subsample()
subsample_X = reg_data[subsample_idx]
subsample_y = reg_target[subsample_idx]
r.fit(subsample_X, subsample_y)
coefs_r[i][0] = r.coef_[0]
coefs_r[i][1] = r.coef_[1]
coefs_r[i][2] = r.coef_[2]
- 可视化结果：
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 5))
ax1 = plt.subplot(311, title ='Coef 0')
ax1.hist(coefs_r[:,0])
ax2 = plt.subplot(312,sharex=ax1, title ='Coef 1')
ax2.hist(coefs_r[:,1])
ax3 = plt.subplot(313,sharex=ax1, title ='Coef 2')
ax3.hist(coefs_r[:,2])
plt.tight_layout()
别被图表的相似宽度欺骗；岭回归的系数要接近零。
- 让我们看一下系数之间的平均差异：
np.var(coefs, axis=0)
array([ 228.91620444, 380.43369673, 297.21196544])
所以，平均而言，线性回归的系数要远高于岭回归的系数。这种差异是系数中的偏差（暂时忽略线性回归系数可能存在的偏差）。那么，岭回归的优势是什么呢？
- 好吧，来看一下我们系数的方差：
np.var(coefs_r, axis=0)
array([ 19.28079241, 15.53491973, 21.54126386])
方差已经显著降低。这就是机器学习中常常讨论的偏差-方差权衡。接下来的内容将介绍如何调整岭回归中的正则化参数，这是这一权衡的核心。
优化岭回归参数
一旦你开始使用岭回归进行预测或了解你正在建模的系统中的关系，你就会开始考虑 alpha 的选择。
例如，使用普通最小二乘 (OLS) 回归可能会显示两个变量之间的关系；然而，当通过 alpha 进行正则化时，这种关系不再显著。这可能是一个是否需要做出决策的问题。
做好准备
通过交叉验证，我们将调整岭回归的 alpha 参数。如果你还记得，在岭回归中，gamma 参数通常在调用 RidgeRegression 时被表示为 alpha，因此，出现的问题是，最优的 alpha 值是什么？创建一个回归数据集，然后我们开始吧：
from sklearn.datasets import make_regression
reg_data, reg_target = make_regression(n_samples=100, n_features=2, effective_rank=1, noise=10)
如何进行…
在 linear_models 模块中，有一个叫做 RidgeCV 的对象，代表岭回归交叉验证。它执行的交叉验证类似于 留一交叉验证 (LOOCV)。
- 在背后，它将为除了一个样本之外的所有样本训练模型。然后它会评估预测这个测试样本的误差：
from sklearn.linear_model import RidgeCV
rcv = RidgeCV(alphas=np.array([.1, .2, .3, .4]))
rcv.fit(reg_data, reg_target)
- 在我们拟合回归之后，alpha 属性将是最佳的 alpha 选择：
rcv.alpha_
0.10000000000000001
- 在之前的例子中，它是第一个选择。我们可能想要集中关注
.1附近的值：
rcv2 = RidgeCV(alphas=np.array([.08, .09, .1, .11, .12]))
rcv2.fit(reg_data, reg_target)
rcv2.alpha_
0.080000000000000002
我们可以继续这个探索，但希望这些机制已经很清晰了。
它是如何工作的…
这些机制可能很清楚，但我们应该再谈谈为什么，并定义最优值的含义。在交叉验证的每个步骤中，模型会对测试样本进行误差评分。默认情况下，本质上是平方误差。
我们可以强制 RidgeCV 对象存储交叉验证值；这将帮助我们可视化它所做的工作：
alphas_to_test = np.linspace(0.01, 1)
rcv3 = RidgeCV(alphas=alphas_to_test, store_cv_values=True)
rcv3.fit(reg_data, reg_target)
如你所见，我们测试了从 0.01 到 1 的一堆点（共 50 个）。由于我们将 store_cv_values 设置为 True，我们可以访问这些值：
rcv3.cv_values_.shape
(100L, 50L)
所以，在初始回归中我们有 100 个值，并测试了 50 个不同的 alpha 值。现在我们可以访问所有 50 个值的误差。因此，我们可以找到最小的均值误差，并选择它作为 alpha：
smallest_idx = rcv3.cv_values_.mean(axis=0).argmin()
alphas_to_test[smallest_idx]
0.030204081632653063
这与 RidgeCV 类的 rcv3 实例找到的最佳值一致：
rcv3.alpha_
0.030204081632653063
也值得可视化正在发生的事情。为此，我们将绘制所有 50 个测试 alpha 的均值：
plt.plot(alphas_to_test, rcv3.cv_values_.mean(axis=0))
还有更多…
如果我们想使用自己的评分函数，也可以这么做。由于我们之前查找了 MAD，我们就用它来评分差异。首先，我们需要定义我们的损失函数。我们将从 sklearn.metrics 导入它：
from sklearn.metrics import mean_absolute_error
在定义损失函数后，我们可以使用 sklearn 中的 make_scorer 函数。这将确保我们的函数被标准化，从而让 scikit 的对象知道如何使用它。此外，因为这是一个损失函数而不是评分函数，所以越低越好，因此需要让 sklearn 翻转符号，将其从最大化问题转变为最小化问题：
from sklearn.metrics import make_scorer
MAD_scorer = make_scorer(mean_absolute_error, greater_is_better=False)
按照之前的方式继续寻找最小的负 MAD 分数：
rcv4 = RidgeCV(alphas=alphas_to_test, store_cv_values=True, scoring=MAD_scorer)
rcv4.fit(reg_data, reg_target)
smallest_idx = rcv4.cv_values_.mean(axis=0).argmin()
查看最低的得分：
rcv4.cv_values_.mean(axis=0)[smallest_idx]
-0.021805192650070034
它发生在 alpha 为 0.01 时：
alphas_to_test[smallest_idx]
0.01
贝叶斯岭回归
此外，scikit-learn 还包含贝叶斯岭回归，它允许轻松估计置信区间。（注意，获取以下贝叶斯岭回归置信区间特别需要 scikit-learn 0.19.0 或更高版本。）
创建一条斜率为3且没有截距的直线，简化起见：
X = np.linspace(0, 5)
y_truth = 3 * X
y_noise = np.random.normal(0, 0.5, len(y_truth)) #normally distributed noise with mean 0 and spread 0.1
y_noisy = (y_truth + y_noise)
导入、实例化并拟合贝叶斯岭回归模型。请注意，一维的X和y变量必须进行重塑：
from sklearn.linear_model import BayesianRidge
br_inst = BayesianRidge().fit(X.reshape(-1, 1), y_noisy)
编写以下代码以获取正则化线性回归的误差估计：
y_pred, y_err = br_inst.predict(X.reshape(-1, 1), return_std=True)
绘制结果。噪声数据是蓝色的点，绿色的线条大致拟合它：
plt.figure(figsize=(7, 5))
plt.scatter(X, y_noisy)
plt.title("Bayesian Ridge Line With Error Bars")
plt.errorbar(X, y_pred, y_err, color='green')
最后，关于贝叶斯岭回归的补充说明，你可以通过交叉验证网格搜索对参数alpha_1、alpha_2、lambda_1和lambda_2进行超参数优化。
使用稀疏性来正则化模型
最小绝对收缩与选择算子（LASSO）方法与岭回归非常相似，也与最小角回归（LARS）类似。它与岭回归的相似之处在于我们通过一定的惩罚量来惩罚回归，而与 LARS 的相似之处在于它可以作为参数选择，通常会导致一个稀疏的系数向量。LASSO 和 LARS 都可以去除数据集中的许多特征，这取决于数据集的特点以及如何应用这些方法，这可能是你想要的，也可能不是。（而岭回归则保留所有特征，这使得你可以用来建模多项式函数或包含相关特征的复杂函数。）
正在准备
为了明确，LASSO 回归并非万能。使用 LASSO 回归可能会带来计算上的后果。正如我们在本配方中看到的，我们将使用一个不可微的损失函数，因此需要特定的、更重要的是、影响性能的解决方法。
如何操作…
执行 LASSO 回归的步骤如下：
- 让我们回到可靠的
make_regression函数，并使用相同的参数创建一个数据集：
import numpy as np
from sklearn.datasets import make_regression
reg_data, reg_target = make_regression(n_samples=200, n_features=500, n_informative=5, noise=5)
- 接下来，我们需要导入
Lasso对象：
from sklearn.linear_model import Lasso
lasso = Lasso()
- LASSO 包含许多参数，但最有趣的参数是
alpha。它用于调整Lasso方法的惩罚项。目前，保持它为 1。顺便说一下，就像岭回归一样，如果这个项为零，LASSO 就等同于线性回归：
lasso.fit(reg_data, reg_target)
- 再次查看有多少系数保持非零：
np.sum(lasso.coef_ != 0)
7
lasso_0 = Lasso(0)
lasso_0.fit(reg_data, reg_target)
np.sum(lasso_0.coef_ != 0)
500
我们的系数没有一个变为零，这正是我们预期的。实际上，如果你运行这个，你可能会收到来自 scikit-learn 的警告，建议你选择LinearRegression。
它是如何工作的…
LASSO 交叉验证 – LASSOCV
选择最合适的 lambda 是一个关键问题。我们可以自己指定 lambda，或者使用交叉验证根据现有数据找到最佳选择：
from sklearn.linear_model import LassoCV
lassocv = LassoCV()
lassocv.fit(reg_data, reg_target)
LASSOCV 将具有作为属性的最合适的 lambda。scikit-learn 在其符号中大多使用 alpha，但文献中使用 lambda：
lassocv.alpha_
0.75182924196508782
系数的数量可以通过常规方式访问：
lassocv.coef_[:5]
array([-0., -0., 0., 0., -0.])
让 LASSOCV 选择最合适的最佳拟合，我们得到15个非零系数：
np.sum(lassocv.coef_ != 0)
15
LASSO 用于特征选择
LASSO 通常可以用于其他方法的特征选择。例如，您可以运行 LASSO 回归来获取合适数量的特征，然后在其他算法中使用这些特征。
为了获取我们想要的特征，创建一个基于非零列的掩码数组，然后过滤掉非零列，保留我们需要的特征：
mask = lassocv.coef_ != 0
new_reg_data = reg_data[:, mask]
new_reg_data.shape
(200L, 15L)
采用更基本的方法进行 LARS 的正则化
借用 Gilbert Strang 对高斯消元法的评估，LARS 是一个你可能最终会考虑的想法，除非它已经在 Efron、Hastie、Johnstone 和 Tibshirani 的工作中被发现[1]。
准备工作
LARS 是一种回归技术，非常适合高维问题，即p >> n，其中p表示列或特征，n是样本的数量。
如何操作…
- 首先，导入必要的对象。我们使用的数据将包含 200 个数据点和 500 个特征。我们还将选择低噪声和少量的信息性特征：
from sklearn.datasets import make_regression
reg_data, reg_target = make_regression(n_samples=200, n_features=500, n_informative=10, noise=2)
- 由于我们使用了 10 个信息性特征，让我们也指定希望 LARS 中有 10 个非零系数。我们可能事先无法确切知道信息性特征的数量，但这对学习来说是有帮助的：
from sklearn.linear_model import Lars
lars = Lars(n_nonzero_coefs=10)
lars.fit(reg_data, reg_target)
- 然后我们可以验证 LARS 返回正确数量的非零系数：
np.sum(lars.coef_ != 0)
10
- 问题是为什么使用更少的特征更有用。为了说明这一点，让我们保留一半的数据并训练两个 LARS 模型，一个具有 12 个非零系数，另一个没有预定数量。我们在这里使用 12 是因为我们可能对重要特征的数量有所了解，但我们可能不确定确切的数量：
train_n = 100
lars_12 = Lars(n_nonzero_coefs=12)
lars_12.fit(reg_data[:train_n], reg_target[:train_n])
lars_500 = Lars() # it's 500 by default
lars_500.fit(reg_data[:train_n], reg_target[:train_n]);
np.mean(np.power(reg_target[train_n:] - lars_12.predict(reg_data[train_n:]), 2))
- 现在，为了查看每个特征如何拟合未知数据，请执行以下操作：
87.115080975821513
np.mean(np.power(reg_target[train_n:] - lars_500.predict(reg_data[train_n:]), 2))
2.1212501492030518e+41
如果你错过了，再看看；测试集上的误差显然非常高。这就是高维数据集的问题所在；给定大量特征，通常不难在训练样本上得到一个拟合良好的模型，但过拟合成为了一个巨大的问题。
它是如何工作的…
LARS 通过反复选择与残差相关的特征来工作。在几何上，相关性实际上是特征和残差之间的最小角度；这也是 LARS 得名的原因。
选择第一个特征后，LARS 将继续沿最小角度方向移动，直到另一个特征与残差的相关性达到相同的程度。然后，LARS 将开始沿这两个特征的组合方向移动。为了直观地理解这一点，考虑以下图表：
所以，我们沿着x1移动，直到x1受到y的拉力与x2受到y的拉力相等。发生这种情况时，我们沿着一条路径移动，这条路径的角度等于x1和x2之间的角度除以二。
还有更多…
就像我们使用交叉验证来调整岭回归一样，我们也可以对 LARS 做同样的操作：
from sklearn.linear_model import LarsCV
lcv = LarsCV()
lcv.fit(reg_data, reg_target)
使用交叉验证可以帮助我们确定使用的非零系数的最佳数量。这里，结果如下：
np.sum(lcv.coef_ != 0)
23
参考文献
- Bradley Efron, Trevor Hastie, Iain Johnstone, 和 Robert Tibshirani, 最小角回归, 《统计年鉴》32(2) 2004: 第 407–499 页, doi:10.1214/009053604000000067, MR2060166。
第五章：线性模型 – 逻辑回归
在本章中，我们将涵盖以下内容：
-
从 UCI 库加载数据
-
使用 pandas 查看 Pima 印第安人糖尿病数据集
-
查看 UCI Pima 印第安人数据集网页
-
使用逻辑回归进行机器学习
-
使用混淆矩阵检查逻辑回归错误
-
改变逻辑回归中的分类阈值
-
接收器操作特征 – ROC 分析
-
绘制没有上下文的 ROC 曲线
-
整合所有内容 – UCI 乳腺癌数据集
介绍
线性回归是一种非常古老的方法，是传统统计学的一部分。机器学习线性回归涉及训练集和测试集。通过这种方式，它可以与其他模型和算法通过交叉验证进行比较。传统线性回归则在整个数据集上进行训练和测试。这仍然是一种常见的做法，可能是因为线性回归往往是欠拟合，而非过拟合。
使用线性方法进行分类 – 逻辑回归
如第一章所示，高效机器学习 – NumPy，逻辑回归是一种分类方法。在某些情况下，它也可以作为回归器使用，因为它计算一个类别的实数概率，然后再进行分类预测。考虑到这一点，我们来探索加利福尼亚大学欧文分校（UCI）提供的 Pima 印第安人糖尿病数据集。
从 UCI 库加载数据
我们将加载的第一个数据集是 Pima 印第安人糖尿病数据集。这需要访问互联网。数据集由 Sigillito V.（1990）提供，存储在 UCI 机器学习库中（archive.ics.uci.edu/ml/machine-learning-databases/pima-indians-diabetes/pima-indians-diabetes.data），约翰·霍普金斯大学应用物理实验室，马里兰州劳雷尔。
如果你是开源老手，首先想到的可能是，这个数据库的许可证/授权是什么？这是一个非常重要的问题。UCI 库有一个使用政策，要求我们在使用数据库时引用它。我们可以使用它，但必须给予他们应有的赞扬并提供引用。
如何操作…
- 打开 IPython 并导入
pandas：
import pandas as pd
- 将 Pima 印第安人糖尿病数据集的网页地址作为字符串输入，如下所示：
data_web_address = "https://archive.ics.uci.edu/ml/machine-learning-databases/pima-indians-diabetes/pima-indians-diabetes.data"
- 将数据的列名输入为列表：
column_names = ['pregnancy_x',
'plasma_con',
'blood_pressure',
'skin_mm',
'insulin',
'bmi',
'pedigree_func',
'age',
'target']
- 将特征名称存储为一个列表。排除
target列，这是column_names中的最后一列，因为它不是特征：
feature_names = column_names[:-1]
- 创建一个 pandas 数据框来存储输入数据：
all_data = pd.read_csv(data_web_address , names=column_names)
使用 pandas 查看 Pima 印第安人糖尿病数据集
如何操作…
- 你可以通过多种方式查看数据。查看数据框的顶部：
all_data.head()
- 这里似乎没有什么问题，除了可能有一个胰岛素水平为零的情况。这可能吗？那
skin_mm变量呢？它可以是零吗？在你的 IPython 中做个注释：
#Is an insulin level of 0 possible? Is a skin_mm of 0 possible?
- 使用
describe()方法获取数据框的粗略概览：
all_data.describe()
- 再次在笔记本中做个注释，记录其他零值的情况：
#The features plasma_con, blood_pressure, skin_mm, insulin, bmi have 0s as values. These values could be physically impossible.
- 绘制
pregnancy_x变量的直方图。将hist()方法中的 bins 参数设置为 50，以获得更多的区间和更高分辨率的图像；否则，图像会难以阅读：
#If within a notebook, include this line to visualize within the notebook.
%matplotlib inline
#The default is bins=10 which is hard to read in the visualization.
all_data.pregnancy_x.hist(bins=50)
- 为数据框中的所有列绘制直方图。在方法中将
figsize设置为元组(15,9)，并再次将 bins 设置为50；否则，图像会很难读取：
all_data.hist(figsize=(15,9),bins=50)
https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/skl-cb-2e/img/5aaa8355-0f01-49a8-aa4a-df447d30b05c.pngblood_pressure和bmi看起来呈正态分布，除了异常的零值。
pedigree_func和plasma_con变量呈偏态正态分布（可能是对数正态分布）。age和pregnancy_x变量在某种程度上呈衰减状态。胰岛素和skin_mm变量看起来除了零值很多外，可能呈正态分布。
- 最后，注意
target变量中的类别不平衡。使用value_counts()pandas 系列方法重新检查这种不平衡：
all_data.target.value_counts()
0 500
1 268
Name: target, dtype: int64
在许多情况下，人的描述是类别零而不是类别一。
查看 UCI Pima Indians 数据集的网页
我们进行了初步的探索性分析，以大致了解数据。现在我们将阅读 UCI Pima Indians 数据集的文档。
怎么做…
查看引用政策
-
这里是关于 UCI Pima Indians 糖尿病数据集的所有信息。首先，滚动到页面底部，查看他们的引用政策。糖尿病数据集的通用 UCI 引用政策可以在以下链接找到：
archive.ics.uci.edu/ml/citation_policy.html。 -
通用政策表示，若要发布使用数据集的材料，请引用 UCI 存储库。
阅读关于缺失值和上下文的说明
- 页面顶部有重要链接和数据集摘要。摘要提到数据集中存在缺失值：
- 在摘要下方，有一个属性描述（这就是我最初为列命名的方式）：
- 那么这些类别变量到底意味着什么呢？目标变量中的零和一分别代表什么？为了弄清楚这一点，点击摘要上方的“数据集描述”链接。滚动到页面的第九点，那里会给出所需的信息：
这意味着 1 代表糖尿病的阳性结果。这是重要信息，并且在数据分析中提供了上下文。
- 最后，有一个免责声明指出：正如一个仓库用户所指出的那样，这不可能是正确的：在某些地方有零值，而这些地方在生物学上是不可能的，例如血压属性。很可能零值表示缺失数据。
因此，我们在数据探索阶段怀疑某些不可能出现的零值时是正确的。许多数据集都有损坏或缺失的值。
使用逻辑回归进行机器学习
你熟悉训练和测试分类器的步骤。使用逻辑回归时，我们将执行以下操作：
-
将数据加载到特征和目标数组
X和y中 -
将数据拆分为训练集和测试集
-
在训练集上训练逻辑回归分类器
-
在测试集上测试分类器的性能
准备中
定义 X、y——特征和目标数组
让我们开始使用 scikit-learn 的逻辑回归进行预测。首先执行必要的导入并设置输入变量 X 和目标变量 y：
import numpy as np
import pandas as pd
X = all_data[feature_names]
y = all_data['target']
如何操作…
提供训练集和测试集
- 导入
train_test_split以为X和y创建测试集和训练集：输入和目标。注意stratify=y，它会将分类变量y进行分层。这意味着y_train和y_test中的零和一的比例是相同的：
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=7,stratify=y)
训练逻辑回归模型
- 现在导入
LogisticRegression并将其拟合到训练数据上：
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression()
lr.fit(X_train,y_train)
- 在测试集上进行预测并将其存储为
y_pred：
y_pred = lr.predict(X_test)
对逻辑回归进行评分
- 使用
accuracy_score检查预测的准确性，这是分类正确的百分比：
from sklearn.metrics import accuracy_score
accuracy_score(y_test,y_pred)
0.74675324675324672
所以，我们得到了一个分数，但这个分数是我们在这些情况下可以使用的最佳度量吗？我们能做得更好吗？也许可以。请看一下下面的结果混淆矩阵。
使用混淆矩阵检查逻辑回归的错误
准备中
导入并查看我们构建的逻辑回归的混淆矩阵：
from sklearn.metrics import confusion_matrix
confusion_matrix(y_test, y_pred,labels = [1,0])
array([[27, 27],
[12, 88]])
我向混淆矩阵传递了三个参数：
-
y_test：测试目标集 -
y_pred：我们的逻辑回归预测 -
labels：指代阳性类别
labels = [1,0] 表示阳性类别为 1，阴性类别为 0。在医学背景下，我们在探索 Pima 印第安糖尿病数据集时发现，类别 1 测试为糖尿病阳性。
这是混淆矩阵，再次以 pandas 数据框形式呈现：
如何操作…
阅读混淆矩阵
这个小的数字数组具有以下含义：
混淆矩阵不仅告诉我们分类过程中发生了什么，还能提供准确率评分。矩阵从左上到右下的对角线元素表示正确的分类。共有 27 + 88 = 115 个正确分类。对角线之外，共犯了 27 + 12 = 39 个错误分类。注意，115 / (115 + 39)再次得到分类器的准确度，大约为 0.75。
让我们再次关注错误。在混淆矩阵中，27 个人被错误地标记为没有糖尿病，尽管他们确实患有糖尿病。在实际情况下，这是比那些被认为患有糖尿病但实际上没有的人更严重的错误。第一类可能会回家后忘记，而第二类则可能会重新测试。
上下文中的一般混淆矩阵
一般的混淆矩阵，其中正类表示识别一种病症（此处为糖尿病），因此具有医学诊断背景：
改变逻辑回归中的分类阈值
准备工作
我们将利用以下事实：在逻辑回归分类中，存在回归过程，用于最小化那些本应被诊断为糖尿病的人却被送回家的次数。通过调用估算器的predict_proba()方法来实现：
y_pred_proba = lr.predict_proba(X_test)
这将得到一个概率数组。查看该数组：
y_pred_proba
array([[ 0.87110309, 0.12889691],
[ 0.83996356, 0.16003644],
[ 0.81821721, 0.18178279],
[ 0.73973464, 0.26026536],
[ 0.80392034, 0.19607966], ...
在第一行，类别0的概率约为 0.87，类别1的概率为 0.13。注意，作为概率，这些数字加起来为 1。由于只有两类，这个结果可以看作是一个回归器，一个关于类别为1（或0）概率的实数。通过直方图可视化类别为1的概率。
取数组的第二列，将其转换为 pandas 序列并绘制直方图：
pd.Series(y_pred_proba[:,1]).hist()
在概率直方图中，相比选择 1 的高概率，低概率的选择较多。例如，选择 1 的概率通常在 0.1 到 0.2 之间。在逻辑回归中，算法默认只有在概率大于 0.5（即一半）时才会选择 1。现在，将其与开始时的目标直方图进行对比：
all_data['target'].hist()
在接下来的步骤中，我们将：
-
调用类方法
y_pred_proba() -
使用具有特定阈值的
binarize函数 -
查看由阈值生成的混淆矩阵
如何操作…
- 要基于阈值选择分类类，请使用
preprocessing模块中的binarize。首先导入它：
from sklearn.preprocessing import binarize
- 查看
y_pred_proba的前两列：
array([[ 0.87110309, 0.12889691],
[ 0.83996356, 0.16003644]
- 然后尝试对
y_pred_proba使用binarize函数，设置阈值为0.5。查看结果：
y_pred_default = binarize(y_pred_proba, threshold=0.5)
y_pred_default
array([[ 1., 0.],
[ 1., 0.],
[ 1., 0.],
[ 1., 0.],
[ 1., 0.],
[ 1., 0.]
binarize函数如果y_pred_proba中的值大于 0.5，则用1填充数组；否则填入0。在第一行，0.87 大于 0.5，而 0.13 不大于。因此，为了二值化，将 0.87 替换为1，将 0.13 替换为0。现在，查看y_pred_default的第一列：
y_pred_default[:,1]
array([ 0., 0., 0., 0., 0., 0 ...
这恢复了默认逻辑回归分类器在阈值0.5下所做的决策。
- 在 NumPy 数组上尝试混淆矩阵函数会得到我们遇到的第一个混淆矩阵（请注意，标签再次选择为
[1,0]）：
confusion_matrix(y_test, y_pred_default[:,1],labels = [1,0])
array([[27, 27],
[12, 88]])
- 尝试不同的阈值，使得类别
1有更高的选择概率。查看其混淆矩阵：
y_pred_low = binarize(y_pred_proba, threshold=0.2)
confusion_matrix(y_test, y_pred_low[:,1],labels=[1,0]) #positive class is 1 again
array([[50, 4],
[48, 52]])
通过更改阈值，我们增加了预测类别1的概率——增加了混淆矩阵第一列中数字的大小。现在，第一列的总和为 50 + 48 = 98。之前，第一列是 27 + 12 = 39，数字要小得多。现在只有四个人被错误分类为没有糖尿病，尽管他们实际上有。请注意，在某些情况下，这是一件好事。
当算法预测为零时，它很可能是正确的。当它预测为一时，它通常不准确。假设你经营一家医院。你可能喜欢这个测试，因为你很少会送走一个看似没有糖尿病的人，尽管他实际上有糖尿病。每当你送走那些确实有糖尿病的人时，他们就无法尽早治疗，导致医院、保险公司和他们自己都面临更高的费用。
你可以衡量测试在预测为零时的准确性。观察混淆矩阵的第二列，[4, 52]。在这种情况下，准确率为 52 / (52 + 4)，大约是 0.93。这叫做负预测值（NPV）。你可以编写一个函数来根据阈值计算 NPV：
from __future__ import division #In Python 2.x
import matplotlib.pyplot as plt
def npv_func(th):
y_pred_low = binarize(y_pred_proba, threshold=th)
second_column = confusion_matrix(y_test, y_pred_low[:,1],labels=[1,0])[:,1]
npv = second_column[1]/second_column.sum()
return npv
npv_func(0.2)
0.9285714285714286
然后绘制它：
ths = np.arange(0,1,0.05)
npvs = []
for th in np.arange(0,1.00,0.05):
npvs.append(npv_func(th))
plt.plot(ths,npvs)
接收者操作特征 – ROC 分析
沿着检查 NPV 的思路，还有一些标准的度量方法用于检查混淆矩阵中的单元格。
准备就绪
敏感性
敏感性，像上一节中的 NPV 一样，是混淆矩阵单元格的数学函数。敏感性是接受测试且实际患有疾病的人群中，正确标记为有疾病的人群的比例，在这个例子中是糖尿病：
从数学角度看，它是正确标记为有疾病（TP）的患者与实际有疾病的所有人（TP + FN）之比。首先，回顾混淆矩阵的单元格。重点关注真实行，对应于所有患有糖尿病的人：
因此：
敏感性另一个名字是真正阳性率（TPR）。
从视觉角度来看
从另一个角度来看，混淆矩阵非常直观。让我们用直方图（下一张图中的左列）来可视化有糖尿病的阳性类别和没有糖尿病的阴性类别。每个类别大致呈正态分布。借助 SciPy，我可以找到最佳拟合的正态分布。右下角注明了阈值被设置为 0.5，这是逻辑回归中的默认设置。观察阈值是如何让我们选择假阴性（FN）和假阳性（FP）的，尽管这种选择并不完美。
如何操作…
在 scikit-learn 中计算 TPR
- scikit-learn 提供了方便的函数，用于根据阳性类别的概率向量
y_pred_proba[:,1]来计算逻辑回归的敏感度或 TPR：
from sklearn.metrics import roc_curve
fpr, tpr, ths = roc_curve(y_test, y_pred_proba[:,1])
这里，给定阳性类别向量，scikit-learn 中的roc_curve函数返回了一个包含三个数组的元组：
-
TPR 数组（记作
tpr） -
FPR 数组（记作
fpr） -
用于计算 TPR 和 FPR 的自定义阈值集（记作
ths）
详细说明假阳性率（FPR），它描述了假警报的发生率。它是错误地认为没有糖尿病的人却被判定为有糖尿病的比例：
这是一个没有糖尿病的人的陈述。从数学角度来看，混淆矩阵的单元格如下：
绘制敏感度图
- 在y轴绘制敏感度，在x轴绘制阈值：
plt.plot(ths,tpr)
- 因此，阈值越低，敏感度越高。查看阈值为
0.1时的混淆矩阵：
y_pred_th= binarize(y_pred_proba, threshold=0.1)
confusion_matrix(y_test, y_pred_th[:,1],labels=[1,0])
array([[54, 0],
[81, 19]])
- 在这种情况下，没有人会误以为自己得了糖尿病。然而，就像我们在计算 NPV 时的情况一样，当测试预测某人患有糖尿病时，它的准确性非常低。这种情况下的最佳情形是当阈值为
0.146时：
y_pred_th = binarize(y_pred_proba, threshold=0.146)
confusion_matrix(y_test, y_pred_th[:,1],labels=[1,0])
array([[54, 0],
[67, 33]])
即便如此，当预测某人患有糖尿病时，测试结果也不起作用。它的有效性为 33 / (33 + 121) = 0.21，或者说 21%的准确率。
还有更多…
非医疗背景下的混淆矩阵
假设你是一位银行工作人员，想要确定一个客户是否值得获得购买房屋的抵押贷款。接下来是一个可能的混淆矩阵，展示了根据银行可用的客户数据是否应该给某人抵押贷款。
任务是对人群进行分类，并决定他们是否应该获得抵押贷款。在这种情况下，可以为每种情况分配一个数字。当混淆矩阵中的每个单元格都有明确的成本时，找到最佳分类器就更容易了：
绘制没有上下文的 ROC 曲线
如何操作…
ROC 曲线是任何分类器的诊断工具，没有任何上下文。没有上下文意味着我们尚不知道哪个错误类型（FP 或 FN）是更不可取的。让我们立即使用概率向量 y_pred_proba[:,1] 绘制它：
from sklearn.metrics import roc_curve
fpr, tpr, ths = roc_curve(y_test, y_pred_proba[:,1])
plt.plot(fpr,tpr)
ROC 曲线是 FPR（假警报）在 x 轴和 TPR（找出所有有病且确实有病的人）在 y 轴上的图示。没有上下文时，它是一个衡量分类器性能的工具。
完美分类器
完美的分类器无论如何都会有 TPR 为 1，假警报率（FAR）：
在前面的图表中，FN 非常小；因此 TPR，TP / (TP + FN)，接近 1. 其 ROC 曲线呈 L 形：
不完美的分类器
在以下图像中，分布重叠，类别之间无法区分：
在不完美的分类器中，FN 和 TN 几乎相等，FP 和 TP 也几乎相等。因此，通过替代，TPR TP/ (TP + FP) 几乎等于假阴性率（FNR）FP/ (FP + TN)。即使您改变阈值，这也是成立的。因此，我们得到的 ROC 曲线是斜率大约为 1 的直线：
AUC – ROC 曲线下面积
L 形完美分类器的面积是 1 x 1 = 1. 坏分类器的面积是 0.5. 为了衡量分类器的性能，scikit-learn 提供了一个便捷的 ROC 曲线下面积（AUC）计算函数：
from sklearn.metrics import auc
auc(fpr,tpr)
0.825185185185
将所有内容整合在一起 – UCI 乳腺癌数据集
如何做到…
数据集感谢 Street, N (1990)，UCI 机器学习库提供 (archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data)，威斯康星大学，麦迪逊，WI：计算机科学系：
- 阅读引用/许可证信息后，从 UCI 加载数据集：
import numpy as np
import pandas as pd
data_web_address = data_web_address = "https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data"
column_names = ['radius',
'texture',
'perimeter',
'area',
'smoothness'
,'compactness',
'concavity',
'concave points',
'symmetry',
'malignant']
feature_names = column_names[:-1]
all_data = pd.read_csv(data_web_address , names=column_names)
- 查看数据类型：
all_data.dtypes
radius int64
texture int64
perimeter int64
area int64
smoothness int64
compactness object
concavity int64
concave points int64
symmetry int64
malignant int64
dtype: object
结果表明，特征紧凑度具有像 ? 这样的特征。目前，我们不使用这个特征。
- 现在，阅读文档，目标变量设置为
2（没有癌症）和4（有癌症）。将变量更改为0代表没有癌症，1代表有癌症：
#changing the state of having cancer to 1, not having cancer to 0
all_data['malignant'] = all_data['malignant'].astype(np.int)
all_data['malignant'] = np.where(all_data['malignant'] == 4, 1,0) #4, and now 1 means malignant
all_data['malignant'].value_counts()
0 458
1 241
Name: malignant, dtype: int64
- 定义
X和y：
X = all_data[[col for col in feature_names if col != 'compactness']]
y = all_data.malignant
- 将
X和y分割成训练集和测试集：
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=7,stratify=y)
#Train and test the Logistic Regression. Use the method #predict_proba().
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression()
lr.fit(X_train,y_train)
y_pred_proba = lr.predict_proba(X_test)
- 绘制 ROC 曲线并计算 AUC 分数：
import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, auc, roc_auc_score
fpr, tpr, ths = roc_curve(y_test, y_pred_proba[:,1])
auc_score = auc(fpr,tpr)
plt.plot(fpr,tpr,label="AUC Score:" + str(auc_score))
plt.xlabel('fpr',fontsize='15')
plt.ylabel('tpr',fontsize='15')
plt.legend(loc='best')
这个分类器表现得相当不错。
未来项目的概要
总体而言，对于未来的分类项目，您可以执行以下操作：
-
加载您能找到的最佳数据，以解决特定问题。
-
确定是否有分类上下文：是 FP 更好还是 FN 更好？
-
使用 ROC-AUC 分数对数据进行无上下文的训练和测试。如果几个算法在此步骤中表现不佳，您可能需要回到第 1 步。
-
如果上下文很重要，可以使用混淆矩阵进行探索。
逻辑回归特别适合执行这些步骤，尽管任何具有 predict_proba() 方法的算法也能非常相似地工作。作为一个练习，如果您有雄心，可以将此过程推广到其他算法甚至是一般回归。这里的主要观点是，并非所有的错误都是相同的，这一点在健康数据集中尤其重要，因为治疗所有患有某种病症的患者非常关键。
关于乳腺癌数据集的最后一点：观察数据由细胞测量值组成。您可以通过自动化计算机查看图片并使用传统程序或机器学习找到这些测量值来收集这些数据。
第六章：使用距离度量构建模型
本章将介绍以下配方：
-
使用 k-means 对数据进行聚类
-
优化质心数量
-
评估聚类的正确性
-
使用 MiniBatch k-means 处理更多数据
-
使用 k-means 聚类对图像进行量化
-
在特征空间中找到最接近的对象
-
使用高斯混合模型进行概率聚类
-
使用 k-means 进行异常值检测
-
使用 KNN 进行回归
介绍
本章我们将讨论聚类。聚类通常与无监督技术一起使用。这些技术假设我们不知道结果变量。这会导致实践中的结果和目标出现模糊，但尽管如此，聚类仍然很有用。正如我们将看到的那样，我们可以在监督环境中使用聚类来定位我们的估计。这或许是聚类如此有效的原因；它可以处理各种情况，且结果通常是，在没有更好术语的情况下，可以说是合理的。
本章我们将涵盖各种应用，从图像处理到回归和异常值检测。聚类与类别的分类相关。你有一个有限的簇或类别集。与分类不同，你事先并不知道这些类别。此外，聚类通常可以通过连续的、概率性的或优化的视角来看待。
不同的解释会导致不同的权衡。我们将讨论如何在这里拟合模型，这样当你遇到聚类问题时，你就有工具可以尝试多种模型。
使用 k-means 对数据进行聚类
在数据集中，我们观察到一些点聚集在一起。使用 k-means，我们将把所有点分到不同的组或簇中。
准备工作
首先，让我们演示一些简单的聚类操作，然后我们将讨论 k-means 如何工作：
import numpy as np
import pandas as pd
from sklearn.datasets import make_blobs
blobs, classes = make_blobs(500, centers=3)
此外，由于我们将进行一些绘图操作，按如下方式导入 matplotlib：
import matplotlib.pyplot as plt
%matplotlib inline #Within an ipython notebook
如何进行……
我们将通过一个简单的示例，演示如何对虚拟数据的簇进行聚类。然后，我们将简要介绍 k-means 如何工作，以寻找最佳的簇数：
- 观察我们的簇，我们可以看到有三个明显的聚类：
f, ax = plt.subplots(figsize=(7.5, 7.5))
rgb = np.array(['r', 'g', 'b'])
ax.scatter(blobs[:, 0], blobs[:, 1], color=rgb[classes])
ax.set_title("Blobs")
现在我们可以使用 k-means 找到这些簇的中心。
- 在第一个示例中，我们假设我们知道有三个中心：
from sklearn.cluster import KMeans
kmean = KMeans(n_clusters=3)
kmean.fit(blobs)
KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
n_clusters=3, n_init=10, n_jobs=1, precompute_distances='auto',
random_state=None, tol=0.0001, verbose=0)
kmean.cluster_centers_
array([[ 3.48939154, -0.92786786],
[-2.05114953, 1.58697731],
[ 1.58182736, -6.80678064]])
f, ax = plt.subplots(figsize=(7.5, 7.5))
ax.scatter(blobs[:, 0], blobs[:, 1], color=rgb[classes])
ax.scatter(kmean.cluster_centers_[:, 0],kmean.cluster_centers_[:, 1], marker='*', s=250,color='black', label='Centers')
ax.set_title("Blobs")
ax.legend(loc='best')
- 以下截图显示了输出结果：
- 其他属性也很有用。例如，
labels_属性将为每个点生成预期的标签：
kmean.labels_[:5]
array([2, 0, 1, 1, 0])
我们可以检查 kmean.labels_ 是否与类别相同，但由于 k-means 在开始时并不知道类别，它无法将样本索引值分配给两个类别：
classes[:5]
array([2, 0, 1, 1, 0])
随时交换类别中的 1 和 0，检查它们是否与 labels_ 匹配。变换函数非常有用，因为它会输出每个点与质心之间的距离：
kmean.transform(blobs)[:5]
array([[ 6.75214231, 9.29599311, 0.71314755], [ 3.50482136, 6.7010513 , 9.68538042], [ 6.07460324, 1.91279125, 7.74069472], [ 6.29191797, 0.90698131, 8.68432547], [ 2.84654338, 6.07653639, 3.64221613]])
它是如何工作的……
k-means 实际上是一个非常简单的算法，旨在最小化群内距离均值的平方和。我们将再次最小化平方和！
它首先设置一个预设的簇数 K，然后在以下步骤之间交替进行：
-
将每个观察值分配给最近的簇
-
通过计算分配给该簇的每个观察值的均值来更新每个质心
这一过程会一直持续，直到满足某个特定的标准。质心很难解释，而且确定我们是否有正确数量的质心也非常困难。理解数据是否标记过是很重要的，因为这将直接影响你可以使用的评估方法。
优化质心数量
在进行 k-means 聚类时，我们实际上无法事先知道正确的簇数，因此找出这一点是一个重要步骤。一旦我们知道（或估计）了质心的数量，问题就开始更像一个分类问题，因为我们可以使用的知识量大大增加了。
准备开始
对无监督技术的模型性能进行评估是一项挑战。因此，sklearn提供了几种在已知真实标签的情况下评估聚类的方法，而在未知真实标签的情况下，方法则非常有限。
我们将从单一簇模型开始并评估其相似度。这更多是为了操作上的需要，因为测量单一簇数量的相似度显然对于找到真实的簇数并没有什么用处。
如何做…
- 为了开始，我们将创建几个簇状数据，用于模拟数据的聚类：
from sklearn.datasets import make_blobs
import numpy as np
blobs, classes = make_blobs(500, centers=3)
from sklearn.cluster import KMeans
kmean = KMeans(n_clusters=3)
kmean.fit(blobs)
KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
n_clusters=3, n_init=10, n_jobs=1, precompute_distances='auto',
random_state=None, tol=0.0001, verbose=0)
- 首先，我们来看轮廓距离。轮廓距离是群内不相似度与最接近的群外不相似度之间的差值与这两个值的最大值之比。可以将其视为衡量簇与簇之间分离度的标准。让我们看看从点到簇中心的距离分布；这对于理解轮廓距离很有帮助：
from sklearn import metrics
silhouette_samples = metrics.silhouette_samples(blobs, kmean.labels_)
np.column_stack((classes[:5], silhouette_samples[:5]))
array([[ 0\. , 0.69568017],
[ 0\. , 0.76789931],
[ 0\. , 0.62470466],
[ 0\. , 0.6266658 ],
[ 2\. , 0.63975981]])
- 以下是部分输出：
- 注意，通常情况下，接近 1 的系数越多（这很好），得分越高。
它是如何工作的…
轮廓系数的平均值通常用于描述整个模型的拟合度：
silhouette_samples.mean()
0.5633513643546264
这是非常常见的，实际上，metrics 模块暴露了一个函数，可以得到我们刚刚得到的值：
metrics.silhouette_score(blobs, kmean.labels_)
0.5633513643546264
import matplotlib.pyplot as plt
%matplotlib inline
blobs, classes = make_blobs(500, centers=10)
silhouette_avgs = []
for k in range(2, 60):
kmean = KMeans(n_clusters=k).fit(blobs)
silhouette_avgs.append(metrics.silhouette_score(blobs, kmean.labels_))
f, ax = plt.subplots(figsize=(7, 5))
ax.plot(silhouette_avgs)
以下是输出：
这张图显示了随着质心数量的增加，轮廓系数的平均值变化。我们可以看到，根据数据生成过程，最佳的簇数是 3；但在这里，看起来它大约是 7 或 8。这就是聚类的现实；通常，我们无法获得正确的簇数。我们只能希望对簇数进行某种程度的估算。
评估簇的正确性
我们稍微谈了一下当真实标签未知时如何评估聚类。然而，我们还没有讨论如何在已知簇的情况下评估 k-means。在很多情况下，这个是不可知的；然而，如果有外部标注，我们将知道真实标签或至少知道某种代理标签。
准备好了
所以，让我们假设一个世界，假设有一个外部代理为我们提供真实标签。
我们将创建一个简单的数据集，以几种方式评估与真实标签的正确性度量，然后进行讨论：
from sklearn import datasets
from sklearn import cluster
blobs, ground_truth = datasets.make_blobs(1000, centers=3,cluster_std=1.75)
如何操作…
- 在讲解度量标准之前，我们先来看一下数据集：
%matplotlib inline
import matplotlib.pyplot as plt
f, ax = plt.subplots(figsize=(7, 5))
colors = ['r', 'g', 'b']
for i in range(3):
p = blobs[ground_truth == i]
ax.scatter(p[:,0], p[:,1], c=colors[i],
label="Cluster {}".format(i))
ax.set_title("Cluster With Ground Truth")
ax.legend()
- 为了拟合一个 k-means 模型，我们将从聚类模块创建一个
KMeans对象：
kmeans = cluster.KMeans(n_clusters=3)
kmeans.fit(blobs)
KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
n_clusters=3, n_init=10, n_jobs=1, precompute_distances='auto',
random_state=None, tol=0.0001, verbose=0)
kmeans.cluster_centers_
array([[ 3.61594791, -6.6125572 ],
[-0.76071938, -2.73916602],
[-3.64641767, -6.23305142]])
- 现在我们已经拟合了模型，让我们看看聚类中心：
f, ax = plt.subplots(figsize=(7, 5))
colors = ['r', 'g', 'b']
for i in range(3):
p = blobs[ground_truth == i]
ax.scatter(p[:,0], p[:,1], c=colors[i], label="Cluster {}".format(i))
以下是输出结果：
- 现在我们可以将聚类性能视为分类练习，其相关的度量标准在这里也同样适用：
for i in range(3):
print (kmeans.labels_ == ground_truth)[ground_truth == i].astype(int).mean()
0.946107784431
0.135135135135
0.0750750750751
显然，我们有一些反向的簇。所以，让我们先解决这个问题，然后再看看准确性：
new_ground_truth = ground_truth.copy()
new_ground_truth[ground_truth == 1] = 2
new_ground_truth[ground_truth == 2] = 1
0.946107784431
0.852852852853
0.891891891892
所以，我们大约有 90% 的时间是正确的。我们将看的第二个相似性度量是互信息得分：
from sklearn import metrics
metrics.normalized_mutual_info_score(ground_truth, kmeans.labels_)
0.66467613668253844
当得分趋向 0 时，标签分配可能不是通过相似的过程生成的；然而，接近 1 的得分意味着两个标签之间有很大的相似性。
例如，让我们看看当互信息得分本身时会发生什么：
metrics.normalized_mutual_info_score(ground_truth, ground_truth)
1.0
从名称上看，我们可以推测这可能是一个未归一化的 mutual_info_score：
metrics.mutual_info_score(ground_truth, kmeans.labels_)
0.72971342940406325
这些非常接近；然而，归一化互信息是互信息除以每个集合真实标签和分配标签的熵乘积的平方根。
还有更多内容…
我们还没有讨论过的一个聚类度量，并且这个度量不依赖于真实标签的是惯性。它目前作为度量并没有得到很好的文档化。然而，它是 k-means 最小化的一个度量。
惯性是每个点与其分配的簇之间的平方差之和。我们可以通过一点 NumPy 来确定这一点：
kmeans.inertia_
4849.9842988128385
使用 MiniBatch k-means 处理更多数据
K-means 是一种不错的方法；然而，它对于大量数据并不理想。这是由于 k-means 的复杂性。话虽如此，我们可以通过 MiniBatch k-means 以更好的算法复杂度获得近似解。
准备好了
MiniBatch k-means 是 k-means 的更快实现。K-means 在计算上非常昂贵；该问题是 NP 难问题。
然而，使用 MiniBatch k-means，我们可以通过数量级加速 k-means。这是通过使用许多叫做 MiniBatches 的子样本来实现的。考虑到子抽样的收敛性，提供了良好的初始条件时，可以实现对常规 k-means 的接近近似。
如何操作…
- 让我们对 MiniBatch 聚类进行一些非常高层次的分析。首先，我们将看一下整体速度差异，然后再看一下估计误差：
import numpy as np
from sklearn.datasets import make_blobs
blobs, labels = make_blobs(int(1e6), 3)
from sklearn.cluster import KMeans, MiniBatchKMeans
kmeans = KMeans(n_clusters=3)
minibatch = MiniBatchKMeans(n_clusters=3)
理解这些指标的目的是为了揭示问题。因此，在确保基准的最高准确性方面非常小心。关于这个主题有很多信息可供参考；如果你真的想深入了解为什么 MiniBatch k-means 在扩展性方面更好，建议查看现有的资料。
- 现在设置完成了，我们可以测量时间差异：
%time kmeans.fit(blobs) #IPython Magic
Wall time: 7.88 s
KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
n_clusters=3, n_init=10, n_jobs=1, precompute_distances='auto',
random_state=None, tol=0.0001, verbose=0)
%time minibatch.fit(blobs)
Wall time: 2.66 s
MiniBatchKMeans(batch_size=100, compute_labels=True, init='k-means++',
init_size=None, max_iter=100, max_no_improvement=10, n_clusters=3,
n_init=3, random_state=None, reassignment_ratio=0.01, tol=0.0,
verbose=0)
- CPU 时间上有很大的差异。聚类性能的差异如下所示：
kmeans.cluster_centers_
array([[-3.74304286, -0.4289715 , -8.69684375],
[-5.73689621, -6.39166391, 6.18598804],
[ 0.63866644, -9.93289824, 3.24425045]])
minibatch.cluster_centers_
array([[-3.72580548, -0.46135647, -8.63339789],
[-5.67140979, -6.33603949, 6.21512625],
[ 0.64819477, -9.87197712, 3.26697532]])
- 查看这两个数组；找到的中心点是按照相同顺序排列的。这是随机的—聚类不必按相同顺序排列。查看第一个聚类中心之间的距离：
from sklearn.metrics import pairwise
pairwise.pairwise_distances(kmeans.cluster_centers_[0].reshape(1, -1), minibatch.cluster_centers_[0].reshape(1, -1))
array([[ 0.07328909]])
- 这似乎非常接近。对角线将包含聚类中心的差异：
np.diag(pairwise.pairwise_distances(kmeans.cluster_centers_, minibatch.cluster_centers_))
array([ 0.07328909, 0.09072807, 0.06571599])
工作原理…
这里的批次很关键。通过批次迭代以找到批次均值；对于下一次迭代，相对于当前迭代更新前一批次均值。有几个选项决定了一般 k-means 的行为和参数，这些参数决定了 MiniBatch k-means 的更新方式。
batch_size参数决定了批处理的大小。仅供娱乐，让我们来运行 MiniBatch；但是这次我们将批处理大小设置为与数据集大小相同：
minibatch = MiniBatchKMeans(batch_size=len(blobs))
%time minibatch.fit(blobs)
Wall time: 1min
MiniBatchKMeans(batch_size=1000000, compute_labels=True, init='k-means++',
init_size=None, max_iter=100, max_no_improvement=10, n_clusters=8,
n_init=3, random_state=None, reassignment_ratio=0.01, tol=0.0,
verbose=0)
显然，这违背了问题的本意，但它确实说明了一个重要的观点。选择不良的初始条件可能会影响模型的收敛效果，特别是聚类模型。对于 MiniBatch k-means，不能保证会达到全局最优解。
MiniBatch k-means 中有许多有力的教训。它利用了许多随机样本的力量，类似于自助法。在创建大数据的算法时，您可以在许多机器上并行处理许多随机样本。
使用 k-means 聚类量化图像
图像处理是一个重要的话题，聚类在其中有一些应用。值得指出的是，Python 中有几个非常好的图像处理库。scikit-image是 scikit-learn 的姊妹项目。如果您想要做一些复杂的事情，不妨看看它。
本章的重点之一是图像也是数据，聚类可以用来尝试猜测图像中某些物体的位置。聚类可以是图像处理流程的一部分。
准备工作
在这个示例中我们将会有一些乐趣。目标是使用一个聚类方法来对图像进行模糊处理。首先，我们将使用 SciPy 来读取图像。图像被转换为一个三维数组；x和y坐标描述了高度和宽度，第三维表示每个像素的 RGB 值。
首先，下载或将一个 .jpg 图像移动到你的 IPython 笔记本所在的文件夹。你可以使用你自己的照片。我使用的是一张名为 headshot.jpg 的自己的照片。
怎么做……
- 现在，让我们用 Python 读取图像：
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import ndimage
img = ndimage.imread("headshot.jpg")
plt.figure(figsize = (10,7))
plt.imshow(img)
看到以下图片：
- 就是我！现在我们有了图像，让我们检查一下它的维度：
img.shape
(379L, 337L, 3L)
要实际对图像进行量化，我们需要将其转换为二维数组，长度为 379 x 337，宽度为 RGB 值。更好的理解方式是将其看作是三维空间中的一组数据点，接着对这些点进行聚类，以减少图像中远离的颜色——这是一种简单的量化方法。
- 首先，让我们重新整理一下我们的数组；它是一个 NumPy 数组，因此操作起来很简单：
x, y, z = img.shape
long_img = img.reshape(x*y, z)
long_img.shape
(127723L, 3L)
- 现在我们可以开始聚类过程。首先，让我们导入聚类模块并创建一个 k-means 对象。我们将传递
n_clusters=5，这样我们就有五个聚类，或者说，五种不同的颜色。这将是一个很好的练习，帮助我们使用轮廓距离，这个在优化质心数量的练习中已经讲解过：
from sklearn import cluster
k_means = cluster.KMeans(n_clusters=5)
k_means.fit(long_img)
centers = k_means.cluster_centers_
centers
array([[ 169.01964615, 123.08399844, 99.6097561 ],
[ 45.79271071, 94.56844879, 120.00911162],
[ 218.74043562, 202.152748 , 184.14355039],
[ 67.51082485, 151.50671141, 201.9408963 ],
[ 169.69235986, 189.63274724, 143.75511521]])
它是如何工作的……
现在我们有了中心点，接下来我们需要的是标签。这将告诉我们哪些点应该与哪些聚类相关联：
labels = k_means.labels_
labels
array([4, 4, 4, ..., 3, 3, 3])
此时，我们只需进行最简单的 NumPy 数组操作，然后稍作调整，就可以得到新的图像：
plt.figure(figsize = (10,7))
plt.imshow(centers[labels].reshape(x, y, z))
以下是结果图像：
聚类将图像分成了几个区域。
在特征空间中寻找最接近的对象
有时，最简单的方法是找出两个对象之间的距离。我们只需要找出一个距离度量，计算成对距离，并将结果与预期的进行比较。
准备工作
scikit-learn 中的一个底层工具是 sklearn.metrics.pairwise。它包含了用于计算矩阵 X 中向量之间或 X 和 Y 中向量之间距离的服务器函数。这对于信息检索很有帮助。例如，给定一组具有属性 X 的客户，我们可能想选择一个参考客户，找到与该客户最相似的客户。
事实上，我们可能想根据相似度来对客户进行排名，而这种相似度是通过距离函数来衡量的。相似度的质量取决于特征空间的选择以及我们对空间进行的任何变换。我们将通过几种不同的距离衡量场景来进行讲解。
如何做……
我们将使用 pairwise_distances 函数来确定对象之间的相似性。请记住，相似性本质上就是通过我们定义的距离函数来衡量的：
- 首先，让我们从 metrics 模块中导入成对距离函数，并创建一个数据集进行操作：
import numpy as np
from sklearn.metrics import pairwise
from sklearn.datasets import make_blobs
points, labels = make_blobs()
- 检查距离的最简单方法是
pairwise_distances：
distances = pairwise.pairwise_distances(points)
distances是一个 N x N 矩阵，对角线上的值为 0。在最简单的情况下，让我们查看每个点与第一个点之间的距离：
np.diag(distances) [:5]
distances[0][:5]
array([ 0\. , 4.24926332, 8.8630893 , 5.01378992, 10.05620093])
- 按照接近程度对点进行排名非常简单，使用 np.argsort：
ranks = np.argsort(distances[0])
ranks[:5]
array([ 0, 63, 6, 21, 17], dtype=int64)
argsort的一个优点是，现在我们可以对points矩阵进行排序，从而获得实际的点：
points[ranks][:5]
array([[-0.15728042, -5.76309092],
[-0.20720885, -5.52734277],
[-0.08686778, -6.42054076],
[ 0.33493582, -6.29824601],
[-0.89842683, -5.78335127]])
sp_points = points[ranks][:5]
- 看到最接近的点是什么样子是很有用的，如下所示。选定的点[0]被涂成绿色。最接近的点被涂成红色（除了选定的点）。
请注意，除了某些保证之外，这个过程按预期工作：
import matplotlib.pyplot as plt
%matplotlib inline
plt.figure(figsize=(10,7))
plt.scatter(points[:,0], points[:,1], label = 'All Points')
plt.scatter(sp_points[:,0],sp_points[:,1],color='red', label='Closest Points')
plt.scatter(points[0,0],points[0,1],color='green', label = 'Chosen Point')
plt.legend()
它是如何工作的…
给定某个距离函数，每个点通过成对的函数进行测量。考虑两个点，表示为 N 维空间中的向量，具有分量p[i]和q[i]；默认情况下是欧几里得距离，其公式如下：
口头上，这个过程是计算两个向量每个分量之间的差值，平方这些差异，求和所有差异，然后求平方根。这个过程看起来非常熟悉，因为我们在研究均方误差时用过类似的东西。如果我们取平方根，我们就得到了相同的结果。事实上，常用的一个度量是均方根偏差（RMSE），它就是应用的距离函数。
在 Python 中，这看起来如下所示：
def euclid_distances(x, y):
return np.power(np.power(x - y, 2).sum(), .5)
euclid_distances(points[0], points[1])
4.249263322917467
在 scikit-learn 中还有其他一些函数，但 scikit-learn 还将使用 SciPy 的距离函数。在写这本书时，scikit-learn 的距离函数支持稀疏矩阵。有关距离函数的更多信息，请查阅 SciPy 文档：
-
cityblock -
cosine -
euclidean -
l1 -
l2 -
manhattan
我们现在可以解决问题了。例如，如果我们站在网格的原点上，而这些线是街道，我们需要走多远才能到达点*(5, 5)*？
pairwise.pairwise_distances([[0, 0], [5, 5]], metric='cityblock')[0]
array([ 0., 10.])
还有更多…
使用成对距离，我们可以找到比特向量之间的相似度。对于 N 维向量p和q，问题就变成了计算汉明距离，其定义如下：
使用以下命令：
X = np.random.binomial(1, .5, size=(2, 4)).astype(np.bool)
X
array([[False, False, False, False],
[False, True, True, True]], dtype=bool)
pairwise.pairwise_distances(X, metric='hamming')
array([[ 0\. , 0.75],
[ 0.75, 0\. ]])
请注意，scikit-learn 的hamming度量返回的是汉明距离除以向量的长度，在此情况下为4。
使用高斯混合模型进行概率聚类
在 k-means 中，我们假设簇的方差是相等的。这导致了一种空间的划分，决定了簇的分配方式；但是，如果方差不相等，并且每个簇点与其有某种概率关联，应该如何处理呢？
准备工作
有一种更具概率性的方式来看待 k-means 聚类。硬性 k-means 聚类等同于应用一个具有协方差矩阵 S 的高斯混合模型（GMM），其中协方差矩阵 S 可以分解为误差和单位矩阵的乘积。这对于每个聚类来说都是相同的协方差结构，导致聚类呈球形。然而，如果我们允许 S 变化，则可以估计 GMM 并用于预测。我们将首先在一维空间中看看这种方法如何运作，然后扩展到更多维度。
如何操作…
- 首先，我们需要创建一些数据。例如，我们可以模拟女性和男性的身高。在整个例子中，我们将使用这个例子。它是一个简单的示例，但希望它能够展示我们在 N 维空间中试图实现的目标，这样稍微更容易进行可视化：
import numpy as np
N = 1000
in_m = 72
in_w = 66
s_m = 2
s_w = s_m
m = np.random.normal(in_m, s_m, N)
w = np.random.normal(in_w, s_w, N)
from matplotlib import pyplot as plt
%matplotlib inline
f, ax = plt.subplots(figsize=(7, 5))
ax.set_title("Histogram of Heights")
ax.hist(m, alpha=.5, label="Men");
ax.hist(w, alpha=.5, label="Women");
ax.legend()
这是输出结果：
- 接下来，我们可能有兴趣对子集进行抽样，拟合分布，然后预测剩余的组：
random_sample = np.random.choice([True, False], size=m.size)
m_test = m[random_sample]
m_train = m[~random_sample]
w_test = w[random_sample]
w_train = w[~random_sample]
- 现在，我们需要根据训练集获取男性和女性身高的经验分布：
from scipy import stats
m_pdf = stats.norm(m_train.mean(), m_train.std())
w_pdf = stats.norm(w_train.mean(), w_train.std())
对于测试集，我们将根据数据点来自某个分布的似然度进行计算，最有可能的分布将被分配到相应的标签。
- 我们当然会查看我们的准确度：
m_pdf.pdf(m[0])
0.19762291119664221
w_pdf.pdf(m[0])
0.00085042279862613103
- 注意似然度的差异。假设我们在男性概率较高时进行猜测，但如果女性的概率更高时，我们会覆盖这些猜测：
guesses_m = np.ones_like(m_test)
guesses_m[m_pdf.pdf(m_test) < w_pdf.pdf(m_test)] = 0
- 显然，问题是我们有多准确。由于如果我们猜对了，
guesses_m将是 1，如果猜错了，则为 0，我们可以取向量的均值来计算准确度：
guesses_m.mean()
0.94176706827309242
- 还不错！现在，为了查看我们在女性组上的表现，我们使用以下命令：
guesses_w = np.ones_like(w_test)
guesses_w[m_pdf.pdf(w_test) > w_pdf.pdf(w_test)] = 0
guesses_w.mean()
0.93775100401606426
- 让我们允许不同组之间的方差有所不同。首先，创建一些新的数据：
s_m = 1
s_w = 4
m = np.random.normal(in_m, s_m, N)
w = np.random.normal(in_w, s_w, N)
- 然后，创建一个训练集：
m_test = m[random_sample]
m_train = m[~random_sample]
w_test = w[random_sample]
w_train = w[~random_sample]
f, ax = plt.subplots(figsize=(7, 5))
ax.set_title("Histogram of Heights")
ax.hist(m_train, alpha=.5, label="Men");
ax.hist(w_train, alpha=.5, label="Women");
ax.legend()
- 现在我们可以创建相同的概率密度函数（PDFs）：
m_pdf = stats.norm(m_train.mean(), m_train.std())
w_pdf = stats.norm(w_train.mean(), w_train.std())
x = np.linspace(50,80,300)
plt.figure(figsize=(8,5))
plt.title('PDF of Heights')
plt.plot(x, m_pdf.pdf(x), 'k', linewidth=2, color='blue', label='Men')
plt.plot(x, w_pdf.pdf(x), 'k', linewidth=2, color='green',label='Women')
- 以下是输出结果：
你可以在一个多维空间中想象这个过程：
class_A = np.random.normal(0, 1, size=(100, 2))
class_B = np.random.normal(4, 1.5, size=(100, 2))
f, ax = plt.subplots(figsize=(8, 5))
plt.title('Random 2D Normal Draws')
ax.scatter(class_A[:,0], class_A[:,1], label='A', c='r')
ax.scatter(class_B[:,0], class_B[:,1], label='B')
它是如何工作的…
好的，现在我们已经看过如何根据分布来分类数据点，接下来我们看看如何在 scikit-learn 中实现这一点：
from sklearn.mixture import GaussianMixture
gmm = GaussianMixture(n_components=2)
X = np.row_stack((class_A, class_B))
y = np.hstack((np.ones(100), np.zeros(100)))
由于我们是认真的数据科学家，我们将创建一个训练集：
train = np.random.choice([True, False], 200)
gmm.fit(X[train])
GaussianMixture(covariance_type='full', init_params='kmeans', max_iter=100,
means_init=None, n_components=2, n_init=1, precisions_init=None,
random_state=None, reg_covar=1e-06, tol=0.001, verbose=0,
verbose_interval=10, warm_start=False, weights_init=None)
拟合和预测的方式与在 scikit-learn 中拟合其他对象时相同：
gmm.fit(X[train])
gmm.predict(X[train])[:5]
array([0, 0, 0, 0, 0], dtype=int64)
在模型拟合之后，还有其他方法值得关注。例如，使用 score_samples，我们可以获得每个标签的每个样本的似然度。
使用 k-means 进行异常值检测
在这个例子中，我们将讨论 k-means 在异常值检测中的辩论和机制。它在隔离某些类型的错误时可能有用，但使用时需要小心。
准备工作
我们将使用 k-means 对一组点进行异常值检测。值得注意的是，在异常值和异常值检测方面有许多不同的观点。一方面，我们可能通过移除异常值去除了一些由数据生成过程产生的点；另一方面，异常值可能是由于测量误差或其他外部因素导致的。
这是我们对这个辩论给出的最大信任。接下来的部分将关注识别异常值；我们将假设移除异常值的选择是合理的。异常值检测的过程是找到簇的质心，然后通过计算与质心的距离来识别潜在的异常值点。
如何操作……
- 首先，我们将生成一个包含 100 个点的单一簇，然后识别出距离质心最远的五个点。这些点可能是异常值：
from sklearn.datasets import make_blobs
X, labels = make_blobs(100, centers=1)
import numpy as np
- k-means 聚类必须具有一个单一的质心，这是很重要的。这个概念类似于用于异常值检测的单类 SVM：
from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=1)
kmeans.fit(X)
KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
n_clusters=1, n_init=10, n_jobs=1, precompute_distances='auto',
random_state=None, tol=0.0001, verbose=0)
- 现在，让我们看看这个图。正在家里一起操作的朋友们，试着猜猜哪些点会被识别为五个异常值之一：
import matplotlib.pyplot as plt
%matplotlib inline
f, ax = plt.subplots(figsize=(8, 5))
ax.set_title("Blob")
ax.scatter(X[:, 0], X[:, 1], label='Points')
ax.scatter(kmeans.cluster_centers_[:, 0],kmeans.cluster_centers_[:, 1], label='Centroid',color='r')
ax.legend()
以下是输出结果：
- 现在，让我们识别出最接近的五个点：
distances = kmeans.transform(X)
# argsort returns an array of indexes which will sort the array in ascending order
# so we reverse it via [::-1] and take the top five with [:5]
sorted_idx = np.argsort(distances.ravel())[::-1][:5]
- 让我们看看哪些点离得最远：
f, ax = plt.subplots(figsize=(7, 5))
ax.set_title("Single Cluster")
ax.scatter(X[:, 0], X[:, 1], label='Points')
ax.scatter(kmeans.cluster_centers_[:, 0],kmeans.cluster_centers_[:, 1],label='Centroid', color='r')
ax.scatter(X[sorted_idx][:, 0], X[sorted_idx][:, 1],label='Extreme Value', edgecolors='g',facecolors='none', s=100)
ax.legend(loc='best')
以下是输出结果：
- 如果我们愿意，去除这些点是很简单的：
new_X = np.delete(X, sorted_idx, axis=0)
此外，随着这些点的移除，质心明显发生了变化：
new_kmeans = KMeans(n_clusters=1)
new_kmeans.fit(new_X)
- 让我们可视化一下旧质心和新质心之间的差异：
f, ax = plt.subplots(figsize=(7, 5))
ax.set_title("Extreme Values Removed")
ax.scatter(new_X[:, 0], new_X[:, 1], label='Pruned Points')
ax.scatter(kmeans.cluster_centers_[:, 0],kmeans.cluster_centers_[:, 1], label='Old Centroid',color='r',s=80, alpha=.5)
ax.scatter(new_kmeans.cluster_centers_[:, 0],new_kmeans.cluster_centers_[:, 1], label='New Centroid',color='m', s=80, alpha=.5)
ax.legend(loc='best')
以下是输出结果：
很明显，质心几乎没有移动，这在移除五个最极端的值时是可以预期的。这个过程可以重复进行，直到我们满意数据能代表这个过程。
它是如何工作的……
正如我们已经看到的，高斯分布和 k-means 聚类之间有着根本的联系。让我们基于质心和样本协方差矩阵创建一个经验高斯分布，并查看每个点的概率——理论上，这就是我们移除的那五个点。这实际上表明我们已经移除了最不可能出现的值。距离和可能性之间的这种关系非常重要，在你的机器学习训练中会经常遇到。使用以下命令创建经验高斯分布：
from scipy import stats
emp_dist = stats.multivariate_normal(kmeans.cluster_centers_.ravel())
lowest_prob_idx = np.argsort(emp_dist.pdf(X))[:5]
np.all(X[sorted_idx] == X[lowest_prob_idx])
True
使用 KNN 进行回归
回归在本书的其他部分已经涉及，但我们也可能想对特征空间中的小范围进行回归。我们可以认为我们的数据集受到多个数据处理过程的影响。如果这是真的，只有在相似的数据点上进行训练才是明智的选择。
准备工作
我们的老朋友——回归，可以在聚类的背景下使用。回归显然是一种监督学习技术，因此我们将使用K-最近邻（KNN）聚类，而不是 K 均值聚类。对于 KNN 回归，我们将使用特征空间中离测试点最近的 K 个点来构建回归模型，而不是像常规回归那样使用整个特征空间。
如何实现…
对于这个示例，我们将使用iris数据集。如果我们想预测每朵花的花瓣宽度，通过鸢尾花种类进行聚类可能会为我们带来更好的结果。KNN 回归不会根据种类进行聚类，但我们会假设相同种类的样本在特征空间中会较为接近，在这个例子中，即花瓣长度：
- 我们将使用
iris数据集进行这个示例：
import numpy as np
from sklearn import datasets
iris = datasets.load_iris()
iris.feature_names
- 我们将尝试基于萼片的长度和宽度预测花瓣长度。我们还将拟合一个常规的线性回归，看看 KNN 回归与它相比表现如何：
X = iris.data[:,:2]
y = iris.data[:,2]
from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(X, y)
print "The MSE is: {:.2}".format(np.power(y - lr.predict(X),2).mean())
The MSE is: 0.41
- 现在，对于 KNN 回归，使用以下代码：
from sklearn.neighbors import KNeighborsRegressor
knnr = KNeighborsRegressor(n_neighbors=10)
knnr.fit(X, y)
print "The MSE is: {:.2}".format(np.power(y - knnr.predict(X),2).mean())
The MSE is: 0.17
- 让我们看看当我们让 KNN 回归使用离测试点最近的 10 个点时，结果会是怎样：
f, ax = plt.subplots(nrows=2, figsize=(7, 10))
ax[0].set_title("Predictions")
ax[0].scatter(X[:, 0], X[:, 1], s=lr.predict(X)*80, label='LR Predictions', color='c', edgecolors='black')
ax[1].scatter(X[:, 0], X[:, 1], s=knnr.predict(X)*80, label='k-NN Predictions', color='m', edgecolors='black')
ax[0].legend()
ax[1].legend()
- 预测结果大部分是接近的，这可能是显而易见的，但让我们看看 Setosa 种类的预测结果与实际值的对比：
setosa_idx = np.where(iris.target_names=='setosa')
setosa_mask = iris.target == setosa_idx[0]
y[setosa_mask][:5]
array([ 1.4, 1.4, 1.3, 1.5, 1.4])
knnr.predict(X)[setosa_mask][:5]
array([ 1.46, 1.47, 1.51, 1.42, 1.48])
lr.predict(X)[setosa_mask][:5]
array([ 1.83762646, 2.1510849 , 1.52707371, 1.48291658, 1.52562087])
- 再次查看图表，我们看到 Setosa 种类（左上角的聚类）在常规线性回归中被大大高估，而 KNN 则非常接近实际值。
它是如何工作的…
KNN 回归非常简单，通过取离测试点最近的K个点的平均值来计算。让我们手动预测一个点：
example_point = X[0]
现在，我们需要找到离我们的example_point最近的 10 个点：
from sklearn.metrics import pairwise
distances_to_example = pairwise.pairwise_distances(X)[0]
ten_closest_points = X[np.argsort(distances_to_example)][:10]
ten_closest_y = y[np.argsort(distances_to_example)][:10]
ten_closest_y.mean()
1.46
我们可以看到，结果非常接近预期。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
