Skip to content

Instantly share code, notes, and snippets.

@Cartman0

Cartman0/histogramからKL.ipynb

Created May 16, 2019 20:41
Show Gist options
  • Save Cartman0/c94c19a1f450c4832e12bfa349b8c2e8 to your computer and use it in GitHub Desktop.
Save Cartman0/c94c19a1f450c4832e12bfa349b8c2e8 to your computer and use it in GitHub Desktop.
Download ZIP
histogramからKL情報量
Display the source blob
Display the rendered blob
Raw
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, 18+1)\nhist1, bins = sp.histogram(sample1, bins, density=True)\nhist1\nhist2, bins = sp.histogram(sample2, bins, density=True)\nhist2",
"execution_count": 45,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 45,
"data": {
"text/plain": "array([0. , 0. , 0. , 0.02, 0.02, 0.08, 0.18, 0.36, 0.36, 0.46, 0.2 ,\n 0.12, 0.14, 0.06, 0. , 0. , 0. , 0. ])"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# 0が無くなるように再構成\n\neps = 1e-9\neps\n\nhist1 = (hist1 + eps)\nhist1 = hist1 / (sp.diff(bins) * hist1.sum())\n\nhist2 = (hist2 + eps)\nhist2 = hist2 / (sp.diff(bins) * hist2.sum())\nhist2",
"execution_count": 58,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 58,
"data": {
"text/plain": "array([1.19999989e-08, 1.19999989e-08, 1.19999989e-08, 2.00000098e-02,\n 2.00000098e-02, 8.00000034e-02, 1.79999993e-01, 3.59999973e-01,\n 3.59999973e-01, 4.59999962e-01, 1.99999990e-01, 1.19999999e-01,\n 1.39999997e-01, 6.00000055e-02, 1.19999989e-08, 1.19999989e-08,\n 1.19999989e-08, 1.19999989e-08])"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true,
"scrolled": true
},
"cell_type": "code",
"source": "hist2.sum()",
"execution_count": 59,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 59,
"data": {
"text/plain": "1.9999999999999996"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true,
"scrolled": 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": 60,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": "<Figure size 432x288 with 1 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEr1JREFUeJzt3XuQXGWZx/HvE0IyC4mhilhRSMqJbhITE0QdLhaWDmhI2FhJYUEVASyoVYOXsN4Aw2pBwsZyvSxu1S6CqV0ktauEy0IZIYKiTHkp0VzkkpANm2JHGHBLiBYywQSCz/4xkziZTDI90z3dM+98P3/1Oec95zz9pvuXd073eTsyE0lSWcY0ugBJUu0Z7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCjW3UiSdPnpzNzc2NOv2g7N69m2OPPbbRZTSUfdDFfrAP9qt3P2zevPn5zHxtf+0aFu7Nzc1s2rSpUacflLa2NlpbWxtdRkPZB13sB/tgv3r3Q0T8ppJ2XpaRpAIZ7pJUIMNdkgrUsGvuksr1yiuv0NHRwZ49expdypCbNGkS27dvr/lxm5qamDp1KkcfffSg9jfcJdVcR0cHEydOpLm5mYhodDlD6sUXX2TixIk1PWZmsmvXLjo6Opg+ffqgjuFlGUk1t2fPHo4//vjig32oRATHH398VX/5GO6ShoTBXp1q+89wl6QCec1d0pBrXnFvTY/X/o+Lanq8arW2tvK1r32NlpaWitrfcccdrFy5ku3bt/OrX/2q4v0GwnDXiFaL0GhvunDgO81aBSuXDGyflS8M/Dwq0ty5c7nrrru47LLLhuwcXpaRVJzdu3ezaNEi3vrWtzJ37lxuu+02AK677jpOOeUU5s6dy7Jly8hMoGvk/elPf5p3v/vdzJ49m40bN/KBD3yAGTNm8IUvfAGA9vZ23vzmN3PJJZdw0kkncd555/HSSy8dcu4f/OAHvPOd7+Ttb387559/Pp2dnYe0mT17NrNmzRrCHjDcJRXovvvu44QTTuCRRx5h69atLFy4EIDly5ezceNGtm7dyp/+9CfuueeeA/uMGzeOn/zkJ3z0ox9lyZIl3HDDDWzdupVbbrmFXbt2AbBjxw6WLVvGo48+ymte8xq+8Y1vHHTe559/ntWrV/PAAw+wZcsWWlpauP766+v3xHsw3CUVZ968eTzwwAN87nOf46c//SmTJk0C4MEHH+S0005j3rx5/PjHP2bbtm0H9lm8ePGBfd/ylrfw+te/nvHjx/PGN76Rp59+GoBp06ZxxhlnAHDxxRfzs5/97KDzPvTQQzz++OOcccYZnHzyyaxdu5bf/Kaieb5qzmvukoozc+ZMNm/ezIYNG7j66qs5++yzueqqq/j4xz/Opk2bmDZtGitXrjzoe+Tjx48HYMyYMQce71/et28fcOjXE3svZybz58/n1ltvHaqnVjFH7pKK8+yzz3LMMcdw8cUXc8UVV7Bly5YDQT558mQ6Ozu58847B3zcp556il/84hcA3HrrrbzrXe86aPvpp5/Oz3/+c3bu3AnASy+9xBNPPFHlsxkcR+6Shly9v7r42GOPceWVVzJmzBiOPvpobrzxRo477jg+8pGPMG/ePJqbmznllFMGfNzZs2ezdu1aLrvsMmbMmMHHPvYxXn311QPbX/va13LLLbewdOlS9u7dC8Dq1auZOXPmQce5++67ufzyy3nuuedYtGgRJ598Mvfff391T7oXw11ScRYsWMCCBQsOWb969WpWr159yPq2trYDj1tbWw/68Y3929rb2xkzZgw33XTTQfu++OKLB+1/1llnsXHjxiPWd+6553Luuef2/0Sq4GUZSSqQ4S5JFWhubmbr1q2NLqNihrskFchwl6QCGe6SVCDDXZIK5FchJQ29lZNqfLzhNcPmQKf8vfLKK/ne977HuHHjeNOb3sS3vvUtjjvuuJrW5Mhdkups/vz5bN26lUcffZSZM2fypS99qebnMNwlFWe4T/l79tlnM3Zs14WT008/nY6Ojpr3geEuqTgjacrfm2++mXPOOafGPWC4SyrQSJny94tf/CJjx47loosuqunzBz9QlVSgkTDl79q1a7nnnnv40Y9+dMhxasGRu6TiDPcpf++77z6+/OUvs379eo455pgB11EJR+6Shl6dv7o43Kf8Xb58OXv37mX+/PlA138KvWebrJbhLqk4w33K3/0j+6FU0WWZiFgYETsiYmdErDhCu/MiIiOism/yS5KGRL/hHhFHATcA5wBzgKURMaePdhOBvwN+WesiJanRSpzy91RgZ2Y+mZkvA+uAJX20+wfgK8CePrZJGmX23yCkwam2/yoJ9xOBp3ssd3SvOyAi3gZMy8x7kDTqNTU1sWvXLgN+kDKTXbt20dTUNOhjVPKBal9fwDzwLxYRY4CvA5f2e6CIZcAygClTphz0IcRI0NnZOeJqrrXh1gefnbev6mO0jVk14H06x59A26wB7jeM+q0WjvRaiAiOPfbYAzf/lCwzh+R76q+++iq7d+8+4k1QR1JJuHcA03osTwWe7bE8EZgLtHU/wdcB6yNicWZu6nmgzFwDrAFoaWnJnp9IjwRtbW2MtJprbbj1waUr7q36GO1N1w54n7ZZq2jdMcD9lg6vmQyrNdxeC40yXPuhkssyG4EZETE9IsYBFwDr92/MzBcyc3JmNmdmM/AQcEiwS5Lqp99wz8x9wHLgfmA7cHtmbouI6yJi8VAXKEkauIpuYsrMDcCGXuuuOUzb1urLkiRVw7llJKlAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SClTRrJDSkFk5qard2wf/K2RS0Ry5S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAlUU7hGxMCJ2RMTOiFjRx/aPRsRjEfFwRPwsIubUvlRJUqX6DfeIOAq4ATgHmAMs7SO8v5OZ8zLzZOArwPU1r1SSVLFKRu6nAjsz88nMfBlYByzp2SAz/9hj8Vgga1eiJGmgxlbQ5kTg6R7LHcBpvRtFxCeAzwDjgLNqUp0kaVAi88iD7Ig4H1iQmR/uXv4gcGpmXn6Y9hd2t7+kj23LgGUAU6ZMece6deuqLL++Ojs7mTBhQqPLaKia98FvH67dseqoc/wJTNj77ID2eezP06s657wTJ1W1f635fuhS734488wzN2dmS3/tKhm5dwDTeixPBY70ql4H3NjXhsxcA6wBaGlpydbW1gpOP3y0tbUx0mqutZr3wcol/bcZhtpmraJ1x7UD2ufSPd+p6pztF7VWtX+t+X7oMlz7oZJr7huBGRExPSLGARcA63s2iIgZPRYXAf9TuxIlSQPV78g9M/dFxHLgfuAo4ObM3BYR1wGbMnM9sDwi3ge8AvwBOOSSjCSpfiq5LENmbgA29Fp3TY/Hn6xxXZKkKniHqiQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKVFG4R8TCiNgRETsjYkUf2z8TEY9HxKMR8aOIeEPtS5UkVarfcI+Io4AbgHOAOcDSiJjTq9mvgZbMPAm4E/hKrQuVJFWukpH7qcDOzHwyM18G1gFLejbIzAcz86XuxYeAqbUtU5I0EJWE+4nA0z2WO7rXHc6HgO9XU5QkqTqRmUduEHE+sCAzP9y9/EHg1My8vI+2FwPLgfdk5t4+ti8DlgFMmTLlHevWrav+GdRRZ2cnEyZMaHQZDdWzDx575oWqjzdvzP9WfYxG6Bx/AhP2Plvfk77+5Pqerx++H7rUux/OPPPMzZnZ0l+7sRUcqwOY1mN5KnDIqzoi3gd8nsMEO0BmrgHWALS0tGRra2sFpx8+2traGGk111rPPrh0xb1VH6+96dqqj9EIbbNW0bqjzrUvrf4/01ry/dBluPZDJZdlNgIzImJ6RIwDLgDW92wQEW8Dvgkszszf1b5MSdJA9BvumbmPrkst9wPbgdszc1tEXBcRi7ubfRWYANwREQ9HxPrDHE6SVAeVXJYhMzcAG3qtu6bH4/fVuC5JUhW8Q1WSCmS4S1KBDHdJKpDhLkkFqugDVakv7U0XNroESYfhyF2SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIK5C8xSSVYOalO53mhPudR1Ry5S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgrkTUzSSFGvG5VUBEfuklQgw12SCmS4S1KBDHdJKlBF4R4RCyNiR0TsjIgVfWx/d0RsiYh9EXFe7cuUJA1Ev+EeEUcBNwDnAHOApRExp1ezp4BLge/UukBJ0sBV8lXIU4GdmfkkQESsA5YAj+9vkJnt3dv+PAQ1SpIGqJLLMicCT/dY7uheJ0kapioZuUcf63IwJ4uIZcAygClTptDW1jaYwzRMZ2fniKv5SB57ZuC/qjPlr+Bfvv1dAObNWlXrkkaMzvEn0DYan3+P139p74fBGq79UEm4dwDTeixPBZ4dzMkycw2wBqClpSVbW1sHc5iGaWtrY6TVfCSXrrh3wPt8dt4+/umxrpdNe9O1tS5pxGibtYrWHaPw+S/9y4CgtPfDYA3XfqjkssxGYEZETI+IccAFwPqhLUuSVI1+wz0z9wHLgfuB7cDtmbktIq6LiMUAEXFKRHQA5wPfjIhtQ1m0JOnIKpo4LDM3ABt6rbumx+ONdF2ukSQNA96hKkkFMtwlqUCGuyQVyHCXpAL5S0yjWHvThQPep23MqlH9/XZppHDkLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQNzFJqtzKSX95PGsVrFwywP0H/utfGhxH7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCeRPTCNa84t6q9m9vqlEhkoYdR+6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAnkTU6P0/EWbQfImJEmH48hdkgpkuEtSgQx3SSqQ4S5JBaoo3CNiYUTsiIidEbGij+3jI+K27u2/jIjmWhcqSapcv+EeEUcBNwDnAHOApRExp1ezDwF/yMy/Br4OfLnWhUqSKlfJyP1UYGdmPpmZLwPrgCW92iwB1nY/vhN4b0RE7cqUJA1EJeF+IvB0j+WO7nV9tsnMfcALwPG1KFCSNHCV3MTU1wg8B9GGiFgGLOte7IyIHRWcfziZDDzf6CIa69P2AWA/wKD6YFWRf9DX+7XwhkoaVRLuHcC0HstTgWcP06YjIsYCk4Df9z5QZq4B1lRS2HAUEZsys6XRdTSSfdDFfrAP9huu/VDJZZmNwIyImB4R44ALgPW92qwHLul+fB7w48w8ZOQuSaqPfkfumbkvIpYD9wNHATdn5raIuA7YlJnrgX8H/iMidtI1Yr9gKIuWJB1ZRROHZeYGYEOvddf0eLwHOL+2pQ1LI/aSUg3ZB13sB/tgv2HZD+HVE0kqj9MPSFKBDPdBiogrIiIjYnKja6m3iPhqRPx3RDwaEXdHxHGNrqle+puKYzSIiGkR8WBEbI+IbRHxyUbX1CgRcVRE/Doi7ml0Lb0Z7oMQEdOA+cBTja6lQX4IzM3Mk4AngKsbXE9dVDgVx2iwD/hsZs4GTgc+MUr7AeCTwPZGF9EXw31wvg5cRR83ao0GmfmD7juRAR6i696H0aCSqTiKl5m/zcwt3Y9fpCvcet+1XryImAosAv6t0bX0xXAfoIhYDDyTmY80upZh4m+B7ze6iDqpZCqOUaV7Bti3Ab9sbCUN8c90DfL+3OhC+uJvqPYhIh4AXtfHps8Dfw+cXd+K6u9IfZCZ3+1u83m6/kT/dj1ra6CKptkYLSJiAvBfwKcy84+NrqeeIuL9wO8yc3NEtDa6nr4Y7n3IzPf1tT4i5gHTgUe6J72cCmyJiFMz8//qWOKQO1wf7BcRlwDvB947iu5GrmQqjlEhIo6mK9i/nZl3NbqeBjgDWBwRfwM0Aa+JiP/MzIsbXNcBfs+9ChHRDrRk5qiaQCoiFgLXA+/JzOcaXU+9dM+b9ATwXuAZuqbmuDAztzW0sDrrns57LfD7zPxUo+tptO6R+xWZ+f5G19KT19w1GP8KTAR+GBEPR8RNjS6oHro/RN4/Fcd24PbRFuzdzgA+CJzV/e//cPcIVsOII3dJKpAjd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KB/h+8wNWwn+XxjQAAAABJRU5ErkJggg==\n"
},
"metadata": {
"needs_background": "light"
}
}
]
},
{
"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": 61,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 61,
"data": {
"text/plain": "0.12159790930303474"
},
"metadata": {}
}
]
},
{
"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\nbins = sp.linspace(-5,5, 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=\"sample 1\")\nplt.bar(bins[:-1], hist2, label=\"sample 2\")\nplt.legend()\nplt.grid(True)\nplt.show()\n\nD_KL(hist1, hist2) ",
"execution_count": 65,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": "<Figure size 432x288 with 1 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAG3xJREFUeJzt3X2QXHWd7/H3J4FkLgYDS6hZQ6IT3EnMkwYdAhSKA5IHDBLLgjJBrFDXdUCJ614NEq4WJDGUgC57y9oopu5GUntXIqLuHTESeeqruCKTIA95uINDHMkYb4kJhQwhiRO/94/uCZ1mwpye6YfpPp9XVSp9zvn9Tn9/Sddnzpxz+ncUEZiZWTqMqnYBZmZWOQ59M7MUceibmaWIQ9/MLEUc+mZmKeLQNzNLEYe+mVmKOPTNzFLEoW9mliInVLuAQhMmTIimpqZql1FWr7zyCm9605uqXUZFeKz1yWMdebZt2/aniDh9sHYjLvSbmprYunVrtcsoq0wmQ2tra7XLqAiPtT55rCOPpN8laefTO2ZmKeLQNzNLEYe+mVmKjLhz+mZW3/7yl7/Q09PDwYMHq11KIuPHj2fXrl3VLuOohoYGJk2axIknnjik/g59M6uonp4eTj75ZJqampBU7XIG9fLLL3PyySdXuwwAIoJ9+/bR09PDlClThrQPn94xs4o6ePAgp512Wk0E/kgjidNOO21YvyU59M2s4hz4QzfcfzuHvplZivicvplVVdPKH5d0f923Lirp/oajtbWVr33ta7S0tCRq/73vfY9Vq1axa9cuHn/88cT9iuHQt7pX6lDJN5ICxmrfrFmz+MEPfsA111xTtvfw6R0zS41XXnmFRYsW8a53vYtZs2bx3e9+F4A1a9Zw9tlnM2vWLNra2ogIIHukvnLlSi644AKmT59OR0cHH/nIR2hubuZLX/oSAN3d3bzjHe9g2bJlvPOd7+Tyyy/nwIEDr3vvn/70p5x33nm8+93v5oorrqC3t/d1baZPn860adPK+C/g0DezFLn//vuZOHEiTz31FNu3b2fhwoUALF++nI6ODrZv386rr77Kfffdd7TPmDFj+NnPfsa1117L4sWLWbduHdu3b+euu+5i3759AHR2dtLW1sbTTz/Nm9/8Zr7xjW8c875/+tOfWLt2LQ8++CBPPPEELS0t3HHHHZUbeB6f3rG6191wZRn3/lIZ922lNnv2bFasWMENN9zApZdeyvve9z4AHnnkEW6//XYOHDjA/v37mTlzJh/60IcA+OAHP3i078yZM3nLW94CwJlnnsmePXs45ZRTmDx5Mueffz4AV111FV//+tdZsWLF0fd97LHH2Llz59E2hw8f5rzzzqvYuPM59M0sNaZOncq2bdvYvHkzN954I/Pnz+cLX/gCn/70p9m6dSuTJ09m1apVx9wHP2bMGABGjRrF2LFjj64fNWoUfX19wOtvoyxcjgjmzZvH3XffXa6hJebTO2aWGnv37uWkk07iqquuYsWKFTzxxBNHA37ChAn09vZy7733Fr3f559/nl/+8pcA3H333bz3ve89Zvu5557LL37xC7q6ugA4cOAAzz777DBHMzQ+0jezqqrkHVDPPPMM119/PaNGjeLEE0/km9/8Jqeccgqf/OQnmT17Nk1NTZx99tlF73f69Ols3LiRa665hubmZj71qU8ds/3000/nrrvuYunSpRw6dAiAtWvXMnXq1GPa/fCHP+Qzn/kML7zwAosWLWLOnDls2bJl6AMegPqvUo8ULS0t4Yeo1I8RMdZV48u479fO6Y+IsVbIcMa6a9cupk+fXtqCymiwuXe6u7u59NJL2b59e8VqGujfUNK2iBj0xv5Ep3ckLZTUKalL0so3aHe5pJDUkrfuxly/TkkLkryfmZmVx6CndySNBtYB84AeoENSe0TsLGh3MvAPwK/y1s0AlgAzgYnAg5KmRsSR0g3BzKx6mpqaKnqUP1xJjvTnAl0RsTsiDgObgMUDtPsycDuQP/3bYmBTRByKiN8CXbn9mZlZFSS5kHsGsCdvuQc4J7+BpLOAyRFxn6QVBX0fK+h7RuEbSGoD2gAaGxvJZDKJiq9Vvb29dT/GfiNirNNWl2/feWMbEWOtkOGMdfz48bz88sulLaiMjhw5MuLqPXjw4JD//ZOE/kDzeB69+itpFPDPwNXF9j26ImI9sB6yF3Lr/WKYL/hV2KqBfjEtkaW+kFusXbt2jZiHkiQxkh6i0q+hoYGzzjprSH2ThH4PMDlveRKwN2/5ZGAWkMl9IeFvgXZJlyXoa2ZmFZQk9DuAZklTgN+TvTB79HvtEfESMKF/WVIGWBERWyW9CnxH0h1kL+Q2A4+Xrnwzq3mlvqV21ciZGqPYqZWvv/56fvSjHzFmzBje/va38+1vf5tTTjmlpDUNeiE3IvqA5cAWYBdwT0TskLQmdzT/Rn13APcAO4H7get8546Z2cDmzZvH9u3befrpp5k6dSpf+cpXSv4eie7Tj4jNETE1It4eEbfk1t0UEe0DtG2NiK15y7fk+k2LiJ+UrnQzs+KM9KmV58+fzwknZE/AnHvuufT09JT838Bz75hZatTS1MobNmzgkksuKfG/gEPfzFJk9uzZPPjgg9xwww38/Oc/Z/z47PWERx55hHPOOYfZs2fz8MMPs2PHjqN9BppaeezYsUenVgZeN7Xyo48+esz75k+tPGfOHDZu3Mjvfve749Z5yy23cMIJJ/Cxj32spOMHT7hmZilSC1Mrb9y4kfvuu4+HHnrodfspBR/pm1lqjPSple+//35uu+022tvbOemkk4quIwkf6ZtZdVXwFsuRPrXy8uXLOXToEPPmzQOyPyzuvPPOIY52YA59M0uNBQsWsGDB6yf7Xbt2LWvXrn3d+kwmc3QKhtbW1mO+hdw/DUJ3dzejRo0aMJzzp0q46KKL6OjoeMP6+n8TKCef3jEzSxGHvpnZMNTa1Mo+vWM2DE0rf3z09edn93F13vJgKvmYwJEmIspyZ0oaDPdphz7SN7OKamhoYN++fcMOrzSKCPbt20dDQ8OQ9+EjfTOrqEmTJtHT08MLL7xQ7VISOXjw4LBCttQaGhqYNGnSkPs79M2sok488USmTJlS7TISy2QyQ567fiTy6R0zsxRx6JuZpYhD38wsRRKFvqSFkjoldUlaOcD2ayU9I+lJSY9KmpFb3yTp1dz6JyWV9vvEZmZWlEEv5EoaDawD5pF95m2HpPaI2JnX7DsRcWeu/WXAHcDC3LbnImJOacs2Gxm6G44+OZTMqNV0N9xcRO+R81g/S48kR/pzga6I2B0Rh4FNwOL8BhHx57zFNwG+AdfMbARKEvpnAHvylnty644h6TpJzwG3A/+Qt2mKpF9L+j+S3jesas3MbFg02LfiJF0BLIiIv88tfxyYGxGfOU77K3Ptl0kaC4yLiH2S3gP8BzCz4DcDJLUBbQCNjY3v2bRp03DHNaL19vYybty4apdRESNirH94siJv0zt2IuMO7U3e4S21e9ZzRPy/VkitjPXCCy/cFhEtg7VL8uWsHmBy3vIk4I0+2ZuAbwJExCHgUO71ttxvAlOBrfkdImI9sB6gpaUl8qcvrUeZTIZ6H2O/ETHWVYsHb1MCmWmrae0s4pz+0to9pz8i/l8rpN7GmuT0TgfQLGmKpDHAEqA9v4Gk5rzFRcBvcutPz10IRtKZQDOwuxSFm5lZ8QY90o+IPknLgS3AaGBDROyQtAbYGhHtwHJJFwN/AV4EluW6XwCskdQHHAGujYj95RiImZkNLtHcOxGxGdhcsO6mvNefPU6/7wPfH06BZmZWOv5GrplZijj0zcxSxKFvZpYiDn0zsxRx6JuZpYhD38wsRRz6ZmYp4tA3M0sRh76ZWYok+kauWdWtGl/tCszqgo/0zcxSxKFvZpYiDn0zsxRx6JuZpYhD38wsRXz3jlmVNK388ZD7dt+6qISVWJokOtKXtFBSp6QuSSsH2H6tpGckPSnpUUkz8rbdmOvXKWlBKYs3M7PiDBr6uWfcrgMuAWYAS/NDPec7ETE7IuYAtwN35PrOIPtM3ZnAQuAb/c/MNTOzyktypD8X6IqI3RFxGNgELM5vEBF/zlt8ExC514uBTRFxKCJ+C3Tl9mdmZlWQ5Jz+GcCevOUe4JzCRpKuAz4HjAEuyuv7WEHfMwbo2wa0ATQ2NpLJZBKUVbt6e3vrfoz9SjbWaauHv48y6x07kUwRdX7+r31Dfq9qf378Ga5dSUJfA6yL162IWAesk3Ql8CVgWRF91wPrAVpaWqK1tTVBWbUrk8lQ72PsV7Kxrlo8eJsqy0xbTWvnzYnbX33wO0N+r+6PtQ65byn4M1y7kpze6QEm5y1PAva+QftNwIeH2NfMzMooSeh3AM2SpkgaQ/bCbHt+A0nNeYuLgN/kXrcDSySNlTQFaAYeH37ZZmY2FIOe3omIPknLgS3AaGBDROyQtAbYGhHtwHJJFwN/AV4ke2qHXLt7gJ1AH3BdRBwp01jMakp3w5XD6P1SyeqwdEn05ayI2AxsLlh3U97rz75B31uAW4ZaoJmZlY6nYTAzSxGHvplZijj0zcxSxKFvZpYiDn0zsxRx6JuZpYhD38wsRRz6ZmYp4tA3M0sRh76ZWYo49M3MUsShb2aWIg59M7MUceibmaWIQ9/MLEUShb6khZI6JXVJWjnA9s9J2inpaUkPSXpb3rYjkp7M/Wkv7GtmZpUz6ENUJI0G1gHzyD7ztkNSe0TszGv2a6AlIg5I+hRwO/DR3LZXI2JOies2M7MhSHKkPxfoiojdEXGY7IPPF+c3iIhHIuJAbvExsg9ANzOzESZJ6J8B7Mlb7smtO55PAD/JW26QtFXSY5I+PIQazcysRJI8I1cDrIsBG0pXAS3A+/NWvzUi9ko6E3hY0jMR8VxBvzagDaCxsZFMJpOk9prV29tb92PsV7KxTls9/H2UWe/YiWQqVWeVPz/+DNeuJKHfA0zOW54E7C1sJOli4IvA+yPiUP/6iNib+3u3pAxwFnBM6EfEemA9QEtLS7S2thY1iFqTyWSo9zH2K9lYVy0evE2VZaatprXz5sq82dKXKvM+x+HPcO1KcnqnA2iWNEXSGGAJcMxdOJLOAr4FXBYRf8xbf6qksbnXE4DzgfwLwGZmVkGDHulHRJ+k5cAWYDSwISJ2SFoDbI2IduCrwDjge5IAno+Iy4DpwLck/ZXsD5hbC+76MTOzCkpyeoeI2AxsLlh3U97ri4/T7z+B2cMp0MzMSsffyDUzSxGHvplZijj0zcxSxKFvZpYiDn0zsxRx6JuZpYhD38wsRRz6ZmYp4tA3M0sRh76ZWYo49M3MUsShb2aWIg59M7MUceibmaWIQ9/MLEUShb6khZI6JXVJWjnA9s9J2inpaUkPSXpb3rZlkn6T+7OslMWbmVlxBg19SaOBdcAlwAxgqaQZBc1+DbRExDuBe4Hbc33/BrgZOAeYC9ws6dTSlW9mZsVIcqQ/F+iKiN0RcRjYBBzzlOqIeCQiDuQWHyP78HSABcADEbE/Il4EHgAWlqZ0MzMrVpLQPwPYk7fck1t3PJ8AfjLEvmZmVkZJnpGrAdbFgA2lq4AW4P3F9JXUBrQBNDY2kslkEpRVu3p7e+t+jP1KNtZpq4e/jzLrHTuRTKXqrPLnx5/h2pUk9HuAyXnLk4C9hY0kXQx8EXh/RBzK69ta0DdT2Dci1gPrAVpaWqK1tbWwSV3JZDLU+xj75Y+1aeWPh7yf7oabS1RR+WSmraa1s0J1Ln2pMu9zHGn9DNeDJKd3OoBmSVMkjQGWAO35DSSdBXwLuCwi/pi3aQswX9KpuQu483PrzMysCgY90o+IPknLyYb1aGBDROyQtAbYGhHtwFeBccD3JAE8HxGXRcR+SV8m+4MDYE1E7C/LSMzMbFBJTu8QEZuBzQXrbsp7ffEb9N0AbBhqgWZmVjqJQt/MRphV44fRt7rXA6y6PA2DmVmKOPTNzFLEoW9mliIOfTOzFHHom5mliEPfzCxFHPpmZini0DczSxGHvplZijj0zcxSxKFvZpYiDn0zsxRx6JuZpYhD38wsRRz6ZmYpkij0JS2U1CmpS9LKAbZfIOkJSX2SLi/YdkTSk7k/7YV9zcyscgZ9iIqk0cA6YB7ZB513SGqPiJ15zZ4HrgZWDLCLVyNiTglqNTOzYUry5Ky5QFdE7AaQtAlYDBwN/Yjozm37axlqNDOzEklyeucMYE/eck9uXVINkrZKekzSh4uqzszMSirJkb4GWBdFvMdbI2KvpDOBhyU9ExHPHfMGUhvQBtDY2Egmkyli97Wnt7e37sfYL3+sd73rN0PeT4bVJaqofHrHTiQzbeTXSQk+e2n9DNeDJKHfA0zOW54E7E36BhGxN/f3bkkZ4CzguYI264H1AC0tLdHa2pp09zUpk8lQ72Psd8xYVy2uai3llpm2mtbOm6tdxuCWDv/B6Kn9DNeBJKd3OoBmSVMkjQGWAInuwpF0qqSxudcTgPPJuxZgZmaVNWjoR0QfsBzYAuwC7omIHZLWSLoMQNLZknqAK4BvSdqR6z4d2CrpKeAR4NaCu37MzKyCkpzeISI2A5sL1t2U97qD7Gmfwn7/CcweZo1mZlYi/kaumVmKOPTNzFLEoW9mliIOfTOzFHHom5mliEPfzCxFHPpmZini0DczSxGHvplZijj0zcxSxKFvZpYiDn0zsxRx6JuZpYhD38wsRRJNrWzWr2nlj4tq//nZfVyd69PdUI6KzKwYiY70JS2U1CmpS9LKAbZfIOkJSX2SLi/YtkzSb3J/lpWqcDMzK96goS9pNLAOuASYASyVNKOg2fPA1cB3Cvr+DXAzcA4wF7hZ0qnDL9vMzIYiyZH+XKArInZHxGFgE3DME64jojsingb+WtB3AfBAROyPiBeBB4CFJajbzMyGIEnonwHsyVvuya1LYjh9zcysxJJcyNUA6yLh/hP1ldQGtAE0NjaSyWQS7r429fb21uwYPz+7r6j2jf/ltT6ZUavLUdKI0Tt2IplpNTDGEnz2avkzXKx6G2uS0O8BJuctTwL2Jtx/D9Ba0DdT2Cgi1gPrAVpaWqK1tbWwSV3JZDLU6hivHsLdO//0TPZj1t1wczlKGjEy01bT2lkDY1z60rB3Ucuf4WLV21iTnN7pAJolTZE0BlgCtCfc/xZgvqRTcxdw5+fWmZlZFQwa+hHRBywnG9a7gHsiYoekNZIuA5B0tqQe4ArgW5J25PruB75M9gdHB7Amt87MzKog0ZezImIzsLlg3U15rzvInroZqO8GYMMwarQRpLvhyqLaZ0atrvvTOma1xNMwmJmliEPfzCxFPPeOWcoUO39Sv+5bF5W4EqsGh75ZyhR7XeY1w7/V06rPp3fMzFLEoW9mliIOfTOzFHHom5mliEPfzCxFHPpmZini0DczSxGHvplZijj0zcxSxKFvZpYiDn0zsxRx6JuZpUii0Je0UFKnpC5JKwfYPlbSd3PbfyWpKbe+SdKrkp7M/bmztOWbmVkxBp1lU9JoYB0wj+yDzjsktUfEzrxmnwBejIi/k7QEuA34aG7bcxExp8R1m5nZECQ50p8LdEXE7og4DGwCFhe0WQxszL2+F/iAJJWuTDMzKwVFxBs3kC4HFkbE3+eWPw6cExHL89psz7XpyS0/B5wDjAN2AM8Cfwa+FBE/H+A92oA2gMbGxvds2rSpBEMbuXp7exk3bly1yxiaPzxZVPPesRMZd2hvmYoZWep+rG957Rf2mv4MF6lWxnrhhRdui4iWwdoleYjKQEfshT8pjtfmD8BbI2KfpPcA/yFpZkT8+ZiGEeuB9QAtLS3R2tqaoKzalclkqNkxrir8Je+NZaatprUzHQ9Gr/uxLn3tISo1/RkuUr2NNcnpnR5gct7yJKDwcOZoG0knAOOB/RFxKCL2AUTENuA5YOpwizYzs6FJcqTfATRLmgL8HlgCFD5vrR1YBvwSuBx4OCJC0ulkw/+IpDOBZmB3yaq3IRvyc1IbSlyImVXUoKEfEX2SlgNbgNHAhojYIWkNsDUi2oF/Bf5NUhewn+wPBoALgDWS+oAjwLURsb8cA7HiDP05qWZWyxI9GD0iNgObC9bdlPf6IHDFAP2+D3x/mDWamVmJ+Bu5ZmYp4tA3M0sRh76ZWYo49M3MUsShb2aWIonu3jEzY9X4115PW13ct7NXvTR4G6sIH+mbmaWIQ9/MLEUc+mZmKeLQNzNLEYe+mVmK+O4dMyu//Dt/iu7rO39KyUf6ZmYp4iP9GjbUOfHB8+KbpZWP9M3MUiRR6EtaKKlTUpeklQNsHyvpu7ntv5LUlLftxtz6TkkLSle6mZkVa9DTO5JGA+uAeWSfhdshqT0iduY1+wTwYkT8naQlwG3ARyXNIPsUrZnAROBBSVMj4kipB5JGfvqVmRUryTn9uUBXROwGkLQJWAzkh/5iYFXu9b3Av0hSbv2miDgE/Db3OMW5ZJ+lazC8uxrMzIqU5PTOGcCevOWe3LoB20REH/AScFrCvmZmViFJjvQ1wLpI2CZJXyS1AW25xV5JnQnqqmUTgD9Vu4jK+G8ea12q4FhXDxQjFVUr/69vS9IoSej3AJPzlicBe4/TpkfSCcB4YH/CvkTEemB9koLrgaStEdFS7ToqwWOtTx5r7UpyeqcDaJY0RdIYshdm2wvatAPLcq8vBx6OiMitX5K7u2cK0Aw8XprSzcysWIMe6UdEn6TlwBZgNLAhInZIWgNsjYh24F+Bf8tdqN1P9gcDuXb3kL3o2wdc5zt3zMyqR9kDcqskSW25U1p1z2OtTx5r7XLom5mliKdhMDNLEYd+lUlaISkkTah2LeUi6auS/q+kpyX9UNIp1a6p1AabqqReSJos6RFJuyTtkPTZatdUbpJGS/q1pPuqXUspOPSrSNJkstNbPF/tWsrsAWBWRLwTeBa4scr1lFTeVCWXADOApbkpSOpRH/D5iJgOnAtcV8dj7fdZYFe1iygVh351/TPwBQb4wlo9iYif5r6pDfAY2e9r1JOjU5VExGGgf6qSuhMRf4iIJ3KvXyYbhnX7LXtJk4BFwP+sdi2l4tCvEkmXAb+PiKeqXUuF/VfgJ9UuosRSOd1Ibjbds4BfVbeSsvofZA/M/lrtQkrFD1EpI0kPAn87wKYvAv8dmF/ZisrnjcYaEf871+aLZE8P/Hsla6uARNON1BNJ44DvA/8YEX+udj3lIOlS4I8RsU1Sa7XrKRWHfhlFxMUDrZc0G5gCPJWdjJRJwBOS5kbE/6tgiSVzvLH2k7QMuBT4QNTffcKJphupF5JOJBv4/x4RP6h2PWV0PnCZpA8CDcCbJf2viLiqynUNi+/THwEkdQMtEVELkzoVTdJC4A7g/RHxQrXrKbXcfFPPAh8Afk926pIrI2JHVQsrg9yU6RuB/RHxj9Wup1JyR/orIuLSatcyXD6nb5XwL8DJwAOSnpR0Z7ULKqXcRer+qUp2AffUY+DnnA98HLgo93/5ZO5I2GqEj/TNzFLER/pmZini0DczSxGHvplZijj0zcxSxKFvZpYiDn0zsxRx6JuZpYhD38wsRf4/TOmaqioyT3kAAAAASUVORK5CYII=\n"
},
"metadata": {
"needs_background": "light"
}
},
{
"output_type": "execute_result",
"execution_count": 65,
"data": {
"text/plain": "0.029529275019493414"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "hist1",
"execution_count": 66,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 66,
"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])"
},
"metadata": {}
}
]
},
{
"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\nbins = sp.linspace(-5,5, 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=\"sample 1\")\nplt.bar(bins[:-1], hist2, label=\"sample 2\")\nplt.legend()\nplt.grid(True)\nplt.show()\n\nD_KL(hist1, hist2) ",
"execution_count": 67,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": "<Figure size 432x288 with 1 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAG1xJREFUeJzt3XuQXGW97vHvk0AyhcHANtRsk4lO0EkMSRR0SFAsHJTEwVDJLotLglihjtsBJeo+ChqOSC5CIWrhOVUnCqm9I2GfDRHxsgeMBBC6vGyBSbjldgaHOMAYT4kJFWwuwYm/88f0xE4zyaye6eme7vV8qlLptdb7rv69k65nVtZa/S5FBGZmlg5jKl2AmZmVj0PfzCxFHPpmZini0DczSxGHvplZijj0zcxSxKFvZpYiDn0zsxRx6JuZpcgxlS6g0KRJk6KxsbHSZYyol19+mTe96U2VLqMsPNba5LGOPlu3bv1zRJw0WLtRF/qNjY1s2bKl0mWMqEwmQ0tLS6XLKAuPtTZ5rKOPpGeTtPPpHTOzFEkU+pJaJXVK6pK04ijtzpcUkprz1l2d69cp6aOlKNrMzIZm0NM7ksYCa4H5QA/QIak9InYWtDse+DzwSN66U4AlwCxgMvCApOkRcbB0QzAzs6SSnNOfC3RFxG4ASRuBxcDOgnZfB74JXJm3bjGwMSIOAL+X1JXb32+HW7iZVae//vWv9PT08Nprr1W6lEQmTpzIrl27Kl3GIXV1dTQ0NHDssccOqX+S0J8CPJ+33APMy28g6TRgakTcI+nKgr4PF/SdUvgGktqANoD6+noymUyi4qtVNput+TH281hr03DGOmHCBOrr65kyZQqSSlvYCDh48CBjx46tdBkARAT79+/nySefJJvNDmkfSUJ/oH+VQ09ekTQG+A5wabF9D62IWAesA2hubo5quFI+HNVyN0ApeKy1aThj3bVrFw0NDVUR+AB/+ctfOP744ytdxiHHH3882WyW5ubmwRsPIEno9wBT85YbgD35NQCzgUzuH/EfgXZJixL0NbMUqpbAH42G+7NLcvdOB9AkaZqkcfRdmG3v3xgR+yNiUkQ0RkQjfadzFkXElly7JZLGS5oGNAGPDqtiMzMbskGP9COiV9JyYDMwFlgfETskrQG2RET7UfrukHQnfRd9e4ErfOeOmeVrXPGzku6v+xsLS7q/4WhpaeHb3/524lMxP/zhD1m1ahW7du3i0UcfHfIpnKNJ9I3ciNgEbCpYd+0R2rYULF8PXD/E+swqZ9XE4trPWA2rFhex//3F7d9q3uzZs/nxj3/MZZddNmLvMeqmYTArtaEeSXbXlbgQq7iXX36ZCy+8kJ6eHg4ePMjXvvY1LrroItasWcPdd9/Nq6++ygc+8AFuueUWJNHS0sKsWbPYtm0bL7zwArfddhs33HAD27Zt46KLLuK6666ju7ub1tZW5s2bx+OPP8706dO57bbbOO644w577/vuu4+VK1dy4MAB3vGOd/D973+fCRMmHNZm5syZI/4z8DQMZpYa9957L5MnT+bJJ59k+/bttLa2ArB8+XI6OjrYvn07r776Kvfcc8+hPuPGjeOXv/wll19+OYsXL2bt2rVs376dW2+9lb179wLQ2dlJW1sbTz31FG9+85v57ne/e9j7/vnPf+a6667jgQce4LHHHqO5uZmbbrqpfAPP49A3s9SYM2cODzzwAF/5ylf41a9+xcSJfafwHnroIebNm8ecOXN48MEH2bFjx6E+H/vYxw71nTVrFm9961sZP348J598Ms8/3/cVpqlTp3LmmWcCcMkll/DrX//6sPd9+OGH2blzJ2eeeSannnoqGzZs4NlnE82PVnI+vWNmqTF9+nS2bt3Kpk2buPrqq1mwYAFf/vKX+exnP8uWLVuYOnUqq1atOuzbwuPGjQNgzJgxjB8//tD6MWPG0NvbC7zxNsrC5Yhg/vz53HHHHSM1tMR8pG9mqbFnzx6OO+44LrnkEq688koee+yxQwE/adIkstksd911V9H7fe655/jtb/tml7njjjv44Ac/eNj2M844g9/85jd0dXUB8Morr/D0008PczRD4yN9q3nddRdXugQ7inLeYrlt2zauuuoqxowZw7HHHsv3vvc9TjjhBD796U8zZ84cGhsbOf3004ve78yZM9mwYQOXXXYZTU1NfOYznzls+0knncStt97K0qVLOXDgAADXXXcd06dPP6zdT37yEz73uc/xwgsvsHDhQk499VQ2b9489AEPQBFvmBWhopqbm8MPUakdo2Ksxd56OUSZGatp6VyZvEMV37I53GkYynGXSqkMNg1Dd3c35513Htu3by9bTQP9DCVtjYhBb+z36R0zsxRx6JuZDUNjY2NZj/KHy6FvZpYiDn0zsxRx6JuZpYhD38wsRXyfvplVVqlvqR1Ft8IWO7XyVVddxd133824ceMOTcp2wgknlLQmH+mbmY0S8+fPZ/v27Tz11FNMnz6dG264oeTv4dA3s9R4+eWXWbhwIe95z3uYPXs2P/jBDwBYs2YNp59+OrNnz6atrY3+L622tLSwYsUKzjrrLGbOnElHRwcf//jHaWpq4pprrgH6vpz1rne9i2XLlvHud7+b888/n1deeeUN733ffffx/ve/n/e+971ccMEFAz7YfMGCBRxzTN8JmDPOOIOenp6S/wwc+maWGtU0tfL69es599xzS/wTSBj6kloldUrqkrRigO2XS9om6QlJv5Z0Sm59o6RXc+ufkHRzqQdgZpZUtUytfP3113PMMcfwiU98oqTjhwQXciWNBdYC84EeoENSe0TszGt2e0TcnGu/CLgJaM1teyYiTi1t2WZmxauGqZU3bNjAPffcwy9+8Ys37KcUkhzpzwW6ImJ3RLwObAQOexBoRLyUt/gmYHTN4mZmxuifWvnee+/lxhtvpL29/Q2PWyyVJLdsTgGez1vuAeYVNpJ0BfBFYBzw4bxN0yQ9DrwEXBMRvxp6uWZWc8p4i+Von1p5+fLlHDhwgPnz5wN9vyxuvrm0Z8UHnVpZ0gXARyPin3PLnwTmRsTnjtD+4lz7ZZLGAxMiYq+k9wE/BWYV/M8ASW1AG0B9ff37Nm7cONxxjWrZbPYND0SuVaNirH98oixvkx0/mQkH9iTv8NbqPes5nH/XiRMn8s53vrPEFY2cgwcPMnbs2CNuf/bZZ7nwwgt55JFHylZTV1cX+/cf/svy7LPPTjS1cpIj/R5gat5yA3C0T/ZG4HsAEXEAOJB7vVXSM8B04LAJ8yNiHbAO+ubTr/j86yNsVMwxXyajYqyrFg/epgSKnk9/6ej5ElGxhjuf/tHmpx9tBptPf8KECYwZM6asY6qrq+O0004bUt8k5/Q7gCZJ0ySNA5YA7fkNJDXlLS4Efpdbf1LuQjCSTgaagN1DqtTMbBSqtqmVBz3Sj4heScuBzcBYYH1E7JC0BtgSEe3AcknnAH8FXgSW5bqfBayR1AscBC6PiH0jMRAzqx4RMSJ3pqTBcJ92mGjunYjYBGwqWHdt3usvHKHfj4AfDadAM6stdXV17N27l7e85S0O/iJFBHv37qWurm7I+/CEa2ZWVg0NDfT09PDCCy9UupREXnvttWGFbKnV1dXR0NAw5P4OfTMrq2OPPZZp06ZVuozEMpnMkC+ajkaee8fMLEUc+mZmKeLQNzNLEYe+mVmKOPTNzFLEoW9mliIOfTOzFHHom5mliEPfzCxFHPpmZini0DczSxGHvplZijj0zcxSxKFvZpYinlrZqkLjip8NuW/36JkK3aziEh3pS2qV1CmpS9KKAbZfLmmbpCck/VrSKXnbrs7165T00VIWb2ZmxRn0SD/3YPO1wHygB+iQ1B4RO/Oa3R4RN+faLwJuAlpz4b8EmAVMBh6QND0iDpZ4HGbVZ9XEYfTdX7o6LFWSHOnPBboiYndEvA5sBBbnN4iIl/IW3wT0P7l3MbAxIg5ExO+Brtz+zMysApKc058CPJ+33APMK2wk6Qrgi8A44MN5fR8u6DtlSJWamdmwJQn9gR5XH29YEbEWWCvpYuAaYFnSvpLagDaA+vp6MplMgrKqVzabrfkx9ivVWL80p3fIfTNjVg/7/ZPIjp9MZkZ53osKf378Ga5eSUK/B5iat9wA7DlK+43A94rpGxHrgHUAzc3N0dLSkqCs6pXJZKj1MfYr1VgvHdbdOyuH/f5JZGaspqWzPO/F0sqe0/dnuHolCf0OoEnSNOAP9F2YvTi/gaSmiPhdbnEh0P+6Hbhd0k30XchtAh4tReGWLt11Fw/eyMwGNWjoR0SvpOXAZmAssD4idkhaA2yJiHZguaRzgL8CL9J3aodcuzuBnUAvcIXv3DEzq5xEX86KiE3ApoJ11+a9/sJR+l4PXD/UAs3MrHQ8DYOZWYo49M3MUsShb2aWIg59M7MUceibmaWIQ9/MLEUc+mZmKeLQNzNLEYe+mVmKOPTNzFLEoW9mliIOfTOzFHHom5mliEPfzCxFHPpmZini0DczSxGHvplZijj0zcxSJFHoS2qV1CmpS9KKAbZ/UdJOSU9J+oWkt+dtOyjpidyf9lIWb2ZmxRn0GbmSxgJrgflAD9AhqT0iduY1exxojohXJH0G+CZwUW7bqxFxaonrNjOzIUhypD8X6IqI3RHxOrARWJzfICIeiohXcosPAw2lLdPMzEohSehPAZ7PW+7JrTuSTwE/z1uuk7RF0sOS/mkINZqZWYkMenoH0ADrYsCG0iVAM/ChvNVvi4g9kk4GHpS0LSKeKejXBrQB1NfXk8lkktRetbLZbM2PsV/Jxjpj9fD3McKy4yeTKVedFf78+DNcvZKEfg8wNW+5AdhT2EjSOcBXgQ9FxIH+9RGxJ/f3bkkZ4DTgsNCPiHXAOoDm5uZoaWkpahDVJpPJUOtj7Feysa5aPHibCsvMWE1L58ryvNnS/eV5nyPwZ7h6JTm90wE0SZomaRywBDjsLhxJpwG3AIsi4k9560+UND73ehJwJpB/AdjMzMpo0CP9iOiVtBzYDIwF1kfEDklrgC0R0Q58C5gA/FASwHMRsQiYCdwi6W/0/YL5RsFdP2ZmVkZJTu8QEZuATQXrrs17fc4R+v0XMGc4BZqZWen4G7lmZini0DczSxGHvplZijj0zcxSxKFvZpYiDn0zsxRx6JuZpYhD38wsRRz6ZmYp4tA3M0sRh76ZWYo49M3MUsShb2aWIg59M7MUceibmaWIQ9/MLEUc+mZmKZIo9CW1SuqU1CVpxQDbvyhpp6SnJP1C0tvzti2T9Lvcn2WlLN7MzIozaOhLGgusBc4FTgGWSjqloNnjQHNEvBu4C/hmru8/ACuBecBcYKWkE0tXvpmZFSPJkf5coCsidkfE68BGYHF+g4h4KCJeyS0+DDTkXn8UuD8i9kXEi8D9QGtpSjczs2IlCf0pwPN5yz25dUfyKeDnQ+xrZmYj6JgEbTTAuhiwoXQJ0Ax8qJi+ktqANoD6+noymUyCsqpXNput+TH2K9lYZ6we/j5GWHb8ZDLlqrPCnx9/hqtXktDvAabmLTcAewobSToH+CrwoYg4kNe3paBvprBvRKwD1gE0NzdHS0tLYZOakslkqPUx9ivZWFctHrxNhWVmrKalc2V53mzp/vK8zxH4M1y9kpze6QCaJE2TNA5YArTnN5B0GnALsCgi/pS3aTOwQNKJuQu4C3LrzMysAgY90o+IXknL6QvrscD6iNghaQ2wJSLagW8BE4AfSgJ4LiIWRcQ+SV+n7xcHwJqI2DciIzEzs0ElOb1DRGwCNhWsuzbv9TlH6bseWD/UAs3MrHT8jVwzsxRJdKRvZqPMqonD6FvZi8BWWT7SNzNLEYe+mVmK+PSOlU3jip8NuW93XQkLMUsxH+mbmaWIQ9/MLEUc+mZmKeJz+lY23XUXV7oEs9Tzkb6ZWYo49M3MUsShb2aWIg59M7MUceibmaWIQ9/MLEUc+mZmKeLQNzNLEYe+mVmKJAp9Sa2SOiV1SVoxwPazJD0mqVfS+QXbDkp6IvenvbCvmZmVz6DTMEgaC6wF5gM9QIek9ojYmdfsOeBS4MoBdvFqRJxaglrNzGyYksy9MxfoiojdAJI2AouBQ6EfEd25bX8bgRrNzKxEkoT+FOD5vOUeYF4R71EnaQvQC3wjIn5a2EBSG9AGUF9fTyaTKWL31Sebzdb8GPsdNtYZqytay0jLjp9MphrGWILPXmo/wzUgSehrgHVRxHu8LSL2SDoZeFDStoh45rCdRawD1gE0NzdHS0tLEbuvPplMhlofY7/DxrpqcUVrGWmZGatp6VxZ6TIGt3T4D0ZP7We4BiS5kNsDTM1bbgD2JH2DiNiT+3s3kAFOK6I+MzMroSSh3wE0SZomaRywBEh0F46kEyWNz72eBJxJ3rUAMzMrr0FDPyJ6geXAZmAXcGdE7JC0RtIiAEmnS+oBLgBukbQj130msEXSk8BD9J3Td+ibmVVIoidnRcQmYFPBumvzXnfQd9qnsN9/AXOGWaOZmZWIv5FrZpYiDn0zsxRx6JuZpYhD38wsRRz6ZmYp4tA3M0sRh76ZWYo49M3MUsShb2aWIg59M7MUceibmaWIQ9/MLEUc+mZmKeLQNzNLEYe+mVmKOPTNzFLEoW9mliKJQl9Sq6ROSV2SVgyw/SxJj0nqlXR+wbZlkn6X+7OsVIWbmVnxBg19SWOBtcC5wCnAUkmnFDR7DrgUuL2g7z8AK4F5wFxgpaQTh1+2mZkNRZIj/blAV0TsjojXgY3A4vwGEdEdEU8Bfyvo+1Hg/ojYFxEvAvcDrSWo28zMhiBJ6E8Bns9b7smtS2I4fc3MrMSOSdBGA6yLhPtP1FdSG9AGUF9fTyaTSbj76pTNZmt+jP0OG+uM1RWtZaRlx08mUw1jLMFnL7Wf4RqQJPR7gKl5yw3AnoT77wFaCvpmChtFxDpgHUBzc3O0tLQUNqkpmUyGWh9jv8PGumrxUdtWu8yM1bR0rqx0GYNbun/Yu0jtZ7gGJDm90wE0SZomaRywBGhPuP/NwAJJJ+Yu4C7IrTMzswoYNPQjohdYTl9Y7wLujIgdktZIWgQg6XRJPcAFwC2SduT67gO+Tt8vjg5gTW6dmZlVQJLTO0TEJmBTwbpr81530HfqZqC+64H1w6jRzMxKJFHom/VrXPGzotp/aU4vl+b6dNeNREVmVgyHvlnarJo4xH7DvwBslefQt6J0111cVPvMmNV011XBHS1mKeEJ18zMUsShb2aWIg59M7MUceibmaWIQ9/MLEUc+mZmKeLQNzNLEYe+mVmKOPTNzFLEoW9mliIOfTOzFHHom5mliEPfzCxFHPpmZini0DczS5FEoS+pVVKnpC5JKwbYPl7SD3LbH5HUmFvfKOlVSU/k/txc2vLNzKwYgz5ERdJYYC0wH+gBOiS1R8TOvGafAl6MiHdKWgLcCFyU2/ZMRJxa4rrNzGwIkhzpzwW6ImJ3RLwObAQWF7RZDGzIvb4L+Igkla5MMzMrBUXE0RtI5wOtEfHPueVPAvMiYnlem+25Nj255WeAecAEYAfwNPAScE1E/GqA92gD2gDq6+vft3HjxhIMbfTKZrNMmDCh0mUMzR+fKKp5dvxkJhzYM0LFjC41P9a3/v0/7FX9GS5StYz17LPP3hoRzYO1S/KM3IGO2At/UxypzR+Bt0XEXknvA34qaVZEvHRYw4h1wDqA5ubmaGlpSVBW9cpkMlTtGFcV/ifv6DIzVtPSmY5n5Nb8WJf+/cHoVf0ZLlKtjTXJ6Z0eYGrecgNQeDhzqI2kY4CJwL6IOBARewEiYivwDDB9uEWbmdnQJAn9DqBJ0jRJ44AlQHtBm3ZgWe71+cCDERGSTspdCEbSyUATsLs0pZuZWbEGPb0TEb2SlgObgbHA+ojYIWkNsCUi2oF/A/5dUhewj75fDABnAWsk9QIHgcsjYt9IDMTMzAaX5Jw+EbEJ2FSw7tq8168BFwzQ70fAj4ZZo5mZlYi/kWtmliIOfTOzFEl0esfMjFUT//56xuribt9dtX/wNlYWDv2UalzxsyH1664rcSFmVlY+vWNmliI+0k+p7rqLK12CmVWAj/TNzFLEoW9mliIOfTOzFHHom5mliEPfzCxFHPpmZini0DczSxGHvplZijj0zcxSxN/IrWJDnT8HPIeOldewPqvfWFjCSsyhX8U8lYKZFStR6EtqBf4XfY9L/NeI+EbB9vHAbcD7gL3ARRHRndt2NfAp+h6X+PmI2Fyy6s2sKgzvAMXTMpfSoOf0cw82XwucC5wCLJV0SkGzTwEvRsQ7ge8AN+b6nkLf83JnAa3Ad/sflG5mZuWX5ELuXKArInZHxOvARqDw6QmLgQ2513cBH5Gk3PqNEXEgIn4PdOX2Z2ZmFZDk9M4U4Pm85R5g3pHaRESvpP3AW3LrHy7oO2XI1Y5SxV6k+tKcXi7N9fF5ebNB5D+xq0iNr90+pH61fPE4SehrgHWRsE2SvkhqA9pyi1lJnQnqqlqfh0nAn2HgH1Bt+e+Hxlr7PNbR57wh9dKNhy1WyVh5e5JGSUK/B5iat9wA7DlCmx5JxwATgX0J+xIR64B1SQquBZK2RERzpesoB4+1Nnms1SvJOf0OoEnSNEnj6Lsw217Qph1Ylnt9PvBgRERu/RJJ4yVNA5qAR0tTupmZFWvQI/3cOfrlwGb6btlcHxE7JK0BtkREO/BvwL9L6qLvCH9Jru8OSXcCO4Fe4IqIODhCYzEzs0Go74DcyklSW+6UVs3zWGuTx1q9HPpmZiniCdfMzFLEoV9hkq6UFJImVbqWkSLpW5L+r6SnJP1E0gmVrqnUJLVK6pTUJWlFpesZKZKmSnpI0i5JOyR9odI1jTRJYyU9LumeStdSCg79CpI0FZgPPFfpWkbY/cDsiHg38DRwdYXrKamEU5XUil7gSxExEzgDuKKGx9rvC8CuShdRKg79yvoO8GUG+MJaLYmI+yKiN7f4MH3f16glSaYqqQkR8ceIeCz3+i/0hWHNfcu+n6QGYCHwr5WupVQc+hUiaRHwh4h4stK1lNl/A35e6SJKbKCpSmo2CPtJagROAx6pbCUj6n/Sd2D2t0oXUiqeT38ESXoA+McBNn0V+B/AgvJWNHKONtaI+M9cm6/Sd3rgP8pZWxkkmm6klkiaAPwI+JeIeKnS9YwESecBf4qIrZJaKl1PqTj0R1BEnDPQeklzgGnAk32TkdIAPCZpbkT8vzKWWDJHGms/ScvomwjlI1F79wknmm6kVkg6lr7A/4+I+HGl6xlBZwKLJH0MqAPeLOn/RMQlFa5rWHyf/iggqRtojohqmNSpaLmH8NwEfCgiXqh0PaWWm2/qaeAjwB/om7rk4ojYUdHCRkBuyvQNwL6I+JdK11MuuSP9KyNiaDO4jSI+p2/l8L+B44H7JT0h6eZKF1RKuYvU/VOV7ALurMXAzzkT+CTw4dy/5RO5I2GrEj7SNzNLER/pm5mliEPfzCxFHPpmZini0DczSxGHvplZijj0zcxSxKFvZpYiDn0zsxT5/1w4gVjkFiFeAAAAAElFTkSuQmCC\n"
},
"metadata": {
"needs_background": "light"
}
},
{
"output_type": "execute_result",
"execution_count": 67,
"data": {
"text/plain": "0.003290250486473783"
},
"metadata": {}
}
]
},
{
"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\n\nbins = sp.linspace(5-8,5+8, 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(10, 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": 94,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": "<Figure size 432x288 with 1 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHtJJREFUeJzt3X10VPW97/H31/AQangS21QJC7DNoYkJCuQgt1oNBRTBxYNiRS3CEot2FfRWe++lB6+A57RSe/VaV/EUVG5tDza19kGWohwr5F711AqiPIWLBMzCCKKgFxMkQOB7/5jJMEye9iSTzOD+vNbKYn57/357fydMPtn5zd57zN0REZFwOCvdBYiISOdR6IuIhIhCX0QkRBT6IiIhotAXEQkRhb6ISIgo9EVEQkShLyISIgp9EZEQ6ZLuAhKde+65PmjQoFj78OHDnH322ekrqBmqKzmqKzmZWhdkbm1hr+utt9464O5fbrWju2fU14gRIzzeunXrPBOpruSoruRkal3umVtb2OsCNniAjNX0johIiCj0RURCRKEvIhIiGfdGroh0jOPHj1NdXU1dXV27ttO7d2+2b9+eoqpSJyx1ZWdnk5eXR9euXds0XqEvEhLV1dX07NmTQYMGYWZt3k5NTQ09e/ZMYWWpEYa63J2DBw9SXV3N4MGD27QNTe+IhERdXR39+vVrV+BLepkZ/fr1a9dfawp9kRBR4J/52vt/qNAXEQmRQHP6ZjYe+AWQBTzh7ksS1t8B/AA4AdQCc9y9Irrux8Ds6Lo73X1N6soXkbYaNP+FlG6vasnEltdXVXHNNdewdevWRutuu+027r77bgoLC1vcxiOPPMI555zDLbfcwh/+8AcWLVrE9u3befPNNxkyZEis3wMPPMCTTz5JVlYWjz76KFdddVWL2/3lL3/JI488wq5du/j4448599xzY+uOHz/OqFGj+Mtf/sItt9zChx9+yFlnncWcOXO46667APjRj37EhAkT+Pa3v93ifjJBq6FvZlnAUmAcUA2sN7NVDaEe9bS7/yrafxLwMDDezAqB6cCFwPnAX83sH9z9RIqfh0inaE9QthaKYfbEE0+02qe+vp4VK1awceNGAIqKivjTn/7E7bffflq/iooKysrK2LZtG3v37mXs2LG8++67ZGVlNbvtSy+9lGuuuYbS0tJG61577TW++c1v0qVLFx566CGGDx9OTU0NI0aMYNy4cRQWFjJv3jy+973vnRGhH2R6ZyRQ6e673f0YUAZMju/g7p/FNc8GPPp4MlDm7kfd/T2gMro9EQmh+vp6Zs6cydChQ5k2bRqff/45AKWlpWzYsAGAnJwcFixYwEUXXcSoUaPYv38/AGvXrmX48OF06RI5Vi0oKDjt6L7Bc889x/Tp0+nevTuDBw/m61//Om+++WaLdQ0bNoz4e37Fe+mll7j66qs577zzGD58OAA9e/akoKCADz74AICBAwdy8OBBPvzww+S/KZ0sSOj3B96Pa1dHl53GzH5gZruAB4E7kxkrIuGwY8cO5syZw+bNm+nVqxePPfZYoz6HDx9m1KhRbNq0icsvv5zHH38cgNdff50RI0a0uo8PPviAAQMGxNp5eXmxcG6LdevWNfoLoKqqirfffptLLrkktmz48OG8/vrrbd5PZwkyp9/UW8XeaIH7UmCpmd0E3AvMDDrWzOYAcwByc3MpLy+PrautrT2tnSlUV3K+KHXdU1zf5n0ls5+O+H717t2bmpqalG4zXmvbrq2tJS8vj6FDh1JTU8O1117Lr371K26//XZOnDjB4cOHqampoVu3blxxxRXU1NRQWFjIunXrqKmpYc+ePQwaNKjRfhrGnjhxgpqaGo4ePcqRI0di/Y4fP05dXV2g5+7u1NbW0r17dwD27dtHr169YttueB5Tp07lgQcewMxiy/v06cPu3bubrC/V3/e6uro2vz6ChH41MCCunQfsbaF/GfCvyYx19+XAcoCSkhKP/61aXl7e5Dxbuqmu5HxR6prVnjn9m4PvpyO+X9u3b+/Qi5da23ZOTg5nnXVWrN+XvvQlunbtSs+ePcnKyuLss8+mZ8+edO3alV69esXGmBk9e/akV69escfxGsZmZWXRs2dPLrjgAg4cOBDrt3//fr72ta8Feu5mRk5OTqzvM888w8SJE2Pt48ePM23aNGbMmMHNN9982tiTJ0/St2/fRvvpiIvGsrOzGTZsWJvGBpneWQ/km9lgM+tG5I3ZVfEdzCw/rjkR2Bl9vAqYbmbdzWwwkA+0PLkmIl9Ye/bs4W9/+xsAv/vd77jssssCjy0oKKCysrLVfpMmTaKsrIyjR4/y3nvvsXPnTkaOjLyVOGbMmKSmehrm8yHyV8Ds2bMpKCjg7rvvbtT33XffpaioKPC206XVI313rzezucAaIqdsrnD3bWZ2P5H7N68C5prZWOA48CmRqR2i/Z4BKoB64Ac6c0ckM7T1bKL2HLkWFBTw1FNPcfvtt5Ofn8/3v//9wGOvvvpqZsyYEWv/+c9/Zt68eXz88cdMnDiRoqIiXnnlFS688EK+853vUFhYSJcuXVi6dClZWVmcPHmSyspKzjnnnEbbfvTRR3nwwQf58MMPGTp0KBMmTGDZsmXs3LmTb3zjG0DkPYXf/va3FBcXc/HFFwPw05/+lAkTJnD8+HEqKyspKSlp0/elMwU6T9/dVwOrE5bdF/f4rhbG/gT4SVsLFJEvhkGDBlFRUdHkusT38RpMmzaNadOmAZEzZPr168fOnTvJz89n6tSpTJ06NdY3ft58wYIFLFiw4LR9VFRUcN1119GjR49G+7/zzju58847T1v22muvMWrUqFj7sssuI/JZJY09//zzTJs2LXZmUSbL/ApFRKKWLFnCvn37yM/Pb71zgqKiIh5++OHA/S+77LLA00/19fXcc889SdeUDgp9ETljDBkypMlz89Pt+uuvT3cJgeneOyIiIaLQFxEJEYW+iEiIKPRFREJEoS8iEiIKfRGRENEpmyJhtah3m4Y1ey3uokOBxi9btow77riDiooKCgoKgMiVui+++CK5ubmMHz+etWvXtnj/+4525MiR0+q49dZbef755/nKV77S5IfAtNexY8cYO3Ysa9eupUuXLo3aqaQjfRHpVJs3b+biiy/mhRciN687evQo+/fvZ+DAgaxYsYJrr702rYEPNKpj1qxZvPTSSx22v27dujFmzBh+//vfN9lOJYW+iHSqLVu2MH/+/Fjob9u2jYKCAsyMlStXMnnyqc9omjp1Kvfeey/f+ta3+OpXv8pf//pXAKZPn84NN9zAJZdcwsCBA2Pbak6y20ms4/LLL2/ynj1t0dw+p0yZwsqVK2P9EtupotAXkU5VUVHBpEmT+Oijjzh06BBbtmyhuLiYY8eOsXv37tM+wWrr1q306dOHV199lcceeywWgps2beKCCy7g73//OytXrmTx4sUt7jOZ7TRVRyo1V3tRURHr16+P9Utsp4rm9EWk07z//vv069ePHj16MG7cONasWcPmzZsZOnQoBw4coE+fPrG+n3/+OYcOHeKHP/whELm/TZ8+fThy5AgHDhxg4cKFABQWFvLpp582u89kt5NYR1Bjx45t8uMS7733XqZPnw7QYu1ZWVl069YtdhfTxHaqKPRFpNNs3ryZ4uJiACZMmMDKlSvZt28fU6ZMoUePHtTV1cX6btu2jREjRsTm1Tdv3kxRURFbt24lPz+f7OxsADZu3MhFF13U7D6T3U5iHUE1TBklir/7Z2u1Hz16NLauqXYqaHpHRDpNw1QOwBVXXMGrr74a+0XQt29fTpw4EQvcrVu3xu5bD8T+Iti0aRN79uyhrq6Ow4cPs3DhwthRPDT+oJRkt5NYRyq1VPvBgwf58pe/TNeuXZtsp4qO9EXCKuAplonaM92wZcsWrrvuOgC6d+9OcXExb7/9dmw65corr+S1115j7NixbNmy5bQPHt+6dStFRUX85je/4eabb6a0tJTPPvuMf/qnf+LSSy+lpqamyQ9KSXY7iXUA3HjjjZSXl3PgwAHy8vJYvHgxs2fPTvr5b9q0qdl9rlu3jgkTJsT6JrZTRaEvIp0m8WyU55577rT23Llzefjhhxk7dmyje9/v3r0biATn448/zs9+9rNG22/qg1Lasp34OiDy0Y6p0NI+n376aR544IFm26mi6R0RyRjDhg1j9OjRnDjR/Keq7tq1q9kPUUnmg1Ja2k6QOtqiuX0eO3aMKVOmxD4rILGdSjrSF5GMcuutt7a4PpkPNm/PdlqrI5X77NatG7fcckuz7VTSkb6ISIgo9EVEQkTTO3LGGTS/5UvuW1O1ZGKKKjnzuDtmlu4ypB3cvV3jdaQvEhLZ2dkcPHiw3aEh6ePuHDx4sF0XbOlIXyQk8vLyqK6u5uOPP27Xdurq6lJ+lWgqhKWu7Oxs8vLy2jxeoS8SEl27dmXw4MHt3k55eTnDhg1LQUWppbqCCTS9Y2bjzWyHmVWa2fwm1t9tZhVmttnMXjGzgXHrTpjZO9GvVaksXkREktPqkb6ZZQFLgXFANbDezFa5e0Vct7eBEnf/3My+DzwI3BBdd8TdL0ZERNIuyJH+SKDS3Xe7+zGgDJgc38Hd17n759HmG0DbJ5xERKTDBAn9/sD7ce3q6LLmzAZejGtnm9kGM3vDzKa0oUYREUkRa+30LTO7HrjK3W+LtmcAI919XhN9vwvMBa5w96PRZee7+14zuwBYC4xx910J4+YAcwByc3NHlJWVxdbV1taSk5PTjqfYMVRXclJZ15YP2nZ3yAbF/U99IHiydbV330Hl9oD9R06142tOtzC8xlKps+oaPXr0W+5e0lq/IGfvVAMD4tp5wN7ETmY2FlhAXOADuPve6L+7zawcGAacFvruvhxYDlBSUuKlpaWxdeXl5cS3M4XqSk4q65rV3ouzbj5VR7J1tXffQd1TXM9DW079eMbXnG5heI2lUqbVFWR6Zz2Qb2aDzawbMB047SwcMxsGLAMmuftHccv7mln36ONzgUuB+DeARUSkE7V6pO/u9WY2F1gDZAEr3H2bmd0PbHD3VcDPgRzgD9FLvPe4+ySgAFhmZieJ/IJZknDWj4iIdKJAF2e5+2pgdcKy++Iej21m3H8Axe0pUCTlFsXNjw9ZDIsmN9+3kadTXo5IZ9K9d0REQkShLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaL76YskoSr7pjaPHVSn0z0l/XSkLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaLQFxEJEYW+iEiIKPRFREJEoS8iEiIKfRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCRGFvohIiCj0RURCRKEvIhIiCn0RkRAJFPpmNt7MdphZpZnNb2L93WZWYWabzewVMxsYt26mme2Mfs1MZfEiIpKcVkPfzLKApcDVQCFwo5kVJnR7Gyhx96HAs8CD0bHnAAuBS4CRwEIz65u68kVEJBlBjvRHApXuvtvdjwFlwOT4Du6+zt0/jzbfAPKij68CXnb3T9z9U+BlYHxqShcRkWQFCf3+wPtx7erosubMBl5s41gREelA5u4tdzC7HrjK3W+LtmcAI919XhN9vwvMBa5w96Nm9l+A7u7+L9H1/x343N0fShg3B5gDkJubO6KsrCy2rra2lpycnHY8xY6hupKTyrq2fHCoXeOLz3ov9ri2+/nkHN3b3pIC2XJycOC+uT1g/5FT7eL+vTugorYJw2sslTqrrtGjR7/l7iWt9esSYFvVwIC4dh7Q6KfEzMYCC4gGftzY0oSx5Ylj3X05sBygpKTES0tPDSkvLye+nSlUV3JSWdes+S+0a3xV9sLY4/IhiyndsbCF3qkzq+7pwH3vKa7noS2nfjyrbi7tgIraJgyvsVTKtLqCTO+sB/LNbLCZdQOmA6viO5jZMGAZMMndP4pbtQa40sz6Rt/AvTK6TERE0qDVI313rzezuUTCOgtY4e7bzOx+YIO7rwJ+DuQAfzAzgD3uPsndPzGzfybyiwPgfnf/pEOeiYiItCrI9A7uvhpYnbDsvrjHY1sYuwJY0dYCRUQkdQKFvkgmqcq+Kd0liJyxdBsGEZEQUeiLiISIQl9EJEQU+iIiIaI3ckU6STJvQJeftfi0i8igfVchizTQkb6ISIgo9EVEQkShLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaLQFxEJEYW+iEiIKPRFREJEoS8iEiIKfRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCRGFvohIiCj0RURCRKEvIhIigULfzMab2Q4zqzSz+U2sv9zMNppZvZlNS1h3wszeiX6tSlXhIiKSvFY/GN3MsoClwDigGlhvZqvcvSKu2x5gFvCjJjZxxN0vTkGtIiLSTq2GPjASqHT33QBmVgZMBmKh7+5V0XUnO6BGERFJkSDTO/2B9+Pa1dFlQWWb2QYze8PMpiRVnYiIpJS5e8sdzK4HrnL326LtGcBId5/XRN9fA8+7+7Nxy853971mdgGwFhjj7rsSxs0B5gDk5uaOKCsri62rra0lJyenjU+v46iu5KS0rn3vpGY7QG3388k5ujdl20uVxLq2nBzcru0V9+/d3pJiQvEaS6HOqmv06NFvuXtJa/2CTO9UAwPi2nlA4J8Sd98b/Xe3mZUDw4BdCX2WA8sBSkpKvLS0NLauvLyc+HamUF3JSaxr0PwX2rytquyFKagoonzIYkp3pG57qZJY16y6p9u1vaqbS9tZ0SlnymssU2RaXUGmd9YD+WY22My6AdOBQGfhmFlfM+sefXwucClx7wWIiEjnajX03b0emAusAbYDz7j7NjO738wmAZjZP5pZNXA9sMzMtkWHFwAbzGwTsA5YknDWj4iIdKIg0zu4+2pgdcKy++Ierycy7ZM47j+A4nbWKCIiKaIrckVEQkShLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaLQFxEJEYW+iEiIKPRFREJEoS8iEiIKfRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCRGFvohIiCj0RURCRKEvIhIigT4uUSTVqrJvSncJIqGk0Bc5A7T/l+ShlNQhZz5N74iIhIhCX0QkRBT6IiIhotAXEQkRhb6ISIgECn0zG29mO8ys0szmN7H+cjPbaGb1ZjYtYd1MM9sZ/ZqZqsJFRCR5rYa+mWUBS4GrgULgRjMrTOi2B5gFPJ0w9hxgIXAJMBJYaGZ921+2iIi0RZAj/ZFApbvvdvdjQBkwOb6Du1e5+2bgZMLYq4CX3f0Td/8UeBkYn4K6RUSkDYKEfn/g/bh2dXRZEO0ZKyIiKRbkilxrYpkH3H6gsWY2B5gDkJubS3l5eWxdbW3tae1MobqS06iuIYvTVku82u7nU54htcRLeV0pfE2cMa+xDJFpdQUJ/WpgQFw7D9gbcPvVQGnC2PLETu6+HFgOUFJS4qWlp4aUl5cT384Uqis5jepaNLnZvp2pfMhiSncsTHcZjaS8rhtTdxuGM+Y1liEyra4g0zvrgXwzG2xm3YDpwKqA218DXGlmfaNv4F4ZXSYiImnQaui7ez0wl0hYbweecfdtZna/mU0CMLN/NLNq4HpgmZlti479BPhnIr841gP3R5eJiEgaBLrLpruvBlYnLLsv7vF6IlM3TY1dAaxoR40iIpIiuiJXRCREFPoiIiGi0BcRCRGFvohIiCj0RURCRKEvIhIiCn0RkRBR6IuIhIhCX0QkRBT6IiIhotAXEQkRhb6ISIgo9EVEQkShLyISIgp9EZEQUeiLiISIQl9EJEQU+iIiIaLQFxEJEYW+iEiIKPRFREJEoS8iEiIKfRGREFHoi4iESJd0FyBnrkHzXwjc957iembF9a/K7oiKRKQ1gULfzMYDvwCygCfcfUnC+u7Ab4ARwEHgBnevMrNBwHZgR7TrG+5+R2pKl3Sryr4pcN/ysxZTlb2wA6uRliTzCzpR1ZKJKaxE0q3V0DezLGApMA6oBtab2Sp3r4jrNhv41N2/bmbTgZ8BN0TX7XL3i1Nct4gkIZlf0I0dSlkdkn5B5vRHApXuvtvdjwFlwOSEPpOBp6KPnwXGmJmlrkwREUmFIKHfH3g/rl0dXdZkH3evJ3Jo0C+6brCZvW1m/9vMvtXOekVEpB3M3VvuYHY9cJW73xZtzwBGuvu8uD7bon2qo+1dRP5CqAVy3P2gmY0A/gJc6O6fJexjDjAHIDc3d0RZWVlsXW1tLTk5Oe1+oqmmuoB97wTuWtv9fHKO7u3AYtpGdQVw3umzs3rtJ6ez6ho9evRb7l7SWr8gb+RWAwPi2nlA4quxoU+1mXUBegOfeOQ3ylEAd38r+svgH4AN8YPdfTmwHKCkpMRLS0tj68rLy4lvZwrVBSxKnOVrXvmQxZTuyLw3clVXADeePqev135yMq2uINM764F8MxtsZt2A6cCqhD6rgJnRx9OAte7uZvbl6BvBmNkFQD6wOzWli4hIslo90nf3ejObC6whcsrmCnffZmb3AxvcfRXwJPBbM6sEPiHyiwHgcuB+M6sHTgB3uPsnHfFERESkdYHO03f31cDqhGX3xT2uA65vYtwfgT+2s0YREUkR3YZBRCREFPoiIiGi0BcRCRGFvohIiCj0RURCRKEvIhIiCn0RkRBR6IuIhIhCX0QkRBT6IiIhotAXEQkRhb6ISIgo9EVEQkShLyISIoFurSxfYIt6p7sCyXSJr5Ehi4N/atqiQ633kU6lI30RkRBR6IuIhIhCX0QkRBT6IiIhotAXEQkRhb6ISIjolE0R6TCD5r/Q5rFVSyamsBJpoCN9EZEQUeiLiISIQl9EJEQCzemb2XjgF0AW8IS7L0lY3x34DTACOAjc4O5V0XU/BmYDJ4A73X1NyqqXiPjL5JO5RF6kg1Vl39SO0bqFQ0do9UjfzLKApcDVQCFwo5kVJnSbDXzq7l8H/ifws+jYQmA6cCEwHngsuj0REUmDINM7I4FKd9/t7seAMiDxUHIy8FT08bPAGDOz6PIydz/q7u8BldHtiYhIGgSZ3ukPvB/XrgYuaa6Pu9eb2SGgX3T5Gwlj+7e52i8q3elSpLFkfy7ipzZ1d89mBQl9a2KZB+wTZCxmNgeYE23WmtmOuNXnAgcC1NnZMrSuH6qupKiu5GVqbXF1LW4qetKms75fA4N0ChL61cCAuHYesLeZPtVm1gXoDXwScCzuvhxY3tTOzWyDu5cEqLNTqa7kqK7kZGpdkLm1qa5ggszprwfyzWywmXUj8sbsqoQ+q4CZ0cfTgLXu7tHl082su5kNBvKBN1NTuoiIJKvVI/3oHP1cYA2RUzZXuPs2M7sf2ODuq4Angd+aWSWRI/zp0bHbzOwZoAKoB37g7ic66LmIiEgrAp2n7+6rgdUJy+6Le1wHXN/M2J8AP2lHjU1O+2QA1ZUc1ZWcTK0LMrc21RWARWZhREQkDHQbBhGREDmjQt/MfmRmbmbnprsWADP7uZn9XzPbbGZ/NrM+aa5nvJntMLNKM5ufzloamNkAM1tnZtvNbJuZ3ZXumuKZWZaZvW1mz6e7lgZm1sfMno2+trab2X9Kd00AZvbD6P/hVjP7nZllp6mOFWb2kZltjVt2jpm9bGY7o//2zZC6Mioj4AwKfTMbAIwD9qS7ljgvA0XuPhR4F/hxugoJeLuMdKgH7nH3AmAU8IMMqavBXcD2dBeR4BfAS+7+DeAiMqA+M+sP3AmUuHsRkZM6pqepnF8Tua1LvPnAK+6eD7wSbXe2X9O4rozJiAZnTOgTuafPf6WJi7vSxd3/3d3ro803iFyHkC5BbpfR6dx9n7tvjD6uIRJgGXFVtpnlAROBJ9JdSwMz6wVcTuSMONz9mLv/v/RWFdMF6BG9FudLNHHNTWdw9/9D5CzBePG3gnkKmNKpRdF0XRmWEcAZEvpmNgn4wN03pbuWFtwKvJjG/Td1u4yMCNcGZjYIGAb8Pb2VxDxC5EDiZLoLiXMB8DHwv6LTTk+Y2dnpLsrdPwD+B5G/tPcBh9z939Nb1Wly3X0fRA40gK+kuZ6mpDsjgAwKfTP7a3SuMPFrMrAAuK+1baShroY+C4hMY6xMR40NZTSxLGP+KjKzHOCPwH92988yoJ5rgI/c/a1015KgCzAc+Fd3HwYcJj1TFaeJzpFPBgYD5wNnm9l301vVmSNDMgLIoM/IdfexTS03s2IiL7RNkRt3kgdsNLOR7v5huuqKq28mcA0wxtN7/mugW16kg5l1JRL4K939T+muJ+pSYJKZTQCygV5m9m/unu4gqwaq3b3hr6FnyYDQB8YC77n7xwBm9ifgm8C/pbWqU/ab2Xnuvs/MzgM+SndBDTIoI4AMOtJvjrtvcfevuPsgdx9E5IdieGcEfmuiHy7z34BJ7v55mssJcruMThe9xfaTwHZ3fzjd9TRw9x+7e170NTWdyK1D0h34RF/X75vZkOiiMUSuaE+3PcAoM/tS9P90DBnwBnOc+FvBzASeS2MtMRmWEcAZEPoZ7pdAT+BlM3vHzH6VrkKibxY13C5jO/CMu29LVz1xLgVmAN+Ofo/eiR5dS/PmASvNbDNwMfDTNNdD9C+PZ4GNwBYi2ZGWK03N7HfA34AhZlZtZrOBJcA4M9tJ5Cy/JS1toxPrypiMaKArckVEQkRH+iIiIaLQFxEJEYW+iEiIKPRFREJEoS8iEiIKfRGREFHoi4iEiEJfRCRE/j+ps1HVQmklYwAAAABJRU5ErkJggg==\n"
},
"metadata": {
"needs_background": "light"
}
},
{
"output_type": "execute_result",
"execution_count": 94,
"data": {
"text/plain": "0.3162509415508368"
},
"metadata": {}
}
]
},
{
"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": [
{
"output_type": "display_data",
"data": {
"text/plain": "<Figure size 432x288 with 1 Axes>",
"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"
},
"metadata": {
"needs_background": "light"
}
},
{
"output_type": "execute_result",
"execution_count": 96,
"data": {
"text/plain": "0.005019163960424132"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "n*p, n*p*(1-p)",
"execution_count": 81,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 81,
"data": {
"text/plain": "(50.0, 25.0)"
},
"metadata": {}
}
]
},
{
"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, 30+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(10000, 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": 110,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": "<Figure size 432x288 with 1 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD8CAYAAABthzNFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X2UVfV97/H31yGAdXgQohPiGMEFWgZUyFBQc7VOMIKkEQ1wxSiaiot0FXKNuVktmEYMa5HG9jakjbqaXKExiBkMbZq5CtIkDGltFBDkaQZphofG8SGEh6CjPAj53j/2j3Fz2Odxzpk5o5/XWmfN3r/93b/z3fvMnO/svc8+P3N3REREUp3V1QmIiEh5UoEQEZFEKhAiIpJIBUJERBKpQIiISCIVCBERSaQCISIiiVQgREQkkQqEiIgk6tHVCeTjwx/+sA8ePLigdd9++23OOeec4iZUJOWcG5R3fsqtMOWcG5R3ft0tt40bN+539/MK6tDdu82jtrbWC9XY2FjwuqVWzrm5l3d+yq0w5Zybe3nn191yA170At9zdYpJREQSqUCIiEgiFQgREUnUrS5Si0jxvfvuu7S2tnL06NGuTqVdv3792LFjR1enkahcc+vduzdmVtQ+VSBEPuBaW1vp06cPgwcPLvobTKHeeust+vTp09VpJCrH3NydAwcOFP3TVTrFJPIBd/ToUQYOHFg2xUHyZ2YMHDiQioqKovarAiEiKg7vA6V4DVUgREQkkQqEiIgkUoEQKZFtrx7u6hREOkSfYhKR0wye+0xR+9v7zU/nFLd27VoWL17M0qVLT2v/5S9/yerVq/n617+edt0jR44wceJE1qxZQ0VFBXfffTdPP/00559/Ptu3b2+Pe/bZZ7n33ns5efIk99xzD3Pnzs26LLV99uzZWftKki4ngC984QvceeedfOxjH+POO+/kjTfe4KyzzmLWrFnce++9HD9+nOuvv541a9bQo0fnvW3rCEJEysLmzZsZPXr0Ge1XX311xuIAsGTJEj772c+2f4rn85//PM8+++xpMSdPnmT27NmsWrWK5uZmfvjDH9Lc3JxxWVL7yy+/nLGvdJJyOmXdunVceeWV9OjRg7/7u79jx44dvPDCCzzyyCM0NzfTs2dPxo8fz/LlyzM+R7GpQIhIWdiyZQuvvvoq48aN4/LLL2ft2rUATJs2jeeeew6AW265hb/6q7/immuu4SMf+Qg/+9nPAFi2bBmTJ09u7+vaa69lwIABp/W/fv16hg4dysUXX0zPnj2ZPn06P/nJTzIuS2p/5plnMvaVTlJOADt27OCSSy6hoqKCQYMG8fGPfxyAPn36MHz4cF599VUAbr75ZpYtW1bAni2cCoSIlIXNmzfTp08f1q1bx6JFi/ja174GwPbt27nsssvap/v3789//Md/8Oijj7Js2TKOHz/O7t27yTYUwKuvvsqFF17YPl9dXd3+5ptuWVL7a6+9lrGvfK1atYqJEyee0b53715eeuklxo0bB8DIkSPZsGFDQc9RKBUIEelyJ06c4MCBA9x///0AXH755ezfv5+jR4/y7rvv0q9fP9555x0OHz7Mfffd175O//792b9/P/3798/6HNE3X5/u1L0D6Zbl216I1atXn1Eg2tramDJlCt/+9rfp27cvABUVFfTs2ZO33nqroOcphAqEiHS55uZmhg4dSs+ePYHoaOKKK66gqamJmpoaAJqamqitrW2/zrB161ZGjhzJ2WefndP3SFVXV/PKK6+0z7e2tvLRj34047Kk9kGDBmXsKx/vvPMOv/vd705b991332XKlCncfvvtfPaznz0t/tixY/Tu3Tvv5ylUTgXCzCaa2U4zazGzMy7Vm1kvM1selq8zs8GxZfNC+04zmxBrv8/Mmsxsu5n90Mw6b6tFpKxs2bKFPXv2cOzYMdra2vjmN7/Jl770JbZt28bll18ORKeXRo0a1b7O1q1bufzyyzn33HM5efJk1iLxR3/0R/zqV79iz549HD9+nPr6em666aaMy5LaJ02alLEvgPHjx+d0yqmxsZG6urr2eXdn5syZDB8+nC9/+cunxR44cIDzzjuPD33oQ9l3aJFk/byUmVUAjwCfAlqBDWbW4O7xS/YzgUPuPtTMpgMPAbeaWQ0wHRgBfBT4mZldAnwE+F9AjbsfMbOnQtz3i7dpIlKIXD+WWkxbtmzh9ttv5+qrr+bIkSN85Stf4corr+RHP/pR+zn4bdu2tU9DVDBGjhwJwA033MBzzz3H9ddfD8Btt93G2rVr2b9/P9XV1Xz9619n5syZPPzww0yYMIGTJ09y9913M2LECAB69OiRdllq+/DhwzPG//73v6elpeWMC9JJOb300ktMnTq1PeY///M/Wbp0KZdddll7MfzGN77BpEmTaGxsZNKkSaXY/ellG3IOuApYHZufB8xLiVkNXBWmewD7AUuNPRUHXAC8AgwI8U8DN2TLRUOOdo1yzq+cc/uHJ/61q1NIK77fmpubuy6RNN5888284jdt2uR33HFHibI5Xbbctm3b5vfdd19OfY0ePdqPHz+eU+wtt9ziL7/8csaYTZs2ndFGiYccPfVmfkpraEuMcfcTwGFgYLp13f1V4P8AvwZeBw67+7/lkIuIyBlGjx5NXV0dJ0+e7OpUGDlyJN/61rdyit20aVNOp4yOHz/OzTffzKWXXtrR9PKSyy15SZfmUy/hp4tJbDezc4HJwBDgd8CPzOwOd3/ijCc3mwXMAqiqqmr/bHS+2traCl631Mo5Nyjv/Mo5t6qzKdvc4vutX79+nfrJmFycPHky75ymTZvGO++8U6KM3lNIbsVwyy23ZH1edy/q71wuBaIVuDA2Xw28liam1cx6AP2AgxnWvR7Y4+6/BTCzfwGuBs4oEO7+PeB7AGPGjPHrrrsuh5TPtHbtWgpdt9TKOTco7/zKObfvLPsJ/7NMc4vvtx07dpTdADjlOCjPKeWcm5kV9e8hl1NMG4BhZjbEzHoSXUxuSIlpAO4K01OBNeHcVwMwPXzKaQgwDFhPdGrpSjP7A4s+PDweKL8x/EREPsCyHkG4+wkzm0N0gbkCWOLuTWa2gOjiRwOwGFhqZi1ERw7Tw7pN4RNKzcAJYLa7nwTWmdkKYFNof4lwlCAiIuUhp68FdPeVwMqUtgdi00eBaWnWXQgsTGifD8zPJ1mRLvdgP3hQX+MtHwy6k1pERBKpQIiISCIVCBERSaQR5URK5LKz9nR1CoV5sF+R+8vtms13v/td/uzP/ozm5maqq6sBGD58OKtWraKqquq0EeO6ypEjR7jxxhv5xS9+kXHkumJJHUmus0eW0xGEiJSFrVu3MmrUKJ55Jhry9NixY/zmN7/hoosuOmPEuK6yZMkSPvOZz2Qcua6YUkeS6+yR5VQgRLpasf9j76a2bdvG3Llz2wtEU1MTw4cPx8zOGDEu3chy06dP59Zbb2XcuHFcdNFF7X2lk28/y5Yt49Offu/LDNONEleIdM+ZOpJcZ44spwIhImWhubmZm266iX379nH48GG2bdvGZZddljhiXNLIchB9K+zFF1/MunXrWLZsWdaxrPPp51QeF110UUm2P13uqSPJdebIcroGISJd7pVXXmHgwIGcffbZfOpTn+LnP/95+3gPqSPGpRtZ7siRI+zfv5/586Pbq2pqajh06FDa58y3n1xHrkt1/fXX88Ybb5zRvnDhwvajoky5x0eS69OnzxnzpaQCISJdbuvWre3jTk+aNInvf//77N+/n5tvvvmMEePSjSy3fft2hg0b1j7i2qZNm7jiiivSPme+/eQ6cl2qU6etMsmWe+pIcp01spxOMYkAg+dmPlctpXXqdBLAH//xH/P888+3F43UEePSjSy3ZcsWfv3rX3P06FHefvtt5s+f3350AGeO8pZvP7mOXFeITLmnjiTXmSPL6QhCRE7XBV8lsm3bNqZMmQJAr169qKmpYdu2be2ndOIjxqUbWe4HP/gBt99+O9dddx1vvvkm999/P5/4xCeA5FHeCunnhhtu4Pnnn28fXjTdyHX5OjWiXtJzpo4k16kjyxU60lBXPDSiXNco5/yKldtFf/l0boHz++bcZ+OTi4reZ7F0txHlchkx7pprrkk74lo+o7xl6mfTpk1+66235tRPPjI9Z+pIcplGluuKEeVERLpULiPG7dq1i2HDhiUuy2eUt0z9jB49mmuvvbboI9ele87UkeQ6e2Q5nWISkW7h7rvvzrg8fn2hI7L1M2PGjKLfsJfuOXv27Mmdd96Zdr7UdAQhIiKJVCBERCRRTgXCzCaa2U4zazGzuQnLe5nZ8rB8nZkNji2bF9p3mtmE0HapmW2OPd40sy8Va6NEJD/RtUzpzkrxGmYtEGZWATwC3AjUALeZWU1K2EzgkLsPBRYBD4V1a4iGHx0BTAQeNbMKd9/p7qPcfRRQC7wD/LhI2ySSt729P9fVKXSZ3r17c+DAARWJbszdOXDgQNEvnudykXos0OLuuwHMrB6YTDTO9CmTgQfD9ArgYTOz0F7v7seAPWHM6rHA87F1xwO73P2/O7IhIlKY6upqWltb+e1vf9vVqbQ7evRop9wpXIhyza137968/fbbRe3Tsv3XYGZTgYnufk+YnwGMc/c5sZjtIaY1zO8CxhEVjRfc/YnQvhhY5e4rYusuATa5+8Npnn8WMAugqqqqtr6+vqANbWtro7KysqB1S62cc4Pyzq9oub2+GQaNKl4c0HZwH5UDzi9qn8VSzq8plHd+3S23urq6je4+pqAOs90oAUwDHovNzwC+kxLTBFTH5ncBA4lOTd0Ra18MTInN9wT2A1W53LShG+W6RjnnV7Tccr1Z7X14o1w5Kuf8ultulPhGuVbgwth8NfBauhgz6wH0Aw7msO6NREcPv8khDxER6US5FIgNwDAzG2JmPYkuOjekxDQAd4XpqcCaULkagOnhU05DgGHA+th6twE/7MgGiIhIaWS9SO3uJ8xsDrAaqACWuHuTmS0gOnRpIDp1tDRchD5IVEQIcU8RXdA+Acx295MAZvYHwKeAL5Rgu0REpINy+qoNd18JrExpeyA2fZToWkXSuguBhQnt7xBdpxARkTKkO6lFRCSRCoSIiCRSgRARkUQqECIikkgFQt63NM60SMeoQIiISCIVCBERSaQCISIiiVQg5H3rgzzGg0gxqECIiEgiFQgREUmkAiEiIolUIEREJJEKhIiIJFKBEBGRRDkVCDObaGY7zazFzOYmLO9lZsvD8nVmNji2bF5o32lmE2Lt/c1shZm9bGY7zOyqYmyQiIgUR9YCYWYVwCNE40fXALeZWU1K2EzgkLsPBRYBD4V1a4hGlxsBTAQeDf0B/D3wrLv/IXAFsKPjmyMiIsWSyxHEWKDF3Xe7+3GgHpicEjMZeDxMrwDGm5mF9np3P+bue4AWYKyZ9QWuJRqqFHc/7u6/6/jmiIhIseRSIC4AXonNt4a2xBh3PwEcJhpONN26FwO/Bf7JzF4ys8fM7JyCtkBERErC3D1zgNk0YIK73xPmZwBj3f2LsZimENMa5ncRHXksAJ539ydC+2Kisa3/G3gB+IS7rzOzvwfedPevJTz/LGAWQFVVVW19fX1BG9rW1kZlZWVB65ZaOecG5Z1fxtxe3wyDRuXWUa6xefTZdnAflQPOL2qfxVLOrymUd37dLbe6urqN7j6moA7dPeMDuApYHZufB8xLiVkNXBWmewD7AUuNPRUHfATYG2u/BngmWy61tbVeqMbGxoLXLbVyzs29vPPLmNv8vrl3lGtsHn02Prmo6H0WSzm/pu7lnV93yw140bO8t6Z75HKKaQMwzMyGmFlPoovODSkxDcBdYXoqsCYk1gBMD59yGgIMA9a7+xvAK2Z2aVhnPNCcY00TEZFO0CNbgLufMLM5RP/9VwBL3L3JzBYQVaYGoovNS82sBThIVEQIcU8RvfmfAGa7+8nQ9ReBZaHo7Ab+tMjbJiIiHZC1QAC4+0qiawfxtgdi00eBaWnWXQgsTGjfDBR2XkxEREpOd1KLiEgiFQgREUmkAiEiIolUIEREJJEKhIiIJFKBEBGRRCoQIiKSSAVCREQSqUCIiEgiFQgREUmkAiEiIolUIEREJJEKhIiIJFKBEHkfGjz3ma5OQd4HVCBERCSRCoR0Pw/26+oMRD4QcioQZjbRzHaaWYuZzU1Y3svMlofl68xscGzZvNC+08wmxNr3mtk2M9tsZi8WY2NERKR4so4oZ2YVwCPAp4BWYIOZNbh7fAzpmcAhdx9qZtOBh4BbzayGaPjREcBHgZ+Z2SWxYUfr3H1/EbdHRIC9vT8HHO7qNKSby+UIYizQ4u673f04UA9MTomZDDweplcA483MQnu9ux9z9z1AS+hPRETKnLl75gCzqcBEd78nzM8Axrn7nFjM9hDTGuZ3AeOAB4EX3P2J0L4YWOXuK8xsD3AIcOC77v69NM8/C5gFUFVVVVtfX1/Qhra1tVFZWVnQuqVWzrlBGeb3+mYYNArIklssLp8+ixIHtB3cR+WA84vaZ86y9Fl2r2mKcs6vu+VWV1e30d3HFNShu2d8ANOAx2LzM4DvpMQ0AdWx+V3AQKJTU3fE2hcDU8L0R8PP84EtwLXZcqmtrfVCNTY2FrxuqZVzbu5lmN/8vu2TGXOLxeXTZ1Hi3L3xyUVF7zNnWfosu9c0RTnn191yA170LO+t6R65nGJqBS6MzVcDr6WLMbMeQD/gYKZ13f3Uz33Aj9GpJxGRspJLgdgADDOzIWbWk+iic0NKTANwV5ieCqwJlasBmB4+5TQEGAasN7NzzKwPgJmdA9wAbO/45oiISLFk/RSTu58wsznAaqACWOLuTWa2gOjQpYHo1NFSM2shOnKYHtZtMrOngGbgBDDb3U+aWRXw4+g6Nj2AJ9392RJsn4iIFChrgQBw95XAypS2B2LTR4muVSStuxBYmNK2G7gi32RFRKTz6E5qERFJpAIhIiKJVCBERCSRCoSIiCRSgRARkUQqECIikkgFQkREEqlAiIhIIhUIERFJpAIhIiKJVCBERCSRCoSIiCRSgRARkUQqECIikkgFQkREEuVUIMxsopntNLMWM5ubsLyXmS0Py9eZ2eDYsnmhfaeZTUhZr8LMXjKzpzu6ISIiUlxZC4SZVQCPADcCNcBtZlaTEjYTOOTuQ4FFwENh3Rqi0eVGABOBR0N/p9wL7OjoRoiISPHlcgQxFmhx993ufhyoByanxEwGHg/TK4DxFo0nOhmod/dj7r4HaAn9YWbVwKeBxzq+GSIiUmy5FIgLgFdi862hLTHG3U8Ah4GBWdb9NvAXwO/zzlpERErO3D1zgNk0YIK73xPmZwBj3f2LsZimENMa5ncRHSksAJ539ydC+2Kisa2PAZPc/c/N7DrgK+7+J2mefxYwC6Cqqqq2vr6+oA1ta2ujsrKyoHVLrZxzgzLM7/XNMGgUkCW3WFw+fRYlDmg7uI/KAecXtc+cZemz7F7TFOWcX3fLra6ubqO7jymoQ3fP+ACuAlbH5ucB81JiVgNXhekewH7AUmNPxQF/TXQ0sRd4A3gHeCJbLrW1tV6oxsbGgtcttXLOzb0M85vft30yY26xuHz6LEqcuzc+uajofeYsS59l95qmKOf8ultuwIue5b013SOXU0wbgGFmNsTMehJddG5IiWkA7grTU4E1IbEGYHr4lNMQYBiw3t3nuXu1uw8O/a1x9ztyrmoiIlJyPbIFuPsJM5tD9N9/BbDE3ZvMbAFRZWoAFgNLzawFOEj0pk+IewpoBk4As939ZIm2RUREiihrgQBw95VE1w7ibQ/Epo8C09KsuxBYmKHvtcDaXPIQEZHOozupRUQkkQqEiIgkUoEQEZFEKhAiIpJIBUJERBKpQIiISCIVCBERSaQCISIiiVQgREQkkQqEiIgkUoEQEZFEKhAiIpJIBUJERBKpQIiISCIVCBERSaQCISIiiXIqEGY20cx2mlmLmc1NWN7LzJaH5evMbHBs2bzQvtPMJoS23ma23sy2mFmTmX29WBskIiLFkbVAmFkF8AhwI1AD3GZmNSlhM4FD7j4UWAQ8FNatIRp+dAQwEXg09HcM+KS7XwGMAiaa2ZXF2SQpJ4PnPtPVKYhIgXI5ghgLtLj7bnc/DtQDk1NiJgOPh+kVwHgzs9Be7+7H3H0P0AKM9UhbiP9QeHgHt0VERIrI3DO/L5vZVGCiu98T5mcA49x9Tixme4hpDfO7gHHAg8AL7v5EaF8MrHL3FeFIYiMwFHjE3f8yzfPPAmYBVFVV1dbX1xe0oW1tbVRWVha0bqmVc27Qsfy2vXqYyy7oV9yEXt8Mg0YBWXKLxeXTZ1HigLaD+6gccH5R+8xZlj7fz79zpdbdcqurq9vo7mMK6tDdMz6AacBjsfkZwHdSYpqA6tj8LmAg0ampO2Lti4EpKev2BxqBkdlyqa2t9UI1NjYWvG6plXNu7h3L76K/fLp4iZwyv2/7ZMbcYnH59FmUOHdvfHJR0fvMWZY+38+/c6XW3XIDXvQs763pHrmcYmoFLozNVwOvpYsxsx5AP+BgLuu6+++AtUTXKESkXD1Y5CNBKXu5FIgNwDAzG2JmPYkuOjekxDQAd4XpqcCaULkagOnhU05DgGHAejM7z8z6A5jZ2cD1wMsd3xwRESmWHtkC3P2Emc0BVgMVwBJ3bzKzBUSHLg1Ep46WmlkL0ZHD9LBuk5k9BTQDJ4DZ7n7SzAYBj4frEGcBT7n706XYQBERKUzWAgHg7iuBlSltD8SmjxJdq0hadyGwMKVtKzA632RFRKTz6E5qERFJpAIhIiKJVCBERCSRCoSU1N7en+vqFESkQCoQIiKSSAVCREQSqUCIiEgiFQgREUmkAiEiIolUIEREJJEKhIiIJFKBEBGRRCoQIiKSSAVCREQSqUCIiEiinAqEmU00s51m1mJmcxOW9zKz5WH5OjMbHFs2L7TvNLMJoe1CM2s0sx1m1mRm9xZrg0REpDiyFogw6tsjwI1ADXCbmdWkhM0EDrn7UGAR8FBYt4ZodLkRRGNOPxr6OwH8b3cfDlwJzE7oU0REulAuRxBjgRZ33+3ux4F6YHJKzGTg8TC9AhhvZhba6939mLvvAVqAse7+urtvAnD3t4AdwAUd3xwRESmWXArEBcArsflWznwzb49x9xPAYWBgLuuG01GjgXW5py0iIqVm7p45wGwaMMHd7wnzM4iOAr4Yi2kKMa1hfhfRkccC4Hl3fyK0LwZWuvs/h/lK4BfAQnf/lzTPPwuYBVBVVVVbX19f0Ia2tbVRWVlZ0LqlVs65QQfze30zDBpV3IRifWbMLZ/nzjU2jz7bDu6jcsD5Re0zZ1n6LOg1LUWeaZTz30R3y62urm6ju48pqEN3z/gArgJWx+bnAfNSYlYDV4XpHsB+wFJjU+I+FOa/nC2HU4/a2lovVGNjY8Hrllo55+bewfzm9y1aHkl9Zswtn+fONTaPPhufXFT0PnOWpc+CXtNS5JlGOf9NdLfcgBc9x/fY1Ecup5g2AMPMbIiZ9SS66NyQEtMA3BWmpwJrQmINwPTwKachwDBgfbg+sRjY4e7fyquiiYhIp+iRLcDdT5jZHKL/9iuAJe7eZGYLiCpTA9Gb/VIzawEOEhURQtxTQDPRJ5dmu/tJM/sfwAxgm5ltDk91v7uvLPYGiohIYbIWCIDwxr0ype2B2PRRYFqadRcCC1PaniM6BSUiImVKd1KLiEgiFQgREUmkAiEiIolUIEREJJEKhIiIJFKBEBGRRCoQIiKSSAVCREQSqUCIiEgiFQgREUmkAiEiIolUIEREJJEKhIiIJFKBEJGu82C/rs5AMlCBEBGRRCoQIiKSKKcCYWYTzWynmbWY2dyE5b3MbHlYvs7MBseWzQvtO81sQqx9iZntM7PtxdgQEREprqwFwswqgEeAG4Ea4DYzq0kJmwkccvehwCLgobBuDdHwoyOAicCjoT+A74c2EREpQ7kcQYwFWtx9t7sfB+qBySkxk4HHw/QKYLyZWWivd/dj7r4HaAn94e7/TjR+tYiIlCFz98wBZlOBie5+T5ifAYxz9zmxmO0hpjXM7wLGAQ8CL7j7E6F9MbDK3VeE+cHA0+4+MsPzzwJmAVRVVdXW19cXtKFtbW1UVlYWtG6plU1ur2+GQaPOaO5Qfmn67JBYnxlzy+e5c43No8+2g/uoHHB+UfvMWZY+C3pNOzHPsvmbSNDdcqurq9vo7mMK6tDdMz6AacBjsfkZwHdSYpqA6tj8LmAg0ampO2Lti4EpsfnBwPZsOZx61NbWeqEaGxsLXrfUyia3+X0TmzuUX5o+OyTWZ8bc8nnuXGPz6LPxyUVF7zNnWfos6DXtxDzL5m8iQXfLDXjRc3yPTX3kcoqpFbgwNl8NvJYuxsx6AP2ITh/lsq6IiJShXArEBmCYmQ0xs55EF50bUmIagLvC9FRgTahcDcD08CmnIcAwYH1xUhcRkVLKWiDc/QQwB1gN7ACecvcmM1tgZjeFsMXAQDNrAb4MzA3rNgFPAc3As8Bsdz8JYGY/BJ4HLjWzVjObWdxNExGRjuiRS5C7rwRWprQ9EJs+SnStImndhcDChPbb8spUREQ6le6kFhGRRCoQIiKSSAVCREQSqUCIiEgiFQgREUmkAiEiIolUIEREJJEKhIiIJFKBEJH3lcFzn+nqFN43VCC6Kf0RiEipqUCIiEgiFYhuam/vz3V1CiJlSX8bxaMCISIiiVQgREQkkQqEiIgkyqlAmNlEM9tpZi1mNjdheS8zWx6WrzOzwbFl80L7TjObkGufIiLStbIWCDOrAB4BbgRqgNvMrCYlbCZwyN2HAouAh8K6NURDlI4AJgKPmllFjn2+b+gjqSLd3IP9ujqDLpHLEcRYoMXdd7v7caAemJwSMxl4PEyvAMabmYX2enc/5u57gJbQXy59vm/oUxUi0h3lUiAuAF6JzbeGtsSYMIb1YWBghnVz6bNrfED/UxCRTtYN3mvM3TMHmE0DJrj7PWF+BjDW3b8Yi2kKMa1hfhfRUcIC4Hl3fyK0LyYa2/qsbH3G+p4FzAqzlwI7C9zWDwP7C1y31Mo5Nyjv/JRbYco5Nyjv/Lpbbhe5+3mFdNa4Jn2mAAAGNklEQVQjh5hW4MLYfDXwWpqYVjPrAfQDDmZZN1ufALj794Dv5ZBnRmb2oruP6Wg/pVDOuUF556fcClPOuUF55/dByi2XU0wbgGFmNsTMehJddG5IiWkA7grTU4E1Hh2aNADTw6echgDDgPU59ikiIl0o6xGEu58wsznAaqACWOLuTWa2AHjR3RuAxcBSM2shOnKYHtZtMrOngGbgBDDb3U8CJPVZ/M0TEZFC5XKKCXdfSXTtIN72QGz6KDAtzboLgYW59FliHT5NVULlnBuUd37KrTDlnBuUd34fmNyyXqQWEZEPJn3VhoiIJOr2BSLcmf2SmT0d5seb2SYz22xmz5nZ0NCe99eBdGJunzez34b2zWZ2T6yPu8zsV+FxV7rnKlJ+nwz5bTezx8Mn0rDIP4R9tNXMPl7q/PLI7TozOxzbdw/E+ij617mY2V4z2xae68XQNsDMfhr2wU/N7NzQ3hX7LZ/8ymHfTTOzJjP7vZmNSYnvtK/pySc3MxtsZkdi++0fY8tqQz8t4bW3Eub3t2b2cvjd+rGZ9Y/FF2ffuXu3fgBfBp4Eng7z/wUMD9N/Dnw/Nv2PYXo6sDxM1wBbgF7AEGAXUNHJuX0eeDhh/QHA7vDz3DB9bin2HdE/C68Al4RlC4CZYXoSsAow4EpgXanzyyO3607t35T1K8JreTHQM7zGNUXIay/w4ZS2vwHmhum5wENduN/yya8c9t1wovub1gJjYu2Jf5dlkttgYHuaftYDV4XXfBVwYwlf1xuAHmH6odjrWrR9162PIMysGvg08Fis2YG+Ybof791fke/XgXRmbulMAH7q7gfd/RDwU6LvtOqwhPwGAsfc/b/C/E+BKWF6MvADj7wA9DezQaXKL8/c0unMr3OJ/249Dtwca++0/VZAful02r5z9x3unnTza5d/TU+G3BKF17avuz/v0Tv1D8i+rzuS37959M0VAC8Q3U8GRdx33bpAAN8G/gL4faztHmClmbUCM4BvhvZ8vw6kM3MDmBIOFVeY2ambCEv5lSSp+e0HPhQ7lJ7KezczdvZXpuSTG8BVZrbFzFaZ2YgsOXeUA/9mZhstussfoMrdXwcIP8/PkkMpX9d88oOu33fpdPa+yyc3gCEWnQL9hZldE8u5tQS55ZLf3URHLKfyKMq+67YFwsz+BNjn7htTFt0HTHL3auCfgG+dWiWhG8/Q3pm5/T9gsLtfDvyM9/7bK3pu6fIL//FMBxaZ2XrgLaJ7VzLl0Sn7Lktum4i+SuAK4DvAv2bJuaM+4e4fJ/om4tlmdm2G2E7bbzH55Kd9V1hurwMfc/fRhFOhZta3hLllzM/Mvkr097DsVFOaPPLOr9sWCOATwE1mtpfoUOmTZvYMcIW7rwsxy4Grw3T7135Y7l8H0im5ufsBdz8W2v8vUJuacxFzS5ffE+HQ+Bp3Hwv8O/CrLHl01r5Lm5u7v+nubWF6JdGRxodLlBvu/lr4uQ/4MdFh+2/C6YVTpxn2hfDO3G9551cm+y6dTt13+eQWTt0cCNMbic7rXxJyq46Flvp1xaIPOPwJcHv4RwqKue+KcQGlqx+Ei21EN/7t572LmTOBfw7Tszn9IvVTYXoEp1/Q2U2RLlLnkdugWPwtwAthegCwh+hC5rlhekAp9l2YPj/87AX8HPhkmP80p19sXd8Z+eWY20d4736escCvQ549wms5hPcuyI3oYD7nAH1i078kunbwt5x+EfhvumK/FZBfl++72PK1nH4hOPHvskxyO4/wHkF0wffVU68f0dcIXcl7F6knlfB1nUj0LRXnpcQXbd8V7Y2mKx+c/kZyC7AtbPxa4OLQ3hv4EdEFm/Wn2sOyrxL9F7CTIn3qIM/c/hpoCu2NwB/G1r875NwC/GmJ993fAjvCfvhSLMaIBnjaFfIf0xn55ZjbnNi+ewG4OrZsEtEnx3YBXy1CPheH59kSnvOroX0gUdH6Vfh56s2iU/dbAfmVw767heg/22PAb4DVsXUS/y67OjeiD0ic2m+bgM/E+hoDbA+5PUwowCXKr4XomsLm8PjHYu873UktIiKJuvM1CBERKSEVCBERSaQCISIiiVQgREQkkQqEiIgkUoEQEZFEKhAiIpJIBUJERBL9f0K5qHmecHfGAAAAAElFTkSuQmCC\n"
},
"metadata": {
"needs_background": "light"
}
},
{
"output_type": "execute_result",
"execution_count": 110,
"data": {
"text/plain": "0.0004201354102465936"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "bins",
"execution_count": 111,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 111,
"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. ])"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "",
"execution_count": null,
"outputs": []
}
],
"metadata": {
"kernelspec": {
"name": "python3",
"display_name": "Python 3",
"language": "python"
},
"hide_input": false,
"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
},
"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"
},
"gist": {
"id": "",
"data": {
"description": "histogramからKL情報量",
"public": true
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment