Python-数值计算-二-

Python 数值计算（二）

原文：annas-archive.org/md5/9e81efca12eeaa9a42c1e05702b8c0c0

译者：飞龙

协议：CC BY-NC-SA 4.0

第四章：无监督学习

无监督学习的目标是发现数据中的隐藏模式或结构，其中不存在用于执行分类或回归方法的目标变量。无监督学习方法通常更具挑战性，因为结果是主观的，并且没有像预测类别或连续变量这样简单的分析目标。这些方法通常作为探索性数据分析的一部分进行。此外，由于没有公认的结果验证机制，因此很难评估无监督学习方法获得的结果。

尽管如此，无监督学习方法在各个领域的应用越来越重要，成为当前的热门话题，许多研究人员正在积极研究这些方法，探索这一新领域。以下是一些良好的应用示例：

• 基因组学：将无监督学习应用于理解基因组范围的生物学见解，从 DNA 入手更好地理解疾病和人类。这些任务本质上更具探索性。

• 搜索引擎：搜索引擎可能会根据其他相似用户的点击历史来选择显示哪些搜索结果给特定个体。

• 知识提取：从原始文本中提取概念的分类法，生成知识图谱，以创建自然语言处理（NLP）领域的语义结构。

• 客户细分：在银行业中，像聚类这样的无监督学习被应用于将相似的客户分组，并根据这些细分，营销部门设计他们的接触策略。例如，年龄较大、低风险的客户将通过定期存款产品进行目标营销，而高风险、年轻的客户则会通过信用卡或共同基金等进行营销，等等。

• 社交网络分析：识别网络中彼此更紧密连接、具有共同特征的凝聚群体。

在本章中，我们将介绍以下技术，使用公开可获得的数据进行无监督学习：

• K-means 聚类

• 主成分分析

• 奇异值分解

• 深度自编码器

K-means 聚类

聚类是将观察值分组的任务，确保同一聚类的成员彼此之间更加相似，而不同聚类的成员之间差异较大。

聚类通常用于探索数据集，以识别其中的潜在模式或创建一组特征。在社交网络中，可以通过聚类来识别社区并建议人们之间缺失的连接。以下是一些示例：

• 在反洗钱措施中，可以通过异常检测识别可疑的活动和个人。

• 在生物学中，聚类用于寻找具有相似表达模式的基因群体。

• 在市场营销分析中，聚类用于寻找相似客户的细分群体，从而为不同的客户群体制定不同的营销策略。

k-means 聚类算法是一个迭代过程，它将聚类中心或质心移动到其组成点的均值位置，并反复将实例分配到最接近的聚类，直到聚类中心数量没有显著变化或达到最大迭代次数。

k-means 的代价函数由属于某个聚类的观测值与其相应质心之间的欧几里得距离（平方范数）决定。理解该公式的一种直观方法是，如果只有一个聚类（k=1），那么所有观测值与该单一均值的距离将被比较。而如果聚类数增加到2（k=2），则计算两个均值，一些观测值被分配到聚类1，另一些观测值根据距离分配到聚类2。随后，通过应用相同的距离度量，但分别对其聚类中心计算距离，从而得到代价函数：

从第一性原理出发的 K-means 工作方法论

k-means 的工作方法论通过以下示例进行了说明，在该示例中考虑了 12 个实例及其X和Y值。任务是从数据中确定最优的聚类。

实例	X	Y
1	7	8
2	2	4
3	6	4
4	3	2
5	6	5
6	5	7
7	3	3
8	1	4
9	5	4
10	7	7
11	7	6
12	2	1

在将数据点绘制到二维图表上后，我们可以看到大致上有两个聚类可能性，其中左下角是第一个聚类，右上角是另一个聚类。但在许多实际案例中，变量（或维度）会非常多，以至于我们无法简单地将其可视化。因此，我们需要一种数学和算法的方法来解决这些问题。

迭代 1：假设从所有12个实例中选择两个中心。这里，我们选择了实例1（X = 7, Y = 8）和实例8（X = 1, Y = 4），因为它们似乎位于两个极端。对于每个实例，我们将计算其相对于两个质心的欧几里得距离，并将其分配到最近的聚类中心。

实例	X	Y	质心 1 距离	质心 2 距离	分配的聚类
1	7	8	7.21	0.00	C2
2	2	4	1.00	6.40	C1
3	6	4	5.00	4.12	C2
4	3	2	2.83	7.21	C1
5	6	5	5.10	3.16	C2
6	5	7	5.00	2.24	C2
7	3	3	2.24	6.40	C1
8	1	4	0.00	7.21	C1
9	5	4	4.00	4.47	C1
10	7	7	6.71	1.00	C2
11	7	6	6.32	2.00	C2
12	2	1	3.16	8.60	C1
聚类中心 1	1	4
聚类中心 2	7	8

两点之间的欧几里得距离 A (X1, Y1) 和 B (X2, Y2) 如下所示：

聚类中心的距离计算是通过计算欧几里得距离来进行的。以下展示了一个样本计算。例如，六号实例与两个聚类中心（聚类中心 1 和聚类中心 2）之间的距离。

以下图表描述了实例到两个聚类中心的分配情况，格式与前面的表格相同：

如果仔细观察前面的图表，我们会发现，除了实例*9 (X =5, Y = 4)*外，所有实例似乎都已被正确分配。然而，在后续阶段，它应该被正确分配。让我们看看下面的步骤中分配是如何演变的。

迭代 2：在此迭代中，新的聚类中心由分配给该聚类或聚类中心的实例计算得出。新的聚类中心是基于所分配点的简单平均值来计算的。

实例	X	Y	分配的聚类
1	7	8	C2
2	2	4	C1
3	6	4	C2
4	3	2	C1
5	6	5	C2
6	5	7	C2
7	3	3	C1
8	1	4	C1
9	5	4	C1
10	7	7	C2
11	7	6	C2
12	2	1	C1
聚类中心 1	2.67	3
聚类中心 2	6.33	6.17

聚类中心 1 和 2 的样本计算如下所示。类似的方法将应用于所有后续的迭代：

更新聚类中心后，我们需要将实例重新分配给最近的聚类中心，这将在迭代 3 中执行。

迭代 3：在此迭代中，新的分配是基于实例与新聚类中心之间的欧几里得距离计算的。如果发生任何变化，将反复计算新的聚类中心，直到分配没有变化或达到迭代次数为止。下表描述了新聚类中心与所有实例之间的距离度量：

实例	X	Y	质心 1 距离	质心 2 距离	先前分配的聚类	新分配的聚类	是否变化？
1	7	8	6.61	1.95	C2	C2	No
2	2	4	1.20	4.84	C1	C1	No
3	6	4	3.48	2.19	C2	C2	No
4	3	2	1.05	5.34	C1	C1	No
5	6	5	3.88	1.22	C2	C2	No
6	5	7	4.63	1.57	C2	C2	No
7	3	3	0.33	4.60	C1	C1	No
8	1	4	1.95	5.75	C1	C1	No
9	5	4	2.54	2.55	C1	C1	No
10	7	7	5.89	1.07	C2	C2	No
11	7	6	5.27	0.69	C2	C2	No
12	2	1	2.11	6.74	C1	C1	No
质心 1	2.67	3
质心 2	6.33	6.17

似乎没有注册到任何变化。因此，我们可以说，解决方案已经收敛。需要注意的一点是，除了实例*9 (X = 5, Y = 4)*外，所有实例都已非常清晰地被分类。根据直觉，它应该分配给质心 2，但经过仔细计算后，该实例实际上更接近聚类 1 而非聚类 2。然而，距离的差异非常小（质心 1 的距离为 2.54，质心 2 的距离为 2.55）。

最优聚类数量与聚类评估

尽管选择聚类的数量更像是一门艺术而非科学，但在选择最优聚类数时，若增加聚类数后，解释能力的提升非常小，则可以选择该数量。在实际应用中，通常业务方应该能够提供他们大致需要的聚类数量。

肘部法则

肘部法则用于确定 k-means 聚类中的最优聚类数量。肘部法则绘制不同k值下由代价函数产生的值。如你所知，当k增加时，平均失真度会降低，每个聚类的构成实例会减少，并且实例会更接近各自的质心。然而，随着k的增加，平均失真度的改善会逐渐减缓。失真度改善下降最明显的k值称为肘部，应该在此停止继续划分数据为更多的聚类。

使用轮廓系数评估聚类：轮廓系数是衡量聚类紧凑性和分离度的指标。较高的值代表更好的聚类质量。轮廓系数对紧凑且分离良好的聚类较高，而对重叠聚类则较低。轮廓系数值变化范围从 -1 到 +1，值越高，聚类质量越好。

轮廓系数是按实例计算的。对于一组实例，它是各个样本得分的平均值。

a 是簇内实例之间的平均距离，b 是该实例与下一个最近簇中实例之间的平均距离。

使用鸢尾花数据集进行 K-means 聚类示例

经典的鸢尾花数据集来自 UCI 机器学习库，用于演示 K-means 聚类。数据下载链接在此：archive.ics.uci.edu/ml/datasets/Iris。鸢尾花数据集包含三种花卉：山鸢尾、变色鸢尾和维吉尼亚鸢尾，以及它们的萼片长度、萼片宽度、花瓣长度和花瓣宽度的相应测量值。我们的任务是根据这些测量值将花卉分组。代码如下：

>>> import os 
""" First change the following directory link to where all input files do exist """ 
>>> os.chdir("D:\\Book writing\\Codes\\Chapter 8") 

K-means algorithm from scikit-learn has been utilized in the following example 

# K-means clustering 
>>> import numpy as np 
>>> import pandas as pd 
>>> import matplotlib.pyplot as plt 
>>> from scipy.spatial.distance import cdist, pdist 

>>> from sklearn.cluster import KMeans 
>>> from sklearn.metrics import silhouette_score

>>> iris = pd.read_csv("iris.csv") 
>>> print (iris.head()) 

以下代码用于将 class 变量作为依赖变量来为图中的颜色创建并应用无监督学习算法，且在给定的 x 变量上进行操作，没有目标变量：

>>> x_iris = iris.drop(['class'],axis=1) 
>>> y_iris = iris["class"] 

作为示例，使用了三个簇，但在现实生活中，我们无法事先知道数据将属于多少个簇，因此需要通过反复试验来测试结果。这里选择的最大迭代次数为 300，当然，这个值也可以调整并相应地检查结果：

>>> k_means_fit = KMeans(n_clusters=3,max_iter=300) 
>>> k_means_fit.fit(x_iris) 

>>> print ("\nK-Means Clustering - Confusion Matrix\n\n",pd.crosstab(y_iris, k_means_fit.labels_,rownames = ["Actuall"],colnames = ["Predicted"]) )      
>>> print ("\nSilhouette-score: %0.3f" % silhouette_score(x_iris, k_means_fit.labels_, metric='euclidean')) 

从之前的混淆矩阵中，我们可以看到所有的山鸢尾花都被正确地聚类，而 50 个变色鸢尾花中的 2 个和 50 个维吉尼亚鸢尾花中的 14 个被错误地分类。

再次强调，在现实生活中的例子中，我们事先并不知道类别名称，因此无法衡量准确性等指标。

以下代码用于执行敏感性分析，检查实际提供更好细分解释的簇的数量：

>>> for k in range(2,10): 
...     k_means_fitk = KMeans(n_clusters=k,max_iter=300) 
...     k_means_fitk.fit(x_iris) 
...     print ("For K value",k,",Silhouette-score: %0.3f" % silhouette_score(x_iris, k_means_fitk.labels_, metric='euclidean')) 

上述结果中的轮廓系数值显示，K 值 2 和 K 值 3 的得分优于其他所有值。作为经验法则，我们需要选择轮廓系数最高的下一个 K 值。在这里，我们可以说 K 值 3 更好。此外，在得出最佳 K 值 之前，我们还需要查看每个簇内的平均变化值和肘部图。

# Avg. within-cluster sum of squares 
>>> K = range(1,10) 

>>> KM = [KMeans(n_clusters=k).fit(x_iris) for k in K] 
>>> centroids = [k.cluster_centers_ for k in KM] 

>>> D_k = [cdist(x_iris, centrds, 'euclidean') for centrds in centroids] 

>>> cIdx = [np.argmin(D,axis=1) for D in D_k] 
>>> dist = [np.min(D,axis=1) for D in D_k] 
>>> avgWithinSS = [sum(d)/x_iris.shape[0] for d in dist] 

# Total with-in sum of square 
>>> wcss = [sum(d**2) for d in dist] 
>>> tss = sum(pdist(x_iris)**2)/x_iris.shape[0] 
>>> bss = tss-wcss 

# elbow curve - Avg. within-cluster sum of squares 
>>> fig = plt.figure() 
>>> ax = fig.add_subplot(111) 
>>> ax.plot(K, avgWithinSS, 'b*-') 
>>> plt.grid(True) 
>>> plt.xlabel('Number of clusters') 
>>> plt.ylabel('Average within-cluster sum of squares') 

从肘部图中看，值为三时，斜率发生了剧烈变化。在这里，我们可以选择最优的 k 值为三。

# elbow curve - percentage of variance explained 
>>> fig = plt.figure() 
>>> ax = fig.add_subplot(111) 
>>> ax.plot(K, bss/tss*100, 'b*-') 
>>> plt.grid(True) 
>>> plt.xlabel('Number of clusters') 
>>> plt.ylabel('Percentage of variance explained')
>>> plt.show()

最后但同样重要的是，解释的总方差百分比值应该大于 80%，以决定最优的聚类数量。即使在这里，k 值为三似乎也能提供一个合理的总方差解释值。因此，我们可以从前述的所有指标（轮廓系数、聚类内平均方差和总方差解释）得出结论，三类聚类是理想的。

使用鸢尾花数据的 k 均值聚类的 R 代码如下：

setwd("D:\\Book writing\\Codes\\Chapter 8")   

iris_data = read.csv("iris.csv")   
x_iris =   iris_data[,!names(iris_data) %in% c("class")]   
y_iris = iris_data$class   

km_fit = kmeans(x_iris,centers   = 3,iter.max = 300 )   

print(paste("K-Means   Clustering- Confusion matrix"))   
table(y_iris,km_fit$cluster)   

mat_avgss = matrix(nrow = 10,   ncol = 2)   

# Average within the cluster   sum of square   
print(paste("Avg. Within   sum of squares"))   
for (i in (1:10)){   
  km_fit =   kmeans(x_iris,centers = i,iter.max = 300 )   
  mean_km =   mean(km_fit$withinss)   
  print(paste("K-Value",i,",Avg.within   sum of squares",round(mean_km, 2)))   
  mat_avgss[i,1] = i   
  mat_avgss[i,2] = mean_km   
}   
 plot(mat_avgss[,1],mat_avgss[,2],type   = 'o',xlab = "K_Value",ylab = "Avg. within sum of square")   
title("Avg. within sum of   squares vs. K-value")   

mat_varexp = matrix(nrow = 10,   ncol = 2)   
# Percentage of Variance   explained   
print(paste("Percent.   variance explained"))   
for (i in (1:10)){   
  km_fit =   kmeans(x_iris,centers = i,iter.max = 300 )   
  var_exp =   km_fit$betweenss/km_fit$totss   
  print(paste("K-Value",i,",Percent   var explained",round(var_exp,4)))   
  mat_varexp[i,1]=i   
  mat_varexp[i,2]=var_exp   
}   

plot(mat_varexp[,1],mat_varexp[,2],type   = 'o',xlab = "K_Value",ylab = "Percent Var explained")   
title("Avg. within sum of   squares vs. K-value") 

主成分分析 - PCA

主成分分析（PCA）是一种具有多种用途的降维技术。PCA 通过将数据投影到低维子空间来减少数据集的维度。例如，可以通过将点投影到一条线上来减少二维数据集。数据集中的每个实例将由单个值表示，而不是一对值。以类似的方式，可以通过将变量投影到平面上将三维数据集减少到二维。PCA 具有以下用途：

• 缓解维度灾难

• 在压缩数据的同时，尽量减少信息丢失。

• 主成分将在监督学习的下一阶段中进一步应用，如随机森林、提升方法等。

• 理解具有数百个维度的数据结构可能很困难，因此，通过将维度减少到二维或三维，可以更容易地可视化观察结果。

主成分分析（PCA）可以通过以下机械支架的示意图来轻松解释，该图已在机械工程课程的机械制图模块中绘制。图的左侧展示了组件的顶部视图、正面视图和侧面视图。而右侧则绘制了等轴测视图，其中使用了一张图像来可视化组件的外观。所以，可以想象，左侧的图像是实际的变量，右侧则是第一个主成分，其中捕获了大部分方差。

最终，通过旋转轴向，三张图像被替换为一张图像。实际上，我们在 PCA 分析中应用了相同的技术。

主成分工作方法在以下示例中进行了说明，实际数据在二维空间中展示，其中使用X和Y轴来绘制数据。主成分是捕捉数据最大变异性的部分。

下图展示了拟合主成分后的效果。第一个主成分涵盖了数据中的最大方差，第二个主成分与第一个主成分正交，正如我们所知，所有主成分彼此正交。我们可以仅用第一个主成分来表示整个数据。事实上，这就是用更少的维度表示数据的优势，不仅可以节省空间，还能抓取数据中的最大方差，这在下一阶段的监督学习中可以得到利用。这就是计算主成分的核心优势。

特征向量和特征值在线性代数、物理学、力学等领域具有重要意义。在学习主成分分析（PCA）时，刷新对特征向量和特征值的基础知识是必要的。特征向量是线性变换作用下仅通过 拉伸/压缩 和/或 翻转 的轴（方向）；而特征值则告诉你压缩发生的倍数。换句话说，线性变换的特征向量是一个非零向量，在应用该线性变换时，其方向保持不变。

更正式地说，A 是从向量空间到 的线性变换，如果 是 的标量倍数，则 是 A 的特征向量。该条件可以写为以下方程：

在前面的方程中， 是一个特征向量，A 是一个方阵，λ 是一个标量，称为特征值。特征向量的方向在被 A 变换后保持不变，只有其大小发生了变化，这一变化由特征值表示。换句话说，将一个矩阵乘以其特征向量等同于对特征向量进行缩放，这是原始矩阵的紧凑表示。下图展示了特征向量和特征值在二维空间中的图形表示：

以下示例描述了如何从方阵及其理解中计算特征向量和特征值。请注意，特征向量和特征值只能针对方阵（行列数相同的矩阵）进行计算。

回想一下方程，即 A 与任何 A 的特征向量的乘积必须等于特征向量与特征值的大小相乘：

特征方程表明矩阵的行列式，即数据矩阵与单位矩阵和特征值的乘积之差为0。

前述矩阵的两个特征值均为*-2*。我们可以使用特征值来替代方程中的特征向量：

将特征值代入前述方程，我们将得到以下公式：

前述方程可以重写为以下方程组：

这个方程表明它可以有多个特征向量的解，我们可以用任何满足前述方程的值进行替代以验证方程。在这里，我们使用了向量*[1 1]*进行验证，似乎已被证明。

PCA 需要单位特征向量进行计算，因此我们需要用范数除以特征向量，或者我们需要对特征向量进行归一化处理。二范数方程如下所示：

输出向量的范数计算如下：

单位特征向量如下所示：

PCA 从基本原理出发的工作方法

PCA 工作方法在以下示例数据中描述，每个实例或数据点有两个维度。这里的目标是将二维数据降维为一维（也称为主成分）：

实例	X	Y
1	0.72	0.13
2	0.18	0.23
3	2.50	2.30
4	0.45	0.16
5	0.04	0.44
6	0.13	0.24
7	0.30	0.03
8	2.65	2.10
9	0.91	0.91
10	0.46	0.32
列均值	0.83	0.69

第一步，在进行任何分析之前，是从所有观察值中减去均值，这样可以去除变量的尺度因素，并使它们在各维度之间更加统一。

X	Y
0.72 - 0.83 = -0.12	0.13 - 0.69 = - 0.55
0.18 - 0.83 = -0.65	0.23 - 0.69 = - 0.46
2.50 - 0.83 = 1.67	2.30 - 0.69 = 1.61
0.45 - 0.83 = -0.38	0.16 - 0.69 = - 0.52
0.04 - 0.83 = -0.80	0.44 - 0.69 = - 0.25
0.13 - 0.83 = -0.71	0.24 - 0.69 = - 0.45
0.30 - 0.83 = -0.53	0.03 - 0.69 = - 0.66
2.65 - 0.83 = 1.82	2.10 - 0.69 = 1.41
0.91 - 0.83 = 0.07	0.91 - 0.69 = 0.23
0.46 - 0.83 = -0.37	0.32 - 0.69 = -0.36

主成分可以通过两种不同的技术进行计算：

• 数据的协方差矩阵

• 奇异值分解

我们将在下一节中介绍奇异值分解技术。在本节中，我们将使用协方差矩阵方法求解特征向量和特征值。

协方差是衡量两个变量共同变化程度的指标，也是衡量两个变量集合之间相关性强度的度量。如果两个变量的协方差为零，我们可以得出结论，说明这两个变量集合之间没有任何相关性。协方差的公式如下：

下面的公式展示了 X 和 Y 变量的样本协方差计算。然而，它是一个 2 x 2 的矩阵，表示整个协方差矩阵（此外，它是一个方阵）。

由于协方差矩阵是方阵，我们可以从中计算特征向量和特征值。你可以参考前面章节中解释的方法。

通过求解上述方程，我们可以获得特征向量和特征值，如下所示：

前述结果可以通过以下 Python 语法获得：

>>> import numpy as np
>>> w, v = np.linalg.eig(np.array([[ 0.91335 ,0.75969 ],[ 0.75969,0.69702]]))
\>>> print ("\nEigen Values\n", w)
>>> print ("\nEigen Vectors\n", v)

一旦获得特征向量和特征值，我们可以将数据投影到主成分上。第一个特征向量具有最大的特征值，是第一个主成分，因为我们希望将原始的 2D 数据压缩成 1D 数据。

从上述结果中，我们可以看到原始 2D 数据的第一个主成分的 1D 投影。此外，特征值 1.5725 表明主成分解释了原始变量 57% 的方差。在多维数据的情况下，一般的经验法则是选择特征值或主成分的值大于某个阈值作为投影的依据。

使用 scikit-learn 对手写数字应用 PCA

PCA 示例已通过来自 scikit-learn 数据集的手写数字示例进行说明，其中手写数字从 0 到 9，并且其相应的 64 个特征（8 x 8 矩阵）表示像素强度。这里的核心思想是将原始的 64 维特征尽可能压缩到更少的维度：

# PCA - Principal Component Analysis 
>>> import matplotlib.pyplot as plt 
>>> from sklearn.decomposition import PCA 
>>> from sklearn.datasets import load_digits 

>>> digits = load_digits() 
>>> X = digits.data 
>>> y = digits.target 

>>> print (digits.data[0].reshape(8,8)) 

使用 plt.show 函数绘制图形：

>>> plt.matshow(digits.images[0])  
>>> plt.show()  

在执行 PCA 之前，建议对输入数据进行缩放，以消除因数据维度不同而可能产生的问题。例如，在对客户数据应用 PCA 时，客户的薪水维度要大于客户的年龄维度。因此，如果我们没有将所有变量缩放到相同的维度，一个变量会解释整个变化，而不是其实际的影响。在下面的代码中，我们已对所有列分别进行了缩放：

>>> from sklearn.preprocessing import scale 
>>> X_scale = scale(X,axis=0)

在下面的代码中，我们使用了两个主成分，以便将性能表示在 2D 图上。在后续部分，我们还应用了 3D。

>>> pca = PCA(n_components=2) 
>>> reduced_X = pca.fit_transform(X_scale) 

>>> zero_x, zero_y = [],[] ; one_x, one_y = [],[] 
>>> two_x,two_y = [],[]; three_x, three_y = [],[] 
>>> four_x,four_y = [],[]; five_x,five_y = [],[] 
>>> six_x,six_y = [],[]; seven_x,seven_y = [],[] 
>>> eight_x,eight_y = [],[]; nine_x,nine_y = [],[] 

在下面的代码部分中，我们将相关的主成分分别附加到每个数字上，以便我们可以创建所有 10 个数字的散点图：

>>> for i in range(len(reduced_X)): 
...     if y[i] == 0: 
...         zero_x.append(reduced_X[i][0]) 
...         zero_y.append(reduced_X[i][1]) 

...     elif y[i] == 1: 
...         one_x.append(reduced_X[i][0]) 
...         one_y.append(reduced_X[i][1]) 

...     elif y[i] == 2: 
...         two_x.append(reduced_X[i][0]) 
...         two_y.append(reduced_X[i][1]) 

...     elif y[i] == 3: 
...         three_x.append(reduced_X[i][0]) 
...         three_y.append(reduced_X[i][1]) 

...     elif y[i] == 4: 
...         four_x.append(reduced_X[i][0]) 
...         four_y.append(reduced_X[i][1]) 

...     elif y[i] == 5: 
...         five_x.append(reduced_X[i][0]) 
...         five_y.append(reduced_X[i][1]) 

...     elif y[i] == 6: 
...         six_x.append(reduced_X[i][0]) 
...         six_y.append(reduced_X[i][1]) 

...     elif y[i] == 7: 
...         seven_x.append(reduced_X[i][0]) 
...         seven_y.append(reduced_X[i][1]) 

...     elif y[i] == 8: 
...         eight_x.append(reduced_X[i][0]) 
...         eight_y.append(reduced_X[i][1]) 

...     elif y[i] == 9: 
...         nine_x.append(reduced_X[i][0]) 
...         nine_y.append(reduced_X[i][1]) 
>>> zero = plt.scatter(zero_x, zero_y, c='r', marker='x',label='zero') 
>>> one = plt.scatter(one_x, one_y, c='g', marker='+') 
>>> two = plt.scatter(two_x, two_y, c='b', marker='s') 

>>> three = plt.scatter(three_x, three_y, c='m', marker='*') 
>>> four = plt.scatter(four_x, four_y, c='c', marker='h') 
>>> five = plt.scatter(five_x, five_y, c='r', marker='D') 

>>> six = plt.scatter(six_x, six_y, c='y', marker='8') 
>>> seven = plt.scatter(seven_x, seven_y, c='k', marker='*') 
>>> eight = plt.scatter(eight_x, eight_y, c='r', marker='x') 

>>> nine = plt.scatter(nine_x, nine_y, c='b', marker='D') 

>>> plt.legend((zero,one,two,three,four,five,six,seven,eight,nine), 
...            ('zero','one','two','three','four','five','six', 'seven','eight','nine'), 
...            scatterpoints=1, 
...            loc='lower left', 
...            ncol=3, 
...            fontsize=10) 

>>> plt.xlabel('PC 1') 
>>> plt.ylabel('PC 2') 

>>> plt.show() 

尽管前面的图看起来有些杂乱，但它确实提供了一些关于数字彼此之间的远近的想法。我们可以得出结论，数字6和8非常相似，而数字4和7则远离中心组，等等。然而，我们也应该尝试使用更多的主成分，因为有时我们可能无法在二维中完全表示每个变化。

在下面的代码中，我们应用了三个主成分，以便可以在 3D 空间中更好地查看数据。这个过程与两个主成分非常相似，除了为每个数字创建一个额外的维度（X、Y 和 Z）。

# 3-Dimensional data 
>>> pca_3d = PCA(n_components=3) 
>>> reduced_X3D = pca_3d.fit_transform(X_scale) 

>>> zero_x, zero_y,zero_z = [],[],[] ; one_x, one_y,one_z = [],[],[] 
>>> two_x,two_y,two_z = [],[],[]; three_x, three_y,three_z = [],[],[] 
>>> four_x,four_y,four_z = [],[],[]; five_x,five_y,five_z = [],[],[] 
>>> six_x,six_y,six_z = [],[],[]; seven_x,seven_y,seven_z = [],[],[] 
>>> eight_x,eight_y,eight_z = [],[],[]; nine_x,nine_y,nine_z = [],[],[] 

>>> for i in range(len(reduced_X3D)): 

...     if y[i]==10: 
...         continue  

...     elif y[i] == 0: 
...         zero_x.append(reduced_X3D[i][0]) 
...         zero_y.append(reduced_X3D[i][1]) 
...         zero_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 1: 
...         one_x.append(reduced_X3D[i][0]) 
...         one_y.append(reduced_X3D[i][1]) 
...         one_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 2: 
...         two_x.append(reduced_X3D[i][0]) 
...         two_y.append(reduced_X3D[i][1]) 
...         two_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 3: 
...         three_x.append(reduced_X3D[i][0]) 
...         three_y.append(reduced_X3D[i][1]) 
...         three_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 4: 
...         four_x.append(reduced_X3D[i][0]) 
...         four_y.append(reduced_X3D[i][1]) 
...         four_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 5: 
...         five_x.append(reduced_X3D[i][0]) 
...         five_y.append(reduced_X3D[i][1]) 
...         five_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 6: 
...         six_x.append(reduced_X3D[i][0]) 
...         six_y.append(reduced_X3D[i][1]) 
...         six_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 7: 
...         seven_x.append(reduced_X3D[i][0]) 
...         seven_y.append(reduced_X3D[i][1]) 
...         seven_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 8: 
...         eight_x.append(reduced_X3D[i][0]) 
...         eight_y.append(reduced_X3D[i][1]) 
...         eight_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 9: 
...         nine_x.append(reduced_X3D[i][0]) 
...         nine_y.append(reduced_X3D[i][1]) 
...         nine_z.append(reduced_X3D[i][2]) 

 # 3- Dimensional plot 
>>> from mpl_toolkits.mplot3d import Axes3D 
>>> fig = plt.figure() 
>>> ax = fig.add_subplot(111, projection='3d') 

>>> ax.scatter(zero_x, zero_y,zero_z, c='r', marker='x',label='zero') 
>>> ax.scatter(one_x, one_y,one_z, c='g', marker='+',label='one') 
>>> ax.scatter(two_x, two_y,two_z, c='b', marker='s',label='two') 

>>> ax.scatter(three_x, three_y,three_z, c='m', marker='*',label='three') 
>>> ax.scatter(four_x, four_y,four_z, c='c', marker='h',label='four') 
>>> ax.scatter(five_x, five_y,five_z, c='r', marker='D',label='five') 

>>> ax.scatter(six_x, six_y,six_z, c='y', marker='8',label='six') 
>>> ax.scatter(seven_x, seven_y,seven_z, c='k', marker='*',label='seven') 
>>> ax.scatter(eight_x, eight_y,eight_z, c='r', marker='x',label='eight') 

>>> ax.scatter(nine_x, nine_y,nine_z, c='b', marker='D',label='nine') 
>>> ax.set_xlabel('PC 1') 
>>> ax.set_ylabel('PC 2') 
>>> ax.set_zlabel('PC 3')

>>> plt.legend(loc='upper left', numpoints=1, ncol=3, fontsize=10, bbox_to_anchor=(0, 0)) 

>>> plt.show()

与 R 图等其他软件绘图相比，matplotlib 图形有一个巨大的优势。它们是交互式的，这意味着我们可以旋转它们并从不同的角度观察它们的外观。我们鼓励读者通过旋转和探索来观察图形。在 3D 图中，我们可以看到类似的情况，并得到更多的解释。数字2位于最左边，数字0位于图形的下方。而数字4位于右上角，数字6似乎更靠近PC 1轴。通过这种方式，我们可以可视化并观察数字如何分布。在使用 4 个主成分时，我们需要使用子图并分别进行可视化。

选择提取多少主成分是无监督学习中的一个开放性问题，但有一些方法可以得到一个大致的估计。我们可以通过两种方式来确定聚类的数量：

• 检查总方差解释是否在边际上开始减少

• 总方差解释大于 80%

以下代码提供了随着主成分数目变化，所解释的总方差。随着主成分数目的增加，能够解释的方差也会增多。但挑战在于尽可能限制主成分的数量，这可以通过限制方差解释的边际增量开始下降来实现。

# Choosing number of Principal Components 
>>> max_pc = 30 

>>> pcs = [] 
>>> totexp_var = [] 

>>> for i in range(max_pc): 
...     pca = PCA(n_components=i+1) 
...     reduced_X = pca.fit_transform(X_scale) 
...     tot_var = pca.explained_variance_ratio_.sum() 
...     pcs.append(i+1) 
...     totexp_var.append(tot_var) 

>>> plt.plot(pcs,totexp_var,'r') 
>>> plt.plot(pcs,totexp_var,'bs') 
>>> plt.xlabel('No. of PCs',fontsize = 13) 
>>> plt.ylabel('Total variance explained',fontsize = 13) 

>>> plt.xticks(pcs,fontsize=13) 
>>> plt.yticks(fontsize=13) 
>>> plt.show() 

从之前的图中可以看到，总方差解释度在 10 个主成分（PCA）时略微减少；而 21 个主成分时，总方差解释度大于 80%。选择哪个值取决于业务和用户的需求。

应用于手写数字数据的 PCA 的 R 代码如下：

# PCA   
digits_data = read.csv("digitsdata.csv")   

remove_cols = c("target")   
x_data =   digits_data[,!(names(digits_data) %in% remove_cols)]   
y_data = digits_data[,c("target")]   

# Normalizing the data   
normalize <- function(x)   {return((x - min(x)) / (max(x) - min(x)))}   
data_norm <-   as.data.frame(lapply(x_data, normalize))   
data_norm <-   replace(data_norm, is.na(data_norm), 0.0)   

# Extracting Principal   Components   
pr_out =prcomp(data_norm)   
pr_components_all = pr_out$x   

# 2- Dimensional PCA   
K_prcomps = 2   
pr_components =   pr_components_all[,1:K_prcomps]   

pr_components_df =   data.frame(pr_components)   
pr_components_df =   cbind(pr_components_df,digits_data$target)   
names(pr_components_df)[K_prcomps+1]   = "target"   

out <- split(   pr_components_df , f = pr_components_df$target )   
zero_df = out$`0`;one_df =   out$`1`;two_df = out$`2`; three_df = out$`3`; four_df = out$`4`   
five_df = out$`5`;six_df =   out$`6`;seven_df = out$`7`;eight_df = out$`8`;nine_df = out$`9`   

library(ggplot2)   
# Plotting 2-dimensional PCA   
ggplot(pr_components_df, aes(x   = PC1, y = PC2, color = factor(target,labels = c("zero","one","two",   "three","four", "five","six","seven","eight","nine"))))   +    
geom_point()+ggtitle("2-D   PCA on Digits Data") +   
labs(color = "Digtis")   

# 3- Dimensional PCA   
# Plotting 3-dimensional PCA   
K_prcomps = 3   

pr_components =   pr_components_all[,1:K_prcomps]   
pr_components_df =   data.frame(pr_components)   
pr_components_df =   cbind(pr_components_df,digits_data$target)   
names(pr_components_df)[K_prcomps+1]   = "target"   

pr_components_df$target =   as.factor(pr_components_df$target)   

out <- split(   pr_components_df , f = pr_components_df$target )   
zero_df = out$`0`;one_df =   out$`1`;two_df = out$`2`; three_df = out$`3`; four_df = out$`4`   
five_df = out$`5`;six_df =   out$`6`;seven_df = out$`7`;eight_df = out$`8`;nine_df = out$`9`   

library(scatterplot3d)   
colors <- c("darkred",   "darkseagreen4", "deeppink4", "greenyellow", "orange",   "navyblue", "red", "tan3", "steelblue1",   "slateblue")   
colors <- colors[as.numeric(pr_components_df$target)]   
s3d =   scatterplot3d(pr_components_df[,1:3], pch = 16, color=colors,   
xlab = "PC1",ylab = "PC2",zlab   = "PC3",col.grid="lightblue",main = "3-D PCA on   Digits Data")   
legend(s3d$xyz.convert(3.1,   0.1, -3.5), pch = 16, yjust=0,   
       legend =   levels(pr_components_df$target),col =colors,cex = 1.1,xjust = 0)   

# Choosing number of Principal   Components   
pr_var =pr_out$sdev ²   
pr_totvar = pr_var/sum(pr_var)   
plot(cumsum(pr_totvar), xlab="Principal   Component", ylab ="Cumilative Prop. of Var.",   ylim=c(0,1),type="b",main = "PCAs vs. Cum prop of Var   Explained") 

奇异值分解 - SVD

许多 PCA 实现使用奇异值分解来计算特征向量和特征值。SVD 通过以下方程给出：

U 的列称为数据矩阵的左奇异向量，V 的列是其右奇异向量， 的对角线条目是其奇异值。左奇异向量是协方差矩阵的特征向量， 的对角元素是协方差矩阵特征值的平方根。

在进行 SVD 之前，了解一些 SVD 的优点和重要点是明智的：

• SVD 可以应用于矩形矩阵；而特征值仅定义于方阵。通过 SVD 方法获得的特征值等价物称为奇异值，得到的与特征向量等价的向量称为奇异向量。然而，由于它们本质上是矩形的，我们需要分别为它们的维度计算左奇异向量和右奇异向量。

• 如果矩阵 A 的特征向量矩阵 P 不可逆，那么 A 没有特征分解。然而，如果 A 是一个 m x n 的实矩阵且 m > n，则 A 可以通过奇异值分解表示。

• U 和 V 都是正交矩阵，这意味着 U^T U = I（其中 I 是 m x m 的单位矩阵）或 V^T V = I（这里 I 是 n x n 的单位矩阵），两个单位矩阵的维度可能不同。

• 是一个具有 m x n 维度的非负对角矩阵。

然后，奇异值和奇异向量的计算通过以下一组方程式完成：

在第一阶段，使用方程计算奇异值/特征值。一旦获得奇异值/特征值，我们将代入计算 V 或右奇异/特征向量：

一旦我们获得了右奇异向量和对角值，我们将代入公式以获得左奇异向量 U，公式如下：

这样，我们将计算原始方程矩阵的奇异值分解。

使用 scikit-learn 对手写数字应用 SVD

可以对相同的手写数字数据应用 SVD，以进行技术的“苹果对苹果”比较。

# SVD 
>>> import matplotlib.pyplot as plt 
>>> from sklearn.datasets import load_digits 

>>> digits = load_digits() 
>>> X = digits.data 
>>> y = digits.target 

在下面的代码中，使用了 15 个奇异向量和 300 次迭代，但我们鼓励读者修改这些值并检查 SVD 的性能。我们使用了两种类型的 SVD 函数，其中randomized_svd函数提供了原始矩阵的分解，而TruncatedSVD函数则可以提供总方差解释比例。在实际应用中，用户可能不需要查看所有的分解，只需使用TruncatedSVD函数来满足实际需求。

>>> from sklearn.utils.extmath import randomized_svd 
>>> U,Sigma,VT = randomized_svd(X,n_components=15,n_iter=300,random_state=42) 

>>> import pandas as pd 
>>> VT_df = pd.DataFrame(VT) 

>>> print ("\nShape of Original Matrix:",X.shape) 
>>> print ("\nShape of Left Singular vector:",U.shape) 
>>> print ("Shape of Singular value:",Sigma.shape) 
>>> print ("Shape of Right Singular vector",VT.shape) 

从之前的截图中，我们可以看到，原始矩阵（维度为 1797 x 64）已经被分解为左奇异向量（1797 x 15）、奇异值（15 维的对角矩阵）和右奇异向量（15 x 64）。我们可以通过按顺序乘以这三种矩阵来恢复原始矩阵。

>>> n_comps = 15 
>>> from sklearn.decomposition import TruncatedSVD 
>>> svd = TruncatedSVD(n_components=n_comps, n_iter=300, random_state=42) 
>>> reduced_X = svd.fit_transform(X) 

>>> print("\nTotal Variance explained for %d singular features are %0.3f"%(n_comps, svd.explained_variance_ratio_.sum())) 

对于 15 个奇异值特征，总方差解释为 83.4%。但是读者需要更改不同的值，以决定最优值。

以下代码说明了奇异值个数变化时总方差解释的变化：

# Choosing number of Singular Values 
>>> max_singfeat = 30 
>>> singfeats = [] 
>>> totexp_var = [] 

>>> for i in range(max_singfeat): 
...     svd = TruncatedSVD(n_components=i+1, n_iter=300, random_state=42) 
...     reduced_X = svd.fit_transform(X) 
...     tot_var = svd.explained_variance_ratio_.sum() 
...     singfeats.append(i+1) 
...     totexp_var.append(tot_var) 

>>> plt.plot(singfeats,totexp_var,'r') 
>>> plt.plot(singfeats,totexp_var,'bs') 
>>> plt.xlabel('No. of Features',fontsize = 13) 
>>> plt.ylabel('Total variance explained',fontsize = 13) 

>>> plt.xticks(pcs,fontsize=13) 
>>> plt.yticks(fontsize=13) 
>>> plt.show()

从之前的图中，我们可以根据需求选择 8 个或 15 个奇异向量。

应用于手写数字数据的 SVD 的 R 代码如下：

#SVD    
library(svd)   

digits_data = read.csv("digitsdata.csv")   

remove_cols = c("target")   
x_data =   digits_data[,!(names(digits_data) %in% remove_cols)]   
y_data = digits_data[,c("target")]   

sv2 <- svd(x_data,nu=15)   

# Computing the square of the   singular values, which can be thought of as the vector of matrix energy   
# in order to pick top singular   values which preserve at least 80% of variance explained   
energy <- sv2$d ^ 2   
tot_varexp = data.frame(cumsum(energy)   / sum(energy))   

names(tot_varexp) = "cum_var_explained"   
tot_varexp$K_value =   1:nrow(tot_varexp)   

plot(tot_varexp[,2],tot_varexp[,1],type   = 'o',xlab = "K_Value",ylab = "Prop. of Var Explained")   
title("SVD - Prop. of Var   explained with K-value")    

深度自编码器

自编码神经网络是一种无监督学习算法，它通过反向传播设置目标值，使其等于输入值 y^((i)) = x^((i))。自编码器试图学习一个函数 hw,b ≈ x，这意味着它试图学习一个近似的恒等函数，使得输出 与 x 类似。

尽管试图复制恒等函数看似一个简单的任务，通过在网络上施加约束（例如限制隐藏单元的数量），我们可以发现数据中的有趣结构。假设输入的图像大小为 10 x 10 像素，具有强度值，总共有 100 个输入值，第二层（隐藏层）中的神经元数量为 50 个，最终输出层有 100 个神经元，因为我们需要将图像传递并映射到自身，在此过程中，我们将迫使网络学习输入的压缩表示，即隐藏单元激活 a^((2))  ε  R¹⁰⁰，通过这个表示，我们必须尝试重建这 100 个像素的输入 x。如果输入数据完全随机，且没有任何相关性等，那么压缩将非常困难；而如果潜在的数据具有某些相关性或可检测的结构，那么这个算法将能够发现这些相关性并紧凑地表示它们。实际上，自编码器通常最终会学习到与 PCA 非常相似的低维表示。

使用编码器-解码器架构的模型构建技术

训练自编码器模型有点复杂，因此提供了详细的说明以便更好地理解。在训练阶段，整个编码器-解码器部分会针对相同的输入进行训练，作为解码器的输出。为了实现所需的输出，特征将在中间层进行压缩，因为我们通过了收敛和发散层。一旦通过多次迭代减少错误值完成足够的训练，就可以使用训练好的编码器部分来为下一阶段的建模创建潜在特征，或进行可视化等。

在下图中，展示了一个样本。输入层和输出层有五个神经元，而中间部分的神经元数量逐渐减少。压缩层只有两个神经元，即我们希望从数据中提取的潜在维度的数量。

下图展示了如何使用训练过的编码器部分，从新的输入数据中创建潜在特征，这些特征可以用于可视化或模型的下一个阶段：

使用 Keras 应用在手写数字上的深度自编码器

深度自编码器使用相同的手写数字数据进行说明，目的是展示这种非线性方法与 PCA 和 SVD 等线性方法的区别。一般来说，非线性方法的表现要好得多，但这些方法类似黑盒模型，我们无法得知其背后的解释。本文中使用了 Keras 软件来构建深度自编码器，因为它们像乐高积木一样工作，用户可以通过不同的模型架构和参数组合进行实验，以便更好地理解：

# Deep Auto Encoders 
>>> import matplotlib.pyplot as plt 
>>> from sklearn.preprocessing import StandardScaler 
>>> from sklearn.datasets import load_digits 

>>> digits = load_digits() 
>>> X = digits.data 
>>> y = digits.target 

>>> print (X.shape) 
>>> print (y.shape) 
>>> x_vars_stdscle = StandardScaler().fit_transform(X) 
>>> print (x_vars_stdscle.shape) 

使用 Keras 的密集神经元模块构建编码器-解码器架构：

>>> from keras.layers import Input,Dense 
>>> from keras.models import Model

本文使用了 NVIDIA GTX 1060 的 GPU，并安装了cuDNN和CNMeM库，以进一步提升性能，使其比常规 GPU 性能提高 4x-5x。这些库利用了 20%的 GPU 内存，剩下的 80%内存用于处理数据。用户需要注意，如果他们使用的是内存较小的 GPU，如 3GB 到 4GB 的显卡，可能无法充分利用这些库。

读者需要考虑一个重要的点，那就是，Keras 代码的语法在 CPU 模式和 GPU 模式下是相同的。

以下几行代码是模型的核心。输入数据有 64 列，我们需要将这些列作为层的输入，因此我们给出了 64 的形状。此外，每一层的神经网络都被赋予了名称，稍后我们会在代码的后续部分解释这样做的原因。在第一隐藏层中，使用了 32 个全连接神经元，这意味着输入层的所有 64 个输入都连接到第一隐藏层的 32 个神经元。整个维度流动是 64, 32, 16, 2, 16, 32, 64。我们将输入压缩到两个神经元，以便在 2D 图上绘制分量；而如果需要绘制 3D 数据（我们将在下一部分讨论），则需要将隐藏层的层数从 2 改为 3。训练完成后，我们需要使用编码器部分来预测输出。

# 2-Dimensional Architecture 

>>> input_layer = Input(shape=(64,),name="input") 

>>> encoded = Dense(32, activation='relu',name="h1encode")(input_layer) 
>>> encoded = Dense(16, activation='relu',name="h2encode")(encoded) 
>>> encoded = Dense(2, activation='relu',name="h3latent_layer")(encoded) 

>>> decoded = Dense(16, activation='relu',name="h4decode")(encoded) 
>>> decoded = Dense(32, activation='relu',name="h5decode")(decoded) 
>>> decoded = Dense(64, activation='sigmoid',name="h6decode")(decoded) 

为了训练模型，我们需要传递架构的起点和终点。在以下代码中，我们提供了输入层为input_layer，输出为decoded，即最后一层（名称为h6decode）：

>>> autoencoder = Model(input_layer, decoded) 

使用了 Adam 优化算法来优化均方误差，因为我们希望在网络的输出层末端重现原始输入：

>>> autoencoder.compile(optimizer="adam", loss="mse") 

网络使用了 100 个周期进行训练，每个批次的批量大小为 256 个观察值。使用了 20%的验证集来检查随机选取的验证数据的准确性，以确保模型的鲁棒性，因为如果只在训练数据上训练，可能会导致过拟合问题，这在高度非线性模型中非常常见：

# Fitting Encoder-Decoder model 
>>> autoencoder.fit(x_vars_stdscle, x_vars_stdscle, epochs=100,batch_size=256, shuffle=True,validation_split= 0.2 ) 

从之前的结果可以看出，模型已经在 1,437 个训练样本上进行了训练，并在 360 个验证样本上进行了验证。从损失值来看，训练损失和验证损失分别从 1.2314 降至 0.9361 和从 1.0451 降至 0.7326。因此，我们正在朝着正确的方向前进。然而，建议读者尝试不同的架构、迭代次数、批次大小等，看看准确度还能进一步提高多少。

一旦编码器-解码器部分训练完成，我们只需要使用编码器部分来压缩输入特征，以获得压缩的潜在特征，这也是降维的核心思想！在以下代码中，我们构建了一个新的模型，包含了训练好的输入层和中间隐藏层（h3latent_layer）。这就是为每一层网络分配名称的原因。

# Extracting Encoder section of the Model for prediction of latent variables 
>>> encoder = Model(autoencoder.input,autoencoder.get_layer("h3latent_layer").output) 

Extracted encoder section of the whole model used for prediction of input variables to generate sparse 2-dimensional representation, which is being performed with the following code 

# Predicting latent variables with extracted Encoder model 
>>> reduced_X = encoder.predict(x_vars_stdscle)  

仅检查降维后的输入变量的维度，可以看到，对于所有观察值，我们可以看到两个维度或两个列向量：

 >>> print (reduced_X.shape) 

以下代码部分非常类似于 2D 主成分分析（PCA）：

>>> zero_x, zero_y = [],[] ; one_x, one_y = [],[] 
>>> two_x,two_y = [],[]; three_x, three_y = [],[] 
>>> four_x,four_y = [],[]; five_x,five_y = [],[] 
>>> six_x,six_y = [],[]; seven_x,seven_y = [],[] 
>>> eight_x,eight_y = [],[]; nine_x,nine_y = [],[] 

# For 2-Dimensional data 
>>> for i in range(len(reduced_X)): 
...     if y[i] == 0: 
...         zero_x.append(reduced_X[i][0]) 
...         zero_y.append(reduced_X[i][1]) 

...     elif y[i] == 1: 
...         one_x.append(reduced_X[i][0]) 
...         one_y.append(reduced_X[i][1]) 

...     elif y[i] == 2: 
...         two_x.append(reduced_X[i][0]) 
...         two_y.append(reduced_X[i][1]) 
 ...     elif y[i] == 3: 
...         three_x.append(reduced_X[i][0]) 
...         three_y.append(reduced_X[i][1]) 

...     elif y[i] == 4: 
...         four_x.append(reduced_X[i][0]) 
...         four_y.append(reduced_X[i][1]) 

...     elif y[i] == 5: 
...         five_x.append(reduced_X[i][0]) 
...         five_y.append(reduced_X[i][1]) 

...     elif y[i] == 6: 
...         six_x.append(reduced_X[i][0]) 
...         six_y.append(reduced_X[i][1]) 

...     elif y[i] == 7: 
...         seven_x.append(reduced_X[i][0]) 
...         seven_y.append(reduced_X[i][1]) 

...     elif y[i] == 8: 
...         eight_x.append(reduced_X[i][0]) 
 ...        eight_y.append(reduced_X[i][1]) 

 ...    elif y[i] == 9: 
 ...        nine_x.append(reduced_X[i][0]) 
 ...        nine_y.append(reduced_X[i][1]) 

>>> zero = plt.scatter(zero_x, zero_y, c='r', marker='x',label='zero') 
>>> one = plt.scatter(one_x, one_y, c='g', marker='+') 
>>> two = plt.scatter(two_x, two_y, c='b', marker='s') 

>>> three = plt.scatter(three_x, three_y, c='m', marker='*') 
>>> four = plt.scatter(four_x, four_y, c='c', marker='h') 
>>> five = plt.scatter(five_x, five_y, c='r', marker='D') 

>>> six = plt.scatter(six_x, six_y, c='y', marker='8') 
>>> seven = plt.scatter(seven_x, seven_y, c='k', marker='*') 
>>> eight = plt.scatter(eight_x, eight_y, c='r', marker='x') 

>>> nine = plt.scatter(nine_x, nine_y, c='b', marker='D') 

>>> plt.legend((zero,one,two,three,four,five,six,seven,eight,nine), 
...  ('zero','one','two','three','four','five','six','seven','eight','nine'), 
...            scatterpoints=1,loc='lower right',ncol=3,fontsize=10) 

>>> plt.xlabel('Latent Feature 1',fontsize = 13) 
>>> plt.ylabel('Latent Feature 2',fontsize = 13) 

>>> plt.show() 

从之前的图中我们可以看到，数据点已经很好地分离，但问题在于视角方向，因为这些特征并没有像主成分那样沿彼此垂直的维度变化。在深度自编码器中，我们需要将视角方向从*(0, 0)*改变，以便可视化这种非线性分类，接下来我们将在 3D 案例中详细展示这一点。

以下是 3D 潜在特征的代码。除了h3latent_layer中的值需要从2替换为3外，其余代码保持不变，因为这是编码器部分的结束，我们将在创建潜在特征时使用它，最终它将用于绘图。

# 3-Dimensional architecture 
>>> input_layer = Input(shape=(64,),name="input") 

>>> encoded = Dense(32, activation='relu',name="h1encode")(input_layer) 
>>> encoded = Dense(16, activation='relu',name="h2encode")(encoded) 
>>> encoded = Dense(3, activation='relu',name="h3latent_layer")(encoded) 

>>> decoded = Dense(16, activation='relu',name="h4decode")(encoded) 
>>> decoded = Dense(32, activation='relu',name="h5decode")(decoded) 
>>> decoded = Dense(64, activation='sigmoid',name="h6decode")(decoded) 

>>> autoencoder = Model(input_layer, decoded) 
autoencoder.compile(optimizer="adam", loss="mse") 

# Fitting Encoder-Decoder model 
>>> autoencoder.fit(x_vars_stdscle, x_vars_stdscle, epochs=100,batch_size=256, shuffle=True,validation_split= 0.2) 

从之前的结果我们可以看到，加入三维度而不是二维度后，得到的损失值比 2D 用例中的值要低。对于两个潜在因子的训练和验证损失，在 100 轮迭代后的损失分别为 0.9061 和 0.7326，而对于三个潜在因子，在 100 轮迭代后的损失为 0.8032 和 0.6464。这表明，增加一个维度可以显著减少误差。这样，读者可以改变不同的参数，确定适合降维的理想架构：

# Extracting Encoder section of the Model for prediction of latent variables 
>>> encoder = Model(autoencoder.input,autoencoder.get_layer("h3latent_layer").output) 

# Predicting latent variables with extracted Encoder model 
>>> reduced_X3D = encoder.predict(x_vars_stdscle) 

>>> zero_x, zero_y,zero_z = [],[],[] ; one_x, one_y,one_z = [],[],[] 
>>> two_x,two_y,two_z = [],[],[]; three_x, three_y,three_z = [],[],[] 
>>> four_x,four_y,four_z = [],[],[]; five_x,five_y,five_z = [],[],[] 
>>> six_x,six_y,six_z = [],[],[]; seven_x,seven_y,seven_z = [],[],[] 
>>> eight_x,eight_y,eight_z = [],[],[]; nine_x,nine_y,nine_z = [],[],[] 

>>> for i in range(len(reduced_X3D)): 

...     if y[i]==10: 
...         continue 

...     elif y[i] == 0: 
...         zero_x.append(reduced_X3D[i][0]) 
...         zero_y.append(reduced_X3D[i][1]) 
...         zero_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 1: 
...         one_x.append(reduced_X3D[i][0]) 
...         one_y.append(reduced_X3D[i][1]) 
...         one_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 2: 
...         two_x.append(reduced_X3D[i][0]) 
...         two_y.append(reduced_X3D[i][1]) 
...         two_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 3: 
...         three_x.append(reduced_X3D[i][0]) 
...         three_y.append(reduced_X3D[i][1]) 
...         three_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 4: 
...         four_x.append(reduced_X3D[i][0]) 
...         four_y.append(reduced_X3D[i][1]) 
...         four_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 5: 
...         five_x.append(reduced_X3D[i][0]) 
...         five_y.append(reduced_X3D[i][1]) 
...         five_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 6: 
...         six_x.append(reduced_X3D[i][0]) 
...         six_y.append(reduced_X3D[i][1]) 
...         six_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 7: 
...         seven_x.append(reduced_X3D[i][0]) 
...         seven_y.append(reduced_X3D[i][1]) 
...         seven_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 8: 
...         eight_x.append(reduced_X3D[i][0]) 
...         eight_y.append(reduced_X3D[i][1]) 
...         eight_z.append(reduced_X3D[i][2]) 

...     elif y[i] == 9: 
...         nine_x.append(reduced_X3D[i][0]) 
...         nine_y.append(reduced_X3D[i][1]) 
...         nine_z.append(reduced_X3D[i][2]) 

 # 3- Dimensional plot 
>>> from mpl_toolkits.mplot3d import Axes3D 
>>> fig = plt.figure() 
>>> ax = fig.add_subplot(111, projection='3d') 

>>> ax.scatter(zero_x, zero_y,zero_z, c='r', marker='x',label='zero') 
>>> ax.scatter(one_x, one_y,one_z, c='g', marker='+',label='one') 
>>> ax.scatter(two_x, two_y,two_z, c='b', marker='s',label='two') 

>>> ax.scatter(three_x, three_y,three_z, c='m', marker='*',label='three') 
>>> ax.scatter(four_x, four_y,four_z, c='c', marker='h',label='four') 
>>> ax.scatter(five_x, five_y,five_z, c='r', marker='D',label='five') 

>>> ax.scatter(six_x, six_y,six_z, c='y', marker='8',label='six') 
>>> ax.scatter(seven_x, seven_y,seven_z, c='k', marker='*',label='seven') 
>>> ax.scatter(eight_x, eight_y,eight_z, c='r', marker='x',label='eight') 
 >>> ax.scatter(nine_x, nine_y,nine_z, c='b', marker='D',label='nine') 

>>> ax.set_xlabel('Latent Feature 1',fontsize = 13) 
>>> ax.set_ylabel('Latent Feature 2',fontsize = 13) 
>>> ax.set_zlabel('Latent Feature 3',fontsize = 13) 

>>> ax.set_xlim3d(0,60) 

>>> plt.legend(loc='upper left', numpoints=1, ncol=3, fontsize=10, bbox_to_anchor=(0, 0)) 

>>> plt.show() 

深度自编码器生成的 3D 图相比于三个 PCA 图提供了更好的分类分离。这里我们得到了更好的数字分离。读者需要注意的一点是，上述图是从*(0, 0, 0)*的旋转视角看到的，因为数据分离并不是在正交平面（如 PCA）上发生的，因此我们需要从原点的视角来看才能看到这种非线性分类。

总结

在本章中，您已经学习了多种无监督学习方法，用于识别数据中的结构和模式，使用了 k-means 聚类、PCA、SVD 和深度自编码器。同时，使用鸢尾花数据解释了 k-means 聚类算法。展示了如何根据各种性能指标选择最优的 k 值。利用来自 scikit-learn 的手写数据对比了线性方法（如 PCA 和 SVD）与非线性技术和深度自编码器之间的差异。详细介绍了 PCA 和 SVD 之间的区别，以便读者理解 SVD，它可以应用于矩形矩阵，即用户数量和产品数量不一定相等的情况。最后，通过可视化，证明了深度自编码器在分离数字方面优于 PCA 和 SVD 等线性无监督学习方法。

在下一章中，我们将讨论各种强化学习方法及其在人工智能等领域的应用。

第五章：强化学习

强化学习（RL）是继监督学习和无监督学习之后的第三大机器学习领域。这些技术在近年来在人工智能应用中获得了很大的关注。在强化学习中，需要做出顺序决策，而不是一次性决策，这在某些情况下使得训练模型变得困难。在本章中，我们将涵盖强化学习中使用的各种技术，并提供实际示例支持。虽然涵盖所有主题超出了本书的范围，但我们确实在这里涵盖了这个主题的最重要基础知识，以激发读者对这一主题产生足够的热情。本章讨论的主题包括：

• 马尔可夫决策过程

• 贝尔曼方程

• 动态规划

• 蒙特卡洛方法

• 时间差分学习

强化学习基础知识

在深入研究强化学习的细节之前，我想介绍一些理解 RL 方法的各种要素所必需的基础知识。这些基础知识将出现在本章的各个部分中，我们将在需要时详细解释：

• 环境： 这是具有状态和状态之间转换机制的任何系统。例如，机器人的环境是其操作的景观或设施。

• 代理： 这是与环境交互的自动化系统。

• 状态： 环境或系统的状态是完全描述环境的变量或特征集。

• 目标或吸收状态或终止状态： 这是提供比任何其他状态更高折现累积奖励的状态。高累积奖励可以防止最佳策略在训练过程中依赖于初始状态。每当代理达到目标时，我们将完成一个回合。

• 动作： 这定义了状态之间的转换。代理负责执行或至少推荐一个动作。执行动作后，代理从环境中收集奖励（或惩罚）。

• 策略： 这定义了在环境的任何状态下要选择和执行的动作。换句话说，策略是代理的行为；它是从状态到动作的映射。策略可以是确定性的，也可以是随机的。

• 最佳策略： 这是通过训练生成的策略。它定义了 Q 学习中的模型，并且会随着任何新的回合不断更新。

• 奖励： 这量化了代理与环境的积极或消极交互。奖励通常是代理到达每个状态时获得的即时收益。

• 回报或值函数： 值函数（也称为回报）是对每个状态未来奖励的预测。这些用于评估状态的好坏，基于这一点，代理将选择/行动以选择下一个最佳状态：

• 回合：这定义了从初始状态到目标状态所需的步骤数。回合也称为试验。

• 视野：这是在最大化奖励过程中所考虑的未来步骤或动作的数量。视野可以是无限的，在这种情况下，未来的奖励会被折扣，以便策略的价值能够收敛。

• 探索与利用：强化学习（RL）是一种试错学习方法。目标是找到最佳策略；同时，要保持警觉，探索一些未知的策略。一个经典的例子就是寻宝：如果我们只是贪婪地去已知的位置（利用），就会忽视其他可能藏有宝藏的地方（探索）。通过探索未知的状态，尽管即时奖励较低且没有失去最大奖励，我们仍然可能实现更大的目标。换句话说，我们是在逃离局部最优解，以便达到全局最优解（即探索），而不是仅仅专注于即时奖励的短期目标（即利用）。以下是几个例子来解释二者的区别：

  • 餐馆选择：通过偶尔探索一些未知的餐馆，我们可能会发现比我们常去的最喜欢的餐馆更好的餐馆：

    • 利用：去你最喜欢的餐馆

    • 探索：尝试一家新的餐馆

  • 油井钻探示例：通过探索新的未开发地点，我们可能会获得比仅仅探索相同地点更有益的新见解：

    • 利用：在已知最佳地点钻探石油

    • 探索：在新地点钻探

• 状态值与状态-动作函数：在动作值中，Q 代表一个智能体在状态S下采取动作A并根据某一策略π(a|s)（即在给定状态下采取某一动作的概率）之后预期获得的回报（累计折扣奖励）。

在状态值中，值是智能体在状态s下根据策略*π(a|s)*行为所预期获得的回报。更具体地说，状态值是基于策略下各动作值的期望：

• 策略内学习与策略外学习的时间差分控制：策略外学习者独立于智能体的行动学习最优策略的值。Q 学习是一个策略外学习者。策略内学习者则学习智能体执行的策略值，包括探索步骤。

• 预测与控制问题：预测是指根据给定策略评估我的表现：也就是说，如果有人给我一个策略，我执行它后能获得多少奖励。而在控制中，问题是找到最佳策略，以便我能够最大化奖励。

• 预测：评估在给定策略下各状态的值。

对于均匀随机策略，所有状态的价值函数是多少？

• 控制： 通过找到最佳策略来优化未来。

最优价值函数是什么，如何在所有可能的策略中找到最优策略？

通常，在强化学习中，我们需要先解决预测问题，之后才能解决控制问题，因为我们需要找出所有策略，才能找出最佳或最优的策略。

• RL 智能体分类： 一个 RL 智能体包括以下一个或多个组件：

  • 策略： 智能体的行为函数（从状态到动作的映射）；策略可以是确定性的或随机的

  • 价值函数： 每个状态的好坏（或）每个状态的未来奖励预期值预测

  • 模型： 智能体对环境的表征。模型预测环境接下来会做什么：

    • 转移： p 预测下一个状态（即动态）：

• • • 奖励： R 预测下一个（即时）奖励

让我们通过基于策略和值的组合以及用以下迷宫示例来解释 RL 智能体分类中的各种可能类别。在以下迷宫中，你既有起点，也有目标；智能体需要尽快到达目标，选择一条路径以获得最大的总奖励和最小的总负奖励。这个问题主要可以通过五种类别的方式来解决：

• 基于价值

• 基于策略

• Actor critic

• 无模型

• 基于模型

类别 1 - 基于价值

价值函数看起来像图像的右侧（折扣未来奖励的总和），其中每个状态都有一个值。假设距离目标一步的状态值为-1；距离目标两步的状态值为-2。类似地，起始点的值为-16。如果智能体卡在错误的位置，值可能达到-24。事实上，智能体确实根据最佳值在网格中移动，以到达目标。例如，智能体处于值为-15 的状态。在这里，它可以选择向北或向南移动，因此由于高奖励，它选择向北移动（-14），而不是向南移动（值为-16）。通过这种方式，智能体选择它在网格中的路径，直到到达目标。

• 价值函数：在所有状态下仅定义值

• 无策略（隐式）： 没有专门的策略；策略根据每个状态的值来选择

类别 2 - 基于策略

以下图中的箭头表示智能体在这些状态中选择的下一个移动方向。例如，智能体首先向东移动，然后向北移动，沿着所有箭头直到目标被达成。这也被称为从状态到动作的映射。一旦我们有了这个映射，智能体只需要读取它并相应地行动。

• 策略：策略或箭头，通过调整这些策略来达到最大可能的未来奖励。顾名思义，只有策略被存储并优化，以最大化奖励。

• 无价值函数：状态没有对应的价值。

第三类 - Actor-Critic

在 Actor-Critic 中，我们有策略和价值函数（或价值基和策略基的结合）。这种方法融合了两者的优点：

• 策略

• 价值函数

第四类 - 无模型

在强化学习中，一个基本的区分是是否基于模型或无模型。在无模型方法中，我们并未显式地建模环境，或者我们不了解完整环境的所有动态。相反，我们直接通过策略或价值函数来获得经验，并了解策略如何影响奖励：

• 策略和/或价值函数

  • 无模型

第五类 - 基于模型

在基于模型的强化学习中，我们首先建立整个环境的动态：

• 策略和/或价值函数

• 模型

经过上述所有类别后，以下维恩图展示了强化学习智能体的分类法的整体框架。如果你拿起任何关于强化学习的论文，这些方法都可以适应这个框架中的任意部分。

顺序决策中的基本类别

顺序决策中有两种基本的类型问题：

• 强化学习（例如，自主直升机等）：

  • 环境最初是未知的

  • 智能体与环境交互并从环境中获得策略、奖励和价值

  • 智能体改善其策略

• 规划（例如，象棋、Atari 游戏等）：

  • 环境模型或完整的环境动态已知

  • 智能体通过其模型进行计算（无需任何外部交互）

  • 智能体改善其策略

  • 这些问题也被称为推理、搜索、自省等问题

尽管前述的两类可以根据具体问题结合在一起，但这基本上是两种设置类型的宏观视角。

马尔可夫决策过程与贝尔曼方程

马尔可夫决策过程（MDP）正式描述了强化学习中的环境。其定义如下：

• 环境是完全可观察的

• 当前状态完全表征过程（这意味着未来状态完全依赖于当前状态，而不是历史状态或历史值）

• 几乎所有的强化学习问题都可以形式化为 MDP（例如，最优控制主要处理连续的 MDP）

MDP 的核心思想： MDP 基于状态的简单马尔可夫性属性工作；例如，S[t+1] 完全依赖于最新的状态 S[t]，而不是任何历史依赖关系。在以下方程中，当前状态捕获了来自历史的所有相关信息，这意味着当前状态是未来的充分统计量：

这个特性可以通过自主直升机示例来直观解释：下一步，直升机将向右、向左、俯仰或滚动，等等，完全取决于直升机当前的位置，而不是五分钟前的位置。

MDP 的建模： 强化学习问题通过 MDP 的五元组(S, A, {P[sa]}, y, R)来建模世界

• S - 状态集（直升机可能的所有朝向）

• A - 动作集（可以拉动控制杆的所有可能位置的集合）

• P[sa] - 状态转移分布（或状态转移概率分布）提供从一个状态到另一个状态的转移及所需的相应概率，供马尔可夫过程使用：

• γ - 折扣因子：

• R - 奖励函数（将状态集映射为实数，可以是正数或负数）：

返回是通过折扣未来奖励计算的，直到到达终止状态为止。

贝尔曼方程在 MDP 中的应用： 贝尔曼方程用于 MDP 的数学公式化，解决这些方程可以获得环境的最优策略。贝尔曼方程也被称为动态规划方程，它是与数学优化方法——动态规划——相关的最优性所必需的条件。贝尔曼方程是线性方程，可以解出整个环境的解。然而，求解这些方程的时间复杂度是 O (n³)，当环境中的状态数非常大时，计算开销会变得非常昂贵；有时，由于环境本身非常大，探索所有状态也不可行。在这种情况下，我们需要寻找其他问题求解方法。

在贝尔曼方程中，价值函数可以分解为两部分：

• 从后续状态你将得到的即时奖励 R[t+1]

• 从那个时间步开始，你将获得的后续状态的折现值 yv(S[t+1])：

MDP 的网格世界示例： 机器人导航任务生活在以下类型的网格世界中。一个障碍物位于单元格（2,2），机器人无法穿越该单元格。我们希望机器人移动到右上角的单元格（4,3），当它到达该位置时，机器人将获得+1 的奖励。机器人应该避免进入单元格（4,2），因为如果进入该单元格，它将获得-1 的奖励。

机器人可以处于以下任何位置：

• 11 个状态 - （除了 (2,2) 这个格子，那里有一个障碍物阻挡了机器人）

• A = {N-北，S-南，E-东，W-西}

在现实世界中，机器人的运动是嘈杂的，机器人可能无法精确地移动到它被要求到达的地方。例如，它的一些轮子可能打滑，部件可能连接松动，执行器可能不正确，等等。当要求它移动 1 米时，机器人可能无法精确地移动 1 米；它可能只会移动 90-105 厘米，等等。

在一个简化的网格世界中，机器人的随机动态可以如下建模。如果我们命令机器人向北移动，机器人有 10% 的概率会被拉向左侧，10% 的概率会被拉向右侧。只有 80% 的概率它才会真正向北移动。当机器人碰到墙壁（包括障碍物）并停留在原地时，什么也不会发生：

这个网格世界示例中的每个状态都由 (x, y) 坐标表示。假设机器人位于状态 (3,1)，如果我们让它向北移动，则状态转移概率矩阵如下：

机器人停留在原地的概率为 0。

如我们所知，所有状态转移概率的总和等于 1：

奖励函数：

对于所有其他状态，有一些小的负奖励值，这意味着在网格中跑来跑去时会扣除机器人的电池或燃料消耗，这就产生了不浪费移动或时间的解决方案，同时尽可能快速地达到奖励 +1 的目标，这鼓励机器人以最少的燃料消耗尽可能快速地达到目标。

当机器人到达 +1 或 -1 状态时，世界结束。到达这些状态后将不再有任何奖励；这些状态可以称为吸收状态。这些是零成本的吸收状态，机器人会永远停留在那里。

MDP 工作模型：

• 在状态 S[0]

• 选择 a[0]

• 到达 S[1] ~ P[s0, a0]

• 选择 a[1]

• 到达 S[2] ~ P[s1, a1]

• 以此类推….

一段时间后，它会获得所有奖励并将其累加：

折扣因子模型化了一种经济应用，其中今天赚到的一美元比明天赚到的一美元更有价值。

机器人需要在一段时间内选择行动（a[0]，a[1]，a[2]，......）以最大化预期回报：

在这个过程中，一个强化学习算法学习一个策略，该策略是每个状态下的行动映射，意味着这是一个推荐的行动，机器人需要根据其所处的状态来采取行动：

网格世界的最优策略： 策略是从状态到行动的映射，这意味着如果你处于某个特定状态，你需要采取这一特定行动。以下策略是最优策略，它最大化了总回报或折扣奖励的期望值。策略总是根据当前状态来决策，而不是之前的状态，这就是马尔可夫性质：

需要注意的一个问题是位置（3,1）：最优策略显示应向左（西）走，而不是向北走，虽然向北走可能涉及的状态数较少；但是，我们也可能进入一个更危险的状态。因此，向左走可能需要更长时间，但可以安全到达目的地而不会陷入负面陷阱。这些类型的结论可以通过计算获得，虽然对人类而言不明显，但计算机在提出这些策略时非常擅长：

定义：V^π, V, π**

V^π = 对于任何给定的策略π，价值函数为 V^π : S -> R，使得 V^π (S) 是从状态 S 开始，执行π后的期望总回报

网格世界的随机策略： 以下是一个随机策略及其价值函数的示例。这个策略是一个相当糟糕的策略，具有负值。对于任何策略，我们都可以写出该策略的价值函数：

简单来说，Bellman 方程说明了当前状态的价值等于即时奖励和折扣因子应用于新状态（S')的期望总回报，这些回报根据它们进入这些状态的概率进行加权。

Bellman 方程用于求解策略的价值函数的闭式解，给定固定策略，如何求解价值函数方程。

Bellman 方程对价值函数施加了一组线性约束。事实证明，通过求解一组线性方程，我们可以在任何状态 S 下求解其价值函数。

Bellman 方程在网格世界问题中的示例：

为单元格 (3,1) 选择的策略是向北移动。然而，我们的系统存在随机性，大约 80%的时间它会朝着指定方向移动，而 20%的时间它会偏离，向左（10%）或向右（10%）偏移。

可以为网格中的所有 11 个 MDP 状态写出类似的方程。我们可以从中获得以下度量值，利用线性方程组的方法来解决所有未知值：

• 11 个方程

• 11 个未知值函数变量

• 11 个约束条件

这是解决一个n个变量与n个方程的问题，我们可以通过使用方程组轻松找到一个解决方案的确切形式，进而得到整个网格的 V (π)的精确解，网格包含了所有的状态。

动态规划

动态规划是一种通过将复杂问题分解为子问题并逐一解决它们来顺序求解问题的方法。一旦子问题解决，它就会将这些子问题的解组合起来解决原始的复杂问题。在强化学习中，动态规划是一种在环境的完美模型作为马尔可夫决策过程（MDP）下计算最优策略的方法论。

动态规划适用于具有以下两个性质的问题。事实上，MDP 满足这两个性质，这使得动态规划非常适合通过求解贝尔曼方程来解决它们：

• 最优子结构

  • 最优性原理适用

  • 最优解可以分解成子问题

• 子问题重叠

  • 子问题重复多次

  • 解可以被缓存和重用

• MDP 满足这两个性质——幸运的是！

  • 贝尔曼方程具有状态值的递归分解

  • 值函数存储并重用解决方案

然而，经典的动态规划算法在强化学习中的应用有限，原因在于它们假设了一个完美的模型，并且计算开销较大。然而，它们依然很重要，因为它们为理解强化学习领域的所有方法提供了必要的基础。

使用动态规划计算最优策略的算法

计算 MDP 最优策略的标准算法利用了动态规划，以下是相关算法，我们将在本章后续部分详细讲解：

• 值迭代算法： 一种迭代算法，其中状态值不断迭代直到达到最优值；随后，最优值被用于确定最优策略

• 策略迭代算法： 一种迭代算法，其中策略评估和策略改进交替进行，以达到最优策略

值迭代算法： 值迭代算法之所以容易计算，是因为它仅对状态值进行迭代计算。首先，我们将计算最优值函数 V，然后将这些值代入最优策略方程，以确定最优策略。为了说明问题的规模，对于 11 个可能的状态，每个状态可以有四个策略（N-北，S-南，E-东，W-西），因此总共有 11⁴ 种可能的策略。值迭代算法包括以下步骤：

1. 初始化 V(S) = 0 对于所有状态 S

2. 对每个 S，更新：

1. 通过反复计算步骤 2，我们最终会收敛到所有状态的最优值：

在算法的步骤 2 中，有两种更新值的方法

• 同步更新 - 通过执行同步更新（或贝尔曼备份操作符），我们将执行右侧计算并替换方程的左侧，如下所示：

• 异步更新 - 一次更新一个状态的值，而不是同时更新所有状态，在这种情况下，状态会按固定顺序更新（先更新状态 1，再更新状态 2，以此类推）。在收敛过程中，异步更新比同步更新稍快。

值迭代在网格世界示例中的说明： 值迭代在网格世界中的应用在下图中进行了说明，解决实际问题的完整代码会在本节末尾提供。在使用贝尔曼方程对 MDP 应用之前的值迭代算法后，我们得到了以下所有状态的最优值 V*（Gamma 值选择为 0.99）：

当我们将这些值代入到我们的策略方程时，我们得到以下的策略网格：

这里，在位置 (3,1) 我们想通过数学证明为什么最优策略建议向左（西）而不是向上（北）移动：

由于墙壁，每当机器人尝试向南（下方）移动时，它会停留在原地，因此我们为当前位置分配了 0.71 的值，概率为 0.1。

同样地，对于北，我们计算了如下的总收益：

因此，向西而非向北移动会是最优的选择，因此最优策略选择了这种方式。

策略迭代算法： 策略迭代是获得 MDP 最优策略的另一种方法，在这种方法中，策略评估和策略改进算法被反复应用，直到解决方案收敛到最优策略。策略迭代算法包括以下步骤：

1. 初始化随机策略 π

2. 重复执行以下操作，直到收敛发生

   • 对当前策略求解贝尔曼方程，以使用线性方程组获得 V^π：

• • 根据新的价值函数更新策略，通过将新值假设为最优值，使用 argmax 公式改进策略：

1. 通过重复这些步骤，值和策略将会收敛到最优值：

策略迭代通常适用于较小的问题。如果一个 MDP 的状态数量非常庞大，策略迭代会在计算上非常昂贵。因此，大型 MDP 通常使用值迭代而不是策略迭代。

如果我们在现实生活中不知道确切的状态转移概率 P[s,a] 该怎么办？

我们需要使用以下简单公式从数据中估计概率：

如果某些状态没有数据可用，导致出现 0/0 问题，我们可以从均匀分布中获取一个默认概率。

使用基本 Python 的值迭代和策略迭代算法的网格世界示例

经典的网格世界示例被用来通过动态规划说明值迭代和策略迭代，以求解 MDP 的贝尔曼方程。在以下网格中，代理从网格的西南角（1,1）位置开始，目标是向东北角（4,3）移动。一旦到达目标，代理将获得 +1 的奖励。在途中，它应该避免进入危险区（4,2），因为这会导致 -1 的负奖励。代理不能从任何方向进入有障碍物（2,2）的位置。目标区和危险区是终止状态，意味着代理会继续移动，直到到达这两个状态之一。其他所有状态的奖励为 -0.02。在这里，任务是为代理在每个状态（共 11 个状态）确定最优策略（移动方向），以使代理的总奖励最大化，或者使代理尽可能快地到达目标。代理可以朝四个方向移动：北、南、东和西。

完整的代码使用 Python 编程语言编写，带有类实现。欲了解更多，请参考 Python 中的面向对象编程，了解类、对象、构造函数等。

导入random包以生成任何 N、E、S、W 方向的动作：

>>> import random,operator

以下argmax函数根据每个状态的值计算给定状态中的最大状态：

>>> def argmax(seq, fn):
...     best = seq[0]; best_score = fn(best)
...     for x in seq:
...         x_score = fn(x)
...     if x_score > best_score:
...         best, best_score = x, x_score
...     return best

要在分量级别上添加两个向量，以下代码已被用于：

>>> def vector_add(a, b):
...     return tuple(map(operator.add, a, b))

方位提供了需要加到代理当前位置的增量值；方位可以作用于x轴或y轴：

>>> orientations = [(1,0), (0, 1), (-1, 0), (0, -1)]

以下函数用于将代理转向正确的方向，因为我们知道在每个命令下，代理大约 80%的时间会按该方向移动，10%的时间会向右移动，10%的时间会向左移动：

>>> def turn_right(orientation):
...     return orientations[orientations.index(orientation)-1]
>>> def turn_left(orientation):
...     return orientations[(orientations.index(orientation)+1) % len(orientations)]
>>> def isnumber(x):
...     return hasattr(x, '__int__')

马尔可夫决策过程在这里定义为一个类。每个 MDP 由初始位置、状态、转移模型、奖励函数和 gamma 值定义。

>>> class MDP:
... def __init__(self, init_pos, actlist, terminals, transitions={}, states=None, gamma=0.99):
...     if not (0 < gamma <= 1):
...         raise ValueError("MDP should have 0 < gamma <= 1 values")
...     if states:
...         self.states = states
...     else:
...         self.states = set()
...         self.init_pos = init_pos
...         self.actlist = actlist
...         self.terminals = terminals
...         self.transitions = transitions
...         self.gamma = gamma
...         self.reward = {}

返回状态的数值奖励：

... def R(self, state):
...     return self.reward[state]

转移模型从某个状态和动作返回每个状态的（概率，结果状态）对的列表：

... def T(self, state, action):
...     if(self.transitions == {}):
...         raise ValueError("Transition model is missing")
...     else:
...         return self.transitions[state][action]

可以在特定状态下执行的动作集：

... def actions(self, state):
...     if state in self.terminals:
...         return [None]
...     else:
...         return self.actlist

GridMDP类用于建模一个二维网格世界，其中每个状态有网格值、终端位置、初始位置和折扣值（gamma）：

>>> class GridMDP(MDP):
... def __init__(self, grid, terminals, init_pos=(0, 0), gamma=0.99):

以下代码用于反转网格，因为我们希望将第 0 行放在底部，而不是顶部：

... grid.reverse()

以下__init__命令是一个构造函数，用于在网格类中初始化参数：

... MDP.__init__(self, init_pos, actlist=orientations,
terminals=terminals, gamma=gamma)
... self.grid = grid
... self.rows = len(grid)
... self.cols = len(grid[0])
... for x in range(self.cols):
...     for y in range(self.rows):
...         self.reward[x, y] = grid[y][x]
...         if grid[y][x] is not None:
...             self.states.add((x, y))

状态转移在 80%的情况下随机朝向期望方向，10%朝左，10%朝右。这是为了模拟机器人在地板上可能滑动的随机性，等等：

... def T(self, state, action):
...     if action is None:
...         return [(0.0, state)]
...     else:
...         return [(0.8, self.go(state, action)),
...                (0.1, self.go(state, turn_right(action))),
...                (0.1, self.go(state, turn_left(action)))]

返回执行指定方向后产生的状态，前提是该状态在有效状态列表中。如果下一个状态不在列表中，例如撞墙，则代理应保持在同一状态：

... def go(self, state, direction):
...     state1 = vector_add(state, direction)
...     return state1 if state1 in self.states else state

将从(x, y)到 v 的映射转换为[[...，v，...]]网格：

... def to_grid(self, mapping):
...     return list(reversed([[mapping.get((x, y), None)
...                         for x in range(self.cols)]
...                         for y in range(self.rows)]))

将方位转换为箭头，以便更好地进行图形表示：

... def to_arrows(self, policy):
...     chars = {(1, 0): '>', (0, 1): '^', (-1, 0): '<', (0, -1):
 'v', None: '.'}
...     return self.to_grid({s: chars[a] for (s, a) in policy.items()})

以下代码用于通过值迭代求解 MDP，并返回最优状态值：

>>> def value_iteration(mdp, epsilon=0.001):
...     STSN = {s: 0 for s in mdp.states}
...     R, T, gamma = mdp.R, mdp.T, mdp.gamma
...     while True:
...         STS = STSN.copy()
...         delta = 0
...         for s in mdp.states:
...             STSN[s] = R(s) + gamma * max([sum([p * STS[s1] for 
...             (p, s1) in T(s,a)]) for a in mdp.actions(s)])
...             delta = max(delta, abs(STSN[s] - STS[s]))
...         if delta < epsilon * (1 - gamma) / gamma:
...             return STS

给定一个 MDP 和一个效用函数STS，确定最佳策略，即从状态到动作的映射：

>>> def best_policy(mdp, STS):
...     pi = {}
...     for s in mdp.states:
...         pi[s] = argmax(mdp.actions(s), lambda a: expected_utility(a, s, STS, mdp))
...     return pi

根据 MDP 和 STS，执行动作a在状态s中的预期效用：

>>> def expected_utility(a, s, STS, mdp):
...     return sum([p * STS[s1] for (p, s1) in mdp.T(s, a)])

以下代码用于通过策略迭代求解 MDP，通过交替执行策略评估和策略改进步骤：

>>> def policy_iteration(mdp):
...     STS = {s: 0 for s in mdp.states}
...     pi = {s: random.choice(mdp.actions(s)) for s in mdp.states}
...     while True:
...         STS = policy_evaluation(pi, STS, mdp)
...         unchanged = True
...         for s in mdp.states:
...             a = argmax(mdp.actions(s),lambda a: expected_utility(a, s, STS, mdp))
...             if a != pi[s]:
...                 pi[s] = a
...                 unchanged = False
...         if unchanged:
...             return pi

以下代码用于返回从每个 MDP 状态到其效用的更新效用映射U，使用近似（修改后的策略迭代）：

>>> def policy_evaluation(pi, STS, mdp, k=20):
...     R, T, gamma = mdp.R, mdp.T, mdp.gamma
 ..     for i in range(k):
...     for s in mdp.states:
...         STS[s] = R(s) + gamma * sum([p * STS[s1] for (p, s1) in T(s, pi[s])])
...     return STS

>>> def print_table(table, header=None, sep=' ', numfmt='{}'):
...     justs = ['rjust' if isnumber(x) else 'ljust' for x in table[0]]
...     if header:
...         table.insert(0, header)
...     table = [[numfmt.format(x) if isnumber(x) else x for x in row]
...             for row in table]
...     sizes = list(map(lambda seq: max(map(len, seq)),
...                      list(zip(*[map(str, row) for row in table]))))
...     for row in table:
...         print(sep.join(getattr(str(x), j)(size) for (j, size, x)
...             in zip(justs, sizes, row)))

以下是一个 4 x 3 网格环境的输入网格，呈现出代理的顺序决策问题：

>>> sequential_decision_environment = GridMDP([[-0.02, -0.02, -0.02, +1],
...                                           [-0.02, None, -0.02, -1],
...                                           [-0.02, -0.02, -0.02, -0.02]],
...                                           terminals=[(3, 2), (3, 1)])

以下代码用于在给定的顺序决策环境中执行价值迭代：

>>> value_iter = best_policy(sequential_decision_environment,value_iteration (sequential_decision_environment, .01))
>>> print("\n Optimal Policy based on Value Iteration\n")
>>> print_table(sequential_decision_environment.to_arrows(value_iter))

策略迭代的代码如下：

>>> policy_iter = policy_iteration(sequential_decision_environment)
>>> print("\n Optimal Policy based on Policy Iteration & Evaluation\n")
>>> print_table(sequential_decision_environment.to_arrows(policy_iter))

从前面的输出结果来看，两个结果表明，价值迭代和策略迭代为智能体提供了相同的最优策略，使其能在最短时间内通过网格到达目标状态。当问题规模足够大时，从计算角度考虑，选择价值迭代更为可取，因为在策略迭代中，每次迭代都需要进行两个步骤：策略评估和策略改进。

蒙特卡洛方法

使用蒙特卡洛（MC）方法，我们首先计算价值函数，并确定最优策略。在这种方法中，我们不假设对环境有完全的了解。蒙特卡洛方法只需要经验，这些经验包括来自与环境实际或模拟互动的状态、动作和回报的样本序列。从实际经验中学习非常具有意义，因为它不依赖于对环境动态的先验知识，但仍能获得最优行为。这与人类或动物从实际经验中学习而非依赖任何数学模型的方式非常相似。令人惊讶的是，在许多情况下，根据所需的概率分布生成样本经验是容易的，但获得这些分布的显式形式却不可行。

蒙特卡洛方法通过对每个回合中的样本回报进行平均来解决强化学习问题。这意味着我们假设经验被划分为多个回合，并且所有回合最终都会终止，无论选择什么动作。只有在每个回合完成后，才会估算价值并改变策略。蒙特卡洛方法是按回合逐步增量的，而不是按步骤增量（这属于在线学习，我们将在时间差学习部分进行相同讨论）。

蒙特卡洛方法通过对每个状态-动作对在整个回合中的回报进行采样和平均。然而，在同一回合内，采取某一动作后的回报依赖于后续状态中采取的动作。由于所有的动作选择都在进行学习，从早期状态的角度来看，问题变得非平稳。为了处理这种非平稳性，我们借用了动态规划中的策略迭代思想，其中，首先计算固定任意策略的价值函数；然后，再改进策略。

蒙特卡洛预测

如我们所知，蒙特卡洛方法预测给定策略的状态值函数。任何状态的值是从该状态开始的预期回报或预期累计未来折扣奖励。这些值在蒙特卡洛方法中通过简单地对访问该状态后的回报进行平均来估算。随着越来越多的值被观察到，平均值应该根据大数法则收敛到预期值。实际上，这是所有蒙特卡洛方法适用的原理。蒙特卡洛策略评估算法包括以下步骤：

1. 初始化：

1. 无限重复：

   • 使用π生成一个回合

   • 对回合中出现的每个状态s：

     • G 回报，跟随第一次出现的s

     • 将G附加到 Returns(s)

     • V(s) 平均值（Returns(s)）

蒙特卡洛预测在网格世界问题上的适用性

下图用于说明目的。然而，实际上，由于并非所有策略都能保证终止，蒙特卡洛方法不能轻易用于解决网格世界类型的问题。如果发现某个策略导致代理停留在同一状态，那么下一个回合将永远不会结束。像**（状态-动作-奖励-状态-动作）**（SARSA，我们将在本章后面讨论的时间差分学习控制方法中讲解）这样的逐步学习方法没有这个问题，因为它们在回合过程中迅速学习到这些策略是差的，并会切换到其他策略。

使用 Python 建模二十一点游戏的蒙特卡洛方法

流行赌场扑克牌游戏“二十一点”的目标是获得一组牌，这些牌的点数总和尽可能大，但不得超过 21 点。所有的面牌（国王、皇后和杰克）都算作 10 点，而 A 牌可以根据玩家的需求算作 1 点或 11 点，只有 A 牌具有这种灵活性选项。其他所有牌按面值计算。游戏开始时，庄家和玩家各发两张牌，其中庄家的其中一张牌是面朝上的，另一张是面朝下的。如果玩家从这两张牌中得到“自然 21 点”（即一张 A 牌和一张 10 点牌），那么玩家获胜，除非庄家也有自然 21 点，在这种情况下，游戏为平局。如果玩家没有自然 21 点，则可以要求继续抽牌，一张一张地抽（称为“要牌”），直到他选择停牌（不再要牌）或超过 21 点（称为“爆牌”）。如果玩家爆牌，则玩家失败；如果玩家选择停牌，则轮到庄家。庄家根据固定策略选择要牌或停牌，无法选择：通常庄家在总点数为 17 点或更高时选择停牌，低于 17 点时选择要牌。如果庄家爆牌，则玩家自动获胜。如果庄家停牌，则游戏结果将是胜利、失败或平局，取决于庄家和玩家的点数哪个更接近 21 点。

黑杰克问题可以被建模为一个周期性有限的马尔可夫决策过程（MDP），其中每一局黑杰克游戏为一个回合。每个回合的奖励分别为+1（胜利）、-1（失败）和 0（平局），这些奖励会在游戏结束时给出，游戏状态中的其余奖励则为 0，不进行折扣（gamma = 1）。因此，终端奖励也就是游戏的回报。我们从一个无限的牌堆中抽取卡牌，以确保没有可追踪的规律。以下是用 Python 语言建模的整个游戏代码。

以下代码片段灵感来自Shangtong Zhang的强化学习 Python 代码，并已获得Richard S. Sutton著作《Reinforcement Learning: An Introduction》学生的许可，在本书中发布（详细信息请见进一步阅读部分）。

以下包用于数组操作和可视化：

>>> from __future__ import print_function 
>>> import numpy as np 
>>> import matplotlib.pyplot as plt 
>>> from mpl_toolkits.mplot3d import Axes3D 

在每个回合中，玩家或庄家可以采取两种可能的行动：抽牌或停牌。这是唯一的两种可能状态：

>>> ACTION_HIT = 0 
>>> ACTION_STAND = 1   
>>> actions = [ACTION_HIT, ACTION_STAND] 

玩家策略使用 21 个值的数组来建模，因为玩家总分超过 21 时将爆掉：

>>> policyPlayer = np.zeros(22) 

>>> for i in range(12, 20): 
...     policyPlayer[i] = ACTION_HIT 

如果玩家的总分为 20 或 21 点，他将选择停牌；否则，玩家将继续抽牌：

>>> policyPlayer[20] = ACTION_STAND 
>>> policyPlayer[21] = ACTION_STAND 

玩家目标策略的函数形式：

>>> def targetPolicyPlayer(usableAcePlayer, playerSum, dealerCard): 
...     return policyPlayer[playerSum] 

玩家行为策略的函数形式：

>>> def behaviorPolicyPlayer(usableAcePlayer, playerSum, dealerCard): 
...     if np.random.binomial(1, 0.5) == 1: 
...         return ACTION_STAND 
...     return ACTION_HIT 

庄家的固定策略是持续抽牌直到达到 17 点，然后在 17 到 21 点之间停牌：

>>> policyDealer = np.zeros(22) 
>>> for i in range(12, 17): 
...     policyDealer[i] = ACTION_HIT 
>>> for i in range(17, 22): 
...     policyDealer[i] = ACTION_STAND 

以下函数用于从牌堆中抽取一张新牌，并进行替换：

>>> def getCard(): 
...     card = np.random.randint(1, 14) 
...     card = min(card, 10) 
...     return card 

让我们开始游戏吧！

>>> def play(policyPlayerFn, initialState=None, initialAction=None): 

1. 玩家总分、玩家轨迹以及玩家是否将王牌算作 11 点的情况：

...     playerSum = 0 
...     playerTrajectory = [] 
...     usableAcePlayer = False 

1. 庄家抽牌状态：

...     dealerCard1 = 0 
...     dealerCard2 = 0 
...     usableAceDealer = False 

...     if initialState is None: 

1. 生成一个随机初始状态：

...         numOfAce = 0 

1. 初始化玩家的手牌：

...         while playerSum < 12: 

1. 如果玩家的总分小于 12 点，始终要抽一张牌：

...             card = getCard() 
...             if card == 1: 
...                 numOfAce += 1 
...                 card = 11 
...                 usableAcePlayer = True 
...             playerSum += card 

1. 如果玩家的总分大于 21 点，他必须至少拥有一张王牌，但也可以有两张王牌。在这种情况下，他将会把王牌当作 1 点来计算，而不是 11 点。如果玩家只有一张王牌，那么他将不再有可用的王牌：

...         if playerSum > 21: 
...             playerSum -= 10 
...             if numOfAce == 1: 
...                 usableAcePlayer = False 

1. 初始化庄家手牌：

...         dealerCard1 = getCard() 
...         dealerCard2 = getCard() 

...     else: 
...         usableAcePlayer = initialState[0] 
...         playerSum = initialState[1] 
...         dealerCard1 = initialState[2] 
...         dealerCard2 = getCard() 

1. 初始化游戏状态：

...     state = [usableAcePlayer, playerSum, dealerCard1] 

1. 初始化庄家的总分：

...     dealerSum = 0 
...     if dealerCard1 == 1 and dealerCard2 != 1: 
...         dealerSum += 11 + dealerCard2 
...         usableAceDealer = True 
...     elif dealerCard1 != 1 and dealerCard2 == 1: 
...         dealerSum += dealerCard1 + 11 
...         usableAceDealer = True 
...     elif dealerCard1 == 1 and dealerCard2 == 1: 
...         dealerSum += 1 + 11 
...         usableAceDealer = True 
...     else: 
...         dealerSum += dealerCard1 + dealerCard2 

1. 游戏从这里开始，因为玩家需要从这里开始抽取额外的牌：

...     while True: 
...         if initialAction is not None: 
...             action = initialAction 
...             initialAction = None 
...         else: 

1. 根据当前玩家的总分选择行动：

...             action = policyPlayerFn(usableAcePlayer, playerSum, dealerCard1) 

1. 跟踪玩家轨迹以便进行重要性采样：

...         playerTrajectory.append([action, (usableAcePlayer, playerSum, dealerCard1)]) 

...         if action == ACTION_STAND: 
...             break 

1. 如果选择抽牌（hit the deck），则抽一张新牌：

...         playerSum += getCard() 

1. 如果玩家的总分大于 21 点，玩家爆掉，游戏结束，玩家获得奖励-1。不过，如果玩家手上有王牌，他可以将王牌用作 11 点来挽救游戏，否则将会失败。

...         if playerSum > 21: 
...             if usableAcePlayer == True: 
...                 playerSum -= 10 
...                 usableAcePlayer = False 
...             else: 
...                 return state, -1, playerTrajectory 

1. 现在轮到庄家行动。他将根据总分来决定是否抽牌：如果他达到 17 点，就停牌，否则继续抽牌。如果庄家也有王牌，他可以使用王牌来达到爆掉的情况，否则庄家会爆掉：

...     while True: 
...         action = policyDealer[dealerSum] 
...         if action == ACTION_STAND: 
...             break 
...         dealerSum += getCard() 
...         if dealerSum > 21: 
...             if usableAceDealer == True: 
...                 dealerSum -= 10 
...                 usableAceDealer = False 
...             else: 
...                 return state, 1, playerTrajectory 

1. 现在我们将玩家的总和与庄家的总和进行比较，决定谁在不爆掉的情况下获胜：

...     if playerSum > dealerSum: 
...         return state, 1, playerTrajectory 
...     elif playerSum == dealerSum: 
...         return state, 0, playerTrajectory 
...     else: 
...         return state, -1, playerTrajectory 

以下代码演示了使用On-Policy的蒙特卡洛采样：

>>> def monteCarloOnPolicy(nEpisodes): 
...     statesUsableAce = np.zeros((10, 10)) 
...     statesUsableAceCount = np.ones((10, 10)) 
...     statesNoUsableAce = np.zeros((10, 10)) 
...     statesNoUsableAceCount = np.ones((10, 10)) 
...     for i in range(0, nEpisodes): 
...         state, reward, _ = play(targetPolicyPlayer) 
...         state[1] -= 12 
...         state[2] -= 1 
...         if state[0]: 
...             statesUsableAceCount[state[1], state[2]] += 1 
...             statesUsableAce[state[1], state[2]] += reward 
...         else: 
...             statesNoUsableAceCount[state[1], state[2]] += 1 
...             statesNoUsableAce[state[1], state[2]] += reward 
...     return statesUsableAce / statesUsableAceCount, statesNoUsableAce / statesNoUsableAceCount 

以下代码讨论了带有探索起始的蒙特卡洛方法，其中每个状态-动作对的所有回报都会被累积并平均，不管观察到时使用的是何种策略：

>>> def monteCarloES(nEpisodes): 
...     stateActionValues = np.zeros((10, 10, 2, 2)) 
...     stateActionPairCount = np.ones((10, 10, 2, 2)) 

行为策略是贪心策略，获取平均回报(s, a)的argmax：

...     def behaviorPolicy(usableAce, playerSum, dealerCard): 
...         usableAce = int(usableAce) 
...         playerSum -= 12 
...         dealerCard -= 1 
...         return np.argmax(stateActionValues[playerSum, dealerCard, usableAce, :] 
                      / stateActionPairCount[playerSum, dealerCard, usableAce, :]) 

游戏会持续若干回合，每一回合随机初始化状态、动作，并更新状态-动作对的值：

...     for episode in range(nEpisodes): 
...         if episode % 1000 == 0: 
...             print('episode:', episode) 
...         initialState = [bool(np.random.choice([0, 1])), 
...                        np.random.choice(range(12, 22)), 
...                        np.random.choice(range(1, 11))] 
...         initialAction = np.random.choice(actions) 
...         _, reward, trajectory = play(behaviorPolicy, initialState, initialAction) 
...         for action, (usableAce, playerSum, dealerCard) in trajectory: 
...             usableAce = int(usableAce) 
...             playerSum -= 12 
...             dealerCard -= 1 

更新状态-动作对的值：

...             stateActionValues[playerSum, dealerCard, usableAce, action] += reward 
...             stateActionPairCount[playerSum, dealerCard, usableAce, action] += 1 
...     return stateActionValues / stateActionPairCount 

打印状态值：

>>> figureIndex = 0 
>>> def prettyPrint(data, tile, zlabel='reward'): 
...     global figureIndex 
...     fig = plt.figure(figureIndex) 
...     figureIndex += 1 
...     fig.suptitle(tile) 
...     ax = fig.add_subplot(111, projection='3d') 
...     x_axis = [] 
...     y_axis = [] 
...     z_axis = [] 
...     for i in range(12, 22): 
...         for j in range(1, 11): 
...             x_axis.append(i) 
...             y_axis.append(j) 
...             z_axis.append(data[i - 12, j - 1]) 
...     ax.scatter(x_axis, y_axis, z_axis,c='red') 
...     ax.set_xlabel('player sum') 
...     ax.set_ylabel('dealer showing') 
...     ax.set_zlabel(zlabel) 

使用或不使用可用 A 牌的On-Policy结果，经过 10,000 和 500,000 次迭代：

>>> def onPolicy(): 
...     statesUsableAce1, statesNoUsableAce1 = monteCarloOnPolicy(10000) 
...     statesUsableAce2, statesNoUsableAce2 = monteCarloOnPolicy(500000) 
...     prettyPrint(statesUsableAce1, 'Usable Ace & 10000 Episodes') 
...     prettyPrint(statesNoUsableAce1, 'No Usable Ace & 10000 Episodes') 
...     prettyPrint(statesUsableAce2, 'Usable Ace & 500000 Episodes') 
...     prettyPrint(statesNoUsableAce2, 'No Usable Ace & 500000 Episodes') 
...     plt.show() 

策略迭代的优化或蒙特卡洛控制：

>>> def MC_ES_optimalPolicy(): 
...     stateActionValues = monteCarloES(500000) 
...     stateValueUsableAce = np.zeros((10, 10)) 
...     stateValueNoUsableAce = np.zeros((10, 10)) 
    # get the optimal policy 
...     actionUsableAce = np.zeros((10, 10), dtype='int') 
...     actionNoUsableAce = np.zeros((10, 10), dtype='int') 
...     for i in range(10): 
...         for j in range(10): 
...             stateValueNoUsableAce[i, j] = np.max(stateActionValues[i, j, 0, :]) 
...             stateValueUsableAce[i, j] = np.max(stateActionValues[i, j, 1, :]) 
...             actionNoUsableAce[i, j] = np.argmax(stateActionValues[i, j, 0, :]) 
...             actionUsableAce[i, j] = np.argmax(stateActionValues[i, j, 1, :]) 
...     prettyPrint(stateValueUsableAce, 'Optimal state value with usable Ace') 
...     prettyPrint(stateValueNoUsableAce, 'Optimal state value with no usable Ace') 
...     prettyPrint(actionUsableAce, 'Optimal policy with usable Ace', 'Action (0 Hit, 1 Stick)') 
...     prettyPrint(actionNoUsableAce, 'Optimal policy with no usable Ace', 'Action (0 Hit, 1 Stick)') 
...     plt.show() 

# Run on-policy function 
>>> onPolicy()

从前面的图表中，我们可以得出结论：即使在低玩家总和的组合下，手中有可用 A 牌也能带来更高的奖励，而对于没有可用 A 牌的玩家，如果回报少于 20，这些值之间的差异就非常明显。

# Run Monte Carlo Control or Explored starts 
>>> MC_ES_optimalPolicy() 

从最优策略和状态值中，我们可以得出结论：如果有可用的 A 牌，我们可以打得比停牌更多，而且奖励的状态值也远高于没有 A 牌时的情况。尽管我们讨论的结果是显而易见的，但我们可以看到持有 A 牌对结果的影响程度。

时序差分学习：

时序差分（TD）学习是强化学习的核心和创新主题。TD 学习结合了蒙特卡洛（MC）和动态规划（DP）的思想。与蒙特卡洛方法类似，TD 方法可以直接从经验中学习，无需环境模型。类似于动态规划，TD 方法基于其他已学习的估计值来更新估计，而不需要等待最终结果，这与 MC 方法不同，后者只有在达到最终结果后才更新估计。

时序差分预测：

TD 和 MC 都使用经验来解决z预测问题。给定某个策略π，两者都更新它们对非终结状态S[t]的估计v，该状态在经验中发生。蒙特卡洛方法会等到访问后的回报已知后，再将该回报作为*V(S[t])*的目标。

上述方法可以被称为常数 - α MC，其中 MC 必须等到回合结束后才能确定对*V(S[t])的增量（只有到那时G[t]*才会知道）。

TD 方法只需要等到下一个时间步。到达t+1时，它们立即形成目标并使用观察到的奖励R[t+1]和估计的V(S[t+1])进行有用的更新。最简单的 TD 方法是TD(0)，其公式为：

MC 更新的目标是 G[t]，而 TD 更新的目标是 R[t+1] + y V(S[t+1])。

在下图中，我们对比了 TD 和 MC 方法。如我们在方程 TD(0) 中所写，我们使用一步的实际数据，然后使用下一个状态的价值函数的估计值。同样，我们也可以使用两步的实际数据，以便更好地理解现实，并估计第三阶段的价值函数。然而，随着步骤的增加，最终需要越来越多的数据来进行参数更新，这将消耗更多的时间。

当我们采取无限步数，直到触及终止点以更新每个回合的参数时，TD 就变成了蒙特卡罗方法。

用于估计 v 的 TD (0) 算法包含以下步骤：

1. 初始化：

1. 重复（对于每个回合）：

   • 初始化 S

   • 重复（对于回合的每一步）：

     • A <- 由 π 给定的动作，针对 S

     • 执行动作 A，观察 R,S'

     • 

     • 

2. 直到 S 是终止状态。

驾驶办公室示例用于 TD 学习

在这个简单的示例中，你每天从家里出发到办公室，试图预测早晨到达办公室需要多长时间。当你离开家时，你记录下这个时间、星期几、天气情况（是否下雨、刮风等）以及你认为相关的其他参数。例如，在星期一早晨，你在早上 8 点整离开，并估计到达办公室需要 40 分钟。到了 8:10，你注意到一位贵宾正在经过，你需要等到完整的车队离开，因此你重新估计从那时起需要 45 分钟，总共需要 55 分钟。15 分钟后，你顺利完成了高速公路部分的行程。现在你进入了一个绕行的道路，你将总行程时间的估计减少到 50 分钟。不幸的是，这时你被一堆牛车堵住了，路太窄无法通过。你不得不跟着那些牛车走，直到你转入一条侧街，在 8:50 到达了办公室。七分钟后，你到达了办公室停车场。状态、时间和预测的顺序如下：

在这个示例中，奖励是每段行程所用的时间，我们使用了折扣因子（gamma，v = 1），因此每个状态的回报就是从该状态到目的地（办公室）的实际时间。每个状态的价值是预计的行驶时间，即前表中的第二列，也被称为每个遇到的状态的当前估计值。

在前面的图示中，蒙特卡洛方法被用来绘制在事件序列中的预测总时间。箭头始终表示由常数-α MC 方法推荐的预测变化。这些是每个阶段的估计值与实际回报（57 分钟）之间的误差。在 MC 方法中，学习仅在结束后进行，因此它需要等待直到 57 分钟过去。然而，实际上，您可以在达到最终结果之前进行估计，并相应地修正您的估计。TD 也遵循相同的原理，在每个阶段它尝试进行预测并相应地修正估计。因此，TD 方法立即学习，不需要等到最终结果。事实上，这也是人类在现实生活中的预测方式。由于这些积极的特性，TD 学习被认为是强化学习中的创新方法。

SARSA 在策略上的 TD 控制

状态-动作-奖励-状态-动作（SARSA）是一个在策略的 TD 控制问题，其中策略将使用策略迭代（GPI）进行优化，仅在 TD 方法用于评估预测策略时。首先，算法学习一个 SARSA 函数。特别地，对于一种策略方法，我们估计当前行为策略 π 下所有状态（s）和动作（a）对应的 q[π] (s, a)，并使用 TD 方法学习 v[π]。

现在，我们考虑从状态-动作对到状态-动作对的转换，并学习状态-动作对的值：

此更新在每次从非终结状态 S[t] 转换后进行。如果 S[t+1] 是终结状态，则 Q (S[t+1,] A[t+1]) 定义为零。此规则使用事件五元组的每个元素（S[t]，A[t]，Rt，St[+1]，A[t+1]），它们构成从一个状态-动作对到下一个的转换。这个五元组使得算法被命名为 SARSA。

与所有策略方法一样，我们不断估计行为策略 π 的 q[π]，并同时将 π 向贪心策略调整，以符合 q[π]。SARSA 计算的算法如下：

1. 初始化：

1. 重复（对于每个回合）：

   • 初始化 S

   • 使用从 Q 中得出的策略从 S 中选择 A（例如，ε-贪心策略）

   • 重复（对于每一步的回合）：

     • 执行动作 A，观察 R, S'

     • 使用从 Q 中得出的策略从 S' 中选择 A'（例如，ε-贪心策略）

     • 

     • 

2. 直到 S 是终结状态

Q-learning - 非策略 TD 控制

Q-learning 是目前在许多强化学习问题中应用最广泛的方法。策略外的 TD 控制算法被称为 Q-learning。在这种情况下，学习到的动作-价值函数 Q 直接逼近最优的动作-价值函数，独立于所遵循的策略。这个逼近简化了算法的分析，并支持早期收敛性的证明。策略仍然有影响，因为它决定了哪些状态-动作对会被访问和更新。然而，对于正确的收敛性而言，只需要所有对都继续被更新。正如我们所知，这是一个最小要求，因为任何能够保证找到最优行为的算法在一般情况下都必须满足这个要求。收敛算法的步骤如下：

1. 初始化：

1. 重复（每个回合）：

   • 初始化 S

   • 重复（每个步骤都进行）：

     • 从 S 中选择 A，使用从 Q 中推导的策略（例如，ε - 贪婪策略）

     • 采取动作 A，观察R,S'

     • 

     • 

2. 直到S为终止状态

TD 控制的策略内和策略外的悬崖行走示例

使用悬崖行走网格世界示例来比较 SARSA 和 Q-learning，突出策略内（SARSA）和策略外（Q-learning）方法之间的区别。这是一个标准的非折扣、回合制任务，具有起始和结束目标状态，并且允许四个方向（北、南、西、东）移动。除了标记为悬崖的区域，所有的转移都使用-1 作为奖励，踩到这个区域会使智能体受到-100 的惩罚，并立即将其送回起始位置。

以下代码片段灵感来自于 Shangtong Zhang 的强化学习 Python 代码，并已获得Richard S. Sutton（《强化学习：导论》一书的著名作者）的学生的许可，在本书中发布（更多详细信息请参阅进一步阅读部分）：

# Cliff-Walking - TD learning - SARSA & Q-learning 
>>> from __future__ import print_function 
>>> import numpy as np 
>>> import matplotlib.pyplot as plt 

# Grid dimensions 
>>> GRID_HEIGHT = 4 
>>> GRID_WIDTH = 12 

# probability for exploration, step size,gamma  
>>> EPSILON = 0.1 
>>> ALPHA = 0.5 
>>> GAMMA = 1 

# all possible actions 
>>> ACTION_UP = 0; ACTION_DOWN = 1;ACTION_LEFT = 2;ACTION_RIGHT = 3 
>>> actions = [ACTION_UP, ACTION_DOWN, ACTION_LEFT, ACTION_RIGHT] 

# initial state action pair values 
>>> stateActionValues = np.zeros((GRID_HEIGHT, GRID_WIDTH, 4)) 
>>> startState = [3, 0] 
>>> goalState = [3, 11] 

# reward for each action in each state 
>>> actionRewards = np.zeros((GRID_HEIGHT, GRID_WIDTH, 4)) 
>>> actionRewards[:, :, :] = -1.0 
>>> actionRewards[2, 1:11, ACTION_DOWN] = -100.0 
>>> actionRewards[3, 0, ACTION_RIGHT] = -100.0 

# set up destinations for each action in each state 
>>> actionDestination = [] 
>>> for i in range(0, GRID_HEIGHT): 
...     actionDestination.append([]) 
...     for j in range(0, GRID_WIDTH): 
...         destinaion = dict() 
...         destinaion[ACTION_UP] = [max(i - 1, 0), j] 
...         destinaion[ACTION_LEFT] = [i, max(j - 1, 0)] 
...         destinaion[ACTION_RIGHT] = [i, min(j + 1, GRID_WIDTH - 1)] 
...         if i == 2 and 1 <= j <= 10: 
...             destinaion[ACTION_DOWN] = startState 
...         else: 
...             destinaion[ACTION_DOWN] = [min(i + 1, GRID_HEIGHT - 1), j] 
...         actionDestination[-1].append(destinaion) 
>>> actionDestination[3][0][ACTION_RIGHT] = startState 

# choose an action based on epsilon greedy algorithm 
>>> def chooseAction(state, stateActionValues): 
...     if np.random.binomial(1, EPSILON) == 1: 
...         return np.random.choice(actions) 
...     else: 
...         return np.argmax(stateActionValues[state[0], state[1], :]) 

# SARSA update 

>>> def sarsa(stateActionValues, expected=False, stepSize=ALPHA): 
...     currentState = startState 
...     currentAction = chooseAction(currentState, stateActionValues) 
...     rewards = 0.0 
...     while currentState != goalState: 

...         newState = actionDestination[currentState[0]][currentState[1]] [currentAction] 

...         newAction = chooseAction(newState, stateActionValues) 
...         reward = actionRewards[currentState[0], currentState[1], currentAction] 
...         rewards += reward 
...         if not expected: 
...             valueTarget = stateActionValues[newState[0], newState[1], newAction] 
...         else: 
...             valueTarget = 0.0 
...             actionValues = stateActionValues[newState[0], newState[1], :] 
...             bestActions = np.argwhere(actionValues == np.max(actionValues)) 
...             for action in actions: 
...                 if action in bestActions: 

...                     valueTarget += ((1.0 - EPSILON) / len(bestActions) + EPSILON / len(actions)) * stateActionValues[newState[0], newState[1], action] 

...                 else: 
...                     valueTarget += EPSILON / len(actions) * stateActionValues[newState[0], newState[1], action] 
...         valueTarget *= GAMMA 
...         stateActionValues[currentState[0], currentState[1], currentAction] += stepSize * (reward+ valueTarget - stateActionValues[currentState[0], currentState[1], currentAction]) 
...         currentState = newState 
...         currentAction = newAction 
...     return rewards 

# Q-learning update 
>>> def qlearning(stateActionValues, stepSize=ALPHA): 
...     currentState = startState 
...     rewards = 0.0 
...     while currentState != goalState: 
...         currentAction = chooseAction(currentState, stateActionValues) 
...         reward = actionRewards[currentState[0], currentState[1], currentAction] 
...         rewards += reward 
...         newState = actionDestination[currentState[0]][currentState[1]] [currentAction] 
...         stateActionValues[currentState[0], currentState[1], currentAction] += stepSize * (reward + GAMMA * np.max(stateActionValues[newState[0], newState[1], :]) - 
...             stateActionValues[currentState[0], currentState[1], currentAction]) 
...         currentState = newState 
...     return rewards 

# print optimal policy 
>>> def printOptimalPolicy(stateActionValues): 
...     optimalPolicy = [] 
...     for i in range(0, GRID_HEIGHT): 
...         optimalPolicy.append([]) 
...         for j in range(0, GRID_WIDTH): 
...             if [i, j] == goalState: 
...                 optimalPolicy[-1].append('G') 
...                 continue 
...             bestAction = np.argmax(stateActionValues[i, j, :]) 
...             if bestAction == ACTION_UP: 
...                 optimalPolicy[-1].append('U') 
...             elif bestAction == ACTION_DOWN: 
...                 optimalPolicy[-1].append('D') 
...             elif bestAction == ACTION_LEFT: 
...                 optimalPolicy[-1].append('L') 
...             elif bestAction == ACTION_RIGHT: 
...                 optimalPolicy[-1].append('R') 
...     for row in optimalPolicy: 
...         print(row) 

>>> def SARSAnQLPlot(): 
    # averaging the reward sums from 10 successive episodes 
...     averageRange = 10 

    # episodes of each run 
...     nEpisodes = 500 

    # perform 20 independent runs 
...     runs = 20 

...     rewardsSarsa = np.zeros(nEpisodes) 
...     rewardsQlearning = np.zeros(nEpisodes) 
...     for run in range(0, runs): 
...         stateActionValuesSarsa = np.copy(stateActionValues) 
...         stateActionValuesQlearning = np.copy(stateActionValues) 
...         for i in range(0, nEpisodes): 
            # cut off the value by -100 to draw the figure more elegantly 
...             rewardsSarsa[i] += max(sarsa(stateActionValuesSarsa), -100) 
...             rewardsQlearning[i] += max(qlearning(stateActionValuesQlearning), -100) 

    # averaging over independent runs 
...     rewardsSarsa /= runs 
...     rewardsQlearning /= runs 

    # averaging over successive episodes 
...     smoothedRewardsSarsa = np.copy(rewardsSarsa) 
...     smoothedRewardsQlearning = np.copy(rewardsQlearning) 
...     for i in range(averageRange, nEpisodes): 
...         smoothedRewardsSarsa[i] = np.mean(rewardsSarsa[i - averageRange: i + 1]) 
...         smoothedRewardsQlearning[i] = np.mean(rewardsQlearning[i - averageRange: i + 1]) 

    # display optimal policy 
...     print('Sarsa Optimal Policy:') 
...     printOptimalPolicy(stateActionValuesSarsa) 
...     print('Q-learning Optimal Policy:') 
...     printOptimalPolicy(stateActionValuesQlearning) 

    # draw reward curves 
...     plt.figure(1) 
...     plt.plot(smoothedRewardsSarsa, label='Sarsa') 
...     plt.plot(smoothedRewardsQlearning, label='Q-learning') 
...     plt.xlabel('Episodes') 
...     plt.ylabel('Sum of rewards during episode') 
...     plt.legend() 

# Sum of Rewards for SARSA versus Qlearning 
>>> SARSAnQLPlot() 

在初步的瞬态之后，Q-learning 学习最优策略的价值，沿着最优路径行走，在此过程中智能体会沿着悬崖边缘行走。不幸的是，这会因为ε-贪婪动作选择导致偶尔从悬崖上掉下去。而 SARSA 则考虑到动作选择，学习通过网格上部的更长且更安全的路径。尽管 Q-learning 学习最优策略的价值，但其在线表现不如 SARSA，后者学习了绕行的最安全策略。即便我们观察下图所显示的奖励总和，SARSA 在整个过程中获得的奖励总和也比 Q-learning 要少负。

进一步阅读

目前有许多经典的强化学习资源可供使用，我们鼓励读者进行学习：

• R.S. Sutton 和 A.G. Barto, 强化学习：导论，MIT Press，剑桥，美国，1998 年

• RL 课程 由 David Silver 在 YouTube 上提供：www.youtube.com/watch?v=2pWv7GOvuf0&list=PL7-jPKtc4r78-wCZcQn5IqyuWhBZ8fOxT

• 机器学习（斯坦福大学）由 Andrew NG 在 YouTube 上提供（第 16-20 讲）：www.youtube.com/watch?v=UzxYlbK2c7E&list=PLA89DCFA6ADACE599

• 强化学习算法，Csaba 著，Morgan & Claypool 出版社

• 人工智能：现代方法 第 3 版，Stuart Russell 和 Peter Norvig 著，Prentice Hall 出版

摘要

在本章中，你学习了各种强化学习技术，如马尔可夫决策过程、贝尔曼方程、动态规划、蒙特卡罗方法、时序差分学习，包括基于策略（SARSA）和非基于策略（Q-learning）的算法，并通过 Python 示例来理解其在实际中的实现。你还了解了 Q-learning 如何在许多实际应用中使用，因为这种方法通过与环境互动并通过试错学习。

最后，如果你希望全职深入学习强化学习，我们为你提供了进一步阅读的资源。我们祝你一切顺利！