{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 1.2 概率论" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "加法率：\n", "\n", "$$\n", "P(X) = \\sum_Y P(X,Y)\n", "$$\n", "\n", "乘法率：\n", "\n", "$$\n", "P(X,Y)=P(Y|X)P(X)=P(Y,X)=P(X|Y)P(Y)\n", "$$\n", "\n", "由乘法率可得 `Bayes` 公式：\n", "\n", "$$\n", "P(Y|X)=\\frac{P(X|Y)P(Y)}{P(X)}\n", "$$\n", "\n", "再由加法率和乘法率，`Bayes` 公式的分母可以写成：\n", "\n", "$$\n", "P(X) = \\sum_Y P(X,Y) = \\sum_Y P(X|Y)P(Y)\n", "$$\n", "\n", "所以 `Bayes` 公式的分母可以看成是一个归一化项，使得条件概率 $P(Y|X)$ 对所有的 $Y$ 求和之后和为 `1`。\n", "\n", "如果\n", "$P(X,Y)=P(X)P(Y)\n", "$，那么 $X, Y$ 是独立的。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.2.1 概率密度函数" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 概率密度函数" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "对于实值变量 $x$，概率密度函数（`probability density function`） $p(x)$ 定义为：当$\\delta x>0, \\delta x\\rightarrow 0$ 时，变量 $x$ 落在区间 $(x,x+\\delta x)$ 范围内的概率为 $p(x)\\delta x$。\n", "\n", "$x$ 落在任意区间 $(a,b)$ 的概率为：\n", "\n", "$$\n", "P(x\\in (a,b)) = \\int_{a}^b p(x) dx \n", "$$\n", "\n", "对于概率密度函数，它需要满足两个条件：\n", "\n", "$$\n", "\\begin{align}\n", "p(x)&\\geq 0 \\\\\n", "\\int_{-\\infty}^\\infty p(x) dx&=1\n", "\\end{align}\n", "$$\n", "\n", "对于随机变量间这样的一个变换 $x=g(y)$，设 $x,y$ 的概率密度函数分别为 $p_x(x)$ 和 $p_y(y)$，那么 $x$ 落在 $(x, x+\\delta x)$ 可以对应于 $y$ 落在 $(y, y+\\delta y)$ 的情况，并有 $p_x(x)\\delta x\\simeq p_y(y)\\delta y$，从而：\n", "\n", "$$\n", "p_y(y)=p_x(x)\\left|\\frac{dy}{dx}\\right|=p_x(g(y))\\left|g'(y)\\right|\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 累积分布函数" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "累积分布函数（`cumulative distribution function`）定义为：\n", "\n", "$$\n", "P(z)=\\int_{-\\infty}^zp(x) dx\n", "$$\n", "\n", "满足：\n", "\n", "$$\n", "P'(x) = p(x)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 联合概率密度分布" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "对于多个变量 $x_1, \\dots, x_D$（用 $\\mathbf x$ 表示），定义联合概率密度为 $p(\\mathbf x)=p(x_1,\\dots,x_D)$，满足当 $\\mathbf x$ 落在一个包含 $\\mathbf x$ 的足够小的空间体积 $\\delta \\mathbf x$ 中时，其概率为 $p(\\mathbf x)\\delta \\mathbf x$。\n", "\n", "它也需要满足：\n", "\n", "$$\n", "\\begin{align}\n", "p(\\mathbf x)&\\geq 0 \\\\\n", "\\int p(\\mathbf x) dx&=1\n", "\\end{align}\n", "$$\n", "\n", "其中积分是对整个空间进行积分。\n", "\n", "对于概率密度函数，我们的加法法则乘法法则仍然适用：\n", "\n", "$$\n", "\\begin{align}\n", "p(x) & = \\int p(x,y) dy \\\\\n", "p(x,y) & = p(y|x)p(x)\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 概率质量函数" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "对于离散的 $x$ 我们有时候将概率密度函数叫做概率质量函数，因为它的质量相当于集中在了某个允许的 $x$ 上。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.2.2 期望和方差" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 期望" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "函数 $f(x)$ 在概率密度函数 $p(x)$ 下的均值叫做 $f(x)$ 的期望（`expectation`）。\n", "\n", "离散分布下，定义为：\n", "\n", "$$\\mathbb E\\left[f\\right] = \\sum_x p(x) f(x)$$\n", "\n", "连续分布下，定义为：\n", "\n", "$$\n", "\\mathbb E\\left[f\\right] = \\int p(x) f(x) dx\n", "$$\n", "\n", "如果给定了 $N$ 个从 $p(x)$ 中随机抽样的点，那么期望可以近似为\n", "\n", "$$\n", "\\mathbb E\\left[f\\right] \\simeq \\frac 1 N \\sum_{n=1}^N f(x_n)\n", "$$\n", "\n", "当 $N\\rightarrow \\infty$ 是，等式成立。\n", "\n", "多元函数可以对其中的一个参数求期望，例如 $\\mathbb E_x\\left[f(x,y)\\right]$ 是函数 $f(x,y)$ 在概率密度 $f(x)$ 上的期望，注意，它是一个关于 $y$ 的函数。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 条件期望 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "我们可以考虑一个函数 $f(x)$ 在条件分布 $p(x|y)$ 下的条件期望（`conditional expectation`），当 $x$ 是离散变量，定义为：\n", "\n", "$$\n", "\\mathbb E_x[f|y]= \\sum_x p(x|y) f(x)\n", "$$\n", "\n", "当 $x$ 是连续变量时，定义为\n", "\n", "$$\n", "\\mathbb E_x[f|y]= \\int p(x|y) f(x) dx\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 方差" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$f(x)$ 的方差（`variance`）定义为：\n", "\n", "$$\n", "\\operatorname{var}[f] = \\mathbb E\\left[(f(x) - \\mathbb E[f(x)])^2\\right]\n", "$$\n", "\n", "它表示 $f(x)$ 偏离其均值 $\\mathbb E[f(x)]$ 的程度。\n", "\n", "对平方进行展开之后，方差可以写成：\n", "\n", "$$\n", "\\operatorname{var}[f] = \\mathbb E[f(x)^2] - \\mathbb E[f(x)]^2\n", "$$\n", "\n", "特别地，我们考虑 $x$ 本身的方差：\n", "\n", "$$\n", "\\operatorname{var}[x] = \\mathbb E[x^2] - \\mathbb E[x]^2\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 协方差" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "对于两个随机变量 $x, y$，其协方差（`covariance`）定义为：\n", "\n", "$$\n", "\\begin{align}\n", "\\operatorname{cov}[x,y] \n", "& = \\mathbb E_{x,y} \\left[(x-\\mathbb E[x])(y-\\mathbb E[y])\\right] \\\\\n", "& = \\mathbb E_{x,y} [xy] - \\mathbb E[x]\\mathbb E[y] \n", "\\end{align}\n", "$$\n", "\n", "当变量 $x,y$ 独立时，协方差为 $0$。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 协方差矩阵" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "对于两个随机向量 $\\mathbf{x, y}$，其协方差为一个矩阵：\n", "\n", "$$\n", "\\begin{align}\n", "\\operatorname{cov}[\\mathbf{x, y}] \n", "& = \\mathbb E_{\\mathbf{x, y}} \\left[(\\mathbf x-\\mathbb E[\\mathbf x])(\\mathbf y^\\top-\\mathbb E[\\mathbf y^\\top])\\right] \\\\\n", "& = \\mathbb E_{\\mathbf{x, y}} [\\mathbf{xy^\\top}] - \\mathbb E[\\mathbf x]\\mathbb E[\\mathbf y^\\top] \n", "\\end{align}\n", "$$\n", "\n", "如果只是考虑随机向量 $\\mathbf x$ 自身分量之间的协方差，那么我们有 $\\operatorname{cov}[\\mathbf x] \\equiv \\operatorname{cov}[\\mathbf{x, x}]$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.2.3 Bayes 概率" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "假设我们有一组模型的参数 $\\mathbf w$，并且做出假定：这组参数服从一定的先验概率分布 $p(\\mathbf w)$。\n", "\n", "$\\mathcal D = \\{t_1, \\dots, t_n\\}$ 是我们观测到的一组数据，这组数据在参数 $\\mathbf w$ 下的条件概率分布为 $p(\\mathcal D|\\mathbf w)$。\n", "\n", "`Bayes` 公式告诉我们：\n", "\n", "$$\n", "p(\\mathbf w|\\mathcal D)=\\frac{p(\\mathcal D|\\mathbf w)p(\\mathbf w)}{p(\\mathcal D)}\n", "$$\n", "\n", "这给了我们一种衡量在观测到数据 $\\mathcal D$ 的情况下，参数 $\\mathbf w$ 的不确定性的方法。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 似然函数" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$p(\\mathcal D|\\mathbf w)$ 可以看成是给定观测数据 $\\mathcal D$ 的情况下关于参数向量 $\\mathbf w$ 的一个函数，通常叫做似然函数（`likelihood function`）。\n", "\n", "似然函数反映了在给定一组参数 $\\mathbf w$ 的情况下，生成这组观测数据的一种可能性。注意它并不是一个关于 $w$ 的概率分布。\n", "\n", "给定似然函数的定义，我们可以将上面的 Bayes 公式表示为：\n", "\n", "$$\n", "\\text{posterior} \\propto \\text{likelihood} \\times \\text{prior} \n", "$$\n", "\n", "这三个量都是 $\\mathbf w$ 的函数。\n", "\n", "对于分母，在给定观测数据 $\\mathcal D$ 的情况下是一个归一化常数，可以写成：\n", "\n", "$$\n", "p(\\mathcal D)=\\int p(\\mathcal D|\\mathbf w)p(\\mathbf w)d\\mathbf w\n", "$$\n", "\n", "在 Bayes 学派和频率学派的眼中，似然函数 $p(\\mathcal D|\\mathbf w)$ 都扮演了一个重要角色。但二者对于似然函数的使用方式是截然不同的。\n", "\n", "在频率学派眼中，$\\mathbf w$ 被看成是一个固定的参数，其值由某些估计量来决定，误差的计算要考虑数据 $\\mathcal D$ 的分布；在 Bayes 学派眼中，数据集 $\\mathcal D$ 是唯一的，参数 $\\mathbf w$的不确定性只来自于 $\\mathbf w$ 的一个概率分布。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 频率估计量：最大似然估计" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "一个最常用的频率估计量就是最大似然（`maximum likelihood`）估计：\n", "\n", "$$\n", "\\mathbf w_{ML} = \\operatorname{argmax}_{\\mathbf w} p(\\mathcal D|\\mathbf w)\n", "$$\n", "\n", "这相当于选择使得观测数据 $\\mathcal D$ 出现的概率最大化的 $\\mathbf w$。\n", "\n", "在机器学习的文献中，似然函数的负对数通常被叫做一个损失函数（`error function`），因为负对数函数是单调递减的函数，因此最大似然就相当于最小化损失函数。\n", "\n", "一种衡量频率估计量误差大小的方法是 `bootstrap`：假设我们有 $N$ 个数据 $\\mathbf X=\\{\\mathbf x_1,\\dots,\\mathbf x_N\\}$，我们从这写数据中有放回的抽样 $N$ 个数据得到一组新的数据 $\\mathbf X_{\\text B}$。\n", "重复进行 $L$ 次这样的操作，得到 $L$ 组这样的抽样数据，然后可以通过这 $L$ 组抽样数据的统计结果估计最大似然估计的误差。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bayes 估计" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`Bayes` 估计的一个重要观点是先验知识的引入，根据后验概率来决定参数 $\\mathbf w$。\n", "\n", "考虑抛硬币的情况，假设我们抛三次，每次都得到正面，那么最大似然估计会得到这枚硬币会 `100%` 得到正面的结论，而 `Bayes` 估计不会得到这么极端的结论。\n", "\n", "两者并没有什么好坏之分，只不过是看问题的角度不同。`Bayes` 估计如果选定的先验不好，也可能得到很差的结果。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.2.4 高斯分布" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "高斯分布（`Gaussian distribution`），又叫正态分布（`normal distribution`）。\n", "\n", "对于实值变量 $x$，高斯分布定义为：\n", "\n", "$$\n", "\\mathcal{N}\\left(x\\left|~\\mu,\\sigma^2\\right.\\right) = \\frac{1}{(2\\pi\\sigma^2)^{1/2}} \\exp\\left\\{-\\frac{1}{2\\sigma^2}(x-\\mu)^2\\right\\}\n", "$$\n", "\n", "参数为均值 $\\mu$ 和方差 $\\sigma^2$。方差的平方根 $\\sigma$ 叫做标准差，方差的倒数 $\\beta = \\frac{1}{\\sigma^2}$ 叫做精度。\n", "\n", "其图像如下所示：" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEDCAYAAAA7jc+ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8zmX+x/HXhRqkKNKuXRSJVEiStEwzJbRLmqKFVstI\nZVCJSkxoOlFNSY0aMfxkKWuIEJUsWUuFbKmsJ1y/Pz7KkuW4t+u+v/f7+Xich0e3c+7zyTRv17m+\n1/X5OO89IiISXflCFyAiIsmloBcRiTgFvYhIxCnoRUQiTkEvIhJxCnoRkYgrELqA3TnndN5TRCQG\n3nu3p9fTLugBdLZf0lX79u1p37596DJE/sC5PWY8oK0bEZHIU9CLiEScgl7kANSsWTN0CSIHzKXb\nfrhzzqdbTSIi6c45t9eHsVrRi4hEnIJeRCTiFPQiIhGnoBcRiTgFvYhIxCnoRUQiTkEvIhJxCnoR\nkYhT0IuIRJyCXkQk4hT0IiIRp6AXEYk4Bb2ISMQp6EVEIk5BLyIScQp6EZGIU9CLiEScgl5EJOIU\n9CIiEaegFxGJOAW9iEjEKehFRCJOQS8iEnEKehGRiFPQi4hEnIJeRCTiFPQiIhGnoBcRiTgFvYhI\nxCnoRUQiTkEvIhJxCnoRkYhT0IuIRJyCXkQk4hT0IiIRp6AXEYk4Bb2ISMQp6EVEIq5A6AJE0tKm\nTTBxIkydCnPmwLJlsH49OAfFi8MJJ0CFCnDhhVC2rL0ukqac9z50Dbtwzvl0q0myxNatMHw4/Pvf\n8MEHUK4cVKlivx53HBQpAtu2wapV8M03MGMGjBtnX3vTTXDHHVC6dNh/B8lazjm893tccWhFL7Jl\nC/TpA507Q7Fi0Lgx9OoFRxyx/6/1Hr78Et58Ey66CKpVg3/8AypWTH7dInmkFb1ktzFj4IEHLNSf\neAJq1Ih9G2bDBnjlFejUCa68Ep55BkqWTGy9InuxrxW9HsZKdvrlF7j7brjtNujQAcaOhYsvjm+v\nvXBh+0vjq69sH//ss6Ffv4SVLBIrregl+3z2GdxwA1SvDt26QdGiyfk+U6dCw4a2ndOzp/1FIJIk\nWtGL/OY//4HLLoP27eG115IX8gDnnQfTpsGvv9rpnKVLk/e9RPZBK3rJDt7bHvzrr8PgwVC+fGq/\nd+fOkJMDQ4fCWWel7ntL1tCpG8lu3kPLljBqFEyeDEcdldrv7xy0aQOlSsEll9i+fa1aqa1BsppW\n9BJt27ZBs2YwfbqdkT/88LD1jBkDN95oxzmvvDJsLRIp+1rRK+glurZtgzvvhIULYcgQOOyw0BWZ\njz+Ga6+Fd9+FmjVDVyMRoYexkn28h4cfhkWLYNiw9Al5sFM477wD118PkyaFrkaygIJeoum552D0\naBg0CA45JHQ1f3TJJbZ9c+21dtxTJIkU9BI9ffvaufVhw6ylQbr685+hRw+oUweWLw9djUSY9ugl\nWkaPhptvtl8z5RjjE0/YscsxY6BQodDVSIbSw1jJDt98AxdcAG+/nVnHF72HW26xY5hvvaWWxxIT\nPYyV6Nu0CerXh1atMivkwYL9tdfsdNDTT4euRiJIK3rJfN7bMcr16+0yUqauiJctg3PPtWcMmfaX\nlQSnFb1EW69eMGUKvPpq5oY8wDHH2Emchg31cFYSSit6yWwzZ9rqd+LE6Ex3+sc/7FLViBGQP3/o\naiRDaEUv0bRxoz3EfO656IQ8QLt2NtawY8fQlUhEaEUvmeuBB+CHHzJ7X35vli2DSpWgf39rcSyy\nH+peKdEzdKjdev3ss+iFPNh+/UsvQaNG9u9YpEjoiiSDaUUvmWflSqhQwYaIXHxx6GqS67bbLOT/\n9a/QlUia04UpiZabboITTrC9+ahbu9Zmz77yClx+eehqJI3pYaxEx8CB1lv+iSdCV5IaxYrZsdHG\njS30RWKgFb1kjjVrbATgO+/YYO9s0rQpbN5soS+yB9q6kWi4/XbrK9+9e+hKUu/nn61J25tvaliJ\n7JFO3UjmGz4cxo2zC1LZ6LDDrPXy3XfD559DwYKhK5IMoj16SX8bNsC998LLL2f3McM6daBcOV2k\nkgOmrRtJf489Zp0d+/ULXUl4S5fa0dIxYyz0RbZL6B69c+5PwMXARcAxQEnAASuApcAEYJz3flOM\nxSroZYc5c6BGDduuOPbY0NWkh5wc26sfPx7y6YdyMQk5XumcO9M51w/4DmgFFAIWAu8Dg4GZ21+7\nH1jknPuPc+7MeIuXLOa9nTZp21Yhv7O77oJff7WwF8mDPK3onXOPAJcA/wQ+8N5v3c/n5wOuAJoB\nE733nfJckFb08pu+faFrV2tBXEDnBnYxdSpcc439xJPOc3ElZeLaunHOtQNmeu8HxPjN6wDneO87\n5PHzFfRixwnPOMP62Zx/fuhq0tNdd9mM2RdeCF2JpIF4g/5k7/3iOAs4xXu/KI+fq6AX+PvfYdUq\nG7Ene7ZqFZx5JowcaW0SJKvpwpRklgULoEoV+PJLOPro0NWkt5wcG4Y+blw0u3hKnsX1MNY5V9U5\n19I599edXjvJOXeTc06NsiXxWra0D4X8/jVpAuvWWVsIkb3Y54reOdcMaAlMBY7b/nId7/0q59xJ\nwELvfUJnnWlFn+VGjbLwmj1btz/zauxYaw8xd67+zLJYPCv684Ey3vsbvPcXYqdoujvnSgIbsfPz\nIomxZQs89BA8/7wC60DUrAkVK+qhrOzV/oL+c+/95t/+wXv/GXAHcC9wIqCltyROr15w5JFw7bWh\nK8k8zz5r/fl/+CF0JZKG9rd1cz1QBOgAXOW9/3Kn37sbeNF7n9ADztq6yVI//ghlysCHH+oESaya\nN7e+QDk5oSuRAOI9XnkKUB5433u/Zbffu9B7PzFhlaKgz1oPPQSbNimk4vHjj3b3YPRo9cHJQjpe\nKelt4UK44AK75XnkkaGryWzdu9vg9OHDQ1ciKZaUUYLOuZHOucKxlyWyXfv28MADCvlEuPdeWLxY\nQS+7iHlF75zbBpzkvV+S0IK0os8uM2fCZZfB/Plw6KGhq4mGwYOhTRvr+KkeQVkjmcPBr3TOdXXO\nveCcu2V7C2ORvHv8cWjdWiGfSFdfDSVLqn2E/C7eFf0PwFvAwUAV4GjgVu/9RzEXpBV99pg8GW64\nAebN07n5RJsyBerVsz/bwtphzQZJeRi7Pehre+9H7/RaReBF4AHv/bQY31dBnw28h0svhQYN4M47\nQ1cTTfXrW+fP1q1DVyIpkKytm5XAjzu/4L2fgfWhvzuO95VsMHIkfP89NGoUupLo6tgRunSxY5eS\n1eIJ+r7AA7u/6L3/BRspKLJn3sOjj8KTT+phYTKVKWO3jJ95JnQlElg8Qd8WOMM51985V/a3F51z\n+dnRAE3kjwYOtL42110XupLoa9cOeve2n54ka8Uc9N77DcClwBJgunPuG+fcBGAuENM0KskCW7fa\nSZunn9Zg61Q4/nh7BvLEE6ErkYAScjPWOVcMqA4UBMZ571fG8V56GBtlb7wBr76qQRmptGaNtUaY\nMMF+lUhSCwRJD5s3W9D07QvVq4euJrt06gQzZsC774auRJIkmRemRPKud2846yyFfAgPPggTJ8K0\nmE49S4bTil5SY/16OO00GDYMzjkndDV5tmnTJp5//nlWrVrF3LlzyZcvH507d6Z8+fKhSztwOTnQ\nv78dbZXI0YpewnvhBbj44owKeYAnn3ySRo0a0a1bN4YNG0alSpWoXr06CxYsCF3agbvzTliyREGf\nhRK+onfOnQPUA7YCK4D+B/JwViv6CPrxRyhdGj7+GE4/PXQ1ebZ582ZKlChB69atefzxxwFYt24d\nRxxxBPfccw/du3cPXGEM3n3XplFNnaqH4RGT6hX9COA+730HYBBwu3Pu3iR8H8kUzz4LdetmVMgD\nbN26lRIlSrB+/frfXytSpAjFixfPzBU92N0F72GATkBnk2RcS7yH7UPDvfdLgeecc7pAla2WLbNZ\nsJ9/HrqSA1a4cGEWL168y2sbNmxgxYoVnHbaabu81rZtW/r27cvKlbv+8FqrVi1GptNWSb581hrh\n4YehTh3dTM4SCV/Re+8Heu8H7PaaruVlq44d4fbb7eJOBLz11lsUKlSIhx56CICNGzdSu3Zt1qxZ\nw/Dhwxk5ciRHHXUUOTk5fP311wwePDhwxXtwxRU25KVv39CVSIrE072ywB5myLbGfkp413s/P8b3\n1R59VCxeDJUrw9y5kZgetXbtWipUqMBTTz1Fw4YNAWjevDlLliyhf//+v39ey5Yt2bZtG127dg1V\n6v5NmAC33gpffQV/0hiJKNjXHn08P7ctds5tAMYDHwHjvffPOOcc0AVoEcd7SxS0awf33x+JkPfe\nc+edd9KhQ4ffQ37t2rW89NJLTJ8+fZfPzc3NxaX7g87q1e1OQ+/ecN99oauRJIsn6KsCVwIXAR2A\nE51z3wOz0bFNmTULRoywEYER0KZNG2699Vbq1q0LwMKFC5kzZw4lS5akbNmyu3zu5MmTf9/aSWtP\nPQVXXQV/+xscckjoaiSJ4mlq9p33/hXvfSPv/clAKaANUBRonKgCJUM9/jj8/e9w2GGhK4lbr169\nqFq16u8hD9CnTx82bdrECSecsMvnzpgxg+XLl1O/fv1Ul3ngKlaEGjUgE4+JygFJxjn6okBT732n\nGL9ee/SZ7pNP7BjfvHlQqFDoauIyePBgunfvzmWXXfb7a+vXr6dQoUI0btyYatWqMXfuXPLnz88v\nv/zCFVdcQefOnalRo0bAqg/AV1/ZNs68eXD44aGrkTgka5TghcAC7/0Pe/i9x733T8X4vgr6TFe7\nNtx4IzRpErqSuKxZs4ZSpUqxcePGP/zeO++8w3XXXcfQoUMZMGAAxx13HEuXLuXuu++mcuXKAaqN\nQ+PGcNRRdkJKMlaygn4qUAmYD4zDHsjOAvIDj3jvr4/xfRX0mWzUKLjnHpg9Gw46KHQ1khdLltg2\nzuzZFviSkZJyM9Z7fx5wBtAZC/e2wHRgCnCqc+4u51wGdn6SmO08InCnkP/2229p2rQpvXr1Clic\neO+pX78+vXv3Jjc3d8dvlCoFDRtqRR9hcZ2O8d4v8N6/7r1v7L0vA5QErgNGA3cA05xza51zfRJQ\nq6S7QYOs5/wNNwA7Ar5ChQoceuihmfGAMsKcc7Rq1Yr+/ftTunTpXQP/0Ufhrbfgm2/CFilJkdQ2\nxc65QkAVoJT3/o08fo22bjLR1q1w9tnW1+Yvf6FPnz40atSIQw45hGOPPZYCumqfVtavX893331H\nkSJF+Omnn+zFtm1ttuxrr4UtTmKSrAtT++W93wiMSeb3kDTx9tt2auOqqwC48cYbWbhwIb169eLo\no4+madOmnH322YGLlNzcXP73v//9fmS0devWO36zRQtrPDd3LpQpE65ISTgNHpH45eZaMLz+up3L\n3snmzZv597//zdNPP02LFi148MEHw9QoeO+pVKkSRx11FO3ataNq1ap//KRnnoFPP9XIwQwU96kb\n51x/4IhYvz+wxnufpw1aBX0GevFFGDLEpkftRW5uLlu2bKFw4cIpLEx2t3btWooVK7b3T1i/3lb1\nQ4ZApUqpK0zipuHgkjwKhujJw1/ckn40SlCSp0cPu1mpkI+OJk1sn/6jj0JXIgmiFb3E7rcRgRMm\nwBlnhK5GEqlPH+ts+dFHGjmYIZKyonfOPe+c+8E5V3qn1xo45/axASiR0qWLTSlSyEdPgwawerW2\nbyIinhYIjwN/AnJ+myDlnCsB3AcM8N5/EeP7akWfCZYvt37mM2bYzUqJngED7Jbzp5/aCEJJa8na\no/fASzuPCfTer/Letwdi6nMjGaRjR7jtNoV8lNWtazNld5qeJZkpnhV9WWAsMGP7r+OAqd77Lc65\n5733MU2Y0oo+A3z9NZx7bmRGBMo+fPihTaCaNUuDxNNcslb0HYDHgMnA5cAo4Gfn3LfAujjeV9Jd\n+/bQrJlCPhvUrg3HHmsPZyVjxbOib7PzcBHn3MHYeMGrsT36j2N8X63o09ns2VCzpo0ILFo0dDWS\nCpMmwU032ZCSggVDVyN7kawV/Qbn3O/HLbz3ud77cVi74ivjeF9JZ23bQqtWCvlsUrUqVKgAL78c\nuhKJUVzn6J1zDwD5vffdtv/zJcBQbEXfIMb31Io+XUVoRKAcoC++gMsvhwULoEiR0NXIHqSsBYJz\nrgBwF/CR9/7LGN9DQZ+OvIdatex8dWPNfs9Kt9xiR2ofeyx0JbIHcQW9c+487/3UOAs433s/JY+f\nq6BPR8OHw8MPw8yZOn2RrebPt22cefPgiFh7HEqyxLtHf45z7v44vvnDQIVYv17SwLZt8MgjdnZe\nIZ+9Tj8d6tWz4TKSUfYb9N773sAK59wE51xj51zx/X2Nc+5I59w9zrnJwLfb30MyVb9+dtqibt3Q\nlUho//iH9cBZtix0JXIA8rxH75wrCTQHbgZysUHgK4Htc8g4HCgOVAQOBt4C/um9X3lABWnrJr3k\n5kLZsvDqq3asUqRFC5sN3LNn6EpkJzHv0TvnemCnapru9noZoDpwHPDbrZmVwPfAeO/9V3EUq6BP\nJz17wvvvq7mV7LBypU0UmzYNTj45dDWyXTxB/yzQAijhvf8xSfXt/j0V9Oli3Trblx0+3M5Ri/ym\nXTtrhfHGG6Erke3ieRjbB9gK1Ex0UZIBunaFSy9VyMsftWhhP+XNmhW6EsmDvByvvBBoAiwDhgMT\nvfdbklaQVvTpYcUK25vXj+eyN126WHuE994LXYmQoAtTzrkjseZlFwIHAfOxrpXTvPfbElOqgj5t\nPPigXZLq3j10JZKuNm6E006D//0PzjsvdDVZLyk3Y51zp2BbOpWx4F+GdbL8xHu/OrZSFfRpYfFi\nqFwZ5syBkiVDVyPpLCfHBpR88EHoSrJeSlogbD9+WQ1rXbzYe39DjO+joA+tQQN7CNu+fehKJN39\n+qtt8eXkWEtjCSZZ3St3VxH4O3AuoK5HmWrKFBg7Flq2DF2JZIKDDoJnnrGHs1u3hq5G9iKuoHfO\nHeuce8Q5Nw8Yhp3Qucp7f1VCqpPU8h6aN7c5oepQKHlVrx4cdpiOWqaxA966cc6dBly7/aMK8DPw\nX+DVvDYu28/7a+smlPfegyeegOnTIX/+0NVIJpkyxVpkfPWVFgmBxL1H75yrDvwFmx5VFliPreD7\nAe9773MTWKyCPoTNm60FrfZaJVYNGtgpnA4dQleSlRIR9CuAEsAcoCPQP5Hhvtv3UtCH0LUrjB4N\nQ4aErkQy1ZIlULEifP45HH986GqyTiKCvjS2kj91+6+HAWuBicBo7/13CSxWQZ9qq1db75KPPrIT\nFCKxevRRWLoUXn89dCVZJ1nn6A8CzgNqAydhTc0Geu8nx1jnb++roE+1Bx+ELVvgxRdDVyKZ7uef\n4YwzrBFepUqhq8kqqTpHfzhQB6gBjPDevxPj+yjoU2nePKhWzS5HHXnk/j9fZH9eftlmGIweDW6P\nuSNJkLKZsdu/WXFgBdDOe/9UDF+voE+lq6+G6tWhdevQlUhUbNkC55xjJ7jq1QtdTdZI1YWp31wP\nOEB359Pd0KF2HO6hh0JXIlFSoID1SGreHDZsCF2NkJwV/THYGfuB3vvlMXy9VvSpsHkzlCtn/4f8\n859DVyNRdMMNcOaZaqWRIinduomXgj5FOneGjz+GwYNDVyJRtWSJPZCdOlWtrlNAQS+7+u4720P9\n5BM49dTQ1UiUdexoMw0GDgxdSeSleo9e0l2rVnDvvQp5Sb4WLWDmTBgxInQlWU0r+mwzdiw0amTH\nKQsXDl2NZIMhQ3YE/sEHh64msrSiF7NlCzzwADz/vEJeUuevf7UeOC+8ELqSrKUVfTbp3h0GDYKR\nI3WRRVJr/nyoWtX64Bx3XOhqIkkPY8UewFasqH42Ek7btrZl2L9/6EoiSVs3AvffD82aKeQlnMce\ngy++0JHeABT02WDgQFtJtWkTuhLJZgULQq9ecN998MsvoavJKtq6ibqffrKBIm+/DTVqhK5GBO64\nAw49VA9nE0x79Nnsvvus3UHv3qErETGrV1v7jUGD4PzzQ1cTGfsK+gKpLkZSaNIkGDAAZs0KXYnI\nDsWL2xHfJk3s1uxBB4WuKPK0Rx9Vublw113QrRscfnjoakR2dfPNcMwx9t+nJJ22bqKqY0eYONEm\n/ejMvKSjxYvhvPPUcylBtEefbb78Ei65xH4sPvHE0NWI7F3XrrZXP2YM5NMGQzx0jj6b/Pqr9bLp\n1EkhL+nvwQdh2za7tS1JoxV91Dz5pPWZHzpUWzaSGRYsgCpV7L/b0qVDV5OxtHWTLT77DC6/HKZP\nh+OPD12NSN717Gl3PcaPh/z5Q1eTkbR1kw02b7Ytmy5dFPKSeZo2hT/9yY5dSsJpRR8VrVrZj8AD\nBmjLRjLT11/bBaoRI6wBnxwQXZiKupEj4T//sRawCnnJVCedBP/8J9xyC3z6qWYmJJBW9Jlu1Sqb\n//r661C7duhqROLXsCEccgjk5ISuJKPoYWxUeQ9168Lpp8Nzz4WuRiQxfvrJtm66dYM6dUJXkzG0\ndRNVOTnwzTfwzjuhKxFJnKJFoW9fqFcPKlWCE04IXVHG04o+U02fDldcYW0OdPZYoqhTJxssPnas\nGp/lgY5XRs3atXD99fDiiwp5ia7WraFYMQ3MSQCt6DON9/Yj7fHHQ48eoasRSa7Vq+Hcc21Iifbr\n90l79FHSrRt8/z306xe6EpHkK17cnkFdfTWULw+nnBK6ooykFX0mGTUKGjSAyZPtzLFItujZ0+bN\nfvwxFCkSupq0pOOVUbBoEVSrZhejLrkkdDUiqeU9NG5sRy/ffVctjfdAD2Mz3bp1tj/5+OMKeclO\nzsG//mXblh07hq4m42hFn+62bbMTNocfbgO+1eJAstmyZdYPp0cPuPba0NWkFT2MzWSPPgrLl1sL\nV4W8ZLtjjrHGfVddZRepzj03dEUZQVs36Swnx/6jHjTIWriKiM2Z7dULrrnGOl7KfmlFn67efx86\ndIAJE6BEidDViKSXunVhyRJb2U+caFubslfao09Hn34KV14J//d/NmJNRPbsoYesPfewYVCwYOhq\ngtLxykwyZw7UqgUvvaSHTSL7s3Ur3HQTbNkC//0vFMjeTQodr8wUixfbzNdnnlHIi+RF/vzw1luw\naRPccYedUpM/UNCni6VLbXDII4/AbbeFrkYkcxx8MLz3nj2Yvf9+u1wlu1DQp4Plyy3k77wTmjUL\nXY1I5ilc2J5pffKJzU9W2O9CQR/a0qVQsybcfLOdmReR2BQtCh98AGPGwMMPK+x3oqAP6dtv4eKL\n4fbboW3b0NWIZL4jjrDmf5MmwX33ac9+OwV9KIsWWcjfe6/ty4tIYhQrZiv7zz6Du++2kzlZTkEf\nwvTpUL267SU2bx66GpHoKVoUhg+3B7T168PGjaErCkpBn2offmiXoXr2tNW8iCTHoYfaDfNDDrHD\nDmvWhK4oGAV9Kr35Jtx6qx0Fq1cvdDUi0Xfwwfb/u6pV4cIL7a5KFsrea2SptHWrDTh+7z0YPRrO\nOit0RSLZI18+6NIFTjzRAr9fPzvplkXUAiHZfvrJjk5u2mRXtIsXD12RSPYaOdLGcbZvH7mtU7VA\nCGXmTLjgAjj1VBgxQiEvElrt2tbtskcPaNIkax7SKuiTwXubBlWrlm3Z9OgBBx0UuioRATjtNLtB\nu369LcTmzg1dUdJp6ybRfv4Z7rnHVvPvvgtly4auSET2xHt45RW7kd6li/WYyuApbtq6SZVRo+Ds\ns+Gww2DKFIW8SDpzzrZvRo+GZ5+F666DFStCV5UUCvpEWLcOmja1VgYvv2wjAAsVCl2ViORF+fI2\n7Of0022h1r9/6IoSTkEfD+9tnmu5crBhg23XXHFF6KpE5EAVLAidO8PAgfDYYzYP4ptvQleVMAr6\nWC1aBFdfDa1bw2uvweuvW48NEclcVavaaMJzz7WPTp0gNzd0VXFT0B+otWutCdl559lNuy++sNM1\nIhINBQtaN9kpU+woZrlytp2TwYdEFPR5tXkzdOsGpUvDqlUW8G3a2BVrEYmeU06BIUOsL9VTT0G1\najB+fOiqYqKg358NG+CFF+zS0+jRNtTglVfguONCVyYiqXD55dZxtlkzO4J56aUwdmxGrfAV9Huz\nZo0N6T71VBg3DgYPtlFl6lMjkn3y5bOGhPPm2a9NmkCNGrbiz4DhJrowtbtZs6B7d7vsdM010LKl\nHb8SEfnNli3Wu6prV/jxRxtK/re/2R2aQHRhan/WrrXz79WqWS+MY4+1a9FvvKGQF5E/KlDAmhVO\nmQJ9+sDHH0OpUra1M2pU2q3ys3dF/+uv1snujTdsEs1ll0GjRnYOXn1pRORArVgBb79tmbJ6NTRs\naB9lyqTk2+9rRZ9dQb92rYX6oEHWTbJMGfsf4sYbbaiwiEgifPGFBX6/fradc8019lGlCuTPn5Rv\nmb1Bv2WLPS0fO9ZG+H3yiT1AqVMH/vpXOOaYxHwfEZE92bbNMmjQIPtYvhyuusru3tSsads9CZI9\nQb9+PcyYAZMmWbhPmGB/kDVr2h/sZZdBkSKJLFdEJO8WL4ahQy2fxo611X7NmnDxxXD++XZPJ19s\nj06jGfSrVsHs2fYj0rRp9rFokd1iO//8HX94Rx6Z9JpFRA7Ytm0wZ47dzRk/3jJs5UqoVAkqV7YW\nDGedZeFfsOB+3y4zg957O7b09dc7PubPtz+Y2bOt/8SZZ1qwV65sH+XK6aaqJNXYsWOpmWXzRiWF\nVq+2TprTptmWz+zZtoAtVcranpctazd2TzrJPkqV+v0vgcwL+vLlLdidg5NP3vEvdcopFu5nnmn7\n6xk8JEAyU/v27Wnfvn3oMiSb5ObCggUW+nPn7rr4/fZbG1F60km4SZP2GvQFUlpwXvXpY8GubpAi\nku0OPnjHAnd3W7fC0qUW+jVq7PUt0jPozzkndAUiIukvf3444QT72Ie03LoJXYOISCbKmD16ERFJ\nLPW6ERGJOAW9iEjEKehFRCJOQS8iEnEKehGRiFPQi4hEnIJeJA+cc1Odcxfu9M/HO+e+dc6l56VD\nkZ0o6EWYSigPAAABGElEQVT2wzl3GlAJWLDTy1cCG7z3W8JUJZJ3CnqR/asOzPfe/7DTaxcB4wPV\nI3JAFPQi+1cdGLfbazWAjwLUInLAFPQi+1ednULdOXc8cCJa0UuGUNCL7INzrgRQGpi108s1ge+9\n94udc62DFCZyABT0Ivt2EeCB/ADOuaLAfcBs55wjXVt9i+xE3StF9sE51wWoBSwEPge2Am8DrwBf\nADne+/nhKhTZPwW9yD445yYDr3nve4WuRSRW2roR2QvnXCGgIjApdC0i8VDQi+xdFWCj935m6EJE\n4qGgF9m7UsDg0EWIxEt79CIiEacVvYhIxCnoRUQiTkEvIhJxCnoRkYhT0IuIRJyCXkQk4hT0IiIR\np6AXEYm4/wcSdOSRJGqxUQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import scipy as sp\n", "import matplotlib.pyplot as plt\n", "\n", "%matplotlib inline\n", "\n", "from scipy.stats import norm\n", "\n", "xx = np.linspace(-3, 3, 200)\n", "\n", "norm_xx = norm.pdf(xx)\n", "\n", "fig, ax = plt.subplots()\n", "\n", "ax.plot(xx, norm_xx, \"r\")\n", "ax.set_ylim(0, 0.5)\n", "ax.set_ylabel(r\"$\\mathcal{N}\\left(x|\\mu,\\sigma^2\\right)$\", fontsize=\"xx-large\")\n", "ax.set_yticks([])\n", "ax.set_yticklabels([])\n", "\n", "ax.set_xticks([0])\n", "ax.set_xticklabels([r\"$\\mu$\"], fontsize=\"xx-large\")\n", "\n", "ax.text(-.1, 0.25, \"$2\\sigma$\", fontsize=\"xx-large\")\n", "\n", "ax.annotate(\"\",\n", " xy=(-1, 0.24), xycoords='data',\n", " xytext=(1, 0.24), textcoords='data',\n", " arrowprops=dict(arrowstyle=\"<->\",\n", " connectionstyle=\"arc3\"), \n", " )\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "下面证明这是一个概率分布，首先，从定义中可以看出：\n", "\n", "$$\\mathcal{N}\\left(x\\left|~\\mu,\\sigma^2\\right.\\right) > 0$$\n", "\n", "然后是积分：\n", "\n", "$$\\int_{-\\infty}^{\\infty} \\mathcal{N}\\left(x\\left|~\\mu,\\sigma^2\\right.\\right) dx = 1$$\n", "\n", "这个结果可以通过计算 $I = \\int_{\\infty}^{\\infty} \\exp (x^2)dx$ 得到（计算 $I^2$，并换成极坐标计算）。\n", "\n", "$x$ 的期望为（令 $y = x + \\mu$，对 $y$ 的积分变成一个奇函数的积分加上 $\\mu$ 乘以一个高斯分布的积分）：\n", "\n", "$$\n", "\\mathbb E[x] = \\int_{-\\infty}^{\\infty} \\mathcal{N}\\left(x\\left|~\\mu,\\sigma^2\\right.\\right)x~dx = \\mu\n", "$$\n", "\n", "其方差为（等式 $\\int_{-\\infty}^{\\infty} \\exp\\left\\{ -\\frac{1}{2\\sigma^2}(x-\\mu)^2 \\right\\} dx = (2\\pi \\sigma^2)^{1/2}$ 两边对 $\\sigma^2$ 求导）：\n", "\n", "$$\n", "\\text var[x] = \\int_{-\\infty}^{\\infty} \\mathcal{N}\\left(x\\left|~\\mu,\\sigma^2\\right.\\right)(x-\\mu)^2dx =\\sigma^2\n", "$$\n", "\n", "因此：\n", "\n", "$$\n", "\\mathbb E[x^2] = \\mathbb E[x]^2 + \\text{var}[x] = \\mu^2 + \\sigma^2\n", "$$\n", "\n", "概率分布的最大值叫做众数（`mode`），高斯分布的众数就是均值 $\\mu$。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 多维高斯分布" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "对于 $D$ 维的向量 $\\mathbf x$，高斯分布定义为：\n", "\n", "$$\n", "\\mathcal{N}\\left(\\mathbf x\\left|~\\mathbf{\\mu, \\Sigma}\\right.\\right) = \\frac{1}{(2\\pi)^{D/2}} \\frac{1}{|\\mathbf\\Sigma|^{1/2}} \\exp \\left\\{-\\frac{1}{2}(\\mathbf x - \\mathbf \\mu)^\\top\\mathbf\\Sigma^{-1}(\\mathbf x - \\mathbf \\mu)\\right\\}\n", "$$\n", "\n", "其中，$D$ 维向量 $\\mathbf \\mu$ 是均值，$D\\times D$ 矩阵 $\\mathbf\\Sigma$ 是方差，$|\\mathbf\\Sigma|$ 是其行列式。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 最大似然估计" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "假设我们现在有 $N$ 组对 $x$ 的观测数据 $\\mathsf x = (x_1,\\dots,x_N)^{\\text T}$，这些数据是独立同分布（`independent and identically distributed, i.i.d.`）的，都服从一个均值 $\\mu$，方差 $\\sigma^2$ 的高斯分布。那么在给定这些参数的情况下，出现这些观测数据的概率，或者从参数的角度来说，似然函数为：\n", "\n", "$$\n", "p(\\mathsf x~|~\\mu, \\sigma^2)=\\prod_{n=1}^N \\mathcal N\\left(x_n \\left|~\\mu, \\sigma^2\\right.\\right)\n", "$$\n", "\n", "通常最大似然的问题经常转化为求最大对数似然的问题：\n", "\n", "$$\n", "\\ln p(\\mathsf x~|~\\mu, \\sigma^2) = -\\frac{1}{2\\sigma^2} \\sum_{n=1}^N(x_n - \\mu)^2 - \\frac N 2 \\ln \\sigma^2 - \\frac N 2 \\ln(2\\pi) \n", "$$\n", "\n", "对 $\\mu$ 最大化，我们得到最大似然解：\n", "\n", "$$\n", "\\mu_{ML} = \\frac 1 N \\sum_{n=1}^N x_n\n", "$$\n", "\n", "即样本均值。\n", "\n", "对 $\\sigma^2$ 最大化，我们得到：\n", "\n", "$$\n", "\\sigma^2_{ML} = \\frac 1 N \\sum_{n=1}^N (x_n-\\mu_{ML})^2\n", "$$\n", "\n", "即样本方差。\n", "\n", "但是这个解不是无偏的，我们可以计算它们的期望：\n", "\n", "$$\n", "\\begin{align}\n", "\\mathbb E[\\mu_{ML}] & = \\frac 1 N \\sum_{n=1}^N \\mathbb E[x_n] = \\mu \\\\\n", "\\mathbb E[\\sigma^2_{ML}] & = \\mathbb E \\left[\\frac{1}{N} \\sum_{n=1}^N(x_n-\\frac 1 N \\sum_{m=1}^N x_m)\\right] \\\\\n", "& = \\frac 1 N \\sum_{i=1}^N \\mathbb E \\left[x_n^2-\\frac 2 N x_n\\sum_{m=1} x_m + \\frac{1}{N^2} \\sum_{m=1}^N\\sum_{l=1}^N x_mx_l\\right] \\\\\n", "& = (\\mu^2 + \\sigma^2) - 2 (\\mu^2+ \\frac 1 N \\sigma^2) + \\mu^2+ \\frac 1 N \\sigma^2 \\\\\n", "& = \\left(\\frac{N-1}{N}\\right)\\sigma^2\n", "\\end{align}\n", "$$\n", "\n", "用到了：\n", "- 当 $m=n$ 时，$\\mathbb E[x_mx_n]=\\mathbb E[x_n^2] = \\mu^2 + \\sigma^2$\n", "- 当 $m\\neq n$ 时，$\\mathbb E[x_mx_n]=\\mathbb E[x_n]\\mathbb E[x_m] = \\mu^2$\n", "\n", "因此，方差的一个无偏估计为：\n", "\n", "$$\n", "\\tilde \\sigma^2 = \\frac{N}{N-1}\\sigma^2_{ML}=\\frac{1}{N-1}\\sum_{n=1}^N(x_n-\\mu_{ML})^2\n", "$$\n", "\n", "随着 $N$ 的增大，方差估计的误差也随之增大。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.2.5 重新理解曲线拟合" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "对于曲线拟合的问题，设训练集输入为 $\\mathsf x=(x_1, \\dots, x_N)^\\top$，对应的目标值为 $\\mathsf t=(t_1, \\dots, t_N)^\\top$。\n", "\n", "我们将我们的不确定性用高斯分布来表示，假设给定 $x$，对应的目标值 $t$ 服从一个均值为 $y(x,\\mathbf w)$ 的高斯分布：\n", "\n", "$$\n", "p(t\\left|~x,\\mathbf w,\\beta\\right.)=\\mathcal N\\left(t\\left|~y(x,\\mathbf w), \\beta^{-1}\\right.\\right)\n", "$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEDCAYAAAABcbKvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmczuX+x/HXd4aRsRxZxillSVIpEsoWKlSWigynplNO\npUWn7XCKKJE5qBAtZKLCVIxjySEtMhnLlF/Z6hghS1FZZqwzZr1+f1xjjmG2e8x9f+975v18POYx\nM/d8l89I99t1fa/FMcYgIiLibUFuFyAiImWDAkdERHxCgSMiIj6hwBEREZ9Q4IiIiE8ocERExCfK\nuV2AP3McR2PGRUSKwRjjnPmaAqcQmqck/mbePAgP74sxc90uRSRPjnNW1gDqUhMRER9R4IiIiE8o\ncEQCUhO3CxDxmAJHJCApcCTwKHBERMQnFDgiIuITChwREfEJBY6IiPiEAkdERHxCgSMiIj6hwBER\nEZ9Q4IiIiE8ocEREimHXrl3MnDnT7TKIjIwkLS3N7TKKRIEjIuKhxMREhg8fTkREhNulEB4ezoAB\nA9wuo0gUOCIiHhoyZAiDBw8mODjY7VK47LLLaNasGTNmzHC7lEIpcEREPJCQkMCuXbu45ppr3C4l\nx4ABA5gwYQLp6elul1IgBY6IiAfeeOMN7r//frfLyKVKlSq0atWKBQsWnNuFfvsNMjNLpqg8KHBE\nRDywbNky2rRp43YZZ2nbti0LFy4s/gUOHoQOHeDzz0uuqDMocEREsi1ZsoTIyEhuvfXWXCO/+vfv\nz7vvvsvu3bs5fPgwl1xySZ7nT5o0iTFjxtC3b1+2bNlCZGQko0ePZtCgQYXee+LEidSrV4+goCAa\nNGhAdHQ0F198MUFBQVx33XX89ttvdO/enZCQEFq0aMFvv/2W6/xWrVoRHx9fvF/85Em4807o0wdu\nu6141ygKY4w+8vmwfzwi/iUmxhiIcbuMUicxMdG8+uqrJjMz01StWtV8//33xhhjMjMzTZUqVczS\npUvNypUrzeWXX57n+ZMmTTIJCQnGGGOioqJMWFiY2bt3rxk+fLgJCwsrUg0xMTHGcRzTtWtXY4wx\nH330kQkKCjJRUVHGGGP2799vGjdunOe5v/76qwkODjZpaWke/d4mM9OYv/zFmL597dclIPu986z3\nVLVwRESAL7/8kr59+xIbG0u5cuVo0sRucrd+/XqSk5Np06YN+/fvp1q1anmen5mZSePGjQH49ddf\nad68ORdeeCGPPvoocXFxRaqhW7duVKxYkdjYWJKSksjIyMAYw7x58wCYP38+vXr1yvPc6tWrY4zh\n8OHDnv3iL7wAu3fD++9DkHcjQYEjIoKdz1K3bl0++ugjwsPDCQkJAWDVqlVcfvnlVKtWjaysLILy\neVN+5plncr5euXIlnTp1AqBOnTpcdtllRaohNDSUbt26kZGRwYIFC5g7dy4NGjRgxYoVJCUlERMT\nQ58+ffI8Nzg4GGMMjuMU/ZeOioI5c2DRIqhYsejnFZMCR0TkNIsXL871pr5q1SratWsHQK1atUhM\nTCzw/LS0NOLj4+nYsWOx7t+nTx+MMUyfPp0ff/yRUaNGkZ6eTlRUFHv27KFFixZ5nnfo0CHKlStH\nzZo1i3ajZcts62bpUqhVq1i1ekqBIyKSLTExkf3799OqVauc104PnDp16nDo0KGzzsvIyGDFihUA\nrF27FiDnGgcPHiQ6OrrINfTo0YPzzjuP+Ph4evbsSc+ePQkJCWHUqFH5dqedqr127dpFu8mGDXDf\nffDvf0MRW18lQYEjIpKtQoUKhISE5IxQmzdvHvv3788JnEaNGhESEsK+fftynTdt2jS6detGSkoK\nixcvpkaNGpQrVw6AN998k+7du+ccu3DhQmrXrs2GDRvyrKFSpUrceuutgO3mq1q1Kp07dyYlJYXw\n8PB8a1+/fn2+rZ9cfvkFevaEN9+E7N/LV8r59G4iIn6sUqVKTJkyhRdeeIGLL76YtWvXUqtWLRo2\nbJhzTNeuXYmLi6Nfv345r3Xo0IHevXszZswYevfuTaVKlRg0aBChoaGEh4fnGmiQlZVFeno68+fP\nz3e1goiICDZu3JgTdPfccw87duygZcuW+dYeFxeXK9jydPiwHfb89NPQt29R/khKlGNHsEleHMcx\n+vMRfzNvHoSHz8OYvB8eS8np3bs3FStWzNUlFhsby+TJk5k/f/45XXvkyJGMGDHiXEsE7HOjpk2b\nsm7dOqpUqZL3QampNmyuugomTQJPBhd4yHEcjDFn3UBdaiIi2Z5//nm++eYbwD4T+eKLLxg4cGCu\nYzp16sTJkyfZtm1bse+zb98+wsLCzqnW082cOZOIiIj8wyYrCx58EKpVg4kTvRo2BVHgiIgABw4c\nYPz48SQlJQEwbNgwBg4cmNOtdbpT3W7F7QGJjo7Od3izp/744w+WLFnCkCFD8j9o6FD4+WeIjgYX\nV7hW4IiIYIc8R0ZGsmnTJv7xj3/QvHlzxo0bl+ex9erVY/Dgwbz11lse32fv3r1UrlyZWiU0FHnc\nuHG88847lC9fPu8DJk+282wWL/bJXJuC6BlOAfQMR/yRnuFIkc2bZwcIrFoF9ev77Lb5PcPRKDUR\nkdIoNhYGDoTPPvNp2BREXWoiIqXNxo122PPHH0Pz5m5Xk0OBIyJSmuzaBd2724mdN93kdjW5KHBE\nREqLAwfgllvguedcmdhZGAWOiEhpcOyYndgZHg5PPOF2NXlS4IiIBLrUVLtjZ8uW8PLLbleTLwWO\niEggy8yEe++F88+Ht95ybRWBotCwaBGRQGUMPPIIJCXBkiWuriJQFAocEZFAZAw8+yz88AN8+SVU\nqOB2RYVS4IiIBKIxY+yunV9/DZUru11NkShwREQCzZtvwvTpEBcH1au7XU2RKXBERALJe+/BK6/A\nypVw4YVuV+MRBY6ISKCYOxeGDYMVK/xmfTRPKHBERALB4sV2Qufnn0Pjxm5XUywKHBERf/fZZ3bH\nziVLoFkzt6spNk38FBHxZytWwF//CgsXQqtWbldzThQ4IiL+Ki4O+vWzz27atnW7mnOmwBER8Uer\nVsFdd8GHH0KnTm5XUyIUOCIByX/Xy5ISsGYN9O4Ns2dD585uV1NiFDgiAcaP12aUkrBmjV35edYs\n6NrV7WpKlAJHRMRfrFplw2bmTLuRWimjYdEiAUnNnFJn5Uro0weio6FLF7er8Qq1cEQCjO1SU+CU\nKsuX27D5+ONSGzagwBEJOHqGU8osXQp33w3z5sFNN7ldjVcpcEQCklKnVFiwAP72N7tsTYcOblfj\ndQockQATFAT6X7cUiI6GgQPtnjbXX+92NT6hQQMiAcYGjlo4Ae3tt+0GasuXw5VXul2NzyhwRAKM\nDRz/3rte8mEMjB1rN09buRIaNHC7Ip9S4IgEGHWpBaisLPjnP+3KzwG4eVpJUOCIBJhy5UD/6waY\n9HR44AH4+WcbNgG0LXRJ0t9akQCjwAkwJ05AeDgEB8MXX0BoqNsVuUbtcpEAExICEOJ2GVIU+/fD\njTdC7dowf36ZDhtQ4IgEnAoVACq4XYYUZts2u4fNbbfBjBlQvrzbFblO7XKRAHPeeQAV3S5DChIf\nD716wahRMGCA29X4DQWOSICpVAmgbHfN+LW5c+Hxx+GDD6BbN7er8SsKHJEAU6UKQBW3y5AznZpj\nM2UKfPklNGvmdkV+R4EjEmCqVgWoijFayNNvpKbCo4/Cxo2wdi3UqeN2RX5JgwZEAowdpZbCkSNu\nVyIA/PGHXeX52DGIi1PYFECBIxKQDrB/v9s1CBs2wHXXQefO9tmNfcAm+VDgiASkX/n1V7drKOM+\n+shulvbKKzBy5Kk1h6QAeoYjEpB2s3u32zWUURkZMGSI3ctGgwM8osARCUjb2LrV7RrKoD/+sLtz\nli8P69aV2TXRikttQJGAtJnNm92uoYxZtQpatID27e220Aobj6mFIxKQNvDdd2hotC9kZcHEifZZ\nzXvvaTLnOVDgiASk3YSEwE8/QePGbtdSih08CP3728/ffAP167tdUUBTl5pIgLrxRvvMWrxk5Upo\n3txuAR0Xp7ApAQockQDVqxfMm+d2FaVQejoMGwb9+sE779iuNK30XCIcY4zbNfgtx3GM/nzEHzmO\nQ0qK4YIL4IcfNLm9xGzbBhERULOmfV5Tu7bbFQUkx3Ewxpz1dFEtHJEAdd559r1xyhS3KykFsrLg\nzTehTRv4619hyRKFjReohVMAtXDEX2X/C5Jt26BdO9i5U6uqFNuePfDAA3D8uN1SQKMwzplaOCKl\nUKNGdvDA+PFuVxKAsrJg6lQ7t+amm+w8G4WNV6mFUwC1cMRfnWrhgG3dtGwJmzbpWU6Rbd8ODz0E\nJ0/a7Z+vvNLtikoVtXBESqkGDeCxx2DgQDsRVAqQmgqjR0Pr1nD77bB6tcLGhxQ4IqXAiy/aRxHT\np7tdiR/7+mu45ho7gfO77+Af/4DgYLerKlPUpVYAdamJvzq9S+2U//4XOnaE5cuhaVOXCvNHe/fC\ns8/ayZuvv24nMGk9IK9Sl5pIKXfllXZk7+2320WNy7yTJ2HcOLt9QIMGsGUL9O6tsHGR1lITKUX6\n9bPvq3fcYZe9qVzZ7YpcYAzExMBzz9mwiY+HSy91uypBXWoFUpea+Ku8utROMQYGDIAdO+z8xdBQ\nHxfnprg4GzQnT8KECdCpk9sVlUn5dakpcAqgwBF/VVDgAGRmwv33w4EDsHAhVKzow+LcsHEjPP+8\nfZA1apRdgkFbPrtGz3BEypDgYHj/fahRA267DY4edbsiL/nhBwgPh1tugVtvhYQEuzSNwsYv6b+K\nSClVrhzMmgVXXAE332y3dCk1Nm2yD6xuvhmuu872Hz7xBFSo4HZlUgAFjkgpFhwMb79tGwDt2tn3\n5YC2Zg307Gl/oRYt7C/0z39qIbkAoVFqIqWc49jJ9RddBO3bw6JFtlEQMDIzbdGvvQa//Wbn1MTE\n2OWyJaBo0EABNGhA/FVhgwbys3ixXRg5KgruvNMLhZWkpCS7J81bb0FYGAwaZIsup38n+7v8Bg3o\nv5xIGdKzJ3z6qZ2ns2sXPPWUH86D/O47u9NmTAx07w6zZ9t9aiTgKXBEypiWLe2jkO7d4eefYeJE\nP1hS7PBh+Phj2/Q6dMiu5LxlC/z5zy4XJiVJXWoFUJea+Kvidqmd7vBhu6xYrVowc6YLj0TS0+1y\nCB98AMuWQefOdsZqly4a1pyHxMREzj//fBy/a5KeTfNwRCSXatVs95oxdgqLT+bqZGXZ1QAGDrSb\n94wcCR062KbWvHl29JnC5iypqanUqVOHV1991e1Szon+y4qUYeedZ3uyrrjCvtd7JXTS0+Grr+Dx\nx+1QuVOf4+Ptx8CBUL26F25cerz33nsYY5gyZQqZmZlul1Ns6lIrgLrUxF+VRJfa6YyBv/8dvv/e\n9m796U/neMHERPj8c/jkE3vBhg3hrrvsas2XXVYiNZcVqampNGrUiKpVq+I4DsOGDeMvf/mL22UV\nyGdrqTmOUx/oYIyZWaIXLkGO4wwDXjXGpBVynAJH/FJJBw7Y0Bk40M6lXLrUw9HH6emwbp19JrNs\nmV1ypmNHOyyuRw+48MISrbUsmTp1Kp988glJSUn069ePadOmsXnzZoJdH+mRP588w3EcpzowGogu\nyet6QQwQ5XYRIv7EceCNN+znwYMLOTgjwwbM+PE2UGrWtGl19Ci8/LJdNXTxYnj4YYXNOdqxYwej\nR48G4LrrrqNFixb8EaAbHpVoC8dxnGnA28aYDSV2US9xHOcfwGFjzIwCjlELR/ySN1o4pyQlwfXX\nw5gxthcMsEPavvnGjqdeswa+/Rbq1bMP/Dt2hBtvtKEjXtOmTRsmTJhAmzPmJJ08eZLx48dz8OBB\nEhISCAoKYuzYsVx99dUuVeqDiZ+O41wO1A+EsMkWBax1HGeWMSbd7WJE/MX55Y8T9cTP3DegPrfN\ne4bQ71fBvn127bK2be1s0TZt7FLU4rqXX36Zxx57jIsuugiAF154gfbt2/Pdd99xqZ9tPFeSXWpP\nAB+U4PW8yhhzDFgH9HK7FhFXZGXZ5Qb+8x+IjIS+faFxY6hdm46zB9Cmxk9MODnQDldOSoLYWPjX\nv2wXmsLGL6SmpjJ58mTef//9nNeee+45UlJSmDx5snuF5aMkVxq4FRhfgtfzhTXAncBctwsR8ZrM\nTNi92+4V89//5v6oUgWuvtp+3HEHvPiiDZ3y5Rm0Du67D4a71zMjhcjMzKRmzZqcOHEi57XKlStT\no0YNtm/f7mJleStS4DiO0x24BrgBuP3U6C7Hcd4HVgFfANWMMT/nc/5TQCjQHBgB9AYc4HxjzKDi\nFu84zjPA08DFwG5gODAWqAP8H3AH8C7QBdgM9DDG/HbaJdYBQ4t7fxG/cugQ/PTT/z62brWft2+3\nywk0bgxXXgmtW8Pf/gZNmhQ4/6VFCzhyxF5CI5n9U2hoKDt37sz1WnJyMvv378/VnZacnMwLL7zA\n7NmzOXDgQK7jb7rpJr788kuf1Fto4DiOcz5wBTAGeBZoAqx3HCcIGxxzgLrA/nzOfxJYZozZ6jjO\nQ0AsNngeA+4FBmUf1wNoCqQAycaYdwqrzRgz0XGcX7AtlJ+MMdGO42RiR8lNM8b85jhOfyDOGNMi\nj0scAOo6jlNez3EkICQm2gABO0t/27b/fWRl2WRo1Mh+7tPHfr7sMqhc2eNbBQVB8+YKnEATHR1N\nxYoVefrppwFISUmhc+fONG7cmGXLlpGUlERERAQjR47k1ltvpaYPB3sUpYXTGfuG3gnIAH7Mfr05\nttWyFrgZOJzP+cHGmK3ZX18ErDfG7HMcZyowC8BxnCrAKGPMtdnff+44zqfGmD1FqG8pNqQ6ZYdj\nOWzrqQ+2ddMbWJDPuYnZx1bDhs9ZnE6nDbSoDzSAER1H8FKnl8469qXYlxj59cizXtfxOv6cjk+8\nmpe+CbWhkp5uAwXs11262OHIjRrZUWIlvM5WACzbJac5fPgwo0ePZsqUKVxyySUADBs2jAsvvJD3\n3nsv57iIiAgSEhJ4+OGHS+S+sbGxxMbGFnpckYdFO44TBWQaYx7N/v4pYIAx5irHccKBp40x7Qq5\nxlfA58aYsWe83hO42xhzT/b3Y4CfjTFFmivjOE4MNlgGYJ/JNMGG25+xc26eM8Z8l8d5IcBJIMwY\nc9YGvBoWLV5jDOzda5+jJCTYlZETEmxzIikJLr3UdoE1apT7IywMHMerw6JPd911dt+zDh28fisp\nRH7Dok8xxtCnTx969uxJ//79ARtAF1xwAd9//z1XXHFFzrFPPvkkjuMwadIkr9RaEsOie2K7wE5p\nD6zO/voAUOBiSNlv7q2BYXn8uA65W0hHAE8a8fOAu4AHsSHzIrb1NACom1fYZKsBZOQVNiIl5sAB\n2LQJNm+2M/B/+MEGTaVKdhGzK66wz1N697Yhc9FFfrGA5e+/20ZV69ZuVyJFMXToUO6991569bID\nb3fs2MGWLVsICwvLFTYA8fHxOV1uvlTUQQPVgTDsQ/ZT2gPPZX+9F/vmfeZ55YAbjDErgFOxvC77\nZzWBW4wx0cD52JbGKWlAlaL/Gvwn+/zWwBvA4uxrvAi8VcB51YHAnLIr/scYOxrs++/tJmLr18OG\nDZCcDM2a2ZFgLVtC//724b2fL1j50Ud2FemQELcrkcJMmzaNNm3acMcdd+S8NnPmTK6++mouvvji\nXMeuX7+e33//nbtyZvX6TlFbOKnYN/AQAMdx+mADaDWAMWab4zhpjuNcaIzZd9p5DwPjswOrJ3DI\nGJOR/bO/A69nf30s+3qnVMQ+XyH7fncC72AD6qyJpcaYE47jLMOOSosxxhx1HOdL4DZsl1p+mgP5\ntX5ECnbkiJ19Hx9vZ95/+63dyaxlS7j2WrusS/PmULduwD0MOXjQrjTw1VduVyKF+eSTT5g7dy5d\nunQhISEBgBMnTlCxYkU6duzI0KFDyczMJDg4mGPHjvH4448ze/ZsKlSo4PNaixQ42W/ojwEvZ48K\nawMcMMbsOO2wz7HDpuec9tpKYD526PF84ITjOOOBZGwwnOpG24FtnZxSAzh9EHkQUB77nCa/lQyi\ngWbGmFPdfB8CDY0x/1fAr3YDsKSAn4v8z759sHKl/YiLg5077djhNm3gwQfttsh16rhdZYl49lm4\n+2646iq3K5GCJCYmcs8995CSksKKFSty/WzOnDnUqlWLSZMm8cgjj1CnTh327dvH5MmTadmypSv1\nFmstNcdx5gMpxpiI017rBDxpjOldjOuFAuuMMU2yv18LRJw5r8dxnBHGmLOH8RRD9jOlTUCr7FUH\n8jpGgwbKsqQkWL7c/jN/+XL7z/4bbrBP0Nu3t62X8uVdKc2bgwamToXJk23DrWpVr9xCPPDVV1/R\nqlUrunbtyoQJEzh+/Dht27alUqVKbpeWr3MaNOA4zr+ARcaYb7K7x7pgVxbIYYyJdRznWcdxGhlj\ntnlSnDEm2XGcfzmO8wK2++79PMLmQvKZ61NM9wHR+YWNlEFZWfb5y9Kldon9zZttsHTubLvHmjb1\ni4f53rRiBYwYAatXK2z8RUxMTE7r5cCBA/Tv35/t27f7deDkp9AWjuM4tYBfgTuMMcscx5kCHDXG\nPJfHsfWAcdghziX6zy/Hcf6JDaI858t4eK3awFSgb0ETPtXCKQNSU20LZuFCu6ZY5crQvTvcdptt\nzZx3ntsV5skbLZy4ODtQbu5cu/iz+IedO3fSsmVLLrnkEurXr0/9+vX9fqvpc9qAzXGcwdjnKH8G\nEowx0wo4tiXQ2hjz5jnUe+Y162CX1JlSQtebAIw1xhTYYlLglFKpqbYFExMDS5bYIcm9etnNwgJk\nSn1JB05cnN2K4MMPbYNO/MuAAQNYtGgRJ0+eZPv27YSFhRV+kot8tuNnaaLAKUWysmx/UXS0bc1c\nfbVdHbl3b7jgArer81hJBs6SJXak9kcfKWz81c6dO2nYsCF333030dH+vr+lAqdYFDilwLZt8P77\nMGuWXVL/3nuhXz87uTKAlVTgvPcePP+8zeDrry+BwsRrhg4dyqOPPkq9evXcLqVQCpxiUOAEqJMn\nYcECmDbNzui/9164/3770L+UONfAMQZGj4YZM2zvYuPGJViclHle3/FTxHW7d9sxvdOn25n9Awfa\nPV40VT6XtDQ76O6HH+xu0QHYoygBqnSP8ZTSzxg7hveuu+zs/pMnYdUq+OILCA9X2Jzh8GG7XE1S\nEnz9tcJGfEuBI4EpMxP+/W87y/++++w43t27YeLEgBlp5mt79thpRVddBfPn27VDvWHXrl3MnDnT\nOxcPUJGRkaSlpbldhusUOBJY0tPtIIAmTeCVV+waLD/9BH//e7E2GSsrNm6Edu3ggQdg0iS75Js3\nJCYmMnz4cCIiIgo/GOjTp493CvEz4eHhDBgwwO0yXKdnOBIYTgVNZKTdK+btt22rJsAWxXTDqlW2\nx/GNN+xIcG8aMmQIgwcPJjiPRIuLi6NWrVpcfvnlOa8lJiaedVyg2717N9HR0TiOwx9//EHDhg15\n4oknaNasGTNmzOCBBx5wu0T3GGP0kc+H/eMRV6WnGzNjhjENGhjTubMxq1e7XZFfKOrfzU8/NaZm\nTWM++8zLBRljtmzZYrp06ZLvzxs2bGhmzZqV67Ubb7zR22X51A8//GDefffdXK/VrFnTJCQkmKNH\nj5omTZqYtLQ0l6rzney/n2e9p6pLTfyTMXZoc9OmtmUzc6YdCNC2rduVBYylS+3jrUWLoGtX79/v\njTfe4P7778/zZ7/99hs7d+6kY8eO3i/EJSdPniQ+Pp4HH3ww57VDhw5x+PBhKlWqRJUqVWjVqhUL\nFuS3433pp8AR/7NmjQ2WkSNh/HiIjbVPu6XIli+3qwd88onvMnrZsmX5bn8cFxfHxRdffNZmYKVJ\nTEwMd999d67XRo4cyUMPPcRF2RON27Zty8KFC90ozy/oGY74j507YcgQGziRkXbCZilfndkb/u//\n7F428+aVzPbQ48aNY9OmTdx8883UrVuXb7/9lu3bt9OsWTOeeuopwD63OHz4MJdcckmucz/55BNi\nYmJYtWoVoaGh3HfffbRt25ZHH3003/tNmjSJ5ORk1q9fz8iRI5k/fz7GGJKSkhg/fnyxf4+JEyfy\n+uuv88svv1CvXj1Gjx7NkCFD2Lt3Ly1btmTRokU89NBDfPHFF1x99dX85z//4QIPxo2npqYSGhrK\n2LFjOXLkCIsWLeKiiy7i888/zzmmVatWjBkzpti/Q8DLq59NH3qG41MnThgzfLgx1asbM2qU/V4K\nlN/fzV9+MaZOHWMWLCiZ+2zYsMEsWrTILFmyxFStWtXMmTPHGGNMenq6qVq1qtm8ebMxxpiVK1ea\nyy+/PN/rXHPNNWbatGlnvX7mM5xJkyaZhIQEY4wxUVFRJiwszOzdu9cMHz7chIWFGWOMWbx4sYmM\njDQTJkwwU6dO9ej3iYmJMY7jmK5duxpjjPnoo49MUFCQiYqKMsYYs3//ftO4cWOPrmmMMbt37zbL\nli0zxhgzffp0M2zYMPPII4+YcuXKmW+//TbnuF9//dUEBweX+uc45PMMx/U3dX/+UOB4WVaWMfPm\nGVO3rjF/+Yt9t5QiyevvZkqKMddea8zYsSV3n3fffdckJyebCRMmnBUOYWFh5u233zbGGDNv3jzT\nunXrPK9x9OhRExwcbP773/+e9bMzrzlhwoScr0eMGGFuueUWY4x9o966das5evSoad68ec4xXbp0\nMbt37y7y73PixAkTGhpqQkJCTGJiopk1a5ZxHCfnPlOnTjVDhgwp8vVO+fjjj01ycvJZr1977bXm\nhRdeyPk+OTnZBAUFmf3793t8j0CSX+Cov0LcsXOn3XfmxRfhgw/sUsUBvqCm24YNgwYN7NSkkvLg\ngw9SsWJFVq9eneuB//79+zlw4ABVqlQBICsri6B8uj/XrFnD+eefzxVXXFHo/Z555pmcr1euXEmn\nTp0AqFOnDpdddhmxsbG5hlW3aNGCzz77rMi/T2hoKN26dSMjI4MFCxYwd+5cGjRowIoVK0hKSiIm\nJqZYc4PWmA4BAAAMgElEQVTS0tKoWLHiWa8fO3aMrKysnO+Dg4MxxuCU0eH8ChzxrfR0GDcOWrWy\nG5ytXw/ZbypSfMuXw5w58M473pmatGbNmlyBExsbS1BQEB06dACgVq1a+c6piYuLo127dh7dLy0t\njfj4+LNGte3du5dq1arlfP+nP/2Jn376yaNr9+nTB2MM06dP58cff2TUqFGkp6cTFRXFnj17aNGi\nhUfXA/v85kwHDx5k586dtD9twMuhQ4coV64cNWvW9PgepYECR3xn/Xq7Bv5XX8G338LQoVrrrASk\npcFjj9l1S2vUKPnr79y5k6SkpFwj0ObMmUOvXr2oW7cuYFsghw4dyvP81atXc8MNNwCwZcuWfIcF\nZ2Rk5GylvHbtWsA+ZAf75j179mySkpI477RdWENCQjh2zLNd4nv06MF5551HfHw8PXv2pGfPnoSE\nhDBq1Ch69erl0bUAfv/9d7Zt23bW61FRUVx11VXceuutOa8lJiZSu3Ztj+9RWihwxPtSU2H4cLjl\nFnjySbse/hmjmaT43n4bGjaEHj28c/1Vq1aRlZWVEyjLli1j/fr1TJ48OeeYRo0aERISwr59+846\n/+DBg1x55ZUYY5gxYwa33357nveZNm0a3bp1IyUlhcWLF1OjRg3KlbMDad9880169OhBlSpVTj1f\nBSAlJYXq1avnus7ChQupXbs2GzZsyPM+lSpVygmB8PBwqlatSufOnUlJSSE8PPys4wu73tq1azl+\n/DhJSUk5r8XHxxMVFcWcOXNyHbt+/fpitaBKCw2LFu/asMHOPrzkErugl5YnLlHJyfCvf9nNTL1l\nzZo19O/fn9dff53Q0FD27NnDypUrzxoy3LVrV+Li4ujXr1+u14cMGcKHH35IXFwcAwcOzHPZG4AO\nHTrQu3dvxowZQ+/evalUqRKDBg0iNDSU8PBwqlWrRsOGDYmPj88559ChQ1x66aW5rpOVlUV6ejrz\n58/nmmuuyfNeERERbNy4Maer75577mHHjh20bNnyrGMLu97hw4cZN24ckZGRlC9fnoyMDI4fP87q\n1avP+jOKi4uje/fuedZUJuQ1kkAfGqV2ztLTjXn5ZWNq1TLmgw/siDQpMaf+bk6dasztt3v3Xk2b\nNjXR0dGFHrdixQrTq1cvj67t6dI2J06cMFdeeWXO961btzY7duzI89iXXnrJo2sXJr/rTZkypUjn\np6ammsaNG5ujR4+WZFl+iXxGqamFIyVv5047abNiRfjuOyjFs8vdZIxdkPONN7x3j6NHj/Ljjz9y\n3XXXFXpsp06deOWVV9i2bRuNGjXySj2hoaE8//zzvPzyy1SoUIH+/fufNdkUYN++fYSFhZXYffO7\nXmJiYs5IvcLMnDmTiIiIIh9fKuWVQvpQC6fYZs2yq0W+9poxmZluV1NqAWbTJjuFyZuNx08//dRU\nq1atyMfv2rXL9OvXz2QVsahOnToVt7QCvfLKKyU61yW/6y1cuND89NNPhZ7/+++/mzvvvLPUT/g8\nBc3DEa86ftw+q4mMtItsDhqkZWm8LCbGbjfgrSkdH374IU8++SSpqalERERw4sSJQs+pV68egwcP\n5q233irSPc584F8S9u7dS+XKlalVq5bXr7d169YitebGjRvHO++8Q/ny5UukpkDl2DCSvDiOY/Tn\nUwQbN9p3vnbtbP+Ot7aSlByO43D99YaxYzWNSfyP4zgYY876p5ACpwAKnEIYA+++C88/b7d2vvde\ntysqMxynCpUqHePAAfuoTMSf5Bc4GjQgxZOcDAMHwrp1EBcHpy03Ir7QnKuuUthIYFEnu3huxw5o\n08YuU/PNNwobV1xN06Zu1yDiGQWOeOazz+yOXgMGwOzZULmy2xWVUVfQpInbNYh4Rl1qUjTGwCuv\nwKRJdmev7LWxxC31adDA7RpEPKPAkcKlpMADD9iutG+/1TYCfqEe2etmigQMdalJwfbtgw4d7Jya\nr79W2PiNMP78Z7drEPGMAkfy9913djuB3r3t8xoNifILdj+v6l7ZikDEm9SlJnlbtAgeesju6NW7\nt9vVyGmOHwdIKfOz1iXwKHAkN2Pg9dfhtddg6VK7M6f4laNHAY4BVV2uRMQzChz5n8xMePppiI2F\nNWugXj23K5I82BbOcbfLEPGYAkeslBS45x44cgRWrYI//cntiiQfKSkAKW6XIeIxDRoQOHQIbr7Z\nLrq5bJnCxs+dPAlw0u0yRDymwCnr9uyB9u3t0OeZMyEkxO2KpBBpaQBpbpch4jEFTln24482bB55\nBMaO1f41ASI9HSDd7TJEPKZnOGXV2rVw550wYQJERLhdjXggIwMgw+0yRDymwCmLvvjChswHH8Bt\nt7ldjXgoMxMg0+0yRDymPpSyZv58Gzbz5ytsApRdaSDL7TJEPKbAKUtmzoTHH7cj0dq3d7saKSa7\nCa0CRwKPutTKiqlTITISVqzQhmkBzrZwtPW5BB4FTlnw+ut2H5vYWGjY0O1q5BzZFo4CRwKPAqe0\nGzsWpk+3WwtoA5VSQYEjgUqBU5q9/DJER9uwufBCt6uREqLAkUClwCmNjIGRI2HuXNuNpp26SiEF\njgQeBU5pYwyMGGGHPcfGQliY2xWJiAAaFl36jBxpw+arrxQ2pZRR40YClFo4pcnIkRATY4c+K2xK\nOaWOBB4FTmkRGQlz5ihsRMRvKXBKg9des6sIxMZC7dpuVyMikicFTqB7802YMsUOfb7gArerERHJ\nlwInkEVFwauv2rC56CK3qxERKZACJ1BFR9tBArGxUL++29WIiBRKgROIFi6EQYNg+XK49FK3qxER\nKRIFTqD54gu7JfTSpdCkidvViIgUmQInkKxZYzdPW7AAWrRwuxoREY9opYFAsWkT9OoFs2ZBu3Zu\nVyMi4jEFTiDYts1uB/3GG3DLLW5XIyJSLAocf7dvnw2Zl16Cvn3drkZEpNgUOP4sKcmGzcMPw4AB\nblcjInJOFDj+KjkZevSALl3guefcrkZE5JwpcPxRejqEh0PDhnadNMdxuyIRkXOmwPE3xvyv+2z6\ndAjSfyIRKR00D8ffDB0KCQl2FYHy5d2uRkSkxChw/Mnrr8OiRbBqFVSq5HY1IiIlSoHjL+bMgfHj\nbdjUqOF2NSIiJU6B4w9iY+GJJ+DLL6FePberERHxCj2RdtvmzdCvH3z8MTRt6nY1IiJeo8Bx06+/\nQvfu9tnNTTe5XY2IiFcpcNxy5Ah062a70u6+2+1qRES8ToHjhrQ0uOsu6NABBg92uxoREZ9Q4Pia\nMfDQQ1C5MkyapFUERKTM0Cg1XxsxArZuhRUrIDjY7WpERHxGgeNLM2ZAdDSsXQuhoW5XIyLiUwoc\nX/n8c7tszcqVEBbmdjUiIj6nwPGFzZvh3nvh3/+Gxo3drkZExBUaNOBt+/bZfW0mTYIbbnC7GhER\n1yhwvOn4cejZ0+7Yqbk2IlLGKXC8JTMT7rkHrrkGnn/e7WpERFynwPGWjAxo3RqmTtVcGxERNGjA\neypUUMtGROQ0auGIiIhPKHBERMQnFDgiIuITChwREfEJBY6IiPiEAkdERHxCgSMiIj6hwBEREZ9Q\n4IiIiE8ocERExCcUOCIi4hMKHBER8QkFjoiI+IQCR0REfEKBIxKQfnS7ABGPKXBEApICRwKPAkdE\nRHxCgSMiIj6hLaYL4TiO2yWI5El/NyXQOMYYt2sQEZEyQF1qIiLiEwocERHxCQWOiIj4hAJHRER8\nQoEjIiI+ocARERGfUOCIiIhPaOKniB9zHOcpIBRoDowAegMOcL4xZpCbtYl4ShM/RfyU4zhPAp8Z\nY7Y6jvMQEIkNnseAh40xtR3H6QE0BVKAZGPMO+5VLFIwtXBE/FewMWZr9tcXAeuNMfscx5kKzHIc\npwowyhhzLYDjOJ87jvOpMWaPWwWLFETPcET8lDFm4mnfdgBis1/fa4z5CegEJJx2zHfALb6qT8RT\nChwRP+c4TgjQGvj6jB/VAQ6f9v0R4DJf1SXiKQWOiB9yHKec4zg3Zn/bJvvzuuyf1XQc517gfODk\naaelAVV8V6WIZxQ4Iv7pYWCp4zgVgZ7AIWNMRvbP/g78BziGHbF2SkUg0adVinhAgwZE/NNKYD4w\nNPvzCcdxxgPJQIwx5rDjODuwXW2n1AC2+7xSkSLSsGiRAOU4TiiwzhjTJPv7tUCEMeZndysTyZsC\nRySAOY4TAVwCpAJHNA9H/JkCR0REfEKDBkRExCcUOCIi4hMKHBER8QkFjoiI+IQCR0REfEKBIyIi\nPqHAERERn1DgiIiITyhwRETEJ/4fzxDyYrBoeKUAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "xx = np.linspace(-0.9, 0.9, 100)\n", "yy = 4 * xx - np.sin(xx * np.pi)\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(xx, yy, color=\"red\")\n", "\n", "ax.set_xlim(-1, 1)\n", "ax.set_ylim(-4, 4)\n", "\n", "ax.set_xticks([0])\n", "ax.set_xticklabels([r'$x_0$'], fontsize=\"xx-large\")\n", "ax.set_yticks([0])\n", "ax.set_yticklabels([r'$y(x_0, \\mathbf{w})$'], fontsize=\"xx-large\")\n", "\n", "xx = np.linspace(-4, 4, 100)\n", "yy = norm.pdf(xx, scale=0.5) / 5\n", "\n", "ax.plot([-1, 0], [0, 0], \"g--\")\n", "ax.plot([0, 0], [-4, 4], \"k\")\n", "ax.plot(yy, xx)\n", "\n", "ax.annotate(\"\",\n", " xy=(0.75, -0.5), xycoords='data',\n", " xytext=(0.75, 0.5), textcoords='data',\n", " arrowprops=dict(arrowstyle=\"<->\",\n", " connectionstyle=\"arc3\"), \n", " )\n", "\n", "ax.text(0.77, -0.2, r'$2\\sigma$', fontsize=\"xx-large\")\n", "ax.text(0.15, -1, r'$p(t|x_0,\\mathbf{w}, \\beta)$', fontsize=\"xx-large\")\n", "ax.text(0.5, 3, r'$y(x, \\mathbf{w})$', fontsize=\"xx-large\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 最大似然" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "设训练集数据是独立同分布的，那么似然函数为：\n", "\n", "$$\n", "p(\\mathsf t \\left|\\mathsf x, \\mathbf w, \\beta\\right.)\n", "= \\sum_{i=1}^N \\mathcal{N}\\left(t_n\\left|~y(x,\\mathbf w), \\beta^{-1}\\right.\\right)\n", "$$\n", "\n", "对数似然为：\n", "\n", "$$\n", "\\ln p(\\mathsf t \\left|\\mathsf x, \\mathbf w, \\beta\\right.) = -\\frac{\\beta}{2} \\sum_{n=1}^N \\{y(x,\\mathbf w) - t_n\\}^2 + \\frac N 2 \\ln \\beta - \\frac N 2 \\ln(2\\pi) \n", "$$\n", "\n", "设系数的最大似然解为 $\\mathbf w_{ML}$，从最大化对数似然的角度来看，求它的问题相当于最小化：\n", "\n", "$$\n", "\\frac{1}{2} \\sum_{n=1}^N \\{y(x,\\mathbf w) - t_n\\}^2\n", "$$\n", "\n", "这就是之前最小化平方误差和的结果。\n", "\n", "因此最小化平方误差和可以看成是高斯噪音假设下的最大似然的结果。\n", "\n", "再对精度 $\\beta$ 求最大似然，我们有（可以理解为照搬之前求 $\\sigma^2$ 的结果）：\n", "\n", "$$\n", "\\frac{1}{\\beta_{ML}} = \\frac{1}{N}\\sum_{i=1}^N \\{y(x,\\mathbf w) - t_n\\}^2\n", "$$\n", "\n", "我们有了最大似然的结果之后，对于一个新的输入 $x$，其输出 $t$ 应当满足：\n", "\n", "$$\n", "p\\left(t\\left|~x,\\mathbf w_{ML}, \\beta_{ML}\\right.\\right) = \\mathcal N\\left(t\\left|~y(x,\\mathbf w_{ML}), \\beta_{ML}^{-1}\\right.\\right)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 最大化后验分布" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "假设我们对系数 $\\mathbf w$ 有一个先验的知识（$M$ 是多项式阶数，加上常数项一共 $M+1$ 维）：\n", "\n", "$$\n", "p(\\mathbf w~|~\\alpha) = \\mathcal{N}(\\mathbf w~|~\\mathbf 0, \\alpha^{-1} I)\n", "= \\left(\\frac{\\alpha}{2\\pi}\\right)^{(M+1)/2} \\exp\\left\\{-\\frac{\\alpha}{2}\\mathbf{w}^\\top\\mathbf{w}\\right\\}\n", "$$\n", "\n", "$\\alpha$ 控制这个模型的先验分布，这一类的参数通常被叫做超参（`hyperparameters`），Bayes 公式告诉我们，后验概率正比与先验概率和似然函数的乘积：\n", "\n", "$$\n", "p(\\mathbf w~|~\\mathsf{x, t}, \\alpha, \\beta) \\propto p(\\mathsf t~|~\\mathsf x, \\mathbf w, \\beta)~p(\\mathbf w~|~\\alpha) \n", "$$\n", "\n", "我们可以通过最大化后验概率（`maximum posterior, MAP`）来决定参数 $\\mathbf w$ 的值，对上式求对数，并去掉跟 $\\mathbf w$ 无关的项，我们相当于要最大化：\n", "\n", "$$\n", "-\\frac{\\beta}{2} \\sum_{n=1}^N \\{y(x,\\mathbf w) - t_n\\}^2 -\\frac{\\alpha}{2}\\mathbf{w}^\\top\\mathbf{w}\n", "$$\n", "\n", "即最小化\n", "\n", "$$\n", "\\frac{\\beta}{2} \\sum_{n=1}^N \\{y(x,\\mathbf w) - t_n\\}^2 +\\frac{\\alpha}{2}\\mathbf{w}^\\top\\mathbf{w}\n", "$$\n", "\n", "因此，`MAP` 的结果相当于给多项式拟合加二范数正则化的结果，其中正则参数 $\\lambda = \\alpha / \\beta$。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.2.6 Bayes 曲线拟合" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "虽然在 `MAP` 中，我们引入了先验分布，但是本质上它还是一个点估计，本质上并不是一个完全的 `Bayes` 估计。\n", "\n", "一个完全的 `Bayes` 估计要求我们对 $\\mathbf w$ 的所有值进行积分。\n", "\n", "对于之前的曲线拟合问题，给定训练集 $\\mathsf x$ 和 $\\mathsf t$，对于一个新的测试样例 $x$，其目标值为 $t$，我们考虑预测的分布 $p(t~|~x,\\mathsf{x, t})$（这里我们假定 $\\beta, \\alpha$ 两个参数已经给定了）\n", "\n", "Bayes 公式给出：\n", "\n", "$$\n", "p(t~|~x,\\mathsf{x, t}) = \\int p(t~|~x,\\mathbf w) p(w~|~\\mathsf{x,t})d\\mathbf w\n", "$$\n", "\n", "其中 $p(t~|~x,\\mathbf w)$ 是之前给定的高斯分布，$p(w~|~\\mathsf{x,t})$ 是训练集上的后验概率（也是一个高斯分布）。\n", "\n", "由于高斯分布的性质，上面的式子本质上也是一个高斯分布，因此可以写成\n", "\n", "$$\n", "p(t~|~x,\\mathsf{x, t})=\\mathcal{N}\\left(t~|~m(x), s^2(x)\\right)\n", "$$\n", "\n", "其中均值和方差分别为\n", "\n", "$$\n", "\\begin{align}\n", "m(x) & = \\beta \\phi(x)^\\text{T} \\mathbf S \\sum_{n=1}^N \\phi(x_n) t_n \\\\\n", "s^2(x) & = \\beta^{-1} + \\phi(x)^\\top \\mathbf S \\phi(x) \n", "\\end{align}\n", "$$\n", "\n", "其中，$\\phi_i(x) = x^i, i = 0, \\dots, M$，矩阵 $\\mathbf S$：\n", "\n", "$$\n", "\\mathbf S^{-1}=\\alpha I +\\beta \\sum_{n=1}^N \\phi(x_n)^\\top\\phi(x_n) \n", "$$\n", "\n", "下图粉红色部分就是 Bayes 估计给出的结果，红色曲线是 `MAP` 给出的结果。" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEOCAYAAACEiBAqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVPX+x/HXOeyLgoCKgIqC4q6V5a1+FdX1Wt1svXVb\n7Gpllpr7WrklprmklVu2uLSY2Z7d9sXq2s30trgrIOACLiiK7DBzfn98RUUZGJgzMMN8no+Hj4Rz\n5szXxPOe810+X80wDIQQQnguvb4bIIQQon5JEAghhIeTIBBCCA8nQSCEEB5OgkAIITycBIEQQng4\n7/puQG1omiZzXoUQohYMw9DO/55bBgFAXa5/mDZtGtOmTauz9xNCeIa6vrdo2gUZAEjXkBBCeDwJ\nAiGE8HASBHZITEys7yYIIRogV7m3aO5Ya0jTNMMd2y2EEPVJ07RKB4vliUAIITycBIEQQng4CQIh\nhPBwEgRCCOHhJAiEEMLDSRAIIYSHkyAQQggPJ0EghBAeToJACCE8nASBEEJ4OAkCIYTwcG67H4Gw\nU3EJnCqAvALIL4DCYigtA4sVyus1aRp46+DjA/5+EBwAwYHQKAh8feq3/UIIp5Oicw2N1QrHc+Ho\nccjJhTIL6Jq68deElw5WA7y9ICwEIppAk8bq+0IIt2Sr6JwEQUNgGJCbBwePwLET6hN+TW/81fHS\n1fuEhUJ0UwhppN5HCOE2JAgaIqsVjhyHjEwoKVNf1wUvHby8oFUkREao3wshXJ4EQUNitUJWNqQf\nVN03dRUA59N10ICWkRDTXAJBCBcnQdAQGIbq+knep/r+6ysAzqfrqpuoTTRENZUuIyFclASBuyso\ngt1pkFfoOgFwPl0HX29IaAOhjeq7NUKI80gQuCvDgH1Z6pfVTf7Mug4RoRDfCnxkhrIQrkKCwB0V\nFMH2FCgqcd2nAFs0TQ0qd2gL4SH13RohBBIE7udQthoLcLcAOJ+uQ7Mw9XQgaxCEqFcSBO7CYoXd\n6WpQ2N1DoJyuqRXKXdtBYEB9t0YIjyVB4A6KimFLMhQXu894QE3oOiTEqicEIUSdkyBwdSfzYGsy\nWCz13RLn0jWIagZtY2SaqRB1TILAlR05BrszGk5XUHV0HUKCoXOcLEITog5JELiq/VmQnuVSIZCW\nlcnkL9dx0GohWvciqU9f2rSIMvdNNA38faF7Avj5mnttIUSl3CYINE17DbgZOGwYRjcb57h/EBgG\npO5XpSJcLAR6f7ia1AH9ISAACguJW7mKr2+/zzlh4O0FPRJkEFmIOmArCFxxPt8KoE99N8KpDEPN\nDHKxEACY/OW6syEAEBBA6oD+TP5ynflvZhhqb4TfdkJuvvnXF0LYxeWCwDCM/wA59d0OpzEM2LkX\njua4XAgAHLRazoZAuYAAMq1OHMS2WOHP3XDilPPeQwhhk8sFQYNmGLAjFY6ddMkQAIjWvaCwsOI3\nCwuJ0p08qGu1qllTx086932EEBdwuTECAE3TWgPrGtQYQfmTgAuHAFQyRpB/kpi1ixl2RQfy/QrJ\nLDnKybJ8yowyNDRCvRsR4RNKXEA0CYGtuSg4gcbewbVvgK5DpzgpSyGEE9gaI3DbimDTpk078/vE\nxEQSExPrrS3VMgxVOfTYCZdfKNamRRTL/n4NY76cTIb/QU75ZRPUIZqthjdxxNCzUSdCvYPx0byx\nGgYnyk5xpDSHX3K3sfLQp2zJS6F9YCt6N+nF7RGJXNa4M7pWgwdPqxV2pEDneLVFphCi1tavX8/6\n9eurPc9VnwhiUU8EXW0cd58nAsOAlH1w6JhLPwnklRWw6vCnvH7oMzKKsri7WW9uCLucq0MuItg7\n0O7rFFtL2HxqJ58f+5kPsr8n31LIo1F38HCLW2juG25/g3QdusSrfZKFEKZwp+mjq4FEIBw4DEw1\nDGPFeee4TxCkZ8L+Qy4bAkdKjvPigXdYlvUBiaGXMLDFrVwfeineujkPi7+d2sXSzPd47+h33Bp+\nNZNjBxIXEGPfi3UdurVXi8+EEA5zmyCwh9sEQeYRtVbABbuDCi1FPLf/LRYceJu7m/2VMTH3Ex/Y\n0mnvl1OaywsH1rDw4FrubHot02IHEeXXtPoXeunQowME2/9UIoSonARBXcs+ATtTax0CzlrdaxgG\n7x39ljGpz9OrcRdmt32ctvZ+QjfBsdITzN73OsuzPmFy7MMMjbqr+qcPby+4uCME+NdNI4VooCQI\n6lJuHvy5p9bdQaav7j15AjLSOLlvNx+lfkxJzhFuCr6UaL+mqvslKBiaNIGwcGjZGmJagbdz5xHs\nzE9jSPJsTpblsarDNLoGx1f9Al8fuKST+q8QolYkCOpKUTFs3uFQFdF+K5fx1t13VFzYVVjI/Ws/\n4M0Bj1b9YsOAtFT4bRP88Rts3wIlxRyLjuC7wCxaNIunV8z/4RMQpM63WiE/D3KOw7Fs2JcBRw5D\ny1bQ7SLofjFcdAmEhNb6z2O7qQYrD61j/N6FPNXqQYbH3FP1DKNAf/VkIIXqhKiVBjd91CWVlcEf\nux0uJV2r1b17U+C7r+D7b9T7X9oLEq+nZOhwJua9ywfZ61nT6SX+ElLpRKyKiosgNQW2/A5ffApz\nk6BDZ7jmerj6Ogg1JxQ0TePBFrdwdejF9Ns5hS+O/5c3OyYR4Wvj+oXFsC1FDSBLCWshTCNBYBbD\nUDepklKHL3Vmde95TwQXrO4tK4OfvocP34XMA9D7RpjyDLTvAJpGZvFR7tw+gQifEH7r+QZhPnbO\ny/fzh05d1K97HoCiIvj1Z1j/Lby8CP5yJdxyJ3TtbsoNOS4ghp96vMyTaUu49Lf+fNRlLt2D2194\nomGomkR7MqB9awkDIUwiXUNm2ZMBh81ZK1DtGEFZGXz5b3hjOTRvDrffDf+XWKFf//dTu7l12xge\njbqDJ1oNqNmirqrknlTv/cn7amzhgYfgiqtNuymvOfwVw1LmsrjdeO5u1rvyk3QdYqOgZaQp7ymE\np5AxAmfKOgop+01dK1A+ayjTaiGqfNZQZAv44fSn8sgWMOBR6Nbjgteuy/6Rh3YnsaTdBO5q9lfT\n2lSB1Qo/rYc3XlOf1Ac9Dr2uMOXSf5zaTd9toxkZcy9jWvar/CQpRSFEjUkQOEtunqqc6ey1Aul7\n4cW5cOIEPD4aLr600tNWZq3jibTFfNRlHr0ad3Fum0CFwIYf4aUXISoahoyE2LYOX3Z/0SFu2DKc\nG8IuZ27ciMqfaHRdDR4HyV4GQthDgsAZiktg83Yoc2KJ5rIyeHMFfLgW+g9UffM2pnYuOrCWOftf\n56tui+gQFOu8NlWmtBQ+ele19W83wcODwd+xef/HS09y67axtPGPYkWHKXhplcwW8vWBnp3BR4a7\nhKiOBIHZrFa1oUpBkfpU7Azpe2HmNDVLZ9wkaNrM5qmzMlbwWtYnfNN9MbEBJu8kVhMnTsCiebBj\nm2rzRT0dulyBpYhbto6muW8YqzpMu3DxmaZBo0C1+lgGj4WokgSB2fakw+Hjzqsh9NknsGwhDBwM\nN99e5U0uKf1V3j7yFd90X2xf2Ya68PNPsOBZuPwq1V3kwNNBoaWIW7eNJdwnhDc6PH1hGOg6REZA\nu1YONlqIhk2CwEyHj6lZQs4IgaIieGGO+kT99LPV9rfP3/8WL2V+wI89lhHpF2F+exyRlwfPPwsp\nyTBtlkNjB4WWIm7bNo4m3o14q1PShd1Eug4JraFZDSqcCuFh3GnPYteWX+i8EMg+CsMfgZISWLqy\n2hvny5kfsPDgWr7tvtj1QgAgOBieSoK774cRj8K/P651N1qAlz8fd5nH0dIchuyZzQUfBKxW2J2h\n/n6EEDUiTwQ1YbHApu1qkNhse1PgiVHQ9w64f0C1/d2rD3/B+NSFrO/xklOrhpomfS9MmwhdusOI\n8eBTec2g6ortnSrL59o/BnND2OXMaDv4wgv4+sClXVShOiFEBdI1ZIYdqc7ZZWzTL/DMFBg2Bq7v\nU+3p3+Vs4t4dk/iuxxI6B8WZ2xZnKshXf87cXJg+G5qEVThsb7G9oyU5/N/vAxkS/Q9GxNxb8T00\nDUIbQdd2MngsxHmka8hRWdmn9xs2OQS++wpmTlU3RjtCYFteCvfumMTazjPdKwQAAoMgaa4qZPfY\nAEhNrnB48pfrzoYAQEAAqQP6M/nLdRXOa+rbhK+6L2Le/jdZe+Triu9hGHAyDw4ccuIfRIiGRYLA\nHvmFartJs8cFvvw3LF4Azy1WlT6rkVl8lL9vHcWC+FFcE3qJuW2pK7quZkINGgpjhsLvm88cqkmx\nvdb+LVjXZT5Dk+fwa+72igetVrUz3Mk8Z/wJhGhwJAiqY7WqYnJmh8CnH8IrS2D+EmhbTS1+IN9S\nyM1bR/FY1B3c1/wGc9tSH67vA1NnwtNPwvfqU/2ZYnvnqqzY3mk9GiWwPGEyt28bx76i854ArAZs\nT4HSMme0XogGRYKgOkUlplQUrWDdB6pg3PMvQes21Z5uGAYP7ZpOt6B4JrYaYG5b6tNFPWHuQlj8\nPHzwDkl9+hK3ctXZMDg9RpDUp6/NS/SNuJoxLe+n79bRnCrLr3iwzKLGddxwHEyIuiSDxdUpKILf\ndoDFpCeC776CJc/DCy9DtH1bRM7et4r3j37Hjz1ext/Lz5x2uJKsTNVN1Pd20hL/emGxvWp2ZTMM\ng8f2zOJwyXE+6DKnYl0iqVQqxBkya6i2zAyCjRvg2elqTMCO7iCAz49tYODuZ9h48Qpi/Js73gZX\ndfQIjBoMff6uSlvXUIm1lMQ/HuXm8Kt4svWDFQ/qmipB0SjIpMYK4Z5k1lB92/KHqhs0Y67dIZBc\nsI/+u55mbedZDTsEQNVRemEZfP05rHylxt05vroP73Z+lkUH1/LV8V8qHrSe3jTImcUBhXBjEgR1\nYV86TJ0Ak5Kgcze7XpJvKeS2bWNJavMYV4Z0d277XEV4hBo3+eFbWLGsxi+P9mvG251m8K+dU0kv\nzKx4sLRM1YcSQlxAgsDZTp6AJ0bDwCFw6V/sftnQPbO5tFEnHo26w4mNc0Fh4bBgKfzwHbz9eo1f\nfk3oJUxo1Z87t0+g0FJ09oBhqHUgh4+Z2FghGgYJAmcqKYEp4+GqRPj7rXa/bGXWOn49tYPF7Sc4\nr22uLLQJzFuktsP8+P0av3xkzL3EB8QwIuW5igesVlUnqrDYpIYK0TBIEDiLYcD8WdAoRG3jaKft\n+amM2/si73aeRZCXB++81bSZGlR/Yzl880WNXqppGq8kPMV3OZt598g3FQ9arWp9gRtOkhDCWSQI\nnOX9NaqEwlPT1RRGO+RbCrlr+xPMbTvc/cpHOENUDMx9Ua2+/vnHGr20sXcwazo9w9DkOaQVHqx4\nsLBYrTwWQgASBM6x5Q94a6WqH3R+yYQqPJ48h0sbdWJAC9sLqDxOmziYNR9mJ6k9GmqgZ+NOTGzV\nn3t3TKLUes4KY6tV1SLKlRIUQoAEgfmOZcP0p2DiFGgRbffL1h75mg0nt7C43Xjz2+TldfapxNsL\n/H0h0F9t+h7or74uL9us6+ClgysV7uzQGSZMgUlj4eCBGr10ZMy9hPk0Zkr6SxUPWA3YnqpKiwvh\n4WRBWXVqsqCsrAzGDIGLLoUBj9j9FgeKDnPx/x7g064LuKxxZwcae5qXrm50wYEQ1lgtpAoMUDf8\nqkozG4baa6GgCE7lQ04unCpQoWDWympHfPw+vLsaFr2m9nG205GS41y0uR+rOkzlr2G9zh7QNWga\nBh2qL/MhREMgK4trqyZBsGyh2mBm1gK7xwWshpW//fk414RezOTYgbVvp5cOBhARCs3DVU1+O9tQ\nJcNQoXA0B44cV4uynLVPsz2WLYStf6iBZD/790H++vhGHt6dxJ89V9PEp/HZA7oOndpCuP3BIoS7\nkiCoLXuDYPNGVT7i1bdq9Gn1+f2rWXv0G37s8fKFm7LbQ9chwA9aRUJEE3Nu/raUh8LBI5CdA2h1\nHwpWq9rcxmqBKTNrtPnMsOS55JTm8manpIoHvL3gsq7gU4v//0K4ESkx4UwnclQITJxSoxDYmpfC\nM/tW8EbH6TUPAV2HkGDo1h56dlabtjszBEDddBsHQ8e2cHl3aN1C3USd/b7n0nUYPxkOZamppTUw\nu+0wNp3aceGUUosVdqXJlFLhsSQIHGUYMCcJet8APXtVf/5pxdYS7t85mdlthxEXYF8VUkDdCAMD\nVAD06KDCoD54e0OrFioQ2kRXHJB2Nj8/mDFPlfP+6Xu7Xxbo5c/rHZ9mWPJcsoqzzx4wDDhxSnV9\nCeGBJAgc9dG7kJ0NDz1Wo5dNS3+ZuIAYHoy0c6qopqmbbbtW0LNT/QXA+XQdYprD5d2gZfO6C4Pw\nCEiaA/NmqnEZO/Vq3IVBUbczcPcMKnQvWq2QnKEGy4XwMBIEjtifAStehslJ4ONj98s25W5nedY6\nXmo/Ec2ePm5dh/AQ6NUVIiNcc1N2Ly+IjYbLukBYSN0EQofOMGwMPDUGTpyw+2WTWw/kUMkxXs36\nqOIBixV27pUuIuFxJAhqy2KB2dOh/0Bo2drulxVbS3hw13QWxI+iuW949S/w0qFDLHSOd4/BTD9f\n6NpOzcTx9lJTNJ3przfAtb1Vddcy+7al9NG9eaPj0zyZtoSMoqyKB08VQFZ25S8UooGSIKitD94B\n3Qtuv7tGL5uR8RrxAS25t1mfqk/UNbXYq2dnNdfd3YSHqieY8FDnPx08PBj8/eHlRXa/pFNQW0bH\n3Meg3TMv7CJK3Q9FUphOeA4JgtrYn6FmrIyfXKOb3G+ndrEs80OWVtclpOvQtAlc0gn83XhrSm9v\n6BSnFmx5OfFHzctL1XT68Xv48Tu7Xza25QMcLc1h1aFPKx6wWmGHdBEJzyFBUFMWi6p7038gxLS0\n+2Ul1lIe3DWd5+JG0sIvwvaJuqb22E1oU7fTMp2paRP1ZBMU4LyuosYhMG0WzH9WBbUdfHRvlidM\nZvzehRVnEQHkF0LmUSc0VAjX00DuNHXog3fUp9sadgnNzFhBK//m9Gt+o+2TdB06xqmN1l1xQNgR\n/n5wcUfnrnfo0EnN3poyAQoL7XpJj0YJDGpxO0OTZ1/YRbT3gOxdIDyCBEFNlC9iGlezLqEtecks\nyXyPl9o/YbtLyEuHbu1UiYiGStchIRbiWzrvyaDv7dAuQe0FYWfXzuTYh9lVkMF7R7+teMAqs4iE\nZ5AgqImF8+Af99aoS8hiWBi4ewYz2wwh2q9Z5Sd56dA9AUIamdRQF9eiKXRLUH37ZtM0GP0EpOyB\nTz6w6yV+ui+vJUxieMo8skvOm4aaXwiZR8xvpxAuRILAXv/5AfZlwD0P1OhlSw6+R4Dux8MtbGxV\n6aWrFcKNgkxopBsJCVaD4X4+5neD+furvSCWvwR7dtn1kstDunFPs78xsrLtLfcelC4i0aBJENij\noABenAujJ4Kvr90vO1B0mOkZr7Ks/ZOVdwnpp58EggNNbKwbCfCDSzqrabJmdxW1bA3Dx0LSU+rv\nzw4z2gxmQ+4Wvjr+S8UD0kUkGjgJAnuseBl6XAIX9azRy4alzGNo1F10CIq98KB+ekzA054Ezufj\nDRd1gEbB5ofB9X2gS3d4YY5dpwd5BbCk3XiG7JlNoaWo4kGZRSQaMAmC6mzdCl98CoNH1OhlHx79\nnp35aTzResCFB3VNrbz1lDGB6nh5qVBs0tj8GUXDx6ktLr/+3K7Tbwy/kosbdeCZjBUVD5TPIpKF\nZqIBkiCozktL4ZEh0MT+1b25ZXkMT5nHsoQn8NPP60rSdYhrKRuhnE/XVRmNcJPrFAUEwJRnYNF8\nOLDfrpc8Hz+aZVkfsD0/teIB6SISZjpwGDZvr+9WALIxTfVy8+CPXWr3LzsNS55LoaWYVztMqnhA\n19SMmfhW5raxITEM2J2udkQzc9ObD96BL/+ttrm0o0Dg4oNrWXPka37osQxdOyeYyoM8qql5bROe\nxTAgZf/Z2WjX1KzL2RGyMU1teXvX6BPqxtxtvHf0W+bEDat4oHxTlzj7p556JE1Taw0iTK5RdPvd\nENEUXlls1+mPRd1JibWUFYfWVTxQXotIylWL2igpVR8sD51eye4ihSQlCExUai1j0O6ZPBc3kjCf\nkIoHfX1U10dDWzHsDJqm6hOZ2U2kaao21Pdfw6Zfqj3dS/NiWfsneGLvYo6UnLdhjWFIF5Goudw8\n1RV0Kr9+9/2uhASBiZ4/8DaRvuEXVhYtnyHk7YQFVA2VpqktMUMbmTebKCQUJkxRtaJyT1Z7eo9G\nCfyr+U2MTX2h4gHDUOWqDx8zp12iYTMMtc/3n7uhtKxG3cx1RYLAJAeKDjN73yqWtJ9Qcc1AeVmF\nwIB6a5vb0jToHKem2Jr1JNWzF1xzHSx41q5P9NNiB/HDid/4NufXigesVkjZJ11EomrW0/th7z0A\nVhdMgNMkCEwyOvV5hkbfVXH/YV2DZmHql6gdXVcb3QT6mxcGg4ZCWip880W1pwZ7B7K43Xge2/Ms\nRZbzpo5arGpgW7qIRGVKSuG3nZBt8sQHJ5AgMMFXx39h86kdTGzVv+IBX1+1x7BwjJeXWoHta/92\noFXy84enkmDxAjh8qNrTb464iq5Bcczd/8aFB0/mqRlOQpyroEiNBxQUuvSTQDkJAgcVW0sYljyX\nF+PHEuDlf/aArkPX+Iazp0B98/FWNZnMGmdplwB33Quzptn1aW1B/GieP7CGtMKDFQ9YrbAnXX36\nEwIgN189CbjoeEBl5C7loPn73yIhsDU3R1x19pu6Dm1jZFzAbP6+6snArHC9519qn+N3VwOQlpVJ\nv5XLuHb5EvqtXEZaVuaZU1v7t2BMy/sZmTL/wutYDRUGQpw4pQaFLZb6bkmNyIKyKqSlpzN50SIO\nZucQrekk9elLmxZRZ45nFGVxyeYH2HTJKtoERJc3Tg1u9kiQqaLOcvwkbE8x55E78wA8NoADT00n\n8X8/kzqgv1qNXFhI3MpVfH37fWf+zoutJXTddC/z40ZWDH44valQG4ho4nibhHs6cQq2JtdsPMDH\nG67o4bw2nUcWlNVQWno6vadO5a3rr2f9gAd46+476P3h6gqfEkelzGdEzD1nQwDUzb9TWwkBZwoL\nUQvzzHgyiIqBx4ZjfXYa+++/T4UAQEAAqQP6M/nLswvK/HRfFrYby4iU5y4sSmc9PXBcWuZ4m4T7\nyc2reQi4EAkCGyYvXUrqPffYvDF8fmwDW/NTGdfynP0JdF3tvuVnf6lqUUtRzSDSpG0vb+zLvuBG\nTHvnnYrfDwgg01rxEb9P2OVcFJzA7H2vX3gdi1W6iDxRfiFs2eO2IQASBDYdLCo6GwLlTt8YiizF\nDEuex8L4sfh7+Z09HhwIkVVsTC/MFd8KGgU6/vSlaay+4v946LPPuHTnzrPfLywkSr9wcHp+/CgW\nHlxLauGBigcMA47nwrETF7xGNFDFJafHBNw3BECCwKZof/8LN0A/fWOYu/8NugXHc0P4FWePlfcR\nS5dQ3dE06BJvyrTScbf9k5kdOrJq1iz8SkrOjBEk9el7wbmt/CMZ17Ifw5PnccFYVfkCojLpImrw\nLBYVAg3g71qCwIakwYOJW7PmbBicvjEMuvYyXjiwhgXxo8+erOvQugX4+1V+MeE83t7Qrb3DXURt\nWkQxcvgETulerBg9hvvXflBhoPh8o1veT2rhAdYd+/HCgxYrJO9zqD3CxRkGbE+FohK3mSJaFZk1\nVIXyWUOZ2TlEnZ41NCJ7Hpc37soTrR88e6K/H1zWRZ4G6tOxk7DDhJlEOcfhoXthxjzo3LXKU785\nvpFH9sxk+6XvEHjuGhI4vb9CnBrYFg1P2kG1n4Cj4wIuMmtIgqA6BUXw2w6wWFmX/SNjU19gy6Vv\nn91wprygnOw2Vv/SD8J+E/5xrv8GXnsJXn1TrUKuwj+3P0FCYGumt3nswoPe3tCrqxQbbGjM+tAB\nLhME0jVkp0JLESNSnmNRu/FnQ0DTVKlkCQHX0DoKQoIdfzJL/CvEt1NhUI3n4kay5OB7JBdU0hVk\nsUByhmNtEa6lqBh2prpF2YiacLkg0DTtBk3TdmmatkfTtAn13Z5yz+5bRc9Gnegd1uvsNzVNdhtz\nJeVrOMzY7GPEBFWUbtufVZ4W49+cCa36Myx57oUDx4YB2SfUAjjh/qxWtVbAzWcIVcalgkDTNB1Y\nBPQBOgP3aprWoX5bBSkF+1l88F3mx408+83yAWKzCqEJc3h7q2qljq4vCA2FkePh2elQVFTlqSNi\n7mFf8WE+yl5/4cEzs4jcq+SAqMTeA2pwuAFyqSAALgOSDcPIMAyjFFgD3FqfDTIMg+F75jKhVX9i\n/JufPeDtBTHNbb9Q1J/gQLWwz9EwuPo6SOgIry6p8jRf3YdF7cYxMmU++ZbCC08oky4it5eTC1lH\n3XrRWFVq9C9F07S7zvl9M03TEk1uTzSw/5yvD5z+Xr35OGUd6YWZjIi55+w3dV2Vl5bKoq6rRVMI\nD3V8d7PhY9X2llt+r/K065pcyhUh3XgmY/mFB6WLyL2VlsEO88cF/ndqJ8dKXGPxYbWdqZqm9QWu\nBb4A2pV/3zCMI5qmtdM07S7DMN51YhsrNW3atDO/T0xMJDEx0Snvsz17J4sTJuCrn9MFFOivbjLC\ntSW0hs15jj3Oh4TCqIlqe8vXVoO/7VlE8+JG0H3TffSP/DsJgbEVD1qtap/jXl1V95VwD4YBu/aa\nPi5QYCnijm3jWdllOteSaOq1z7V+/XrWr19f7XnVTh/VNC0ceBDVb38l8F/gm9O/NgMPGIZRSeGV\nmtM07S/ANMMwbjj99UTAMAxj9nnn1cv0UUB9wuyeAI2D6+b9hWPyC1VteEcf6WdMgibhMHRUlafN\n3/8WXxz/L192W1hxy1JQg9kRodApzrG2iLpzKFstDjS5S2jS3qWkFh3g7e6z3WP6qGEYxwzDmGcY\nRm9gPjAHaAasAnKBu6p6fQ1tAuI1TWutaZovcA/wiYnXd1xIIwkBdxIUYM54wbCx8N1XsLXqWUTD\nov9JZvFR3jv67YUHDUPNQc92je4AUY3iEqeEQHLBPl7KfJ95cSNMva4javqvY6NhGF8ahjHKMIxO\nQGtMHMymr9O6AAAgAElEQVQ1DMMCPA58BWwH1hiGsbPqV9UhXaaLuqXICAhr7Nj6gpBQGDEO5kyH\nYtuziHx0b5a0n8Do1AXklRVceEL5LCIpV+3aDEN15Znc82AYBiNSnmNCq/5E+zUz9dqOqFEQGIax\n7ryvjxuGYWpcGobxhWEYCYZhtDMM41kzr+2wiCZqfEC4F02DhDaOry+4+jqIT4Dly6o+LfRiEkMv\nISnjtcpPKA8D4boOZcOpAtOD4JNjP5JWdN7kExcg017sYTXU00Cbep3AJBzh7aVq/zg6i2jEOPj6\nc9i+pcrT5rYdzvJDn7Ajf++FBw1D7WZ1+JhjbRHOUVwCKftN7xIqtBQxMmU+i9qNqzj5xAVIENjD\nMKBZuFQXdXeNg6FVC8fGC0KbwPBxaqFZFV1EkX4RTG79MI9XtuIY1E0mOUPddITrMAz1tOaEySjP\n7lvFpY06cX2Ty87uj/3KQvpNmEBaerrp71cTEgT2iq28HLFwM61aQJA/OPJgkHg9tI2Hla9UedqQ\nqH9wvPQk7xz5uvITLFbYYX4/tHDA4WOQm2/630lq4QEWH3yX5+JGkJaVSe8PV/PWXbfzwwP381Zi\nIr2nTq3XMJAgqE6AH/TsLNtPNhSaBp3iHZ9FNHI8fPEp7Nhm8xRv3ZvF7ScwJvV5csvyKj8prwAy\njzjWFmGOklJIMX+WkGEYjEh+jvEt7qXlrzs4MGks/97wH/LuvJO7169X2+Decw+Tly419X1rQoKg\nOpqmpiCKhsPfF9q1diwMmoTBsDGnZxEV2zztypDu/C2sF0+n23h6sFpVDZuCSkpTiLq1J90pVUU/\nObqe7j/tZNxTH8O7b/NbeBh3JiXR4r33eOe669RJAQFkVlPTypkkCIRnah7u+JTSa3tDq1h4/dUq\nT5vddhivH/6MrXkplZ9gNWBbaoOtY+MWjp2AnFOmdwkV7E8hcvxTjP+zMdqkJFj4Cps6dWV7ZCSn\ngoLOnlhYSFQVq9adTYJAeK6EWMc2jdE0GDkBPvsEdm23eVoz3zCejh3E0OQ5lQ8cgxo0TjtY+7aI\n2iuzqAFis4N44wYsQx9kT8+2hCx9G7p0ByCpT1/iVq6quA3umjUkDR5s7vvXgOxQJjxbTi5sS3as\nS+CbL+DNFfDyG+Bb+ViSxbBw2f8GMCrmXvpF3lT5dXQNuraHUNnoqE7tTleDxGbeU955k9K1r3Pb\nnaW8eue7tPCLqHA4LSuTyV+uIxMrUc0iSBo8mDaxsea9vw2yVaUQtiRnwKFjtf9EaBgwaRy0iYOB\ntj/Vbczdxu3bxrHj0rWE+ti42ft4q/2vpTBd3Th5CrbsMW9swDBgxTKM9d9y90ONuDq+D8Ni/mn7\nfNmqUggXEdfSsVXHmgajJ8K/P4I9u2ye1qtxF/4efiVT06tYmVxmgZ3OmccuzlNeEdbMEFi2EDb8\nyPuT7iY1uITBUXeac20nkyAQQtcdX3UcHgGDR8CzT0Npqc3TZrV5nLePfMWfeXsqP6F81fGh7Nq3\nRdgnPdPcmk9vvw4bfyZ37jyGZy9nafuJeOvu8WQnQSAEQKMgiG7u2JTS3jdCZAt4s5LNaU6L8A1l\nRpvHGLJnNlZbZbqsVlXiQKaUOk9+IRw8bN7TwBefwsfvwZwXmXRsNX3Dr6JX4y7mXLsOSBAIUS42\nCvwcqAGjaTD6Cfj4fUjebfO0gS1uo8yw8Pqhf9u+ltUKW1NkSqkzGIa5O45t/Fl1Cc15kd/8j7P2\nyNfMbDvEnGvXEQkCIcrputo0xpEuooim8NjwKruIdE1ncbvxTNy7iJzSXNvXKimBPbLXsekOHDZv\nE/p96TBrGkyfg7VVawbveZZZbYcS7uNeOxhKEAhxruBAxwvT9fk7NG0Gq1faPKVn407c2fQ6Ju5d\nZPs6VgOO5sCR47Vvi6iosFiNDZjxpJWfp2aLDRwMXbvzatZH+Gje9I+82fFr1zEJAiHO16qFKkNR\nW5oGY56ED9+FFBuDwsDMtkP597EN/OfEH7avZbWqee6F9Vd+oMEo32zGjBCwWmHmNOh+Mdx8O0dK\njjMp7SWWtJ+ArrnfbdX9WiyEs2ma411ETZvBo4/D7OlQVvnMlBDvYF5oN4ZBe2ZSbK2iq8Jqha3J\npm+g7nEyj6pBYjOsXgUnc2D4WABGpsznwci+dAtuZ87165gEgRCVCQpwvIvohr5q/4LVq2yeckfE\ntcQHxDBn3+tVX6u4RBVFE7VTVKyK+5nxNLDtT3h/DUyZCT4+fH5sAxtztzE19hHHr11PJAiEsMWM\nLqJxT6mbxt7KC85pmsbiduN54cAadhek276W1VCb3mcdrX17PNWZ/YdNCIFTuZA0GcY+Cc2ak28p\nZEjyHJa2n0igl/tuYytBIIQtZnQRNYuER4aqHc1sdBG19I9kcuzDPLp7lu2idHB6fcE+OJVf+/Z4\noqyjkFcIjs4WNQyY+wxceTVceQ0AU9OW8X8h3flb2F8cb2c9kiAQoipBAdAy0rEuor/fCo0bwztv\n2Dzl8ei7ybcWsuLQuqqvZTXUeEGJ7dXL4hyFxZBqUpfQug/h4H54dBgAv53axZuHP2d+3CjHr13P\nJAiEqE6rFo4vNBs3CdauhrTUSk/x0rx4uf1TPLF3MUdKqpkuWlYG22SxWbXOLBwz4f9TWiq8thSm\nzgQ/P8qsZQzcPYM5ccNp6tvE8evXMwkCIapjxkKz5pHw8OAqZxFd1CiBf0XexKiU+VVfywDyCyB5\nX+3b4wn2H4ICE6bdlpXBzKnwyBC1ERHwwsE1hHk35oHmNkqKuxkJAiHsERx4uhaRA2HQ93YIDIS1\nb9k8ZVrsIH7O3cqXx/9b9bWshlpolimDx5XKK4AMkxaOvbEcwsLh77cBkFZ4kFkZK3mp/RNojuxw\n50IkCISwV2wU+DjaRTQZ3nkTMtIqPSXIK4Bl7Z9g0O6Ztje8L2e1Quo+OFFFmQpPZLGe7jozoZbQ\n7p2qmNy4SaBq+fPonlmMbdmP+MCWjl/fRUgQCGEvXYdObR17KmgRBQ89qmYRWSyVnvK3sL/wt7Be\njE9dWP31rIa66ZnRBdJQpOyrshS43YqLYdZUeHy0qiEFvJr1ETlluYxt2c/x67sQCQIhaqJxMERG\nONhFdAf4+8O7q22eMi9uJJ8d38A3xzdWfz2LFf7cbc7Nz91ln67NZMbTwIpl0KoNXN8HgH1Fh3gy\nbQkrEqa4zT4D9pIgEKKm2saAlwOb3uu6Wmi2ehVkpFd6Soh3MC+3f5JH9szkVJkd6wZKS+HPPZ5d\nhqKoWO3uZsa4wNY/4avP1M5zp7uEBu6ewaiY++gSHO/49V2MBIEQNeXlBR3bOra2ICoGBgyCOba7\niG4Iv4LrQnsyfq8dXUQGqjDd9hTP3OayvB6TGSFQWKhKS4+aqEqEoLqEjpfmMr7lA45f3wVJEAhR\nG00aQ0SoY11Et/1DbVL/wTs2T3kubiSfHvuJb3N+rf56VgNO5qlqpZ4WBnsy1BOBGZYthC7d4KpE\n4GyX0MoODa9LqJwEgRC11a6VY08Fug7jp6jpiQcqXxMQ6tOIZe2fZODuZ+zrIrJa1R4Gew/Uvl3u\nJuuo+jObMS6weSNs+AGGjQFo8F1C5SQIhKgtb29IaONYGETHQP+BVc4iuin8ShJDL2Zs6gv2XdNq\nVesLMjJr3y53kZunZgmZ0SWUlwdzZ6ipoo0aAw2/S6icBIEQjogIhSaN1BqB2rr9bvDS1foCG56P\nH8NXORtZl/2jfde0WmHfIbW6tqEqLjk9LmBSN9iS5+HSv8BllwOQXLCPJ9OWsKrj1AbbJVROgkAI\nR7WPdWysQNfhyelqxfGuHZWeEuIdzOsdpjFoz0wOlxyz77pWq9qWsSGGQZkF/tit/muGXzbAb5tg\nyEgASq1l9Ns5hSmtB9I5KM6c93BhEgRCOMrXB+IdHC9oHgkjxsGMyVBQUOkpV4VexEORt/DQrqSq\ny1Wf60wYZNW+ba6mfIZQsWMb0KdlZdJv5TL6Ll3A8elPkvXwYAgMAmB6xiuE+TTm8ei7zWixeuLz\n9oKwEPXL2ws0oLTyulN1TYJACDM0D4dGgeofd21d21vNVllku+jctNhBHCk9zkuZ79t/XasV0rMg\n7aD7zyYyDLVWIK/AoT9LWlYmvT9czVt338Fdx4+y+rrruCp5G2lZmfznxB+8mvUxKxKmOF5LSNfU\nB4UObeGKHtC1nfp1eXeIiYQAP8eubxLN7k8WLkTTNMMd2y0auKJi2LTdsYHLgnx45AEY9Dhcc12l\np+wuSOf/fn+En3q8QoegWPuvresqsNq1cmxMo74YhpomeuS4w4PD/VYu46277+CW//2P55Yupfur\nr1IA3PXuGjYlfMOL8WPoG3G1Y+3VdQgPgYRYxxYgmkhTi+Mu+MuXJwIhzOLvB22jHesiCgyCSUnw\n/Gw4UnnffkJgLDPaDOb+nZOr3vT+fFYrHD6mFp25214GhgEp+00JAYCDVgthJSUsff55HpwwgYKA\nAAgI4MewDfRp8hdzQiCmuVp46CIhUBUJAiHMFNUMAh3cu7ZjZ7jzHpg5zeaU0kEtbqe1fyTjUl+s\n2bWtVsjJhd92us8uZ4ahpogeyjYtwKJ1LxYtWMCaa6/lP926qW/uW0ehXzbPxY907OK6Dq1bQJto\nt3nykiAQwkyadrr8hIM3gHv/pW6Aa1638TYayxOm8Omxn3jvyLc1u7bVgIJC2Lxd9bW7MsOAXWlw\n6JipTzHPhYbR67ffeOr++9U3ju1ET17MO+2TCPIKqP2FdR2im6pd7dyIjBEI4QwZmWoevyM3ryOH\n4NH+8PRs6Naj0lM25+7gpq0j+e/Fy4kLiKn5e+i6GjOIjKh9O53FYlEltnPzzFsrAJBzHB66l8xR\nExm/dw/7jCL+jF7HpFb9GdfBgYVjugbhoeqDgIs+CdgaI5AgEMIZDAM2bVObpzti4waYNxNeefNM\nAbTzLTqwlhWH1rHholfx96rFLBQXHNSkqBi27IGiEnNnOhkGTJ0A0S3h0WEYhsEDO6fgq/uwvMMU\nx64dFAAXd3RsjMjJZLBYiLqkaY7vcwzQ60r4203wzBSb4wVDo++ibUA0o1MX1O49rFbIPgG/boOT\npxxorEmOn4TNO1SImv2B7/N1cGC/qvyKKiHxZ34yi9qNd+y63l7Qrb1Lh0BV3LPVQriD4ECIbuZ4\nGDz4KJQUw1srKj2saRqvJkzi25xNrMj6pHbvYRhq8HjLHlW9tKweFjpZrJCcoWY12Qg9hxzYryqL\nTkoCPz825W7nybQlvNvpWQK9HBjg1zXoEq/WC7gpCQIhnCk22rF9jkEVt5v8DHz8viqDUIkQ72A+\n6jKPCXsX8cvJrbV/L6uhppj+slUVrqurLticXNi09fSgsBPes6xMPVU98BC0jSerOJs7to/nlfZP\n1Wwtxvl0HVpHQUgj05paHyQIhHAmM/Y5BrVn7hPT1M0s+2ilp3QMasNrCZP4x/aJZBZXfo5dDEN9\nIk/dDxu3qrn7zgqEgiL1FLItBYpLnbe+4Y3XIDgY7vgnRZZi7tg+nkda3MZtTRNrf00NtZq8ZaRZ\nraw3MlgsRF1Izjg9D97Bn9s3lsN//wPPvwS+vpWe8kzGctZl/8T6Hi/VbvD4fLquauW0jITIcMef\ncAwDTuVDRhacyHXOE8C5tv0Jk8fDq29hhIXz0O7pnCorYG3nWeiaA5+Fvb3g0i5u1SUks4aEqE8W\nixqMdXQRl2HA1Inq0+24SZVOUzQMg3/ueJJA3Y8VHaY6Xi+nnK6pLTEbB6lSFU1CwL/yMKq03fmF\nanP5Q8dUsbW6WN2cnwcD74cho+CqRF448DbLs9QMq2DvwNpfV9ehc5wqIOdGJAiEqG8nTsHWPY5/\nAi4ogKEPwS13qL0MKpFvKeSa3x/lloirmBL7iGPvVxldBwz136AANTDu73e6qqambvylZWrmT16B\nCgFQ36+rf7uGAdOfhOBGMOZJPs7+gcf2zOLni16jTUB07a+rayoI28ea1tS6YisIGvZuC0K4ktBG\n0Cwcjjg4IBoYCDPmweMPQ5s46HHJBacEeQXwadf5XPH7w7T0a86DLW6xebm0rEwmf7mOg1YL0boX\nSX360qZFVNVtKP80b7WofZJP5qkA0DTVd24AhlX9t758/B7s3weLl/PLya0M3D2Dz7o+71gIgOoa\ni2tlThtdhDwRCFGXLBY1AGtGHfrNG2HmVFiyAiIrL2mwuyCda35/lFUdp9En7PILjpeXY04d0B8C\nAqCwkLiVq/j69vuqDwNXtnsnjB8Oi5ezJ8zK1b8P4rUOk/l7+P85dl1dgx4doFGQOe2sY7KgTAhX\n4OVlTi0igJ694J4H4Kkxqi+8EgmBsbzfZQ4P7JzKptztFxyf/OW6syEAEBBA6oD+TP5ynePtqy+n\nTsHTT8CoCRxqGsCNW0bwTNshJoSArvYQcNMQqIoEgRB1rUlj1UVkRhjcdR907gbTnrC5COzKkO68\nljCJm7eOZktecoVjB62WsyFQLiCATKsTFnTVBcOAOdOh15VkX9GTv/45lAcj+/Jwi1sdv7a/L8S6\n8VNSFSQIhKgP8S3NqeujaTB8LHh5w/xZNgdi+0ZczcJ2Y7lhy3B25aef+X607gWFhRVPLiwkSneR\nmkM1tXolHD1CzsAH6b1lKLdGXMNTrR9y/Lq6pmYJuWgxOUdJEAhRH7y8zKlFBGrl8ZRnIGWPWmdg\nw93NejOr7VB6bxlKauEBAJL69CVu5aqzYXB6jCCpT1/H21XXNvwAH77HqWnTuGHXWK4N7cmMNoNN\n2G5SV08CgQ6Up3ZxMlgsRH3ak65KOpixqOpYNgx5CB4cBDfcbPO0lzM/ICnjNb7stpBOQW3PzBrK\ntFqIsnfWkKtJS4WRj3EyKYk+pcu4pFEHFrUbb84aiuBAVVW0ATwNyDoCIVyRWQvNymWkwajBMGoC\nXHWtzdPePPQZY1Nf4NOuC+jZuJM5711fTp6AwQPIue8ermr2ITeGX8GctsPNCQFdh56dIMDBXedc\nhMwaEsIVeXmpvmczuogAWreBWQvguVmw6Rebp/WLvIllCU9y09aR/HDif+a8d30oLoLJ4zlxeU96\nRqyhX/MbzQ2BttENJgSqIkEgRH1rHKz2Ojarln1CR5g+B2ZMhi1/2Dzt1ohreLvTDP6xfSKrD39h\nznvXJYsFZkzmWIgPXS7awPhW/2Ji6wHmhICGWjEd1czxa7kBCQIhXEGbaHOLl3XrAZNnwJTxsPVP\nm6dd3+Qyvuu+lKfSljIl7SWsRh3U/zGDYcCLczmUs4/uf93Dkg4TeTTqDvOur+lqML8BjAvYQ4JA\nCFeg62pzE7O6iEAtOHvyaZg8Dn633f3TNTiejRev4JucX7lr+0ROllW+OM2VWF5/lczfvqP3Pwr4\n/JKl3BJxjXkX13U1vdfegnoNgASBEC4gLT2dftOnce2br9Jv1TLSsjLNufBll8PUWWql7a//tXla\nM98wvuu+lBa+EVy8uR+bc3eY8/5OcOKNJWSuW8noR2L59opVdA2ON+/iGqq6amSEedd0AzJrSIh6\nlpaeTu+pU0m95x7n1fvZ9idMGgcjx0PiX6s89d0j3zAkeTbjW/6L0S3vw0tzncVlW1+bRtBnn/PZ\npH4M6THUsf0EKuOlw2Vd3WqPgZqQ6aNCuKh+EybwVmJixVIPhYXcv/YD3hzwqHlvlLwbnhwNd94D\n/+xXZf/33sIDDNz9DHmWApYnTKaLmZ+6ayGrOJtvljzOVf/JIHvOTHrG2Z4aW2u6Dh3bQEQT86/t\nImT6qBAu6mBRUeX1fsweuG2XAItfg68+g+dnV7lBfduAGL7tvoRHWtzGtX8OZlzqC+SU5prbHjuU\nWEt5cd/bfDz9Nv62+RiRi991UghoEB7SoEOgKi4TBJqm/UPTtG2aplk0Tbu4vtsjRF2J9vevvN6P\nv7/5s1aaRcLCV+DgAXhilFqMZYOmaTwSdTtbe75Nblk+Cb/+gwX7V1NgKTK3TZWwGBbWHP6K7j//\ng24vvM4DObE0f+l9/CNjnPOGXl5uudGMWVwmCICtwO3AD/XdECHqUtLgwcStWVOx3s+aNSSNGw0+\nTtg7KigYnn0e4trBoH/BrgvLU58r0i+CZQlP8n2Ppfx48ndif7mFqWnLOFJy3PSm5ZUVsOjAWjr8\nehev71jFf98MIDGwC0ELXoVGjU1/P0A9DXSKU7ureSiXGyPQNO17YIxhGL9VcY6MEYgGJS09nclL\nl5JZVESUvz9JgwfTJjZWbfL+xy7nbfD+0/dqFfL9A9TYgR2L2nYXpLNg/9u8feRLrg69iPua3UDf\n8KtqvQdwqbWMb3I28s6Rr/nk2E9cG3oJU05dRrfnVqD1vQ0eeNi8xXbn03WIagpxLZ1zfRfjNoPF\nEgRCnGf/IUjPdN5m7wcPwKxpasbMhCkQZV/3y6myfD7O/oG3jnzBTyf+4KJGCVwTcjHdg9vRJSiO\naL+mNPIKqrDSt9BSRGZJNjvz09ian8JPJ/9gw8k/6RzUln82681dYdcS9cHn8P4amDgFel3pnD9z\nuUB/uKST84LGxbhEEGia9jXQ/NxvoXY1fcowjHWnz5EgEOJchgFb9qh9gZ31c2+xwLur4e3X4d5/\nqacDH/unUOZbCvn55BZ+Ovk7W/NT2Z6/l8zio1iwEqD7oaFRbC2h1CijhW8EHQJb0zmoLfHWVnz7\ny36OlflyaW4uSVu34BfaBMZPhuaRzvmzltN1FQKBDb+WUDmXCAJ72BsEU6dOPfN1YmIiiYmJddA6\nIepJaRls2mbOXsdV2Z8Bi5+HA/tgyEi4/P8cGrDOtxRSZC3GMMBP9yHYK/DME0L5fslH7r6LJ99/\nn4H//jcvxLfjobGTaRPl4Abz1dF1aNeqwS8cW79+PevXrz/z9dNPP+1WQTDWMAyba+LliUB4pNw8\n+HO388YLzrVxAyx5AQID4f4H4YqrTO8+efiVRTQO9GX8++/zxWWX8dTDD5MVFGT++onzaaeninau\n37UR9cHlnwg0TbsNWAhEACeAPwzDuNHGuRIEwjM5e7zgXFYr/LQe3lwOxcVw0y3Q+0YId/BT9P4M\n+OJTct5fwze9ejHz/vv5o127M4evXfE63z00xLH3qIqvD1zaxSNnCbl8ENSEBIHwWIYB21IgJ9d5\n4wWVveeW3+HLf8OP30P7DqqG0UU9oW189WMJxUWQkgybf4FfNsChLOh9I+NKi5n3yEPOX1F9Ll2D\nizqqXcc8kASBEA1FmQU2b4fikrp/78JC+O1XtenNH/+DzEyIaQnNmkNYOPj6qfOKCiHnOBzOUue0\nbAWXXKYCpMcl4O19ZowgdUB/59VYOpeuQ1yMx+wxUBkJAiEakoIi+N+OuukiqkpxEWSkQ/YROH7s\n7Jabvr6qCymiqdo1zbfyks51tl+yrkF4KHRs6zF7DFRGgkCIhib7BOzcW/9h4A78/dTew16eNy5w\nLik6J0RDExEKMc09ZjFUrXnp0K29x4dAVeQnSAh3FhsFoY3M3dmsISmvIxTgV98tcWkSBEK4M02D\nTm1V14cH931XStehdRSEhdR3S1yeBIEQ7s7LS3V9eOC8eJvK9xdo6eQyFQ2EBIEQDYGfrwoDGS9Q\nT0aBAdChjTwl2Ul+aoRoKIIDoUu8jBf4eEso1pD8nxKiIWnSWO205alh4OUFPTo4Z0OfBkyCQIiG\npnk4tInxvE/Eug7d28sMoVrwsJ8UITxETHNo6UFrDHRNdYs1CqrvlrglD/kpEcIDtY5S2zA29DDQ\nNegYp7rFRK008J8QITyYpkHbGGjRgMNA1yAhVq2yFrXWQH86hBCACoO4GIhugGGga2pgvFl4fbfE\n7TWwnwwhxAU0Ddq2hFaRDWc2ka6pdQLNJQTMINVHhfAkmUchdV/dbHfpLLoOneOkdEQtSBlqIYTi\nzuWry8tpNJbZQbUhQSCEOOtUPmzZo3Y7cweapvYa7p4g6wQcIEEghKiopBS27oGCYtd+OtB19QTQ\nOQ68ZcWwIyQIhBAXsloh9QAcOuqa4wa6BjGRat8FKSDnMAkCIYRtx0+qcQOLFVzh35auqfGATnFq\n4x1hCgkCIUTVSstgTzocz63friJdg/Am0L6VdAWZTIJACGGf4ydVIJRa6jYQdF1VDe0QC6FSLsIZ\nJAiEEPazWiHrKKRlqq4iZwaCrqv+/zbR0CKi4a2AdiESBEKImrOcDoR9WSoMLCYGgpeubvqtWqh6\nSF4SAM4mQSCEqD3DgJxcyDyixhB0rXah4KWr2UlNGkNUMwhrLLOB6pAEgRDCHFYrnDilxhJOnILC\nYhUUugace48x1E1fAwL81eyfJo1V/798+q8XEgRCCOcwDLVCuaQUysrU15qmZvz4+oC3l3zqdxG2\ngkBi2Q7r16+v7yYI4bo0Tc32CQqAkEbqE39II/W1j7eEQBVc5d4iQWAHV/nLEkI0LK5yb5EgEEII\nDydBIIQQHs5tB4vruw1CCOGOGsysISGEEOaRriEhhPBwEgRCCOHhJAiqoGnaDZqm7dI0bY+maRPq\nuz1CiIZB07TXNE07rGnalvpuC0gQ2KRpmg4sAvoAnYF7NU3rUL+tEkI0ECtQ9xaXIEFg22VAsmEY\nGYZhlAJrgFvruU1CiAbAMIz/ADn13Y5yEgS2RQP7z/n6wOnvCSFEgyJBIIQQHk6CwLaDQKtzvo45\n/T0hhGhQJAhs2wTEa5rWWtM0X+Ae4JN6bpMQouE4fwOHeiNBYINhGBbgceArYDuwxjCMnfXbKiFE\nQ6Bp2mrgZ6C9pmn7NE17sF7bIyUmhBDCs8kTgRBCeDgJAiGE8HASBEII4eEkCIQQwsNJEAghhIeT\nIBBCCA8nQSCEEB5OgkAIITycBIEQQng4CQIhhPBwEgRCCOHhvOu7AUI0BJqmtQOuBsKAMqAI6Ahs\nNAzjrfpsmxDVkSAQwkGapmnADYZhLNQ0rSmQBrQEBgHH67VxQthBqo8K4SBN03QgwDCMfE3TbgUe\nMW/5Qy0AAACiSURBVAzj5vpulxD2kicCIRxkGIYVyD/95Y3AlwCapvkAumEYxfXVNiHsIYPFQjhI\n07SemqbNP72TXV9g4+lD/ZEPW8INSNeQEA7SNO164CHgF+AU0BnYAewwDGNjVa8VwhVIEAghhIeT\nriEhhPBwEgRCCOHhJAiEEMLDSRAIIYSHkyAQQggPJ0EghBAeToJACCE8nASBEEJ4OAkCIYTwcP8P\nol2Ma66SdT8AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def phi(x, M):\n", " return x[:,None] ** np.arange(M + 1)\n", "\n", "N = 10\n", "\n", "# 生成 0，1 之间等距的 N 个 数\n", "x_tr = np.linspace(0, 1, N)\n", "\n", "# 计算 t\n", "t_tr = np.sin(2 * np.pi * x_tr) + 0.25 * np.random.randn(N)\n", "\n", "# 加正则项的解\n", "M = 9\n", "alpha = 5e-3\n", "beta = 11.1\n", "lam = alpha / beta\n", "\n", "phi_x_tr = phi(x_tr, M)\n", "A_0 = phi_x_tr.T.dot(phi_x_tr) + lam * np.eye(M+1)\n", "y_0 = t_tr.dot(phi_x_tr)\n", "\n", "# 求解 Aw=y\n", "coeff = np.linalg.solve(A_0, y_0)[::-1]\n", "\n", "f = np.poly1d(coeff)\n", "\n", "# 绘图\n", "\n", "xx = np.linspace(0, 1, 500)\n", "\n", "# Bayes估计的均值和标准差\n", "S = np.linalg.inv(A_0 * beta)\n", "\n", "m_xx = beta * phi(xx, M).dot(S).dot(y_0)\n", "s_xx = np.sqrt(1 / beta + phi(xx, M).dot(S).dot(phi(xx, M).T).diagonal())\n", "\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(x_tr, t_tr, 'co')\n", "ax.plot(xx, np.sin(2 * np.pi * xx), 'g')\n", "ax.plot(xx, f(xx), 'r')\n", "ax.fill_between(xx, m_xx-s_xx, m_xx+s_xx, color=\"pink\")\n", "ax.set_xlim(-0.1, 1.1)\n", "ax.set_ylim(-1.5, 1.5)\n", "ax.set_xticks([0, 1])\n", "ax.set_yticks([-1, 0, 1])\n", "ax.set_xlabel(\"$x$\", fontsize=\"x-large\")\n", "ax.set_ylabel(\"$t$\", fontsize=\"x-large\")\n", "\n", "plt.show()" ] } ], "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 }