{
"cells": [
{
"cell_type": "markdown",
"id": "eb072a11",
"metadata": {},
"source": [
"# 第 2 章 概率编程\n",
"\n",
"\n",
"\n",
"上一章对贝叶斯统计有了基本了解，本章将学习如何使用计算工具构建概率模型。我们将学习使用 `PyMC3` 进行概率编程。其基本思想是使用代码指定模型，然后以或多或少自动的方式求解它们。\n",
"\n",
"选择概率编程的背后原因是：许多模型无法得到封闭形式解，因此只能使用依托于计算机的数值方法来求解。学习概率编程的另一个原因是，现代贝叶斯统计主要是通过编写代码来完成的，既然已经了解 `Python` ，为什么还要用另一种方式呢？概率编程提供了一种构建和求解复杂模型的有效方法，使我们可以更多地关注模型设计、评估和解释，而更少地关注数学或计算细节。\n",
"\n",
"在本章以及本书的其余部分中，我们将使用 `PyMC3` 和 `ArviZ` ， `PyMC3` 是一个非常灵活的概率编程 `Python` 库，`ArviZ` 是一个新的 `Python` 库，它将帮助我们解释概率模型的结果。了解 `PyMC3` 和 `ArviZ` 还将帮助我们以更实际的方式学习先进的贝叶斯概念。\n",
"\n",
"本章涵盖以下主题：\n",
"\n",
"- 概率编程\n",
"- 推断引擎\n",
"- `PyMC3` 指南\n",
"- 重温抛硬币问题\n",
"- 模型检查和诊断\n",
"- 高斯模型和学生 $\\boldsymbol{t}$ 模型\n",
"- 分组比较和有效容量\n",
"- 分层模型和收缩\n",
" \n",
"---\n",
"\n",
"## 2.1 概率编程\n",
"\n",
"贝叶斯统计在概念上非常简单：我们有已知量和未知量，通过贝叶斯定理以前者为条件推断后者。如果幸运的话，该过程将减少未知量的不确定性。通常，我们把已知量称为**数据**并将其视为确切值，将未知量称为参数并将其视为概率分布。\n",
"\n",
"在更正式的术语中，我们为未知量分配先验概率分布；然后利用贝叶斯定理将其先验分布转化为后验分布。尽管概念简单，但贝叶斯定理中，分母项的全概率公式常常无法得到封闭形式的解析表达式，使得获得真实后验的难度很大。多年来，该问题一直是阻碍贝叶斯方法广泛采用的主要原因之一。\n",
"\n",
"随着计算时代的到来，数值方法使得解决任意后验分布的推断问题成为可能。这极大改变了贝叶斯数据分析的应用。我们可以把这些数值方法看作**通用推断引擎**，或者 `PyMC3` 核心开发者 `Thomas Wiecki` 所称的**推断按钮**。自动化推断过程的可能性导致了 **概率编程语言（PPL）** 的发展，它导致了模型创建和推断之间的解耦，从而让用户专注于模型构建和优化，而不用过多考虑如何得到后验分布的问题。\n",
"\n",
"在概率编程语言的框架中，用户只需要寥寥数行代码描述概率模型，后面的推断过程就能自动完成了。概率编程使得人们能够更快速地构建复杂的概率模型并减少出错的可能。可以预见，这将给数据科学和其他学科带来极大的影响。\n",
"\n",
"我认为，概率编程语言对科学计算的影响可以与 60 多年前 `Fortran` 语言的问世相对比。尽管如今 `Fortran` 语言风光不再，不过在当时 `Fortran` 语言被认为是相当革命性的。当时的科学家们第一次从计算细节中解放出来，开始用一种更自然的方式将注意力放在构建数值化方法、模型和仿真系统上。类似地，概率编程将处理概率和推断的过程对用户隐藏起来，从而使得用户更多去关注模型构建和结果分析。\n",
"\n",
"在本章中，我们将学习如何使用 `PyMC3` 定义和求解模型。\n",
"\n",
"我们将把贝叶斯推断看作一个黑盒，只需要运行一次，它就能够为我们提供了抽取自后验的适当样本。该方法是随机的，所以每次运行黑盒子时，得到的样本都会有所不同。不过理论上，如果推断过程能够按照预期工作，即使样本不同，我们也应当能够从其中任何一个样本中得到相同的结论。\n",
"\n",
"至于当按下 “推断驱动按钮” 时，引擎中会发生什么，以及如何检查样本是否确实值得信任等问题，将在第 8 章中解释。\n",
"\n",
"## 2.2 用 `PyMC3` 做后验推断\n",
"\n",
" `PyMC3` 是一个 `Python` 库，用于概率编程。撰写本文时的最后一个版本是 3.6 。 `PyMC3` 提供了非常简单直观的语法，易于阅读，与统计文献中用于描述概率模型的语法非常接近。 `PyMC3` 的基本代码是用 `Python` 编写的，计算要求高的部分是用 `Numpy` 和 `Theano` 编写的。\n",
"\n",
"`Theano` 是为深度学习而开发的一个 `Python` 库，允许我们高效地定义、优化和计算涉及多维数组的数学表达式。 `PyMC3` 使用 `Theano` 的主要原因是，有些采样方法需要计算梯度，而 `Theano` 知道如何使用自动微分来计算梯度。此外，`Theano` 将 `Python` 代码编译成 C 代码，因此 `PyMC3` 非常快。这是关于 `Theano` 的所有信息，我们必须使用 `PyMC3` 。\n",
"\n",
"如果您还想了解更多，请阅读[官方 `Theano` 教程](http://deeplearning.net/software/theano/tutorial/index.html#tutorial)。\n",
"\n",
"> 你可能听说过 Theano 已经不再开发了，但这没什么好担心的。PyMC 开发人员将接管 Theano 维护，确保 Theano 在未来几年内继续为 `PyMC3` 服务。与此同时，PyMC 开发人员正在迅速行动，以创建 `PyMC3` 的继任者。这可能会基于 TensorFlow 作为后端，尽管其他选项也在分析中。有关这方面的更多信息，请访问[博客](https://medium.com/@pymc_devs/theano-tensorflow-and-the-future-of-pymc-6c9987bb19d5) 。\n",
"\n",
"重新回顾抛硬币问题，这次使用 `PyMC3` 。首先获取数据，此处使用手动构造的数据。由于数据是我们通过代码生成的，所以知道真实参数值，以下代码中用 `theta_real` 变量表示。显然，在实际数据中，通常并不知道参数的真实值，而是要将其估计出来。"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "27867052",
"metadata": {},
"outputs": [
{
"ename": "AttributeError",
"evalue": "module 'numpy.distutils.__config__' has no attribute 'blas_opt_info'",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNoSectionError\u001b[0m Traceback (most recent call last)",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/theano/configparser.py:238\u001b[0m, in \u001b[0;36mTheanoConfigParser.fetch_val_for_key\u001b[0;34m(self, key, delete_key)\u001b[0m\n\u001b[1;32m 237\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m--> 238\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_theano_cfg\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[43msection\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43moption\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 239\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m InterpolationError:\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/configparser.py:781\u001b[0m, in \u001b[0;36mRawConfigParser.get\u001b[0;34m(self, section, option, raw, vars, fallback)\u001b[0m\n\u001b[1;32m 780\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m--> 781\u001b[0m d \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_unify_values\u001b[49m\u001b[43m(\u001b[49m\u001b[43msection\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mvars\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 782\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m NoSectionError:\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/configparser.py:1149\u001b[0m, in \u001b[0;36mRawConfigParser._unify_values\u001b[0;34m(self, section, vars)\u001b[0m\n\u001b[1;32m 1148\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m section \u001b[38;5;241m!=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdefault_section:\n\u001b[0;32m-> 1149\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m NoSectionError(section) \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;28mNone\u001b[39m\n\u001b[1;32m 1150\u001b[0m \u001b[38;5;66;03m# Update with the entry specific variables\u001b[39;00m\n",
"\u001b[0;31mNoSectionError\u001b[0m: No section: 'blas'",
"\nDuring handling of the above exception, another exception occurred:\n",
"\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/theano/configparser.py:354\u001b[0m, in \u001b[0;36mConfigParam.__get__\u001b[0;34m(self, cls, type_, delete_key)\u001b[0m\n\u001b[1;32m 353\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m--> 354\u001b[0m val_str \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mcls\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfetch_val_for_key\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mname\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdelete_key\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdelete_key\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 355\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mis_default \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mFalse\u001b[39;00m\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/theano/configparser.py:242\u001b[0m, in \u001b[0;36mTheanoConfigParser.fetch_val_for_key\u001b[0;34m(self, key, delete_key)\u001b[0m\n\u001b[1;32m 241\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m (NoOptionError, NoSectionError):\n\u001b[0;32m--> 242\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m(key)\n",
"\u001b[0;31mKeyError\u001b[0m: 'blas__ldflags'",
"\nDuring handling of the above exception, another exception occurred:\n",
"\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)",
"Input \u001b[0;32mIn [1]\u001b[0m, in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mpandas\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mpd\u001b[39;00m\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mseaborn\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01msns\u001b[39;00m\n\u001b[0;32m----> 6\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mpymc3\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mpm\u001b[39;00m\n\u001b[1;32m 7\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01marviz\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01maz\u001b[39;00m\n\u001b[1;32m 9\u001b[0m az\u001b[38;5;241m.\u001b[39mstyle\u001b[38;5;241m.\u001b[39muse(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124marviz-darkgrid\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/pymc3/__init__.py:23\u001b[0m, in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 20\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mplatform\u001b[39;00m\n\u001b[1;32m 22\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01msemver\u001b[39;00m\n\u001b[0;32m---> 23\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\n\u001b[1;32m 25\u001b[0m _log \u001b[38;5;241m=\u001b[39m logging\u001b[38;5;241m.\u001b[39mgetLogger(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpymc3\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 27\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m logging\u001b[38;5;241m.\u001b[39mroot\u001b[38;5;241m.\u001b[39mhandlers:\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/theano/__init__.py:83\u001b[0m, in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 75\u001b[0m \u001b[38;5;66;03m# This is the api version for ops that generate C code. External ops\u001b[39;00m\n\u001b[1;32m 76\u001b[0m \u001b[38;5;66;03m# might need manual changes if this number goes up. An undefined\u001b[39;00m\n\u001b[1;32m 77\u001b[0m \u001b[38;5;66;03m# __api_version__ can be understood to mean api version 0.\u001b[39;00m\n\u001b[1;32m 78\u001b[0m \u001b[38;5;66;03m#\u001b[39;00m\n\u001b[1;32m 79\u001b[0m \u001b[38;5;66;03m# This number is not tied to the release version and should change\u001b[39;00m\n\u001b[1;32m 80\u001b[0m \u001b[38;5;66;03m# very rarely.\u001b[39;00m\n\u001b[1;32m 81\u001b[0m __api_version__ \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[0;32m---> 83\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m scalar, tensor\n\u001b[1;32m 84\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mcompile\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m (\n\u001b[1;32m 85\u001b[0m In,\n\u001b[1;32m 86\u001b[0m Mode,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 93\u001b[0m shared,\n\u001b[1;32m 94\u001b[0m )\n\u001b[1;32m 95\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mcompile\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mfunction\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m function, function_dump\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/theano/tensor/__init__.py:20\u001b[0m, in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mcompile\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m SpecifyShape, specify_shape\n\u001b[1;32m 10\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mgradient\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m (\n\u001b[1;32m 11\u001b[0m Lop,\n\u001b[1;32m 12\u001b[0m Rop,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 18\u001b[0m verify_grad,\n\u001b[1;32m 19\u001b[0m )\n\u001b[0;32m---> 20\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m nnet \u001b[38;5;66;03m# used for softmax, sigmoid, etc.\u001b[39;00m\n\u001b[1;32m 21\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m sharedvar \u001b[38;5;66;03m# adds shared-variable constructors\u001b[39;00m\n\u001b[1;32m 22\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m (\n\u001b[1;32m 23\u001b[0m blas,\n\u001b[1;32m 24\u001b[0m blas_c,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 29\u001b[0m xlogx,\n\u001b[1;32m 30\u001b[0m )\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/theano/tensor/nnet/__init__.py:3\u001b[0m, in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mwarnings\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01m.\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m opt\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mabstract_conv\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m conv2d \u001b[38;5;28;01mas\u001b[39;00m abstract_conv2d\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mabstract_conv\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m conv2d_grad_wrt_inputs, conv3d, separable_conv2d\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/theano/tensor/nnet/opt.py:32\u001b[0m, in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 24\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mnnet\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mblocksparse\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m (\n\u001b[1;32m 25\u001b[0m SparseBlockGemv,\n\u001b[1;32m 26\u001b[0m SparseBlockOuter,\n\u001b[1;32m 27\u001b[0m sparse_block_gemv_inplace,\n\u001b[1;32m 28\u001b[0m sparse_block_outer_inplace,\n\u001b[1;32m 29\u001b[0m )\n\u001b[1;32m 31\u001b[0m \u001b[38;5;66;03m# Cpu implementation\u001b[39;00m\n\u001b[0;32m---> 32\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mnnet\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mconv\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m ConvOp, conv2d\n\u001b[1;32m 33\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mnnet\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mcorr\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m CorrMM, CorrMM_gradInputs, CorrMM_gradWeights\n\u001b[1;32m 34\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mnnet\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mcorr3d\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m Corr3dMM, Corr3dMMGradInputs, Corr3dMMGradWeights\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/theano/tensor/nnet/conv.py:20\u001b[0m, in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 18\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mgraph\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mbasic\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m Apply\n\u001b[1;32m 19\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mgraph\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mop\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m OpenMPOp\n\u001b[0;32m---> 20\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m blas\n\u001b[1;32m 21\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mbasic\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m (\n\u001b[1;32m 22\u001b[0m NotScalarConstantError,\n\u001b[1;32m 23\u001b[0m as_tensor_variable,\n\u001b[1;32m 24\u001b[0m get_scalar_constant_value,\n\u001b[1;32m 25\u001b[0m patternbroadcast,\n\u001b[1;32m 26\u001b[0m )\n\u001b[1;32m 27\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mnnet\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mabstract_conv\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m get_conv_output_shape, get_conv_shape_1axis\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/theano/tensor/blas.py:163\u001b[0m, in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 161\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mscalar\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;28mbool\u001b[39m \u001b[38;5;28;01mas\u001b[39;00m bool_t\n\u001b[1;32m 162\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m basic \u001b[38;5;28;01mas\u001b[39;00m tt\n\u001b[0;32m--> 163\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mblas_headers\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m blas_header_text, blas_header_version\n\u001b[1;32m 164\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mopt\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m in2out, local_dimshuffle_lift\n\u001b[1;32m 165\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtheano\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtensor\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtype\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m values_eq_approx_remove_inf_nan\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/theano/tensor/blas_headers.py:1016\u001b[0m, in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 997\u001b[0m header \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m textwrap\u001b[38;5;241m.\u001b[39mdedent(\n\u001b[1;32m 998\u001b[0m \u001b[38;5;124;03m\"\"\"\\\u001b[39;00m\n\u001b[1;32m 999\u001b[0m \u001b[38;5;124;03m static float sdot_(int* Nx, float* x, int* Sx, float* y, int* Sy)\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 1010\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[1;32m 1011\u001b[0m )\n\u001b[1;32m 1013\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m header \u001b[38;5;241m+\u001b[39m blas_code\n\u001b[0;32m-> 1016\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[43mconfig\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mblas__ldflags\u001b[49m:\n\u001b[1;32m 1017\u001b[0m _logger\u001b[38;5;241m.\u001b[39mwarning(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mUsing NumPy C-API based implementation for BLAS functions.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 1020\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mmkl_threads_text\u001b[39m():\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/theano/configparser.py:358\u001b[0m, in \u001b[0;36mConfigParam.__get__\u001b[0;34m(self, cls, type_, delete_key)\u001b[0m\n\u001b[1;32m 356\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m:\n\u001b[1;32m 357\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m callable(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdefault):\n\u001b[0;32m--> 358\u001b[0m val_str \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdefault\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 359\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 360\u001b[0m val_str \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdefault\n",
"File \u001b[0;32m/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/theano/link/c/cmodule.py:2621\u001b[0m, in \u001b[0;36mdefault_blas_ldflags\u001b[0;34m()\u001b[0m\n\u001b[1;32m 2617\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m 2618\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mhasattr\u001b[39m(numpy\u001b[38;5;241m.\u001b[39mdistutils, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m__config__\u001b[39m\u001b[38;5;124m\"\u001b[39m) \u001b[38;5;129;01mand\u001b[39;00m numpy\u001b[38;5;241m.\u001b[39mdistutils\u001b[38;5;241m.\u001b[39m__config__:\n\u001b[1;32m 2619\u001b[0m \u001b[38;5;66;03m# If the old private interface is available use it as it\u001b[39;00m\n\u001b[1;32m 2620\u001b[0m \u001b[38;5;66;03m# don't print information to the user.\u001b[39;00m\n\u001b[0;32m-> 2621\u001b[0m blas_info \u001b[38;5;241m=\u001b[39m \u001b[43mnumpy\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdistutils\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m__config__\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mblas_opt_info\u001b[49m\n\u001b[1;32m 2622\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 2623\u001b[0m \u001b[38;5;66;03m# We do this import only here, as in some setup, if we\u001b[39;00m\n\u001b[1;32m 2624\u001b[0m \u001b[38;5;66;03m# just import theano and exit, with the import at global\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 2630\u001b[0m \u001b[38;5;66;03m# This happen with Python 2.7.3 |EPD 7.3-1 and numpy 1.8.1\u001b[39;00m\n\u001b[1;32m 2631\u001b[0m \u001b[38;5;66;03m# isort: off\u001b[39;00m\n\u001b[1;32m 2632\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mnumpy\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mdistutils\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01msystem_info\u001b[39;00m \u001b[38;5;66;03m# noqa\u001b[39;00m\n",
"\u001b[0;31mAttributeError\u001b[0m: module 'numpy.distutils.__config__' has no attribute 'blas_opt_info'"
]
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"import scipy.stats as stats\n",
"import numpy as np\n",
"import pandas as pd\n",
"import seaborn as sns\n",
"import pymc3 as pm\n",
"import arviz as az\n",
"\n",
"az.style.use('arviz-darkgrid')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0eec17a3",
"metadata": {},
"outputs": [],
"source": [
"np.random.seed(123)\n",
"trials = 4\n",
"theta_real = 0.35 # unknown value in a real experiment\n",
"data = stats.bernoulli.rvs(p=theta_real, size=trials)"
]
},
{
"cell_type": "markdown",
"id": "933ce50c",
"metadata": {},
"source": [
"### 2.2.1 建立模型\n",
"\n",
"现在有了数据，需要进一步指定模型。回想一下，模型可通过指定似然函数和先验分布完成。对于似然，可用参数为 $n=1$ 和 $p=\\theta$ 的二项分布或者伯努利分布来描述。对于先验，用参数为 $\\alpha=\\beta=1$ 的贝塔分布描述，该贝塔分布与 [0,1] 区间内的均匀分布等价。用数学表达式描述如下： \n",
"\n",
"$$\n",
"\\theta \\sim \\operatorname{Beta}(\\alpha, \\beta) \n",
"$$\n",
"\n",
"$$\n",
"y \\sim \\operatorname{Bern}(p=\\theta) \\tag{式2.1}\n",
"$$\n",
"\n",
"\n",
"该统计模型与 `PyMC3` 的语法几乎一一对应。"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "aa027580",
"metadata": {},
"outputs": [],
"source": [
"with pm.Model() as our_first_model:\n",
" θ = pm.Beta('θ', alpha=1., beta=1.)\n",
" y = pm.Bernoulli('y', p=θ, observed=data)\n",
" trace = pm.sample(1000, random_seed=123)"
]
},
{
"cell_type": "markdown",
"id": "c8b2aafa",
"metadata": {},
"source": [
"第 1 行代码构建了一个模型的容器， `PyMC3` 使用 `with` 语法将所有位于该语法块内的代码都指向同一个模型，你可以把它看作是简化模型描述的“语法糖”，此处将模型命名为 `our_first_model` 。\n",
"\n",
"第 2 行代码指定了先验，可以看到语法与数学表示法很接近。\n",
"\n",
"```{tip}\n",
"我们把随机变量命名为 θ ，其名称与贝塔函数的第 1 个参数名保持一样，这是一个能避免混淆好习惯。 在后期我们可以通过变量名从后验样本中提取信息。注意此处变量指随机变量，可以被视为从某个分布中生成样本的方法，而不是某个确切值。\n",
"```\n",
"\n",
"第 3 行代码用跟先验相同的语法描述了似然，唯一不同的是： 我们用 `observed` 参数向似然传递了观测数据 `data`。`data` 可以是 `Python` 语言的一个 `list` 或者 `Numpy` 的一个 `array` 或者 `Pandas` 的一个 `DataFrame` 。\n",
"\n",
"### 2.2.2 执行推断\n",
"\n",
"最后一行可以视为按下了*推断按钮*。我们要求从后验中做 1000 次采样， 并且将其存在 `trace` 对象中。在这行代码背后，`PyMC3` 将执行数百个贝叶斯推断任务！如果您运行该代码，您将收到如下消息：\n",
"\n",
"```\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [θ]\n",
"100%|██████████| 3000/3000 [00:00<00:00, 3695.42it/s]\n",
"```\n",
"\n",
"第一行和第二行告诉我们，`PyMC3` 已经自动分配了 NUTS 采样器（一种非常适用于连续型变量的 `MCMC` 推断引擎），并使用了一种方法来初始化该采样器。第三行说明 `PyMC3` 将同时运行两个链，因此我们可以从后验中同时获得两个独立样本。链的确切数量根据计算机中处理器的数量计算，您也可以使用 `sample` 函数的 `chains` 参数来更改链的数量。\n",
"\n",
"下一行告诉我们哪些变量是由哪个采样器采样的。对于此特定情况，此行不会添加新信息。因为 `NUTS` 是用来对我们拥有的唯一随机变量 $θ$ 进行采样的。但情况并不总是如此，因为 `PyMC3` 可以将不同采样器分配给不同的随机变量。这一般由 `PyMC3` 根据变量属性自动完成，以确保为每个变量使用尽可能好的采样器。用户可以使用 `sample` 函数的 `step` 参数为随机变量手动分配采样器。\n",
"\n",
"最后一行是进度条，其中有几个度量指标指示了采样器的工作速度，包括每秒迭代的次数。如果运行代码，您将看到进度条更新得非常快。在示例中，我们看到的是采样器完成了其工作的最后一个阶段。数字是 `3000/3000`，其中第一个数字是运行的采样器编号（从 1 开始），最后一个数字是样本总数。您会注意到，我们要求了 1,000 个样本，但 `PyMC3` 实际计算 3,000 个样本。我们每条链会有 500 个样本用于自动优化调整采样算法（本例中为 NUTS）。默认情况下，这些样本将被丢弃。同时，两条链各有 1000 个有效的抽取，因此总共生成了 3000 个样本。优化调整阶段帮助 `PyMC3` 从后验中提供更加可靠样本。我们也可以使用 `sample` 函数的 `tune` 参数来更改用于调优的步骤数量。\n",
"\n",
"### 2.2.3 分析和诊断后验分布\n",
"\n",
"#### 2.2.3.1 诊断工具和方法\n",
"\n",
"现在我们根据有限数量的样本对后验做出了（近似）推断，接下来要做的是检查推断是否合理。可以通过可视化的或者定量的手段做一些检查，尝试从中发现样本中的问题。需要强调的是：诊断并不能证明得到的分布是否正确，但能够提供证明样本合理性的证据。\n",
"\n",
"**（1）可视化诊断工具**\n",
"\n",
"`ArviZ` 中的 `plot_trace` 函数非常适合执行可视化诊断任务："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a14a3adf",
"metadata": {},
"outputs": [],
"source": [
"az.plot_trace(trace)"
]
},
{
"cell_type": "markdown",
"id": "f52460e9",
"metadata": {},
"source": [
"\n",
"\n",
"通过使用 `az.plot_trace`，我们可以为每个非观测变量绘制两个子图，由于本例中只有唯一的非观测变量 $\\theta$ ，因此只需要绘制了 $\\theta$ 的后验分布图和样本生成轨迹图。\n",
"\n",
"> $\\mathbb{y}$ 表示观测变量，其值已知，无需采样。\n",
"\n",
"上图中，左边为样本经过核密度估计（ `KDE` ）后的图，类似于直方图的平滑版本；右图为采样过程中每一步的采样值，被称为轨迹图。轨迹图中可以直观地反映在后端获得的采样结果。你可以尝试将 `PyMC3` 的结果与上一章通过解析方式获得的结果进行比较。\n",
"\n",
"**（2）生成关于后验分布的摘要信息**\n",
"\n",
"`ArviZ` 还提供了其他几个绘图来帮助解释轨迹图，我们将在后面逐步看到它们。此外 `az.summary` 以 `DataFrame` 形式提供了关于后验分布的摘要数据。"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "68d9e5da",
"metadata": {},
"outputs": [],
"source": [
"az.summary(trace)"
]
},
{
"cell_type": "markdown",
"id": "1761910e",
"metadata": {},
"source": [
"\n",
"\n",
"通过上图所示的摘要数据，可以得到均值 $\\mu$ 、标准差 $\\sigma$ 和 ` 94% HPDI` 区间。正如在第一章中所讨论，可以使用这些统计特征来解释和汇报贝叶斯推断的结果。后两个指标与诊断样本相关。详细信息请参阅第 8 章。\n",
"\n",
"**（3）绘制后验并标识关键摘要数据**\n",
"\n",
"另一种直观总结后验的方法是使用 `ArviZ` 附带的 `plot_posterior` 函数。默认情况下，`plot_posterior` 函数会为离散型随机变量显示直方图，为连续型随机变量显示 `KDE` 。\n",
"\n",
"图中还标识了后验分布的均值（你也可以使用 `point_estimate` 参数指定想要的点估计值，如：中位数、众数等），`94% HPDI` 作为图底部的一条黑线。可以使用 `credible_interval` 参数为 `HPD` 设置不同的间隔值。\n",
"\n",
"> **注意：** 此类型的图由 `John K.Kruschke` 在其著作《 `Doing Bayesian Data Analysis` 》中引入。"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0d000a61",
"metadata": {},
"outputs": [],
"source": [
"az.plot_posterior(trace)"
]
},
{
"cell_type": "markdown",
"id": "cc27cf5b",
"metadata": {},
"source": [
"\n",
"\n",
"#### 2.2.3.2 诊断出问题怎么办？\n",
"\n",
"如果通过样本发现了问题，解决办法有如下几种：\n",
"\n",
"**（1）增加样本数量**\n",
"\n",
"从样本链（轨迹）的前面部分去掉一定数量的样本，被称为老化（`Burn-in`）。在实践中，`MCMC` 方法通常需要经过一段时间采样之后，才能够趋近于真正的目标分布。老化在无限多次采样中并不是必须的，因为这部分并没有包含在马尔科夫理论中。事实上，去掉前面部分样本只不过是在有限次采样中提升结果的一个小技巧。\n",
"\n",
"注意不要被 “数学对象” 和 “数学对象的近似” 弄糊涂了，球体、高斯分布以及马尔科夫链等数学对象只存在于柏拉图式想象世界中，并不存在于不完美但却真实的世界中。\n",
"\n",
"**（2）重新参数化你的模型**\n",
"\n",
"也就是说换一种不同但却等价的方式描述模型。\n",
"\n",
"**（3）转换数据**\n",
"\n",
"这么做有可能得到更高效的采样。转换数据的时候需要注意对结果在转换后的空间内进行解释，或者再重新转换回去，然后再解释结果。\n",
"\n",
"#### 2.2.3.3 收敛性问题\n",
"\n",
"在 `az.plot_trace` 图中，我们需要观察什么呢？首先，`KDE` 图看起来应该是光滑曲线。通常随着数据的增加，根据中心极限定理，参数分布会趋近于高斯分布。当然，这并不一定是正确的。右侧图看起来应该像白噪声，也就是说有很好的混合度，通常不希望看到任何可以识别的模式，相反希望看到曲线在某个值附近震荡。对于多峰分布或者离散分布，我们希望曲线不要在某个值或区域停留过多时间，我们希望看到采样值在多个区间自由移动。此外，我们希望迹表现出稳定的相似性，也就是说，前 10% 看起来跟后 50% 或者 10% 差不多。再次强调，我们不希望看到规律的模式，相反期望看到的是噪声。下图展示了一些迹呈现较好混合度（右侧）与较差混合度（左侧）的对比。\n",
"\n",
"![](\"https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210504202103_98.webp\") \n",
"\n",
"如果迹的前面部分跟其他部分看起来不太一样，那就意味着需要进行老化处理，如果其他部分没有呈现稳定的相似性或者可以看到某种模式，那这意味着需要更多的采样，或者需要更换采样方法或者参数化方法。对于一些复杂的模型，我们可能需要结合使用前面所有的策略。\n",
"\n",
" `PyMC3` 可以实现并行地运行一个模型多次，因而对同一个参数可以得到多条并行的迹。这可以通过在采样函数中指定 `njobs` 参数实现。此时使用 `plot_trace` 函数，便可在同一幅图中得到同一个参数的所有迹。由于每组迹都是相互独立的，所有的迹看起来都应该差不多。除了检查收敛性之外，这些并行的迹也可以用于推断，我们可以将这些并行的迹组合起来提升采样大小而不是扔掉多余的迹："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2fe9ad47",
"metadata": {},
"outputs": [],
"source": [
"with our_first_model:\n",
" step = pm.Metropolis()\n",
" multi_trace = pm.sample(1000, step=step)\n",
"az.plot_trace(multi_trace, lines={'theta':theta_real});"
]
},
{
"cell_type": "markdown",
"id": "e9e7ae50",
"metadata": {},
"source": [
"一种定量检测收敛性的方法是 **Gelman-Rubin 检验** （ `az.summary` 输出的表格中的 $\\hat R$ 值）。该检验的思想是比较不同迹之间的差异和迹内部的差异，因此，需要多组迹来执行检验。理想状态下，我们希望得到 $\\hat R=1$ 。但根据经验，值低于 1.1 也可以认为是收敛的，更高的值则意味着不收敛。\n",
"\n",
"#### 2.2.3.4 自相关问题\n",
"\n",
"最理想的采样应该不会是自相关的，也就是说，某点的值应该与其他点的值相互独立。在实际中，从 `MCMC` 方法（特别是 `Metropolis-Hastings`）中得到的采样值通常是自相关的。由于参数间的相互依赖关系，模型会导致更多自相关采样。 `PyMC3` 有一个很方便的函数用来描述自相关。"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4ae37725",
"metadata": {},
"outputs": [],
"source": [
"az.plot_autocorr(trace)"
]
},
{
"cell_type": "markdown",
"id": "eb8e13ca",
"metadata": {},
"source": [
"\n",
"\n",
"该图显示了采样值与相邻连续点之间的平均相关性。理想状态下，不应看到自相关性，实际中希望看到自相关性降低到较低水平。参数越自相关，要达到指定精度的采样次数就需要越多，也就是说，自相关性不利于降低采样次数。\n",
"\n",
"#### 2.2.3.5 有效后验采样次数\n",
"\n",
"有自相关性的采样要比没有自相关性的采样包含的信息量更少，给定采样大小和采样的自相关性之后，可以尝试估计出该采样在采样次数为多少时，没有自相关性且包含信息量不变，该值称为有效采样次数（ `az.summary` 输出表格中的 `eff_n` 值）。理想情况下，两个值是一模一样的；二者越接近，采样效率越高。有效采样次数可以作为一个参考，如果想要估计出分布的均值，需要的最小采样数至少为 100；如果想要估计出依赖于尾部分布的量，比如可信区间的边界，那么可能需要 1000 到 10000 次采样。\n",
"\n",
"提高采样效率的一个方法是换一个更好的采样方法；另一个办法是转换数据或者对模型重新设计参数，此外，还有一个常用的办法是对采样链压缩。所谓压缩其实就是每隔 k 个观测值取一个，在 `Python` 中我们称为切片。压缩会降低自相关性，但代价是同时降低了样本量。因此，实际使用中通常更倾向于增加样本量而不是切片。\n",
"\n",
"### 2.2.4 基于后验推断的决策\n",
"\n",
"有时候，仅描述后验还不够，还需要根据推断结果做决策。也就是说，我们不得不将一个连续的估计（后验分布）简化为一个二分法，如：是或否、健康或生病、污染或安全等。\n",
"\n",
"回到抛硬币问题上，我们需要回答硬币是否公平。一枚公平的硬币 $\\theta$ 值应当恰好是 0.5，严格来说，出现这种情况的概率为 0，因而实际中会对定义稍稍放松，即假如一枚硬币的 $\\theta$ 值在 0.5 左右，就可以认为判定该硬币是公平的。此处“左右”的具体含义依赖于具体问题，并没有一个满足所有问题的普适准则。因此决策是偏主观的，我们的任务就是根据目标做出最可能的决策。\n",
"\n",
"直观上，一个明智的做法是将 `HPD` 区间与我们感兴趣的值进行比较，在抛硬币的例子中，该值是 0.5。前面的图中可以看出 `HPD` 的范围是 0.02～0.71 ，包含 0.5 ，根据后验分布来看，硬币似乎倾向于反面朝上，但我们并不能完全排除这枚硬币的公平的。如果想要一个更清晰的决定，将需要收集更多数据来减少后验数据的扩散，或者需要找出如何定义一个更有信息量的先验。\n",
"\n",
"**（1）`实际等效区间 ROPE` 与决策**\n",
"\n",
"严格地说，恰好观察到 0.5 的机会为零。此外，实践中我们通常不关心确切结果，而是在一定范围内的结果。因此，可以放宽公平的定义，例如我们可以说区间 [0.45，0.55] 中的任何值实际上等于 0.5。我们称这个区间为实际等效区间（`ROPE`）。一旦定义了 `ROPE`，将其与最高后验密度（`HPD`）比较，至少可以得到三个场景：\n",
"\n",
"- `ROPE` 与 `HPD` 区间没有重叠，因此我们可以说硬币是不公平的。\n",
"- `ROPE` 包含整个 `HPD` 区间，我们可以认为硬币是公平的。\n",
"- `ROPE` 与 `HPD` 区间部分重叠，此时我们不能判断硬币是否公平。\n",
"\n",
"当然，如果选择区间 [0,1] 作为 `ROPE`，那么不管结果怎样我们都会说这枚硬币是公平的，不过恐怕没人会同意这个对 `ROPE` 的定义。\n",
"\n",
"> `ROPE` 是根据背景或专业知识选择的任意区间，假设在该区间内的任何值都具有实际等效性。\n",
"\n",
"plot_posterior 函数可以用来刻画 `ROPE`。从图中可以看到，`ROPE` 是一段较宽的半透明的红色线段，同时上面有两个数值表示 `ROPE` 的两个端点。"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e3ec1786",
"metadata": {},
"outputs": [],
"source": [
"az.plot_posterior(trace, ROPE=[0.45, .55])"
]
},
{
"cell_type": "markdown",
"id": "e3b3147f",
"metadata": {},
"source": [
"\n",
"\n",
"还可以给 `plot_posterior` 传递一个参考值 `ref_val`，例如 0.5，用来和后验进行对比。从图中可以看出我们会得到一个橙色的垂直线以及大于该值和小于该值的后验比例。"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "58932e15",
"metadata": {},
"outputs": [],
"source": [
"az.plot_posterior(trace,ref_val=0.5)"
]
},
{
"cell_type": "markdown",
"id": "5ef056e8",
"metadata": {},
"source": [
"\n",
"\n",
"关于如何使用 `ROPE` 的更多细节，可以阅读 Kruschke 的《`Doing Bayesian Data Analysis`》一书的第 12 章。该章还讨论了在贝叶斯框架下如何做假设检验，以及一些（贝叶斯或者非贝叶斯的）假设检验方面的警告。\n",
"\n",
"**（2）损失函数、点估计与决策**\n",
"\n",
"如果你认为 `ROPE` 准则听起来有点笨拙，而你想要更正式的东西，那么损失函数就是你想要的！要做出好的决策，重要的是参数的估计值有尽可能高的精度，但也要考虑犯错的代价。成本/收益的权衡在数学上可以使用损失函数形式化。损失函数或其逆函数的名称在不同的领域中各不相同，可以找到成本函数、目标函数、适应度函数、效用函数等名称。无论名称如何，关键思想都是使用一个函数来捕获参数的真实值和估计值的差别。损失函数值越大，估计就越差（根据损失函数）。\n",
"\n",
"损失函数的一些常见示例包括：\n",
"\n",
"- 二次损失 $(\\theta-\\hat{\\theta})^{2}$ \n",
"- 绝对损失 $|\\theta-\\hat{\\theta}|$ \n",
"- 0-1 损失 $I(\\theta \\neq \\hat{\\theta})$ ，其中 $I$ 为指示函数\n",
"\n",
"实践中通常手头没有真实参数 $\\theta$ 的值，仅有一个后验分布形式的估计。因此，我们能做的就是找出最小化期望损失函数的 $\\hat \\theta$ 值。期望损失函数是指整个后验分布的平均损失函数。在下面代码块中，我们有两个损失函数：绝对损失（ `lossf_a`）和平方损失（ `lossf_b`）。我们将尝试超过 200 个网格的 $\\hat \\theta$ 值，然后绘制其曲线，还将包括最小化每个损失函数的 $\\hat \\theta$ 值："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8132678c",
"metadata": {},
"outputs": [],
"source": [
"grid = np.linspace(0, 1, 200)\n",
"θ_pos = trace['θ']\n",
"lossf_a = [np.mean(abs(i - θ_pos)) for i in grid]\n",
"lossf_b = [np.mean((i - θ_pos)**2) for i in grid]\n",
"for lossf, c in zip([lossf_a, lossf_b], ['C0', 'C1']):\n",
" mini = np.argmin(lossf)\n",
" plt.plot(grid, lossf, c)\n",
" plt.plot(grid[mini], lossf[mini], 'o', color=c)\n",
" plt.annotate('{:.2f}'.format(grid[mini]),\n",
" (grid[mini], lossf[mini] + 0.03), color=c)\n",
" plt.yticks([])\n",
" plt.xlabel(r' $\\hat \\theta $ ')"
]
},
{
"cell_type": "markdown",
"id": "2a633a4c",
"metadata": {},
"source": [
"\n",
"\n",
"正如所看到的，结果显示 `lossf_a` 的最优解是 $\\hat \\theta=0.32$ ，`lossf_b` 的最优解是 $\\hat \\theta=0.33$ 。该结果中比较有趣的是，前一个值等于后验的中位数，后一个值等于后验的平均值。通过计算 `np.Means(θ_pos)`、`np.Medium(θ_pos)` 可以检查这个情况。这似乎在暗示：**不同的损失函数与不同的点估计有关**。\n",
"\n",
"OK，如果想要形式化的结果并给出点估计，必须决定想要哪个损失函数，或者反过来，如果选择一个给定点估计，就隐含地（甚至可能是无意识地）决定了一个损失函数。显式选择损失函数的好处是可以根据问题调整函数，而不是使用一些可能不适合特定情况的预定义规则。\n",
"\n",
"在许多问题中，做出决定的成本是不对称的；例如，决定给五岁以下的儿童接种疫苗还是不接种疫苗。一个糟糕的决定可能会造成数千人生命损失，并产生健康危机，而这场危机本可以通过接种一种廉价而安全的疫苗来避免。因此，如果问题需要的话，可以构造一个不对称损失函数。\n",
"\n",
"同样重要的是要注意，由于后验分布是数字采样的形式，因此理论上可以计算任意复杂的损失函数，而不需要受数学形式上的限制。\n",
"\n",
"以下只是一个愚蠢的例子："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "94377017",
"metadata": {},
"outputs": [],
"source": [
"lossf = []\n",
"for i in grid:\n",
" if i < 0.5:\n",
" f = np.mean(np.pi * θ_pos / np.abs(i - θ_pos))\n",
" else:\n",
" f = np.mean(1 / (i - θ_pos))\n",
" lossf.append(f)\n",
"mini = np.argmin(lossf)\n",
"plt.plot(grid, lossf)\n",
"plt.plot(grid[mini], lossf[mini], 'o')\n",
"plt.annotate('{:.2f}'.format(grid[mini]),\n",
" (grid[mini] + 0.01, lossf[mini] + 0.1))\n",
"plt.yticks([])\n",
"plt.xlabel(r' $\\hat \\theta $ ')"
]
},
{
"cell_type": "markdown",
"id": "d6f9534b",
"metadata": {},
"source": [
"\n",
"\n",
"话虽如此，我想澄清一点。这并不是说，每次人们使用点估计时，他们都是真的在考虑损失函数。事实上，在我或多或少熟悉的许多科学领域，损失函数并不是很常见。人们经常选择中位数，只是因为它对异常值比平均值更可靠，或者仅仅因为它是一个简单而熟悉的概念，或者因为他们认为它们的可观测性真的是某种程度上某个过程的平均值，比如分子相互弹跳，或者基因与自己和环境相互作用。\n",
"\n",
"我们刚刚看到对损失函数的简单而浅显的介绍。如果你想了解更多这方面知识，可以尝试阅读决策理论，这是一个研究正式决策的领域。\n",
"\n",
"## 2.3 高斯模型与更为稳健的学生 $t$ 模型\n",
"\n",
"前面我们用`贝塔--二项分布`模型介绍了贝叶斯思想，该模型比较简单。另外还有一个非常简单的模型是高斯分布或者叫高斯分布。从数学角度来看，高斯分布受欢迎的原因是其处理起来非常简单，例如：高斯分布的均值计算，其共轭先验还是高斯分布等优点。此外，许多现象可以用高斯分布来近似；本质上来说，每当测量某种均值时，只要采样样本量足够大，观测值的分布就会呈现高斯分布。至于这种近似什么时候是对的，什么时候是错的，可以了解下中心极限定理。\n",
"\n",
"此处举一个例子：身高是受基因和许多环境因素影响的，因而我们观测到的成年人的身高符合高斯分布。不过事实上，我们得到身高观测数据其实是一个双峰分布，男人和女人的身高分布重叠在了一起。\n",
"\n",
"总的来说，高斯分布用起来很简单，而且自然界中随处可见，这也是为什么许多统计方法都基于高斯分布。学习如何构建这类模型非常重要，此外，学会如何放宽高斯分布的假设也同等重要，这一点在贝叶斯框架中利用 `PyMC3` 之类的现代计算工具很容易处理。\n",
"\n",
"### 2.3.1 高斯推断\n",
"\n",
"下面的例子与核磁共振实验有关，核磁共振是一种研究分子和生物的技术。下面这组数据，可能来自一群人身高的测量值、回家的平均时间、从超市买回来橙子的重量、大壁虎的伴侣个数或者任何可以用高斯分布近似的测量值。在这个例子中，我们有 48 个测量值："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1f783f19",
"metadata": {},
"outputs": [],
"source": [
"data = np.loadtxt('../data/chemical_shifts.csv')\n",
"az.plot_kde(data, rug=True)\n",
"plt.yticks([0], alpha=0"
]
},
{
"cell_type": "markdown",
"id": "a6a66a9b",
"metadata": {},
"source": [
"除了两个远离平均值的数据点外，该数据集的 `KDE` 图显示出类似高斯的分布：\n",
"\n",
"\n",
"\n",
"暂且先不考虑偏离均值的那两个点，假设以上分布就是高斯分布。由于不知道均值和方差，需要先对这两个参数设置先验。然后，顺理成章地得到如下模型：\n",
"\n",
"$$\n",
"\\mu \\sim U(l, h) \n",
"$$\n",
"\n",
"$$\n",
"\\sigma \\sim\\left|\\mathcal{N}\\left(0, \\sigma_{\\sigma}\\right)\\right| \n",
"$$\n",
"\n",
"$$\n",
"y \\sim \\mathcal{N}(\\mu, \\sigma) \\tag{式2.2} \n",
"$$\n",
"\n",
"\n",
"\n",
"\n",
"其中， $\\mu$ 来自上下界分别为 $l$ 和 $h$ 的均匀分布， $\\sigma$ 来自标准差为 $\\sigma_\\sigma$ 的半高斯分布。半高斯分布和普通高斯分布很像，不过只包含正数，看起来就好像将普通的高斯分布沿着均值对折了。通过从高斯分布中采样，然后取绝对值，可以获取半高斯分布的样本。最后，在我们的模型中，数据 $y$ 来自参数为 $\\mu$ 和 $\\sigma$ 的高斯分布，可以用 Kruschke 图表示法将其画出来：\n",
"\n",
"\n",
"\n",
"如果不知道 $\\mu$ 和 $\\sigma$ 的值，可以通过先验来表示该未知信息。例如，可以将均匀分布的上下界分别设为（$l = 40, h = 75$），该范围要比数据本身的范围稍大一些。或者，可以根据先验知识设得更广一些，比如知道这类观测值不可能小于 $0$ 或者大于 $100$ ，可以将均匀先验参数设为（$l = 0,h = 100$）。对于半高斯分布而言，可以把 $\\sigma_\\sigma$ 的值设为 $10$ ，该值相对于数据分布而言算是较大的。 利用 `PyMC3` ，可以将模型编码为："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5a78fbe7",
"metadata": {},
"outputs": [],
"source": [
"with pm.Model() as model_g:\n",
" μ = pm.Uniform('μ', lower=40, upper=70)\n",
" σ = pm.HalfNormal('σ', sd=10)\n",
" y = pm.Normal('y', mu=μ, sd=σ, observed=data)\n",
" trace_g = pm.sample(1000)\n",
"az.plot_trace(trace_g)"
]
},
{
"cell_type": "markdown",
"id": "cc486b07",
"metadata": {},
"source": [
"\n",
"\n",
"您可能已经注意到了，使用 `ArviZ` 的 `plot_trace` 函数绘制的图中，每个未知参数都占用一行。对于这个模型，后验是一个二维分布，所以图中分别显示了每个参数的边缘分布。我们可以使用 `ArviZ` 中的 `plot_joint` 绘制联合分布的函数来查看二维后验分布的形态，以及 $μ$ 和 $σ$ 的边缘分布："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "aa229c6d",
"metadata": {},
"outputs": [],
"source": [
"az.plot_joint(trace_g, kind='kde', fill_last=False)"
]
},
{
"cell_type": "markdown",
"id": "ecdb0c85",
"metadata": {},
"source": [
"\n",
"\n",
"如果要访问存储在轨迹对象中的任何参数值，可以使用相关参数的名称做为轨迹索引。因此，可以获得一个 NumPy 数组。尝试执行 `trace_g['σ']` 或 `az.plot_kde(trace_g['σ'])`。\n",
"\n",
"> 顺便说一下，使用 Jupyter Notebook/lab，您可以直接在代码单元格中写入 Latex 的公式代码 \\sigma，然后按 Tab 键来获取字符。\n",
"\n",
"我们将打印该摘要以供以后使用："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f6d0ae93",
"metadata": {},
"outputs": [],
"source": [
"az.summary(trace_g)"
]
},
{
"cell_type": "markdown",
"id": "6e00fca6",
"metadata": {},
"source": [
"\n",
"\n",
"---\n",
"\n",
"现在已经有了后验分布，可以用它来生成模拟数据，并检查模拟数据与观测数据的一致性。在第 1 章中曾将这种比较称为`后验预测检查`，因为我们要使用后验分布计算预测分布，并使用预测分布来检查模型。\n",
"\n",
"使用 `PyMC3` 的 `sample_posterior_predictive` 函数可以非常容易获得后验预测分布的样本。以下代码将从后端生成 100 个预测，每个预测的大小与数据相同。请注意，我们必须将轨迹和模型传递给 `sample_posterior_predictive` 函数，其他参数是可选的："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5e4cdac7",
"metadata": {},
"outputs": [],
"source": [
"y_pred_g = pm.sample_posterior_predictive(trace_g, 100, model_g)"
]
},
{
"cell_type": "markdown",
"id": "be984f15",
"metadata": {},
"source": [
"`y_pred_g` 变量是一个 “key-value” 字典，`keys` 是模型中的可观测变量名称，`values` 是形状为（`samples`、`size`）的数组，在本例中为 $(100，len(data))$ 。本例中我们将一个字典，因为该模型仅有一个可观测变量。可以使用 `plot_ppc` 函数进行可视化后验预测检查："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bc482544",
"metadata": {},
"outputs": [],
"source": [
"data_ppc = az.from_PyMC3(trace=trace_g, posterior_predictive=y_pred_g)\n",
"ax = az.plot_ppc(data_ppc, figsize=(12, 6), mean=False)\n",
"ax[0].legend(fontsize=15)"
]
},
{
"cell_type": "markdown",
"id": "2f3e57c2",
"metadata": {},
"source": [
"\n",
"\n",
"上图中，黑色的线是观测数据的 `KDE` ，半透明（青色）线反映了我们对预测结果的不确定性。图中可以看到模拟数据的平均值稍微向右偏移，并且模拟数据的方差似乎比实际数据方差更大。我相信，这是观测数据中最右边那两个脱离群体的观察点所导致的。此时，我们能不能用这张图自信地说，模型有问题，需要改变？像往常一样，模型的解释和评估取决于上下文。根据测量经验和数据使用方式，我个人认为该模型已经是一个足够合理的数据表示方法了，对我的大多数分析都够用。但在下一节中，我们会学习如何改进 `model_g` 并获得与数据更接近的预测。\n",
"\n",
"### 2.3.2 更稳健的推断\n",
"\n",
"对于 `model_g` 模型，您可能会有异议，那就是假设观测数据呈高斯分布，但实际数据分布的尾部有两个特殊的数据点，使得高斯假设有点牵强。对于高斯分布，其尾部应随远离平均值而迅速下降，因此，当尾部出现这两个点时，高斯分布的直接反应就是均值向这些点移动，同时可能增加标准差。可以想象这些点对决定高斯分布的参数具有过大的影响。\n",
"\n",
"怎么处理呢？ \n",
"\n",
"有两个最常用的做法，一是将两个点作为异常样本剔除，二是调整模型。\n",
"\n",
"#### （1）剔除异常值\n",
"\n",
"第一种做法是将两个点视为异常并将其剔除，因为这些点有可能是仪器异常或者人为疏忽导致的，又或是处理数据时代码出了问题。但更多时候，我们希望有一些可以遵循的处理规则，来自动消除异常点，其中两个常见的规则如下：\n",
"\n",
"- 低于下四分位数 1.5 倍或高于上四分位数 1.5 倍的数据点视为异常值\n",
"\n",
"- 低于或大于数据标准差两倍的数据点视为异常值\n",
"\n",
"\n",
"#### （2）调整模型\n",
"\n",
"除了改变原始数据外，还可以修改模型。按照贝叶斯思想，**我们更倾向于通过使用不同的先验或似然将假设编码到模型中，而不是通过剔除异常值这种方法**。\n",
"\n",
"一个解决异常值的常用方法是将高斯分布替换成 $\\boldsymbol{t}$ 分布。 $\\boldsymbol{t}$ 分布有 3 个参数：`均值 mean `、`尺度 scale` 、 `自由度 df`（通常用 $\\nu$ 表示，取值范围为 $[0,∞]$ )。根据 Kruschke 的命名方式，我们将 $\\nu$ 称为高斯参数，该参数决定了 $\\boldsymbol{t}$ 分布与高斯分布的相似程度。\n",
"\n",
"- 当 $\\boldsymbol{\\nu \\leq 1}$ 时\n",
" - 该分布没有准确定义的均值（当然，实际从 $\\boldsymbol{t}$ 分布得到的采样终究是一些数字，总是可以算出经验性的均值，不过理论上没有一个准确定义的均值）。 直观地说，这可以理解为：$\\boldsymbol{t}$ 分布的尾部太重了，以至于随时从曲线中获得一个样本，都会影响均值的结果，即该数字永远不会收敛于某个固定值。\n",
"\n",
"- 当 $\\boldsymbol{\\nu=1}$ 时\n",
" - $\\boldsymbol{t}$ 分布表现为重尾分布，此时的分布在不同领域也被称做**柯西分布**或者**洛伦兹分布**，后者常用于物理学。尾部更重的意思是：相比高斯分布，更有可能观测到偏离均值的点，即该分布不像高斯分布那样聚集在均值附近。例如：柯西分布中 95% 的值位于 $-12.7$ 和 $+12.7$ 之间，范围很宽;\n",
"\n",
"- 当 $\\boldsymbol{\\nu \\leq 2}$ 时\n",
" - 分布的方差没有明确定义，因此，需要注意 $\\boldsymbol{t}$ 分布的尺度参数与标准差不是同一个概念。对于的分布，方差并没有明确定义，因而也没有明确定义的标准差。 \n",
"\n",
"- 当 $\\boldsymbol{\\nu \\to \\infty}$ 时\n",
" - $\\boldsymbol{t}$ 分布就是高斯分布， `尺度 scale`趋近于高斯分布的标准差；\n",
"\n",
"\n",
"你可以尝试多次运行下面的代码，查看是否有稳定的均值（ 其中 `df` 为自由度）。将参数 `df` 换成一个更大100在尝试一下。"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ad3143e3",
"metadata": {},
"outputs": [],
"source": [
"np.mean(stats.t(loc=0, scale=1, df=1).rvs(100))"
]
},
{
"cell_type": "markdown",
"id": "3de65d79",
"metadata": {},
"source": [
"下面的代码绘制了自由度分别为 1、2、5、30 时的 $\\boldsymbol{t}$ 分布曲线："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b6f0b419",
"metadata": {},
"outputs": [],
"source": [
"x_values = np.linspace(-10, 10, 200)\n",
"for df in [1, 2, 5, 30]:\n",
" distri = stats.t(df)\n",
" x_pdf = distri.pdf(x_values)\n",
" plt.plot(x_values, x_pdf, label=r\\'\\ $\\\\nu\\ $ = {}\\'.format(df))\n",
"x_pdf = stats.norm.pdf(x_values)\\\n",
"plt.plot(x_values, x_pdf, label=r\\'\\ $\\\\nu = \\\\infty\\ $\\') plt.xlabel(\\'x\\')\n",
"plt.ylabel(\\'p(x)\\', rotation=0)\n",
"plt.legend(loc=0, fontsize=14)\n",
"plt.xlim(-7, 7)"
]
},
{
"cell_type": "markdown",
"id": "39918b24",
"metadata": {},
"source": [
"\n",
"\n",
"> 图 2.12 不同 $\\boldsymbol{\\nu}$ 值的学生 $\\boldsymbol{t}$ 分布示意图\n",
"\n",
"利用 $\\boldsymbol{t}$ 分布可以将模型调整为如下形式：\n",
"\n",
"$$\n",
"\\mu \\sim U(l, h) \n",
"$$\n",
"\n",
"$$\n",
"\\sigma \\sim\\left|\\mathcal{N}\\left(0, \\sigma_{\\sigma}\\right)\\right| \n",
"$$\n",
"\n",
"$$\n",
"\\nu \\sim \\operatorname{Exp}(\\lambda) \n",
"$$\n",
"\n",
"$$\n",
"y \\sim \\mathcal{T}(\\mu, \\sigma, \\nu) \\tag{式2.3}\n",
"$$\n",
" \n",
"\n",
"此模型与高斯模型的主要区别是：似然调整为 $\\boldsymbol{t}$ 分布，由于 $\\boldsymbol{t}$ 分布多了一个新的参数 $\\nu$，需要为其增加一个先验。此处计划采用均值为 $30$ 的指数分布。通过上图可以看出，当 $\\nu = 30$ 时， $\\boldsymbol{t}$ 分布看起来与高斯分布很相似。 从图中也可以看出， $\\nu$ 在 $30$ 附近是一个比较适中的值，既可以调大也可以调小，因此属于弱信息先验。我们的模型可以表示如下：\n",
"\n",
"\n",
"\n",
"> 图 2.13 学生 $\\boldsymbol{t}$ 分布作为似然的概率模型可视图\n",
"\n",
"同样， `PyMC3` 让我们只需要几行代码便可修改模型。唯一需要注意的是， `PyMC3` 中指数分布的参数采用的是均值的倒数。"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fe0465d6",
"metadata": {},
"outputs": [],
"source": [
"with pm.Model() as model_t:\n",
" μ = pm.Uniform('μ', 40, 75)\n",
" σ = pm.HalfNormal('σ', sd=10)\n",
" ν = pm.Exponential('ν', 1/30)\n",
" y = pm.StudentT('y', mu=μ, sd=σ, nu=ν, observed=data)\n",
" trace_t = pm.sample(1000)\n",
"\n",
"az.plot_trace(trace_t)"
]
},
{
"cell_type": "markdown",
"id": "c43da4e7",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"> 图 2.14 学生 $\\boldsymbol{t}$ 分布作为似然的后验分布与轨迹图\n",
"\n",
"\n",
"将 `model_g` 的轨迹图（图2.9）与 `model_t` 的跟踪（图2.14）进行比较。现在，打印 `model_t` 的摘要，并将其与 `model_g` 的摘要进行比较。在继续阅读之前，请花点时间找出两个结果之间的差异。你注意到什么有趣的事情了吗？"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "12182551",
"metadata": {},
"outputs": [],
"source": [
"az.summary(trace_t)"
]
},
{
"cell_type": "markdown",
"id": "3fcd040b",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"可以看到，高斯分布模型和 $\\boldsymbol{t}$ 分布模型对 $\\mu$ 的估计比较接近，只相差 $0.5$ 左右，而 $σ$ 的估计则从 3.5 变成了 2.1，这是因为 $\\boldsymbol{t}$ 分布对于偏离均值的点赋予的权重较小。此外还可以看到， $\\nu$ 的值接近 4.5，也就是说，该分布并不太像高斯分布，而是更接近重尾分布。\n",
"\n",
"接下来对 $\\boldsymbol{t}$ 分布模型做后验检查，并将其与高斯分布对比："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6fab1dec",
"metadata": {},
"outputs": [],
"source": [
"y_ppc_t = pm.sample_posterior_predictive(\n",
" trace_t, 100, model_t, random_seed=123)\n",
"y_pred_t = az.from_pymc3(trace=trace_t, posterior_predictive=y_ppc_t)\n",
"az.plot_ppc(y_pred_t, figsize=(12, 6), mean=False)\n",
"ax[0].legend(fontsize=15)\n",
"plt.xlim(40, 70)"
]
},
{
"cell_type": "markdown",
"id": "303f4e8d",
"metadata": {},
"source": [
"\n",
"\n",
"可以看到，使用 $\\boldsymbol{t}$ 分布之后，从分布的峰值和形状来看，模型的预测值与观测数据更吻合了。留意预测值远离观测值中心的部分，那是因为 $\\boldsymbol{t}$ 分布希望看到在偏离数据中心的两个方向上都有数据。\n",
"\n",
"在新模型中 $\\boldsymbol{t}$ 分布的估计值更稳健，因为异常点降低了自由度 $\\nu$ 的值，而不是将均值拉向异常点方向并且增加方差。也就是说，均值和尺度的估计给占多数的数据点的权值要大于异常点（再次强调：`尺度 scale`不是标准差）。同样要记住，虽然尺度不是标准差，但是它确实与数据的分散程度有关，其值越低，则分布越集中。此外，对应自由度 $\\nu$ 大于等于 2 的情况，尺度值倾向于接近剔除异常值后的标准差。因此，对于大的 $\\nu$ 值，可以近似将 $\\boldsymbol{t}$ 分布的尺度视为剔除异常值后的数据标准差。\n",
"\n",
"## 2.4 组间比较问题\n",
"\n",
"统计分析中一个常见任务是对不同组进行比较，例如：想知道病人对某种药的反应如何、引入某种交通法规后车祸数量是否会降低、学生对不同教学方式的表现如何等。\n",
"\n",
"传统频率主义通常将此类问题统一归入`假设检验`的框架下，其目的是得到统计学意义上的显著性。但仅依赖统计显著性可能会带来很多问题：一方面，统计显著性并非实际显著性；另一方面，即使作用不明显，只要收集尽可能多数据，假设都会被看做具有显著性，亦即显著性被人为放大了。此外，统计显著性的核心思想往往需要计算 `p 值` ，但其与显著性之间并非本质联系，而是一种文化联系，人们习惯于这样做，并非应当这样做，而是因为在大多数统计学入门课程中要求这么做。已有很多文章和研究记录表明，`p 值` 会被错误地使用和解释，即使每天与统计打交道的科学家也难以避免。\n",
"\n",
"在贝叶斯框架下，不做假设检验，而是走一条不同的技术路线。贝叶斯方法将重点放在估计 **效应大小 （`effect size`）** 上，也就是量化两组之间的差异。从效果大小角度思考的好处是，我们不会去回答是或否的问题，比如：它起作用了吗？有什么效果吗？或者更细微的问题，比如：它工作得有多好？影响有多大？\n",
"\n",
"```{note}\n",
" 效应大小是量化两个组之间差异大小的一种方法。\n",
"```\n",
"\n",
"在比较不同组的数据时，人们往往会将其分为一个实验组和一个对照组（也可能超过一个），例如当测试一个新药时，希望将使用新药的组（实验组）和不使用新药的组（对照组）进行对比，以便知道对于治疗某种疾病，使用药物对比不用药物（通常使用安慰剂）的作用有多大。此外，还有另一个有趣的问题是：与治疗某种疾病最常用的某种药相比，新药效果如何？此时，对照组不再是使用安慰剂的组，而是使用其他药物的组。\n",
"\n",
"从统计学上讲，使用假对照组是一种不错的撒谎方式，例如：假设某家乳制品公司想要卖某种含过量糖分的酸奶给小朋友，同时告诉家长乳酸有利于免疫系统。为证明该说法，他们可以和牛奶或者水做对照，而不是选择更便宜、糖更少、市场更小的其他酸奶品种。人们的关注点有时会趋向于主题，而忽略了参照物的选择。因此，下次再听到有人说某种东西更好、更快、更强时，记得问一下他对比的基准是什么。\n",
"\n",
"要在组间做比较，必须决定使用哪些特征进行比较。一个常见特征是每组的平均值。**与频率主义不同，贝叶斯方法致力于获得组间均值差异的后验分布，而不仅仅获得均值差异的点估计。**\n",
"\n",
"为帮助解释后验，会使用三个工具：\n",
"\n",
"- 绘制带有参照值的后验图\n",
"- 计算 `Cohen's d` 值\n",
"- 计算`优势概率（The probability of superiority）`\n",
"\n",
"在之前的案例中，已经看到了如何使用 `az.plot_posterior` 函数绘制带参考值的后验。主要的新概念是 `Cohen's d` 和 `优势概率` 。\n",
"\n",
"\n",
"### 2.4.1 Cohen's d\n",
"\n",
"`Cohen's d` 是一种用来描述组间差异性大小（效应大小）的常见方式：\n",
"\n",
"```{math}\n",
"\\text{Cohen's d} = \\frac{\\mu_{2}-\\mu_{1}}{\\sqrt{\\frac{\\sigma_{2}^{2}+\\sigma_{1}^{2}}{2}}} \\tag{式2.4} \\label{式2.4} \n",
"```\n",
"\n",
"根据该表达式，效应大小是在合并两组标准差的情况下，各组均值相对于合并标准差的差异。由于在贝叶斯框架下，两个组的均值和标准差经过推断后，均表现为后验分布，因此，`Cohen's d` 也表现为概率分布，而不是某个确切值。如果只需要或只想要一个确切值，可以计算 `Cohen's d` 后验的平均值，得到一个 `Cohen's d Value` 。需要说明的是，通常在计算合并标准差时，会显式考虑每组的样本量，但式（ 2.4 ）省略了样本量，主要是因为从后验得到标准差中，已经体现了其不确定性。\n",
"\n",
"`Cohen‘s d` 通过使用各组标准差来引入每组的差异性。这点很重要，因为在相同绝对差的情况下，标准差为 $0.1$ 的组显然应当比标准差为 $10$ 的组效应变化更明显。一组数据相比另一组数据变化了 $x$ 个单位，可能是每个点都整体变化了 $x$ 个单位，也可能是其中一半的数据没有变化而另外一半数据变化了 $2x$ ，还可能是其他组合方法。因此，将组内方差包含在公式中，将有助于将绝对差放在合理上下文中做比较。这种类似对绝对差做标准化的方法，有助于我们正确理解不同组之间差异的重要性。\n",
"\n",
"```{tip}\n",
"根据 `Cohen's d` 得到的效应值可以看做是 `Z-score`，因而 `Cohen's d` 为 $0.5$ 可以解释为一组数据相比另一组数据的差别是 $0.5$ 倍的标准差。\n",
"```\n",
"即使将均值的差异做了标准化，我们可能仍然需要基于给定问题的上下文来校准自己，以便能够判断给定值是大、中、小等等。还好，这种校准可以通过足够的练习获得。举例来说，如果我们对相同类型的问题执行多个分析，如果 `Cohen‘s d =1` 是个经验值，那么当分析遇到 `Cohen's d =2 ` 时，意味着有重要的事情发生。如果您还没有这样的经验，可以向[领域专家](http://rpsychologist.com/d3/cohend)咨询。该网页可以用于了解 不同 `Cohen's d` 值都长什么样，此外，还可以看到描述效应值的其他方式。\n",
"\n",
"### 2.4.2 优势概率\n",
"\n",
"优势概率是表示效应值的另一种方式，描述的是从一组数据中取出的一个点大于从另外一组中取出的点的概率。假设两个组中数据的分布都是高斯分布，我们可以通过以下表达式从 `Cohen's d` 中得到优势概率：\n",
"\n",
"```{math}\n",
"ps=\\Phi\\left(\\frac{\\delta}{\\sqrt{2}}\\right) \\tag{式2.5} \\label{式2.5} \n",
"```\n",
"\n",
"此处， $\\Phi $ 是累积高斯分布， $\\delta $ 是 `Cohen's d`。我们可以计算优势概率的点估计，也可以计算值的整个后验分布。如果同意高斯假设，可以使用该公式从 `Cohen's d` 中计算得到优势概率。如果没有封闭形式，则当有后验样本时，可通过后验样本计算它。这是马尔可夫链蒙特卡罗 (`MCMC`) 方法的一个优点：一旦从后验获得样本，就可以从它计算出很多量来。\n",
"\n",
"### 2.4.3 小费数据集\n",
"\n",
"接下来，使用 `seaborn` 中的 `tips` 数据集来讨论前面提到的想法。\n",
"\n",
"我们希望知道星期几对于餐馆小费数量的影响。该例中实际并没有明确的实验组和对照组之分，只是观察性的实验。如果愿意，可以任选一天（例如星期四）作为实验组或对照组。需注意的一点是：对于观察性实验，仅能得到相关性，无法得出因果关系。\n",
"\n",
"事实上，如何从数据中得出因果关系是非常热门的研究问题，我们会在第 4 章重新讨论该问题。现在，首先用一行代码将数据导入成 `Pandas` 中的 `DataFrame`，`tail` 函数返回数据中的最后一条（当然你也可以用 `head` 函数返回第一条数据）。"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3e438923",
"metadata": {},
"outputs": [],
"source": [
"tips = pd.read_csv('../data/tips.csv')\n",
"tips.tail()"
]
},
{
"cell_type": "markdown",
"id": "460a2a12",
"metadata": {},
"source": [
"\n",
"\n",
"对于该数据集，我们只关心 `day` 和 `tip` 列，利用 `seaborn` 中的 `violinplot` 函数可以将其画出来："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "04d5a81e",
"metadata": {},
"outputs": [],
"source": [
"sns.violinplot(x='day', y='tip', data=tips)t"
]
},
{
"cell_type": "markdown",
"id": "e2a906d6",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"> 图 2.16 不同日期对应的小费分布\n",
"\n",
"\n",
"把问题简化下，我们创建两个变量：变量 `y`表示 `tips`；变量 `idx` 表示分类变量的编码。即用数字 0、1、2、3 表示星期四、星期五、星期六和星期天。 groups 则包含其中每一组的数量。"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "15310e0c",
"metadata": {},
"outputs": [],
"source": [
"tip = tips['tip'].values\n",
"idx = pd.Categorical(tips['day'],\n",
" categories=['Thur', 'Fri', 'Sat', 'Sun']).codes\n",
"groups = len(np.unique(idx))"
]
},
{
"cell_type": "markdown",
"id": "bee60d20",
"metadata": {},
"source": [
"该模型几乎与 `model_g` 相同；唯一区别是， $\\mu$ 和 $\\sigma$ 为向量而不是标量。对于这种情况，`PyMC3` 语法非常有用：可以用向量化方式编写模型，而不用编写 `for` 循环。这意味着对于先验和似然，我们可以使用 `idx` 变量正确地索引到其 `mean` 和 `sd` 变量："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b42bee6e",
"metadata": {},
"outputs": [],
"source": [
"with pm.Model() as comparing_groups:\n",
" μ = pm.Normal('μ', mu=0, sd=10, shape=groups)\n",
" σ = pm.HalfNormal('σ', sd=10, shape=groups)\n",
" y = pm.Normal('y', mu=μ[idx], sd=σ[idx], observed=tip)\n",
" trace_cg = pm.sample(5000)\n",
"az.plot_trace(trace_cg)"
]
},
{
"cell_type": "markdown",
"id": "d95fa2af",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"> 图 2.17 周四到周日小费的均值和方差分布情况（后验）\n",
"\n",
"\n",
"下面代码是在不重复比较情况下绘制差异的一种方式，。与绘制完全比较的矩阵不同，我们在此只绘制了上三角部分："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c030ac42",
"metadata": {},
"outputs": [],
"source": [
"dist = stats.norm()\n",
"_, ax = plt.subplots(3, 2, figsize=(14, 8), constrained_layout=True)\n",
"comparisons = [(i, j) for i in range(4) for j in range(i+1, 4)]\n",
"pos = [(k, l) for k in range(3) for l in (0, 1)]\n",
"for (i, j), (k, l) in zip(comparisons, pos):\n",
" means_diff = trace_cg['μ'][:, i] - trace_cg['μ'][:, j]\n",
" d_cohen = (means_diff / np.sqrt((trace_cg['σ'][:, i]**2 +\n",
"trace_cg['σ'][:, j]**2) / 2)).mean()\n",
" ps = dist.cdf(d_cohen/(2**0.5))\n",
" az.plot_posterior(means_diff, ref_val=0, ax=ax[k, l])\n",
" ax[k, l].set_title(f'$\\mu_{i}-\\mu_{j}$')\n",
" ax[k, l].plot(\n",
" 0, label=f\"Cohen's d = {d_cohen:.2f}\\nProb sup = {ps:.2f}\",\n",
"alpha=0)\n",
" ax[k, l].legend()"
]
},
{
"cell_type": "markdown",
"id": "43510983",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"> 图 2.18 不同日期之间的效应大小\n",
"\n",
"\n",
"\n",
"一种解释该结果的方式是将参考值与 `HPDI` 区间进行比较。其中，只有星期四与星期天的组间对比情况下， `94% HPDI` 没有包含参考值 $0$ ；其他的所有情况，均没法得出组件差别显著大于 $0$ 的结论（根据 `HPDI` 与参考值的重叠性准则）。但是即便如此，平均下来 $0.5$ 美元的小费差别是否足够大了呢？该差别是否大到让人们牺牲星期日陪家人或者朋友的时间去工作呢？是否大到就应该这 4 天都给相同的小费而且男服务员和女服务员的小费一模一样呢？诸如此类的问题很难用统计学来回答，只能从统计学中找到些启发。\n",
"\n",
"## 2.5 分层模型\n",
"\n",
"### 2.5.1 什么是分层模型？\n",
"假设想要分析一个城市的水质，然后将城市分成了多个相邻（或者水文学上）的区域。可以自然想到用如下两种方法进行分析：\n",
"\n",
"- 分别对每个区域单独进行估计；\n",
"\n",
"- 将所有数据都混合在一起，把整个城市看做一个整体进行估计。\n",
"\n",
"两种方式都合理，具体使用哪种取决于想知道什么。如果想了解具体细节，可以采用第一种方式。采用第二种方式则可以将数据都聚在一起，得到一个更大的样本集，从而得出更准确的估计，但是会掩盖细节的信息。两种方式都有合理性，不过或许还可以找到些中间方案。\n",
"\n",
"我们可以在对相邻区域的水质进行评估同时，对整个城市的水质进行评估。这类模型被称为 `分层模型（hierarchical model）` 或者 \n",
"`多层模型(multilevel model)`，名称来自于对数据采用了一种层次化的建模方式。\n",
"\n",
"那么如何构建分层模型呢？\n",
"\n",
"简单说，就是在先验之上使用一个共享先验。也就是说，我们不再固定先验的参数值，而是从数据中将参数值估计出来。这类更高层的先验通常被称为 `超先验（hyper-prior）` ，其参数称为 `超参数`。 \"超\"在希腊语中是\"在某某之上\"的意思。当然，还可以在超先验之上再增加先验，做到尽可能分层。问题是这么做会使模型变得相当复杂而且难以理解。**原则上，除非问题确实需要更复杂的结构，增加分层对于做推断并没有更大帮助，相反，会使分析人员陷入 `超参数` 和 `超先验` 的混乱中而无法做出有意义的解释，并降低模型的可解释性**。\n",
"\n",
"为了更好地解释分层模型中的主要概念，以本节开头提到的水质模型为例，使用一些人工数据来讲解。假设我们从同一个城市的 $3$ 个不同水域得到了含铅量的采样值：其中高于世界卫生组织标准的值标记为 $0$ ；低于标准的值标记为 $1$ 。（注：此处做了简化，实际中会使用铅含量的连续值，并分成更多组）。\n",
"\n",
"我们通过以下代码合成数据："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ef9d9396",
"metadata": {},
"outputs": [],
"source": [
"N_samples = [30, 30, 30]\n",
"G_samples = [18, 18, 18]\n",
"group_idx = np.repeat(np.arange(len(N_samples)), N_samples)\n",
"data = []\n",
"for i in range(0, len(N_samples)):\n",
" data.extend(np.repeat([1, 0], [G_samples[i], N_samples[i]-G_samples[i]]))"
]
},
{
"cell_type": "markdown",
"id": "346b7043",
"metadata": {},
"source": [
"上述代码模拟了一个实验，其中测量了三组数据，每组包含一定数量样本。我们将每组的样本总数存储在 `N_samples` 中，将其中水质好的样本数量存储在 `G_samples` 中。代码的其余部分用来生成了对应的数据集，其中填满了 $0$ 和 $1$ 。\n",
"\n",
"该模型的数学表示法见式（ 2.6 ），基本与抛硬币模型相同，只是有两个重要特征：\n",
"\n",
"- 模型定义了两个影响贝塔先验的超先验，即贝塔先验的两个参数 $\\alpha,\\beta$ 不再是确切的常数值。\n",
"\n",
"- 模型并没有把超先验直接定义在贝塔分布的参数 $\\alpha$ 和 $\\beta$ 上，而是间接地定义在贝塔分布的均值 $\\mu$ 和精度 $\\kappa$ 上（ 用于描述贝塔分布的另外一种参数形式，其中精度可以粗略地理解为标准差的倒数）， $\\kappa$ 值越大，贝塔分布越集中。\n",
"\n",
"$$\n",
"\\mu \\sim \\operatorname{Beta}\\left(\\alpha_{\\mu}, \\beta_{\\mu}\\right) \n",
"$$\n",
"\n",
"$$\n",
"\\kappa \\sim\\left|\\operatorname{Normal}\\left(0, \\sigma_{\\kappa}\\right)\\right| \n",
"$$\n",
"\n",
"$$\n",
"\\alpha =\\mu * \\kappa \n",
"$$\n",
"\n",
"$$\n",
"\\beta =(1-\\mu) * \\kappa \n",
"$$\n",
"\n",
"$$\n",
"\\theta_{i} \\sim \\operatorname{Beta}\\left(\\alpha_{i}, \\beta_{i}\\right) \n",
"$$\n",
"\n",
"$$\n",
"y_{i} \\sim \\operatorname{Bern}\\left(\\theta_{i}\\right) \\tag{式2.6} \n",
"$$\n",
"\n",
"\n",
"注意，式中使用下标索引 $i$ 来指示模型中不同组的贝塔分布参数具有不同的值。也就是说，并非所有参数都在组间共享值。使用 Kruschke 图可以很明显看到新模型中出现了一个额外级别。\n",
"\n",
"\n",
"\n",
"![](\"https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210504181501_31.webp\") \n",
"\n",
"> **图 2.19 三个区域水质概率模型的 Kruschke 图**\n",
"\n",
"\n",
"注意上图中，使用符号 $=$ 而不是 $∼$ 来定义 $\\alpha_i$ 和 $\\beta_i$ 。这是因为一旦知道 $\\mu$ 和 $\\kappa$ 的值， $\\alpha_i$ 和 $\\beta_i$ 的值确定了。此类变量被称为 `确定性变量` ，与 `随机变量`（如 $\\mu,\\kappa,\\theta$ 等）相对。在 PyMC3 中对两者有比较明确的说明。\n",
"\n",
"模型对均值 $\\mu$ 和精度 $\\kappa$ 的参数化，在数学上最终就是对 $\\alpha$ 和 $\\beta$ 的参数化。那为什么要舍近求远呢？这有两个原因：\n",
"\n",
"- 首先，均值和精度虽然数学上等价，但在数值上更适合采样器。本书将在第 8 章中对其做解释。\n",
"\n",
"- 其次是处于教学目的。这个例子表明模型的表达方式可能不止一种，但在数学上等价时在实现上可能会有差异。采样器的效率是一个方面，而另一个方面可能是模型的可解释性。对于某些特定问题或特定受众，报告贝塔分布的均值和精度比报告 $\\alpha$ 和 $\\beta$ 参数更好。\n",
"\n",
"让我们在 `PyMC3` 中实现并求解该模型，这样我们就可以继续讨论层次模型："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c82ac3aa",
"metadata": {},
"outputs": [],
"source": [
"with pm.Model() as model_h:\n",
" μ = pm.Beta('μ', 1., 1.)\n",
" κ = pm.HalfNormal('κ', 10)\n",
" θ = pm.Beta('θ', alpha=μ*κ, beta=(1.0-μ)*κ, shape=len(N_samples))\n",
" y = pm.Bernoulli('y', p=θ[group_idx], observed=data)\n",
" \n",
" trace_h = pm.sample(2000)\n",
"az.plot_trace(trace_h)"
]
},
{
"cell_type": "markdown",
"id": "6eee8822",
"metadata": {},
"source": [
"\n",
"\n",
"![](\"https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210504181604_78.webp\") \n",
"\n",
"> **图 2.20 三个区域水质的分层模型推断结果**\n",
"\n",
"\n",
"### 2.5.2 收缩\n",
"\n",
"现在做个简单实验，使用 `az.summary(trace_h)` 打印并保存摘要数据。然后，希望对合成数据进行微调后再重新运行模型两次，并始终保留摘要记录。也就是总共运行三次：\n",
"\n",
"- 第一次运行将 `G_Samples` 的所有元素值设置为 18\n",
"- 第二次运行将 `G_Samples` 的所有元素值设置为 3\n",
"- 第三次运行将 `G_Samples` 的其中一个元素设置为 18，其他元素设置为 3\n",
"\n",
"继续之前先想想该实验会是什么结果。重点关注每次实验中 $\\theta_i$ 的均值。根据前两次模型的运行结果，能猜出第 3 种情况下的结果吗？如果将结果汇总在表格中，会得到类似如下的值（注意：由于 `NUTS` 采样方法的随机性，结果可能会有小幅波动）：\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"表中第一行，可以看到对于 30 个样本中有 18 个正样本的情况，估计值 $\\theta$ 的均值为 0.6 ；注意 $θ$ 是一个向量，因为有 3 组实验，每组都有一个均值。\n",
"\n",
"表中第二行，可以看到对于 30 个样本中有 3 个正样本的情况，估计值 $\\theta$ 的均值为 0.11 。\n",
"\n",
"表中第三行，结果有点意外， $\\theta$ 的均值并非前面两组中均值的组合（比如 `0.6, 0.11, 0.11`），而是（`0.53, 0.14, 0.14`）。\n",
"\n",
"为什么呢？是模型收敛问题还是模型选型出了问题？其实都不是，这是因为估计结果趋向于整体的均值。事实上，这正是模型预期的结果：在设置了超先验后，直接从数据中估计贝塔先验，每个组的估计值会都受到其他组估计值的影响，同时也影响着其他组的估计值。换句话说，所有组都通过超先验共享了部分信息，从而出现一种称为 `收缩（shrinkage`）的现象。 收缩现象的效果相当于对数据做了 `部分池化（pooling）` ，既不是对数据分组的建模，也不是将数据作为一个大组的建模，而是介于二者之间。\n",
"\n",
"`收缩` 有利于更稳健的推断。这和前面讨论过的学生 $\\boldsymbol{t}$ 分布与异常点的关系很像，那个例子中，使用重尾分布后的模型对偏离均值的异常点表现得更健壮。引入超先验后，在更高层次上进行推断，从而得到一个更 `\"保守\"` 的模型，更少地受到组内极限值的影响。举例来说，假设在某个相邻区域得到了一组不同数量的采样值；那么采样数量越小就越容易得到错误结果；极限情况下，假设在该区域只有一个采样值，而你恰好是从该区域的某个铅管中得到的采样值（或相反，恰好从 PVC 管道中得到的采样值），那么可能导致这片区域水质的整体高估或者低估。在分层模型中，估计出错的情况可通过其他组提供的信息进行改善。当然，通过增加采样值能达到类似效果，但考虑到样本采集成本问题，大多数情况下这并不是最佳候选方案。\n",
"\n",
"显然，收缩的程度取决于数据，数量大的组会对数量较小的组造成更大的影响。如果大多数组都比较相似，而其中某组不太一样，相似组之间会共享这种相似性，从而强化共同的估计值，并拉近表现不太一样的那一组的估计值，前面的例子中已经体现了这一点。\n",
"\n",
"此外，超先验对收缩的程度有影响。如果对所有组的整体分布有一些可信的先验信息，那么将其加入到模型中将有助于收缩程度调整到比较合理的值。\n",
"\n",
"理论上完全可以只用两个组来构建分层模型，不过通常更倾向于使用多个组。直观原因是，收缩其实是将每个组看成一个数据点，而我们是在组这个层级上估计标准差，除非我们对估计值有很强的先验，否则通常不会太相信点数较少的估计值，。\n",
"\n",
"你可能对估计到的先验分布比较感兴趣，以下是将其表示出来的一种方式："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "61b383b0",
"metadata": {},
"outputs": [],
"source": [
"x = np.linspace(0, 1, 100)\n",
"for i in np.random.randint(0, len(trace_h), size=100):\n",
" u = trace_h['μ'][i]\n",
" k = trace_h['κ'][i]\n",
" pdf = stats.beta(u*k, (1.0-u)*k).pdf(x)\n",
" plt.plot(x, pdf, 'C1', alpha=0.2)\n",
"\n",
"u_mean = trace_h['μ'].mean()\n",
"k_mean = trace_h['κ'].mean()\n",
"dist = stats.beta(u_mean*k_mean, (1.0-u_mean)*k_mean)\n",
"pdf = dist.pdf(x)\n",
"mode = x[np.argmax(pdf)]\n",
"mean = dist.moment(1)\n",
"plt.plot(x, pdf, lw=3, label=f'mode = {mode:.2f}\\nmean = {mean:.2f}')\n",
"plt.yticks([])\n",
"plt.legend()\n",
"plt.xlabel(' $ θ_{prior} $ ')\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"id": "16dd3417",
"metadata": {},
"source": [
"\n",
"\n",
"![](\"https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210504181857_6e.webp\") \n",
"\n",
"图2.21\n",
"\n",
"\n",
"### 2.5.3 另一个例子\n",
"\n",
"我们再一次以化学漂移数据为例。这些数据来自一组蛋白质分子。准确地说，该化学漂移数据来自蛋白质的 $^{13}C_\\alpha$ 原子核，因为这是我们能够观测的有限类型原子核。蛋白质是由若干基础类型的 `氨基酸`（总共存在20种基础类型） 序列组合而成。每种氨基酸在序列中出现零次或多次，序列可由几个、数百个、数千个氨基酸组成。每种氨基酸都有且只有一种 $^{13}C_\\alpha$ ，所以可以很有把握地将每一种化学漂移与蛋白质中特定氨基酸联系在一起。此外，这 20 种氨基酸中的每一种都有不同化学性质，对蛋白质的生物学特性有贡献。有些带净电荷，有些是中性的；有些喜欢被水分子包围，而另一些喜欢与同类型氨基酸作伴等。关键点在于：这些氨基酸相似但不相等，因此，根据氨基酸类型的定义，将任何与化学漂移相关的推论建立在 20 个基团上可能是合理甚至自然的。你可以通过这个精彩的视频了解更多关于蛋白质的知识：https://www.youtube.com/watch?v=wvTv8TqWC48 。\n",
"\n",
"在下面的代码块中，我们将数据加载到 `Pandas` 的 `DataFrame` 中。可以看到四列：第一列是蛋白质的 ID ；第二列包括氨基酸的名称，使用标准的三字母代码；而其余列对应于理论计算的化学漂移和实验测量的化学漂移。这个例子的目的是比较理论计算值与实验测量值间的差异以了解其吻合程度即差异原因。出于这个原因，我们计算了系列差异 `diff`："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e4a94aad",
"metadata": {},
"outputs": [],
"source": [
"cs_data = pd.read_csv('../data/chemical_shifts_theo_exp.csv')\n",
"diff = cs_data.theo.values - cs_data.exp.values\n",
"idx = pd.Categorical(cs_data['aa']).codes\n",
"groups = len(np.unique(idx))"
]
},
{
"cell_type": "markdown",
"id": "c5916999",
"metadata": {},
"source": [
"为了解分层模型和非分层模型间的区别，我们构建两个模型。第一个与之前 `comparing_groups` 模型基本相同："
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "16f4ec5b",
"metadata": {},
"outputs": [],
"source": [
"with pm.Model() as cs_nh:\n",
" μ = pm.Normal('μ', mu=0, sd=10, shape=groups)\n",
" σ = pm.HalfNormal('σ', sd=10, shape=groups)\n",
" y = pm.Normal('y', mu=μ[idx], sd=σ[idx], observed=diff)\n",
" trace_cs_nh = pm.sample(1000)"
]
},
{
"cell_type": "markdown",
"id": "36c34490",
"metadata": {},
"source": [
"现在，构建模型的分层版本。我们添加两个超先验：一个用于 $\\mu$ 的均值，另一个用于 $\\mu$ 的标准差（处于教学目的，此处未对 $\\sigma$ 设置超先验）。"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7fa6ea49",
"metadata": {},
"outputs": [],
"source": [
"with pm.Model() as cs_h:\n",
" # hyper_priors\n",
" μ_μ = pm.Normal('μ_μ', mu=0, sd=10)\n",
" σ_μ = pm.HalfNormal('σ_μ', 10)\n",
" # priors\n",
" μ = pm.Normal('μ', mu=μ_μ, sd=σ_μ, shape=groups)\n",
" σ = pm.HalfNormal('σ', sd=10, shape=groups)\n",
" y = pm.Normal('y', mu=μ[idx], sd=σ[idx], observed=diff)\n",
" trace_cs_h = pm.sample(1000)"
]
},
{
"cell_type": "markdown",
"id": "093517f1",
"metadata": {},
"source": [
"我们使用 `ArviZ` 的 `plot_forest` 函数比较结果。当想要比较不同模型（如本例）的参数值时，该函数很有用。注意，此处将几个参数传递给 `plot_forest` 以获得想要的图，例如：设置 `combined=True` 以合并所有链的结果。"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1c3d2b84",
"metadata": {},
"outputs": [],
"source": [
"_, axes = az.plot_forest([trace_cs_nh, trace_cs_h],\n",
" model_names=['n_h', 'h'],\n",
" var_names='μ', combined=False, colors='cycle')\n",
"y_lims = axes[0].get_ylim()\n",
"axes[0].vlines(trace_cs_h['μ_μ'].mean(), *y_lims)"
]
},
{
"cell_type": "markdown",
"id": "6b20897e",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"图 2.22\n",
"\n",
"\n",
"OK，在图 2.22 中能看到什么呢？ 我们有一个 40 个估计平均值的曲线图（共 20 个氨基酸 ，2 种模型），也有 94% HDP可信区间和二分位数范围（分布中心周边 50% )。垂直的黑色线是根据分层模型得到的全局平均值。该值接近于零，表明理论值能够如预期的一样复现实验值。\n",
"\n",
"此图最相关的部分是，分层模型的估计值被拉向 `部分池化` 的均值，或者说，相对于 `未池化的估计`，分层模型的估计值被收缩了。您还会注意到，这种影响对于远离均值的组（如 13) 来说更明显，并且其不确定性比非分层模型更小。这些估计是部分池化的，因为我们对每组都有一个估计，但单个组的估计通过超先验相互制约。因此，我们得到了一种中间情况，一种是只有一个基团，所有的化学漂移都在一起，另一个是 20 个独立基团，每个氨基酸一个。这就是，层次分明的模型之美。\n",
"\n",
"可以肯定地说，分层模型是一个非常棒的想法。在接下来的章节中，我们将继续构建分层模型，并学习如何使用它构建更好的模型。在第 5 章“模型比较”中，我们还将讨论分层模型与统计学和机器学习中普遍存在的 `过拟合/欠拟合` 问题之间的相关性。在第 8 章“推断引擎”中，我们将讨论从分层模型中采样时可能会出现的一些技术问题，以及如何诊断和修复这些问题。\n",
"\n",
"## 2.6 总结\n",
"\n",
"虽然贝叶斯统计概念上很简单，但全概率公式经常无法获得解析解。多年来，这个问题严重阻碍了贝叶斯方法的广泛应用。幸运的是，数学、统计学、物理学和计算机科学以数值方法的方式拯救了人类，这些方法至少在原理上能够解决任何推断问题。自动化推断过程的这种可能性导致了概率编程语言的发展，允许模型定义和推断过程之间的明确分离。\n",
"\n",
"`PyMC3` 是一个 `Python` 库，用于概率编程，具有非常简单、直观和易于阅读的语法，也接近于描述概率模型的统计语法。本章回顾第 1 章“概率思维”中的抛硬币模型，此次没有解析地推导后验，而是采用了 `PyMC3` 构建模型。在 `PyMC3` 中，只需编写一行代码即可将概率分布添加到模型中。概率分布可以组合，可用于先验（未观测变量）或似然（观测变量）。如果将数据传递给一个分布，那它就成为一种似然。`PyMC3` 允许我们从后验分布中获得样本， 而且可以用一行代码就能实现。如果一切正常，这些样本将代表正确的后验分布，从而代表模型和数据的逻辑结果。\n",
"\n",
"可以使用 `ArviZ` 探索 `PyMC3` 生成的后验分布，`ArviZ` 是一个 `Python` 库，它与 `PyMC3` 协同工作，可以用来帮助我们解释和可视化后验分布。一种使用后验来帮助我们做出推断性决策的方法是将 `ROPE` 与 `HPD` 间隔进行比较。我们还简要提到了 `损失函数` 的概念，这是一种在存在不确定性的情况下，量化决策权衡和成本的形式化方法。我们也了解到 `损失函数` 和 `点估计` 是密切相关的。\n",
"\n",
"到目前为止，讨论仅限于简单的单参数模型。在 `PyMC3` 中，将参数推广到任意数量是非常简单的，我们举例说明了如何使用高斯模型和学生 $\\boldsymbol{t}$ 模型来实现这一点。高斯分布是学生 $\\boldsymbol{t}$ 分布的特例，我们展示了在存在异常值的情况下，如何使用后者进行更为稳健的推断。在下一章中，我们将研究如何将这些模型用作线性回归模型的一部分。\n",
"\n",
"我们使用高斯模型来比较组间常见的数据分析任务。虽然人们经常被假设检验背景限制，但贝叶斯方法采取了另一种路线，将组间比较任务转换为推断效应大小的问题，这是一种更丰富、更有成效的方法。我们还探索了解释和报告效应大小的不同方式。\n",
"\n",
"我们把从本书中最重要的概念之一 -- **分层模型** -- 留到了最后。我们可以在每次识别数据中的子组时建立分层模型。在这种情况下，我们可以构建一个模型在组间做部分池化，而不是将子组作为独立对象对待，或者忽略子组将整体作为一个对象来对待。这种部分池化的主要影响是，每个子组的估计值将被其余子组估计值拉偏，这种效果被称为 `收缩`。 收缩是一个非常有用的技巧，可以通过使推断更保守和通过更多信息来帮助改进推断。在接下来的章节中，我们将看到更多分层模型的示例。每个例子都有助于从不同角度更好地理解它们。\n",
"\n",
"\n",
"## 2.7 练习\n",
"\n",
"（1）对于本章的第一个模型，将高斯分布的先验均值修改为一个经验性均值，用几个对应的标准差多跑几遍，观察推断过程对这些变化的鲁棒性/敏感性如何。你觉得用一个没有限制上下界的高斯分布对有上下界的数据建模的效果怎样？记住我们说过数据不可能大于 100 或者小于 0。\n",
"\n",
"（2）利用第一个例子中的数据，分别在包含和不包含异常值情况下，计算出经验均值和标准差。将结果与使用高斯分布和 $\\boldsymbol{t}$ 分布的贝叶斯估计进行比较，增加更多异常值并重复该过程。\n",
"\n",
"（3）修改小费例子中的模型，使其对于异常点更鲁棒。分别尝 试对所有组使用一个共享的 $\\nu$ 和单独为每个组设置一个 $\\nu$ ，最后对这 3 个模型进行后验预测检查。\n",
"\n",
"（4）直接从后验中计算出优势概率（先不要计算 `Cohen's d`），你可以用 `sample_ppc()` 函数从每个组中获取一个采样值。对比这样做与基于高斯假设的计算是否不同？并对结果做出解释。\n",
"\n",
"（5）重复水质对比的例子，不过不用多层模型，而是使用一个均匀分布（比如 ）。比较两种模型的结果。\n",
"\n",
"（6）在小费的例子中，对一个星期中的不同天使用部分池化操作，构建一个多层模型，将结果与不使用多层模型的结果进行对比。\n",
"\n",
"（7）重复本章中的所有例子，用 `findMAP()` 函数的返回值来初始化采样。看是否能得到相同的推断结果。同时看一下 `find_MAP ()` 函数对退化过程的数量以及推断的速度有什么影响。\n",
"（8）对所有模型进行诊断测试并采取相应措施，比如，如果有 必要，增加采样次数。\n",
"\n",
"（9）对本章中的至少一个模型使用你自己的数据并运行。牢记第 1 章中提到的构建模型的 3 个步骤。\n",
"\n",
"\\begin{align*} \\tag{式2.7} \\label{式2.7} \n",
"\\alpha &\\sim \\mathcal{N}(\\mu_{\\alpha},\\sigma_{\\alpha})\\\\\n",
"\\beta &\\sim \\mathcal{N}(\\mu_{\\beta},\\sigma_{\\beta})\\\\\n",
"\\epsilon &\\sim |\\mathcal{N}(0,\\sigma_{\\epsilon})|\\\\\n",
"\\nu &\\sim \\text{Exp}(\\lambda)\\\\\\\\\n",
"y &\\sim \\mathcal{T}(\\alpha + x \\beta, \\epsilon,\\nu)\n",
"\\end{align}\n",
"\n",
"\n",
"\\begin{align*} \\tag{式2.8} \\label{式2.8} \n",
"\\beta_0 &\\sim \\mathcal{N}(\\mu_{\\beta_0},\\sigma_{\\beta_0})\\\\\n",
"\\beta_1 &\\sim \\mathcal{N}(\\mu_{\\beta_1},\\sigma_{\\beta_1})\\\\\n",
"&\\text{...}\\\\\n",
"\\beta_n &\\sim \\mathcal{N}(\\mu_{\\beta_n},\\sigma_{\\beta_n})\\\\\n",
"\\epsilon &\\sim |\\mathcal{N}(0,\\sigma_{\\epsilon})|\\\\\\\\\n",
"y &\\sim \\mathcal{N}(\\beta_0 + \\beta_1x^1 + \\beta_2x^2 + ... + \\beta_nx^n,\\epsilon)\n",
"\\end{align*}\n",
"\n",
"---\n",
"\n",
"\\begin{align*} \\tag{式2.9} \\label{式2.9} \n",
"\\beta_0 &\\sim \\mathcal{N}(\\mu_{\\beta_0},\\sigma_{\\beta_0})\\\\\n",
"\\beta_1 &\\sim \\mathcal{N}(\\mu_{\\beta_1},\\sigma_{\\beta_1})\\\\\n",
"&\\text{...}\\\\\n",
"\\beta_n &\\sim \\mathcal{N}(\\mu_{\\beta_n},\\sigma_{\\beta_n})\\\\\n",
"\\epsilon &\\sim |\\mathcal{N}(0,\\sigma_{\\epsilon})|\\\\\\\\\n",
"y &\\sim \\mathcal{N}(\\beta_0 + \\beta_1x_1 + \\beta_2x_2 + ... + \\beta_nx_n,\\epsilon)\n",
"\\end{align*}\n",
"\n",
"\n",
"---\n",
"\n",
"\\begin{align*} \\tag{式2.10} \\label{式2.10} \n",
"\\mu_{\\alpha} &\\sim \\mathcal{N}(\\mu_{\\mu_{\\alpha}},\\sigma_{\\mu_{\\alpha}})\\\\\n",
"\\sigma_{\\alpha} &\\sim |\\mathcal{N}(0,\\sigma_{\\sigma_{\\alpha}})|\\\\\n",
"\\end{align*}\n",
"\n",
"---\n",
"\n",
"\\begin{align*} \\tag{式2.11} \\label{式2.11} \n",
"\\mu_{\\beta} &\\sim \\mathcal{N}(\\mu_{\\mu_{\\beta}},\\sigma_{\\mu_{\\beta}})\\\\\n",
"\\sigma_{\\beta} &\\sim |\\mathcal{N}(0,\\sigma_{\\sigma_{\\beta}})|\\\\\n",
"\\alpha &\\sim \\mathcal{N}(\\mu_{\\alpha},\\sigma_{\\alpha})\\\\\n",
"\\beta &\\sim \\mathcal{N}(\\mu_{\\beta},\\sigma_{\\beta})\\\\\n",
"\\epsilon &\\sim |\\mathcal{N}(0,\\sigma_{\\epsilon})|\\\\\n",
"\\nu &\\sim \\text{Exp}(\\lambda)\\\\\n",
"y &\\sim \\mathcal{T}(\\alpha + x \\beta,\\epsilon,\\nu)\n",
"\\end{align*}"
]
}
],
"metadata": {
"jupytext": {
"formats": "ipynb,md:myst",
"text_representation": {
"extension": ".md",
"format_name": "myst",
"format_version": 0.13,
"jupytext_version": "1.11.5"
}
},
"kernelspec": {
"display_name": "Python 3",
"language": "python3",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.12"
},
"source_map": [
13,
68,
80,
85,
102,
107,
147,
149,
163,
165,
179,
181,
213,
218,
226,
228,
262,
264,
270,
272,
290,
303,
317,
332,
352,
356,
385,
392,
398,
400,
410,
412,
422,
424,
428,
432,
476,
478,
482,
493,
526,
535,
546,
548,
560,
567,
637,
640,
646,
648,
659,
664,
668,
675,
686,
702,
736,
743,
797,
806,
847,
866,
881,
886,
890,
896,
900,
910,
914,
920
]
},
"nbformat": 4,
"nbformat_minor": 5
} |