{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3ca15dd5",
   "metadata": {},
   "source": [
    "# 第 5 章 模型比较与模型平均\n",
    "\n",
    "<style>p{text-indent:2em;2}</style>\n",
    "\n",
    "模型应该被设计成帮助我们理解特定问题的近似值，而不是被设计成真实世界的翻版。从这个意义上讲，所有模型都是错误的。即使在有先验的情况下，模型也都是错误的，只不过不同模型的错误可能有所不同，而其中一些比其他模型更好地描述了给定问题。\n",
    "\n",
    "此前的章节将注意力集中在推断问题上，即如何从数据中学习参数的分布。本章将重点讨论一个互补问题：如何对解释同一数据的多个模型进行比较？这是对数据进行分析必须解决的关键问题之一。 \n",
    "\n",
    "**为什么你的模型是比较好的那个？**\n",
    "\n",
    "本章将讨论以下内容：\n",
    "\n",
    "- 后验预测检查\n",
    "- 奥卡姆剃刀\n",
    "- 过拟合和欠拟合\n",
    "- 信息准则\n",
    "- 贝叶斯因子\n",
    "- 正则化先验\n",
    "\n",
    "---\n",
    "\n",
    "## 5.1 后验预测分布\n",
    "\n",
    "`第 1 章`介绍了 **后验预测检查** 的概念，实际上它作为一种度量，可以被用来评估 *根据同一数据拟合出的不同模型* 对该数据的**解释程度** 。\n",
    "\n",
    "如前所述，所有的模型都是错误的，因此后验预测检查的目的并非判定模型是否正确或错误，而是希望通过后验预测检查更好地把握模型的局限性，以便做出适当改进。通常模型不会再现数据中的所有问题，但这本身并不是问题，因为构建模型都是为了特定目的。后验预测检查本身就是在特定目的背景下评估模型的一种方式，因此在考虑多个模型的时候，可以使用后验预测检查来对其进行比较。\n",
    "\n",
    "让我们读取并绘制一个简单的数据集："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "77178000",
   "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<cell line: 6>\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<module>\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<module>\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<module>\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<module>\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<module>\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<module>\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<module>\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<module>\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": "c4e17070",
   "metadata": {},
   "outputs": [],
   "source": [
    "dummy_data = np.loadtxt('../data/dummy.csv')\n",
    "x_1 = dummy_data[:, 0]\n",
    "y_1 = dummy_data[:, 1]\n",
    "order = 2\n",
    "x_1p = np.vstack([x_1**i for i in range(1, order+1)])\n",
    "x_1s = (x_1p - x_1p.mean(axis=1, keepdims=True)) / x_1p.std(axis=1, keepdims=True)\n",
    "y_1s = (y_1 - y_1.mean()) / y_1.std()\n",
    "plt.scatter(x_1s[0], y_1s)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51c46e62",
   "metadata": {},
   "source": [
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/bayesian_stat_20210510111807b3.webp)\n",
    "\n",
    "> **图 5.1 一个简单数据集的散点图**\n",
    "</center>\n",
    "\n",
    "现在用两个不同的模型来拟合该数据，第一个是线性模型 `model_l`，第二个是二阶多项式模型 `model_p` ："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1caa5ab0",
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model_l:\n",
    "    α = pm.Normal('α', mu=0, sd=1)\n",
    "    β = pm.Normal('β', mu=0, sd=10)\n",
    "    ϵ = pm.HalfNormal('ϵ', 5)\n",
    "    μ = α + β * x_1s[0]\n",
    "    y_pred = pm.Normal('y_pred', mu=μ, sd=ϵ, observed=y_1s)\n",
    "    \n",
    "    trace_l = pm.sample(2000)\n",
    "\n",
    "with pm.Model() as model_p:\n",
    "    α = pm.Normal('α', mu=0, sd=1)\n",
    "    β = pm.Normal('β', mu=0, sd=10, shape=order)\n",
    "    ϵ = pm.HalfNormal('ϵ', 5)\n",
    "    μ = α + pm.math.dot(β, x_1s)\n",
    "    y_pred = pm.Normal('y_pred', mu=μ, sd=ϵ, observed=y_1s)\n",
    "    \n",
    "    trace_p = pm.sample(2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61feb5be",
   "metadata": {},
   "source": [
    "现在绘制两个模型的平均拟合曲线："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dede5ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_new = np.linspace(x_1s[0].min(), x_1s[0].max(), 100)\n",
    "α_l_post = trace_l['α'].mean()\n",
    "β_l_post = trace_l['β'].mean(axis=0)\n",
    "y_l_post = α_l_post + β_l_post *x_new\n",
    "plt.plot(x_new, y_l_post, 'C1', label='linear model')\n",
    "\n",
    "α_p_post = trace_p['α'].mean()\n",
    "β_p_post = trace_p['β'].mean(axis=0)\n",
    "idx = np.argsort(x_1s[0])\n",
    "y_p_post = α_p_post + np.dot(β_p_post, x_1s)\n",
    "plt.plot(x_1s[0][idx], y_p_post[idx], 'C2', label=f'model order {order}')\n",
    "\n",
    "α_p_post = trace_p['α'].mean()\n",
    "β_p_post = trace_p['β'].mean(axis=0)\n",
    "x_new_p = np.vstack([x_new**i for i in range(1, order+1)])\n",
    "y_p_post = α_p_post + np.dot(β_p_post, x_new_p)\n",
    "plt.scatter(x_1s[0], y_1s, c='C0', marker='.')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "520aff97",
   "metadata": {},
   "source": [
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/bayesian_stat_20210510112129a4.webp)\n",
    "\n",
    "> **图 5.2  数据集的散点图以及两种拟合结果**\n",
    "</center>\n",
    "\n",
    "图中二阶模型似乎做得更好，但线性模型也并没有那么糟糕。此时可以使用 PyMC3 来获得两个模型的后验预测样本，并执行检查："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bb9d538",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_l = pm.sample_posterior_predictive(trace_l, 2000, model=model_l)['y_pred']\n",
    "y_p = pm.sample_posterior_predictive(trace_p, 2000, model=model_p)['y_pred']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b2adab3",
   "metadata": {},
   "source": [
    "正如已经看到的，后验预测检查通常使用可视化方式来做，如下例所示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cc9d47b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 3))\n",
    "data = [y_1s, y_l, y_p]\n",
    "labels = ['data', 'linear model', 'order 2']\n",
    "for i, d in enumerate(data):\n",
    "    mean = d.mean()\n",
    "    err = np.percentile(d, [25, 75])\n",
    "    plt.errorbar(mean, -i, xerr=[[-err[0]], [err[1]]], fmt='o')\n",
    "    plt.text(mean, -i+0.2, labels[i], ha='center', fontsize=14)\n",
    "    plt.ylim([-i-0.5, 0.5])\n",
    "    plt.yticks([])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66f6fcaf",
   "metadata": {},
   "source": [
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/bayesian_stat_20210510112319ed.webp)\n",
    "\n",
    "**图 5.3 不同模型的均值和四分位数区间**。 圆点表示均值，线的起点和终点分别为 25% 和 75 % 分位点。\n",
    "</center>\n",
    "\n",
    "图 5.3 显示了数据、线性模型和二次多项式模型的均值和四分位数范围。该图对 *各模型的后验预测样本* 做了平均，而且两个模型的均值都复现得很好，分位数范围也不是很差。不过在实际问题中，有一些小差异可能是值得注意的。为此，可以尝试做更多的曲线图来分析后验预测分布，以发现问题。例如，绘制均值和四分位数间相对于数据真实值的离散度。\n",
    "\n",
    "下图是一个示例："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f06ae712",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 3), constrained_layout=True)\n",
    "\n",
    "def iqr(x, a=0):\n",
    "    return np.subtract(*np.percentile(x, [75, 25], axis=a))\n",
    "\n",
    "for idx, func in enumerate([np.mean, iqr]):\n",
    "    T_obs = func(y_1s)\n",
    "    ax[idx].axvline(T_obs, 0, 1, color='k', ls='--')\n",
    "    for d_sim, c in zip([y_l, y_p], ['C1', 'C2']):\n",
    "    T_sim = func(d_sim, 1)\n",
    "    p_value = np.mean(T_sim >= T_obs)\n",
    "    az.plot_kde(T_sim, plot_kwargs={'color': c},\n",
    "    label=f'p-value {p_value:.2f}', ax=ax[idx])\n",
    "    ax[idx].set_title(func.__name__)\n",
    "    ax[idx].set_yticks([])\n",
    "    ax[idx].legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da2cc90a",
   "metadata": {},
   "source": [
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/bayesian_stat_2021051011245094.webp)\n",
    "\n",
    "> **图 5.4 两种模型的均值参数和四分位范围参数的分布图和 p-value 值**\n",
    "</center>\n",
    "\n",
    "图 5.4 中，黑色虚线表示根据真实计算的数据均值和四分位数（来自真实数据，为确切值而非分布）。图中不同颜色的曲线分别表示两种模型的统计量的分布，橙色曲线为线性模型根据后验预测样本计算得出的均值分布（左图）或四分位数范围分布（右图）；绿色曲线表示二次多项式模型的同一统计量的分布。图 5.4 还包括 `p-value` 值，该值来自于样本，是预测数据与实际观测数据的比较和计算。\n",
    "\n",
    "对于两个预测数据集，我们分别计算了其平均值和四分位数范围，然后计算了两个统计量等于或大于根据实际数据统计量的比例。**一般而言，如果真实数据和预测结果一致，预期 `p-value` 值在 0.5 左右，否则将处于有偏的后验预测分布**。\n",
    "\n",
    "贝叶斯 `p-value` 与频率派的 `p-value` 名字相似，定义基本上也相同：\n",
    "\n",
    "$$\n",
    "\\text{Bayesian p-value}\\triangleq p\\left(T\\_{sim} \\geq T_{o b s} \\mid y \\right) \\tag{式 5.1}\n",
    "$$\n",
    "\n",
    "贝叶斯 `p-value` 可以被解释为：从模拟数据中获得与观测数据相同或更高统计量值的概率。式中 $T$ 几乎可以代表数据的任意统计量。在图 5.4 中，统计量是左侧的`均值`和右侧的`四分位数范围`。通常 $T$ 应该在最初定义推断任务时就选择好。\n",
    "\n",
    "贝叶斯 `p-value` 与频率主义的内涵完全不同，因为预测并计算统计量所用到的样本，来自于对后验预测分布的采样，而非来自一个参数已经确定的模型。\n",
    "\n",
    "需要注意：贝叶斯的 `p-value` 不需要频率主义的任何零假设作为条件；事实上，我们拥有基于观测数据的完整后验分布。此外，贝叶斯也没有使用类似置信度的任何预定义阈值来声明统计显著性，也没有执行假设检验。贝叶斯方法只是试图计算一个数字来评估后验预测分布与数据集的拟合度而已。\n",
    "\n",
    "无论使用曲线图还是数据摘要（如贝叶斯 `p-value` ），或是两者组合，后验预测检查都是非常灵活的。该概念可让分析师思考采用不同方法来探索后验预测分布，并使用合适的方法来讲述一个数据驱动的故事，**包括但不限于模型比较**。\n",
    "\n",
    "> **小结：**\n",
    ">\n",
    "> 后验预测分布指在获得 **参数的后验分布** 后，通过采样和边缘化得出的 **结果变量的预测分布** 。由于模型参数是不确定的随机变量，因此根据其得到的预测结果也是不确定的，需要用概率分布形式来表示。这两种分布在概念和内涵上有非常显著的区别，但很容易让初学者困惑。建议初学者在阅读和使用时，突出两者的修饰定语 “参数的” 和 “结果变量的” ，以帮助自己有效区别二者。\n",
    ">\n",
    "> 本节介绍了如何利用统计推断获得的后验分布，来比较不同模型对数据的解释程度。其中涉及的主要环节或手段包括：\n",
    "> \n",
    "> **① 利用后验分布获取预测样本** ：蒙特卡洛方法是用于处理复杂概率分布的常用方法，因此，当后验分布无法获得封闭形式解时，对其采样并做蒙特卡洛积分来计算后验预测样本，就成为一种自然的选择。`pymc3` 提供了 `pm.sample_posterior_predictive()` 函数，对后验抽样并边缘化获得后验预测样本；\n",
    "> \n",
    "> **② 绘制预测曲线和真实数据散点图**：绘制后验预测均值曲线和 `HDPI` 区间（ 因为预测结果是随机量 ），并与实际数据进行可视化的分析比较；\n",
    "> \n",
    "> **③ 计算预测的统计量，并与真实数据做比较**：计算真实数据和后验预测样本的均值、方差、分位数等统计特征，而后进行可视化地分析比较，需清楚真实数据的统计量是确定值，而后验预测样本的统计量是随机量；\n",
    "> \n",
    "> **④ 计算模型的贝叶斯 `p-value` 评价指标**，为不同的模型分别计算一些统计量的贝叶斯 `p-value` 值，用于量化和比较模型对真实数据的解释程度。\n",
    "\n",
    "\n",
    "在接下来几节中，我们探索一些其他模型比较的方法。\n",
    "\n",
    "## 5.2 好模型应当具备的特征\n",
    "\n",
    "假如对同一个问题（或数据）有两个模型，二者对数据解释得同样好，应该选哪个模型呢？有一个基本准则叫做**奥卡姆剃刀**，即**如果对同一现象有两种不同假说，应选用比较简单的那一种**。\n",
    "\n",
    "关于奥卡姆剃刀的论证很多，其中一种说法与波普尔的可证伪性有关，还有一种说法是从实用角度提出的，因为简单模型相比复杂模型更容易理解，此外还有一种论证是基于贝叶斯统计的。这里不深入讨论该准则的论证细节，只将该准则当做一个有用而合理的常识。\n",
    "\n",
    "在比较模型时，既要考虑**简约型**，也需要同时考虑**准确性**（ 即模型对数据拟合得如何 ）。之前章节已出现过一些度量准确性的指标，如： $R^2$ 系数可视为线性回归中可解释方差的比例。但**如果有两个模型，其中一个对数据的解释比另一个更准确，是否应该选更准确率的模型呢？**直觉上，似乎最好选择准确度高且简单的模型。但如果简单模型准确度最差，该怎么办？如何才能平衡这两种要素呢？\n",
    "\n",
    "为简化问题，此处引入一个示例来帮助理解准确性与简约性之间的平衡。为了更形象些，该例使用一系列逐渐复杂的多项式来拟合同一个简单数据集，并且未采用贝叶斯方法，而是采用频率主义的最小二乘估计（ `MSE` ）来建模。\n",
    "\n",
    "> 具有贝叶斯基础的人大部分都知道： `MSE` 可转化为带均匀先验的贝叶斯模型，因此，此处将其理解成贝叶斯方法也没问题。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62c6b88e",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([4., 5., 6., 9., 12, 14.])\n",
    "y = np.array([4.2, 6., 6., 9., 10, 10.])\n",
    "plt.figure(figsize=(10, 5))\n",
    "order = [0, 1, 2, 5]\n",
    "plt.plot(x, y, 'o')\n",
    "for i in order:\n",
    "    x_n = np.linspace(x.min(), x.max(), 100)\n",
    "    coeffs = np.polyfit(x, y, deg=i)\n",
    "    ffit = np.polyval(coeffs, x_n)\n",
    "    p = np.poly1d(coeffs)\n",
    "    yhat = p(x)\n",
    "    ybar = np.mean(y)\n",
    "    ssreg = np.sum((yhat-ybar)**2)\n",
    "    sstot = np.sum((y - ybar)**2)\n",
    "    r2 = ssreg / sstot\n",
    "    plt.plot(x_n, ffit, label=f'order {i}, $R^2$= {r2:.2f}')\n",
    "plt.legend(loc=2)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y', rotation=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d160a04d",
   "metadata": {},
   "source": [
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210512104113_44.webp)\n",
    "\n",
    "图 5.5 \n",
    "</center>\n",
    "\n",
    "### 5.2.1 参数太多会导致过拟合\n",
    "\n",
    "从图 5.5 可看出，模型复杂度增加时，对应的 $R^2$ 系数在上升。当多项式为 5 阶时，模型完美拟合了数据（ $R^2$ 趋近于 $1$ 表示更好地拟合了数据）。前面章节中讨论过，用多项式去解决实际问题并非特别好的办法。为什么 $5$ 阶多项式能完美拟合所有数据呢？原因是模型中参数数量与样本数量相同，都是 $6$ 。**此时模型只是用另一种方式对数据进行了编码，并没有从数据中学到任何内容，只是记住了全部数据而已。** 此外，如果使用这几种模型做预测，$5$ 阶多项式模型对数据的预测看起来也会非常奇怪。\n",
    "\n",
    "假设收集了更多数据点。例如，收集到点 $\\{(10，9)，(7，7)\\}$ （ 参见图 5.5 ）。与 $1$ 阶或 $2$ 阶模型相比，$5$ 阶模型对新点的解释效果如何？不是很好 ！！！  $5$ 阶模型没有在数据中学习到任何有趣的模式，只是记住了一些表面的东西，因此其在泛化到未来数据方面做得非常差：\n",
    "\n",
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210512104753_be.webp)\n",
    "\n",
    "> **图 5.5 相同样本支撑下的不同阶次多项式模型**\n",
    "</center>\n",
    "\n",
    "当一个模型与训练数据集非常吻合，但在测试数据集上却非常差时，被称为 `过拟合` 。过拟合是统计学和机器学习中一个普遍问题。\n",
    "\n",
    "描述过拟合问题的一个有效方法是**将数据集视为由 `信号` 和 `噪声` 两部分组成**。信号是想要从数据中了解到的东西，如果使用某个数据集，那必然是因为我们认为该数据集中有一个信号，否则训练毫无意义；而噪声是数据中无用的部分，往往是测量误差、数据生成方式、数据损坏等因素的产物。当某个模型过于灵活，以至于能够学到噪声而掩盖了信号时，该模型就会变得过拟合。\n",
    "\n",
    "避免过拟合是奥卡姆剃刀的确切理由之一。上例表明，如果仅关注模型对数据的拟合程度，很容易被过拟合误导，因为理论上通过增加模型参数数量总是能够提高数据拟合的准确率。\n",
    "\n",
    "### 5.2.2 参数太少会导致欠拟合\n",
    "\n",
    "继续关注该例，不过重点放在 $0$ 阶模型上。在 $0$ 阶模型中，所有 $\\beta$ 参数都为 $0$ ，因而变量 $x$ 和 $y$ 间的线性关系变成了只描述结果变量的一个高斯模型。对于 $0$ 阶模型来说，预测变量对模型不再有任何影响，模型只捕捉到结果变量的均值。换句话说，模型认为数据能够通过结果变量的均值以及一些高斯噪声来解释。我们称这种模型是欠拟合的，因为它实在太简单了，以至于不能从数据中获取有意义的模式。通常，一个参数很少的模型容易出现欠拟合。\n",
    "\n",
    "### 5.2.3 简约性与准确性之间的平衡\n",
    "\n",
    "经常与奥卡姆剃刀准则一起提到的是爱因斯坦的一句名言 **“事情应该尽可能简单，但不必过于简单”** 。我们在建模时需要保持某种平衡。理想状态下，模型既不过拟合也不欠拟合，因此，通常需要优化或者调整模型来权衡二者。\n",
    "\n",
    "机器学习领域中，通常从`方差（variance）`和 `偏差（bias）` 两个角度来讨论和权衡二者：\n",
    "\n",
    "- 高偏差（ `bias` ）是模型适应数据的能力不足导致的。高偏差可能使模型无法捕捉数据中一些关键模式，导致欠拟合。\n",
    "- 高方差（ `variance`）是模型对数据中细节过于敏感导致的。高方差会使模型捕捉到数据中的噪声，导致过拟合。\n",
    "\n",
    "图 5.5 中，$0$ 阶模型具有较高偏差（和较低的方差），因为它偏向于在变量 $y$ 的平均值处返回一条平坦直线，而与 $x$ 值无关。$5$ 阶模型具有较高的方差（和较低的偏差），你可以采用差别很大的方式设置六个点，会发现曲线将完美拟合其中的大多数点。\n",
    "\n",
    "具有高偏差的模型具有更多偏见或惯性，而具有高方差的模型是思想更开放的模型。太有偏见存在的问题是没有能力容纳新证据；太开放的问题是最终会相信荒唐的东西。总体来说，如果提升其中一个方面，就会导致另外一方面的下降，就像天平的两端，这也是它们被称为 `偏差-方差平衡` 的原因。我们最希望得到二者平衡的模型。如何做到呢？下面是一些经验方法。\n",
    "\n",
    "处理 `方差` 较大的问题：\n",
    "\n",
    "- 减少预测变量的数量\n",
    "\n",
    "- 使用更简单的模型\n",
    "\n",
    "- 增加训练数据集\n",
    "\n",
    "- 使用正则化\n",
    "\n",
    "- 增加噪声（加入随机因子，例如采用 bagging 和 boosting 方法）\n",
    "  \n",
    "处理 `偏差` 较大的问题：\n",
    "\n",
    "- 增加预测变量的数量\n",
    "\n",
    "- 使用更复杂的模型\n",
    "\n",
    "- 去掉正则化\n",
    "\n",
    "## 5.3 频率主义的模型比较\n",
    "\n",
    "在上例中，很容易看出 $0$ 阶模型非常简单，而 $5$ 阶模型相对数据过于复杂，但其他两个模型呢？要回答该问题，需要一种原则性的方式，在考虑准确性同时，兼顾考虑简约性。要做到这一点，需要引入几个新概念：\n",
    "\n",
    "- **样本内精度**：基于拟合模型的样本数据测量得到的模型精度。\n",
    "- **样本外精度**：用拟合模型的样本数据以外的数据测量得到的模型精度（也称为 `预测精度`）。\n",
    "\n",
    "对于数据和模型的任意组合，平均而言，样本内精度将优于样本外精度。通常样本内精确度会让我们对模型过于自信；似乎采用样本外精度比样本内精度更可取，但也存在问题。因此，合理做法是放弃一部分样本数据，仅将其用于对模型的测试。但对大多数情况来说，仅将花大成本得到的数据用于测试，似乎过于奢侈了。为避免该问题，人们花了很多精力来思考**使用样本内数据估计样本外精度的方法**。其中两种常见方法包括：\n",
    "\n",
    "- **交叉验证**：这是一种经验性策略，将数据分为多个子集，并轮流将其中一个子集作为测试集，剩余子集作为训练集进行评估。\n",
    "- **信息准则**：通过一些封闭形式的表达式来量化模型适用程度，可以认为表达式的值能够体现与交叉验证近似的结果。\n",
    "\n",
    "### 5.3.1 交叉验证（ `Cross-validation` ） \n",
    "\n",
    "交叉验证是一种简单有效的解决方案，可在不遗漏数据的情况下评估模型。此过程的示意见下图。通常把数据分成大致相等的 $K$ 份，使用其中 $K-1$ 份训练模型 $A_1$，剩下的 1 份用来测试模型；然后，从训练集中重新选择不同的 $K-1$ 份用于训练模型 $A_2$，并用剩余的 1 份测试模型；如此直到完成所有 $K$ 轮，得到模型 $A_K$；然后对结果 $A$ 求平均。\n",
    "\n",
    "上述交叉验证过程被称为 `K-折交叉验证` 。当 $K$ 与样本数量相同时（即 $ K = N$ 时），就是常称的 `留一法交叉验证 （LOO-CV）`。在执行留一法交叉验证时，如果数据数量太多，有时会出现轮数少于数据总数的情况。\n",
    "\n",
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210512112629_82.webp)\n",
    "\n",
    "图 5.7\n",
    "</center>\n",
    "\n",
    "交叉验证是机器学习从业者的谋生之本，有关更多细节，可以阅读 `Sebastian Raschka` 的《`Python Machine Learning`》一书，或 `Jack Vanderplas` 的《`Python Data Science Handbook`》。\n",
    "\n",
    "交叉验证简单而强大，不过对某些模型或者量很大的数据而言，交叉验证的计算量可能超出可接受范围。因此，许多人尝试提出了一些更容易计算的量，来得到近似交叉验证的效果，或者应用到不能直接使用交叉验证的情况，其中比较出名的是**信息准则**。\n",
    "\n",
    "### 5.3.2 信息准则 \n",
    "\n",
    "信息准则是一系列用于比较模型对数据拟合程度的方法，此类方法引入了一个惩罚项来平衡模型的复杂度。换句话说，信息准则形式化地表达了在本章开始建立的直觉，用一种合适的方式平衡模型的准确性和简约性。这些衡量方式的推导过程与信息论相关，超出了本书范围，我们只从实用的角度去理解这些概念。\n",
    "\n",
    "#### **（ 1 ）`log 似然`与离差**\n",
    "\n",
    "一种衡量模型对数据的拟合程度的方法是计算模型预测结果与真实数据之间的均方差：\n",
    "\n",
    "$$\n",
    "\\frac{1}{N} \\sum_{i=1}^{N}\\left(y_{i}-\\mathrm{E}\\left(y_{i} \\mid \\theta\\right)\\right)^{2}  \\tag{式 5.2} \n",
    "$$\n",
    "\n",
    "其中，$E(y_i|\\theta)$ 是根据估计的参数值计算得到的预测值。\n",
    "\n",
    "可以看到基本上就是观察值和预测值之间平均差值，求平方是为保证误差为正，不会相互抵消。相比其他的度量指标（比如绝对值误差），平方度量更强调较大的误差。\n",
    "\n",
    "一种更通用的方法是计算 `log 似然`：\n",
    "\n",
    "$$\n",
    "\\sum_{i=1}^{N} \\log p\\left(y_{i} \\mid \\theta\\right)  \\tag{式 5.3} \n",
    "$$\n",
    "\n",
    "当似然服从正态分布时，已经证明 `log 似然`与二次均方误差成正比。由于历史原因，实践中人们通常不直接使用 `log 似然`，而是使用一个称作 `离差（deviance）` 的量：\n",
    "\n",
    "$$\n",
    "-2 \\sum_{i=1}^{N} \\log p\\left(y_{i} \\mid \\theta\\right)  \\tag{式 5.4}  \n",
    "$$\n",
    "\n",
    "离差在贝叶斯方法和非贝叶斯方法中类似，区别在于：贝叶斯框架中 $θ$ 来自后验的采样。而在非贝叶斯方法中，$θ$ 是一个点估计。在使用离差时，需注意以下两点：\n",
    "\n",
    "- 离差越小，`log 似然`值越大，模型的预测结果与数据越吻合。因此我们**希望离差越小越好**。\n",
    "- 离差衡量的是样本内的模型精度，因而复杂模型通常会比简单模型的离差小，此时**需要给复杂模型加入惩罚项**。\n",
    "\n",
    "下面我们将学习几个不同的信息准则方法，**它们的共同点是都使用了离差和正则项，区别在于离差和惩罚项的计算方式不同**。\n",
    "\n",
    "#### **（ 2 ）`AIC` 信息准则**\n",
    "\n",
    "`AIC 信息准则`（Akaike Information Criterion）是一个广泛应用的信息准则，其定义如下：\n",
    "\n",
    "$$\n",
    "\\text{AIC} = -2\\sum_{i=1}^{n} \\log p\\left(y_{i} \\mid \\hat{\\theta}_{m l e}\\right)+2 \\text{pAIC}  \\tag{式 5.5} \n",
    "$$\n",
    "\n",
    "其中，$pAIC$ 表示参数的个数， $\\hat{\\theta}_{m l e}$ 为 $\\theta$ 的最大似然估计。最大似然估计在非贝叶斯方法中经常用到，等价于贝叶斯方法中基于均匀先验的最大后验估计。注意这里 $\\hat{\\theta}_{mle}$ 是点估计而不是分布。\n",
    "\n",
    "同样，此处 −2 也是出于历史原因。从实用角度来看，上式中的第 1 项考虑的是模型对数据的拟合效果，第 2 项衡量的是模型复杂度。因此，如果两个模型对数据的解释能力相同，但是其中一个比另一个的参数更多的话，`AIC` 会告诉我们应该选择参数更少的那个。\n",
    "\n",
    "`AIC` 对非贝叶斯方法来说很有用，但对贝叶斯方法可能会有些问题。原因是 `AIC` 没有使用后验，因而将估计中的不确定信息丢失了，此外将均匀分布作为先验，对使用非均匀先验的模型来说不太合适。因为在使用非均匀先验时，不能简单地计算模型中参数的个数，合理使用非均匀先验实际上相当于对模型已经使用了正则，并且会降低过拟合的可能，也就是说带正则模型的有效参数个数可能比真实参数个数要少。类似情况在多层模型中也会出现，毕竟多层模型可视为从数据中学习先验的有效方式。\n",
    "\n",
    "#### **（ 3 ） `WAIC` 通用信息准则**\n",
    "\n",
    "`通用信息准则（Widely Available Information Criterion， WAIC）` 是 `AIC` 的完全贝叶斯版本。与 `AIC` 一样， `WAIC` 有两个项：一项衡量模型对数据的拟合效果；另外一项衡量模型的复杂程度。\n",
    "\n",
    "$$\n",
    "\\text{ WAIC }=-2 \\times lppd + 2 \\times \\text{pWAIC} \\tag{5.5}\n",
    "$$\n",
    "\n",
    "如果您想更好地理解这两个术语是什么，请阅读后面的 `深入 WAIC` 部分。从应用角度看，只需要知道我们更喜欢较低的值。\n",
    "\n",
    "#### **（ 4 ）`Pareto` 平滑重要性采样留一交叉验证**\n",
    "\n",
    "`Pareto 平滑重要性采样留一交叉验证` 是一种用于近似留一法交叉验证结果但又不实际执行 `K` 次迭代的方法。该方法不是一个信息准则，但提供的结果与 `WAIC` 非常相似，并且在某些条件下， `WAIC` 和 `LOO` 都是渐近收敛的。该方法主要思想是通过对似然适当重新加权来近似留一法交叉验证，在统计学中可以通过重要性采样来实现。但普通的重要性采样结果不稳定，为引入了称为 `Pareto 平滑重要性采样 (PSIS)` 的新方法，用来计算更可靠的留一法估计值。\n",
    "\n",
    "该方法与 `AIC` 和 `WAIC` 类似，其结果数值越低，模型估计预测的精度就越高。因此，通常更倾向于选择数值较低的模型。\n",
    "\n",
    "#### **（ 5 ）`DIC` 与 `BIC` 准则**\n",
    "\n",
    "另一种常见的信息准则是 `差分信息准则（DIC）` 。但无论在理论上还是在实践上， `WAIC` 都被证明比 `DIC` 更有效，因此推荐使用 `WAIC` 而不是 `DIC`。\n",
    "\n",
    "另一个信息准则是 `贝叶斯信息准则（BIC）`，它类似于 Logistic 回归。 `BIC` 的提出是为了纠正 `AIC` 的一些问题，作者建议采用贝叶斯纠正。但 `BIC` 并不是真正的贝叶斯，实际上它与 `AIC` 非常相似。它假设平坦的先验，并使用最大似然估计。更重要的是， `BIC` 不同于 `AIC` 和 `WAIC` ，而更多涉及 `贝叶斯因子（Bayesian Factor）` 的概念，这点将在本章后面讨论。\n",
    "\n",
    "### 5.3.3 使用 `PyMC3` 做模型比较\n",
    "\n",
    "采用 `ArviZ` 进行模型比较比想像中容易得多！"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3713fc16",
   "metadata": {},
   "outputs": [],
   "source": [
    "waic_l = az.waic(trace_l)\n",
    "waic_l"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3b07a93",
   "metadata": {},
   "source": [
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210525103006_9f.webp)\n",
    "\n",
    "如果你想计算 `LOO` 而不是 `WAIC` ，需要使用 `az.loo` 。对于 `WAIC` 和 `LOO` ，`PyMC3` 报告了四个值（见上表）：\n",
    "- 一个点估计值\n",
    "\n",
    "- 点估计的标准差（假设正态分布，样本量较少时不太可靠）\n",
    "\n",
    "- 有效参数的数量\n",
    "\n",
    "- 警告数量\n",
    "\n",
    "> **注意：**\n",
    ">\n",
    ">在计算 `WAIC` 或 `LOO` 时，可能会收到一些警告消息，指出计算的结果可能不可靠。此警告是根据经验确定的阈值提出的（请参阅相关文献资料）。虽然这不一定是错误，但可能表明这些度量计算存在问题。 `WAIC` 和 `LOO` 相对较新，或许需要开发更好的方法来获得其可靠性。\n",
    "> \n",
    "> 无论如何，如果出现警告的情况，首先应当确保有足够样本，并且是一个混合良好、可靠的样本（ 参见第 8 章 ）。如果仍然接收到警告， `LOO` 方法的提出者建议使用更健壮的模型，如使用学生 $t$ 分布而不是高斯分布。如果上述建议都不起作用，可能需要考虑换一种验证方法，例如直接执行 `K-折交叉验证`。 \n",
    "\n",
    "*`WAIC` 和 `LOO` 只能帮助你在一组给定的模型中进行选择，但不能帮助你决定模型是否真的是解决特定问题的好方法*。因此， **`WAIC` 和 `LOO` 应该得到后验预测检查以及任何其他信息和测试的补充**，这些信息和测试可以帮助我们根据待解决的特定问题和领域知识来设置模型和数据。\n",
    "\n",
    "由于 `WAIC` 和 `LOO` 总是以多个模型的相对方式进行解释，`ArviZ` 提供了两个辅助函数来简化比较。第一个是 `az.compare` ："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64f0c6eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmp_df = az.compare({'model_l':trace_l, 'model_p':trace_p}, method='BB-pseudo-BMA')\n",
    "cmp_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd41999f",
   "metadata": {},
   "source": [
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210525103242_82.webp)\n",
    "\n",
    "这里有很多列：\n",
    "\n",
    "- 第 1 列为 `WAIC` 的值。表格默认以该列升序排列，索引列则反映了该排序。\n",
    "- 第 2 列是估计的有效参数个数。一般来说，参数越多的模型数据拟合越灵活，但也更可能导致过拟合。因此，可以将 `pwaic` 解释为惩罚性术语，也可将其解释为度量每个模型在拟合数据方面的灵活性。\n",
    "- 第 3 列是 `WAIC` 的相对值，以排名最高的模型 `WAIC` 值为基准，列出各模型 `WAIC` 值与基准值间的相对差，第一个模型的值始终为 0 。\n",
    "- 第 4 列为权重。在比较模型时，有时并不想选择量值指示最好的模型，而是希望通过平均若干模型来进行预测，并且通过加权平均，赋予不同模型适当的权重（见 5.5 节）。比较常用的方法是基于 `WAIC` 值给每个模型赋予一个 `Akaike 权重` 。在给定数据时，这些权重解释为每个模型的概率。此方法存在的一个问题是：由于该权重的计算基于 `WAIC` 的点估计，因此不确定性被忽略了。\n",
    "- 第 5 列记录了 `WAIC` 值的标准差。标准差可用于评估 `WAIC` 估计的不确定度。\n",
    "- 第 6 列记录了第 2 列相对值的标准差。由于 `WAIC` 的不确定性在不同模型之间相关程度不同，不同模型应当拥有不同的值。\n",
    "- 第 7 列名为 `WARNING` 。值 1 表示 `WAIC` 的计算可能不可靠。\n",
    "\n",
    "我们还可以通过使用 `az.plot_compare` 函数可视化上述信息。该函数接受 `az.compare` 的输出，并以 `Richard McElreath` 的《统计反思》一书中使用的样式生成汇总图："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "032e34f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "az.plot_compare(cmp_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2b4039b",
   "metadata": {},
   "source": [
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210525103532_b8.webp)\n",
    "\n",
    "> **图 5.8 用于比较模型的 WAIC 信息准则图示**。 空圆圈表示 `WAIC` 值，水平直线表示其标准差，黑色圆圈表示样本内离差。 `WAIC` 值最小（最左侧）的模型通常被视为相对最优模型，除最优模型外，其他模型均绘制了一个以三角形为中心的灰色直线，代表该模型与最优模型之间的相对值。\n",
    "\n",
    "在此详细描述一下图 5.8：\n",
    "\n",
    "- 空圆圈代表 `WAIC` 的值，与之相关的黑色误差条是 `WAIC` 的标准差。\n",
    "\n",
    "- 最低的 `WAIC` 值用一条垂直的灰色虚线表示，以便于与其他 `WAIC` 值进行比较。\n",
    "\n",
    "- 实心黑圆圈是每个模型的`样本内离差`，对于 `WAIC` 来说，它与相应的 `WAIC` 值相差 $2 \\times \\text{pWAIC}$ 。\n",
    "\n",
    "- 除最佳模型外，其他模型会有一个三角形，表示该模型和最佳模型之间的 `WAIC` 相对值，配套的有一个灰色误差条，表示 `WAIC` 相对值的标准差。\n",
    "\n",
    "依据信息准则来选择模型，可能是一种最简单的决策方式。此时只需要选择信息准则值较低的模型，而忽略其他模型即可。如果遵循这种方式，前例中的二次多项式模型可能是最佳选择。请注意，标准差不重叠给做出此选择提供了信心。相反，如果标准差是重叠的，则应该提供一个更微妙的答案。\n",
    "\n",
    "## 5.4 贝叶斯框架下的模型比较\n",
    "\n",
    "> **背景：**\n",
    ">\n",
    ">零假设显著性检验（ `NHST` ）是频率主义进行数据分析的主要工具。但在统计学领域，`NHST` 受到了广泛批评。越来越多的统计学者提倡使用贝叶斯方法检验研究假设，在实证研究中也有越来越多的学者使用贝叶斯因子进行数据分析。\n",
    " \n",
    "### 5.4.1 基本概念\n",
    "\n",
    "在贝叶斯世界中，评估和比较模型的一种常见方式是 `贝叶斯因子（Bayes factor, BF）` 。 为理解什么是贝叶斯因子，让我们重温一遍贝叶斯定理：\n",
    "\n",
    "$$\n",
    "p(\\theta \\mid \\mathcal{D})=\\frac{p(\\mathcal{D} \\mid \\theta) p(\\theta)}{p(\\mathcal{D})} \\tag{式 5.6} \n",
    "$$\n",
    "\n",
    "这里，$\\mathcal{D}$ 表示数据，由于通常模型是默认的，所以被隐藏了。这里为了区分不同的模型，可以显式地将模型 $M_k$ 加入到条件概率公式中：\n",
    "\n",
    "$$\n",
    "p\\left(\\theta \\mid \\mathcal{D}, M_{k}\\right)=\\frac{p\\left(\\mathcal{D} \\mid \\theta, M_{k}\\right) p\\left(\\theta \\mid M_{k}\\right)}{p\\left(\\mathcal{D} \\mid M_{k}\\right)}\\tag{式 5.7}  \n",
    "$$\n",
    "\n",
    "第一章中曾经介绍过，分母中的术语为边缘似然（或证据），可视为一个归一化常数。在进行单模型推断时，通常不需要真实计算它，而是基于一个常数因子来计算后验（如： `MCMC` 和 `变分推断` 方法都巧妙地规避了边缘似然的计算）。但对于模型比较和模型平均来说，边缘似然是不得不计算的重要量。\n",
    "\n",
    "前面曾经提到过，**模型选择**的主要目标是从 $k$ 个模型中选择一个最好的模型，理论上我们应当选择边缘似然 $p(\\mathcal{D}|M_k)$ 最大的那个（`边缘似然`之所以又被称为`证据`，就是因为它可以证明哪个模型更能解释数据集）。一般来说， $p(\\mathcal{D}|M_k)$ 的绝对数值本身并不能告诉我们太多信息，重要的是在不同模型之间的相对值。因此，实践中经常计算两个边缘似然的比率，这个比率就是著名的**贝叶斯因子**：\n",
    "\n",
    "$$\n",
    "B F=\\frac{p\\left(\\mathcal{D} \\mid M_{0}\\right)}{p\\left(\\mathcal{D} \\mid M_{1}\\right)} \\tag{式 5.8} \n",
    "$$\n",
    "\n",
    "当 $BF（M_0,M_1） > 1$ 时，表明模型 $0$ 比模型 $1$ 更好地解释了数据。\n",
    "\n",
    "一些研究者为了简化对 $BF$ 的解释，设计了带有范围的表格。例如，下面的列表显示了 “支持模式 $0$ 而不支持模式 $1$ ” 的证据强度：\n",
    "\n",
    "|  BF 值   | 证据强度 |\n",
    "| :-----: | :------: |\n",
    "|   1-3   |   初级   |\n",
    "|  3-10   |   中等   |\n",
    "|  10-30  |    强    |\n",
    "| 30-100  |  非常强  |\n",
    "| 大于 100 |   极端   |\n",
    "\n",
    "不过需要注意的是：这些规则只是一些约定，最终结果始终应放在任务的上下文中并伴随足够细节，以便其他人可以检查是否同意我们的结论。\n",
    "\n",
    "现在对问题做更深入的探讨，如果假设所有模型具有相同的先验概率，则使用 $p(\\mathcal{D}|M_k)$ 来比较模型完全没有问题。但如果先验不同的时候怎么办呢？也就是说，我们**根据专业知识对不同的模型赋予有差异的信念时怎么办？**\n",
    "\n",
    "此时就要对贝叶斯定理和贝叶斯因子做进一步的推广，必须计算模型的后验赔率：\n",
    "\n",
    "$$\n",
    "\\underbrace{\\frac{p\\left(M_{0} \\mid \\mathcal{D}\\right)}{p\\left(M_{1} \\mid \\mathcal{D}\\right)}}_{\\text {posterior odds }}=\\underbrace{\\frac{p\\left(\\mathcal{D}\\mid M_{0}\\right)}{p\\left(\\mathcal{D} \\mid M_{1}\\right)}}_{\\text {Bayes factors}} \\underbrace{\\frac{p\\left(M_{0}\\right)}{p\\left(M_{1}\\right)}}_{\\text{prior odds} } \\tag{式 5.9}  \n",
    "$$\n",
    "\n",
    "### 5.4.2 一些讨论\n",
    "\n",
    "现在简要讨论有关边缘似然的一些关键事实。通过仔细检查定义，可以理解边缘似然的性质和应用效果：\n",
    "\n",
    "$$\n",
    "p\\left(\\mathcal{D} \\mid M_{k}\\right)=\\int_{\\theta_{k}} p\\left(\\mathcal{D} \\mid \\theta_{k}, M_{k}\\right) p\\left(\\theta_{k}, M_{k}\\right) d \\theta_{k} \\tag{式 5.10}  \n",
    "$$\n",
    "\n",
    "- **好处**：参数多的模型比参数少的模型具有更大惩罚。贝叶斯因子内置奥卡姆剃刀，因为参数数量越多，先验分布相对于似然就越宽。结合贝叶斯因子公式，越宽的先验积分（质量）越大，而越聚集的先验积分（质量）越小，从而间接实现了对参数数量的惩罚。\n",
    "- **缺点**：计算边缘似然是艰巨的任务，因为要计算高维参数空间上的多变量函数积分，需要使用复杂方法进行数值求解。\n",
    "- **尴尬之处**：边缘似然对先验取值的依赖过于敏感。\n",
    "\n",
    "使用边缘似然来比较模型是一个好主意，因为复杂模型的惩罚已经包括在内。但同时先验信息的变化会影响边缘似然的计算，其中关键词是 ”敏感“ 。也就是说先验的微小变化，可能会对边缘似然的值产生巨大影响。在前例中，标准差为 $100$ 的正态先验与标准差为 $1000$ 的正态先验变化很小，但相应的贝叶斯因子受这些变化却产生了较大变化。\n",
    "\n",
    "另一个相关的批评是，贝叶斯因子可以被用作进行假设检验的贝叶斯统计方法。这种批评本身没错，但许多文章指出，推断方法比假设检验方法（无论是否为贝叶斯方法）更适合于大多数问题。\n",
    "\n",
    "### 5.4.3 贝叶斯因子的计算\n",
    "\n",
    "贝叶斯因子的计算可以被视为一个分层模型，其中超参数是分配给每个模型并从`类别分布`中采样的索引值（ `index` ）。换句话说，我们同时对多个相互竞争的模型进行推断，并使用在模型间跳跃的离散变量。我们花在每个模型上的采样时间与 $p(M_k|\\mathcal{D})$ 成正比。然后，应用公式 5.9 求出贝叶斯因子。\n",
    "\n",
    "为举例说明贝叶斯因子的计算，我们再来一次抛硬币：\n",
    "\n",
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210525130312_77.webp)\n",
    "\n",
    "图 5.10\n",
    "\n",
    "</center>\n",
    "\n",
    "让我们创建一些数据，以便在示例中使用："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bf4435c",
   "metadata": {},
   "outputs": [],
   "source": [
    "coins = 30\n",
    "heads = 9\n",
    "y_d = np.repeat([0, 1], [coins-heads, heads])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e0fe8c5",
   "metadata": {},
   "source": [
    "现在，来看一下 `PyMC3` 模型。为在之前的代码之间切换，我们使用了 `pm.math.switch` 函数。如果此函数的第一个参数的计算结果为 `true`，则返回第二个参数，否则返回第三个参数。请注意，还使用 `pm.math.eq` 函数来检查 `model_index` 变量是否等于 0 ："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b50be23",
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model_BF:\n",
    "    p = np.array([0.5, 0.5])\n",
    "    model_index = pm.Categorical('model_index', p=p)\n",
    "    m_0 = (4, 8)\n",
    "    m_1 = (8, 4)\n",
    "    m = pm.math.switch(pm.math.eq(model_index, 0), m_0, m_1)\n",
    "    # a priori\n",
    "    θ = pm.Beta('θ', m[0], m[1])\n",
    "    # likelihood\n",
    "    y = pm.Bernoulli('y', θ, observed=y_d)\n",
    "\n",
    "    trace_BF = pm.sample(5000)\n",
    "\n",
    "az.plot_trace(trace_BF)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c74fcf6",
   "metadata": {},
   "source": [
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210525130426_94.webp)\n",
    "\n",
    "图 5.11\n",
    "\n",
    "</center>\n",
    "\n",
    "现在，需要通过计算 `model_index` 变量来计算贝叶斯因子。请注意，我们已经包括了每个模型的先验值："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91708aa3",
   "metadata": {},
   "outputs": [],
   "source": [
    "pM1 = trace_BF['model_index'].mean()\n",
    "pM0 = 1 - pM1\n",
    "BF = (pM0 / pM1) * (p[1] / p[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bab4b82",
   "metadata": {},
   "source": [
    "结果，我们得到的值为$≈11$，这意味着模型 $0$ 比模型 $1$ 高出一个数量级。这非常有意义，因为数据的正面值比预期的 $\\theta=0.5$ 要少，两个模型之间的唯一区别是模型 $0$ 的先验更兼容 $\\theta<0.5$ （背面比正面多），模型 $1$ 更兼容 $\\theta>0.5$ （正面比背面多）。\n",
    "\n",
    "### 5.4.4 计算贝叶斯因子时的常见问题\n",
    "\n",
    "用我们定义的方式计算贝叶斯因子会有一些问题，比如当其中一个模型比另一个模型更好时，根据定义，我们会对更好的这个模型采样次数更多，这可能会导致我们对另外一个模型欠采样。另外，第一个问题是：即使某些参数没有用于拟合数据，也会更新。也就是说，当模型 $0$ 被选择时，模型 $1$ 中的参数也会更新，不过由于这部分参数并没有用于解释数据，值受限于先验。如果先验太模糊，有可能当我们选到模型 $1$ 时，参数值距离上一次被接受的值太远了，因而该步被拒绝，从而导致采样会出现问题。\n",
    "\n",
    "> 此处关于模型平均时的参数更新问题，可以参考 `附录 G` 。\n",
    "\n",
    "如果遇到此类问题，可以对模型进行两个调整以改进采样：\n",
    "\n",
    "- 理想情况下，如果两个模型都访问相同次数，我们会得到一个更好的采样，因此我们对模型的先验做出调整（前一个模型中的 $p$ 值），从而向原来访问频次较低的模型倾斜。这个过程对贝叶斯因子的计算不会有多大影响，因为我们在计算过程中包含了先验。\n",
    "\n",
    "- 根据 Kruschke 以及其他人的建议，可以使用伪先验，其思想很简单：当没被选择的模型的参数出现自由漂移时，可以尝试手动限制它们，不过是在这个模型没被使用的时候。您可以找到一个在 `Kruschke` 书中使用的模型示例，我将该模型移植到了 [`PyMC3`](https:/​/​github.​com/​aloctavodia/​Doing_​bayesian_​data_​analysis) 中。\n",
    "\n",
    "### 5.4.5 用序贯蒙特卡罗方法计算贝叶斯因子\n",
    "\n",
    "另一种计算贝叶斯因子的方法是使用 `序贯蒙特卡罗 (SMC) 采样方法`。我们将在 `第 8 章-推断引擎` 中学习此方法的详细信息。现在只需要知道这个采样器计算的边缘似然估计是一个副产品，可以直接使用它来计算贝叶斯因子。要在 `PyMC3` 中使用 `SMC`，需将 `pm.SMC()` 传递给 `sample` 的 `step` 参数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba487bff",
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model_BF_0:\n",
    "    θ = pm.Beta('θ', 4, 8)\n",
    "    y = pm.Bernoulli('y', θ, observed=y_d)\n",
    "    trace_BF_0 = pm.sample(2500, step=pm.SMC())\n",
    "\n",
    "with pm.Model() as model_BF_1:\n",
    "    θ = pm.Beta('θ', 8, 4)\n",
    "    y = pm.Bernoulli('y', θ, observed=y_d)\n",
    "    trace_BF_1 = pm.sample(2500, step=pm.SMC())\n",
    "    \n",
    "model_BF_0.marginal_likelihood / model_BF_1.marginal_likelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3036f77c",
   "metadata": {},
   "source": [
    "根据 `SMC 方法`，贝叶斯因子也在 11 左右，如果你想用 `PyMC3` 计算贝叶斯因子，我强烈推荐使用 `SMC 方法`。本书中提出的另一种方法在计算上更加繁琐，需要更多手动调整，主要是因为模型间的跳跃需要用户通过反复试验进行更多调整。从这点上来说，`SMC` 是一种自动化程度更高的方法。\n",
    "\n",
    "### 5.4.6 贝叶斯因子与信息准则\n",
    "\n",
    "如果对贝叶斯因子求对数，可以将两个边缘似然的比值转换成求差，这样比较边缘似然就与比较信息准则类似了。但衡量模型的数据拟合程度项以及惩罚项去哪儿了呢？前者包含在了似然部分，而后者包含在对先验取平均的部分。参数越多，先验空间相比似然就越大，因而平均之后似然就会较低，而且参数越多，先验就会越分散，因而在计算证据的时候惩罚越大。这也是为什么人们说贝叶斯理论会很自然地惩罚更复杂的模型，或者称贝叶斯理论自带奥卡姆剃刀。\n",
    "\n",
    "此前说过，贝叶斯因子对先验过于敏感。这在执行推断时会导致本来不相关的差异，在计算贝叶斯因子时被证明为非常重要。现在我们来看一个示例，它将有助于阐明贝叶斯因子在做什么，信息准则在做什么，以及它们如何在相似的情况下专注于两个不同的方面。回到抛硬币示例的数据定义，现在设置 $300$ 枚硬币和 $90$ 个正面；这与以前的比例相同，但数据多了 $10$ 倍。然后，分别运行每个模型："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d2890fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "traces = []\n",
    "waics = []\n",
    "for coins, heads in [(30, 9), (300, 90)]:\n",
    "    y_d = np.repeat([0, 1], [coins-heads, heads])\n",
    "    for priors in [(4, 8), (8, 4)]:\n",
    "        with pm.Model() as model:\n",
    "            θ = pm.Beta('θ', *priors)\n",
    "            y = pm.Bernoulli('y', θ, observed=y_d)\n",
    "            trace = pm.sample(2000)\n",
    "            traces.append(trace)\n",
    "            waics.append(az.waic(trace))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "795e102a",
   "metadata": {},
   "source": [
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210525130903_77.webp)\n",
    "\n",
    "图 5.12\n",
    "\n",
    "</center>\n",
    "\n",
    "通过增加数据，我们几乎完全克服了先验，现在两个模型都做出了类似的预测。用 $30$ 枚硬币和 $9$ 个正面作为数据，可以看到的 $BF \\approx 11$ ，如果用 $300$ 个硬币和 $90$ 个正面的数据重复计算，我们会看到 $BF \\approx 25$ 。贝叶斯因子表明模型 $0$ 比模型 $1$ 更受青睐。当增加数据时，模型之间的决定变得更加清晰。这完全有道理，因为现在我们更确定模型 $1$ 有一个与数据不一致的先验。\n",
    "\n",
    "还要注意，随着数据量增加，两个模型的 $\\theta$ 值趋于一致；实际上，两个模型的值都大约是 $0.3$。因此，如果决定用 $\\theta$ 来预测新的结果，将与计算 $\\theta$ 的分布的模型几乎没有什么不同。\n",
    "\n",
    "现在，比较一下 `WAIC` 告诉我们的内容（参见图 5.13）。模型 $0$ 的 `WAIC` 是 $368.4$ ，模型 $1$ 的是 $368.6$ ，直觉上差别不大。比实际差异更重要的是，如果重新计算数据的信息准则，也就是 $30$ 枚硬币和 $9$ 个正面，你会得到模型 $0$ 的 $38.1$ 和模型 $1$ 的 $39.4$ 。也就是说，在增加数据时，相对差异变得越小，$\\theta$ 的估计值越相近，与信息准则估计出的预测准确度的值就越相似。如果你用 `LOO` 代替 `WAIC` ，会发现本质上是一样的："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad529573",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, sharey=True)\n",
    "labels = model_names\n",
    "indices = [0, 0, 1, 1]\n",
    "for i, (ind, d) in enumerate(zip(indices, waics)):\n",
    "    mean = d.waic\n",
    "    ax[ind].errorbar(mean, -i, xerr=d.waic_se, fmt='o')\n",
    "    ax[ind].text(mean, -i+0.2, labels[i], ha='center')\n",
    "ax[0].set_xlim(30, 50)\n",
    "ax[1].set_xlim(330, 400)\n",
    "plt.ylim([-i-0.5, 0.5])\n",
    "plt.yticks([])\n",
    "plt.subplots_adjust(wspace=0.05)\n",
    "fig.text(0.5, 0, 'Deviance', ha='center', fontsize=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ba11323",
   "metadata": {},
   "source": [
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210525131100_4d.webp)\n",
    "\n",
    "图 5.13\n",
    "\n",
    "</center>\n",
    "\n",
    "**贝叶斯因子关注的是哪个模型更好，而 `WAIC` （和 `LOO` ) 关注的是哪个模型能给出更好的预测。** 如果检查公式 5.5 和 5.11，你就会看到这些不同。 `WAIC` 和其他信息准则一样，以这样或那样的方式使用对数似然，先验并不直接作为计算的一部分。先验只间接参与，辅助我们估计。取而代之的是，贝叶斯因子直接使用先验，因为我们需要对先验值的整个范围内的似然进行平均。\n",
    "\n",
    "\n",
    "## 5.5 模型平均\n",
    "\n",
    "> 注：关于本节更详细的信息，参见 `附录 G` 。\n",
    "\n",
    "通过交叉验证、信息准则、贝叶斯因子等方法对众多模型进行比较，然后做出模型选择的策略，实践起来比较简单。但大家心里应当清楚，**当我们在选择模型的时候，也正在抛弃有关模型的不确定性信息** 。这种情况有些类似于计算了完整的后验分布，但最后只保留了后验均值。这会造成我们对模型过于自信。\n",
    "\n",
    "一种变通的方案是在执行模型选择时，报告和讨论不同模型的信息准则值、标准差等统计量以及后验预测检查情况。将所有这些数字和检查放在问题上下文中很重要，只有这样，相关人士才能更好地感受到模型可能存在的局限性和缺点。在学术界中，这种作法比较常见，在很多论文、演示文稿的讨论部分中，能够找到大量的相关案例。\n",
    "\n",
    "除了上述对各模型均做出报告和讨论的方法外，还有一种做法是充分利用模型比较中的出现的不确定性，执行模型平均。其原理可以理解为：假设有 $S$ 个模型（或者 $S$ 种估计方法），每个模型得出了一个预测值，可记为 $\\tilde\\mu_i$， 其中 $i \\in \\{1...S\\}$ ，那么模型平均方法得到的最终估计值为：\n",
    "\n",
    "$$\n",
    "\\tilde \\mu = \\tilde W_1\\tilde\\mu_1 + \\tilde W_2\\tilde\\mu_2 + ... + \\tilde W_s\\tilde\\mu_S\n",
    "$$\n",
    "\n",
    "\n",
    "从公式中不难看出，模型平均中比较重要的是权重，该权重是数据依赖的。模型平均是处理模型不确定性和提高预测效果的方法，其核心研究内容是权重选择以及不确定性的衡量。\n",
    "\n",
    "> **注解：**\n",
    "> \n",
    "> - 之前讨论的重点是参数的不确定性问题，此处为模型的不确定性。\n",
    "> - 能够解释数据集的模型很多，**模型比较**试图从其中选择一个最优的，而**模型平均**则认为所有模型的加权平均可能是最好的估计。\n",
    "> - 从广义上理解，模型平均是集成方法的一种。\n",
    "> - 有必要将**模型平均**与**模型组合**区分开，模型组合通常对预测变量的不同空间采用不同的模型，常见方法是**决策树**。 在模型平均中整个数据集由单一的模型生成，而在模型组合中，数据集里的不同数据点可能由不同模型生成。\n",
    "\n",
    "下面介绍其中几种比较常用的模型平均方法：\n",
    "\n",
    " ### 5.5.1 基于信息准则值计算权重\n",
    "\n",
    "本方法使用每个模型的加权平均值来生成 `元模型（meta-model）` 和 `元预测（meta-predictions）` 。基于某些信息准则值（如 `WAIC`）计算不同模型的权重，公式如下：\n",
    "\n",
    "$$\n",
    "w_{i}=\\frac{e^{\\frac{1}{2} d E_{i}}}{\\sum_{j=1}^{M} e^{-\\frac{1}{2} d E_{j}}} \\tag{式 5.11}  \n",
    "$$\n",
    "\n",
    "这里 $dE_i$ 是第 $i$ 个模型相对于最佳模型（`WAIC`值最小的模型）的 `WAIC` 相对差值。除 `WAIC` 外，此处也可以使用其他信息准则值，如 `AIC` 或 `LOO` 等。此公式是根据 `WAIC` 值计算各模型相对概率的启发式方法。分母为归一化因子，第 4 章中有过类似的表达式。\n",
    "\n",
    "采用式 5.11 的权重对模型进行加权平均称为 `伪贝叶斯模型平均` 。真正的贝叶斯模型平均应当使用边缘似然，而非 `WAIC` 或 `LOO` 。不过尽管边缘似然在理论上很有吸引力，但在模型比较和模型平均中使用很少。大多还是选择 `WAIC` 或 `LOO` 而非边缘似然。\n",
    "\n",
    ">**注解：** 边缘似然是模型解释数据的能力体现，理论上边缘似然越大赋予的权重也应当越大。\n",
    "\n",
    "使用 `PyMC3` ，可以通过将 `method=‘pseudo-BMA’` （伪贝叶斯模型平均）参数传递给 `az.compare` 函数来计算式 5.7 的权重。其问题是未考虑计算 $E_i$ 时的不确定性。通过高斯近似可以计算每一个 $E_i$ 的标准差。这也是函数 `az.waic`、`az.loo` 和 `az.compare` 在传递 `method=‘pseudo-BMA’` 参数时返回的误差值。此外，还可以使用 `贝叶斯自举（Bayesian bootstrapping）法` 来估计不确定性。这是一种比高斯近似更可靠的方法。通过将 `method=‘BB-pseudo-BMA’` 传递给 `az.compare` 函数即可实现。\n",
    "\n",
    "### 5.5.2 基于后验预测分布计算权重\n",
    "\n",
    "另一种计算平均模型权重的方法被称为 `预测分布堆叠（stacking of predictive distributions）` 。这在 `PyMC3` 中通过将 `method=‘stacking’` 传递给 `az.compare` 实现。其基本思想是通过最小化元模型和真实生成模型之间的差异，将多个模型组合到一个元模型中。当使用对数打分规则时，这等价于：\n",
    "\n",
    "$$\n",
    "\\max _{n} \\frac{1}{n} \\sum_{i=1}^{n} \\log \\sum_{k=1}^{K} w_{k} p\\left(y_{i} \\mid y_{-i}, M_{k}\\right) \\tag{式 5.12}  \n",
    "$$\n",
    "\n",
    "这里，$n$ 是数据点的数量，$k$ 是模型的数量。为了强制实施方案，我们将 $w$ 约束为 $w_k \\geq 0$ 并且 $\\sum w_k =1$。量 $p(y_i|y_{-i},M_k)$ 是模型 $M_k$ 的留一预测分布。根据留一法，计算需要拟合每个模型 $n$ 次，每次遗留一个数据点。幸运的是，`PyMC3` 可以使用 `WAIC` 或 `LOO` 来近似留一预测分布。\n",
    "\n",
    "### 5.5.3 其他模型平均方法\n",
    "\n",
    "**（1）经典贝叶斯模型平均法**\n",
    "\n",
    "在上述基于预测能力的模型平均方法之外，还有一种最为经典的模型平均方法，该方法显式构建包括所有感兴趣模型作为子模型的元模型。可以构建这样一个模型：我们对每个子模型的参数进行推断，同时计算每个模型的相对概率（有关此方面的示例，请参阅 `贝叶斯因子` 一节 和 `附录 G` 的`第 2.2 节`）。\n",
    "\n",
    "**（2）连续贝叶斯模型平均**\n",
    "\n",
    "除了对离散的模型求平均之外，有时候还可以将模型视为连续的。一个简单的示例就是，假设我们有一个抛硬币问题以及两个不同的模型，其中之一的先验偏向正面朝上，另一个偏向于反面朝上。则我们可以按照离散的模型平均方式，用两个模型去拟合并用 `dIC` 权重求平均；但也可以构建一个分层模型估计先验分布，而此时构建的不再是两个离散的模型了，而是一个连续模型，其中包含两个离散的模型作为特例。\n",
    "\n",
    "哪种方法更好呢？还是要具体问题具体分析，最终使用哪一个取决于实际问题是更适合用离散模型还是连续模型去描述。\n",
    "\n",
    "- 我们是否真的有很好的理由考虑离散模型，或者我们的问题更好地表示为连续模型？\n",
    "- 对于问题来说，挑出一个模型很重要，因为我们是从相互竞争的解释角度思考的，或者平均是更好的想法，因为我们对预测更感兴趣，或者真的可以将流程生成过程视为子流程的平均吗？\n",
    "\n",
    "所有这些问题都不是由统计数据来回答的，而是由领域知识背景下的统计数据来提供信息的。\n",
    "\n",
    "### 5.5.4 模型平均后验预测示例\n",
    "\n",
    "以下只是如何从 `PyMC3` 获得加权后验预测样本的一个模拟示例。在这里，我们使用的是 `pm.sample_posterior_predictive_w` 函数（注意函数名称末尾的 `w` ）。`pm.sample_posterior_predictive` 和`pm.sample_posterior_predictive_w` 之间的区别在于，后者接受多个模型的轨迹和模型，以及相应权重列表（默认值权重相同）。你可以通过 `az.compare` 或其他来源获取这些权重："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0805837d",
   "metadata": {},
   "outputs": [],
   "source": [
    "w = 0.5\n",
    "y_lp = pm.sample_posterior_predictive_w([trace_l, trace_p],\n",
    "                                        samples=1000,\n",
    "                                        models=[model_l, model_p],\n",
    "                                        weights=[w, 1-w])\n",
    "_, ax = plt.subplots(figsize=(10, 6))\n",
    "az.plot_kde(y_l, plot_kwargs={'color': 'C1'},\n",
    "            label='linear model', ax=ax)\n",
    "az.plot_kde(y_p, plot_kwargs={'color': 'C2'},\n",
    "            label='order 2 model', ax=ax)\n",
    "az.plot_kde(y_lp['y_pred'], plot_kwargs={'color': 'C3'},\n",
    "           label='weighted model', ax=ax)\n",
    "plt.plot(y_1s, np.zeros_like(y_1s), '|', label='observed data')\n",
    "plt.yticks([])\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "266c492c",
   "metadata": {},
   "source": [
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210524232628_ac.webp)\n",
    "\n",
    "> **图 5.9 对一阶和二阶多项式模型做模型平均**。 橙色为一阶模型的预测曲线，绿色为二阶模型的预测曲线，中间的粉色为模型平均后的预测曲线。\n",
    "\n",
    "</center>\n",
    "\n",
    "前面提到这是一个虚拟示例，因为与线性模型相比，二次多项式模型的 `WAIC` 值非常低，第一个模型的权重基本上是 1 ，而后者权重基本上是 0 ，为生成图 5.9 ，我假设了这两个模型具有相同的权重 $0.5$ 。\n",
    "\n",
    "## 5.6 其他\n",
    "\n",
    "### 5.6.1 正则先验\n",
    "\n",
    "使用强信息和弱信息先验是在模型中引入偏差的一种方式，如果操作得当，这可能是一个非常好的方法，因为偏差可以防止过拟合，从而有助于模型做出泛化性能更好的预测。在不影响模型建模能力的情况下，添加偏差以减少泛化误差的想法称为正则化。这种正则化通常采用对模型参数数量的较大值实施惩罚的形式。\n",
    "\n",
    "正则先验是一种减少模型所能表示信息的方法，能够降低了型捕获噪声而不是信号的机会。\n",
    "\n",
    "正则化思想如此强大和有用，以至于本书中已经出现了很多次。在非贝叶斯统计中，正则化思想表现为对最小二乘法的两种修正： 岭回归 和套索回归 。从贝叶斯观点来看：\n",
    "\n",
    "- **岭回归**：可解释为『对线性模型的 $\\beta$ 系数采用了标准差趋近于 0 的正态分布作为先验』，使该系数趋向于零。\n",
    "\n",
    "- **套索回归**：可解释为『对线性模型的  $\\beta$ 系数采用了拉普拉斯先验分布』。拉普拉斯分布看起来类似于高斯分布，但其一阶导数在零处没有定义，因为它在零处有一个非常尖锐的峰值（参见图 5.14）。与正态分布相比，拉普拉斯分布使其概率质量更接近于零。使用这种先验的出发点是提供 `正则化` 的同时实现 `变量选择`。其思路是，由于峰值为零，预计先验会导致稀疏性，也就是说，我们创建了一个具有许多参数的模型，先验将自动使大多数参数为零，只保留对模型输出有贡献的相关变量。不幸的是，贝叶斯套索不是这样工作的，基本上是为了有很多参数，拉普拉斯先验迫使非零参数变小。幸运的是，并不是所有东西都丢失了 -- 有一些贝叶斯模型可以用来诱导稀疏性和执行变量选择。\n",
    "\n",
    "值得注意的是，经典版本的岭回归和套索回归对应于点估计，而贝叶斯版本则给出了完整的后验分布结果：\n",
    "\n",
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210524233016_7c.webp)\n",
    "\n",
    "图 5.14\n",
    "</center>\n",
    "\n",
    "### 5.6.2 深入探讨 `WAIC` \n",
    "\n",
    "如果展开公式 5.5，会得到以下结果：\n",
    "\n",
    "$$\n",
    "\\text{WAIC}=-2 \\sum_{i}^{n} \\log \\left(\\frac{1}{S} \\sum_{s=1}^{S} p\\left(y_{i} \\mid \\theta^{s}\\right)\\right)+2 \\sum_{i}^{n}\\left(\\text{V}_{s=1}^{S}\\left(\\log p\\left(y_{i} \\mid \\theta^{s}\\right)\\right)\\right. \\tag{式 5.14}  \n",
    "$$\n",
    "\n",
    "该表达式中的两项看起来非常相似。第一项是式 5.5 中的`对数点预测密度（lppd）`，计算的是后验样本集 $S$ 的平均似然。我们对每个数据点都先求平均似然，然后取对数，最后对所有数据点求和。请将这一项与公式 5.3 和 5.4 进行比较。其实该项就是考虑了后验的样本内离差（deviance）。因此，如果我们认为计算对数似然是衡量模型适合性的好方法，那么在贝叶斯方法中，从后验计算对数似然就顺理成章。观测数据的 lddp 是对未来数据 lppd 的高估（此处意指样本内离差通常小于样本外离差），因此引入第二项来修正这种过高的估计。第二项计算后验样本的对数似然方差，我们对每个数据点执行此方差计算，然后对所有数据点进行汇总。为什么方差会给出惩罚条件？这与贝叶斯因子内置奥卡姆剃须刀的原理相似。有效参数越多，后验分布越大。当向模型添加结构时（如具有信息性/正则化的先验或分层依赖），与非正则化的模型相比，我们约束了后验，进而减少了有效参数的数量。\n",
    "\n",
    "### 5.6.3 熵与最大熵原理\n",
    "\n",
    "从数学上讲，熵定义为：\n",
    "\n",
    "$$\n",
    "H(p)=-\\sum_{i} p_i\\text{log} (p_i) \\tag{式 5.15}  \n",
    "$$\n",
    "\n",
    "直观地说，分布越分散，其熵越大。通过运行以下代码并查看图 5.15，可以看到这一点："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4975bddc",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(912)\n",
    "x = range(0, 10)\n",
    "q = stats.binom(10, 0.75)\n",
    "r = stats.randint(0, 10)\n",
    "true_distribution = [list(q.rvs(200)).count(i) / 200 for i in x]\n",
    "q_pmf = q.pmf(x)\n",
    "r_pmf = r.pmf(x)\n",
    "_, ax = plt.subplots(1, 3, figsize=(12, 4), sharey=True,\n",
    "                     constrained_layout=True)\n",
    "for idx, (dist, label) in enumerate(zip([true_distribution, q_pmf, r_pmf], ['true_distribution', 'q', 'r'])):\n",
    "    ax[idx].vlines(x, 0, dist, label=f'entropy ={stats.entropy(dist):.2f}')\n",
    "    ax[idx].set_title(label)\n",
    "    ax[idx].set_xticks(x)\n",
    "    ax[idx].legend(loc=2, handlelength=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4039357",
   "metadata": {},
   "source": [
    "<center>\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210525131621_8e.webp)\n",
    "\n",
    "图 5.15\n",
    "</center>\n",
    "\n",
    "如图所示，图中的分布 $r$ 是三种分布中较广的一个，也是熵最大的一个。建议使用代码并探索熵是如何变化的（参见练习 10）。\n",
    "\n",
    "#### （1） 熵与方差\n",
    "\n",
    "熵和方差概念上相关，以至于很多人将其声明为度量 ”数据分布的方差“ 的一种特殊形式。不过尽管两个概念相关，但本质上并不相同。在某些情况下，熵增加意味着方差增加，例如高斯分布。然而，也存在方差增加但熵不变的示例。例如混合了两个高斯分布的混合模型，当增加两个高斯分布的众数之间距离的时候，导致大部分点到平均值的距离增加了，即方差增大了。但此时，熵受到的影响极小，因为随着众数之间距离的增加，众数之间的点的概率越来越小，因此它们对总熵的贡献逐步可以忽略不计。从熵的角度来看，如果从两个重叠的高斯开始，将一个相对于另一个移动，则在某个点上，将有两个分离的高斯。\n",
    "\n",
    "#### （2） 最大熵原理\n",
    "\n",
    "熵与信息及其不确定性也有关。事实上，更分散或更平坦的先验是弱信息先验。这不仅直观上是正确的，而且有熵的理论支撑。事实上，在贝叶斯学派中有一个群体，在用熵来证明弱信息先验或正则化先验是合理的。这就是`最大熵原理`：**我们总是期望在问题定义的约束下找到具有最大可能熵的分布**。这是一个可用数学方法解决的优化问题，但本书不讲解细节。下面仅列出一些常见约束条件下的最大熵分布：\n",
    "\n",
    "- 无约束：均匀分布（连续或离散，取决于变量类型）\n",
    "- 正均值：指数分布\n",
    "- 给定方差：正态分布\n",
    "- 只有 2 个非定序输出和 1 个常数均值：二项分布，如果有罕见事件，则为泊松分布（泊松可视为 $p$ 很小，$n$ 很大时的二项分布）\n",
    "\n",
    "有趣的是，许多传统的广义线性模型（如在第 4 章中看到的模型），都是在给定模型约束下，使用最大熵分布来定义的。\n",
    "\n",
    "### 5.6.4 关于 $KL$ 散度\n",
    "\n",
    "现在简单谈谈 `Kullback-Leibler(KL) 散度`，或简称 `KL 散度`。这是在阅读统计学、机器学习、信息论或统计力学文献时经常遇到的概念。你或许会说，`KL 散度`、`熵`、`边缘似然`等概念反复出现的原因很简单，因为所有这些学科都在讨论同一组问题，只是观点略有不同。`KL 散度` 非常有用，因为**它是衡量两个分布接近程度的一种方法**，其定义如下：\n",
    "\n",
    "$$\n",
    "D_{K L}(p \\| q)=\\sum_{i} p_{i} \\log \\frac{p_{i}}{q_{i}} \\tag{式 5.16} \n",
    "$$\n",
    "\n",
    "上式可读为 $q$ 到 $p$ 的 `Kullback-Leibler 散度`（两者顺序不能相反，因为 `KL 散度` 不符合交换率），其中 $p$ 和 $q$ 是两个概率分布。对于连续变量应该计算积分而非求和，但主要思想相同。\n",
    "\n",
    "可以将 $D_{KL}({p||q})$ 散度解释为 **”通过使用概率分布 $q$ 来近似真实分布 $p$ 而引入的额外熵或不确定性“**。事实上，`KL 散度` 是两个熵之间的差值：\n",
    "\n",
    "$$\n",
    "D_{K L}(p \\| q)=\\underbrace{\\sum_{i} p_{i} \\log p_{i}}_{\\text {entropy of p }}-\\underbrace{\\sum_{i} p_{i} \\log q_{i}}_{\\text {crossentropy of p,q }}=\\sum_{i} p_{i}\\left(\\log p_{i}-\\log q_{i}\\right) \\tag{式 5.17} \n",
    "$$\n",
    "\n",
    "利用对数性质，可以重新排列式 5.17 以恢复式 5.16。从式 5.17 的角度来看，也可以将 $D_{KL}({p||q})$ 理解为 $p$ 相对于 $q$ 的相对熵（这一次顺着念）。\n",
    "\n",
    "作为一个简单示例，我们可以使用 $KL$ 散度来评估哪个分布（ $q$ 或 $r$ ）更接近真实分布。使用 `Scipy`，可以计算 $D_{KL}({真实分布||q})$  和 $D_{KL}({真实分布||r})$ ："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b40219be",
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.entropy(true_distribution, q_pmf), stats.entropy(true_distribution,r_pmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6af58387",
   "metadata": {},
   "source": [
    "如果运行上段代码，您将获得 $\\approx 0.0096,\\approx 0.7394$ 。因此可以判定，$q$ 比 $r$ 更接近于真实分布，因为它引入的额外不确定性更小。我希望您同意我的观点，即这个数值结果与您检查图 5.15 所预期的一致。\n",
    "\n",
    "您可能很想将 $KL$ 散度描述为距离，但它是不对称的，因此不是真实距离。如果运行下面代码，将获得 $\\approx 2.7,\\approx 0.7$ 。由此可见，结果数字是不同的。在此例中，可以看到 $r$ 是 $q$ 的更好近似，但反之可能不成立："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8df5b470",
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.entropy(r_pmf, q_pmf), stats.entropy(q_pmf, r_pmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "998e2d6c",
   "metadata": {},
   "source": [
    " $D_{KL}({p||q})$ 表示 $q$ 有多像 $p$ 。也可从惊喜的角度来考虑，如果在预期 $p$ 的时候突然看到了 $q$ ，我们会有多惊讶。对一个事件的惊讶程度取决于用于判断该事件的信息。我在一个非常干旱的城市长大，每年可能会有一两场真正的暴风雨。然后。我搬到另一省份去上大学，我真的很震惊，至少在雨季，平均每周有一场真正的暴风雨！我的一些同学来自布宜诺斯艾利斯，这是阿根廷最潮湿多雨的省份之一。但对他们来说，降雨频率或多或少是意料之中的。更重要的是，他们可能因为空气不够潮湿，而认为天气可能会下多一点雨。\n",
    "\n",
    "我们也可以使用 $KL$ 散度来比较模型，因为它将给出哪个模型更接近真实分布的后验。但问题是我们并不知道真实分布。因此，KL 散度不能直接适用，但可用它作为论据来修正离差（式 5.3）。如果假设真实分布存在（如下式所示），则其应当独立于任何模型和常数，并以同样方式影响 $KL$ 散度，而与用于近似真实分布的后验分布无关。因此，可以使用离差（依赖于每个模型的部分）来估计我们离真实分布相对有多近，即使我们不知道它。对于公式 5.17 ，通过使用一些代数，可以得到：\n",
    "\n",
    "$$\n",
    "D_{K L}(p \\| q)-D_{K L}(p \\| r) =\\left(\\sum_{i} p_{i} \\log p_{i}-\\sum_{i} p_{i} \\log q_{i}\\right)-\\left(\\sum_{i} p_{i} \\log p_{i}-\\sum_{i} p_{i} \\log r_{i}\\right) \\tag{式 5.18}  \n",
    "$$\n",
    "\n",
    "$$\n",
    "=\\sum_{i} p_{i} \\log q_{i}-\\sum_{i} p_{i} \\log r_{i}\n",
    "$$\n",
    "\n",
    "即使不知道 $p$，我们也可以得出结论，具有更大对数似然（或离差）的分布就是在 $KL$ 散度中更接近真实分布的分布。实践中对数似然（或离差）是从有限样本拟合的模型中获得的。因此，还必须增加一个惩罚项，以纠正对离差的高估，这就引出了 `WAIC` 等信息准则。\n",
    "\n",
    "## 5.7 总结\n",
    "\n",
    "**（1）后验预测检查是一个通用的概念和实践工具。** 它可以帮助我们直观地了解模型捕获数据的能力，以及模型捕获我们感兴趣问题的能力。我们既可以只用一个模型进行后验预测检查，也可以用多个模型进行后验预测检查，因此可以把它作为模型比较的一种最基础方法。后验预测检查大多是通过可视化完成的，但类似 `贝叶斯 p-value` 的数字总结也很有帮助。\n",
    "\n",
    "**（2）好的模型在复杂性和准确性之间有好的平衡。** 我们用多项式回归的经典示例来说明这一特征。我们讨论了两种估计样本外准确度的方法：**交叉验证法和信息准则法**，并集中讨论了后者。从实践角度来看，信息准则是一系列平衡`数据拟合度`和`模型复杂度`的方法。在众多可用信息准则中， `WAIC` 是贝叶斯模型中最有用的。另一个有用的方法是 `PSIS-LOO-CV` （或 `LOO` ），它在实践中提供了与 `WAIC` 非常相似的结果。 `WAIC` 和 `LOO` 既可用于模型选择，也可用于计算模型平均的权重。模型平均不是选择单个最佳模型，而是通过对所有可用模型进行加权平均来组合所有可用模型。\n",
    "\n",
    "**（3）`贝叶斯因子`是两个模型的边缘似然之比，是贝叶斯框架下进行`模型比较`和`模型平均` 的重要测度**。  贝叶斯因子的计算很有挑战性。本章介绍了使用 `PyMC3` 计算贝叶斯因子的两种方法：一种是直接尝试使用离散 `index` 估计每个模型的相对概率的分层模型，另一种是称为 `序贯蒙特卡罗（Sequential Monte Carlo，SMC）采样方法` ，建议使用后者。\n",
    "\n",
    "**（4）贝叶斯因子对先验非常敏感而且有计算上的挑战性，同时也存在健壮性问题**。我们比较了贝叶斯因子和信息准则，通过示例表明，两者解决的是相关但不相同的问题：贝叶斯因子侧重于确定正确的模型，而信息准则侧重于实现`最佳预测`或`更低的泛化损失`。所有这些方法都或多或少存在问题，但 `WAIC` 和 `LOO` 在实践中相对要健壮得多。\n",
    "\n",
    "## 5.8 习题\n",
    "\n",
    "![](https://gitee.com/XiShanSnow/imagebed/raw/master/images/articles/spatialPresent_20210525214039_d1.webp)"
   ]
  }
 ],
 "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,
   44,
   56,
   67,
   78,
   96,
   100,
   119,
   130,
   133,
   137,
   148,
   161,
   178,
   234,
   254,
   423,
   426,
   449,
   452,
   468,
   470,
   570,
   574,
   578,
   593,
   605,
   609,
   629,
   641,
   651,
   663,
   679,
   693,
   779,
   795,
   850,
   865,
   911,
   913,
   919,
   921
  ]
 },
 "nbformat": 4,
 "nbformat_minor": 5
}