R-贝叶斯模型学习指南-全-

R 贝叶斯模型学习指南（全）

原文：annas-archive.org/md5/2c729bafbaf8d63831aa1a39dabebfc5

译者：飞龙

协议：CC BY-NC-SA 4.0

前言

贝叶斯推理为使用机器学习模型从数据中学习模式并用于预测未来观察结果时处理各种不确定性提供了一个统一的框架。然而，由于涉及到的数学处理水平，对于数据科学从业者来说，学习和实现贝叶斯模型并不容易。此外，将贝叶斯方法应用于现实世界问题需要大量的计算资源。随着云计算和高性能计算技术的最新进展以及计算资源的易于获取，贝叶斯建模现在在实用应用中变得更加可行。因此，对于所有数据科学家和数据工程师来说，了解贝叶斯方法并在他们的项目中应用它们以实现更好的结果将是有益的。

本书涵盖的内容

本书全面介绍了贝叶斯机器学习模型及其实现它们的 R 包。它从对概率论和 R 编程基础知识的介绍开始，对于那些对该主题不熟悉的人来说。然后，本书涵盖了使用贝叶斯推理和 R 实现的一些最重要的机器学习方法，包括监督学习和无监督学习。每一章都以非常简单的方式对方法进行理论描述。然后，讨论相关的 R 包，并给出一些使用 UCI 机器学习仓库数据集的示例。每一章结束时，还有一些简单的练习，帮助你获得本章讨论的概念和 R 包的动手经验。章节中涵盖的尖端主题包括使用线性模型和广义线性模型的贝叶斯回归、使用逻辑回归的贝叶斯分类、使用朴素贝叶斯模型对文本数据进行分类，以及使用潜在狄利克雷分配的贝叶斯混合模型和主题建模。

最后两章致力于该领域的最新发展。其中一章讨论了深度学习，它使用一类处于人工智能前沿的神经网络模型。本书以使用 Hadoop 和 Spark 等框架在大型数据上应用贝叶斯方法作为结尾。

第一章, 介绍概率论，涵盖了概率论的基础概念，特别是那些对学习贝叶斯推理所必需的方面，这些方面以简单和连贯的方式呈现给你。

第二章，R 环境，向您介绍了 R 环境。阅读完本章后，您将学习如何将数据导入 R，选择数据子集进行分析，并使用函数和控制结构编写简单的 R 程序。此外，您还将熟悉 R 的图形功能以及一些高级功能，如循环函数。

第三章，贝叶斯推理介绍，向您介绍了贝叶斯统计框架。本章包括贝叶斯定理的描述，如先验概率和后验概率等概念，以及估计后验分布的不同方法，如 MAP 估计、蒙特卡洛模拟和变分估计。

第四章，使用贝叶斯推理进行机器学习，概述了机器学习是什么以及其一些高级任务。本章还讨论了贝叶斯推理在机器学习中的重要性，特别是在如何帮助避免重要问题，如模型过拟合以及如何选择最佳模型。

第五章，贝叶斯回归模型，介绍了在贝叶斯框架下最常见的一种监督式机器学习任务，即回归建模。它通过一个示例展示了如何使用贝叶斯回归模型获得更紧的预测置信区间。

第六章，贝叶斯分类模型，介绍了如何使用贝叶斯框架进行另一种常见的机器学习任务——分类。本章讨论了分类的两种贝叶斯模型，即朴素贝叶斯和贝叶斯逻辑回归，以及评估分类器性能的一些重要指标。

第七章，无监督学习的贝叶斯模型，介绍了无监督和半监督机器学习的概念及其贝叶斯处理方法。本章讨论了两种最重要的贝叶斯无监督模型，即贝叶斯混合模型和 LDA。

第八章，贝叶斯神经网络，介绍了机器学习的一个重要类别，即神经网络及其贝叶斯实现。神经网络模型受到人脑架构的启发，并且继续是活跃的研究和开发领域。本章还讨论了深度学习，这是神经网络领域最新的进展之一，它以显著的准确性解决了计算机视觉和自然语言处理中的许多问题。

第九章，大数据规模下的贝叶斯建模，涵盖了执行大规模贝叶斯机器学习的各种框架，如 Hadoop、Spark 和 R 的本地并行化框架。本章还讨论了如何在云服务上设置实例，例如亚马逊网络服务和微软 Azure，并在其上运行 R 程序。

你需要这本书的什么

要学习本书中提供的示例和尝试练习，你需要安装 R 编程环境的最新版本和 RStudio IDE。除此之外，你还需要分别安装本书每一章中提到的特定 R 包。

这本书面向的对象

这本书旨在帮助数据科学家分析大型数据集以生成见解，以及帮助基于机器学习开发平台、解决方案或应用的数据工程师。尽管许多数据科学实践者对机器学习技术非常熟悉，并且对 R 也很熟悉，但他们可能不了解贝叶斯推断及其优点。因此，这本书对经验丰富的数据科学家和数据工程师来说也很有帮助，他们可以学习贝叶斯方法并将其纳入他们的项目中以获得更好的结果。使用这本书不需要在 R 或概率理论方面有先前的经验。

术语约定

在这本书中，你会找到许多文本样式，用于区分不同类型的信息。以下是一些这些样式的示例及其含义的解释。

文本中的代码词、数据库表名、文件夹名、文件名、文件扩展名、路径名、虚拟 URL、用户输入和 Twitter 用户名应如下所示：“第一个函数是gibbs_met。”

代码块应如下设置：

myMean ←function(x){
    s ←sum(x)
    l ←length(x)
    mean ←s/l
    mean
}
>x ←c(10,20,30,40,50)
>myMean(x)
[1]  30

任何命令行输入或输出都应如下编写：

setwd("directory path")

新术语和重要词汇以粗体显示。你会在屏幕上看到的词，例如在菜单或对话框中，在文本中显示如下：“你也可以通过点击 RStudio 的菜单栏中的会话 | 设置工作目录来设置此选项。”

注意

警告或重要注意事项将以如下框的形式出现。

提示

小技巧和窍门看起来像这样。

读者反馈

读者反馈始终欢迎。请告诉我们你对这本书的看法——你喜欢什么或不喜欢什么。读者反馈对我们来说很重要，因为它帮助我们开发出你真正能从中获得最大收益的标题。

要发送给我们一般性的反馈，只需发送电子邮件至<feedback@packtpub.com>，并在邮件主题中提及书的标题。

如果你在某个领域有专业知识，并且对撰写或参与一本书感兴趣，请参阅我们的作者指南，网址为www.packtpub.com/authors。

客户支持

现在你已经是 Packt 图书的骄傲拥有者，我们有许多事情可以帮助你从购买中获得最大收益。

下载示例代码

您可以从您在www.packtpub.com的账户下载示例代码文件，适用于您购买的所有 Packt Publishing 书籍。如果您在其他地方购买了这本书，您可以访问www.packtpub.com/support并注册，以便将文件直接通过电子邮件发送给您。

勘误

尽管我们已经尽一切努力确保内容的准确性，但错误仍然可能发生。如果您在我们的书中发现错误——可能是文本或代码中的错误——如果您能向我们报告这一点，我们将不胜感激。通过这样做，您可以节省其他读者的挫败感，并帮助我们改进本书的后续版本。如果您发现任何勘误，请通过访问www.packtpub.com/submit-errata，选择您的书籍，点击勘误提交表单链接，并输入您的勘误详情来报告它们。一旦您的勘误得到验证，您的提交将被接受，勘误将被上传到我们的网站或添加到该标题的勘误部分下的现有勘误列表中。

要查看之前提交的勘误表，请访问www.packtpub.com/books/content/support，并在搜索字段中输入书籍名称。所需信息将出现在勘误部分。

盗版

互联网上版权材料的盗版是所有媒体持续存在的问题。在 Packt，我们非常重视我们版权和许可证的保护。如果您在互联网上发现任何形式的非法复制我们的作品，请立即提供位置地址或网站名称，以便我们可以寻求补救措施。

请通过以下链接联系我们 <copyright@packtpub.com>，并提供涉嫌盗版材料的链接。

我们感谢您在保护我们的作者和我们为您提供有价值内容的能力方面的帮助。

问题

如果您对本书的任何方面有问题，您可以通过以下链接 <questions@packtpub.com> 联系我们，我们将尽力解决问题。

第一章：介绍概率论

贝叶斯推理是一种在现实世界问题中，在存在不确定性的情况下，从数据中学习变量之间关系的方法。它是概率理论的一个框架。任何对贝叶斯推理感兴趣的读者都应该对概率理论有很好的了解，以便理解和使用贝叶斯推理。本章概述了概率论，这将足以理解本书中其余章节的内容。

是皮埃尔-西蒙·拉普拉斯首先提出了一个具有数学严谨性的概率的正式定义。这个定义被称为古典定义，它陈述如下：

	偶然性的理论在于将同一类事件减少到一定数量的可能情况，也就是说，这些情况在存在方面我们可能同样不确定，并确定有利于所求概率的事件的可能情况的数量。这个数量与所有可能情况数量的比是这种概率的度量，因此它只是一个分子是有利情况数量，分母是所有可能情况数量的分数。
	--皮埃尔-西蒙·拉普拉斯，《关于概率的哲学论文》

这个定义的含义是，如果一个随机实验可以导致相互排斥且等可能的结果，那么事件的概率由以下公式给出：

这里，是事件发生的次数。

为了说明这个概念，让我们以掷骰子的简单例子来说明。如果一个骰子是公平的，那么当掷骰子时，所有面出现的概率都是相等的。那么，每个面出现的概率是 1/6。然而，当掷骰子 100 次时，由于随机波动，所有面不会以 1/6 的比例出现。每个面出现概率的估计是出现次数除以掷骰子的次数。由于分母非常大，这个比率将接近 1/6。

从长远来看，这个经典定义将不确定事件的概率视为其发生的相对频率。这也被称为概率的频率主义方法。尽管这种方法适用于大量问题，但在某些情况下，这种方法不能被使用。例如，考虑以下问题：帕塔利普特拉是古代城市的名字还是国王的名字？在这种情况下，我们对各种可能的答案有一定的信念，但这并不是基于实验结果的计数（在梵文中Putra意为儿子，因此有些人可能认为帕塔利普特拉是印度一位古代国王的名字，但它是一个城市）。

另一个例子是，2016 年美国民主党赢得选举的机会有多大？有些人可能认为是一半，有些人可能认为是三分之二。在这种情况下，概率被定义为一个人对不确定事件结果的信念程度。这被称为概率的主观定义。

经典或频率主义概率定义的一个局限性是它不能处理主观概率。正如我们将在本书后面看到的那样，贝叶斯推理是处理频率主义和主观概率解释的自然框架。

概率分布

在经典和贝叶斯方法中，概率分布函数是核心量，它捕捉了在不确定性存在时变量之间关系的所有信息。概率分布将概率值分配给随机实验结果的可测子集。涉及的变量可以是离散的或连续的，单变量的或多变量的。尽管人们使用略有不同的术语，但不同类型随机变量常用的概率分布如下：

概率质量函数（pmf）用于离散数值随机变量
分类分布用于分类随机变量
概率密度函数（pdf）用于连续随机变量

其中一个著名的分布函数是正态分布或高斯分布，它以德国著名数学家和物理学家卡尔·弗里德里希·高斯的名字命名。由于其形状，它也被称为钟形曲线。这种分布的数学形式如下：

这里，是均值或位置参数，是标准差或尺度参数（称为方差）。以下图表显示了不同位置和尺度参数值下的分布情况：

可以看到，随着均值的改变，分布峰值的定位也会改变。同样，当标准差改变时，分布的宽度也会改变。

许多自然数据集遵循正态分布，因为根据中心极限定理，任何可以表示为独立随机变量均值的随机变量都将具有正态分布。这不受该随机变量分布形式的影响，只要它们具有有限的均值和方差，并且都来自相同的原始分布。正态分布也在数据科学家中非常受欢迎，因为在许多统计推断中，如果基础分布是正态的，就可以推导出理论结果。

现在，让我们来看看正态分布的多维版本。如果随机变量是一个 N 维向量，x表示为：

然后，相应的正态分布由以下公式给出：

在这里，对应均值（也称为位置），而是一个 N x N 协方差矩阵（也称为尺度）。

为了更好地理解多维正态分布，让我们考虑二维的情况。在这种情况下，和协方差矩阵由以下给出：

在这里，和分别是沿和方向的方差，而是和之间的相关系数。以下图像展示了、和的二维正态分布图：

如果，那么二维正态分布将简化为两个一维正态分布的乘积，因为在这种情况下将变为对角线。以下和相同值但带有和的 2D 正态分布投影说明了这种情况：

在第一种情况下，x 和 y 之间的高度相关性迫使大多数数据点沿着 45 度线分布，使得分布更加各向异性；而在第二种情况下，当相关性为零时，分布更加各向同性。

在这里，我们将简要回顾一些在贝叶斯推理中使用的其他著名分布。

条件概率

通常，人们会对在问题中其他随机变量保持固定的情况下，寻找一组随机变量发生概率感兴趣。作为一个人口健康研究的例子，人们可能会感兴趣地了解在 40-50 岁年龄范围内，患有高血压和糖尿病的人患心脏病的概率是多少。这些问题可以使用条件概率来建模，条件概率定义为在另一个事件已经发生的情况下，某个事件发生的概率。更正式地说，如果我们取变量 A 和 B，这个定义可以重写如下：

同样：

下面的维恩图更清楚地解释了这个概念：

在贝叶斯推理中，我们感兴趣的是对应于多元分布的条件概率。如果表示整个随机变量集，那么在固定在某个值的情况下，的条件概率由的联合概率与的联合概率的比值给出：

在二维正态分布的情况下，感兴趣的条件概率如下：

可以证明（本章“练习”部分的第 2 题）右侧可以简化，从而得到的表达式，再次以正态分布的形式出现，均值为，方差为。

贝叶斯定理

从条件概率和的定义中，可以很容易地证明以下内容：

托马斯·贝叶斯牧师（1701–1761）使用了这个规则，并提出了他著名的贝叶斯定理，如果代表在观察B之前对随机变量A的值的初始信念（或先验概率），那么在考虑到B之后，其后验概率或信念将根据前面的方程更新。因此，贝叶斯推理本质上是在对系统进行一些观察后更新关于不确定系统的信念。从这个意义上说，这也是我们人类了解世界的方式。例如，在我们访问一个新城市之前，我们将在阅读书籍或网络后对该地有一定的先验知识。

然而，在我们到达那里不久之后，这种信念将根据我们对该地的初始经验进行更新。随着我们越来越多地探索这个新城市，我们不断地更新信念。我们将在第三章介绍贝叶斯推理中更详细地描述贝叶斯推理。

边际分布

在许多情况下，我们只对随机变量子集的概率分布感兴趣。例如，在前一节提到的冠心病问题中，如果我们想根据年龄推断人群中患有冠心病的人的概率，我们需要积分掉其他随机变量的影响，如血压和糖尿病。这被称为边际化：

或者：

注意，边际分布与条件分布非常不同。在条件概率中，我们是在寻找其他随机变量的值固定（条件）在给定值的情况下，随机变量子集的概率。在边际分布的情况下，我们通过从联合分布中积分它们（在某种意义上是平均它们的影响）来消除随机变量子集的影响。例如，在二维正态分布的情况下，关于一个变量的边际化将导致另一个变量的一维正态分布，如下所示：

这个积分的详细情况作为练习给出（本章节的练习部分中的练习 3）。

期望和协方差

已知一组随机变量的分布，在现实应用中，人们通常感兴趣的通常是估计这些随机变量的平均值以及它们之间的相关性。这些是通过以下表达式正式计算的：

例如，在二维正态分布的情况下，如果我们对寻找变量和之间的相关性感兴趣，它可以通过以下公式从联合分布中形式化计算得出：

二项分布

二项分布是一种离散分布，它给出了在n个独立试验中发生正面的概率，其中每个试验有两种可能的结果，正面或反面，正面的概率为p。每个试验被称为伯努利试验。二项分布的函数形式由以下公式给出：

在这里，表示在n次试验中发生k次正面的概率。二项分布的均值由np给出，方差由np(1-p)给出。请查看以下图表：

前面的图表显示了n的两个值；100 和 1000，对于p = 0.7的二项分布。正如你所看到的，当n变得很大时，二项分布变得尖锐峰值。可以证明，在大的n极限下，二项分布可以用具有均值np和方差np(1-p)的正态分布来近似。这是许多离散分布的一个特征，在大的n极限下，它们可以被某些连续分布近似。

Beta 分布

由表示的 Beta 分布是的幂的函数，其反射由以下公式给出：

在这里，是确定分布函数形状的参数，而是给定的 Beta 函数，它是伽玛函数比值的函数：。

Beta 分布是贝叶斯推理中一个非常重要的分布。它是二项分布、伯努利分布、负二项分布和几何分布的共轭先验概率分布（将在下一章中更精确地定义）。它用于模拟百分比和比例的随机行为。例如，Beta 分布已被用于模拟群体遗传学中的等位基因频率、项目管理中的时间分配、岩石中的矿物比例以及 HIV 传播概率的不均匀性。

伽玛分布

由表示的伽马分布是贝叶斯推断中常用的另一种常见分布。它用于模拟等待时间，如生存率。伽马分布的特殊情况是众所周知的指数分布和卡方分布。

在贝叶斯推断中，伽马分布被用作一维正态分布或如指数分布或泊松分布的速率（）的方差的共轭先验分布。

伽马分布的数学形式如下：

这里，和分别是形状和速率参数（两者都取大于零的值）。还有一种以尺度参数表示的形式，这在计量经济学中很常见。另一个相关的分布是逆伽马分布，它是根据伽马分布分布的变量的倒数分布。它主要在贝叶斯推断中用作一维正态分布方差的共轭先验分布。

狄利克雷分布

狄利克雷分布是贝塔分布的多变量类似物。它在贝叶斯推断中常用作多项分布和分类分布的共轭先验分布。主要原因在于，在狄利克雷-多项分布上实施推断技术，如吉布斯抽样，相对容易。

阶数为的狄利克雷分布定义在开维简单形上，如下所示：

这里，、和。

威沙特分布

威沙特分布是伽马分布的多变量推广。它定义在对称非负矩阵值随机变量上。在贝叶斯推断中，它被用作正态分布协方差矩阵的逆（或精度矩阵）的共轭先验分布。当我们讨论伽马分布时，我们说它被用作一维正态分布方差的共轭分布。

威沙特分布的数学定义如下：

在这里，表示维度为的矩阵的行列式，而是自由度。

Wishart 分布的一个特例是当对应于具有自由度的著名卡方分布函数时。

维基百科给出了超过 100 种常用的分布列表，这些分布通常被统计学家使用（本章参考文献部分的参考 1）。感兴趣的读者应参考这篇文章。

练习

通过使用条件概率的定义，证明任何 N 个随机变量的多元联合分布具有以下平凡的分解：
双变量正态分布由以下公式给出：

这里：

通过使用条件概率的定义，证明条件分布可以写成形式的正态分布，其中和。
通过对练习 2 中的表达式进行显式积分，证明二项正态分布的边缘化将导致单变量正态分布。

在以下表中，展示了 15 种不同鸢尾花花瓣和萼片尺寸的测量数据集（来自鸢尾花数据集，UCI 机器学习数据集仓库）。所有单位均为厘米：

萼片长度	萼片宽度	花瓣长度	花瓣宽度	花朵类别
5.1	3.5	1.4	0.2	爱丽丝·塞托萨
4.9	3	1.4	0.2	爱丽丝·塞托萨
4.7	3.2	1.3	0.2	爱丽丝·塞托萨
4.6	3.1	1.5	0.2	爱丽丝·塞托萨
5	3.6	1.4	0.2	爱丽丝·塞托萨
7	3.2	4.7	1.4	爱丽丝·弗洛里斯卡
6.4	3.2	4.5	1.5	爱丽丝·弗洛里斯卡
6.9	3.1	4.9	1.5	爱丽丝·弗洛里斯卡
5.5	2.3	4	1.3	爱丽丝·弗洛里斯卡
6.5	2.8	4.6	1.5	爱丽丝·弗洛里斯卡
6.3	3.3	6	2.5	爱丽丝·维吉妮卡
5.8	2.7	5.1	1.9	爱丽丝·维吉妮卡
7.1	3	5.9	2.1	爱丽丝·维吉妮卡
6.3	2.9	5.6	1.8	爱丽丝·维吉妮卡
6.5	3	5.8	2.2	爱丽丝·维吉妮卡

回答以下问题：

找到花瓣长度超过 5 厘米且花瓣宽度小于 3 厘米的花朵的概率是多少？
在花瓣宽度等于 0.2 厘米的条件下，找到花瓣长度小于 1.5 厘米的花朵的概率是多少？
在花朵类别为爱丽丝·弗洛里斯卡的条件下，找到萼片长度小于 6 厘米且花瓣宽度小于 1.5 厘米的花朵的概率是多少？

参考文献

概率分布列表
Feller W. 《概率论及其应用导论》。第 1 卷。威利概率与数学统计系列。1968 年。ISBN-10: 0471257087
Jayes E.T. 《概率论：科学的逻辑》。剑桥大学出版社。2003 年。ISBN-10: 0521592712
Radziwill N.M. 《使用 R 的统计学（更简单的方式）：应用统计学非正式文本》。Lapis Lucera。2015 年。ISBN-10: 0692339426

摘要

总结本章内容，我们讨论了概率论的基本要素；特别是那些对学习贝叶斯推断所必需的方面。由于篇幅限制，我们没有涵盖该主题的许多基本方面。关于这个主题有一些优秀的书籍，例如，威廉·费勒（参考 2，本章“参考文献”部分），E. T. 杰恩斯（参考 3，本章“参考文献”部分），以及 M. Radziwill（参考 4，本章“参考文献”部分）的书籍。鼓励读者阅读这些书籍，以获得对概率论更深入的理解以及它在现实生活中的应用。

在下一章中，我们将介绍 R 编程语言，这是数据分析中最受欢迎的开源框架，尤其是在贝叶斯推断方面。

第二章：R 环境

R 目前是统计计算中最受欢迎的编程环境之一。它是从贝尔实验室开发的 S 编程语言演变而来的开源语言。R 的主要创造者是新西兰奥克兰大学的两位学者，罗伯特·甘特曼和罗斯·伊哈卡。

R 流行的主要原因，除了在 GNU 通用公共许可证下的免费软件之外，还有以下几点：

R 非常易于使用。它是一种解释型语言，同时也可以用于过程式编程。
R 支持函数式和面向对象范式。它具有非常强大的图形和数据可视化能力。
通过其类似 LaTeX 的文档支持，R 可以用来制作高质量的文档。
作为开源软件，R 拥有大量贡献的包，使得几乎所有的统计建模都可以在这个环境中实现。

本章旨在提供一个基本的 R 介绍，以便任何不熟悉该语言的读者可以通过阅读本章来跟随本书的其余部分。不可能在一章中详细描述 R 语言，感兴趣的读者应参考专门为 R 编程编写的书籍。我推荐对那些主要对使用 R 进行数据分析和建模感兴趣的用户的《R 编程艺术》（本章参考文献部分的第 1 条参考文献）和《R 食谱》（本章参考文献部分的第 2 条参考文献）。对于那些对学习 R 的高级功能感兴趣的读者，例如编写复杂程序或 R 包，《高级 R》（本章参考文献部分的第 3 条参考文献）是一本优秀的书籍。

设置 R 环境和包

R 是在 GNU 开源许可证下的免费软件。R 附带一个基本包，并且还有大量用户贡献的高级分析和建模包。它还有一个名为 RStudio 的基于图形用户界面的编辑器。在本节中，我们将学习如何下载 R，在您的计算机上设置 R 环境，并编写一个简单的 R 程序。

安装 R 和 RStudio

综合 R 档案网络（CRAN）托管了 R 的所有版本和贡献的包。R for Windows 可以通过从cran.r-project.org下载基础包的二进制文件来安装；标准安装应该足够。对于 Linux 和 Mac OS X，网页提供了下载和安装软件的说明。在撰写本书时，最新版本是 3.1.2。需要安装的各种包可以从包页面单独安装。用户可以使用以下命令从 R 命令提示符安装任何包：

install.packages("package name")

安装包后，需要使用以下命令加载包才能使用：

library("package name")

一个非常实用的 R 的集成开发环境（IDE）是 RStudio。它可以从www.rstudio.com/免费下载。RStudio 在 Windows、Linux 和 Mac 平台上运行。它既有桌面版本，也有服务器版本，可以通过远程服务器上的浏览器界面编写 R 程序。安装 R 和 RStudio 后，将默认工作目录设置为所选目录很有用。RStudio 将包含 R 代码的文件读入和写入工作目录。要找出当前目录，请使用 R 命令getwd()。要将工作目录更改为你偏好的目录，请使用以下命令：

setwd("directory path")

你也可以通过点击 RStudio 菜单栏的会话 | 设置工作目录来设置此选项。

你的第一个 R 程序

让我们编写一个简单的程序来添加两个整数x和y，得到它们的和z。在 RStudio 的命令提示符中，输入以下命令并按Enter键：

>x <-2
>y <-3
>z <-x+y
>print(z)
[1]  5

现在，你可以为x和y分配不同的值，并打印z以查看z如何变化。除了print(z)，你也可以简单地输入z来打印其值。

R 中的数据管理

在我们开始使用 R 进行任何严肃的编程之前，我们需要学习如何将数据导入 R 环境以及 R 支持哪些数据类型。通常，对于特定的分析，我们可能不会使用整个数据集。因此，我们还需要学习如何为任何分析选择数据子集。本节将涵盖这些方面。

R 中的数据类型

R 有五种基本数据类型，如下所示：

整数型
数值型（实数）
复数型
字符型
逻辑型（True/False）

R 中数字的默认表示形式是双精度实数（数值型）。如果你想要显式地表示整型，需要添加后缀L。例如，在命令提示符中简单地输入1将存储为数值对象。要存储1为整型，需要输入1L。命令class(x)将给出对象x的类别（类型）。因此，在命令提示符中输入class(1)将得到答案numeric，而输入class(1L)将得到答案integer。

R 还有一个特殊的数字Inf，表示无穷大。数字NaN（非数字）用于表示未定义的值，例如 0/0。缺失值使用符号NA表示。

R 中的数据结构

R 中的数据结构可以分为同质型（所有元素包含相同的数据类型）或异质型（元素包含不同的数据类型）。此外，这些数据结构根据维数的不同而具有不同的结构：

同质型：
- 原子向量：一维
- 矩阵：二维
- 数组：N 维
异质型：
- 列表：一维
- 数据框：二维

R 中最基本的对象是向量。要在 R 提示符中创建一个大小为 10 的空整型向量，请输入以下命令：

>v <-vector("integer",10)
>v
[1]  0000000000

您可以使用以下命令将值 m 赋给向量的 n 个组件：

> v[5] <-1
> v
[1]  0000100000

读者应注意，与许多编程语言不同，R 中的数组索引从 1 开始，而不是从 0 开始。

与向量只能包含相同类型的对象不同，列表虽然与向量相似，但可以包含不同类型的对象。以下命令将创建一个包含整数、实数和字符的列表：

> l <-list(1L, 2L, 3, 4, "a", "b")
> str(l)
List of 6
$: int 1
$: int 2
$: num 3
$: num 4
$: chr "a"
$: chr "b"

在这里，我们使用了 R 中的 str() 函数，该函数显示了任何 R 对象的结构。

R 有一个特殊的函数 c()，用于将多个基本数据组合成一个向量或列表。例如，c(1,3,6,2,-1) 将生成一个包含数字 1,2,3,6,-1 的向量：

> c(1, 3, 6, 2, -1)
[1]  1 3 6 2 -1

矩阵是向量在二维上的推广。考虑以下命令：

>m <-matrix(c(1:9),nrow=3,ncol=3)

此命令将生成一个 3 x 3 大小的矩阵 m，包含从 1 到 9 的数字。

R 中用于存储数据的常见数据结构是数据框。数据框，就像列表一样，可以包含不同类型的数据（数值、整数、布尔或字符）。它本质上是一个长度相等的向量的列表。因此，它具有与矩阵相同的二维结构。数据框的长度（使用 length( ) 查找）是底层列表的长度，即数据框中的列数。有简单的命令 nrow( ) 和 ncol( ) 用于查找数据框的行数和列数。数据框的其他两个属性是 rownames( ) 和 colnames( )，可以用来设置或查找行或列的名称。

将数据导入 R

以表格形式存在的数据可以很容易地使用 read.table(…) 函数加载到 R 中。它有几个参数使导入非常灵活。其中一些有用的参数如下：

file: 文件名或完整的 URL
header: 一个逻辑值，指示文件是否包含包含变量名称的标题行
sep: 表示列分隔符的字符
row.names: 行名称的向量
col.names: 变量名称的向量
skip: 在读取数据之前要跳过的数据文件中的行数
nrows: 数据集中的行数
stringsASFactors: 一个逻辑值，指示字符变量是否可以编码为因子

对于小型数据集，可以使用 read.table("filename.txt") 而不指定其他参数；其余的 R 会自动处理。另一个有用的函数是 read.csv()，用于仅读取 CSV 文件。

除了从文本文件加载数据外，还可以通过连接外部数据库通过各种接口将数据导入 R。其中一种流行的接口是开放数据库连接（ODBC）。R 中的RODBC包通过 ODBC 接口提供对不同数据库的访问。此包包含用于与数据库连接和执行各种操作的不同函数。RODBC 包中的一些重要函数如下：

odbcConnect(dsn, uid="user_name", pwd="password"): 用于打开一个连接到已注册数据源名称 dsn 的 ODBC 数据库。
sqlFetch(channel, sqtable): 用于从 ODBC 数据库读取表到数据框。
sqlQuery(channel, query): 用于向 ODBC 数据库提交查询并返回结果。
sqlSave(channel, mydf, tablename = sqtable, append = FALSE): 用于将数据框写入或更新到 ODBC 数据库中的表（append = TRUE）。
close(channel): 用于关闭连接。在这里，channel是odbcConnect返回的连接句柄。

切割和切块数据集

在数据分析中，通常需要切割整个数据框以选择一些变量或观测值。这被称为子集化。R 有一些强大且快速的方法来完成这项工作。

要提取 R 对象的子集，可以使用以下三个运算符：

单方括号 [ ]: 返回与原始对象相同类的对象。单方括号运算符可以用于选择对象的一个或多个元素。以下是一些示例：

>x <-c(10,20,30,40,50)
>x[1:3]
[1]  10 20 30

>x[x >25]
[1]  30 40 50

>f <-x >30
>x[f]
[1]  40 50

>m <-matrix(c(1:9),nrow=3,ncol=3)
>m[1 ,] #select the entire first row
[1]  1 4 7

>m[  ,2] #select the entire second column
[1]  4 5 6

双方括号 [[ ]]: 用于从列表或数据框中提取单个元素。返回的对象不必与初始对象类型相同。以下是一些示例：
```
>y <-list("a", "b", "c", "d", "e")

>y[1]
[[1]]
[1]  "a"

>class(y[1])
[1]  "list"

>y[[1]]
[1]  "a"

>class(y[[1]])
[1]  "character"
```
美元符号 $: 用于通过名称提取列表或数据框的元素。以下是一些示例：
```
>z <-list(John = 12 ,Mary = 18,Alice = 24 ,Bob = 17 ,Tom = 21)

>z$Bob
[1] 17
```
使用负索引值: 这用于删除特定的索引或列——对应索引有一个带负号的子集。例如，要从前面的列表中删除 Mary 和 Bob，请使用以下代码：
```
> y <-z[c(-2, -4)]
> y
```

向量化操作

在 R 中，许多操作，如涉及向量和矩阵的算术运算，可以使用向量化操作非常高效地完成。例如，如果你正在添加两个向量x和y，它们的元素是并行添加的。这也使得代码更加简洁，更容易理解。例如，在代码中添加两个向量不需要for( )循环：

>x <-c(1,2,3,4,5)

>y <-c(10,20,30,40,50)

>z <-x+y

>z
[1]  11 22 33 44 55

>w <-x*y

>w
[1]  10 40 90 160 250

向量化操作的另一个非常有用的例子是在矩阵的情况下。如果X和Y是两个矩阵，则可以在 R 中以向量化的形式执行以下操作：

>X*Y  ## Element-wise multiplication
>X/Y  ## Element-wise division
>X  %*%  Y  ## Standard matrix multiplication

编写 R 程序

尽管在 R 中，许多数据分析可以通过命令提示符以交互式方式进行，但通常对于更复杂的任务，需要编写 R 脚本。如介绍中所述，R 具有函数式和面向对象编程语言的视角。在本节中，将描述 R 编程的一些标准语法。

控制结构

控制结构用于控制程序执行的流程。标准控制结构如下：

if 和 else: 测试条件
for: 循环一组语句固定次数
while: 当条件为真时循环一组语句
repeat: 执行无限循环
break: 用于中断循环的执行
next: 跳过循环的一次迭代
return: 退出函数

函数

如果想用 R 进行更严肃的编程，了解如何编写函数是至关重要的。它们使语言更加强大和优雅。R 有许多内置函数，例如 mean()、sort()、sin()、plot() 等，这些函数都是使用 R 命令编写的。

函数的定义如下：

>fname<-function(arg1,arg2,   ){
        R Expressions
   }

这里，fname 是函数的名称；arg1、arg2 等是传递给函数的参数。请注意，与在其他语言中不同，R 中的函数不需要以返回语句结束。默认情况下，函数体内执行的最后一个语句由函数返回。

一旦定义了函数，只需通过输入函数名并传递参数值来执行它：

>fname(arg1,arg2,…)

R 中函数的重要特性如下：

函数是一等公民
函数可以作为参数传递给其他函数
可以在另一个函数内部定义一个函数（嵌套）
函数的参数可以通过位置或名称进行匹配

让我们考虑一个函数的简单示例，给定一个输入向量 x，计算其平均值。要编写此函数，请通过菜单栏中的文件 | 新建文件 | R 脚本 在 RStudio 中打开一个新的 R 脚本窗口。在此 R 脚本中，输入以下代码行：

myMean <-function(x){
    s <-sum(x)
    l <-length(x)
    mean <-s/l
    mean
}

选择整个代码，然后使用 Ctrl + Enter 键来执行脚本。这完成了 myMean 函数的定义。要在命令提示符中使用此函数，请输入以下内容：

>x <-c(10,20,30,40,50)
>myMean(x)

这将生成以下结果：

>myMean(x)
[1]  30

作用域规则

在编程语言中，理解所有变量的作用域非常重要，以避免执行过程中的错误。在函数中，变量的作用域有两种类型：词法作用域和动态作用域。在词法作用域的情况下，函数中的变量值是在函数定义的环境中进行查找的。通常，这是全局环境。在动态作用域的情况下，变量的值是在函数被调用的环境中查找的（调用环境）。

R 使用词汇作用域，这使得在函数内部编写函数成为可能。以下示例说明了这一点：

>x <-0.1
>f <-function(y){
          x*y
    }
>g <-function(y){
          x<-5
          x-f(y)
    }
>g(10)
[1]  4

答案是4，因为在评估函数f时，x的值是从全局环境中获取的，即0.1，而在评估函数g时，x的值是从g的局部环境中获取的，即5。

词汇作用域有一些缺点。由于变量的值是从定义函数的环境中进行查找的，因此所有函数都必须携带指向它们各自定义环境的指针。此外，所有对象都必须在程序执行期间存储在内存中。

循环函数

通常，我们有一个包含一些对象的列表，并且想要将函数应用于列表的每个元素。例如，我们有一个包含来自n个参与者的m个问题的调查结果列表。我们希望找到每个问题的平均响应（假设所有问题都有作为数值的响应）。可以使用 R 中的mean()函数通过一个for循环遍历问题集并找到n个用户的平均值。在这种情况下，循环函数很有用，可以更紧凑地进行此类计算。这些在其他语言（如 Java）中类似于迭代器。

以下是在 R 中的标准循环函数：

lapply: 遍历列表并对每个元素评估函数
sapply: 与lapply相同，但输出形式更简单
mapply: sapply的多变量版本
apply: 在数组边缘应用函数
tapply: 将函数应用于锯齿形数组的每个单元格

lapply

lapply()函数的使用方式如下：

>lapply(X,FUN,   )

在这里，X是一个包含数据的列表或向量。FUN是需要应用于列表或向量每个元素的函数的名称。最后一个参数表示可选参数。使用lapply的结果始终是一个列表，无论输入类型如何。

例如，考虑四家公司的季度收入（以十亿美元计，非真实数据）。我们希望计算四家公司的年收入平均值，如下所示：

>X<-list(HP=c(12.5,14.3,16.1,15.4),IBM=c(22,24.5,23.7,26.2),Dell=c(8.9,9.7,10.8,11.5),Oracle=c(20.5,22.7,21.8,24.4)  )
>lapply(X,mean)
$HP
[1]  14.575

$IBM
[1]  24.1

$Dell
[1]  10.225

$Oracle
[1]  22.35

sapply

sapply()函数与lapply()类似，但增加了将输出简化为所需形式的选择。例如，sapply()可以在之前的数据集中如下使用：

> sapply(X,mean,simplify="array")
     HP       IBM       Dell      Oracle
   14.575    24.100    10.225     22.350

mapply

lapply()和sapply()函数只能有一个参数。如果您想要使用多个变量参数应用函数，则mapply()很有用。以下是它的用法：

>mapply(FUN,L1,L2,   ,Ln,SIMPLIFY=TRUE)

在这里，是需要应用函数FUN的列表。例如，考虑以下列表生成命令：

>rep(x=10,times=5)
[1] 10 10 10 10 10

在这里，rep函数重复x的值五次。假设我们想要创建一个列表，其中数字 10 出现 1 次，数字 20 出现 2 次，依此类推，我们可以使用mapply如下所示：

>mapply(rep,x=c(10,20,30,40,50),times=1:5)

apply

apply() 函数用于将函数应用于数组或矩阵的边缘。函数的形式如下：

>apply(X,MARGIN,FUN,   )

在这里，MARGIN 是一个向量，它给出了函数将被应用到的子索引。例如，在矩阵的情况下，1 表示行，2 表示列，而 c(1,2) 表示行和列。以下是一个示例来说明：

>Y <-matrix(1:9,nrow=3,ncol=3)
>Y
          [,1]        [,2]          [,3]
[1,]        1           4             7
[2,]        2           5             8
[1,]        3           6             9
>apply(Y,1,sum) #sum along the row
[1]  12 15 18
>apply(Y,2,sum) #sum along the column
[1]  6 15 24

tapply

tapply() 函数用于在向量的子集中应用函数。函数描述如下：

>tapply(X,INDEX,FUN,SIMPLIFY=TRUE)

让我们考虑之前提到的五家公司的季度营收示例：

>X<-X(HP=c(12.5,14.3,16.1,15.4),IBM=c(22,24.5,23.7,26.2),Dell=c(8.9,9.7,10.8,11.5),Oracle=c(20.5,22.7,21.8,24.4)  )

使用 lapply()，我们找到了每家公司的平均年营收。假设我们想找到所有四家公司按季度平均的营收，我们可以使用 tapply() 如下；这里我们使用函数 c 而不是列表来创建 X：

>X<-c(HP=c(12.5,14.3,16.1,15.4),IBM=c(22,24.5,23.7,26.2),Dell=c(8.9,9.7,10.8,11.5),Oracle=c(20.5,22.7,21.8,24.4)  )

>f<-factor(rep(c("Q1","Q2","Q3","Q4"),times=4) ) 
>f
[1]  Q1 Q2 Q3 Q4 Q1 Q2 Q3 Q4 Q1 Q2 Q3 Q4 Q1 Q2 Q3 Q4
Levels Q1 Q2 Q3 Q4

>tapply(X,f,mean,simplify=TRUE)
Q1           Q2         Q3        Q4
15.97      17.80      18.10     19.37

通过使用以季度值作为级别的因子列表，我们可以使用 tapply() 对每个季度应用 mean() 函数。

数据可视化

R 的一个强大功能是其生成高质量绘图和可视化数据的函数。R 中的绘图函数可以分为三组：

高级绘图函数用于创建新绘图、添加坐标轴、标签和标题。
低级绘图函数用于向现有绘图添加更多信息。这包括添加额外的点、线和标签。
交互式绘图函数用于交互式地向现有绘图添加信息或从现有绘图中提取信息。

R 的基础包本身包含几个绘图函数。对于更高级的图形应用，可以使用如 ggplot2、grid 或 lattice 等包。特别是，ggplot2 对于生成视觉上吸引人的、多层图形非常有用。它基于 图形语法 的概念。由于篇幅限制，我们在这本书中不涵盖这些包。感兴趣的读者应参考本章 参考文献 部分的第 4 个参考文献 Hadley Wickham 的书籍。

高级绘图函数

让我们从 R 中最基本的绘图函数开始，如下所示：

plot( )：这是 R 中最常见的绘图函数。它是一个泛型函数，其输出取决于第一个参数的类型。
plot(x, y)：这会产生 y 与 x 的散点图。
plot(x)：如果 x 是实值向量，输出将是 x 的值与其在 X 轴上的索引的绘图。如果 x 是复数，则将实部与虚部进行绘图。
plot(f, y)：在这里，f 是因子对象，y 是数值向量。该函数为 f 的每个级别生成 y 的箱线图。
plot(y ~ expr)：在这里，y 是任何对象，expr 是由 + 分隔的对象名称列表。该函数将 y 与 expr 中命名的每个对象进行绘图。

R 中有两个用于可视化多元数据的实用函数：

pairs(X): 如果X是一个包含数值数据的 data frame，那么这个函数将生成由X的列定义的变量的成对散点图矩阵。
coplot(y ~ x | z): 如果y和x是数值向量，而z是一个因子对象，那么这个函数将为z的每个水平绘制y与x的关系图。

对于绘制数据的分布，可以使用以下函数：

hist(x): 这条命令生成数值向量x的直方图。
qqplot(x, y): 这条命令将x的量数与y的量数进行比较，以比较它们的相应分布。
qqnorm(x): 这条命令将数值向量x与预期的正态顺序得分进行比较。

低级绘图命令

要向图表添加点和线，可以使用以下命令：

points(x, y): 这条命令将点(x, y)添加到当前图表。
lines(x, y): 这条命令向当前图表添加连接线。
abline(a, b): 这条命令向当前图表添加斜率为b、截距为a的直线。
polygon(x, y, …): 这条命令用于绘制由有序顶点(x, y, …)定义的多边形。

要向图表添加文本，请使用以下函数：

text(x, y, labels): 这条命令在当前图表的点(x, y)处添加文本。
legend(x, y, legend): 这条命令在当前图表的点(x, y)处添加图例。
title(main, sub): 这条命令在当前图表的顶部添加一个标题main，使用大字体，并在底部添加一个副标题sub，使用小字体。
axis(side, …): 这条命令在由第一个参数给出的侧边添加一个轴到当前图表。side可以取从 1 到 4 的值，从底部开始按顺时针方向计数。

以下示例展示了如何绘制散点图并添加趋势线。为此，我们将使用由 R.A. Fisher 创建的著名 Iris 数据集，该数据集在 R 中本身可用：

data(iris)
str(iris)
plot(iris$Petal.Width, iris$Petal.Length, col = "blue", xlab = "X", ylab = "Y")
title(main = "Plot of Iris Data", sub = "Petal Length (Y) Vs Petal Width (X)")
fitlm <- lm(iris$Petal.Length ~ iris$Petal.Width)
abline(fitlm[1], fitlm[2], col = "red")

交互式图形函数

R 中有函数允许用户以交互方式使用鼠标从图表中添加或提取信息：

locator (n , type): 该命令等待用户使用左键在当前图表上选择n个位置。在这里，type可以是n、p、l或o，用于在这些位置绘制点或线。例如，要将图例Outlier放置在异常点附近，可以使用以下代码：
```
>text(locator(1),"Outlier" ,adj=0")
```
identify(x, y, label): 这允许用户通过在附近放置标签来突出显示使用左键选择的任何点x和y。

抽样

通常，我们会对从总体中抽取的代表性数据集感兴趣，用于某些分析或实验设计。这在贝叶斯推断中尤其如此，正如我们将在后面的章节中看到的那样，样本是从后验分布中抽取的。因此，在本章中学习如何从某些已知分布中抽取N个点将是有用的。

在我们使用任何特定的抽样方法之前，读者应该注意，R，像任何其他计算机程序一样，使用伪随机数生成器进行抽样。提供起始种子数以获得可重复的结果是有用的。这可以通过使用 set.seed(n) 命令并使用整数 n 作为种子来完成。

从区间进行随机均匀抽样

要生成在区间 [a, b] 内均匀分布的 n 个随机数（数值），可以使用 runif() 函数：

>runif(5,1,10)  #generates 5 random numbers between 1 and 10
[1]  7.416    9.846    3.093   2.656   1.561

如果没有提供任何参数，runif() 将生成介于 0 和 1 之间的均匀随机数。

如果我们想要生成在区间内均匀分布的随机整数，可以使用 sample() 函数：

>sample(1:100,10,replace=T)   #generates 10 random integers between 1 and 100
[1]  24 51 46 87 30 86 50 45 53 62

选项 replace=T 表示允许重复。

从正态分布中进行抽样

通常，我们可能想要生成符合特定分布的数据，例如正态分布。在单变量分布的情况下，R 有几个内置函数用于此目的。对于从正态分布中抽样数据，要使用的函数是 rnorm()。例如，考虑以下代码：

>rnorm(5,mean=0,sd=1)
[1]  0.759  -1.676   0.569  0.928 -0.609

这将生成五个随机数，这些数遵循均值为 0、标准差为 1 的正态分布。

类似地，可以使用 rbinom() 函数从二项分布中进行抽样，rpois() 用于从泊松分布中进行抽样，rbeta() 用于从 Beta 分布中进行抽样，以及 rgamma() 用于从 Gamma 分布中进行抽样，以提及其他一些分布。

练习

对于本章的以下练习，我们使用来自 UCI 机器学习存储库的 Auto MPG 数据集（本章参考文献部分的第 5 和 6 条）。数据集可以从 archive.ics.uci.edu/ml/datasets.html 下载。该数据集包含 1970-1982 年间在美国测量的汽车油耗。除了消耗值外，还有属性变量，如气缸数、排量、马力、重量、加速度、年份、产地和汽车名称：

使用 read.table() 函数将数据集加载到 R 中。
为每个汽车名称生成 mpg 值的箱线图。
编写一个函数，该函数将计算给定列名的缩放值（减去平均值并除以标准差）。
使用 lapply() 函数计算所有变量的缩放值。
使用 coplot() 函数为每个汽车名称生成 mpg 与加速度的散点图。使用图例注释图形。

参考文献

Matloff N. R 编程艺术 – 统计软件设计之旅. No Starch Press. 2011. ISBN-10: 1593273843
Teetor P. R 烹饪书. O'Reilly Media. 2011. ISBN-10: 0596809158
Wickham H. 高级 R. Chapman & Hall/CRC The R Series. 2015. ISBN-10: 1466586966
Wickham H. ggplot2: 数据分析的优雅图形（使用 R!）. Springer. 2010. ISBN-10: 0387981403
Auto MPG 数据集，UCI 机器学习仓库，archive.ics.uci.edu/ml/datasets/Auto+MPG
Quinlan R. "Combining Instance-Based and Model-Based Learning". In: Tenth International Conference of Machine Learning. 236-243. University of Massachusetts, Amherst. Morgan Kaufmann. 1993

小贴士

下载示例代码

您可以从www.packtpub.com下载您购买的所有 Packt Publishing 书籍的示例代码文件。如果您在其他地方购买了这本书，您可以访问www.packtpub.com/support并注册，以便将文件直接通过电子邮件发送给您。

摘要

在本章中，您被介绍了 R 环境。阅读完本章后，您学习了如何将数据导入 R，选择数据子集进行分析，以及使用函数和控制结构编写简单的 R 程序。此外，您现在应该熟悉 R 的图形功能以及一些高级功能，例如循环函数。在下一章中，我们将开始介绍本书的核心主题，贝叶斯推理。

第三章。介绍贝叶斯推理

在第一章中，介绍概率论，我们学习了贝叶斯定理作为两个随机变量，如 A 和 B 的条件概率之间的关系。这个定理是贝叶斯推理中根据观察更新信念或模型参数值的基础。在本章中，将给出贝叶斯推理的更正式处理。首先，让我们尝试理解在贝叶斯方法中如何处理现实世界问题中的不确定性。

贝叶斯不确定性视图

经典的或频率派统计学通常认为，任何产生含噪声数据的物理过程都可以用具有固定参数值的随机模型来建模。参数值通过如最大似然估计等程序从观察数据中学习。基本思想是在参数空间中搜索以找到最大化观察到的数据概率的参数值。既没有以正式方式处理从数据中估计模型参数的不确定性，也没有处理解释研究现象的模型本身的不确定性。另一方面，贝叶斯方法使用概率来处理所有不确定性的来源。因此，解释观察数据集的模型及其参数都不是固定的，而是被视为不确定变量。贝叶斯推理提供了一个框架来学习模型参数的整个分布，而不仅仅是最大化观察数据概率的值。学习可以来自观察数据提供的证据和专家的领域知识。还有一个框架来选择最适合解释给定数据集的模型家族中的最佳模型。

一旦我们有了模型参数的分布，我们就可以通过联合概率分布的边缘化来消除使用学习模型预测的随机变量的未来值中参数估计不确定性的影响。这正如在第一章中解释的，介绍概率论。

再次考虑 N 个随机变量的联合概率分布，正如在第一章中讨论的，介绍概率论：

这次，我们在概率分布的参数中添加了一个额外的项 m，以明确表示参数是由模型 m 生成的。然后，根据贝叶斯定理，给定观察数据和模型 m 的条件下，模型参数的概率分布如下：

正式来说，方程的左侧的项被称为后验概率分布。右侧分子中出现的第二个项被称为先验概率分布。它代表了在观察任何数据之前，例如从领域知识中得到的对模型参数的先验信念。先验分布也可以有参数，这些参数被称为超参数。项是模型 m 解释观察数据的似然。由于，它可以被视为一个归一化常数。前面的方程可以重写为以下迭代形式：

这里，代表在时间步 n 获得的观测值，是直到时间步 n - 1 更新的边缘参数分布，而是在时间步 n 观察到观测值后更新的模型参数分布。

将贝叶斯定理以这种迭代形式表述对于在线学习是有用的，并且它暗示以下内容：

模型参数可以在获得更多数据和证据的迭代过程中学习。
使用迄今为止看到的数据估计的后验分布，在获得下一组观测值时可以被视为先验模型。
即使没有数据可用，也可以仅基于领域知识创建的先验分布进行预测。

为了使这些观点更加清晰，让我们通过一个简单的说明性例子来说明。考虑这样一个案例，即一个人试图估计某个地区男性身高的分布。用于此例的数据是从人口中随机抽取的 M 个志愿者的身高测量值（以厘米为单位）。我们假设身高服从均值为和方差的正态分布：

如前所述，在经典统计学中，人们试图从观测数据中估计和的值。除了每个参数的最佳估计值外，还可以确定估计的误差项。另一方面，在贝叶斯方法中，和也被视为随机变量。为了简化，我们假设是一个已知的常数。此外，我们假设的先验分布是一个具有(超)参数和的正态分布。在这种情况下，的后验分布的表达式如下：

为了方便起见，我们在这里使用符号来表示。将乘积中的项展开并完成指数中的平方是一个简单的练习。这作为本章末尾的练习给出。后验分布的结果表达式如下：

在这里，代表样本均值。尽管前面的表达式看起来很复杂，但它有一个非常简单的解释。后验分布也是一个正态分布，具有以下均值：

方差如下：

后验均值是先验均值 和样本均值 *的加权平均值。随着样本大小 M 的增加，样本均值的权重增加，而先验的权重减少。同样，后验精度（方差的倒数）是先验精度和样本均值精度的和：

随着 M 的增加，来自观测（证据）的精度贡献超过了来自先验知识的贡献。

让我们举一个具体的例子，其中我们考虑具有人口均值 5.5 和人口标准差 0.5 的年龄分布。我们使用以下 R 脚本来从这个总体中抽取 100 人：

>set.seed(100)
>age_samples <- rnorm(10000,mean = 5.5,sd=0.5)

我们可以使用以下 R 函数来计算后验分布：

>age_mean <- function(n){
  mu0 <- 5
  sd0 <- 1
  mus <- mean(age_samples[1:n])
  sds <- sd(age_samples[1:n])
  mu_n <- (sd0²/(sd0² + sds²/n)) * mus + (sds²/n/(sd0² + sds²/n)) * mu0
  mu_n
}
>samp <- c(25,50,100,200,400,500,1000,2000,5000,10000)
>mu <- sapply(samp,age_mean,simplify = "array")
>plot(samp,mu,type="b",col="blue",ylim=c(5.3,5.7),xlab="no of samples",ylab="estimate of mean")
>abline(5.5,0)

可以看到，随着样本数量的增加，估计的均值渐近地接近总体均值。初始的低值是由于先验的影响，在这种情况下，先验是 5.0。

这个关于先验知识和观察证据如何贡献于整体模型参数估计的简单直观图景在任何贝叶斯推理中都成立。它们如何结合的精确数学表达式会有所不同。因此，一个人可以仅使用先验信息，无论是来自领域知识还是过去收集的数据，就开始使用模型进行预测。此外，随着新观察数据的到来，模型可以使用贝叶斯方案进行更新。

选择合适的先验分布

在前面的简单例子中，我们看到了如果似然函数具有正态分布的形式，并且当先验分布被选为正态分布时，后验分布也变成了正态分布。此外，我们还可以得到后验均值的封闭形式解析表达式。由于后验是通过将先验和似然函数相乘，并通过在参数变量上的积分进行归一化得到的，因此先验分布的形式对后验有显著影响。本节将给出有关不同类型先验分布的更多细节，以及在某些特定情境下如何选择先验分布的指导方针。

有不同的方式以正式的方式对先验分布进行分类。其中一种方法是基于先验提供的信息量。在这个方案中，先验分布被分类为 信息性先验、弱信息性先验、最少信息性先验 和 非信息性先验。对每个这些类别的详细讨论超出了本书的范围，感兴趣的读者应参考相关书籍（本章“参考文献”部分的参考 1 和 2）。在这里，我们采取更多实践者的方法，并说明了在实践中最常用的先验分布的一些重要类别。

非信息先验

让我们从我们没有关于模型参数的任何先验知识的情况开始。在这种情况下，我们希望通过数学表达式表达对模型参数的完全无知。这是通过所谓的非信息先验实现的。例如，在单个随机变量 x 可以取和之间任何值的情形下，其均值的非信息先验如下：

在这里，参数值的完全无知通过参数空间中的均匀分布函数来捕捉。请注意，均匀分布不是一个合适的分布函数，因为其在域上的积分不等于 1；因此，它不是可归一化的。然而，只要乘以似然函数，就可以使用不适当的先验分布函数；结果的后验可以归一化。

如果感兴趣的参数是方差，那么根据定义，它只能取非负值。在这种情况下，我们变换变量，使得变换后的变量在到的范围内具有均匀概率：

使用简单的微分微积分很容易证明，在原始变量中对应的非信息分布函数如下：

在实际应用中另一个著名的非信息先验是杰弗里斯先验，它是以英国统计学家哈罗德·杰弗里斯的名字命名的。这个先验在的重参数化下是不变的，并且定义为与费舍尔信息矩阵的行列式的平方根成比例：

在这里，讨论一下费舍尔信息矩阵是值得的。如果X是一个按照分布的随机变量，我们可能想知道观察X携带了多少关于未知参数的信息。这正是费舍尔信息矩阵提供的内容。它定义为得分（似然函数对数的一阶导数）的第二矩：

让我们通过一个简单的二维问题来理解费舍尔信息矩阵和杰弗里斯先验。这个例子是由加州大学的 D. Wittman 教授给出的（本章“参考文献”部分的第 3 个参考文献）。让我们考虑两种食品：面包和热狗。

假设它们通常成对产生（热狗和面包对），但偶尔热狗也会在单独的过程中独立生产。有两个可观测量，如热狗的数量 () 和面包的数量 ()，以及两个模型参数，如成对的产量 () 和单独热狗的产量 ()。我们假设这两种食品产品数量的测量不确定性分别按照正态分布分布，方差分别为和。在这种情况下，该问题的费舍尔信息矩阵如下：

在这种情况下，费舍尔信息矩阵的逆对应于协方差矩阵：

我们在本章的练习部分包含了一个问题，用于计算费舍尔信息矩阵和杰弗里先验。请读者尝试解决这个问题，以了解如何从观察结果中计算杰弗里先验。

主观先验

与经典（频率主义）统计学相比，贝叶斯统计学的关键优势之一是框架允许人们捕捉关于任何随机变量的主观信念。通常，人们会对随机变量的最小值、最大值、平均值以及最可能或峰值有直观的感觉。例如，如果一个人对热带国家冬季每小时温度的分布感兴趣，那么熟悉热带气候或气候学专家将相信，在冬季，温度可以低至 15°C，高至 27°C，最可能的温度值是 23°C。这可以通过三角分布作为先验分布来捕捉，如图所示。

三角分布有三个参数，分别对应最小值 (a)、最可能值 (b) 和最大值 (c)。该分布的均值和方差如下给出：

一个人还可以使用 PERT 分布来表示对随机变量的最小值、最大值和最可能值的信念。PERT 分布是一个重新参数化的 Beta 分布，如下所示：

在这里：

PERT 分布常用于项目完成时间分析，其名称来源于项目评估和审查技术。另一个经常使用三角形和 PERT 分布的领域是风险管理。

通常，人们也会对随机变量值的相对概率有所信念。例如，当研究日本或某些欧洲国家等人口中老年人比年轻人多的年龄分布时，专家可以为不同年龄在人口中的概率给出相对权重。这可以通过包含以下详细信息的相对分布来捕捉：

在这里，min 和 max 分别代表最小值和最大值，{values} 代表可能观察到的值的集合，而 {weights} 代表它们的相对权重。例如，在人口年龄分布问题中，这些可能如下所示：

权重不需要总和为 1。

共轭先验

如果先验分布和后验分布属于同一分布族，那么它们被称为共轭分布，相应的先验被称为似然函数的共轭先验。共轭先验对于获取后验分布的解析封闭形式表达式非常有帮助。在我们考虑的简单例子中，我们看到了当噪声按照正态分布分布时，选择均值正态先验会导致后验正态分布。以下表格给出了我们将在这本书的后续章节中使用的某些著名共轭对的例子：

似然函数	模型参数	共轭先验
二项分布	(概率)	Beta 分布
泊松分布	(速率)	伽马分布
分类变量	(概率，类别数量)	Dirichlet 分布
单变量正态分布（已知方差）	(均值)	正态分布
单变量正态分布（已知均值）	(方差)	逆伽马分布

层次先验

有时，为超参数本身定义先验分布是有用的。这与贝叶斯观点一致，即所有参数都应通过使用概率被视为不确定。这些分布被称为超先验分布。在理论上，人们可以将此扩展到许多层次，作为一个层次模型。这是获取最优先验分布的一种方式。例如：

是具有超参数的先验分布。我们可以通过第二组方程定义的先验分布，如下所示：

在这里，是超参数的超先验分布，由超超参数参数化。人们可以用相同的方式为定义一个先验分布，并永远继续这个过程。将此类模型形式化的实际原因是，在某个层次级别上，人们可以为超参数定义一个均匀先验，反映对参数分布的完全无知，并有效地截断层次结构。在实际情况下，通常在第二级进行此操作。这对应于前一个示例中使用的均匀分布。

我想通过强调一个重要观点来结束本节。尽管先验分布在贝叶斯推理中起着重要作用，但只要选择的先验是合理的，并与迄今为止看到的领域知识和证据一致，就无需过分担心。原因如下：首先，随着我们获得更多证据，先验的重要性会减弱。其次，当我们使用贝叶斯模型进行预测时，我们将使用后验分布对参数估计的不确定性进行平均。这种平均是贝叶斯推理的关键成分，它消除了选择正确先验时许多的模糊性。

后验分布的估计

到目前为止，我们讨论了贝叶斯推理背后的基本概念以及如何选择先验分布。由于在使用模型进行预测之前，需要计算模型参数的后验分布，因此我们在这节中讨论这个任务。尽管贝叶斯规则看起来非常简单，但在实际可用的方式下计算后验分布通常非常具有挑战性。这主要是因为当有 N 个参数时，计算归一化常数涉及到 N-维积分。即使使用共轭先验，这种计算也可能非常难以分析或数值跟踪。这就是为什么直到最近几十年，人们不使用贝叶斯推理进行多元建模的主要原因之一。在本节中，我们将探讨在实践中使用的各种计算后验分布的近似方法。

最大后验概率估计

最大后验概率（MAP）估计是一种点估计，对应于取后验分布的最大值或众数。尽管取点估计不能捕捉参数估计中的变异性，但与最大似然估计相比，它在一定程度上考虑了先验分布的影响。MAP 估计也被称为穷人的贝叶斯推理。

从贝叶斯规则，我们有：

这里，为了方便起见，我们使用了 X 表示 N-维向量。最后一个关系成立是因为贝叶斯规则右侧分母与无关。将此与以下最大似然估计进行比较：

MAP 估计和 ML 估计之间的区别在于，ML 找到似然函数的众数，而 MAP 找到似然函数与先验的乘积的众数。

拉普拉斯近似

我们看到，最大后验概率估计（MAP）只是找到后验分布的最大值。拉普拉斯近似更进一步，并计算最大值周围的局部曲率，直到二次项。这相当于假设后验分布在大数据量相对于参数数量时近似为高斯（正态）分布：M >> N。

这里，A 是通过取后验分布对数的导数得到的 N x N 赫 essian 矩阵：

使用以下条件概率的定义，可以简单地评估前面的表达式在：

我们可以从拉普拉斯近似中得到一个关于 P(X|m) 的表达式，其形式如下：

在大量样本的极限情况下，可以证明这个表达式可以简化为以下形式：

术语被称为贝叶斯信息准则（BIC），可用于模型选择或模型比较。这是统计模型的一个拟合优度术语。另一个常用的类似标准是赤池信息准则（AIC），其定义为。

现在我们将讨论如何使用 BIC 来比较不同的模型以进行模型选择。在贝叶斯框架中，使用贝叶斯因子比较两个模型，如和。贝叶斯因子的定义是后验优势与先验优势之比，由以下给出：

在这里，后验优势是给定数据的两个模型的后验概率之比，先验优势是两个模型的先验概率之比，如前一个方程所示。如果，数据更倾向于模型；如果，数据更倾向于模型。

在现实中，计算贝叶斯因子是困难的，因为很难得到精确的先验概率。可以证明，在大的 N 极限下，可以被视为的粗略近似。

蒙特卡洛模拟

我们迄今为止讨论的两个近似，即最大似然估计（MAP）和拉普拉斯近似，当后验概率在最大值附近有一个非常尖锐的峰值时是有用的。在现实生活中的许多情况下，后验概率通常会有长尾。例如，在电子商务中，用户购买产品的概率在所有产品空间中具有长尾。因此，在许多实际情况下，MAP 和拉普拉斯近似都无法给出良好的结果。另一种方法是直接从后验分布中进行采样。蒙特卡洛模拟是一种用于从后验分布中采样的技术，是实际应用中贝叶斯推理的得力工具。在本节中，我们将向读者介绍马尔可夫链蒙特卡洛（MCMC）模拟，并讨论实践中常用的两种 MCMC 方法。

如前所述，让是我们通过后验分布从数据中感兴趣估计的参数集。考虑参数是离散的情况，其中每个参数有K个可能的值，即。设置一个具有状态和转移概率矩阵的马尔可夫过程。MCMC 模拟背后的基本思想是，可以选择转移概率，使得马尔可夫链的稳态分布将对应于我们感兴趣的后续分布。一旦这样做，在马尔可夫链达到稳态后，从马尔可夫链输出中采样将给出的样本，这些样本遵循后验分布。

现在，问题是如何设置马尔可夫过程，使其稳态分布对应于感兴趣的后续分布。为此有两种众所周知的方法。一种是 Metropolis-Hastings 算法，另一种是 Gibbs 抽样。我们将在下面详细讨论这两种方法。

Metropolis-Hasting 算法

Metropolis-Hasting 算法是第一个被提出的用于 MCMC（本章“参考文献”部分的第 4 个参考文献）的主要算法之一。它有一个非常简单的概念——类似于优化中的爬山算法：

让表示时间步长t时系统的状态。
要在时间步长t + 1时将系统移动到另一个状态，通过从提议分布中采样生成一个候选状态。提议分布被选择得易于从中采样。
以以下概率接受提议的移动：
如果被接受， = ；如果不接受，。
继续这个过程，直到分布收敛到稳态。

在这里，是我们想要模拟的后验分布。在特定条件下，前面的更新规则将保证，在长时间极限下，马尔可夫过程将趋近于的稳态分布。

Metropolis-Hasting 算法背后的直觉很简单。提议分布给出了从当前状态在下一个时间步过渡到状态的条件概率。因此，是系统当前处于状态并且会在下一个时间步过渡到状态的概率。同样，是系统当前处于状态并且会在下一个时间步过渡到状态的概率。如果这两个概率的比率大于 1，则接受移动。或者，仅以比率给出的概率接受移动。因此，Metropolis-Hasting 算法类似于一种爬山算法，其中接受所有向上的移动，并且偶尔以较小的概率接受向下的移动。向下的移动有助于系统不陷入局部最小值。

让我们重新回顾一下在引言部分讨论的估计人群中身高均值和方差后验分布的例子。这次我们将使用 Metropolis-Hasting 算法来估计后验分布。以下 R 代码行完成这项工作：

>set.seed(100)
>mu_t <- 5.5
>sd_t <- 0.5
>age_samples <- rnorm(10000,mean = mu_t,sd = sd_t)

>#function to compute log likelihood
>loglikelihood <- function(x,mu,sigma){
    singlell <- dnorm(x,mean = mu,sd = sigma,log = T)
    sumll <- sum(singlell)
    sumll
    }

>#function to compute prior distribution for mean on log scale
>d_prior_mu <- function(mu){
  dnorm(mu,0,10,log=T)
  }

>#function to compute prior distribution for std dev on log scale
>d_prior_sigma <- function(sigma){
  dunif(sigma,0,5,log=T)
  }

>#function to compute posterior distribution on log scale
>d_posterior <- function(x,mu,sigma){
  loglikelihood(x,mu,sigma) + d_prior_mu(mu) + d_prior_sigma(sigma)
   }

>#function to make transition moves
  tran_move <- function(x,dist = .1){
  x + rnorm(1,0,dist)
  }

>num_iter <- 10000
>posterior <- array(dim = c(2,num_iter))
>accepted <- array(dim=num_iter - 1)
>theta_posterior <-array(dim=c(2,num_iter))

>values_initial <- list(mu = runif(1,4,8),sigma = runif(1,1,5))
>theta_posterior[1,1] <- values_initial$mu
>theta_posterior[2,1] <- values_initial$sigma

>for (t in 2:num_iter){
   #proposed next values for parameters
    theta_proposed <- c(tran_move(theta_posterior[1,t-1]),tran_move(theta_posterior[2,t-1]))
    p_proposed <- d_posterior(age_samples,mu = theta_proposed[1],sigma = theta_proposed[2])
    p_prev <-d_posterior(age_samples,mu = theta_posterior[1,t-1],sigma = theta_posterior[2,t-1])
    eps <- exp(p_proposed - p_prev)

    # proposal is accepted if posterior density is higher w/ theta_proposed
    # if posterior density is not higher, it is accepted with probability eps
    accept <- rbinom(1,1,prob = min(eps,1))
    accepted[t - 1] <- accept
    if (accept == 1){
      theta_posterior[,t] <- theta_proposed
    } else {
      theta_posterior[,t] <- theta_posterior[,t-1]
    }
}

要绘制得到的后验分布，我们使用 R 中的 sm 包：

>library(sm)
x <- cbind(c(theta_posterior[1,1:num_iter]),c(theta_posterior[2,1:num_iter]))
xlim <- c(min(x[,1]),max(x[,1]))
ylim <- c(min(x[,2]),max(x[,2]))
zlim <- c(0,max(1))

sm.density(x,
           xlab = "mu",ylab="sigma",
           zlab = " ",zlim = zlim,
           xlim = xlim ,ylim = ylim,col="white")
title("Posterior density")

得到的后验分布将类似于以下图示：

虽然 Metropolis-Hasting 算法对于任何贝叶斯推断问题来说实现起来都很简单，但在实践中，它在许多情况下可能并不非常高效。主要原因在于，除非一个人仔细选择一个提议分布，否则会有太多的拒绝，并且需要大量的更新才能达到稳态。当参数数量较高时，这种情况尤为明显。有各种基本 Metropolis-Hasting 算法的修改版本，试图克服这些困难。在下一节讨论各种 R 包的 Metropolis-Hasting 算法时，我们将简要介绍这些修改。

Metropolis-Hasting 算法的 R 包

R 中存在多个用于 Metropolis-Hasting 算法的 MCMC 模拟的贡献包，在此我们将介绍一些流行的包。

由 Charles J. Geyer 和 Leif T. Johnson 贡献的 mcmc 包是 R 中用于 MCMC 模拟的流行包之一。它包含一个名为 metrop 的函数，用于运行基本的 Metropolis-Hasting 算法。metrop 函数使用多元正态分布作为提议分布。

有时，进行变量变换以改善 MCMC 收敛速度是有用的。mcmc 包有一个名为 morph 的函数用于此目的。结合这两个函数，morph.metrop 函数首先对变量进行变换，然后在变换后的密度上进行 Metropolis 采样，并将结果转换回原始变量。

除了 mcmc 包之外，R 中还有两个有用的包：由 Corey Chivers 贡献的 MHadaptive 包和 Gopi Goswami 的 进化蒙特卡洛（EMC）算法包。由于篇幅限制，我们不会在本书中讨论这两个包。有兴趣的读者请从 C-RAN 项目网站下载这些包并尝试使用它们。

吉布斯采样

如前所述，Metropolis-Hasting 算法由于过多的拒绝而存在收敛性差的缺点，如果不选择一个好的提议分布。为了避免这个问题，两位物理学家 Stuart Geman 和 Donald Geman 提出了一种新的算法（本章“参考文献”部分的第 5 个参考文献）。这个算法被称为吉布斯采样，并以著名的物理学家 J W Gibbs 命名。目前，吉布斯采样是贝叶斯推理中 MCMC 的主要工具。

设为我们希望估计的模型参数集合：

从一个初始状态开始。
在每个时间步，逐个更新组件，通过从基于其他组件最近值的分布中进行抽样：
经过 N 步后，参数的所有组件都将被更新。
继续进行步骤 2，直到马尔可夫过程收敛到稳态。

吉布斯采样是一个非常有效的算法，因为没有拒绝。然而，要使用吉布斯采样，必须知道后验分布的条件分布形式。

吉布斯采样的 R 包

不幸的是，在 R 中，贡献的通用吉布斯采样包并不多。gibbs.met 包提供了两个通用函数，用于以朴素的方式为用户定义的目标分布执行 MCMC。第一个函数是 gibbs_met。该函数使用 Metropolis 算法对每个一维分布进行采样，提议分布为正态分布。第二个函数 met_gaussian 使用独立正态分布更新整个状态，其中心位于前一个状态。gibbs.met 包对于通用 MCMC 在中等维度问题上是很有用的。

在本章的练习部分，我们将讨论一个涉及通过使用 Metropolis-Hasting 算法和 Gibbs 抽样从二维正态分布中采样的问题，以使这些概念更加清晰。读者可以使用这些提到的包来解决这个练习。

除了通用的 MCMC 包之外，R 中还有几个包专门用于解决特定类型的机器学习问题。GibbsACOV包可用于单因素混合效应方差分析（ANOVA）和方差分析协方差（ANCOVA）模型。lda包执行主题（LDA）模型的折叠 Gibbs 抽样方法。stocc包通过 Gibbs 抽样拟合空间占用模型。binomlogit包实现了二项 Logit 模型的效率 MCMC。Bmk是一个用于进行 MCMC 输出的诊断的包。Bayesian Output Analysis Program（BOA）是另一个类似的包。RBugs是知名OpenBUGS MCMC 包的接口。ggmcmc包是分析 MCMC 模拟的图形工具。MCMCglm是一个用于广义线性混合模型的包，BoomSpikeSlab是一个用于进行 Spike 和 Slab 回归的 MCMC 的包。最后，SamplerCompare是一个用于比较各种 MCMC 包性能的包（更多是一个框架）。

变分近似

在变分近似方案中，假设后验分布可以被近似为因式分解形式：

注意，因式分解形式也是一种条件分布，因此每个都可以通过条件变量X依赖于其他。换句话说，这不是一个简单的因式分解，使得每个参数相互独立。这种因式分解的优势在于可以选择更易于分析的分布函数形式。实际上，可以通过调整函数，使其尽可能接近真实的后验。这在数学上被表述为一个变分法问题，如这里所述。

让我们使用一些度量来计算两个概率分布之间的距离，例如和，其中。概率分布之间距离的一个标准度量是 Kullback-Leibler 散度，或简称 KL 散度。它被定义为以下：

它被称为散度而不是距离的原因是在Q和P上不对称。可以使用关系并将前面的表达式重写为log P(X)的方程：

这里：

注意，在ln P(X)的方程中，左侧对Q没有依赖。因此，相对于Q最大化将最小化，因为它们的和是一个与Q无关的项。通过选择对Q进行解析处理的函数，可以在实践中进行这种最大化。这将导致后验分布的一个近似以及ln P(X)的下界，即证据的对数或边缘似然，因为。

因此，变分近似一次给出了两个量。后验分布可以用来对未来观测进行预测（如下一节所述），而证据的下界可以用于模型选择。

如何在实际中实现 KL 散度的最小化？不深入数学细节，我们在这里给出解的最终表达式：

在这里，意味着对联合分布的对数期望是在所有参数上取的，除了。因此，KL 散度的最小化导致了一组耦合方程；每个都需要自洽地求解以获得最终解。尽管变分近似在数学上看起来非常复杂，但它有一个非常简单、直观的解释。每个参数的后验分布是通过对所有其他变量的联合分布的对数进行平均得到的。这在物理学中类似于平均场理论，其中，如果有N个相互作用的带电粒子，系统可以通过说每个粒子处于一个由所有其他粒子产生的场的平均值所构成的恒定外部场来近似。

我们将通过提及一些用于变分近似的 R 包来结束本节。VBmix包可用于贝叶斯混合模型中的变分近似。一个类似的包是vbdm，用于贝叶斯离散混合模型。vbsr包用于 Spike 回归正则化线性模型中的变分推理。

未来观测的预测

一旦我们使用已经描述的一些方法从数据中推断出后验分布，就可以用它来预测未来的观测值。给定观察到的数据X和参数的后验分布，观察到值Y的概率由以下公式给出：

注意，在这个表达式中，似然函数是通过后验给出的参数分布进行平均的。这实际上是贝叶斯推理的核心优势。这种贝叶斯平均消除了估计参数值的不确定性，并使预测更加稳健。

练习

通过展开指数中的平方，对每个i收集所有类似的幂次项，并再次使其成为完全平方，推导出后验均值的方程。注意，指数的乘积可以写成求和项的指数。
对于这个练习，我们使用的是来自 UCI 机器学习仓库（archive.ics.uci.edu/ml/datasets/Smartphone-Based+Recognition+of+Human+Activities+and+Postural+Transitions）对应的智能手机基于人类活动和姿势转换识别的数据集。它包含从智能手机加速度计获取的加速度值。原始数据集包含加速度的x、y和z分量以及相应的时间戳值。在这个练习中，我们只使用了加速度的两个水平分量x和y。在这个练习中，我们假设加速度遵循正态分布。我们还假设加速度均值具有正态先验分布，均值的超参数在区间(-0.5, 0.5)内均匀分布，且已知方差等于 1。通过方程给出的表达式找到后验均值。
编写一个 R 函数来计算 Fisher 信息矩阵。使用本节练习 1 中提到的数据集，获取这个问题的 Fisher 信息矩阵。
使用 R 中的mcmc包设置这个问题的 MCMC 模拟。绘制模拟数据的直方图。
使用 Gibbs 抽样设置 MCMC 模拟。将结果与 Metropolis 算法的结果进行比较。

参考文献

Berger J.O. 《统计决策理论及贝叶斯分析》。Springer 统计学系列。1993 年。ISBN-10: 0387960988
Jayes E.T. 《概率论：科学的逻辑》。剑桥大学出版社。2003 年。ISBN-10: 052159271
Wittman D. "Fisher 矩阵入门". 加州大学戴维斯分校物理系 (www.physics.ucdavis.edu/~dwittman/Fisher-matrix-guide.pdf)
Metropolis N, Rosenbluth A.W., Rosenbluth M.N., Teller A.H., Teller E. "快速计算机器的物态方程计算". 化学物理杂志 21 (6): 1087–1092. 1953
Geman S., Geman D. "随机松弛、吉布斯分布以及图像的贝叶斯恢复". IEEE 交易杂志：模式分析与机器智能 6 (6): 721-741. 1984

摘要

在本章中，我们介绍了贝叶斯推理的基本原理。从贝叶斯统计与经典统计在处理不确定性方面的不同之处开始，我们深入讨论了贝叶斯定理的各个组成部分。首先，我们学习了不同类型的先验分布以及如何为你的问题选择正确的分布。然后，我们学习了使用如 MAP 估计、拉普拉斯近似和 MCMC 模拟等技术来估计后验分布。一旦读者理解了本章内容，他们就能将贝叶斯原理应用于他们的数据分析问题。在我们开始讨论具体的贝叶斯机器学习问题之前，下一章我们将回顾机器学习的一般知识。

第四章：使用贝叶斯推理的机器学习

既然我们已经了解了贝叶斯推理和 R，现在是时候将两者都用于机器学习了。在本章中，我们将概述不同的机器学习技术，并在随后的章节中详细讨论每一个。机器学习是计算机科学和统计学交叉的领域，也是人工智能或 AI 的一个子分支。这个名字本质上来源于 AI 早期的工作，当时研究人员试图开发学习机器，这些机器能够仅从数据中自动学习输入和输出变量之间的关系。一旦机器在一个给定问题的数据集上进行了训练，它就可以作为一个黑盒来预测新输入变量的输出变量的值。

将机器的这个学习过程设置在数学框架中是有用的。设X和Y是两个随机变量，我们寻求一个学习机器，它能从数据中学习这两个变量之间的关系，并预测给定X值的Y值。系统完全由联合概率分布P(X,Y)来表征，然而，这个分布的形式是未知的。学习的目标是找到一个函数f(X)，它将X映射到Y，使得预测包含尽可能小的误差。为了实现这一点，选择一个损失函数L(Y, f(X))，并找到一个f(X)，它最小化给定X和Y的联合分布的期望或平均损失，如下所示：

在统计决策理论中，这被称为经验风险最小化。典型的损失函数是平方损失函数(使用贝叶斯推理的机器学习)，如果Y是一个连续变量，以及Hinge 损失函数(使用贝叶斯推理的机器学习)，如果Y是一个具有值的二进制离散变量。第一种情况通常被称为回归，第二种情况称为二元分类，正如我们将在本章后面看到的。

这里描述的数学框架被称为监督学习，其中机器被呈现一个包含对应于对(Y, X)对的地面真实值的训练数据集。让我们再次考虑平方损失函数的情况。在这里，学习任务是要找到一个f(X)，它最小化以下内容：

由于目标是预测给定X值的Y值，我们在积分内使用了P(X, Y)的分解来使用条件分布P(Y|X)。可以证明，先前损失函数的最小化导致以下解：

前述方程的意义是，对于任何输入值X，Y的最佳预测是条件概率分布P(Y|X)在X处的均值或期望，用E表示。

在第三章《介绍贝叶斯推理》中，我们提到了最大似然估计（MLE）作为学习任何分布P(X)参数的方法。一般来说，如果基础分布是正态分布，MLE 与平方损失函数的最小化相同。

注意，在经验风险最小化中，我们正在学习参数，E[(Y|X)]，给定X值的条件分布的均值。我们将使用一个特定的机器学习任务，线性回归，来解释贝叶斯推理相对于经典学习方法的优势。然而，在这之前，我们将简要解释一些机器学习的更一般方面。

监督机器学习模型有两种类型，即生成模型和判别模型。在生成模型的情况下，算法试图从数据中学习X和Y的联合概率，即P(X,Y)，并使用它来估计均值P(Y|X)。在判别模型的情况下，算法试图直接学习期望的函数，即P(Y|X)的均值，并且不尝试对X变量进行建模。

在训练数据中对目标变量的值进行标记是手动完成的。这使得当需要使用非常大的数据集，如文本分析中的情况时，监督学习变得非常昂贵。然而，非常常见的是，监督学习方法产生最准确的结果。

如果没有足够的训练数据可供学习，一个人仍然可以通过无监督学习来使用机器学习。在这里，学习主要是通过发现数据集中变量之间关联模式的过程。聚类具有相似特征的点是一个经典的例子。

强化学习是机器学习的第三种类型，学习发生在动态环境中，机器需要根据其当前状态执行某些动作。与每个动作相关联的是奖励。机器需要学习在每个状态下需要采取什么动作，以便最大化总奖励。这通常是机器人学习在现实生活环境中执行任务，如驾驶车辆的方式。

为什么机器学习要使用贝叶斯推理？

我们已经在上一章讨论了贝叶斯统计相对于经典统计的优势。在本章中，我们将更详细地探讨我们上一章学到的贝叶斯推理的一些概念在机器学习中的应用。为此，我们选取一个简单的机器学习任务，即线性回归。让我们考虑一个学习任务，其中我们有一个包含 N 对点的数据集 D ，目标是使用线性回归构建一个机器学习模型，该模型可以用来预测的值，给定新的值。

在线性回归中，首先，我们假设 Y 的形式如下：

在这里，F(X) 是一个函数，它捕捉了 X 和 Y 之间真实关系的函数，是一个误差项，它捕捉了数据中固有的噪声。假设这种噪声由均值为 0、方差为的正态分布来表征。这意味着，如果我们有一个无限大的训练数据集，我们可以从数据中学习到 F(X) 的形式，即使在这种情况下，我们也只能预测到加上一个噪声项的 Y。在实践中，我们只有有限大小的训练数据集 D；因此，我们只能学习到 F(X) 的一个近似形式，用表示。

注意，我们在这里讨论了两种类型的误差。一种是误差项，这是由于数据中固有的噪声造成的，我们对此无能为力。第二种误差是在学习 F(X) 时，通过函数从数据集 D 中近似地学习。

通常，，它是输入变量X和输出变量Y之间近似映射的函数，是一个包含一组参数的X的函数。当是参数的线性函数时，我们说学习模型是线性的。一个普遍的误解是，只有当是X的线性函数时，线性回归才对应这种情况。参数线性而不是X的原因是，在最小化损失函数的过程中，实际上是通过最小化参数值来寻找最佳的。因此，在上是线性的函数会导致一个可以更容易地解析和数值解决的线性优化问题。因此，线性回归对应以下：

这是一个在M个基函数上的扩展。在这里，每个基函数是X的函数，没有任何未知参数。在机器学习中，这些被称为特征函数或模型特征。对于线性回归问题，因此损失函数可以写成以下形式：

这里，是参数向量的转置，B(X)是由基函数组成的向量。从数据集中学习意味着通过一些优化方案（如梯度下降）最小化损失函数来估计的值。

选择尽可能多的基函数来捕捉数据中的有趣模式是很重要的。然而，选择更多的基函数或特征会导致模型过拟合，即它甚至开始拟合数据中包含的噪声。过拟合会导致对新输入数据的预测效果变差。因此，选择最佳特征的最佳数量以最大化任何机器学习模型的预测准确性是很重要的。在基于经典统计学的机器学习中，这通过所谓的偏差-方差权衡和模型正则化来实现。而在基于贝叶斯推理的机器学习中，可以通过贝叶斯模型平均来最大化预测模型的准确性，无需施加模型正则化或偏差-方差权衡。我们将在以下各节中学习这些概念。

模型过拟合与偏差-方差权衡

在上一节中提到的预期损失，在采用平方损失函数的线性回归情况下，可以写成三个项的和，如下所示：

在这里，偏差是真实模型 F(X) 与对数据集集合取平均值的差异。偏差是衡量集合中所有数据集的平均预测与真实回归函数 F(X) 差异的程度。方差由给出。它是衡量给定数据集的解在所有数据集的平均值周围变化的程度的度量。因此，方差是衡量函数对特定数据集 D 选择敏感性的度量。第三项噪声，如前所述，是观察值与真实回归函数之间的差异的期望，在所有 X 和 Y 的值上。将这些放在一起，我们可以写出以下公式：

机器学习的目标是学习从数据中学习函数，以最小化预期的损失 E[L]。可以通过在模型中保持更多的基函数来不断减少偏差，从而增加模型复杂度。然而，由于每个模型参数都是从给定的数据集中学习得到的，因此模型越复杂，其参数估计对数据集的敏感性就越高。这导致复杂模型具有更大的方差。因此，在任何监督机器学习任务中，模型偏差和模型复杂度之间存在权衡。必须选择一个最优复杂度的模型，以最小化未见数据集上的预测误差。在经典或频率论方法中，这是通过将标记数据分为三组来实现的。第一组是训练集，第二组是验证集，第三组是测试集。使用训练集训练的不同复杂度的模型将使用验证数据集进行评估，以选择最优复杂度的模型。然后，最终使用测试集来估计预测误差。

选择最优复杂度的模型

有不同的方法来选择具有适当复杂性的模型，以便在未见数据上的预测误差更小。让我们在线性回归模型的背景下讨论这些方法中的每一个。

子集选择

在子集选择方法中，只选择整个变量集中对模型显著的那部分变量。这不仅通过减少模型方差来提高模型的预测精度，而且从解释的角度来看也是有益的。有几种不同的子集选择方法，但以下两种是最常用的方法：

向前选择法：在向前选择法中，从一个没有变量（仅包含截距）的模型开始，通过使用贪婪算法，逐个添加其他变量。在每一步中，选择最能改善拟合的变量添加到模型中。
向后选择法：在向后选择法中，从一个完整的模型开始，依次删除对拟合影响最小的变量。在每一步中，选择 Z 分数最小的变量进行删除。在统计学中，随机变量的 Z 分数是衡量一个元素与其均值之间标准差的一个度量。Z 分数的值较小（通常小于 2）表明变量的影响更有可能是偶然的，并且不具有统计学意义。

模型正则化

在这种方法中，向损失函数中添加一个惩罚项，该惩罚项不允许参数的大小在最小化过程中变得非常大。主要有两种实现方式：

岭回归：这种简单的正则化类型是，额外的项与由给出的参数向量的幅度成正比。具有正则化项的线性回归的损失函数可以写成以下形式：具有较大幅度的参数将对损失做出更大的贡献。因此，前面损失函数的最小化通常会产生具有较小值的参数并减少过拟合。的最优值是从验证集中找到的。
Lasso：在 Lasso 中，也添加了一个类似于岭回归的惩罚项，但这个项是每个参数模的加权和，而不是其平方：

虽然这看起来是一个简单的变化，但 Lasso 与岭回归有一些非常重要的区别。首先，项的存在使得损失函数在参数上是非线性的。与岭回归中的线性规划问题相比，相应的最小化问题称为二次规划问题，对于后者，存在闭式解。由于惩罚项的特殊形式，当系数在最小化过程中缩小，其中一些最终会变成零。因此，从某种意义上说，Lasso 也是一个子集选择问题。

关于各种子集选择和模型正则化方法的详细讨论，可以在 Trevor Hastie 等人所著的书籍中找到（本章“参考文献”部分的第 1 条参考文献）。

贝叶斯平均

到目前为止，我们已经了解到，仅仅最小化损失函数（或等价地，在正态分布的情况下最大化对数似然函数）对于开发给定问题的机器学习模型是不够的。我们必须担心模型过度拟合训练数据，这将在新的数据集上导致更大的预测误差。贝叶斯方法的主要优势在于，原则上可以摆脱这个问题，而无需使用显式的正则化和为训练和验证使用不同的数据集。这被称为贝叶斯模型平均，将在下面进行讨论。这是本章主要问题之一，“为什么机器学习的贝叶斯推理？”的答案之一。

对于这个问题，我们将对线性回归问题进行完整的贝叶斯处理。由于我们只想解释贝叶斯推理如何避免过拟合问题，我们将跳过所有的数学推导，只在此处陈述重要的结果。更多细节，感兴趣的读者可以参考 Christopher M. Bishop 的书籍（本章“参考文献”部分的第 2 条参考文献）。

线性回归方程，其中具有零均值和方差（等价于精度），可以转换为具有Y具有均值f(X)和精度的概率分布形式。因此，线性回归等价于估计正态分布的均值：

由于，其中基础函数集B(X)是已知的，并且我们假设这里的噪声参数也是一个已知的常数，因此只需要将作为一个不确定变量来完全贝叶斯处理。

贝叶斯推理的第一步是计算参数向量的后验分布。为此，我们假设的先验分布是一个M维正态分布（因为有M个分量），均值为，协方差矩阵为。正如我们在第三章中看到的，介绍贝叶斯推理，这对应于对先验取共轭分布：

相应的后验分布由以下给出：

这里，和。

这里，B是一个由基础向量B在不同X值上堆叠而成的N x M矩阵，如图所示：

现在我们已经得到了的后验分布作为一个封闭形式的解析表达式，我们可以用它来预测Y的新值。为了得到Y预测分布的解析封闭形式，我们假设和。这对应于一个具有零均值和各向同性协方差矩阵的先验，其特征是一个精度参数。预测分布或预测新值X = x为y的概率由以下给出：

这个等式是本节的核心主题。在经典或频率派方法中，一个人从训练数据集中估计参数的特定值，并通过简单地使用来预测 y 的概率。除非使用正则化，否则这不会解决模型的过拟合问题。在贝叶斯推理中，我们通过使用从数据中学习到的参数变量的后验概率分布来积分出参数变量。这种平均将消除使用正则化或通过偏差-方差权衡将参数保持在最佳水平的必要性。这可以从线性回归问题的封闭形式表达式 P(y|x) 中看出，在我们将和的表达式代入并进行积分后。由于两者都是正态分布，积分可以解析地进行，从而得到以下简单的 P(y|x) 表达式：

这里，。

这个等式表明预测分布的方差由两个项组成。一个项是来自数据固有的噪声和第二个项来自与从数据估计模型参数的不确定性。可以证明，随着训练数据大小 N 变得非常大，第二个项会减小，并在极限时变为零。

这里展示的例子说明了贝叶斯推理的力量。由于可以通过贝叶斯平均处理参数估计中的不确定性，因此不需要保留单独的验证数据，所有数据都可以用于训练。因此，对问题的全面贝叶斯处理可以避免过拟合问题。贝叶斯推理的另一个主要优点，我们将在本节中不深入探讨，是处理机器学习模型中的潜在变量。在下一节中，我们将对各种常见的机器学习任务进行高层次概述。

常见机器学习任务的概述

本节是以下章节的序言，我们将详细讨论不同的机器学习技术。从高层次来看，机器学习试图解决的任务只有少数几个。然而，对于每个这样的任务，都有几种方法和算法可供选择。

任何机器学习中的典型任务通常是以下之一：

分类
回归
聚类
关联规则
预测
维度降低
密度估计

在分类中，目标是把新的数据点分配到预定的类别之一。通常，这是一个监督学习或半监督学习问题。用于分类的著名机器学习算法包括逻辑回归、支持向量机（SVM）、决策树、朴素贝叶斯、神经网络、Adaboost 和随机森林。在这里，朴素贝叶斯是基于贝叶斯推理的方法。其他算法，如逻辑回归和神经网络，也已在贝叶斯框架中实现。

回归可能是最常见的机器学习问题。它用于确定一组输入变量（通常是连续变量）与一个连续的输出（因变量）之间的关系。我们在上一节中详细讨论了线性回归的最简单例子。更复杂的回归例子包括广义线性回归、样条回归、使用神经网络的非线性回归、支持向量回归和贝叶斯网络。回归的贝叶斯表述包括贝叶斯网络和贝叶斯线性回归。

聚类是无监督学习的经典例子。在这里，目标是根据数据的某些特征将数据集中的相似项分组。簇的数量事先并不知道。因此，聚类更像是模式检测问题。著名的聚类算法包括 K 均值聚类、层次聚类和潜在狄利克雷分配（LDA）。在这里，LDA 被表述为一个贝叶斯推理问题。其他使用贝叶斯推理的聚类方法包括贝叶斯混合模型。

关联规则挖掘是一种无监督方法，它在大规模数据交易中寻找共同出现的项。基于关联规则挖掘的市场篮子分析，是寻找在超市中一起销售的物品。Apriori 算法和频繁模式匹配算法是用于关联规则挖掘的两种主要方法。

预测与回归类似，不同之处在于数据是时间序列，其中存在具有不同时间戳值的观测值，目标是基于当前和过去值预测未来的值。为此，可以使用 ARIMA、神经网络和动态贝叶斯网络等方法。

机器学习中的一个基本问题被称为维度诅咒。由于机器学习模型中可能存在大量特征，为了估计模型参数，通常需要进行大量维度的搜索和优化。在更高维度的空间中，数据通常非常稀疏。这可能会使寻找最佳参数变得非常低效。为了避免这个问题，人们试图将这个高维空间投影到一个包含少数重要变量的低维空间中。然后可以使用这些低维变量作为特征。降维的两个著名例子是主成分分析和自组织映射。

通常，从少量观察数据中直接估计总体概率分布，而不使用任何参数模型，以进行推断。这被称为密度估计。密度估计的最简单形式是直方图，尽管它对于许多实际应用来说并不充分。更复杂的密度估计包括核密度估计（KDE）和矢量量化。

参考文献

Friedman J., Hastie T., and Tibshirani R. 《统计学习的要素 – 数据挖掘、推理和预测》. Springer Series in Statistics. 2009
Bishop C.M. 《模式识别与机器学习（信息科学和统计学）》. Springer. 2006. ISBN-10: 0387310738

摘要

在本章中，我们概述了机器学习是什么以及它的一些高级任务。我们还讨论了贝叶斯推理在机器学习中的重要性，特别是在它如何帮助避免重要问题，例如模型过拟合以及如何选择最佳模型。在接下来的章节中，我们将详细了解一些贝叶斯机器学习方法。

第五章：贝叶斯回归模型

在上一章中，我们详细介绍了贝叶斯线性回归的理论。在本章中，我们将通过一个示例问题来说明它如何应用于实际情况。为此，我们将使用 R 中的广义线性模型（GLM）包。首先，我们将向读者简要介绍 GLM 的概念。

广义线性回归

回想一下，在线性回归中，我们假设依赖变量Y和自变量X之间的以下函数形式：

这里，是一组基函数，是参数向量。通常，假设，因此代表一个截距或偏差项。此外，假设是一个按照均值为零、方差的正态分布分布的噪声项。我们还展示了这导致以下方程：

可以将前面的方程推广，不仅包括噪声的正态分布，还包括指数族中的任何分布（参考本章“参考文献”部分的第 1 条参考文献）。这是通过定义以下方程来实现的：

这里，g被称为连接函数。众所周知，如逻辑回归、对数线性模型、泊松回归等模型，都是广义线性模型（GLM）的特殊情况。例如，在普通线性回归的情况下，连接函数将是。对于逻辑回归，它是，即逻辑函数的逆，而对于泊松回归，它是。

在 GLM 的贝叶斯公式中，与普通线性回归不同，没有封闭形式的解析解。需要指定回归系数的先验概率。然后，通常通过蒙特卡洛模拟获得它们的后验概率。

arm 包

在本章中，为了说明贝叶斯回归模型，我们将使用 R 的arm包。这个包是由 Andrew Gelman 及其同事开发的，可以从CRAN.R-project.org/package=arm网站下载。

arm 包具有bayesglm函数，该函数实现了具有独立正态、t 或 Cauchy 先验分布的贝叶斯广义线性模型，用于模型系数。我们将使用此函数构建贝叶斯回归模型。

能源效率数据集

我们将使用来自 UCI 机器学习仓库的能源效率数据集来展示贝叶斯回归（本章参考文献部分的参考 2）。数据集可以从网站archive.ics.uci.edu/ml/datasets/Energy+efficiency下载。该数据集包含具有不同建筑参数的建筑能源效率测量值。测量了两个能源效率参数：供暖负荷（Y1）和冷却负荷（Y2）。

使用的建筑参数包括：相对紧凑度（X1）、表面积（X2）、墙体面积（X3）、屋顶面积（X4）、整体高度（X5）、朝向（X6）、玻璃面积（X7）和玻璃面积分布（X8）。我们将尝试使用普通回归和贝叶斯回归，利用 arm 包的glm函数，将供暖负荷作为所有建筑参数的函数进行预测。我们将表明，对于同一数据集，贝叶斯回归给出了显著更小的预测区间。

建筑参数与能源效率的回归

在本节中，我们将对建筑的能源效率指标进行线性回归，将供暖负荷（Y1）作为建筑参数的函数。进行初步的描述性分析以找出哪些建筑变量具有统计学意义将是有用的。为此，我们将首先创建Y1和所有X变量的双变量图。我们还将计算Y1和所有X变量之间的 Spearman 相关系数。执行这些任务的 R 脚本如下：

>library(ggplot2)
>library(gridExtra)

>df <- read.csv("ENB2012_data.csv",header = T)
>df <- df[,c(1:9)]
>str(df)
>df[,6] <- as.numeric(df[,6])
>df[,8] <- as.numeric(df[,8])

>attach(df)
>bp1 <- ggplot(data = df,aes(x = X1,y = Y1)) + geom_point()
>bp2 <- ggplot(data = df,aes(x = X2,y = Y1)) + geom_point()
>bp3 <- ggplot(data = df,aes(x = X3,y = Y1)) + geom_point()
>bp4 <- ggplot(data = df,aes(x = X4,y = Y1)) + geom_point()
>bp5 <- ggplot(data = df,aes(x = X5,y = Y1)) + geom_point()
>bp6 <- ggplot(data = df,aes(x = X6,y = Y1)) + geom_point()
>bp7 <- ggplot(data = df,aes(x = X7,y = Y1)) + geom_point()
>bp8 <- ggplot(data = df,aes(x = X8,y = Y1)) + geom_point()
>grid.arrange(bp1,bp2,bp3,bp4,bp5,bp6,bp7,bp8,nrow = 2,ncol = 4)
>detach(df)

>cor.val <- cor(df[,1:8],df[,9],method = "spearman")
>cor.val           
[,1]
X1  0.622134697
X2 -0.622134697
X3  0.471457650
X4 -0.804027000
X5  0.861282577
X6 -0.004163071
X7  0.322860320
X8  0.068343464

从 b 图和相关系数值中，我们可以得出结论，变量 X6 和 X8 对 Y1 没有显著影响，因此可以从模型中删除。

普通回归

在我们查看贝叶斯线性回归之前，让我们先进行普通线性回归。以下 R 代码使用基础 R 中的lm函数在训练数据上拟合线性回归模型，并在测试数据集上预测Y1的值：

>#Removing X6 and X8 since they don't have significant correlation with Y1
>df <- df[,c(1,2,3,4,5,7,9)]
>str(df)

>Splitting data set to Train and Test set in the ratio 80:20
>set.seed(123)
>samp <- sample.int(nrow(df),as.integer(nrow(df)*0.2),replace = F)
>dfTest <- df[samp,]
>dfTrain <- df[-samp,]
>xtest <- dfTest[,1:6]
>ytest <- dfTest[,7]

>library(arm)
>attach(dfTrain)

>#Ordinary Multivariate Regression
>fit.ols <- lm(Y1 ~ X1 + X2 + X3 + X4 + X5 + X7,data = dfTrain)
>summary(fit.ols)
>fit.coeff <- fit.ols$coefficients
>ypred.ols <- predict.lm(fit.ols,xtest,interval = "prediction",se.fit = T)
>ypred.ols$fit
>yout.ols <- as.data.frame(cbind(ytest,ypred.ols$fit))
>ols.upr <- yout.ols$upr
>ols.lwr <- yout.ols$lwr

贝叶斯回归

要执行贝叶斯线性回归，我们使用 arm 包的bayesglm()函数。正如我们在引言中所描述的，对于 GLM，如果我们选择gaussian（与正态分布相同）作为家族，并且选择identity作为连接函数，那么 GLM 就等同于普通线性回归。因此，如果我们使用bayesglm()函数并带有gaussian家族和identity连接函数，那么我们就是在执行贝叶斯线性回归。

对于贝叶斯模型，我们需要指定一个先验分布。对于高斯分布，默认设置是 prior.mean = 0，prior.scale = NULL，和 prior.df = Inf。以下 R 代码可用于贝叶斯线性回归：

>fit.bayes <- bayesglm(Y1 ~ X1 + X2 + X3 + X4 + X5 + X7,family=gaussian(link=identity),data=dfTrain,prior.df = Inf,prior.mean = 0,prior.scale = NULL,maxit = 10000)
>ypred.bayes <- predict.glm(fit.bayes,newdata = xtest,se.fit = T)
>ypred.bayes$fit

为了比较普通回归和贝叶斯回归的结果，我们在一个图表上绘制测试数据的预测值以及两种方法的预测误差。为此，我们将使用ggplot2包：

>library(ggplot2)
>library(gridExtra)
>yout.ols <- as.data.frame(cbind(ytest,ypred.ols$fit))
>ols.upr <- yout.ols$upr
>ols.lwr <- yout.ols$lwr

>p.ols <- ggplot(data = yout.ols,aes(x = yout.ols$ytest,y = yout.ols$fit)) + geom_point() + ggtitle("Ordinary Regression Prediction on Test Data") + labs(x = "Y-Test",y = "Y-Pred")
>p.ols + geom_errorbar(ymin = ols.lwr,ymax = ols.upr)yout.bayes <- as.data.frame(cbind(ytest,ypred.bayes$fit))
>names(yout.bayes) <- c("ytest","fit")
>critval <- 1.96 #approx for 95% CI
>bayes.upr <- ypred.bayes$fit + critval * ypred.bayes$se.fit
>bayes.lwr <- ypred.bayes$fit - critval * ypred.bayes$se.fit

>p.bayes <- ggplot(data = yout.bayes,aes(x = yout.bayes$ytest,y = yout.bayes$fit)) + geom_point() + ggtitle("Bayesian Regression Prediction on Test Data") + labs(x = "Y-Test",y = "Y-Pred")
>p.bayes + geom_errorbar(ymin = bayes.lwr,ymax = bayes.upr)

>p1 <-  p.ols + geom_errorbar(ymin = ols.lwr,ymax = ols.upr)
>p2 <-  p.bayes + geom_errorbar(ymin = bayes.lwr,ymax = bayes.upr)

>grid.arrange(p1,p2,ncol = 2)

可以看到，与普通回归相比，贝叶斯方法给出了更加紧凑的 95%置信预测区间。这是因为，在贝叶斯方法中，我们计算参数的分布。预测是通过从后验分布中抽取一系列值并平均得到最终预测和置信区间来进行的。

后验分布的模拟

如果想找出模型参数的后验分布，arm 包中的sim( )函数就变得很有用。以下 R 脚本将模拟参数的后验分布并生成一系列直方图：

>posterior.bayes <- as.data.frame(coef(sim(fit.bayes)))
>attach(posterior.bayes)

>h1 <- ggplot(data = posterior.bayes,aes(x = X1)) + geom_histogram() + ggtitle("Histogram X1")
>h2 <- ggplot(data = posterior.bayes,aes(x = X2)) + geom_histogram() + ggtitle("Histogram X2")
>h3 <- ggplot(data = posterior.bayes,aes(x = X3)) + geom_histogram() + ggtitle("Histogram X3")
>h4 <- ggplot(data = posterior.bayes,aes(x = X4)) + geom_histogram() + ggtitle("Histogram X4")
>h5 <- ggplot(data = posterior.bayes,aes(x = X5)) + geom_histogram() + ggtitle("Histogram X5")
>h7 <- ggplot(data = posterior.bayes,aes(x = X7)) + geom_histogram() + ggtitle("Histogram X7")
>grid.arrange(h1,h2,h3,h4,h5,h7,nrow = 2,ncol = 3)

>detach(posterior.bayes)

练习

使用来自 UCI 机器学习仓库的多变量数据集 Auto MPG（本章参考文献部分的第 3 个参考文献）。数据集可以从网站archive.ics.uci.edu/ml/datasets/Auto+MPG下载。该数据集描述了在美国城市行驶的汽车的每加仑燃油消耗量（mpg）。从包含数据集的文件夹中下载两个文件：auto-mpg.data和auto-mpg.names。auto-mpg.data文件包含数据，格式为空格分隔。auto-mpg.names文件包含关于数据集的多个细节，包括每列的变量名。使用 OLS 和贝叶斯 GLM 构建燃油效率回归模型，作为排量（disp）、马力（hp）、重量（wt）和加速度（accel）的函数。使用 OLS 模型和贝叶斯 GLM 模型（使用bayesglm函数）预测测试数据集中的 mpg 值。计算 OLS 和贝叶斯 GLM 的均方根误差（RMSE）值，并比较两种方法的准确性和预测区间。

参考文献

Friedman J., Hastie T., and Tibshirani R. 统计学习的要素 – 数据挖掘、推理与预测. 施普林格统计学系列. 2009
Tsanas A. and Xifara A. "使用统计机器学习工具准确估计住宅建筑的能源性能". 能源与建筑. 第 49 卷，第 560-567 页. 2012
Quinlan R. "Combining Instance-based and Model-based Learning". In: Tenth International Conference of Machine Learning. 236-243. University of Massachusetts, Amherst. Morgan Kaufmann. 1993. Original dataset is from StatLib library maintained by Carnegie Mellon University.

摘要

在本章中，我们展示了如何使用能源效率数据集和 arm 包中的bayesglm函数，通过更紧的置信区间来提高贝叶斯回归在预测中的实用性。我们还学习了如何使用同一 R 包中的sim函数来模拟后验分布。在下一章中，我们将学习贝叶斯分类。

第六章。贝叶斯分类模型

我们在第四章，使用贝叶斯推理进行机器学习中介绍了分类机器学习任务，并指出分类的目标是将数据记录分配到预定的类别之一。分类是机器学习中最被研究的任务之一，并且有几种已建立的先进方法。这些包括逻辑回归模型、支持向量机、随机森林模型和神经网络模型。在有足够标记训练数据的情况下，这些模型可以在许多实际问题中实现超过 95% 的准确率。

那么，显然的问题是，为什么你需要使用贝叶斯方法进行分类？对此有两个答案。一是，通常很难获得大量标记数据用于训练。当给定问题中有数百或数千个特征时，通常需要大量训练数据来避免过拟合。贝叶斯方法可以通过贝叶斯平均来克服这个问题，因此只需要少量到中等大小的训练数据。其次，大多数方法，如 SVM 或 NN，就像黑盒机器。它们会给出非常准确的结果，但很少能洞察到哪些变量对例子很重要。在许多实际问题上，例如疾病的诊断，识别主要原因非常重要。因此，黑盒方法是不够的。贝叶斯方法有一个固有的特性，称为自动相关性确定（ARD），可以通过它来识别问题中的重要变量。

在本章中，我们将讨论两种贝叶斯分类模型。第一个是流行的朴素贝叶斯文本分类方法。第二个是贝叶斯逻辑回归模型。在我们讨论这些模型之前，让我们回顾一下在分类任务中常用的一些性能指标。

分类性能指标

为了更容易理解这些概念，让我们以二元分类为例，这里的任务是判断一个输入特征向量属于两种状态之一：-1 或 1。假设 1 是正类，-1 是负类。预测输出只包含 -1 或 1，但可能存在两种类型的错误。测试集中的某些 -1 可能被预测为 1。这被称为假阳性或第一类错误。同样，测试集中的某些 1 可能被预测为 -1。这被称为假阴性或第二类错误。在二元分类的情况下，这两种类型的错误可以用下面的混淆矩阵来表示。

混淆矩阵	预测类别
正类	负类
---	---
实际类别	正类
负类	FP

从混淆矩阵中，我们可以推导出以下性能指标：

精确度：这给出了输出预测为正的准确答案的百分比
召回率：这给出了测试数据集中正确预测的正类别的百分比
F 分数：这是精确度和召回率的几何平均值
真阳性率：这与召回率相同
假阳性率：这给出了将负类别分类为正类别的百分比

此外，Tpr 被称为 灵敏度，而 1 - Fpr 被称为分类器的 特异性。Tpr 与 Fpr（灵敏度与 1 - 特异性）的图像称为 ROC 曲线（代表 接收者操作特征 曲线）。这用于找到决定预测输出（通常是一个分数或概率）属于类别 1 或 -1 的最佳阈值（分类器的操作点）。

通常，阈值被取为 ROC 曲线的膨胀点，该曲线给出了最佳性能且错误预测最少。ROC 曲线下方的面积或 AUC 是分类器性能的另一个衡量标准。对于纯随机模型，ROC 曲线将沿着对角线是一条直线。相应的 AUC 值将为 0.5。AUC 大于 0.8 的分类器将被认为是好的，尽管这很大程度上取决于要解决的问题。

朴素贝叶斯分类器

Naïve Bayes 的名字来源于模型中的基本假设，即特定特征的概率在给定类别标签的情况下与其他任何特征独立。这暗示了以下内容：

使用这个假设和贝叶斯定理，可以证明，在特征给定的情况下，类别的概率如下：

这里，是通过求和所有 k 值的分子得到的归一化项。它也被称为贝叶斯证据或配分函数 Z。分类器选择一个类别标签作为目标类别，该类别最大化后验类别概率：

朴素贝叶斯分类器是文档分类的基线分类器。其中一个原因是，在给定类别标签的情况下，每个特征（单词或 m-gram）与其他特征相互独立的基本假设通常适用于文本。另一个原因是，当存在大量文档时，朴素贝叶斯分类器的扩展性很好。

朴素贝叶斯有两种实现方式。在伯努利朴素贝叶斯中，特征是二进制变量，表示一个特征（m-gram）是否存在于文档中。在多项式朴素贝叶斯中，特征是文档中 m-gram 的频率。为了避免频率为零时的问题，通过对特征向量添加 1 来进行拉普拉斯平滑。让我们详细看看多项式朴素贝叶斯。

令为特征在训练数据中类别中出现的次数。那么，给定类别标签观察到特征向量的似然函数如下：

在这里，是在类别中观察到特征的概率。

使用贝叶斯规则，给定特征向量 X 观察到类别的后验概率如下：

对等式两边取对数并忽略常数项 Z，我们得到以下结果：

因此，通过对后验分布取对数，我们将问题转化为一个线性回归模型，其中作为需要从数据中确定的系数。这可以很容易地解决。通常，人们使用 TF-IDF（词频乘以逆频率）而不是词频，并将文档长度归一化以提高模型性能。

R 包 e1071 (统计学部门的各种函数) 由 T.U. Wien 编写，其中包含了对朴素贝叶斯的 R 实现。对于本章，我们将使用来自 UCI 机器学习仓库的 SMS 炎言数据集（本章“参考文献”部分的第 1 个参考文献）。该数据集包含从英国论坛 Grumbletext 收集的 425 条短信垃圾信息，消费者可以在此提交垃圾短信。数据集还包含来自新加坡国立大学维护的 SMS 语料库的 3375 条正常（非垃圾）短信。

该数据集可以从 UCI 机器学习仓库下载（archive.ics.uci.edu/ml/datasets/SMS+Spam+Collection）。假设我们已经将其保存为 R 的工作目录中的文件 SMSSpamCollection.txt（实际上，您需要将其在 Excel 中打开并保存为制表符分隔的文件，以便 R 正确读取）。

然后，将文件读入 tm（文本挖掘）包的命令如下：

>spamdata <-read.table("SMSSpamCollection.txt",sep="\t",stringsAsFactors = default.stringsAsFactors())

我们将首先将因变量 y 和自变量 x 分离，并将数据集按 80:20 的比例分为训练集和测试集，使用以下 R 命令：

>samp<-sample.int(nrow(spamdata),as.integer(nrow(spamdata)*0.2),replace=F)
>spamTest <-spamdata[samp,]
>spamTrain <-spamdata[-samp,]
>ytrain<-as.factor(spamTrain[,1])
>ytest<-as.factor(spamTest[,1])
>xtrain<-as.vector(spamTrain[,2])
>xtest<-as.vector(spamTest[,2])

由于我们处理的是文本文档，在我们可以使用数据为任何机器学习模型之前，我们需要进行一些标准的预处理。我们可以使用 R 中的 tm 包来完成这个目的。在下一节中，我们将对此进行详细描述。

使用 tm 包进行文本处理

tm 包提供了数据导入、语料库处理、预处理、元数据管理和创建词-文档矩阵的方法。数据可以导入到 tm 包中，无论是从目录、每个组件都是一个文档的向量，还是从数据框中。tm 的基本数据结构是一个抽象的文本文档集合，称为 Corpus。它有两种实现方式；一种是将数据存储在内存中，称为 VCorpus（易失性语料库），另一种是将数据存储在硬盘上，称为 PCorpus（永久性语料库）。

我们可以使用以下 R 命令创建我们的 SMS 垃圾邮件数据集的语料库；在此之前，您需要使用 R 中的 install.packages("packagename") 命令安装 tm 包和 SnowballC 包：

>library(tm)
>library(SnowballC)
>xtrain <- VCorpus(VectorSource(xtrain))

首先，我们需要进行一些基本的文本处理，例如删除额外的空白字符，将所有单词转换为小写，删除停用词，以及进行词干提取。这可以通过使用 tm 包中的以下函数来实现：

>#remove extra white space
>xtrain <- tm_map(xtrain,stripWhitespace)
>#remove punctuation
>xtrain <- tm_map(xtrain,removePunctuation)
>#remove numbers
>xtrain <- tm_map(xtrain,removeNumbers)
>#changing to lower case
>xtrain <- tm_map(xtrain,content_transformer(tolower))
>#removing stop words
>xtrain <- tm_map(xtrain,removeWords,stopwords("english"))
>#stemming the document
>xtrain <- tm_map(xtrain,stemDocument)

最后，数据被转换成机器学习模型可以消费的形式。这就是所谓的文档-词矩阵形式，其中每个文档（在这种情况下为 SMS）是一行，所有文档中出现的术语是列，每个单元格中的条目表示每个单词在一个文档中出现的次数：

>#creating Document-Term Matrix
>xtrain <- as.data.frame.matrix(DocumentTermMatrix(xtrain))

同样的处理过程也应用于 xtest 数据集。我们将 y 转换为因子，将 xtrain 转换为数据框的原因是为了匹配 e1071 包中朴素贝叶斯分类器的输入格式。

模型训练和预测

您需要首先从 CRAN 安装 e1071 包。可以使用 naiveBayes() 函数来训练朴素贝叶斯模型。该函数可以通过两种方法调用。以下为第一种方法：

>naiveBayes(formula,data,laplace=0,…,subset,na.action=na.pass)

其中 formula 代表独立变量的线性组合，用于预测以下类别：

>class ~ x1+x2+…

此外，data代表一个数据框或由分类和数值变量组成的列联表。

如果我们有一个向量y作为类别标签和依赖变量作为数据框x，那么我们可以使用调用函数的第二种方法，如下所示：

>naiveBayes(x,y,laplace=0,…)

在我们的示例中，我们将使用调用方法的第二种方法。一旦我们有一个训练好的模型，它是一个名为naiveBayes的 R 对象，我们就可以预测新实例的类别，如下所示：

>predict(object,newdata,type=c(class,raw),threshold=0.001,eps=0,…)

因此，我们可以在我们的训练数据集上训练朴素贝叶斯模型，并使用以下命令在测试数据集上进行评分：

>#Training the Naive Bayes Model
>nbmodel <- naiveBayes(xtrain,ytrain,laplace=3)
>#Prediction using trained model
>ypred.nb <- predict(nbmodel,xtest,type = "class",threshold = 0.075)
>#Converting classes to 0 and 1 for plotting ROC
>fconvert <- function(x){
     if(x == "spam"){ y <- 1}
   else {y <- 0}
   y
}

>ytest1 <- sapply(ytest,fconvert,simplify = "array")
>ypred1 <- sapply(ypred.nb,fconvert,simplify = "array")
>roc(ytest1,ypred1,plot = T)

在这里，展示了该模型和数据的 ROC 曲线。这是使用 CRAN 中的 pROC 包生成的：

>#Confusion matrix
>confmat <- table(ytest,ypred.nb)
>confmat
pred.nb
ytest  ham spam
  ham  143  139
  spam   9   35

从 ROC 曲线和混淆矩阵中，可以选择分类器最佳阈值和精确度、召回率指标。请注意，这里所示示例仅用于说明目的。模型需要进一步调整以提高准确性。

我们还可以打印出两个类别中最频繁出现的单词（模型特征）及其由模型生成的后验概率。这将使我们对模型练习有更直观的感觉。以下 R 代码完成这项工作：

>tab <- nbmodel$tables
>fham <- function(x){
  y <- x[1,1]
  y
}
>hamvec <- sapply(tab,fham,simplify = "array")
>hamvec <- sort(hamvec,decreasing = T)

>fspam <- function(x){
  y <- x[2,1]
  y
}
>spamvec <- sapply(tab,fspam,simplify = "array")
>spamvec <- sort(spamvec,decreasing = T)
>prb <- cbind(spamvec,hamvec)
>print.table(prb)

输出表格如下：

| word | Prob(word | spam) | Prob(word | ham) |
| --- | --- | --- |
| call | 0.6994 | 0.4084 |
| free | 0.4294 | 0.3996 |
| now | 0.3865 | 0.3120 |
| repli | 0.2761 | 0.3094 |
| text | 0.2638 | 0.2840 |
| spam | 0.2270 | 0.2726 |
| txt | 0.2270 | 0.2594 |
| get | 0.2209 | 0.2182 |
| stop | 0.2086 | 0.2025 |

表格显示，例如，给定一个文档是垃圾邮件，该文档中出现单词call的概率为 0.6994，而该单词出现在正常文档中的概率仅为 0.4084。

贝叶斯逻辑回归模型

逻辑回归的名称来源于回归的依赖变量是逻辑函数。它是响应变量为二元变量（例如，欺诈或非欺诈、点击或未点击等）的问题中广泛使用的模型之一。

逻辑函数由以下方程定义：

它具有特定的特征，即当y从变化到时，函数值从 0 变化到 1。因此，逻辑函数非常适合模拟任何二元响应，因为输入信号的变化。

逻辑函数的逆称为logit。它定义如下：

在逻辑回归中，y被视为解释变量X的线性函数。因此，逻辑回归模型可以定义为以下：

在这里，是基函数集合，是模型参数，如第四章中所述，使用贝叶斯推理进行机器学习。从第五章中广义线性模型（GLM）的定义，贝叶斯回归模型，可以立即识别出逻辑回归是具有logit函数作为连接函数的 GLM 的特殊情况。

与线性回归相比，贝叶斯处理逻辑回归更为复杂。在这里，似然函数由逻辑函数的乘积组成；每个数据点对应一个。为了计算后验概率，必须将这个函数与先验概率相乘（以得到贝叶斯公式的分母）。一种方法是使用拉普拉斯近似，如第三章中所述，介绍贝叶斯推理。读者可能还记得，在拉普拉斯近似中，后验概率被近似为关于后验概率最大值的正态（高斯）分布。这是通过首先找到最大后验概率（MAP）解，然后计算 MAP 解周围负对数似然函数的二阶导数来实现的。感兴趣的读者可以在 D.J.C. MacKay 的论文中找到逻辑回归拉普拉斯近似的详细内容（本章参考文献部分的第 2 条参考文献）。

除了使用解析近似，Polson 和 Scott 最近提出了一种使用数据增强策略的完全贝叶斯处理方法（本章参考文献部分的第 3 条参考文献）。作者们将他们的方法实现在了 R 包：BayesLogit 中。我们将使用这个包在本章中展示贝叶斯逻辑回归。

BayesLogit R 包

该包可以从 CRAN 网站cran.r-project.org/web/packages/BayesLogit/index.html下载。该包包含一个logit函数，可用于执行贝叶斯逻辑回归。调用此函数的语法如下：

>logit(Y,X,n=rep(1,length(Y) ),m0=rep(0,ncol(X) ),P0=matrix(0,nrow=ncol(X),ncol=ncol(X) ),samp=1000,burn=500)

在这里，Y 是一个包含响应值的 N-维向量；X 是一个包含独立变量值的 N x P 维矩阵；n 是一个 N-维向量，是一个 P-维先验均值，是一个 P x P 维先验精度。其他两个参数与 MCMC 模拟参数相关。保存的 MCMC 模拟次数用 samp 表示，运行开始前丢弃的 MCMC 模拟次数用 burn 表示。

数据集

为了说明贝叶斯逻辑回归，我们使用了 UCI 机器学习仓库中的 Parkinsons 数据集(archive.ics.uci.edu/ml/datasets/Parkinsons)。Little 等人使用该数据集通过分析声音障碍来检测帕金森病（本章“参考文献”部分的第 4 条参考文献）。数据集包含 31 个人的声音测量数据，其中 23 人患有帕金森病。有 195 行对应于单个个体的多次测量。测量可以分组为以下集合：

声音基频
震动
Shimmer
噪声与音调成分的比率
非线性动力学复杂性度量
信号分形尺度指数
基频变化的非线性度量

总共有 22 个数值属性。

准备训练和测试数据集

在我们能够训练贝叶斯逻辑模型之前，我们需要对数据进行一些预处理。数据集包含来自同一个体的多次测量。在这里，我们取所有观察值，每个来自一组被采样个体的集合，以创建训练和测试集。此外，我们需要将因变量（类别标签 Y）从自变量（X）中分离出来。以下 R 代码完成这项工作：

>#install.packages("BayesLogit") #One time installation of package
>library(BayesLogit)
>PDdata <- read.table("parkinsons.csv",sep=",",header=TRUE,row.names = 1)
>rnames <- row.names(PDdata)
>cnames <- colnames(PDdata,do.NULL = TRUE,prefix = "col")
>colnames(PDdata)[17] <- "y"
>PDdata$y <- as.factor(PDdata$y)

>rnames.strip <- substr(rnames,10,12)
>PDdata1 <- cbind(PDdata,rnames.strip)
>rnames.unique <- unique(rnames.strip)
>set.seed(123)
>samp <- sample(rnames.unique,as.integer(length(rnames.unique)*0.2),replace=F)
>PDtest <- PDdata1[PDdata1$rnames.strip %in% samp,-24]  # -24 to remove last column
>PDtrain <- PDdata1[!(PDdata1$rnames.strip %in% samp),-24] # -24 to remove last column
>xtrain <- PDtrain[,-17]
>ytrain <- PDtrain[,17]
>xtest <- PDtest[,-17]
>ytest<- PDtest[,17]

使用贝叶斯逻辑模型

我们可以使用 xtrain 和 ytrain 通过 logit( ) 函数训练贝叶斯逻辑回归模型：

>blmodel <- logit(ytrain,xtrain,n=rep(1,length(ytrain)),m0 = rep(0,ncol(xtrain)),P0 = matrix(0,nrow=ncol(xtrain),ncol=ncol(xtrain)),samp = 1000,burn = 500)

summary( ) 函数将给出拟合模型的概述：

>summary(blmodel)

为了预测新数据集的 Y 值，我们需要编写一个自定义脚本，如下所示：

>psi <- blmodel$beta %*% t(xtrain)  # samp x n
>p   <- exp(psi) / (1 + exp(psi) )  # samp x n
>ypred.bayes <- colMeans(p)

预测误差可以通过与 ytest 中存在的实际 Y 值进行比较来计算：

>table(ypred.bayes,ytest)

可以使用 pROC 包绘制 ROC 曲线，如下所示：

>roc(ytrain,ypred.bayes,plot = T)

ROC 曲线的 AUC 为 0.942，表明分类精度良好。再次强调，这里展示的模型是为了说明目的，并未调整以获得最佳性能。

练习

在这个练习中，我们将使用来自 UCI 机器学习仓库的 DBWorld 电子邮件数据集来比较朴素贝叶斯和 BayesLogit 方法的相对性能。该数据集包含 64 封来自 DBWorld 通讯的电子邮件，任务是将电子邮件分类为会议公告或其他一切。该数据集的参考是 Michele Filannino 教授的课程（本章参考文献部分的第 5 条）。数据集可以从 UCI 网站下载，网址为archive.ics.uci.edu/ml/datasets/DBWorld+e-mails#。

使用这两种方法都需要对数据集进行一些预处理。数据集采用 ARFF 格式。您需要下载外国的 R 包(cran.r-project.org/web/packages/foreign/index.html)，并使用其中的read.arff()方法将文件读入 R 数据框。

参考文献

Almeida T.A., Gómez Hidalgo J.M. 和 Yamakami A. "对短信垃圾邮件过滤研究的新贡献：新的收集和结果". 在：2011 年 ACM 文档工程研讨会（DOCENG'11）. 加利福尼亚州山景城，美国. 2011
MacKay D.J.C. "将证据框架应用于分类网络". 神经计算 4(5)
"使用 Pólya-Gamma 潜在变量进行逻辑模型的贝叶斯推理". 美国统计学会杂志. 第 108 卷，第 504 期，第 1339 页. 2013
Costello D.A.E., Little M.A., McSharry P.E., Moroz I.M. 和 Roberts S.J. "利用非线性回溯和分形尺度特性进行语音障碍检测". 生物医学工程在线. 2007
Filannino M. "使用非常小语料库的 DBWorld 电子邮件分类". 机器学习课程项目. 曼彻斯特大学. 2011

摘要

在本章中，我们讨论了使用贝叶斯推理进行分类任务的多种优点。我们回顾了一些用于分类任务的常见性能指标。我们还学习了两种基本且流行的分类方法，即朴素贝叶斯和逻辑回归，这两种方法都采用了贝叶斯方法实现。在了解了某些重要的贝叶斯监督机器学习技术之后，在下一章中，我们将讨论一些无监督的贝叶斯模型。

第七章。无监督学习的贝叶斯模型

我们在前两章中讨论的机器学习模型具有一个共同特征：它们需要包含真实值的训练数据。这意味着包含谓词或因变量的真实值的数据集通常需要手动标记。这种使用标记数据训练算法的机器学习被称为监督学习。这种类型的机器学习在预测准确性方面表现出色。实际上，这是大多数使用机器学习的工业系统中事实上的方法。然而，这种方法的一个缺点是，当人们想要用大数据集训练模型时，获取标记数据会变得困难。这在大数据时代尤为重要，因为组织从各种日志、交易和与消费者的互动中获得了大量数据；组织希望从这些数据中获得洞察力，并对消费者的兴趣做出预测。

在无监督学习方法中，学习过程中不需要标记数据。学习过程是通过识别数据集中存在的占主导地位的模式和相关性来进行的。无监督学习的常见例子包括聚类、关联规则挖掘、密度估计和降维。在聚类中，使用一种合适的算法识别数据中的自然分组，该算法利用数据点之间的某种距离度量。在关联规则挖掘中，从事务数据集中识别出在事务中经常一起出现的项目。在降维技术，如主成分分析中，包含大量变量（维度）的原始数据集被投影到一个较低维度的空间，其中数据中包含最大信息。尽管无监督学习不需要标记的训练数据，但人们需要大量的数据来学习所有感兴趣的模式，而且通常学习过程计算量更大。

在许多实际情况下，创建少量标记数据是可行的。第三种学习方法，即半监督学习，是一种利用这个小规模标记数据集，并通过合适的算法将标签传播到其余未标记的训练数据的方法。在本章中，我们将介绍无监督学习的贝叶斯方法。我们将详细讨论两个重要模型：用于聚类的高斯混合模型和用于主题模型的潜在狄利克雷分配。

贝叶斯混合模型

通常，混合模型对应于使用概率分布的混合来表示数据。最常见的混合模型如下所示：

在这里，是X的概率分布，其参数为，而代表混合中第k个成分的权重，使得。如果潜在的概率分布是正态（高斯）分布，那么混合模型被称为高斯混合模型（GMM）。因此，GMM 的数学表示如下：

在这里，我们使用了与前面章节相同的符号，其中X代表一个N-维数据向量，它表示每个观测值，数据集中有M个这样的观测值。

当簇之间存在重叠时，这种混合模型适合于聚类。GMM 的一个应用领域是计算机视觉。如果想在视频中跟踪移动对象，减去背景图像是有用的。这被称为背景减法或前景检测。GMM 用于此目的，其中每个像素的强度使用高斯分布的混合来建模（本章“参考文献”部分的参考 1）。

学习 GMM 的任务对应于学习所有成分的模型参数和混合权重。学习 GMM 的标准方法是通过使用最大似然方法。对于一个由M个观测值组成的数据集，似然函数的对数由以下给出：

与单个高斯模型不同，在 GMM 中，无法以直接的方式对参数进行对数似然函数的最大化。这是因为在这种情况下，由于难以计算和的对数，没有该导数的闭式表达式。因此，使用所谓的期望最大化（EM）算法来最大化对数似然函数。EM 算法是一个迭代算法，其中每个迭代包括两个计算：期望和最大化。EM 算法的步骤如下：

初始化参数、和，并评估对数似然函数的初始值。
在期望步骤中，使用当前参数值和从对数似然函数评估混合成分。
在最大化步骤中，使用步骤 2 中计算的的值，通过最大化对数似然来估计新的参数值和。
使用步骤 2 和 3 中估计的、和的值，计算对数似然函数的新值。
重复步骤 2-4，直到对数似然函数收敛。

在贝叶斯处理 GMM 中，通过引入一个潜在变量Z，简化了对对数似然函数的最大化。设Z为一个具有一个元素1和其余K – 1个元素为0的K-维二元随机变量。使用Z，可以写出X和Z的联合分布如下：

这里：

然后：

因此：

然后：

在这个问题中引入潜在变量Z的优势在于，对数似然的表达式被简化，其中对数直接作用于正态分布，就像单个高斯模型的情况一样。因此，最大化P(X, Z)是直接的。然而，仍然存在的问题是我们不知道Z!的值。所以，技巧是使用类似于 EM 的迭代算法，在 E 步骤中估计Z的期望值，在 M 步骤中，使用上一步估计的Z的最后一个值，找到高斯分布的参数值。贝叶斯版本的 EM 算法对于 GMM 的步骤如下：

初始化参数、和，并评估对数似然函数的初始值。
在期望步骤中，使用这些值来计算期望值。
在最大化步骤中，固定，通过最大化来估计和。
计算新的似然函数。
重复步骤 2-4，直到收敛。

在 Christopher M. Bishop 所著的书中（本章“参考文献”部分的第 2 条参考文献）可以找到对贝叶斯版本的 EM 算法和 GMM 的更详细处理。在这里，我们省略了贝叶斯 GMM 的理论处理，转而查看其在bgmm包中的 R 实现。

贝叶斯混合模型的 bgmm 包

bgmm 包是由 Przemyslaw Biecek 和 Ewa Szczurek 开发的，用于建模基因表达数据（本章“参考文献”部分的第 3 个参考文献）。可以从 CRAN 网站下载：cran.r-project.org/web/packages/bgmm/index.html。该包不仅包含 GMM 的无监督版本，还包含完全监督和半监督的实现。bgmm 包中包含以下不同的模型：

完全监督的 GMM: 这是指训练集中所有记录都有的标记数据。这包括以下内容：
- supervised( ) 函数
半监督的 GMM: 这是指训练集中所有记录的小子集可用的标记数据。这包括以下内容
- semisupervised( ) 函数
部分监督的 GMM: 这是指所有记录的小子集可用的标记数据，但这些标签是不确定的。标签的值以某种概率给出。包中包含两个用于部分监督 GMM 的函数。这包括以下内容::
- belief( ) 函数：标签的不确定性以对其组件的概率分布来表示。对于前 m 个观测值，输入一个维度为 m x k 的信念矩阵 B，其中矩阵条目表示 i^(th) 条记录具有 j^(th) 标签的概率。
- soft( ) 函数：在此方法中，定义了一个维度为 M x k 的可信度矩阵，跨越了大小为 M 的训练集中的所有记录。矩阵元素被解释为 i^(th) 条记录具有 j^(th) 标签的先验概率的权重。如果没有关于任何记录标签的特定信息，它们可以被赋予相等的权重。为了实施目的，对矩阵元素施加了一个约束：。
无监督的 GMM: 这是指任何记录都没有标记数据。这包括以下内容：
- unsupervised( ) 函数

传递给这些函数的典型参数如下：

X: 这是一个包含未标记 X 数据的数据框。
knowns: 这是一个包含标记 X 数据的数据框。
B: 这是一个信念矩阵，指定了标记记录的信念分布。B 的行数应与 knowns 相同。
P: 这是一个先验概率（可信度）的权重矩阵。
class: 这是一个包含标记记录的类别或标签的向量。
k: 这是矩阵 B 的组成部分或列数。
init.params: 这些是模型参数估计的初始值。

belief( )和soft( )函数之间的区别在于，在前一种情况下，输入是一个包含每个可能标签的先验概率值的矩阵，而在第二种情况下，输入是一个包含每个先验的权重而不是先验概率本身的矩阵。更多细节，读者请参阅 Przemyslaw Biecek 等人撰写的论文（本章“参考文献”部分的第 3 条参考文献）。

现在，让我们通过一个小示例来展示如何使用 bgmm。我们将使用来自 UCI 机器学习仓库的 ADL 数据集。这个数据集包含了 16 个志愿者佩戴在手腕上的加速度计的加速度数据。数据集和元数据详情可以在archive.ics.uci.edu/ml/datasets/Dataset+for+ADL+Recognition+with+Wrist-worn+Accelerometer找到。关于生成此数据集的 ADL 监控系统的研究工作发表在 Bruno B.等人撰写的两篇论文中（本章“参考文献”部分的第 4 和第 5 条参考文献）。

对于 bgmm 的示例，我们将在数据集目录中仅使用一个文件夹，即Brush_teeth。首先，我们将进行少量预处理，将来自不同志愿者的数据合并到一个文件中。下面的 R 脚本完成这项工作：

>#Set working directory to folder containing files (provide the correct path)
>setwd("C:/…/ADL_Dataset/HMP_Dataset/Brush_teeth")
>flist <- list.files(path = "C:/../ADL_Dataset/HMP_Dataset/Brush_teeth",pattern = "*.txt")
>all.data <- lapply(flist,read.table,sep = " ",header = FALSE)
>combined.data <- as.data.frame(do.call(rbind,all.data))
>combined.data.XZ <- combined.data[,c(1,3)]

最后一步是从加速度中选择X和Z分量来创建一个二维数据集。

以下 R 脚本调用bgmm函数并执行聚类。数据的一个简单散点图表明数据集中可能有四个聚类，选择k = 4就足够了：

>modelbgmm <- unsupervised(combined.data.XZ,k=4)
>summary(modelbgmm)
>plot.mModel(modelbgmm)

bgmm 生成的聚类可以在以下图中看到；有四个聚类，其中心由四个彩色点表示，其各自的高斯密度由椭圆表示：

使用贝叶斯推理进行主题建模

我们在第六章中看到了文本文档的监督学习（分类），即使用朴素贝叶斯模型。通常，一篇大型的文本文档，如新闻文章或短篇小说，可以包含不同的主题作为子章节。对于分类、摘要、压缩等目的，对这种文档内部统计相关性进行建模是有用的。在前一节中学习的高斯混合模型更适合数值数据，如图像，而不是文档。这是因为文档中的单词很少遵循正态分布。一个更合适的选择是多项式分布。

混合模型向文档的强大扩展是 T. Hofmann 关于概率语义索引的工作（本章“参考文献”部分的第 6 条参考文献）以及 David Blei 等人关于潜在狄利克雷分配的工作（本章“参考文献”部分的第 7 条参考文献）。在这些工作中，文档被描述为主题混合，每个主题由单词分布描述。LDA 是一个用于文本文档的生成无监督模型。LDA 的任务是从数据中学习主题分布、单词分布和混合系数的参数。下一节将简要介绍 LDA。强烈建议读者阅读 David Blei 等人的论文，以理解他们的方法。

Latent Dirichlet allocation

在 LDA 中，假设单词是文档的基本单元。一个单词是称为词汇表的一个集合中的一个元素，由索引。在这里，V表示词汇表的大小。一个单词可以通过一个单位基向量来表示，其所有分量都是零，除了对应于单词的那个分量，其值为 1。例如，词汇表中的第n个单词由一个大小为V的向量描述，其第n个分量和所有其他分量对于。同样，一个文档是由表示的N个单词的集合，一个语料库是由表示的M个文档的集合（注意，在这里文档用粗体w表示，而单词则没有粗体 w）。

如前所述，LDA 是一个语料库的生成概率模型，其中文档被表示为潜在主题上的随机混合，每个主题由单词分布表征。在 LDA 模型中生成语料库中的每个文档w时，执行以下步骤：

根据参数定义的泊松分布选择与文档大小相对应的N值：
选择参数的值，该参数表征了由参数定义的狄利克雷分布的主题分布：
对于每个N个单词
1. 根据步骤 2 中绘制的参数定义的多项式分布选择一个主题：
2. 从由描述并由条件化的多项式概率分布中选择一个单词：

给定 N、和，主题混合、主题集合 z 和单词集合 w 的联合分布如下：

注意，在这种情况下，只有 w 是可观察的（文档），而和 z 都被视为潜在（隐藏）变量。

LDA 中的贝叶斯推断问题是给定一个文档，估计潜在变量和 z 的后验密度：

如同许多贝叶斯模型一样，这种情况在分析上难以处理，因此必须使用近似技术，如 MCMC 或变分贝叶斯，来估计后验。

LDA 的 R 包

R 中有两个主要包可以用于在文档上执行 LDA。一个是 Bettina Grün 和 Kurt Hornik 开发的topicmodels包，另一个是 Jonathan Chang 开发的lda包。在这里，我们将描述这两个包。

主题模型包

主题模型包是 LDA 和相关主题模型（CTM）论文作者开发的 C 和 C++代码的接口（本章参考文献部分的 7、8 和 9）。该包中的主要函数 LDA 用于拟合 LDA 模型。可以通过以下方式调用：

>LDA(X,K,method = "Gibbs",control = NULL,model = NULL,...)

在这里，X 是一个文档-词矩阵，可以使用 tm 包生成，而 K 是主题的数量。method 是用于拟合的方法。支持两种方法：Gibbs 和 VEM。

让我们通过这个包构建 LDA 模型的小例子。所使用的数据集是来自 UCI 机器学习仓库的Reuter_50_50数据集（本章参考文献部分的 10 和 11）。数据集可以从archive.ics.uci.edu/ml/datasets/Reuter_50_50下载。对于这个练习，我们只将使用一个目录中的文档，即在C50train目录下的AlanCrosby。所需的预处理可以使用以下 R 脚本来完成；在尝试此练习之前，读者应已安装 tm 和 topicmodels 包：

>library(topicmodels)
>library(tm)
>#creation of training corpus from reuters dataset
>dirsourcetrain <- DirSource(directory = "C:/…/C50/C50train/AaronPressman")
>xtrain <- VCorpus(dirsourcetrain)
>#remove extra white space
>xtrain <- tm_map(xtrain,stripWhitespace)
>#changing to lower case
>xtrain <- tm_map(xtrain,content_transformer(tolower))
>#removing stop words
>xtrain <- tm_map(xtrain,removeWords,stopwords("english"))
>#stemming the document
>xtrain <- tm_map(xtrain,stemDocument)
>#creating Document-Term Matrix
>xtrain <-  as.data.frame.matrix(DocumentTermMatrix(xtrain))

可以使用相同的步骤从 /…/C50/C50test/ 目录创建测试数据集。

一旦我们有了文档-词矩阵 xtrain 和 xtest，可以使用以下 R 脚本来构建和测试 LDA 模型：

>#training lda model
>ldamodel <- LDA(xtrain,10,method = "VEM")
>#computation of perplexity, on training data (only with VEM method)
>perp <- perplexity(ldamodel)
>perp
[1] 407.3006

大约 100 的困惑度值表示拟合良好。在这种情况下，我们需要添加更多训练数据或更改K的值以提高困惑度。

现在让我们使用训练好的 LDA 模型来预测测试数据集上的主题：

>#extracting topics from test data)
>postprob <- posterior(ldamodel,xtest)
>postprob$topics

在这里，测试集只包含一个文件，即 42764newsML.txt。它在其由 LDA 模型产生的 10 个主题中的分布情况如下所示。

lda 包

lda 包是由 Jonathan Chang 开发的，他为后验估计实现了折叠吉布斯抽样方法。该包可以从 CRAN 网站下载，网址为 cran.r-project.org/web/packages/lda/index.html。

包中的主要函数 lda.collapsed.gibbs.sampler 使用折叠吉布斯抽样的方法来拟合三个不同的模型。这些是潜在狄利克雷分配（LDA）、监督 LDA（sLDA）和混合成员随机块模型（MMSB）。这些函数接受输入文档并返回潜在参数的点估计。这些函数可以在 R 中使用如下：

>lda.collapsed.gibbs.sampler(documents,K,vocab,num.iterations,alpha,eta,initial = NULL,burnin = NULL,compute.log.likelihood = FALSE,trace = 0L,freeze.topics = FALSE)

在这里，documents 代表一个包含文档的列表，列表的长度等于 D，而 K 是主题的数量；vocab 是一个指定词汇的字符向量；alpha 和 eta 是超参数的值。

练习

对于 Reuter_50_50 数据集，使用 lda 包中的 lda.collapsed.gibbs.sampler 函数拟合 LDA 模型，并与 topicmodels 包的性能进行比较。请注意，您需要使用 topicmodels 包中的 dtm2ldaformat( ) 函数将文档-词矩阵转换为 lda 格式，以便使用 lda 包。

参考文献

Bouwmans, T., El Baf F. 和 "Vachon B. 使用高斯混合模型进行背景建模以检测前景 – 一篇综述" (PDF). 近期计算机科学专利 1: 219-237. 2008
Bishop C.M. 模式识别与机器学习. Springer. 2006
Biecek P., Szczurek E., Tiuryn J., 和 Vingron M. "R 包 bgmm：不确定知识混合建模". 统计软件杂志. 第 47 卷，第 3 期. 2012
Bruno B., Mastrogiovanni F., Sgorbissa A., Vernazza T. 和 Zaccaria R. "基于加速度数据的人行行为识别算法分析". 在：IEEE 国际机器人与自动化会议（ICRA），第 1602-1607 页. 2013
Bruno B., Mastrogiovanni F., Sgorbissa A., Vernazza T. 和 Zaccaria R. "人体运动建模与识别：一种计算方法". 在：IEEE 国际自动化科学和工程会议（CASE），第 156-161 页. 2012
Hofmann T. "概率潜在语义索引". 在：第二十二届国际 SIGIR 会议. 1999
Blei D.M., Jordan M.I., 和 Ng A.Y. "潜在狄利克雷分配". 机器学习研究杂志 3. 993-1022. 2003
Blei D.M. 和 Lafferty J.D. "科学相关主题模型"。应用统计学年鉴。1(1)，17-35。2007
Phan X.H.，Nguyen L.M. 和 Horguchi S. "从大规模数据集中学习对短文本和网页进行分类，并使用隐藏主题"。在第 17 届国际万维网会议（WWW 2008）中。第 91-100 页。北京，中国。2008

摘要

在本章中，我们讨论了无监督和半监督机器学习背后的概念及其贝叶斯处理方法。我们学习了两个重要的贝叶斯无监督模型：贝叶斯混合模型和 LDA。我们详细讨论了贝叶斯混合模型的 bgmm 包，以及用于主题建模的 topicmodels 和 lda 包。由于无监督学习的主题非常广泛，我们只能在本章中涵盖一些贝叶斯方法，仅为了给读者一个对该主题的初步了解。我们没有涵盖使用项目标注和特征标注的半监督方法。对此感兴趣的读者应参考该主题的更专业书籍。在下一章中，我们将学习另一类重要的模型，即神经网络。

第八章. 贝叶斯神经网络

正如其名所示，人工神经网络是受生物大脑的架构和认知能力启发的统计模型。神经网络模型通常具有分层架构，每一层包含大量神经元，不同层之间的神经元是相互连接的。第一层被称为输入层，最后一层被称为输出层，中间的其他层被称为隐藏层。每个神经元的状态由与其相连的所有神经元状态的非线性函数确定。每个连接都有一个权重，该权重由包含一组输入和输出对的训练数据确定。这种神经元及其连接的分层架构存在于人脑的新皮层区域，并被认为是负责诸如感官感知和语言理解等高级功能的原因。

神经网络的第一种计算模型是由沃伦·麦克洛克和沃尔特·皮茨在 1943 年提出的。大约在同一时间，心理学家唐纳德·赫布基于神经元的兴奋和适应机制提出了一个学习假设，即赫布规则。这个假设可以总结为“一起放电的神经元，一起连接”。尽管有几位研究人员试图实现神经网络计算模型，但直到 1958 年，弗兰克·罗森布拉特才首次使用一个两层神经网络感知器创建了一个用于模式识别的算法。

1970-2010 年间，神经网络的研究和应用经历了停滞和快速发展的时期。神经网络历史上的里程碑包括 1975 年保罗·沃伯斯发明的反向传播算法，这是一种用于学习多层神经网络（也称为深度学习网络）的快速学习算法，由杰弗里·辛顿在 2006 年提出，以及在上个十年后半期使用 GPGPUs 实现处理神经网络所需的更大计算能力。

现在，神经网络模型及其应用再次在人工智能领域占据中心舞台，应用于计算机视觉、语音识别和自然语言理解。这就是本书专门用一章来探讨这个主题的原因。当我们深入了解后续章节时，贝叶斯推理在神经网络模型中的重要性将变得清晰。

两层神经网络

让我们看看双层神经网络的正式定义。我们遵循 David MacKay（本章“参考文献”部分的第 1、2 和 3 条参考文献）使用的符号和描述。NN 的输入由给出。输入值首先乘以一组权重，以产生加权线性组合，然后使用非线性函数转换，以产生隐藏层神经元的状态值：

在第二层执行类似的操作以产生最终的输出值：

函数通常取为sigmoid函数或。用于多类分类的另一个常见函数是softmax，其定义如下：

这是一个归一化的指数函数。

所有这些都是高度非线性的函数，具有输出值随着输入值急剧增加的性质。这种非线性特性使得神经网络比标准的线性或广义线性模型具有更多的计算灵活性。在这里，被称为偏置参数。权重与偏置一起形成权重向量w。

双层神经网络的示意图结构如下所示：

神经网络中的学习对应于寻找权重向量如w的值，使得对于给定的由真实值输入和目标（输出）组成的训练数据集，网络预测目标值的误差最小。对于回归问题，这是通过最小化误差函数来实现的：

对于分类任务，在神经网络训练中，人们使用交叉熵而不是平方误差，其定义如下：

为了避免过拟合，通常也在目标函数中包含一个正则化项。正则化函数的形式通常是，它对大的w值给予惩罚，从而降低过拟合的可能性。结果的目标函数如下：

在这里，和是自由参数，其最佳值可以通过交叉验证实验找到。

要最小化关于 w 的 M(w)，可以使用 Rumelhart、Hinton 和 Williams 的经典论文中描述的逆传播算法（本章“参考文献”部分的第 3 个参考文献）。在逆传播的每个输入/输出对中，预测输出的值通过从输入层的前向传递来计算。误差，即预测输出和实际输出之间的差异，被反向传播，并且在每个节点，权重都被重新调整，以便误差最小。

贝叶斯神经网络处理

为了将神经网络学习置于贝叶斯框架中，考虑回归案例中的误差函数。它可以被视为观察给定数据集的条件权重 w 的高斯噪声项。这正是可以写成以下形式的似然函数：

这里，是由给出的噪声项的方差，而代表一个概率模型。正则化项可以被认为是参数先验概率分布的对数：

这里，是权重先验分布的方差。可以使用贝叶斯定理轻松证明，目标函数 M(w) 然后对应于参数 w 的后验分布：

在神经网络的情况下，我们感兴趣的是的局部极大值。然后，后验被近似为围绕每个极大值的高斯分布，如下所示：

在这里，A 是 M(w) 关于 w 的二阶导数的矩阵，它代表协方差矩阵的逆。它也被称为海森矩阵。

超参数和的值是通过证据框架找到的。在这里，概率被用作证据，以从数据 D 中找到和的最佳值。这是通过以下贝叶斯规则完成的：

通过使用证据框架和后验的高斯近似（本章“参考文献”部分的第 2 和第 5 条参考文献），可以证明的最佳值满足以下条件：

此外，的最佳值也满足以下条件：

在这些方程中，是给出的确定参数的数量，其中k是w的长度。

brnn R 包

brnn包是由 Paulino Perez Rodriguez 和 Daniel Gianola 开发的，它实现了上一节中描述的两层贝叶斯正则化神经网络。包中的主要函数是brnn( )，可以通过以下命令调用：

>brnn(x,y,neurons,normalize,epochs,…,Monte_Carlo,…)

在这里，x是一个n x p矩阵，其中n是数据点的数量，p是变量的数量；y是一个包含目标值的n维向量。网络隐藏层中的神经元数量可以通过变量neurons指定。如果指示函数normalize为TRUE，则将归一化输入和输出，这是默认选项。模型训练期间的最大迭代次数由epochs指定。如果指示二进制变量Monte_Carlo为真，则使用 MCMC 方法估计 Hessian 矩阵逆的迹。

让我们尝试一个例子，使用我们在第五章中使用的 Auto MPG 数据集，贝叶斯回归模型。下面的 R 代码将导入数据，创建训练集和测试集，使用训练数据训练神经网络模型，并对测试集进行预测：

>install.packages("brnn")  #one time installation
>library(brnn)
>mpgdataall <- read.csv("C:/…/auto-mpg.csv")#give the correct full path
>mpgdata <- mpgdataall[,c(1,3,5,6)]
>#Fitting Bayesian NN Model
>ytrain <- mpgdata[1:100,1]
>xtrain <- as.matrix(mpgdata[1:100,2:4])
>mpg_brnn <- brnn(xtrain,ytrain,neurons=2,normalize = TRUE,epochs = 1000,Monte_Carlo = TRUE)
>summary(mpg_brnn)
A Bayesian regularized neural network
3 - 2 - 1 with 10 weights,biases and connection strengths
Inputs and output were normalized
Training finished because Changes in F= beta*SCE + alpha*Ew in last 3 iterations less than 0.001

>#Prediction using trained model
>ytest <- mpgdata[101:150,1]
>xtest <- as.matrix(mpgdata[101:150,2:4])
>ypred_brnn <- predict.brnn(mpg_brnn,xtest)
>plot(ytest,ypred_brnn)
>err <-ytest-ypred
>summary(err)

深度信念网络和深度学习

在过去十年中，神经网络研究的一些开创性进展为机器学习开辟了一个新的前沿，通常被称为深度学习（本章“参考文献”部分的第 5 和第 7 条参考文献）。深度学习的一般定义是，一类机器学习技术，其中在层次监督架构中利用许多信息处理阶段，用于无监督特征学习和模式分析/分类。深度学习的本质是计算观测数据的层次特征或表示，其中高级特征或因素由低级特征定义（本章“参考文献”部分的第 8 条参考文献）。尽管深度学习有许多类似定义和架构，但所有这些定义中都包含两个共同元素：多层非线性信息处理和在每个层次上从上一层次学习到的特征中进行监督或无监督学习特征表示。深度学习的最初工作基于多层神经网络模型。最近，还使用了许多其他形式的模型，例如深度核机和深度 Q 网络。

即使在几十年前，研究人员也尝试过多层神经网络。然而，有两个原因限制了使用这种架构进行学习的任何进展。第一个原因是网络参数的学习是一个非凸优化问题。从随机初始条件开始，在最小化误差的过程中会陷入局部最小值。第二个原因是相关的计算需求巨大。对于第一个问题的突破发生在 Geoffrey Hinton 开发了一种快速算法来学习一种称为深度信念网络（DBN）的特殊类别的神经网络时。我们将在后面的章节中更详细地描述 DBN。通过使用通用图形处理单元（GPGPU）的计算能力的提升，满足了高计算能力需求。深度学习之所以在实用应用中如此受欢迎，是因为在自动语音识别和计算机视觉中实现了显著的准确性提升。例如，自动语音识别的词错误率在经过多年的研究后，达到了大约 40%的饱和状态。

然而，使用深度学习，词错误率在短短几年内大幅降低，接近 10%。另一个众所周知的例子是深度卷积神经网络在 2012 年 ImageNet 大规模视觉识别挑战赛中的表现，其最低错误率为 15.3%，而当时最先进的方法的最低错误率为 26.2%（本章“参考文献”部分的第 7 条参考文献）。

在本章中，我们将描述一类深度学习模型，称为深度信念网络。感兴趣的读者可能希望阅读李登和董宇所著的书籍（本章“参考文献”部分的第 9 条参考文献），以详细了解深度学习的各种方法和应用。我们将在本章的其余部分遵循他们的符号。我们还将使用 R 包darch来展示 DBN 的使用。

受限玻尔兹曼机

受限玻尔兹曼机（RBM）是一个两层网络（二部图），其中一层是可见层（v），另一层是隐藏层（h）。可见层中的所有节点和隐藏层中的所有节点通过无向边连接，同一层中的节点之间没有连接：

RBM 的特征是所有可见单元和所有隐藏单元的状态的联合分布。

在这里，被称为能量函数，而是名为配分函数的归一化常数，来自统计物理学的术语。

RBM 主要有两种类型。在第一种类型中，v和h都是伯努利随机变量。在第二种类型中，h是伯努利随机变量，而v是高斯随机变量。对于伯努利 RBM，能量函数如下所示：

在这里，代表节点和之间边的权重；和分别是可见层和隐藏层的偏置参数。对于这个能量函数，条件概率的精确表达式可以推导如下：

在这里，是逻辑函数。

如果输入变量是连续的，可以使用高斯 RBM；其能量函数如下所示：

此外，在这种情况下，和的条件概率将变为以下形式：

这是一个均值为和方差为 1 的正态分布。

现在我们已经描述了 RBM 的基本架构，那么它是如何被训练的呢？如果我们尝试使用标准方法，即对对数似然函数求梯度，我们得到以下更新规则：

在这里，是使用数据集计算得到的的期望，而是使用模型计算得到的相同期望。然而，由于难以计算，因此不能使用这个精确的表达式来更新权重。

第一个突破是在解决这个问题的同时，训练深度神经网络时，当 Hinton 及其团队提出了一种名为对比散度（CD）的算法（本章“参考文献”部分的第 7 条）时出现的。算法的精髓将在下一段中描述。

策略是使用 Gibbs 采样从之前提到的条件分布中生成的和的值来近似。实现这一策略的一种方案如下：

从数据集中初始化。
通过从条件分布中采样来找到。
通过从条件分布中采样来找到。
通过从条件分布中采样来找到。

一旦我们找到了和的值，可以使用，即的第i个分量与的第j个分量的乘积，作为的近似。这被称为CD-1 算法。可以将此推广到使用 Gibbs 采样的k步的值，这被称为CD-k 算法。可以很容易地看到 RBM 与贝叶斯推理之间的联系。由于 CD 算法类似于后验密度估计，可以说 RBM 是使用贝叶斯推理方法进行训练的。

虽然对比散度算法看起来很简单，但在训练 RBMs 时需要非常小心，否则模型可能导致过拟合。对在实用应用中使用 RBMs 感兴趣的读者应参考本章参考文献部分的第 10 条参考文献中的技术报告，其中对此进行了详细讨论。

深度信念网络

可以堆叠多个 RBMs，一个叠在另一个上面，使得层中隐藏单元的值成为第n层中可见单元的值，依此类推。所得到的网络称为深度信念网络。它是早期深度学习网络预训练中使用的几种主要架构之一。预训练神经网络的想法如下：在标准的三层（输入-隐藏-输出）神经网络中，可以从权重的随机初始值开始，使用反向传播算法找到对数似然函数的良好最小值。然而，当层数增加时，直接应用反向传播算法不起作用，因为从输出层开始，当我们计算深层层的梯度值时，它们的幅度变得非常小。这被称为梯度消失问题。因此，网络将陷入某些较差的局部最小值。如果我们从良好最小值附近开始，反向传播仍然有效。为了实现这一点，DNN 通常以无监督的方式预训练，使用 DBN。不是从权重的随机值开始，而是以无监督的方式训练 DBN，并使用 DBN 的权重作为相应监督 DNN 的初始权重。观察到使用 DBN 预训练的 DNN 表现更好（本章参考文献部分的第 8 条参考文献）。

DBN 的逐层预训练过程如下。从第一个 RBM 开始，使用可见层中的输入数据和使用 CD 算法（或其最新的更好变体）对其进行训练。然后，在这个 RBM 的上面堆叠第二个 RBM。对于这个 RBM，使用中采样的值作为可见层的值。继续这个过程，直到达到所需的层数。顶层隐藏单元的输出也可以用作训练监督模型的输入。为此，在 DBN 的顶部添加一个具有所需类别数作为输出节点数的传统神经网络层。这个 NN 的输入将是 DBN 顶层的输出。这被称为DBN-DNN 架构。在这里，DBN 的作用是从输入数据自动生成高度有效的特征（DBN 顶层输出）供顶层监督 NN 使用。

用于二分类任务的五层 DBN-DNN 架构如图所示：

最后一个层使用监督方式使用反向传播算法对两个类别和进行训练。我们将使用 darch R 包通过这样的 DBN-DNN 来展示训练和分类。

darch R 包

由 Martin Drees 编写的 darch 包是 R 包之一，使用它可以开始使用 R 进行深度学习。它实现了上一节中描述的 DBN（参考文献 5 和 7 在本章“参考文献”部分的引用）。该包可以从cran.r-project.org/web/packages/darch/index.html下载。

darch 包中的主要类实现了深度架构，并提供了使用对比散度训练以及使用反向传播、弹性反向传播和共轭梯度微调进行微调的能力。该类的新的实例是通过newDArch构造函数创建的。它使用以下参数调用：一个包含每层节点数的向量，批量大小，一个布尔变量，用于指示是否使用ff包来计算权重和输出，以及生成权重矩阵的函数名称。让我们创建一个具有两个输入单元、四个隐藏单元和一个输出单元的网络：

install.packages("darch") #one time
>library(darch)
>darch <- newDArch(c(2,4,1),batchSize = 2,genWeightFunc = generateWeights)
INFO [2015-07-19 18:50:29] Constructing a darch with 3 layers.
INFO [2015-07-19 18:50:29] Generating RBMs.
INFO [2015-07-19 18:50:29] Construct new RBM instance with 2 visible and 4 hidden units.
INFO [2015-07-19 18:50:29] Construct new RBM instance with 4 visible and 1 hidden units.

让我们使用一个玩具数据集来训练 DBN。我们之所以这样做，是因为训练任何真实示例将花费很长时间：如果不是几天，可能就是几个小时。让我们创建一个包含两列和四行的输入数据集：

>inputs <- matrix(c(0,0,0,1,1,0,1,1),ncol=2,byrow=TRUE)
>outputs <- matrix(c(0,1,1,0),nrow=4)

现在，让我们使用输入数据对 DBN 进行预训练：

>darch <- preTrainDArch(darch,inputs,maxEpoch=1000)

我们可以使用getLayerWeights( )函数查看任何层的权重。让我们看看隐藏层看起来如何：

>getLayerWeights(darch,index=1)
[[1]]
          [,1]        [,2]       [,3]       [,4]
[1,]   8.167022    0.4874743  -7.563470  -6.951426
[2,]   2.024671  -10.7012389   1.313231   1.070006
[3,]  -5.391781    5.5878931   3.254914   3.000914

现在，让我们为监督学习进行一次反向传播。为此，我们首先需要将层函数设置为sigmoidUnitDerivatives：

>layers <- getLayers(darch)
>for(i in length(layers):1){
     layers[[i]][[2]] <- sigmoidUnitDerivative
    }
>setLayers(darch) <- layers
>rm(layers)

最后，以下两行执行反向传播：

>setFineTuneFunction(darch) <- backpropagation
>darch <- fineTuneDArch(darch,inputs,outputs,maxEpoch=1000)

我们可以通过运行darch来查看 DBN 在训练数据本身上的预测质量，如下所示：

>darch <- getExecuteFunction(darch)(darch,inputs)
>outputs_darch <- getExecOutputs(darch)
>outputs_darch[[2]]
        [,1]
[1,] 9.998474e-01
[2,] 4.921130e-05
[3,] 9.997649e-01
[4,] 3.796699e-05

与实际输出相比，DBN 对第一行和第二行输入的输出预测是错误的。由于这个例子只是为了说明如何使用 darch 包，所以我们在这里并不担心 50%的准确率。

R 中的其他深度学习包

尽管 R 语言中存在其他深度学习包，例如deepnet和RcppDL，但与 C++的Cuda和 Python 的Theano等语言的库相比，R 语言在深度学习方面还没有好的本地库。唯一可用的包是 Java 基础的深度学习开源项目 H2O 的包装器。这个 R 包，h2o，允许在 R 中通过其 REST API 运行 H2O。对深度学习项目和应用程序感兴趣的用户应使用 R 中的 h2o 包来使用 H2O。需要在你机器上安装 H2O 才能使用 h2o。我们将在讨论大数据和名为 Spark 的分布式计算平台时，在下一章中介绍 H2O。

练习

对于 Auto MPG 数据集，比较使用普通回归、贝叶斯 GLM 和贝叶斯神经网络进行预测模型的性能。

参考文献

MacKay D. J. C. 信息论、推理与学习算法. 剑桥大学出版社. 2003. ISBN-10: 0521642981
MacKay D. J. C. "应用于分类网络的证据框架". 神经计算. 第 4 卷，第 3 期，698-714. 1992
MacKay D. J. C. "可能的网络与合理的预测——监督神经网络实用贝叶斯方法的综述". 网络：神经系统的计算
Hinton G. E., Rumelhart D. E., 和 Williams R. J. "通过反向传播错误学习表示". 自然. 第 323 卷，533-536. 1986
MacKay D. J. C. "贝叶斯插值". 神经计算. 第 4 卷，第 3 期，415-447. 1992
Hinton G. E., Krizhevsky A., 和 Sutskever I. "使用深度卷积神经网络进行 ImageNet 分类". 神经信息处理系统（NIPS）的进展. 2012
Hinton G., Osindero S., 和 Teh Y. "深度信念网的快速学习算法". 神经计算. 第 18 卷:1527–1554. 2006
Hinton G. 和 Salakhutdinov R. "使用神经网络降低数据维度". 科学. 第 313 卷，第 5786 期:504–507. 2006
李登辉和杜东宇. 深度学习：方法与应用（信号处理的基础与趋势）. Now Publishers Inc. 第 7 卷，第 3-4 期. 2014. ISBN-13: 978-1601988140
Hinton G. "训练受限玻尔兹曼机的实用指南". UTML 技术报告 2010-003. 多伦多大学. 2010

摘要

在本章中，我们学习了一个重要的机器学习模型类别，即神经网络，以及它们的贝叶斯实现。这些模型受到人脑架构的启发，并且继续是活跃的研究和开发领域。我们还学习了一种最新的神经网络进展，称为深度学习。它可以用于解决许多问题，例如涉及高度认知元素的计算机视觉和自然语言处理。使用深度学习的人工智能系统能够在语音识别和图像分类等任务中达到与人类智能相当的高精度。通过本章，我们已经涵盖了贝叶斯机器学习模型的重要类别。在下一章中，我们将探讨不同的方面：大规模机器学习及其在贝叶斯模型中的某些应用。

第九章。大数据规模下的贝叶斯建模

当我们在第三章中学习了贝叶斯推理原理，即《介绍贝叶斯推理》时，我们发现随着训练数据量的增加，数据对参数估计的贡献超过了先验分布的贡献。此外，参数估计的不确定性也会降低。因此，你可能想知道为什么在大规模数据分析中需要贝叶斯建模。为了回答这个问题，让我们看看这样一个问题，即为电子商务产品构建推荐系统。

在典型的电子商务商店中，可能会有数百万用户和数万种产品。然而，每个用户在其一生中可能只购买过商店中所有产品的一小部分（不到 10%）。让我们假设电子商务商店正在收集每个销售产品的用户反馈，以 1 到 5 分的评分尺度进行。然后，商店可以创建一个用户-产品评分矩阵来捕捉所有用户的评分。在这个矩阵中，行对应于用户，列对应于产品。每个单元格的值将是用户（对应于行）对产品（对应于列）给出的评分。现在，很容易看出，尽管这个矩阵的整体大小很大，但只有不到 10%的条目会有值，因为每个用户只从商店购买了不到 10%的产品。所以，这是一个高度稀疏的数据集。每当存在一个机器学习任务，尽管整体数据量很大，但数据高度稀疏时，可能会发生过拟合，并且应该依赖于贝叶斯方法（参考本章“参考文献”部分的第 1 条）。此外，许多模型，如贝叶斯网络、潜在狄利克雷分配和深度信念网络，都是建立在贝叶斯推理范式之上的。

当这些模型在大数据集上训练，例如来自路透社的文本语料库时，那么潜在的问题是大规模贝叶斯建模。实际上，贝叶斯建模是计算密集型的，因为我们必须估计参数的整个后验分布，并且还要对预测进行模型平均。大数据集的存在将使情况变得更糟。那么，我们能够使用哪些计算框架在 R 中进行大规模贝叶斯学习呢？在接下来的两个部分中，我们将讨论这个领域的一些最新发展。

使用 Hadoop 进行分布式计算

在过去十年中，当两位来自 Google 的研究工程师开发了一种名为 MapReduce 框架的计算范式以及一个相关的分布式文件系统 Google 文件系统（本章 参考文献 部分的参考 2）时，分布式计算取得了巨大的进步。后来，Yahoo 开发了一个名为 Hadoop 的开源分布式文件系统版本，这成为了大数据计算的一个标志。Hadoop 通过将数据分布到多台计算机并在每个节点上从磁盘本地执行计算，非常适合处理无法适应单个大型计算机内存的大量数据。一个例子是从日志文件中提取相关信息，通常一个月的数据量在千兆字节级别。

要使用 Hadoop，必须使用 MapReduce 框架编写程序以并行化计算。Map 操作将数据分割成多个键值对，并发送到不同的节点。在每个节点上，对每个键值对进行计算。然后，有一个洗牌操作，将具有相同键值的所有键值对聚集在一起。之后，Reduce 操作将前一步计算中对应相同键的所有结果求和。通常，这些 MapReduce 操作可以使用称为 Pig 的高级语言编写。也可以使用 RHadoop 包在 R 中编写 MapReduce 程序，我们将在下一节中描述。

从 R 使用 RHadoop 操作 Hadoop

RHadoop 是一系列开源包的集合，R 用户可以使用这些包管理和分析存储在 Hadoop 分布式文件系统（HDFS）中的数据。在后台，RHadoop 将这些操作转换为 Java 中的 MapReduce 操作并在 HDFS 上运行。

RHadoop 中的各种包及其用途如下：

rhdfs：使用此包，用户可以从 R 连接到 HDFS 并执行基本操作，如读取、写入和修改文件。
rhbase：这是连接到 HBASE 数据库并读取、写入和修改表的包。
plyrmr：使用此包，R 用户可以执行常见的数据操作任务，如数据集的切片和切块。这与 plyr 或 reshape2 等包的功能类似。
rmr2：使用此包，用户可以在 R 中编写 MapReduce 函数并在 HDFS 中执行它们。

与本书中讨论的其他包不同，与 RHadoop 相关的包不可从 CRAN 获取。它们可以从 GitHub 仓库 github.com/RevolutionAnalytics 下载，并从本地驱动器安装。

下面是一个使用 rmr2 包编写的 MapReduce 代码示例，用于统计语料库中的单词数量（本章 参考文献 部分的参考 3）：

第一步是加载 rmr2 库：

>library(rmr2)
>LOCAL <- T #to execute rmr2 locally

第二步涉及编写 Map 函数。这个函数将文本文档中的每一行分割成单词。每个单词被视为一个标记。该函数输出键值对，其中每个不同的单词是一个键，值 = 1：
```
>#map function
>map.wc <- function(k,lines){ 
       words.list <- strsplit(lines,'\\s+^' )
        words <- unlist(words.list)         
        return(keyval(words,1))
    }
```
第三步涉及编写 Reduce 函数。这个函数将来自不同 Mapper 的相同键进行分组并求和它们的值。由于在这种情况下，每个单词都是一个键，且 值 = 1，Reduce 的输出将是单词的计数：
```
>#reduce function
>reduce.wc<-function(word,counts){
          return(keyval(word,sum(counts) ))
}
```

第四步涉及编写一个结合 Map 和 Reduce 函数的词频统计函数，并在名为 hdfs.data 的文件上执行此函数，该文件存储在包含输入文本的 HDFS 中：

>#word count function
>wordcount<-function(input,output=NULL){
            mapreduce(input = input,output = output,input.format = "text",map = map.wc,reduce = reduce.wc,combine = T)
}
>out<-wordcount(hdfs.data,hdfs.out)

第五步涉及从 HDFS 获取输出文件并打印前五行：

>results<-from.dfs(out)
>results.df<-as.data.frame(results,stringAsFactors=F)
>colnames(results.df)<-c('word^' ,^' count^')
>head(results.df)

Spark – 内存分布式计算

Hadoop 的问题之一是在 MapReduce 操作之后，结果文件被写入硬盘。因此，当进行大量数据处理操作时，硬盘上会有许多读写操作，这使得 Hadoop 的处理速度非常慢。此外，网络延迟，即在不同节点之间洗牌数据所需的时间，也加剧了这个问题。另一个缺点是，无法从存储在 HDFS 中的文件中进行实时查询。对于机器学习问题，在训练阶段，MapReduce 不会在迭代中持久化。所有这些都使得 Hadoop 不是一个理想的机器学习平台。

2009 年，加州大学伯克利分校的 AMP 实验室发明了这种解决方案。这是罗马尼亚出生的计算机科学家 Matei Zaharia 博士的博士论文 Resilient Distributed Datasets: A Fault-Tolerant Abstraction for In-Memory Cluster Computing（本章“参考文献”部分的第 4 条参考文献）的成果，该论文催生了 Spark 项目，最终成为 Apache 下的一个完全开源项目。Spark 是一个内存分布式计算框架，解决了之前提到的许多 Hadoop 问题。此外，它支持比 MapReduce 更多的操作类型。Spark 可以用于处理迭代算法、交互式数据挖掘和流式应用。它基于一个称为 Resilient Distributed Datasets（RDD）的抽象。类似于 HDFS，它也是容错的。

Spark 是用一种名为 Scala 的语言编写的。它具有从 Java 和 Python 使用的接口，从最近的 1.4.0 版本开始；它还支持 R。这被称为 SparkR，我们将在下一节中描述。Spark 中可用的四个库类别是 SQL 和 DataFrames、Spark Streaming、MLib（机器学习）和 GraphX（图算法）。目前，SparkR 仅支持 SQL 和 DataFrames；其他肯定在路线图中。Spark 可以从 Apache 项目页面spark.apache.org/downloads.html下载。从 1.4.0 版本开始，SparkR 包含在 Spark 中，无需单独下载。

SparkR

与 RHadoop 类似，SparkR 是一个 R 包，允许 R 用户通过RDD类使用 Spark API。例如，使用 SparkR，用户可以从 RStudio 上运行 Spark 作业。SparkR 可以从 RStudio 中调用。为了启用此功能，请在 R 启动时初始化环境的.Rprofile文件中包含以下行：

Sys.setenv(SPARK_HOME/.../spark-1.5.0-bin-hadoop2.6")
#provide the correct path where spark downloaded folder is kept for SPARK_HOME 
.libPaths(c(file.path(Sys.getenv("SPARK_HOME"),""R",""lib"),".libPaths()))

完成这些后，启动 RStudio 并输入以下命令以开始使用 SparkR：

>library(SparkR)
>sc <- sparkR.init(master="local")

如前所述，当本章撰写时，SparkR 支持 R 的有限功能。这主要包括数据切片和切块以及汇总统计函数。当前版本不支持使用贡献的 R 包；然而，计划在未来的版本中实现。在机器学习方面，当前 SparkR 支持glm()函数。我们将在下一节中做一个示例。

使用 SparkR 进行线性回归

在以下示例中，我们将说明如何使用 SparkR 进行机器学习。为此，我们将使用与第五章中线性回归相同的能源效率测量数据集，贝叶斯回归模型：

>library(SparkR)
>sc <- sparkR.init(master="local")
>sqlContext <- sparkRSQL.init(sc)

#Importing data
>df <- read.csv("/Users/harikoduvely/Projects/Book/Data/ENB2012_data.csv",header = T)
>#Excluding variable Y2,X6,X8 and removing records from 768 containing mainly null values
>df <- df[1:768,c(1,2,3,4,5,7,9)]
>#Converting to a Spark R Dataframe
>dfsr <- createDataFrame(sqlContext,df) 
>model <- glm(Y1 ~ X1 + X2 + X3 + X4 + X5 + X7,data = dfsr,family = "gaussian")
 > summary(model)

云计算集群

为了使用 Hadoop 和相关 R 包处理大型数据集，需要一个计算机集群。在当今世界，使用亚马逊、微软和其他人提供的云计算服务很容易。只需支付使用的 CPU 和存储量。无需在基础设施上进行前期投资。前四大云计算服务是亚马逊的 AWS、微软的 Azure、谷歌的 Compute Cloud 和 IBM 的 Bluemix。在本节中，我们将讨论在 AWS 上运行 R 程序。特别是，你将学习如何创建 AWS 实例；在该实例中安装 R、RStudio 和其他包；开发和运行机器学习模型。

亚马逊网络服务

广为人知的 AWS，即亚马逊网络服务，始于 2002 年亚马逊的一个内部项目，旨在满足支持其电子商务业务的动态计算需求。这作为一个基础设施即服务项目发展起来，2006 年亚马逊向世界推出了两项服务，简单存储服务（S3）和弹性计算云（EC2）。从那时起，AWS 以惊人的速度增长。如今，他们拥有超过 40 种不同类型的服务，使用数百万台服务器。

在 AWS 上创建和运行计算实例

学习如何设置 AWS 账户并开始使用 EC2 的最佳地方是亚马逊 Kindle 商店中免费提供的电子书，名为 Amazon Elastic Compute Cloud (EC2) 用户指南（参考本章节 参考文献 部分的第 6 条）。

这里，我们仅总结涉及此过程的基本步骤：

创建 AWS 账户。
登录 AWS 管理控制台 (aws.amazon.com/console/)。
点击 EC2 服务。
选择 Amazon Machine Instance (AMI)。
选择实例类型。
创建公私钥对。
配置实例。
添加存储。
标记实例。
配置安全组（指定谁可以访问实例的策略）。
审查并启动实例。

使用 SSH（从 Linux/Ubuntu）、Putty（从 Windows）或使用配置安全时提供的私钥和启动时给出的 IP 地址通过浏览器登录到您的实例。这里，我们假设您启动的实例是一个 Linux 实例。

安装 R 和 RStudio

要安装 R 和 RStudio，您需要成为认证用户。因此，创建一个新用户并授予用户管理员权限（sudo）。之后，从 Ubuntu shell 执行以下步骤：

编辑 /etc/apt/sources.list 文件。

在末尾添加以下行：

deb http://cran.rstudio.com/bin/linux/ubuntu trusty .

获取运行存储库的密钥：

sudo apt-key adv  --keyserver keyserver.ubuntu.com –recv-keys 51716619E084DAB9

更新包列表：
```
sudo apt-get update
```
安装 R 的最新版本：
```
sudo apt-get install r-base-core
```
安装 gdebi 以从本地磁盘安装 Debian 软件包：
```
sudo apt-get install gdebi-core
```

下载 RStudio 软件包：

wget http://download2.rstudio.org/r-studio-server-0.99.446-amd64.deb

安装 RStudio：

sudo gdebi r-studio-server-0.99.446-amd64.deb

安装成功完成后，运行在您的 AWS 实例上的 RStudio 可以通过浏览器访问。为此，请打开浏览器并输入 URL <your.aws.ip.no>:8787。

如果你能够使用运行在 AWS 实例上的 RStudio，那么你可以从 RStudio 安装其他包，例如 rhdfs、rmr2 等，并在 R 中构建任何机器学习模型，然后在 AWS 云上运行它们。

除了 R 和 RStudio，AWS 还支持 Spark（因此也支持 SparkR）。在下一节中，您将学习如何在 EC2 集群上运行 Spark。

在 EC2 上运行 Spark

您可以使用位于您本地机器 Spark 的 ec2 目录中的 spark-ec2 脚本来在 Amazon EC2 上启动和管理 Spark 集群。要在 EC2 上启动 Spark 集群，请按照以下步骤操作：

在您的本地机器 Spark 文件夹中的 ec2 目录下。
运行以下命令：
```
./spark-ec2 -k <keypair> -i <key-file> -s <num-slaves> launch <cluster-name>
```
在这里，<keypair> 是您用于启动本章“在 AWS 上创建和运行计算实例”部分中提到的 EC2 服务的密钥对名称。<key-file> 是您本地机器上私钥已下载并保存的路径。工作节点数由 <num-slaves> 指定。
要在集群上运行您的程序，请首先使用以下命令通过 SSH 连接到集群：
```
./spark-ec2 -k <keypair> -i <key-file> login <cluster-name>
```
登录到集群后，您可以使用与在本地机器上使用相同的方式使用 Spark。

在 Spark 文档和 AWS 文档中可以找到有关如何在 EC2 上使用 Spark 的更多详细信息（章节“参考文献”部分的第 5、6 和 7 条）。

Microsoft Azure

Microsoft Azure 对 R 和 Spark 提供了全面支持。Microsoft 收购了 Revolution Analytics 公司，该公司开始构建并支持 R 的企业版。除此之外，Azure 还有一个机器学习服务，其中包含一些贝叶斯机器学习模型的 API。有关如何在 Azure 上启动实例以及如何使用其机器学习服务的视频教程可以在 Microsoft Virtual Academy 网站上找到（章节“参考文献”部分的第 8 条）。

IBM Bluemix

Bluemix 通过其实例上可用的完整 R 库集对 R 完全支持。IBM 在其路线图计划中也包含了将 Spark 集成到其云服务中。更多详细信息可以在他们的文档页面上找到（章节“参考文献”部分的第 9 条）。

其他用于大规模机器学习的 R 包

除了 RHadoop 和 SparkR 之外，还有几个专门为大规模机器学习构建的本地 R 包。在这里，我们简要概述它们。感兴趣的读者应参考 CRAN 任务视图：使用 R 进行高性能和并行计算（章节“参考文献”部分的第 10 条）。

虽然 R 是单线程的，但存在几个用于 R 的并行计算包。其中一些知名的包是 Rmpi（流行的消息传递接口的 R 版本）、multicore、snow（用于构建 R 集群）和 foreach。从 R 2.14.0 版本开始，一个新的名为 parallel 的包开始随基础 R 一起发货。我们将在下面讨论其一些功能。

并行 R 包

并行包建立在多核和 snow 包之上。它适用于在多个数据集上运行单个程序，例如 K 折交叉验证。它可以用于在单个机器上的多个 CPU/核心或跨多台机器进行并行化。对于跨机器集群的并行化，它通过 Rmpi 包调用 MPI（消息传递接口）。

我们将通过计算列表 1:100000 的数字平方的简单示例来说明并行包的使用。此示例在 Windows 上无法工作，因为相应的 R 不支持多核包。它可以在任何 Linux 或 OS X 平台上进行测试。

执行此操作的顺序方法是使用lapply函数，如下所示：

>nsquare <- function(n){return(n*n)}
>range <- c(1:100000)
>system.time(lapply(range,nsquare))

使用并行包中的mclapply函数，可以在更短的时间内完成这个计算：

>library(parallel) #included in R core packages, no separate installation required
>numCores<-detectCores( )  #to find the number of cores in the machine
>system.time(mclapply(range,nsquare,mc.cores=numCores))

如果数据集非常大，需要使用计算机集群，我们可以使用parLapply函数在集群上运行程序。这需要 Rmpi 包：

>install.packages(Rmpi)#one time
>library(Rmpi)
>numNodes<-4 #number of workers nodes
>cl<-makeCluster(numNodes,type="MPI")
>system.time(parLapply(cl,range,nsquare))
>stopCluster(cl)
>mpi.exit( )

The foreach R package

这是在 R 中的一种新的循环结构，可以在多核或集群上并行执行。它有两个重要的运算符：%do%用于重复执行任务，%dopar%用于并行执行任务。

例如，我们可以在前一个章节中讨论的平方函数中使用 foreach 包的单行命令来实现：

>install.packages(foreach)#one time
>install.packages(doParallel)#one time
>library(foreach)
>library(doParallel)
>system.time(foreach(i=1:100000)   %do%  i²) #for executing sequentially
>system.time(foreach(i=1:100000)   %dopar%  i²) #for executing in parallel

我们还将使用foreach函数的例子来演示快速排序：

>qsort<- function(x) {
  n <- length(x)
  if (n == 0) {
    x
  } else {
    p <- sample(n,1)
    smaller <- foreach(y=x[-p],.combine=c) %:% when(y <= x[p]) %do% y
    larger  <- foreach(y=x[-p],.combine=c) %:% when(y >  x[p]) %do% y
    c(qsort(smaller),x[p],qsort(larger))
  }
}
qsort(runif(12))

这些包仍在进行大量开发。它们还没有被广泛用于贝叶斯建模。它们很容易用于贝叶斯推断应用，如蒙特卡洛模拟。

练习

回顾第六章中的分类问题，贝叶斯分类模型。使用 SparkR 的glm()函数重复相同的问题。
使用 SparkR 回顾本章中我们做的线性回归问题。在创建 AWS 实例后，使用 AWS 上的 RStudio 服务器重复这个问题。

参考文献

"变分贝叶斯概率矩阵分解算法的 MapReduce 实现"。在：IEEE 大数据会议。第 145-152 页。2013 年
Dean J. and Ghemawat S. "MapReduce: Simplified Data Processing on Large Clusters". Communications of the ACM 51 (1). 107-113
github.com/jeffreybreen/tutorial-rmr2-airline/blob/master/R/1-wordcount.R
Chowdhury M., Das T., Dave A., Franklin M.J., Ma J., McCauley M., Shenker S., Stoica I.和 Zaharia M. "弹性分布式数据集：内存集群计算的容错抽象"。NSDI 2012。2012 年
亚马逊弹性计算云（EC2）用户指南，亚马逊网络服务 Kindle 电子书，更新于 2014 年 4 月 9 日
Spark 文档在 AWS 上的说明，请参阅spark.apache.org/docs/latest/ec2-scripts.html
AWS 上的 Spark 文档，请参阅aws.amazon.com/elasticmapreduce/details/spark/
微软虚拟学院网站，请参阅www.microsoftvirtualacademy.com/training-courses/getting-started-with-microsoft-azure-machine-learning
IBM Bluemix 教程，请参阅www.ibm.com/developerworks/cloud/bluemix/quick-start-bluemix.html
CRAN Task View for contributed packages in R at cran.r-project.org/web/views/HighPerformanceComputing.html

摘要

在本书的最后一章中，我们介绍了各种实现大规模机器学习的框架。这些框架对贝叶斯学习也非常有用。例如，为了从后验分布中进行模拟，可以在机器集群上运行 Gibbs 抽样。我们学习了如何使用 RHadoop 包从 R 连接到 Hadoop，以及如何使用 SparkR 与 Spark 一起使用 R。我们还讨论了如何在 AWS 等云服务中设置集群，以及如何在它们上运行 Spark。还涵盖了某些本地并行化框架，如 parallel 和 foreach 函数。

本书的主要目的是向读者介绍使用 R 进行贝叶斯建模的领域。读者应该已经对贝叶斯机器学习模型背后的理论和概念有了很好的理解。由于示例主要是为了说明目的，我敦促读者将这些技术应用于实际问题，以更深入地理解贝叶斯推断的主题。

posted @ 2025-09-03 10:22 绝不原创的飞龙阅读(43) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟

掠夺·扩张·投机·博弈

R-贝叶斯模型学习指南-全-

R 贝叶斯模型学习指南（全）

前言

本书涵盖的内容

你需要这本书的什么

这本书面向的对象

术语约定

注意

提示

读者反馈

客户支持

下载示例代码

勘误

盗版

问题

第一章：介绍概率论

概率分布

条件概率

贝叶斯定理

边际分布

期望和协方差

二项分布

Beta 分布

伽玛分布

狄利克雷分布

威沙特分布

练习

参考文献

摘要

第二章：R 环境

设置 R 环境和包

安装 R 和 RStudio

你的第一个 R 程序

R 中的数据管理

R 中的数据类型

R 中的数据结构

将数据导入 R

切割和切块数据集

向量化操作

编写 R 程序

控制结构

函数

作用域规则

循环函数

lapply

sapply

mapply

apply

tapply

数据可视化

高级绘图函数

低级绘图命令

交互式图形函数

抽样

从区间进行随机均匀抽样

从正态分布中进行抽样

练习

参考文献

小贴士

摘要

第三章。介绍贝叶斯推理

贝叶斯不确定性视图

选择合适的先验分布

非信息先验

主观先验

共轭先验

层次先验

后验分布的估计

最大后验概率估计

拉普拉斯近似

蒙特卡洛模拟

Metropolis-Hasting 算法

Metropolis-Hasting 算法的 R 包

吉布斯采样

吉布斯采样的 R 包

变分近似

未来观测的预测

练习