Archive/2023/Notebooks/Lecture5.ipynb · hcp4715/PyBayesian

加入 Gitee

与超过 1200万开发者一起发现、参与优秀开源项目，私有仓库也完全免费：）

文件

克隆/下载

Lecture5.ipynb 69.16 KB

{"cells":[{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"83CE419F2111489993E63BA9A17A2710","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":" # Lecture 5: Conjugate Families  \n \n ## Instructor: Dr. Hu Chuan-Peng  \n"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"F9D74954E247457C8BCCD7DB645F4502","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"## 回顾对先验的选择"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"8DC0D7793C8E407CBD7F9032896D6C0F","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"在lec3 与 lec4 中，我们使用了Beta分布来反映我们对参数$\\pi$的先验认识。  \n\n我们知道若先验可以用$Beta(\\alpha,\\beta)$描述，收集到的数据可以用$Bin(n, \\pi)$描述，后验分布就可以用$Beta(\\alpha+y, \\beta+n-y)$描述  \n\n在lec3中，我们介绍过共轭先验的概念(conjugate prior):  \n\n* 如果后验分布与先验分布属于同类，先验分布被称为似然函数的共轭先验  \n\n事实上，Beta-Binomial分布这类先验-数据组合，被称为共轭家族(conjugate family)，利用这些家族来获得后验，既有计算上的便利性，由于我们明确知道后验的分布类型，解释起来也很容易。"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"F207F3443B504C219F4AC7DADB819982","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**非共轭先验会带来什么**  \n\n让我们再回到纣王支持率的例子，考虑一个非共轭先验的状况  \n\n假如此时的先验分布$f(\\pi)$不是$beta$分布，而是以下形式：  \n$$  \n\\begin{equation}  \nf(\\pi)=e-e^\\pi\\; \\text{ for } \\pi \\in [0,1].  \n\\tag{5.1}  \n\\end{equation}  \n$$  \n\n![](https://www.bayesrulesbook.com/bookdown_files/figure-html/non-conjugate-1.png)  \n"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"E3886F1CC1EE47DDBAAD2275CB4BCCDD","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"\n> 在有共轭先验$Beta(45,55)$的情况下，后验分布可以简洁地写为$Beta(75,75)$  \n> 但在非共轭先验的情况下，该后验分布的结果变得很繁琐。   \n> 并且在非共轭先验的情况下，我们很难从这个后验表达式中获得类似的直觉。   \n\n\n假设我们在50人的投票结果中观察到10人投了支持票，那么似然函数仍是一个二项分布，可以写成：  \n$$  \nL(\\pi | y=10) = \\left(\\!\\begin{array}{c} 50 \\\\ 10 \\end{array}\\!\\right) \\pi^{10} (1-\\pi)^{40} \\; \\; \\text{ for } \\pi \\in [0,1]  \n$$  \n\n后验可以写成：  \n$$  \nf(\\pi | y = 10) \\propto f(\\pi) L(\\pi | y = 10) = (e-e^\\pi) \\cdot \\binom{50}{10} \\pi^{10} (1-\\pi)^{40}.  \n$$  \n\n加入归一化常数后：  \n\n$$  \n\\begin{equation}  \nf(\\pi|y=10)= \\frac{(e-e^\\pi)  \\pi^{10} (1-\\pi)^{40}}{\\int_0^1(e-e^\\pi)  \\pi^{10} (1-\\pi)^{40}d\\pi}  \\; \\; \\text{ for } \\pi \\in [0,1].  \n\\tag{5.2}  \n\\end{equation}  \n$$"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"C10A1B47DB19447BB4F1402304B08E62","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"## Gamma-Poisson conjugate family"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"C27F7E9428AA4A7DAE1766A06BD73671","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"### The Poisson data model  \n\n**例子：亚运会中国每天获得金牌的数量**：🤔我们该使用哪种分布来描述下面这个例子？  \n\n假设一位体育爱好者每天都会查询中国在亚运会期间每天获得金牌的数量。  \n\n她觉得平均下来，中国每天会获得12枚金牌，并且这个数字大概在5-20之间波动  \n\n我们假设中国每天会获得金牌的**平均数量**设为$\\lambda$  \n> 请注意一些新的希腊字母：λ = lambda  \n\n<iframe src=\"https://tiyu.baidu.com/major/home/%E6%9D%AD%E5%B7%9E%E4%BA%9A%E8%BF%90%E4%BC%9A/tab/%E5%A5%96%E7%89%8C%E6%A6%9C\"></iframe>"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"47F86AF2B3DB44F7BDA59D64851ADE2C","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**🤔我们该使用哪种分布来描述下面这个例子？**  \n\n需要注意的是，通常我们会：  \n1. 为$\\lambda$选择一个合适的先验  \n2. 接着收集数据，并且选择一个合适的数据模型  \n3. 结合先验与数据，更新我们对$\\lambda$的信念  \n\n在这一节中，为了方便探讨共轭先验，我们先从似然函数出发，然后再选取先验，最后讨论后验和共轭先验的关系。  \n- 我们需要先假定似然函数$L(\\lambda|y)$是已知的  \n- 再去选取一个合适的先验分布$p(\\lambda)$  \n- 最后使得后验分布与先验分布具有相同的数学形式。"},{"cell_type":"markdown","metadata":{"id":"797880D6806B41AF9C9042E2FDAEB8E7","notebookId":"65251b1ab0f4556a3142003a","runtime":{"status":"default","execution_status":null,"is_visible":false},"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"}},"source":"假设$\\lambda$是已知的，接下来我们对中国每天会获得金牌的数量$(Y_1,Y_2,\\ldots,Y_n)$进行估计  \n\n$$  \nf(Y_n|\\lambda)  \n$$  \n\n> * 注意，在这里，我们并不能使用二项模型来描述事件的分布情况  \n> * 在二项模型中，我们需要有事件重复的总次数n，和事件成功发生的概率p  \n>$$ f(y|\\pi) = P(Y=y | \\pi) = \\binom{n}{y} \\pi^y (1-\\pi)^{n-y}  $$  \n> * 在这里，我们已知的是$\\lambda$，即每天会获得金牌的平均数量"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"1BC0FF53004F41569EE381D24EF5BAFD","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**Poisson分布**  \n\n\n我们可以使用 **泊松分布(Poisson distribution)** 来表示在给定$\\lambda$下，中国每天会获得金牌在不同数量下的概率  \n\n**Poisson分布**的概率质量函数(pmf)可以表示为：  \n\n$$  \nf(y) =  \\frac{\\lambda^{y}e^{-\\lambda}}{y!}  \n$$  \n\n* Poisson分布只有一个参数$\\lambda$，表示事情发生平均发生率(event rate)或事件发生的期望次数。  \n* y是数据，指在一定时间间隔中事件发生的次数。  \n* $f(y)$表示在某单位时间内，某事件y发生的平均次数的概率。  \n* 例如，已知中国每天平均会获得12枚金牌($\\lambda=12$), 那么中国队明天获得5枚金牌(y=5)的概率为 $f(y)=0.1$小于明天获得13枚金牌(y=13)的概率为 $f(y)=0.3$。  \n"},{"cell_type":"markdown","metadata":{"id":"C728D4D4AAAE43659E61C9A3FC148FD3","notebookId":"65251b1ab0f4556a3142003a","runtime":{"status":"default","execution_status":null,"is_visible":false},"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"}},"source":" poisson分布对数据做出的假设：  \n> * 事件发生比率($\\lambda$)是常数  \n> * 事件的发生是相互独立的  \n\n在这个例子中，中国每天会获得金牌的平均数量是固定的，并且前一天获得金牌的数量不会影响后一天获得金牌的数量"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"8FD26889C30C4D0DAEAF585B4D6324E4","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**Poisson分布图示**  \n\n下图展示了不同$\\lambda$下，事件发生y次的可能性分布"},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"skip"},"id":"407018403F2C4609BC76F35A85162502","notebookId":"65251b1ab0f4556a3142003a","trusted":true},"source":"# @title setup\nimport ipywidgets\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport scipy.stats as st\nfrom scipy.stats import gamma\nfrom ipywidgets import interact\nfrom matplotlib.lines import Line2D\nimport seaborn as sns\n","outputs":[],"execution_count":1},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"fragment"},"trusted":true,"id":"5C0B2B6C5EDC4ACEA6B1B000E19D8941","notebookId":"65251b1ab0f4556a3142003a"},"source":"y = np.linspace(0, 25, 26)  # 假设明天获得金牌的平均数量 0到25次\nλ = 12                      # 假设每天平均获得12枚金牌\nst.poisson.pmf(y, λ)","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 640x480 with 1 Axes>","text/html":"<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAj0AAAHICAYAAAClJls2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8PklEQVR4nO3de3zO9eP/8ee12WZsJuzgTEgyhjmNQqWNfBwqH1I+RJRCtE5UTqkWbqQipz70SU1HUZTMnMJEshzKaU4LY2ZtNppt1/X7w9f7Z40Y297Xrvfjfrtdt891va/34bmWPD+v9+v9ftscDodDAAAALs7N7AAAAADFgdIDAAAsgdIDAAAsgdIDAAAsgdIDAAAsgdIDAAAsgdIDAAAsgdIDAAAsgdIDAAAsgdIDAP/HbrcrODhYb7zxhtlRDK1bt9aLL75odgzAJVB6AEiSDh06pGHDhum2225TmTJlVKZMGd1xxx0aOnSoduzYYXa8YrFo0SIlJiZq2LBhxrLHHntMNpvtqq9jx44V+DiHDx++6v4+/fTTPOu+9NJLmjlzppKSkm765wOszsaztwAsW7ZMvXv3VqlSpfToo48qJCREbm5u2rNnjxYvXqwjR47o0KFDqlmzptlRi1STJk3UqlUrzZkzx1gWFxenhISEPOs5HA4NGTJEtWrV0u7duwt8nMOHD6t27drq06eP7r///jzf3XXXXXn+OdvtdlWtWlWDBw/Wa6+9VuBjAfj/SpkdAIC5EhIS9PDDD6tmzZqKjY1V5cqV83w/adIkvf/++3Jzc+2B4e3bt+vXX3/V1KlT8ywPCwtTWFhYnmUbNmzQuXPn9Oijj97UMZs1a6a+ffv+4zpubm7q2bOnPvroI02YMEE2m+2mjglYmWv/VwzANU2ePFmZmZlasGBBvsIjSaVKldIzzzyj6tWrG8t27Nihxx57TLfeeqtKly6toKAgDRw4UCkpKXm2HT9+vGw2m/bt26e+ffvKz89P/v7+GjNmjBwOhxITE9W9e3eVK1dOQUFB+QrH2rVrZbPZ9Pnnn2vChAmqWrWqfH191bNnT6WlpSkrK0sjR45UQECAfHx8NGDAAGVlZeXZx4IFC3TPPfcoICBAXl5euuOOOzRr1qx8P+eSJUvk6empdu3aXfOfWXR0tGw2mx555JFrrnstmZmZunDhwj+uc9999+nIkSOKj4+/6eMBVsZID2Bxy5YtU926ddWqVavr3iYmJkYHDx7UgAEDFBQUpN27d2vu3LnavXu3Nm/enG80onfv3mrQoIHeeustLV++XK+//roqVKigOXPm6J577tGkSZP0ySef6Pnnn1eLFi3yFY+oqCh5e3tr1KhROnDggN577z15eHjIzc1NqampGj9+vDZv3qwPP/xQtWvX1tixY41tZ82apYYNG6pbt24qVaqUvv32Wz399NOy2+0aOnSosd6mTZsUHBwsDw+Pf/zZs7Oz9fnnn6tNmzaqVavWdf8zu5IJEybohRdekM1mU2hoqN544w2Fh4fnWy80NFSStHHjRjVt2vSmjglYmgOAZaWlpTkkOXr06JHvu9TUVEdycrLxOnfunPHd5e8vWbRokUOSY/369caycePGOSQ5nnjiCWNZTk6Oo1q1ag6bzeZ466238hzP29vb0b9/f2PZmjVrHJIcwcHBjgsXLhjL+/Tp47DZbI7OnTvnyRAWFuaoWbNmnmVXyhoREeG49dZb8yyrVq2a46GHHsq37t99++23DkmO999//5rrXs2RI0cc4eHhjlmzZjm++eYbx/Tp0x01atRwuLm5OZYtW3bFbTw9PR1PPfXUDR8TgMPB6S3AwtLT0yVJPj4++b7r0KGD/P39jdfMmTON77y9vY33f/31l06fPq3WrVtLkn755Zd8+xo0aJDx3t3dXc2bN5fD4dDjjz9uLC9fvrzq16+vgwcP5tu+X79+eUZgWrVqJYfDoYEDB+ZZr1WrVkpMTFROTs4Vs6alpen06dNq3769Dh48qLS0NOO7lJQU3XLLLfmO/XfR0dHy8PBQr169rrnu1dSoUUM//PCDhgwZoq5du2rEiBHavn27/P399dxzz11xm1tuuUWnT5++4WMCYE4PYGm+vr6SpIyMjHzfzZkzRzExMfr444/zfXfmzBmNGDFCgYGB8vb2lr+/v2rXri1JeYrEJTVq1Mjz2c/PT6VLl1alSpXyLU9NTb2u7SXlmWd0abndbs+TYePGjerYsaPKli2r8uXLy9/fXy+//PIVszqucTFrRkaGli5dqoiICFWsWPEf1y2oChUqaMCAAdq7d6/++OOPfN87HA4mMQM3iTk9gIX5+fmpcuXK2rVrV77vLs3xOXz4cL7vevXqpU2bNumFF15QkyZN5OPjI7vdrk6dOslut+db393d/bqWSVcuHldb91r7SEhI0L333qvbb79d06ZNU/Xq1eXp6anvvvtOb7/9dp6sFStWvGLhutySJUsK5aqtq7lU4s6cOaNq1arl+e7PP//MVxIBFAylB7C4Ll266IMPPtCWLVvUsmXLa66fmpqq2NhYTZgwIc+E4f379xdlzBvy7bffKisrS998802e0aI1a9bkW/f222/XoUOH/nF/n3zyiXx8fNStW7dCzyrJOLXn7++fZ/mxY8d04cIFNWjQoEiOC1gFp7cAi3vxxRdVpkwZDRw4UCdPnsz3/d9HXi6Nrvx9+fTp04ss4426Uta0tDQtWLAg37phYWHatWtXvkveL0lOTtaqVav0wAMPqEyZMjeVKzk5Od+yY8eOaf78+WrcuHG+Wwds27ZNktSmTZubOi5gdYz0ABZXr149RUdHq0+fPqpfv75xR2aHw6FDhw4pOjpabm5uxumWcuXKqV27dpo8ebKys7NVtWpVrVy58pqjJGYIDw+Xp6enunbtqieffFIZGRmaN2+eAgICdOLEiTzrdu/eXRMnTtS6deuueNn4Z599ppycnH88tTV+/HhNmDBBa9asUYcOHa663osvvmiceqtSpYoOHz6sOXPmKDMzU++8806+9WNiYlSjRg0uVwduEqUHgLp3766dO3dq6tSpWrlypebPny+bzaaaNWuqS5cuGjJkiEJCQoz1o6OjNXz4cM2cOVMOh0Ph4eH6/vvvVaVKFRN/ivzq16+vL7/8Uq+++qqef/55BQUF6amnnpK/v3++K79CQ0PVuHFjff7551csPZ988okCAgLUsWPHqx4vIyNDNptNQUFB/5grPDxcs2fP1syZM5Wamqry5curXbt2evXVV9WsWbM869rtdn311Vd6/PHHmcgM3CSevQUA/2fhwoUaOnSojh49qvLlyxd4+5YtW6pmzZr64osvCi3TkiVL9MgjjyghIeGKd8wGcP0oPQDwf+x2uxo3bqw+ffrolVdeKdC26enp8vf3V3x8fKFOOA4LC9Ndd92lyZMnF9o+Aaui9AAAAEvg6i0AAGAJTlF6oqKi1KJFC/n6+iogIEA9evTQ3r17r7ndF198odtvv12lS5dWo0aN9N133xVDWgAAUBI5RelZt26dhg4dqs2bNysmJkbZ2dkKDw9XZmbmVbfZtGmT+vTpo8cff1zbt29Xjx491KNHjyveWRYAAMAp5/QkJycrICBA69atU7t27a64Tu/evZWZmally5YZy1q3bq0mTZpo9uzZxRUVAACUEE4x0vN3lx4CWKFChauuExcXl+9+GREREYqLi7vi+llZWUpPTzdeaWlpSk5OvuYDBgEAgGtwutJjt9s1cuRItW3bVsHBwVddLykpSYGBgXmWBQYGKikp6YrrR0VFyc/Pz3iVL19eAQEBOnv2bKHmBwAAzsnpSs/QoUO1a9cuffrpp4W639GjRystLc14JSYmFur+AQCAc3Oqx1AMGzZMy5Yt0/r1643n/FxNUFBQvocjnjx58qq3f/fy8pKXl1ehZQUAACWLU4z0OBwODRs2TF9//bVWr16t2rVrX3ObsLAwxcbG5lkWExOjsLCwoooJAABKMKcY6Rk6dKiio6O1dOlS+fr6GvNy/Pz85O3tLUnq16+fqlatqqioKEnSiBEj1L59e02dOlVdunTRp59+qp9//llz58417ecAAADOyylGembNmqW0tDR16NBBlStXNl6fffaZsc7Ro0d14sQJ43ObNm0UHR2tuXPnKiQkRF9++aWWLFnyj5OfAQCAdTnlfXqKQ3p6uvz8/JSWlqZy5cqZHQcAABQxpxjpAQAAKGqUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAmUHgAAYAlOUXrWr1+vrl27qkqVKrLZbFqyZMk/rr927VrZbLZ8r6SkpOIJDAAAShynKD2ZmZkKCQnRzJkzC7Td3r17deLECeMVEBBQRAkBAEBJV8rsAJLUuXNnde7cucDbBQQEqHz58oUfCAAAuBynGOm5UU2aNFHlypV13333aePGjWbHAQAATswpRnoKqnLlypo9e7aaN2+urKwsffDBB+rQoYN++uknNWvW7IrbZGVlKSsry/icnp5eXHEBAIATKJGlp379+qpfv77xuU2bNkpISNDbb7+thQsXXnGbqKgoTZgwobgiAgAAJ1OiT29drmXLljpw4MBVvx89erTS0tKMV2JiYjGmAwAAZiuRIz1XEh8fr8qVK1/1ey8vL3l5eRVjIgAA4EycovRkZGTkGaU5dOiQ4uPjVaFCBdWoUUOjR4/WsWPH9NFHH0mSpk+frtq1a6thw4b666+/9MEHH2j16tVauXKlWT8CAABwck5Ren7++WfdfffdxufIyEhJUv/+/fXhhx/qxIkTOnr0qPH9hQsX9Nxzz+nYsWMqU6aMGjdurFWrVuXZBwAAwOVsDofDYXYIM6Snp8vPz09paWkqV66c2XEAAEARc5mJzAAAAP+E0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyB0gMAACyh1M1snJ2draSkJJ07d07+/v6qUKFCYeUCAAAoVAUe6Tl79qxmzZql9u3bq1y5cqpVq5YaNGggf39/1axZU4MHD9bWrVuLIisAAMANK1DpmTZtmmrVqqUFCxaoY8eOWrJkieLj47Vv3z7FxcVp3LhxysnJUXh4uDp16qT9+/cXVW4AAIACKVDp2bp1q9avX68tW7ZozJgxioiIUKNGjVS3bl21bNlSAwcO1IIFC5SUlKQePXroxx9/vK79rl+/Xl27dlWVKlVks9m0ZMmSa26zdu1aNWvWTF5eXqpbt64+/PDDgvwoAADAYgo0p2fRokXXtZ6Xl5eGDBly3fvNzMxUSEiIBg4cqAcffPCa6x86dEhdunTRkCFD9Mknnyg2NlaDBg1S5cqVFRERcd3HBQAA1mFzOBwOs0Nczmaz6euvv1aPHj2uus5LL72k5cuXa9euXcayhx9+WH/++adWrFhxXcdJT0+Xn5+f0tLSVK5cuZuNDQAAnNxNXb11SW5urvbs2aNdu3YZr6+//rowdn1FcXFx6tixY55lERERGjly5FW3ycrKUlZWlvE5PT29qOIBAAAnVODSc/DgQe3cuTNPwdm/f7+ys7Pl6empBg0aqFGjRkWR1ZCUlKTAwMA8ywIDA5Wenq7z58/L29s73zZRUVGaMGFCkeYCAADOq0Clp2/fvlq0aJFsNpvKlCmjzMxMdenSRWPHjlWjRo1Ur149ubu7F1XWmzJ69GhFRkYan9PT01W9enUTEwEAgOJUoKu3vvzyS7377rvKyMjQ8ePHNWzYMK1cuVJbt25VzZo1i63wBAUF6eTJk3mWnTx5UuXKlbviKI90cXJ1uXLl8rwAAIB1FKj0PPvss+rXr59Kly4tHx8fvfPOO9q4caPWrFmjhg0bXvck4psVFham2NjYPMtiYmIUFhZWLMcHAAAlT4FKT1RUlHx9ffMsCw0N1ZYtWzRixAj17t1bjzzyiJKTkwsUIiMjQ/Hx8YqPj5d08ZL0+Ph4HT16VNLFU1P9+vUz1h8yZIgOHjyoF198UXv27NH777+vzz//XM8++2yBjgsAAKyjUB44arPZNGLECP3222/KysrS7bffXqDtf/75ZzVt2lRNmzaVJEVGRqpp06YaO3asJOnEiRNGAZKk2rVra/ny5YqJiVFISIimTp2qDz74gHv0AACAqyqS+/QsX75cXbp0KezdFiru0wMAgLUUaKTn8tGWf3Kp8Bw7dqzgiQAAAIpAgUpPixYt9OSTT/7jU9TT0tI0b948BQcH66uvvrrpgAAAAIWhQPfp6dKli3x8fHTfffepdOnSCg0NVZUqVVS6dGmlpqbqt99+0+7du9WsWTNNnjxZ999/f1HlBgAAKJACzenx9PRUYmKifH195e/vrz59+iglJUXnz59XpUqV1LRpU0VERCg4OLgoMxcK5vQAAGAtBRrpqVKliuLj4xUREaHz58/rzTffVEBAQFFlAwAAKDQFmtPz3HPPqWvXrrrrrrtks9n0ySefaOvWrTp//nxR5QMAACgUBb5kfceOHfr22281ZswY3XrrrTp8+LBsNpvq1q2rkJAQNWnSRCEhIercuXNRZS4UnN4CAMBabvg+PfXq1VNcXJzKli2rHTt2GHdUjo+P165du3T27NnCzlqoKD0AAFhLkdyc0OFwyGazFfZuCxWlBwAAaymUx1D8nbMXHgAAYD1FUnoAAACcDaUHAABYAqUHAABYAqUHAABYAqUHAABYAqUHAABYAqUHAABYAqUHAABYAqUHAABYAqUHAABYAqUHAABYAqUHAABYAqUHAABYAqUHAABYAqUHAABYAqUHAABYAqUHAABYQimzA8B5HT9+XJs2bdKpU6fk6+urBg0aqFmzZnJzoysDAEoeSg/yWbduncaPH6+1a9fm+y4oKEjDhg3T8OHDVa5cueIPBwDADbI5HA6H2SHMkJ6eLj8/P6WlpfGX9/85f/68IiMjNXv27GuuW61aNc2bN0+dOnUqhmQAANw8Sg+lR5L0559/qmvXrtqwYYOxrF69eurVq5fq1q2r1NRUrVmzRsuXL5fdbpck2Ww2TZw4US+//LJsNptZ0QEAuC6UHkqPzp07p3vuuUc//fSTJMnb21tTp07VE088IXd39zzrHjhwQEOHDtXKlSuNZcOHD9c777xD8QEAODVmpFqcw+HQf/7zH6PwVKpUST/++KOeeuqpfIVHkurWrasVK1YoKirKKDnvvfeeXnjhhWLNDQBAQVF6LG769OlavHixJMnX11erVq1SaGjoP25js9k0atQoLViwwCg+U6dO1bvvvlvkeQEAuFGc3rLw6a1ff/1VLVq0UHZ2tiRp+fLluv/++wu0j3nz5umJJ56QJLm7u2vdunVq27ZtoWcFAOBmMdJjUbm5uRo8eLBReF544YUCFx5JGjx4sEaNGmXss3fv3kpOTi7UrAAAFAZKj0XNmTNHW7dulSQ1aNBAEydOvOF9TZw4Ue3bt5ckHTt2TP/5z3+MK7wAAHAWlB4LSktL07hx44zPc+fOlZeX1w3vr1SpUoqOjlZAQIAk6Ycffriue/0AAFCcKD0WFBUVpdOnT0uSHnnkEd155503vc8qVaro448/Nj6PGjVKiYmJN71fAAAKC6XHYlJSUjRjxgxJkpeXl954441C2/d9992nQYMGSZLOnj2rp59+WhadJw8AcEKUHouZOXOmMjMzJUmDBg1SrVq1CnX/kydPVlBQkCRp2bJl+uKLLwp1/wAA3CguWbfQJeuZmZmqWbOmUlJS5O7urgMHDhR66ZGkr776Sj179pR08bTXvn37VLZs2UI/DgAABcFIj4V88MEHSklJkST16dOnSAqPJD344IPq0qWLJOn48eOaMmVKkRwHAICCYKTHIiM9Fy5cUN26dY3JxTt37lRwcHCRHW/v3r0KDg5WTk6OvL29tW/fPlWrVq3IjgcAwLUw0mMRixcvNgpP165di7TwSFL9+vU1dOhQSdL58+f18ssvF+nxAAC4FkqPRXzwwQfG+2effbZYjjl27FhVqFBBkrRw4UJt27atWI4LAMCVUHosICEhQbGxsZIuPiW9Q4cOxXLcChUq5LkJ4pgxY4rluAAAXAmlxwL++9//Gu8HDRpkPBm9ODz55JOqWbOmJOn777/Xxo0bi+3YAABcjtLj4rKzs7VgwQJJFx8X0b9//2I9vpeXl8aOHWt8fuWVV7hhIQDAFJQeF/fdd98pKSlJ0sUJzJduHFic+vXrp3r16kmS1q1bZ5xqAwCgOFF6XNzlp7YGDx5sSoZSpUppwoQJxudx48Yx2gMAKHbcp8eF79OTmpqqwMBAZWdnq2rVqjpy5Ijc3d1NyWK329W4cWPt3r1bkrRmzZpim1ANAIDESI9L+/rrr5WdnS1J6tWrl2mFR5Lc3Nzy3KvnzTffNC0LAMCaKD0u7NNPPzXe9+7d28QkF/Xq1Uu33nqrJCkmJkZbt241OREAwEooPS4qOTlZq1evliTVqlVLLVu2NDnRxbk9o0aNMj4z2gMAKE6UHhf11VdfKTc3V9LFUZ7ivDfPP+nXr5+qVq0qSVqyZIl27dplciIAgFVQelyUs53ausTLy0vPP/+88fmtt94yMQ0AwEq4essFr946fvy4qlWrJofDodtuu0179uxxmpEeScrMzFStWrV0+vRpubm5ad++fapTp47ZsQAALo6RHhf0zTffGPfB6dWrl1MVHkkqW7asRo4cKenipexTp041NxAAwBIoPS7om2++Md4/8MADJia5uqFDh6ps2bKSpA8//FApKSkmJwIAuDpKj4vJyMgwHvNQtWpVNW3a1OREV1a+fHk9/vjjkqTz589r9uzZJicCALg6So+LWblypS5cuCBJ6tatm9Od2rrcyJEj5eZ28V/BGTNmKCsry+REAABXRulxMZef2urWrZuJSa6tdu3aevDBByVJSUlJio6ONjkRAMCVcfWWC129lZOTo6CgIKWkpMjHx0enT5+Wl5eX2bH+0ebNmxUWFiZJCg4O1o4dO5x6dAoAUHIx0uNC4uLijAnBnTp1cvrCI0mtW7dWmzZtJEm7du3SypUrTU4EAHBVlB4XUpJObV3uueeeM95z+ToAoKhwesuFTm/dcccd+v333+Xm5qZTp06pYsWKZke6Lrm5ubrtttt08OBBSdKvv/6qxo0bm5wKAOBqGOlxEUePHtXvv/8uSWrVqlWJKTyS5O7ubtysUJKmTZtmXhgAgMui9LiIH374wXjfqVMnE5PcmAEDBqh8+fKSpOjoaJ04ccLcQAAAl0PpcRErVqww3pfE0uPj46MhQ4ZIkrKzszVjxgyTEwEAXA1zelxgTk92drYqVaqk9PR0VahQQadOnZK7u7vZsQrs2LFjqlWrlnJyclShQgUlJiaqTJkyZscCALgIRnpcwE8//aT09HRJUnh4eIksPNLFx2b07t1bknTmzBktXLjQ5EQAAFdC6XEBJX0+z+WeffZZ4/306dNlt9tNTAMAcCWUHhdw+Xye8PBwE5PcvNDQULVr106StGfPnjyFDgCAm0HpKeGSk5O1bds2SVJISIgqV65scqKbd/loz9tvv21iEgCAK6H0lHBr1qzRpbnoERERJqcpHF27dtWtt94qSYqJidGuXbtMTgQAcAWUnhJuzZo1xvt77rnHxCSFx93dXSNGjDA+T58+3bwwAACXwSXrJfyS9dtvv1179+5VqVKllJqaKh8fH7MjFYqzZ8+qWrVqSk9Pl5eXl44ePaqAgACzYwEASjBGekqw48ePa+/evZKkli1bukzhkSRfX18NHjxYkpSVlaVZs2aZnAgAUNJRekqwtWvXGu/vvvtu84IUkeHDh8vN7eK/ou+//77++usvkxMBAEoypyo9M2fOVK1atVS6dGm1atVKW7Zsueq6H374oWw2W55X6dKlizGt+S6fz+OKpadmzZp66KGHJEmnTp3SokWLTE4EACjJnKb0fPbZZ4qMjNS4ceP0yy+/KCQkRBERETp16tRVtylXrpxOnDhhvI4cOVKMic13qfR4enqqTZs2JqcpGpGRkcb7t99+WxadggYAKAROU3qmTZumwYMHa8CAAbrjjjs0e/ZslSlTRvPnz7/qNjabTUFBQcYrMDCwGBObKzExUQkJCZKk1q1by9vb2+RERaN169Zq3bq1JGnnzp1avXq1yYkAACWVU5SeCxcuaNu2berYsaOxzM3NTR07dlRcXNxVt8vIyFDNmjVVvXp1de/eXbt3777qullZWUpPT8/zKslc/dTW5bhZIQCgMDhF6Tl9+rRyc3PzjdQEBgYqKSnpitvUr19f8+fP19KlS/Xxxx/LbrerTZs2+uOPP664flRUlPz8/IxX9erVC/3nKE5WKj0PPvigatSoIUlavny5ccUaAAAF4RSl50aEhYWpX79+atKkidq3b6/FixfL399fc+bMueL6o0ePVlpamvFKTEws5sSF61LpKV26tHH6x1WVKlVKw4cPNz5zs0IAwI1witJTqVIlubu76+TJk3mWnzx5UkFBQde1Dw8PDzVt2lQHDhy44vdeXl4qV65cnldJdejQIWPSdps2beTl5WVyoqI3aNAglS1bVpL0v//9TykpKSYnAgCUNE5Rejw9PRUaGqrY2Fhjmd1uV2xsrMLCwq5rH7m5udq5c6dLPHDzWqx0auuS8uXLa+DAgZKk8+fPa+7cuSYnAgCUNE5ReqSLlybPmzdP//vf//T777/rqaeeUmZmpgYMGCBJ6tevn0aPHm2s/9prr2nlypU6ePCgfvnlF/Xt21dHjhzRoEGDzPoRio0VS48kjRgxQjabTZI0Y8YMXbhwweREAICSpJTZAS7p3bu3kpOTNXbsWCUlJalJkyZasWKFMbn56NGjxt15JSk1NVWDBw9WUlKSbrnlFoWGhmrTpk264447zPoRioXD4TBKT5kyZdSiRQuTExWfOnXqqHv37lqyZImOHz+uL774Qo8++qjZsQAAJQQPHC1hDxzdv3+/brvtNklSeHi4fvjhB5MTFa/169erffv2kqTQ0FBt3brVGP0BAOCfOM3pLVyfy09tdejQwbwgJrnrrrvUrFkzSdK2bdu0YcMGkxMBAEoKSk8JY9X5PJfYbDZuVggAuCGc3ipBp7ccDocqV66skydPysfHR2fOnJGHh4fZsYrdhQsXVKtWLZ04cUI2m0379+9XnTp1zI4FAHByjPSUIHv27DHuZXTXXXdZsvBIF29xMGzYMEkXiyA3KwQAXA9KTwmydu1a470VT21d7sknn1SZMmUkSf/973+VnJxsciIAgLOj9JQgVp/Pc7mKFStq8ODBki7erPDdd981OREAwNkxp6eEzOlxOBwKDAxUcnKy/Pz8lJKSInd3d7NjmSoxMVG33nqrcnJyVL58eR09elS+vr5mxwIAOClGekqI3bt3G6dw2rVrZ/nCI0nVq1dX3759JUl//vnnVR82CwCAROkpMTi1dWUvvviicXPCadOmKSsry+REAABnRekpISg9V9agQQP16NFDknTixAktXLjQ3EAAAKfFnJ4SMKfHbrfL399fZ86cUYUKFZScnJznOWRWt2XLFrVq1UqSVK9ePf3++++c/gMA5MPfnCXAjh07dObMGUlS+/btKTx/07JlS91zzz2SLj6bbPHixSYnAgA4I/72LAE4tXVto0aNMt5HRUXJogOYAIB/QOkpASg919axY0eFhoZKkrZv366VK1eanAgA4GwoPU4uNzdX69evlyT5+/urYcOGJidyTjabLc9oz8SJExntAQDkQelxctu3b1daWpokqUOHDsbl2cjvgQce0B133CFJ2rhxo1atWmVyIgCAM6H0ODlObV0/d3d3jR071vg8fvx4RnsAAAZKj5Oj9BRMz549jdGeTZs2KSYmxuREAABnQelxYtnZ2frxxx8lSUFBQapfv77JiZyfu7u7xo0bZ3xmtAcAcAmlx4lt27ZNGRkZki6O8jCf5/r07NnTmPAdFxfHlVwAAEmUHqfGqa0b4+bmxmgPACAfSo8To/TcuIceekjBwcGSpM2bN+u7774zOREAwGyUHid14cIFbdy4UZJUrVo11alTx+REJYubm5vGjx9vfB49erRyc3PNCwQAMB2lx0lt2bJF586dk8R8nhv14IMPqmXLlpKknTt3Kjo62uREAAAzUXqcVGxsrPH+0sM0UTA2m01vvfWW8XnMmDHKysoyMREAwEyUHid1eem59957TUxSst19993q1KmTJOnIkSOaNWuWyYkAAGaxOSx6WUt6err8/PyUlpamcuXKmR0nj8zMTN1yyy3Kzs5WvXr1tG/fPrMjlWjx8fFq2rSpJKlixYpKSEiQn5+fyakAAMWNkR4ntGHDBmVnZ0vi1FZhaNKkiR599FFJUkpKiqZMmWJyIgCAGSg9TohTW4Vv4sSJ8vDwkCRNmzZNiYmJJicCABQ3So8Turz0cH+ewlG7dm09/fTTkqTz58/rpZdeMjkRAKC4MafHyeb0nDlzRpUqVZLD4VBISIji4+PNjuQyUlNTVa9ePaWkpEiSfvzxR915550mpwIAFBdGepzM2rVrjUcmcGqrcN1yyy164403jM/PPPMMNywEAAuh9DgZ5vMUrUGDBikkJESStH37ds2fP9/kRACA4sLpLSc7vdWgQQPt2bNHpUqV0pkzZ+Tr62t2JJezfv16tW/fXpJUqVIl7d+/X+XLlzc3FACgyDHS40SOHTumPXv2SJJatmxJ4Ski7dq1U69evSRJp0+fzvNEdgCA66L0OJHVq1cb7zm1VbSmTJkib29vSdKMGTO0detWkxMBAIoapceJXF56uClh0apRo4YxwmO32zV48GDjhpAAANdE6XESDodDq1atkiR5e3srLCzM5ESuLzIy0pjU/Ouvv+rtt982OREAoChRepzEb7/9pj/++EPSxTknXl5eJidyfR4eHpo3b55sNpskafz48UpISDA5FQCgqFB6nMSKFSuM9507dzYxibW0aNFCzzzzjKSLd2oeMmSILHpBIwC4PEqPk/j++++N9506dTIxifW8/vrrqlGjhiRp1apVWrBggcmJAABFgdLjBDIyMvTjjz9KkmrVqqXbbrvN5ETW4uPjo/fff9/4PGLECB08eNDERACAokDpcQJr167VhQsXJF0c5bk0xwTFp0uXLhowYICkiyW0X79+PKICAFwMpccJXD6fh1Nb5nnnnXdUu3ZtSdLGjRs1adIkkxMBAAoTj6FwgsdQ1K1bVwkJCfLw8FBKSgp3YjbRxo0b1a5dO9ntdpUqVUqbN29WaGio2bEAAIWAkR6THThwwLhM+s4776TwmKxt27YaNWqUJCknJ0d9+/ZVZmamyakAAIWB0mOyZcuWGe8jIiJMTIJLxo0bp2bNmkmS9uzZoyeeeILL2AHABVB6TLZ06VLjfdeuXU1Mgks8PT0VHR0tHx8fSVJ0dHSeq7sAACUTpcdEKSkpxqXqdevWVYMGDUxOhEvq16+f5349zz77rDZv3mxiIgDAzaL0mGj58uXGZdHdu3fnUnUn07NnT0VGRkqSsrOz9e9//1vJyckmpwIA3ChKj4kuP7XVvXt3E5Pgat566y3dddddkqQ//vhDvXr1Mu6pBAAoWbhk3aRL1v/66y9VqlRJmZmZqlSpkpKSkuTu7l7sOXBtJ06cULNmzZSUlCRJ6t+/vxYsWMDIHACUMIz0mCQ2Nta4FLpr164UHidWuXJlLV68WKVLl5Yk/e9//9Prr79ucioAQEFRekyyZMkS4z2ntpxfWFiYFi5caHweO3asPv74YxMTAQAKitJjguzsbGM+T+nSpXXfffeZnAjXo2fPnpoyZYrxeeDAgYqNjTUxEQCgICg9JoiNjTWuAuratavKlCljciJcr+eee05PPfWUpIvltVu3bsZtBwAAzo3SY4Lo6Gjj/SOPPGJiEhSUzWbTu+++q27dukmSzp07p/vvv597+ABACcDVW8V89da5c+cUGBiojIwM+fn56eTJk/Ly8iq246NwZGVlqUePHlqxYoUkyc/PT6tWrVLz5s1NTgYAuBpGeorZ8uXLlZGRIUl66KGHKDwllJeXlxYvXqx7771XkpSWlqbw8HD99NNPJicDAFwNpaeYcWrLdXh7e2vp0qVq166dJCk1NVX33HOPMfoDAHAulJ5ilJqaqu+++06SFBQUpA4dOpgbCDetbNmyWrZsme6++25JF09fdu3alcvZAcAJUXqK0aJFi4xHGDz88MPckNBF+Pr66vvvv1fPnj0lSTk5OfrPf/6jKVOmyKJT5gDAKVF6ionD4dCsWbOMz4899ph5YVDovLy89Omnn+rpp582lr344ovq16+fzp07Z2IyAMAllJ5ismHDBu3atUuS1KZNG4WEhJicCIXN3d1dM2bM0GuvvWYs+/jjj3XnnXfq8OHD5gUDAEii9BSby0d5Lt3cDq7HZrNpzJgx+vLLL1W2bFlJ0vbt29W8eXMmOAOAybhPTzHcp+fkyZOqXr26srOzValSJSUmJhoPr4Tr2r17t3r06KEDBw4Yy4YNG6ZJkyZxF24AMAEjPcVg/vz5ys7OliQ9/vjjFB6LaNiwobZu3aouXboYy2bMmKHQ0FBt27bNxGQAYE2UniKWlZWlmTNnSrp46uPJJ580ORGKU/ny5fXtt9/qvffeM8runj171Lp1a7300kvGjSoBAEWP0lPEPvroIx07dkyS1K1bN9WuXdvkRChuNptNw4YN0/bt2xUaGirp4mXtkydPVoMGDbR48WIubQeAYkDpKULZ2dmaNGmS8fmVV14xMQ3MdvvttysuLk7jx4+Xp6enJOmPP/7QQw89pE6dOmn79u0mJwQA10bpKULz589XQkKCJKljx45q0aKFyYlgNg8PD40bN067du1SRESEsXzlypVq1qyZHn74Ye3bt8/EhADgurh6q4iu3srMzFS9evV04sQJSdLmzZvVqlWrQj8OSi6Hw6HFixcrMjJSR48eNZa7u7vr0Ucf1fPPP69GjRqZmBAAXAsjPUXkzTffNArPAw88QOFBPjabTQ899JD27dun6dOny9/fX5KUm5urjz76SI0bN1anTp20atUq5vwAQCFgpKcIRnr27t2rRo0aKTs7Wx4eHtq5c6fq169fqMeA68nIyND06dM1bdo0paam5vmuXr16GjRokPr166egoCCTEgJAycZITyHLyclR//79jfvyvPDCCxQeXBcfHx+9+uqrOnr0qN59913VqlXL+G7//v166aWXVK1aNXXv3l2LFi3S2bNnzQsLACUQIz2FPNIzfvx4TZgwQZJUt25d/frrr9x9FzckJydHX3/9tWbNmqU1a9bk+97Ly0udO3fWAw88oPDwcEaAAOAaKD2FWHq++eYbde/eXZLk5uamjRs3qnXr1oWyb1hbQkKC5s+frwULFhhzxf4uJCRE4eHhioiIUNu2bbnzNwD8DaWnkErPhg0bFBERoXPnzkm6OJF59OjRN71f4HK5ubnasGGDPv/8c3311Vc6efLkFdfz8PBQs2bNFBYWZryqV69ezGkBwLlQegqh9Pzwww/q2bOn8UiB3r17a9GiRbLZbIURFbii3NxcxcXFacWKFfrhhx+0bdu2f7zKy9/fX40aNVLjxo2N/61fv758fX2LMTUAmMepSs/MmTM1ZcoUJSUlKSQkRO+9955atmx51fW/+OILjRkzRocPH1a9evU0adIk3X///dd1rMIoPTk5OYqKitKECROUm5srSYqIiNDSpUvl5eV1Q/sEbtTp06cVExOjmJgYbdq0SXv37r2u7fz9/VWnTh3VrVtXderUUZ06dVSzZk1VrVpVVapUkbe3dxEnB4Di4TSl57PPPlO/fv00e/ZstWrVStOnT9cXX3yhvXv3KiAgIN/6mzZtUrt27RQVFaV//etfio6O1qRJk/TLL78oODj4mse7mdKTnZ2txYsXa/z48dqzZ4+x/IEHHtAnn3zCXxJwCmfOnNHmzZsVFxenn376STt27Ljq6bB/Ur58eaMABQUFqWLFiqpQoYIqVqyY7+Xr6ysfHx95eHgUwU8EADfHaUpPq1at1KJFC82YMUOSZLfbVb16dQ0fPlyjRo3Kt37v3r2VmZmpZcuWGctat26tJk2aaPbs2dc83rVKj8PhUG5urs6ePavTp08rMTFR8fHx2rZtm3744QelpKQY67q7u+uVV17RuHHj5ObGXQDgvE6ePKmdO3carwMHDighIUHHjx8v1ON4eXnJx8fHKEGXvy9btqxKly4tLy+vK748PT3zLfPw8JC7u7tKlSold3f3PK+/L7vaOm5ubrLZbAV6SbqudQCUDKXMDiBJFy5c0LZt2/JM/HVzc1PHjh0VFxd3xW3i4uIUGRmZZ1lERISWLFlyxfWzsrKUlZVlfE5PT5ckVa5cWQ6HQ3a7Xbm5ubLb7bLb7dedvU2bNpo+fTrP1UKJEBgYqMDAQHXs2DHP8nPnzungwYNKSEhQQkKC/vjjDx07dkzHjx/X8ePHdezYsTx/fq7l0p+3y//PgSu7nuJ0s/svjIxkKJx9wHxpaWk3tJ1TlJ7Tp08rNzdXgYGBeZYHBgbmOX10uaSkpCuun5SUdMX1L829+btLV1sVhI+Pj+6//349+eSTuvvuu/lDhBKvTJkyCg4OvuqpYYfDodTUVJ06dUopKSnG68yZM3nenz17VhkZGcrIyMjz/kb+nJUkDoeDR4UAJYBTlJ7iMHr06DwjQ+np6apevboaNGggDw8Pubm5yc3NTe7u7sZ7Nzc3lSlTRv7+/goICFBwcLCaNGmi4OBgJirDUmw2mypUqKAKFSrc0Pa5ubnKzMw0StBff/2lCxcuGCNC13rl5OQoNzdXubm5ed5f6fPV1rHb7UY5uZ6XpAKt/0/7uBmFUabIUHj7QMnmFKWnUqVKcnd3zzfJ8uTJk1e9y2xQUFCB1r80N+DvNm/eXCRPWQfw/7m7u6tcuXL8WQNgKqeYdevp6anQ0FDFxsYay+x2u2JjYxUWFnbFbcLCwvKsL0kxMTFXXR8AAFibU4z0SFJkZKT69++v5s2bq2XLlpo+fboyMzM1YMAASVK/fv1UtWpVRUVFSZJGjBih9u3ba+rUqerSpYs+/fRT/fzzz5o7d66ZPwYAAHBSTlN6evfureTkZI0dO1ZJSUlq0qSJVqxYYUxWPnr0aJ7Lwdu0aaPo6Gi9+uqrevnll1WvXj0tWbLkuu7RAwAArMdp7tNT3IrqKesAAMA5Oc1IT3G71PUu3a8HAACUHL6+vgW+ZYxlS8+lm6bx5GkAAEqeU6dOyd/fv0DbWLb0XLrfyNGjR+Xn52dyGhS1S/dlSkxM5HSmBfD7thZ+39Zy6fft6elZ4G0tW3ouTYr28/PjD4mFcK8Ya+H3bS38vq3lRp6G4BT36QEAAChqlB4AAGAJli09Xl5eGjduHM/Qsgh+39bC79ta+H1by838vi17nx4AAGAtlh3pAQAA1kLpAQAAlkDpAQAAlkDpAQAAlmDZ0jNz5kzVqlVLpUuXVqtWrbRlyxazI6EIREVFqUWLFvL19VVAQIB69OihvXv3mh0LxeStt96SzWbTyJEjzY6CInLs2DH17dtXFStWlLe3txo1aqSff/7Z7FgoArm5uRozZoxq164tb29v1alTRxMnTlRBrseyZOn57LPPFBkZqXHjxumXX35RSEiIIiIidOrUKbOjoZCtW7dOQ4cO1ebNmxUTE6Ps7GyFh4crMzPT7GgoYlu3btWcOXPUuHFjs6OgiKSmpqpt27by8PDQ999/r99++01Tp07VLbfcYnY0FIFJkyZp1qxZmjFjhn7//XdNmjRJkydP1nvvvXfd+7DkJeutWrVSixYtNGPGDEmS3W5X9erVNXz4cI0aNcrkdChKycnJCggI0Lp169SuXTuz46CIZGRkqFmzZnr//ff1+uuvq0mTJpo+fbrZsVDIRo0apY0bN+rHH380OwqKwb/+9S8FBgbqv//9r7HsoYcekre3tz7++OPr2oflRnouXLigbdu2qWPHjsYyNzc3dezYUXFxcSYmQ3FIS0uT9P8fOAvXNHToUHXp0iXPn3O4nm+++UbNmzfXv//9bwUEBKhp06aaN2+e2bFQRNq0aaPY2Fjt27dPkvTrr79qw4YN6ty583Xvw3IPHD19+rRyc3MVGBiYZ3lgYKD27NljUioUB7vdrpEjR6pt27YKDg42Ow6KyKeffqpffvlFW7duNTsKitjBgwc1a9YsRUZG6uWXX9bWrVv1zDPPyNPTU/379zc7HgrZqFGjlJ6erttvv13u7u7Kzc3VG2+8oUcfffS692G50gPrGjp0qHbt2qUNGzaYHQVFJDExUSNGjFBMTIxKly5tdhwUMbvdrubNm+vNN9+UJDVt2lS7du3S7NmzKT0u6PPPP9cnn3yi6OhoNWzYUPHx8Ro5cqSqVKly3b9vy5WeSpUqyd3dXSdPnsyz/OTJkwoKCjIpFYrasGHDtGzZMq1fv17VqlUzOw6KyLZt23Tq1Ck1a9bMWJabm6v169drxowZysrKkru7u4kJUZgqV66sO+64I8+yBg0a6KuvvjIpEYrSCy+8oFGjRunhhx+WJDVq1EhHjhxRVFTUdZcey83p8fT0VGhoqGJjY41ldrtdsbGxCgsLMzEZioLD4dCwYcP09ddfa/Xq1apdu7bZkVCE7r33Xu3cuVPx8fHGq3nz5nr00UcVHx9P4XExbdu2zXcLin379qlmzZomJUJROnfunNzc8tYWd3d32e32696H5UZ6JCkyMlL9+/dX8+bN1bJlS02fPl2ZmZkaMGCA2dFQyIYOHaro6GgtXbpUvr6+SkpKkiT5+fnJ29vb5HQobL6+vvnma5UtW1YVK1ZkHpcLevbZZ9WmTRu9+eab6tWrl7Zs2aK5c+dq7ty5ZkdDEejataveeOMN1ahRQw0bNtT27ds1bdo0DRw48Lr3YclL1iVpxowZmjJlipKSktSkSRO9++67atWqldmxUMhsNtsVly9YsECPPfZY8YaBKTp06MAl6y5s2bJlGj16tPbv36/atWsrMjJSgwcPNjsWisDZs2c1ZswYff311zp16pSqVKmiPn36aOzYsfL09LyufVi29AAAAGux3JweAABgTZQeAABgCZQeAABgCZQeAABgCZQeAABgCZQeAABgCZQeAABgCZQeAABgCZQeAC5n1KhR8vLy0iOPPGJ2FABOhDsyA3A5aWlpWrhwoYYPH679+/erbt26ZkcC4AQY6QHgcvz8/PT444/Lzc1NO3fuNDsOACdB6QHgknJyclSmTBnt2rXL7CgAnASlB4BLevXVV5WRkUHpAWBgTg8Al7Nt2za1adNG9913nw4dOqTdu3ebHQmAE6D0AHApdrtdLVu2VPv27dWqVSv17dtXmZmZ8vDwMDsaAJNxeguAS3nvvfd0+vRpvfbaa2rUqJGys7O1Z88es2MBcAKUHgAu49ixYxozZoxmzpypsmXLql69evLy8mJeDwBJlB4ALuSZZ55R586d1aVLF0lSqVKl1KBBA0oPAElSKbMDAEBhWLZsmVavXq3ff/89z/JGjRpRegBIYiIzAACwCE5vAQAAS6D0AAAAS6D0AAAAS6D0AAAAS6D0AAAAS6D0AAAAS6D0AAAAS6D0AAAAS6D0AAAAS6D0AAAAS6D0AAAAS6D0AAAAS/h/NnoZ+zuIHh8AAAAASUVORK5CYII=\">"},"metadata":{}}],"execution_count":2},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"A51E37A95B5345C6BA3F1A0A76D8CABF","notebookId":"65251b1ab0f4556a3142003a","trusted":true},"source":"y = np.linspace(0, 25, 26)        # 假设明天获得金牌的平均数量 0到25次\n\nλs = [5,10,15]                      # 生成三种每天获得枚金牌平均数量 lambda\n\nfig, axes = plt.subplots(1,3,figsize=(15,5),sharex=True,sharey=True)   # 生成三个子图\n\nfor i,λ in enumerate(λs):\n    y_p = st.poisson.pmf(y,λ)     # 在给定λ下，y发生的概率\n    axes[i].stem(y,y_p,           # 画图\n                 linefmt='black',\n                 bottom=-1)\n    axes[i].set_title(f\"poisson($\\lambda=${λ})\")  # 设置标题\n    axes[i].set_xlabel('$y$')                     # 设置x轴标题\n\n\naxes[0].set_ylabel('$f(y)$')                      # 设置y轴标题\nplt.ylim(0, 0.5)                                  # 设置y轴范围\nplt.yticks(np.arange(0, 0.5, 0.1))                # 设置y轴刻度\nsns.despine()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 1500x500 with 3 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/3FAF373AE01F4B8EB9669C1ECEE11F14/s2326h5d7a.png\">"},"metadata":{}}],"execution_count":3},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"E5F26085AD824826BB4C4B8BAD948C37","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**计算Poisson的似然函数**  \n\n\n🤔如果我们一共获得了四天的数据，此时的似然函数怎么写？  \n\n我们将四天内中国获得的金牌数量$(Y_1, Y_2, \\ldots, Y_n)$组合为一个向量$\\vec{y} = (y_1,y_2,\\ldots,y_n)$  \n\n在给定的$\\lambda$下，$\\vec{y}$发生的可能性可以表示为：  \n\n$$  \nf(\\vec{y} | \\lambda) = \\prod_{i=1}^n f(y_i | \\lambda) = f(y_1 | \\lambda) \\cdot f(y_2 | \\lambda) \\cdot \\cdots \\cdot f(y_n | \\lambda)  \n$$  \n\n> 注意：由于每天获得金牌的数量是互不影响的，所以我们可以直接相乘：$P(A \\cap B) = P(A)P(B)$  \n"},{"cell_type":"markdown","metadata":{"id":"886B263892CC43189EA9B9B8405FF561","notebookId":"65251b1ab0f4556a3142003a","runtime":{"status":"default","execution_status":null,"is_visible":false},"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"}},"source":"我们可以把公式代入并进行整理，最后可得到  \n\n$$  \n\\begin{split}  \nf(\\vec{y} | \\lambda) =\\frac{\\lambda^{\\sum y_i}e^{-n\\lambda}}{\\prod_{i=1}^n y_i!} \\\\  \n\\end{split}  \n$$  \n\n其中n为数据的数量(在这个例子中是收集多少天的数据)，当n很大时，$\\prod y_i!$的计算将变得麻烦，同样的我们可以暂时忽略分母  \n\n$$  \nL(\\lambda | \\vec{y}) = \\frac{\\lambda^{\\sum y_i}e^{-n\\lambda}}{\\prod_{i=1}^n y_i!} \\propto \\lambda^{\\sum y_i}e^{-n\\lambda} \\;\\; \\text{ for } \\lambda > 0  \n$$"},{"cell_type":"markdown","metadata":{"id":"B0597CD0EE994B21AA75F1A53E5FDD87","notebookId":"65251b1ab0f4556a3142003a","runtime":{"status":"default","execution_status":null,"is_visible":false},"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"}},"source":"**具体计算**  \n\n$$  \n\\begin{split}  \nf(\\vec{y} | \\lambda)  \n& = \\frac{\\lambda^{y_1}e^{-\\lambda}}{y_1!} \\cdot \\frac{\\lambda^{y_2}e^{-\\lambda}}{y_2!} \\cdots \\frac{\\lambda^{y_n}e^{-\\lambda}}{y_n!} \\\\  \n& = \\frac{\\left[\\lambda^{y_1}\\lambda^{y_2} \\cdots \\lambda^{y_n}\\right] \\left[e^{-\\lambda}e^{-\\lambda} \\cdots e^{-\\lambda}\\right]}{y_1! y_2! \\cdots y_n!} \\\\  \n& =\\frac{\\lambda^{\\sum y_i}e^{-n\\lambda}}{\\prod_{i=1}^n y_i!} \\\\  \n\\end{split}  \n$$"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"140797C9D58143F28DF58E856EEF774B","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**代入具体数据**  \n\n假设四天内中国获得金牌的数量分别是20，13，13，8，我们可以将其组合成一个向量$\\vec{y} = (y_1,y_2,y_3,y_4) = (20，13，13，8)$  \n\n代入公式：  \n$$  \nL(\\lambda | \\vec{y}) = \\frac{\\lambda^{\\sum y_i}e^{-n\\lambda}}{\\prod_{i=1}^n y_i!} \\;\\;\\;\\;\\;\\;\\;\\;\\;\\sum_{i=1}^4y_i = 20 + 13 + 13 + 8 = 54  \n$$  \n\n可得：  \n$$  \nL(\\lambda | \\vec{y}) = \\frac{\\lambda^{54}e^{-4\\lambda}}{20!\\times13!\\times13!\\times8!} \\propto \\lambda^{54}e^{-4\\lambda} \\;\\;\\;\\; \\text{ for } \\lambda > 0  \n$$  \n\n此时，输入不同的参数$\\lambda$，就可以计算得到对应参数的似然值。"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"ABDD78BB34B24D32BAFC931AA13BA557","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"首先，我们尝试绘制数据的频率分布图，以确定数据的形态。  \n\n其次，我们假设中国每天平均获得12枚金牌($\\lambda=12$)，并尝试计算对应数据的似然值。"},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"fragment"},"trusted":true,"id":"F7DB7053FE674477B0909182E87E5D90","notebookId":"65251b1ab0f4556a3142003a"},"source":"# 绘制数据的分布\n\n#------------------------------------------------\n#                你可以更改以下数据，体会不同数据形成的分布，例如 observed_data = [3,3,2,2,2,1,...]\n#------------------------------------------------\nobserved_data = [20,13,13,8]   # 假设连续四天中国获得金牌的数量\n\nplt.hist(observed_data)\nplt.title('Histogram of Observations')\nplt.xlabel('Number of Gold Medal')\nplt.ylabel('Frequency')\nplt.show()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 640x480 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/1E65914BB68140E4B7763D517C0DAB09/s23291dovs.png\">"},"metadata":{}}],"execution_count":17},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"C1B5866FF6CD4EF68F8B20A52EC8A4CF","notebookId":"65251b1ab0f4556a3142003a","trusted":true},"source":"# 计算数据的似然 \n\n#------------------------------------------------\n#                你可以尝试不同的lambda 参数值\n#------------------------------------------------\nλ = 12 # 假设每天获得金牌的平均数量为12\n\nprint('The likelihood of observing the data given the parameter lambda is: ', st.poisson.pmf(observed_data,λ))","outputs":[{"output_type":"stream","name":"stdout","text":"The likelihood of observing the data given the parameter lambda is:  [0.00968203 0.10557038 0.10557038 0.06552328]\n"}],"execution_count":18},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"CD0779EB0A424B799C38AB95E2D12CA4","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**Poisson似然函数图示**  \n\n如何画出Poisson的似然函数图？  \n\n- 我们需要横轴为不同的$\\lambda$，纵轴为在不同的$\\lambda$下，$\\vec{y}$发生的似然$L(\\lambda | \\vec{y})$："},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"fragment"},"trusted":true,"id":"20D5F657F7404D92AC4A2C3F94159BB1","notebookId":"65251b1ab0f4556a3142003a"},"source":"# 定义似然函数，由于该似然函数不能通过调用已有的函数获得，因此我们根据得到的公式定义一个\n\ndef poisson_likelihood(y):\n    lambdas = np.arange(0, 25, 0.1)\n    \n    # Calculate the Poisson likelihood for each lambda\n    likelihood_values = np.exp(-len(y) * lambdas) * np.power(lambdas, np.sum(y)) / np.prod([np.math.factorial(val) for val in y])\n\n    return lambdas, likelihood_values\n\n#传入需要的参数vec y\nobserved_data = [20,13,13,8]\nlambdas, likelihood_values = poisson_likelihood(observed_data)\nlikelihood_values /= np.sum(likelihood_values)\n\nplt.plot(lambdas, likelihood_values, color='black', lw=2)\nplt.xlabel('$\\lambda$')\nplt.ylabel('$L(\\lambda | Y=y)$')\nsns.despine()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 640x480 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/8593B03A2A20436981E5AD8108C72BDB/s2326hee52.png\">"},"metadata":{}}],"execution_count":6},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"9881119C9F2248488FB9EF7AC6D19F0D","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"###  Potential priors"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"A76CD1F5723B40A5A44788B3C2665F6A","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**🤔思考：**  \n\n在之前的例子中，我们主要关注泊松分布相关的似然函数，而忽略了参数$\\lambda$的先验分布。  \n$$  \nL(\\lambda | \\vec{y}) = \\frac{\\lambda^{\\sum y_i}e^{-n\\lambda}}{\\prod_{i=1}^n y_i!} \\propto \\lambda^{11}e^{-4\\lambda}  \n$$  \n\n\n请根据公式的形态猜测一下，以下哪种先验分布可以作为poisson似然的共轭先验？？  \n\n1. Gamma模型：$f(\\lambda) \\propto \\lambda^{s - 1} e^{-r \\lambda}$  \n\n2. Weibull模型：$f(\\lambda) \\propto \\lambda^{s - 1} e^{(-r \\lambda)^s}$  \n\n3. \"F\"模型：$f(\\lambda) \\propto \\lambda^{\\frac{s}{2} - 1}\\left( 1 + \\lambda\\right)^{-s}$  \n\n"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"85CC3737FC894EC780887552897962E5","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"### Gamma prior"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"26F55FBE4B9F4B94A7CCDAD1756486BA","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"考虑到 $\\lambda$ 是一个可以取任何正值的连续变量。  \n\n我们可以使用Gamma分布对它进行模拟，以此作为$\\lambda$的先验分布。  \n\n$$  \n\\lambda \\sim \\text{Gamma}(s, r)  \n$$  \n\nGamma分布由两个参数指定：  \n- 形状参数（shape parameter）或称为尺度形状参数（scale shape parameter）, 一般用 s 或者 $\\alpha$ 表示。它是一个正实数，控制着Gamma分布的形状。较大的s值使得分布更加陡峭，而较小的s值使得分布更加平坦。  \n- 尺度参数（scale parameter）或称为逆尺度参数（inverse scale parameter），一般用 r 或者 $\\beta$ 表示。它也是一个正实数，控制着Gamma分布的尺度。较大的r值使得分布的尺度变大，而较小的r值使得分布的尺度变小。  \n\n\n![Image Name](https://cdn.kesci.com/upload/s233qulcqf.png?imageView2/0/w/640/h/640)  \n> source: https://www.mbbpxt.co/gamma%E5%88%86%E5%B8%83/  \n\n其pdf为：  \n\n$$  \nf(\\lambda) = \\frac{r^s}{\\Gamma(s)} \\lambda^{s-1} e^{-r\\lambda} \\;\\; \\text{ for } \\lambda > 0  \n$$"},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"skip"},"id":"79192B773A5A4387B0AA9D6953C7195C","notebookId":"65251b1ab0f4556a3142003a","trusted":true},"source":"# @title interactive gamma\nimport ipywidgets\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport scipy.stats as st\nfrom scipy.stats import gamma\nfrom ipywidgets import interact\nimport seaborn as sns\n\n\nx = np.linspace(0, 20, 1000)\ndef plot_beta(a, b):\n    y = gamma.pdf(x, a=a, scale=1/b)\n    plt.plot(x, y, color='black', lw=2)\n    plt.title(f\"Gamma({a}, {b})\")\n    plt.xlabel('$\\lambda$')\n    plt.ylabel('$f(\\lambda)$')\n    plt.gca().yaxis.set_major_formatter('{x:1.1f}')\n    plt.xlim(0, 8)\n    plt.xticks(np.arange(0, 9, 2))\n    plt.yticks(np.arange(0, 2.5, 0.5))\n    sns.despine()","outputs":[],"execution_count":21},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"fragment"},"id":"80BC5C759F9F4F61A629239310D9D513","notebookId":"65251b1ab0f4556a3142003a","trusted":true},"source":"ipywidgets.interact(plot_beta, a=(0, 10), b=(0, 10))","outputs":[{"output_type":"execute_result","data":{"text/plain":"<function __main__.plot_beta(a, b)>"},"metadata":{},"execution_count":23},{"output_type":"display_data","data":{"text/plain":"<Figure size 640x480 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/CEE649D338CB4EDCA5DB83BCB1779A35/s2335zcem6.png\">"},"metadata":{}}],"execution_count":23},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"6A04DF616CBF4C868EB88BE2839482C3","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"已知$\\lambda$的平均值为12，我们选用$\\lambda \\sim \\text{Gamma}(24,2)$作为$\\lambda$的先验"},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"fragment"},"trusted":true,"id":"CD10F35202054F248B87BC981787C3E1","notebookId":"65251b1ab0f4556a3142003a"},"source":"x = np.arange(0, 15, 0.1)\nprior_y = st.gamma.pdf(x, a=10, scale=1/2)\nplt.plot(x, prior_y, color='black', lw=2)\nplt.title(f\"Gamma(10,2)\")\nplt.xlabel('$\\lambda$')\nplt.ylabel('$f(\\lambda)$')\nplt.gca().yaxis.set_major_formatter('{x:1.1f}')\nplt.xticks(np.arange(0, 16, 5))\nplt.yticks(np.arange(0, 0.3, 0.1))\nsns.despine()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 640x480 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/B5F73221BDD34C34906D231451626029/s2326io68d.png\">"},"metadata":{}}],"execution_count":10},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"D4CB988ACF564E5A94E98B357E4B86A3","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"### Gamma-Poisson conjugacy  \n\n由于Gamma先验分布和Poisson似然来自一组共轭家族，既二者组成的后验分布同样也是Gamma分布，具有计算上的便利性，因此我们可以直接代入公式来计算后验：  \n\n- Gamma-Poisson 家族的公式表达：  \n\n$$\\begin{align*}  \n先验：& \\lambda \\sim \\text{Gamma}(s, r) \\\\  \n似然：& Y_i | \\lambda \\stackrel{ind}{\\sim} \\text{Pois}(\\lambda) \\\\  \n后验：& \\lambda|\\vec{y} \\; \\sim \\; \\text{Gamma}\\left(s + \\sum y_i, \\; r + n\\right)  \n\\end{align*}  \n$$  \n"},{"cell_type":"markdown","metadata":{"id":"476CF5036BCA469E8F9D20D2F732A83F","notebookId":"65251b1ab0f4556a3142003a","runtime":{"status":"default","execution_status":null,"is_visible":false},"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"}},"source":"**代入当前数据可得：**  \n\n> $先验：\\lambda \\sim \\text{Gamma}(24, 2)$  \n\n> $\\vec{y} = (y_1,y_2,y_3,y_4) = (20,13,13,8); \\;\\;\\;\\; $  \n\n> $似然：L(\\lambda | \\vec{y}) = \\frac{\\lambda^{54}e^{-48\\lambda}}{20!\\times13!\\times13!\\times8!} \\propto \\lambda^{54}e^{-48\\lambda}$  \n\n> $s + \\sum y_i = 24 + 54; \\;\\;\\;\\; r + n = 2 + 4 \\\\  \n后验：\\lambda|\\vec{y} \\; \\sim \\; \\text{Gamma}(78,6)$  \n\n\n我们可以将这三者画出来：  \n\n在收集了四天的数据后，我们对$\\lambda$的信念发生了更新，"},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"9C407BB8A2E04C59B3470817B065ACC8","notebookId":"65251b1ab0f4556a3142003a","trusted":true},"source":"x = np.arange(0, 25, 0.1)\nprior_y = st.gamma.pdf(x, a=24, scale=1/2) / np.sum(st.gamma.pdf(x, a=24, scale=1/2))\nposterior = st.gamma.pdf(x, a=78, scale=1/6) / np.sum(st.gamma.pdf(x, a=78, scale=1/6))\n\nobserved_data = [20,13,13,8]\nlambdas, likelihood_values = poisson_likelihood(observed_data)\n\nlikelihood_values /= likelihood_values.sum()\n\nplt.plot(x, prior_y, color='black')\nplt.fill_between(x, prior_y, color=\"#f0e442\", alpha=0.5)\nplt.plot(lambdas, likelihood_values, color='black')\nplt.fill_between(lambdas, likelihood_values, color=\"#0071b2\", alpha=0.5)\nplt.plot(x, posterior, color='black')\nplt.xticks(np.arange(0, 16, 5))\nplt.fill_between(x, posterior, color=\"#009e74\", alpha=0.5)\n\n# 设置 x 和 y 轴标签\nplt.xlabel('$\\lambda$')\nplt.ylabel('density')\n\n# 创建自定义图例\ncustom_lines = [Line2D([], [], color=\"#f0e442\", lw=5),\n                Line2D([], [], color=\"#0071b2\", lw=5),\n                Line2D([], [], color=\"#009e74\", lw=5)]\n        \n# 将图例放置在子图外部的右上角\nplt.legend(custom_lines, ['prior', 'likelihood', 'posterior'], loc='upper right', bbox_to_anchor=(1, 1))\n\n# 移除图的上、右边框线\nsns.despine()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 640x480 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/270707D794574D52A8D2C199C8A13EA0/s2326ijr57.png\">"},"metadata":{}}],"execution_count":11},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"27D65062551E48A0A7A1A5814B31EFF4","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**扩展**  \n\n除了 Gamma分布外，Poisson分布的共轭先验还包括 Negative Binomial分布:  \n\n$$  \n\\begin{align*}  \nGamma: & f(\\lambda) \\propto \\lambda^{s - 1} e^{-r \\lambda} \\\\  \nNegative-binomial: & P(X = k) = (k + r - 1) C (k) * p^k * (1 - p)^r  \n\\end{align*}  \n$$  \n\n然而，Weibull模型和\"F\"模型不是Poisson分布的共轭先验：  \n\n$$  \n\\begin{align*}  \nWeibull模型：& f(\\lambda) \\propto \\lambda^{s - 1} e^{(-r \\lambda)^s} \\\\  \nF模型：& f(\\lambda) \\propto \\lambda^{\\frac{s}{2} - 1}\\left( 1 + \\lambda\\right)^{-s}  \n\\end{align*}  \n$$"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"ADF44E0833C0467FA8D3231C53A1215F","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"## Normal-Normal conjugate family  \n\n我们已经学习了两个共轭族：Beta-Binomial 和 Gamma-Poisson。  \n\n但还有更多的共轭族存在！最常见的就是：Normal-Normal。  \n\n**来看一个经典的例子：身高分布**  \n\n国家卫健委：2020年我国18-44岁成年人平均身高为165厘米。随着生成质量的提高，我们有理由相信，这一数据有所提高，所以我们打算新收集一万名被试(n=10000)的身高数据，并使用正态-正态贝叶斯模型来估计**根据数据更新后**的平均身高。  \n\n![](https://th.bing.com/th/id/OIP.fdUphNKKkEvn9msoc8VRFwHaE8?pid=ImgDet&rs=1)  \n\n> source: https://baijiahao.baidu.com/s?id=1736757271631341254&wfr=spider&for=pc"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"1C79B769B2114BFBB0CF68BE057B0A63","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"### The Normal data model  \n\n同样，我们的贝叶斯分析也是从对平均身高(μ)的先验开始。  \n然而，要为μ指定一个适当的先验模型结构，需要考虑数据(Yi)的特征。  \n\n由于身高 Yi 是连续数据，因此有很多潜在可能的模型：贝塔模型、指数模型、伽马模型、正态分布模型、F 模型等。  \n- 从这个列表中，我们可以立即排除贝塔模型。因为它假设 Yi ∈ [0，1]，而身高往往在 165cm 左右。  \n- 在剩下的选项中，正态模型是相当合理的。这也是我们直觉上的第一反应，毕竟身高通常围绕某个总体平均值（μ）呈现为对称或正态分布。"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"2C205A057BF44432AEA99BEDBB6A1EB0","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**正态模型介绍**  \n\n假设 Y 是一个连续随机变量，可以取 -∞ 和 ∞ 之间的任意值，即 Y∈ (-∞, ∞)。  \n\n那么，Y 的变异性可以很好地用均值参数 μ∈ (-∞, ∞) 和标准差参数 σ > 0 的正态模型来表示：  \n\n$$  \nY \\sim N(\\mu, \\sigma^2).  \n$$  \n\n正态模型的连续 pdf 形式:  \n\n$$  \n\\begin{equation}  \nf(y) = \\frac{1}{\\sqrt{2\\pi\\sigma^2}} \\exp\\bigg[{-\\frac{(y-\\mu)^2}{2\\sigma^2}}\\bigg] \\;\\; \\text{ for } y \\in (-\\infty,\\infty)  \n\\tag{5.11}  \n\\end{equation}  \n$$  \n\n正太模型的基本特征：  \n\n$$  \n\\begin{split}  \nE(Y) & = \\text{ Mode}(Y) = \\mu \\\\  \n\\text{Var}(Y) & = \\sigma^2 \\\\  \n\\text{SD}(Y) & = \\sigma. \\\\  \n\\end{split}  \n$$  \n\n> 此外，大家都知道的，大约 95% 的 Y 值都在 μ 的 2 个标准差范围内。"},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"DACFED0C13884EA39316F073DFEA3ABD","notebookId":"65251b1ab0f4556a3142003a","trusted":true},"source":"import seaborn as sns\nimport matplotlib.pyplot as plt\n\n# 创建数据\nimport numpy as np\nnp.random.seed(0)\ndata1 = np.random.normal(loc=0, scale=1, size=1000)\ndata2 = np.random.normal(loc=5, scale=2, size=1000)\ndata3 = np.random.normal(loc=10, scale=5, size=1000)\n\n# 绘制图形\nfig, axes = plt.subplots(1, 3, figsize=(15, 5))\n\n# 绘制第一个图\nsns.histplot(data=data1, kde=True, ax=axes[0], stat = 'density')\naxes[0].set_title('Mean=0, Std=1')\naxes[0].set_ylim([0, 0.5])\n\n# 绘制第二个图\nsns.histplot(data=data2, kde=True, ax=axes[1], stat = 'density')\naxes[1].set_title('Mean=5, Std=2')\naxes[1].set_ylim([0, 0.5])\n\n# 绘制第三个图\nsns.histplot(data=data3, kde=True, ax=axes[2], stat = 'density')\naxes[2].set_title('Mean=10, Std=5')\naxes[2].set_ylim([0, 0.5])\n\n# 显示图形\nplt.tight_layout()\n# 移除图的上、右边框线\nsns.despine()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 1500x500 with 3 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/F842A10623C94892A79EBCC2254CB572/s2326jew71.png\">"},"metadata":{}}],"execution_count":12},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"685E88CE7DF147AC8B688BE5E3A5E711","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"上图展示了各种平均值和标准差参数值 μ 和 σ 下的正态模型。  \n\n- 无论参数如何变化，正态模型都呈钟形，并围绕 μ 对称  \n- 因此，当 μ 变大时，模型也会随之向右移动。  \n- 此外，σ 控制着正态模型的变异性--σ 越大，模型越分散。  \n- 最后，虽然正态变量 Y 的取值范围可以从-∞到∞，但如果 Y 值与均值 μ 的差值超过 3 个标准差 σ，则正态模型赋予 Y 值的可信度可以忽略不计。"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"7B4D10BDE5FF42B1BC819C6D89FC3CD2","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**Normal data model**  \n\n回到我们的例子，我们假设收集了 n = 10000 名被试的身高（Y1, Y2, ... , Yn）  \n- 这些被试的身高显然是彼此独立的，  \n- 并且围绕平均身高 µ 的正态分布，标准偏差为 σ。  \n- 此外为了保持对 µ 的关注，我们将在整个分析过程中的标准偏差固定为 σ = 5 cm。3 σ 的范围表明，大多数人的身高在 150 cm 到 180 cm 之间。  \n- 因此，10000 名被试的身高的均值为 168cm，标准差为 5cm。 显然，相比于 2020 年的数据 165cm，现在的数据 168cm 体现了大众身高的增长。  \n\n接着我们可以定义出正态模型的似然函数：  \n\n$$  \nL(\\mu |\\vec{y}) \\propto \\prod_{i=1}^{n} \\exp\\bigg[{-\\frac{(y_i-\\mu)^2}{2\\sigma^2}}\\bigg] =  \\exp\\bigg[{-\\frac{\\sum_{i=1}^n(y_i-\\mu)^2}{2\\sigma^2}}\\bigg]  .  \n$$  \n\n"},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"fragment"},"trusted":true,"id":"9989C239AE6844D787DF91469D5D5BEC","notebookId":"65251b1ab0f4556a3142003a"},"source":"# 模拟10000名被试身高的数据并绘图\n\n# 创建数据\nnp.random.seed(0)\ndata = np.random.normal(loc=168, scale=5, size=10000)\n\n# 绘图\nsns.histplot(data=data, kde=True, stat = 'density')\n\nplt.xlabel('Height (cm) for 10000 participants')\n# 显示图形\nplt.tight_layout()\n# 移除图的上、右边框线\nsns.despine()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 640x480 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/B4D57DB263DD4A039F6E101C3862D897/s2326kwjfd.png\">"},"metadata":{}}],"execution_count":13},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"F8F28D92354346008848E0EB4ADB11E6","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"### Normal prior  \n\n根据2020年之前的数据，我们可以用 165cm 作为我们正态先验模型的均值，即 µ = 165。  \n\n具体来说，我们假设 μ 本身是围绕某个均值 θ 和标准差 τ 的正态分布 (其中， θ 和 τ是超参)：  \n\n$$  \n\\begin{align*}  \n\\mu & \\sim N(\\theta, \\tau^2) \\\\  \nf(\\mu) & = \\frac{1}{\\sqrt{2\\pi\\tau^2}} \\exp\\bigg[{-\\frac{(\\mu - \\theta)^2}{2\\tau^2}}\\bigg] \\;\\; \\text{ for } \\mu \\in (-\\infty,\\infty)  \n\\tag{5.14}  \n\\end{align*}  \n$$  \n\n注意👁️，正态先验假设 μ∈（-∞，∞）与我们对数据模型的假设一致，我们还将在后面证明这是一个共轭先验。因为两个模型都正比于该项：  \n\n$$  \n\\exp\\bigg[{-\\frac{(\\mu - \\blacksquare)^2}{2\\blacksquare^2}}\\bigg]  \n$$  \n\n> 先验模型和数据模型的差异仅表现为均值的差异，即先验模型的均值为 165cm，而数据模型的均值为 168cm。"},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"fragment"},"trusted":true,"id":"832B32A8DCA042A6A0A78BA552079BB4","notebookId":"65251b1ab0f4556a3142003a"},"source":"# 绘制先验\n\n# 创建数据\nnp.random.seed(0)\ndata = np.random.normal(loc=165, scale=5, size=1000)\n\n# 绘图\nsns.histplot(data=data, kde=True, stat = 'density')\n\nplt.xlabel('Height (cm) of prior')\n# 显示图形\nplt.tight_layout()\n# 移除图的上、右边框线\nsns.despine()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 640x480 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/570F27F402E949DBA2EA527D5409BED0/s2326kzk4w.png\">"},"metadata":{}}],"execution_count":14},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"91CFE79150FD49F4800F00CB7962C1CD","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"### Normal-Normal conjugacy  \n\n为了得到 μ 的后验模型，我们必须将先验模型和数据信息结合起来。  \n\n> 需要强调的是，我们选择的正态先验模型属于正态-正态的共轭族！因此，μ 的后验模型也将是正态的，并根据先验数据和观测数据更新参数。  \n\n以下公式展示了μ的后验模型：  \n\n$$  \n\\begin{align*}  \nY_i | \\mu & \\stackrel{ind}{\\sim} N(\\mu, \\sigma^2) \\\\  \n\\mu & \\sim N(\\theta, \\tau^2) \\\\  \n\\mu|\\vec{y} \\; & \\sim \\;  N\\bigg(\\theta\\frac{\\sigma^2}{n\\tau^2+\\sigma^2} + \\bar{y}\\frac{n\\tau^2}{n\\tau^2+\\sigma^2}, \\; \\frac{\\tau^2\\sigma^2}{n\\tau^2+\\sigma^2}\\bigg)  \n\\tag{5.15}  \n\\end{align*}  \n$$  \n\n首先，后验均值是先验均值 E(μ) = θ 和样本均值 y 的加权平均值；  \n\n$$  \n\\text{posterior mean} \\sim \\theta\\frac{\\sigma^2}{n\\tau^2+\\sigma^2} + \\bar{y}\\frac{n\\tau^2}{n\\tau^2+\\sigma^2}  \n$$  \n\n其次，后验方差受先验变异性 τ 和数据变异性 σ 的影响：  \n\n$$  \n\\text{posterior variance}~\\frac{\\tau^2\\sigma^2}{n\\tau^2+\\sigma^2}  \n$$  \n\n两者都受到样本量n的影响。  \n- 首先，随着n的增加，后验均值对先验均值的权重减少，对样本均值y的权重增加: $\\frac{\\sigma^2}{n\\tau^2+\\sigma^2} \\to 0 \\;\\; \\text{ and } \\;\\;\\frac{n\\tau^2}{n\\tau^2+\\sigma^2} \\to 1 .$  \n- 此外，随着 n 的增大，后验方差会减小 $\\frac{\\tau^2\\sigma^2}{n\\tau^2+\\sigma^2} \\to 0$  \n- 也就是说，我们掌握的数据越多，我们对 μ 的后验确定性就越高，也就越与数据的结果相似。"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"1C261AE222894E79AC078DD66BB96E13","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"结合之前的先验模型和数据模型，现在我们可以来计算后验并绘制它。  \n\n- 先验模型假设全国成人的平均身高为 165cm，标准偏差为 5cm。  \n- 数据模型为 10000 名被试的身高数据，平均身高为 168cm，标准偏差为 5cm。"},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"3FBE3B9404E1448AA357D44CC1938488","notebookId":"65251b1ab0f4556a3142003a","trusted":true},"source":"# 定义身高范围\nx = np.linspace(150, 200, 10000)\n\n# 定义先验分布\nprior_mean = 165\nprior_std = 5\nprior_y = st.norm.pdf(x, loc = prior_mean, scale = prior_std) / np.sum(st.norm.pdf(x, prior_mean, prior_std))\n\n# 生成似然\nlikelihood_mean = 168\nlikelihood_std = 5\nlikelihood_values = st.norm.pdf(x, loc = likelihood_mean, scale = likelihood_std) / np.sum(st.norm.pdf(x, likelihood_mean, likelihood_std))\n\n# 计算后验分布\nposterior_mean = (prior_mean * prior_std**2 + likelihood_mean * likelihood_std**2) / (prior_std**2 + likelihood_std**2)\nposterior_std = np.sqrt((prior_std**2 * likelihood_std**2) / (prior_std**2 + likelihood_std**2))\nposterior = st.norm.pdf(x, loc = posterior_mean, scale = posterior_std) / np.sum(st.norm.pdf(x, posterior_mean, posterior_std))\n\n\nplt.plot(x, prior_y, color='#f0e442', label=\"prior\")\nplt.fill_between(x, prior_y, color=\"#f0e442\", alpha=0.5)\nplt.plot(x, likelihood_values, color='#0071b2', label=\"likelihood\")\nplt.fill_between(x, likelihood_values, color=\"#0071b2\", alpha=0.5)\nplt.plot(x, posterior, color='#009e74', label=\"posterior\")\nplt.fill_between(x, posterior, color=\"#009e74\", alpha=0.5)\n\n# 设置 x 和 y 轴标签\nplt.xlabel('$\\mu$ for height')\nplt.ylabel('density')\nplt.legend()\n\n# 移除图的上、右边框线\nsns.despine()","outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 640x480 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/61DB0475CDD549A2B738DD8AD4518234/s233jeqkw0.png\">"},"metadata":{}}],"execution_count":25},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"BD4C13C0503C42DC8C6E5AE75CD6DF42","notebookId":"65251b1ab0f4556a3142003a","trusted":true},"source":"print(f\"后验分布的均值：{posterior_mean}。 介于先验和似然之间。\")","outputs":[{"output_type":"stream","name":"stdout","text":"后验分布的均值：166.5。 介于先验和似然之间。\n"}],"execution_count":26},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"263381EC101D4ECA828C521F6CAE96A7","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**扩展阅读：证明 Normal-Normal conjugacy**  \n\n为完整起见，我们在此证明正态-正态模型产生后验模型。如果你不喜欢推导，也没关系。请随意跳到下一节。  \n\n让我们直接进入正题。μ 的后验 pdf 与正态先验 pdf 和似然函数 的乘积成正比。对于 μ∈ (-∞, ∞)：  \n$$  \nf(\\mu|\\vec{y}) \\propto f(\\mu)L(\\mu|\\vec{y}) \\propto \\exp\\bigg[{\\frac{-(\\mu - \\theta)^2}{2\\tau^2}}\\bigg] \\cdot \\exp\\bigg[{-\\frac{(\\bar{y}-\\mu)^2}{2\\sigma^2/n}}\\bigg]  .  \n$$  \n\n接下来，我们可以展开指数中的平方，包括第一个指数分子中的$\\theta^2$ 和第二个指数分子中的 $y^2$：  \n$$  \n\\begin{split}  \nf(\\mu|\\vec{y})  \n& \\propto  \\exp\\Bigg[{\\frac{-\\mu^2+2\\mu\\theta-\\theta^2}{2\\tau^2}}\\Bigg]\\exp\\Bigg[{\\frac{-\\mu^2+2\\mu\\bar{y}-\\bar{y}^2}{2\\sigma^2/n}}\\Bigg] \\\\  \n& \\propto  \\exp\\Bigg[{\\frac{-\\mu^2+2\\mu\\theta}{2\\tau^2}}\\Bigg]\\exp\\Bigg[{\\frac{-\\mu^2+2\\mu\\bar{y}}{2\\sigma^2/n}}\\Bigg]. \\\\  \n\\end{split}  \n$$  \n\n我们给出指数的共同分母，并将它们合并成一个指数：  \n$$  \n\\begin{split}  \nf(\\mu|\\vec{y})  \n& \\propto  \\exp\\Bigg[{\\frac{(-\\mu^2+2\\mu\\theta)\\sigma^2/n}{2\\tau^2\\sigma^2/n}}\\Bigg]\\exp\\Bigg[{\\frac{(-\\mu^2+2\\mu\\bar{y})\\tau^2}{2\\tau^2\\sigma^2/n}}\\Bigg] \\\\  \n& \\propto  \\exp\\Bigg[{\\frac{(-\\mu^2+2\\mu\\theta)\\sigma^2 +(-\\mu^2+2\\mu\\bar{y})n\\tau^2}{2\\tau^2\\sigma^2}}\\Bigg]. \\\\  \n\\end{split}  \n$$  \n\n现在，让我们把 μ 项合并起来，重新排列，这样 $\\mu^2$ 就是：  \n$$  \n\\begin{split}  \nf(\\mu|\\vec{y})  \n& \\propto  \\exp\\Bigg[{\\frac{-\\mu^2(n\\tau^2+\\sigma^2)+2\\mu(\\theta\\sigma^2+ \\bar{y}n\\tau^2) }{2\\tau^2\\sigma^2}}\\Bigg] \\\\  \n& \\propto  \\exp\\Bigg[{\\frac{-\\mu^2+2\\mu\\left(\\frac{\\theta\\sigma^2 + \\bar{y}n\\tau^2}{n\\tau^2+\\sigma^2}\\right) }{2(\\tau^2\\sigma^2) /(n\\tau^2+\\sigma^2)}}\\Bigg]. \\\\  \n\\end{split}  \n$$  \n\n仔细观察，你会发现我们可以带回一些不依赖于 μ 的常数来完成分子中的平方：  \n$$  \n\\begin{split}  \nf(\\mu|\\vec{y})  \n& \\propto  \\exp\\Bigg[{\\frac{-\\bigg(\\mu - \\frac{\\theta\\sigma^2 + \\bar{y}n\\tau^2}{n\\tau^2+\\sigma^2}\\bigg)^2 }{2(\\tau^2\\sigma^2) /(n\\tau^2+\\sigma^2)}}\\Bigg]. \\\\  \n\\end{split}  \n$$  \n\n这看起来仍然很混乱，但是一旦我们完成了正方形，我们实际上就得到了 μ 的正态 pdf 内核，即 exp部分。通过确定缺失的部分 ∎，我们可以得出以下结论  \n$$  \n\\mu|\\vec{y} \\;  \\sim  \\; N\\left(\\frac{\\theta\\sigma^2+ \\bar{y}n\\tau^2}{n\\tau^2+\\sigma^2}, \\;{\\frac{\\tau^2\\sigma^2}{n\\tau^2+\\sigma^2}} \\right)  \n$$  \n\n我们可以将后验均值重组为先验均值 μ 和数据均值 y 的加权平均值：  \n$$  \n\\frac{\\theta\\sigma^2+ \\bar{y}n\\tau^2}{n\\tau^2+\\sigma^2} = \\theta\\frac{\\sigma^2}{n\\tau^2+\\sigma^2} + \\bar{y}\\frac{n\\tau^2}{n\\tau^2+\\sigma^2} .  \n$$  \n\n"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"53D9D311257748CE862CCB95889B56D6","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**Critiques of conjugate family models**  \n\n最后，我们补充一下共轭族模型的**缺点**  \n\n1. 共轭先验模型可能并不能适应你的先验信念。  \n  - 例如，正态模型总是围绕均值 μ 对称。因此，如果你的先验认识不是对称的，那么正态先验模型可能就不适合作为先验。  \n2. 共轭族模型并不允许有一个完全平坦的先验(uniform prior)。  \n  - 虽然我们可以通过设置 α = β = 1 来使得 Beta 先验变得更加平坦，但无论是 Normal 还是 Gamma 先验（或任何具有无限支持的适当模型）都**无法调整为完全平坦**。  \n- 我们能做的最好的办法就是调整先验，使其具有非常高的方差，这样它们就几乎是平的了。"},{"cell_type":"markdown","metadata":{"jupyter":{},"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"92FA9F4DF3D44B819F47CB20648AAD23","runtime":{"status":"default","execution_status":null,"is_visible":false},"notebookId":"65251b1ab0f4556a3142003a"},"source":"**练习**  \n\n在先前的示例中，我们将身高的标准差固定为5，这显然不符实实时。  \n\n在接下来的练习中，你可以尝试将标准差设置为不同的值，并观察结果的变化。  \n\n你练习的目标：  \n- 如何调整似然，使得后验均值接近更先验的均值？  \n- 如何调整似然，使得后验均值更接近后验的均值？"},{"cell_type":"code","metadata":{"jupyter":{},"collapsed":false,"scrolled":false,"tags":[],"slideshow":{"slide_type":"slide"},"id":"42D89B281B0A45249D050C0BC2E3A751","notebookId":"65251b1ab0f4556a3142003a","trusted":true},"source":"# 定义身高范围\nx = np.linspace(150, 200, 10000)\n\n# 定义先验分布\nprior_mean = 165\nprior_std = 5\n\n# 生成似然\n#===========================================================================\n#                            请调整似然，观察不同似然参数对于先验的影响。\n#                            请修改 ... 中的值。\n#===========================================================================\nlikelihood_mean = ...\nlikelihood_std = ...\nlikelihood_values = st.norm.pdf(x, loc = likelihood_mean, scale = likelihood_std) / np.sum(st.norm.pdf(x, likelihood_mean, likelihood_std))\n\n#===========================================================================\n#                            进阶：你可以尝试修改先验为不同分布，体验非共轭先验带来的影响 \n#                            例如，将 norm 修改为 gamma，并修改对应的参数。(这里不做任何修改也可以运行)\n#===========================================================================\nprior_y = st.norm.pdf(x, loc = prior_mean, scale = prior_std) / np.sum(st.norm.pdf(x, prior_mean, prior_std))\n\n# 计算后验分布\nposterior_mean = (prior_mean * prior_std**2 + likelihood_mean * likelihood_std**2) / (prior_std**2 + likelihood_std**2)\nposterior_std = np.sqrt((prior_std**2 * likelihood_std**2) / (prior_std**2 + likelihood_std**2))\nposterior = st.norm.pdf(x, loc = posterior_mean, scale = posterior_std) / np.sum(st.norm.pdf(x, posterior_mean, posterior_std))\n\n\nplt.plot(x, prior_y, color='#f0e442', label=\"prior\")\nplt.fill_between(x, prior_y, color=\"#f0e442\", alpha=0.5)\nplt.plot(x, likelihood_values, color='#0071b2', label=\"likelihood\")\nplt.fill_between(x, likelihood_values, color=\"#0071b2\", alpha=0.5)\nplt.plot(x, posterior, color='#009e74', label=\"posterior\")\nplt.fill_between(x, posterior, color=\"#009e74\", alpha=0.5)\n\n# 设置 x 和 y 轴标签\nplt.xlabel('$\\mu$ for height')\nplt.ylabel('density')\nplt.legend()\n\n# 移除图的上、右边框线\nsns.despine()","outputs":[],"execution_count":null}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python","nbconvert_exporter":"python","file_extension":".py","version":"3.5.2","pygments_lexer":"ipython3"}},"nbformat":4,"nbformat_minor":2}

一键复制编辑原始数据按行查看历史

提交于 2024-08-25 14:20 +08:00 . move 2023 folder

马建仓 AI 助手

尝试更多

代码解读

代码找茬

代码优化

https://gitee.com/hcp4715/PyBayesian.git

git@gitee.com:hcp4715/PyBayesian.git

hcp4715

PyBayesian

master

hcp4715/PyBayesian

简介

发行版 (2)

贡献者

语言

近期动态

hcp4715/PyBayesian .gitee-modal { width: 500px !important; }

简介

发行版 (2)

贡献者

语言

近期动态

搜索帮助

hcp4715/PyBayesian