{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 2.1 二元变量" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 伯努利分布" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "考虑二元随机变量 $x\\in \\{0,1\\}$（抛硬币，正面为 1，反面为 0），其概率分布由参数 $\\mu$ 决定：\n", "\n", "$$\n", "p(x=1)=\\mu\n", "$$\n", "\n", "其中 $(0 \\leq\\mu \\leq 1)$，并且有 $p(x=0)=1-\\mu$。这就是伯努利分布（`Bernoulli distribution`），其概率分布可以写成：\n", "\n", "$$\n", "\\text{Bern}(x|\\mu) = \\mu^{x}(1-\\mu)^{1-x}\n", "$$\n", "\n", "均值和方差为：\n", "\n", "$$\n", "\\begin{align}\n", "\\mathbb E[x] & = \\mu \\\\\n", "\\text{var}[x] & = \\mu(1-\\mu)\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "### 伯努利分布的最大似然估计" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "考虑一组 $x$ 的观测数据 $\\mathcal D = \\{x_1, \\dots, x_N\\}$，在独立同分布的假设下，其似然函数为\n", "\n", "$$\n", "p(\\mathcal D|\\mu) = \\prod_{n=1}^N p(x_n|\\mu) = \\prod_{n=1}^N \\mu^{x_n}(1-\\mu)^{1-x_n}\n", "$$\n", "\n", "对数似然为\n", "\n", "$$\n", "\\ln p(\\mathcal D|\\mu) = \\sum_{n=1}^N \\ln p(x_n|\\mu) \n", "= \\sum_{n=1}^N \\left\\{x_n\\ln \\mu + (1-x_n)\\ln (1-\\mu)\\right\\}\n", "$$\n", "\n", "对数似然值只依赖于 $\\sum_{n=1}^N x_n$ 的取值，而事实上 $\\sum_{n=1}^N x_i$ 就是伯努利分布的一个充分统计量，它可以提供参数 $\\mu$ 的全部信息。\n", "\n", "对 $\\mu$ 最大化对数似然，我们很容易得到\n", "\n", "$$\n", "\\mu_{ML} = \\frac 1 N \\sum_{n=1}^N {x_n}\n", "$$\n", "\n", "即最大似然估计值为样本均值（`sample mean`），若样本中 $x=1$ 的数目为 $m$ 则：\n", "\n", "$$\n", "\\mu_{ML} = \\frac{m}{N}\n", "$$\n", "\n", "考虑抛三次硬币出现了三次正面的情况，此时 $N=m=3, \\mu_{ML} = 1$。在这种情况下，最大似然估计会得到每次都是正面的结果，这显然违背了我们的正常认知。事实上，这是一种过拟合的典型表现。\n", "\n", "为了解决这个问题，我们可以考虑引入先验知识。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 二项分布" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "给定数据总数 $N$，$x=1$ 的总次数 $m$ 满足一定的分布，这个分布叫做二项分布（`binomial distribution`）。\n", "\n", "从伯努利分布的似然函数中可以看出它应该正比于 $\\mu^{m}(1-\\mu)^{N-m}$，事实上它可以写成：\n", "\n", "$$\n", "\\text{Bin}(m~|~N,\\mu) = \\begin{pmatrix} N \\\\m \\end{pmatrix} \\mu^{m}(1-\\mu)^{N-m}\n", "$$\n", "\n", "其中\n", "\n", "$$\n", "\\begin{pmatrix} N \\\\m \\end{pmatrix} \\equiv \\frac{N!}{(N-m)!m!}\n", "$$\n", "\n", "是组合数。\n", "\n", "验证它是一个概率分布，二项式定理给出：\n", "\n", "$$\n", "\\sum_{m=0}^N \\begin{pmatrix} N \\\\m \\end{pmatrix} \\mu^{m}(1-\\mu)^{N-m} = (\\mu + 1 - \\mu)^N = 1\n", "$$\n", "\n", "下图是 $N = 10, \\mu=0.25$ 的分布的一个直方图示意。" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAAEgCAYAAAB7MvKtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE01JREFUeJzt3X+07XVd5/Hni3tBQPyBWqBwjVQULEnUgFFv3hTzaina\n1OCvsHSZM0ljK0eJGUauuVyBLrNmWCklFv6CyknCMRVNT2kJAQpIXkxQpgsp/koNy7pX3vPH93to\ne7jn7PNj7+8+n3uej7XOuvv73Z/zfb/3uWe99vd89vdHqgpJUlv2m3UDkqSVM7wlqUGGtyQ1yPCW\npAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4a0mJTkoyZuTbF2wfkuS85O8NMlbkhw1mw5X1kuSByc5J8lv\nJ/lQkkcveP4VSV6e5CFJnpTkD6fdv9a3zbNuQFqpJC8BHgT8DPCuBU+/HTirqj6e5CTgIuA/DNzi\ninpJsh/wcuClVVVJfg64LMnDquqr/bC7A68CXg/8PfAfB3kFWrfc89ZEJblbkv+Z5NYkO5Nk5Lnn\nJrkhyXuTPGu1Narq/Ko6A/inBbUfCvxIVX28H3c58LAkP7DaWqu1wl6OBh4LPKBfvogurJ85MqaA\nHwJOBB5cVVdPq3e1wfDWRFXVvwKvo9vrPIKRAKqqdwEXVtXTq+o9Uyh/HLBrwbpb6IJxaCvp5Xbg\nSOD+AFW1p1936OigqtpZVVdW1Xcn365a47SJpmEr8KfAAcArgNGgvmOKdb8f+JcF6/4ZOGyKNdfc\nS1XdCtxvfjnJFuA+wCdGxyV5MRDgeOAPquqKCfeshhjemoYfBc6h29O8KcnWqvpYkmOAz44OTHJv\n4DfpQmkpb6yq68aMORT4twXr/g24x1LflOQ04GFV9T/65RcDh1fVa9bQ56p66b0UeN/8lEvv/wLX\nV9W/9uF+RZKHVtXty9ie9kGGt6ZhU3XXGt6V5CLglcDHgB8H3jk6sKq+AbxwQnX/aS/rDga+Pub7\nfhp428jyzwLnjQ5YRZ+r6iXJDwNPAZ60oP7VI493JdkDPIO7fmCrDcI5b01UksOAL42sej3w1CQP\nBw6tqm9NsfwXgUMWrLs7cNti35BkE/B44M/75QOBk4CPzqCXe9D9xfK0qvr6yPq7Jzmr723UkWvs\nUQ1zz1uT9iTgsvmFqro+yQeAs4C7HCGR5FDgDUxm2uRyRgItyWbgB4Drl/ieRwOfr6pv9stb+/Hf\nTnKf+RBdRZ+r6eUc4PSq+mJ/lM6zq+oi4Bi6v14upPtrZhPdHPlNY3rRPszw1qQ9qD+qZNTrgDng\nfy8cXFX/yNqmTe4M03464QtJTuw/zHsy8MmqugEgyWOAe1bVR0a+/2Tg8yPLzwX+Engi8A/00xwr\n7XOlvSR5Jd1nBMf0nw08kG7vHeBaujeF+aNXngj8P+C9y+1H+x7DWxOR5FHAfwZOSbJfVf36/HNV\n9ZdJ3g9M5OiIJC8AngpsAd6Y5GPAf6uq3cBpwFn9STGPBp438q3PpQu+R46sexLw3SQ/D+wPfKRf\n99Cq+vAaWx3Xy/cBH0lyLPBaYNPI8/PHdVNVe5K8J8nrgd10R6w8uaoWfiCqDSTj7mGZZDvwW3S/\nWG+pqnMXPH8K8Ot0h4DdAbxiwZ6NtG4keVFVXdA/PhC4le7Ikt2z7UxamSX3vPu5tfPo/rS8Fbgy\nyaVVtXNk2Ier6k/78Y+gO6b3IVPqV1qre408fjxwjcGtFo072uQE4Maqurn/Bb8YOGV0QFV9e2Tx\nEOCrSOtQkp+kP6qk9wjgkhm1I63JuDnvI/jeU3xvobu2wvdI8kzgN+hO7/2JiXUnTVBVvW/B8htn\n1Yu0VuP2vJeeEJ8fVHVJVR0LPJ3umhaSpCkat+d9K90n+vO20O1971V/CvTmJPetqq+NPpdkWW8E\nkqTvVVV3Ob9g3J73VcDRSY5KcgBwKnDp6ID+IvLpHz+qL/S1u2ypW7+qr7PPPnvV37uWr1nV9TVv\njLq+Zusu52sxS+55V3d86enAB+kOFbygqnb2F8Onqs6nuyj8aUl2013G8tlj3hAkSWs09iSdqno/\n8P4F684fefw6ujPoJEkDaeLCVNu2bdtQdWdZ29e8MWr7mtuvO/YMy4kVSmqoWpK0r0hCreIDS0nS\nOmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1yPCWpAZ5G7SB9Jd/mSqPo5c2DsN7UNMM1+m/OUhaP5w2\nkaQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQgw1uSGmR4S1KDDG9J\napDhLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSg7wB8T5uiLvWg3eul4ZmeG8I0w5W71wv\nDc1pE0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQgw1uSGmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1\nyPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1JDTK8JalBhrckNcjwlqQGGd6S1CDDW5IaZHhLUoMM\nb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQgw1uSGmR4S1KDDG9JatDY8E6yPckNST6X5Iy9PP+8\nJNcmuS7JXyU5bjqtSpLmLRneSTYB5wHbgYcDz0ly7IJhnwd+rKqOA14D/O40GpUk/btxe94nADdW\n1c1VtRu4GDhldEBVfaKqvtkvXgEcOfk2JUmjxoX3EcCukeVb+nWLeRHwZ2ttSpK0tM1jnq/lbijJ\njwMvBB632JgdO3bc+Xjbtm1s27ZtuZuXpA1hbm6Oubm5seNStXg+JzkJ2FFV2/vlM4E7qurcBeOO\nA/4E2F5VNy6yrVqq1r4uCSt4L1xNBfb2851+3cVrS1q7JFRVFq4fN21yFXB0kqOSHACcCly6YMMP\npAvu5y8W3JKkyVpy2qSq9iQ5HfggsAm4oKp2JnlJ//z5wKuAQ4E3dXt57K6qE6bbtiRtbEtOm0y0\nkNMmOG0iaaVWO20iSVqHDG9JapDhLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ\n4S1JDTK8JalBhrckNcjwlqQGGd6S1KBx97Dcp/Q3i5g6r20tado2VHh3pn9jAkmaNqdNJKlBhrck\nNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQgw1uSGmR4S1KD\nNuAlYTWUIa6f7rXTtVEZ3pqyaYar107XxuW0iSQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1J\nDTK8JalBhrckNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQg\nw1uSGmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1JDTK8\nJalBY8M7yfYkNyT5XJIz9vL8MUk+keQ7SV4+nTYlSaM2L/Vkkk3AecDJwK3AlUkuraqdI8O+Bvwy\n8MypdSlJ+h7j9rxPAG6sqpurajdwMXDK6ICq+kpVXQXsnlKPkqQFxoX3EcCukeVb+nWSpBlactoE\nqEkW27Fjx52Pt23bxrZt2ya5eUlq3tzcHHNzc2PHpWrxfE5yErCjqrb3y2cCd1TVuXsZezZwe1W9\nYZFt1VK1hpCECb8f7a0Ke3ud0689q7qzrL33utK+JAlVlYXrx02bXAUcneSoJAcApwKXLlZjjT1K\nkpZpyWmTqtqT5HTgg8Am4IKq2pnkJf3z5yc5HLgSuCdwR5KXAQ+vqtun3LskbVhLTptMtJDTJlOu\n7bSJtC9a7bSJJGkdMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1JDTK8JalBhrck\nNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQgw1uSGmR4S1KD\nNs+6AWnSkky9RlVNvYa0FMNb+6hphuv03xykcZw2kaQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0a\n/FDB8857E5dd9vGp1jjhhOM466wzplpDkmYpQ51skKSqimc96zQuueRgYOuUKl3LiSdex+WXf2Bv\nPTDd438BstcTOKZfe1Z1Z1l7fdWVpiEJVXWXkwtmdJLOY4HnTWnb9wWum9K2JWl9cM5bkhpkeEtS\ngwxvSWqQ4S1JDTK8JalBhrckNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI\n8JakBhnektQgw1uSGmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1aPOsG5D2FUkGqVNVg9TR+mZ4SxM1\n7WAd5g1C65/TJpLUIMNbkhpkeEtSgwxvSWqQ4S1JDTK8JalBhrckNWhseCfZnuSGJJ9LcsYiY/5X\n//y1SY6ffJuSpFFLhneSTcB5wHbg4cBzkhy7YMzTgIdU1dHALwJvmnybc5Pf5LquO8vas6o7y9qz\nqju52kmm/jUpc3NzE9vWRq47bs/7BODGqrq5qnYDFwOnLBjzDOBCgKq6Arh3ksMm2+bcZDe37uvO\nsvas6s6y9qzqTrp2reDr7BWOn5x9LURnVXfc6fFHALtGlm8BTlzGmCOB29bcnaR1bbV75K9+9atX\nNH5v13NZTe1J1F0vxoX3cjtf+FNc9Ps2bYKDDvpN9t//j5e5afjOdz7LgQdevayxe/Z8mc2bD132\ntiWt1UoDbkf/tVxLhfRKak+y7uxlqXeWJCcBO6pqe798JnBHVZ07MubNwFxVXdwv3wA8oapuW7Ct\n9fsWJknrWFXd5Z1k3J73VcDRSY4C/gE4FXjOgjGXAqcDF/dh/42Fwb1YcUnS6iwZ3lW1J8npwAeB\nTcAFVbUzyUv658+vqj9L8rQkNwLfBn5h6l1L0ga35LSJJGl9WvdnWC7nJKEp1HxrktuSfHqIegtq\nb0ny0SR/m+T6JP91oLoHJrkiyTVJPpPkN4aoO1J/U5JPJXnvwHVvTnJdX/tvBqx77yTvTrKz/3mf\nNFDdh/Wvdf7rmwP+jp3Z/15/Osm7ktxtoLov62ten+RlU651l+xIcp8kH0ryd0kuS3LviRSrqnX7\nRTdVcyNwFLA/cA1w7AB1twLHA5+ewWs+HHhk//gQ4LNDvOa+3sH9v5uBy4HHD/i6fxV4J3DpwD/v\nLwD3mcH/84XAC0d+3veaQQ/7AV8EtgxQ6yjg88Dd+uU/BF4wQN0fBj4NHNjnyYeAB0+x3l2yA3gd\n8Mr+8RnAOZOotd73vJdzktDEVdXHgH+cdp1Fan+pqq7pH98O7AQeMFDtf+4fHkD3i/71IeomORJ4\nGvAWZnN81qA1k9wL2FpVb4Xus6Wq+uaQPfROBm6qql1jR67dt4DdwMFJNgMHA7cOUPcY4Iqq+k5V\nfRf4C+Cnp1Vskey480TG/t9nTqLWeg/vvZ0AdMSMehlcf5TP8cAVA9XbL8k1dCdYfbSqPjNEXeCN\nwCuAOwaqN6qADye5KsmLB6r5g8BXkvx+kk8m+b0kBw9Ue9SzgXcNUaiqvg68Afh7uiPXvlFVHx6g\n9PXA1n7q4mDgJ+lOIhzSYfXvR+DdBkzkDPT1Ht4b9tPUJIcA7wZe1u+BT11V3VFVj6T75f6xJNum\nXTPJTwFfrqpPMZu97sdV1fHAU4GXJtk6QM3NwKOA36mqR9EdpfVrA9S9U5IDgKcDyz9bbm31Hgz8\nCt30yQOAQ5I8b9p1q+oG4FzgMuD9wKeYzU7CfD8Tu97Aeg/vW4EtI8tb6Pa+92lJ9gf+D/COqrpk\n6Pr9n/DvAx4zQLnHAs9I8gXgIuCJSd42QF0AquqL/b9fAd5DN1U3bbcAt1TVlf3yu+nCfEhPBa7u\nX/cQHgP8dVV9rar2AH9C938/dVX11qp6TFU9AfgG3edIQ7otyeEASe4PfHkSG13v4X3nSUL9nsKp\ndCcF7bPSXbDhAuAzVfVbA9a93/yn4EkOAp5Mt5cyVVX136tqS1X9IN2f8R+pqtOmXRcgycFJ7tE/\nvjvwE3Qfbk1VVX0J2JXkof2qk4G/nXbdBZ5D92Y5lBuAk5Ic1P+OnwwMMi2X5Pv7fx8IPIuBpopG\nXAq8oH/8AmAiO2TjzrCcqVrkJKFp101yEfAE4L5JdgGvqqrfn3bd3uOA5wPXJZkPzzOr6gNTrnt/\n4MIk+9G9qb+9qv58yjX3ZsipssOA9/QXONoMvLOqLhuo9i8D7+x3Sm5iwJPb+jeqk4Gh5vipqmv7\nv6iuopu2+CTwuwOVf3eS+9J9YPpLVfWtaRUayY77zWcHcA7wR0leBNwM/KeJ1OoPX5EkNWS9T5tI\nkvbC8JakBhnektQgw1uSGmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUoHV9YSppkpIc\nRnflwicCr6W7RRbAU+iu+fwIuh2ap1TVs2fSpLRMXphKG0aS/wK8Gfg7uvsIXtCv/zLw8qp6e798\nG/BDVfXVmTUrjWF4a8Po7x15T+DKqpq/OP4WupsSzF/z+UF0Nw04fHadSuM5560No79D0Ml0dxCf\n9+QFy88HLkpyYJK7DdmftBKGtzaahWF9Mt39DeedCrwDeBGw/4B9SStieGujeQjdnZkWW76a7n6S\nNw1142dpNZzzlqQGuectSQ0yvCWpQYa3JDXI8JakBhnektQgw1uSGmR4S1KDDG9JapDhLUkN+v+f\nQhSzz2scsAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import scipy as sp\n", "\n", "%matplotlib inline\n", "\n", "from scipy.stats import binom\n", "\n", "n = 10\n", "mu = 0.25\n", "\n", "yy = binom.rvs(n, mu, size=1000)\n", "\n", "fig, ax = plt.subplots()\n", "\n", "ax.hist(yy, bins=range(11), normed=True, rwidth=0.8)\n", "\n", "ax.set_xlabel(\"$m$\", fontsize='x-large')\n", "\n", "ax.set_xlim(0, 11)\n", "\n", "ax.set_yticks(np.arange(0, 0.31, 0.1))\n", "\n", "ax.set_xticks(np.arange(0.5, 10.6, 1))\n", "ax.set_xticklabels(range(11))\n", "\n", "ax.set_title(r'$N = 10, \\mu=0.25$', fontsize='x-large')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "其均值和方差分别为：\n", "\n", "$$\n", "\\begin{align}\n", "\\mathbb E[m] & = N\\mu \\\\\n", "\\text{var}[m] & = N\\mu(1-\\mu)\n", "\\end{align}\n", "$$\n", "\n", "计算均值时考虑下式对 $\\mu$ 的导数，方差考虑对 $\\mu$ 的二阶导数： \n", "\n", "$$\n", "\\sum_{m=0}^N \\begin{pmatrix} N \\\\m \\end{pmatrix} \\mu^{m}(1-\\mu)^{N-m} = 1\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.1.1 beta 分布" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "之前看到，当数据量较少时，最大似然的结果很可能会过拟合。为了减少这样的情况，从 Bayes 概率的观点出发，我们引入一个关于 $\\mu$ 的先验分布 $p(\\mu)$。\n", "\n", "我们观察到似然函数是一系列 $\\mu^x(1-\\mu)^{1-x}$ 形式的乘积，如果我们选择一个正比于 $\\mu$ 的某个幂次和 $1-\\mu$ 的某个幂次的先验分布，那么对 $\\mu$ 来说，后验分布应当满足同样的形式。这样的性质叫做共轭性（`conjugacy`）。\n", "\n", "在这里，我们引入的是 $0-1$ 间的 `beta` 分布：\n", "\n", "$$\n", "\\text{Beta}(\\mu~|~a,b)=\\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)} \\mu^{a-1} (1-\\mu)^{b-1}\n", "$$\n", "\n", "其中：\n", "\n", "$$\n", "\\Gamma(x) \\equiv \\int_0^{\\infty} u^{x-1}e^{-u} du\n", "$$\n", "\n", "满足如下性质：\n", "\n", "$$\n", "\\begin{align}\n", "\\Gamma(x+1) & = \\int_0^{\\infty} u^{x}e^{-u} du = \\left[-e^{-u}u^x\\right]_0^{\\infty} + x \\int_0^{\\infty} u^{x-1}e^{-u} du = 0 + x\\Gamma(x) = x\\Gamma(x) \\\\\n", "\\Gamma(1) & = \\int_0^{\\infty} e^{-u} du = \\left[-e^{-u}\\right]_0^{\\infty} = 1\n", "\\end{align}\n", "$$\n", "\n", "验证它是一个概率分布：\n", "\n", "由定义\n", "\n", "$$\n", "\\Gamma(a)\\Gamma(b) = \\int_0^\\infty x^{a-1}e^{-x} dx \\int_0^\\infty y^{b-1}e^{-y} dy\n", "$$\n", "\n", "令 $t = y + x, dt = dy$，则有：\n", "\n", "$$\n", "\\Gamma(a)\\Gamma(b) = \\int_0^\\infty x^{a-1} \\left\\{\\int_x^\\infty (t-x)^{b-1}e^{-t} dt \\right\\} dx \n", "$$\n", "\n", "交换积分次序，原来 $t$ 是从 $x$ 积分到 $\\infty$，现在 $x$ 是从 $0$ 积分到 $t$：\n", "\n", "$$\n", "\\Gamma(a)\\Gamma(b) = \\int_0^\\infty \\int_0^t x^{a-1} (t-x)^{b-1}e^{-t} dxdt \n", "$$\n", "\n", "令 $x = t\\mu, dx = td\\mu$，则有\n", "\n", "$$\n", "\\begin{align}\n", "\\Gamma(a)\\Gamma(b) & = \\int_0^\\infty e^{-t} t^{a-1} t^{b-1} tdt \\int_0^1 \\mu^{a-1} (1-\\mu)^{b-1} d\\mu \\\\\n", "& = \\Gamma(a + b) \\int_0^1 \\mu^{a-1} (1-\\mu)^{b-1} d\\mu\n", "\\end{align}\n", "$$\n", "\n", "于是我们有：\n", "\n", "$$\n", "\\int_0^1 \\text{Beta}(\\mu~|~a,b) d\\mu = 1\n", "$$\n", "\n", "其均值和方差为\n", "\n", "$$\n", "\\begin{align}\n", "\\mathbb E[\\mu] & = \\frac{a}{a+b} \\\\\n", "\\text{var}[\\mu] & = \\frac{ab}{(a+b)^2(a+b+1)}\n", "\\end{align}\n", "$$\n", "\n", "求均值：\n", "\n", "从归一化我们知道：\n", "\n", "$$\n", "\\int_0^1 \\mu^{a-1} (1-\\mu)^{b-1} d\\mu = \\frac{\\Gamma(a)\\Gamma(b)}{\\Gamma(a + b)}\n", "$$\n", "\n", "从而，利用 $\\Gamma(x+1) = x\\Gamma(x)$：\n", "\n", "$$\n", "\\mathbb E[\\mu] = \\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)} \\int \\mu^{a+1-1} (1-\\mu)^{b-1} d\\mu \n", "= \\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)} \\frac{\\Gamma(a+1)\\Gamma(b)}{\\Gamma(a+b+1)} = \\frac{a}{a+b}\n", "$$\n", "\n", "类似可以得到：\n", "\n", "$$\n", "\\mathbb E[\\mu^2] = \\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)} \\int \\mu^{a+2-1} (1-\\mu)^{b-1} d\\mu \n", "= \\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)} \\frac{\\Gamma(a+2)\\Gamma(b)}{\\Gamma(a+b+2)} = \\frac{a(a+1)}{(a+b)(a+b+1)}\n", "$$\n", "\n", "从而可以计算出方差。\n", "\n", "$a$ 和 $b$ 被叫做超参，因为它们控制分布的参数 $\\mu$。" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAG/CAYAAABIVpOQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VOXZx/HvYQmrEHaQTQUUBRHEBRA1RMAFqUpdUUSt\nCu6CtvJardhqxYq1bhSLYAVRirgrIKgERRYLFRBkRwKIbErYQkISzvvHQyAgTGaSc+aZM+f3ua5z\nxQmTMzcm3LnnWe7HcV0XERERETmyMrYDEBEREUlkKpZEREREIlCxJCIiIhKBiiURERGRCFQsiYiI\niESgYklEREQkAhVLIiIiIhGUi+ZJjuNcBlQFmgFbXdcd5mtUIiIeUg4TkdJwimtK6ThOKrARSAVy\nga3A6a7rZvofnohI6SiHiUhpFTsN57puFtDedd0c11RW5QDH98hERDygHCYipRXVmiXXdRcDOI7T\nGchwXXeNn0GJiHhJOUxESiOqNUsAjuP0Aq4CHvAvHBERfyiHiUhJFbtm6ZAnO05V4Fug25HemTmO\no1N5RULGdd3ATGlFymHKXyLhFE0OK3YaznGcHo7jfL3/hruAzcCVEV4Ut1073HnzzH/rCsX12GOP\nWY9BVxyun3/GrVHjwOMgiCWHWf//q8vKpfwV3ita0axZKgAy9icdB2gMLIz4FRUqQE5O1EGISEDk\n5Jh/38ESew4TESmi2DVLrutOdhznBMdx7gGaAk+6rjsl4hdVrAi5uR6FKCIJIzfX/PsOkBLlMBGR\nIqJa4O3G2sBNI0uhk5aWZjsEiYdgjizFnsMkVJS/pDj+HHeikaXQUbIJiQCOLIkUR/lLiuNPsaSR\nJZHkFNCRJRGR0tDIkohETyNLIhJCGlkSkehpZElEQkgjSyISPY0siUgIaWRJRKKnkSURCSH/RpZU\nLIkkH40siUgIaRpORKKXk6NiSURCR9NwIhK93FxNw4lI6GhkSUSip5ElEQkhjSyJSPQ0siQiIaSR\nJRGJnkaWRCSENLIkItFT6wARCSGNLIlI9NQ6QERCSCNLh1m3bh39+vXj5Zdf5tZbb2XNmjURn79n\nzx769+/PV199FZ8ARWzSyFLCU04S8V45X+4a4JGlPn368MQTT9C5c2dmz57Nddddx6xZs4743Fde\neYXVq1czYcIEevfuHedIRSzQyFJCU04S8Yc/xVJAR5aWL1/OggUL6Ny5MwAdOnRg2bJlZGZm0rRp\n0189v1+/fgCMHz8+rnGKWKORpYSmnCTiDx13UsTChQtp3LjxIZ9r1KgRM2fOtBSRSILRyJKIhJB/\na5YCOA23efNmKlWqdMjnKleuzKZNmyxFJJJgNLIkIiHk35oln0eW9u3bx5NPPkmdOnUoKChg5syZ\n3HjjjVx44YUHnpOVlcXAgQNxXTfivQYMGECbNm3Ytm0bKSkph/xZSkoKO3fu9OXvIBI4GlnyhB/5\nS0T8E9gF3nfeeSdNmjShf//+5OXl8fvf/57nnnvukOekpqYyatSoqO95zDHH/Opz2dnZ1KxZs9Tx\niiQFNaX0hB/5S0T8E8gF3gsXLuT1119ny5YtgFmY3axZM+rWrVuq+zZo0IBdu3Yd8rndu3dTr169\nUt1XJGnouJNS8yt/iYh/Ajmy9Nlnn9GxY0eqVq0KwKeffkrXrl3Jzs6mfPnylC9fHoBt27bxwAMP\nRD2M3aFDB9avX3/g8/n5+WRmZtK6dWvf/i4igaKRpVLzK3+JiH/8HVlyXXAcz29fs2ZN6tevD0BO\nTg5vvfUWgwYNYuzYsfTt2/fA82rUqBHTMHbjxo05/vjjmTNnDmeffTZTp07l9NNPp2XLlgDMnTuX\nHTt2kJ6e/quvLS6hiSQFjSyVml/563DKSSLecbz8B+U4jnvgfuXLQ3a2+eix3Nxc7rjjDrp168ae\nPXvYsWMH69evp127dlx//fWluvfSpUt54oknOPPMM5k3bx5/+ctfDvRYGjhwIFu2bGHMmDEAvP76\n60yaNIkJEybQpk0bzj33XIYOHXrgnaFI0mnYEL75xnwEHMfBdV3v3xFZcEj+8pGf+Us5SSQ20eYw\n/4qlqlXhp5/gCIumRSSgateGpUvNR1QsiUiwRZvD/OmzBIE+8kREjkJrlkQkhPwrlgJ65ImIRKCm\nlCISQv6OLKlYEkke+flm00Y5f/aFiIgkKn9HljQNJ5I8CnfC+bDDVUQkkWlkSUSio6NORCSktMBb\nRKKjxd0iElJa4L3funXr6NevHy+//DK33nora9asKfZr9uzZQ//+/fnqq6/8D1DENjWkTGjKRyL+\n8W+lZsBGlvr06cMTTzxB586dmT17Ntdddx2zZs066vNfeeUVVq9ezYQJE+jdu3ccIxWxRCNLCUv5\nSMRfgR1ZGjZsGNWrV2fSpEmlvtfy5ctZsGABnTt3BqBDhw4sW7aMzMzMo35Nv379ePrppzlGTTcl\nLDSy5Bkv8xcoH4n4LbBrlvr27YvjOHTp0qXU91q4cCGNGzc+5HONGjVi5syZpb63SNLQyJJnvMxf\nIuK/wI4sTZs2jU6dOlHRg+S9efNmKlWqdMjnKleuzKZNm0p9b5GkoZElz3iZv0TEf4FdszRlyhQq\nVKjA+PHjmTNnDhdeeCHdu3cHICsri4EDBxZ76vaAAQNo06YN27ZtIyUl5ZA/S0lJYefOnb7FLxI4\nGlnyjJf5S0T851+x5PPI0pQpUxg+fDhpaWm0atWKXr16sWzZMgBSU1MZNWpU1Pc60jx/dnY2NWvW\n9CxekcDTyJJnvMxfIuK/QDalzMzMJDc3l7S0NAA2btzIli1bSny/Bg0asGvXrkM+t3v3burVq1ea\nMEWSi0aWPOF1/hIR//k7suTTNNy8efMOJBqAyZMn07NnzwOPt23bxgMPPBD1MHaHDh1Yv379gc/n\n5+eTmZlJ69atPY9dJLA0suQJr/OXiPjP3zVLu3f7cuvU1FRq1KgBmHdlH374IV988cWBP69Ro0ZM\nw9iNGzfm+OOPZ86cOZx99tlMnTqV008/nZYtWx54zty5c9mxYwfp6em/+vrikppIUtDIkie8zl+H\nUz4S8V4gD9JNT0+nXLlyjB49mscff5yPP/6Yhg0bluqeo0eP5sUXX+T555/nrbfeYuzYsYf8+Ztv\nvslrr7124PHrr7/Otddey7p16xgwYAD33XcfeXl5pYpBJKHl5GhkyQN+5C/lIxF/OV6+C3Ecxz1w\nv5degiVL4OWXPbu/iFg0dChs3Gg+7uc4Dq7rOhaj8swh+UtEQiHaHBbYppQiEmeahhORkApsU0oR\niTMt8BaRkNLIkohERyNLIhJSGlkSkehoZElEQkojSyISHY0siUhIRdVnyXGc3kAD4CzgPdd1xxX7\nRRpZEkkuAR5ZKlEOExHZr9hiyXGc5kAt13WfdRynNrDCcZw5ruv+EPELfTzuREQsCOjIUolzmIjI\nftFMw7UC/gDguu5WYCXQvtiv8rEppYhYENyRpZLlMBGR/aIpliYCFwM4juNghrJXFvtVGlkSSS4B\nHVmipDlMRGS/YqfhXNfNAxbtf9gDmOu67vxi76yRJZHkEtCRpRLnMBGR/aI+SNdxnFTgJuCGSM8b\nPHiw+Y+dO0nbsYO0EocmIgklJ4eM778nY/p025GUSDQ57ED+AtLS0khLS/M7LBGJo4yMDDIyMmL+\nuqjOhts/dP0UMMR13SzHcZq6rpt5hOcdPFtp2zY44QTzUUSC7/TTYeRIaNfuwKeCcjZcNDlMZ8OJ\nhI/XZ8PdA7wNVHQc5yzguGK/Qq0DRJJLTk4gp+H2iz2HiYjsF03rgM7Ac0Bh5eUCTYq9c+GaJdcF\nJ+HfeIpIcQK6wLvEOUxEZL9oFnjPAMrGfOeyZc2VlwcpKSWJTUQSSXAXeJcsh4mI7OffcSegI09E\nkklAR5ZERErL32JJ65ZEkkdAR5ZERErL/5ElFUsiyUEjSyISUv6PLGkaTiT48vPNZo1yUbdmExFJ\nGhpZEpHi5eZqVElEQksjSyJSPK1XEpEQ08iSiBRP65VEJMQ0siQixdPIkoiEmEaWRKR4GlkSkRBT\nU0oRKZ4WeItIiKkppYgUL9iH6IqIlIq/xVKVKrBrl68vISJxsHs3VK5sOwoRESv8LZZq14atW319\nCRGJgy1boE4d21GIiFjhb7FUty5s3uzrS4hIHGzebP49i4iEkL/FUp065h2piASbRpZEJMT8H1lS\nsSQSfFu2aGRJRELL/5ElTcOJBN/mzRpZEpHQ0jSciBRP03AiEmLxKZZc19eXERGfaYG3iISYv8VS\npUqQkgI7dvj6MiLiM40siUiI+VssgRZ5iwRdXh7s3Ak1atiORETECv+LJS3yFgm2rVuhVi0o43+6\nEBFJRBpZEpHINAUnIiEXn5ElFUsiwaXF3SIScpqGE5HINLIkIiGnaTgRiUwNKUUk5DSyJCKR6agT\nEQk5jSyJSGQaWRKRkNMCbxGJTGuWRCTkNA0nIpFpGk5EQi5+I0s6H04kmDQNJyIh53+xVLGiubZv\n9/2lRMQHGlkSkZCLz/kFWuQtEkx798KuXZCaajsSERFr4lMsaZG3SDDpXDgRkTgWS1rkLRI8moIT\nEdE0nIhEoMXdIiIaWRKRCNRjSUREI0siEsHmzZqGE5HQ0wJvETk6jSyJiGgaTkQi0MiSiIim4UQk\nAo0siYhoZElEIlCxJCISp2Kpfn3Ytg327InLy4mIR1avhqZNbUchImJVfIqlcuWgeXNYtiwuLyci\nHti2zRx10rix7UhERKyK3xkGJ58M338ft5cTkVJasgRatgTHsR2JiIhV8S2WliyJ28uJSCktWQKn\nnGI7ChER6+JXLJ1yikaWRILk++/NmxwRkZDTyJKIHNmSJSqWRESIZ7F04olmZ01eXtxeUkRKQdNw\nIiJAPIulihXNrpqVK+P2kiJSQrt3w8aNcPzxtiMREbEupmLJcZy2juMMLfGraSpOJBiWLYMWLUzb\njyRR6vwlIqEVdSZ0HGcg0BnYXuJXU7EkEgxJtl7Jk/wlIqEV9ciS67p/Bz4o1atpR5xIMCTZTjhP\n8peIhFasa5ZK151OI0siwZCci7vVXVNESiTWYskt1au1bGnWQuzbV6rbiIjPkmwabr/S5S8RCa1Y\nV28W+85s8ODBB/47LS2NtLS0g39YrRrUrAmZmdplI5Ko9u6FH34w7T4Ok5GRQUZGRvxj8kbx+avI\n0S5p+y8RSR4Z+69YOa4b/Zstx3H6Ammu6958lD93i71f9+5w//1wySWxxCki8fL993D55bB8ebFP\ndRwH13UDMb3lSf4SkaQSbQ6L75olMOsgFi8u9W1ExCeLFyfjFBxozZKIlFDUxZLjOHcDtwBpjuM8\n5jhOtRK9YseO8NVXJfpSEYmDr74y/06TiGf5S0RCKaZpuGJvFs0w9pYtptnd1q1J1fBOJGm0agWv\nvw5nnFHsU4M0DVccTcOJhI9f03ClV6cONG0Kc+fG/aVFpBg//WSudu1sRyIikjDiXywBdO0Kn31m\n5aVFJILPP4e0NChb1nYkIiIJw06xdMEFJimLSGL5/HPz71NERA6I/5olgF27oEED2LQJKlf27PVF\npBRc10yRT50KJ50U1ZdozZKIBFnirlkCqFoV2raFr7+28vIicgQrV5ru+kdoRikiEmZ2iiUwQ/1a\ntySSOD77zPy7dJJioEhExDN2iyWtWxJJHFqvJCJyRHbWLIE5f6puXVi6FOrX9ywGESmB3Fw49lhY\nuBAaNoz6y7RmSUSCLLHXLAGkpMBvfwujR1sLQUT2++ADs44whkJJRCQs7BVLAL/7HYwcaXbhiIg9\nI0eaf48iIvIrdouljh2hTBmYMcNqGCKhlplpOupfcYXtSEREEpLdYslx4NZbzbtaEbHjtdegd2+o\nVMl2JCIiCcneAu9Cmzebvi5r10I1HQQuElcFBXDCCQfXLMVIC7xFJMgSf4F3obp1zXblt96yHYlI\n+Hz+OdSuXaJCSUQkLOwXSwB33w3PPgt5ebYjEQkP14UhQ+DOO21HIiKS0BKjWOrSxZxJNWqU7UhE\nwmPKFNiwAfr2tR2JiEhCs79mqdDcuXDZZbBihQ7XFfHbvn3Qvj088ojpd1ZCWrMkIkEWnDVLhc44\nA845B55/3nYkIsnvP/8xjWF79bIdiYhIwkuckSWA5cuhUydYsgTq1PEsLhEpIicHWreGESPMFHgp\naGRJRIIseCNLYFoI3HIL3H67unqL+OXhh6Fdu1IXSiIiYZFYxRLAX/4CP/xgGuWJiLc++wzGj4fh\nw21HIiISGIk1DVdo8WJIS4NZs6B589LfT0Tg559NP6VRo6BbN09uqWk4EQmyYE7DFWrVCh59FK65\nBnbtsh2NSPDl58ONN8KVV3pWKImIhEVijiyBWbN0222wbh189JHZuSMisSv8t7R+vfm3VL68Z7fW\nyJKIBFmwR5bAHLI7fLgpkn73O9MXRkRi99hjsGABTJjgaaEkIhIWiVssAZQrZ/rBrF5t3hnn59uO\nSCQ4XBcefxzGjYNPPoGqVW1HJCISSIldLIHp5v3pp/Djj3DFFZCdbTsikcSXnw/9+8OHH8JXX5kD\nq0VEpEQSv1gC8474o4+gRg3TG2btWtsRiSSuX36Byy83LTgyMqBePdsRiYgEWjCKJTBrLV5/3Zxj\ndeaZ8P77tiMSSTxff20aTp54Inz8MRxzjO2IREQCL3F3w0UyezZce63ZAv3001Czpv+vKZLIsrPh\niSdMD6VXX4VLL43Ly2o3nIgEWfB3w0XSoYPZ3VOxounJNHq0dstJOLmuWbzdurWZdvv227gVSiIi\nYRHMkaWi5s6Fu+6C3Fx48km45BLTdkAk2c2cac5527gRXngBunePewgaWRKRIIs2hwW/WALz7vr9\n9+GRR6BKFfj9783OuXLl4h+LiJ9cF6ZOhaFDYdkyGDwY+vSx9rOuYklEgixcxVKhggKzVXroUPjp\nJ7j9dujbFxo0sBeTiBe2bYM33zx4AO4DD8B110GFClbDUrEkIkEWzmKpqG++gREjTNfic881v1h6\n9lRjPgmO3FzTY2zcOJg4ES66CG69FS64IGGmmlUsiUiQqVgqtHMnvPee6QQ+YwakpcFvfgM9ekD9\n+rajEznUtm0webIZIf30Uzj1VLPz86qroHZt29H9ioolEQkyFUtH8ssvMGmS+UU0ZQo0bmzaD3Tp\nAuecY5peisTTrl0waxZMn27WIi1ZAuefb0ZBL70Ujj3WdoQRqVgSkSBTsVSc/Hyzk27qVPOLas4c\nOO4405agQwfT+PLkk3XwqHinoABWrDA/d3PmmH5hS5aYJpLnnQddu0KnTtbXIcVCxZKIBJmKpVjl\n5cH8+eaX2Jw5MG8erFkDp5wCbdqYPjatW0PLlmZEKkHWjEgCcl3YtMkUQosXw3ffHbzq1oX27eHs\ns811xhmmX1hAqVgSkSBTseSFXbtg4ULzS27RIvOLb+lSsw6qeXNo0cJ8POEEOP54MzLVqFGgRgak\nhPLzzeHOa9aYZpCrV8PKleZasQLKljUjk6ecYtYdtW4NbdtCaqrtyD2lYklEgkzFkp+2bze/EFes\nML8cC39ZZmbChg3m+JXGjU3h1LChaV3QoIFZUF6vnrnq1IGUFNt/Ezlcfj5s3WpGhgqvDRtMK4oN\nG2D9eli3zny+Xj1o2tQUyscff7CAbtEiIRdj+0HFkogEmYolWwoKzC/WH380v1Q3bDh4bdpkui1v\n2gQ//2zaGNSuDbVqmY81a5pF5kWv6tUPXsccA9Wqma+rVElTgUfiumbL/c6dB6/t2w9e27ZBVpb5\n+Msv5vr5Z3Nt3Wr+rGZNM11WWNwee6y5GjQwRXDDhubSejYVSyISaCqWEt2+feYXdtFf1IW/uLOy\nDv5CL/wlv2PHwWvnTjMCUqXKr6/KlU0hVXhVrHjwqlDBXCkph17ly5urXLmDV9myBz+WKXPwo+Mc\n/Ajmo+OYIgXMx8Jr376DV0HBoVd+vrny8syVnw979x565eaaKyfn0Cs7G/bsMR+zs2H3bnPt2mU+\nliljCspq1Q5ehQVnjRpmKiw11RSpNWseLFYLH5cta+/nImBULIlIkKlYSnZ5eYcWCYVFQ9FComiB\nUVh0FC1E8vIOfiwsXPLzTTGTl3ewsCla8BQWQXCwKCpUtIAqLKgKi6zCgqvwKl/+4MfCQq1CBfPf\nRQu6wiKvaNFXWBBWrnywQKxSxRRIVapoejOOVCyJSJCpWBIR36lYEpEgizaHlYlHMCIiIiJBpWJJ\nREREJAIVSyIiIiIRqFgSERERiaBcNE9yHOdy4BRgH/Cj67pjfI1KRMRDymEiUhrF7oZzHKc68IXr\nuu33P54F9HRdd+sRnqvdJCGVkZFBWlqa7TAkzoKwGy7aHKb8FV7KX+Hl5W6484DvizxeAHQpaWCS\nnDIyMmyHIHI0ymESkfKXFCeaYqkRkFXkcRbQwp9wREQ8pxwmIqUSTbGUCuQUebwXqOpPOCIinlMO\nE5FSiWaB906gVpHHlYBNR3uyo8NdQ+vxxx+3HYLIkUSdw5S/wkv5SyKJplhaBZxR5HFt4H9HemKi\nL/QUkVCKKocpf4nI0UQzDfcl0L7I49OBz/0JR0TEc8phIlIqUR2k6zhOH6Apprha5bruWL8DExHx\ninKYiJRGVMWSiIiISFhF1cG7OOqOGw7RfJ8dx1nFwa3av3ddd3R8o5R4chynLXCD67oP2o6lpJS/\nwkM5TIqKJX+VemQplg7fElwxdEG+DfgU2OC6bn78I5V4cRxnINAZ2O667s224ykJ5a/wUA6TomLN\nX14cpKvuuOEQ7fd5r+u6a5Vkkp/run8HPrAdRykpf4WHcpgcEGv+8mIaTt1xwyHa7/OZjuNUAKoB\ny13X/TAewYk1Qd9ur/wVHsphcrio85cXxZK644ZDtN/nz13XfQ/AcZz5juN86bpu1hGeJ8kh6DtE\nlL/CQzlMDhd1/vJiGm4nh1ZnlYBfPLivJJZov89FhzW3AWk+xiT2BX1kSfkrPJTD5HBR5y8viqVV\nmI64hWoDP3pwX0ksxX6fHce5ARhX5FNVAM37J7egjywpf4WHcpgcLq4jS+qOGw5H/T47jtPMMYdq\nrQGG7/9cZaAO8EV8w5Q4C/rIkvJXeCiHyeGizl+eNKVUd9xwONr32XGc/wG/c133W8dxrsckmKbA\nONd151gLWHzlOM7dwNVAY+DfwHOu6+6wGlQJKH+Fh3KYFIo1f6mDt4iIiEgEXkzDiYiIiCQtFUsi\nIiIiEahYEhEREYlAxZKIiIhIBCqWRERERCJQsSQiIiISgYolERERkQhULImIiIhEoGJJREREJAIV\nSyIiIiIRqFiSEnMc50bHcZ4s8vg2x3EetRmTiEi0lMMkWiqWpDR6AfOKPL4KWGApFhGRWCmHSVR0\nkK6UiOM4ZYFNQDPXdbc7jlMR2Aw0dF13p93oREQiUw6TWGhkSUqqPbDadd3t+x+fCywCdjuOU9Ne\nWCIiUVEOk6ipWJKS6gqsLvK4N/AlkA7UtxKRiEj0lMMkalFNwzmOcxlQFWgGbHVdd5jfgUlicxzn\nc6AAeBMoD+QAFwDf6OdDEo1ymBxOOUxiUWyx5DhOKrARSAVyga3A6a7rZvofniSi/XP7PwL1XdfN\nsx2PSCTKYXI45TCJVbHTcK7rZgHtXdfNcU1lVQ5wfI9MEllnYL6SjASBcpgcgXKYxCSqNUuu6y4G\ncBynM5Dhuu4aP4OShHcq8L7tIESipRwmh1EOk5hE3TrAcZxemB4Uj7quu9LXqEREPKYcJiIlFVOf\nJcdxqgLfAt2O9M7McRw1bRIJGdd1AzOlFSmHKX+JhFM0OaxccU9wHKcH8LDruue4rrvLcZzNwJXA\n0KO8aMyBSvANHjyYwYMH2w5D4sxxEr9OiiWHKX+Fk/JXeEWbw6JZs1QAZOy/qQM0BhaWNDARkThT\nDhORUil2ZMl13cmO45zgOM49QFPgSdd1p/gfmohI6SmHSdLbtQtmzICdO6FGDTj2WDj5ZAjAyG9Q\nFFssAahBlxQnLS3NdggiR6UcJpEEMn+5Lrz/Pjz7LMyfD2ecAbVqwbZtsGqVKZr69YM+faBqVdvR\nBp6nB+k6juNqzl8kPBzHCdQC70iUvyQwFi6E++6DLVvgySehe3eoVOngn+/bB198AS+/bJ47fjy0\nb28v3gQWbQ7T2XAiIiJB8fbbcMEFcPXVZkTpsssOLZQAypSBrl3hvfdgyBC46CJ46SUzGiUlEtU0\nnIiIiFjkuvDMM/DiizBlCrRrF93XXXWVee5ll8Hu3fDQQ/7GmaRULImIiCS6xx+Hd9+FWbOgUaPY\nvrZ5c5g6Fc45B+rWhZtv9ifGJKZiSUREJJGNHAljxsDMmVCvXsnuceyxMHkypKVBnTpw6aWehpjs\ntMBbREpMC7xFfDZpkhkJ+vJLOPHE0t9v9mz4zW/g22+hYcPS3y/gos1hKpZEpMRULIn4aPVqOPts\n+OAD6NTJu/s+/riZzps0KfS9mLQbTkREJKjy8qB3b/jjH70tlAAefhh++QWGD/f2vklMI0siUmIa\nWRLxyR//aKbKPvnEn9GfpUuhc2f473/h+OO9v39AaBpORHynYknEB9Onw3XXmWKppAu6o/GXv8Di\nxTBunH+vkeBULImI71QsiXgsJwfatIGhQ81CbD9lZ5tF4xMmQIcO/r5WgtKaJRERkaB5+mlo1cr/\nQgmgcmV44gl44AF19y6GiiUREZFEsGKF6dD9wgvxe80+fUxn73ffjd9rBpCm4USkxDQNJ+IR1zUH\n4l58MQwcGN/XnjoV7rnHrF8qWza+r21ZtDlMHbwPs2rVKkaMGMGePXv4/vvvGTJkCO11WrOIBMCW\nLVt48cUXqVatGj/++CM9e/YkPT3ddlgSjY8+gg0b4N574//aXbtC9eqmn1OvXvF//QBQsVTEvn37\nePbZZ3n55ZdxHIcxY8bQvXt3li1bRu3atW2HJyIS0Z///GeeffZZUlJS2LdvHz169KBt27bUrFnT\ndmgSSX6+OeB26FAoZ+HXsuPAoEHw1FNwxRWhb1R5JFqzVMSKFSuYOXMmGzZsAOC6665j9+7dvP/+\n+5YjExF3OvrvAAAgAElEQVQp3rRp0yhfvjwAZcqU4bTTTmPt2rWWo5JijRwJDRrAJZfYi+Gyy2DH\nDpg2zV4MCUzFUhFVq1Zl/fr1/PTTTwCUK1eOqlWrsm3bNsuRiYgUr6CggBtuuIHs7GxycnJYtmwZ\nrVu3th2WRLJrlzl+5Jln7I7olCljRreGDLEXQwLTAu8I1q1bR9OmTfnyyy/p3Lmz7XBEEo4WeCeW\nKVOmcPnll9O0aVM6d+7MQw89RPPmzW2HJZH8+c+wfDm88YbtSGDvXmjWzKxdOv1029HERdI3pdy3\nbx9PPvkkderUoaCggJkzZ3LjjTdy4YUXHnhOVlYWAwcOpLiYBgwYQJs2bX71+UGDBrF48WI++ugj\nz+MXSQYqlkrGr/yVnZ3N/fffz4wZM1i9ejXPPfccd9xxh69/FymFrCxo3hzmzDFFSiIYMgRWroRX\nX7UdSVwkfbHUv39/mjRpwsMPP0xeXh7Vq1dnzZo11K1b15P7L1q0iD59+vD5559rcaTIUahYKhm/\n8te1117L3/72N+rWrcvAgQN55ZVXmDFjBh07dvQocvHU4MGwdi2MGmU7koM2bYKWLeGHHyA11XY0\nvkvq1gELFy7k9ddfZ8uWLQAsX76cZs2aeVYo7dy5k0GDBjFx4kQVSiLiKb/y15dffskJJ5xAkyZN\nABg2bBh169Zl/PjxKpYSUVYWvPSSGVVKJPXqwYUXmmnBu++2HU3CCGSx9Nlnn9GxY0eqVq0KwKef\nfkrXrl3Jzs6mfPnyB3aDbNu2jQceeCDmabhBgwbx0ksv0aBBA1zXZdy4cVx33XX+/YVEJDT8yl9b\ntmyhYcOGh/zZVVddxejRo/35i0jp/OMf0LNn4ky/FdW/vymU7rpLbQT2C2SxVLNmTerXrw9ATk4O\nb731FoMGDWLs2LH07dv3wPNq1KjBqBiHN//2t7/RqFEjli5dytKlS1m7di0NGjTwNH4RCS+/8ld6\nejrXXHMNN910E1WqVAFg4sSJ9OnTx9u/gJRe4ajS7Nm2Izmy8883vZ++/hq0uQkI6Jql3Nxc7rjj\nDrp168aePXvYsWMH69evp127dlx//fUlvu+SJUto06YNBQUFBz7nOA6LFy+mZcuWXoQuklS0Zil2\nfuUvgHnz5jF8+HAaNmzInj17SEtL4+KLL/YocvHMkCGwaFFi7IA7mn/8A+bOTewYPZD0C7xFxD4V\nSyIxysmB44+HTz+FI+zCThi//GLizMxM6oXe0eYwNaUUERGJl9GjoV27xC6UAGrWNGfGvfOO7UgS\ngoolERGReCgoMOe/PfSQ7Uiic8MNMGaM7SgSgoolERGReHj/fTNic955tiOJziWXmLVVOl9QxZKI\niEhcDB0Kf/hDcLbjV6gAV14Jb75pOxLrVCyJiIj4bfZs0x37sstsRxKbwqm4kG9+ULEkIiLit+ef\nh3vvhbJlbUcSm3POgexsmD/fdiRWqXXAfgsWLODdd9/lmGOOYeHChdxxxx06IkCkGGodkHg+++wz\n/vvf//J///d/tkORQuvXm91va9ZAtWq2o4ndI49Abi4884ztSDyn1gExuv3227n66qt58MEH6dOn\nDxdffDG//PKL7bBERKK2b98+HnzwQfLy8myHIkW9/DLceGMwCyWAq64yLQSS4M1ESQW2WBo2bBjV\nq1dn0qRJntwvPz+fZcuWAdCkSRN27NjBypUrPbm3iEhRXuevQuPHj6d27dqe3lNKKTsbXn0V7rnH\ndiQl16aNmT789lvbkVgT2GKpb9++OI5Dly5dPLnfvHnz6NWrFwCZmZlUrFiRk046yZN7i4gU5XX+\nAtizZw8///wzjRs3LvbwXYmjsWOhY8fEPDA3Wo5jdsVNmGA7EmsCWyxNmzaNTp06UbFiRc/vPWrU\nKJ555hmqV6/u+b1FRPzIX6NHj9ahuYnGdc0U3N13246k9K68Et5+O7RTceVsB1BSU6ZMoUKFCowf\nP545c+Zw4YUX0r17dwCysrIYOHBgse+uBgwYQJsiLednzpzJlClTKFeuHLfccouv8YtIeHmdvzZt\n2kTFihWpFtQ1Mclq5kwzDde1q+1ISu/00yE/H777LvGPavFBoIul4cOHk5aWRqtWrejVq9eBNUep\nqamMGjUq5nt26tSJTp06MXnyZDp27Mi0adOoUaOG16GLSMh5nb/GjBnDgAEDDjx2gtL0MNm9/DLc\neSeUCewkzkFFp+JCWCwF8juYmZlJbm4uaWlpAGzcuJEtW7Z4dv+LLrqIzMxMnn/+ec/uKSIC3uev\nxYsX06JFC8ru79/juq7WLCWCTZtg0iS46SbbkXincCouhAI5sjRv3rwDiQZg8uTJ9OzZ88Djbdu2\n8cADD0Q9jD179mx69erFnDlzaNy4MQDly5dnx44dvsQvIuHldf6aMWMGa9asYfbs2QDMmjWLVatW\nkZKSol5LNo0YYbbcp6bajsQ7Z50FO3fCsmUQsg1QgSyWUlNTD0yPbdy4kQ8//JAvvvjiwJ/XqFEj\npmHslJQUqlatSqVKlQBYunQpWVlZ9O7d29vARST0vM5f/fr1O+TxrFmz6NKliwolm/Lz4ZVX4KOP\nbEfiLceBnj3N30vFUuJLT09n8uTJjB49mlmzZvHxxx/TsGHDEt/v9NNPZ8iQIQwbNoy9e/eydOlS\nPvjgA8444wwPoxYR8T5/FdqyZQt/+tOfmD9/Pps2bSIvL48nnnjCg4glZpMmQcOG0Lat7Ui817Mn\nDBkCDz5oO5K40nEnIlJiOu5E5Ah69DBTcMm0XqlQTg7UqwerV0OtWrajKTUddyIiIhJvmZkwezZc\nfbXtSPxRsSKkp8PEibYjiSsVSyIiIl559VW4/nqoXNl2JP75zW/gww9tRxFXmoYTkRLTNJxIEXl5\n0LQpTJ0KrVrZjsY/mzfDiSea9ggVKtiOplQ0DSciIhJPH31kzoBL5kIJoG5dOOUUmD7ddiRxo2JJ\nRETECyNGwO23244iPnr2DNVUXFTTcI7j9AYaAGcB77muO+4oz9MwtkiIBGUaLpocpvwlpZKZac5P\nW78e9vfsS2oLFsBvfwsrV9qOpFSizWHF9llyHKc5UMt13Wcdx6kNrHAcZ47ruj94EaiIiJ+UwyQu\nRo2C3r3DUSiBOR9uzx5YsQJatLAdje+imYZrBfwBwHXdrcBKoL2fQYmIeEg5TPxVUGCKpdtusx1J\n/DgOXHSRacAZAtEUSxOBiwEcc5R1A0yyEREJAuUw8dfkyXDssWa0JUxULB3kum6e67qL9j/sAcx1\nXXe+v2GJiHhDOUx8N2JEuEaVCnXrBjNmmOm4JBf12XCO46QCNwE3RHre4MGDD/x3WlraIadri0iw\nZWRkkJGRYTuMEokmhyl/Scx++slsoR8zxnYk8Zeaas6/mz7djDIFQElzWLS74RzgKWCI67pZjuM0\ndV038wjP024SkRAJ0G64YnOY8peUyJAhsGqVGV0Ko7/+1TSnfP5525GUiNdNKe8B3gYqOo5zFnBc\nKWITEYk35TDxnuvCyJHwu9/ZjsSeiy8OxbqlaFoHdAaeAworLxdo4mdQIiJeUQ4T33z5pTnu4+yz\nbUdiz2mnwY4dZnStWTPb0fim2GLJdd0ZQNk4xCIi4jnlMPHNyJFw661mG31YlSkD3bub8/CSuFjS\ncSciIiKxysoyx33cEHHPUzh062aKpSSmYklERCRWb71lRlRq17YdiX1du8K0aZCfbzsS36hYEhER\niVXhFJxAgwbQqBHMnWs7Et+oWBIREYnFggWweTNccIHtSBJHkk/FqVgSERGJxciRcPPNUFb7Bg4o\nXOSdpKJqShn1zdTUTSRUgtKUMhrKXxKVnJyDU07HHWc7msSRnQ316sGGDXDMMbajiZrXTSlFRETk\n/fehXTsVSoerXNn0mwrocUjFUbEkIiISrbB37I6kWzeYMsV2FL7QNJwt2dkwfz4sWwbLl0Nmpjlf\nZ/Nm2L0bcnOhoABSUqBiRaheHerWNcOcxx9vmn+dfDKccorpICtigabhJFTWrIEzzoD1601elkPN\nm2f6Ti1ZYjuSqEWbw1QsxcuOHaYPxdSpMGOGKZBOPhlatYITTzQFUL16UKcOVKliCqCyZWHvXlM4\nZWXBli2wcSOsXm1ayy9ebD6eeCJ06gTnngtpaXDssbb/thISKpYkVB57DH75BV580XYkiamgwLyp\nX7gQGja0HU1UVCwlgu3bzfz2+PHw1VdmPrdbN1PQnHaaNyNCe/bAd9/B11+b15g+3Sw+vPhiuPxy\n85phbsUvvlKxJKFRUGDe1H74IbRtazuaxPXb35rfPX362I4kKiqWbPrmGxg+HN57D84/H665Bi69\nND47BPLzzet/8gm8846Z7rv6arjpJmjd2v/Xl1BRsSSh8emn8PDDZqpJjm7YMPM76N//th1JVFQs\nxZvrwsSJ8NRT8OOP0L+/6cNRt67dmBYtgjffhDFjTJfV22+H3r3NVJ9IKalYktC4+mozK3DnnbYj\nSWzLlpnjT9auDcSshoqleJo6FQYNMqM6gwbBVVdBuXK2ozpUQYGJ85//NGum+vaF+++HJk1sRyYB\npmJJQmHrVmje3CzwTk21HU1ic11o3Bi++MKsp01w6rMUD4sXm66ld90F//d/ZnfbddclXqEEZrH4\nRRfBBx+YYeQyZUyvkBtuMKNPIiJyZG+8YZZSqFAqnuOYY2A+/9x2JJ5SsVQSu3fDQw+ZIdnf/MYU\nTVdeGYghR8A0Uxs61Oyka93a/GBffbWKJhGRw7mueivFSsWSMH06nHqqWZf03Xdw991QvrztqEom\nNdVMG65eDWedZeaZ+/QxQ80iImIWK+fkmDfHEp0LLjCtcvbtsx2JZ1QsRSsnBx580CyOfvFFMyxb\nv77tqLxRpYr5u61YYZpdtm8PAwea3k4iImH26qtmVCkoMweJoGFD0zNw/nzbkXhGxVI0Vq0yTR/X\nrIEFC6BHD9sR+eOYY2DwYPj+e9i1C1q2NC0QCgpsRyYiEn87d8KECWZDjMTmggvMIu8koWKpOB98\nAB07wi23wNtvQ+3atiPyX7168K9/weTJ8NZbZopuzhzbUYmIxNf48XDeeabtisQmPd1MxSUJtQ44\nGteFJ56AESNMkXT22bYjssN1YexY+P3vzWL2p5/WjhA5QK0DJKl17GgaUfbsaTuS4Nm61Szr2Lo1\nodf1qnVAaeTkmC31H39sRlTCWiiBmacvPBjRccxZdu+9ZzsqERF/LV5sDji/+GLbkQRT7drmeJi5\nc21H4gkVS4fLyjK9k/LzISNDw6+FUlPN+qW33jI76K65xrxjEBFJRiNGmOUXidg3LyiSaCpOxVJR\nGzea7aHt2pmioFIl2xElnvPOMzscmjQxLRQ0yiQiySYnx+x4Vm+l0klPT5pF3lqzVCgz06zev+km\n+OMftU00Gl9/bf5/de4Mzz8P1arZjkjiTGuWJCmNHQujR5vDc6Xktm+HRo3MLESFCrajOSKtWYpF\nZiZ06QL33AOPPKJCKVrnnAPffmsW77VtCzNn2o5IRKT0RoyA226zHUXwVa8Op5wCs2fbjqTUVCwV\nFkr33WcuiU3VqqbNwHPPQa9e8OST6sskIsG1fDksXWp2/0rpJclUXLiLpU2bzBEf996rQqm0LrvM\nHNA7dar5f7phg+2IRERiN2KEaUKZkmI7kuTQpYuKpUDLyoILL4Qbb4T777cdTXJo2NAcntilC5xx\nBnz2me2IRESil5sLr78Ot95qO5LkUbhcY/du25GUSjiLpT17TJOx8883a5TEO2XLwp/+ZHaS3Hgj\nPP54Uh2mKCJJ7N134bTToEUL25EkjypVzA7zr7+2HUmphK9Y2rfP/BJv1Miss9Fibn+kp8P//meG\nX3v0gF9+sR2RiEhkw4dDv362o0g+XboEvt9S+Iqlhx+Gn36C116DMuH768dV/fpmKq5VK2jf3gzF\niogkoiVLzOLuyy6zHUnySYLmlOGqFl591Zwg/f77ULGi7WjCoXx5GDrUnCnXvTu8+abtiEREfu2V\nV0zH7gQ+xyywOnSARYtgxw7bkZRYeJpSfv01XHEFfPUVnHSS7WjC6bvvzPfgsstM8aRjBAJPTSkl\nKezZA40bm3PMjjvOdjTJKT0dHnjALMtIIGpKWdSPP8LVV5upNxVK9px6KnzzDSxcaP7BbNtmOyIR\nEfjPf+Css1Qo+SngU3HJXyzl5sKVV8IddyRcRRtKNWvCpElw8slmaHb5ctsRiUjYDRsGd91lO4rk\nFvB+S8k/DXfPPbBundkSqgXdieXVV805fG++ac7lk8DRNJwE3n//a2YeVq40rU/EH3v3Qu3asGaN\nedOcIDQNB/DOO/DJJ/Dvf6tQSkS33mqGv6+/3iyuFBGJt2HDoH9/FUp+S0mBTp1g+nTbkZRI8o4s\nrV5tpnk++QTOPNN2NBLJypVmirRnT7PwW0krMDSyJIH288/QvDmsWGFGPcRfTz9t1hC/8ILtSA4I\n98hSXh5ce63pqaRCKfE1bw6zZpmdKFdeGfi2+CISEK+9Zt6kqVCKjwAfqpucxdKf/2x++HU4bnDU\nrAlTpkC1amYh4KZNtiMSkWRWUAD//CfceaftSMKjXTtYvz6Q+T35iqWZM83C4VGjdJRJ0KSkmPVl\nl1wCHTuajroiIn6YNMm8STv7bNuRhEe5cuZM1gC2EEiuYmnnTujTx5zvU7++7WikJBwHBg82h/Gm\npZkmoiIiXnvhBbj3Xr2pjreA9ltKrgXe/fub9UojR9qLQbwzZYrZKffPf5q1TJJwtMBbAun77027\nkjVroEIF29GEy3ffQa9eZlF9Aog2hyXPeRNffAETJ5pvhCSH7t1h6lS49FKzg0Jr0ETECy+9BP36\nqVCyoVUr2L4d1q6FJk1sRxO15JiG27XL9Ox55RWoXt12NOKltm3NuX7Dh8NDD8G+fbYjEpEgy8qC\nceNMsSTxV6aMWWIRsKm45CiWHn4YzjsPLr7YdiTih6ZNYcYMs37pppvMVKuISEmMHGl+VzRoYDuS\n8ApgC4Hgr1maPdvMfy5alFAt1MUH2dmmf1ZeHkyYAFWq2I4o9LRmSQIlPx+aNTOnO5xxhu1owmv5\nclMwrVtnfYF9OJpS5uWZodS//12FUhhUrmzO+GvQwPxD27rVdkQiEiQTJsBxx6lQsq1FC1MkJcgi\n72gEu1h6/nmoVw+uucZ2JBIv5cqZYfT0dDj3XLNIUESkOK4Lzz4LAwfajkQcx+xG/Owz25FELaZi\nyXGcto7jDPUrmJisXQtDhphDENUnI1wcB556Cm6/HTp3VvNKiUpC5S+Jv6+/Nou7e/a0HYkAdO0K\nn39uO4qoRd06wHGcgUBnYLt/4cTg/vvNVvLmzW1HIrYMGAC1apnjUT74QJ145agSLn9J/D37rMkZ\nZYI9oZI00tNNU9CCgkAcnh71T43run8HPvAxluhNnQoLFsDvf287ErHtxhvNtFzPnubnQuQIEip/\nSfwtX2521PbtazsSKXTssWb96bff2o4kKrGW2Pbnu/LyzIjSc89BxYq2o5FE0KOH2d1yww3w9tu2\no5HEZT9/iR3PPAN33aUdtIkmQOuWYi2W7O+rffll0/VT885S1LnnmuNR7r8fRoywHY0kJvv5S+Jv\nwwbzZuruu21HIocL0LqlWI87Kfad2eDBgw/8d1paGmlpaTG+RARbtsCTT5rmhFrULYc77TSYPh26\ndYNffjEdv8VTGRkZZGRk2A6jpOzmL7HjH/8wB6zXrm07Ejnc+eeb8z9zcuI2U1TSHBZTU0rHcfoC\naa7r3nyUP/e3qdtdd0H58uaHX+RofvzRnCt36aVmx6QKa98EqSml9fwl8ZeVZZpQ/u9/5iQASTwd\nO5pBkPR0Ky/vV1NKe0lx6VIYPx4efdRaCBIQDRvCl1+as4f69TO7LUS0Zil8hg0zaxpVKCWugKxb\nirpYchznbuAWIM1xnMccx6nmX1hH8NBD5qpVK64vKwFVq5aZC1+1Cq67DvbutR2RWGQ9f0n87dpl\nGhcPGmQ7Eomke3ez3jTBBeNsuIwMuPlm03xQO+AkFjk5pljas8cs8tRuGE8FaRquOJqGSzJDh8J/\n/wv/+Y/tSCSSvXuhTh1YudJ8jLPkORvOdeEPf4C//lWFksSuYkXTTqBePbjwQrOGQUSSW3a2KZYe\necR2JFKclBRIS0v4qbjEL5bee8/0VtL5b1JS5crBa6+ZwzPT0mDTJtsRiYif/vUv6NQJTj3VdiQS\njQBMxSX2NFxBgflhHzoULrnEu/tKOLku/PnPMHas6fatRZ+lpmk4STg5OWYH3Ecfwemn245GorFi\nhXkju3593HcvJ8c03JgxZqHuxRfbjkSSgePAY4+ZFhTnnqsDeEWS0fDhZhRZhVJwNG8OFSrA4sW2\nIzmqWJtSxk9uLgweDG+8oT454q377oPUVNPX4+OPoX172xGJiBd27TK91RJ8SkcO4zgHp+Jat7Yd\nzREl7sjSyJFwyinQubPtSCQZ9e0L//ynGbWcPt12NCLihRdegC5doE0b25FIrLp3h08/tR3FUSXm\nmqXcXDMs9847cNZZpb+fyNF88QVce60pznXeYMy0ZkkSRlYWtGgBX38NJ55oOxqJ1fbt0KgRbN4M\nlSrF7WWDvWZp5EjzzkCFkvgtPR0++QRuu82skRORYHr2WfOGR4VSMFWvDm3bmr6KCSjx1izl5sJT\nT5lRJZF4OPNMczTKhReaA3jvu892RCISiw0bzNEm//uf7UikNHr0MG9eE3BTV+KNLGlUSWw4+WT4\n6iuTcB991LQZEJFg+NOf4NZb1Q4k6AqLpQTMv4m1ZmnvXrNWacIEFUtix+bN5l3NWWfBSy9B2bK2\nI0poWrMk1i1aZA5jXbbM7HKV4HJdU/B++ql5AxsHwVyz9MYb0LKlCiWxp25dMyW3bJk5Uy4313ZE\nIhLJH/4ADz+sQikZOI5pQP3JJ7Yj+ZXEKZYKCkx/jIcfth2JhF21ajBxovmZvPRS2LnTdkQiciRT\np8Ly5XDHHbYjEa8UTsUlmMQplt55B2rXhvPPtx2JiDmAd/x4OP5407dl82bbEYlIUXv3wr33wt//\nbg5jleSQng5z55pWAgkkMYol14W//tWMKqlbtySKsmXhlVfMsHDnzvDDD7YjEpFCL71k1reoP1py\nqVLF5NupU21HcojEKJYmTzYFU48etiMROZTjmMN377vP/AP+9lvbEYnIxo3mDfbzz+sNdjLq0cMc\nRZVAEmM3XHo63HIL3HCDZ7GIeG7CBLjzTnjrLbP7RrQbTuy46SazGeNvf7Mdifhh7VpzEPLGjVDO\n33aQwdkNN28erFwJ11xjOxKRyK68Et5+2+ySe+st29GIhNO0aeaYokcftR2J+KVJEzjuONP7LkHY\nL5aefdZMcZQvbzsSkeKdfz58/jk89BA880xCNk8TSVo5OdCvn1mvdMwxtqMRP11xBbz3nu0oDrA7\nDZeZaYbafvjBbNcWCYr1603zyrQ0+Mc/Qtu8UtNwEld/+hMsXqzjsMLg++/NEVRr1/q6Li0Y03DP\nP2/WKqlQkqBp1AhmzDD/oHv1gt27bUckktwWL4Z//hNeeMF2JBIPJ58MlSubNgIJwF6xtH07/Pvf\npk+GSBBVrw6TJpnOwV26wKZNtiMSSU55edC3Lzz5JDRsaDsaiQfHSaipOHvF0muvQffu0LixtRBE\nSi0lxRT9PXpAhw7m3a+IeOupp6BOHbjtNtuRSDwlULFkZ81SQQGceKI5C65jR89eX8SqMWPggQdg\n7Fjo1s12NHGhNUviu3nzTGPY//1Po0phs2+f2Rk3ZQqccoovL5HYa5Y++QRq1TLvxEWSRZ8+phdT\nnz4wfLjtaESCLzsbbrwRnntOhVIYlSkDV18N48bZjsTSyNIFF5iF3ddf79lriySMlSvNEQzdu5vW\nGD43VbNJI0viq9tvN5sn3nhDnbrDau5cuPZaWLHCl5+BxB1Z+u47WLIErroq7i8tEhfNm8OsWebn\nvEcP2LbNdkQiwTNuHGRkmFFaFUrh1b69ac3y3/9aDSP+xdJLL0H//jolWpJbaipMnGi2v3boAMuW\n2Y5IJDhWrjQ7pf/zHzWfDDvHgd69zVpQm2HEdRpu+3bTwnzJEqhf37PXFUloI0fC//2f2QGaZIdF\naxpOPLdrF3TqZN5U33mn7WgkEaxYAeeea5oBe7ysITGn4UaPNh05VShJmPzud/D+++aYhieeMDs8\nROTXXBduvhnOPBPuuMN2NJIoWrQwu+KmTbMWQvyKJdeFYcP0TkHCqVMn+OYbMzXXq5cZZRWRQ/31\nr7BunfldoXVKUpTlqbj4FUvTpplFWueeG7eXFEkoxx5rFqw2amTeOX/3ne2IRBLH22+b40zefRcq\nVLAdjSSa666DDz6w9kYzfsXSsGFw1116tyDhlpJiNjn86U+Qnm4aWYqE3Zdfmt8Pn3xi3lSIHK5e\nPdOO5Y03rLx8fBZ4b9gArVtDZqZ2NogUWrQIrrwSzjkHXnzRHBoZMFrgLaX2/ffmbMWxY6FrV9vR\nSCKbNs3skly40LOBl8Ra4P3aa6avkgolkYNatzYN13Jz4ayzTPEkEiYrVhxs3qpCSYqTlgZ795o+\ndnHmf7G0bx+8+qrpxCoih6pa9eCZcl26mAZ8Gt2QMFi92pzmMHgw3HCD7WgkCBzH7Cq2cJyU/9Nw\nU6bAoEHmEEQRObply0xb/yZNzBuMOnVsR1QsTcNJiaxebdbsPfSQWgRIbH7+GZo1g1WrzBmzpZQ4\n03AjRmhUSSQaJ50Es2dDy5Zw2mlmsatIslm40OyKVqEkJVGrFlx2WdxHl/wdWdq0yST+zEyoVs2z\n1xFJetOnw003mXfff/87VK9uO6Ij0siSxGTGDPjtb82Ghquvth2NBFXhpoDVq6FKlVLdKjFGll5/\nHa64QoWSSKzOP9+8Ay9fHk49FSZPth2RSOmMHm0aso4Zo0JJSueUU+C88+Bf/4rbS/o3suS65hDR\nUYJy+AMAAAhSSURBVKNM92IRKZmpU81U9jnnwD/+AbVr247oAI0sSbEKCsy61ffegw8/NL/oREpr\n/nxz1uaqVVCxYolvY39kqXBrX8eOvr2ESCh062baCtSrZ9oNjBql8+UkGH780ex4+/ZbmDNHhZJ4\np21baNfOtCaKA/+KpddeMwciqmO3SOlVqWJ60UycCK+8YhbIzp9vOyqRo/voI2jf3vRP+vRTT3Yu\niRzikUfgqadg927fX8qfabjdu6FxY/NuWK3rRbxV2Lvs0UfNrpC//MWMOlmgaTj5la1bYcAAs5h7\nzBjo3Nl2RJLMeveG444zhzCXgN1puHfeMdNvKpREvFemjFnDtGyZ6YrfqpVJFHF4dyVyVPv2mUXc\np55q1tUtWqRCSfw3dKhZ6L10qa8v40+x9NprcMstvtxaRPZLTTVTc7Nnm51zJ55oTm3fu9d2ZBI2\nM2bA2WebA9M/+ACee67UW7pFonLssfDHP8Ldd/t6+oH3xdIPP5h3FD17en5rETmC5s1h3DjzS+qj\nj8zj4cPNmXMifvr6a7jwQnNcyYABZmPPWWfZjkrC5p57YPNmM7LpE++LpTFjzJENKSme31pEIjjj\nDLMAfPx4s0X7+OPh6achK8t2ZJJM8vNhwgTTC+yGG+DKK2H5crN2RBt6xIZy5eCNN+DBB307Ws37\nBd7Nm8PYsXp3IWLbggVmPv+TT8wbmLvuMuubPKQF3iGybJl5M/z669C0Kdx7r2k6XL687chEjHfe\ngYED4Ztvot70Em0Oi6pYchzncuAUYB/wo+u6Y47yPNc98USz0ErvMEQSw4YNZgHkv/4FJ5wAffua\nDsoeHKESlGIpmhymYukwrmuOlXj/fXOtWwfXXw833mjOLhRJRI89Bp99BpMmRXV6iGfFkuM41YEv\nXNdtv//xLKCn67pbj/Bc133iCbPYSkIlIyODtLQ022FIJHl55tiUf//bJJP0dLjqKrj00hIfSRSE\nYinaHBb6Ysl1YeVKmDkTpk0zPyNly5r2FJdfbo6XKFfOdpS+UP5KIvv2mcXe06eb5QjNmkV8upet\nA84Dvi/yeAHQ5ajPvv76KG4pySYjI8N2CFKc8uXNxot33oE1a8wvwbFjoWFDSEsz65vmzDFFVXKJ\nLYeFwS+/mO/16NFmnUe3blCnjmkgOXGi2dn2xRfm5+SFF0xhnaSFEih/JZUyZcyuzLvvNkdEvfuu\nJyceRPPT3wgoukI0C2hx1Gcfd1zpIhIR/9WoATfdZK7duyEjw3RZ7tfPnLXUrp252rY1LQmaNTNr\nAII5vR5bDgsa14XsbPN93LEDtm83i/p//tlcmzbBxo1mOnbtWsjMNOe1nXQStGhh+iINGGC+1+qN\nJ8nijjvM8ToDBpgGvg8+CN27m5/xEuSxaIqlVCCnyOO9QNWYX0lEElOVKuZAyh49zOOsLJg71ywQ\n/+ILc7zKqlWwc6cZfahbF6pWhQoV7MYdvehzmBctT440lVf0gPGinzv82rfPXAUFh175+ebKyzN9\ntHJzIScH9uwxHytWNN/HatXMWrTUVHO8SK1apsg97TS46CJo0sQszq5ZM6iFr0j0zj8f5s0zU8ov\nvAB/+IP5uW/eHCpViukA3mjWLN0NHOe67oP7Hz8DbHJdd+gRnhviCX+RcArAmqWocpjyl0g4RZPD\nohlZWgWcUeRxbeCIjQwSPWmKSChFlcOUv0TkaKJZ4P0l0L7I49OBz/0JR0TEc8phIlIq0fZZ6gM0\nxRRXq1zXHet3YCIiXlEOE5HS8LSDt4iIiEiy8aRxRrQdviXYouyCvIqDW7V/77qufycbinWO47QF\nbihcPB1Eyl/hoRwmRcWSv0o9shRLh28Jrhi6IN8GfApscF03P/6RSrw4jjMQ6Axsd133ZtvxlITy\nV3goh0lRseavaBZ4F0fdccMh2u/zXtd11yrJJD/Xdf8OfGA7jlJS/goP5TA5INb85cU0XHJ3x5VC\n0X6fz3QcpwJQDVjuuu6H8QhOrAn6dnvlr/BQDpPDRZ2/vCiW1OE7HKL9Pn/uuu57AI7jzHcc50vX\ndbOO8DxJDkHfIaL8FR7KYXK4qPOXF9NwOzm0OqsE/OLBfSWxRPt9LjqsuQ1I8zEmsS/oI0vKX+Gh\nHCaHizp/eVEsrcJ0xC1UG/jRg/tKYin2++w4zg3AuCKfqgJo3j+5BX1kSfkrPJTD5HBxHVlSd9xw\nOOr32XGcZo7jOMAaYPj+z1UG6gBfxDdMibOgjywpf4WHcpgcLur85UlTSnXHDYejfZ8dx/kf8DvX\ndb91HOd6TIJpCoxzXXeOtYDFV/sPqL0aaAz8G3jOdd0dVoMqAeWv8FAOk0Kx5i918BYRERGJwItp\nOBEREZGkpWJJREREJAIVSyIiIiIRqFgSERERiUDFkoiIiEgEKpZEREREIlCxJCIiIhKBiiURERGR\nCFQsiYiIiESgYklEREQkAhVLUmKO49zoOM6TRR7f5jjOozZjEhGJlnKYREvFkpRGL2BekcdXAQss\nxSIiEivlMImKDtKVEnEcpyywCWjmuu52x3EqApuBhq7r7rQbnYhIZMphEguNLElJtQdWu667ff/j\nc4FFwG7HcWraC0tEJCrKYRI1FUtSUl2B1UUe9wa+BNKB+lYiEhGJnnKYRE3TcFIijuN8DhQAbwLl\ngRzgAuAb13WH2YxNRKQ4ymESCxVLErP9c/s/AvVd182zHY+ISCyUwyRWmoaTkugMzFeSEZGAUg6T\nmKhYkpI4FXjfdhAiIiWkHCYx0TSciIiISAQaWRIRERGJQMWSiIiISAQqlkREREQiULEkIiIiEoGK\nJREREZEIVCyJiIiIRKBiSURERCQCFUsiIiIiEfw/8RaZIdan3cIAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from scipy.stats import beta\n", "\n", "fig, axes = plt.subplots(2, 2,figsize=(10, 7))\n", "\n", "axes = axes.flatten()\n", "\n", "A = (0.1, 1, 2, 8)\n", "B = (0.1, 1, 3, 4)\n", "\n", "xx = np.linspace(0, 1, 100)\n", "\n", "for a, b, ax in zip(A, B, axes):\n", " yy = beta.pdf(xx, a, b)\n", " ax.plot(xx, yy, 'r')\n", " ax.set_ylim(0, 3)\n", " \n", " ax.set_xticks([0, 0.5, 1])\n", " ax.set_xticklabels([\"$0$\", \"$0.5$\", \"$1$\"], fontsize=\"large\")\n", " ax.set_yticks([0, 1, 2, 3])\n", " ax.set_yticklabels([\"$0$\", \"$1$\", \"$2$\", \"$3$\"], fontsize=\"large\")\n", " ax.set_xlabel(\"$\\mu$\", fontsize=\"x-large\")\n", " \n", " ax.text(0.1, 2.5, r\"$a={}$\".format(a), fontsize=\"x-large\")\n", " ax.text(0.1, 2.2, r\"$b={}$\".format(b), fontsize=\"x-large\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "有了先验分布，我们的后验分布为\n", "\n", "$$\n", "p(\\mu~|~m,l,a,b) \\propto \\mu^{m+a+1} (1-\\mu)^{l+b-1}\n", "$$\n", "\n", "其中 $l = N - m$。\n", "\n", "由共轭性，我们知道后验分布还是一个 `beta` 分布，从而有\n", "\n", "$$\n", "p(\\mu~|~m,l,a,b) \\sim \\text{Beta}(\\mu~|~a+m, b+l)\n", "$$\n", "\n", "至此，我们看到，如果我们观测到了一组 $m$ 个 $x = 1$ 和 $l$ 个 $x= 0$ 的数据，那么超参由 $a, b$ 变成 $a + m, b + l$。因此，超参 $a, b$ 可以看成是 $x=1$ 和 $x=0$ 的有效观测次数（注意：$a, b$ 可以不是整数）。\n", "\n", "更进一步，当有新的数据到来时，后验概率可以看成新的先验概率，因此这个过程可以序列化进行。下图表示观测到一个新的 $x=1$ 数据时，先验到后验的变化。" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAACkCAYAAABl7HEFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VNXWBvB3B0JCFUFK6KEJ5IIKAhYwQWmCBYQQkVCk\nFwFFUBGCIBIF0e9KUUBEvdKGdhEQRECaCNKRLiUYpEa8SDMMSfb3x0pCBpLMJFPOmZn39zx5IJnJ\nOTvJrDnr7LK20lqDiIiIiLIWYHQDiIiIiMyOCRMRERGRHUyYiIiIiOxgwkRERERkBxMmIiIiIjuY\nMBERERHZYTdhUkqVV0qtV0odVEodUEoN8kTDiMyKMUF0G+OB/IWyV4dJKVUaQGmt9V6lVCEAuwC0\n0Vof9kQDicyGMUF0G+OB/IXdHiat9Xmt9d7U/18DcBhAGXc3jMisGBNEtzEeyF/kaA6TUqoSgIcA\n/OKOxhB5G8YE0W2MB/JlDidMqV2tiwAMTr2LIPJrjAmi2xgP5OvyOvIkpVQggMUAZmutl97xGDej\nI1PRWit3n4MxQd7E3TGRXTykPs6YIFPJTUw4MulbAfgawCWt9WuZPK65ga/zRo8ejdGjRxvdDK+n\nlPLExYEx4QGMCddwd0zYi4fU5zAmXIAx4Rq5jQlHhuQeBxANoIlSak/qR8sct5DIdzAmiG5jPJBf\nsDskp7X+CSxwSZSOMUF0G+OB/AVf5CYRERFhdBOITIUxQWSLMWEsu3OY7B6AY9NkIp6Yw+RAGxgT\nZBqMCSJb7pzDREREROTXmDARERER2cGEiYiIiMgOJkxEREREdjBhIiIiIrKDCRMRERGRHV6dMEVE\nRKB3795GN4PI5bp164ZmzZoBkO0QqlWrlv7YnZ87e3x3ncMZd7aPyNuEhoYiNjbW6GbYd/kysGwZ\nMGoU0K4d8PDDQJUqQNmy8m+DBkCHDsDYscDatcDNm0a32DAObb5rVkuXLkXevF79IxBlSbboAoYN\nG4ZBgwZl+pgzx854jMzOYTRnf0YiR82ePRtdunRBSkqKy465c+dOFChQwGXHc6m//gIsFmD+fGD3\nbuDRR4GGDSUxqlQJKF4cCAqS5OjPP4G4OGDPHkmqDh4EWrcGevYEmjQB/ChOvTLbsFqtyJcvH4oW\nLeqyYxGZTVqhv4IFC6JgwYKZPubMsTMeI7NzGI2FDskbpV1Tihcv7vSxbt26hcDAQBe0KtXRo8CH\nHwKLFwMtWwJDhwLNmgHBwVl/T9WqwCOPAB07yucJCZJsDR4M5M0rSVSbNn6ROJliSC4iIgI9evTA\nW2+9hRIlSuCee+5Br169kJiYmP54z549ERMTg5CQEJQvXz7967169Uo/zq1bt/DWW2+hXLlyCAoK\nQlhYGObNm2dzroCAAEyePBkvvfQSihYtis6dO3vuByXKgbQeFnvDY3/99RcaNWqEJk2a4MqVKwCA\n+fPn48EHH0T+/PkRGhqK119/HTdu3MjyGFmdY9myZahRowYKFSqEJ598EidPnrR5fOXKlahXrx6C\ng4NRqlQpDBgw4K7zTJw4EZUrV0ZQUBCqVq2KTz755K72R0VFoVChQihdujRiYmKYLJENe9cIR977\nZ86ciZo1ayJ//vwoXrw4wsPDcebMGWzYsAFdunQBINeHgIAAdO/ePf37Jk+ejBo1aiB//vyoXr06\nYmNjkZycnP54pUqVEBMTg/79++O+++5D48aN078+bty49OddvXoVffr0QcmSJREcHIz69etjzZo1\n6Y+fOnUKAQEBmDt3Llq1aoVChQphxIgRrvkFxsUBnToBjRsDFSpI4jRvHvDss9knS5kpUQJ45RXg\n11+Bd98FRo+WnqYDB1zTVjNLu9PM7Yccwjnh4eG6SJEiunfv3vrIkSN6+fLlumTJknrQoEHpjxcu\nXFj369dPHz58WB84cEBrrXVERITu1atX+nGGDh2qixcvrhctWqSPHTumY2NjdUBAgF63bl36c5RS\nunjx4nrKlCn65MmT+tixY063n8wj9fXo9OvamQ9XxETXrl11s2bNtNZav/POO7pq1arpj2X8/Pff\nf9c1atTQkZGR2mq1aq21/vLLL/W9996rZ8+erePi4vSmTZt0nTp1dOfOnW2O37Rp00yPmfZ5wYIF\n9dNPP613796t9+3bpx966CEdHh6e/px9+/bpPHny6CFDhuijR4/qVatW6QoVKticZ8qUKTp//vz6\n888/18ePH9fTpk3TwcHB+osvvkh/Tps2bXS1atX0+vXr9cGDB3V0dLQuUqRI+s9PzvGFmLB3jbD3\n3r9z506dN29e/c033+j4+Hi9f/9+/cUXX+g//vhDW61WPXXqVK2U0hcuXNAXLlzQV65c0VpLHFSs\nWFEvXbpUnzp1Sq9cuVJXqFBBx8TEpLetYsWKukiRInrMmDH62LFj+vDhw1prrStVqqTHjRuX/rz2\n7dvr0NBQ/cMPP+gjR47owYMH63z58ukjR45orbWOi4vTSildrlw5PWfOHB0XF6fj4uKc+r3p69e1\nHj5c62LFtB4zRuvUn8ulbt3SeupUre+7T+vx47VOSnL9OVwstzFheCBoLcEQGhqqU1JS0r82Y8YM\nHRQUpK9fv67Dw8P1/ffff9f3ZUyYrl+/roOCgvRnn31m85y2bdvqJ598Mv1zpZTu2bOn020mc/KF\ni4PWjiVM+/bt02XKlNGvvPKKzfdWrFhRT58+3eZrGzdu1Eopffny5fTj20uY8ubNq//888/0r82f\nP18HBATomzdvaq21jo6O1g0bNrQ5z7fffqsDAgJ0fHy81lrrcuXK6TfffNPmOa+99pquXLmy1lrr\nY8eOaaWUXrt2bfrjVqtVly1blgmTi/hCTNi7RuTLly/b9/4lS5boe+65Jz0RutM333yjU/e7S3f9\n+nVdoEABvXr1apuvf/3117po0aLpn1esWNEmltJkTJjSXuerVq2yeU7dunV19+7dtda3E6b33nsv\n29+Fw9at0zo0VOuOHbU+c8Y1x8xOXJzWjRtr3aKF1pcuuf98TshtTJhiSA4AGjRoYDPJ87HHHoPV\nasWJEyeglEK9evWy/f7jx4/DarXiiSeesPn6E088gYMHD951LiJvlpCQgPDwcHTq1AmTJ0+2+Xp8\nfDxee+01FC5cOP2jVatWUErh+PHjDp+jTJkyNvMwypQpA601Ll68CAA4dOhQpvGmtcahQ4dw5coV\nnDlzJtPnnDp1ComJiTh06BAAifc0gYGBqF+/vuO/DPIL2V0jbt26le17f/PmzVG5cmWEhoaiY8eO\n+Pzzz3Hp0qVsz3fw4EH8888/eOGFF2xiqW/fvrhy5Ur69yul7F5T0l7nHrk+JSYCAwcCXbsCU6cC\nc+cCZco4d0xHVKoE/PgjUKuWrKz77Tf3n9PDTDPpW5K+rLlyQqrZJrcS5VTRokXxwAMPYOnSpRg8\neDDKli0LAOmrfCZNmoQmTZrc9X1pz3PEnYsh0i5WGVcS2Yvb3HLXccl7OfOaKFiwIHbu3IktW7Zg\n7dq1mDZtGt544w2sW7cOdevWzfR70l7nixYtQvXq1e96/N5777U5fm5k9jM5dX06ehSIjARq1JA5\nRhna6BF58wIffyxJU3i4lCvwoZsf0/Qw7dixw+aN+Oeff0ZQUBCqVKniUKBUrVoVQUFB2Lhxo83X\nN27ciNq1a7u8vURGypcvH5YsWYLatWsjPDwc8fHxAIBSpUqhfPnyOHLkCCpXrnzXR1BQkMvaEBYW\nhk2bNtl8bePGjVBKISwsDEWKFEG5cuUyjcnKlSsjODgYtWrVAgBs2bIl/XGr1YodO3a4rJ3kG7K7\nRjjy3h8QEIDGjRtjzJgx2LVrF0JCQtInhqfdHGS81oSFhSE4OBgnTpzINJYCAhy/fIaFhaW3KaNN\nmza57vq0ZAnQqJFMyLZYPJ8sZdSzJzB9upQf+OUX49rhYqbpYbp06RIGDBiAwYMH48SJExg1ahR6\n9eqVXscis6RJ3x4fR4ECBTBo0CDExMSgRIkSqFOnDhYtWoRly5Zh7dq1Hv1ZiFwhuxsFrTXy5MmD\nBQsWoFOnTggPD8ePP/6I0NBQjBs3Dj169EDRokXx/PPPIzAwEIcPH8b333+PadOmuax9w4YNQ926\ndTFkyBD07t0bp06dwsCBAxEdHY1y5coBAIYPH47XX38d1apVS2/jtGnT8OmnnwKQG53nnnsOAwYM\nwPTp01GyZEl88MEHuHbtmsvaSb4hu2uEvff+b7/9FnFxcWjcuDFKlCiBXbt24fTp0+kJe2hoaPrz\nHn/8cRQoUACFChXC22+/jbfffhsA0LRpUyQlJWH//v3Yu3cvPvjgAwBZx2nGr1epUgWRkZHo378/\npk+fjgoVKuCzzz7DoUOHMH/+fOd+MVpLUcmZM4GVK83To/Pcc8CXX8pKvFWrADvTaryBKRImpRQi\nIyNRuHBhNGrUCFarFR06dMCECRPSH8+siN2dXx83bhwCAgLw6quvIiEhAdWqVcOcOXMyHZogMrOM\nr+07X+cZP8+TJw/mzp2Lrl27Ijw8HOvXr0d0dDQKFy6M8ePHIzY2Fnnz5kXlypXRrl27TI/hyOcZ\nv56mdu3aWLZsGWJiYvDpp5+iSJEiiIyMxMSJE9Of069fP1y/fh2xsbHo378/KlSogPHjx+Pll19O\nf86sWbPQr18/PPPMMyhQoAB69eqFtm3b4uzZs878CsmH2LtG2HvvL1asGCZNmoTY2FhcvXoVFSpU\nQExMTPrrsH79+hg8eDD69u2LhIQEdO3aFbNmzcLIkSMREhKCKVOmYOjQocifPz/uv/9+dOvWzaZt\nWbU5o5kzZ2LYsGGIjo7GlStXUKdOHaxYscJmuC/HxVpv3gS6dwdOnAC2bwdKl87Z97tb69bAjBmS\nNG3eLJXDvZhydq5A6soCp47RpEkTVKtWDTNmzHDqOERKKWitDa2g5oqYIHIVX4gJXiMy8fffQNu2\nQNGiwJw5QP78Rrcoa9OmAR99JMNzxYoZ3Zpcx4Qp5jBlHFojIiLKiNeIOyQkSLHImjWBhQvNnSwB\nQN++MkT34otAUpLRrck1UyRMWXX/ExER8RqRwdmzUrG7dWtgyhQgTx6jW+SY8eOBgABg+HCjW5Jr\nphiSI3IVXxh+IHIlxoQPOX0aePJJoEcP4K23jG5Nzl26BNStK/WhnnnGsGbkNiaYMHlacjJw5ox8\nJCQAV68C//wjj+XNCxQoIGO8pUrJnj/33GNse70MLw5EthgTPuLMGalt1K8f8PrrRrcm937+WeZe\n7doFpK6m9TQmTGZktQI7d8oLZOdOYP9+Wc1QvLi8UEqUkIQoOFh2er51C7hxQ7Lw8+eB+HgZm65V\nS7Ly+vWBJ54w7EXmDXhxILLFmPAB589LstSzJzBsmNGtcd577wGbNgGrV8u1z8OYMJnFxYvA0qXA\nihXAhg1A1apSTOzhh4EHHgDuv9/x3aG1lkA5cADYvRvYtk2WZhYvDrRqJVn64497zxi2B/DiQGSL\nMeHl/vpLkqWoKGDkSKNb4xpJSXLt6toV6N/f46dnwmSkmzeB//4X+OorSWpatgSefx5o3lySG1dK\nSQH27pWEbMkSSdA6dZJaHDVruvZcXogXByJbjAkvdu0a8NRTkjCNH29Ib4zbHD0qSdPOnbIPnQcx\nYTLC2bOySmHmTOk96t5dEqXU6uQecfgw8PXX8lGjBvDqq1IkLAdl+30JLw5EthgTXspqlffy8uWB\nzz/3rWQpTWws8NNPwHffefTnY8LkSadOyR964UIgOhoYNAioVs3YNlmtwOLFUhzs2jVgxAigY0eZ\nSO5HeHEgssWY8EIpKUDnzsD168CiRb77Pm61ypYpI0ZIjSYP8erClV7jwgXZ2LBePZmwfewYMHmy\n8ckSAOTLJwnSjh3S6zVjBlC7tgzb8Y2KiMh7vP02EBcHzJvnu8kSINet6dOBoUPlRt/kmDA54uZN\nYMIEICwMCAwEjhwBxo0D7rvP6JbdTSmgaVNZgfDxx8CYMTL+vWeP0S0jIiJ7pk2TObHLlpm/grcr\nPPaY1JYaN87oltjFITl7Nm6Usu5VqgD/93/m6E3KieRk4IsvgJgY6YEaOxYoXNjoVrkNhx+IbDEm\nvMjq1bJybMsWr9+oNkfOnZMRkW3bZGW5m3FIztWuXJFEKToaeP99YPly70uWACk50Ls3cPCg/Ez/\n+pcEJRERmcfBgzJvadEi/0qWACAkBBgyRIYiTYwJU2Y2b5ZVb0lJUgOpTRvvX6Fw333ArFmyoq9P\nH6kWe/260a0iIqKEBFkR99FHUrfPH736qhR5/uUXo1uSJSZMGSUlAaNGAR06yGTumTN9b2uSZs2A\nfftkgl29evJ/IiIyhtUKtG8vhSk7dza6NcYpUAB4912pZG7S4VsmTGnOn5fJ0lu3ygRpAzcGdLt7\n7gG++UaqxjZtKjU+TPoCJSLyaYMGyXuyF0x6druuXaUY85o1RrckU0yYAOkGfPhhoEkTmd9TurTR\nLfKM6GgpGvbvf8seRYmJRreIiMh/TJ8uK5pnz/bbYsM28uSRld0xMaa8iedfaNYsmaM0fTrwzjv+\n96K9/34ZM75yRRLG8+eNbhERke/bskWmgHz7LVCkiNGtMY/ISOCff6T6t8n4WXaQQUoK8MYbwAcf\nSIbfurXRLTJOoULAggXA008DDRsC+/cb3SIiIt919qzMlf3yS+9cfe1OAQHSeTF2rOl6mfwzYUpM\nlDLsW7fKR40aRrfIeErJ3c748bLZ49q1RreIiMj3WK3Si9KvH9CqldGtMae2bWXUY906o1tiw/8S\npr//lp4UQCaWFS9ubHvM5sUXZU+6Tp2A+fONbg0RkW8ZOlSuOyavOWSogABg+HDTTYT3r4QpIUHm\n6dSqJXv0BAcb3SJzatxYepiGDpUy/URE5Lx584CVK4H//Mf/5svmVMeOstH91q1GtySd//zFzp6V\nPdVatZLNafPkMbpF5la7tsztmjABmDjR6NYQEXm3gwelhMDixUDRoka3xvwCA4HXXpNinibhH3vJ\nnTkjPUvdurEbNKf++EM2Rnz5ZekiNTnum0VkizFhAlevAvXrA2+9Jdchcsy1a0ClSrKS24XbxeQ2\nJnw/YTp7FoiIAHr0AN580+jWeKezZyVp6tZNAt7EeHEgssWYMJjWwEsvAQULyu4RlDPDh0vCOWWK\nyw7JhCkzFy/KMFyXLl7RO2JqaUOa/ftLN6lJ8eJAZIsxYbDPPpO5oNu2AfnzG90a73P2rGwaHxfn\nsq3KchsTvjuH6fJloHlzqXXBZMl5ZcoAP/4ITJokW6kQEVH2du2Sci0LFzJZyq0yZYAWLaRmlcHs\nJkxKqVlKqQtKKe+pZnjjhhSibNIEGD3a6Nb4jvLlpRTD6NFS6NJPeWVMELkRYyITf/8tN+xTpgDV\nqxvdGu82cCAwdaoUnDaQIz1MXwJo6e6GuExSkrxIq1aV2fXK0J5o31O1KrBqFfDKK9Lj5J+8KyaI\n3I8xkZHWQK9e0jMSFWV0a7zfo4/K9jHff29oM+wmTFrrzQD+54G2OE9roG9fIDlZJtexzoV71Kkj\nPUwvvgjs3Wt0azzOq2KCyAMYE3eYNg04dgz4+GOjW+IblAIGDJD5YAbyrYxi3Dhgzx4ZLw4MNLo1\nvi0iQrqan31WSg8QEZHcRI4aBVgsLI7sSlFRwM8/A6dPG9YE30mY5s6VycgrVshmsuR+HTrI2HLr\n1rLsk4hMIykJ+OEHo1vhZ65elffFf/+b85ZcrWBBqf5tYGmGvK44yOgME6sjIiIQERHhisM6butW\nYPBgmVMTEuLZc/u7YcOA48eB6GhgyRKPV1DfsGEDNmzY4NFzOsLwmCC/lJwMTJq0AXPmbMDhw+Yq\nKO3zMaG1lF154gnZi5Ncr08f2Qs2JgbI63j64qrrhEN1mJRSlQAs11rXzuQxY+trnD4NNGwovUut\nWxvXDn9mtcrkxgYNgPHjDW2Kp2rOmDomyK+kpMhIhcUCLFok94wdOsgIRmgoY8JjvvoK+PBDYMcO\noEABo1vjux59FBgxAnjmmVwfwm11mJRS8wD8DKC6Uuq0Uurl3DTQLW7cANq0kUKKTJaMky+fvFMv\nWiRDoz7O1DFBfkFrYPt2YMgQoGJFoF8/oFQpYONGYPduKcgfGuq59vh9TBw+LL3tFguTJXd7+WXD\najJ5b6VvrYHOneX/33zD8gFmsH+/bKGyejVQt64hTWBVY/JVWsuaFotFFqnmyye9SFFRQFhY1t/H\nmHCzf/6RUY6BA6WUALnX33/LXcKxY0CJErk6hP9tjTJ5MvDFF9IXzYzePBYuBN54A9i5Eyhe3OOn\n58WBfM2BA8D8+ZIkJSXdTpIeeMCx+0TGhJv17w9cuiR/JN64e0aXLnJT/uqrufp2/0qYtmwBXnhB\nJntXruzZc5N9Q4fKu/zKlR6vhcWLA/mCo0elJ8liub3wKioKePjhnF+TGRNutHix3CDu3u2yfc7I\nAevXS7K0b1+uvt1/EqaEBMksp03jvCWzSkqSobnmzYGRIz16al4cyFudPCkJ0vz58jYXGSlJ0iOP\nOHffwZhwk1OnZKHLihXyL3lOSgpQqZL87uvUyfG3+0fClJwMtGoF1KsHxMZ65pyUO2fPyu3w7NmS\nPHkILw7kTeLjZajNYpH/t2snSVKjRq6r0MGYcINbt6R8QLt20qNOnjd8uCROuViZ7R8JU2ys7CXz\n4485qsFABlm7FujaVbqrS5XyyCl5cSCzO3NGpvpZLMBvv8nsgqgoKZ7vjrc1xoQbDB8uw0ErVnAL\nLqMcPCjlbH7/Pcd3F76fMG3ZItn8zp1AuXLuPx+5xogR8jdbtcojbyy8OJAZXbwoVTcsFllM+uyz\nshVj06bu38WJMeFiP/wAdO8uSxZzuUqLXOShh4CJE4GnnsrRt7mtDpMpXL4slVNnzGCy5G3GjAGu\nXQM++sjolhB51KVLsotDs2ayS8aWLcDrrwPnzgFffy0Fi7nlpZc5fx7o1k1K2TBZMl6nTsC8eR47\nnfl7mLQGXnpJlqhPmeK+85D7nDoF1K8vd2YPPeTWU/Fumox0+TKwdKn0JP38s4wYdOggUy+Nqn7C\nmHCR5GT5gz7+uNwIkvHi4+Wacu6cFCZzUG5jwvwTgWbPBn79VYZ1yDtVqgR88okkvrt2sW4W+ZSr\nV4HlyyVJWr9e1jh06SLzlLgPuA/54ANZATxqlNEtoTQVKgA1awJr1nhk1by5e5h+/11WWq1ZAzz4\noHvOQZ7z0kvAffcBkya57RS8myZPuHFD5vsuWCBvT40aycTt5583XzkexoQLbN4sXYU7dwJlyxrd\nGspoyhTgl19kmNRBvjfpOyVFJnK1aCEbI5H3+9//pDzxzJlSo8kNeHEgd0lMlEW6FousYWjQQCZu\nt2kDFCtmdOuyxphw0p9/yrDPjBky8YzM5fx5oEYNGZbLn9+hb/G9Sd+TJwNWq2xoSL7h3ntlO5se\nPWSyB5HJWa3Ad9/JEFtIiIwsh4fLNlZpi6XMnCyRk1JSpDRKx45MlsyqdGkpZr16tdtPZc4ept9+\nAx57TLY+qVbNtccm4/XrJ7frbthxmnfT5KykJCn1ZrHIBO6aNWW4rX17SZq8DWPCCRMmAN9+C2zY\nwCWNZvbpp7LKYvZsh57uO0NyyclSQTUqChg0yHXHJfO4dk3K2U+e7PKJerw4UG4kJwObNkmStGQJ\nEBoqb0GRkUD58ka3zjmMiVxKq/23Y4f3vwh83blzQK1aMjwXFGT36b6zSm7KFKna+corRreE3KVQ\nIZnH1LWrbNJrtlmy5BdSUuSm1GKRopIhIZIk/fKLJEzkxxISZILarFlMlrxBSAjwr3/J7hJuXC1n\nrh6mkydlJiWH4vxD375ya//55y47JO+mKTtaA9u3S5K0cCFQtKgsfoqKkuKSvogxkUNpe5bWrQu8\n/77RrSFHffIJsHevQ1M9vH9ITmspiduiBSd6+4srV4CwMOA//wGaNHHJIXlxoDtpLe+jFot85Msn\nnQcdOsjLz9cxJnLo3XeBdevkg3uWeo/4eElyz5+3+3fz/iG5r7+WZeevvWZ0S8hTihQBpk4FeveW\n4qQOLgklskdrGe1NS5KSk6UX6b//lcoWytD0gUxrzRpg+nSpt8RkybtUqABUrChzz8LD3XIKc/Qw\nXbwI1K4txU3q1nXuWOR9oqKAKlWA2FinD8W7af929KgkSPPny9qCDh2kN6lePf9NkhgTDoqPlykh\n8+cDERFGt4ZyY+xY4K+/gP/7v2yf5t1Dcp07A6VKya7D5H/On5eEef16mbjnBF4c/M+JE5IkLVgg\n916RkZKDP/IIEGDeSnMew5hwwM2bsjq7fXtOCfFmv/4q5fZPnsz2Dsl7E6Z166T628GD3HjJn332\nGTBnjqztduIqx4uDf4iPlwTJYpH/t2snSVKjRrLIlm5jTDigTx+p6L1okf92RfoCrWW0Im3sPQve\nWen75k2gf3+px8Nkyb/16SMVA91QzJJ8w5kzshDmscdk5P7oUVnEdOaM1K0LD2eyRLkwc6bcqH35\nJZMlb6eU9DB9+617Dm9oD1NsrBQ9cdMPR15m927ZfuDQIaB48VwdgnfTvuXiRbnpt1iA/fuB556T\neUnNmrHwsqMYE9nYtg149llJmGrWNLo15Ao//ggMHy65RRa8b0ju1Cng4YdlNUKlSk61gXzIwIHS\n8zhjRq6+nRcH73fpklTbtljk7aF1axlua9HCoSK+dAfGRBbOnZNJ3lOnSiZOvsFqlTnRR47Iv5nw\nvoTphRekX33kSKfOTz7m8mW501u2DKhfP8ffzouDd7p8WfZts1ik+nbz5rK6rVUrVptwFmMiEzdv\nykq4p58GRo0yujXkah06AC1byvzoTHhXwrR6NTBggBRKCQ526vzkg776SiaBb92a4wngvDh4j6tX\nJS+2WICNG4Enn5T3uWef5ZRGV2JM3EFruZBeuyYvPi6l9D3/+Y/cgS1ZkunD3pMwWa2y8erEicAz\nzzh1bvJRKSnA448DvXpleYeQFV4czO36dWDlSrlOrVkjq9pefFHmaRYpYnTrfBNj4g4TJ8qK3J9+\nAgoWNLpzjRf4AAAMM0lEQVQ15A4JCUDVqjIJMpNxfO+p9D1liuxsyWSJshIQICsnn31W6qLwSurV\nEhOlJq3FAnz/PdCwofQkzZgBFCtmdOvIryxfLkUNt21jsuTLSpQAatUCNm8GmjZ12WE928N08aJs\n3rR5M1CjhlPnJT/w8svywp8wweFv4d20OVit0oNkscg16sEHZeJ2u3byJyXPYUyk2rNHVg4sXy5Z\nO/m2MWNk3D+TgtjeMSTXu7dMTvj4Y6fOSX7i3Dmp/L1tG1CtmkPfwouDcZKSZEWvxSLTB2rWlCSp\nfXsgJMTo1vkvxgSAP/4AHn0U+Pe/JWsn37d9u9x0Hzx410PmT5h+/VWKpxw5Atx7r1PnJD/y/vvA\njh1ZTt67Ey8OnpWcLCVsLBb5E4WGSpIUGQmUL2906whgTODyZaBxY6BrV2DoUGPaQJ6XkiJlBXbt\nko15MzB3wqS1rBNu00ZWxxE5KjFRhm+/+sqhDTH9/uLgASkpsvTfYpGikiEhkiR16CAJE5mLX8fE\nzZtSOuBf/5Iy8azk7V+io2WPwN69bb5s7q1RVq6ULtE7Gk1kV3Aw8MEHwJAhcqUmQ2gthXOHDAEq\nVgT69QNKlpRyALt3A2++yWSJTCY5WS6YxYrJRG8mS/6nZUtZceIi7u9hSkqSMgITJnBlHOWO1jL/\nYMAAoHPnbJ/q13fTLqY1sHev9CRZLLI6N60nKSzM6NaRo/wyJrSWfUqPHpULJkvE+6eLF4H775cy\nA3lvFwUwb1mBWbNkHLF1a7efinyUUrLS4aWXZAYxSz+7jdZST9ZiARYskPudqKjbm3/zJp1MT2vp\n8ty5E1i3jsmSPytZUrZe275ddu12knt7mK5dA6pXl3K+Dz/s1HmI8MILshz4zTezfIpf3k27wJEj\nt3uSrl2TXqSoKAlbJkneze9iYswYYOFCGS/O5Sbe5EPefFNuskePTv+SOSd9jx0LHD4MzJ3r1DmI\nAEj3eqNG8m8WFQ/97uLghBMnpBfJYpGe68hISZIeeYS7RfgSv4qJ2FjZFmPDBqB0afefj8xv7Vrg\nnXeALVvSv2S+hCkhQQqxbN8OVK7s1DmI0vXtCxQuDHz4YaYP+9XFIRfi428nSfHxUpImKkry0Dx5\njG4duYPfxMTYscDs2cD69UCZMu49F3mPxESplnv6NFC0KAAzJkyDB8tY8qRJTh2fyEZaMcs9e+6q\nrQH40cUhB86elREKiwX47TegbVtJkiIibOZBko/y+ZjQGhg5UqqlrlvHniW6W8uWQJ8+8uYHsyVM\np04B9eoBhw7JhG8iV3r7bRlDmjnzrod8/uLgoIsXpUaSxSI1Y59/XpKkpk2BwEBDm0Ye5tMxkZIC\nDBwIbN0KrF7NfXcocx9+CPz+u+xlC7MlTN26SZnfsWOdOjZRpv73P1lMkMmehD59cbDj0iWptm2x\nyAKh1q0lSWrRgguF/JnPxkRiItClC3Dhgiwsuuce1x6ffMfu3UCnTjKnGmZKmA4dkr7+Y8f4Aib3\nGT9esoKFC22+7LMXhyxcviwjERaLVN9u0UKSpFatWH2BhE/GxJ9/yqrZ0qVlkndwsOuOTb4nJUV6\nH3/9FShb1kSVvt95R/brYbJE7jRwoKx62LPH6JZ43NWrwJw5wHPPSdXtpUvlRvvMGZnQ3a4dkyXy\nYYcOSXmRxx4D5s9nskT2BQQATZrI7uBOcG0P0549Mg5w/DhQoIBTxyWya/JkmbewYkX6l3zybhrA\njRvAd99JT9KaNbKqLSpK5ibx3oSy41MxsWCBVPz/6CO5SyBy1Kefyqr9r74yyZDcM8/ImMDAgU4d\nk8ghN28C1apJFvHoowB86+KQmAh8/738eKtWAQ0aSJLUtm2WZaiI7uITMZGYKCMX330HLF4M1K3r\nusaRfzh6FGjWDPj9d6iAAIO3Rtm2TcYHFy922SGJshUUBMTEAKNGSbeLD7Ba5UexWIDly4EHH5Qk\nadIkLgAiP7V3r/Qm1aghoxiptXSIcqR6dZnLdOJErg/hujlMo0fLcm8uxyFP6tZNhoB/+snoluTa\nrVvADz8APXoAISHA++8D9evLVI3166VWJ5Ml8jv//CP1lZo3l94li4XJEuWeUrIgbePGXB/CbsKk\nlGqplDqilDqmlMp8E6+ff5bNqLp3z3VDiHIlMFDeVN95x2OndCgm7EhOvp0MlS0rHWVhYXIz/dNP\nMqodEuLqlhO5hytiIp3WUkQsLExWW+/ZIz1M3NSQnBURIdvm5FK2CZNSKg+AKQBaAqgFoKNSquZd\nTxw9GhgxAsiXL9cN8XcbnPgj+r0uXaRY6ubNbj+VwzGRiZSU28lQuXLA668DoaHAL7/Ix5AhUr6M\nBGPCOzgTEza0lgUcjzwCjBsHfP659CqVLeviFnsvxoSTIiLkTjWX7PUwNQBwXGt9Smt9C8B8AM/f\n9azffgO6ds11I4iB4JTAQEnYx4zxxNkci4lUWt9OhipUAPr1A0qWlF7h3btlI+3QUE802/swJrxG\njmLiLpcuAVOnypZHI0fKncSuXcBTT7mrvV6LMeGkKlWc+nZ7k77LAjid4fM/ADS861nDh7N3iYzV\nuTPw3nueOJNDMbF7t9wcL1gg0/qioqQCQliYJ5pI5FGOXSfS3LwJHDgAbNokQbF1K/D005I0hYdz\n6I3cJ20e05w5ufp2ewmTY+tAu3XL1cmJXCYwUBYd9Orl7jM5FBPt20uS9N//Ag88wGsA+TTHrhON\nGslO0OfOAVWrSimQnj2lWn/hwm5uIlEqJxKmbOswKaUeATBaa90y9fPhAFK01uMzPMc827ITAW6t\nOcOYIG/EmCCy5fLClUqpvACOAngKwFkA2wF01Fofzm0jibwZY4LIFmOC/EW2Q3Ja6ySl1CsAVgPI\nA+ALBgH5M8YEkS3GBPkLp7dGISIiIvJ1Dlf6dqQwmVJqUurj+5RSD7mumb7D3u9RKRWhlPpbKbUn\n9WOkEe00M6XULKXUBaXU/mye4/bXImPCNRgTzmNM+BbGhPPcEhNaa7sfkG7W4wAqAQgEsBdAzTue\n0wrAytT/NwSwzZFj+9OHg7/HCADLjG6rmT8ANAbwEID9WTzu9tciY8Kjv0fGhP3fI2PCRz4YEy77\nPbo8JhztYXKkMNlzAL4GAK31LwCKKqVKOXh8f+FogTcuQs+G1nozgP9l8xRPvBYZE67BmHABxoRP\nYUy4gDtiwtGEKbPCZHfWq8/sOeUcPL6/cOT3qAE8ltpFuFIpVctjrfMdnngtMiZcgzHhGYwJ78GY\n8IwcvxbtFa5M4+jM8DszXs4ot+XI72M3gPJa6xtKqacBLAVQ3b3N8knufi0yJlyDMeE5jAnvwJjw\nnBy9Fh3tYToDIOO2oOUh2Vh2zymX+jW6ze7vUWt9VWt9I/X/qwAEKqWKea6JPsETr0XGhGswJjyD\nMeE9GBOekePXoqMJ004A1ZRSlZRS+QBEAVh2x3OWAegCpFd+vay1vuDg8f2F3d+jUqqUUrKRhlKq\nAaT0w1+eb6pX88RrkTHhGowJz2BMeA/GhGfk+LXo0JCczqIwmVKqT+rj07XWK5VSrZRSxwFcB/Cy\nEz+IT3Lk9wigPYB+SqkkADcAvGhYg01KKTUPQDiA+5RSpwG8A1lN4rHXImPCNRgTrsGY8B2MCddw\nR0ywcCURERGRHQ4XriQiIiLyV0yYiIiIiOxgwkRERERkBxMmIiIiIjuYMBERERHZwYSJiIiIyA4m\nTERERER2MGEiIiIisoMJExEREZEdTJgMppTqopQal+HzXkqpGCPbRGQkxgSRLcaEOTBhMt4LAHZl\n+DwSwD6D2kJkBowJIluMCRPgXnIGUkrlAXABQBWt9d9KqWAAFwGU1VpfNbZ1RJ7HmCCyxZgwD/Yw\nGasegJNa679TP28M4ACA60qpYsY1i8gwjAkiW4wJk2DCZKymAE5m+PwlAJsAPAmgtCEtIjIWY4LI\nFmPCJDgkZyCl1DoAyQDmAggEkAjgKQDbtdafGtk2IiMwJohsMSbMgwmTQVLHoc8AKK21vmV0e4iM\nxpggssWYMBcOyRmnEYC9DAKidIwJIluMCRNhwmSc2gCWGt0IIhNhTBDZYkyYCIfkiIiIiOxgDxMR\nERGRHUyYiIiIiOxgwkRERERkBxMmIiIiIjuYMBERERHZwYSJiIiIyA4mTERERER2MGEiIiIisuP/\nAZ0lT8mMg+24AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "xx = np.linspace(0, 1, 100)\n", "\n", "fig, axes = plt.subplots(1, 3, figsize=(10, 2))\n", "\n", "axes = axes.flatten()\n", "\n", "axes[0].plot(xx, beta.pdf(xx, 2, 2), 'r')\n", "axes[0].set_ylim(0, 2)\n", "axes[0].text(0.1, 1.6, \"prior\", fontsize=\"x-large\")\n", "axes[0].set_xlabel(\"$\\mu$\", fontsize=\"x-large\")\n", "axes[0].set_xticks([0, 0.5, 1])\n", "axes[0].set_yticks([0, 1, 2])\n", "\n", "axes[1].plot(xx, xx, 'b')\n", "axes[1].set_ylim(0, 2)\n", "axes[1].text(0.1, 1.6, \"likelihood\", fontsize=\"x-large\")\n", "axes[1].set_xlabel(\"$\\mu$\", fontsize=\"x-large\")\n", "axes[1].set_xticks([0, 0.5, 1])\n", "axes[1].set_yticks([0, 1, 2])\n", "\n", "axes[2].plot(xx, beta.pdf(xx, 3, 2), 'r')\n", "axes[2].set_ylim(0, 2)\n", "axes[2].text(0.1, 1.6, \"posterior\", fontsize=\"x-large\")\n", "axes[2].set_xlabel(\"$\\mu$\", fontsize=\"x-large\")\n", "axes[2].set_xticks([0, 0.5, 1])\n", "axes[2].set_yticks([0, 1, 2])\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "如果我们的目的是预测下一次实验的结果，那么我们有\n", "\n", "$$\n", "p(x=1|\\mathcal D) = \\int_{0}^1 p(x=1|\\mu) p(\\mu|\\mathcal D) d\\mu\n", "\\int_0^1 \\mu p(\\mu|\\mathcal D) = \\mathbb E[\\mu|\\mathcal D]\n", "$$\n", "\n", "而我们知道后验概率的分布，所以可以计算出其均值：\n", "\n", "$$\n", "p(x=1|\\mathcal D) = \\frac{m+a}{m+a+l+b} = \\frac{m+a}{N+a+b}\n", "$$\n", "\n", "当 $m, l \\to \\infty$ 时，有\n", "\n", "$$\n", "p(x=1|\\mathcal D) = \\frac{m+a}{m+a+l+b} = \\frac{m+a}{N+a+b} \\to \\frac{m}{N}\n", "$$\n", "\n", "即趋向于最大似然的解。当 $N$ 有限时，我们的解总在先验均值和最大似然解之间。\n", "\n", "另外我们观察到，随着 $(a,b)$ 的增加，函数分布变得越来越尖，即方差越来越小。事实上，随着观测数据的增加，我们对于后验分布的不确定性也在不断减小。\n", "\n", "另外考虑条件期望和方差的性质，我们有：\n", "\n", "$$\n", "\\mathbb E_\\mathbf\\theta[\\mathbf\\theta] = \\mathbb E_\\mathcal{D}[\\mathbb E_\\mathbf\\theta[\\mathbf\\theta|\\mathcal D]]\n", "$$\n", "\n", "和\n", "\n", "$$\n", "\\text{var}_\\mathbf\\theta[\\mathbf\\theta] = \\mathbb E_\\mathcal{D}[\\text{var}_\\mathbf\\theta[\\mathbf\\theta|\\mathcal D]] + \\text{var}_\\mathcal{D}[\\mathbb E_\\mathbf\\theta[\\mathbf\\theta|\\mathcal D]]\n", "$$\n", "\n", "从而我们知道，从平均意义上来说，后验分布的方差要比先验分布的方差要小。后验分布在数据分布上的均值等于先验分布的均值" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.6" } }, "nbformat": 4, "nbformat_minor": 0 }