机器学习口袋参考（四）
原文：
annas-archive.org/md5/dd771df8c19ec4613e3638bc1f862b92译者：飞龙
第十六章：解释回归模型
大多数用于解释分类模型的技术同样适用于回归模型。在本章中，我将展示如何使用 SHAP 库解释回归模型。
我们将解释波士顿房屋数据集的 XGBoost 模型：
>>> import xgboost as xgb
>>> xgr = xgb.XGBRegressor(
... random_state=42, base_score=0.5
... )
>>> xgr.fit(bos_X_train, bos_y_train)
Shapley
我非常喜欢 Shapley 因为它对模型不可知。这个库还为我们提供了全局对模型的洞察，并帮助解释单个预测。如果你有一个黑盒模型，我发现它非常有用。
我们首先来看看索引 5 的预测。我们的模型预测值为 27.26：
>>> sample_idx = 5
>>> xgr.predict(bos_X.iloc[[sample_idx]])
array([27.269186], dtype=float32)
要使用模型，我们必须从我们的模型创建一个TreeExplainer，并估算我们样本的 SHAP 值。如果我们想在 Jupyter 上使用交互界面，我们还需要调用initjs函数：
>>> import shap
>>> shap.initjs()
>>> exp = shap.TreeExplainer(xgr)
>>> vals = exp.shap_values(bos_X)
有了解释器和 SHAP 值，我们可以创建一个力图来解释预测（见图 16-1）。这告诉我们基础预测值为 23，人口状态（LSTAT）和财产税率（TAX）将价格推高，而房间数（RM）将价格推低：
>>> shap.force_plot(
... exp.expected_value,
... vals[sample_idx],
... bos_X.iloc[sample_idx],
... )
图 16-1. 回归的力图。由于人口状态和税率，预期值从 23 增至 27。
我们也可以查看所有样本的力图，以获得整体行为的感觉。如果我们在 Jupyter 上使用交互式 JavaScript 模式，我们可以将鼠标悬停在样本上，查看影响结果的特征（见图 16-2）：
>>> shap.force_plot(
... exp.expected_value, vals, bos_X
... )
图 16-2. 所有样本的回归力图。
从样本的力图中，我们看到 LSTAT 特征有很大的影响。为了可视化 LSTAT 如何影响结果，我们可以创建一个依赖图。库将自动选择一个特征进行着色（您可以提供interaction_index参数来设置您自己的）。
从 LSTAT 的依赖图中（见图 16-3），我们可以看到随着 LSTAT 的增加（低社会地位人口的百分比），SHAP 值下降（推动目标向下）。非常低的 LSTAT 值会提升 SHAP。通过查看 TAX（财产税率）的颜色，我们可以看到随着税率的降低（更蓝色），SHAP 值上升：
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> shap.dependence_plot("LSTAT", vals, bos_X)
>>> fig.savefig(
... "images/mlpr_1603.png",
... bbox_inches="tight",
... dpi=300,
... )
图 16-3. LSTAT 的依赖图。随着 LSTAT 的增加，预测值下降。
这里是另一个依赖图，在图 16-4 中展示了 DIS（到就业中心的距离）。看起来这个特征的影响很小，除非它非常小：
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> shap.dependence_plot(
... "DIS", vals, bos_X, interaction_index="RM"
... )
>>> fig.savefig(
... "images/mlpr_1604.png",
... bbox_inches="tight",
... dpi=300,
... )
图 16-4. DIS 的依赖图。除非 DIS 非常小，否则 SHAP 保持相对平稳。
最后，我们将使用总结图来查看特征的全局效果（见 Figure 16-5）。顶部的特征对模型影响最大。从这个视角可以看出，RM（房间数）的大值显著提升了目标值，而中等和较小值则略微降低了它：
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> shap.summary_plot(vals, bos_X)
>>> fig.savefig(
... "images/mlpr_1605.png",
... bbox_inches="tight",
... dpi=300,
... )
图 16-5. 总结图。最重要的特征位于顶部。
SHAP 库是您工具箱中的一个很好的工具。它帮助理解特征的全局影响，同时也有助于解释单个预测。
第十七章：降维
有许多技术将特征分解为较小的子集。这对于探索性数据分析、可视化、制作预测模型或聚类很有用。
在本章中，我们将使用各种技术探索泰坦尼克号数据集。我们将查看 PCA、UMAP、t-SNE 和 PHATE。
这是数据：
>>> ti_df = tweak_titanic(orig_df)
>>> std_cols = "pclass,age,sibsp,fare".split(",")
>>> X_train, X_test, y_train, y_test = get_train_test_X_y(
... ti_df, "survived", std_cols=std_cols
... )
>>> X = pd.concat([X_train, X_test])
>>> y = pd.concat([y_train, y_test])
主成分分析（PCA）
主成分分析（PCA）接受一个行（样本）和列（特征）的矩阵（X）。PCA 返回一个新的矩阵，其列是原始列的线性组合。这些线性组合最大化方差。
每列与其他列正交（成直角）。列按方差递减的顺序排序。
Scikit-learn 有这个模型的实现。在运行算法之前最好标准化数据。调用.fit方法后，您将可以访问一个.explained_variance_ratio_属性，列出每列中方差的百分比。
PCA 对于在二维（或三维）中可视化数据很有用。它还用作预处理步骤，以过滤数据中的随机噪声。它适合于找到全局结构，但不适合找到局部结构，并且在处理线性数据时表现良好。
在这个例子中，我们将在泰坦尼克号特征上运行 PCA。PCA 类是 scikit-learn 中的一个转换器；您调用.fit方法来教它如何获取主成分，然后调用.transform将矩阵转换为主成分矩阵：
>>> from sklearn.decomposition import PCA
>>> from sklearn.preprocessing import (
... StandardScaler,
... )
>>> pca = PCA(random_state=42)
>>> X_pca = pca.fit_transform(
... StandardScaler().fit_transform(X)
... )
>>> pca.explained_variance_ratio_
array([0.23917891, 0.21623078, 0.19265028,
0.10460882, 0.08170342, 0.07229959,
0.05133752, 0.04199068])
>>> pca.components_[0]
arrayarray([-0.63368693, 0.39682566,
0.00614498, 0.11488415, 0.58075352,
-0.19046812, -0.21190808, -0.09631388])
实例参数：
n_components=None
生成的组件数量。如果为None，则返回与列数相同的组件数量。可以是一个浮点数（0, 1），那么将创建所需数量的组件以获得该方差比例。
copy=True
如果为True，将在.fit时改变数据。
whiten=False
转换后的白化数据以确保无关联的组件。
svd_solver='auto'
如果n_components小于最小维度的 80％，则运行'randomized' SVD（更快但是近似）。否则运行'full'。
tol=0.0
对奇异值的容差。
iterated_power='auto'
'randomized' svd_solver的迭代次数。
random_state=None
'randomized' svd_solver的随机状态。
属性：
components_
主成分（原始特征的线性组合权重列）。
explained_variance_
每个组件的方差量。
explained_variance_ratio_
每个组件的方差量归一化（总和为 1）。
singular_values_
每个组件的奇异值。
mean_
每个特征的均值。
n_components_
当n_components是浮点数时，这是组件的大小。
noise_variance_
估计的噪声协方差。
绘制解释方差比例的累积和被称为屏幕图（见图 17-1）。它将显示组件中存储了多少信息。您可以使用肘方法来确定使用多少个组件：
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> ax.plot(pca.explained_variance_ratio_)
>>> ax.set(
... xlabel="Component",
... ylabel="Percent of Explained variance",
... title="Scree Plot",
... ylim=(0, 1),
... )
>>> fig.savefig(
... "images/mlpr_1701.png",
... dpi=300,
... bbox_inches="tight",
... )
图 17-1. PCA 屏幕图。
查看数据的另一种方法是使用累计图（见图 17-2）。我们的原始数据有 8 列，但从图中看来，如果仅使用 4 个 PCA 成分，我们可以保留大约 90% 的方差：
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> ax.plot(
... np.cumsum(pca.explained_variance_ratio_)
... )
>>> ax.set(
... xlabel="Component",
... ylabel="Percent of Explained variance",
... title="Cumulative Variance",
... ylim=(0, 1),
... )
>>> fig.savefig("images/mlpr_1702.png", dpi=300)
图 17-2. PCA 累计解释方差。
特征对成分的影响有多大？使用 matplotlib 的 imshow 函数将成分沿 x 轴和原始特征沿 y 轴绘制出来（见图 17-3）。颜色越深，原始列对成分的贡献越大。
看起来第一个成分受到 pclass、age 和 fare 列的影响很大。（使用光谱色图（cmap）强调非零值，并提供 vmin 和 vmax 为色条图例添加限制。）
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> plt.imshow(
... pca.components_.T,
... cmap="Spectral",
... vmin=-1,
... vmax=1,
... )
>>> plt.yticks(range(len(X.columns)), X.columns)
>>> plt.xticks(range(8), range(1, 9))
>>> plt.xlabel("Principal Component")
>>> plt.ylabel("Contribution")
>>> plt.title(
... "Contribution of Features to Components"
... )
>>> plt.colorbar()
>>> fig.savefig("images/mlpr_1703.png", dpi=300)
图 17-3. 主成分分析中的 PCA 特征。
另一种查看数据的方式是使用条形图（见图 17-4）。每个成分显示了来自原始数据的贡献：
>>> fig, ax = plt.subplots(figsize=(8, 4))
>>> pd.DataFrame(
... pca.components_, columns=X.columns
... ).plot(kind="bar", ax=ax).legend(
... bbox_to_anchor=(1, 1)
... )
>>> fig.savefig("images/mlpr_1704.png", dpi=300)
图 17-4. 主成分分析中的 PCA 特征。
如果我们有很多特征，可能希望通过仅显示满足最小权重要求的特征来限制上述图形。以下是找出前两个成分中具有至少 0.5 绝对值的所有特征的代码：
>>> comps = pd.DataFrame(
... pca.components_, columns=X.columns
... )
>>> min_val = 0.5
>>> num_components = 2
>>> pca_cols = set()
>>> for i in range(num_components):
... parts = comps.iloc[i][
... comps.iloc[i].abs() > min_val
... ]
... pca_cols.update(set(parts.index))
>>> pca_cols
{'fare', 'parch', 'pclass', 'sibsp'}
PCA 常用于以两个成分可视化高维数据集。在这里，我们在 2D 中可视化了 Titanic 的特征。它们根据生存状态进行了着色。有时可视化中可能会出现聚类。在这种情况下，似乎没有幸存者的聚类现象（见图 17-5）。
我们使用 Yellowbrick 生成此可视化：
>>> from yellowbrick.features.pca import (
... PCADecomposition,
... )
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> colors = ["rg"[j] for j in y]
>>> pca_viz = PCADecomposition(color=colors)
>>> pca_viz.fit_transform(X, y)
>>> pca_viz.poof()
>>> fig.savefig("images/mlpr_1705.png", dpi=300)
图 17-5. Yellowbrick PCA 绘图。
如果要根据一列着色散点图并添加图例（而不是色条图），则需要在 pandas 或 matplotlib 中循环每个颜色并单独绘制该组（或使用 seaborn）。下面我们还将纵横比设置为我们查看的成分解释方差的比率（见图 17-6）。因为第二个成分仅具有第一个成分的 90%，所以它会稍微短一些。
这是 seaborn 版本：
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> pca_df = pd.DataFrame(
... X_pca,
... columns=[
... f"PC{i+1}"
... for i in range(X_pca.shape[1])
... ],
... )
>>> pca_df["status"] = [
... ("deceased", "survived")[i] for i in y
... ]
>>> evr = pca.explained_variance_ratio_
>>> ax.set_aspect(evr[1] / evr[0])
>>> sns.scatterplot(
... x="PC1",
... y="PC2",
... hue="status",
... data=pca_df,
... alpha=0.5,
... ax=ax,
... )
>>> fig.savefig(
... "images/mlpr_1706.png",
... dpi=300,
... bbox_inches="tight",
... )
图 17-6. 带有图例和相对纵横比的 Seaborn PCA。
下面，我们通过在散点图上显示一个载荷图来增强可视化效果。这种图被称为双标图，因为它包含散点图和载荷（见图 17-7）。载荷显示特征的强度和它们的相关性。如果它们的角度接近，则它们可能相关。如果角度为 90 度，则它们可能不相关。最后，如果它们之间的角度接近 180 度，则它们具有负相关性：
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> pca_df = pd.DataFrame(
... X_pca,
... columns=[
... f"PC{i+1}"
... for i in range(X_pca.shape[1])
... ],
... )
>>> pca_df["status"] = [
... ("deceased", "survived")[i] for i in y
... ]
>>> evr = pca.explained_variance_ratio_
>>> x_idx = 0 # x_pc
>>> y_idx = 1 # y_pc
>>> ax.set_aspect(evr[y_idx] / evr[x_idx])
>>> x_col = pca_df.columns[x_idx]
>>> y_col = pca_df.columns[y_idx]
>>> sns.scatterplot(
... x=x_col,
... y=y_col,
... hue="status",
... data=pca_df,
... alpha=0.5,
... ax=ax,
... )
>>> scale = 8
>>> comps = pd.DataFrame(
... pca.components_, columns=X.columns
... )
>>> for idx, s in comps.T.iterrows():
... plt.arrow(
... 0,
... 0,
... s[x_idx] * scale,
... s[y_idx] * scale,
... color="k",
... )
... plt.text(
... s[x_idx] * scale,
... s[y_idx] * scale,
... idx,
... weight="bold",
... )
>>> fig.savefig(
... "images/mlpr_1707.png",
... dpi=300,
... bbox_inches="tight",
... )
图 17-7. Seaborn 双图绘图和加载图。
根据先前的树模型，我们知道年龄、票价和性别对乘客是否生存很重要。第一个主成分受 pclass、年龄和票价的影响，而第四个主成分受性别的影响。让我们将这些组件相互绘制。
同样，这个图是根据组件方差比例调整了绘图的纵横比（见图 17-8）。
此图似乎更准确地区分了幸存者：
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> pca_df = pd.DataFrame(
... X_pca,
... columns=[
... f"PC{i+1}"
... for i in range(X_pca.shape[1])
... ],
... )
>>> pca_df["status"] = [
... ("deceased", "survived")[i] for i in y
... ]
>>> evr = pca.explained_variance_ratio_
>>> ax.set_aspect(evr[3] / evr[0])
>>> sns.scatterplot(
... x="PC1",
... y="PC4",
... hue="status",
... data=pca_df,
... alpha=0.5,
... ax=ax,
... )
>>> fig.savefig(
... "images/mlpr_1708.png",
... dpi=300,
... bbox_inches="tight",
... )
图 17-8. PCA 组件 1 与 4 的图示。
Matplotlib 可以创建漂亮的图表，但在交互图上不太实用。在执行 PCA 时，查看散点图通常很有用。我已经包含了一个使用 Bokeh 库 的函数，用于与散点图进行交互（见图 17-9）。它在 Jupyter 中运行良好：
>>> from bokeh.io import output_notebook
>>> from bokeh import models, palettes, transform
>>> from bokeh.plotting import figure, show
>>>
>>> def bokeh_scatter(
... x,
... y,
... data,
... hue=None,
... label_cols=None,
... size=None,
... legend=None,
... alpha=0.5,
... ):
... """
... x - x column name to plot
... y - y column name to plot
... data - pandas DataFrame
... hue - column name to color by (numeric)
... legend - column name to label by
... label_cols - columns to use in tooltip
... (None all in DataFrame)
... size - size of points in screen space unigs
... alpha - transparency
... """
... output_notebook()
... circle_kwargs = {}
... if legend:
... circle_kwargs["legend"] = legend
... if size:
... circle_kwargs["size"] = size
... if hue:
... color_seq = data[hue]
... mapper = models.LinearColorMapper(
... palette=palettes.viridis(256),
... low=min(color_seq),
... high=max(color_seq),
... )
... circle_kwargs[
... "fill_color"
... ] = transform.transform(hue, mapper)
... ds = models.ColumnDataSource(data)
... if label_cols is None:
... label_cols = data.columns
... tool_tips = sorted(
... [
... (x, "@{}".format(x))
... for x in label_cols
... ],
... key=lambda tup: tup[0],
... )
... hover = models.HoverTool(
... tooltips=tool_tips
... )
... fig = figure(
... tools=[
... hover,
... "pan",
... "zoom_in",
... "zoom_out",
... "reset",
... ],
... toolbar_location="below",
... )
...
... fig.circle(
... x,
... y,
... source=ds,
... alpha=alpha,
... **circle_kwargs
... )
... show(fig)
... return fig
>>> res = bokeh_scatter(
... "PC1",
... "PC2",
... data=pca_df.assign(
... surv=y.reset_index(drop=True)
... ),
... hue="surv",
... size=10,
... legend="surv",
... )
图 17-9. 带有工具提示的 Bokeh 散点图。
Yellowbrick 也可以在三维中绘图（见图 17-10）：
>>> from yellowbrick.features.pca import (
... PCADecomposition,
... )
>>> colors = ["rg"[j] for j in y]
>>> pca3_viz = PCADecomposition(
... proj_dim=3, color=colors
... )
>>> pca3_viz.fit_transform(X, y)
>>> pca3_viz.finalize()
>>> fig = plt.gcf()
>>> plt.tight_layout()
>>> fig.savefig(
... "images/mlpr_1710.png",
... dpi=300,
... bbox_inches="tight",
... )
图 17-10. Yellowbrick 3D PCA。
scprep 库（这是 PHATE 库的依赖项，我们稍后会讨论）具有一个有用的绘图函数。rotate_scatter3d 函数可以生成一个在 Jupyter 中动画显示的图形（见图 17-11）。这使得理解 3D 图形变得更容易。
您可以使用此库可视化任何 3D 数据，而不仅限于 PHATE：
>>> import scprep
>>> scprep.plot.rotate_scatter3d(
... X_pca[:, :3],
... c=y,
... cmap="Spectral",
... figsize=(8, 6),
... label_prefix="Principal Component",
... )
图 17-11. scprep 3D PCA 动画。
如果您在 Jupyter 中将 matplotlib 的单元格魔术模式更改为 notebook，您可以从 matplotlib 获取交互式 3D 绘图（见图 17-12）。
>>> from mpl_toolkits.mplot3d import Axes3D
>>> fig = plt.figure(figsize=(6, 4))
>>> ax = fig.add_subplot(111, projection="3d")
>>> ax.scatter(
... xs=X_pca[:, 0],
... ys=X_pca[:, 1],
... zs=X_pca[:, 2],
... c=y,
... cmap="viridis",
... )
>>> ax.set_xlabel("PC 1")
>>> ax.set_ylabel("PC 2")
>>> ax.set_zlabel("PC 3")
图 17-12. Matplotlib 在笔记本模式下的交互式 3D PCA。
警告
请注意，从：
% matplotlib inline
转到：
% matplotlib notebook
有时可能会导致 Jupyter 停止响应。请谨慎使用。
UMAP
统一流形逼近和投影 (UMAP) 是一种使用流形学习的降维技术。因此，它倾向于将相似的项目在拓扑上保持在一起。它试图同时保留全局和局部结构，与偏好局部结构的 t-SNE 相反（详见“t-SNE”）。
Python 实现不支持多核。
特征归一化是将值放在相同尺度上的一个好主意。
UMAP 对超参数（n_neighbors、min_dist、n_components 或 metric）非常敏感。以下是一些示例：
>>> import umap
>>> u = umap.UMAP(random_state=42)
>>> X_umap = u.fit_transform(
... StandardScaler().fit_transform(X)
... )
>>> X_umap.shape
(1309, 2)
实例参数：
n_neighbors=15
本地邻域大小。较大意味着使用全局视图，较小意味着更局部。
n_components=2
嵌入的维度数。
metric='euclidean'
用于距离的度量。可以是一个接受两个 1D 数组并返回浮点数的函数。
n_epochs=None
训练时期的数量。默认为 200 或 500（取决于数据大小）。
learning_rate=1.0
嵌入优化的学习率。
init='spectral'
初始化类型。谱嵌入是默认值。可以是'random'或 numpy 数组的位置。
min_dist=0.1
0 到 1 之间。嵌入点之间的最小距离。较小意味着更多聚集，较大意味着分散。
spread=1.0
确定嵌入点的距离。
set_op_mix_ratio=1.0
0 到 1 之间：模糊联合（1）或模糊交集（0）。
local_connectivity=1.0
本地连接的邻居数。增加此值会创建更多的局部连接。
repulsion_strength=1.0
排斥强度。较高的值给负样本更多权重。
negative_sample_rate=5
负样本每个正样本。更高的值具有更多的排斥力，更多的优化成本和更好的准确性。
transform_queue_size=4.0
用于最近邻搜索的侵略性。较高的值是较低的性能但更好的准确性。
a=None
参数来控制嵌入。如果等于None，UMAP 则从min_dist和spread中确定这些值。
b=None
参数来控制嵌入。如果等于None，UMAP 则从min_dist和spread中确定这些值。
random_state=None
随机种子。
metric_kwds=None
如果函数用于metric，则用于额外参数的指标字典。也可以使用minkowski（和其他指标）来参数化。
angular_rp_forest=False
使用角度随机投影。
target_n_neighbors=-1
简易设置的邻居数。
target_metric='categorical'
用于使用监督降维。也可以是'L1'或'L2'。还支持一个接受X中两个数组作为输入并返回它们之间距离值的函数。
target_metric_kwds=None
如果target_metric的函数被使用，使用的指标字典。
target_weight=0.5
权重因子。介于 0.0 和 1.0 之间，其中 0 表示仅基于数据，而 1 表示仅基于目标。
transform_seed=42
变换操作的随机种子。
verbose=False
冗余性。
属性：
embedding_
嵌入结果
让我们可视化泰坦尼克号数据集上 UMAP 的默认结果（参见图 17-13）：
>>> fig, ax = plt.subplots(figsize=(8, 4))
>>> pd.DataFrame(X_umap).plot(
... kind="scatter",
... x=0,
... y=1,
... ax=ax,
... c=y,
... alpha=0.2,
... cmap="Spectral",
... )
>>> fig.savefig("images/mlpr_1713.png", dpi=300)
图 17-13. UMAP 结果。
要调整 UMAP 结果，首先关注n_neighbors和min_dist超参数。以下是更改这些值的示例（见图 17-14 和 17-15）：
>>> X_std = StandardScaler().fit_transform(X)
>>> fig, axes = plt.subplots(2, 2, figsize=(6, 4))
>>> axes = axes.reshape(4)
>>> for i, n in enumerate([2, 5, 10, 50]):
... ax = axes[i]
... u = umap.UMAP(
... random_state=42, n_neighbors=n
... )
... X_umap = u.fit_transform(X_std)
...
... pd.DataFrame(X_umap).plot(
... kind="scatter",
... x=0,
... y=1,
... ax=ax,
... c=y,
... cmap="Spectral",
... alpha=0.5,
... )
... ax.set_title(f"nn={n}")
>>> plt.tight_layout()
>>> fig.savefig("images/mlpr_1714.png", dpi=300)
图 17-14. 调整 UMAP 结果n_neighbors。
>>> fig, axes = plt.subplots(2, 2, figsize=(6, 4))
>>> axes = axes.reshape(4)
>>> for i, n in enumerate([0, 0.33, 0.66, 0.99]):
... ax = axes[i]
... u = umap.UMAP(random_state=42, min_dist=n)
... X_umap = u.fit_transform(X_std)
... pd.DataFrame(X_umap).plot(
... kind="scatter",
... x=0,
... y=1,
... ax=ax,
... c=y,
... cmap="Spectral",
... alpha=0.5,
... )
... ax.set_title(f"min_dist={n}")
>>> plt.tight_layout()
>>> fig.savefig("images/mlpr_1715.png", dpi=300)
图 17-15. 调整 UMAP 结果min_dist。
有时在 UMAP 之前执行 PCA 以减少维数并加快计算速度。
t-SNE
t-分布随机邻域嵌入（t-SNE）技术是一种可视化和降维技术。它使用输入的分布和低维嵌入，并最小化它们之间的联合概率。由于计算量大，可能无法在大数据集上使用这种技术。
t-SNE 的一个特征是对超参数非常敏感。此外，虽然它能够很好地保留局部聚类，但全局信息并未保留。最后，这不是一个确定性算法，可能不会收敛。
在使用此技术之前标准化数据是一个好主意：
>>> from sklearn.manifold import TSNE
>>> X_std = StandardScaler().fit_transform(X)
>>> ts = TSNE()
>>> X_tsne = ts.fit_transform(X_std)
实例参数：
n_components=2
嵌入的维度数。
perplexity=30.0
建议的取值范围为 5 到 50。较小的数值倾向于形成更紧密的聚类。
early_exaggeration=12.0
控制簇紧密度和它们之间的间距。较大的值意味着较大的间距。
learning_rate=200.0
通常在 10 到 1000 之间。如果数据看起来像一个球，则降低它。如果数据看起来压缩，则增加它。
n_iter=1000
迭代次数。
n_iter_without_progress=300
如果在这些迭代次数之后没有进展，则中止。
min_grad_norm=1e-07
如果梯度范数低于此值，则优化停止。
metric='euclidean'
来自scipy.spatial.distance.pdist、pairwise.PAIRWISE_DISTANCE_METRIC或函数的距离度量。
init='random'
嵌入初始化。
verbose=0
冗长性。
random_state=None
随机种子。
method='barnes_hut'
梯度计算算法。
angle=0.5
用于梯度计算。小于 0.2 会增加运行时间。大于 0.8 会增加错误。
属性：
embedding_
嵌入向量。
kl_divergence_
Kullback-Leibler 散度。
n_iter_
迭代次数。
这里展示了使用 matplotlib 进行 t-SNE 结果的可视化（见图 17-16）：
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> colors = ["rg"[j] for j in y]
>>> scat = ax.scatter(
... X_tsne[:, 0],
... X_tsne[:, 1],
... c=colors,
... alpha=0.5,
... )
>>> ax.set_xlabel("Embedding 1")
>>> ax.set_ylabel("Embedding 2")
>>> fig.savefig("images/mlpr_1716.png", dpi=300)
图 17-16。使用 matplotlib 进行 t-SNE 结果。
改变perplexity的值可能会对绘图产生重大影响（见图 17-17）。
>>> fig, axes = plt.subplots(2, 2, figsize=(6, 4))
>>> axes = axes.reshape(4)
>>> for i, n in enumerate((2, 30, 50, 100)):
... ax = axes[i]
... t = TSNE(random_state=42, perplexity=n)
... X_tsne = t.fit_transform(X)
... pd.DataFrame(X_tsne).plot(
... kind="scatter",
... x=0,
... y=1,
... ax=ax,
... c=y,
... cmap="Spectral",
... alpha=0.5,
... )
... ax.set_title(f"perplexity={n}")
... plt.tight_layout()
... fig.savefig("images/mlpr_1717.png", dpi=300)
图 17-17。改变t-SNE的perplexity。
PHATE
通过 Affinity-based Trajectory Embedding（PHATE）进行热扩散的潜力是高维数据可视化的工具。它倾向于同时保留全局结构（如 PCA）和局部结构（如 t-SNE）。
PHATE 首先编码局部信息（接近的点应该保持接近）。它使用“扩散”来发现全局数据，然后减少维度：
>>> import phate
>>> p = phate.PHATE(random_state=42)
>>> X_phate = p.fit_transform(X)
>>> X_phate.shape
实例参数：
n_components=2
维度数。
knn=5
核的邻居数。如果嵌入是断开的或数据集大于 10 万个样本，则增加。
decay=40
核的衰减率。降低此值会增加图的连通性。
n_landmark=2000
用于标记的地标点。
t='auto'
扩散力度。对数据进行平滑处理。如果嵌入缺乏结构，请增加；如果结构紧密而紧凑，请减少。
gamma=1
对数潜力（在 -1 到 1 之间）。如果嵌入集中在单个点周围，请尝试将其设置为 0。
n_pca=100
用于邻域计算的主成分数。
knn_dist='euclidean'
KNN 距离度量。
mds_dist='euclidean'
多维缩放（MDS）度量。
mds='metric'
MDS 算法用于降维。
n_jobs=1
要使用的 CPU 数量。
random_state=None
随机种子。
verbose=1
冗余性。
属性（请注意这些后面没有 _）：
X
输入数据
embedding
嵌入空间
diff_op
扩散算子
graph
基于输入构建的 KNN 图
这是使用 PHATE 的一个示例（见 Figure 17-18）：
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> phate.plot.scatter2d(p, c=y, ax=ax, alpha=0.5)
>>> fig.savefig("images/mlpr_1718.png", dpi=300)
图 17-18. PHATE 结果。
如上所述的实例参数中，有一些参数可以调整以改变模型的行为。以下是调整 knn 参数的示例（见 Figure 17-19）。请注意，如果使用 .set_params 方法，它将加快计算速度，因为它使用预计算的图和扩散算子：
>>> fig, axes = plt.subplots(2, 2, figsize=(6, 4))
>>> axes = axes.reshape(4)
>>> p = phate.PHATE(random_state=42, n_jobs=-1)
>>> for i, n in enumerate((2, 5, 20, 100)):
... ax = axes[i]
... p.set_params(knn=n)
... X_phate = p.fit_transform(X)
... pd.DataFrame(X_phate).plot(
... kind="scatter",
... x=0,
... y=1,
... ax=ax,
... c=y,
... cmap="Spectral",
... alpha=0.5,
... )
... ax.set_title(f"knn={n}")
... plt.tight_layout()
... fig.savefig("images/mlpr_1719.png", dpi=300)
图 17-19. 更改 PHATE 的 knn 参数。
第十八章：聚类
聚类是一种无监督机器学习技术，用于将群体分成几个组。它是无监督的，因为我们没有给模型任何标签；它只是检查特征并确定哪些样本相似并属于一个簇。在本章中，我们将研究 K 均值和层次聚类方法。我们还将再次使用各种技术探索泰坦尼克号数据集。
K 均值
K 均值算法需要用户选择簇的数量或“k”。然后它随机选择 k 个质心，并根据从质心到每个样本的距离度量将每个样本分配给一个簇。分配完成后，它根据分配给一个标签的每个样本的中心重新计算质心。然后它根据新的质心再次将样本分配到簇中。经过几次迭代后，它应该会收敛。
因为聚类使用距离度量来确定哪些样本相似，所以行为可能会根据数据的规模而变化。您可以标准化数据并使所有特征处于相同的比例。有些人建议，如果规模提示某些特征更重要，SME 可能会建议不要标准化。我们将在这个例子中对数据进行标准化。
在这个例子中，我们将对泰坦尼克号乘客进行聚类。我们将从两个簇开始，看看聚类是否能够分开生存（我们不会将生存数据泄漏到聚类中，只使用X，而不是y）。
无监督算法具有.fit方法和.predict方法。我们只将X传递给.fit：
>>> from sklearn.cluster import KMeans
>>> X_std = preprocessing.StandardScaler().fit_transform(
... X
... )
>>> km = KMeans(2, random_state=42)
>>> km.fit(X_std)
KMeans(algorithm='auto', copy_x=True,
init='k-means', max_iter=300,
n_clusters=2, n_init=10, n_jobs=1,
precompute_distances='auto',
random_state=42, tol=0.0001, verbose=0)
在模型训练后，我们可以调用.predict方法将新样本分配给一个簇：
>>> X_km = km.predict(X)
>>> X_km
array([1, 1, 1, ..., 1, 1, 1], dtype=int32)
实例参数：
n_clusters=8
创建的簇的数量。
init='kmeans++'
初始化方法。
n_init=10
使用不同质心运行算法的次数。最佳得分将获胜。
max_iter=300
运行的迭代次数。
tol=0.0001
收敛的公差。
precompute_distances='auto'
预计算距离（需要更多内存但更快）。如果n_samples * n_clusters小于或等于 1200 万，auto将预先计算。
verbose=0
冗余性。
random_state=None
随机种子。
copy_x=True
在计算之前复制数据。
n_jobs=1
要使用的 CPU 数量。
algorithm='auto'
K 均值算法。'full'适用于稀疏数据，但'elkan'更高效。'auto'在密集数据中使用'elkan'。
属性：
cluster_centers_
中心的坐标
labels_
样本的标签
inertia_
到聚类质心的平方距离之和
n_iter_
迭代次数
如果事先不知道需要多少个簇，可以以一系列大小运行算法并评估各种指标。这可能有些棘手。
您可以使用 .inertia_ 计算自己的肘部图。寻找曲线弯曲的位置，因为那可能是选择聚类数量的一个良好选择。在这种情况下，曲线很平滑，但在八个之后似乎没有太多改善（见 Figure 18-1）。
对于没有肘部的图，我们有几个选择。我们可以使用其他指标，其中一些如下所示。我们还可以检查聚类的可视化，看看聚类是否可见。我们可以向数据添加特征，看看是否有帮助进行聚类。
这里是肘部图的代码：
>>> inertias = []
>>> sizes = range(2, 12)
>>> for k in sizes:
... k2 = KMeans(random_state=42, n_clusters=k)
... k2.fit(X)
... inertias.append(k2.inertia_)
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> pd.Series(inertias, index=sizes).plot(ax=ax)
>>> ax.set_xlabel("K")
>>> ax.set_ylabel("Inertia")
>>> fig.savefig("images/mlpr_1801.png", dpi=300)
图 18-1. 看起来相当平滑的肘部图。
当地面真实标签未知时，Scikit-learn 具有其他聚类指标。我们也可以计算并绘制这些指标。轮廓系数 是介于 -1 和 1 之间的值。得分越高越好。1 表示紧密的聚类，0 表示重叠的聚类。根据这个度量，两个聚类给我们带来了最佳分数。
Calinski-Harabasz 指数 是介于类间离散度和类内离散度之间的比率。分数越高越好。对于这个指标，两个聚类给出了最佳分数。
Davis-Bouldin 指数 是每个聚类与最接近的聚类之间相似性的平均值。分数从 0 开始。0 表示更好的聚类。
在这里，我们将绘制惯性、轮廓系数、Calinski-Harabasz 指数和 Davies-Bouldin 指数在一系列聚类大小上的情况，以查看数据是否有明确的聚类大小（见 Figure 18-2）。大多数这些指标都同意选择两个聚类：
>>> from sklearn import metrics
>>> inertias = []
>>> sils = []
>>> chs = []
>>> dbs = []
>>> sizes = range(2, 12)
>>> for k in sizes:
... k2 = KMeans(random_state=42, n_clusters=k)
... k2.fit(X_std)
... inertias.append(k2.inertia_)
... sils.append(
... metrics.silhouette_score(X, k2.labels_)
... )
... chs.append(
... metrics.calinski_harabasz_score(
... X, k2.labels_
... )
... )
... dbs.append(
... metrics.davies_bouldin_score(
... X, k2.labels_
... )
... )
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> (
... pd.DataFrame(
... {
... "inertia": inertias,
... "silhouette": sils,
... "calinski": chs,
... "davis": dbs,
... "k": sizes,
... }
... )
... .set_index("k")
... .plot(ax=ax, subplots=True, layout=(2, 2))
... )
>>> fig.savefig("images/mlpr_1802.png", dpi=300)
图 18-2. 聚类指标。这些指标大多数同意选择两个聚类。
另一种确定聚类的技术是可视化每个聚类的轮廓分数。Yellowbrick 有一个此类的可视化工具（见 Figure 18-3）。
在这个图中，垂直的虚线是平均分数。一种解释方法是确保每个聚类都突出于平均水平之上，并且聚类分数看起来还不错。确保使用相同的 x 轴限制 (ax.set_xlim)。我会从这些图中选择两个聚类：
>>> from yellowbrick.cluster.silhouette import (
... SilhouetteVisualizer,
... )
>>> fig, axes = plt.subplots(2, 2, figsize=(12, 8))
>>> axes = axes.reshape(4)
>>> for i, k in enumerate(range(2, 6)):
... ax = axes[i]
... sil = SilhouetteVisualizer(
... KMeans(n_clusters=k, random_state=42),
... ax=ax,
... )
... sil.fit(X_std)
... sil.finalize()
... ax.set_xlim(-0.2, 0.8)
>>> plt.tight_layout()
>>> fig.savefig("images/mlpr_1803.png", dpi=300)
图 18-3. Yellowbrick 轮廓可视化器
聚合（层次）聚类
聚合聚类是另一种方法。你从每个样本单独的聚类开始。然后将“最近的”聚类组合起来。重复此过程，同时跟踪最近的大小。
当您完成这些操作时，您将得到一棵dendrogram，或者一棵跟踪聚类创建时间和距离度量的树。您可以使用 scipy 库可视化这棵树。
我们可以使用 scipy 创建一棵树状图（见 Figure 18-4）。如您所见，如果样本很多，则叶节点很难阅读：
>>> from scipy.cluster import hierarchy
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> dend = hierarchy.dendrogram(
... hierarchy.linkage(X_std, method="ward")
... )
>>> fig.savefig("images/mlpr_1804.png", dpi=300)
图 18-4. Scipy 层次聚类树状图
一旦您有了树状图，您就拥有了所有的聚类（从样本数为一到样本数为止）。高度表示加入时相似聚类的相似程度。为了找出数据中有多少个聚类，您需要在最高的线交叉处“切割”一条水平线。
在这种情况下，当您执行该切割时，看起来有三个聚类。
前一个图表太嘈杂了，其中包含了所有的样本。您还可以使用truncate_mode参数将叶子节点合并为单个节点（参见图 18-5）：
>>> from scipy.cluster import hierarchy
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> dend = hierarchy.dendrogram(
... hierarchy.linkage(X_std, method="ward"),
... truncate_mode="lastp",
... p=20,
... show_contracted=True,
... )
>>> fig.savefig("images/mlpr_1805.png", dpi=300)
图 18-5. 截断的层次聚类树状图。如果我们在最大的垂直线上切割，我们将得到三个聚类。
一旦我们知道需要多少个聚类，我们可以使用 scikit-learn 创建一个模型：
>>> from sklearn.cluster import (
... AgglomerativeClustering,
... )
>>> ag = AgglomerativeClustering(
... n_clusters=4,
... affinity="euclidean",
... linkage="ward",
... )
>>> ag.fit(X)
注意
fastcluster 包 提供了一个优化的凝聚聚类包，如果 scikit-learn 的实现速度太慢的话。
理解聚类
在 Titanic 数据集上使用 K-means，我们将得到两个聚类。我们可以使用 pandas 中的分组功能来检查聚类之间的差异。下面的代码检查每个特征的均值和方差。看起来 pclass 的均值变化相当大。
我将生存数据重新放回，看看聚类是否与此相关：
>>> km = KMeans(n_clusters=2)
>>> km.fit(X_std)
>>> labels = km.predict(X_std)
>>> (
... X.assign(cluster=labels, survived=y)
... .groupby("cluster")
... .agg(["mean", "var"])
... .T
... )
cluster 0 1
pclass mean 0.526538 -1.423831
var 0.266089 0.136175
age mean -0.280471 0.921668
var 0.653027 1.145303
sibsp mean -0.010464 -0.107849
var 1.163848 0.303881
parch mean 0.387540 0.378453
var 0.829570 0.540587
fare mean -0.349335 0.886400
var 0.056321 2.225399
sex_male mean 0.678986 0.552486
var 0.218194 0.247930
embarked_Q mean 0.123548 0.016575
var 0.108398 0.016345
embarked_S mean 0.741288 0.585635
var 0.191983 0.243339
survived mean 0.596685 0.299894
var 0.241319 0.210180
注意
在 Jupyter 中，您可以将以下代码添加到 DataFrame 中，并突出显示每行的高低值。这对于直观地查看哪些值在上述聚类摘要中显著是有用的：
.style.background_gradient(cmap='RdBu', axis=1)
在图 18-6 中，我们绘制了每个聚类的均值条形图：
>>> fig, ax = plt.subplots(figsize=(6, 4))
... (
... X.assign(cluster=labels, survived=y)
... .groupby("cluster")
... .mean()
... .T.plot.bar(ax=ax)
... )
>>> fig.savefig(
... "images/mlpr_1806.png",
... dpi=300,
... bbox_inches="tight",
... )
图 18-6. 每个聚类的均值
我还喜欢绘制 PCA 组件，但是按聚类标签着色（见图 18-7）。在这里，我们使用 Seaborn 来执行此操作。将hue的值更改为深入研究聚类中显著的特征也很有趣。
>>> fig, ax = plt.subplots(figsize=(6, 4))
>>> sns.scatterplot(
... "PC1",
... "PC2",
... data=X.assign(
... PC1=X_pca[:, 0],
... PC2=X_pca[:, 1],
... cluster=labels,
... ),
... hue="cluster",
... alpha=0.5,
... ax=ax,
... )
>>> fig.savefig(
... "images/mlpr_1807.png",
... dpi=300,
... bbox_inches="tight",
... )
图 18-7. 聚类的 PCA 图
如果我们想检查单个特征，可以使用 pandas 的.describe方法：
>>> (
... X.assign(cluster=label)
... .groupby("cluster")
... .age.describe()
... .T
... )
cluster 0 1
count 362.000000 947.000000
mean 0.921668 -0.280471
std 1.070188 0.808101
min -2.160126 -2.218578
25% 0.184415 -0.672870
50% 0.867467 -0.283195
75% 1.665179 0.106480
max 4.003228 3.535618
我们还可以创建一个替代模型来解释这些聚类。在这里，我们使用决策树来解释它们。这还显示了 pclass（其均值差异很大）非常重要：
>>> dt = tree.DecisionTreeClassifier()
>>> dt.fit(X, labels)
>>> for col, val in sorted(
... zip(X.columns, dt.feature_importances_),
... key=lambda col_val: col_val[1],
... reverse=True,
... ):
... print(f"{col:10}{val:10.3f}")
pclass 0.902
age 0.074
sex_male 0.016
embarked_S 0.003
fare 0.003
parch 0.003
sibsp 0.000
embarked_Q 0.000
我们可以通过图 18-8 来可视化决策。它显示 pclass 是第一个特征，用于做出决策：
>>> dot_data = StringIO()
>>> tree.export_graphviz(
... dt,
... out_file=dot_data,
... feature_names=X.columns,
... class_names=["0", "1"],
... max_depth=2,
... filled=True,
... )
>>> g = pydotplus.graph_from_dot_data(
... dot_data.getvalue()
... )
>>> g.write_png("images/mlpr_1808.png")
图 18-8. 解释聚类的决策树
第十九章：管道
Scikit-learn 使用管道的概念。使用 Pipeline 类，您可以将转换器和模型链在一起，并像对待整个 scikit-learn 模型一样对待整个过程。甚至可以插入自定义逻辑。
分类管道
这是在管道内使用 tweak_titanic 函数的示例：
>>> from sklearn.base import (
... BaseEstimator,
... TransformerMixin,
... )
>>> from sklearn.pipeline import Pipeline
>>> def tweak_titanic(df):
... df = df.drop(
... columns=[
... "name",
... "ticket",
... "home.dest",
... "boat",
... "body",
... "cabin",
... ]
... ).pipe(pd.get_dummies, drop_first=True)
... return df
>>> class TitanicTransformer(
... BaseEstimator, TransformerMixin
... ):
... def transform(self, X):
... # assumes X is output
... # from reading Excel file
... X = tweak_titanic(X)
... X = X.drop(column="survived")
... return X
...
... def fit(self, X, y):
... return self
>>> pipe = Pipeline(
... [
... ("titan", TitanicTransformer()),
... ("impute", impute.IterativeImputer()),
... (
... "std",
... preprocessing.StandardScaler(),
... ),
... ("rf", RandomForestClassifier()),
... ]
... )
有了管道，我们可以对其调用 .fit 和 .score：
>>> from sklearn.model_selection import (
... train_test_split,
... )
>>> X_train2, X_test2, y_train2, y_test2 = train_test_split(
... orig_df,
... orig_df.survived,
... test_size=0.3,
... random_state=42,
... )
>>> pipe.fit(X_train2, y_train2)
>>> pipe.score(X_test2, y_test2)
0.7913486005089059
管道可以在网格搜索中使用。我们的 param_grid 需要将参数前缀设为管道阶段的名称，后跟两个下划线。在下面的示例中，我们为随机森林阶段添加了一些参数：
>>> params = {
... "rf__max_features": [0.4, "auto"],
... "rf__n_estimators": [15, 200],
... }
>>> grid = model_selection.GridSearchCV(
... pipe, cv=3, param_grid=params
... )
>>> grid.fit(orig_df, orig_df.survived)
现在我们可以提取最佳参数并训练最终模型。（在这种情况下，随机森林在网格搜索后没有改善。）
>>> grid.best_params_
{'rf__max_features': 0.4, 'rf__n_estimators': 15}
>>> pipe.set_params(**grid.best_params_)
>>> pipe.fit(X_train2, y_train2)
>>> pipe.score(X_test2, y_test2)
0.7913486005089059
我们可以在使用 scikit-learn 模型的管道中使用：
>>> metrics.roc_auc_score(
... y_test2, pipe.predict(X_test2)
... )
0.7813688715131023
回归管道
这是在波士顿数据集上执行线性回归的管道示例：
>>> from sklearn.pipeline import Pipeline
>>> reg_pipe = Pipeline(
... [
... (
... "std",
... preprocessing.StandardScaler(),
... ),
... ("lr", LinearRegression()),
... ]
... )
>>> reg_pipe.fit(bos_X_train, bos_y_train)
>>> reg_pipe.score(bos_X_test, bos_y_test)
0.7112260057484934
如果我们想要从管道中提取部分来检查它们的属性，我们可以使用 .named_steps 属性进行操作：
>>> reg_pipe.named_steps["lr"].intercept_
23.01581920903956
>>> reg_pipe.named_steps["lr"].coef_
array([-1.10834602, 0.80843998, 0.34313466,
0.81386426, -1.79804295, 2.913858 ,
-0.29893918, -2.94251148, 2.09419303,
-1.44706731, -2.05232232, 1.02375187,
-3.88579002])_
我们也可以在度量计算中使用管道：
>>> from sklearn import metrics
>>> metrics.mean_squared_error(
... bos_y_test, reg_pipe.predict(bos_X_test)
... )
21.517444231177205
PCA 管道
Scikit-learn 管道也可以用于 PCA。
这里我们对泰坦尼克号数据集进行标准化并对其执行 PCA：
>>> pca_pipe = Pipeline(
... [
... (
... "std",
... preprocessing.StandardScaler(),
... ),
... ("pca", PCA()),
... ]
... )
>>> X_pca = pca_pipe.fit_transform(X)
使用 .named_steps 属性，我们可以从管道的 PCA 部分提取属性：
>>> pca_pipe.named_steps[
... "pca"
... ].explained_variance_ratio_
array([0.23917891, 0.21623078, 0.19265028,
0.10460882, 0.08170342, 0.07229959,
0.05133752, 0.04199068])
>>> pca_pipe.named_steps["pca"].components_[0]
array([-0.63368693, 0.39682566, 0.00614498,
0.11488415, 0.58075352, -0.19046812,
-0.21190808, -0.09631388])
