Cartman0/histogramからKL.ipynb

## histogramからKL.ipynb
{
  "cells": [
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "# histogramからKL"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "参考：https://paper.hatenadiary.jp/entry/2018/10/13/164539\n\n厳密な確率密度関数を知らなくても密度化したヒストグラムから\nKL情報量を近似して求められるのでは\n\n"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## 離散分布のKL情報量"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "離散分布のKL情報量\n\n$$\nD_{KL}(p||q) = \\sum_{i} p(i) \\log{(\\frac{p(i)}{q(i)})}\n$$"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "\np, qをヒストグラム密度にする．\n\n"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## histogramの密度化"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "密度化するとビン幅を考慮して合計が1になる．\nこれでヒストグラムを生成するデータ間のデータ数（標本の大きさ）が異なっても考慮できる．\n\nbinの数をnとして，bin幅のベクトルを定義する．\n\n$$\n\\vec{b} = (bin_1 - bin_0, \\cdots, bin_{n}-bin_{n-1})\n=(b_1 , \\cdots, b_{n})\n$$\n\n合計して1より以下を満たす．\n\n$$\n\\sum_i b_i p(i) =1\n$$\n"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "頻度を，bin幅を考慮して総頻度で割れば密度化できるので，\nbinのi番目の密度は，\n\n$$\np(b_i) \n= \\frac{\\text{counts}(b_i)}{\\sum{\\text{counts}(b_i)* b_i}}\n= \\frac{\\text{counts}(b_i)}{\\sum{\\text{counts}(b_i)* (bin_{i} - bin_{i-1})}}\n$$"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## 実装"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "参考：https://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram.html"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "### 標準正規分布のヒストグラムからKLを計算"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "import scipy as sp\nimport matplotlib.pyplot as plt\nimport scipy.stats as stats\n\nsample1 = stats.norm(loc=0, scale=1).rvs(100)\nsample2 = stats.norm(loc=0, scale=1).rvs(100)\n\nbins = sp.linspace(-4.5, 4.5, 25+1)\nhist1, bins = sp.histogram(sample1, bins, density=True)\nhist1\nhist2, bins = sp.histogram(sample2, bins, density=True)\nhist2\n\ndef create_hists(data1, data2, bins, eps=1e-9):\n    hist1, bins = sp.histogram(sample1, bins, density=True)\n    hist2, bins = sp.histogram(sample2, bins, density=True)\n    hist1 = (hist1 + eps)\n    hist1 = hist1 / (sp.diff(bins) * hist1.sum())\n    hist2 = (hist2 + eps)\n    hist2 = hist2 / (sp.diff(bins) * hist2.sum())\n    return hist1, hist2",
      "execution_count": 15,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# 0が無くなるように再構成\n\nhist1, hist2 = create_hists(sample1, sample2, bins)\nhist1",
      "execution_count": 16,
      "outputs": [
        {
          "data": {
            "text/plain": "array([9.99999991e-10, 9.99999991e-10, 9.99999991e-10, 9.99999991e-10,\n       9.99999991e-10, 2.77777785e-02, 2.77777785e-02, 1.38888889e-01,\n       8.33333336e-02, 8.33333336e-02, 4.16666664e-01, 5.27777774e-01,\n       4.44444441e-01, 3.61111109e-01, 2.77777776e-01, 1.66666666e-01,\n       1.11111111e-01, 5.55555561e-02, 2.77777785e-02, 9.99999991e-10,\n       9.99999991e-10, 2.77777785e-02, 9.99999991e-10, 9.99999991e-10,\n       9.99999991e-10])"
          },
          "execution_count": 16,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "scrolled": true,
        "trusted": true
      },
      "cell_type": "code",
      "source": "hist2.sum()",
      "execution_count": 17,
      "outputs": [
        {
          "data": {
            "text/plain": "2.7777777777777772"
          },
          "execution_count": 17,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "scrolled": true,
        "trusted": true
      },
      "cell_type": "code",
      "source": "%matplotlib inline\n\nplt.bar(bins[:-1], hist1, label=\"sample 1\")\nplt.bar(bins[:-1], hist2, label=\"sample 2\")\nplt.legend()\nplt.grid(True)\nplt.show()",
      "execution_count": 18,
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFKRJREFUeJzt3X2QXVWZ7/Hvk9e+mBiqwMoAnZqOmqQSEwRpXiws7aCBYKaSKgemiGBBzWhAJ4zjAApXC5PcWIp6cW6VqJPrRVJ3RoIyWjdiLjgoXb6UYF6Gl4RMMEW1kxanBuIUY4MJRp/7RzrcTqeT3qdzXrpXfz9/nb33Oms/a6f7l9X7nLNOZCaSpLJMaHUBkqT6M9wlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBZrUqhOffvrp2dHR0arT1+yll17iNa95TavLaKnxfg3G+/jBazAaxr99+/YXMvN1w7VrWbh3dHSwbdu2Vp2+Zt3d3XR1dbW6jJYa79dgvI8fvAajYfwR8Ysq7bwtI0kFMtwlqUCGuyQVqGX33CWV63e/+x29vb0cOHCg1aXU1YwZM9i9e3dTztXW1kZ7ezuTJ08e0fMNd0l119vby/Tp0+no6CAiWl1O3fzmN79h+vTpDT9PZrJ//356e3uZPXv2iPrwtoykujtw4ACnnXZaUcHeTBHBaaeddlJ/+RjukhrCYD85J3v9DHdJKpD33CU1XMet361rfz2fWVbX/k5WV1cXn//85+ns7KzU/pvf/CZr1qxh9+7d/OxnP6v8vFoY7hpz6h0UxzPaAkTlWLhwId/61re4/vrrG3YOb8tIKs5LL73EsmXLePOb38zChQu57777AFi3bh3nn38+CxcuZNWqVWQmcHjm/ZGPfIS3v/3tzJ8/n61bt/Ke97yHOXPm8IlPfAKAnp4ezjvvPK699lrOPvtsrrjiCl5++eVjzv29732Pt771rbzlLW/hyiuvpK+v75g28+fPZ968eQ28Aoa7pAI9+OCDnHnmmTzxxBPs3LmTpUuXArB69Wq2bt3Kzp07+e1vf8sDDzzw6nOmTJnCD3/4Q2644QZWrFjBXXfdxc6dO7nnnnvYv38/AD//+c9ZtWoVTz75JK997Wv50pe+dNR5X3jhBdavX8/DDz/Mjh076Ozs5M4772zewAcw3CUVZ9GiRTz88MN87GMf40c/+hEzZswA4JFHHuHCCy9k0aJF/OAHP2DXrl2vPmf58uWvPvdNb3oTZ5xxBlOnTuX1r389+/btA6C9vZ2LL74YgGuuuYYf//jHR5330Ucf5emnn+biiy/mnHPOYePGjfziF5XW+ao777lLKs7cuXPZvn07W7Zs4bbbbuPSSy/lox/9KB/60IfYtm0bs2bNYs2aNUe9j3zq1KkATJgw4dXHR7YPHToEHPv2xMHbmcmSJUu49957GzW0ypy5SyrOc889xymnnMI111zDzTffzI4dO14N8tNPP52+vj7uv//+mvvdt28fP/3pTwG49957edvb3nbU8Ysuuoif/OQn7N27F4CXX36ZZ5555iRHMzLO3CU1XLPfefTUU09xyy23MGHCBCZPnsyXv/xlTj31VD7wgQ+waNEiOjo6OP/882vud968eWzcuJHrr7+eOXPm8MEPfvCo46973eu45557WLlyJQcPHgRg/fr1zJ0796h23/72t7nxxht5/vnnWbZsGeeccw4PPfTQyAc8hDjyavEJG0UsBf4HMBH4amZ+ZtDx64DPAb/s3/XFzPzqifrs7OxMv6xjbBkt16BVb4UcLeNvparXYPfu3cyfP7/xBTVRT08P7373u3n66aebds6hrmNEbM/MYd8YP+zMPSImAncBS4BeYGtEbM7MwSO8LzNXVy9bktQoVe65XwDszcxnM/MVYBOworFlSdLo0tHRwWOPPdbqMiqrEu5nAfsGbPf27xvsTyPiyYi4PyJm1aU6SdKIDHvPPSKuBC7LzPf3b78PuCAzbxzQ5jSgLzMPRsQNwJ9l5iVD9LUKWAUwc+bM8zZt2lS/kTRYX18f06ZNa3UZLTVarsFTv3yxKedZdNaMo7ZHy/hbqeo1mDFjBm984xubUFFz/f73v2fixIlNO9/evXt58cWjf94XL15cn3vuHJ6pD5yJtwPPDWyQmfsHbP5P4I6hOsrMDcAGOPyC6lh6ccoX00bPNbiuWS+oXt111PZoGX8r1fKCajO+1KLZmvVlHUe0tbVx7rnnjui5VW7LbAXmRMTsiJgCXAVsHtggIs4YsLkcaM73UEmShjTszD0zD0XEauAhDr8V8u7M3BUR64BtmbkZ+KuIWA4cAn4NXNfAmiWNNWtmDN+mpv6ac2uuqlqX/L3lllv4zne+w5QpU3jDG97A1772NU499dS61lTpE6qZuSUz52bmGzLzU/37bu8PdjLztsx8U2a+OTMXZ+a/1LVKSSrIkiVL2LlzJ08++SRz587l05/+dN3P4fIDkooz2pf8vfTSS5k06fCNk4suuoje3t66XwPDXVJxxtKSv3fffTeXX355na+A4S6pQGNlyd9PfepTTJo0iauvvrqu4wcXDpNUoLGw5O/GjRt54IEH+P73v39MP/XgzF1ScUb7kr8PPvggd9xxB5s3b+aUU06puY4qnLlLarwmv3VxtC/5u3r1ag4ePMiSJUuAw/8pfOUrXxnhaIdmuEsqzmWXXcZll112zP7169ezfv36Y/Z3d3e/+rirq+uoT+EeOdbT08OECROGDOGBz7/kkkvYunXrCes7MrNvJG/LSFKBDHdJqqDEJX8lqWZVvuVNx3ey189wl1R3bW1t7N+/34Afocxk//79tLW1jbgPX1CVVHft7e309vby/PPPt7qUujpw4MBJBW4t2traaG9vH/HzDXdJdTd58mRmz57d6jLqrru7e8Trqzebt2UkqUCGuyQVyHCXpAJ5z10aoKftvf9/Y82gg/PWwpoV9TnRKPsmIZXHmbskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekArn8gFpvzYyamvc0ZzltaUxz5i5JBTLcJalAlcI9IpZGxJ6I2BsRt56g3RURkRHRWb8SJUm1GjbcI2IicBdwObAAWBkRC4ZoNx34K+CxehcpSapNlZn7BcDezHw2M18BNgFDLWr934DPAgfqWJ8kaQSqhPtZwL4B2739+14VEecCszLzgTrWJkkaoSpvhYwh9uWrByMmAF8Arhu2o4hVwCqAmTNn0t3dXanI0aCvr29M1dsIDbsG89bWv88G6Jt6Jt31qnWM/iyN99+DsTT+KuHeC8wasN0OPDdgezqwEOiOCIA/AjZHxPLM3Dawo8zcAGwA6OzszK6urpFX3mTd3d2MpXoboWHXoF5fXddg3fPW0rXnk/XpbOXY/Jq98f57MJbGX+W2zFZgTkTMjogpwFXA5iMHM/PFzDw9MzsyswN4FDgm2CVJzTNsuGfmIWA18BCwG/hGZu6KiHURsbzRBUqSaldp+YHM3AJsGbTv9uO07Tr5siRJJ8NPqEpSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCVQr3iFgaEXsiYm9E3DrE8Rsi4qmIeDwifhwRC+pfqiSpqmHDPSImAncBlwMLgJVDhPfXM3NRZp4DfBa4s+6VSpIqqzJzvwDYm5nPZuYrwCZgxcAGmfmfAzZfA2T9SpQk1WpShTZnAfsGbPcCFw5uFBF/CfwNMAW4pC7VSZJGJDJPPMmOiCuByzLz/f3b7wMuyMwbj9P+vf3trx3i2CpgFcDMmTPP27Rp00mW3zx9fX1Mmzat1WW0VMOuwa8er3+fDdA39UymHXyuPp2dcU59+mmy8f57MBrGv3jx4u2Z2Tlcuyoz915g1oDtduBEP+GbgC8PdSAzNwAbADo7O7Orq6vC6UeH7u5uxlK9jdCwa7BmxfBtRoHueWvp2vPJ+nS28sX69NNk4/33YCyNv8o9963AnIiYHRFTgKuAzQMbRMScAZvLgJ/Xr0RJUq2Gnbln5qGIWA08BEwE7s7MXRGxDtiWmZuB1RHxLuB3wH8Ax9ySkSQ1T5XbMmTmFmDLoH23D3j84TrXJUk6CX5CVZIKZLhLUoEMd0kqkOEuSQWq9IKqpDpbM6OBfY/N99Crvpy5S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBVoUqsLUKHWzGh1BdK45sxdkgpkuEtSgQx3SSqQ4S5JBaoU7hGxNCL2RMTeiLh1iON/ExFPR8STEfH9iPjj+pcqSapq2HCPiInAXcDlwAJgZUQsGNTsn4HOzDwbuB/4bL0LlSRVV2XmfgGwNzOfzcxXgE3AioENMvORzHy5f/NRoL2+ZUqSahGZeeIGEVcASzPz/f3b7wMuzMzVx2n/ReDfMnP9EMdWAasAZs6ced6mTZtOsvzm6evrY9q0aa0uo6Vquga/eryxxbRA39QzmXbwuVaXMbwzzmlY1+P992A0jH/x4sXbM7NzuHZVPsQUQ+wb8n+EiLgG6ATeMdTxzNwAbADo7OzMrq6uCqcfHbq7uxlL9TZCTddgzYrh24wx3fPW0rXnk60uY3grX2xY1+P992Asjb9KuPcCswZstwPHTF8i4l3Ax4F3ZObB+pQnSRqJKvfctwJzImJ2REwBrgI2D2wQEecCfwcsz8x/r3+ZkqRaDBvumXkIWA08BOwGvpGZuyJiXUQs72/2OWAa8M2IeDwiNh+nO0lSE1RaOCwztwBbBu27fcDjd9W5LknSSfATqpJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SClTpa/ZUpo5bv1tT+5sWHeK6is/paRtJRZLqxZm7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAK5/IBUmFqXlej5zLIGVaJWqjRzj4ilEbEnIvZGxK1DHH97ROyIiEMRcUX9y5Qk1WLYcI+IicBdwOXAAmBlRCwY1OxfgeuAr9e7QElS7arclrkA2JuZzwJExCZgBfD0kQaZ2dN/7A8NqFGSVKMqt2XOAvYN2O7t3ydJGqWqzNxjiH05kpNFxCpgFcDMmTPp7u4eSTct0dfXN6bqreKmRYdqaj/zv1R/TveEtSMpaVTrm3om3fNG/7hu+kNt/661/FyX+HtQi7E0/irh3gvMGrDdDjw3kpNl5gZgA0BnZ2d2dXWNpJuW6O7uZizVW0XVL9444qZFh/jvT1V7g1VP2ydHUtKo1j1vLV17Rv+4rjtQ20tfPVd3VW5b4u9BLcbS+KvcltkKzImI2RExBbgK2NzYsiRJJ2PYcM/MQ8Bq4CFgN/CNzNwVEesiYjlARJwfEb3AlcDfRcSuRhYtSTqxSn9jZ+YWYMugfbcPeLyVw7drJEmjgMsPSFKBXH5gHOtpe29N7bsnrC3yhdLS1Prvypoa2s5bC2tW1ND3i7XVorpx5i5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUoEmtLkDH13Hrd2tqX/O33kuNtmZGTc07Dny9QYVAz2eWNazv0ciZuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSpQpXCPiKURsSci9kbErUMcnxoR9/UffywiOupdqCSpumHDPSImAncBlwMLgJURsWBQs78A/iMz3wh8Abij3oVKkqqrMnO/ANibmc9m5ivAJmDFoDYrgI39j+8H3hkRUb8yJUm1qBLuZwH7Bmz39u8bsk1mHgJeBE6rR4GSpNpVWX5gqBl4jqANEbEKWNW/2RcReyqcf7Q4HXih1UWcSOP/VPrIqL8GjTXexw+NvwZ/0rCeoz43i0fDz8AfV2lUJdx7gVkDttuB547TpjciJgEzgF8P7igzNwAbqhQ22kTEtszsbHUdrTTer8F4Hz94DcbS+KvcltkKzImI2RExBbgK2DyozWbg2v7HVwA/yMxjZu6SpOYYduaemYciYjXwEDARuDszd0XEOmBbZm4G/hfwvyNiL4dn7Fc1smhJ0olVWvI3M7cAWwbtu33A4wPAlfUtbdQZk7eT6my8X4PxPn7wGoyZ8Yd3TySpPC4/IEkFMtxHICJujoiMiNNbXUuzRcTnIuJfIuLJiPh2RJza6pqaYbglOEoWEbMi4pGI2B0RuyLiw62uqVUiYmJE/HNEPNDqWoZjuNcoImYBS4B/bXUtLfJPwMLMPBt4BritxfU0XMUlOEp2CLgpM+cDFwF/Oc7GP9CHgd2tLqIKw712XwA+yhAf0hoPMvN7/Z9CBniUw597KF2VJTiKlZm/yswd/Y9/w+FwG/wp9eJFRDuwDPhqq2upwnCvQUQsB36ZmU+0upZR4s+B/9vqIpqgyhIc40L/iq/nAo+1tpKW+FsOT+z+0OpCqqj0VsjxJCIeBv5oiEMfB/4rcGlzK2q+E12DzPw//W0+zuE/1/+hmbW1SKXlNUoXEdOAfwT+OjP/s9X1NFNE/Anw75m5PSK6Wl1PFYb7IJn5rqH2R8QiYDbwRP+Cl+3Ajoi4IDP/rYklNtzxrsEREXEthxcBeec4+SRylSU4ihYRkzkc7P+Qmd9qdT0tcDGwPCLeDbQBr42Iv8/Ma1pc13H5PvcRiogeoDMzW72IUFNFxFLgTuAdmfl8q+tphv71kp4B3gn8ksNLcrw3M3e1tLAm6V++eyPw68z861bX02r9M/ebM7Nxq5zVgffcVasvAtOBf4qIxyPiK60uqNH6X0A+sgTHbuAb4yXY+10MvA+4pP/f/PH+GaxGMWfuklQgZ+6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAv0/JdmWl4Uv5uUAAAAASUVORK5CYII=\n",
            "text/plain": "<Figure size 432x288 with 1 Axes>"
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "def D_KL(p, q):\n    return p.dot(sp.log(p/q))\n    \nD_KL(hist1, hist2)    ",
      "execution_count": 19,
      "outputs": [
        {
          "data": {
            "text/plain": "1.1641780255612815"
          },
          "execution_count": 19,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# sample 1000\nimport scipy as sp\nimport matplotlib.pyplot as plt\nimport scipy.stats as stats\n\nN = 1000\nsample1 = stats.norm(loc=0, scale=1).rvs(N)\nsample2 = stats.norm(loc=0, scale=1).rvs(N)\n\nhist1, hist2 = create_hists(sample1, sample2, bins, eps=1e-9)\n\nplt.bar(bins[:-1], hist1, label=\"sample 1\")\nplt.bar(bins[:-1], hist2, label=\"sample 2\")\nplt.legend()\nplt.grid(True)\nplt.savefig(\"hist_1000.png\")\nplt.show()\n\nD_KL(hist1, hist2) ",
      "execution_count": 22,
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAG9NJREFUeJzt3X9wXWW97/H3p4UmA8WClMmxtJrCSWtLq0VDwcGLAWkpp9jMOHBtFQaVY0CpeuCAlqvTX5YR1EEZLWLHU+Gee6QC/rgRKhWEff1xBNJigf641VAjjfWM2DJIgBZTv/ePvcrdDWmzdrJ/Jevzmumw11rPs9b3WWR/8+RZaz1LEYGZmWXDqGoHYGZmleOkb2aWIU76ZmYZ4qRvZpYhTvpmZhnipG9mliFO+mZmGeKkb2aWIU76ZmYZclS1A+hr/Pjx0djYWO0wivLSSy9x7LHHVjuMqsl6+8HnIOvth+qfg02bNv0lIk4aqFzNJf3GxkY2btxY7TCKksvlaGlpqXYYVZP19oPPQdbbD9U/B5L+kKach3fMzDLESd/MLEOc9M3MMqTmxvTNbGT729/+Rnd3N/v27at2KCU1btw4tm/fXvbj1NfXM3HiRI4++uhB1XfSN7OK6u7u5rjjjqOxsRFJ1Q6nZF588UWOO+64sh4jItizZw/d3d1Mnjx5UPvw8I6ZVdS+ffs48cQTR1TCrxRJnHjiiUP6K8lJ38wqzgl/8IZ67pz0zcwyxGP6ZlZVjUvuL+n+um6aX9L9DUVLSwtf+cpXaG5uTlX+nnvuYfny5Wzfvp3HH388db1iOOlbdi0fV7p9TV0By1sPc5wXSnccG9FmzJjBD37wA6688sqyHcPDO2aWGS+99BLz58/n7W9/OzNmzOB73/seACtXruSMM85gxowZtLW1ERFAvqd+zTXXcM455zBt2jQ6Ojp4//vfT1NTE5///OcB6Orq4q1vfStXXnklb3vb27j44ot5+eWXX3fsn/70p7zrXe/iHe94B5dccgk9PT2vKzNt2jSmTp1axjPgnr6NJKXsuduI9MADDzBhwgTuvz8/pPTCC/m/whYvXszSpUsBuOyyy7jvvvt43/veB8CYMWP4+c9/zq233kprayubNm3ijW98I6eeeirXXHMNADt27ODrX/86c+bM4aMf/Si33XYb11133WvH/ctf/sKqVat46KGHOPbYY7n55pu55ZZbXjtmJbmnb2aZMXPmTB566CE++9nP8otf/IJx4/IdhUceeYQzzzyTmTNn8vDDD7N169bX6ixYsOC1uqeddhpvetObqKur45RTTmHXrl0ATJo0ibPOOguASy+9lF/+8peHHPfRRx9l27ZtnH322cyaNYs777yTP/wh1fxoJeeevpllxpQpU9i0aRPr16/nhhtuYO7cuXzmM5/hE5/4BBs3bmTSpEksX778kPvg6+rqABg1atRrnw8u9/b2Aq+/jbLvckQwZ84c7rrrrnI1LTX39M0sM3bv3s0xxxzDpZdeynXXXccTTzzxWoIfP348PT093HvvvUXv99lnn+Wxxx4D4K677uLd7373IdvPOussfvWrX9HZ2QnAyy+/zG9/+9shtmZw3NM3s6qq5C2WTz/9NNdffz2jRo3i6KOP5pvf/CbHH388H/vYx5g5cyaNjY2cccYZRe932rRp3HXXXVx77bU0NTXx8Y9//JDtJ510EnfccQeLFi1i//79AKxatYopU6YcUu6HP/whn/zkJ3nuueeYP38+s2bNYsOGDYNvcD908Cr1EQtJ84BbgdHAtyPipsOUuxi4BzgjIjYm624ArgAOAJ+KiCO2oLm5OfwSleGlZtpfxQu5uakraNmxrP+NGbhls5ifge3btzNt2rTyBlRBXV1dXHTRRfz6178u+9w7B/V3DiVtiogBb+wfsKcvaTSwGpgDdAMdktojYlufcscBnwIeK1g3HVgInAZMAB6SNCUiDgzYKjMzK7k0Y/qzgc6I2BkRrwLrgP6eQvkC8CWgcCagVmBdROyPiN8Dncn+zMxGhMbGRrZs2VLtMFJLk/RPBnYVLHcn614j6XRgUkTcV2xdMzOrnDQXcvub0u21CwGSRgFfBT5cbN2CfbQBbQANDQ3kcrkUYdWOnp6eYRdzKdVM+6euqNqhe+omkDvc8Wvh3JRZMT8D48aN48UXXyxvQFVw4MCBirVr3759g/7OpUn63cCkguWJwO6C5eOAGUAuuTf1H4B2SQtS1AUgItYAayB/IbcmLgoWoWYuZFZJzbT/cHPfVMARL+Qu8oXcQtu3b6/YBc9KqsRLVA6qr6/n9NNPH1TdNMM7HUCTpMmSxpC/MNt+cGNEvBAR4yOiMSIagUeBBcndO+3AQkl1kiYDTcDjg4rUzMyGbMCefkT0SloMbCB/y+baiNgqaSWwMSLaj1B3q6S7gW1AL3C179wxs0OU+lbbGrpFttipla+//np+/OMfM2bMGE499VS+853vcPzxx5c0plRP5EbE+oiYEhGnRsSNybql/SX8iGg5eI9+snxjUm9qRPykdKGbmY0sc+bMYcuWLTz11FNMmTKFL37xiyU/hqdhMLPMqPWplefOnctRR+UHYM466yy6u7tLfg6c9M0sMw5Orfzkk0+yZcsW5s2bB+SnVu7o6GDLli288sor3Hff/7/7/ODUyldddRWtra2sXr2aLVu2cMcdd7Bnzx4gP7XyRz7yEZ566ine8IY3cNtttx1y3MKplZ944gmam5u55ZZbjhjr2rVrufDCC0t8Bpz0zSxDhsvUyjfeeCNHHXUUH/rQh0rafvCEa2aWIcNhauU777yT++67j5/97Gev208puKdvZplR61MrP/DAA9x88820t7dzzDHHFB1HGu7pm1l1VfAWy1qfWnnx4sXs37+fOXPmAPlfFrfffvsgW9s/J30zy4wLLriACy644HXrV61axapVq163vnCqg5aWlkOeOj64rauri1GjRvG1r33tdU/kFtY/77zz6OjoOGJ8B/8SKCcP75iZZYiTvpnZEIzEqZXNzEoqzRv7rH9DPXdO+mZWUfX19ezZs8eJfxAigj179lBfXz/offhCrplV1MSJE+nu7ua5556rdigltW/fviEl47Tq6+uZOHHioOs76ZtZRR199NFMnjy52mGUXC6XG/Qc95Xk4R0zswxxT9+s3IqZL76G5oK3kSlVT1/SPEk7JHVKWtLP9qskPS1ps6RfSpqerG+U9EqyfrOk0j5aZmZmRRmwpy9pNLAamEP+nbcdktojYltBse9GxO1J+QXALcC8ZNszETGrtGGbmdlgpOnpzwY6I2JnRLwKrAMOeQN1RPy1YPFYwPdimZnVoDRJ/2RgV8Fyd7LuEJKulvQM8CXgUwWbJkv6jaT/I+m/DSlaMzMbEg30gISkS4ALIuKfk+XLgNkR8cnDlP9gUv5ySXXA2IjYI+mdwI+A0/r8ZYCkNqANoKGh4Z3r1q0barsqqqenh7Fjx1Y7jKqpmfb/aXPVDt1TN4Gx+3cPfUdvGp4joTXzM1BF1T4H55577qaIGPAN7Gnu3ukGJhUsTwSO9NO9DvgmQETsB/YnnzclfwlMATYWVoiINcAagObm5iicyW44yOVyDLeYS6lm2r+8deAyZZKbuoKWHcuGvqNFw/PunZr5Gaii4XIO0gzvdABNkiZLGgMsBNoLC0hqKlicD/wuWX9SciEYSacATcDOUgRuZmbFG7CnHxG9khYDG4DRwNqI2CppJbAxItqBxZLOB/4GPA9cnlQ/B1gpqRc4AFwVEXvL0RAzMxtYqoezImI9sL7PuqUFnz99mHrfB74/lADNzKx0PA2DmVmGOOmbmWWIk76ZWYY46ZuZZYiTvplZhjjpm5lliJO+mVmGOOmbmWWIk76ZWYY46ZuZZYjfkWtWQxqX3F9U+a6b5pcpEhup3NM3M8sQJ30zswxx0jczyxCP6VtNK2aMu6u+jIFUSFf9B4usMTzftGXV456+mVmGpOrpS5oH3Er+zVnfjoib+my/Cria/NuxeoC2iNiWbLsBuCLZ9qmI2FC68G2kK77na2ZHMmBPP3nH7WrgQmA6sEjS9D7FvhsRMyNiFvAl4Jak7nTy79Q9DZgH3HbwnblmZlZ5aYZ3ZgOdEbEzIl4F1gGthQUi4q8Fi8cCkXxuBdZFxP6I+D3QmezPzMyqIM3wzsnAroLlbuDMvoUkXQ1cC4wBziuo+2ifuicPKlIzMxuyNElf/ayL162IWA2slvRB4PPA5WnrSmoD2gAaGhrI5XIpwqodPT09wy7mUipr+6euKM9+S6ynbgK5asRaIz93Wf8OwPA5B2mSfjcwqWB5IrD7COXXAd8spm5ErAHWADQ3N0dLS0uKsGpHLpdjuMVcSmVt//LWgcvUgNzUFbTsWFb5Ay+qjVs2s/4dgOFzDtKM6XcATZImSxpD/sJse2EBSU0Fi/OB3yWf24GFkuokTQaagMeHHraZmQ3GgD39iOiVtBjYQP6WzbURsVXSSmBjRLQDiyWdD/wNeJ780A5JubuBbUAvcHVEHChTW8zMbACp7tOPiPXA+j7rlhZ8/vQR6t4I3DjYAM3MrHT8RK6ZWYY46ZuZZYiTvplZhjjpm5lliJO+mVmGOOmbmWWIk76ZWYY46ZuZZYiTvplZhjjpm5lliJO+mVmGOOmbmWWIk76ZWYY46ZuZZYiTvplZhjjpm5llSKqkL2mepB2SOiUt6Wf7tZK2SXpK0s8kvaVg2wFJm5N/7X3rmplZ5Qz45ixJo4HVwBzyLzrvkNQeEdsKiv0GaI6IlyV9HPgS8IFk2ysRMavEcZuZ2SCk6enPBjojYmdEvAqsA1oLC0TEIxHxcrL4KDCxtGGamVkpKCKOXEC6GJgXEf+cLF8GnBkRiw9T/hvAf0XEqmS5F9hM/sXoN0XEj/qp0wa0ATQ0NLxz3bp1g29RFfT09DB27Nhqh1E1ZW3/nzaXZ78l1lM3gbH7d1f+wG+qjT+is/4dgOqfg3PPPXdTRDQPVC7Ni9HVz7p+f1NIuhRoBt5TsPrNEbFb0inAw5KejohnDtlZxBpgDUBzc3O0tLSkCKt25HI5hlvMpVTW9i9vHbhMDchNXUHLjmWVP/CiFyp/zH5k/TsAw+ccpBne6QYmFSxPBF7XpZF0PvA5YEFE7D+4PiJ2J//dCeSA04cQr5mZDUGapN8BNEmaLGkMsBA45C4cSacD3yKf8P9csP4ESXXJ5/HA2UDhBWAzM6ugAYd3IqJX0mJgAzAaWBsRWyWtBDZGRDvwZWAscI8kgGcjYgEwDfiWpL+T/wVzU5+7fszMrILSjOkTEeuB9X3WLS34fP5h6v0nMHMoAZqZWemkSvpmVpsal9xfVPmum+aXKRIbLjwNg5lZhjjpm5lliJO+mVmGOOmbmWWIk76ZWYY46ZuZZYiTvplZhjjpm5lliJO+mVmGOOmbmWWIk76ZWYY46ZuZZYiTvplZhjjpm5llSKqkL2mepB2SOiUt6Wf7tZK2SXpK0s8kvaVg2+WSfpf8u7yUwZuZWXEGTPqSRgOrgQuB6cAiSdP7FPsN0BwRbwPuBb6U1H0jsAw4E5gNLJN0QunCNzOzYqTp6c8GOiNiZ0S8CqwDWgsLRMQjEfFysvgo+ZenA1wAPBgReyPieeBBYF5pQjczs2KlSfonA7sKlruTdYdzBfCTQdY1M7MySvO6RPWzLvotKF0KNAPvKaaupDagDaChoYFcLpcirNrR09Mz7GIupaLa/6fNxe186oqi46mGnroJ5KoQ6x38rqjy5fo5zfp3AIbPOUiT9LuBSQXLE4HdfQtJOh/4HPCeiNhfULelT91c37oRsQZYA9Dc3BwtLS19i9S0XC7HcIu5lIpq//LWgcsMQ7mpK2jZsazaYQxs0Qtl2W3WvwMwfM5BmuGdDqBJ0mRJY4CFQHthAUmnA98CFkTEnws2bQDmSjohuYA7N1lnZmZVMGBPPyJ6JS0mn6xHA2sjYquklcDGiGgHvgyMBe6RBPBsRCyIiL2SvkD+FwfAyojYW5aWmJnZgNIM7xAR64H1fdYtLfh8/hHqrgXWDjZAMzMrHT+Ra2aWIU76ZmYZ4qRvZpYhTvpmZhnipG9mliFO+mZmGeKkb2aWIU76ZmYZ4qRvZpYhTvpmZhnipG9mliGp5t4xsxFi+bgiypZnGmarLvf0zcwyxEnfzCxDnPTNzDLESd/MLENSJX1J8yTtkNQpaUk/28+R9ISkXkkX99l2QNLm5F9737pmZlY5A969I2k0sBqYQ/5F5x2S2iNiW0GxZ4EPA9f1s4tXImJWCWI1M7MhSnPL5mygMyJ2AkhaB7QCryX9iOhKtv29DDGamVmJpEn6JwO7Cpa7gTOLOEa9pI1AL3BTRPyobwFJbUAbQENDA7lcrojdV19PT8+wi7mUimr/1BVljaVaeuomkBtpbSviZzrr3wEYPucgTdJXP+uiiGO8OSJ2SzoFeFjS0xHxzCE7i1gDrAFobm6OlpaWInZffblcjuEWcykV1f7lrWWNpVpyU1fQsmNZtcMorUXpH87K+ncAhs85SHMhtxuYVLA8Edid9gARsTv5704gB5xeRHxmZlZCaZJ+B9AkabKkMcBCINVdOJJOkFSXfB4PnE3BtQAzM6usAZN+RPQCi4ENwHbg7ojYKmmlpAUAks6Q1A1cAnxL0tak+jRgo6QngUfIj+k76ZuZVUmqCdciYj2wvs+6pQWfO8gP+/St95/AzCHGaGZmJeIncs3MMsRJ38wsQzyfvg3Z0398gQ8vuT9V2a76MgdjZkfknr6ZWYa4p29DNnPU7+mqH2EPJpmNUO7pm5lliJO+mVmGOOmbmWWIk76ZWYY46ZuZZYiTvplZhjjpm5lliJO+mVmGOOmbmWWIk76ZWYakSvqS5knaIalT0pJ+tp8j6QlJvZIu7rPtckm/S/5dXqrAzcyseAMmfUmjgdXAhcB0YJGk6X2KPQt8GPhun7pvBJYBZwKzgWWSThh62GZmNhhpevqzgc6I2BkRrwLrgNbCAhHRFRFPAX/vU/cC4MGI2BsRzwMPAvNKELeZmQ1CmqR/MrCrYLk7WZfGUOqamVmJpZlaWf2si5T7T1VXUhvQBtDQ0EAul0u5+9rQ09Mz7GIupZ66CeSmrqh2GFU1Is9BET/TWf8OwPA5B2mSfjcwqWB5IrA75f67gZY+dXN9C0XEGmANQHNzc7S0tPQtUtNyuRzDLeZSyt31NVp2ZHs+/dzUFSPvHCx6IXXRrH8HYPicgzTDOx1Ak6TJksYAC4H2lPvfAMyVdEJyAXduss7MzKpgwKQfEb3AYvLJejtwd0RslbRS0gIASWdI6gYuAb4laWtSdy/wBfK/ODqAlck6MzOrglSvS4yI9cD6PuuWFnzuID9001/dtcDaIcRoZmYl4idyzcwyxEnfzCxDUg3vWLY0Lrm/qPJ3vL1MgZhZybmnb2aWIU76ZmYZ4qRvZpYhTvpmZhnipG9mliG+e8fM+lXMXVz/OrP3kEm2rHa5p29mliFO+mZmGeKkb2aWIR7Tt9fpqv9gUeVzjLCXhxhQ3M9BbpR/BoYL9/TNzDLESd/MLEOc9M3MMiRV0pc0T9IOSZ2SlvSzvU7S95Ltj0lqTNY3SnpF0ubk3+2lDd/MzIox4IVcSaOB1cAc8i8675DUHhHbCopdATwfEf8oaSFwM/CBZNszETGrxHFbMZaPq3YEZlYj0vT0ZwOdEbEzIl4F1gGtfcq0Ancmn+8F3itJpQvTzMxKIc0tmycDuwqWu4EzD1cmInolvQCcmGybLOk3wF+Bz0fEL/oeQFIb0AbQ0NBALpcrpg1V19PTU9sxTy3v7XQ9dRPIlfkYtS7r56CnbkJtfwcqoObzQCJN0u+vxx4py/wJeHNE7JH0TuBHkk6LiL8eUjBiDbAGoLm5OVpaWlKEVTtyuRw1HfPyvn+YlVZu6gpadiwr6zFqXdbPQW7qClpa/nu1w6iqms8DiTTDO93ApILlicDuw5WRdBQwDtgbEfsjYg9ARGwCngGmDDVoMzMbnDRJvwNokjRZ0hhgIdDep0w7cHny+WLg4YgISSclF4KRdArQBOwsTehmZlasAYd3kjH6xcAGYDSwNiK2SloJbIyIduDfgH+X1AnsJf+LAeAcYKWkXuAAcFVE7C1HQ8zMbGCp5t6JiPXA+j7rlhZ83gdc0k+97wPfH2KMZmZWIn4i18wsQzzLppmVRDFv2uq6aX4ZI7EjcdI3s5IobkruF8oWhx2Zh3fMzDLESd/MLEM8vGNmFVfM+D/4GkApuadvZpYhTvpmZhnipG9mliEe0x+u/GIUG8aKu70TfItn6binb2aWIU76ZmYZ4qRvZpYhTvpmZhnipG9mliGp7t6RNA+4lfxLVL4dETf12V4H/E/gncAe4AMR0ZVsuwG4gvxLVD4VERtKFr2ZZYKf4C2dAXv6yesOVwMXAtOBRZKm9yl2BfB8RPwj8FXg5qTudPJv0ToNmAfcdvD1iWZmVnlpevqzgc6I2AkgaR3QCmwrKNMKLE8+3wt8Q5KS9esiYj/w++R1irOBX5cm/JGj6J5MfZkCMatBRd/Xv7yYstl6BiDNmP7JwK6C5e5kXb9lIqKX/JMUJ6asa2ZmFZKmp69+1kXKMmnqIqkNaEsWeyTtSBFXLRkP/KWSB+zvxFbPNRVvf+3J+jkYxu1fUbJvU7XPwVvSFEqT9LuBSQXLE4HdhynTLekoYBywN2VdImINsCZNwLVI0saIaK52HNWS9faDz0HW2w/D5xykGd7pAJokTZY0hvyF2fY+ZdqBy5PPFwMPR0Qk6xdKqpM0GWgCHi9N6GZmVqwBe/oR0StpMbCB/C2bayNiq6SVwMaIaAf+Dfj35ELtXvK/GEjK3U3+om8vcHVEHChTW8zMbADKd8htKCS1JUNUmZT19oPPQdbbD8PnHDjpm5lliKdhMDPLECf9EpJ0naSQNL7asVSapC9L+r+SnpL0Q0nHVzumSpA0T9IOSZ2SllQ7nkqTNEnSI5K2S9oq6dPVjqkaJI2W9BtJ91U7loE46ZeIpEnAHODZasdSJQ8CMyLibcBvgRuqHE/ZpZyiZKTrBf41IqYBZwFXZ/AcAHwa2F7tINJw0i+drwKfoZ+Hz7IgIn6aPI0N8Cj5ZzJGutemKImIV4GDU5RkRkT8KSKeSD6/SD7xZeqpe0kTgfnAt6sdSxpO+iUgaQHwx4h4stqx1IiPAj+pdhAV4GlGCkhqBE4HHqtuJBX3NfIdvr9XO5A0/GL0lCQ9BPxDP5s+B/wPYG5lI6q8I52DiPjfSZnPkf+T/z8qGVuVpJpmJAskjQW+D/xLRPy12vFUiqSLgD9HxCZJLdWOJw0n/ZQi4vz+1kuaCUwGnsxPLMpE4AlJsyPivyoYYtkd7hwcJOly4CLgvZGNe4FTTTMy0kk6mnzC/4+I+EG146mws4EFkv4JqAfeIOl/RcSlVY7rsHyffolJ6gKaI2J4Tj41SMmLdm4B3hMRz1U7nkpI5pn6LfBe4I/kpyz5YERsrWpgFZRMoX4nsDci/qXa8VRT0tO/LiIuqnYsR+IxfSuVbwDHAQ9K2izp9moHVG7JheuDU5RsB+7OUsJPnA1cBpyX/H/fnPR6rUa5p29mliHu6ZuZZYiTvplZhjjpm5lliJO+mVmGOOmbmWWIk76ZWYY46ZuZZYiTvplZhvw/HyKhBLLiTq0AAAAASUVORK5CYII=\n",
            "text/plain": "<Figure size 432x288 with 1 Axes>"
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "text/plain": "0.05993882587548316"
          },
          "execution_count": 22,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "hist1",
      "execution_count": 66,
      "outputs": [
        {
          "data": {
            "text/plain": "array([9.99999990e-10, 9.99999990e-10, 9.99999990e-10, 4.00000096e-03,\n       1.20000009e-02, 4.40000006e-02, 1.10000000e-01, 1.97999999e-01,\n       2.49999998e-01, 4.23999997e-01, 3.71999997e-01, 3.07999998e-01,\n       1.59999999e-01, 8.60000001e-02, 2.80000007e-02, 4.00000096e-03,\n       9.99999990e-10, 9.99999990e-10, 9.99999990e-10, 9.99999990e-10])"
          },
          "execution_count": 66,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# sample 10000\nimport scipy as sp\nimport matplotlib.pyplot as plt\nimport scipy.stats as stats\n\nN = 10000\nsample1 = stats.norm(loc=0, scale=1).rvs(N)\nsample2 = stats.norm(loc=0, scale=1).rvs(N)\n\nhist1, hist2 = create_hists(sample1, sample2, bins, eps=1e-9)\n\nplt.bar(bins[:-1], hist1, label=\"sample 1\")\nplt.bar(bins[:-1], hist2, label=\"sample 2\")\nplt.legend()\nplt.grid(True)\nplt.savefig(\"hist_10000.png\")\nplt.show()\n\nD_KL(hist1, hist2) ",
      "execution_count": 23,
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAG2lJREFUeJzt3X9wXWW97/H3p4Umg8XCoUyObaopmtaWVlsNBQcPBqWl3DLtjANXijjo9RhQqh4RtFwd+sMyB38McmasYsdT4Z57bAX8cSNWKgj7+OOIpMUC/XGroUaI9YzYMkiAFlO/94/scndD2qyd7OyV3efzmumw11rPs/Z3LbI/efLsvddSRGBmZmkYk3cBZmZWPQ59M7OEOPTNzBLi0DczS4hD38wsIQ59M7OEOPTNzBLi0DczS4hD38wsISfkXUB/EydOjKamprzLKMvzzz/Pq171qrzLyE3qxw8+B6kfP+R/DrZu3frniDh9sHajLvSbmprYsmVL3mWUpVAo0NramncZuUn9+MHnIPXjh/zPgaTfZ2nn6R0zs4Q49M3MEuLQNzNLSKY5fUkLgX8BxgLfiIibj9LuEuAu4KyI2FJcdwPwQeAQ8LGI2FyJws2sNv31r3+lu7ubAwcO5F1KRU2YMIFdu3aN+PPU19fT2NjIiSeeOKT+g4a+pLHAWmA+0A10SGqPiJ392p0MfAz4Vcm6mcBlwJnAJOB+SdMi4tCQqjWzmtfd3c3JJ59MU1MTkvIup2Kee+45Tj755BF9johg3759dHd3M3Xq1CHtI8v0zjygMyL2RMRLwEZgyQDtPgd8ASj99b0E2BgRByPid0BncX9mlqgDBw5w2mmnHVeBXy2SOO2004b1V1KW0J8MPFWy3F1cV1rIXGBKRNxTbl8zS48Df+iGe+6yzOkP9Awv32NR0hjgy8D7y+1bso82oA2goaGBQqGQoazRo6enp+ZqrqTUjx98Dso5/gkTJvDcc8+NbEE5OHToUNWO68CBA0P+ecsS+t3AlJLlRmBvyfLJwCygUPwN9PdAu6TFGfoCEBHrgHUALS0tUWtf8sj7Sxl5S/34weegnOPftWvXEXPfTct/WNFaum5eVNH9ZTXQnH5raytf+tKXaGlpybSPu+66i5UrV7Jr1y4efvjho/arr69n7ty5Q6ozS+h3AM2SpgJ/oO+N2csPb4yIZ4GJh5clFYDrImKLpBeBb0m6hb43cpuBh4dUqVmFVTJsPjm7l/cfZX95hZDVnlmzZvHd736Xq666asSeY9A5/YjoBZYBm4FdwJ0RsUPS6uJo/lh9dwB3AjuBe4Fr/MkdM8vL888/z6JFi3jzm9/MrFmz+Pa3vw3A6tWrOeuss5g1axZtbW1E9M1Ct7a28olPfILzzjuPGTNm0NHRwbvf/W6am5v57Gc/C0BXVxdvfOMbueqqq3jTm97EJZdcwgsvvPCK5/7xj3/M2972Nt7ylrdw6aWX0tPT84o2M2bMYPr06SN4BjJ+OSsiNkXEtIh4fUTcVFx3Y0S0D9C29fBn9IvLNxX7TY+IH1WudDOz8tx7771MmjSJRx99lO3bt7Nw4UIAli1bRkdHB9u3b+fFF1/knnv+/2dSxo0bx09/+lOuvvpqlixZwtq1a9m+fTu33347+/btA2D37t184AMf4LHHHuPVr341X/3qV4943j//+c+sWbOG+++/n0ceeYSWlhZuueWW6h14CX8j18ySMXv2bO6//34+/elP87Of/YwJEyYA8OCDD3L22Wcze/ZsHnjgAXbs2PFyn8WLF7/c98wzz+Q1r3kNdXV1nHHGGTz1VN+HE6dMmcI555wDwBVXXMHPf/7zI573oYceYufOnZx77rnMmTOHO+64g9//PtP10Spu1F1l06xauuovH7xRRoUxq+iqX3GUrc9W7HlseKZNm8bWrVvZtGkTN9xwAwsWLOBTn/oUH/nIR9iyZQtTpkxh5cqVR3wOvq6uDoAxY8a8/Pjwcm9vL/DKj1H2X44I5s+fz4YNG0bq0DLzSN/MkrF3715OOukkrrjiCq677joeeeSRlwN+4sSJ9PT0cPfdd5e93yeffJJf/arvYgQbNmzg7W9/+xHbzznnHH7xi1/Q2dkJwAsvvMBvfvObYR7N0Hikb2a5quanmx5//HGuv/56xowZw4knnsjXvvY1TjnlFD70oQ8xe/ZsmpqaOOuss8re74wZM9iwYQPXXnstzc3NfPjDHz5i++mnn87tt9/O0qVLOXjwIABr1qxh2rRpR7T73ve+x0c/+lGefvppFi1axJw5c9i8ubKXK9Phd6lHi5aWlvBNVGpLzR7/ygkV21Vh+ipadx9lemfl8T+9U+7n9GfMmDGyBVVRV1cXF198Mb/85S9H/No7hw10DiVtjYhBvxDg6R0zs4Q49M3MhqGpqYnt27fnXUZmntM3G2nlTCMlMBVk+fJI38wsIQ59M7OEeHrHjh8V/DSO2fHKoW9m+ar0L+tR9L5IuZdWvv766/nBD37AuHHjeP3rX883v/lNTjnllIrW5OkdM7NRYv78+Wzfvp3HHnuMadOm8c///M8Vfw6HvpklY7RfWnnBggWccELfBMw555xDd3d3xc+BQ9/MklFLl1Zev349F110UYXPgEPfzBJSK5dWvummmzjhhBN473vfW9HjB7+Ra2YJqYVLK99xxx3cc889/OQnP3nFfioh00hf0kJJuyV1Slo+wParJT0uaZukn0uaWVzfJOnF4vptkm6r9AGYmWU12i+tfO+99/L5z3+e9vZ2TjrppLLryGLQkb6kscBaYD7QDXRIao+InSXNvhURtxXbLwZuARYWtz0REXMqW7aZHTeq+BHL0X5p5WXLlnHw4EHmz58P9P2yuO22yo6Vs0zvzAM6I2IPgKSNwBL6bnYOQET8paT9q4DRdb1mMzPgwgsv5MILL3zF+jVr1rBmzZpXrC8UCi8/bm1tPeLy0Ye3dXV1MWbMGG699dZXXFq5tP873/lOOjo6jlnf4b8ERlKW6Z3JwFMly93FdUeQdI2kJ4AvAB8r2TRV0q8l/YekfxhWtWZmNiyD3kRF0qXAhRHxj8Xl9wHzIuKjR2l/ebH9lZLqgPERsU/SW4HvA2f2+8sASW1AG0BDQ8NbN27cONzjqqqenh7Gjx+fdxm5GTXH/8dtuT11T90kxh/cO/wdvaY2Z0LL+RmYMGECb3jDG0a4ouo7dOgQY8eOrcpzdXZ28uyzR06LnX/++ZluopJleqcbmFKy3Agc66d7I/A1gIg4CBwsPt5a/EtgGnDErbEiYh2wDvrunFVrd2Gq2TtHVcioOf6VS3J76mPeOascS0fPJQTKUe6ds8aPHz8in0zJ03PPPVeVO2dFBPX19cydO3dI/bNM73QAzZKmShoHXAa0lzaQ1FyyuAj4bXH96cU3gpF0BtAM7BlSpWZ2XKivr2ffvn2Mtlu11oKIYN++fdTX1w95H4OO9COiV9IyYDMwFlgfETskrQa2REQ7sEzSBcBfgWeAK4vdzwNWS+oFDgFXR8T+IVdrZjWvsbGR7u5unn766bxLqagDBw4MK4yzqq+vp7Gxccj9M305KyI2AZv6rbux5PHHj9LvO8B3hlydmR13TjzxRKZOnZp3GRVXKBSGPOVSTb4Mg5lZQhz6ZmYJceibmSXEoW9mlhCHvplZQhz6ZmYJceibmSXEoW9mlhCHvplZQhz6ZmYJceibmSXEoW9mlpBMF1wzsypZOaHM9rV5/X3Lj0f6ZmYJceibmSXEoW9mlhDP6dvoVu4ct5kdU6aRvqSFknZL6pS0fIDtV0t6XNI2ST+XNLNk2w3FfrslXVjJ4s3MrDyDhn7xxuZrgYuAmcDS0lAv+lZEzI6IOcAXgFuKfWfSdyP1M4GFwFcP3yjdzMyqL8tIfx7QGRF7IuIlYCOwpLRBRPylZPFVwOHb3C8BNkbEwYj4HdBZ3J+ZmeUgy5z+ZOCpkuVu4Oz+jSRdA1wLjAPeWdL3oX59Jw+pUjMzG7Ysoa8B1sUrVkSsBdZKuhz4LHBl1r6S2oA2gIaGBgqFQoayRo+enp6aq7mSRvT4p68amf1WWE/dJAp51DpKfu5Sfw1A7ZyDLKHfDUwpWW4E9h6j/Ubga+X0jYh1wDqAlpaWaG1tzVDW6FEoFKi1mitpRI9/5ZLB24wChemraN29ovpPvHR0fCM39dcA1M45yDKn3wE0S5oqaRx9b8y2lzaQ1FyyuAj4bfFxO3CZpDpJU4Fm4OHhl21mZkMx6Eg/InolLQM2A2OB9RGxQ9JqYEtEtAPLJF0A/BV4hr6pHYrt7gR2Ar3ANRFxaISOxczMBpHpy1kRsQnY1G/djSWPP36MvjcBNw21QDMzqxxfhsHMLCEOfTOzhDj0zcwS4tA3M0uIQ9/MLCEOfTOzhDj0zcwS4tA3M0uIQ9/MLCEOfTOzhDj0zcwS4tA3M0uIQ9/MLCEOfTOzhDj0zcwS4tA3M0uIQ9/MLCEOfTOzhGQKfUkLJe2W1Clp+QDbr5W0U9Jjkn4i6XUl2w5J2lb8196/r5mZVc+g98iVNBZYC8wHuoEOSe0RsbOk2a+Bloh4QdKHgS8A7yluezEi5lS4bjMzG4IsI/15QGdE7ImIl4CNwJLSBhHxYES8UFx8CGisbJlmZlYJiohjN5AuARZGxD8Wl98HnB0Ry47S/ivAf0XEmuJyL7AN6AVujojvD9CnDWgDaGhoeOvGjRuHfkQ56OnpYfz48XmXkZsRPf4/bhuZ/VZYT90kxh/cW/0nfs3o+CM69dcA5H8Ozj///K0R0TJYu0GndwANsG7A3xSSrgBagHeUrH5tROyVdAbwgKTHI+KJI3YWsQ5YB9DS0hKtra0Zyho9CoUCtVZzJY3o8a9cMnibUaAwfRWtu1dU/4mXPlv95xxA6q8BqJ1zkGV6pxuYUrLcCLxiSCPpAuAzwOKIOHh4fUTsLf53D1AA5g6jXjMzG4Ysod8BNEuaKmkccBlwxKdwJM0Fvk5f4P+pZP2pkuqKjycC5wKlbwCbmVkVDTq9ExG9kpYBm4GxwPqI2CFpNbAlItqBLwLjgbskATwZEYuBGcDXJf2Nvl8wN/f71I+ZmVVRljl9ImITsKnfuhtLHl9wlH7/CcweToFmZlY5/kaumVlCHPpmZglx6JuZJcShb2aWEIe+mVlCMn16x8xGp6blPyyrfdfNi0aoEqsVHumbmSXEI32rqrJHpvUjVIhZojzSNzNLiEf6VlVd9ZfnXYJZ0hz6ZjWs/F+io+NSzJYfT++YmSXEoW9mlhCHvplZQhz6ZmYJceibmSXEoW9mlpBMoS9poaTdkjolLR9g+7WSdkp6TNJPJL2uZNuVkn5b/HdlJYs3M7PyDBr6ksYCa4GLgJnAUkkz+zX7NdASEW8C7ga+UOz7d8AK4GxgHrBC0qmVK9/MzMqRZaQ/D+iMiD0R8RKwEVhS2iAiHoyIF4qLDwGNxccXAvdFxP6IeAa4D1hYmdLNzKxcWb6ROxl4qmS5m76R+9F8EPjRMfpO7t9BUhvQBtDQ0EChUMhQ1ujR09NTczVXUlnHP33ViNaSl566SRRq4dhG6Oc09dcA1M45yBL6GmBdDNhQugJoAd5RTt+IWAesA2hpaYnW1tYMZY0ehUKBWqu5kso6/pVLBm9TgwrTV9G6e0XeZQxu6chchiH11wDUzjnIMr3TDUwpWW4E9vZvJOkC4DPA4og4WE5fMzOrjiyh3wE0S5oqaRxwGdBe2kDSXODr9AX+n0o2bQYWSDq1+AbuguI6MzPLwaDTOxHRK2kZfWE9FlgfETskrQa2REQ78EVgPHCXJIAnI2JxROyX9Dn6fnEArI6I/SNyJGZmNqhMl1aOiE3Apn7rbix5fMEx+q4H1g+1QDMzqxx/I9fMLCEOfTOzhDj0zcwS4tA3M0uIQ9/MLCEOfTOzhDj0zcwS4tA3M0uIQ9/MLCEOfTOzhDj0zcwS4tA3M0uIQ9/MLCEOfTOzhDj0zcwS4tA3M0uIQ9/MLCGZQl/SQkm7JXVKWj7A9vMkPSKpV9Il/bYdkrSt+K+9f18zM6ueQW+XKGkssBaYD3QDHZLaI2JnSbMngfcD1w2wixcjYk4FajUzs2HKco/ceUBnROwBkLQRWAK8HPoR0VXc9rcRqNHMzCokS+hPBp4qWe4Gzi7jOeolbQF6gZsj4vv9G0hqA9oAGhoaKBQKZew+fz09PTVXcyWVdfzTV41oLXnpqZtEoRaObYR+TlN/DUDtnIMsoa8B1kUZz/HaiNgr6QzgAUmPR8QTR+wsYh2wDqClpSVaW1vL2H3+CoUCtVZzJRU23ErrlhV5l5GrwvRVtO6ugXOw9NkR2W3qrwGonXOQ5Y3cbmBKyXIjsDfrE0TE3uJ/9wAFYG4Z9ZmZWQVlCf0OoFnSVEnjgMuATJ/CkXSqpLri44nAuZS8F2BmZtU1aOhHRC+wDNgM7ALujIgdklZLWgwg6SxJ3cClwNcl7Sh2nwFskfQo8CB9c/oOfTOznGSZ0yciNgGb+q27seRxB33TPv37/Scwe5g1mlmFNC3/Yea2XTcvGsFKLC/+Rq6ZWUIyjfTN7PjQVX95Ga1H5pM+li+P9M3MEuLQNzNLiEPfzCwhDn0zs4Q49M3MEuLQNzNLiEPfzCwhDn0zs4Q49M3MEuLQNzNLiEPfzCwhDn0zs4Q49M3MEuLQNzNLiEPfzCwhmUJf0kJJuyV1Slo+wPbzJD0iqVfSJf22XSnpt8V/V1aqcDMzK9+goS9pLLAWuAiYCSyVNLNfsyeB9wPf6tf374AVwNnAPGCFpFOHX7aZmQ1FlpH+PKAzIvZExEvARmBJaYOI6IqIx4C/9et7IXBfROyPiGeA+4CFFajbzMyGIEvoTwaeKlnuLq7LYjh9zcyswrLcI1cDrIuM+8/UV1Ib0AbQ0NBAoVDIuPvRoaenp+ZqrqSeukkUpq/Ku4xcHZfnoIyf6dRfA1A75yBL6HcDU0qWG4G9GfffDbT261vo3ygi1gHrAFpaWqK1tbV/k1GtUChQazVXUmHDrbTuXpF3GbkqTF91/J2DpdlvjJ76awBq5xxkmd7pAJolTZU0DrgMaM+4/83AAkmnFt/AXVBcZ2ZmORg09COiF1hGX1jvAu6MiB2SVktaDCDpLEndwKXA1yXtKPbdD3yOvl8cHcDq4jozM8tBlukdImITsKnfuhtLHnfQN3UzUN/1wPph1GhmZhXib+SamSUk00jfzNLTtPyHmdt+cnbvEZ/YsNHLI30zs4Q49M3MEuLpHXullRPKa3+8fSnJ7Djmkb6ZWUIc+mZmCXHom5klxHP6ZjagrvrLM7ctjPH7OrXCI30zs4Q49M3MEuLQNzNLiEPfzCwhDn0zs4Q49M3MEuLQNzNLiEPfzCwhmUJf0kJJuyV1Slo+wPY6Sd8ubv+VpKbi+iZJL0raVvx3W2XLNzOzcgz6jVxJY4G1wHygG+iQ1B4RO0uafRB4JiLeIOky4PPAe4rbnoiIORWu28zMhiDLSH8e0BkReyLiJWAjsKRfmyXAHcXHdwPvkqTKlWlmZpWQJfQnA0+VLHcX1w3YJiJ6gWeB04rbpkr6taT/kPQPw6zXzMyGIcsF1wYasUfGNn8EXhsR+yS9Ffi+pDMj4i9HdJbagDaAhoYGCoVChrJGj56enpqr+ZjKvClKT90kConfSCX1c9BTN+n4eg0MQa3kQJbQ7wamlCw3AnuP0qZb0gnABGB/RARwECAitkp6ApgGbCntHBHrgHUALS0t0draWv6R5KhQKFBrNR/Tyv6zd8dWmL6K1t0rRqiY2pD6OShMX0Vr63/Pu4xc1UoOZJne6QCaJU2VNA64DGjv16YduLL4+BLggYgISacX3whG0hlAM7CnMqWbmVm5Bh3pR0SvpGXAZmAssD4idkhaDWyJiHbgX4F/k9QJ7KfvFwPAecBqSb3AIeDqiNg/EgdiR9e0/Idlte+qH6FCzCx3mW6iEhGbgE391t1Y8vgAcOkA/b4DfGeYNZqZWYX4zlkJKOcOSGZDVc5flF03LxrBSuxYfBkGM7OEeKRvZhVR3l+Uz45YHXZsHumbmSXEoW9mlhCHvplZQhz6ZmYJceibmSXEoW9mlhCHvplZQvw5/Vq1ckLeFZhZDfJI38wsIR7pm1nVlX3lV1+rp2I80jczS4hH+mZWdeVf+dXX6qkUj/TNzBLikb6ZjX7lflptpf8yOJpMoS9pIfAv9N0u8RsRcXO/7XXA/wLeCuwD3hMRXcVtNwAfpO92iR+LiM0Vq/544o9gmlkVDDq9U7yx+VrgImAmsFTSzH7NPgg8ExFvAL4MfL7YdyZ998s9E1gIfPXwjdLNzKz6soz05wGdEbEHQNJGYAmws6TNEmBl8fHdwFckqbh+Y0QcBH5XvHH6POCXlSl/FPPI3cxGoSyhPxl4qmS5Gzj7aG0iolfSs8BpxfUP9es7ecjV5u1oQT59FaxcUt1azOzoRnLQVePvF2QJfQ2wLjK2ydIXSW1AW3GxR9LuDHWNIp+YCPw57yryk/rxg89BQse/aqBYAyDvc/C6LI2yhH43MKVkuRHYe5Q23ZJOACYA+zP2JSLWAeuyFDwaSdoSES1515GX1I8ffA5SP36onXOQ5XP6HUCzpKmSxtH3xmx7vzbtwJXFx5cAD0REFNdfJqlO0lSgGXi4MqWbmVm5Bh3pF+folwGb6fvI5vqI2CFpNbAlItqBfwX+rfhG7X76fjFQbHcnfW/69gLXRMShEToWMzMbhPoG5DYcktqKU1RJSv34wecg9eOH2jkHDn0zs4T42jtmZglx6FeQpOskhaSJeddSbZK+KOn/SnpM0vcknZJ3TdUgaaGk3ZI6JS3Pu55qkzRF0oOSdknaIenjedeUB0ljJf1a0j151zIYh36FSJoCzAeezLuWnNwHzIqINwG/AW7IuZ4Rl/ESJce7XuCTETEDOAe4JsFzAPBxYFfeRWTh0K+cLwOfYoAvn6UgIn4cEb3FxYfo+07G8e7lS5RExEvA4UuUJCMi/hgRjxQfP0df8NXut+6HQFIjsAj4Rt61ZOHQrwBJi4E/RMSjedcySvwP4Ed5F1EFA12iJKnAKyWpCZgL/CrfSqruVvoGfH/Lu5AsfD39jCTdD/z9AJs+A/xPYEF1K6q+Y52DiPg/xTafoe9P/n+vZm05yXSZkRRIGg98B/iniPhL3vVUi6SLgT9FxFZJrXnXk4VDP6OIuGCg9ZJmA1OBR/suLEoj8IikeRHxX1UsccQd7RwcJulK4GLgXZHGZ4EzXWbkeCfpRPoC/98j4rt511Nl5wKLJf03oB54taT/HRFX5FzXUflz+hUmqQtoiYg0Lj5VVLzRzi3AOyLi6bzrqYbidaZ+A7wL+AN9lyy5PCJ25FpYFRUvoX4HsD8i/invevJUHOlfFxEX513LsXhO3yrlK8DJwH2Stkm6Le+CRlrxjevDlyjZBdyZUuAXnQu8D3hn8f/7tuKo10Ypj/TNzBLikb6ZWUIc+mZmCXHom5klxKFvZpYQh76ZWUIc+mZmCXHom5klxKFvZpaQ/wfT44pH2tD13gAAAABJRU5ErkJggg==\n",
            "text/plain": "<Figure size 432x288 with 1 Axes>"
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "text/plain": "0.008646225250418648"
          },
          "execution_count": 23,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "### 二項分布と正規分布"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# sample 10000\nimport scipy as sp\nimport matplotlib.pyplot as plt\nimport scipy.stats as stats\n\nN = 10000\nn = 10\np = 1/2\nsample1 = stats.binom(n, p).rvs(N)\nsample2 = stats.norm(loc=n*p, scale=sp.sqrt(n*p*(1-p))).rvs(N)\n\nbins = sp.linspace(5-8,5+8, 20+1)\nhist1, hist2 = create_hists(sample1, sample2, bins, eps=1e-9)\n\nplt.bar(bins[:-1], hist1, label=\"bin(10, 1/2)\")\nplt.bar(bins[:-1], hist2, label=\"$N(np, np(1-p))$\")\nplt.legend()\nplt.grid(True)\nplt.savefig(\"Bin(10,0.5)_hist.png\")\nplt.show()\n\nD_KL(hist1, hist2) ",
      "execution_count": 25,
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHvpJREFUeJzt3Xt4VPW97/H313AJNWxUtKkSDuDTbJqYIJcUaLU2FFBED6BCjVrFisVeEE/Rp4cWK+BuK7WtVSvdlSqnthubtvYij6JsEXOOdrcVRLkEDhKQYgSUiwcJ9+D3/DGTYZhMkjXJJDO4Pq/n4WF+a/1+a30nDJ+s+c1aa8zdERGRcDgt0wWIiEjHUeiLiISIQl9EJEQU+iIiIaLQFxEJEYW+iEiIKPRFREJEoS8iEiIKfRGREOmU6QISnX322d63b99Y+8CBA5x++umZK6gJqis1qis12VoXZG9tYa/rtdde2+3u57TY0d2z6s+QIUM83ksvveTZSHWlRnWlJlvrcs/e2sJeF7DSA2SspndEREJEoS8iEiIKfRGREMm6D3JFpH0cO3aM2tpaDh8+3Kbt9OjRgw0bNqSpqvQJS125ubkUFBTQuXPnVo1X6IuERG1tLd27d6dv376YWau3s3//frp3757GytIjDHW5O3v27KG2tpZ+/fq1ahua3hEJicOHD9OzZ882Bb5klpnRs2fPNr1bU+iLhIgC/9TX1n/DQKFvZmPMbKOZ1ZjZzCTrv2pma83sDTN7xcyK49Z9Ozpuo5ld1qZqRUSkTVqc0zezHGA+MBqoBVaY2WJ3Xx/X7Ul3/0W0/zjgAWBMNPwrgAuA84BlZvav7n48zc9DRFLUd+azad3e1nlXNL9+61auvPJK1q1b12jdrbfeyowZMyguLk4y8oQHH3yQs846i5tuuok//OEPzJkzhw0bNvDqq6/Sv3//WL/77ruPxx9/nJycHB5++GEuu6z5481HHnmEBx98kM2bN7Nr1y7OPvvs2Lpjx44xfPhw/vKXv3DTTTexc+dOTjvtNKZOncodd9wBwF133cXYsWP5whe+0Ox+skGQD3KHAjXuvgXAzCqB8UAs9N39g7j+pwMN37Y+Hqh09yPAW2ZWE93e39JQu0iHa0tQthSKYfbYY4+12Ke+vp6FCxeyatUqAEpKSvjTn/7EbbfddlK/9evXU1lZSXV1Ndu3b2fUqFG8+eab5OTkNLntiy66iCuvvJLy8vJG61555RU++9nP0qlTJ37yk58wePBg9u/fz5AhQxg9ejTFxcXcfvvtfOUrXzklQj/I9E4v4O24dm102UnM7Btmthm4H5ieylgRCYf6+nomT57MgAEDmDhxIgcPHgSgvLyclStXApCXl8esWbO48MILGT58OO+++y4Ay5cvZ/DgwXTqFDlWLSoqOunovsHTTz9NRUUFXbt2pV+/fnzyk5/k1VdfbbauQYMGEX/Pr3jPP/88l19+Oeeeey6DBw8GoHv37hQVFfHOO+8A0KdPH/bs2cPOnTtT/6F0sCBH+sk+NfBGC9znA/PN7HrgbmBy0LFmNhWYCpCfn09VVVVsXV1d3UntbKG6UvNRqevO0vpW7yuV/bTHz6tHjx7s378/rduM19K26+rq2LhxIz/72c945JFH+PrXv85Pf/pTpk+fzvHjxzlw4AD79+/nwIEDXHjhhcycOZPvfve7PPLII3zrW99i+fLlXHDBBY320zD2+PHj7N+/n7feeotPf/rTsX75+fnU1NRQUlLS4nNwd+rq6ujatWts2bJly5gxY8ZJ+/3nP//JqlWrKC4uji0vLS1l2bJljB8/vlF96f65Hz58uNWvjyChXwv0jmsXANub6V8J/HsqY919AbAAoKyszOPfYlVVVSV9y5Vpqis1H5W6bm7L9M4NwffTHj+vDRs2tOt57C1tOy8vj969ezN69GgAvvzlL/Pwww/TvXt3cnJyOP300+nevTtdunRh0qRJmBmf+cxneOGFF+jevTt79+5l4MCBjfbTMDYnJ4fu3bvTuXNnunXrFuvXuXNnPvaxjwV67mZGXl5erO/27ds555xzyM/Pj/Wpq6tj8uTJPPTQQ/TqdWLiolevXrz//vuN9tMe1w/k5uYyaNCgVo0NMr2zAig0s35m1oXIB7OL4zuYWWFc8wpgU/TxYqDCzLqaWT+gEGj+fZaIfGQlnm6Y7PTDzp07x5bn5ORQXx95d9WtW7dA56cXFBTw9tsnZpVra2s577zzWlXvc889d9KHwMeOHeOaa67hhhtu4Oqrrz6p7+HDh+nWrVur9tORWgx9d68HpgFLgQ3A79292szujZ6pAzDNzKrN7A1gBpGpHdy9Gvg9kQ99nwe+oTN3RMJr27Zt/O1vkfM4fvvb33LxxRcHHltUVERNTU2L/caNG0dlZSVHjhzhrbfeYtOmTQwdOhSAkSNHxubhg2iYz4fI1M+UKVMoKipixowZjfq++eabgaaQMi3QbRjcfQmwJGHZPXGP72hm7PeB77e2QBFpH609m6gt0xVFRUU88cQT3HbbbRQWFvK1r30t8NjLL7+cG2+8Mdb+85//zO23386uXbu44oorKCkp4cUXX+SCCy7gi1/8IsXFxXTq1In58+eTk5PDhx9+SE1NDWeddVajbT/88MPcf//97Ny5kwEDBjB27FgeffRRNm3axKc+9SkA/vrXv/Kb3/yG0tJSBg4cCMAPfvADxo4dy7Fjx6ipqaGsrKxVP5eOpHvviEiH6Nu3L+vXr0+6LvHkjQYTJ05k4sSJQOQMmZ49e7Jp0yYKCwu56qqruOqqq2J94z8snTVrFrNmzTppH+vXr+eaa65JOgUzffp0pk+fftKyV155heHDh8faF198MZHvKmnsmWeeYeLEibEzi7JZ9lcoIhI1b948duzYQWFhYcudE5SUlPDAAw8E7n/xxRcHnn6qr6/nzjvvTLmmTFDoi8gpo3///knPzc+0SZMmZbqEwHTDNRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCRGFvohIiOiUTZGwmtOjVcOavBZ3zr5A4x999FG++tWvsn79eoqKioDIlbrPPfcc+fn5jBkzhuXLlzd7//v2dujQoZPquOWWW3jmmWf4+Mc/nvRLYNrq6NGjjBo1iuXLl9OpU6dG7XTSkb6IdKg1a9YwcOBAnn02csfSI0eO8O6779KnTx8WLlzI1VdfndHABxrVcfPNN/P888+32/66dOnCyJEj+d3vfpe0nU4KfRHpUGvXrmXmzJmx0K+urqaoqAgzY9GiRSfdj/6qq67i7rvv5nOf+xyf+MQnWLZsGQAVFRVce+21DBs2jD59+sS21ZRUt5NYxyWXXJL0nj2t0dQ+J0yYwKJFi2L9EtvpotAXkQ61fv16xo0bx3vvvce+fftYu3YtpaWlHD16lC1btpz0DVbr1q3jjDPO4OWXX+bnP/95LARXr17N+eefzz/+8Q8WLVrE3Llzm91nKttJVkc6NVV7SUkJK1asiPVLbKeL5vRFpMO8/fbb9OzZk27dujF69GiWLl3KmjVrGDBgALt37+aMM86I9T148CD79u3jm9/8JhC5v80ZZ5zBoUOH2L17N7NnzwaguLiY999/v8l9prqdxDqCGjVqVNKvS7z77rupqKgAaLb2nJwcunTpEruLaWI7XRT6ItJh1qxZQ2lpKQBjx45l0aJF7NixgwkTJjT6kpTq6mqGDBkSm1dfs2YNJSUlrFu3jsLCQnJzcwFYtWoVF154YZP7THU7Qb+sJVHDlFGi+Lt/tlT7kSNHYuuStdNB0zsi0mEapnIAPv/5z/Pyyy/HfhGceeaZHD9+PBa469ati923Hoi9I1i9ejXbtm3j8OHDHDhwgNmzZ8eO4qHxF6Wkup3EOtKpudr37NnDOeecQ+fOnZO200VH+iJhFfAUy0RtmW5Yu3Yt11xzDQBdu3altLSU119/PTadcumll/LKK68watQo1q5dy7Bhw2Jj161bR0lJCb/+9a+54YYbKC8v54MPPuA73/kOF110Efv370/6RSmpbiexDoDrrruOqqoqdu/eTUFBAXPnzmXKlCkpP//Vq1c3uc+XXnqJsWPHxvomttNFoS8iHSbxbJSnn376pPa0adN44IEHGDVqVKN732/ZsgWIBOcvf/lLfvjDHzbafrIvSmnNduLrgMhXO6ZDc/t88sknue+++5psp4umd0QkawwaNIgRI0Zw/HjTX6W9efPmJr9EJZUvSmluO0HqaI2m9nn06FEmTJgQ+66AxHY66UhfRLLKLbfc0uz6VL7YvC3baamOdO6zS5cu3HTTTU2200lH+iIiIaIjfZEO0ndm81eNxruztJ6b4/pvnXdFe5QkIaQjfZEQcfdMlyBt1NZ/Q4W+SEjk5uayZ88eBf8pzN3Zs2dPmy7Y0vSOSEgUFBRQW1vLrl272rSdw4cPp/0q0XQIS125ubkUFBS0enyg0DezMcBDQA7wmLvPS1g/A7gVqAd2Abe4+z+j644Da6Ndt7n7uFZXKyKt1rlzZ/r169fm7VRVVTFo0KA0VJReqiuYFkPfzHKA+cBooBZYYWaL3X19XLfXgTJ3P2hmXwPuB66Nrjvk7gMRSZNUPhBNRh+KSpgFmdMfCtS4+xZ3PwpUAuPjO7j7S+5+MNr8O9D69x4iItJugoR+L+DtuHZtdFlTpgDPxbVzzWylmf3dzCa0okYREUkTa+mTfDObBFzm7rdG2zcCQ9399iR9vwRMAz7v7keiy85z9+1mdj6wHBjp7psTxk0FpgLk5+cPqaysjK2rq6sjLy+vDU+xfaiu1KSzrrXvtO5GYQ1Ke534bthU62rrvoPK7wbvHjrRjq8508LwGkunjqprxIgRr7l7WUv9gnyQWwv0jmsXANsTO5nZKGAWcYEP4O7bo39vMbMqYBBwUui7+wJgAUBZWZmXl5fH1lVVVRHfzhaqKzXprOvmts7p33CijlTrauu+g7qztJ6frD3x3zO+5kwLw2ssnbKtriDTOyuAQjPrZ2ZdgApgcXwHMxsEPAqMc/f34pafaWZdo4/PBi4C4j8AFhGRDtTikb6715vZNGApkVM2F7p7tZndC6x098XAj4A84A9mBidOzSwCHjWzD4n8gpmXcNaPiIh0oEDn6bv7EmBJwrJ74h6PamLcfwGlbSlQJO3mxM2P958Lc8Y33beRJ9NejkhH0m0YRERCRLdhEEnB1tzrWz2272G9S5DM05G+iEiIKPRFREJEoS8iEiIKfRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCRGFvohIiCj0RURCRKEvIhIiCn0RkRBR6IuIhIhCX0QkRBT6IiIhotAXEQkRhb6ISIgo9EVEQkShLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaLQFxEJkUChb2ZjzGyjmdWY2cwk62eY2XozW2NmL5pZn7h1k81sU/TP5HQWLyIiqWkx9M0sB5gPXA4UA9eZWXFCt9eBMncfADwF3B8dexYwGxgGDAVmm9mZ6StfRERSEeRIfyhQ4+5b3P0oUAmMj+/g7i+5+8Fo8+9AQfTxZcAL7r7X3d8HXgDGpKd0ERFJlbl78x3MJgJj3P3WaPtGYJi7T2ui/yPATnf/npndBeS6+/ei674LHHL3HyeMmQpMBcjPzx9SWVkZW1dXV0deXl5rn1+7UV2pSWdda9/Z16bxpae9FXtc1/U88o5sb2tJgaz9sF/gvvnd4N1DJ9qlvXq0Q0WtE4bXWDp1VF0jRox4zd3LWurXKcC2LMmypL8pzOxLQBnw+VTGuvsCYAFAWVmZl5eXx9ZVVVUR384Wqis16azr5pnPtmn81tzZscdV/edSvnF2M73T5+bDTwbue2dpPT9Ze+K/59YbytuhotYJw2ssnbKtriDTO7VA77h2AdDo0MjMRgGzgHHufiSVsSIi0jGChP4KoNDM+plZF6ACWBzfwcwGAY8SCfz34lYtBS41szOjH+BeGl0mIiIZ0OL0jrvXm9k0ImGdAyx092ozuxdY6e6LgR8BecAfzAxgm7uPc/e9ZvZvRH5xANzr7nvb5ZmIiEiLgszp4+5LgCUJy+6JezyqmbELgYWtLVBERNJHV+SKiISIQl9EJEQU+iIiIRJoTl8km2zNvT7TJbRKKnVXnTb3pOsJoG0XpIk00JG+iEiIKPRFREJEoS8iEiIKfRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCRGFvohIiCj0RURCRKEvIhIiCn0RkRBR6IuIhIhCX0QkRBT6IiIhotAXEQkRhb6ISIgo9EVEQkShLyISIgp9EZEQUeiLiIRIoNA3szFmttHMasxsZpL1l5jZKjOrN7OJCeuOm9kb0T+L01W4iIikrlNLHcwsB5gPjAZqgRVmttjd18d12wbcDNyVZBOH3H1gGmoVEZE2ajH0gaFAjbtvATCzSmA8EAt9d98aXfdhO9QoIiJpYu7efIfIdM0Yd7812r4RGObu05L0/RXwjLs/FbesHngDqAfmuftfkoybCkwFyM/PH1JZWRlbV1dXR15eXurPrJ2prtSkta4db6RnO0Bd1/PIO7I9bdtLl0Z1nZs9b5ZD8RpLo46qa8SIEa+5e1lL/YIc6VuSZc3/pjjZf3P37WZ2PrDczNa6++aTNua+AFgAUFZW5uXl5bF1VVVVxLezhepKTVrrmjM+PdsBqvrPpXzj7LRtL10a1XXdvswVkyAUr7E0yra6gnyQWwv0jmsXAIEPjdx9e/TvLUAVMCiF+kREJI2CHOmvAArNrB/wDlABXB9k42Z2JnDQ3Y+Y2dnARcD9rS1WPjr6zny21WO35qaxEJGQaTH03b3ezKYBS4EcYKG7V5vZvcBKd19sZp8G/gycCfx3M5vr7hcARcCj0Q94TyMyp7++iV2JSBPa8ksSYOu8K9JUiZzqghzp4+5LgCUJy+6Je7yCyLRP4rj/AkrbWKOIiKSJrsgVEQkRhb6ISIgo9EVEQkShLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaLQFxEJEYW+iEiIKPRFREJEoS8iEiIKfRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCRGFvohIiCj0RURCRKEvIhIigb4YXUQya2vu9W3cwr601CGnPoW+ZETbQ0xEWkPTOyIiIaLQFxEJEYW+iEiIBAp9MxtjZhvNrMbMZiZZf4mZrTKzejObmLBuspltiv6ZnK7CRUQkdS2GvpnlAPOBy4Fi4DozK07otg24GXgyYexZwGxgGDAUmG1mZ7a9bBERaY0gR/pDgRp33+LuR4FKYHx8B3ff6u5rgA8Txl4GvODue939feAFYEwa6hYRkVYwd2++Q2S6Zoy73xpt3wgMc/dpSfr+CnjG3Z+Ktu8Cct39e9H2d4FD7v7jhHFTgakA+fn5QyorK2Pr6urqyMvLa/UTbC+qKzWN6trxRuaKiVPX9TzyjmzPdBmNpL2ucwembVOnzGssS3RUXSNGjHjN3cta6hfkPH1Lsqz53xQpjnX3BcACgLKyMi8vL4+tq6qqIr6dLVRXahrVNWd8k307UlX/uZRvnJ3pMhpJe13Xpe/irFPmNZYlsq2uINM7tUDvuHYBEPQQpC1jRUQkzYKE/gqg0Mz6mVkXoAJYHHD7S4FLzezM6Ae4l0aXiYhIBrQY+u5eD0wjEtYbgN+7e7WZ3Wtm4wDM7NNmVgtMAh41s+ro2L3AvxH5xbECuDe6TEREMiDQvXfcfQmwJGHZPXGPVxCZukk2diGwsA01iohImuiKXBGREFHoi4iEiEJfRCREFPoiIiGi0BcRCRGFvohIiCj0RURCRKEvIhIiCn0RkRBR6IuIhIhCX0QkRBT6IiIhotAXEQkRhb6ISIgo9EVEQkShLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaLQFxEJEYW+iEiIKPRFREKkU6YLkFPYnB7B+/afC3PGt18t0qy+M59t9dit865IYyWSaQp9kRDYmnt9G0bvS1sdknmBpnfMbIyZbTSzGjObmWR9VzP7XXT9P8ysb3R5XzM7ZGZvRP/8Ir3li4hIKlo80jezHGA+MBqoBVaY2WJ3Xx/XbQrwvrt/0swqgB8C10bXbXb3gWmuW0REWiHIkf5QoMbdt7j7UaASSJycHQ88EX38FDDSzCx9ZYqISDoECf1ewNtx7drosqR93L2eyCRgz+i6fmb2upn9bzP7XBvrFRGRNjB3b76D2STgMne/Ndq+ERjq7rfH9amO9qmNtjcTeYdQB+S5+x4zGwL8BbjA3T9I2MdUYCpAfn7+kMrKyti6uro68vLy2vxE0011ATveCNy1rut55B3Z3o7FtI7qCuDck2dn9dpPTUfVNWLEiNfcvaylfkHO3qkFese1C4DEV2NDn1oz6wT0APZ65DfKEQB3fy36y+BfgZXxg919AbAAoKyszMvLy2PrqqqqiG9nC9VFSqdgVvWfS/nG2e1YTOuorgCuO/nsHb32U5NtdQWZ3lkBFJpZPzPrAlQAixP6LAYmRx9PBJa7u5vZOdEPgjGz84FCYEt6ShcRkVS1eKTv7vVmNg1YCuQAC9292szuBVa6+2LgceA3ZlYD7CXyiwHgEuBeM6sHjgNfdfe97fFERESkZYEuznL3JcCShGX3xD0+DExKMu6PwB/bWKOIiKSJ7r0jIhIiCn0RkRBR6IuIhIhCX0QkRBT6IiIhotAXEQkRhb6ISIgo9EVEQkShLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaLQFxEJEYW+iEiIKPRFREIk0JeoyEfYnB6ZrkCyXeJrpP/c4N+PPGdfy32kQ+lIX0QkRBT6IiIhotAXEQkRhb6ISIjog1wRaTd9Zz7b6rFb512RxkqkgY70RURCRKEvIhIiCn0RkRAJNKdvZmOAh4Ac4DF3n5ewvivwa2AIsAe41t23Rtd9G5gCHAemu/vStFUvEfEXz6Ry4YxIO9uae30bRuvCrvbQ4pG+meUA84HLgWLgOjMrTug2BXjf3T8J/BT4YXRsMVABXACMAX4e3Z6IiGRAkOmdoUCNu29x96NAJZB4KDkeeCL6+ClgpJlZdHmlux9x97eAmuj2REQkA4JM7/QC3o5r1wLDmurj7vVmtg/oGV3+94SxvVpd7UeV7n8j0liq/y/ipzZ1z58mBQl9S7LMA/YJMhYzmwpMjTbrzGxj3Oqzgd0B6uxoWVrXN1VXSlRX6rK1tri65iaLnozpqJ9XnyCdgoR+LdA7rl0AbG+iT62ZdQJ6AHsDjsXdFwALku3czFa6e1mAOjuU6kqN6kpNttYF2Vub6gomyJz+CqDQzPqZWRciH8wuTuizGJgcfTwRWO7uHl1eYWZdzawfUAi8mp7SRUQkVS0e6Ufn6KcBS4mcsrnQ3avN7F5gpbsvBh4HfmNmNUSO8CuiY6vN7PfAeqAe+Ia7H2+n5yIiIi0IdJ6+uy8BliQsuyfu8WFgUhNjvw98vw01Jp32yQKqKzWqKzXZWhdkb22qKwCLzMKIiEgY6DYMIiIhckqFvpndZWZuZmdnuhYAM/uRmf1fM1tjZn82szMyXM8YM9toZjVmNjOTtTQws95m9pKZbTCzajO7I9M1xTOzHDN73cyeyXQtDczsDDN7Kvra2mBmn8l0TQBm9s3ov+E6M/utmeVmqI6FZvaema2LW3aWmb1gZpuif5+ZJXVlVUbAKRT6ZtYbGA1sy3QtcV4AStx9APAm8O1MFRLwdhmZUA/c6e5FwHDgG1lSV4M7gA2ZLiLBQ8Dz7v4p4EKyoD4z6wVMB8rcvYTISR0VGSrnV0Ru6xJvJvCiuxcCL0bbHe1XNK4razKiwSkT+kTu6fMtklzclSnu/p/uXh9t/p3IdQiZEuR2GR3O3Xe4+6ro4/1EAiwrrso2swLgCuCxTNfSwMz+BbiEyBlxuPtRd/9/ma0qphPQLXotzsdIcs1NR3D3/0PkLMF48beCeQKY0KFFkbyuLMsI4BQJfTMbB7zj7qszXUszbgGey+D+k90uIyvCtYGZ9QUGAf/IbCUxDxI5kPgw04XEOR/YBfyv6LTTY2Z2eqaLcvd3gB8Teae9A9jn7v+Z2apOku/uOyByoAF8PMP1JJPpjACyKPTNbFl0rjDxz3hgFnBPS9vIQF0NfWYRmcZYlIkaG8pIsixr3hWZWR7wR+B/uPsHWVDPlcB77v5apmtJ0AkYDPy7uw8CDpCZqYqTROfIxwP9gPOA083sS5mt6tSRJRkBZNF35Lr7qGTLzayUyAttdeTGnRQAq8xsqLvvzFRdcfVNBq4ERnpmz38NdMuLTDCzzkQCf5G7/ynT9URdBIwzs7FALvAvZvYf7p7pIKsFat294d3QU2RB6AOjgLfcfReAmf0J+CzwHxmt6oR3zexcd99hZucC72W6oAZZlBFAFh3pN8Xd17r7x929r7v3JfKfYnBHBH5Lol8u8z+Bce5+MMPlBLldRoeL3mL7cWCDuz+Q6XoauPu33b0g+pqqIHLrkEwHPtHX9dtm1j+6aCSRK9ozbRsw3Mw+Fv03HUkWfMAcJ/5WMJOBpzNYS0yWZQRwCoR+lnsE6A68YGZvmNkvMlVI9MOihttlbAB+7+7VmaonzkXAjcAXoj+jN6JH19K024FFZrYGGAj8IMP1EH3n8RSwClhLJDsycqWpmf0W+BvQ38xqzWwKMA8YbWabiJzlN6+5bXRgXVmTEQ10Ra6ISIjoSF9EJEQU+iIiIaLQFxEJEYW+iEiIKPRFREJEoS8iEiIKfRGREFHoi4iEyP8HqmlCIY5XYI4AAAAASUVORK5CYII=\n",
            "text/plain": "<Figure size 432x288 with 1 Axes>"
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "text/plain": "0.31268925739647135"
          },
          "execution_count": 25,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# sample 10000\nimport scipy as sp\nimport matplotlib.pyplot as plt\nimport scipy.stats as stats\n\nN = 10000\nn = 100\np = 1/2\nsample1 = stats.binom(n, p).rvs(N)\nsample2 = stats.norm(loc=n*p, scale=sp.sqrt(n*p*(1-p))).rvs(N)\n\n\nbins = sp.linspace(n*p-20,n*p+20, 20+1)\nhist1, bins = sp.histogram(sample1, bins, density=True)\nhist2, bins = sp.histogram(sample2, bins, density=True)\n\neps = 1e-9\n\nhist1 = (hist1 + eps)\nhist1 = hist1 / (sp.diff(bins) * hist1.sum())\nhist2 = (hist2 + eps)\nhist2 = hist2 / (sp.diff(bins) * hist2.sum())\n\nplt.bar(bins[:-1], hist1, label=\"$bin(100, 1/2)$\")\nplt.bar(bins[:-1], hist2, label=\"$N(np, np(1-p))$\")\nplt.legend()\nplt.grid(True)\nplt.show()\n\nD_KL(hist1, hist2) ",
      "execution_count": 96,
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X90VeWd7/H31/AjtqGhRc2o0QIFKQgKkoJ3bL1J8Qfae8UWGEO5ohUX09XSa6cz917aWxVcdUZ6q7Rr6UxLC16ljLE6q2NWRbnT4pmlrkoB+WXCsBqQOwS0XgKlHiFAwvf+cTbx5OScnH1OTpKT7M9rrSz2fvbz7P09T8L37POcvfdj7o6IiETDef0dgIiI9B0lfRGRCFHSFxGJECV9EZEIUdIXEYkQJX0RkQhR0hcRiRAlfRGRCFHSFxGJkCH9HUCqCy64wEePHl2QfX3wwQd89KMfLci+Ck2x5aeYY4Pijk+x5WegxLZt27Yj7n5h1kbuXlQ/06dP90J55ZVXCravQlNs+Snm2NyLOz7Flp+BEhuw1UPkWA3viIhEiJK+iEiEKOmLiERIqC9yzWw28COgBPiZuz+Ssn048DQwHWgB7nD3A2Y2FPgZcE1wrKfd/e8KGL+IhHDmzBnKysrYs2dPf4eSVnl5uWILqbS0lMrKSoYOHZpX+6xJ38xKgCeAG4FmYIuZ1bt7Y1K1xcAxdx9nZrXASuAOYD4w3N2nmNlHgEYze8bdD+QVrYjkpbm5mYqKCiorKzGz/g6ni/fff58RI0b0dxhpFVNs7k5LSwvNzc2MGTMmr32EGd6ZATS5+353Pw3UAXNS6swBngqWnwdmWeIvy4GPmtkQ4HzgNPCnvCIVkby1trZSXl5elAlfwjMzRo0aRWtra977CJP0LwUOJq03B2Vp67h7G3AcGEXiDeAD4B3g34EfuPvRvKMVkbwp4Q8OPf09hhnTT3eE1DkWM9WZAbQDlwAfB141s1+7+/5Ojc2WAEsAKioqiMViIcLKLh6PF2xfhabY8lPMsUHxxldeXk57ezvvv/9+f4eSlmLLTWtrK7FYLK+/tzBJvxm4LGm9EjicoU5zMJRTDhwFvgy87O5ngPfM7HWgCuiU9N19NbAaoKqqyqurq3N6EZnEYjEKta9CU2z5KebYoHjj27NnDyUlJUUzNp2qmMbNUxVjbKWlpUybNi2vv7cwSX8LMN7MxgCHgFoSyTxZPXAX8FtgHrDJ3d3M/h34vJn9HPgIcC3ww5wiFClmy8s7r09YAcuDr7yWH+/7eESyyJr03b3NzJYCG0lcsrnW3RvM7CESt/3WA2uAdWbWROIMvzZo/gTwJPAWiSGgJ919Vy+8DhHJwehlLxZ0fwce+UKoerFYjDVr1rBu3bpO5Zs3b+bVV19lxYoVGduePHmS2bNns2nTJkpKSrjnnnv41a9+xUUXXcRbb73VUe/ll1/mvvvuo729nXvvvZdly5aF2pbq3P4vuOACGhsbO237y7/8SxYtWsTll1/OokWLePfddznvvPNYsmQJ9913H6dPn+aGG25g06ZNDBlSXI84C3VzlrtvcPcr3P1T7v5wUPZAkPBx91Z3n+/u49x9xrkxe3ePB+VXuvskd/9fvfdSRKTY7dixg2nTpnUpnzlzZrcJH2Dt2rV86UtfoqSkBIC7776bl19+uVOd9vZ2vv71r/PSSy/R2NjIM88805Gwu9uWTrr9n7N582auvfZahgwZwqOPPsqePXt44403eOKJJ2hsbGTYsGHMmjWLZ599ttvX1B90R66I9JmdO3dy6NAhZs6cydixYzu+hFy0aBGvvfYaAF/84hf57ne/y+c+9zn+7M/+jF//+tcArF+/njlzPrxa/Prrr+cTn/hEp/3/7ne/Y9y4cYwdO5Zhw4ZRW1vLCy+8kHVbOun2D4nvR6644gpKSkq4+OKLueaaawAYMWIEEydO5NChQwDcfvvtrF+/Ps+e6j1K+iLLy9P/SMHt2LGDESNGsHnzZn784x9z//33A9DY2MiUKVMAeOuttxg5ciSvvvoqf//3f8/69es5ffo0+/fvJ9tj1w8dOsRll3143UllZWVHEu5uWy5eeuklZs+e3aX8wIEDbN++nZkzZwIwefJktmzZkvP+e5uSvoj0iba2NlpaWvjOd74DwNSpUzly5Aitra2cOXOG8vJyTpw4wfHjx/mrv/qrjjYjR47kyJEjjBw5MusxEk8Y7uzcde3dbcvFxo0buyT9eDzO3Llz+eEPf8jHPvYxAEpKShg2bFjRXe6ppC8ifaKxsZFx48YxbNgwAN58802uvvpqGhoa+PSnPw1AQ0MD06dP7xi337VrF5MnT+b8888PdRdqZWUlBw9+eC9pc3Mzl1xySdZtYZ04cYI//vGPndqdOXOGuXPnsnDhQr70pS91qn/q1ClKS0tzOkZvU9IXkT6xc+dO3n77bU6dOkU8HmfFihV885vfZPfu3Vx55ZVAYmhn6tSpHW127drFVVddxcc//nHa29uzJv7PfOYz/P73v+ftt9/m9OnT1NXVcdttt2XdNmvWrFBDPa+88go1NTUd6+7O4sWLmThxIt/61rc61W1paeHCCy/M+8FovaW4riUSkT4R9hLLQtq5cycLFy7kz//8zzl58iT3338/1157Lc899xxXXXUVALt37+4YE4fEm8DkyZMBuOmmm3jttde44YYbAFiwYAGxWIwjR45QWVnJihUrWLx4MY8//jg333wz7e3t3HPPPR1vKEOGDEm77ezZszQ1NXX50jbd/rdv3868efM66rz++uusW7eOKVOmdLxZ/e3f/i233norr7zyCrfeemvvdWielPRFpE/84Ac/AOB73/tep/JHH320Y9z7scce67Rt//4Pb95funQpjz32WEfSf+aZZ9Ie59Zbb82YbNNta2xsZO7cuZx//vmdys/tP/mO3GuuuYZVq1Z11PnsZz+b9rsCgH/8x3/k7/6u+J4kr6QvIgPCtGnTqKmpob29vWPMvxAmT57c5c0mkzfffDNUvdOnT3P77bczYcKEnoTWK5T0RWTAuOeee/o7hFCGDRvGokWL+juMtPRFrohIhCjpi4hEiJK+iEiEKOmLiESIkr6ISIQo6YuIRIiSvohIhOg6fZEoKvSjo0NODfmTn/yEr371qzQ2NjJx4kQAJk6cyHPPPcenPvWpTjNj9ZfUGbq+9rWvsXHjxi4zdBVK6ixbvT3rVqgzfTObbWZ7zazJzLrML2Zmw83s2WD7ZjMbHZQvNLMdST9nzWxqansRiYZdu3YxdepUXnwxMV3jqVOn+MMf/sDll1/eZWas/pIax8KFCzPOoFUIqbNs9fasW1mTvpmVkJjr9hZgErDAzCalVFsMHHP3ccAqYCWAu69396nuPhW4Ezjg7jsK+QJEZODYvXs3y5Yt60j6DQ0NTJw4ETPrMjNWphm0amtrueOOO5g5cyaf/OQnO/aVSa77SY3juuuuSzuDVj4yHTN1lq3enHUrzJn+DKDJ3fe7+2mgDpiTUmcO8FSw/Dwwy7rOTrAASP+EJJGIGr3sxbQ/g1VjYyO33XYb7733HsePH2f37t1MmTIl7cxY6WbQgsTTOseOHcvmzZtZv3591rl1c9lP2Bm68pUp9tRZtnpz1q0wA0aXAgeT1puBmZnquHubmR0HRgFHkurcQdc3CxGJiIMHDzJq1CjOP/98brzxRjZu3NjxvPyWlpZOM2NlmkHr5MmTHDlyhAcffBCASZMmcezYsYzHzHU/YWfoSnXDDTfw7rvvdil/+OGHOz41dBd78ixbI0aM6LJeSGGSfrr5xFKfJdptHTObCZxw97TfgpjZEmAJQEVFRcdkyT0Vj8cLtq9CU2z56ZXYJmQ4UwxznJS28eGXEDtXFqL9X09py3DoEMfOQXl5Oe3t7R2PMC5sGiHUlIBvvPEGEydO5P3336e6uppf/OIXvPvuu9x0000MGzaMkydPduxn27ZtXH311Zw4cQKALVu2MG7cODZv3szYsWM5c+YMZ86c4bXXXuPKK6/MePxc99PW1tYpDoD29nbi8Thnz57NeJxf/vKXWftm27Zt3cZ+btrITOvJWltbicVief1/CJP0m4HLktYrgcMZ6jSb2RCgHDiatL2WboZ23H01sBqgqqrKq6urQ4SVXSwWo1D7KjTFlp9eiW15hg+gC0JckZLSNjZhBdV7Hwzd/u4MQzkHFlZnP3YO9uzZQ0lJScHPGs8Js999+/Yxbdo0RowYwS233MK3vvUtTpw4wcyZMykpKeHs2bMMHTqU0tJS3n77baqqqjr2u3fvXv7iL/6CnTt3cujQIYYOHUp7ezsrV67k+9//fke9WbNm8fTTT3PppZcC5Lyfyy+/vFMckEjaZWVlnHfeeT3qv6ampoyxt7S0cNFFF3V8d5C6nqq0tJRp06bl9f8hTNLfAow3szHAIRIJ/MspdeqBu4DfAvOATR7MLGBm5wHzgetzikxEek/ISywLaffu3cydOxeA4cOHM2XKFLZv387IkSN5//33O82MlWkGraeffpqFCxdSXV3Nn/70J77zne9w3XXXAaSdASuf/aTO0PWVr3yF119/vcsMXbk6N3NYumOmzrLVm7NuZU36wRj9UmAjUAKsdfcGM3sI2Oru9cAaYJ2ZNZE4w69N2sX1QLO770/dt4hER+rVKC+88EKn9eSZsTLNoLVz505++tOfsnLlyi77TzcDVj77SZ2h68knnyzIJ6Tujpk6y1ZvzroV6sp/d98AbEgpeyBpuZXE2Xy6tjHg2vxDFJEoCDMz1r59+xg/fnzabbnMgNXdfnprhq5Mx0ydZau3Z93SHbkiUjSyzYx16NChghwn2356Y4auTMdMnWWrt2fdUtKXgS/TIwX6YdxapNjpgWsiIhGipC8SEcEFdTLA9fT3qKQvEgGlpaUcP35ciX+Ac3daWlo67iHIh8b0RSKgsrKSnTt3Eo/H+zuUtFpbW3uUyHpTscVWWlpKZWVl3u2V9EUiYOjQocTjcaqqqvo7lLRisRjTpk3r7zDSKubY8qHhHRGRCFHSFxGJEA3viPSjA6Wpj7E6R/cYSO/Qmb6ISIQo6YuIRIiSvohIhCjpi4hEiJK+iEiEKOmLiESIkr6ISISESvpmNtvM9ppZk5ktS7N9uJk9G2zfbGajk7ZdZWa/NbMGM9ttZsXzEAsRkYjJmvTNrAR4ArgFmAQsMLNJKdUWA8fcfRywClgZtB0C/Bz4qrtfCVQDZwoWvYiI5CTMmf4MoMnd97v7aaAOmJNSZw7wVLD8PDDLzAy4Cdjl7jsB3L3F3dsLE7qIiOQqTNK/FDiYtN4clKWt4+5tJO4hHwVcAbiZbTSzN83sv/c8ZBERyZdlm1TBzOYDN7v7vcH6ncAMd/9GUp2GoE5zsL6PxCeErwBfBz4DnAB+A3zX3X+TcowlwBKAioqK6XV1dQV5cfF4nLKysoLsq9AUW37SxvbOjvSVL54abqc9aZ/SNj78EspOHc67fU7HztGA+70WiYESW01NzTZ3z/rs7DAPXGsGLktarwQOZ6jTHIzjlwNHg/J/dfcjAGa2AbiGRPLv4O6rgdUAVVVVXl1dHSKs7GKxGIXaV6EptvykjW156mhjYEHIh5b1pH1K29iEFVTvfTDv9jkdO0cD7vdaJAZbbGGGd7YA481sjJkNA2qB+pQ69cBdwfI8YJMnPkJsBK4ys48Ebwb/EWjMKUIRESmYrGf67t5mZktJJPASYK27N5jZQ8BWd68H1gDrzKyJxBl+bdD2mJk9RuKNw4EN7v5iL70WERHJItTz9N19A7AhpeyBpOVWYH6Gtj8ncdmmiIj0M92RKyISIUr6IiIRoqQvIhIhSvoiIhGipC8iEiFK+iIiEaKkLyISIUr6IiIRoqQvIhIhSvoiIhGipC8iEiGhnr0jIsVp9LL0zy888MgX+jgSGSh0pi8iEiFK+iIiEaKkLyISIUr6IiIRoqQvIhIhoZK+mc02s71m1mRmy9JsH25mzwbbN5vZ6KB8tJmdNLMdwc+PCxu+iIjkIuslm2ZWAjwB3Ag0A1vMrN7dkyc4Xwwcc/dxZlYLrATuCLbtc/epBY5bRETyEOZMfwbQ5O773f00UAfMSakzB3gqWH4emGVmVrgwRUSkEMIk/UuBg0nrzUFZ2jru3gYcB0YF28aY2XYz+1cz+1wP4xURkR4wd+++gtl84GZ3vzdYvxOY4e7fSKrTENRpDtb3kfiEEAfK3L3FzKYD/wxc6e5/SjnGEmAJQEVFxfS6urqCvLh4PE5ZWVlB9lVoii0/aWN7Z0f6yheHHFXsSfuUtvHhl1B26nDe7XM6NrD70PG05VMuLe9SNuB+r0VioMRWU1Ozzd2rsrUJ8xiGZuCypPVK4HCGOs1mNgQoB4564h3lFIC7bwveDK4AtiY3dvfVwGqAqqoqr66uDhFWdrFYjELtq9AUW5LlXRPUh9s6J7W0sS1PHW0MLEifELseowftU9rGJqygeu+DebfP6djA3Zkew7CwukuZ/ubyM9hiCzO8swUYb2ZjzGwYUAvUp9SpB+4KlucBm9zdzezC4ItgzGwsMB7Yn1OEIiJSMFnP9N29zcyWAhuBEmCtuzeY2UPAVnevB9YA68ysCThK4o0B4HrgITNrA9qBr7r70d54ISIikl2op2y6+wZgQ0rZA0nLrcD8NO3+CfinHsYoIiIFojtyRUQiRElfRCRClPRFRCJESV9EJEI0XaLIAHag9MsZtoS8R0EiR2f6IiIRoqQvIhIhSvoiIhGipC8iEiFK+iIiEaKkLyISIUr6IiIRoqQvIhIhSvoiIhGipC8iEiFK+iIiEaKkLyISIUr6IiIREirpm9lsM9trZk1mtizN9uFm9mywfbOZjU7ZfrmZxc3sbwoTtoiI5CNr0jezEuAJ4BZgErDAzCalVFsMHHP3ccAqYGXK9lXASz0PV0REeiLMmf4MoMnd97v7aaAOmJNSZw7wVLD8PDDLzAzAzG4H9gMNhQlZRETyZe7efQWzecBsd783WL8TmOnuS5PqvBXUaQ7W9wEzgZPAr4Ebgb8B4u7+gzTHWAIsAaioqJheV1dXgJcG8XicsrKyguyr0BRbknd2ZN528dROq2ljy9Q+pW3Oxw/TPqVtfPgllJ06nHf7nI6dY3v9zeVnoMRWU1Ozzd2rsrUJM3OWpSlLfafIVGcFsMrd48GJf1ruvhpYDVBVVeXV1dUhwsouFotRqH0VmmJLsjz1g2OSBZ1ngEobW6b2C0LOHtWT9iltYxNWUL33wbzb53TsHNvrby4/gy22MEm/Gbgsab0SOJyhTrOZDQHKgaMkzvbnmdn3gZHAWTNrdffHc4pSREQKIkzS3wKMN7MxwCGgFkidmLMeuAv4LTAP2OSJcaPPnatgZstJDO8o4YuI9JOsSd/d28xsKbARKAHWunuDmT0EbHX3emANsM7Mmkic4df2ZtAiIpKfMGf6uPsGYENK2QNJy63A/Cz7WJ5HfCIiUkC6I1dEJEKU9EVEIkRJX0QkQpT0RUQiRElfRCRClPRFRCJESV9EJEJCXacvIoPP6GUvZtx24JEv9GEk0pd0pi8iEiFK+iIiEaKkLyISIUr6IiIRoqQvIhIhSvoiIhGipC8iEiFK+iIiEaKkLyISIaHuyDWz2cCPSEyX+DN3fyRl+3DgaWA60ALc4e4HzGwGsPpcNWC5u/+yUMFLEVlenqH8eN/GISLdynqmb2YlwBPALcAkYIGZTUqpthg45u7jgFXAyqD8LaDK3acCs4GfmJke/SAi0k/CJOAZQJO77wcwszpgDtCYVGcOsDxYfh543MzM3U8k1SkFvMcRi0hBHCj9cjdb9QltsDL37vOwmc0DZrv7vcH6ncBMd1+aVOetoE5zsL4vqHPEzGYCa4FPAnemG94xsyXAEoCKiorpdXV1BXlx8XicsrKyguyr0AZdbO/sSF9+8dT826Zpnza2nhy7p+1T2saHX0LZqcN5t8/p2Dm279J3OfR7bxt0/x/6SHJsNTU129y9KlubMGf6lqYs9Z0iYx133wxcaWYTgafM7CV3b+1U0X01wdh/VVWVV1dXhwgru1gsRqH2VWiDLrblc9KXLwhxxpipbZr2aWPrybF72j6lbWzCCqr3Pph3+5yOnWP7Ln2XQ7/3tkH3/6GP5BNbmKt3moHLktYrgcOZ6gRj9uXA0eQK7r4H+ACYnFOEIiJSMGGS/hZgvJmNMbNhQC1Qn1KnHrgrWJ4HbHJ3D9oMATCzTwITgAMFiVxERHKWdXjH3dvMbCmwkcQlm2vdvcHMHgK2uns9sAZYZ2ZNJM7wa4PmnwWWmdkZ4CzwNXc/0hsvREREsgt1+aS7bwA2pJQ9kLTcCsxP024dsK6HMYqISIHojlwRkQhR0hcRiRAlfRGRCFHSFxGJECV9EZEIUdIXEYkQJX0RkQhR0hcRiRAlfRGRCFHSFxGJECV9EZEIUdIXEYkQJX0RkQhR0hcRiRAlfRGRCFHSFxGJECV9EZEICZX0zWy2me01syYzW5Zm+3AzezbYvtnMRgflN5rZNjPbHfz7+cKGLyIiucia9M2sBHgCuAWYBCwws0kp1RYDx9x9HLAKWBmUHwH+s7tPITFxuqZOFBHpR2HO9GcATe6+391PA3XAnJQ6c4CnguXngVlmZu6+3d0PB+UNQKmZDS9E4CIikrswSf9S4GDSenNQlraOu7cBx4FRKXXmAtvd/VR+oYqISE+Zu3dfwWw+cLO73xus3wnMcPdvJNVpCOo0B+v7gjotwfqVQD1wk7vvS3OMJcASgIqKiul1dXWFeG3E43HKysoKsq9CG3SxvbMjffnFU/Nvm6Z92th6cuyetk9pGx9+CWWnDufdPqdj59i+S9/l0O+9bdD9f+gjybHV1NRsc/eqbG2GhNhvM3BZ0nolcDhDnWYzGwKUA0cBzKwS+CWwKF3CB3D31cBqgKqqKq+urg4RVnaxWIxC7avQBl1sy1NH/AILjuffNk37tLH15Ng9bZ/SNjZhBdV7H8y7fU7HzrF9l77Lod9726D7/9BH8oktzPDOFmC8mY0xs2FALYmz9mT1JL6oBZgHbHJ3N7ORwIvAt9399ZwiExGRgst6pu/ubWa2FNgIlABr3b3BzB4Ctrp7PbAGWGdmTSTO8GuD5kuBccD9ZnZ/UHaTu79X6BciIn1r9LIX05YfeOQLfRyJ5CLM8A7uvgHYkFL2QNJyKzA/TbvvAd/rYYwiIlIguiNXRCRClPRFRCJESV9EJEKU9EVEIkRJX0QkQpT0RUQiRElfRCRCQl2nLxGwvPzD5QkrPrxFf3nf3o4vA8eB0i9n2KK/mWKmM30RkQhR0hcRiRAlfRGRCFHSFxGJECV9EZEIUdIXEYkQJX0RkQhR0hcRiRAlfRGRCAmV9M1stpntNbMmM1uWZvtwM3s22L7ZzEYH5aPM7BUzi5vZ44UNXUREcpU16ZtZCfAEcAswCVhgZpNSqi0Gjrn7OGAVsDIobwXuB/6mYBGLiEjewpzpzwCa3H2/u58G6oA5KXXmAE8Fy88Ds8zM3P0Dd3+NRPIXEZF+FibpXwocTFpvDsrS1nH3NhJPXBpViABFRKRwzN27r2A2H7jZ3e8N1u8EZrj7N5LqNAR1moP1fUGdlmD9bqDK3ZdmOMYSYAlARUXF9Lq6up6+LgDi8ThlZWUF2VehFV1s7+zoWIwPv4SyU4cTKxdPzbl9J2HaZ2qbpn3afuvJsXvaPqVtzn3Xh7F36bsc+r2nx86m6P4/JBkosdXU1Gxz96psbcI8WrkZuCxpvRI4nKFOs5kNAcqBo2GCBnD31cBqgKqqKq+urg7btFuxWIxC7avQii625R+O2MUmrKB674OJlQUhH5O7PHXEj/DtM7VN0z5tv/Xk2D1tn9I2577rw9i79F0O/d7TY2dTdP8fkgy22MIM72wBxpvZGDMbBtQC9Sl16oG7guV5wCbP9hFCRET6XNYzfXdvM7OlwEagBFjr7g1m9hCw1d3rgTXAOjNrInGGX3uuvZkdAD4GDDOz24Gb3L2x8C9FRESyCTVzlrtvADaklD2QtNwKzM/QdnQP4hORQWr0shc7lv96Sht3B+sHHvlCf4UUCbojV0QkQpT0RUQiRBOji0i/SJ5YPXbeCg6UBlc9aWL1XqUzfRGRCNGZ/iCR/KVYKn0xJiLnKOkPEskflbvSx2URSdDwjohIhCjpi4hEiJK+iEiEaExfRAakTBcv6MKF7ulMX0QkQpT0RUQiRMM7xWR5eYZyXXIpIoWhM30RkQjRmb6IDEiZb0jUJ+Pu6ExfRCRClPRFRCIk1PCOmc0GfkRiusSfufsjKduHA08D04EW4A53PxBs+zawGGgH/qu7byxY9CIieYjyAwqzJn0zKwGeAG4EmoEtZlafMs/tYuCYu48zs1pgJXCHmU0iMV/ulcAlwK/N7Ap3by/0CykKuvpGRIpcmDP9GUCTu+8HMLM6YA6QnPTnAMuD5eeBx83MgvI6dz8FvB1MnD4D+G1hwhcRyV1Pn0o7kO8GDpP0LwUOJq03AzMz1XH3NjM7DowKyt9IaXtp3tH2tkxn6qCzdRHpMJCvHDJ3776C2XzgZne/N1i/E5jh7t9IqtMQ1GkO1veROKN/CPitu/88KF8DbHD3f0o5xhJgSbA6AdhbgNcGcAFwpED7KjTFlp9ijg2KOz7Flp+BEtsn3f3CbA3CnOk3A5clrVcChzPUaTazIUA5cDRkW9x9NbA6RCw5MbOt7l5V6P0WgmLLTzHHBsUdn2LLz2CLLcwlm1uA8WY2xsyGkfhitj6lTj1wV7A8D9jkiY8Q9UCtmQ03szHAeOB3uQQoIiKFk/VMPxijXwpsJHHJ5lp3bzCzh4Ct7l4PrAHWBV/UHiXxxkBQ7xckvvRtA74+aK/cEREZAEJdp+/uG4ANKWUPJC23AvMztH0YeLgHMfZEwYeMCkix5aeYY4Pijk+x5WdQxZb1i1wRERk89BgGEZEIGRRJ38xKzex3ZrbTzBrMbEVQPsbMNpvZ783s2eCL6GKJ7X+b2dtmtiP4mdrXsSXFWGJm283sV8F6v/dbN7EVU78dMLPdQRxbg7JPmNm/BH33L2b28SKKbbmZHUrqu1v7KbYRLJF5AAADj0lEQVSRZva8mf2bme0xs/9QLP3WTXz93ndmNiHp+DvM7E9m9s1c+25QJH3gFPB5d78amArMNrNrSTwOYpW7jweOkXhcRLHEBvDf3H1q8LOjH2I75z5gT9J6MfTbOamxQfH0G0BNEMe5y+aWAb8J+u43wXp/SY0NEr/Xc323IWPL3vUj4GV3/zRwNYnfbzH1W7r4oJ/7zt33njs+ieecnQB+SY59NyiSvifEg9WhwY8DnyfxWAiAp4Dbiyi2omBmlcAXgJ8F60YR9Fu62AaIOST6DPqx74qVmX0MuJ7EFX+4+2l3/yNF0m/dxFdsZgH73P3/kmPfDYqkDx3DADuA94B/AfYBf3T3tqBKvz0CIjU2d98cbHrYzHaZ2SpLPKm0P/wQ+O/A2WB9FEXSb3SN7Zxi6DdIvHn/HzPbFtxVDlDh7u8ABP9eVESxASwN+m5tPw2hjAX+H/BkMGz3MzP7KMXTb5nig/7vu2S1wDPBck59N2iSvru3Bx97Kkk8AmJiump9G1Vw0JTYzGwy8G3g08BngE8A/6Ov4zKz/wS85+7bkovTVO3zfssQGxRBvyW5zt2vAW4Bvm5m1/djLKnSxfYPwKdIDDO+AzzaD3ENAa4B/sHdpwEf0L9DOakyxVcMfQdA8B3bbcBz+bQfNEn/nOCjWAy4FhhpicdCQIZHQPSlpNhmu/s7wdDPKeBJEm9Ufe064DYzOwDUkRjW+SHF0W9dYjOznxdJvwHg7oeDf98jMbY6A/iDmV0MEPz7XrHE5u5/CE5AzgI/pX/6rhloTvq0+zyJJFsU/ZYpviLpu3NuAd509z8E6zn13aBI+mZ2oZmNDJbPB24g8eXLKyQeCwGJx0S8UCSx/VvSL8lIjMG91dexufu33b3S3UeT+Li4yd0XUgT9liG2/1IM/RYc/6NmNuLcMnBTEEvyI0n6628ubWzn+i7wRfrnb+5d4KCZTQiKZpG4Y7/f+w0yx1cMfZdkAR8O7UCufefuA/4HuArYDuwi8ct4ICgfS+JZP00kPgoNL6LYNgG7g7KfA2X93IfVwK+Kpd+6ia0o+i3oo53BTwPwP4PyUSSuoPh98O8niii2dUHf7QoSxcX91HdTga1BHP8MfLwY+i1LfMXSdx8hMTtheVJZTn2nO3JFRCJkUAzviIhIOEr6IiIRoqQvIhIhSvoiIhGipC8iEiFK+iIiEaKkLyISIUr6IiIR8v8BYXvDNfgm2qQAAAAASUVORK5CYII=\n",
            "text/plain": "<Figure size 432x288 with 1 Axes>"
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "text/plain": "0.005019163960424132"
          },
          "execution_count": 96,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "n*p, n*p*(1-p)",
      "execution_count": 81,
      "outputs": [
        {
          "data": {
            "text/plain": "(50.0, 25.0)"
          },
          "execution_count": 81,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# sample 10000\n# n = 10000\nimport scipy as sp\nimport matplotlib.pyplot as plt\nimport scipy.stats as stats\n\nN = 10000\nn = 10000\np = 1/2\nsample1 = stats.binom(n, p).rvs(N)\nsample2 = stats.norm(loc=n*p, scale=sp.sqrt(n*p*(1-p))).rvs(N)\n\n\nbins = sp.linspace(n*p-200,n*p+200, 200+1)\nhist1, hist2 = create_hists(sample1, sample2, bins, eps=1e-9)\n\nplt.bar(bins[:-1], hist1, label=\"$Bin(10000, 1/2)$\")\nplt.bar(bins[:-1], hist2, label=\"$N(np, np(1-p))$\")\nplt.legend()\nplt.grid(True)\nplt.savefig(\"Bin(10000,0.5)_hist.png\")\nplt.show()\n\nD_KL(hist1, hist2) ",
      "execution_count": 28,
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt0lPW97/H3l2CIW26KitZ4ABdoQfDSUK/HC4VdKFpRkVW8u4VFu472eLRnnUIvgnZh625ru8+urtqq1SqKbvfeRzaCtN0Ej7YKKGCupYbLKaFYBS0SMQTC9/zx/BImk5lkZjIzySSf11pZ+T2/+T2/5/s8mcx3nuvP3B0REZF+3R2AiIj0DEoIIiICKCGIiEighCAiIoASgoiIBEoIIiICKCGIiEighCAiIoASgoiIBP27O4B0HH/88T5y5MiM5v3kk0845phjshtQFiiu9PXU2BRXehRXejKN6+23397t7iek1NjdC+anrKzMM1VeXp7xvLmkuNLXU2NTXOlRXOnJNC7gLU/xM1aHjEREBNA5BBERCZQQREQEKLCTyiKSfQcPHqS+vp7Gxsa8LnfIkCHU1tbmdZmpKNS4SkpKKC0t5aijjsp4GUoIIn1cfX09gwYNYuTIkZhZ3pa7b98+Bg0alLflpaoQ43J39uzZQ319PaNGjcp4GTpkJNLHNTY2MmzYsLwmA8kuM2PYsGFd3stTQhARJYNeIBt/QyUEEREBlBBERCRQQhDpzKIh3R2BSF7oKiMRaWPk/Jez2t/2H1yRUrtHH32URYsWMXz4cBoaGhg/fjwvvPACb731FqtWreK+++5LOu+nn37KtGnTWL16NUVFRdx+++0sX76cE088kaqqqjZtX3nlFe666y6am5uZO3cu8+fPb1N/8OBB5s2b164+vn1nr8XrKCaAr371q9xwww0sWrSI9957j379+jFv3jzuuusumpqamDZtGq+++ir9++fuY1t7CCLSI1RUVPDAAw+wadMm/vSnP1FVVUVFRQUXXXRRh8kA4IknnuDaa6+lqKgIgNtuu41XXnmlXbvm5mbuuOMOVq5cSU1NDc899xw1NTVt6tevX5+wPrZ9R30lkyymFmvXrmX06NH8+Mc/pra2ljfffJOHH36YmpoaiouLueyyy3j++edT2ZQZU0IQkR6hsrKSc889F4C6ujrcndNPP51Zs2bx+uuvA3DNNdfwne98h0suuYSTTjqJ3/3udwAsWbKEGTNmtPZ16aWXctxxx7Vbxrp16xg9ejSnnXYaxcXFzJ49m5deeint+o76SiZZTAC1tbWcfvrpnHLKKXzuc58DYNCgQYwdO5adO3cCcOWVV7JkyZJ0N2talBBEpEeorq7mlltuYezYsZSVlfGrX/2KwYMHU1VVxYQJEwCoqqpi6NChvPbaazzyyCMsWbKEpqYmtm7dSiqPxt+5cyennnpq63RpaSk7d+5Mu76jvjKxcuVKpk2b1qZu+/btbNy4kfPPPx+AcePGsX79+oz6T5USgkg6dII5J3bs2MGJJ55IRUUFtbW1PPzww3zve9+jsbGRgwcPMmTIEPbv38/evXu5++67ATh06BBDhw5l9+7dDB06NKXlRE+DbsvM0q7vqK9MrFq1qk1CaGhoYObMmfz0pz9l8ODBABQVFVFcXMy+ffsyWkYqlBBEpNtVVFQwbty41umzzz6b999/n+rq6tb66upqysrKWs8TVFRUMH78eI4++uiU79AtLS1lx44drdP19fV85jOfSbu+o77StX//fv72t7+1znvw4EFmzpzJjTfeyLXXXtum7YEDBygpKUl7GalSQhCRbldZWcnYsWOB6Jv3U089xZQpU6isrOSss84CosNF55xzTus8FRUVnHXWWRx77LE0NzenlBQ+//nP8+6777Jt2zaamppYunQpV111Vdr1HfUFMHny5JQPH5WXlzNp0qTWdZ8zZw5jx47lnnvuadNuz549nHDCCV16eF1ndNmpSBaNnP9yypdZ9lTdEX9lZSWvvvoqL7/8MmbGBRdcwI9+9CO++93vth5Dr6ysbC1DlCDGjx8PwBe/+EVef/11pkyZAsD111/PmjVr2L17N6Wlpdx3333MmTOH/v3787Of/YypU6fS3NzM7bffzplnngnQWn/w4EHmzp3brj6+fbK+Dh8+TF1dXbsTyMliWrlyJddddx0Av//973n66aeZMGFCa/J74IEHmD59Oq+99hrTp0/P1Z8gkurQaj3hR0No5k9Pjcu9G2JbODhxOU55ebmP+ObyPASUns62V01NTX4CifPxxx9nra8NGzb4TTfdlJW+uhpXZWWl33333Sm3P/fcc72pqanTdl/+8pf9j3/8Y4dtEv0t0RCaInmkE83d7txzz2XSpEk0Nzd3dyiMHz+ehx56KOX2GzZs6PQwUFNTE1dccQVnnHFGV8PrkA4ZiUivcPvtt3d3CDlTXFzMDTfckPPlaA9BREQAJQQREQmUEEREBFBCEMlczMnkbD8hVKQ7KCGIiAighCAiIkFKCcHMppnZZjOrM7N2I0CY2QAzez68vtbMRsa8tiDUbzazqTH1d5tZtZlVmdlzZpa7B3SI5JPuS5AC1el9CGZWBDwM/D1QD6w3s2XuHjsSxBzgI3cfbWazgQeBr5jZOGA2cCbwGeB3ZnY6cBLw34Fx7v6pmb0Q2j2ZvVUTKWCLhsCivd237Kz2l9p6PProo3zta1+jpqam9blGY8eOZeXKlQwfPrzNiGjdJZ2R2bKhqamJKVOmsHr16nbTuRg5LZU9hPOAOnff6u5NwFJgRlybGcBTofwiMNmi58DOAJa6+wF33wbUhf4gSkZHm1l/4O+Av3RtVUSkkFVUVHDOOefw8svRCfoDBw7w17/+lREjRrQbEa27pDoyW7YUFxczefLk1pHS4qezLZWEcAqwI2a6PtQlbOPuh4C9wLBk87r7TuBHwJ+BXcBed/9NJisgIr1DZWUl8+fPb00I1dXVjB07FjNrNyJaspHTZs+ezVe+8hXOP/98RowY0dpXMun2k+rIbJlItsyrr766zUhp8dPZlMo+R6IRH+JHhkjWJmG9mR1LtPcwCvgb8C9mdpO7P9Nu4WbzgHkAw4cPZ82aNSmE3F5DQ0PG8+aS4kpf3mM74z5oWV6i8hn38Y3Dh2hoaOAbE5pZ0y+mTTaW2UWdba8hQ4a0GXRlUFaWekSyAV2am5vbvFZdXc2kSZNYtGgR9fX1rFu3js9+9rPs2bOHLVu2MGzYsNb2FRUVTJw4kRUrVrBs2TKefPJJzj//fDZu3MiVV17JY489xhtvvMGCBQu49NJLk8aWqJ+JEycm7OeCCy5oFwdE2/fw4cNdHrgmWewjRoxg3bp1rdurZTrR8hobG7v0v5FKQqgHTo2ZLqX94Z2WNvXhENAQ4MMO5p0CbHP3DwDM7N+Ai4B2CcHdfwH8AmDixIl++eWXpxBye2vWrCHTeXNJcaUv77EtmgHX701eXjSD2xqf5clpx/Dj1z9he8nCI22yscwu6mx71dbWMmhQttPAEcn63rdvX+trO3bs4Pjjj+fEE09k6tSp/OEPf+Ddd9+lrKyMAwcOcOyxx7a23b9/P/v27WP+/PkUFRVx1FFHccIJJ9C/f38+/PBDFi9eTElJCRMnTuTjjz9Ouvxk/TQ1NSXsJz6OFgMHDqRfv35JlzNlyhTee++9dvWLFy9u3dv49NNPO4x9wIAB7N+/v3UQnQEDBiTctiUlJa3jUmcilYSwHhhjZqOAnUQnf+OfsrQMuBV4A7gOWO3ubmbLgGfN7CGik8pjgHXAYeACM/s74FNgMvBWxmsh0hN154nhAlNRUdE6bvL06dNZsmQJu3bt4uqrr243IlqykdOqqqoYM2ZM64hiGzZs4Oyzz066zGT91NTUJOwnnZHZYrUchupIZ7HHj5SWq5HTOj2HEM4J3AmsAmqBF9y92szuN7OrQrPHgWFmVgfcA8wP81YDLwA1wCvAHe7e7O5riU4+bwAqQxy/yOqaiUjBqKysbE0Il112Ga+99lprkogfES3ZyGnvvPMOf/7zn2lsbOSTTz5h4cKFreMvQ/tRzJL1U1VVlbCfdEZmS1dHscePlJbLkdNSum7J3VcAK+Lq7o0pNwKzksy7GFicoH4hsDCdYEUkD7phr6ayspKZM2cC0eGQCRMmsHHjRoYOHQq0HREt2chpv/71r7nxxhu5/PLL+fjjj/nWt77FxRdfDJBwFLNk/Tz22GNJ+0l1ZLZ0vfPOO0mXWV5e3maktPjpbNJ4CCLS7eKvmnnppZfaTN9555089NBDTJkypd3gM1u3bgWiD9Vf/vKXPPjgg+36r6mpYebMmRx99NGtdcn6qaqq4oknnkjYT2wcAM8991yqq9ihjmJ/9tln+f73v590Opv06AoR6fFSGRFty5YtjBkzJuFr6Yxitm3btqT95GpktmSxNzU1cfXVV7eOlBY/nW3aQxCRgtDZiGix5we6YvPmzfTrl/y7ci5GZksWe3FxMbfcckvS6WzTHoJIIrl8HpGedSQ9lBKCiIgASggiPUZ3DrLjHv/wASk02fgbKiFI31AAh2m2l8Tf75kfJSUl7NmzR0mhgLk7e/bs6fLNajqpLNLHlZaWUl9fzwcffJDX5TY2NubkbtuuKtS4SkpKKC0t7dIylBBE+rijjjqKUaNG5X25a9as6dJzd3KlL8elQ0YiWdBdh3tEskkJQfqebjif0J0njEVSpYQgkkXaU5BCpoQg0tMVwBVS0jsoIYiICKCEICIigRKCiIgASgjS2+h4u0jGlBBERARQQhDJC12OKoVACUFERAAlBBERCZQQRLqRHmkhPYkSgkigD2fp65QQREQEUEIQEZFACUF6LR0CEkmPEoL0OkoEIplRQhAREUAJQUREAiUE6dN0eEnkCCUEEREBlBCkl+t1ewB6vLfkkBKCiIgASggiIhIoIYgUqF53OEy6nRKCSNAdg9ho4BzpSZQQREQESDEhmNk0M9tsZnVmNj/B6wPM7Pnw+lozGxnz2oJQv9nMpsbUDzWzF83sj2ZWa2YXZmOFREQkM50mBDMrAh4GvgSMA643s3FxzeYAH7n7aOAnwINh3nHAbOBMYBrwSOgP4J+AV9z9s8DZQG3XV0f6qr54PF2HmyTbUtlDOA+oc/et7t4ELAVmxLWZATwVyi8Ck83MQv1Sdz/g7tuAOuA8MxsMXAo8DuDuTe7+t66vjoiIZCqVhHAKsCNmuj7UJWzj7oeAvcCwDuY9DfgA+JWZbTSzx8zsmIzWQITMvy0X6rfsvrhHJLln7t5xA7NZwFR3nxumbwbOc/evx7SpDm3qw/QWoj2L+4E33P2ZUP84sAL4f8CbwMXuvtbM/gn42N2/m2D584B5AMOHDy9bunRpRiva0NDAwIEDM5o3lxRX+hLGtmsTnHwO7NpE5eFRTDhlCJU79wIw4ZQhR16PbZut8q5NUVyDRjNwX130Wkt9jpZZuXMvE/ptO1Kf7vbqARRXejKNa9KkSW+7+8SUGrt7hz/AhcCqmOkFwIK4NquAC0O5P7AbsPi2Le2Ak4DtMfWXAC93FktZWZlnqry8PON5c0lxpS9hbAsHt/4e8c3l7u4+4pvLW8utr+eivHCw+8LBUVyhnOtljvjm8rb1Heipf0vFlZ5M4wLe8k4+W1t+UjlktB4YY2ajzKyY6CTxsrg2y4BbQ/k6YHUIZBkwO1yFNAoYA6xz9/eAHWZ2RphnMlCTUgYTSUOhHhIS6Q79O2vg7ofM7E6ib/dFwBPuXm1m9xNlnmVEJ4efNrM64EOipEFo9wLRh/0h4A53bw5dfx1YEpLMVuAfsrxuIiKShk4TAoC7ryA69h9bd29MuRGYlWTexcDiBPWbgNSOa4mISM7pTmXpdfrCYaK+sI6Sf0oIUlB0uWV72iaSLUoI0uNl+oGnb9Ei6VFCEBERQAlBpHfREJvSBUoIUrB07Fwku5QQREQEUEIQEZFACUFERAAlBBERCZQQREQEUEIQKXjbS27QFVeSFUoIIiICKCGIiEighCAiIoASgoiIBEoIIiICKCGIiEighCAiIoASghQwDYDTMd2bIOlSQhDppWITppKDpEIJQUREACUEKQA6NCSSH0oIUlCUHERyRwlBpBdQopRsUEIQERFACUFERAIlBBERAZQQREQkUEIQERFACUFERAIlBBERAZQQpIfSs3dE8k8JQaQP0I1rkgolBBERAZQQREQkUEIQERFACUF6KB3zFsm/lBKCmU0zs81mVmdm8xO8PsDMng+vrzWzkTGvLQj1m81satx8RWa20cyWd3VFRESkazpNCGZWBDwMfAkYB1xvZuPims0BPnL30cBPgAfDvOOA2cCZwDTgkdBfi7uA2q6uhPQSi4Z0dwQifVoqewjnAXXuvtXdm4ClwIy4NjOAp0L5RWCymVmoX+ruB9x9G1AX+sPMSoErgMe6vhoikqqk93goIfd55u4dNzC7Dpjm7nPD9M3A+e5+Z0ybqtCmPkxvAc4HFgFvuvszof5xYKW7v2hmLwLfBwYB/9Pdr0yy/HnAPIDhw4eXLV26NKMVbWhoYODAgRnNm0uKK8auTXDyOZ2WGwaNjmJL1GbXpmg6UTmFvjMqh+U0DBrNwH11eV1mKuvZ0NAQxZXOcvJA7/30ZBrXpEmT3nb3iSk1dvcOf4BZwGMx0zcD/xzXphoojZneAgwjOtR0U0z948BM4ErgkVB3ObC8szjcnbKyMs9UeXl5xvPmkuKKsXBwSuXW2BK1WTg4eTnN5aRcDsspLy/P+zITluPatsaVznLyQO/99GQaF/CWp/D56u4pHTKqB06NmS4F/pKsjZn1B4YAH3Yw78XAVWa2negQ1BfM7JkUYhERkRxJJSGsB8aY2SgzKyY6Sbwsrs0y4NZQvg5YHTLTMmB2uAppFDAGWOfuC9y91N1Hhv5Wu/tNWVgfERHJUP/OGrj7ITO7E1gFFAFPuHu1md1PtCuyjOhQ0NNmVke0ZzA7zFttZi8ANcAh4A53b87RuoiISBd0mhAA3H0FsCKu7t6YciPRuYZE8y4GFnfQ9xpgTSpxiIhI7uhOZRERAZQQREQkUEIQERFACUFERAIlBBERAZQQREQkUEIQERFACUFERAIlBBERAZQQRKQTScdPkF5HCUG6lwZlEekxlBBERARQQhARkUAJQUREACUEEREJlBBERARQQpBc01VEIgVDCUHyT0lCpEdSQhAREUAJQUREAiUEEREBlBBEJJGY8zzbS27oxkAkn5QQREQEUEKQfNGVRQVPTz3t/ZQQREQEUEIQEZFACUFERAAlBMmSlI4v6zyCSI+mhCAiKdHlp72fEoKIiABKCCIiEighiEjqwnkg3ZPQOykhiIgIoIQgIiKBEoKIiABKCCKSodbzCLq/pNdQQpDc0QeFSEFJKSGY2TQz22xmdWY2P8HrA8zs+fD6WjMbGfPaglC/2cymhrpTzazczGrNrNrM7srWColI7ukmtd6p04RgZkXAw8CXgHHA9WY2Lq7ZHOAjdx8N/AR4MMw7DpgNnAlMAx4J/R0CvuHuY4ELgDsS9CkiInmUyh7CeUCdu2919yZgKTAjrs0M4KlQfhGYbGYW6pe6+wF33wbUAee5+y533wDg7vuAWuCUrq+OiIhkKpWEcAqwI2a6nvYf3q1t3P0QsBcYlsq84fDSucDa1MMWEZFsM3fvuIHZLGCqu88N0zcTfcv/ekyb6tCmPkxvIdqzuB94w92fCfWPAyvc/V/D9EDgVWCxu/9bkuXPA+YBDB8+vGzp0qUZrWhDQwMDBw7MaN5c6i1xVe7cy4RT4k4i79oU/T75nPblk8850ibNcsOg0VFsidrkaJkdlsNyGgaNZuC+urwuM5X1bGhoiOLK8jIrD4+K/uaxbdPQW977+ZJpXJMmTXrb3Sem1NjdO/wBLgRWxUwvABbEtVkFXBjK/YHdgMW3jWt3VJi+p7MYWn7Kyso8U+Xl5RnPm0u9Ja4R31zevnLh4OgnUTm2TZrl1tgStcnRMjssh+WUl5fnfZkJy3FtW+PK8jJb/+axbdPQW977+ZJpXMBbnuJnbCqHjNYDY8xslJkVE50kXhbXZhlwayhfB6wOgSwDZoerkEYBY4B14fzC40Ctuz+UUuaSHq3NVSe63FSkIPXvrIG7HzKzO4m+zRcBT7h7tZndT5R5lhF9uD9tZnXAh0RJg9DuBaCG6MqiO9y92cz+K3AzUGlmYT+Ub7n7imyvoIiIpKbThAAQPqhXxNXdG1NuBGYlmXcxsDiu7nWiQ0oi0ouMnP8y239wRXeHIRnSncqSMT0CWTqkQ4cFRwlBRDKiu5V7HyUEEREBlBBEJIu011DYlBAkY/rnF+ldlBBERARQQhARkUAJQUREACUESYWuJ5cu0j0rhUEJQUREACUEEREJlBBERARQQpBkWgZDEZE+QwlB0qMTzJIB3cRYGJQQRCS/9KWix1JCEJEeR5epdg8lBBERAZQQRKQbaU+gZ1FCEBERQAlBYk/w6WSfSJ+mhNCX6ANfRDqghCAiIoASgoh0I92w1rMoIfRFOnQkIgkoIYiICKCEICIigRJCH9HuBiAdNpICo5vYck8JoRfSP44UokTvW72X80sJQUR6nJarj2KvQkp6RZL2drNGCaG3SPJPoW9YUih0CWr3U0Lo5RJ90xLptbS30CVKCCIiAighiEhPl+ABjDoUmhtKCIVGu8QikiNKCL2QzhdIb5fKezzTvYjKnXszmq83UELoQdJ+A2tvQSR79P+khNAj6I0oklsx/2M6/5BcSgnBzKaZ2WYzqzOz+QleH2Bmz4fX15rZyJjXFoT6zWY2NdU+e7tO35R6A4t0SbLDSrk83FToOk0IZlYEPAx8CRgHXG9m4+KazQE+cvfRwE+AB8O844DZwJnANOARMytKsc8+oS8frxTpdto7byOVPYTzgDp33+ruTcBSYEZcmxnAU6H8IjDZzCzUL3X3A+6+DagL/aXSp4iI5FEqCeEUYEfMdH2oS9jG3Q8Be4FhHcybSp85lXCXMIVvC23mSzJAfbLdzZb6kfNfTrjb2pVdXBHJjgn9th35n0v2mZDkfz9Zm0I5BGXu3nEDs1nAVHefG6ZvBs5z96/HtKkOberD9BaivYD7gTfc/ZlQ/ziwgigRddhnTN/zgHlh8gxgc4brejywO8N5c0lxpa+nxqa40qO40pNpXCPc/YRUGvZPoU09cGrMdCnwlyRt6s2sPzAE+LCTeTvrEwB3/wXwixTi7JCZveXuE7vaT7YprvT11NgUV3oUV3ryEVcqh4zWA2PMbJSZFROdJF4W12YZcGsoXwes9mjXYxkwO1yFNAoYA6xLsU8REcmjTvcQ3P2Qmd0JrAKKgCfcvdrM7gfecvdlwOPA02ZWR7RnMDvMW21mLwA1wCHgDndvBkjUZ/ZXT0REUpXKISPcfQXRsf/Yuntjyo3ArCTzLgYWp9JnjnX5sFOOKK709dTYFFd6FFd6ch5XpyeVRUSkb9CjK0REBOgFCSHc+bzRzJaH6clmtsHMNpnZ62Y2OtSn/XiNPMV1m5l9EOo3mdncmD5uNbN3w8+tyZbVxbi+EOKqMrOnwlViWOR/h+1SYWaf6yFxXW5me2O2170xfWT9cShmtt3MKsOy3gp1x5nZb8P6/9bMjg31edtmacaVt22WJK5ZZlZtZofNbGJc+7w82iaduMxspJl9GrO9fh7zWlnopy78rS1Hsf3QzP4Y3kf/bmZDY9rnbpu5e0H/APcAzwLLw/SfgLGh/N+AJ2PKPw/l2cDzoTwOeAcYAIwCtgBFeYzrNuBnCeY/Dtgafh8bysdmMy6iLwQ7gNPDa/cDc0J5OrASMOACYG0Pievylm0aN39R+NudBhSHv+m4LMS1HTg+ru4fgfmhPB94MN/bLM248rbNksQ1lugeojXAxJj6hP97PSCukUBVkn7WAReGv/FK4Es5eo99Eegfyg/G/C1zus0Keg/BzEqBK4DHYqodGBzKQzhyf0O6j9fIV1zJTAV+6+4fuvtHwG+JngeVzbiGAQfc/U9h+rfAzFCeAfzaI28CQ83s5B4QVzL5fBxK7HvpKeDqmPq8bLM040omL9vM3WvdPdENpd36aJsO4koo/C0Hu/sbHn06/5rOt3Gmsf3Go6c+ALxJdK8W5HibFXRCAH4K/C/gcEzdXGCFmdUDNwM/CPXpPl4jX3EBzAy7hi+aWcsNe/mIazdwVMzu8nUcuWEwn48dSScugAvN7B0zW2lmZ3YSb1c58Bsze9uiu+YBhrv7LoDw+8ROYshFbOnEBfnbZoniSqa7t1dHRll0CPNVM7skJt76LMeVSmy3E+2NtMSQs21WsAnBzK4E3nf3t+NeuhuY7u6lwK+Ah1pmSdCNd1Cfr7j+Axjp7mcBv+PIN7ycxxW+5cwGfmJm64B9RPeLdLT87o5rA9Gt+GcD/wz8n07i7aqL3f1zRE/mvcPMLu2gbV62WQZx5XOb9YbttQv4L+5+LuFQppkNzlFcHcZmZt8meu8vaalKEkNWYivYhABcDFxlZtuJdo++YGYvA2e7+9rQ5nngolBufYyGpf54jZzH5e573P1AqP8lUBYfbw7jeibs/l7i7ucB/xd4t5Pld2tc7v6xuzeE8gqiPYnjcxAXYRl/Cb/fB/6daNf8r+HwQcthhPdD83xts7Tiyuc2SxJXMt29vZK1PeDue0L5baJj86eHuEpjmubyPYZFFx9cCdwYviRBrrdZJidBetoP4aQZ0Y12uzlyMnIO8K+hfAdtTyq/EMpn0vYkzVaycFI5jbhOjml/DfBmKB8HbCM6CXlsKB+XzbhC+cTwewDwn8AXwvQVtD1Buq6HxHUSR+6fOQ/4c4ixf/jbjeLISbUzuxjPMcCgmPIfiI79/5C2J2//MZ/bLIO48rLNksUV8/oa2p68Tfi/1wPiOoHwGUB0knZny9+L6LE7F3DkpPL0HL3HphE94eGEuPY53WZd/ifuCT+0/SC5BqgMG2QNcFqoLwH+hegkzLqW+vDat4m+BWwmC1cNpBnX94HqUF8OfDZm/ttDvHXAP+Qorh8CtWHd/0dMGyMaxGhLiHtiD4nrzpjt9SZwUcxr04mu5toCfDsK6jsiAAAAi0lEQVQL8ZwWlvNOWOa3Q/0woiT1bvjd8mGRl22WQVx52WYdxHUN0TfYA8BfgVWd/e91Z1xEFzC0bK8NwJdj+poIVIW4fkZItDmIrY7onMCm8PPzfGwz3aksIiJAYZ9DEBGRLFJCEBERQAlBREQCJQQREQGUEEREJFBCEBERQAlBREQCJQQREQHg/wMw/wza6XZpngAAAABJRU5ErkJggg==\n",
            "text/plain": "<Figure size 432x288 with 1 Axes>"
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        },
        {
          "data": {
            "text/plain": "0.016068499410022707"
          },
          "execution_count": 28,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "bins",
      "execution_count": 111,
      "outputs": [
        {
          "data": {
            "text/plain": "array([4800.        , 4813.33333333, 4826.66666667, 4840.        ,\n       4853.33333333, 4866.66666667, 4880.        , 4893.33333333,\n       4906.66666667, 4920.        , 4933.33333333, 4946.66666667,\n       4960.        , 4973.33333333, 4986.66666667, 5000.        ,\n       5013.33333333, 5026.66666667, 5040.        , 5053.33333333,\n       5066.66666667, 5080.        , 5093.33333333, 5106.66666667,\n       5120.        , 5133.33333333, 5146.66666667, 5160.        ,\n       5173.33333333, 5186.66666667, 5200.        ])"
          },
          "execution_count": 111,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    }
  ],
  "metadata": {
    "hide_input": false,
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "language_info": {
      "name": "python",
      "version": "3.7.0",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "toc": {
      "nav_menu": {},
      "number_sections": true,
      "sideBar": true,
      "skip_h1_title": false,
      "base_numbering": 1,
      "title_cell": "Table of Contents",
      "title_sidebar": "Contents",
      "toc_cell": false,
      "toc_position": {},
      "toc_section_display": true,
      "toc_window_display": true
    },
    "gist": {
      "id": "",
      "data": {
        "description": "histogramからKL",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}