00_maximum_likelihood_estimationa_and_3d_visualization/maximum_likelihood_estimation_and_3d_visualization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><strong>A Journey Into Math For Machine Learning</strong></center>\n",
    "<center><strong>机器学习之数学之旅</strong></center>\n",
    "<center><strong>图解极大似然估计$(maximum \\ likelihood \\ estimation \\ with \\ 3D \\ visualization)$</strong></center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<script>requirejs.config({paths: { 'plotly': ['https://cdn.plot.ly/plotly-latest.min']},});if(!window.Plotly) {{require(['plotly'],function(plotly) {window.Plotly=plotly;});}}</script>"
      ],
      "text/vnd.plotly.v1+html": [
       "<script>requirejs.config({paths: { 'plotly': ['https://cdn.plot.ly/plotly-latest.min']},});if(!window.Plotly) {{require(['plotly'],function(plotly) {window.Plotly=plotly;});}}</script>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "import plotly.graph_objs as go\n",
    "from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\n",
    "from IPython.display import Image \n",
    "init_notebook_mode(connected=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_likelihood(observations, mu, sigma):\n",
    "    # 定义y轴取值\n",
    "    plt.ylim(-0.02,1)\n",
    "    # 定义一个画图范围\n",
    "    x_locs =  np.linspace(-10, 10, 500)\n",
    "    # 画出推断的概率分布的概率密度函数\n",
    "    plt.plot(x_locs, stats.norm.pdf(x_locs, loc=mu, scale=sigma), label=\"inference\")\n",
    "    for obs in observations:\n",
    "        plt.axvline(x=obs, ymin=0, ymax=stats.norm.pdf(obs, loc=mu, scale=sigma)+0.01, c=\"g\")\n",
    "    plt.axvline(x=obs, ymin=0, ymax=stats.norm.pdf(obs, loc=mu, scale=sigma)+0.01, c=\"g\", label=\"probabilities\")\n",
    "    # 画出观测数据的概率\n",
    "    plt.scatter(x=observations, y=[0 for _ in range(len(observations))], c=\"r\", marker=\"o\", label=\"obsevations\")\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    plt.title(\"mean={} sigma={}\".format(str(mu), str(sigma)))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**极大似然估计**是神经网络和很多复杂模型得以优化求解的理论基础, 我们今天来学习并试着深入理解极大似然估计的原理和推导, 最后我们对极大似然估计进行3D可视化, 建立一种直观的认识.\n",
    "   \n",
    "要理解**极大似然估计**是什么, 首先要明白**概率密度(质量)函数**是什么, 如果你不知道的话, 那就简短解释一下：   \n",
    "概率密度函数用来描述某个随机变量取某个值的时候，取值点所对应的的**概率**$(probability)$的函数.   \n",
    "如下图, 我们现在有一个**概率分布**, 属于**正态分布**: $$X \\sim N(\\mu,\\sigma^2), \\quad  f(x;\\mu,\\sigma)=\\frac{1}{\\sigma\\sqrt{2\\pi}} \\, \\exp \\left( -\\frac{(x- \\mu)^2}{2\\sigma^2} \\right)$$\n",
    "\n",
    "其中$\\mu$是均值, $\\sigma$是标准差. 如果你不熟悉正态分布, 我们简单回顾一下$\\mu$指的是均值, 在下图中, 均值是$0$则正态分布的概率在均值处概率最高, 以均值为中心两边是对称的,  $\\sigma$是标准差, 标准差控制着概率分布偏离均值的程度, 标准差越大概率分布越扁平, 越小的话, 概率分布越集中于均值.   \n",
    "   \n",
    "我们另有一个**数据点**, 是一个**随机变量**, 取值$2.5$, 我们将$x=2.5$代入$f(x;\\mu=5,\\sigma=2)$得出下图出中绿色直线的长度, 也就是得到了$P(x=2.5 \\mid \\mu=5, \\sigma=2)$    \n",
    "意义为$x=2.5$在上面定义的正态分布中的概率, 也就是给定一个概率分布, 随机变量在这个概率分布中出现的可能性, 而$f(x;\\mu,\\sigma)=\\frac{1}{\\sigma\\sqrt{2\\pi}} \\, \\exp \\left( -\\frac{(x- \\mu)^2}{2\\sigma^2} \\right)$就是概率密度函数.   \n",
    "概率质量函数是离散的, 概率密度函数是连续的, 意义相同, 为了可视化的方便, 今天用概率密度函数来讲解极大似然估计."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW9//HXJwmEQDBhE4EgQWWHsIQtIBpUKlIVrbivqEW06PW2LmDdt59evdZqsbR1r9QIWhDRXlCEiig7yCoS1oQqsiSQECAk+f7+mEkcQpbJOsnx/Xw8hpzlO2c+c2Z4z5nvmfmOOecQERFvCQt1ASIiUv0U7iIiHqRwFxHxIIW7iIgHKdxFRDxI4S4i4kEKd5FKMLP1ZpYc6jpESqNwl58V83nWzPb5L8+amVV0O865Hs65BTVQYrUws3vNbJ2ZZZnZNjO7N9Q1Se2KCHUBIrVsHHAJ0BtwwKfANmBKKIuqAQbcAKwBTgfmmlmacy4ltGVJbdGRu1SJmW33HyWuMbNDZvaambU2s3/5jxo/M7NmAe0Hm9lXZpZpZt8Edm2Y2Vgz2+i/3lYzuy1gXbKZpZvZ78zsRzP73szGVqLkG4H/dc6lO+d2Af8L3FTKfWtpZrP9te43s4VmFhZwv8/zT0eZ2VtmluGv/z4zS6/CPppuZj+Y2QEz+8LMelT0Tjrn/sc5t9I5l+ec2wR8CAyt6Hak/lK4S3W4DBgBdAYuAv4FPAC0wvccuwvAzNoBHwNPAs2Be4APzKyVfzs/AhcCJwFjgT+YWb+A2zkFiAHaAbcAkwtD0cwm+kO4xEvANnoA3wTMf+NfVpLfAen++9Haf59KGq/jESAeOM2/H66r7D7y+xfQCTgZWAlMLVxRgftJwHUMGAasL+V+ihc553TRpdIXYDtwbcD8B8CfA+bvBGb6p+8H/l7s+nOAG0vZ9kzgv/zTycBhICJg/Y/A4ArWmw90DZjvhC+wrYS2j+M74j2jlPt9nn96K3B+wLpbgfTK7KMSbifWX19MFR6jx/C9iEWG+vmiS+1ddOQu1WF3wPThEuaj/dMdgMuLHWmeCbQBMLMLzGyxvwskExgFtAzY1j7nXF7AfE7AtoOVje+dQaGTgGznT8FingNS8fVXbzWziaVssy2QFjCfVkKboPaRmYWb2TNmtsXMDuJ7YYDj90PQzGwCvr73XzrnjlZmG1I/KdylNqXhO3KPDbg0cc49Y2aR+I5onwdaO+digU/wnRgsl5k9YGbZpV0Cmq7HdzK1UG9K6a5wzmU5537nnDsNuBj4rZmdW0LT74G4gPn2wdRcimuA0cB5+Lqg4v3LDSp0PzGzm4GJwLnOuXTkZ0XhLrXpHeAiMzvff4TayH+iNA5oCEQCe4A8M7sA+EWwG3bOPe2ciy7tEtD0bXwh3c7M2uLrV3+zpG2a2YVmdoa/z/oAvi6dghKaTgMmmVkz/3mFCcHWXYKmwFFgH9AYeLoy99PMrvVfd4RzbmsV6pF6SuEutcY5l4bvqPQBfCGeBtwLhDnnsvCdVJwGZOA7gp1VA2X8BfgIWAusw3eC9y+ltO0EfIavK+dr4BXn3PwS2j2O78TrNn/79/EFdGW8DewAdgEbgMWV3M6TQAtgWcCRvdc+7illsJK7GkWksszsduAq59zZoa5Ffr505C5SRWbWxsyGmlmYmXXB19UzI9R1yc9bueFuZq/7vzSyrpT1ZmYvmVmq/0sa/UpqJ+JhDfF17WQBn+P7+OQrIa1IfvbK7ZYxs7Pw9Tm+7ZzrWcL6Ufg+pzsKGAT80Tk3qAZqFRGRIJV75O6c+wLYX0aT0fiC3znnFgOxZtamugoUEZGKq46Bw9px/Jc20v3Lvi/e0MzG4Ru4iaioqMT27Sv3ceCCggLCwure6QLVVTGqq+Lqam2qq2KqUtd333231znXqtyGwXyNFd8XKdaVsm42cGbA/Dygf3nbTExMdJU1f/78Sl+3JqmuilFdFVdXa1NdFVOVuoDlrpaGH9jF8d/Ii/MvExGREKmOcJ8F3OD/1Mxg4IBz7oQuGRERqT3l9rmb2bv4RuRr6R+j+hGgAYBzbgq+8T9G4RtgKQffUK0iIhJC5Ya7c+7qctY74DfVVpGIhNyxY8dIT0/nyJEjQbWPiYlh48aNNVxVxdXnuho1akRcXBwNGjSo1G3oZ/ZE5ATp6ek0bdqU+Ph4LIifmM3KyqJp06a1UFnF1Ne6nHPs27eP9PR0OnbsWKnbqHufERKRkDty5AgtWrQIKtil+pkZLVq0CPqdU0kU7iJSIgV7aFV1/yvcRUQ8SOEuInXSkCFDym2zcOFCevToQZ8+fTh8+HAtVFV/KNxFpE766quvym0zdepUJk2axOrVq4mKiiq3vXOOgoKSfkzLexTuIlInRUf7fjVwwYIFJCcnM2bMGLp27cq1116Lc45XX32VadOm8dBDD3HttdcC8NxzzzFgwAASEhJ45JFHANi+fTtdunThhhtuoGfPnqSlpTF37lySkpLo168fl19+OdnZvp+fjY+P55FHHqFfv3706tWLb7/9FoDs7GzGjh1Lr169SEhI4IMPPgAodTt1gT4KKSJleuyj9Wz4z8Ey2+Tn5xMeHh70Nru3PYlHLuoRdPtVq1axfv162rZty9ChQ1m0aBG33norX375JRdeeCFjxoxh7ty5bN68maVLl+Kc4+KLL2bRokV07dqVzZs389ZbbzF48GD27t3Lk08+yWeffUaTJk149tlneeGFF3j44YcBaNmyJStXruSVV17h+eef59VXX+WJJ54gJiaGtWvXApCRkVHudkJN4S4idd7AgQOJi4sDoE+fPmzfvp0zzzzzuDZz585l7ty59O3bF/AdbW/ZsoWuXbvSoUMHBg8eDMDixYvZsGEDQ4cOBSA3N5ekpKSi7fzqV78CIDExkX/+858AfPbZZ6SkpBS1adasGbNnzy5zO6GmcBeRMgVzhF3TXxaKjIwsmg4PDycvL++ENs45Jk2axG233XZcXfv27aNJkybHtRsxYgTvvvtumbdV2u0Eu51QU5+7iHjC+eefz+uvv17U771r1y727NlzQrvBgwezaNEiUlNTATh06BDfffddmdseMWIEkydPLprPyMio1HZqk8JdRDzhF7/4Bddccw1JSUn06tWLMWPGkJWVdUK7Vq1a8eabb3L11VeTkJBAUlJS0YnT0jz44INkZGTQs2dPevfuzfz58yu1nVoVzKDvNXHRj3XUHtVVMXW1Ludqr7YNGzZUqP3BgwdrqJKqqe91lfQ4UIs/1iEiInWMwl1ExIMU7iIiHqRwFxHxIIW7iIgHKdxFRDxI4S4inpScnMzKlSuDbv/mm28yYcKEEtcVDj+8fft2evbsCcDy5cu56667AN/gZoGjWE6ZMoW33367sqVXCw0/ICL1VkUHLKuskoYf7t+/P/379wd84R4dHV30IjB+/Pgar6k8OnIXkTpp+/btRUP8duvWjTFjxpCTk0N8fDz3338//fr1Y/r06axevZrBgweTkJDApZdeSkZGRtE2UlJS6NOnDz179mTp0qUALF26lKSkJPr27cuQIUPYtGlTUfu0tDSSk5Pp1KkTjz32WNHywuGHAy1YsIALL7yQ7du3M2XKFP7whz/Qp08fFi5cyKOPPsrzzz8PwJYtWxg5ciSJiYkMGzas6Fus06dPL/rG61lnnVXt+09H7iJSprv/725W/7C6zDYVPYLuc0ofXhz5YrntNm3axGuvvcbQoUO5+eabeeWVVwBo0aJFUZdLQkICL7/8MmeffTYPP/wwjz32GC++6Nv24cOHWb16NV988QU333wz69ato2vXrixcuJCIiAg+++wzHnjggaLx2ZcuXcq6deto3LgxAwYM4Je//GXR0Xlp4uPjGT9+PNHR0dxzzz0AzJs3r2j9uHHjmDJlCp06dWLJkiXccccdfPjhhzz++OPMmTOHdu3akZmZGfS+C5bCXUTqrPbt2xcNqXvdddfx0ksvAXDllVcCcODAATIzMzn77LMBuPHGG7n88suLrj9mzBgAzjrrLA4ePEhmZiZZWVnceOONbN68GTPj2LFjRe1HjBhBixYtAN/Qv19++WW54V6W7Oxsvvrqq+NqOnr0KABDhw7lpptu4oorrigaZrg6KdxFpEzBHGHX1JC/ZlbifOAQvhW9/kMPPcTw4cOZMWMG27dvJzk5udzbq6yCggJiY2NZvfr4dz5ZWVlMmTKFJUuW8PHHH5OYmMiKFSuKXliqg/rcRaTO2rlzJ19//TUA//jHP074gY6YmBiaNWvGwoULAfj73/9edBQPFP3YxpdffklMTAwxMTEcOHCAdu3aAb5PyAT69NNP2b9/P4cPH2bmzJlF7xrK07Rp0xJHoDzppJPo2LEj06dPB3wDNX7zzTeAry9+0KBBPP7447Rq1Yq0tLSgbitYCncRqbO6dOnC5MmT6datGxkZGdx+++0ntHnrrbe49957SUhIYPXq1cf9zF1kZCR9+/Zl/PjxvPbaawDcd999TJo0ib59+57wYxwDBw7ksssuIyEhgcsuuyzoLpmLLrqIGTNmFJ1QDTR16lRee+01evfuTY8ePfjwww8BuPfee+nVqxc9e/ZkyJAh9O7du0L7plzBDB1ZExcN+Vt7VFfF1NW6nPt5Dfm7bds216NHjyptQ0P+ioiIpyjcRaROio+PZ926daEuo95SuIuIeJDCXUTEgxTuIiIepHAXEfGgoMLdzEaa2SYzSzWziSWsP9XM5pvZKjNbY2ajqr9UEfm5CxxytzbVxSF9y1Pu8ANmFg5MBkYA6cAyM5vlnNsQ0OxBYJpz7s9m1h34BIivgXpFRGpdXRzStzzBHLkPBFKdc1udc7lACjC6WBsHnOSfjgH+U30likhdFzFtGsTHQ1iY7+/UqdWy3RdeeIGePXvSs2fPopEe8/LyThgGGGDixIl0796dhISEotEZ9+7dy2WXXcaAAQMYMGAAixYtoqCggPj4+ONGYuzUqRO7d+/mo48+YtCgQfTt25fzzjuP3bt3lzukb2lDDicnJ3P//fczcOBAOnfuXPTN1fXr15OcnEyfPn1ISEhg8+bN1bKvijPfF57KaGA2BhjpnLvVP389MMg5NyGgTRtgLtAMaAKc55xbUcK2xgHjAFq3bp2YkpJSqaKzs7NLHF851FRXxaiuiqut2mJiYjjjjDOCahsxbRqN7rwTO3y4aJmLiuLIyy+Td8UVla5h1apV3H777cybNw/nHOeccw5/+9vfGDZsGHPnzmXw4MHccccdRWO+jxgxghUrVmBmZGZmEhsby9ixYxk3bhxJSUmkpaVx6aWXsnz5cu677z4SEhK47rrrWLZsGU888QSzZs0iIyOD2NhYzIy33nqLTZs28fTTT/P0008THR1d9MtLgfNJSUk899xznHnmmTz55JNkZWXx7LPPMmrUKPr06cPTTz/NnDlzmDx5MrNmzeKee+4hMTGRq6++mtzcXPLz84mKiipxH6SmpnLgwIHjlg0fPnyFc678cRHK+worMAZ4NWD+euBPxdr8FvidfzoJ2ACElbVdDT9Qe1RXxdTVupyro8MPdOjgHJx46dChSjW8+OKL7qGHHiqaf/DBB90f//hH1759+6Jl8+bNc6NHj3bHjh1zCQkJbuzYse6DDz5wR48edc4517JlS9e7d++iS9u2bV1WVpZbtGiRO//8851zzt19993ur3/9q3POuTVr1rgRI0a4nj17us6dOxe1eeSRR9xzzz1XdLuF85mZmcfVk5qa6vr27eucc+7ss892X375pXPOuR9++MGdfvrpzjnnpk6d6rp27eqeeeYZ991335W5D2p6+IFdQPuA+Tj/skC3ANP8LxZfA42AlkFsW0Tqu507K7a8ikoaljciIoKlS5cyZswYZs+ezciRIwHfkLuLFy9m9erVrF69ml27dhEdHU1SUhKpqans2bOHmTNnFo2nfueddzJhwgTWrl3LX/7yF44cOVKlWiMjIwEIDw8vGqTsmmuuISUlhaioKEaNGsXnn39epdsoTTDhvgzoZGYdzawhcBUwq1ibncC5AGbWDV+476nOQkWkjjr11IotD9KwYcOYOXMmOTk5HDp0iBkzZjBs2LAShwHOzs7mwIEDjBo1ij/84Q9Fw+qec845vPzyy0XbLBxX3cy49NJL+e1vf0u3bt2KxlEPHA74rbfeKrpeaUP6ljfkcEm2bt1Kx44dueuuuxg9ejRr1qyp7C4qU7nh7pzLAyYAc4CN+D4Vs97MHjezi/3Nfgf82sy+Ad4FbvK/fRARr3vqKVzxPuPGjeGpp6q02X79+nHTTTcxcOBABg0axK233kqzZs1KHAY4KyuLCy+8kISEBM4880xeeOEFAJ577jmWL19OQkIC3bt3Z8qUKUXbv/LKK3nnnXeKftUJ4NFHH+Xyyy8nMTGRli1/6nwoa0jfsoYcLsm0adMYNGgQffr0Yd26ddxwww1V2k+lCqbvpiYu6nOvPaqrYupqXc7V0T5351zOq6/6+tjNfH/feadG6qqon/OQv/qZPRGpsrwrroBbbgl1GRJAww+IiHiQwl1ESuR02iykqrr/Fe4icoJGjRqxb98+BXyIOOfYt28fjRo1qvQ21OcuIieIi4sjPT2dPXuC+0TzkSNHqhRENaU+19WoUSPi4uIqfRsKdxE5QYMGDejYsWPQ7RcsWEDfvn1rsKLK+TnXpW4ZEREPUriLiHiQwl1ExIMU7iIiHqRwFxHxIIW7iIgHKdxFRDxI4S4i4kEKdxERD1K4i4h4kMJdRMSDFO4iIh6kcBcR8SCFu4iIByncRUQ8SOEuIuJBCncREQ9SuIuIeJDCXUTEgxTuIiIepHAXEfEghbuIiAcp3EVEPEjhLiLiQQp3EREPUriLiHiQwl1ExIOCCnczG2lmm8ws1cwmltLmCjPbYGbrzewf1VumiIhURER5DcwsHJgMjADSgWVmNss5tyGgTSdgEjDUOZdhZifXVMEiIlK+YI7cBwKpzrmtzrlcIAUYXazNr4HJzrkMAOfcj9VbpoiIVEQw4d4OSAuYT/cvC9QZ6Gxmi8xssZmNrK4CRUSk4sw5V3YDszHASOfcrf7564FBzrkJAW1mA8eAK4A44Augl3Mus9i2xgHjAFq3bp2YkpJSqaKzs7OJjo6u1HVrkuqqGNVVcXW1NtVVMVWpa/jw4Succ/3LbeicK/MCJAFzAuYnAZOKtZkCjA2YnwcMKGu7iYmJrrLmz59f6evWJNVVMaqr4upqbaqrYqpSF7DclZPbzrmgumWWAZ3MrKOZNQSuAmYVazMTSAYws5b4umm2BrFtERGpAeWGu3MuD5gAzAE2AtOcc+vN7HEzu9jfbA6wz8w2APOBe51z+2qqaBERKVu5H4UEcM59AnxSbNnDAdMO+K3/IiIiIaZvqIqIeJDCXUTEgxTuIiIepHAXEfEghbuIiAcp3EVEPEjhLiLiQQp3EREPUriLiHiQwl1ExIMU7iIiHqRwFxHxIIW7iIgHKdxFRDxI4S4i4kEKdxERD1K4i4h4kMJdRMSDFO4iIh6kcBcR8SCFu4iIByncRUQ8SOEuIuJBCncREQ9SuIuIeJDCXUTEgxTuIiIepHAXEfEghbuIiAcp3EVEPEjhLiLiQQp3EREPUriLiHiQwl1ExIOCCnczG2lmm8ws1cwmltHuMjNzZta/+koUEZGKKjfczSwcmAxcAHQHrjaz7iW0awr8F7CkuosUEZGKCebIfSCQ6pzb6pzLBVKA0SW0ewJ4FjhSjfWJiEglmHOu7AZmY4CRzrlb/fPXA4OccxMC2vQDfu+cu8zMFgD3OOeWl7CtccA4gNatWyempKRUqujs7Gyio6Mrdd2apLoqRnVVXF2tTXVVTFXqGj58+ArnXPld3865Mi/AGODVgPnrgT8FzIcBC4B4//wCoH95201MTHSVNX/+/EpftyapropRXRVXV2tTXRVTlbqA5a6cfHXOBdUtswtoHzAf519WqCnQE1hgZtuBwcAsnVQVEQmdYMJ9GdDJzDqaWUPgKmBW4Urn3AHnXEvnXLxzLh5YDFzsSuiWERGR2lFuuDvn8oAJwBxgIzDNObfezB43s4trukAREam4iGAaOec+AT4ptuzhUtomV70sERGpCn1DVUTEgxTuIiIepHAXEfEghbuIiAcp3EVEPEjhLiLiQQp3EREPUriLiHiQwl1ExIMU7iIiHqRwFxHxIIW7iIgHKdxFRDxI4S4i4kEKdxERD1K4i4h4kMJdRMSDFO4iIh6kcBcR8SCFu4iIByncRUQ8SOEuIuJBCncREQ9SuIuIeJDCXUTEgxTuIiIepHAXEfEghbuIiAcp3EVEPEjhLiLiQQp3EREPUriLiHiQwl1ExIMU7iIiHhRUuJvZSDPbZGapZjaxhPW/NbMNZrbGzOaZWYfqL1VERIJVbribWTgwGbgA6A5cbWbdizVbBfR3ziUA7wP/U92FiohI8CKCaDMQSHXObQUwsxRgNLChsIFzbn5A+8XAddVZpEhNyC9wbN2TzTfpB9i6J5sfs46yL/so+/Yf4Z0dy4iJakhcsyhObd6YvqfG0rFlE8ws1GWLBMWcc2U3MBsDjHTO3eqfvx4Y5JybUEr7PwE/OOeeLGHdOGAcQOvWrRNTUlIqVXR2djbR0dGVum5NUl0VE4q6juY71uzJZ+XuPL7Zk09Onm95mEFMQ+OkSKOgIB8snEPHHBlHHIX/Q5o2hJ4twhnYJoKeLcNpEFb7Qa/HsmK8WNfw4cNXOOf6l9cumCP3oJnZdUB/4OyS1jvn/gr8FaB///4uOTm5UrezYMECKnvdmqS6KqY269q+9xB/X7yD91ekc+DwMZo1bsCo3m0Y1LE5vdvHcnqraML9YR1YV25eATv2HWLFjgyWbtvPvG9/5OvvjxIT1YCrBrbnhqR42sVG1cp9KF5bXaK6KqY26gom3HcB7QPm4/zLjmNm5wG/B852zh2tnvJEqmbb3kO8NG8zH67eRZgZI3uewjUDT2Vgx+ZEhJf/eYKGEWF0at2UTq2bctXAU8nNK2DRlr1MW5bG377YyqsLt3FJn3b894hOxDVrXAv3SCQ4wYT7MqCTmXXEF+pXAdcENjCzvsBf8HXf/FjtVYpU0IGcY/zvp5t4Z/EOGkaE8ethp3HLsI6c3LRRlbbbMCKM4V1OZniXk0nbn8Mbi7bzzpIdfPTNf7hucAfuHtGJkxo1qKZ7IVJ55Ya7cy7PzCYAc4Bw4HXn3HozexxY7pybBTwHRAPT/SecdjrnLq7BukVKVFDgeH9FOs/837dk5uRy3eAO3HlOJ1o1jaz222rfvDEPX9SdW4d15KV5m3njq23MXvMfHr6oO7/s1UYnXyWkgupzd859AnxSbNnDAdPnVXNdIhWWnpHDvdPX8PXWffTv0IzHRg+kR9uYGr/dtrFRPHNZAtcMOpUHZqxlwj9WMb1zOs9c1os2MbXXHy8SSN9QlXrPOce05WmMfHEha9IzeeZXvZg+PqlWgj1QQlwsH/7mTB65qDtLt+1n5IsLmb3mP7Vag0ihav20jEhtO3jkGBM/WMMna39gUMfmPH95b9o3D92JzfAwY+zQjiR3OZn/fm81E/6xis83/sgTl/SkSaT+u0nt0bNN6q31/znAb6auJC3jMBMv6Mq4YacRFoLPnpekY8smvD8+iZc/T+XlzzezZtcBplzXjzNObhrq0uRnQt0yUi+lLN3Jpa98xZFjBbw3bjDjzz69zgR7oYjwMP57RGfeuWUQGYdyufhPi5j1jbpppHYo3KVeOZZfwIMz1zLxn2sZ1LE5H991Jv3jm4e6rDINOaMlH981jG5tTuKud1fx2EfrycsvCHVZ4nEKd6k3MnNyufH1pbyzeCe3nX0ab44dSIvo6v+IY004JaYRKeMGM3ZoPG8s2s7YN5dxIOdYqMsSD1O4S72Q+mMWl0xexPLtGTx/eW8mXdCtaLiA+qJBeBiPXNSDZy/rxeKt+7j0lUVs3ZMd6rLEoxTuUuct2PQjl07+iuyjebw7bhBjEuNCXVKVXDngVKbeOpjMw8e4ZPIiFm7eE+qSxIMU7lJnOed4/ctt3PzmMuKaN2bmb4aS2KFu968Ha2DH5nz4m6G0jY3ipjeW8eaibZQ3QqtIRSjcpU7KzStg0j/X8vjsDZzXrTXvj0/y3MBc7Zs35v3bhzC8y8k8+tEGfj9zHcd0olWqicJd6pz9h3K5/rUlpCxL4zfDT2fKdYme/QJQdGQEf70+kTuST+cfS3Zyw2tLyTiUG+qyxAMU7lKnfLfbd+J0VVomL17Zh3vP71rnPr9e3cLCjPtGduWFK3qzYkcGl7yyiNQfdaJVqkbhLnXG59/u5levfEVObj4p4wZzSd92oS6pVv2qXxzvjhvMoaN5XPrKIv79nU60SuUp3CXknHP87Yut3PLWcjq0aMysCUPpd2qzUJcVEokdmjHzN0NpFxvF2DeW8oZOtEolKdwlpI7m5XPf+2t46pONjOxxCtPHJ9G2Fn+2ri6Ka9aYD24fwrndWvPYRxt4YIZOtErFefMsldQLB486rnt1Ccu2Z3DXuZ24+9xOnu9fD1aTyAj+cl0iz83dxJ8XbGHb3mz+fG1iqMuSekThLiHxTVomj319mOy8I7x8dV8u6t021CXVOWFhxv0ju9Lp5GgmfrCWS15ZxLiu6qKR4KhbRmqVc45/LNnJ5VO+BuD98UMU7OX46URrPo8tPszMVSf8Pr3ICXTkLrXmyLF8Hpq5jukr0jmrcyuuiDtEr7ja/bWk+iqxQzOONv8dP+zN5u73XmLFjgwevLAbkRHhoS5N6igduUut2LInm8v+/BXTV6Rz17mdeOOmAUQ3VP96RaQd3EFB2G7GnXUaf1+8gyumfE16Rk6oy5I6SuEuNco5R8rSnVz40pf8J/Mwr9/Un9+O6FzvRnSsSx4Y1Y0p1yWydc8hLvjjQmau2qWPS8oJ1C0jNSYzJ5eJH6zl/9b/wNAzWvDCFX1ofVKjUJflCSN7nkKPtifx3++t5u73VvPZxt08dUkvYho3CHVpUkco3KVG/N+6H3jow3Vk5uQy6YKu/LoO/b6pV7Rv3pj3bktiyr+38IdPvysa6/7MTi1DXZo0er6oAAANoUlEQVTUAeqWkWq1J+sod0xdwfh3VtAqOpIZdwzltjr4+6ZeER5m/Gb4Gcy4YyhNIsO57rUl/G7aN+zX4GM/ezpyl2qRX+BIWbaT5+ZsIic3n3vP78K4s06jQbiOH2pDr7gYPr5rGC9/vpm//Hsrn3+7mwdGdWNMYhxmemH9OVK4S5V9vWUfj8/ewMbvDzKwY3OevrQXZ5wcHeqyfnYaNQjn3vO7MrpPOyb9cy33vr+Gd5fu5Pe/7OaZHzmR4CncpdI2/ZDFC59uYs763bSLjeKVa/txQc9TdKQYYp1bN2X6bUlMX5HG83O/47I/f80FPU/h/pFdiW/ZJNTlSS1RuEuFfbc7iz/O28wna7+nScMIfjeiM78+6zQaNdAXauqKsDDjygGncmFCW/62cCt//WIrn27YzSV923F78umc3krvrLxO4S5Bcc6xZNt+3ly0nTkbfqBxg3DuSD6dW888jWZNGoa6PClFk8gI7j6vM9cMOpVX5m8hZdlOPliZzqiebbjt7NNIiIsNdYlSQxTuUqZDR/OYveY/vLFoO9/+kEVs4wYK9Xro5KaNePTiHkw45wzeWLSNt7/awcdrv6d3XAzXDDqVi3q3pXFDxYGX6NGUE+TlF/Bl6l5mrtrFnPW7OXwsn66nNOWZX/VidJ92RDVU90t91TI6knvP78ptZ5/OjJW7mLpkB/d/sJYnZ29kVK82XNi7DUmntSBCn3Kq9xTuAsDBI8dY+N1e5n27mwWb9rD/UC4xUQ24tF87Lu3bjv4dmulEqYec1KgBNw6J54akDizfkcG7S3by8drveW95Gs2bNOT8HqcwvEsrkk5vQdNG+tZrfaRw/5k6kHOMFTv3s2x7Bsu27Wd1WiZ5BY7Yxg04u3MrRvVqQ3KXVhp10OPMjAHxzRkQ35wjx/L593d7mL3mez5cvYt3l+4kIszod2ozhpzRgn6nNqN3+1hiohT29UFQ4W5mI4E/AuHAq865Z4qtjwTeBhKBfcCVzrnt1VuqVEZBgWPf4QI+/3Y3G7/P4tsfstj4/UG27MnGOYgIM3rFxfDrs07j3K4n06d9rN6S/0w1ahDO+T1O4fwep5CbV8CKHRks3LyHLzbv4Y/zNlM4NtnprZrQu30sXVo3pXPrppxxcjQFGriszik33M0sHJgMjADSgWVmNss5tyGg2S1AhnPuDDO7CngWuLImCpafHM3LZ/+hXPZl57I3+yh7/X/TM3JI23+YtP05pGccJje/AP69HIB2sVF0a9OU0b3b0j++OX3ax6oPvT6YOhUOHoTISIiPh6eegmuvrbGbaxgRRtLpLUg6vQX3jexK1pFjrE0/wKq0TFbtzGDh5r38c+VPPxrSMBziV/+btrFRtIuNol0z398WTSJp1qQBzZs0pFnjhvq4bC0K5sh9IJDqnNsKYGYpwGggMNxHA4/6p98H/mRm5mpgHNIfDx5hS2Y+J+3M4Ket+yYK513gtH/iuGX81LBwE8XXuePW/bQNSmhX+Hftj3kc27C7xNssvGKBg2P5BeTmFZBX4Iqmj+U78vILfPP5vuVH8/I5dDSfQ0fzOJSbR3bh9NE8so/kkXU0r8R9FBPVgFObN6Zrm6aM6NGaw3vSueisRLqc0pST1H9a/0ydCuPGwV3+J9OOHb55qNGAD9S0UQOGnNGSIWf8NChZZk4um3/MZvPubOav3AhNmrAr4zCr0zLJzDlW4naaNAwntnFDmjaKoHHDcBo39P1tEhlBVMNwmjQMJzIinIhwo0F4GA3CjYiwMBpEhNEgzIjwL2sQHkaY+bqVDAgzwwz/5adl6/bmE7F5L2EGmL8dvjZh/vbBCf58UzDbPHC05t/pBBPu7YC0gPl0YFBpbZxzeWZ2AGgB7K2OIgPNWLWL/7f4CCz+qro3XT1WLq/yJhr6n8ANI8JoEhlBdGQETSIjiIlqQLvYRjRp6Jtv0aQhLaIjaRl9/N/oyOMf1gULdjMgXl8/r7d+/3vIKfajHDk5vuW1FO4liW3csKi/vu3hrSQn9y9al300j+8zD7P/UC4ZObnsO5RLxqFc9h86RkZOLtlH8zicm8+h3Dz2Zh/lUK5vPvtoHrl5BRRUZ/YtX1KNG6seN3RvyOgavo1aPaFqZuOAcQCtW7dmwYIFFd5Gs5wCbu/uiIryjQte+CL506ulFf1b/BU0cLZwXbnLik1YKcsADh8+TOOoqGL1BNbom4owiAjzXcLDjPDCeaPoaOR4BUCu/1LMYd8law9kAdtObEF2dnal9nVNU11BuvNOAPKyH8SFhbHg+ed/WldH6ixrnzXCd/TXLgxo6r+cIMx/8b2zLHCOvALId5BfAHnOke+fL1zunDvu3XEBFL0bL3znnZNzmEb+/5MFP71h968P7hWkIv0PwTaNtSM1/xxzzpV5AZKAOQHzk4BJxdrMAZL80xH4jtitrO0mJia6ypo/f36lr1uTVFfFqK4gdejgHLiYibgmjzZyzp9NrkOHUFdWpM7tMz8v1gUsd+XktnMuqPHclwGdzKyjmTUErgJmFWszC7jRPz0G+NxfhIhU1VNPQePGxy9r3Ni3XKQU5XbLOF8f+gR8R+fhwOvOufVm9ji+V5BZwGvA380sFdiP7wVARKpDYb/6+ut9fzt0qPFPy0j9Z6E6wDazPcCOSl69JTVwsrYaqK6KUV0VV1drU10VU5W6OjjnWpXXKGThXhVmttw517/8lrVLdVWM6qq4ulqb6qqY2qhLX0UUEfEghbuIiAfV13D/a6gLKIXqqhjVVXF1tTbVVTE1Xle97HMXEZGy1dcjdxERKYPCXUTEg+psuJvZ5Wa23swKzKx/sXWTzCzVzDaZ2fmlXL+jmS3xt3vP/+3a6q7xPTNb7b9sN7PVpbTbbmZr/e2qPrJY+XU9ama7AmobVUq7kf59mGpmE2uhrufM7FszW2NmM8ysxF9nrq39Vd79N7NI/2Oc6n8uxddULQG32d7M5pvZBv/z/79KaJNsZgcCHt+Ha7ou/+2W+biYz0v+/bXGzPrVQk1dAvbDajM7aGZ3F2tTa/vLzF43sx/NbF3AsuZm9qmZbfb/bVbKdW/0t9lsZjeW1KZCghmjIBQXoBvQBVgA9A9Y3h34BogEOgJbgPASrj8NuMo/PQW4vYbr/V/g4VLWbQda1uK+exS4p5w24f59dxrQ0L9Pu9dwXb8AIvzTzwLPhmp/BXP/gTuAKf7pq4D3auGxawP08083Bb4roa5kYHZtPZ+CfVyAUcC/8I2VNxhYUsv1hQM/4PuST0j2F3AW0A9YF7Dsf4CJ/umJJT3vgebAVv/fZv7pZlWppc4euTvnNjrnNpWwajSQ4pw76pzbBqTiG3O+iPmGVTwH39jyAG8Bl9RUrf7buwJ4t6ZuowYUjdPvnMsFCsfprzHOubnOucJB6BcDcTV5e+UI5v6PxvfcAd9z6Vw7ccjOauWc+945t9I/nQVsxDeoYn0wGnjb+SwGYs2sTS3e/rnAFudcZb/5XmXOuS/wDcESKPB5VFoWnQ986pzb75zLAD4FRlalljob7mUoaXz54k/+FkBmQJCU1KY6DQN2O+c2l7LeAXPNbIV/2OPaMMH/1vj1Ut4GBrMfa9LN+I7ySlIb+yuY+3/c7xQAhb9TUCv83UB9gZIGJE8ys2/M7F9m1qOWSirvcQn1c+oqSj/ACsX+KtTaOfe9f/oHoHUJbap934X0B7LN7DPglBJW/d4592Ft11OSIGu8mrKP2s90zu0ys5OBT83sW/8rfI3UBfwZeALff8Yn8HUZ3VyV26uOugr3l5n9HsgDppaymWrfX/WNmUUDHwB3O+cOFlu9El/XQ7b/fMpMoFMtlFVnHxf/ObWL8Q1JXlyo9tcJnHPOzGrl8+chDXfn3HmVuNouoH3AfJx/WaB9+N4SRviPuEpqUy01mlkE8Ct8Pw5e2jZ2+f/+aGYz8HUJVOk/RbD7zsz+BswuYVUw+7Ha6zKzm4ALgXOdv7OxhG1U+/4qQTD3v7BNuv9xjsH33KpRZtYAX7BPdc79s/j6wLB3zn1iZq+YWUvnXI0OkBXE41Ijz6kgXQCsdM7tLr4iVPsrwG4za+Oc+97fTfVjCW124Ts3UCgO3/nGSquP3TKzgKv8n2ToiO8VeGlgA39ozMc3tjz4xpqvqXcC5wHfOufSS1ppZk3MrGnhNL6TiutKaltdivVzXlrK7QUzTn911zUSuA+42DmXU0qb2tpfdfJ3Cvx9+q8BG51zL5TS5pTCvn8zG4jv/3GNvugE+bjMAm7wf2pmMHAgoDuippX67jkU+6uYwOdRaVk0B/iFmTXzd6P+wr+s8mrjDHJlLvhCKR04Cuzm+F+D+j2+TzpsAi4IWP4J0NY/fRq+0E8FpgORNVTnm8D4YsvaAp8E1PGN/7IeX/dETe+7vwNrgTX+J1ab4nX550fh+zTGllqqKxVfv+Jq/2VK8bpqc3+VdP+Bx/G9+IDvF+Km++teCpxWC/voTHzdaWsC9tMoYHzh8wyY4N833+A7MT2kFuoq8XEpVpcBk/37cy0Bn3Kr4dqa4AvrmIBlIdlf+F5gvgeO+fPrFnznaeYBm4HPgOb+tv2BVwOue7P/uZYKjK1qLRp+QETEg+pjt4yIiJRD4S4i4kEKdxERD1K4i4h4kMJdRMSDFO4iIh6kcBcR8aD/D8eGUrIYLmm/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw_likelihood([2.5], mu=0, sigma=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**概率**和**似然**的区别:   \n",
    "**概率**，是在已知一些概率分布参数的情况下，预测观测的结果；   \n",
    "**似然**，则是用于在已知某些观测所得到的结果时，对观测结果所属于的的概率分布的参数进行估值。\n",
    "   \n",
    "**极大似然估计**的目的在于找到一个最符合当前观测数据的概率分布.   \n",
    "1. 我们先理解似然函数是什么\n",
    "例如下面两图中: **红色圆点**指的是观测到的随机变量, **蓝色的线**是概率密度函数的图像, **绿色的直线**的**长度**所对应的$y$轴的值指的是是观测的数据出现在当前概率分布中的可能性, 也就是概率, 概率是介于$[0，1]$之间的实数, 我们用$P(x \\mid \\mu,\\sigma)$来表示绿色线的长度, 也就是概率.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4VNX9x/H3Nws7RAQEIUBQEYEQtrAEVGILghalCLjhAmrRUrT+bG2l7rhUK3UtiChWrFSFUhFxQUQimwiIyBJkEaIEZZE1EQIkOb8/ZpKGkGQmIWQml8/reebJ3Dvnnvudw+STy507Z8w5h4iIeEtEqAsQEZHyp3AXEfEghbuIiAcp3EVEPEjhLiLiQQp3EREPUriLlIGZrTWz5FDXIVIchbucUszsITM7amaZBW5nlbYf51xb51zKSSixXJjZ3Wa2xswyzGyLmd0d6pqkYkWFugCREHjbOXddqIs4yQy4AVgFnA18bGZbnXNvhbYsqSg6cpcTYmZp/qPEVWb2s5lNMrOGZvah/6jxEzOrW6B9dzNbbGb7zOzrgqc2zGy4ma3zb7fZzG4t8FiymaWb2R/MbKeZ/Whmw0/yc6tvZrP8te4xswVmFlHgeff2369uZpPNbK+//j+ZWfoJjNE0M9tuZvvNbL6ZtS1t7c65vznnVjjnsp1z64F3gZ4nPipSWSjcpTwMAvoA5wKXAR8CfwEa4HuN3QFgZk2A94FHgdOBPwLTzayBv5+dQH+gDjAceMbMOhXYTyMgBmgC3AyMywtFM7vHH8JF3grVe5k/rNea2W9LeF5/ANL9z6Oh/zkVNV/Hg0AccJZ/HIr6X0FQY+T3IdASOANYAUzJe6CUzzNvGwMuANaW8FzFa5xzuulW5huQBgwtsDwdeLHA8u3ADP/9PwP/KrT9bODGYvqeAfzefz8ZOAREFXh8J9C9lPW2ARoDkUAP4EfgmmLajsF3xHtOMc+7t//+ZqBvgcduAdLLMkZF7Oc0fH9QYk7g3+hh4GugaqhfL7pV3E1H7lIedhS4f6iI5Vr++82BIYWONM8HzgQws0vMbIn/qHofcClQv0Bfu51z2QWWDxboOyjOuVTn3A/OuRzn3GLgOWBwMc2fAjbhO1+92czuKaZdY2BrgeWtRbQJaozMLNLMnjCzb83sAL4/DHDsOATNzEbhO/f+K+fc4bL0IZWTwl0q0lZ8R+6nFbjVdM49YWZV8R3RjgUaOudOAz7A98ZgQGb2l0JXwBxzK2FTV9w+nHMZzrk/OOfOAi4H7jKzXxbR9EcgtsBy02BqLsa1wACgN75TUHH+9Qale55mdhNwD/BL51w6ckpRuEtFegPf+e6+/iPUav43SmOBKkBVYBeQbWaXABcH27Fz7nHnXK3ibnntzGyAmdU1n674znW/W1SfZtbfzM7xn7PeD+QAuUU0nQqM9vfbBBgVbN1FqA0cBnYDNYDHy/g8h/q37eOc23wC9UglpXCXCuOc24rvqPQv+EJ8K3A3EOGcy8AXtFOBvfiOYGeehDKuxneqJQN4HXjSOTe5mLYtgU+ATOBzYLxzbl4R7cbge+N1i7/9f/AFdFm8DnwHbANSgSVl7OdRoB6wrMCR/YQy9iWVkDmnL+sQKU/+K3Cuds71CnUtcurSkbvICTKzM82sp5lFmFkrfJdQvhPquuTUFjDczexV/4dG1hTzuJnZ82a2yf8hjU5FtRPxsCrAS/hO9XyK7xz++JBWJKe8gKdlzOxCfOccX3fOxRfx+KX4rtO9FOgGPOec63YSahURkSAFPHJ3zs0H9pTQZAC+4HfOuSXAaWZ2ZnkVKCIipVceE4c14dgPbaT71/1YuKGZjQBGAFSvXr1z06Zluxw4NzeXiAi9XRCIxikwjVFwNE6BVdQYbdiw4SfnXINA7Sp0Vkjn3ERgIkBiYqJbvnx5mfpJSUkhOTm5HCvzJo1TYBqj4GicAquoMTKz74JpVx5/ZrZx7CfyYv3rREQkRMoj3GcCN/ivmukO7HfOHXdKRkREKk7A0zJm9ia+Gfnq++eofhCIBnDOTcA3/8el+D71dxDfVK0iIhJCAcPdOXdNgMcd8Ltyq0hEQu7o0aOkp6eTlZUFQExMDOvWrQtxVeGtvMeoWrVqxMbGEh0dXabt9TV7InKc9PR0ateuTVxcHGZGRkYGtWvXDnVZYa08x8g5x+7du0lPT6dFixZl6kPXNonIcbKysqhXrx6+CTGlopkZ9erVy/+fU1ko3EWkSAr20DrR8Ve4i4h4kMJdRMJSjx49ArZZsGABbdu2pUOHDhw6dKgCqqo8FO4iEpYWL14csM2UKVMYPXo0K1eupHr16gHbO+fIzS3qy7S8R+EuImGpVi3ftwbmfax/8ODBnHfeeQwdOhTnHK+88gpTp07l/vvvZ+jQoQA89dRTdOnShYSEBB588EEA0tLSaNWqFTfccAPx8fFs3bqVjz/+mKSkJDp16sSQIUPIzPR9/WxcXBwPPvggnTp1ol27dnzzzTcAZGZmMnz4cNq1a0dCQgLTp08HOKafG264Ib+fcKBLIUWkRA+/t5bVW/cSGRlZbn22aVyHBy9rG3T7r776irVr19K4cWN69uzJokWLuOWWW1i4cCH9+/dn8ODBfPzxx2zcuJGlS5finOPyyy9n/vz5NGvWjI0bNzJ58mS6d+/OTz/9xKOPPsonn3xCzZo1efLJJ3n66ad54IEHAKhfvz4rVqxg/PjxjB07lldeeYVHHnmEmJgYVq9eDcDevXuP62fMmDHH9BNqCncRCXtdu3YlNjYWgA4dOpCWlsb5559/TJuPP/6Yjz/+mI4dOwK+o+2NGzfSrFkzmjdvTvfu3QFYsmQJqamp9OzZE4AjR46QlJSU388VV1wBQOfOnfnvf/8LwCeffMJbb72V36Zu3brMmjXrmH6ysrLy74cDhbuIlOjBy9qG/ENMVatWzb8fGRlJdnb2cW2cc4wePZpbb731mPVpaWnUrFnzmHZ9+vThzTffLHFfxe2nuH5CPUaF6Zy7iHhC3759efXVV/PPe2/bto2dO3ce16579+4sWrSITZs2AfDzzz+zYcOGEvvu06cP48aNy1/eu3dvmfqpSAp3EfGEiy++mGuvvZakpCTatWvH4MGDycjIOK5dgwYNeO2117jmmmtISEggKSkp/43T4tx3333s3buX+Ph42rdvz7x5847rp3fv3gH7qUgBv0P1ZNGXdZx8GqfANEZFW7duHa1bt85fDrdTDuHoZIxR4X8HADP70jmXGGhbHbmLiHiQwl1ExIMU7iIiHqRwFxHxIIW7iIgHKdxFRDxI4S4inpScnExpLrd+7bXXGDVqVJGP5U0/nJaWRnx8PADLly/njjvuAHyX1H7xxRf57SdMmMDrr79e1tLLhaYfEJFKKycnp1wnNCtOUdMPJyYmkpjou9w8JSWF6OhoevfuDcBtt9120msKREfuIhKW0tLS8qf4bd26NYMHD+bgwYPExcXx5z//mU6dOjFt2jRWrlxJ9+7dSUhIYODAgezduze/j3/961906NCB+Ph4li5dCsDSpUtJSkqiY8eO9OjRg/Xr1+e337p1K8nJybRs2ZKHH344f33e9MMFpaSk0L9/f9LS0pgwYQLjxo2jQ4cOLFiwgIceeoixY8cC8O2339KvXz86d+7MBRdckP8p1mnTpuV/4vXCCy8s9/HTkbuIlOjOj+7ky21flusRcodGHXi237MB261fv55JkybRs2dPbrrpJsaPHw9AvXr1WLFiBQAJCQm88MIL9OrViwceeICHH36YZ5/19X3w4EFWrlzJ/Pnzuemmm1izZg3nnXceCxYsICoqik8++YS//OUv+fOzL126lDVr1lCjRg26dOnCr371q/yj8+LExcVx2223ER0dzb333gvA3Llz8x8fMWIEEyZMoGXLlnzxxReMHDmSTz/9lDFjxjB79myaNGnCvn37Sj+IASjcRSRsNW3aNH8a3euuu47nn38egKuuugqA/fv3s2/fPnr16gXAjTfeyJAhQ/K3v+aaawC48MILOXDgAPv27SMjI4Mbb7yRjRs3YmYcPXo0v32fPn2oV68e4Jv6d+HChQHDvSSZmZksXrz4mJoOHz4MQM+ePRk2bBhXXnll/jTD5UnhLiIlerbfsyGbW8bMilwuOIVvabe///77ueiii3jnnXdIS0s7Zm6h4vZXVrm5uZx22mmsXLnyuMcmTJjAF198wfvvv0/nzp358ssv8/+wlAedcxeRsPX999/z+eefA/Dvf//7uC/oiImJoW7duixYsADwnWPPO4oHePvttwFYuHAhMTExxMTEsH//fpo0aQL4rpApaM6cOezZs4dDhw4xY8aMoL98o3bt2kXOQFmnTh1atGjBtGnTAN8c8F9//TXgOxffrVs3xowZQ4MGDdi6dWtQ+wqWwl1EwlarVq0YN24crVu3Zu/evfz2t789rs3kyZO5++67SUhIYOXKlcd8zV21atXo2LEjt912G5MmTQLgT3/6E6NHj6Zjx47HfRlH165dGTRoEAkJCQwaNCjoUzKXXXYZs2bNyn9DtaApU6YwadIk2rdvT9u2bXn33XcBuPvuu2nXrh3x8fH06NGD9u3bl2psAtGUvx6mcQpMY1S0cJjyNy0tjf79+7NmzZoK3W9ZacpfERE56RTuIhKW4uLiKs1RezhSuIuIeJDCXUTEgxTuIiIepHAXEfGgoMLdzPqZ2Xoz22Rm9xTxeDMzm2dmX5nZKjO7tPxLFZFTXcEpdytSSkrKMTNDhsOUvoEEnH7AzCKBcUAfIB1YZmYznXOpBZrdB0x1zr1oZm2AD4C4k1CviEiFS0lJoVatWvnzuofDlL6BBHPk3hXY5Jzb7Jw7ArwFDCjUxgF1/PdjgB/Kr0QRCXtTpkBcHERE+H5OmVIu3T799NPEx8cTHx+fP9Njdnb2cdMAA9xzzz20adOGhIQE/vjHPwKwa9cuBg0aRJcuXejSpQuLFi0iNzeXuLi4Y2ZibNmyJTt27OC9996jW7dudOzYkd69e7Njx478KX2feeaZIqf0zZtyOCkp6Zgph5OTk/nzn/9M165dOffcc/M/ubp27Vq6du1Khw4dSEhIYOPGjeUyVoUFM3FYE6DgpAfpQLdCbR4CPjaz24GaQO+iOjKzEcAIgIYNG5KSklLKcn0yMzPLvO2pROMUmMaoaDExMcfMlZKTk1Pk3CkAUVOnUu3227FDh3wrvvsO95vfkJWVRfaVV5a5hq+++opJkyYxd+5cnHP84he/IDExkfXr1/PCCy8wYcIERo4cyTPPPMPQoUOZPn06X375JWaWP/vjyJEjufXWW0lKSmLr1q0MHDiQ5cuXc8kll/Dmm29y3XXXsWzZMmJjY6lRowbt27dnzpw5mBmTJ0/m0Ucf5fHHH2f48OHUqlUr/5uXPvjgA6Kjo8nIyOC6667jqaeeIikpib/+9a/ce++9PPnkk+Tk5HDw4EHmzp3L7NmzeeCBB5g5cybPP/88I0aM4KqrruLIkSMljm1WVlbZX5/OuRJvwGDglQLL1wP/KNTmLuAP/vtJQCoQUVK/nTt3dmU1b968Mm97KtE4BaYxKlpqauoxywcOHCi+cfPmzsHxt+bNT6iGZ5991t1///35y/fdd5977rnnXNOmTfPXzZ071w0YMMAdPXrUJSQkuOHDh7vp06e7w4cPO+eca9CggWvfvn3+rXHjxi4jI8MtWrTI9e3b1znn3J133ukmTpzonHNu1apVrk+fPi4+Pt6de+65+W0efPBB99RTT+XvN2953759+fUcOHDAbdq0yXXs2NE551yvXr3cwoULnXPObd++3Z199tnOOeemTJni2rRp45544gm3YcOGEseg8L+Dc84By12A3HbOBXVaZhvQtMByrH9dQTcDU/1/LD4HqgH1y/bnRkQqle+/L936E1TUtLxRUVEsXbqUwYMHM2vWLPr16wf4ptxdsmQJK1euZOXKlWzbto1atWqRlJTEpk2b2LVrFzNmzMifT/32229n1KhRrF69mpdeeomsrKwTqrVq1aoAREZG5k9Sdu211zJz5kyqV6/OpZdeyqeffnpC+yhOMOG+DGhpZi3MrApwNTCzUJvvgV8CmFlrfOG+qzwLFZEw1axZ6dYH6YILLmDGjBkcPHiQn3/+mXfeeYcLLrigyGmAMzMz2b9/P5deeinPPPNM/rS6F198MS+88EJ+n3nzqpsZAwcO5K677qJ169b586gXnA548uTJ+dsVN6VvoCmHi7J582bOOuss7rjjDgYMGMCqVavKOkQlChjuzrlsYBQwG1iH76qYtWY2xswu9zf7A/AbM/saeBMY5v/vg4h43WOPQY0ax66rUcO3/gR06tSJYcOG0bVrV7p168Ytt9xC3bp1i5wGOCMjg/79+5OQkMD555/P008/DcDzzz/P8uXLSUhIoE2bNkyYMCG//6uuuoo33ngj/1udAB566CGGDBlC586dqV//fycfLrvsMt55550ip/TNm3I4KSnpuCmHizJ16lTi4+Pp0KEDa9as4YYbbjihcSpWMOduTsZN59xPPo1TYBqjopXqnLtzzr3xhu8cu5nv5xtvnLTawlXAMSqDEznnrq/ZE5ETN3So7yZhQ9MPiIh4kMJdRIrk9LZZSJ3o+CvcReQ41apVY/fu3Qr4EHHOsXv3bqpVq1bmPnTOXUSOExsbS3p6Ort2+a5ozsrKOqGgORWU9xhVq1aN2NjYMm+vcBeR40RHR9OiRYv85ZSUFDp27BjCisJfuI2RTsuIiHiQwl1ExIMU7iIiHqRwFxHxIIW7iIgHKdxFRDxI4S4i4kEKdxERD1K4i4h4kMJdRMSDFO4iIh6kcBcR8SCFu4iIByncRUQ8SOEuIuJBCncREQ9SuIuIeJDCXUTEgxTuIiIepHAXEfEghbuIiAcp3EVEPEjhLiLiQQp3EREPUriLiHiQwl1ExIMU7iIiHhRUuJtZPzNbb2abzOyeYtpcaWapZrbWzP5dvmWKiEhpRAVqYGaRwDigD5AOLDOzmc651AJtWgKjgZ7Oub1mdsbJKlhERAIL5si9K7DJObfZOXcEeAsYUKjNb4Bxzrm9AM65neVbpoiIlEYw4d4E2FpgOd2/rqBzgXPNbJGZLTGzfuVVoIiIlF7A0zKl6KclkAzEAvPNrJ1zbl/BRmY2AhgB0LBhQ1JSUsq0s8zMzDJveyrROAWmMQqOximwcBujYMJ9G9C0wHKsf11B6cAXzrmjwBYz24Av7JcVbOScmwhMBEhMTHTJycllKjolJYWybnsq0TgFpjEKjsYpsHAbo2BOyywDWppZCzOrAlwNzCzUZga+o3bMrD6+0zSby7FOEREphYDh7pzLBkYBs4F1wFTn3FozG2Nml/ubzQZ2m1kqMA+42zm3+2QVLSIiJQvqnLtz7gPgg0LrHihw3wF3+W8iIhJi+oSqiIgHKdxFRDxI4S4i4kEKdxERD1K4i4h4kMJdRMSDFO4iIh6kcBcR8SCFu4iIByncRUQ8SOEuIuJBCncREQ9SuIuIeJDCXUTEgxTuIiIepHAXEfEghbuIiAcp3EVEPEjhLiLiQQp3EREPUriLiHiQwl1ExIMU7iIiHqRwFxHxIIW7iIgHKdxFRDxI4S4i4kEKdxERD1K4i4h4kMJdRMSDFO4iIh6kcBcR8SCFu4iIByncRUQ8KKhwN7N+ZrbezDaZ2T0ltBtkZs7MEsuvRBERKa2A4W5mkcA44BKgDXCNmbUpol1t4PfAF+VdpIiIlE4wR+5dgU3Ouc3OuSPAW8CAIto9AjwJZJVjfSIiUgZRQbRpAmwtsJwOdCvYwMw6AU2dc++b2d3FdWRmI4ARAA0bNiQlJaXUBQNkZmaWedtTicYpMI1RcDROgYXbGAUT7iUyswjgaWBYoLbOuYnARIDExESXnJxcpn2mpKRQ1m1PJRqnwDRGwdE4BRZuYxTMaZltQNMCy7H+dXlqA/FAipmlAd2BmXpTVUQkdIIJ92VASzNrYWZVgKuBmXkPOuf2O+fqO+finHNxwBLgcufc8pNSsYiIBBQw3J1z2cAoYDawDpjqnFtrZmPM7PKTXaCIiJReUOfcnXMfAB8UWvdAMW2TT7wsERE5EfqEqoiIByncRUQ8SOEuIuJBCncREQ9SuIuIeJDCXUTEgxTuIiIepHAXEfEghbuIiAcp3EVEPEjhLiLiQQp3EREPUriLiHiQwl1ExIMU7iIiHqRwFxHxIIW7iIgHKdxFRDxI4S4i4kEKdxERD1K4i4h4kMJdRMSDFO4iIh6kcBcR8SCFu4iIByncRUQ8SOEuIuJBCncREQ9SuIuIeJDCXUTEgxTuIiIepHAXEfEghbuIiAcp3EVEPCiocDezfma23sw2mdk9RTx+l5mlmtkqM5trZs3Lv1QREQlWwHA3s0hgHHAJ0Aa4xszaFGr2FZDonEsA/gP8rbwLFRGR4EUF0aYrsMk5txnAzN4CBgCpeQ2cc/MKtF8CXFeeRYpIxdiVcZjV2/bxzfYMvt99kP2HjpJ5OJvM/VnM2P4VjWKq0/KMWrQ+sw7nNapNRISFumQpRjDh3gTYWmA5HehWQvubgQ+LesDMRgAjABo2bEhKSkpwVRaSmZlZ5m1PJRqnwDRGsP3nXBb/kM3Xu3L47kBu/vo6VYxdkROJioBm3EL6+h/Zk+XIcXmPQ9t6kXQ7M4p29SOJPMWDPtxeS8GEe9DM7DogEehV1OPOuYnARIDExESXnJxcpv2kpKRQ1m1PJRqnwE7VMcrNdcxZt4N/LtrCks172Bs9kTNqV+Xuvo/TJe50WjWqTUz1aJJf851hfSiuFsnJyWTn5JK2+yBfb93H/I27mL9hF5//eJgzalflmq7NGNYjjro1q4T42YVGuL2Wggn3bUDTAsux/nXHMLPewL1AL+fc4fIpT0TKk3OO2Wt38OwnG/hmewaxdatzd99W/HvzXqIjI/jdReeUuH1UZATnnFGLc86oxaDOsRzJzmXe+p28vWwrz83dyMsLNnN9UnNGJp9DTPXoCnpWUpRgwn0Z0NLMWuAL9auBaws2MLOOwEtAP+fcznKvUkRO2MYdGTw4cy2Lv93NWQ1q8sxV7bksoTFRkRFM+65sV0VXiYqgb9tG9G3biPXbMxifsomJ8zfzn+Xp/LFvK65MbHrKn64JlYDh7pzLNrNRwGwgEnjVObfWzMYAy51zM4GngFrANDMD+N45d/lJrFtEgpR1NIen52zg1YVb2F/1ZTq1qssHN04q99Bt1ag2z13dkd9ccBYPv7eW0f9dzb+/+J6nhiRwXqM65bovCSyoc+7OuQ+ADwqte6DA/d7lXJeIlIPV6fv5v6kr2bQzk2u6NmXB3j0ctr0n9Wg6vkkMU29NYubXP/DIrFQuf2ER/9fnXEZceJaO4iuQPqEq4kG5uY4X5m5k4PhFZGZl86+bu/LXKxKIiqyYcDUzBnRowuw7L+SXrc/gyY++4cqXPueHfYcqZP+icBfxnD0/H2HYa8v4+5wNXNruTGbfeSEXtGwQklrq1arK+KGdeO7qDqzfnsGvnl/AZxt2haSWU43CXcRDVm7dx2UvLGTJt7t5fGA7nru6AzE1QnvVSt5R/MxRPWlYpxrD/rmUp+dsICfXhbQur1O4i3jE1GVbGTJhMWYw/bc9uLZbM/wXOISFsxrU4p2RPRnUKZbn527k5snLOJB1NNRleZbCXaSSy811/PWDdfxp+iq6n1WPWbefT7vYmFCXVaTqVSIZO6Q9jw2MZ+HGnxg0fjFb9xwMdVmepHAXqcQOHsnmtje+5KX5m7m+e3P+OawLp9UI/0+IDu3WnNdv6srOjMMMGLeIZWl7Ql2S5yjcRSqp7fuzGDLhcz5Zt4OHLmvDmAFtiYqsPL/SPc6pzzsje3Ba9WiufXkJ05ZvDbyRBK3yvBJEJN+abfsZMG4haT/9zKQbuzCsZ4uwOr8erLzz8F1bnM7d/1nFkx99Q67eaC0XCneRSmZO6g6GTPicqIgIpo/swUXnnRHqkk5ITI1oXhvelWu7NePFlG+5/c2vyDqaE+qyKr1ynRVSRE6uyYvTePi9tbRrEsPLNyZyRu1qoS6pXERHRvDYr+OJq1eDv374DT/sP8TLNyRSv1bVUJdWaenIXaQSyM11PDIrlds//D3769xF0xbTPBPsecyMEReezYtDO7HuxwMMHL+ITTszQl1WpaVwFwlzWUdzGDllBZMWbuH0034gx3awZtfXoS7rpOkXfyZvj0ji0JFcBo5fzKJNP4W6pEpJ4S4SxnZnHuaal5cwO3U7D/RvQ1y9mqEuqUK0b3oaM37Xg8Yx1bnx1aVMXaYraUpL4S4SpjbvyuSKFxeT+sMBXhzamZvObxHqkipUbN0aTPttEkln1+NP01fxN11JUyoKd5EwtDxtD1e8uJjMrGzeGtGdfvGNQl1SSNSpFs2rw7pwbbdmjNeVNKWiq2VEwsy05Vu59501xJ5endeGdaVZvRqhLimk8q6kaVGvJo9/uE5X0gRJR+4iYSI7J5dHZqVy939W0aVFXf772x6nfLDnMTN+c+FZvDi0M+t+PMCvxy1i4w5dSVMShbtIGNh/8CjDX1vGpIVbGNYjjsnDu1aKOWIqWr/4Rrw9Iomso7lc8aKupCmJwl0kxDbtzOTX4xexZPNunhzUjocur1xzxFS0wlfSTFq4Bef0RmthegWJhNCMr7Zx+T8WkpF1lDd/052rujQLdUmVQt6VNL847wwemZXKyCkryNDc8MdQuIuEQNbRHEb/dxV3vr2S+MYxzLr9AhLjTg91WZVKnWrRvHR9Z+69tDUfp+7g8n8sYt2PB0JdVthQuItUsG93ZfLrcYt4c+lWmrZ4m0bN36JRjLemEqgoeW+0vvmb7vx8OJuB4xfxryXf6TQNCneRCpOb6/jnoi386vkF7Mw4zGvDuxBZ9XtW7fDuVAIVpWuL03n/jgvo2qIe989Yw7B/LmPHgaxQlxVSCneRCvD97oNc8/ISHn4vlaSz6vHh7y8guVXlnqo33DSoXZXJw7vwyIC2fLFlN32fnc97X/9wyh7F60NMIifR0ZxcJi9O4+k5G4gw42+DEhiSGFspv1ijMjAzrk+Ko8c59blr6tfc/uZXTF+sEKcLAAAL9ElEQVSRzpjL40+5zwzoyF3kJFmyeTe/en4Bj76/jm4tTqd7p/dYvPvvCvYKcHaDWky/LYkH+rdh2ZY99HnmM8bN28Th7FNn6gIduYuUs7Sffubvczbw3tc/0OS06nRKeJc6p1fn6+1rQl3aKSUqMoKbzm/BJe0aMea9VJ6avZ63ln3PHy9uxWUJjYmI8PYfWR25i5STH/YdYvR/V/HLpz/jk9Qd3PGLc/jkrl7sObKBr7frTdNQOTOmOi9e15nXb+pK7arR/P6tlVw+biHzN+zy9Pl4HbmLnKBNOzOZtHAL01ek45zj+u7NGXnR2Z77pqTK7sJzG3D+OfV59+ttjJ29gRteXUp8kzrc1utsLok/k0iPHckr3EXKIDfXsfjb3by6aAuffrOTKlERDOrUhN9ddA6xdU+tN+4qk4gIY2DHWC5tdyYzvtrGS59tZtS/v6J5vfVc3705V3SK5fSa3pjTR+EuUgrf7f6Z6V+mM33FNrbtO0T9WlX4v97nMrR7M01BW4lUjYrkqi7NGNy5KXNStzNx/mYefX8df/toPRe3bciQxKb0OLse0ZV4jh+Fu0gJnHN8sz2DOak7mJO6g9Xb9hNhcEHLBtxzyXn0adOQatGRoS5TyigywugXfyb94s9k/fYM3lr2Pf9dsY1Zq37ktBrR9GndkH7xjeh5Tv1K9++scBcpwDnHd7sPsnTLHpZs2c0Xm/ewbd8hzKBj09MYfcl5DOjQRNMFeFCrRrV58LK2/LnfeXy2YRcfrdnOR2u2M+3LdKpERdAlri49zq5Pj7Pr0aZxHapGhXfYBxXuZtYPeA6IBF5xzj1R6PGqwOtAZ2A3cJVzLq18SxUpX9k5uWz/OZcPVv9I6g8HSP3xAGu27WdnxmEA6tWsQtcWp1On0RSanl6dGlWiWHcIbo15NsSVy8lULTqSvm0b0bdtIw5n5/D5t7tZuPEnFm76iadmrwcgOtI4r1EdEmJjaB97Guc2qk1WdnhdeRPwhJKZRQLjgEuANsA1ZtamULObgb3OuXOAZ4Any7vQMpkyBeLiICLC93PKlFBXVLSRIyEqCsx8P0eODF0t5TlmwfY1ZQrUr+97/ma++yXtN4h+c3IdOzOySP3hAJ9t2MX0L9OZ8Nm33DdjNddP+oJeT82j1f0fcc+CQ4ycsoIX565n24Kl9Fy9gEcaH2TO/13I8vt686J9w6ElU9jwzjhWvv0cK99+rnTjMmUKLFkCn332v1ve9gUfC6bPKVNg0SLYv9+3TaBxKqmWgvsLZv3q1eH7+3MSVY2KJLnVGdzXvw0f3Xkhy+7tzYtDO3Hz+WdRu1oUM1f+wJ+mr+LX4xbx+3kHi/4C7xDlUDBH7l2BTc65zQBm9hYwAEgt0GYA8JD//n+Af5iZuZNwEenOA1l8uy+HmO/3AlB4B/l7nD0bnnwKsqtD41aQDTwwFn6OxF188bHbFNdH/vKxKwK2L9yiiFHIX/Xcc/DBYlxs/P8e/GAx7vYH4I7bS9hH4BpX78rGrd9ZQg2FtpmXAv+YAFFnwFln4DB4/CU4EAW9koPaZ77P5uNefBmqNIZzGvvWPfkK7I+GCy/43zbzF8C4iVCvFdSDXDOORkZzZOy/OLq7Ckc7J3I0x3EkJ5ejObkc/Xo1R2fP58g5fTnUuiqZVaqTOWsLP3/7Lpl16pJ5OJvMw9kcOHSUon7PYqpH07xeDdo1iaF/wpkcSvuGXz9xPxNbpxLlcnnmI+C9GlBvom+DESPgysPHdvLdd771AEOHHr+TPFOmFL/9TTf5/lGHHg2uzylTYPhwGJr9v3W7d/v6CVRHUbXk7W/RIpg8OfD6I0eCe84e16B2VS5pdyaXtDsT8F019d2eg2zYkcGSFauP/2BU3rgfPOhbDva1Uw4sUP6a2WCgn3PuFv/y9UA359yoAm3W+Nuk+5e/9bcp9juwEhMT3fLly0td8EuffctfP/ym1NuJN0Tm5hCdc5To3BxqHMmi5pGD1D58iJqRjloXXUitqlHUrBpF3RrRNKhdlfq1qub/rF+7KrWqHns8k/LCCyTfcQfJw/zLr/kfaN7c9/O77/Ify9+mYJu0tOKLjYsreXsoer9F9Vmgr5WNoMP2UtRRRC3520ZGQk5OwPUPxY8l+Y9/DG5fp6iUlBSSk5OPXekf9+OcwDia2ZfOucRA7Sr0DVUzGwGMAGjYsCEpKSml7qPuwVx+28ZRvXrxb2gZwKZNBZYL/QFr2TKv1bHbHFNryXUE+rhD4e2Lbb9+ffHbtGpVuhoKrcg6dIjq1auXXEPBleu+OX4c8saudeuiayjuia1bV+w+rU2Bs3qpqcc+hiMKR5T5f3ZoT2QERBpERUDEihUA/OPwuwCMqjrgfxs3yzx2R4d9t8zdkAmkFVFL5hlnkDJ2LPsOjgcgZezxp8TyHstzTJuSXsO33x5w+yL3W1SfBfrKzv2BfWc1Dr6OImop/DwDrc+MjSVl7Njg9nWKyszMPD7T/ONe5Ov1ZI+jc67EG5AEzC6wPBoYXajNbCDJfz8K+An//wqKu3Xu3NmV1bx58wI3at7cOXC/7+e7Ofy35s3LvN+TIjLSOXAdbvXd8uuMjDzhroMap4L8Y3ZcLWUZs2DHv0C7DrcG8W/lb99rmO92ov+u855/vvj+Cu2r1PsMtH1pnkeBvmLuKcNzL27c/K+/QOvnjR0bnr8/YaTI37dyfr065xyw3AXIbedcUHPLLANamlkLM6sCXA3MLNRmJnCj//5g4FN/EaHz2GNQowYrG/n+GwtAjRq+9eHEf/5tS13frfD6CuUfs2NqKeuYBTv+jz0GVaqwspFvn/lto6OL3q+/32OcyL9rkybF91fUvkqzz5K2r1LF9xyD7fOxx45vn9dPMM+9uHEbMaJ068Pt9yfclffrtRQCnpZxzmWb2Sh8R+eRwKvOubVmNgbfX5CZwCTgX2a2CdiD7w9AaOW9WbHgZjh82HeO67HHwu/NoPF5/2V/0fcjMtL3izV+fLGbnDR5Y7P2et/xxYmMWbDjn7ecMgzfu95AvXq+N5qL2m95/7uefjpMnFhyf3mP5Ql2n4VrLbx9aZ5HaccpUC0F99ezZ+D1Var4xincfn/CXQhzKOAbqidtx2a7gCLeaQhKfXynfqRkGqfANEbB0TgFVlFj1Nw51yBQo5CF+4kws+UuiHeLT3Uap8A0RsHROAUWbmNUeWfFERGRYincRUQ8qLKG+8RQF1BJaJwC0xgFR+MUWFiNUaU85y4iIiWrrEfuIiJSAoW7iIgHVapwN7MhZrbWzHLNLLHQY6PNbJOZrTezvqGqMZyY2UNmts3MVvpvl4a6pnBiZv38r5dNZnZPqOsJR2aWZmar/a+f0s/051Fm9qqZ7fRPmpi37nQzm2NmG/0/65bUx8lWqcIdWANcAcwvuNI/v/zVQFugHzDePw+9wDPOuQ7+2wehLiZcBPk9BeJzkf/1EzbXcIeB1/BlTUH3AHOdcy2Buf7lkKlU4e6cW+ecW1/EQwOAt5xzh51zW4BN+OahFylO/vcUOOeOAHnfUyASkHNuPr6pVgoaAEz2358M/LpCiyqkUoV7CZoAWwssp/vXCYwys1X+/0aG9L+JYUavmeA44GMz+9I/ZbcUr6Fz7kf//e1Aw1AWE3ZfkG1mnwCNinjoXufcuxVdT7grabzwzUb2CL5f0EeAvwM3VVx14gHnO+e2mdkZwBwz+8Z/1ColcM45MwvpdeZhF+7Oud5l2Gwb0LTAcqx/necFO15m9jIw6ySXU5mcsq+Z0nDObfP/3Glm7+A7naVwL9oOMzvTOfejmZ0J7AxlMV45LTMTuNrMqppZC6AlsDTENYWc/wWWZyC+N6TFJ5jvKTilmVlNM6uddx+4GL2GSlLwey1uBEJ6piHsjtxLYmYDgReABsD7ZrbSOdfXP7/8VHxf2p0N/M45lxPKWsPE38ysA77TMmnAraEtJ3wU9z0FIS4r3DQE3jHfdylGAf92zn0U2pLCg5m9CSQD9c0sHXgQeAKYamY345vO/MrQVajpB0REPMkrp2VERKQAhbuIiAcp3EVEPEjhLiLiQQp3EREPUriLiHiQwl1ExIP+H2iorF7WNZpuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8leWd///Xh+x72LcEghUVjCEIAhEXnIFKrcpYcbdu7VDraMdfW6dS6lKrHR0dtTo61J/2IVaqg2O1VG1FrRkRq4gaZBNZjCZhD9lD9uv7x31yPIQk5yQEEm7ez8fjPHLu+1z3dV/nysk797nOfa7bnHOIiIi/9OvtBoiISM9TuIuI+JDCXUTEhxTuIiI+pHAXEfEhhbuIiA8p3EW6wczWmdmM3m6HSEcU7nJUMc99ZlYauN1nZtbVepxzJzrn8g9BE3uUmcWa2QYzK+7ttsjhFd3bDRA5zOYB/wRMABzwBvAFsLA3G3UI3QLsBlJ6uyFyeOnIXQ6KmRWa2S1m9qmZ1ZjZU2Y21Mz+YmZVZvammfUPKT/NzN4zs3IzWx06tGFm1waOMqvMbKuZ/SDksRlmVmxmPzGzXWa23cyu7UaTrwb+0zlX7JwrAf4TuKaD5zbIzF4JtHWvmS03s34hz3tm4H6CmS0ys7JA+/8t9Ei5G330gpntMLMKM3vHzE7sxvPEzMYAVwL/3p3t5cimcJeecCEwCzgOOA/4C/BzYDDea+xHAGY2EngVuBsYAPwUeNHMBgfq2QWcC6QC1wIPmdnJIfsZBqQBI4HvAY+1hqKZ3RoI4XZvIXWcCKwOWV4dWNeenwDFgecxNPCc2puv4w4gCzgm0A9XdrePAv4CjAWGAB8Di1sf6MLzBHg0sI99HTw/8TGFu/SER51zOwNHwsuBD5xznzjn6oCXgImBclcCrznnXnPOtTjn3gBWAecAOOdedc5tcZ7/A5YBp4fspxG4yznX6Jx7DagGjg9se69zLr2jW0gdyUBFyHIFkNzBuHsjMBwYHdjnctf+ZEwXA792zpU554qBRw6ij3DO/c45V+WcqwfuBCaYWVpXnqeZXQBEOedeaqctchRQuEtP2Blyf187y8mB+6OBi9ocaZ6GF6CY2bfM7P3AEEg5XugPCqmr1DnXFLJcG1J3pKrx3hm0SgWqOwjt+4HNwLLAMNGtHdQ5AigKWS5qp0xEfWRmUWZ2r5ltMbNKoDBQJrQfOmVmScB/sP+7ATnKKNzlcCoCft/maDPJOXevmcUBLwIPAEMDR6GvARGdyWJmPzez6o5uIUXX4X2Y2mpCYN0BAkfPP3HOHQOcD/zYzP6xnaLbgYyQ5cxI2tyBy4E5wEy8IaiswHqDiJ/n2MB2y81sB/BHYHhgHD8LOSoo3OVwehY4z8zODhyhxgc+KM0AYoE4vDM7mszsW8A3I63YOfdr51xyR7eQos/ghfRIMxuBN67+dHt1mtm5ZnZsYMimAmgGWtopugSYb2b9A58r3Bhpu9uRAtQDpUAi8OtuPM+1eP9gcgO37+O9U8il/XcV4kMKdzlsnHNFeEelP8cL8SK8U/X6Oeeq8IYRlgBleEewSw9BM34L/BlYgxeCrwbWtWcs8CbeUM7fgcedc2+3U+4uvA9evwiU/1+8gO6OZ4AvgRJgPfB+VytwzjU553a03oC9QEtgubmb7ZIjjOliHSI9y8x+CFzqnDuzt9siRy8duYscJDMbbmbTzayfmR2PN9Sjs1SkV4UNdzP7XeBLI2s7eNzM7BEz2xz4ksbJ7ZUT8bFYvKGdKuBvwJ+Ax3u1RXLUCzssY2Zn4I05PuOcy27n8XOAm/BOW5sK/MY5N/UQtFVERCIU9sjdOfcO3gcyHZmDF/zOOfc+kG5mw3uqgSIi0nU9MXHYSPY/vao4sG5724JmNg9v4iYSEhImZWZ273TglpYW+vXTxwXhqJ/CUx9FRv0U3uHqo88//3yPc25wuHKHdVZI59wTwBMAkydPdqtWrepWPfn5+cyYMaMHW+ZP6qfw1EeRUT+Fd7j6yMy+jKRcT/ybKWH/b+RlBNaJiEgv6YlwXwpcFThrZhpQ4Zw7YEhGREQOn7DDMmb2HDADGBSYo/oOIAbAObcQb/6Pc/AmWKrFm6pVRER6Udhwd85dFuZxB/xLj7VIRHpdY2MjxcXF1NXVAZCWlsaGDRt6uVV9W0/3UXx8PBkZGcTExHRre11mT0QOUFxcTEpKCllZWZgZVVVVpKToSn2d6ck+cs5RWlpKcXExY8aM6VYdOrdJRA5QV1fHwIEDaf8aJnKomRkDBw4MvnPqDoW7iLRLwd67Drb/Fe4iIj6kcBeRPunUU08NW2b58uWceOKJ5Obmsm+frgMeSuEuIn3Se++9F7bM4sWLmT9/PgUFBSQkJIQt75yjpaW9i2n5j8JdRPqk5GTvqoGtX+ufO3cuJ5xwAldccQXOOZ588kmWLFnCbbfdxhVXXAHA/fffzymnnEJOTg533HEHAIWFhRx//PFcddVVZGdnU1RUxLJly8jLy+Pkk0/moosuorrau/xsVlYWd9xxByeffDInnXQSn332GQDV1dVce+21nHTSSeTk5PDiiy8C7FfPVVddFaynL9CpkCLSqV/+eR1risqIiorqsTrHj0jljvNOjLj8J598wrp16xgxYgTTp09nxYoVfP/73+fdd9/l3HPPZe7cuSxbtoxNmzaxcuVKnHOcf/75vPPOO4waNYpNmzaxaNEipk2bxp49e7j77rt58803SUpK4r777uPBBx/k9ttvB2DQoEF8/PHHPP744zzwwAM8+eST/OpXvyItLY01a9YAUFZWdkA9d91113719DaFu4j0eVOmTCEjIwOA3NxcCgsLOe200/Yrs2zZMpYtW8bEiRMB72h706ZNjBo1itGjRzNt2jQA3n//fdavX8/06dMBaGhoIC8vL1jPd77zHQAmTZrEH//4RwDefPNNnn/++WCZ/v3788orr+xXT11dXfB+X6BwF5FO3XHeib3+Jaa4uLjg/aioKJqamg4o45xj/vz5/OAHP9hvfWFhIUlJSfuVmzVrFs8991yn++poPx3V09t91JbG3EXEF84++2x+97vfBce9S0pK2LVr1wHlpk2bxooVK9i8eTMANTU1fP75553WPWvWLB577LHgcllZWbfqOZwU7iLiC9/85je5/PLLycvL46STTmLu3LlUVVUdUG7w4ME8/fTTXHbZZeTk5JCXlxf84LQjv/jFLygrKyM7O5sJEybw9ttvH1DPzJkzw9ZzOIW9huqhoot1HHrqp/DUR+3bsGED48aNCy73tSGHvuhQ9FHb3wOAmX3knJscblsduYuI+JDCXUTEhxTuIiI+pHAXEfEhhbuIiA8p3EVEfEjhLiK+NGPGDLpyuvXTTz/NjTfe2O5jrdMPFxYWkp2dDcCqVav40Y9+BHin1H7wwQfB8gsXLuSZZ57pbtN7hKYfEJEjVnNzc49OaNaR9qYfnjx5MpMne6eb5+fnExMTw8yZMwG4/vrrD3mbwtGRu4j0SYWFhcEpfseNG8fcuXOpra0lKyuLn/3sZ5x88sm88MILFBQUMG3aNHJycrjgggsoKysL1vH73/+e3NxcsrOzWblyJQArV64kLy+PiRMncuqpp7Jx48Zg+aKiImbMmMHYsWP55S9/GVzfOv1wqPz8fM4991wKCwtZuHAhjz32GLm5uSxfvpw777yTBx54AIAtW7Ywe/ZsJk2axOmnnx78FusLL7wQ/MbrGWec0eP9pyN3EenUzX+9mY9KPurRI+TcYbk8PPvhsOU2btzIU089xfTp07nuuut4/PHHARg4cCAff/wxADk5OTz66KOceeaZ3H777fzyl7/k4Ye9umtraykoKOCdd97huuuuY+3atZxwwgksX76c6Oho3nzzTX7+858H52dfuXIla9euJTExkVNOOYVvf/vbwaPzjmRlZXH99dcTExPDggULAHjrrbeCj8+bN4+FCxcyduxYPvjgA2644Qb+9re/cdddd/H6668zcuRIysvLu96JYSjcRaTPyszMDE6je+WVV/LII48AcMkllwBQUVFBeXk5Z555JgBXX301F110UXD7yy67DIAzzjiDyspKysvLqaqq4uqrr2bTpk2YGY2NjcHys2bNYuDAgYA39e+7774bNtw7U11dzXvvvbdfm+rr6wGYPn0611xzDRdffHFwmuGepHAXkU49PPvhXptbxszaXQ6dwrer2992222cddZZvPTSSxQWFu43t1BH++uulpYW0tPTKSgoOOCxhQsX8sEHH/Dqq68yadIkPvroo+A/lp6gMXcR6bO++uor/v73vwPwhz/84YALdKSlpdG/f3+WL18OeGPsrUfxAP/zP/8DwLvvvktaWhppaWlUVFQwcuRIwDtDJtQbb7zB3r172bdvHy+//HLEF99ISUlpdwbK1NRUxowZwwsvvAB4c8CvXr0a8Mbip06dyl133cXgwYMpKiqKaF+RUriLSJ91/PHH89hjjzFu3DjKysr44Q9/eECZRYsWccstt5CTk0NBQcF+l7mLj49n4sSJXH/99Tz11FMA/Nu//Rvz589n4sSJB1yMY8qUKVx44YXk5ORw4YUXRjwkc9555/HKK68EP1ANtXjxYp566ikmTJjAiSeeyJ/+9CcAbrnlFk466SSys7M59dRTmTBhQpf6JhxN+etj6qfw1Eft6wtT/hYWFnLuueeydu3aw7rf7tKUvyIicsgp3EWkT8rKyjpijtr7IoW7iIgPKdxFRHxI4S4i4kMKdxERH4oo3M1stpltNLPNZnZrO4+PMrO3zewTM/vUzM7p+aaKyNEudMrdwyk/P3+/mSH7wpS+4YSdfsDMooDHgFlAMfChmS11zq0PKfYLYIlz7r/NbDzwGpB1CNorInLY5efnk5ycHJzXvS9M6RtOJEfuU4DNzrmtzrkG4HlgTpsyDkgN3E8DtvVcE0Wkz1u8GLKyoF8/7+fixT1S7YMPPkh2djbZ2dnBmR6bmpoOmAYY4NZbb2X8+PHk5OTw05/+FIDdu3dz4YUXcsopp3DKKaewYsUKWlpayMrK2m8mxrFjx7Jz507+/Oc/M3XqVCZOnMjMmTPZuXNncErfhx56qN0pfVunHM7Ly9tvyuEZM2bws5/9jClTpnDccccFv7m6bt06pkyZQm5uLjk5OWzatKlH+qqtSCYOGwmETnpQDExtU+ZOYJmZ3QQkATPbq8jM5gHzAIYOHUp+fn4Xm+uprq7u9rZHE/VTeOqj9qWlpe03V0pzc3O7c6cARC9ZQvxNN2H79nkrvvwS98//TF1dHU0XX9ztNnzyySc89dRTvPXWWzjn+Id/+AcmT57Mxo0befTRR1m4cCE33HADDz30EFdccQUvvvgiH330EWYWnP3xhhtu4Ac/+AF5eXkUFRVxwQUXsGrVKr71rW/x3HPPceWVV/Lhhx+SkZFBYmIiEyZM4I033sDMWLRoEXfffTe//vWvufbaa0lOTg5eeem1114jJiaGqqoqrrzySu6//37y8vL493//dxYsWMB9991Hc3MztbW1vPXWW7z++uvcfvvtLF26lEceeYR58+ZxySWX0NDQ0Gnf1tXVdf/16Zzr9AbMBZ4MWf4u8F9tyvwY+Engfh6wHujXWb2TJk1y3fX22293e9ujifopPPVR+9avX7/fcmVlZceFR492Dg68jR59UG14+OGH3W233RZc/sUvfuF+85vfuMzMzOC6t956y82ZM8c1Nja6nJwcd+2117oXX3zR1dfXO+ecGzx4sJswYULwNmLECFdVVeVWrFjhzj77bOecczfffLN74oknnHPOffrpp27WrFkuOzvbHXfcccEyd9xxh7v//vuD+21dLi8vD7ansrLSbd682U2cONE559yZZ57p3n33Xeecczt27HDf+MY3nHPOLV682I0fP97de++97vPPP++0D9r+HpxzDljlwuS2cy6iYZkSIDNkOSOwLtT3gCWBfxZ/B+KBQd37dyMiR5Svvura+oPU3rS80dHRrFy5krlz5/LKK68we/ZswJty9/3336egoICCggJKSkpITk4mLy+PzZs3s3v3bl5++eXgfOo33XQTN954I2vWrOG3v/0tdXV1B9XWuLg4AKKiooKTlF1++eUsXbqUhIQEzjnnHP72t78d1D46Ekm4fwiMNbMxZhYLXAosbVPmK+AfAcxsHF647+7JhopIHzVqVNfWR+j000/n5Zdfpra2lpqaGl566SVOP/30dqcBrq6upqKignPOOYeHHnooOK3uN7/5TR599NFgna3zqpsZF1xwAT/+8Y8ZN25ccB710OmAFy1aFNyuoyl9w0053J6tW7dyzDHH8KMf/Yg5c+bw6aefdreLOhU23J1zTcCNwOvABryzYtaZ2V1mdn6g2E+Afzaz1cBzwDWBtw8i4nf33AOJifuvS0z01h+Ek08+mWuuuYYpU6YwdepUvv/979O/f/92pwGuqqri3HPPJScnh9NOO40HH3wQgEceeYRVq1aRk5PD+PHjWbhwYbD+Sy65hGeffTZ4VSeAO++8k4suuohJkyYxaNDXgw/nnXceL730UrtT+rZOOZyXl3fAlMPtWbJkCdnZ2eTm5rJ27Vquuuqqg+qnDkUydnMobhpzP/TUT+Gpj9rXpTF355x79llvjN3M+/nss4esbX1V2D7qhoMZc9dl9kTk4F1xhXeTPkPTD4iI+JDCXUTa5fSxWa862P5XuIvIAeLj4yktLVXA9xLnHKWlpcTHx3e7Do25i8gBMjIyKC4uZvdu74zmurq6gwqao0FP91F8fDwZGRnd3l7hLiIHiImJYcyYMcHl/Px8Jk6c2Ist6vv6Wh9pWEZExIcU7iIiPqRwFxHxIYW7iIgPKdxFRHxI4S4i4kMKdxERH1K4i4j4kMJdRMSHFO4iIj6kcBcR8SGFu4iIDyncRUR8SOEuIuJDCncRER9SuIuI+JDCXUTEhxTuIiI+pHAXEfEhhbuIiA8p3EVEfEjhLiLiQwp3EREfUriLiPiQwl1ExIcU7iIiPqRwFxHxoYjC3cxmm9lGM9tsZrd2UOZiM1tvZuvM7A8920wREemK6HAFzCwKeAyYBRQDH5rZUufc+pAyY4H5wHTnXJmZDTlUDRYRkfAiOXKfAmx2zm11zjUAzwNz2pT5Z+Ax51wZgHNuV882U0REuiKScB8JFIUsFwfWhToOOM7MVpjZ+2Y2u6caKCIiXRd2WKYL9YwFZgAZwDtmdpJzrjy0kJnNA+YBDB06lPz8/G7trLq6utvbHk3UT+GpjyKjfgqvr/VRJOFeAmSGLGcE1oUqBj5wzjUCX5jZ53hh/2FoIefcE8ATAJMnT3YzZszoVqPz8/Pp7rZHE/VTeOqjyKifwutrfRTJsMyHwFgzG2NmscClwNI2ZV7GO2rHzAbhDdNs7cF2iohIF4QNd+dcE3Aj8DqwAVjinFtnZneZ2fmBYq8DpWa2HngbuMU5V3qoGi0iIp2LaMzdOfca8FqbdbeH3HfAjwM3ERHpZfqGqoiIDyncRUR8SOEuIuJDCncRER9SuIuI+JDCXUTEhxTuIiI+pHAXEfEhhbuIiA8p3EVEfEjhLiLiQwp3EREfUriLiPiQwl1ExIcU7iIiPqRwFxHxIYW7iIgPKdxFRHxI4S4i4kMKdxERH1K4i4j4kMJdRMSHFO4iIj6kcBcR8SGFu4iIDyncRUR8SOEuIuJDCncRER9SuIuI+JDCXUTEhxTuIiI+pHAXEfEhhbuIiA8p3EVEfCiicDez2Wa20cw2m9mtnZS70MycmU3uuSaKiEhXhQ13M4sCHgO+BYwHLjOz8e2USwH+FfigpxspIiJdE8mR+xRgs3Nuq3OuAXgemNNOuV8B9wF1Pdg+ERHphugIyowEikKWi4GpoQXM7GQg0zn3qpnd0lFFZjYPmAcwdOhQ8vPzu9xggOrq6m5vezRRP4WnPoqM+im8vtZHkYR7p8ysH/AgcE24ss65J4AnACZPnuxmzJjRrX3m5+fT3W2PJuqn8NRHkVE/hdfX+iiSYZkSIDNkOSOwrlUKkA3km1khMA1Yqg9VRUR6TyTh/iEw1szGmFkscCmwtPVB51yFc26Qcy7LOZcFvA+c75xbdUhaLCIiYYUNd+dcE3Aj8DqwAVjinFtnZneZ2fmHuoEiItJ1EY25O+deA15rs+72DsrOOPhmiYjIwdA3VEVEfEjhLiLiQwp3EREfUriLiPiQwl1ExIcU7iIiPqRwFxHxIYW7iIgPKdxFRHxI4S4i4kMKdxERH1K4i4j4kMJdRMSHFO4iIj6kcBcR8SGFu4iIDyncRUR8SOEuIuJDCncRER9SuIuI+JDCXUTEhxTuIiI+pHAXEfEhhbuIiA8p3EVEfEjhLiLiQwp3EREfUriLiPiQwl1ExIcU7iIiPqRwFxHxIYW7iIgPKdxFRHxI4S4i4kMRhbuZzTazjWa22cxubefxH5vZejP71MzeMrPRPd9UERGJVNhwN7Mo4DHgW8B44DIzG9+m2CfAZOdcDvC/wH/0dENFRCRykRy5TwE2O+e2OucagOeBOaEFnHNvO+dqA4vvAxk920wREemK6AjKjASKQpaLgamdlP8e8Jf2HjCzecA8gKFDh5Kfnx9ZK9uorq7u9rZHE/VTeOqjyKifwutrfRRJuEfMzK4EJgNntve4c+4J4AmAyZMnuxkzZnRrP/n5+XR326OJ+ik89VFk1E/h9bU+iiTcS4DMkOWMwLr9mNlMYAFwpnOuvmeaJyIi3RHJmPuHwFgzG2NmscClwNLQAmY2EfgtcL5zblfPN1NERLoibLg755qAG4HXgQ3AEufcOjO7y8zODxS7H0gGXjCzAjNb2kF1IiJyGEQ05u6cew14rc2620Puz+zhdomIyEHQN1RFRHxI4S4i4kMKdxERH1K4i4j4kMJdRMSHFO4iIj6kcBcR8SGFu4iIDyncRUR8SOEuIuJDCncRER9SuIuI+JDCXUTEhxTuIiI+1KOX2RPpTc45ahuaKattoKa+mZqGJmrqm7z79U3UNTXT3OL2u23a2sAGthATZcRG9yMmqvVmJMREkRIfQ0p8NKmBnynx0URH6ZhI+j6FuxwRquoaKSnfx7byfZSU7aOkvI7dVfWU1tSzt6aB0uoGSmvqqWts6Xrln3/WpeKp8dEMTokL3OIZnOzdH5ISR0b/BDIHJDI0NZ6oftb1toj0EIW79BnltQ1s2V3DF3tq2Lq7mq27aygsraGkfB9VdU37lY2JMgYnxzEwOY4BSbEcOySZQYH7/RNjSI6LISkuiqS4aBJjo0iOiyY+JoqofkaUGVFR3s8V7y7n9NPP4CfL/j9anOPOM+6nsbmF+qYW6hqbqaxr5P7359PQ7GhoaqahyXHWsFvYXVXP7qp61pZUsLuqnur6A9s3Mj2BjP6JZA5IIGtgEscOSeYbg5PJHJCo4JdDTuEuh11VXSOf7ajis+2VrN9exec7q9i6u5qy2sZgmeh+xqgBiWQNSmLKmAGMTE9gRHoCI/snMDI9gV+v+BlmxsOzHz6otsRGGQmxUazf8ykAw9LiDyhTtnyTd6cfWCzcef6JB5SpbWhiZ2U9RXtrKS7bR1FZLUV7aykq28eydTsprWkI2Wc/xgxqDfskjhuWwrjhqWQNTFLoS49RuMsh09Li+HJvLRu2VwaD/LMdlRSX7QuWSUuI4fihKUQPeIaWpI+ZPPw0Hpj1IJkDEonpZGx79c7Vh+MpRCwxNpoxg6IZMyjpgMdu/uvN1De1UFHbQPm+Rk4f/FO27K5m/fZK/rJ2Oy3OK5cQE8Xxw1IYPyKV8cNTGTc8lROGpZAUpz9T6Tq9aqTH7Kmup+CrclYXl1NQVM7qonIqA8Mp/QyOGZxMbmY6l00ZxbjhKZwwLJXhafGYGTOenk9l9Tb2NnzOMYOTe/mZ9KyCHQX7Lc8/Z1zw/k2v/Svl+xo4f/TPWb+9kg3bK3ll9Tb+8MFXAJjBmEFJ5IxMY0JmOjkZ6Zw4IpX4mKjD+hzkyKNwl27Z19DM2m0VrC4q55Oicgq+Kqek3Dsij+pnHDc0hW/njKCg8mEGJsfy5Jz/UiC1Y80u7x3IRZMzg+ucc5SU72PD9irWb6tkTUkF720p5eWCbYA3ZHXc0BQmZKYzISONnIx0jhuarLN4ZD8KdwmrucWxZXc1BV+VU1DsBfnGnVU0B8YTRqYnkJuZzjWnZjEhM53skakkxnovrRlPb6G4GgV7F5gZGf0TyeifyKzxQ4Prd1TU8cNXbmJPdT0Dov6FVz/dxnMrvSP8+Jh+ZI/wju4nZKYzMTOdjP4JmGkM/2ilcJcD7Kys45PW4ZWvyllTUhE8GyQlPpoJGen88MxvBIIkjSEpB34IKT1vWFo8ZY2b+KJmM5PGD2D19x6isLSWTwPDYJ8WV/Ds+1/y1LtfADAwKTZwdO/9nnIz00lPjO3lZyGHi8L9KFdd38Sa4orgGHlBUTk7KusA7+3/+BGpXDBxJLmBI8JjBiXRT2d09KrqhmoKdhRgZowZlMSYQUnMyR0JQGNzCxt3VO33+3x74y5c4EPbrIGJwd/lhMx0xg/X+L1fKdyPIq1/+KuLvT/81UUVbNpVFTxbY/TARKaMGcDL264gITaK1dcX6A//CBMT1Y/skWlkj0zjymmjAe/U0zXFFRQEfu9/3/r1+H1MlDFueKoX+Bnp5I5KZ8xA/QP3A4W7T7W0OHbUtPDyJyWBt+zlrNtWSX2T9w3O/okx5GSkMzt7GLmjvD/sAUneW/Zn7t1GZaPGyf0iJT6GU48dxKnHDuLmv95M3mC4Ne9eCorKKCiqoKCojBc/KuaZv38ZKO8NvbUe4edmpvfyM5DuULj7gHOOXVX1rA6Mu7YemXunIRaQEBNF9shUvjttdHAMNnOAPmw7GrWeljksLZ7ZacOZnT0cgH/9y79Svq+Rb2fOD35o/t//tyX4ofnAeGNaycdMyExjQkY6J2WkBT80l75Jv50jTHOLo7C0hnXbKlm/rZJ12yrYsL2SPdXeNyCj+hnHB05DjK/ZwcUzpzJ2iE6Tk861fins4lMyufgU77TMfQ3NXPfyv7C7qp4BVZewuricV9dsB7zvLRw3NIXcTO+8+3E5b91FAAALpklEQVTDUzlheCrJ+sJVn6HfRB9WU9/Epl3VbNjuhfj6bZVs2F7FvsZmwBsvHTskhbOOH8L4EankZKQxfngaCbHecEp+finjhqf25lOQI1hCbBTbaz+DKLghN54ZM2awp7o+8HlNOQXFFfxl7Q6e/7AouM2oAYmcEJhOwbulkNk/UWP4vUDh3gdU1TWyeVc1m3ZWs2lXFZsC91u/FASQEhfNuBGpXDolk/HDUzlxRBrHDkkmNlpH5HL4DEqO4x/HDeUfx3nn3zvn2FZRx2eBb9du2F7Fhu2VvLFhZ/AMnaTYKE4YnsrYIcnBydOOHZLMyPQEhf4hpHA/TBqbWygu20fhHm+mwy9La9m6p4bNO6vYVlEXLBcb3Y9jByczOas/lw8dxbFDkhk3LFVj5NInmXmzX45MTwgGPnhDOht3VgUCv5LPtlexbP3O/Y7y46L7ccxgb/K01tA/ZnASowYkkhIf0xtPx1cU7j3EOcee6ga2V3hzjheX7ePL0tpgkJeU7wt+OAXe0UzrjIdjh6bw4Op/Ij4mik9/WKCZAeWIlxAbRW47Z9rsrWlgy+5qtuyqZvOuarbsrg6O5buv/zzonxjDqAGJZA5IZFTILXNAIsPT4vUZUgQU7hFoam6htKbBm8O7up6dFXXeRSPK64Jhvq2ijoam/S8UkRLvzRI4ITOdObkjGD0wiayBiYwemMSg5Nj9jsQXfFAMdSjYxdcGJMUyIGkAp2QN2G99XWMzW3fX8GVpDV/trQ3e1pZU8Ne1O2gKOTDqZzA4JY5haQkMS41jeFoCw9LiGZ4Wz7DUeIanJTAkNe6oP5U3onA3s9nAb4Ao4Enn3L1tHo8DngEmAaXAJc65wp5tas+pa2ymvLaR8n0N3s/aRioC9/dU17OnuiF4MYbd1fWU1Tbsd1QB3gtsaKr3gsoemcbZJw5jRHoCw9PivXnH0xNIT4zRUIpIBOJjorypjkcceAJAU3MLOyrrvMAvrWVbRR07KvaxvaKOrbtreG9L6QEXcwHvc6oBybEMTIplYHIcg5JjGZjkXdBlYHIsg5LjSEuIIS0hhtSEGFLion31GUDYcDezKOAxYBZQDHxoZkudc+tDin0PKHPOHWtmlwL3AZccigZX7GtkW3ULBUXl1NY3UV3fRG1Dc+BnE9WB62XWNjRRvXELtes3UtVilKcMoHzQMMqJPuAIO1RcdL/gJdRGD0xkUlb/4GXUBqfEMSjZu5zasLT4Tucb75IbboC0Cu9+dDTMmwePP94zdXfV4sVQWQnOQVYW3HMPXHFF9+t6/32or4c7O6lr8WJYsQIGNcHG/4OfDoLf/Kbj/UZab0+0M/SxVl3ZZ2fbQ9eeR1f7KdLnGcn6qjVQUnJw/dxN0VH9ghOpnfqN9stU1zexo6KOHRXeu+mdlXWUhlx+sWhvLQVF5eytadhveDSUGSTHRXthHx9DasLX91PiY0iOiyIxLpqk2CgSYr2frcuJsdHsqu0gVxYvhgUL4KuvYNSog3+9RiiSI/cpwGbn3FYAM3semAOEhvsc4M7A/f8F/svMzLm2x7sH7w8ffMV97+6Dd1d0WCYpNorElkaSd+8hKSqB5KZajtn5BeklG0g7czppkyeQnhBLemIM6QkxpCV6/73TE2NJio06vEfbN9wA//3fcGtgubnZW4bDH/CLF3v/WH4U+LV9+aW3DF1/MbbWdXF953UtXgzXXQeXhxx5lZbCtde2v99I643U3r0d1wf7P9Yq0n22bWvo9tdd5/0DvaIxsjoXL/b65Io2/XTddeHb0V5bWve3YgUsWhR+fUPDwfXzIZYcF82xgbNxOtPS4qjY10hpjfcOvWJfI5X7Gr2fdU1UBpYr67x1hXtqA481UtvQ3GndSTFw8TltVrb2e22tt3ywr9cuiCTcRwJFIcvFwNSOyjjnmsysAhgI7OmJRoaaOW4Ildu/4JSJJ5EUG01SXOAW610vMyEmyntrlZUFX37JzbO97R7+a6CCNX+EwsKeblb3PfFEx+sPd7gvWPD1i7BVba23vqsvxEjrWrDAC462Ghvb329PthG8o9GO6mu9355I9tleW1u195w7q3PBAq9P2qsnkufeUb898YR3QBHp+u72cx/Rr5/RPymW/kmxHDuka9u2tDjqmpqpqW9mX0MzNQ3eCEFNfTO1DU18unbdgRv19Ou1CyzcwbWZzQVmO+e+H1j+LjDVOXdjSJm1gTLFgeUtgTJ72tQ1D5gHMHTo0EnPP/98txpdXV1NcnKYq/V89BEAN9d6Aflw4g1fPzZpUrf2e0gE2nlu9S8AeCX57q8fO8h2RtRPh6otkfZ/SLnNLds4tt+Izn9XPfx7rd65k+Ti4vbrC2h9rFXE+2zT1va2j/h5dLWfwrSl7fMMt/7u435NcnFxZPs6SrX79xbo9/+q/xMAN8bN+fqxbvbjWWed9ZFzbnLYgs65Tm9AHvB6yPJ8YH6bMq8DeYH70XhH7NZZvZMmTXLd9fbbb4cvNHq0c+DOvMa7OQK30aO7vd9DIirKOXBpt3q3YDujog666oj6KVSgzw5oS3f6LNL+DymXdmsEv6se/r2+/cgjHdfXZl9d3me47bvyPLraT2Hast/rLIL1bz/wQN/8++lD2v17OwQ5BKxyYXLbOUcknwh+CIw1szFmFgtcCixtU2YpcHXg/lzgb4FG9J577oHExP3XJSZ+/UFWXxE6vhvJ+kOpJ/ss0rruuQdi27mARExM+/vt6d/ryJEd19fevrqyz862j431nmOkdd5zz4HlW+uJ5Ll31G/z5nVtfV/7++nrejGHwo65O28M/Ua8o/Mo4HfOuXVmdhfef5ClwFPA781sM7AX7x9A72odz1r+Pe/T/tGjD9un1F0SHFcPfIgaFdV7Z8u09s2673rHFwfTZ5H2f+ty/jVA4MPCgQM7Pgukp3+vAwZ448ud1df6WKtI99m2rW2378rz6Go/hWtL6P6mTw+/PjbW66e+9vfT1/ViDoUdcz9kOzbbDXzZzc0HcQg+rPUh9VN46qPIqJ/CO1x9NNo5NzhcoV4L94NhZqtcJB8oHOXUT+GpjyKjfgqvr/WRJmgQEfEhhbuIiA8dqeHewTd/pA31U3jqo8ion8LrU310RI65i4hI547UI3cREemEwl1ExIeOqHA3s4vMbJ2ZtZjZ5DaPzTezzWa20czO7q029iVmdqeZlZhZQeDWds66o5qZzQ68Xjab2a3htzj6mFmhma0JvH5W9XZ7+goz+52Z7QrMq9W6boCZvWFmmwI/+/dmG4+ocAfWAt8B3gldaWbj8b4VeyIwG3g8MA+9wEPOudzA7bXebkxfEXKdgm8B44HLAq8jOdBZgddPnzmHuw94Gi9rQt0KvOWcGwu8xdcTefeKIyrcnXMbnHMb23loDvC8c67eOfcFsBlvHnqRjgSvU+CcawBar1MgEpZz7h28qVZCzQEWBe4vAv7psDaqjSMq3DvR3pzzI3upLX3NjWb2aeBtZK++Texj9JqJjAOWmdlHgSm7pWNDnXPbA/d3AEN7szF97gLZZvYmMKydhxY45/50uNvT13XWX3izkf0K7w/0V8B/AtcdvtaJD5zmnCsxsyHAG2b2WeCoVTrhnHNm1qvnmfe5cHfOzezGZiVAZshyRmCd70XaX2b2/wOvHOLmHEmO2tdMVzjnSgI/d5nZS3jDWQr39u00s+HOue1mNhzY1ZuN8cuwzFLgUjOLM7MxwFhgZS+3qdcFXmCtLsD7QFo8kVyn4KhmZklmltJ6H/gmeg11JvS6FlcDvTrS0OeO3DtjZhcAjwKDgVfNrMA5d3ZgfvkleBftbgL+xTnX+dVsjw7/YWa5eMMyhcAPerc5fUdH1yno5Wb1NUOBlwIXjI8G/uCc+2vnmxwdzOw5YAYwyMyKgTuAe4ElZvY9vOnML+69Fmr6ARERX/LLsIyIiIRQuIuI+JDCXUTEhxTuIiI+pHAXEfEhhbuIiA8p3EVEfOj/AcEADm3+sdNiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "obs_mu, obs_sigma = 0, 4\n",
    "observations = np.random.normal(obs_mu, obs_sigma, 20)\n",
    "draw_likelihood(observations, mu=5, sigma=2)\n",
    "draw_likelihood(observations, mu=0, sigma=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. 我们假设有一组观测到的数据, 一共有$N$个, $X_{observation} = \\{ x_1, x_2, ... , x_N\\}$   \n",
    "我们推断这一组数据属于一个概率分布, 是一个正态分布, 它的概率密度函数为$f(x;\\mu,\\sigma)$, 我们将$X_{observation}$里的数据点带入到$f(x;\\mu,\\sigma)$里, 得到每个数据点在我们假设的概率分布中的出现的可能性, 注意我们接下来简写$X_{observation}$为$X$:\n",
    "$$P(x_1 \\mid \\mu,\\sigma), \\,P(x_2 \\mid \\mu,\\sigma)\\, , ... \\, ,\\,P(x_N \\mid \\mu,\\sigma)$$\n",
    "3. 那么这一组数据$X_{observation} = \\{ x_1, x_2, ... , x_N\\}$在假设的概率分布中的出现的可能性就是他们概率的乘积:\n",
    "$$L(\\mu ,\\sigma \\mid X)=P(X \\mid \\mu ,\\sigma)=\\prod _{i=1}^{N}P(x_{i}\\mid \\mu,\\sigma)$$\n",
    "上式中, 我们用$L(\\mu ,\\sigma \\mid X)$来表示似然函数, 已知的观测到的数据点$X$, 用似然函数$L(\\mu ,\\sigma \\mid X)$来估计参数$\\mu ,\\sigma$的可能性, 由此可见，似然函数也是一种条件概率函数，但我们关注的变量改变了.   \n",
    "   \n",
    "4. 从下面两图可得, 上图图中所假设的概率分布求出的似然函数取值显然比下面的小, 因为有很接近于$0$的数值, 这让似然函数乘积的结果变得非常小. 我们可以得出结论, 下图的概率分布参数更符合观测到的数据点的概率分布,而**最大似然估计**的目的就是找到一个最符合当前数据的分布的参数."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们做个试验:\n",
    "1. 先从均值为$0$, 标准差为$4$正态分布中随机抽取$200$个数据点, 作为我们观测到的数据;\n",
    "之后我们定义一个估计参数的取值范围, 均值在$(-0.5, 0.5)$之间, 标准差在$(3.5, 4.5)$之间;   \n",
    "2. 然后我们将所有的数据点和我们估计的每一组参数代入似然函数, 也就是$L(\\mu ,\\sigma \\mid X)=\\prod _{i=1}^{N}P(x_{i}\\mid \\mu,\\sigma)$中, 求得每一组参数的似然值, 3.下图可见, 似然函数的值约在$\\mu=0, \\ \\sigma=4$时取得极大值(因为数据点太少有些误差);\n",
    "4. 我们发现似然函数的图像是凸函数, 我们就可以用很多优化的方法求它的极大值了.\n",
    "5. 找到似然函数极大值的过程就是极大似然估计的过程."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 首先定义观测的数据分布, 我们定一个一个均值为0, 标准差为4的正态分布, \n",
    "# 并从中随机抽200个数据点作为观测到的数值\n",
    "obs_mu, obs_sigma = 0, 4\n",
    "observations = np.random.normal(obs_mu, obs_sigma, 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "def likelihood(observations, infer_mu, infer_sigma):\n",
    "    # 定义似然函数, observations为观测到的数据点\n",
    "    # infer_mu, infer_sigma为推断的均值和标准差\n",
    "    product_ = 1\n",
    "    for obs in observations:\n",
    "        # 代入每一个数据点到我们假设的概率密度函数内, 并求它们的积\n",
    "        product_ *= stats.norm.pdf(obs, loc=infer_mu, scale=infer_sigma)\n",
    "    return product_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 我们定义一个参数的取值范围, 均值在(-0.5, 0.5)之间\n",
    "# 标准差在(3.5, 4.5)之间\n",
    "all_infer_mu = [i/10 for i in range(-5, 5)]\n",
    "all_infer_sigma = [i/10 for i in range(35, 45)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 求得每一个参数组合的似然值\n",
    "mle = []\n",
    "for infer_mu in all_infer_mu:\n",
    "    temp_lis = []\n",
    "    for infer_sigma in all_infer_sigma:\n",
    "        temp_lis.append(likelihood(observations, infer_mu, infer_sigma))\n",
    "    mle.append(temp_lis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.plotly.v1+json": {
       "data": [
        {
         "type": "surface",
         "x": [
          -0.5,
          -0.4,
          -0.3,
          -0.2,
          -0.1,
          0,
          0.1,
          0.2,
          0.3,
          0.4
         ],
         "y": [
          3.5,
          3.6,
          3.7,
          3.8,
          3.9,
          4,
          4.1,
          4.2,
          4.3,
          4.4
         ],
         "z": [
          [
           1.460825681631413e-249,
           1.0373120935883478e-248,
           4.682481547217632e-248,
           1.419451486742629e-247,
           3.03023899775314e-247,
           4.7476206183476214e-247,
           5.6591073457077994e-247,
           5.296104474699062e-247,
           4.000180909092284e-247,
           2.4982970345070237e-247
          ],
          [
           4.476409474673629e-249,
           2.9894940695276932e-248,
           1.275416701641574e-247,
           3.670235624951869e-247,
           7.467305457483984e-247,
           1.1189830134533918e-246,
           1.2798339607468397e-246,
           1.1526214309907195e-246,
           8.40009226031549e-247,
           5.074219710609222e-247
          ],
          [
           1.1650806670689444e-248,
           7.383549620589285e-248,
           3.0017958929467434e-247,
           8.262581349916492e-247,
           1.613409129439688e-246,
           2.327470419846075e-246,
           2.5697351244185605e-246,
           2.2396342461337547e-246,
           1.583114519199049e-246,
           9.294570610141049e-247
          ],
          [
           2.5755913849581866e-248,
           1.5628307442373045e-247,
           6.104684901647834e-247,
           1.6195186743637765e-246,
           3.0564596192454286e-246,
           4.2722635540242456e-246,
           4.580913843902494e-246,
           3.8853287510968715e-246,
           2.6777156728008803e-246,
           1.5354096355898816e-246
          ],
          [
           4.83607146694548e-248,
           2.8349012866247653e-247,
           1.0727496457009362e-246,
           2.7637870756392706e-246,
           5.076756973695156e-246,
           6.920621412681512e-246,
           7.25011540178241e-246,
           6.017813259175593e-246,
           4.06481179228702e-246,
           2.2874629048579498e-246
          ],
          [
           7.712643101726141e-248,
           4.407001339848886e-247,
           1.6288704370282895e-246,
           4.1064979727080796e-246,
           7.393458592634441e-246,
           9.893394746319798e-246,
           1.0187479491605617e-245,
           8.321658792255388e-246,
           5.537838366397335e-246,
           3.0733974441126817e-246
          ],
          [
           1.0447407591638561e-247,
           5.87121076092514e-247,
           2.137114232999409e-246,
           5.312354564801568e-246,
           9.440662565341089e-246,
           1.2481269857029873e-245,
           1.2709182264773008e-245,
           1.0274042676853264e-245,
           6.771173818490415e-246,
           3.7240736313842286e-246
          ],
          [
           1.2020114932797596e-247,
           6.703342433119301e-247,
           2.42282455330681e-246,
           5.983439501146988e-246,
           1.0569412905956128e-245,
           1.3895858986836018e-245,
           1.40765912878501e-245,
           1.1324866440591987e-245,
           7.430388161429455e-246,
           4.0696087659692717e-246
          ],
          [
           1.1746364697984904e-247,
           6.558951005692979e-247,
           2.3733906647249843e-246,
           5.867636000603116e-246,
           1.0375112463729928e-245,
           1.3652909555156847e-245,
           1.3842241650741372e-245,
           1.1145128373283335e-245,
           7.317839732550941e-246,
           4.0107161668927876e-246
          ],
          [
           9.749741907642282e-248,
           5.499923069240653e-247,
           2.0089518483418072e-246,
           5.009839752459894e-246,
           8.929524248058134e-246,
           1.1837996790456628e-245,
           1.2084936832477347e-245,
           9.792584422798124e-246,
           6.468121041765527e-246,
           3.56472526918372e-246
          ]
         ]
        }
       ],
       "layout": {
        "scene": {
         "xaxis": {
          "title": "mean"
         },
         "yaxis": {
          "title": "sigma"
         },
         "zaxis": {
          "title": "likelihood"
         }
        },
        "title": "Likelihood"
       }
      },
      "text/html": [
       "<div id=\"cc3ffd98-82f6-466c-a8a5-0440824d0070\" style=\"height: 525px; width: 100%;\" class=\"plotly-graph-div\"></div><script type=\"text/javascript\">require([\"plotly\"], function(Plotly) { window.PLOTLYENV=window.PLOTLYENV || {};window.PLOTLYENV.BASE_URL=\"https://plot.ly\";Plotly.newPlot(\"cc3ffd98-82f6-466c-a8a5-0440824d0070\", [{\"type\": \"surface\", \"x\": [-0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4], \"y\": [3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4], \"z\": [[1.460825681631413e-249, 1.0373120935883478e-248, 4.682481547217632e-248, 1.419451486742629e-247, 3.03023899775314e-247, 4.7476206183476214e-247, 5.6591073457077994e-247, 5.296104474699062e-247, 4.000180909092284e-247, 2.4982970345070237e-247], [4.476409474673629e-249, 2.9894940695276932e-248, 1.275416701641574e-247, 3.670235624951869e-247, 7.467305457483984e-247, 1.1189830134533918e-246, 1.2798339607468397e-246, 1.1526214309907195e-246, 8.40009226031549e-247, 5.074219710609222e-247], [1.1650806670689444e-248, 7.383549620589285e-248, 3.0017958929467434e-247, 8.262581349916492e-247, 1.613409129439688e-246, 2.327470419846075e-246, 2.5697351244185605e-246, 2.2396342461337547e-246, 1.583114519199049e-246, 9.294570610141049e-247], [2.5755913849581866e-248, 1.5628307442373045e-247, 6.104684901647834e-247, 1.6195186743637765e-246, 3.0564596192454286e-246, 4.2722635540242456e-246, 4.580913843902494e-246, 3.8853287510968715e-246, 2.6777156728008803e-246, 1.5354096355898816e-246], [4.83607146694548e-248, 2.8349012866247653e-247, 1.0727496457009362e-246, 2.7637870756392706e-246, 5.076756973695156e-246, 6.920621412681512e-246, 7.25011540178241e-246, 6.017813259175593e-246, 4.06481179228702e-246, 2.2874629048579498e-246], [7.712643101726141e-248, 4.407001339848886e-247, 1.6288704370282895e-246, 4.1064979727080796e-246, 7.393458592634441e-246, 9.893394746319798e-246, 1.0187479491605617e-245, 8.321658792255388e-246, 5.537838366397335e-246, 3.0733974441126817e-246], [1.0447407591638561e-247, 5.87121076092514e-247, 2.137114232999409e-246, 5.312354564801568e-246, 9.440662565341089e-246, 1.2481269857029873e-245, 1.2709182264773008e-245, 1.0274042676853264e-245, 6.771173818490415e-246, 3.7240736313842286e-246], [1.2020114932797596e-247, 6.703342433119301e-247, 2.42282455330681e-246, 5.983439501146988e-246, 1.0569412905956128e-245, 1.3895858986836018e-245, 1.40765912878501e-245, 1.1324866440591987e-245, 7.430388161429455e-246, 4.0696087659692717e-246], [1.1746364697984904e-247, 6.558951005692979e-247, 2.3733906647249843e-246, 5.867636000603116e-246, 1.0375112463729928e-245, 1.3652909555156847e-245, 1.3842241650741372e-245, 1.1145128373283335e-245, 7.317839732550941e-246, 4.0107161668927876e-246], [9.749741907642282e-248, 5.499923069240653e-247, 2.0089518483418072e-246, 5.009839752459894e-246, 8.929524248058134e-246, 1.1837996790456628e-245, 1.2084936832477347e-245, 9.792584422798124e-246, 6.468121041765527e-246, 3.56472526918372e-246]]}], {\"title\": \"Likelihood\", \"scene\": {\"xaxis\": {\"title\": \"mean\"}, \"yaxis\": {\"title\": \"sigma\"}, \"zaxis\": {\"title\": \"likelihood\"}}}, {\"showLink\": true, \"linkText\": \"Export to plot.ly\"})});</script>"
      ],
      "text/vnd.plotly.v1+html": [
       "<div id=\"cc3ffd98-82f6-466c-a8a5-0440824d0070\" style=\"height: 525px; width: 100%;\" class=\"plotly-graph-div\"></div><script type=\"text/javascript\">require([\"plotly\"], function(Plotly) { window.PLOTLYENV=window.PLOTLYENV || {};window.PLOTLYENV.BASE_URL=\"https://plot.ly\";Plotly.newPlot(\"cc3ffd98-82f6-466c-a8a5-0440824d0070\", [{\"type\": \"surface\", \"x\": [-0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4], \"y\": [3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4], \"z\": [[1.460825681631413e-249, 1.0373120935883478e-248, 4.682481547217632e-248, 1.419451486742629e-247, 3.03023899775314e-247, 4.7476206183476214e-247, 5.6591073457077994e-247, 5.296104474699062e-247, 4.000180909092284e-247, 2.4982970345070237e-247], [4.476409474673629e-249, 2.9894940695276932e-248, 1.275416701641574e-247, 3.670235624951869e-247, 7.467305457483984e-247, 1.1189830134533918e-246, 1.2798339607468397e-246, 1.1526214309907195e-246, 8.40009226031549e-247, 5.074219710609222e-247], [1.1650806670689444e-248, 7.383549620589285e-248, 3.0017958929467434e-247, 8.262581349916492e-247, 1.613409129439688e-246, 2.327470419846075e-246, 2.5697351244185605e-246, 2.2396342461337547e-246, 1.583114519199049e-246, 9.294570610141049e-247], [2.5755913849581866e-248, 1.5628307442373045e-247, 6.104684901647834e-247, 1.6195186743637765e-246, 3.0564596192454286e-246, 4.2722635540242456e-246, 4.580913843902494e-246, 3.8853287510968715e-246, 2.6777156728008803e-246, 1.5354096355898816e-246], [4.83607146694548e-248, 2.8349012866247653e-247, 1.0727496457009362e-246, 2.7637870756392706e-246, 5.076756973695156e-246, 6.920621412681512e-246, 7.25011540178241e-246, 6.017813259175593e-246, 4.06481179228702e-246, 2.2874629048579498e-246], [7.712643101726141e-248, 4.407001339848886e-247, 1.6288704370282895e-246, 4.1064979727080796e-246, 7.393458592634441e-246, 9.893394746319798e-246, 1.0187479491605617e-245, 8.321658792255388e-246, 5.537838366397335e-246, 3.0733974441126817e-246], [1.0447407591638561e-247, 5.87121076092514e-247, 2.137114232999409e-246, 5.312354564801568e-246, 9.440662565341089e-246, 1.2481269857029873e-245, 1.2709182264773008e-245, 1.0274042676853264e-245, 6.771173818490415e-246, 3.7240736313842286e-246], [1.2020114932797596e-247, 6.703342433119301e-247, 2.42282455330681e-246, 5.983439501146988e-246, 1.0569412905956128e-245, 1.3895858986836018e-245, 1.40765912878501e-245, 1.1324866440591987e-245, 7.430388161429455e-246, 4.0696087659692717e-246], [1.1746364697984904e-247, 6.558951005692979e-247, 2.3733906647249843e-246, 5.867636000603116e-246, 1.0375112463729928e-245, 1.3652909555156847e-245, 1.3842241650741372e-245, 1.1145128373283335e-245, 7.317839732550941e-246, 4.0107161668927876e-246], [9.749741907642282e-248, 5.499923069240653e-247, 2.0089518483418072e-246, 5.009839752459894e-246, 8.929524248058134e-246, 1.1837996790456628e-245, 1.2084936832477347e-245, 9.792584422798124e-246, 6.468121041765527e-246, 3.56472526918372e-246]]}], {\"title\": \"Likelihood\", \"scene\": {\"xaxis\": {\"title\": \"mean\"}, \"yaxis\": {\"title\": \"sigma\"}, \"zaxis\": {\"title\": \"likelihood\"}}}, {\"showLink\": true, \"linkText\": \"Export to plot.ly\"})});</script>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 进行似然函数的3D可视化\n",
    "data = [go.Surface(x=all_infer_mu, y=all_infer_sigma, z=mle)]\n",
    "layout = go.Layout(title=\"Likelihood\", scene={\"xaxis\": {'title': \"mean\"}, \"yaxis\": {\"title\": \"sigma\"},\"zaxis\": {\"title\": \"likelihood\"}})\n",
    "fig = go.Figure(data=data, layout=layout)\n",
    "iplot(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们看到上图中的最高点$\\mu \\approx 0, \\ \\sigma \\approx 4$, 是产生观测数据的真实概率分布的参数, 但是似然函数输出的值极小, 最高的值仅有$1*10^{-245}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对数似然值$(log \\ likelihood)$: \n",
    "我们对似然函数取$\\log$, 就得到了对数似然函数:   \n",
    "$$\\mathcal {L}(\\mu ,\\sigma \\mid X)=\\sum _{i=1}^{N}\\log P(x_{i}\\mid \\mu,\\sigma)$$\n",
    "为什么要对似然函数取对数?  \n",
    "1. 首先原本的似然函数是很多条件概率的乘积, 我们在找极大值的时候需要求似然函数的导数, 而乘积的导数不方便计算, 取对数可以吧乘除变成加减;\n",
    "2. 对似然函数取对数, 原本函数的极大值的位置没有改变;\n",
    "3. 如果观测到的数据点比较多, 原始似然函数的乘积可能非常接近于0, 甚至超出计算机的储存位数限制, 这样就全变成0了, 取对数可以把接近于$0$的数变成很大的负数, 也就是把原本似然函数的取值范围从$0$到$1$扩展到了$-\\infty$到$0$, 方便了计算.   \n",
    "   \n",
    "我们下面对$log \\ likelihood$进行3D可视化, 可看到原始的似然函数和对数似然函数的最高点是一样的位置:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_likelihood(observations, infer_mu, infer_sigma):\n",
    "    sum_ = 0\n",
    "    for obs in observations:\n",
    "        sum_ += stats.norm.logpdf(obs, loc=infer_mu, scale=infer_sigma)\n",
    "    return sum_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "mle = []\n",
    "for infer_mu in all_infer_mu:\n",
    "    temp_lis = []\n",
    "    for infer_sigma in all_infer_sigma:\n",
    "        temp_lis.append(log_likelihood(observations, infer_mu, infer_sigma))\n",
    "    mle.append(temp_lis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.plotly.v1+json": {
       "data": [
        {
         "type": "surface",
         "x": [
          -0.5,
          -0.4,
          -0.3,
          -0.2,
          -0.1,
          0,
          0.1,
          0.2,
          0.3,
          0.4
         ],
         "y": [
          3.5,
          3.6,
          3.7,
          3.8,
          3.9,
          4,
          4.1,
          4.2,
          4.3,
          4.4
         ],
         "z": [
          [
           -572.9646863442881,
           -571.0044702204175,
           -569.4972748479921,
           -568.3882474494818,
           -567.6298764759722,
           -567.180874399439,
           -567.0052518025246,
           -567.0715464239934,
           -567.3521783821593,
           -567.8229086559649
          ],
          [
           -571.8448668869746,
           -569.9459988969535,
           -568.4952450194975,
           -567.4382621065395,
           -566.727983750654,
           -566.3235126274345,
           -566.1892025251844,
           -566.2938940230815,
           -566.6102752803693,
           -567.1143452080269
          ],
          [
           -570.8883127357834,
           -569.041848561144,
           -567.6393072289874,
           -566.6267809187221,
           -565.9575834645206,
           -565.5911508554285,
           -565.4921300473679,
           -565.6296203069767,
           -565.9765387551076,
           -566.5090875452121
          ],
          [
           -570.095023890715,
           -568.292019212989,
           -566.9294614764602,
           -565.9538038860285,
           -565.318675617572,
           -564.9837890834228,
           -564.9140343690765,
           -565.0787252756792,
           -565.450968806374,
           -566.0071356675214
          ],
          [
           -569.4650003517685,
           -567.6965108524885,
           -566.3657077619174,
           -565.41933100846,
           -564.811260209808,
           -564.501427311417,
           -564.4549154903077,
           -564.6412089291887,
           -565.0335654341679,
           -565.6084895749547
          ],
          [
           -568.998242118945,
           -567.2553234796419,
           -565.9480460853588,
           -565.023362286016,
           -564.4353372412288,
           -564.1440655394114,
           -564.1147734110646,
           -564.3170712675054,
           -564.7243286384894,
           -565.3131492675119
          ],
          [
           -568.6947491922439,
           -566.9684570944495,
           -565.6764764467837,
           -564.7658977186967,
           -564.1909067118343,
           -563.9117037674058,
           -563.8936081313446,
           -564.1063122906297,
           -564.5232584193385,
           -565.1211147451934
          ],
          [
           -568.5545215716648,
           -566.8359116969119,
           -565.5509988461928,
           -564.6469373065023,
           -564.0779686216246,
           -563.8043419954001,
           -563.7914196511487,
           -564.0089319985611,
           -564.430354776717,
           -565.0323860079986
          ],
          [
           -568.5775592572089,
           -566.8576872870285,
           -565.5716132835861,
           -564.6664810494316,
           -564.0965229705998,
           -563.8219802233942,
           -563.8082079704772,
           -564.0249303912991,
           -564.4456177106218,
           -565.0469630559279
          ],
          [
           -568.763862248875,
           -567.0337838647996,
           -565.7383197589627,
           -564.8245289474871,
           -564.2465697587592,
           -563.9646184513889,
           -563.9439730893296,
           -564.1543074688458,
           -564.5690472210554,
           -565.164845888981
          ]
         ]
        }
       ],
       "layout": {
        "scene": {
         "xaxis": {
          "title": "mean"
         },
         "yaxis": {
          "title": "sigma"
         },
         "zaxis": {
          "title": "likelihood"
         }
        },
        "title": "Log Likelihood"
       }
      },
      "text/html": [
       "<div id=\"950446c6-e1c9-4748-b730-7eeaa58d9f4b\" style=\"height: 525px; width: 100%;\" class=\"plotly-graph-div\"></div><script type=\"text/javascript\">require([\"plotly\"], function(Plotly) { window.PLOTLYENV=window.PLOTLYENV || {};window.PLOTLYENV.BASE_URL=\"https://plot.ly\";Plotly.newPlot(\"950446c6-e1c9-4748-b730-7eeaa58d9f4b\", [{\"type\": \"surface\", \"x\": [-0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4], \"y\": [3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4], \"z\": [[-572.9646863442881, -571.0044702204175, -569.4972748479921, -568.3882474494818, -567.6298764759722, -567.180874399439, -567.0052518025246, -567.0715464239934, -567.3521783821593, -567.8229086559649], [-571.8448668869746, -569.9459988969535, -568.4952450194975, -567.4382621065395, -566.727983750654, -566.3235126274345, -566.1892025251844, -566.2938940230815, -566.6102752803693, -567.1143452080269], [-570.8883127357834, -569.041848561144, -567.6393072289874, -566.6267809187221, -565.9575834645206, -565.5911508554285, -565.4921300473679, -565.6296203069767, -565.9765387551076, -566.5090875452121], [-570.095023890715, -568.292019212989, -566.9294614764602, -565.9538038860285, -565.318675617572, -564.9837890834228, -564.9140343690765, -565.0787252756792, -565.450968806374, -566.0071356675214], [-569.4650003517685, -567.6965108524885, -566.3657077619174, -565.41933100846, -564.811260209808, -564.501427311417, -564.4549154903077, -564.6412089291887, -565.0335654341679, -565.6084895749547], [-568.998242118945, -567.2553234796419, -565.9480460853588, -565.023362286016, -564.4353372412288, -564.1440655394114, -564.1147734110646, -564.3170712675054, -564.7243286384894, -565.3131492675119], [-568.6947491922439, -566.9684570944495, -565.6764764467837, -564.7658977186967, -564.1909067118343, -563.9117037674058, -563.8936081313446, -564.1063122906297, -564.5232584193385, -565.1211147451934], [-568.5545215716648, -566.8359116969119, -565.5509988461928, -564.6469373065023, -564.0779686216246, -563.8043419954001, -563.7914196511487, -564.0089319985611, -564.430354776717, -565.0323860079986], [-568.5775592572089, -566.8576872870285, -565.5716132835861, -564.6664810494316, -564.0965229705998, -563.8219802233942, -563.8082079704772, -564.0249303912991, -564.4456177106218, -565.0469630559279], [-568.763862248875, -567.0337838647996, -565.7383197589627, -564.8245289474871, -564.2465697587592, -563.9646184513889, -563.9439730893296, -564.1543074688458, -564.5690472210554, -565.164845888981]]}], {\"title\": \"Log Likelihood\", \"scene\": {\"xaxis\": {\"title\": \"mean\"}, \"yaxis\": {\"title\": \"sigma\"}, \"zaxis\": {\"title\": \"likelihood\"}}}, {\"showLink\": true, \"linkText\": \"Export to plot.ly\"})});</script>"
      ],
      "text/vnd.plotly.v1+html": [
       "<div id=\"950446c6-e1c9-4748-b730-7eeaa58d9f4b\" style=\"height: 525px; width: 100%;\" class=\"plotly-graph-div\"></div><script type=\"text/javascript\">require([\"plotly\"], function(Plotly) { window.PLOTLYENV=window.PLOTLYENV || {};window.PLOTLYENV.BASE_URL=\"https://plot.ly\";Plotly.newPlot(\"950446c6-e1c9-4748-b730-7eeaa58d9f4b\", [{\"type\": \"surface\", \"x\": [-0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4], \"y\": [3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4], \"z\": [[-572.9646863442881, -571.0044702204175, -569.4972748479921, -568.3882474494818, -567.6298764759722, -567.180874399439, -567.0052518025246, -567.0715464239934, -567.3521783821593, -567.8229086559649], [-571.8448668869746, -569.9459988969535, -568.4952450194975, -567.4382621065395, -566.727983750654, -566.3235126274345, -566.1892025251844, -566.2938940230815, -566.6102752803693, -567.1143452080269], [-570.8883127357834, -569.041848561144, -567.6393072289874, -566.6267809187221, -565.9575834645206, -565.5911508554285, -565.4921300473679, -565.6296203069767, -565.9765387551076, -566.5090875452121], [-570.095023890715, -568.292019212989, -566.9294614764602, -565.9538038860285, -565.318675617572, -564.9837890834228, -564.9140343690765, -565.0787252756792, -565.450968806374, -566.0071356675214], [-569.4650003517685, -567.6965108524885, -566.3657077619174, -565.41933100846, -564.811260209808, -564.501427311417, -564.4549154903077, -564.6412089291887, -565.0335654341679, -565.6084895749547], [-568.998242118945, -567.2553234796419, -565.9480460853588, -565.023362286016, -564.4353372412288, -564.1440655394114, -564.1147734110646, -564.3170712675054, -564.7243286384894, -565.3131492675119], [-568.6947491922439, -566.9684570944495, -565.6764764467837, -564.7658977186967, -564.1909067118343, -563.9117037674058, -563.8936081313446, -564.1063122906297, -564.5232584193385, -565.1211147451934], [-568.5545215716648, -566.8359116969119, -565.5509988461928, -564.6469373065023, -564.0779686216246, -563.8043419954001, -563.7914196511487, -564.0089319985611, -564.430354776717, -565.0323860079986], [-568.5775592572089, -566.8576872870285, -565.5716132835861, -564.6664810494316, -564.0965229705998, -563.8219802233942, -563.8082079704772, -564.0249303912991, -564.4456177106218, -565.0469630559279], [-568.763862248875, -567.0337838647996, -565.7383197589627, -564.8245289474871, -564.2465697587592, -563.9646184513889, -563.9439730893296, -564.1543074688458, -564.5690472210554, -565.164845888981]]}], {\"title\": \"Log Likelihood\", \"scene\": {\"xaxis\": {\"title\": \"mean\"}, \"yaxis\": {\"title\": \"sigma\"}, \"zaxis\": {\"title\": \"likelihood\"}}}, {\"showLink\": true, \"linkText\": \"Export to plot.ly\"})});</script>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = [go.Surface(x=all_infer_mu, y=all_infer_sigma, z=mle)]\n",
    "layout = go.Layout(title=\"Log Likelihood\", scene={\"xaxis\": {'title': \"mean\"}, \"yaxis\": {\"title\": \"sigma\"},\"zaxis\": {\"title\": \"likelihood\"}})\n",
    "fig = go.Figure(data=data, layout=layout)\n",
    "iplot(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}