{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Metropolis-Hastings算法（后面简称MH算法）是MCMC的代表算法，下面先对MH做一个简单的推导，然后写出MH算法流程，最后再对细节做讨论分析，本部分内容主要参考了[刘建平：《MCMC(三)MCMC采样和M-H采样》](https://www.cnblogs.com/pinard/p/6638955.html)\n",
    "\n",
    "### 一.算法推导\n",
    "\n",
    "根据上一节的讨论，我们知道只要找到一个概率转移函数$p(x\\rightarrow x')$，使得它与我们的目标分布$p(x)$满足如下的一个细致平衡方程：   \n",
    "\n",
    "$$\n",
    "p(x)p(x\\rightarrow x')=p(x')p(x'\\rightarrow x)\n",
    "$$  \n",
    "\n",
    "那么，就能保证对马尔科夫链$p(x\\rightarrow x')$随机游走采样的经验分布近似于我们的目标分布$p(x)$，但是这样的$p(x\\rightarrow x')$没有那么容易构造，这里借鉴接收-拒绝采样的思路，首先我们造一个容易采样的分布$q(x\\rightarrow x')$，称为建议分布，利用它来采样，显然这时它们很难满足细致平衡方程： \n",
    "\n",
    "$$\n",
    "p(x)q(x\\rightarrow x')\\neq p(x')q(x'\\rightarrow x)\n",
    "$$\n",
    "\n",
    "然后，我们再造一个接收概率$\\alpha(x\\rightarrow x')$来调节：   \n",
    "\n",
    "$$\n",
    "p(x)q(x\\rightarrow x')\\alpha(x\\rightarrow x')= p(x')q(x'\\rightarrow x)\\alpha(x'\\rightarrow x)\n",
    "$$  \n",
    "\n",
    "很Nice的思路！将一个难的问题拆解为两个相对容易解决的问题，那么$\\alpha(x\\rightarrow x')$如果构造勒，这就很简单了，可以令：   \n",
    "\n",
    "$$\n",
    "\\alpha(x\\rightarrow x')=p(x')q(x'\\rightarrow x)\\\\\n",
    "\\alpha(x'\\rightarrow x)=p(x)q(x\\rightarrow x')\n",
    "$$  \n",
    "\n",
    "简直是完美，$\\alpha(x\\rightarrow x')$也不用费解心机去造了，我们只要找一个好抽样的$q(x\\rightarrow x')$就可以啦，所以捋一下现在的采样流程：\n",
    "\n",
    "####  采样流程\n",
    "\n",
    "输入：抽样的目标分布密度函数$p(x)$，正整数$m,n,m<n$  \n",
    "输出：$p(x)$的随机样本$x_{m+1},x_{m+2},...,x_n$\n",
    ">（1）任意选择一个初始值$x_0$；    \n",
    "\n",
    ">（2）对$i=1,2,...,n$循环执行：   \n",
    "\n",
    ">>（a）设状态$x_{i-1}=x$，按照建议分布$q(x\\rightarrow x')$随机抽取一个候选状态$x'$；    \n",
    "\n",
    ">>（b）计算接收概率：   \n",
    "\n",
    "$$\n",
    "\\alpha(x\\rightarrow x')=p(x')q(x'\\rightarrow x)\n",
    "$$    \n",
    "\n",
    ">>（c）从区间$(0,1)$中按均匀分布随机抽取一个数$u$，如果$u\\leq \\alpha(x\\rightarrow x')$，则取状态$x_i=x'$，否则取$x_i=x$   \n",
    "\n",
    ">（3）得到样本集合$\\{x_{m+1},x_{m+2},...,x_n\\}$\n",
    "\n",
    "但是呢，这里有个问题，那就是$\\alpha(x\\rightarrow x')=p(x')q(x'\\rightarrow x)$可能会低哟，这样会导致很多采样样本被拒绝，效率低下，那么有没有可能找到一个方法提高$\\alpha(x\\rightarrow x')$呢？有的，那就是MH算法...."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 二.MH算法\n",
    "MH算法的思想很朴实，那就是等比例的提升$\\alpha(x\\rightarrow x')$与$\\alpha(x'\\rightarrow x)$，但又不能超过1，这样既能保证细致平衡方程成立，又能提高接受率，举一例子来直观理解一下，比如我们在$(0,1)$上做均匀采样，接受率为0.1和接受率为1其实效果是一样的，采样结果最终都符合$(0,1)$上的均匀分布，但前者会拒绝掉90%的样本，那MH是怎么调整的勒，很简单，它分两种情况：   \n",
    "\n",
    ">（1）若$p(x)q(x\\rightarrow x')>p(x')q(x'\\rightarrow x)$，那么，对等式两边同除以$p(x)q(x\\rightarrow x')$：\n",
    "\n",
    "$$\n",
    "1\\cdot\\alpha(x\\rightarrow x')=\\frac{p(x')q(x'\\rightarrow x)}{p(x)q(x\\rightarrow x')}\\alpha(x'\\rightarrow x)\n",
    "$$   \n",
    "\n",
    ">所以，我们可以令：   \n",
    "\n",
    "$$\n",
    "\\alpha(x\\rightarrow x')=\\frac{p(x')q(x'\\rightarrow x)}{p(x)q(x\\rightarrow x')}\\\\\n",
    "\\alpha(x'\\rightarrow x)=1\n",
    "$$\n",
    "\n",
    ">（2）反过来，若$p(x')q(x'\\rightarrow x)>p(x)q(x\\rightarrow x')$，那么，对等式两边同除以$p(x')q(x'\\rightarrow x)$：\n",
    "\n",
    "$$\n",
    "\\frac{p(x)q(x\\rightarrow x')}{p(x')q(x'\\rightarrow x)}\\cdot\\alpha(x\\rightarrow x')=1\\cdot\\alpha(x'\\rightarrow x)\n",
    "$$   \n",
    "\n",
    ">所以，我们可以令：   \n",
    "\n",
    "$$\n",
    "\\alpha(x\\rightarrow x')=1\\\\\n",
    "\\alpha(x'\\rightarrow x)=\\frac{p(x)q(x\\rightarrow x')}{p(x')q(x'\\rightarrow x)}\n",
    "$$\n",
    "\n",
    "而对于这两种情况，我们其实是可以综合一下，那就是：    \n",
    "\n",
    "$$\n",
    "\\alpha(x\\rightarrow x')=min\\{1,\\frac{p(x')q(x'\\rightarrow x)}{p(x)q(x\\rightarrow x')}\\}\\\\\n",
    "\\alpha(x'\\rightarrow x)=min\\{1,\\frac{p(x)q(x\\rightarrow x')}{p(x')q(x'\\rightarrow x)}\\}\n",
    "$$   \n",
    "\n",
    "到这里，$\\alpha(x\\rightarrow x')$的接收率一定比之前高（$p(x)q(x\\rightarrow x')$大部分情况是远小于1的正数），即：    \n",
    "\n",
    "$$\n",
    "min\\{1,\\frac{p(x')q(x'\\rightarrow x)}{p(x)q(x\\rightarrow x')}\\}>p(x')q(x'\\rightarrow x)\n",
    "$$\n",
    "\n",
    "所以,MH的算法流程就可以梳理出来了\n",
    "\n",
    "#### MH采样流程\n",
    "\n",
    "输入：抽样的目标分布密度函数$p(x)$  \n",
    "输出：$p(x)$的随机样本$x_{m+1},x_{m+2},...,x_n$\n",
    ">（1）任意选择一个初始值$x_0$；    \n",
    "\n",
    ">（2）对$i=1,2,...,n$循环执行：   \n",
    "\n",
    ">>（a）设状态$x_{i-1}=x$，按照建议分布$q(x\\rightarrow x')$随机抽取一个候选状态$x'$；    \n",
    "\n",
    ">>（b）计算接收概率：   \n",
    "\n",
    "$$\n",
    "\\alpha(x\\rightarrow x')=min\\{1,\\frac{p(x')q(x'\\rightarrow x)}{p(x)q(x\\rightarrow x')}\\}\n",
    "$$    \n",
    "\n",
    ">>（c）从区间$(0,1)$中按均匀分布随机抽取一个数$u$，如果$u\\leq \\alpha(x\\rightarrow x')$，则取状态$x_i=x'$，否则取$x_i=x$   \n",
    "\n",
    ">（3）得到样本集合$\\{x_{m+1},x_{m+2},...,x_n\\}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 三.细致平衡方程校验\n",
    "\n",
    "下面我们还是对MH算法做个检验，看它是否满足细致平衡方程（虽然它就是从该方程推导来的...），此时的转移概率为建议分布和接收概率之积：    \n",
    "\n",
    "$$\n",
    "p(x\\rightarrow x')=q(x\\rightarrow x')\\alpha(x \\rightarrow x')\n",
    "$$   \n",
    "\n",
    "那么，我们需要检验的就是下面的等式是否成立：       \n",
    "\n",
    "$$\n",
    "p(x)p(x\\rightarrow x')=p(x')p(x'\\rightarrow x)\n",
    "$$   \n",
    "\n",
    "证明一下：    \n",
    "\n",
    "$$\n",
    "p(x)p(x\\rightarrow x')=p(x)q(x\\rightarrow x')min\\{1,\\frac{p(x')q(x'\\rightarrow x)}{p(x)q(x\\rightarrow x')}\\}\\\\\n",
    "=min\\{p(x)q(x\\rightarrow x'),p(x')p(x'\\rightarrow x)\\}\\\\\n",
    "=p(x')q(x'\\rightarrow x)min\\{1,\\frac{p(x)q(x\\rightarrow x')}{p(x')q(x'\\rightarrow x)}\\}\\\\\n",
    "=p(x')q(x'\\rightarrow x)\n",
    "$$  \n",
    "\n",
    "好勒，最关键的一步得到了证明，接下来就到了放心的造$q(x\\rightarrow x')$的时候了"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 四.$q(x\\rightarrow x')如何选择$   \n",
    "\n",
    "$q(x\\rightarrow x')$比较省事儿的一种方法就是选择对称的建议分布，即：   \n",
    "\n",
    "$$\n",
    "q(x\\rightarrow x')=q(x'\\rightarrow x)\n",
    "$$   \n",
    "\n",
    "这时接收概率的计算也变得简单了：    \n",
    "\n",
    "$$\n",
    "\\alpha(x\\rightarrow x')=min\\{1,\\frac{p(x')}{p(x)}\\}\n",
    "$$    \n",
    "\n",
    "\n",
    "$q(x\\rightarrow x')$通常可以选择：   \n",
    "\n",
    "（1）以$x$为均值，其协方差为常数矩阵的高斯分布；    \n",
    "\n",
    "（2）或者选择：   \n",
    "\n",
    "$$\n",
    "q(x\\rightarrow x')\\propto exp(-\\frac{(x-x')^2}{2})\n",
    "$$  \n",
    "\n",
    "接下来，就用一个案例来实操一下"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 五.案例\n",
    "\n",
    "我还是选择前面的那个分布，哈哈哈哈哈~~~    \n",
    "\n",
    "$$\n",
    "p(x)=\\frac{1}{1179}[(x-2)^2+(x-5)^3+100cos(x)+106],0<x<10\n",
    "$$  \n",
    "\n",
    "概率密度函数如下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1be1ef9c320>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8FdXZwPHfkx1CCGSBhIQQlrCELUCIuKBUAUFkad3Ava+vVltarXbRLtra2tb2rVqttdq6tipSFUVAcUHEDSRACIQ1hCUbJCwJkZD9ef/IpU1jIDeQ3Ln35vl+PveTuTNnZp6B3PtkzjlzjqgqxhhjTIDTARhjjPEOlhCMMcYAlhCMMca4WEIwxhgDWEIwxhjjYgnBGGMMYAnBGGOMiyUEY4wxgJsJQUSmich2EckVkbtb2H6niGwRkWwR+UBE+jXZdoOI7HS9bmiyfpyIbHId81ERkfa5JGOMMadDWntSWUQCgR3AFKAAWAvMU9UtTcp8DVijqpUichswSVWvEpEoIBNIBxRYB4xT1SMi8gVwO7AaWAY8qqpvnyqWmJgYTU5OPr0rNcaYTmrdunUHVTW2tXJBbhwrA8hV1TwAEVkAzAb+nRBU9cMm5VcD17qWLwbeU9XDrn3fA6aJyEqgu6p+7lr/AjAHOGVCSE5OJjMz042QjTHGnCAie90p506VUQKQ3+R9gWvdydzEf77YT7ZvgmvZ3WMaY4zpYO7cIbRUt99iPZOIXEtj9dAFrezblmPeAtwCkJSU1FqsxhhjTpM7dwgFQN8m7xOBouaFRGQy8FNglqpWt7JvgWv5lMcEUNWnVDVdVdNjY1utAjPGGHOa3EkIa4EUEekvIiHAXGBx0wIiMgZ4ksZkUNJk03Jgqoj0FJGewFRguaoWAxUiMsHVu+h64M12uB5jjDGnqdUqI1WtE5H5NH65BwLPqGqOiNwPZKrqYuAPQDfgX67eo/tUdZaqHhaRX9GYVADuP9HADNwGPAd0obHN4ZQNysYYYzpWq91OvUl6erpaLyNjjGkbEVmnqumtlbMnlY0xxgDu9TIyxvi4ypo6Vucd4sDRao5V13Gsup74yDBGJESS0rsbwYH2t6GxhGCM36qtb+CNDYW8s3k/n+QepLquocVyoUEBTBoSy03nDWB8ck9sFJnOyxKCMX5GVflwewm/XrqVvNJjJPTowtVnJTFlWG/6x4YTHhpEl+BA8g9XsqmwnA37yngjq5DlOQcYlRjJHZNTuHBob6cvwzjAGpWN8SMlR6v4wavZrNpRyoCYcH5yyTAuGtar1b/6j9fU89r6Ap7+ZDe7Dx7jsrGJ3DszlcguwR6K3HQkdxuVLSEY4yc2F5Zz8wuZlFXWctfUwVx/djIhQW1rG6ipa+CxFTv5y8pd9IoI5Y9XjOacQTEdFLHxFOtlZEwnsjS7mMv/+hkCvHrb2fzvxAFtTgYAIUEB3DV1CK/fdg5dQwK5/pkveHVdQes7Gr9gCcEYH7dwbT7feWk9qfHdeXP+eQzvE3nGxxzdtwdvfOdczhoQxQ/+tZG/rMzFl2oTzOmxhGCMD3tn837ufj2biSkxvHTzBGIjQtvt2BFhwTx7Ywaz0/rw+3e286slWy0p+DnrZWSMj/ps10G+9/IGRvftwZPXjSMsOLDdzxESFMDDV6YRFR7CM5/uJiIsiO9PGdzu5zHewRKCMT5o2/6j3PLCOpJjuvLsjePpGtJxH+WAAOHeS1P5sqqOP32wk5iIUK6b0K/1HY3PsYRgjI85Vl3Ht19cT9eQQF74n7Po0TWkw88pIvz2GyM5UlnDvW9uJqprCDNGxXf4eY1nWRuCMT5EVfnZG5vZc/AYf5o7hrjIMI+dOygwgMfmjWVcUk++vzCLzYXlHju38QxLCMb4kIWZ+SzaUMjtFw3m7IHRHj9/l5BAnrxuHDHhIdz6z3WUVdZ4PAbTcSwhGOMjdhyo4L7FOZw7KJr5Fw5yLI7obqH85dpxlByt5vYFWTQ0WM8jf2EJwRgfUN+g/OjVbLqGBPHwVWkEBjg7AF1a3x7cOzOVj3aU8sgHOx2Nxd+VVlTz3pYDVNbUdfi5LCEY4wNeXLOXrPwyfn7pMHpFeK7d4FSuOSuJb4xN4LEVO1mdd8jpcPzWh9tLuPmFTAqPHO/wc7mVEERkmohsF5FcEbm7he3ni8h6EakTkcubrP+aiGQ1eVWJyBzXtudEZHeTbWntd1nG+I/i8uP8/p3tTEyJYU5agtPh/JuI8KvZI0iK6spdCzdSUVXrdEh+afWuQ0SHhzCoV7cOP1erCUFEAoHHgelAKjBPRFKbFdsH3Ai81HSlqn6oqmmqmgZcCFQC7zYp8sMT21U16/Qvwxj/dd+bOdQ1NPDAnJFeN1dBeGgQD12ZRnH5ce5/a4vT4fgdVWV13iEmDIj2yP+9O3cIGUCuquapag2wAJjdtICq7lHVbKDlGTgaXQ68raqVpx2tMZ3M8pz9vLvlAHdMHkxSdFenw2nRuH49+fakQfxrXQHLc/Y7HY5fyT98nKLyKiYMiPLI+dxJCAlAfpP3Ba51bTUXeLnZugdEJFtEHhaR9huExRg/UFPXwG+WbSWlVzduOq+/0+Gc0vcuSmF4n+7c8/omDh+zrqjt5fO8gwAe62LsTkJo6T6lTf3MRCQeGAksb7L6HmAoMB6IAn58kn1vEZFMEcksLS1ty2mN8WkvrtnL3kOV/OSSYV4/53FIUAAPXZlGRVUtv15qVUftZXXeYWK6hTAwtuPbD8C9hFAA9G3yPhEoauN5rgQWqeq/W51UtVgbVQPP0lg19RWq+pSqpqtqemxsbBtPa4xvKj9ey6Mf7OTcQdFMGuIbv/dD4iL41vkDeX19IZ/sPOh0OD7vRPvBWR5qPwD3EsJaIEVE+otICI1VP4vbeJ55NKsuct01II1XOgfY3MZjGuO3/rIyl7LjtdwzfZjXNSSfyvwLB9E/JpyfLNrE8Zp6p8PxaXsPVVJcXsWEAZ57Ir3VhKCqdcB8Gqt7tgILVTVHRO4XkVkAIjJeRAqAK4AnRSTnxP4ikkzjHcZHzQ79oohsAjYBMcCvz/xyjPF9BUcqefbTPXw9LYERCWc+2Y0nhQUH8sDXR7DvcCWPrrAH1s7EiWc7zvZgQnBrtFNVXQYsa7bu3ibLa2msSmpp3z200Aitqhe2JVBjOotH3m/8Ir3r4iEOR3J6zhkYw+XjEnlqVR5z0hIYEhfhdEg+aXXeIWK6hTIwNtxj5/TulipjOpl9hypZtKGQa85KIqFHF6fDOW0/uWQY3UKD+MXiHJtl7TSoKp/nHWLCgCiPVhlaQjDGi/xlZS6BAcKtFwx0OpQzEhUewg+mDubzvEO8vdmeTWirPYcqOXC02qPtB2AJwRivkX+4klfXFTBvfF96d/eO8YrOxNVn9WNYfHceWLrVGpjb6PNdje0HlhCM6aSe+GgXASLcOsm37w5OCAwQfjlrOIVlx3liZa7T4fiUj3eWEtc9zKPtB2AJwRivUFR2nH9l5nPl+ETiI3237aC5jP5RzE7rw19X5ZF/2EatcUddfQOf5h7k/MExHu9ybAnBGC/w5Ee7ALhtknMT33SUe6YPI1CE372zzelQfEJ2YTlHq+qYmOL5BxItIRjjsLLKGhZmFjA7LcGnexadTFxkGLecP4Cl2cWs23vE6XC83qodpYjAeYNiPH5uSwjGOOzFNfs4XlvP/0707gHszsS3LhhAr4hQfrVki3VDbcXHOw8yKiGSnuEhHj+3JQRjHFRdV89zn+1hYkoMQ+O6Ox1Oh+kaEsQPLx5CVn4Zb2UXOx2O1yo/XktWfhnnD3Zm/CpLCMY4aHFWEaUV1dw8cYDToXS4y8YmkhrfnQff3kZVrXVDbcnnuw5S36COtB+AJQRjHKOqPP3JbobGRTAxxfP1xZ4WECD87NJhFJYd5/nP9jgdjldatfMg3UKDGJPUw5HzW0IwxiEf7zzItv0V3HRef58a0fRMnDMwhq8NieXxD3Mpq7SJdJpSVVbtKOXsgdGOzX9hCcEYhzz9yW5iI0KZldbH6VA86sfTh1JRXcfjH9rDak3tOVRJwZHjjrUfgCUEYxyx5+AxPtpRyjVnJREaFOh0OB41NK47l41N5PnP9lJwxB5WO2HVjsYZIc93sPrQEoIxDnhxzV6CAoR5GUlOh+KIO6cMRgQeeneH06F4jQ+2ldA/Jpx+0Z4drqIpSwjGeFhVbT0LMwu4eHicXwxidzr69OjCjecmsyirkJyicqfDcVxFVS2f7zrIlNTejsZhCcEYD1u8sYjy47VcO6Gf06E46tuTBtE9LJjfv7Pd6VAc9/HOg9TWK5OH+UBCEJFpIrJdRHJF5O4Wtp8vIutFpE5ELm+2rV5EslyvxU3W9xeRNSKyU0Recc3XbIzf++fqvaT06saEAVFOh+KoyC7BfHvSQD7aUfrv4Z47q/e3HKBn12DGOtTd9IRWE4KIBAKPA9OBVGCeiKQ2K7YPuBF4qYVDHFfVNNdrVpP1DwIPq2oKcAS46TTiN8anbMwvI7ugnOvO7tdpupqeyg3nJBMfGcaD72zrtENa1NU3sGJ7CV8b2osgh7qbnuDO2TOAXFXNU9UaYAEwu2kBVd2jqtlAgzsnlcZPwoXAq65VzwNz3I7aGB/1j9V7CQ8J5OtjvjLNeKcUFhzIHZNTyMovY3nOAafDcUTm3iOUVdYyxeHqInAvISQA+U3eF7jWuStMRDJFZLWInPjSjwbKVLWutWOKyC2u/TNLS0vbcFpjvEv58VqWZBcxe0wCEWHBTofjNS4bm8jA2HD+sHwbdfVu/U3pV97fcoCQwAAmOvj8wQnuJISW7mvbcm+XpKrpwNXAIyIysC3HVNWnVDVdVdNjY53/BzPmdC3eWERVbQNzx/d1OhSvEhQYwA8vHsqu0mO8tr7A6XA8SlV5b+sBzh4YTbfQIKfDcSshFABNf4MTgSJ3T6CqRa6fecBKYAxwEOghIif+Bdp0TGN80cK1+QyL787IhEinQ/E6Fw/vzZikHjz83s5ONfDdrtIv2XuokskOdzc9wZ2EsBZIcfUKCgHmAotb2QcAEekpIqGu5RjgXGCLNrYefQic6JF0A/BmW4M3xlfkFJWzqbCcq9ITrTG5BSLCj6cNZf/Rqk418N27WxrbTSYP6+VwJI1aTQiuev75wHJgK7BQVXNE5H4RmQUgIuNFpAC4AnhSRHJcuw8DMkVkI40J4HequsW17cfAnSKSS2ObwtPteWHGeJOFa/MJCQpgjjUmn9SEAdFMGhLLX1buoryy1ulwPGLZpmJGJUZ6zTzablVaqeoyYFmzdfc2WV5LY7VP8/0+A0ae5Jh5NPZgMsavVdXWs2hDIdNHxNGjqz1ucyo/ungoMx77mCc+2sXd04c6HU6H2n3wGJsLj/KzGcOcDuXf7EllYzrYO5v3c7SqjqvSrTG5Nal9ujN7dB+e/XQ3+8urnA6nQy3Z2NhsOmNUvMOR/IclBGM62Ctr80mK6sqEAdFOh+IT7po6hAZVHnnfvwe+eyu7iIzkKK+pLgJLCMZ0qPzDlXyed4jLxyUSEGCNye7oG9WVayf0Y2FmPjsPVDgdTofYvr+CHQe+5NLR3nN3AJYQjOlQb2woBLAnk9vouxemEB4SxIN+OvDdWxuLCBCYPsISgjGdgqry+oZCJgyIom9UV6fD8SlR4SHcOmkg7289wBe7DzsdTrtSVZZkF3HOwBhiI0KdDue/WEIwpoOs31fG7oPHuGzsVzrgGTf8z7n9iesexm+WbfWrge82Fx5lz6FKZnpZdRFYQjCmw7y2voAuwYFMH+l9H3xf0CUkkDunDCYrv4y3N+93Opx2s3hjIcGBwsXD45wO5SssIRjTAapq61mysYhpI+K8YowaX3XZuESG9I7gwXe2UV3n+0Na1NQ1sGhDIZOG9PLKZ1IsIRjTAT7YWsLRqjq+MdYak89EYIDw0xnD2Huokhc+2+t0OGdsxbYDHPyyhnkZ3vlMiiUEYzrAa+sLiOsexjkDY5wOxeedPziWSUNieXTFTg4fq3E6nDPy8hf5xHUP44LB3jF2UXOWEIxpZwe/rOajHaXMGZNAoD170C5+eskwKmvq+ZMPP6xWcKSSVTtLuTI90Wt/LywhGNPOlm0qpr5BmTOmj9Oh+I2U3hFcnZHEP9fsI7fENx9W+1dm41wPV3jxECaWEIxpZ29sKGRoXARD47o7HYpfuWNyCl1DAvnVEt/rhlrfoPwrM5/zBsV49TMplhCMaUf7DlWyfl8Zs9OsMbm9RXcL5Y7Jg/loRynvbfGt+ZdX7SilqLyKeRlJTodySpYQjGlHb2Y1DlUxK82qizrC9Wf3Y3Dvbty/ZItPzaz24pp9RIeHMHmYd8yMdjKWEIxpJ6rKG1mFZCRHkdDDe0aw9CfBgQH8ctYICo4c54mVu5wOxy25JV/ywbYDzMtIIiTIu79y3YpORKaJyHYRyRWRu1vYfr6IrBeROhG5vMn6NBH5XERyRCRbRK5qsu05EdktIlmuV1r7XJIxzsgpOsqu0mPMtsbkDnX2wGhmju7DEx/tYt+hSqfDadXfP84jODCAG89NdjqUVrWaEEQkEHgcmA6kAvNEJLVZsX3AjcBLzdZXAter6nBgGvCIiPRosv2HqprmemWd5jUY4xXezGockmCGDVXR4X56yTCCAoRfvJXj1Q3MJRVVvL6+kCvGJRLTzbsGsmuJO3cIGUCuquapag2wAJjdtICq7lHVbKCh2fodqrrTtVwElACx7RK5MV6kvkFZvLGICwZ755AE/iYuMow7pwxmxbYSlmQXOx3OST336R5qGxq4eeIAp0NxizsJIQHIb/K+wLWuTUQkAwgBmlb8PeCqSnpYRLw/fRpzEmv3HObA0WprTPagG89JZlRiJL9YnMMRL3yC+cvqOv6xei/TR8SRHBPudDhucSchtPRIXZvu0UQkHvgH8E1VPXEXcQ8wFBgPRAE/Psm+t4hIpohklpaWtuW0xnjMWxuL6BIcyORh3jkkgT8KCgzgd98YRfnxWn69dKvT4XzFgi/2UVFVx7fOH+h0KG5zJyEUAE0frUsEitw9gYh0B5YCP1PV1SfWq2qxNqoGnqWxauorVPUpVU1X1fTYWKttMt6ntr6BtzfvZ3Jqb7qG2MimnpTapzvfumAAr60v4OOd3vMH47HqOp5alceEAVGM7tuj9R28hDsJYS2QIiL9RSQEmAssdufgrvKLgBdU9V/NtsW7fgowB9jclsCN8Raf7TrE4WM1zBxljclO+O6FKQyIDefu1zZRfrzW6XAA+NvHeZRUVPPDi4c6HUqbtJoQVLUOmA8sB7YCC1U1R0TuF5FZACIyXkQKgCuAJ0Ukx7X7lcD5wI0tdC99UUQ2AZuAGODX7XplxnjIWxuLiAgL4oIhdgfrhLDgQP54xWj2H63i529sdrzXUcnRKp78KI8ZI+MZ16+no7G0lVv3t6q6DFjWbN29TZbX0liV1Hy/fwL/PMkxL2xTpMZ4oeq6epZv3s/FI+IIDQp0OpxOa0xST+64KIU/vreDSUNi+YaD05Y+9N4O6hoa+NG0IY7FcLq8+7E5Y7zcR9tLqaiuY+Zo613ktG9/bRAZyVHc+2aOYw+sbd9fwcLMfK4/O5l+0b7Rs6gpSwjGnIG3souJCg/hnIHRTofS6QUGCA9dNRoR+O6CDR6fclNV+fXSLXQLDeK7Fw7y6LnbiyUEY05TZU0d7285wPQRcQQH2kfJGyT27MrvLxvFxvwyfrrIs+0JC9bm8/HOg9w1dYjPPpxov8XGnKYV20o4XlvPDOtd5FWmj4znexel8Oq6Ap75dI9Hzrnn4DF+tWQLE1NiuG5CP4+csyNYQjDmNC3NLiamWyhn9bfqIm9zx0UpXDy8Nw8s3cKqHR37fEJdfQN3LswiODCAP1w+mgAvnR7THZYQjDkNx6rrWLGthEtGxnnt/LidWUCA8NCVaQzuHcF3XlzPxvyyDjvXEyt3sX5fGb+eM4K4yLAOO48nWEIw5jR8sK2E6roGG9nUi4WHBvHMjePpER7MdU+vYVNBebuf4+1NxTz8/g5mje7jFz3NLCEYcxqWZhfRKyKU9OQop0Mxp9CnRxdevnkC3bsEc+3Ta9hc2H5J4ZOdB7l9QRZjknry4GWj2u24TrKEYEwbfVldx4fbS7lkZLxVF/mAxJ5defnmCXQLDeLqv63mw20lZ3zMDfuOcMs/MhkQG84zN4ynS4h/PJRoCcGYNvpg6wFq6hq41HoX+Yy+UV1ZcMsEEnp25ZvPreWhd7dT33B6XVLf2VzM9c98QWxEKC/8TwaRXYPbOVrnWEIwpo2WZBcT1z2MsUm+NU5NZ9c3qiuLvn0OV4xL5NEVudzwzBfsPFDh9v5VtfX8dNEmbv3nevrHhPPi/55Fr+6+3YjcnI3Va0wbVFTV8tH2Uq6d0M+nuxd2VmHBgfzhitGM69eTXy3ZwtRHVjF7dB9unzyY/ieZxObL6joWZxXx9Cd57Co9xi3nD+AHU4cQEuR/f09bQjCmDd7feoCa+gZmjIpzOhRzBuZmJDF1eBxPrtrF85/t4Y2sIgbEhpPerycjE3tQW9dAWWUNBUeO807Ofipr6hnSO4LnvjmeSUP8dxIkSwjGtMHS7P3ER4Yxpq9VF/m6qPAQ7pk+jJvO689r6wrJ3HOYd7ccYGFmAQAiENU1hEtHxTM3I4kxfXvQOH2L/7KEYIybKqpqWbWzlGvPsuoif9IrIozbJg0EBtLQoOw/WkXXkEAiwoI7XS8ySwjGuOmDrSXU1Fl1kT8LCBD69OjidBiO8b9WEWM6yNJNjb2LrLrI+Cu3EoKITBOR7SKSKyJ3t7D9fBFZLyJ1InJ5s203iMhO1+uGJuvHicgm1zEfFX+vnDM+raKqlo92lDJ9ZJxVFxm/1WpCEJFA4HFgOpAKzBOR1GbF9gE3Ai812zcKuA84C8gA7hORE39ePQHcAqS4XtNO+yqM6WD/ri6ysYuMH3PnDiEDyFXVPFWtARYAs5sWUNU9qpoNNDTb92LgPVU9rKpHgPeAaSISD3RX1c+1cQaLF4A5Z3oxxnSUE9VF9jCa8WfuJIQEIL/J+wLXOnecbN8E13KrxxSRW0QkU0QyS0s7dlxzY1ryZXWdVReZTsGdhNDSJ8DdQUBOtq/bx1TVp1Q1XVXTY2Nj3TytMe3nxNhFVl1k/J07CaEA6NvkfSJQ5ObxT7ZvgWv5dI5pjEctzS6md/dQqy4yfs+dhLAWSBGR/iISAswFFrt5/OXAVBHp6WpMngosV9VioEJEJrh6F10PvHka8RvTob6srmPljlKmj4i36iLj91pNCKpaB8yn8ct9K7BQVXNE5H4RmQUgIuNFpAC4AnhSRHJc+x4GfkVjUlkL3O9aB3Ab8HcgF9gFvN2uV2ZMO1ixrbF30SVWXWQ6AbeeVFbVZcCyZuvubbK8lv+uAmpa7hngmRbWZwIj2hKsMZ62LLu4cWa0flZdZPyfPalszEkcq67jw+0lTB9hvYtM52AJwZiTWLGthGqrLjKdiCUEY05i2aZiYiNCSU+OcjoUYzzCEoIxLais+U91UWcbAtl0Xjb8tRuqaus5cLSK/eVVlB+vZWCvbvSPDrd6ZT+2YlsJVbUNTB9h1UWm87CEcBKqyqe5h3jusz2s2HaAhmbPUUeEBjEyMZK5GUnMGBlvf0X6mWWbionpFkpGf6suMp2HJYQWfJp7kF8szmFnyZdEh4dw88QBDOrVjbjIMLqFBrHzwJdkF5bxae4hvvfyBh55fwfzvzaI2WkJlhj8QGVNHSu2lXDFuL72/2k6FUsITTQ0KH/+MJeH399B/5hw/njFaGaMiicsOPC/yo1J6smV4/vS0KC8vXk/j63YyZ0LN/La+gIeuWoMsRGhDl2BaQ8nqotmjLLqItO5WEJwKaus4Y5Xsli5vZQ5aX34zTdG0jXk1P88AQHCjFHxTB8Rx8LMfO5bnMMlj37MY/PGMGFAtIciN+3tRHXReOtdZDoZ62VEYxXBDc+u5bPcQ/xqzggeviqt1WTQVECAMDcjiTe+cy4RoUFc/bfVvLhmbwdGbDrKieoi611kOqNOnxDq6huY/9IGNhWU8djVY7huQj9OdzbPYfHdWfzd85g0pBc/XbSZf6y2pOBrrLrIdGadOiGoKj9ZtIkV20q4f/YILh4ed8bH7BYaxBPXjmXysF78/A1LCr7GqotMZ9apE8JfP8pjYWYB371wENdO6Nduxw0NCuTxa/6TFBauzW99J+M4qy4ynV2nTQhbio7y0HvbuWRkHHdOGdzuxz+RFCamxPCTRZtYk3eo3c9h2teJ6iIbu8h0Vp0yIdTUNXDnwiwiu4TwwJyRp91m0JrQoED+fPVYkqK7ctuL68k/XNkh5zHtY8nGxrGL7GE001l1yoTw2IqdbNtfwW+/MZKe4SEdeq7ILsH8/fp06uobuPmFTI5V13Xo+czp+dI11LU9dW46M7cSgohME5HtIpIrIne3sD1URF5xbV8jIsmu9deISFaTV4OIpLm2rXQd88S2Xu15YSeTXVDGX1bu4rKxiUxJ7e2JUzIgtht/vnosOw5U8OPXslHV1ncyHvX+lgNU1zVwqfUuMp1YqwlBRAKBx4HpQCowT0RSmxW7CTiiqoOAh4EHAVT1RVVNU9U04Dpgj6pmNdnvmhPbVbWkHa7nlBoalJ+9sZnYbqHcO7P5JXSs8wfHctfUISzJLub19YUePbdp3ZLsYuK6hzE2yWZGM52XO3cIGUCuquapag2wAJjdrMxs4HnX8qvARfLVivl5wMtnEuyZWrqpmOyCcn548RAiuwR7/Py3XjCQjP5R3PvmZvYeOubx85uWlR+vZdWOUmaMircRbE2n5k5CSACa9psscK1rsYyq1gHlQPOxG67iqwnhWVd10c9bSCDtqqaugT8s387QuAjmjGkevmcEBggPX5VGQIBw+4IsausbHInD/Lf3thygpt6qi4xxJyG09EXdvBL8lGVE5CygUlU3N9l+jaqOBCa6Xte1eHKRW0QkU0QyS0tL3Qi3ZS+t2cu+w5X8ePpQRxsNE3p04TdfH0lWfhl/XpHrWBzmP5ZkF5HYswtpfXs4HYoxjnInIRQAfZu8TwSKTlZGRIJhAYwzAAAUp0lEQVSASOBwk+1zaXZ3oKqFrp8VwEs0Vk19hao+parpqpoeGxvrRrhfVVFVy6Mrcjl7QDSTBp/eMdrTzNF9mJPWh8c/zGXb/qNOh9OplVXW8MnOg8wYFd9h3Y+N8RXuJIS1QIqI9BeREBq/3Bc3K7MYuMG1fDmwQl1daUQkALiCxrYHXOuCRCTGtRwMXApspoP8bVUeh4/VcM8lQ73mQ3/vzOF07xLM3a9tor757DvGY5bn7KeuQbl0ZB+nQzHGca0mBFebwHxgObAVWKiqOSJyv4jMchV7GogWkVzgTqBp19TzgQJVzWuyLhRYLiLZQBZQCPztjK/mJMqO1zI7rQ+jEr2nSiAqPIR7L00lK7+Mf3y+x+lwOq3FG4tIju7KiITuTodijOPEl/rEp6ena2Zm5mnt29CgXteDRFW54dm1rNtzmHfvvICEHl2cDqlTKamoYsJvPmD+1wZx59QhTodjTIcRkXWqmt5auU7zpLK3JQMAEeGBOSNoULjvzQ6rMTMnsTS7mAaFWWlWXWQMdKKE4K36RnXljskpvL+1hA+3dfizeaaJtzYWMSy+O4N6RTgdijFewRKCF/jmuf0ZEBvO/Uu2UF1X73Q4nUL+4UrW7ytj5mh79sCYEywheIGQoADumzmc3QeP8cwne5wOp1N4K7ux5/TMUVZdZMwJlhC8xAWDY5mS2pvHVuxkf3mV0+H4vcVZRYxN6kHfqK5Oh2KM17CE4EV+PiOVugblt29vdToUv7bzQAXb9lcwa7TdHRjTlCUEL5IU3ZVbJg7gzawisvLLnA7Hby3eWESAwCU2dpEx/8USgpe5ddJAYrqF8MDSLTZvQgdQVRZtKOTcQTH0ighzOhxjvIolBC/TLTSI708ZzNo9R1iec8DpcPzOur1HKDhynK87NOKtMd7MEoIXuiq9Lym9uvG7t7dSU2dDZLenRRsK6RIcyMXD45wOxRivYwnBCwUFBvCTS4ax51AlL67Z63Q4fqOmroEl2cVMHd6b8NAgp8MxxutYQvBSk4bEcu6gaP70wU6OVtU6HY5fWLm9hPLjtY5NkGSMt7OE4KVEhHumD6OsspanPsprfQfTqjeyCokOD2HioBinQzHGK1lC8GIjEiK5dFQ8T3+ym5Kj9rDamSg/Xsv7W0uYOboPQYH2a29MS+yT4eV+MHUItfUNPLpip9Oh+LR3NhdTU9dgvYuMOQVLCF4uOSacuRl9WfBFPnsOHnM6HJ/12rpCBsSGMyox0ulQjPFalhB8wPcuSiE4MID/e3e706H4pD0Hj/HFnsNcPi7Ra6ZQNcYbuZUQRGSaiGwXkVwRubuF7aEi8opr+xoRSXatTxaR4yKS5Xr9tck+40Rkk2ufR8U+qSfVKyKMm87rz5LsYjYXljsdjs95fX0BAQLfGJPodCjGeLVWE4KIBAKPA9OBVGCeiKQ2K3YTcERVBwEPAw822bZLVdNcr1ubrH8CuAVIcb2mnf5l+L+bzx9AZJdg/mh3CW3S0KC8tr6Q81JiiYu0oSqMORV37hAygFxVzVPVGmABMLtZmdnA867lV4GLTvUXv4jEA91V9XNtHLDnBWBOm6PvRCK7BHPrBQP5cHspa/ccdjocn/F53iEKy45zxTi7OzCmNe4khAQgv8n7Ate6Fsuoah1QDkS7tvUXkQ0i8pGITGxSvqCVY5pmbjwnmdiIUP7wznYb+M5Nr64rICIsiCmpvZ0OxRiv505CaOkv/ebfRicrUwwkqeoY4E7gJRHp7uYxGw8scouIZIpIZmlpqRvh+q8uIYF878JBfLHnMKt2HnQ6HK93tKqWtzcXM2t0H8KCA50Oxxiv505CKAD6NnmfCBSdrIyIBAGRwGFVrVbVQwCqug7YBQx2lW96D9/SMXHt95SqpqtqemxsrBvh+rerxieR2LMLf1i+ze4SWrEsu5iq2gYut+oiY9ziTkJYC6SISH8RCQHmAoublVkM3OBavhxYoaoqIrGuRmlEZACNjcd5qloMVIjIBFdbw/XAm+1wPX4vJCiA708ezObCo7yzeb/T4Xi1VzLzGdSrG2l9ezgdijE+odWE4GoTmA8sB7YCC1U1R0TuF5FZrmJPA9Eikktj1dCJrqnnA9kispHGxuZbVfVEi+htwN+BXBrvHN5up2vye3PGJDCoVzf+793t1DfYXUJLthYfZcO+MuaO72vPHhjjJrfGAFbVZcCyZuvubbJcBVzRwn6vAa+d5JiZwIi2BGsaBQYId00ZzG0vrmfRhkKrEmnBgi/2ERIYwGVj7d/GGHfZk8o+atqIOEYkdOeR93fYJDrNHK+pZ9GGQqaPjKNneIjT4RjjMywh+CgR4QdTh1Bw5DivrN3ndDheZdmmYo5W1TEvI8npUIzxKZYQfNgFg2MZn9yTx1bkcrym3ulwvMbLX+xjQEw4Z/WPcjoUY3yKJQQfJiL88OKhlFRU8/zne5wOxyvsOFBB5t4jzMtIssZkY9rIEoKPy+gfxaQhsTyxchflx22qzZfWuBqTraHdmDazhOAHfjB1COXHa/nbqs491eax6jpeW1fAtBFxRFljsjFtZgnBD5yYavOZT3dTWlHtdDiOeX1DIRXVddxwTrLToRjjkywh+Im7pg6huq6Bxz/MdToUR6gqz3+2h5EJkYxNsieTjTkdlhD8RP+YcK5MT+SlNfvIP1zpdDge92nuIXJLvuTGc5KtMdmY02QJwY/cftFgROCh93Y4HYrHPffZHqLDQ7h0dLzToRjjsywh+JG4yDC+eW5/3sgqJKeo80y1mX+4kg+2HWBeRhKhQTbMtTGnyxKCn7lt0kAiuwTz4DudZ6rNFz7fQ4AI10ywJ5ONOROWEPxMZJdg5n9tEKt2lPJprv9PolNRVcuCtflMGxFHfGQXp8MxxqdZQvBD107oR0KPLvz27a00+Pnw2C+t2UdFVR3fOn+A06EY4/MsIfihsOBA7praOInO4o0tTkTnF6rr6nn6k92cMzCaUYnW1dSYM2UJwU/NSUtgREJ3fv/ONr8d+O6NDYWUVFRz6wUDnQ7FGL9gCcFPBQQIP5+RSlF5FX//2P+GtGhoUJ5clcfwPt2ZmBLjdDjG+AW3EoKITBOR7SKSKyJ3t7A9VERecW1fIyLJrvVTRGSdiGxy/bywyT4rXcfMcr16tddFmUZnDYhm2vA4nvhoFweOVjkdTrt6b+sB8kqP8a0LBtqDaMa0k1YTgogEAo8D04FUYJ6IpDYrdhNwRFUHAQ8DD7rWHwRmqupI4AbgH832u0ZV01yvkjO4DnMS91wylLp65f+W+083VFXliZW76BvVhUtGxDkdjjF+w507hAwgV1XzVLUGWADMblZmNvC8a/lV4CIREVXdoKonWjVzgDARCW2PwI17+kWHc+O5yby6voBNBf7xsNrKHaVk5Zdx6wUDCQq0Wk9j2os7n6YEIL/J+wLXuhbLqGodUA5ENytzGbBBVZsOx/msq7ro53KS+34RuUVEMkUks7S01I1wTXPzLxxEdHgI9y7e7PPdUFWVh97dQd+oLlwxrq/T4RjjV9xJCC19UTf/VjllGREZTmM10reabL/GVZU00fW6rqWTq+pTqpququmxsbFuhGua6x4WzN3Th7FhXxmvritwOpwz8u6WA2wqLOd7F6YQEmR3B8a0J3c+UQVA0z/FEoHmndv/XUZEgoBI4LDrfSKwCLheVXed2EFVC10/K4CXaKyaMh3ksrEJjE/uye/e2UZZZY3T4ZyWhobGu4MBMeF8fUzzm1RjzJlyJyGsBVJEpL+IhABzgcXNyiymsdEY4HJghaqqiPQAlgL3qOqnJwqLSJCIxLiWg4FLgc1ndinmVESE+2ePoPx4LX/w0QbmpZuK2X6ggtsnp1jbgTEdoNVPlatNYD6wHNgKLFTVHBG5X0RmuYo9DUSLSC5wJ3Cia+p8YBDw82bdS0OB5SKSDWQBhcDf2vPCzFcNi+/O9Wf346Uv9rExv8zpcNqktr6Bh9/fwZDeEcwc1cfpcIzxS6LqO42M6enpmpmZ6XQYPu1oVS2T//gRUeEhLJ5/ns/Uwz/zyW7uX7KFv1+fzuTU3k6HY4xPEZF1qpreWjnf+DYw7aZ7WDAPfH0k2/ZX8MTKXa3v4AUOfVnNw+/vYGJKDBcNs+cXjekolhA6oSmpvZk1ug9//nAn2/dXOB1Oqx56bweVNfXcNzPVnko2pgNZQuik7puZSkRYMD96dSN19Q1Oh3NSW4qO8vIX+7j+7H4M6hXhdDjG+DVLCJ1UdLdQfjlrOBsLynlylXcOfqeq/PKtHHp0DeGOiwY7HY4xfs8SQid26ah4ZoyK56H3drB+3xGnw/mKBWvzWbP7MD+8eAiRXYOdDscYv2cJoRMTEX7z9ZHEdQ/jey9v4GhVrdMh/VvBkUp+vWQL5w6K5qp0G6LCGE+whNDJRXYJ5tF5Yygur+Inr2/CG7ohNzQoP3o1G4AHLxtFQIA1JBvjCZYQDOP69eTOKYNZkl3MgrX5re/QwV78Yh+f7TrEzy5NJbFnV6fDMabTsIRgALj1goFMTInh3jc3szrvkGNx5JZU8NtlW5mYEsPc8VZVZIwnWUIwAAQGCH++eix9o7py6z/XsffQMY/HUH68lptfWEfXkCD+cPloe+bAGA+zhGD+LbJLMM/cMB6A/3lurUcbmesblNsXbKDgSCV/vXYscZFhHju3MaaRJQTzX5JjwnnimnHsPVTJ/z6fSWVNnUfO+9B721m5vZT7Zg4nPTnKI+c0xvw3SwjmK84eGM1DV6WRuecwNz67lmPVHZsUXvh8D49/uIt5GX255qykDj2XMebkLCGYFs0a3YdH5o4hc89hvtmBSeHZT3dz75s5TEntzS9njbB2A2McZAnBnNSs0X3409wxrNt3hKv/tprCsuPtevy/f5zHL9/awrThcTx+9VifGYrbGH9ln0BzSjNH9+HJa8eRV3qMSx/9mE92HjzjY1bV1vPzNzbz66VbmTEynseuHmPJwBgv4NanUESmich2EckVkbtb2B4qIq+4tq8RkeQm2+5xrd8uIhe7e0zjPSan9ubN+efSKyKM655Zwx/f3X7ajc3b91cw+8+f8o/Ve7l5Yn/+NDeNYJsO0xiv0OqMaSISCOwApgAFNM6xPE9VtzQp821glKreKiJzga+r6lUikgq8DGQAfYD3gRPDVp7ymC2xGdOcVVlTx8/e2Mzr6wvpFRHKnVMGc/m4RLfmN95fXsXTn+Txwud7iQgL4o9XpnHB4FgPRG2McXfGtCA3jpUB5KpqnuvAC4DZQNMv79nAL1zLrwJ/lsbWwdnAAlWtBna75lzOcJVr7ZjGy3QNCeKhK9O4OiOJ3yzbyt2vb+KxFblcPDyOKam9GZ/c87+Sw4GjVWzML+P9rQdYtKGQ+gZl5ug+/HTGMHpF2HMGxngbdxJCAtB0gJsC4KyTlVHVOhEpB6Jd61c32zfBtdzaMY2XSk+O4rXbzmF5zn5eWZvPP9fs5ZlPdxMYIHQLDaJbaBA19Q2UVlQDEBoUwLyMJG6eOIC+UTY2kTHeyp2E0FI/wOb1TCcrc7L1LdUxtFh3JSK3ALcAJCVZH3VvISJMGxHPtBHxHKuuY9WOUnKKjvJldR0VVXWIwPA+3RmVGElqfCRdQgKdDtkY0wp3EkIB0HSUsUSg6CRlCkQkCIgEDreyb2vHBEBVnwKegsY2BDfiNR4WHhrE9JHxTB8Z73Qoxpgz4E73jrVAioj0F5EQYC6wuFmZxcANruXLgRXa2Fq9GJjr6oXUH0gBvnDzmMYYYzyo1TsEV5vAfGA5EAg8o6o5InI/kKmqi4GngX+4Go0P0/gFj6vcQhobi+uA76hqPUBLx2z/yzPGGOOuVrudehPrdmqMMW3nbrdTeyLIGGMMYAnBGGOMiyUEY4wxgCUEY4wxLpYQjDHGAD7Wy0hESoG9p7l7DHDmYzf7FrvmzsGu2f+d6fX2U9VWR5P0qYRwJkQk051uV/7ErrlzsGv2f566XqsyMsYYA1hCMMYY49KZEsJTTgfgALvmzsGu2f955Ho7TRuCMcaYU+tMdwjGGGNOoVMkBBGZJiLbRSRXRO52Op6OJCJ9ReRDEdkqIjkicrvTMXmKiASKyAYRWeJ0LJ4gIj1E5FUR2eb6/z7b6Zg6moh83/V7vVlEXhYRv5uLVUSeEZESEdncZF2UiLwnIjtdP3t2xLn9PiGISCDwODAdSAXmiUiqs1F1qDrgLlUdBkwAvuPn19vU7cBWp4PwoD8B76jqUGA0fn7tIpIAfA9IV9URNA6dP9fZqDrEc8C0ZuvuBj5Q1RTgA9f7duf3CQHIAHJVNU9Va4AFwGyHY+owqlqsqutdyxU0fkkknHov3yciicAM4O9Ox+IJItIdOJ/GuUhQ1RpVLXM2Ko8IArq4ZmbsyklmWvRlqrqKxnllmpoNPO9afh6Y0xHn7gwJIQHIb/K+gE7wBQkgIsnAGGCNs5F4xCPAj4AGpwPxkAFAKfCsq5rs7yIS7nRQHUlVC4H/A/YBxUC5qr7rbFQe01tVi6Hxjz6gV0ecpDMkBGlhnd93rRKRbsBrwB2qetTpeDqSiFwKlKjqOqdj8aAgYCzwhKqOAY7RQdUI3sJVbz4b6A/0AcJF5Fpno/IvnSEhFAB9m7xPxA9vM5sSkWAak8GLqvq60/F4wLnALBHZQ2OV4IUi8k9nQ+pwBUCBqp64+3uVxgThzyYDu1W1VFVrgdeBcxyOyVMOiEg8gOtnSUecpDMkhLVAioj0F5EQGhuhFjscU4cREaGxXnmrqj7kdDyeoKr3qGqiqibT+P+7QlX9+i9HVd0P5IvIENeqi2icu9yf7QMmiEhX1+/5Rfh5Q3oTi4EbXMs3AG92xEmCOuKg3kRV60RkPrCcxl4Jz6hqjsNhdaRzgeuATSKS5Vr3E1Vd5mBMpmN8F3jR9YdOHvBNh+PpUKq6RkReBdbT2JtuA374xLKIvAxMAmJEpAC4D/gdsFBEbqIxMV7RIee2J5WNMcZA56gyMsYY4wZLCMYYYwBLCMYYY1wsIRhjjAEsIRhjjHGxhGCMMQawhGCMMcbFEoIxxhgA/h/lQqeXkaivGgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1be1cf41e80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def func(x):\n",
    "    #注意，超出范围的概率要设置为0\n",
    "    w=np.ones_like(x)\n",
    "    w=np.where(x<0,0,w)\n",
    "    w=np.where(x>10,0,w)\n",
    "    return w*((x-2)*(x-2)+(x-5)*(x-5)*(x-5)+100*np.cos(x)+106)/1179\n",
    "x=np.linspace(0,10,100)\n",
    "plt.plot(x,func(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们的$q(x\\rightarrow x')$可以任意选择（但是要包含到区间$(0,10)$），比如均值为$x$，方差为1的建议分布：   \n",
    "\n",
    "$$\n",
    "q(x,x')=\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(x'-x)^2}{2}}\n",
    "$$  \n",
    "\n",
    "接下来采样看看效果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#采样的样本量\n",
    "nums=10000\n",
    "count=0\n",
    "points=[]\n",
    "#采样x0\n",
    "point=np.random.randn()\n",
    "points.append(point)\n",
    "while count<nums:\n",
    "    #按照q(x,x')，采下一个点：均值为0，方差为1的标准正态分布采样+point偏置\n",
    "    new_point=np.random.randn()+point\n",
    "    #alpha(x,x')\n",
    "    alpha=min(1.,func(new_point)/(1e-12+func(point)))\n",
    "    #从(0,1)均匀采样一个u\n",
    "    u=np.random.random()\n",
    "    #判断是否接收新点还是旧点\n",
    "    if u<alpha:\n",
    "        points.append(new_point)\n",
    "        point=new_point\n",
    "    else:\n",
    "        points.append(point)\n",
    "    count+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xd8VFX+//HXZ0pCC6GGjqEjggICIoqrqAiigAiK/iyga9tVV3HXshYQRbH35bsKKBaaKIiIIgooNiBABEINNQklgUAghUxm5vz+SFxjDGRCZuZO+Twfjzyc3Ll37nsk+czJueecK8YYlFJKRQeb1QGUUkoFjxZ9pZSKIlr0lVIqimjRV0qpKKJFXymloogWfaWUiiJa9JVSKopo0VdKqSiiRV8ppaKIw+oAZTVo0MAkJiZaHUMppcLK6tWrDxpjGla0X8gV/cTERJKSkqyOoZRSYUVEdvuyn3bvKKVUFNGir5RSUUSLvlJKRREt+kopFUW06CulVBTRoq+UUlFEi75SSkURLfpKKRVFtOgrpVQUCbkZuUqpIBoXX+pxjnU5VNBoS18ppaKItvSViialW/YqKmnRVyqcafeMqiTt3lFKqSiiRV8ppaKIdu8oFSm0q0f5QIu+UqFEC7cKMO3eUUqpKKJFXymlooh27ygVbnSsvaoCbekrpVQU8anoi8gAEdkiIqki8nA5z48RkY0isk5EvhWR00o9d7OIbCv5utmf4ZVSSlVOhUVfROzAW8BAoBNwnYh0KrPbWqCHMeZMYA7wfMmx9YCxwDlAL2CsiNT1X3yllFKV4UtLvxeQaozZYYxxATOBIaV3MMYsNcbkl3z7C9C85PFlwGJjTLYx5jCwGBjgn+hKKaUqy5ei3wxIK/V9esm2E7kV+PIUj1VKKRVAvozekXK2mXJ3FLkB6AH8pTLHisjtwO0ALVu29CGSUsrvdGJYVPClpZ8OtCj1fXNgb9mdROQS4FFgsDGmsDLHGmPeNsb0MMb0aNiwoa/ZlVJKVZIvRX8V0E5EWolIDDASmF96BxHpBvyX4oKfWeqpRUB/EalbcgG3f8k2pZRSFqiwe8cY4xaRuyku1nZgqjEmRUTGA0nGmPnAC0At4GMRAdhjjBlsjMkWkaco/uAAGG+MyQ7IO1FKKVUhn2bkGmMWAgvLbHui1ONLTnLsVGDqqQZUKmJUts9c+9hVAOiMXKWUiiK69o5SESDH1OR7bxeWec4i2bQl95lvyC/04PYa2ibU4oymtTmjWTyXmzjqyzGr4yoLadFXKozt8Sbwsns4n3vPxYOdeHLpZdtMvfa9qRFrB2DbgVy+3niAmavSeIo3udy2ghsdiznbts3i9MoKWvSVCkOHTByvuq9mhqcfDjyMsi/icvsKzpLtOMQLw8f+YX9jDKmZuXz0+r/5xHMB81znc6ktifHO92giOrYimmjRVyrMJHnb83fXvRyiNiPtS7nXMZcEOXLSY0SEdo3iGOd8n385ZvG+pz+vuYdxaeHz/NMxmxvti7FLuXMuVYTRC7lKhQljDFN+2MlI12NUExfzYx7jaee7FRb8smpKIXc5PufrmIfoZktlnHsUtxc9wDFTPUDJVSjRoq9UGPB6DQ99so6nFmzkIlsy82Meo5NtT5Ves6Utk/edExnveJdl3rMY7hpLmreBnxKrUKVFX6kQ5zXCw5+uY3ZSOvf2a8vbzpeJl/yKD/SBCNzkWMw053PsM/UZ4nqaZG8bv7y2Ck1a9JUKYV4j/Nt9a3HBv7gdY/p3QMpbxrCscfG/f/ngfPsG5sY8QS0p4EbXw6xLr1yXkQofWvSVstpJCvRT7huY6enH3Re15f5L2gU0RhvbPmbEPE285HHjlJWk7C1nFnAlP0xU6NGir1SImuW+kHc9Axlt/5IH+rdHfGriV00zOcQM5wRqxti5YfIKtuw/yUQu/QAIS1r0lQpBa7xtedw9mvNt63nU8VFQCv5vWtiymH5bb5x2G7e8t4qsY4UVH6TChhZ9pUJMpqnDXa77aCSHecP5RvFkqyBLbFCTKTf35FBeIbd/kMTxIk/QM6jA0KKvVAjxGOFu1z0cpQZvO1+iruRalqVL83hevbYra/cc4aFP1mGMTt6KBFr0lQoh73oGstKcznjHe5xuS6v4gAAb0LkJ/7qsA58l7+U/y7ZbHSeiTfhiI498uj7g59Gir1SISPU25Xn3NVxqS2K4/Xur4/zP3y5sw5VnNeWlr7ewwtvR6jgRyRjDgnX7OJLvCvi5tOgrFQKKjJ0xRXdRk0KecU72bSx+kIgIzw7rQst6NbjXdTeHTJzVkSLO9qw89uUc5/x2gZ8RrUVfqRAwyTOYdaYNE5xTaChHrY7zp+GYtWIdvHl9dw4TxwNFd+E1IfSpFAF+2JYFQN+2DQN+Li36Sllsl7cRb7qHcoXtZy63r7Q6zgl1bhbPY44PWebtyjueQVbHiSg/pB6kZb0atKxfI+Dn0qKvlMXGu28khiKecH5gdZQK3WhfzGW2lbzkHsFmbwur40SEIo+XX3ZkB6VrB7ToK2WpbzcdYIm3O/c5Pq30EslWEIFnnFOII58Hiu7EZexWRwp7yWlHyC1007etFn2lIte4eI6PbcCT739BW0nnZvsiqxP5rL4cY4JzCimmFW+6h1odJ+wt33YQm0CfNlr0lYpo73gGscc04knHNJwSXjNeB9iTuMq2nLc8Q1nvbWV1nLD2w7YsujSvQ3wNZ1DOp0VfKQtkmjr8xz2YgbYVnGdPsTrOKRnnfJ8G5BR387iDv1REJMgpKCI57UjQunZAi75SlnjDPZQiHDzkmGl1lFMWL3lMcE5lq2nBO8t3WB0nLP28/RBeQ9Au4oLeGF2poNvtTWCGpx/X2peRaDtgdZwqucS+hss9K3jtWxuXd2lCqwY1/7jU8rhy1uRX//NDahY1Yux0b1k3aOfUlr5SQfaKezgOPNzr+NTqKH4xzjmNWIeNf3+6XhdlqwRjDN9vPUjv1vWJcQSvFGvRVyqQysxs3ehtyWfePoy2f0WjMBii6YsEOcIjA0/n5x2HmLM63eo4YWN7Vh57svO5qGNCUM+rRV+pIHrRfQ1xFHCnY4HVUfxqZM8W9Eqsx4SFmzhsalkdJyws2VzctddPi75SkelXb2uWeLtzh+Nz4iWvcgeH+K0JbTbhqaGdOXbczfPukVbHCQtLNmfSsXEczepUD+p5tegrFSRvuocSTy432RdbHSUgOjSO45bzEpnpuZC13jbFG0P8w8oqOQVFrNp1OOitfNCir1RQbPa2YLG3B6Psi4iTAqvjBMw/LmlPAkd4vOgWPLoS5wkt35aFx2u4+HQt+kpFpLfcQ6hJAaMdX1kdJaBqxTp41PkRG0wrpnsutjpOyFqyKZO6NZx0bRG8oZq/0XH6SgXYDm9jFnh7c4d9AXUq25fvD0HuWrnS9jMzbRfxgvtaBtlXUE+OBfX8oc7jNSzbmsWFHRKw24L/15AWfaUCbJJnMDG4udWx0Ooop6aSHxoi8KRjGgNcE3nJPYIJzql/fp0onrSVnHaE7DyXJf35oN07SgXUXlOPuZ7zuc6+JDTuiBUk7WwZ3GT/mumefqR4T7M6TkhZsvkAdptwQfvA3yWrPFr0lQqgae7LMAh/DddWfhXc5/iEuuTyZNFN6ETd3327KZMep9UlvnpwVtUsS4u+UgGSW+hmuqcfA20raS4HrY4TdPGSz78cs1hpTudz77lWxwkJuw7msXn/MS5Ne82yoaxa9JUKkI+T0jhGzahs5f/mGvsyzpCdPFt0Pfkm1uo4lvtyw34ABlp4L2Sfir6IDBCRLSKSKiIPl/P8BSKyRkTcIjK8zHMeEUku+Zrvr+BKhTKP1zD1x530kC10tW23Oo5l7GIY53yffdTnv+4rrI5juS837KNrizo0k0OWZaiw6IuIHXgLGAh0Aq4TkU5ldtsDjAKml/MSBcaYriVfg6uYV6mw8HXKftKyC6K6lf+bnrYtXGH7mf/zXEmGqW91HMukZeezLj2HgZ0bW5rDl5Z+LyDVGLPDGOMCZgJDSu9gjNlljFkH6O1zlAIm/7CTlvVqcKktyeooIeERZ3F7cGLRdRYnCYITLD3x1W9dO52bWJHqf3wp+s2AtFLfp5ds81U1EUkSkV9ERO+irCLer090Y/Xuw4w+Ogm76LAVgGZyiDvsC/jc24dV3g5Wx7HEwg376NysNi3r17A0hy9Fv7wpY5X5SW5pjOkBXA+8KiJt/nQCkdtLPhiSsrKyKvHSSoWeae7+1KSA4fbvrY4SUu50fE5jDjG+6Ea83uj6MNyXU8DaPUcsb+WDb0U/HWhR6vvmwF5fT2CM2Vvy3x3AMqBbOfu8bYzpYYzp0bChNRMWlPKH7DwXC7y9GWZfHtELq52KGlLIw86ZrDetmbOmkjdbCfPVOn/v2rG2Px98K/qrgHYi0kpEYoCRgE+jcESkrojEljxuAJwHbDzVsEqFulmr0nARw40RunxyVQ2x/Ug32cYLi7aQW+i2Ok7QfLl+Px0bx9G6ofU3mKmw6Btj3MDdwCJgEzDbGJMiIuNFZDCAiPQUkXRgBPBfEUkpOfx0IElEfgWWAhONMVr0VUTyeA0f/rKb3rYU2tsyrI4TkkRgrPN9so4V8tbSVKvjBEX64XxW7spmUBfru3bAxwXXjDELgYVltj1R6vEqirt9yh73E9ClihmVCgtLN2eScaSAR53ayj+ZrrbtDOvejCnLd3Jdz5aWX9gMtHlrixsAQ7tVZvxL4OiMXKX85P1fdtOodiyX2lZbHSXkPTSgIw678MzCTVZHCShjDJ+uyeCcVvVoUe8EH25Bvl6hRV8pP9h5MI/vt2Zxfa/TcIrH6jghr1HtavztwjZ8lbKfn1LDdF0iH4r12rQj7DiYx9Xd/9QRYhkt+kr5wcxXxmDHw8jll1gdJWz8tW9rWtSrzrjPU3B7InNe56dr0qnmtDGwi/Wjdn6jRV+pKnK5vczx/IWLbWtoJEesjhM2qjntPDaoE1sP5PLhL7utjuN3hcbB57/u47IzGhNXzZpllMujRV+pKvpm0wEOEc919iVWRwk7/Ts1om+7Bry8eCvZeS6r4/jVEm83cgqKGBZCXTugRV+pKpuxcg/NyOIC2zqro4QdEWHslZ3Id3l48estVsfxq088fUmIi+W8NqG1yJwWfaWqYM+hfJZvO8g1jmW6zs4papsQx03nJjJj5R42ZETGvXMzTR2WebsytFszHPbQKrOhlUapMDMraQ82gWvs31kdJazdd2k76teM5bF5GyJiXZ7pnn64cXBdr5ZWR/kTLfpKnaIij5fZSelc1CGBJpJtdZywVruak0cHdSQ57QizktIqPiCEuYydj9wXc6EtmVYNalod50+06Ct1ipZuziTrWCEjQ7A1F46Gdm3GOa3q8dxXm8P6ou4ib0+yqMvN9q+tjlIuLfpKnaLZSek0jIvlog66Mqw/iAhPDe1M7nE3z3252eo4p+x9d39aygH+YvvV6ijl0qKv1CnIPHacpVsyGdY99C7UhZ1SM1vbN4rjlvNbMSspjdW7w6/LLGVvDqtMR26yL8YWohf29adVqVMwb20GHq9hxNktKt5ZVco/Lm5HszrVeeiT9RS6w2tJiw9+3k01Chnx24X9ELwPgBZ9pSrJGMPspHS6t6xD2wTr10ePNDVjHTx9VWdSM3OZtGy77wdaXGCz81zMS87gKvuPxEueJRl84dPSykqp3yWnHSE1M5eJw3TVcL8rKdgXAYPPWsZbS1MZ1KUJ7axN5ZPJy3dQ6PZyi/NLq6OclBZ9pSppdlLxIlqDzgyNm2KEhAC0rp/YPJjvPS/w8Kvv8HGMWNdHXvq9jSt/8tjhPBfTftrF5V2a0G5raN9AR4u+UpVQ4PKw4Ne9XN6lSUgtohWJGshRHnd+yANFdzHN05/RjkV/3ulEHzY+FGp/mvrjTvJcHu7t1w62Bvx0VaJ9+kpVwlcp+zhW6NYLuEEyzLacfrY1THRfR6q3qdVxypVjavLej7sY2LkxHRrHWR2nQlr0laqET9dk0Lxudc5pVc/qKFFBBCY636EGhdxf9DeKjN1/L+6nC79T3QM4Vujmnn7hcOVBu3eU8tm+nAJ+SD3IPf3aYbOJ1XHCWyUKbYLk8IxzMncV3c8b7qGMCWCsysoxNZnqGUD/To3o1LT2yXcOkWGb2tJXykfz1u7FGBgWIje4jiYD7asYZlvOW56hrNlz2Oo4//Oyezh5VOe+S9pbHcVn2tJXygfGGD5Zk06P0+qS+NsiWiHScosW45zvsaKwI/dMX8uCe86nbs2YU3gR//2bbfS25APPpdxg/4ZOTa/02+sGmrb0lfLB+owcUjNzQ+4uSNGkthQwKeZVso4Vcv/sZEuXYDbGMLZoFHXI5QHHx5blOBXa0lfKB5+uySDGoWPzrXambSePX96Jx+dt4D/LUrnbl4MC8BfZvOQMVpmOPOd4O6Rn35ZHW/pKVcDl9vJZcgb9OzUivrqOzbfaDee0ZEjXpry8eCvLPZ2Dfv6jx4t4ZuFmzpLU39fYCSNa9JWqwLItmRzOL+Jq7doJCSLCM1d1oV1CHH8ruo/N3uDNmTDG8Mgn68nOc/GU892QXUnzZLToK1WBuWszaFArlr7tGlgdRZWoGevg3dE9qcFxRrkeZJ8JzryJD1fs4Yv1+3igf3vOtO0Myjn9TYu+UieRk1/Et5syGXxWU103P8Q0rVOdd2OeJ5fqjHY9yFFTPaDnS/GexlMLNvKX9g2584I2AT1XIOlPsVIn8cX6fbg8Xq7SsfkhqZNtD5Ocr5JqmjLK9RA5JjD3pD1mqnN30b3UreHk5WvOCuvJeTp6R6mTmLc2g7aSTud3WkL4/p5HtL72DbzJG9xTdA/Xuh7jg5hnaShH/fb6uaYao1wPkmYa8uHIbtSvFeu317aCtvSVOoG07HxW7srmKvuPiBb8kDbAvoopzhfYbRpxjWss6cY/119yC93c7HqIZNOWN5xv0Lt1fb+8rpW06Ct1AvPWFq+LPsT+o8VJlC8usK/ng5iJHDS1ubJwAos93av0ekePFzFq6kqSTVvedL7OQPsqPyW1lhZ9pcphjGFucga9WtWjuRy0Oo7yUQ/bVubGjKWJHOK2on/yeNEojpvKz634KfUgA19dztq0I7zufLPigh+C98I9ES36SpVjXXoOO7Ly9AJuGGpr28vcmCe4zb6ADzz9ubjwRd53X+pT8c8pKGLc/BSun7yCWIeNj+88l0H2FUFIHTx6IVepcsxdm0GM3cblXZrAQqvTqMqKFTePOqdzkS2Zl9wjeMI9mtfdVzHC/h1n27bR1ZZKAzlKoXGQQ03WbzrAp2szWLzxAC63l9HnJfLgZR2pHuPH9ftDhBZ9pcpwe7wsWLeXi09P0GUXwlwf+0bOtT3JCtORt9xD+a/nSrye4g6OWFwUUrJS57Qk6tWM4fpeLRl+dnM6Nwv9bppTpUVfqTJ+SD3IwVwXQ7pq104kEIHespneMRPJN7FsMIkke9uSZeKJlzziyaP5jZM4r00DYhyR3+OtRV+pMj5L3kvtag4u6tjQ6ijKz2pIIb1kC71sW/74RIcEawJZIPI/1pSqhHyXm0Up+xl0ZhNiHZHXn6uUT0VfRAaIyBYRSRWRh8t5/gIRWSMibhEZXua5m0VkW8nXzf4KrlQgLN54gHyXR7t2ok3pIZdhMOyyKirs3hERO/AWcCmQDqwSkfnGmI2ldtsDjAL+WebYesBYoAdggNUlx4bOTS6VKmXu2gyaxlejV2JwVm1UYSYCPhB8aen3AlKNMTuMMS5gJjCk9A7GmF3GmHWAt8yxlwGLjTHZJYV+MTDAD7mV8ruDuYUs33aQwV2bhfWCWkqdjC9FvxmQVur79JJtvvDpWBG5XUSSRCQpKyvLx5dWyr++mHgDHq/hqp+vtjqKUgHjy+id8po8vt4uxqdjjTFvA28D9OjRI/xuRaMiwlzPeXSU3XSwpVsdRZ1MBHSxWMmXln46UPp+ZM2BvT6+flWOVSpodh3MI9m0Y6gurqYinC9FfxXQTkRaiUgMMBKY7+PrLwL6i0hdEakL9C/ZplRI+Sx5L4KXwfafrI6iVEBVWPSNMW7gboqL9SZgtjEmRUTGi8hgABHpKSLpwAjgvyKSUnJsNvAUxR8cq4DxJduUChnGGD5LzqCXbKap6I+nimw+zcg1xiykzLJTxpgnSj1eRXHXTXnHTgWmViGjUgG1PiOHHQfzuN3xg9VRlAo4XYahskpfRBqXY10O5Te/rag50L7S6ihKBZwW/dK0oEcdt8fL57/u46KODYnfnv/7EzpCREUoXXtHRbWfth/iYG4hQ3XZBRUltOiX4TY2Dpg6ZB47jjE6ZSDSzUvOIC7WwUUdo2eVRRXdor57x+3xsnRLFjNX7uHX45M4RBwGG0z4lga1Yji9SW16JdbjunNa0qBWrNVxlR8VuDws2rCfK85sSjWnrqipokPUFn2P1zD1h51M+WEn+48eJyEulovta2hENglyBNeAl9i07ygpe4/y8jdbeWNpKld3b85t3sa0tu33/UQn6hvWawaWW7zpAHkuD0P1PrgqikRl0d97pID7ZyWzYmc257dtwJNDzqBfxwScT5Vac+X8Vv97uD0rl8nLd/LJmnQ+dj/PGMfH3GFfgLYNw9u8tRk0ia/GOa10RU0VPaKu6C9K2c+Dc9bh9nh5ccRZXN29GSInX1GxTcNaPDusC2Mubc+4Z5/mefd1LPN05aXsfFrUq/HnA3TkR8g7lFvId1uzuK1va11RU0WVqLqQ+/mve7nrw9Uk1q/BF/f2ZfjZzSss+KU1jIvlTefrvOScxEZzGpc/v4CVT/TSIh+GFqzbh8drGNqtqdVRlAqqqGnpf52yn/tnJdMjsR7TRveiesypdc6IwNX25fSULYwu+hc3uh5hkvNV+vk5rwqseckZdGwcR8fGta2OokJRBDfkoqKl/93WLO6evpbOzeKZOqrnKRf80lraMpkdM572ks7tRWP4LDnDD0lVwJS6Fd6uJ9qzds8Rrjr4X6tTKRV0Ed/S356Vy10frqZtQi2mje5FrViH3z7F68sxpsdM4K+uB7hvlgOHzcagM5v45bVV4MzznqcraqqoFdEt/UK3h3tnrCXGYWPKqB7E13D6/RxxUsC0mOc4u2Vd7p+dzMqdukpjKDMG5nnOp7dtE010RU0VhSK66D//1RZS9h7lheFn0SS+esDOU02KeGf/CJp70vnrf79hm1fHfYeqNaYdu0xjrrLpipoqOkVs0V+6OZMpP+zk5nNP49JOjQJ+vrqSyzTnRGJxMcr1IFlGLxCGorme84nFpStqqqgVeX364+LJMTX5Z+GLdGzcnEcuP/1/2wOthe0g78a8wHDXWO4tuocPnc9gF12/J1S4jJ0Fnt5cZltFnBRYHUcpS0RkS/8V99UcJo6Xsu+m2oR6QR1+1dm2i6cc7/Kz9wxecQ8P2nlVxZZ6u3GEOK6ya9eOil4RV/S3eJvzgedSrrd/yxm23ZZkGOH4nmvtS3nTcxVLPF0tyaD+bK7nfBpwhL629VZHUcoyEVX0jTE86b6JWhTwgONjS7M86XiP02UX9xf9jQxT39IsCnJMTZZ4uzHY/hMO8VodRynLRFTRX5Syn5+8nRnjmENdya3ai5WazHMqqkkRk5yv4cbOP4vuxGt8XO6hiudV5VvgOQcXToZp146KchFT9I8XeXj6i010kD38P/s3VscBINF2gMcdH/Cz9wze81xmdZyoNtfTl3aSzhmyy+ooSlkqYor+oTwX9WvFMtbxfkj9+X6tfRn9bGt4zj2SVK8u7mWFnQfzSDIdGGZfTiXW11MqIkVM0W9Wpzrz/taHPvaNVkf5AxGY6JxMdQoZU3QXRaacdX+0SyegPl2Tjg0vw+zLrY6ilOUipugDlVomuVKqWJQT5AgTnFNZZ9rwtucKP4dTJ+P1Gj5dk0Ff2zoayRGr4yhluYgq+qFskH0FA20reM19FTu9ja2OEzV+2XGIjCMFXK2tfKWASJyRG8LGOafxQ2EXHnXfwkfGBO4vE/U/c1anE1fNQX+TVP4O2qWmooy29IOokRzhIccMfvJ25pM1uv5+oOUWuvlyw36uOLMp1aTI6jhKhQQt+kF2vX0JZ8sWnv5iI4dyC62OE9EWrt9HQZGH4Wc3tzqKUiFDi36Q2cQw0TmZvPwCnn32ce1eCKA5q9Np1aAm3VvWsTqKUiFDi74F2tkyuNW+kDmev7Da287qOBFp18E8Vu7MZvjZzfXaiVKlaNG3yD2OuTThEE8UjcLj6xINymezk9KwCdq1o1QZWvQtUlMKedT5ISmmFdM9F1sdJ6K4PV7mrE7nwg4JNKpdzeo4SoUULfoWGmRbQR/bBl5wX8shE2d1nIjx3dYsMo8Vck2PFlZHUSrkaNG3kAiMd7xHPrG86L7W6jgRY9aqNBrUiuHi0xOsjqJUyNGib7G2tr3cZP+amZ4LSfGeZnWcsJd1rJAlmzMZ1r05Trv+eCtVlv5WhIB/OOZSl1zGF92I0VvqVsnctem4vUa7dpQ6AS36ISBe8hjj+JgVphNfentZHSdsGWOYtSqNs0+rS9uEWlbHUSokadEPESPtS+koe3jGfT3HjdPqOGFp5c5stmflcW1PbeUrdSI+FX0RGSAiW0QkVUQeLuf5WBGZVfL8ChFJLNmeKCIFIpJc8vV//o0fORzi5QnH+6SbBKZ4Lrc6TliavnIPcdUcXHmm3qxGqROpsOiLiB14CxgIdAKuE5FOZXa7FThsjGkLvAI8V+q57caYriVfd/opd0TqY99If9sq/uMeTKbR5RkqIzvPxZfr93N19+ZUjynnRjVKKcC3ln4vINUYs8MY4wJmAkPK7DMEmFbyeA5wsejc91PyiGMGhTh5xT3C6ihhZc7qNFweL9ef09LqKEqFNF+KfjMgrdT36SXbyt3HGOMGcoD6Jc+1EpG1IvKdiPStYt6I18q2n5vsXzPLcyGb9h21Ok5Y8HoN01fsoWdiXdo30kluSp2ML0W/vBZ72YGFJ9pnH9DSGNMNGANMF5HafzqByO0ikiQiSVlZWT5Eimz3OuYSRz4TvtiE0TGcFfp5xyF2HcrXVr5SPvCl6Kcg4XF3AAAM+0lEQVQDpYdDNAf2nmgfEXEA8UC2MabQGHMIwBizGtgOtC97AmPM28aYHsaYHg0bNqz8u4gwdSSPfzg+5YfUgyzboh+CFfloxW7q1HAysHMTq6MoFfJ8KfqrgHYi0kpEYoCRwPwy+8wHbi55PBxYYowxItKw5EIwItIaaAfs8E/0yHaDfTGtGtRkwsJNuD1eq+OErANHj/N1ygGGd29ONadewFWqIhXeI9cY4xaRu4FFgB2YaoxJEZHxQJIxZj4wBfhARFKBbIo/GAAuAMaLiBvwAHcaY7ID8UYiTYx4eDjnKe4oGsOMsSO40fHN70+Oy7EuWIj5aMUePMZwQ29dwkIpX/h0Y3RjzEJgYZltT5R6fBz403ATY8wnwCdVzBi1+tuSOEc28op7OEPsP1JbCqyOFFIK3R6mr9jDRR0SSGxQ0+o4SoUFnZEbwkTgMedHZFOb/7jLjpJVC9fv42BuIaP6JFodRamwoUU/xHWx7WSY7XumegaS5m1gdZyQ8t5Pu2ndsCbnt9X/L0r5Sot+GPiXczY2vDznvs7qKCFj7Z7D/Jp2hFF9ErHZdB6gUr7yqU9fWauJZHO7fQGve65mtPcrzrY6UAiY9tMuasU6GNa91D1wx5VZukIveCv1J9rSDxN3OBaQwGGeKroBrze6J2xlHj3OF+v3MaJHc2rFartFqcrQoh8makoh/3LMItm04/N1ZefGRZepP+7C4zUVX8AdF//7l1IK0KIfVq62L6ez7OS5LzdT4PJYHccSx44X8dGK3Qzs3ITT6uswTaUqS4t+GLGJ4XHnB+zNOc7k5dE5sXnmyjSOHXdz+wWtizdoa16pStGiH2bOsW1mwBmNmfTddvbnHLc6TlAVebxM/XEnvVvX46wWdayOo1RY0qIfhh4ddDpur2Hil5usjhJUn/+6l305x7njgjZWR1EqbGnRD0Mt6tXg9r6tmZe8l6Rd0bGUkTGGt7/fQftGtbiwg67EqtSp0qIfpv52URsa167GuM9T8ETBEM4lmzPZvP8Yt/Vtjd6UTalTp0U/TNWIcfDI5R3ZkHGUj5PSKj4gjBljePWbbbSoV52h3cretE0pVRla9MPY4LOa0jOxLi8s2kJOfpHVcQJm6ZZM1mfkcPdFbXHa9UdWqarQ36AwJiKMG3wGh/NdvPD1ZqvjBMRvrfzmdav/vuSCDtNU6pRp0Q9zZzSN5+Y+iXy0Yg+/ph2xOo7fLduSxbr0HO7pp618pfxBf4siwJhL29OwViyPzdsQURd1i1v5W//YyldKVYkW/QgQV83Jo4NOZ31GDtNX7rE6jt98uymTX9Nz+Lv25SvlN/qbFCEGn9WUPm3q8/xXm8k8Gv4zdd0eL89+uYnWDWoy/Gxt5SvlL1r0I4SI8PTQzhS6vTzxWYrVcapsVlIa27PyeGhgR23lK+VH+tsUQVo3rMX9l7Tnq5T9fLl+n9VxTlluoZtXFm+jp2ym/+wOOkpHKT/Soh9hbuvbis7NavP4ZykcyXdZHeeUvP3ddg7mFvJv50fo5Ful/EuLfjgqPU69TCvYYbfx3NVncjjfxdNfhN+CbPtzjvPO8p1ccWYTutm2Wx1HqYijRT8CndE0njv/0po5q9P5dtMBq+NUylPPPY2nqIAHN19rdRSlIpLeYDQSlG7tl9wM/N6L27FkcxYPzlnHV/ddQMO4WIvC+W7plky+8J7LGMfHtLRl/vFJ7ddXyi+0pR+hYh12XhvZlWOFbh76ZB3GhPakrQKXh8fnbaCNZHCH/XOr4ygVsbToR7D2jeJ4ZGBHlmzO5KMVoT1p6/Ul20g/XMAE5xRixW11HKUilhb9CHfzuYn0bdeAp7/YyJb9x6yOU67N+4/yzvc7GHF2c3rbInPhOKVChRb9CGezCS+NOIu4okPc+dosjo5tbHWkPzhe5OEfr31EHe8RHtkwyOo4SkU8LfpRIKF2Nf4T8zpppiEPFN2FN4QWZXt24Sa2mJa86JxEPQnNv0SUiiRa9CPNCcbv97Rt4d+O6Sz29mDSd6Ex/v3bTQeY9vNubrUv5EL7OqvjKBUVdMhmFBlt/4pfvW148Wtol1CL/mdY19WTefQ4/5qzjk5NavNg9kzLcigVbbSlH0VE4FnnZM5sXod7ZqwlaVe2JTkKXB5u+2A1BS4Pr1/XTUfrKBVEWvSjTA0p5N1RPWlWpzq3vLeKrQeC24/u8Rrum7WWdelHeG1kV9om1Arq+ZWKdlr0o1C9mjFMu6UX1Zx2bp66krTs/KCd+9mFm1iUcoDHB3WytHtJqWilRT9KtahXg/dG9yLf5eHqST8FZQz/299vZ/IPOxnVJ5Fbzm8V8PMppf5Mi34U69S0Nh977kOO7eOaVxeyevfhgJzHGMMLizbzzMLNDOrShMev6BSQ8yilKqajdyKZD4uUtbdlMCdmHDcW/ZsbJq/ghRFncsWZTf0WweM1PDZvAzNW7uG6Xi15emhn7DbRBdSUsoi29BUtbAf5OGYcHRrHcff0tTw451fyXVUfUbM/5zij3l3JjJV7uPuitjxzVUnBV0pZxqeWvogMAF4D7MBkY8zEMs/HAu8DZwOHgGuNMbtKnnsEuBXwAPcaYxb5Lb06NeW0shvKUT7OvJJX7MOZlHQlSbsP8/TQzvRp06DSL2+MYe7aDMbNT6HoeC7POD7k+p+XwM/+CK+UqooKi76I2IG3gEuBdGCViMw3xmwstdutwGFjTFsRGQk8B1wrIp2AkcAZQFPgGxFpb4zx+PuNqKpziocHnbM437aef2bdyfXv5HGuLYUxt/+Vnon1Kjze6zUs3ZLJ5OU7+XnHIXqcVpcX991Doi28buSiVCTzpaXfC0g1xuwAEJGZwBCgdNEfAowreTwHeFNEpGT7TGNMIbBTRFJLXk/bfCGsj30jS2wPMMPTj7fcQxjxfz/TLqEWF3ZoyIUdEmjTsBY1Y+3UjHFw4NhxUjKOsmFvDvPWZrDrUD6Na1dj7JWduOncROzjteArFUp8KfrNgLRS36cD55xoH2OMW0RygPol238pc2yzU06rgqaaFDHasYiR9qXM9lzI4oNnMy2zI+8s31nu/iLQrUUdHujfgQGdG+O06+UipUKRL0W/vCtvZZdpPNE+vhyLiNwO3F7yba6IbPEh14k0AA5W4fhwFOD3PKfk6+R2AXMDF6I0/TeODtH3np+Uqrzn03zZyZeinw60KPV9c2DvCfZJFxEHEA9k+3gsxpi3gbd9CVwREUkyxvTwx2uFi2h7z9H2fkHfc7QIxnv25W/wVUA7EWklIjEUX5idX2af+cDNJY+HA0tM8U1Z5wMjRSRWRFoB7YCV/omulFKqsips6Zf00d8NLKJ4yOZUY0yKiIwHkowx84EpwAclF2qzKf5goGS/2RRf9HUDf9eRO0opZR2fxukbYxYCC8tse6LU4+PAiBMcOwGYUIWMleWXbqIwE23vOdreL+h7jhYBf89S3AujlFIqGui4OqWUiiIRU/RFZICIbBGRVBF52Oo8gSYiLURkqYhsEpEUEfmH1ZmCRUTsIrJWRBZYnSUYRKSOiMwRkc0l/97nWp0p0ETk/pKf6w0iMkNEqlmdyd9EZKqIZIrIhlLb6onIYhHZVvLfuv4+b0QU/VJLRQwEOgHXlSwBEcncwAPGmNOB3sDfo+A9/+YfwCarQwTRa8BXxpiOwFlE+HsXkWbAvUAPY0xnigeQjLQ2VUC8Bwwos+1h4FtjTDvg25Lv/Soiij6lloowxriA35aKiFjGmH3GmDUlj49RXAgifraziDQHBgGTrc4SDCJSG7iA4hFyGGNcxpgj1qYKCgdQvWTeTw3Kmd8T7owx31M82rG0IcC0ksfTgKH+Pm+kFP3yloqI+AL4GxFJBLoBK6xNEhSvAg8CXquDBElrIAt4t6RLa7KI1LQ6VCAZYzKAF4E9wD4gxxjztbWpgqaRMWYfFDfsgAR/nyBSir5Pyz1EIhGpBXwC3GeMOWp1nkASkSuATGPMaquzBJED6A5MMsZ0A/IIwJ/8oaSkH3sI0Iri1XlrisgN1qaKHJFS9H1a7iHSiIiT4oL/kTHmU6vzBMF5wGAR2UVxF14/EfnQ2kgBlw6kG2N++ytuDsUfApHsEmCnMSbLGFMEfAr0sThTsBwQkSYAJf/N9PcJIqXo+7JUREQpWbp6CrDJGPOy1XmCwRjziDGmuTEmkeJ/4yXGmIhuARpj9gNpItKhZNPF/HFZ80i0B+gtIjVKfs4vJsIvXpdSekmbm4HP/H2CiLhH7omWirA4VqCdB9wIrBeR5JJt/y6ZPa0iyz3ARyUNmh3AaIvzBJQxZoWIzAHWUDxKbS0RODtXRGYAFwINRCQdGAtMBGaLyK0Uf/iVu9JBlc6rM3KVUip6REr3jlJKKR9o0VdKqSiiRV8ppaKIFn2llIoiWvSVUiqKaNFXSqkookVfKaWiiBZ9pZSKIv8fTF3LTiwchcoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1be1cf41ba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#看看效果\n",
    "plt.plot(x,func(x))\n",
    "plt.hist(points,normed=True,bins=100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 六.小结\n",
    "这一节的主要内容是将细致平衡方程成功落地到MH算法中，但是对于实际问题中建议分布$q(x\\rightarrow x')$的采样往往比较麻烦，因为实际数据中，随机变量$x$往往都是高维数据，比如对于文本数据，用单词做随机变量的话，能高达上万维，在那么高维的空间造一个容易采样的建议分布极其困难，那我们借鉴一下优化中常使用的坐标轮换法（比如SVM中所使用的SMO）思路，采样有没有可能只需在单个维度上进行呢？然后遍历随机变量的所有维度完成一次采样，可以的，那就是下一节将要介绍的单分量MH算法...."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}