Skip to content

Instantly share code, notes, and snippets.

@sshojiro

sshojiro/ab-test.ipynb

Last active May 9, 2019 17:02
Show Gist options
  • Save sshojiro/7fed028defbb3fc30abd66462aee00c5 to your computer and use it in GitHub Desktop.
Save sshojiro/7fed028defbb3fc30abd66462aee00c5 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
ab-test.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
bayesian-gaussian-mixture-model.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
bayesian-least-squares.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
challenger_causal_analysis.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
recieving-messages.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# メールの受信数の推定\n",
"\n",
"PyMC and Poisson distribution\n",
"\n",
"[Pythonで体験するベイズ推論](https://www.amazon.co.jp/Python%E3%81%A7%E4%BD%93%E9%A8%93%E3%81%99%E3%82%8B%E3%83%99%E3%82%A4%E3%82%BA%E6%8E%A8%E8%AB%96-PyMC%E3%81%AB%E3%82%88%E3%82%8BMCMC%E5%85%A5%E9%96%80-%E3%82%AD%E3%83%A3%E3%83%A1%E3%83%AD%E3%83%B3-%E3%83%87%E3%83%93%E3%83%83%E3%83%89%E3%82%BD%E3%83%B3-%E3%83%94%E3%83%AD%E3%83%B3/dp/4627077912)より"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"('txtdata.csv', <http.client.HTTPMessage at 0x245cf6a8a20>)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from urllib.request import urlretrieve\n",
"# urlretrieve('https://git.io/vXTVC','txtdata.csv')"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pymc as pm\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<BarContainer object of 74 artists>"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEGlJREFUeJzt3V2MXHd9xvHv07wUCNAkZJNaMVsHyQrhonFgSYNSIXCAAkVJLghKhJBVufINRURFgNNKrZC4CELiRWpFZRHASDQvhKSOIgpETqKKXoQ4L0CCSQ1pCFaMbWii8CIBhl8v5hgWs2Zmd2c8Z/77/UirM+fsGc/jnfGz//mfOcepKiRJs++Pph1AkjQeFrokNcJCl6RGWOiS1AgLXZIaYaFLUiMsdElqhIUuSY2w0CWpESefyAc766yzasOGDSfyISVp5j3wwAM/rKq5Yfud0ELfsGEDe/bsOZEPKUkzL8n3RtnPKRdJaoSFLkmNsNAlqREWuiQ1wkKXpEZY6JLUCAtdkhphoUtSIyx0SWrECT1TVLPplR+++ze373/v5ikmkfSHOEKXpEYMLfQk5yd5eNHXs0muTXJmkruS7OuWZ5yIwJKkpQ0t9Kp6rKo2VdUm4BXAz4Dbge3A7qraCOzu1iVJU7LcKZfLgO9W1feAK4Cd3fadwJXjDCZJWp7lHhS9Grixu31OVR0AqKoDSc5e6g5JtgHbAObn51eaU1oTPACt1Rh5hJ7kVOBy4PPLeYCq2lFVC1W1MDc39PrskqQVWs6Uy5uAB6vqYLd+MMk6gG55aNzhJEmjW06hX8Nvp1sA7gC2dLe3ALvGFUqStHwjFXqS5wGvB25btPl64PVJ9nXfu3788SRJoxrpoGhV/Qx40THbfsTgUy+SpB7wTFFJaoSFLkmNsNAlqREWuiQ1wkKXpEZY6JLUCAtdkhphoUtSIyx0SWqEhS5JjbDQJakRFrokNcJCl6RGWOiS1AgLXZIaYaFLUiMsdElqhIUuSY2w0CWpERa6JDVipEJPcnqSW5N8O8neJK9KcmaSu5Ls65ZnTDqsJOn4Rh2hfxz4UlW9FLgQ2AtsB3ZX1UZgd7cuSZqSoYWe5IXAq4EbAKrqF1X1DHAFsLPbbSdw5aRCSpKGG2WE/hLgMPDpJA8l+WSS04BzquoAQLc8e6k7J9mWZE+SPYcPHx5bcEnS7xql0E8GXg58oqouAn7KMqZXqmpHVS1U1cLc3NwKY0qShhml0PcD+6vqvm79VgYFfzDJOoBueWgyESVJoxha6FX1A+D7Sc7vNl0GfAu4A9jSbdsC7JpIQknSSE4ecb93AZ9LcirwOPA3DH4Z3JJkK/AkcNVkIkqSRjFSoVfVw8DCEt+6bLxxJEkr5ZmiktQIC12SGmGhS1IjRj0oqgl55Yfv/s3t+9+7eYpJJM06R+iS1AgLXZIaYaFLUiMsdElqhIUuSY2w0CWpERa6JDXCQpekRljoktQIC12SGuGp/9IM89IRWswRuiQ1wkKXpEZY6JLUCAtdkhox0kHRJE8APwZ+BRypqoUkZwI3AxuAJ4C3VdXTk4kpSRpmOSP011bVpqo6+p9Fbwd2V9VGYHe3LkmaktVMuVwB7Oxu7wSuXH0cSdJKjVroBXwlyQNJtnXbzqmqAwDd8uxJBJQkjWbUE4suraqnkpwN3JXk26M+QPcLYBvA/Pz8CiJKkkYx0gi9qp7qloeA24GLgYNJ1gF0y0PHue+OqlqoqoW5ubnxpJYk/Z6hhZ7ktCQvOHobeAPwCHAHsKXbbQuwa1IhJUnDjTLlcg5we5Kj+/97VX0pyf3ALUm2Ak8CV00uprQ2ea0WLcfQQq+qx4ELl9j+I+CySYSSJC2fZ4pKUiMsdElqhIUuSY2w0CWpERa6JDXCQpekRljoktQIC12SGmGhS1IjLHRJaoSFLkmNsNAlqREWuiQ1wkKXpEZY6JLUCAtdkhphoUtSIyx0SWqEhS5JjbDQJakRIxd6kpOSPJTkzm79vCT3JdmX5OYkp04upiRpmOWM0N8N7F20/iHgo1W1EXga2DrOYJKk5Rmp0JOsB/4a+GS3HmAzcGu3y07gykkElCSNZtQR+seA9wG/7tZfBDxTVUe69f3AuWPOJklahqGFnuQtwKGqemDx5iV2rePcf1uSPUn2HD58eIUxJUnDjDJCvxS4PMkTwE0Mplo+Bpye5ORun/XAU0vduap2VNVCVS3Mzc2NIbIkaSlDC72qrquq9VW1AbgauLuq3g7cA7y1220LsGtiKSVJQ508fJfjej9wU5IPAg8BN4wnkiSdWK/88N2/s37/ezdPKcnqLKvQq+pe4N7u9uPAxeOPJElaCc8UlaRGWOiS1AgLXZIasZqDopLUrMUHSmflIKkjdElqhIUuSY2w0CWpERa6JDXCQpekRqzpT7nM4lFsqe/8dzU9jtAlqREWuiQ1wkKXpEZY6JLUCAtdkhphoUtSIyx0SWqEhS5JjVjTJxZpZTxxROonR+iS1IihhZ7kOUm+luTrSR5N8oFu+3lJ7kuyL8nNSU6dfFxJ0vGMMuXyc2BzVf0kySnAV5P8J/D3wEer6qYk/wZsBT4xwaySVmnxdBnM7pSZ035LGzpCr4GfdKundF8FbAZu7bbvBK6cSEJJ0khGmkNPclKSh4FDwF3Ad4FnqupIt8t+4NzJRJQkjWKkT7lU1a+ATUlOB24HLlhqt6Xum2QbsA1gfn5+hTEl6ficghlY1qdcquoZ4F7gEuD0JEd/IawHnjrOfXZU1UJVLczNza0mqyTpDxg6Qk8yB/yyqp5J8lzgdcCHgHuAtwI3AVuAXZMMqgFHIpKOZ5Qpl3XAziQnMRjR31JVdyb5FnBTkg8CDwE3TDCnJGmIoYVeVd8ALlpi++PAxZMIJUlaPs8UlaRGWOiS1AgLXZIa4dUW1Ut+mkfj1MolD4ZxhC5JjbDQJakRFrokNcJCl6RGWOiS1AgLXZIaYaFLUiMsdElqxMyeWDTsxJO1ciKBJB3lCF2SGmGhS1IjZnbKRW3x2i1rh9Ohk+MIXZIaYaFLUiOcchmjY99KHqsvby2d3miXz21/TOO5cIQuSY0YOkJP8mLgs8CfAr8GdlTVx5OcCdwMbACeAN5WVU9PLqqk5Rr2rlFtGWWEfgR4T1VdAFwCvDPJy4DtwO6q2gjs7tYlSVMytNCr6kBVPdjd/jGwFzgXuALY2e22E7hyUiElScMt66Bokg3ARcB9wDlVdQAGpZ/k7OPcZxuwDWB+fn41WbWGebCvP3wu+mvkg6JJng98Abi2qp4d9X5VtaOqFqpqYW5ubiUZJUkjGKnQk5zCoMw/V1W3dZsPJlnXfX8dcGgyESVJoxjlUy4BbgD2VtVHFn3rDmALcH233DWRhGvcct/ezspp1b5tH5j2z2FWXi8azShz6JcC7wC+meThbts/MCjyW5JsBZ4ErppMREnSKIYWelV9Fchxvn3ZeONIklbKU//1e6Y9DaB+8fUwOzz1X5IaYaFLUiOcclETnBbor2HXk/GTNuPjCF2SGmGhS1IjnHKRjsNpHM0aR+iS1AgLXZIa4ZTLMvgWfG1b7XV1juVrSOPmCF2SGmGhS1IjnHJZxCmVtcXnW6vRx9ePI3RJasTMjND7+NtQo5nGqd2eTr62TaMv+tBRjtAlqREWuiQ1YmamXLQyo0w99OGtomaXr5/+cIQuSY0YWuhJPpXkUJJHFm07M8ldSfZ1yzMmG1OSNMwoUy6fAf4F+OyibduB3VV1fZLt3fr7xx9vuoa9lVzJW81J/JnSOPkanF1DR+hV9V/A/x2z+QpgZ3d7J3DlmHNJkpZppXPo51TVAYBuefb4IkmSVmLin3JJsg3YBjA/Pz/ph5OW5IlGWgtWOkI/mGQdQLc8dLwdq2pHVS1U1cLc3NwKH06SNMxKC/0OYEt3ewuwazxxJEkrNXTKJcmNwGuAs5LsB/4ZuB64JclW4EngqkmGVL/5qQipH4YWelVdc5xvXTbmLJKkVfBMUUlqRDPXchn2/zdKUuscoUtSIyx0SWpEM1Muo/DTGCfGWpn+Wqv/K476yxG6JDXCQpekRqypKZcWrdW34JP4e6/Vn6Xa4QhdkhrhCF2SRjAL7+AcoUtSIyx0SWqEUy6SemcWpjf6yBG6JDXCQpekRljoktQIC12SGmGhS1Ij/JSLNCZ+MmNy/NmOxhG6JDViVYWe5I1JHkvynSTbxxVKkrR8Ky70JCcB/wq8CXgZcE2Sl40rmCRpeVYzQr8Y+E5VPV5VvwBuAq4YTyxJ0nKtptDPBb6/aH1/t02SNAWpqpXdMbkK+Kuq+ttu/R3AxVX1rmP22wZs61bPBx5beVwAzgJ+uMo/Y9LMOB6zkBFmI6cZx2NaGf+squaG7bSajy3uB168aH098NSxO1XVDmDHKh7ndyTZU1UL4/rzJsGM4zELGWE2cppxPPqecTVTLvcDG5Ocl+RU4GrgjvHEkiQt14pH6FV1JMnfAV8GTgI+VVWPji2ZJGlZVnWmaFV9EfjimLKMamzTNxNkxvGYhYwwGznNOB69zrjig6KSpH7x1H9JasTMFHpfLzOQ5FNJDiV5ZNG2M5PclWRftzxjyhlfnOSeJHuTPJrk3X3LmeQ5Sb6W5Otdxg90289Lcl+X8ebuAPxUJTkpyUNJ7uxjxiRPJPlmkoeT7Om29ea5XpTz9CS3Jvl299p8VZ9yJjm/+xke/Xo2ybV9ynismSj0nl9m4DPAG4/Zth3YXVUbgd3d+jQdAd5TVRcAlwDv7H5+fcr5c2BzVV0IbALemOQS4EPAR7uMTwNbp5jxqHcDexet9zHja6tq06KP2PXpuT7q48CXquqlwIUMfqa9yVlVj3U/w03AK4CfAbf3KePvqarefwGvAr68aP064Lpp51qUZwPwyKL1x4B13e11wGPTznhM3l3A6/uaE3ge8CDwFwxO4jh5qdfBlLKtZ/CPeDNwJ5AeZnwCOOuYbb16roEXAv9LdxyvrzkX5XoD8N99zlhVszFCZ/YuM3BOVR0A6JZnTznPbyTZAFwE3EfPcnZTGQ8Dh4C7gO8Cz1TVkW6XPjzvHwPeB/y6W38R/ctYwFeSPNCdqQ09e66BlwCHgU9301efTHIa/ct51NXAjd3tvmacmULPEtv8eM4yJXk+8AXg2qp6dtp5jlVVv6rB29v1DC7+dsFSu53YVL+V5C3Aoap6YPHmJXad9mvz0qp6OYMpyncmefWU8yzlZODlwCeq6iLgp/Rp6mKR7pjI5cDnp51lmFkp9JEuM9AjB5OsA+iWh6achySnMCjzz1XVbd3m3uUEqKpngHsZzPefnuTo+RLTft4vBS5P8gSDq4tuZjBi71NGquqpbnmIwZzvxfTvud4P7K+q+7r1WxkUfN9ywuAX44NVdbBb72NGYHYKfdYuM3AHsKW7vYXBnPXUJAlwA7C3qj6y6Fu9yZlkLsnp3e3nAq9jcJDsHuCt3W5TzVhV11XV+qrawOA1eHdVvZ0eZUxyWpIXHL3NYO73EXr0XANU1Q+A7yc5v9t0GfAtepazcw2/nW6BfmYcmPYk/jIOSrwZ+B8G86r/OO08i3LdCBwAfslg1LGVwbzqbmBftzxzyhn/ksE0wDeAh7uvN/cpJ/DnwENdxkeAf+q2vwT4GvAdBm95/3jaz3mX6zXAnX3L2GX5evf16NF/K316rhdl3QTs6Z7z/wDO6FtOBgfofwT8yaJtvcq4+MszRSWpEbMy5SJJGsJCl6RGWOiS1AgLXZIaYaFLUiMsdElqhIUuSY2w0CWpEf8Pl/EiriTJfEYAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"count_data = np.loadtxt('data/txtdata.csv')\n",
"n_count_data = len(count_data)\n",
"plt.bar(np.arange(n_count_dta),count_data,color=\"#348ABD\")"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"alpha = 1/ count_data.mean()\n",
"l1 = pm.Exponential(\"lambda_1\", alpha)\n",
"l2 = pm.Exponential(\"lambda_2\", alpha)\n",
"tau = pm.DiscreteUniform(\"tau\", lower=0, upper=n_count_data)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# Model definition\n",
"@pm.deterministic\n",
"def lambda_func(tau=tau, lambda_1=l1, lambda_2=l2):\n",
" out=np.zeros(n_count_data)\n",
" out[:tau] = lambda_1\n",
" out[tau:] = lambda_2\n",
" return out "
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" [-----------------100%-----------------] 40000 of 40000 complete in 15.8 sec"
]
}
],
"source": [
"observation = pm.Poisson(\"obs\", lambda_func, \n",
" value=count_data, observed=True)\n",
"model = pm.Model([observation, l1, l2, tau])\n",
"mcmc = pm.MCMC(model)\n",
"sample_end = 40000\n",
"sample_begin = 10000\n",
"mcmc.sample(sample_end, sample_begin)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"l1_post = mcmc.trace('lambda_1')[:]\n",
"l2_post = mcmc.trace('lambda_2')[:]\n",
"tau_post = mcmc.trace('tau')[:]"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADwJJREFUeJzt3X+s3Xddx/Hni3YDw6+BvdNlrbZqJzSGsHGtJKD8EEy3mdY/iHSROHWhCXEoP0RLMBPnP/thJBqn2MDCD5FZEKFxJQNxSGIc9A62ua5WrrPYa6crP5UQGdW3f5wzPNyee8/33p7bc/rp85Gc7Hy/388995Wz+3n1e7/n+/3eVBWSpHY9YdIBJElry6KXpMZZ9JLUOItekhpn0UtS4yx6SWrcyKJPcnuSR5M8uMT2JPmDJPNJHkhyxfhjSpJWq8se/buAHctsvxLY2n/sAf74zGNJksZlZNFX1aeALy8zZBfwnuq5B7goySXjCihJOjPrx/AalwLHB5YX+useWTwwyR56e/08+clPft6znvWsMXx7STp/3HvvvV+sqpmVfM04ij5D1g29r0JV7QP2AczOztbc3NwYvr0knT+SfGGlXzOOs24WgE0DyxuBE2N4XUnSGIyj6A8AP98/++b5wNeq6rTDNpKkyRh56CbJ+4EXAxuSLAC/BVwAUFVvBw4CVwHzwDeAX1yrsJKklRtZ9FV1zYjtBfzy2BJJksbKK2MlqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcSP/wpRWb/PeO5fcduymq89iEknnM/foJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXGdij7JjiRHk8wn2Ttk+/cluTvJ55I8kOSq8UeVJK3GyKJPsg64DbgS2AZck2TbomG/CeyvqsuB3cAfjTuoJGl1uuzRbwfmq+rhqnoMuAPYtWhMAU/rP386cGJ8ESVJZ6JL0V8KHB9YXuivG/RW4FVJFoCDwGuHvVCSPUnmksydPHlyFXElSSvVpegzZF0tWr4GeFdVbQSuAt6b5LTXrqp9VTVbVbMzMzMrTytJWrEuRb8AbBpY3sjph2auA/YDVNXfA08CNowjoCTpzHQp+kPA1iRbklxI78PWA4vG/CvwkwBJnk2v6D02I0lTYGTRV9Up4HrgLuAIvbNrDie5McnO/rA3Aq9Ocj/wfuAXqmrx4R1J0gSs7zKoqg7S+5B1cN0NA88fAl4w3miSpHHwylhJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktS49ZMOcL7avPfOZbcfu+nqs5REUuvco5ekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXGeR38GRp0LL0nTwD16SWqcRS9JjbPoJalxnYo+yY4kR5PMJ9m7xJifTfJQksNJ/my8MSVJqzXyw9gk64DbgJcDC8ChJAeq6qGBMVuBNwMvqKqvJLl4rQJLklamyx79dmC+qh6uqseAO4Bdi8a8Gritqr4CUFWPjjemJGm1uhT9pcDxgeWF/rpBlwGXJfm7JPck2THshZLsSTKXZO7kyZOrSyxJWpEuRZ8h62rR8npgK/Bi4BrgHUkuOu2LqvZV1WxVzc7MzKw0qyRpFboU/QKwaWB5I3BiyJiPVNW3qupfgKP0il+SNGFdiv4QsDXJliQXAruBA4vGfBh4CUCSDfQO5Tw8zqCSpNUZWfRVdQq4HrgLOALsr6rDSW5MsrM/7C7gS0keAu4G3lRVX1qr0JKk7jrd66aqDgIHF627YeB5AW/oPyRJU8QrYyWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZ1KvokO5IcTTKfZO8y416RpJLMji+iJOlMjCz6JOuA24ArgW3ANUm2DRn3VOBXgE+PO6QkafW67NFvB+ar6uGqegy4A9g1ZNzvALcA/z3GfJKkM9Sl6C8Fjg8sL/TXfVuSy4FNVfVXy71Qkj1J5pLMnTx5csVhJUkr16XoM2RdfXtj8gTgbcAbR71QVe2rqtmqmp2ZmemeUpK0al2KfgHYNLC8ETgxsPxU4EeATyY5BjwfOOAHspI0HboU/SFga5ItSS4EdgMHHt9YVV+rqg1VtbmqNgP3ADuram5NEkuSVmRk0VfVKeB64C7gCLC/qg4nuTHJzrUOKEk6M+u7DKqqg8DBRetuWGLsi888liRpXLwyVpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxnU6j15n3+a9dy67/dhNV5+lJJLOde7RS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY3zb8aOMOpvt0rStHOPXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDWuU9En2ZHkaJL5JHuHbH9DkoeSPJDkE0m+f/xRJUmrMbLok6wDbgOuBLYB1yTZtmjY54DZqnoO8EHglnEHlSStTpc9+u3AfFU9XFWPAXcAuwYHVNXdVfWN/uI9wMbxxpQkrVaXor8UOD6wvNBft5TrgI8O25BkT5K5JHMnT57snlKStGpdij5D1tXQgcmrgFng1mHbq2pfVc1W1ezMzEz3lJKkVetyr5sFYNPA8kbgxOJBSV4GvAV4UVV9czzxJElnqsse/SFga5ItSS4EdgMHBgckuRz4E2BnVT06/piSpNUaWfRVdQq4HrgLOALsr6rDSW5MsrM/7FbgKcAHktyX5MASLydJOss63aa4qg4CBxetu2Hg+cvGnEuSNCZeGStJjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUuPWTDqDV2bz3ziW3Hbvp6rOYRNK0c49ekhpn0UtS4yx6SWrceX+Mfrlj3ZLUAvfoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcZ2ujE2yA/h9YB3wjqq6adH2JwLvAZ4HfAl4ZVUdG29UdTXqal/vbimdX0YWfZJ1wG3Ay4EF4FCSA1X10MCw64CvVNUPJdkN3Ay8ci0Cr4a3OZB0PuuyR78dmK+qhwGS3AHsAgaLfhfw1v7zDwJ/mCRVVWPMqjHxXvbS+aVL0V8KHB9YXgB+bKkxVXUqydeA7wa+ODgoyR5gT3/x60mOrib0CBsWf99zwNRkzs2dh05N5o7Otbxg5rPhXMsL8MMr/YIuRZ8h6xbvqXcZQ1XtA/Z1+J6rlmSuqmbX8nuMm5nX3rmWF8x8NpxreaGXeaVf0+WsmwVg08DyRuDEUmOSrAeeDnx5pWEkSePXpegPAVuTbElyIbAbOLBozAHg2v7zVwB/4/F5SZoOIw/d9I+5Xw/cRe/0ytur6nCSG4G5qjoAvBN4b5J5envyu9cy9AhremhojZh57Z1recHMZ8O5lhdWkTnueEtS27wyVpIaZ9FLUuPO6aJPcnuSR5M8uGj9a5McTXI4yS2TyjfMsMxJnpvkniT3JZlLsn2SGQcl2ZTk7iRH+u/nr/bXPzPJx5N8vv/fZ0w66+OWyXxrkn9M8kCSv0xy0aSzwtJ5B7b/WpJKsmFSGRdbLvO0zr9lfi6mcv4leVKSzyS5v5/3t/vrtyT5dH/u/Xn/JJnlVdU5+wB+ArgCeHBg3UuAvwae2F++eNI5O2T+GHBl//lVwCcnnXMg2yXAFf3nTwX+CdgG3ALs7a/fC9w86awdMv8UsL6//uZpybxU3v7yJnonQnwB2DDprB3e46mdf8tknsr5R+/6pKf0n18AfBp4PrAf2N1f/3bgNaNe65zeo6+qT3H6+fqvAW6qqm/2xzx61oMtY4nMBTyt//zpnH6dwsRU1SNV9dn+8/8CjtC7EnoX8O7+sHcDPzOZhKdbKnNVfayqTvWH3UPvmpCJW+Y9Bngb8OsMuQBxkpbJPLXzb5nMUzn/qufr/cUL+o8CXkrvVjPQce6d00W/hMuAH+//avO3SX500oE6eB1wa5LjwO8Cb55wnqGSbAYup7dn8T1V9Qj0JhBw8eSSLW1R5kG/BHz0bOcZZTBvkp3Av1XV/RMNNcKi9/icmH+LMk/t/EuyLsl9wKPAx4F/Br46sMOywP/vFCypxaJfDzyD3q84bwL2Jxl2i4Zp8hrg9VW1CXg9vesSpkqSpwB/Abyuqv5z0nm6WCpzkrcAp4D3TSrbMIN56eV7C3DDREONMOQ9nvr5NyTz1M6/qvqfqnouvd8+twPPHjZs1Ou0WPQLwIf6v/Z8BvhfejcummbXAh/qP/8Avf+hUyPJBfQmxvuq6vGc/5Hkkv72S+jtcUyNJTKT5Frgp4Gfq/5BzmkwJO8PAluA+5McozfRP5vkeyeX8jst8R5P9fxbIvNUzz+Aqvoq8El6/4Be1L/VDAy/Jc1pWiz6D9M7hkWSy4ALmf67050AXtR//lLg8xPM8h36e2PvBI5U1e8NbBq87cW1wEfOdralLJU5vT+g8xvAzqr6xqTyLTYsb1X9Q1VdXFWbq2ozvQK9oqr+fYJRv22Zn4upnX/LZJ7K+Zdk5vEzw5J8F/Ayep8r3E3vVjPQde5N+pPlM/xU+v3AI8C36E2E6+j9YP0p8CDwWeClk87ZIfMLgXuB++kdM3zepHMO5H0hvV8NHwDu6z+uoncb6k/QmxSfAJ456awdMs/Tu5324+vePumsy+VdNOYY03XWzVLv8dTOv2UyT+X8A54DfK6f90Hghv76HwA+0/95/gD9M5yWe3gLBElqXIuHbiRJAyx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1Lj/AzbSMd16yV2zAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADutJREFUeJzt3X2MZXddx/H3h90uGJ4K7lSb7uquuhU2htAyriSgPIhkW8yufxCzG4lVGzYhFuVBdAmmYv2ntEYSYxU30vAgUgsibGRJQSySGAudQlu7rStjLe641S7yoIRIWfz6xz3Fy+ydmTOzd/be+fF+JTc953d+c+eT2zmfPXPuPWdSVUiS2vW4SQeQJK0vi16SGmfRS1LjLHpJapxFL0mNs+glqXErFn2Sm5M8kuS+JbYnye8nmU9yb5LLxx9TkrRWfY7o3wHsXWb7FcCu7nEI+KNzjyVJGpcVi76qPgl8cZkp+4F31cAdwIVJLh5XQEnSudk8hue4BDg5tL7QjT28eGKSQwyO+nniE5/4nGc84xlj+PaS9J3jrrvu+kJVzazma8ZR9BkxNvK+ClV1BDgCMDs7W3Nzc2P49pL0nSPJ51f7NeP41M0CsH1ofRtwagzPK0kag3EU/VHg57tP3zwX+EpVnXXaRpI0GSueuknyXuCFwNYkC8BvARcAVNXbgGPAlcA88DXgF9crrCRp9VYs+qo6uML2An55bIkkSWPllbGS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjehV9kr1JTiSZT3J4xPbvS3J7ks8muTfJleOPKklaixWLPskm4CbgCmA3cDDJ7kXTfhO4taouAw4AfzjuoJKktelzRL8HmK+qB6vqUeAWYP+iOQU8pVt+KnBqfBElSeeiT9FfApwcWl/oxoa9GXhFkgXgGPDqUU+U5FCSuSRzp0+fXkNcSdJq9Sn6jBirResHgXdU1TbgSuDdSc567qo6UlWzVTU7MzOz+rSSpFXrU/QLwPah9W2cfWrmauBWgKr6e+AJwNZxBJQknZs+RX8nsCvJziRbGLzZenTRnH8FfhIgyTMZFL3nZiRpCqxY9FV1BrgGuA14gMGna44nuS7Jvm7a64FXJrkHeC/wC1W1+PSOJGkCNveZVFXHGLzJOjx27dDy/cDzxhtNkjQOXhkrSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXG9ij7J3iQnkswnObzEnJ9Ncn+S40n+bLwxJUlrtXmlCUk2ATcBPwUsAHcmOVpV9w/N2QW8EXheVX0pyUXrFViStDp9juj3APNV9WBVPQrcAuxfNOeVwE1V9SWAqnpkvDElSWvVp+gvAU4OrS90Y8MuBS5N8ndJ7kiyd9QTJTmUZC7J3OnTp9eWWJK0Kn2KPiPGatH6ZmAX8ELgIPAnSS4864uqjlTVbFXNzszMrDarJGkN+hT9ArB9aH0bcGrEnA9V1Teq6l+AEwyKX5I0YX2K/k5gV5KdSbYAB4Cji+Z8EHgRQJKtDE7lPDjOoJKktVnxUzdVdSbJNcBtwCbg5qo6nuQ6YK6qjnbbXprkfuCbwBuq6j/XM7g0TjsOf7jXvIeuf9k6J5HGb8WiB6iqY8CxRWPXDi0X8LruIUmaIl4ZK0mN63VEL2mg7yke8DSPpodH9JLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1LheRZ9kb5ITSeaTHF5m3suTVJLZ8UWUJJ2LFYs+ySbgJuAKYDdwMMnuEfOeDPwK8Klxh5QkrV2fI/o9wHxVPVhVjwK3APtHzPsd4Abgf8aYT5J0jvoU/SXAyaH1hW7sW5JcBmyvqr9a7omSHEoyl2Tu9OnTqw4rSVq9PkWfEWP1rY3J44C3Aq9f6Ymq6khVzVbV7MzMTP+UkqQ161P0C8D2ofVtwKmh9ScDPwJ8IslDwHOBo74hK0nToU/R3wnsSrIzyRbgAHD0sY1V9ZWq2lpVO6pqB3AHsK+q5tYlsSRpVVYs+qo6A1wD3AY8ANxaVceTXJdk33oHlCSdm819JlXVMeDYorFrl5j7wnOPJUkal15FL21EOw5/eNIRpKngLRAkqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGea8baZ30vdfOQ9e/bJ2T6DudR/SS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjetV9En2JjmRZD7J4RHbX5fk/iT3Jvl4ku8ff1RJ0lqsWPRJNgE3AVcAu4GDSXYvmvZZYLaqngW8H7hh3EElSWvT54h+DzBfVQ9W1aPALcD+4QlVdXtVfa1bvQPYNt6YkqS16lP0lwAnh9YXurGlXA18ZNSGJIeSzCWZO336dP+UkqQ161P0GTFWIycmrwBmgRtHba+qI1U1W1WzMzMz/VNKktZsc485C8D2ofVtwKnFk5K8BHgT8IKq+vp44kmSzlWfI/o7gV1JdibZAhwAjg5PSHIZ8MfAvqp6ZPwxJUlrtWLRV9UZ4BrgNuAB4NaqOp7kuiT7umk3Ak8C3pfk7iRHl3g6SdJ51ufUDVV1DDi2aOzaoeWXjDmXJGlMvDJWkhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1LheF0xJ02TH4Q9POoK0oXhEL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOD91I03Yaj5F9ND1L1vHJGqVR/SS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOP/ClKbGav7SkqT+LHppA+n7j6F/clDDPHUjSY3rVfRJ9iY5kWQ+yeER2x+f5M+77Z9KsmPcQSVJa7Ni0SfZBNwEXAHsBg4m2b1o2tXAl6rqh4C3Am8Zd1BJ0tr0OUe/B5ivqgcBktwC7AfuH5qzH3hzt/x+4A+SpKpqjFm1AfkG62R4Ll/D+hT9JcDJofUF4MeWmlNVZ5J8Bfhu4AvDk5IcAg51q19NcmItoVewdfH33QDMvP42Wl44D5kz/t+9N9rrvNHyAvzwar+gT9FnxNjiI/U+c6iqI8CRHt9zzZLMVdXsen6PcTPz+ttoecHM58NGywuDzKv9mj5vxi4A24fWtwGnlpqTZDPwVOCLqw0jSRq/PkV/J7Aryc4kW4ADwNFFc44CV3XLLwf+xvPzkjQdVjx1051zvwa4DdgE3FxVx5NcB8xV1VHg7cC7k8wzOJI/sJ6hV7Cup4bWiZnX30bLC2Y+HzZaXlhD5njgLUlt88pYSWqcRS9JjdvQRZ/k5iSPJLlv0firu1s2HE9yw6TyjTIqc5JnJ7kjyd1J5pLsmWTGYUm2J7k9yQPd6/mr3fjTk3wsyee6/z5t0lkfs0zmG5P8Y5J7k/xlkgsnnRWWzju0/deSVJKtk8q42HKZp3X/W+bnYir3vyRPSPLpJPd0eX+7G9/Z3Wrmc92tZ7as+GRVtWEfwE8AlwP3DY29CPhr4PHd+kWTztkj80eBK7rlK4FPTDrnULaLgcu75ScD/8TgVhg3AIe78cPAWyadtUfmlwKbu/G3TEvmpfJ269sZfBDi88DWSWft8RpP7f63TOap3P8YXJ/0pG75AuBTwHOBW4ED3fjbgFet9Fwb+oi+qj7J2Z/XfxVwfVV9vZvzyHkPtowlMhfwlG75qZx9ncLEVNXDVfWZbvm/gQcYXAm9H3hnN+2dwM9MJuHZlspcVR+tqjPdtDsYXBMyccu8xjC4d9SvM+ICxElaJvPU7n/LZJ7K/a8GvtqtXtA9Cngxg1vNQM99b0MX/RIuBX68+9Xmb5P86KQD9fAa4MYkJ4HfBd444TwjdXclvYzBkcX3VNXDMNiBgIsml2xpizIP+yXgI+c7z0qG8ybZB/xbVd0z0VArWPQab4j9b1Hmqd3/kmxKcjfwCPAx4J+BLw8dsCzw/wcFS2qx6DcDT2PwK84bgFuTjLpFwzR5FfDaqtoOvJbBdQlTJcmTgL8AXlNV/zXpPH0slTnJm4AzwHsmlW2U4bwM8r0JuHaioVYw4jWe+v1vROap3f+q6ptV9WwGv33uAZ45atpKz9Ni0S8AH+h+7fk08L8Mblw0za4CPtAtv4/B/9CpkeQCBjvGe6rqsZz/keTibvvFDI44psYSmUlyFfDTwM9Vd5JzGozI+4PATuCeJA8x2NE/k+R7J5fy2y3xGk/1/rdE5qne/wCq6svAJxj8A3phd6sZGH1LmrO0WPQfZHAOiySXAluY/rvTnQJe0C2/GPjcBLN8m+5o7O3AA1X1e0Obhm97cRXwofOdbSlLZU6yF/gNYF9VfW1S+RYblbeq/qGqLqqqHVW1g0GBXl5V/z7BqN+yzM/F1O5/y2Seyv0vycxjnwxL8l3ASxi8r3A7g1vNQN99b9LvLJ/ju9LvBR4GvsFgR7iawQ/WnwL3AZ8BXjzpnD0yPx+4C7iHwTnD50w651De5zP41fBe4O7ucSWD21B/nMFO8XHg6ZPO2iPzPIPbaT829rZJZ10u76I5DzFdn7pZ6jWe2v1vmcxTuf8BzwI+2+W9D7i2G/8B4NPdz/P76D7htNzDWyBIUuNaPHUjSRpi0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TG/R+SwSkeynf+5QAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"for proba_var in [l1_post, l2_post]:\n",
" plt.hist(proba_var,normed=True)\n",
" plt.xlim([15, 30])\n",
" plt.ylim([0,1])\n",
" plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(30000,)"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"l1_post.shape"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFMlJREFUeJzt3X+s3Xd93/Hna07J6FqE01yiECezoaZbwlpTvBAJgbJmEAMVCe1YnU3EY1QmKJHabdpwNmlBoZHCCgOlhaAAboJGEzIoxQLT1I0Y2TRS7BSTH0Dmmx8lF3uJiVkBgVIlvPfH+dz24M+5P3zPvT7XzvMhfXW+5/39fM75fPI93Je/P84hVYUkScP+zqQHIElafQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdU6Z9ACW6vTTT6/169dPehiSdEK55557vlNVUwu1O2HDYf369ezbt2/Sw5CkE0qSv1xMO08rSZI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6J+w3pCVpktbv+PxE3vfR699wXN7HIwdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1FgyHJDuTPJHk/qHaJ5Psb8ujSfa3+vokPxra9uGhPi9Pcl+S6SQ3JEmrn5ZkT5ID7XHtSkxUkrR4izlyuBnYMlyoqt+oqk1VtQn4NPBHQ5sfmt1WVVcM1W8EtgMb2zL7mjuAO6tqI3Bney5JmqAFw6Gq7gKOjNrW/vX/z4Fb53uNJGcCz6uqL1dVAR8HLm2bLwFuaeu3DNUlSRMy7jWHVwGPV9WBodqGJF9N8qUkr2q1s4CZoTYzrQZwRlUdAmiPLxhzTJKkMY37w3uX8ZNHDYeAc6rqySQvB/44yXlARvStY32zJNsZnJrinHPOWcJwJUmLseQjhySnAL8GfHK2VlVPVdWTbf0e4CHgJQyOFNYNdV8HHGzrj7fTTrOnn56Y6z2r6qaq2lxVm6emppY6dEnSAsY5rfRPgW9W1d+cLkoylWRNW38RgwvPD7fTRd9PckG7TnE58NnWbRewra1vG6pLkiZkMbey3gp8GfiFJDNJ3tY2baW/EP1q4N4kXwM+BVxRVbMXs98BfBSYZnBE8YVWvx54TZIDwGvac0nSBC14zaGqLpuj/q9G1D7N4NbWUe33AS8dUX8SuGihcUiSjh+/IS1J6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6iwYDkl2Jnkiyf1DtXcl+XaS/W15/dC2q5NMJ3kwycVD9S2tNp1kx1B9Q5I/T3IgySeTPGc5JyhJOnaLOXK4Gdgyov7+qtrUlt0ASc4FtgLntT4fSrImyRrgg8DrgHOBy1pbgPe019oIfBd42zgTkiSNb8FwqKq7gCOLfL1LgNuq6qmqegSYBs5vy3RVPVxVfw3cBlySJMCvAJ9q/W8BLj3GOUiSltk41xyuSnJvO+20ttXOAh4bajPTanPVfw74f1X19FH1kZJsT7Ivyb7Dhw+PMXRJ0nyWGg43Ai8GNgGHgPe1eka0rSXUR6qqm6pqc1VtnpqaOrYRS5IW7ZSldKqqx2fXk3wE+Fx7OgOcPdR0HXCwrY+qfwd4fpJT2tHDcHtJ0oQs6cghyZlDT98EzN7JtAvYmuTUJBuAjcBXgL3AxnZn0nMYXLTeVVUFfBH4Z63/NuCzSxmTJGn5LHjkkORW4ELg9CQzwDXAhUk2MTgF9CjwdoCqeiDJ7cDXgaeBK6vqmfY6VwF3AGuAnVX1QHuLdwK3Jfkd4KvAx5ZtdpKkJVkwHKrqshHlOf+AV9V1wHUj6ruB3SPqDzO4m0mStEr4DWlJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1FgyHJDuTPJHk/qHa7yb5ZpJ7k3wmyfNbfX2SHyXZ35YPD/V5eZL7kkwnuSFJWv20JHuSHGiPa1diopKkxVvMkcPNwJajanuAl1bVLwL/B7h6aNtDVbWpLVcM1W8EtgMb2zL7mjuAO6tqI3Bney5JmqAFw6Gq7gKOHFX706p6uj29G1g332skORN4XlV9uaoK+Dhwadt8CXBLW79lqC5JmpDluObwr4EvDD3fkOSrSb6U5FWtdhYwM9RmptUAzqiqQwDt8QXLMCZJ0hhOGadzkv8EPA18opUOAedU1ZNJXg78cZLzgIzoXkt4v+0MTk1xzjnnLG3QkqQFLfnIIck24FeBf9lOFVFVT1XVk239HuAh4CUMjhSGTz2tAw629cfbaafZ009PzPWeVXVTVW2uqs1TU1NLHbokaQFLCockW4B3Am+sqh8O1aeSrGnrL2Jw4fnhdrro+0kuaHcpXQ58tnXbBWxr69uG6pKkCVnwtFKSW4ELgdOTzADXMLg76VRgT7sj9e52Z9KrgWuTPA08A1xRVbMXs9/B4M6n5zK4RjF7neJ64PYkbwO+Bbx5WWYmSVqyBcOhqi4bUf7YHG0/DXx6jm37gJeOqD8JXLTQOCRJx4/fkJYkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdRYVDkl2Jnkiyf1DtdOS7ElyoD2ubfUkuSHJdJJ7k/zyUJ9trf2BJNuG6i9Pcl/rc0OSLOckJUnH5pRFtrsZ+H3g40O1HcCdVXV9kh3t+TuB1wEb2/IK4EbgFUlOA64BNgMF3JNkV1V9t7XZDtwN7Aa2AF8Yb2rSs8/6HZ+fyPs+ev0bJvK+WjmLOnKoqruAI0eVLwFuaeu3AJcO1T9eA3cDz09yJnAxsKeqjrRA2ANsadueV1VfrqpiEECXIkmamHGuOZxRVYcA2uMLWv0s4LGhdjOtNl99ZkRdkjQhK3FBetT1glpCvX/hZHuSfUn2HT58eIwhSpLmM044PN5OCdEen2j1GeDsoXbrgIML1NeNqHeq6qaq2lxVm6empsYYuiRpPuOEwy5g9o6jbcBnh+qXt7uWLgD+qp12ugN4bZK17c6m1wJ3tG3fT3JBu0vp8qHXkiRNwKLuVkpyK3AhcHqSGQZ3HV0P3J7kbcC3gDe35ruB1wPTwA+BtwJU1ZEk7wb2tnbXVtXsRe53MLgj6rkM7lLyTiVJmqBFhUNVXTbHpotGtC3gyjleZyewc0R9H/DSxYxFkrTy/Ia0JKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKljOEiSOoaDJKmz5HBI8gtJ9g8t30vy20neleTbQ/XXD/W5Osl0kgeTXDxU39Jq00l2jDspSdJ4Tllqx6p6ENgEkGQN8G3gM8BbgfdX1XuH2yc5F9gKnAe8EPizJC9pmz8IvAaYAfYm2VVVX1/q2CRJ41lyOBzlIuChqvrLJHO1uQS4raqeAh5JMg2c37ZNV9XDAElua20NB0makOW65rAVuHXo+VVJ7k2yM8naVjsLeGyozUyrzVWXJE3I2OGQ5DnAG4H/3ko3Ai9mcMrpEPC+2aYjutc89VHvtT3JviT7Dh8+PNa4JUlzW44jh9cBf1FVjwNU1eNV9UxV/Rj4CH976mgGOHuo3zrg4Dz1TlXdVFWbq2rz1NTUMgxdkjTKcoTDZQydUkpy5tC2NwH3t/VdwNYkpybZAGwEvgLsBTYm2dCOQra2tpKkCRnrgnSSn2Zwl9Hbh8r/JckmBqeGHp3dVlUPJLmdwYXmp4Erq+qZ9jpXAXcAa4CdVfXAOOOSJI1nrHCoqh8CP3dU7S3ztL8OuG5EfTewe5yxSJKWj9+QliR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1xg6HJI8muS/J/iT7Wu20JHuSHGiPa1s9SW5IMp3k3iS/PPQ621r7A0m2jTsuSdLSLdeRwz+pqk1Vtbk93wHcWVUbgTvbc4DXARvbsh24EQZhAlwDvAI4H7hmNlAkScffSp1WugS4pa3fAlw6VP94DdwNPD/JmcDFwJ6qOlJV3wX2AFtWaGySpAUsRzgU8KdJ7kmyvdXOqKpDAO3xBa1+FvDYUN+ZVpur/hOSbE+yL8m+w4cPL8PQJUmjnLIMr/HKqjqY5AXAniTfnKdtRtRqnvpPFqpuAm4C2Lx5c7ddkrQ8xj5yqKqD7fEJ4DMMrhk83k4X0R6faM1ngLOHuq8DDs5TlyRNwFjhkOTvJfnZ2XXgtcD9wC5g9o6jbcBn2/ou4PJ219IFwF+10053AK9NsrZdiH5tq0mSJmDc00pnAJ9JMvtaf1hVf5JkL3B7krcB3wLe3NrvBl4PTAM/BN4KUFVHkrwb2NvaXVtVR8YcmyRpicYKh6p6GPilEfUngYtG1Au4co7X2gnsHGc8kqTl4TekJUkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1FlyOCQ5O8kXk3wjyQNJfqvV35Xk20n2t+X1Q32uTjKd5MEkFw/Vt7TadJId401JkjSucf4/pJ8G/l1V/UWSnwXuSbKnbXt/Vb13uHGSc4GtwHnAC4E/S/KStvmDwGuAGWBvkl1V9fUxxiZJGsOSw6GqDgGH2vr3k3wDOGueLpcAt1XVU8AjSaaB89u26ap6GCDJba2t4SBJE7Is1xySrAdeBvx5K12V5N4kO5OsbbWzgMeGus202lx1SdKEjB0OSX4G+DTw21X1PeBG4MXAJgZHFu+bbTqie81TH/Ve25PsS7Lv8OHD4w5dkjSHscIhyU8xCIZPVNUfAVTV41X1TFX9GPgIf3vqaAY4e6j7OuDgPPVOVd1UVZuravPU1NQ4Q5ckzWOcu5UCfAz4RlX916H6mUPN3gTc39Z3AVuTnJpkA7AR+AqwF9iYZEOS5zC4aL1rqeOSJI1vnLuVXgm8Bbgvyf5W+4/AZUk2MTg19CjwdoCqeiDJ7QwuND8NXFlVzwAkuQq4A1gD7KyqB8YYlyRpTOPcrfS/GH29YPc8fa4DrhtR3z1fP0nS8eU3pCVJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQZ51dZpQWt3/H5ib33o9e/YWLvLZ3oPHKQJHUMB0lSx3CQJHUMB0lSx3CQJHUMB0lSZ9WEQ5ItSR5MMp1kx6THI0nPZqsiHJKsAT4IvA44F7gsybmTHZUkPXutinAAzgemq+rhqvpr4DbgkgmPSZKetVbLN6TPAh4bej4DvGKl3mxS39r1G7uSThSrJRwyolZdo2Q7sL09/UGSB5f4fqcD31li3yXLe1bkZScylxWw7PNYof/ei3Gy7BNY5Fwm+N/6WJwU+yXvGXsef38xjVZLOMwAZw89XwccPLpRVd0E3DTumyXZV1Wbx32d1eBkmcvJMg9wLqvVyTKX4zWP1XLNYS+wMcmGJM8BtgK7JjwmSXrWWhVHDlX1dJKrgDuANcDOqnpgwsOSpGetVREOAFW1G9h9nN5u7FNTq8jJMpeTZR7gXFark2Uux2Ueqequ+0qSnuVWyzUHSdIqctKFQ5I1Sb6a5HPt+Sfaz3Lcn2Rnkp+ao9+2JAfasu34jnq0MebyTJL9bVkVF/ZHzOVjSb6W5N4kn0ryM3P0u7r9pMqDSS4+vqMebSlzSbI+yY+G9suHj//IuzH9xDyG6r+X5Afz9Fv1+2SoPudcVuM+gZGfr5uTPDI0zk1z9FvWv2EnXTgAvwV8Y+j5J4B/APwj4LnAbx7dIclpwDUMvnh3PnBNkrUrP9QFHfNcmh9V1aa2vHGFx7hYR8/l31TVL1XVLwLfAq46ukP7CZWtwHnAFuBD7adWJu2Y59I8NLRfrljxUS7s6HmQZDPw/Lk6nED7ZMG5NKttn8CIuQD/fmic+4/usBJ/w06qcEiyDngD8NHZWlXtrgb4CoPvUBztYmBPVR2pqu8Cexh88CdmjLmsOnPM5XttWxgE3aiLX5cAt1XVU1X1CDDN4IM/MWPMZVUZNY/2R/53gf8wT9cTYp8sci6rzqi5LNKy/w07qcIB+ACDD8OPj97QTsG8BfiTEf1G/XzHWSsxwGOw1LkA/N0k+5LcneTSFRzjYo2cS5I/AP4vg6Oh3xvR74TZL4uYC8CGdrrgS0letbLDXNCoeVwF7KqqQ/P0O1H2yWLmAqtrn8Dc/7u/rp22fH+SU0f0W/b9ctKEQ5JfBZ6oqnvmaPIh4K6q+p+juo+oTexff2POBeCc9g3KfwF8IMmLV2KcizHfXKrqrcALGRxC/8ao7iNqq3K/LGIuhxjsl5cB/xb4wyTPW8nxzmXUPJK8EHgzcwfb3zQdUVtV++QY5rJq9gnM+/m6msE/Ov4xcBrwzlHdR9TG2i8nTTgArwTemORRBr/q+itJ/htAkmuAKQYfgFEW9fMdx9E4c6GqDrbHh4H/Abxshcc7nznnAlBVzwCfBH59RN8TZr/A/HNpp2GebOv3AA8BLzkegx6hmwfwAPDzwHSr/3SS6RF9V/0+YZFzWWX7BOb4fFXVoXY2+SngDxh9Gm/590tVnXQLcCHwubb+m8D/Bp47T/vTgEeAtW15BDht0vNY4lzWAqe29dOBA8C5k57H8FwY/Cvn51stwHuB945ofx7wNeBUYAPwMLBm0vNY4lymZscOvAj49mr4jA1/vo6q/2CO9qt+nxzDXFblPjl6LsCZQ5+vDwDXj2i/7H/DTqYjh7l8GDgD+HK7Dew/w+BOhiQfBaiqI8C7GfzG017g2lZbbRacC/APgX1JvgZ8kcEH6euTGe6cAtyS5D7gPuBM4FqAJG9Mci1ADX5C5Xbg6wyur1xZg3+dryaLmgvwauDetl8+BVyxSj9jnRNwn8zpBN0nnxj6fJ0O/A6s/N8wvyEtSeo8G44cJEnHyHCQJHUMB0lSx3CQJHUMB0lSx3CQJHUMB0lSx3CQJHX+P64FguEWQPBcAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(tau_post)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"30000\n"
]
}
],
"source": [
"N = tau_post.shape[0]\n",
"print(N)\n",
"expec = [0] * n_count_data\n",
"for day in range(n_count_data):\n",
" ix = day < tau_post\n",
" expec[day] = (l1_post[ix].sum()+l2_post[~ix].sum()) / N\n",
"\n",
"# l1_post[ix].sum()"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<BarContainer object of 74 artists>"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEtJJREFUeJzt3XuwnHV9x/H3l3BTtIaQQxqJ6cGZDOIfJeiZiEPLIFEq4gBjgYGhTtpJJzOOdWCqo8FiFcc/YrWCtB1rBpA4Y7mIYhC8MQdoa/9Aw00ugUZohEhMIhLjZUZu3/6xT8JycsI+u2c3++zvvF8zZ3afy57nk93NJ7/97T6byEwkSaPvgGEHkCT1h4UuSYWw0CWpEBa6JBXCQpekQljoklQIC12SCmGhS1IhLHRJKsSB+/Ng8+fPz/Hx8f15SEkaeXffffcvM3Os0377tdDHx8fZsGHD/jykJI28iPhZnf2ccpGkQljoklQIC12SCmGhS1IhLHRJKoSFLkmFsNAlqRAWuiQVwkKXpELs1zNFNZrGV9+65/rmNacPMYmkV+IIXZIK0bHQI+KYiLiv7WdXRFwUEfMi4raI2FRdHr4/AkuSptex0DPz0cxcmplLgbcCvwduAlYDk5m5BJisliVJQ9LtlMty4LHM/BlwJrCuWr8OOKufwSRJ3en2TdHzgGur6wsycytAZm6NiCOnu0FErAJWASxevLjXnNKs4BvQmonaI/SIOBg4A/h6NwfIzLWZOZGZE2NjHb+fXZLUo26mXE4D7snMbdXytohYCFBdbu93OElSfd0U+vm8NN0CcDOworq+Aljfr1CSpO7VKvSIeDXwLuCbbavXAO+KiE3VtjX9jydJqqvWm6KZ+XvgiCnrnqb1qRdJUgN4pqgkFcJCl6RCWOiSVAgLXZIKYaFLUiEsdEkqhIUuSYWw0CWpEBa6JBXCQpekQljoklQIC12SCmGhS1IhLHRJKoSFLkmFsNAlqRAWuiQVwkKXpEJY6JJUCAtdkgpRq9AjYm5E3BgRj0TExoh4e0TMi4jbImJTdXn4oMNKkvat7gj9i8D3MvNNwHHARmA1MJmZS4DJalmSNCQdCz0i/gg4CbgKIDOfzcydwJnAumq3dcBZgwopSeqszgj9jcAO4CsRcW9EXBkRhwELMnMrQHV55HQ3johVEbEhIjbs2LGjb8ElSS9Xp9APBN4CfCkzjwd+RxfTK5m5NjMnMnNibGysx5iSpE7qFPoWYEtm3lUt30ir4LdFxEKA6nL7YCJKkuroWOiZ+QvgyYg4plq1HHgYuBlYUa1bAawfSEJJUi0H1tzvQ8DXIuJg4HHgb2j9Y3BDRKwEngDOGUxESVIdtQo9M+8DJqbZtLy/cSRJvfJMUUkqhIUuSYWw0CWpEHXfFNWAjK++dc/1zWtOH2ISSaPOEbokFcJCl6RCWOiSVAgLXZIKYaFLUiEsdEkqhIUuSYWw0CWpEBa6JBXCQpekQnjqvzTC/OoItXOELkmFsNAlqRAWuiQVwkKXpELUelM0IjYDvwFeAJ7PzImImAdcD4wDm4FzM/OZwcSUJHXSzQj9HZm5NDN3/2fRq4HJzFwCTFbLkqQhmcmUy5nAuur6OuCsmceRJPWqbqEn8IOIuDsiVlXrFmTmVoDq8shBBJQk1VP3xKITM/OpiDgSuC0iHql7gOofgFUAixcv7iGiJKmOWiP0zHyqutwO3AQsA7ZFxEKA6nL7Pm67NjMnMnNibGysP6klSXvpWOgRcVhEvHb3deBU4EHgZmBFtdsKYP2gQkqSOqsz5bIAuCkidu//H5n5vYj4MXBDRKwEngDOGVxMaXbyu1rUjY6FnpmPA8dNs/5pYPkgQkmSuueZopJUCAtdkgphoUtSISx0SSqEhS5JhbDQJakQFrokFcJCl6RCWOiSVAgLXZIKYaFLUiEsdEkqhIUuSYWw0CWpEBa6JBXCQpekQljoklQIC12SCmGhS1IhLHRJKkTtQo+IORFxb0TcUi0fHRF3RcSmiLg+Ig4eXExJUifdjNAvBDa2LX8WuCwzlwDPACv7GUyS1J1ahR4Ri4DTgSur5QBOAW6sdlkHnDWIgJKkeuqO0C8HPgq8WC0fAezMzOer5S3AUX3OJknqQsdCj4j3Atsz8+721dPsmvu4/aqI2BARG3bs2NFjTElSJ3VG6CcCZ0TEZuA6WlMtlwNzI+LAap9FwFPT3Tgz12bmRGZOjI2N9SGyJGk6HQs9My/OzEWZOQ6cB9yemRcAdwBnV7utANYPLKUkqaMDO++yTx8DrouIzwD3Alf1J5Ik7V/jq2992fLmNacPKcnMdFXomXkncGd1/XFgWf8jSZJ64ZmiklQIC12SCmGhS1IhZvKmqCQVq/2N0lF5k9QRuiQVwkKXpEJY6JJUCAtdkgphoUtSIWb1p1xG8V1sqen8ezU8jtAlqRAWuiQVwkKXpEJY6JJUCAtdkgphoUtSISx0SSqEhS5JhZjVJxapN544IjWTI3RJKkTHQo+IQyPiRxFxf0Q8FBGXVuuPjoi7ImJTRFwfEQcPPq4kaV/qTLn8ATglM38bEQcBP4yI7wJ/D1yWmddFxL8DK4EvDTCrpBlqny6D0Z0yc9pveh1H6Nny22rxoOongVOAG6v164CzBpJQklRLrTn0iJgTEfcB24HbgMeAnZn5fLXLFuCowUSUJNVR61MumfkCsDQi5gI3AcdOt9t0t42IVcAqgMWLF/cYU5L2zSmYlq4+5ZKZO4E7gROAuRGx+x+ERcBT+7jN2sycyMyJsbGxmWSVJL2CjiP0iBgDnsvMnRHxKuCdwGeBO4CzgeuAFcD6QQZViyMRSftSZ8plIbAuIubQGtHfkJm3RMTDwHUR8RngXuCqAeaUJHXQsdAz8yfA8dOsfxxYNohQkqTueaaoJBXCQpekQljoklQIv21RjeSnedRPpXzlQSeO0CWpEBa6JBXCQpekQljoklQIC12SCmGhS1IhLHRJKoSFLkmFGNkTizqdeDJbTiSQpN0coUtSIUZ2hC6V6Ijf7eT0R/6bAzLhisf46w0PvbTxisf22r/T9lfcv+ZtutVtpn4co6v7aZp9ptvetz/HsmVwwgm9374LFroawe9uafn4nVfzlw/e3lqYhE+1b5zce/9O219x/5q36dbLjjGA3z/dMTod81NTV0x23t7pd9b2iU9Y6NKss2ULZzz8n3z1+NP55z//K+7/5Kkcd+kP9my+/5On7nWTTttfaf+6t+lWN5l7zTD1GDM95nTbu71v9+nQQ3u/bZcsdKkprriCAzJZ+7b38etXvRbmzWtd7jZv3l436bT9FfeveZtudZW5xwxTjzHTY063vdv7tgks9D6a+smaqZoyleD0RgPt2gVf/jLfOeZEtrxuQc+/xse2OYbxWPgpF6kJrrwSdu1i7bL3DTuJRljHEXpEvAH4KvDHwIvA2sz8YkTMA64HxoHNwLmZ+czgokqFeu45uPxyOPlkHli4pK+/utOrRpWlzgj9eeDDmXkscALwwYh4M7AamMzMJbTeA149uJhSwW64AZ58Ej7ykWEn0YjrWOiZuTUz76mu/wbYCBwFnAmsq3ZbB5w1qJBSsTLhc5+DY4+F004bdhqNuK7eFI2IceB44C5gQWZuhVbpR8SR+7jNKmAVwOLFi2eSVbNYsW/2TU7C/ffDVVfBAaPxllaxj0UBaj+DIuI1wDeAizJzV93bZebazJzIzImxsbFeMkplevpp+MAH4PWvhwsuGHYaFaDWCD0iDqJV5l/LzG9Wq7dFxMJqdL4Q2D6okFJxnn0Wzj4bnngCbr8dDjlk2IlUgDqfcgngKmBjZn6hbdPNwApgTXW5fiAJZ7luX96OyrdMzuqX7Zmtkfmdd3LRez/Mt769E75961Duh1F5vqieOiP0E4H3Aw9ExH3Vuo/TKvIbImIl8ARwzmAiSoX5/Ofh6qvhkkv41nP75zs+NDt0LPTM/CEQ+9i8vL9xpAJkwosvwgsvtH5+9SvYsgV+/nN44AG49FI499zW5ce/O+y0KsjonPr/61/D4YfvWXw827b90967v2x7nX2m2d6tvY45VQ8Zpm7fH3/uF9u2HxDdZ6ij0+8cxDH3i+z0JABOOgmuucZPtajvRqfQDzkELrlkz+K/TG7ac/3C5XufXde+vc4+023v1tRjTtVLhqnb98efu9Mx+3G/dXuMfj9WA3XAATBnzkuXc+fCokWtn6OOgvnzIfb1olfq3egU+qGHwqc/vWfxsmdfGjVc+Om9Rw3t2+vsM932bk095lS9ZJi6fX/8uTsdsx/3W7fH6PdjJZVodApdegVOCzRXp++T8ZM2/TMak3iSpI4sdEkqhFMu0j44jaNR4whdkgphoUtSIZxy6YIvwWe3mX6vzlQ+h9RvjtAlqRAWuiQVwimXNk6pzC4+3pqJJj5/HKFLUiFGZoTexH8NVc8wTu32dPLZbRh90YSOcoQuSYWw0CWpECMz5aLe1Jl6aMJLRY0unz/N4QhdkgrRsdAj4uqI2B4RD7atmxcRt0XEpury8Ff6HZKkwasz5XIN8K/AV9vWrQYmM3NNRKyulj/W/3jD1emlZC8vNQfxO6V+8jk4ujqO0DPzv4BfTVl9JrCuur4OOKvPuSRJXep1Dn1BZm4FqC6P7F8kSVIvBv4pl4hYBawCWLx48aAPJ03LE400G/Q6Qt8WEQsBqsvt+9oxM9dm5kRmToyNjfV4OElSJ70W+s3Aiur6CmB9f+JIknrVccolIq4FTgbmR8QW4JPAGuCGiFgJPAGcM8iQajY/FSE1Q8dCz8zz97FpeZ+zSJJmwDNFJakQxXyXS6f/v1GSSucIXZIKYaFLUiGKmXKpw09j7B+zZfprtv6vOGouR+iSVAgLXZIKMaumXEo0W1+CD+LPPVvvS5XDEbokFcIRuiTVMAqv4ByhS1IhLHRJKoRTLpIaZxSmN5rIEbokFcJCl6RCWOiSVAgLXZIKYaFLUiH8lIvUJ34yY3C8b+txhC5JhZhRoUfEuyPi0Yj4aUSs7lcoSVL3ei70iJgD/BtwGvBm4PyIeHO/gkmSujOTEfoy4KeZ+XhmPgtcB5zZn1iSpG7NpNCPAp5sW95SrZMkDUFkZm83jDgH+IvM/Ntq+f3Assz80JT9VgGrqsVjgEd7jwvAfOCXM/wdg2bG/hiFjDAaOc3YH8PK+CeZOdZpp5l8bHEL8Ia25UXAU1N3ysy1wNoZHOdlImJDZk706/cNghn7YxQywmjkNGN/ND3jTKZcfgwsiYijI+Jg4Dzg5v7EkiR1q+cRemY+HxF/B3wfmANcnZkP9S2ZJKkrMzpTNDO/A3ynT1nq6tv0zQCZsT9GISOMRk4z9kejM/b8pqgkqVk89V+SCjEyhd7UrxmIiKsjYntEPNi2bl5E3BYRm6rLw4ec8Q0RcUdEbIyIhyLiwqbljIhDI+JHEXF/lfHSav3REXFXlfH66g34oYqIORFxb0Tc0sSMEbE5Ih6IiPsiYkO1rjGPdVvOuRFxY0Q8Uj03396knBFxTHUf7v7ZFREXNSnjVCNR6A3/moFrgHdPWbcamMzMJcBktTxMzwMfzsxjgROAD1b3X5Ny/gE4JTOPA5YC746IE4DPApdVGZ8BVg4x424XAhvblpuY8R2ZubTtI3ZNeqx3+yLwvcx8E3Acrfu0MTkz89HqPlwKvBX4PXBTkzLuJTMb/wO8Hfh+2/LFwMXDztWWZxx4sG35UWBhdX0h8OiwM07Jux54V1NzAq8G7gHeRuskjgOnex4MKdsiWn+JTwFuAaKBGTcD86esa9RjDfwR8H9U7+M1NWdbrlOB/2lyxswcjRE6o/c1AwsycytAdXnkkPPsERHjwPHAXTQsZzWVcR+wHbgNeAzYmZnPV7s04XG/HPgo8GK1fATNy5jADyLi7upMbWjYYw28EdgBfKWavroyIg6jeTl3Ow+4trre1IwjU+gxzTo/ntOliHgN8A3goszcNew8U2XmC9l6ebuI1pe/HTvdbvs31Usi4r3A9sy8u331NLsO+7l5Yma+hdYU5Qcj4qQh55nOgcBbgC9l5vHA72jS1EWb6j2RM4CvDztLJ6NS6LW+ZqBBtkXEQoDqcvuQ8xARB9Eq869l5jer1Y3LCZCZO4E7ac33z42I3edLDPtxPxE4IyI20/p20VNojdiblJHMfKq63E5rzncZzXustwBbMvOuavlGWgXftJzQ+ofxnszcVi03MSMwOoU+al8zcDOworq+gtac9dBERABXARsz8wttmxqTMyLGImJudf1VwDtpvUl2B3B2tdtQM2bmxZm5KDPHaT0Hb8/MC2hQxog4LCJeu/s6rbnfB2nQYw2Qmb8AnoyIY6pVy4GHaVjOyvm8NN0CzczYMuxJ/C7elHgP8L+05lX/Ydh52nJdC2wFnqM16lhJa151EthUXc4bcsY/ozUN8BPgvurnPU3KCfwpcG+V8UHgH6v1bwR+BPyU1kveQ4b9mFe5TgZuaVrGKsv91c9Du/+uNOmxbsu6FNhQPebfAg5vWk5ab9A/DbyubV2jMrb/eKaoJBViVKZcJEkdWOiSVAgLXZIKYaFLUiEsdEkqhIUuSYWw0CWpEBa6JBXi/wFB22lbEtvEGwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# plt.hist(np.arange(n_count_data),)\n",
"plt.plot(np.arange(n_count_data),expec,'-r')\n",
"plt.bar(np.arange(n_count_data),count_data)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([44, 45, 45, ..., 45, 45, 45])"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tau_post"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment