ABC python - Approximate Bayesian Computation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resources\n",
    "\n",
    "https://www.youtube.com/watch?v=nKCT-Cdk0xY\n",
    "\n",
    "https://academic.oup.com/bioinformatics/article/36/10/3286/5727756"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generative model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_pairs = 9\n",
    "n_odds = 5\n",
    "n_socks_types = n_pairs + n_odds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "paired_socks = sorted(list(range(1, n_pairs+1))*2)\n",
    "odds_socks = list(range(n_pairs+1, n_pairs+n_odds+1))\n",
    "socks = paired_socks + odds_socks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[8, 2, 1]"
      ]
     },
     "execution_count": 156,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import random\n",
    "random.sample(socks, 3) # sample without replacement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[5, 4, 1, 2, 9, 10, 5, 6, 6, 4, 8]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def pick_socks(n_pairs=1, n_odd=1, n_pick=1):\n",
    "    paired_socks = sorted(list(range(1, n_pairs+1))*2)\n",
    "    odds_socks = list(range(n_pairs+1, n_pairs+n_odds+1))\n",
    "    socks = paired_socks + odds_socks\n",
    "    assert n_pick <= n_odd + 2*n_pairs\n",
    "    return random.sample(socks, n_pick)\n",
    "\n",
    "pick_socks(n_pairs=9, n_odd=5, n_pick=11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prior distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import nbinom, beta\n",
    "import pylab as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  8.,  29.,  30., 103., 199., 190., 313., 298., 202., 227., 178.,\n",
       "         87.,  57.,  42.,  11.,  11.,  10.,   2.,   1.,   2.]),\n",
       " array([10.  , 12.65, 15.3 , 17.95, 20.6 , 23.25, 25.9 , 28.55, 31.2 ,\n",
       "        33.85, 36.5 , 39.15, 41.8 , 44.45, 47.1 , 49.75, 52.4 , 55.05,\n",
       "        57.7 , 60.35, 63.  ]),\n",
       " <BarContainer object of 20 artists>)"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPxUlEQVR4nO3df6zddX3H8edrFH9MjQV713Rtt8u0m8FkFHLDMJoFISqiGZo4gtm0MST1D0wwcZnVf9RlJJhMmSYbSRUmLio2KqNRorKKcf4heKsILZVYsYQ2hV4V/DEzkuJ7f5wPeqi3vT/Ovb33fu7zkZyc7/fz/XzPeX/S09f95nO+3+9JVSFJ6ssfLHUBkqSFZ7hLUocMd0nqkOEuSR0y3CWpQ2uWugCAdevW1fj4+FKXIUkryt69e39SVWPTbVsW4T4+Ps7k5ORSlyFJK0qSh0+2zWkZSeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nq0LK4QlUrx/iOL89730PXv34BK5F0Kh65S1KHDHdJ6pDhLkkdMtwlqUMzhnuS5yS5J8n3k+xP8sHWfk6Su5McTPK5JM9q7c9u6wfb9vFFHoMk6QSzOXJ/Erikqs4DtgKXJbkI+BBwQ1W9BHgcuLr1vxp4vLXf0PpJkk6jGcO9Bn7VVs9sjwIuAT7f2m8B3tiWr2jrtO2XJslCFSxJmtms5tyTnJHkXuAYcCfwI+CJqjreuhwGNrbljcAjAG37z4EXTfOa25NMJpmcmpoaaRCSpGeaVbhX1VNVtRXYBFwIvHTUN66qnVU1UVUTY2PT/gSgJGme5nSFalU9keQu4OXA2iRr2tH5JuBI63YE2AwcTrIGeCHw0wWsWSuUV7dKp89szpYZS7K2LT8XeDVwALgLeHPrtg24vS3vbuu07V+vqlrAmiVJM5jNkfsG4JYkZzD4Y7Crqr6U5AHg1iT/DHwPuKn1vwn4zyQHgZ8BVy1C3ZKkU5gx3KvqPuD8adofYjD/fmL7/wF/uyDVSZLmxStUJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHZox3JNsTnJXkgeS7E9ybWv/QJIjSe5tj8uH9nlvkoNJHkzy2sUcgCTp962ZRZ/jwLur6rtJXgDsTXJn23ZDVf3LcOck5wJXAS8D/hj47yR/XlVPLWThkqSTm/HIvaqOVtV32/IvgQPAxlPscgVwa1U9WVU/Bg4CFy5EsZKk2ZnTnHuSceB84O7W9M4k9yW5OclZrW0j8MjQboeZ5o9Bku1JJpNMTk1Nzb1ySdJJzTrckzwf+ALwrqr6BXAj8GJgK3AU+PBc3riqdlbVRFVNjI2NzWVXSdIMZhXuSc5kEOyfrqovAlTVY1X1VFX9Bvg4v5t6OQJsHtp9U2uTJJ0mszlbJsBNwIGq+shQ+4ahbm8C9rXl3cBVSZ6d5BxgC3DPwpUsSZrJbM6WeQXwVuD+JPe2tvcBb0myFSjgEPAOgKran2QX8ACDM22u8UwZSTq9Zgz3qvoWkGk23XGKfa4DrhuhLknSCLxCVZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHZnNvGWlFG9/x5ZH2P3T96xeoEun08chdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQzOGe5LNSe5K8kCS/Umube1nJ7kzyQ/b81mtPUk+luRgkvuSXLDYg5AkPdNsjtyPA++uqnOBi4BrkpwL7AD2VNUWYE9bB3gdsKU9tgM3LnjVkqRTmvF+7lV1FDjaln+Z5ACwEbgCuLh1uwX4BvCe1v6pqirg20nWJtnQXkeal1HvyS6tNnOac08yDpwP3A2sHwrsR4H1bXkj8MjQbodb24mvtT3JZJLJqampudYtSTqFWf8SU5LnA18A3lVVv0jy221VVUlqLm9cVTuBnQATExNz2ne185eFJM1kVkfuSc5kEOyfrqovtubHkmxo2zcAx1r7EWDz0O6bWpsk6TSZzdkyAW4CDlTVR4Y27Qa2teVtwO1D7W9rZ81cBPzc+XZJOr1mMy3zCuCtwP1J7m1t7wOuB3YluRp4GLiybbsDuBw4CPwaePtCFixJmtlszpb5FpCTbL50mv4FXDNiXZKkEXiFqiR1aNZny6gfnjMu9c8jd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOuR57tIMRrkuwDtwaql45C5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOjRjuCe5OcmxJPuG2j6Q5EiSe9vj8qFt701yMMmDSV67WIVLkk5uNkfunwQum6b9hqra2h53ACQ5F7gKeFnb59+TnLFQxUqSZmfGcK+qbwI/m+XrXQHcWlVPVtWPgYPAhSPUJ0mah1Hm3N+Z5L42bXNWa9sIPDLU53BrkySdRvMN9xuBFwNbgaPAh+f6Akm2J5lMMjk1NTXPMiRJ05lXuFfVY1X1VFX9Bvg4v5t6OQJsHuq6qbVN9xo7q2qiqibGxsbmU4Yk6STmFe5JNgytvgl4+kya3cBVSZ6d5BxgC3DPaCVKkuZqzUwdknwWuBhYl+Qw8H7g4iRbgQIOAe8AqKr9SXYBDwDHgWuq6qlFqVySdFIzhntVvWWa5ptO0f864LpRipIkjcYrVCWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR2aMdyT3JzkWJJ9Q21nJ7kzyQ/b81mtPUk+luRgkvuSXLCYxUuSpjebI/dPAped0LYD2FNVW4A9bR3gdcCW9tgO3LgwZUqS5mLNTB2q6ptJxk9ovgK4uC3fAnwDeE9r/1RVFfDtJGuTbKiqowtWcSfGd3x5qUuQ1LH5zrmvHwrsR4H1bXkj8MhQv8OtTZJ0Go38hWo7Sq+57pdke5LJJJNTU1OjliFJGjLfcH8syQaA9nystR8BNg/129Tafk9V7ayqiaqaGBsbm2cZkqTpzDfcdwPb2vI24Pah9re1s2YuAn7ufLsknX4zfqGa5LMMvjxdl+Qw8H7gemBXkquBh4ErW/c7gMuBg8CvgbcvQs2SpBnM5myZt5xk06XT9C3gmlGLkiSNxitUJalDhrskdchwl6QOzTjnLmn+RrkS+dD1r1/ASrTaeOQuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdcjfUJWWKX9/VaPwyF2SOjTSkXuSQ8AvgaeA41U1keRs4HPAOHAIuLKqHh+tTEnSXCzEkfurqmprVU209R3AnqraAuxp65Kk02gxpmWuAG5py7cAb1yE95AkncKo4V7A15LsTbK9ta2vqqNt+VFg/XQ7JtmeZDLJ5NTU1IhlSJKGjXq2zCur6kiSPwLuTPKD4Y1VVUlquh2raiewE2BiYmLaPpKk+RnpyL2qjrTnY8BtwIXAY0k2ALTnY6MWKUmam3mHe5LnJXnB08vAa4B9wG5gW+u2Dbh91CIlSXMzyrTMeuC2JE+/zmeq6itJvgPsSnI18DBw5ehlSpLmYt7hXlUPAedN0/5T4NJRipI0Gq9ulVeoSlKHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIX+JaQSjnEssSYtp1Ye7AS2pR07LSFKHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHVr195aR9Ewr9X5L/rD3M3nkLkkdWvFH7iv1KEOSFpNH7pLUoUUL9ySXJXkwycEkOxbrfSRJv29RpmWSnAH8G/Bq4DDwnSS7q+qBxXg/SRplinbUL2OX8r1PZrHm3C8EDlbVQwBJbgWuAAx3SctOj9/dLVa4bwQeGVo/DPzVcIck24HtbfVXSR6c53utA34yz31XktUwztUwRlgd41wNY4QFGGc+NNL7/+nJNizZ2TJVtRPYOerrJJmsqokFKGlZWw3jXA1jhNUxztUwRlje41ysL1SPAJuH1je1NknSabBY4f4dYEuSc5I8C7gK2L1I7yVJOsGiTMtU1fEk7wS+CpwB3FxV+xfjvViAqZ0VYjWMczWMEVbHOFfDGGEZjzNVtdQ1SJIWmFeoSlKHDHdJ6tCKCvckNyc5lmTfUNvZSe5M8sP2fNZS1jiqJJuT3JXkgST7k1zb2nsb53OS3JPk+22cH2zt5yS5u9224nPtC/kVLckZSb6X5EttvccxHkpyf5J7k0y2tt4+s2uTfD7JD5IcSPLy5TzGFRXuwCeBy05o2wHsqaotwJ62vpIdB95dVecCFwHXJDmX/sb5JHBJVZ0HbAUuS3IR8CHghqp6CfA4cPXSlbhgrgUODK33OEaAV1XV1qHzvnv7zH4U+EpVvRQ4j8G/6fIdY1WtqAcwDuwbWn8Q2NCWNwAPLnWNCzze2xnco6fbcQJ/CHyXwVXMPwHWtPaXA19d6vpGHNsmBv/pLwG+BKS3MbZxHALWndDWzWcWeCHwY9pJKCthjCvtyH0666vqaFt+FFi/lMUspCTjwPnA3XQ4zjZdcS9wDLgT+BHwRFUdb10OM7iVxUr2r8A/Ar9p6y+ivzECFPC1JHvbrUWgr8/sOcAU8B9tiu0TSZ7HMh5jD+H+WzX489nFuZ1Jng98AXhXVf1ieFsv46yqp6pqK4Oj2wuBly5tRQsryRuAY1W1d6lrOQ1eWVUXAK9jMJX418MbO/jMrgEuAG6sqvOB/+WEKZjlNsYewv2xJBsA2vOxJa5nZEnOZBDsn66qL7bm7sb5tKp6AriLwRTF2iRPX1y30m9b8Qrgb5IcAm5lMDXzUfoaIwBVdaQ9HwNuY/DHuqfP7GHgcFXd3dY/zyDsl+0Yewj33cC2tryNwRz1ipUkwE3Agar6yNCm3sY5lmRtW34ug+8VDjAI+Te3bit6nFX13qraVFXjDG7B8fWq+js6GiNAkuclecHTy8BrgH109JmtqkeBR5L8RWu6lMEtzJftGFfUFapJPgtczOA2m48B7wf+C9gF/AnwMHBlVf1siUocWZJXAv8D3M/v5mnfx2Devadx/iVwC4PbU/wBsKuq/inJnzE4yj0b+B7w91X15NJVujCSXAz8Q1W9obcxtvHc1lbXAJ+pquuSvIi+PrNbgU8AzwIeAt5O++yyDMe4osJdkjQ7PUzLSJJOYLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDv0/O9c+VfuWU+YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_socks_prior = {'r': 30, 'p': 0.5}\n",
    "r = nbinom.rvs(n_socks_prior['r'], n_socks_prior['p'], size=2000)  # get random values\n",
    "\n",
    "plt.hist(r, bins=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  3.,   1.,   6.,   4.,   6.,  10.,   6.,  17.,  19.,  39.,  43.,\n",
       "         50.,  60.,  78., 104., 105., 123., 150., 115.,  61.]),\n",
       " array([0.55348782, 0.57566924, 0.59785066, 0.62003208, 0.64221349,\n",
       "        0.66439491, 0.68657633, 0.70875774, 0.73093916, 0.75312058,\n",
       "        0.77530199, 0.79748341, 0.81966483, 0.84184625, 0.86402766,\n",
       "        0.88620908, 0.9083905 , 0.93057191, 0.95275333, 0.97493475,\n",
       "        0.99711616]),\n",
       " <BarContainer object of 20 artists>)"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPj0lEQVR4nO3df4xlZX3H8fenbJFaq4vuSOkuONu6aqmpkUwojYkS6Y8VjEtbQyBtXe3GjZVaW211sX/QNDHBtCnVxJqsQkFjQUpt2RSsJQghbYS6iCA/RLfIj10XdxShP2yq2G//uMc6zM7u3Lnnzr27D+9XMplznnPOPd882fnsM8+955lUFZKktvzQtAuQJI2f4S5JDTLcJalBhrskNchwl6QGrZl2AQDr1q2r2dnZaZchSUeV22+//RtVNbPUsSMi3GdnZ9m9e/e0y5Cko0qShw51zGkZSWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lq0BHxhKok9TW747qRr33w4rPHWMmRwZG7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1aNlwT3JZkgNJ7l7i2DuTVJJ13X6SfCDJniR3JTl1NYqWJB3eMCP3y4HNixuTnAT8EvDwgubXAJu6r+3Ah/qXKElaqWXDvapuAR5b4tAlwLuAWtC2BfhoDdwKrE1y4lgqlSQNbaQ59yRbgH1VdeeiQ+uBRxbs7+3aJEkTtOJVIZM8E3gPgymZkSXZzmDqhpNPPrnPS0mSFhll5P5TwEbgziQPAhuAzyf5cWAfcNKCczd0bQepqp1VNVdVczMzMyOUIUk6lBWHe1V9saqeX1WzVTXLYOrl1Kp6FNgFvKH71MzpwBNVtX+8JUuSljPMRyGvBD4LvDjJ3iTbDnP69cADwB7gw8Bbx1KlJGlFlp1zr6rzlzk+u2C7gAv6lyVJ6sMnVCWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNWjZP7MnSZMyu+O6aZfQjGH+QPZlSQ4kuXtB258m+VKSu5L8XZK1C45dmGRPkvuT/PIq1S1JOoxhpmUuBzYvarsBeGlV/SzwZeBCgCSnAOcBP9Nd85dJjhlbtZKkoSwb7lV1C/DYorZ/qqonu91bgQ3d9hbgqqr6n6r6KrAHOG2M9UqShjCON1R/C/hUt70eeGTBsb1d20GSbE+yO8nu+fn5MZQhSfq+XuGe5I+AJ4GPr/TaqtpZVXNVNTczM9OnDEnSIiN/WibJG4HXAmdWVXXN+4CTFpy2oWuTJE3QSCP3JJuBdwGvq6pvLzi0CzgvyTOSbAQ2Af/av0xJ0kosO3JPciVwBrAuyV7gIgafjnkGcEMSgFur6i1VdU+Sq4F7GUzXXFBV31ut4iVJS1s23Kvq/CWaLz3M+e8F3tunKElSPy4/IEkNcvkBSU97fZc9ePDis8dUyfg4cpekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQq0JKGqu+KyxqPAx3SU9hOLfBaRlJapDhLkkNWjbck1yW5ECSuxe0PTfJDUm+0n0/vmtPkg8k2ZPkriSnrmbxkqSlDTNyvxzYvKhtB3BjVW0Cbuz2AV4DbOq+tgMfGk+ZkqSVWDbcq+oW4LFFzVuAK7rtK4BzFrR/tAZuBdYmOXFMtUqShjTqnPsJVbW/234UOKHbXg88suC8vV3bQZJsT7I7ye75+fkRy5AkLaX3G6pVVUCNcN3OqpqrqrmZmZm+ZUiSFhg13L/+/emW7vuBrn0fcNKC8zZ0bZKkCRo13HcBW7vtrcC1C9rf0H1q5nTgiQXTN5KkCVn2CdUkVwJnAOuS7AUuAi4Grk6yDXgIOLc7/XrgLGAP8G3gTatQsyRpGcuGe1Wdf4hDZy5xbgEX9C1KktSPT6hKUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWpQr3BP8vtJ7klyd5IrkxyXZGOS25LsSfKJJMeOq1hJ0nBGDvck64HfBeaq6qXAMcB5wPuAS6rqhcC3gG3jKFSSNLy+0zJrgB9JsgZ4JrAfeDVwTXf8CuCcnveQJK3QyOFeVfuAPwMeZhDqTwC3A49X1ZPdaXuB9Utdn2R7kt1Jds/Pz49ahiRpCX2mZY4HtgAbgZ8AfhTYPOz1VbWzquaqam5mZmbUMiRJS+gzLfMLwFerar6qvgt8EngFsLabpgHYAOzrWaMkaYXWLH/KIT0MnJ7kmcB/A2cCu4GbgNcDVwFbgWv7FilpZWZ3XDftEjRlfebcb2PwxunngS92r7UTeDfwjiR7gOcBl46hTknSCvQZuVNVFwEXLWp+ADitz+tKkvrxCVVJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQb3CPcnaJNck+VKS+5L8fJLnJrkhyVe678ePq1hJ0nD6jtzfD/xjVb0EeBlwH7ADuLGqNgE3dvuSpAlaM+qFSZ4DvBJ4I0BVfQf4TpItwBndaVcANwPv7lOkJB3JZndcN/K1D1589hgr+YGRwx3YCMwDf5XkZcDtwNuBE6pqf3fOo8AJS12cZDuwHeDkk0/uUYbUpj6BIfWZllkDnAp8qKpeDvwXi6ZgqqqAWuriqtpZVXNVNTczM9OjDEnSYn3CfS+wt6pu6/avYRD2X09yIkD3/UC/EiVJKzVyuFfVo8AjSV7cNZ0J3AvsArZ2bVuBa3tVKElasT5z7gBvAz6e5FjgAeBNDP7DuDrJNuAh4Nye95AkrVCvcK+qLwBzSxw6s8/rSpL68QlVSWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoP6PqEq6TBc2VHT4shdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CAfYpKW4YNIOhr1HrknOSbJHUn+odvfmOS2JHuSfKL7+6qSpAkax7TM24H7Fuy/D7ikql4IfAvYNoZ7SJJWoFe4J9kAnA18pNsP8Grgmu6UK4Bz+txDkrRyfUfufwG8C/jfbv95wONV9WS3vxdY3/MekqQVGjnck7wWOFBVt494/fYku5Psnp+fH7UMSdIS+nxa5hXA65KcBRwHPBt4P7A2yZpu9L4B2LfUxVW1E9gJMDc3Vz3qkA7LT7vo6WjkkXtVXVhVG6pqFjgP+ExV/TpwE/D67rStwLW9q5QkrchqPMT0buAdSfYwmIO/dBXuIUk6jLE8xFRVNwM3d9sPAKeN43UlSaNx+QFJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQSOHe5KTktyU5N4k9yR5e9f+3CQ3JPlK9/348ZUrSRpGn5H7k8A7q+oU4HTggiSnADuAG6tqE3Bjty9JmqCRw72q9lfV57vt/wDuA9YDW4ArutOuAM7pWaMkaYXGMueeZBZ4OXAbcEJV7e8OPQqccIhrtifZnWT3/Pz8OMqQJHV6h3uSZwF/C/xeVf37wmNVVUAtdV1V7ayquaqam5mZ6VuGJGmBNX0uTvLDDIL941X1ya7560lOrKr9SU4EDvQtUprdcd20S5COKn0+LRPgUuC+qvrzBYd2AVu77a3AtaOXJ0kaRZ+R+yuA3wS+mOQLXdt7gIuBq5NsAx4Czu1VoZrh6FuanJHDvar+GcghDp856utKkvrzCVVJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg3qtLaOjU58nRR+8+OwxViJptThyl6QGGe6S1CDDXZIaZLhLUoMMd0lqkJ+W6WFanzqZ5rrorskuHR0cuUtSgxy5T4kjYEmryZG7JDXoqB+59x0B+8SlpBat2sg9yeYk9yfZk2THat1HknSwVRm5JzkG+CDwi8Be4HNJdlXVvatxvz6c+5bUotUauZ8G7KmqB6rqO8BVwJZVupckaZHVmnNfDzyyYH8v8HMLT0iyHdje7f5nkvtXqZbVtg74xrSLOMLYJwezT57K/ujkff+/OUqfvOBQB6b2hmpV7QR2Tuv+45Jkd1XNTbuOI4l9cjD75Knsj4ONu09Wa1pmH3DSgv0NXZskaQJWK9w/B2xKsjHJscB5wK5VupckaZFVmZapqieT/A7waeAY4LKqumc17nUEOOqnllaBfXIw++Sp7I+DjbVPUlXjfD1J0hHA5QckqUGGuyQ1yHAf0jDLKSQ5N8m9Se5J8teTrnHSluuTJJck+UL39eUkj0+hzIkZoj9OTnJTkjuS3JXkrGnUOUlD9MkLktzY9cfNSTZMo85JSXJZkgNJ7j7E8ST5QNdfdyU5deSbVZVfy3wxeFP434CfBI4F7gROWXTOJuAO4Phu//nTrnvafbLo/LcxeGN96rVP8d/ITuC3u+1TgAenXfcR0Cd/A2zttl8NfGzada9yn7wSOBW4+xDHzwI+BQQ4Hbht1Hs5ch/OMMspvBn4YFV9C6CqDky4xklb6RIT5wNXTqSy6RimPwp4drf9HOBrE6xvGobpk1OAz3TbNy1xvClVdQvw2GFO2QJ8tAZuBdYmOXGUexnuw1lqOYX1i855EfCiJP+S5NYkmydW3XQM0yfA4FdvYCM/+CFu0TD98cfAbyTZC1zP4LeZlg3TJ3cCv9pt/wrwY0meN4HajlRD/1wtx3AfnzUMpmbOYDBK/XCStdMs6AhyHnBNVX1v2oVM2fnA5VW1gcGv3x9L8nT/GfwD4FVJ7gBexeBJ9qf7v5OxOOr/WMeEDLOcwl4G82PfBb6a5MsMwv5zkylx4layxMR5wAWrXtF0DdMf24DNAFX12STHMVgsqtUpvGX7pKq+RjdyT/Is4Neq6vFJFXgEGtvSLU/3UcOwhllO4e8ZjNpJso7BNM0DE6xx0oZaYiLJS4Djgc9OuL5JG6Y/HgbOBEjy08BxwPxEq5ysZfskyboFv71cCFw24RqPNLuAN3SfmjkdeKKq9o/yQob7EKrqSeD7yyncB1xdVfck+ZMkr+tO+zTwzST3Mnhj6A+r6pvTqXj1DdknMPiBvqq6jwK0asj+eCfw5iR3Mnhz+Y0t98uQfXIGcH/3m+4JwHunUuyEJLmSwUDnxUn2JtmW5C1J3tKdcj2DQeEe4MPAW0e+V8P/tiTpacuRuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDfo/sWi/HgPWy4kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prop_pairs_prior = {'alpha': 15, 'beta': 2}\n",
    "r = beta.rvs(prop_pairs_prior['alpha'], prop_pairs_prior['beta'], size=1000)\n",
    "\n",
    "plt.hist(r, bins=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perform ABC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[26, 36, 36, 42, 35, 36, 37, 41, 42, 54, 46, 41, 55, 39]"
      ]
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import math\n",
    "from collections import Counter\n",
    "from statistics import mean\n",
    "\n",
    "def simulation(niter=10):\n",
    "    n_socks_prior = {'r': 30, 'p': 0.5}\n",
    "    prop_pairs_prior = {'alpha': 15, 'beta': 2}\n",
    "    n_pick = 11\n",
    "    \n",
    "    posterior = []\n",
    "    \n",
    "    for i in range(niter):\n",
    "        # sample one value from each of the priors\n",
    "        n_socks = nbinom.rvs(n_socks_prior['r'], n_socks_prior['p'], size=1)[0]\n",
    "        prop_pairs = beta.rvs(prop_pairs_prior['alpha'], prop_pairs_prior['beta'], size=1)[0]\n",
    "        \n",
    "        n_pairs = int(round(math.floor(n_socks/2)*prop_pairs))\n",
    "        n_odd = n_socks - n_pairs * 2\n",
    "        n_sock_types = n_pairs + n_odd\n",
    "        \n",
    "        if n_pick > n_odd + 2 * n_pairs:\n",
    "            continue\n",
    "        \n",
    "        # run experiment\n",
    "        picked_socks = pick_socks(n_pairs=n_pairs, n_odd=n_odd, n_pick=n_pick)\n",
    "\n",
    "        # compute metric: count the number of unique socks\n",
    "        counts = Counter(picked_socks)\n",
    "        unique = len([x for x in counts if counts[x] == 1])\n",
    "        pairs = len([x for x in counts if counts[x] == 2])\n",
    "        \n",
    "        assert pairs*2 + unique == n_pick\n",
    "        \n",
    "        # accept or reject experiment\n",
    "        if unique == 11:\n",
    "            posterior.append(n_socks)\n",
    "    \n",
    "    return posterior\n",
    "\n",
    "posterior = simulation(100)\n",
    "posterior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## visualize posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "33\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([  1.,  11.,  41.,  74., 159., 249., 357., 252., 551., 616., 683.,\n",
       "        599., 510., 478., 210., 311., 254., 171., 128.,  81.,  58.,  28.,\n",
       "          8.,  20.,   8.,   7.,   4.,   0.,   0.,   1.]),\n",
       " array([13.        , 14.86666667, 16.73333333, 18.6       , 20.46666667,\n",
       "        22.33333333, 24.2       , 26.06666667, 27.93333333, 29.8       ,\n",
       "        31.66666667, 33.53333333, 35.4       , 37.26666667, 39.13333333,\n",
       "        41.        , 42.86666667, 44.73333333, 46.6       , 48.46666667,\n",
       "        50.33333333, 52.2       , 54.06666667, 55.93333333, 57.8       ,\n",
       "        59.66666667, 61.53333333, 63.4       , 65.26666667, 67.13333333,\n",
       "        69.        ]),\n",
       " <BarContainer object of 30 artists>)"
      ]
     },
     "execution_count": 168,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAARpUlEQVR4nO3db4zdVZ3H8fdHEHXVWIHZpmnrFmMj4cFScMLWaIzSaASM5YGyGFcqadJ9wG40unGrT8xudpP6RMRkQ9KIbtmoiCihUaKSqnF3E1inyPqvEitbtm2AjsgflahBv/vgnsql2zJ3Zu7McA/vV3Jzz+/8zm/uOczlM6fn/n6/m6pCktSX5610ByRJ42e4S1KHDHdJ6pDhLkkdMtwlqUOnr3QHAM4+++zasGHDSndDkibK/v37f15VUyfb96wI9w0bNjAzM7PS3ZCkiZLk/lPtm3NZJsmrk9wz9Hg8yfuTnJnkjiQ/bc8vb+2T5JNJDib5fpILxzkYSdLc5gz3qrq3qjZV1SbgNcATwK3ATmBfVW0E9rVtgEuAje2xA7h+CfotSXoG8/1AdQvws6q6H9gK7Gn1e4DLW3krcGMN3AmsSrJmHJ2VJI1mvuF+JfD5Vl5dVQ+08oPA6lZeCxweOuZIq5MkLZORwz3JGcDbgS+euK8GN6iZ101qkuxIMpNkZnZ2dj6HSpLmMJ+Z+yXA3VX1UNt+6PhyS3s+1uqPAuuHjlvX6p6mqnZX1XRVTU9NnfRMHknSAs0n3N/FU0syAHuBba28DbhtqP6qdtbMZuCxoeUbSdIyGOk89yQvBt4M/PVQ9S7g5iTbgfuBK1r97cClwEEGZ9ZcPbbeSpJGMlK4V9WvgbNOqHuYwdkzJ7Yt4Jqx9E6StCDPiitU1Z8NO786UrtDuy5b4p5Iz03eOEySOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkN+hqnkZ9btRJa0sZ+6S1KGRwj3JqiS3JPlJkgNJXpvkzCR3JPlpe355a5skn0xyMMn3k1y4tEOQJJ1o1Jn7dcDXqupc4HzgALAT2FdVG4F9bRvgEmBje+wArh9rjyVJc5oz3JO8DHgDcANAVf2uqh4FtgJ7WrM9wOWtvBW4sQbuBFYlWTPmfkuSnsEoH6ieA8wCn0lyPrAfeB+wuqoeaG0eBFa38lrg8NDxR1rdA0N1JNnBYGbPK17xioX2XxNu1A9oD+26bIl7IvVllGWZ04ELgeur6gLg1zy1BANAVRVQ83nhqtpdVdNVNT01NTWfQyVJcxgl3I8AR6rqrrZ9C4Owf+j4ckt7Ptb2HwXWDx2/rtVJkpbJnOFeVQ8Ch5O8ulVtAX4M7AW2tbptwG2tvBe4qp01sxl4bGj5RpK0DEa9iOlvgc8mOQO4D7iawR+Gm5NsB+4HrmhtbwcuBQ4CT7S2kqRlNFK4V9U9wPRJdm05SdsCrllctyRJi+EVqpLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CG/Q1V+L6rUIWfuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtSh0a6t0ySQ8Avgd8DT1bVdJIzgS8AG4BDwBVV9UiSANcBlwJPAO+tqrvH33U9l4x6/5tDuy5b4p5Ik2E+M/c3VdWmqppu2zuBfVW1EdjXtgEuATa2xw7g+nF1VpI0msUsy2wF9rTyHuDyofoba+BOYFWSNYt4HUnSPI0a7gV8I8n+JDta3eqqeqCVHwRWt/Ja4PDQsUda3dMk2ZFkJsnM7OzsArouSTqVUe/n/vqqOprkT4E7kvxkeGdVVZKazwtX1W5gN8D09PS8jpUkPbORZu5VdbQ9HwNuBS4CHjq+3NKej7XmR4H1Q4eva3WSpGUyZ7gneXGSlx4vA28BfgjsBba1ZtuA21p5L3BVBjYDjw0t30iSlsEoyzKrgVsHZzhyOvC5qvpaku8CNyfZDtwPXNHa387gNMiDDE6FvHrsvZYkPaM5w72q7gPOP0n9w8CWk9QXcM1YeidJWhCvUJWkDo16tow0EbySVRpw5i5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1aORwT3Jaku8l+UrbPifJXUkOJvlCkjNa/Qva9sG2f8MS9V2SdArzmbm/DzgwtP0x4NqqehXwCLC91W8HHmn117Z2kqRlNFK4J1kHXAZ8qm0HuBi4pTXZA1zeylvbNm3/ltZekrRMRp25fwL4EPCHtn0W8GhVPdm2jwBrW3ktcBig7X+stZckLZM5wz3J24BjVbV/nC+cZEeSmSQzs7Oz4/zRkvScN8rM/XXA25McAm5isBxzHbAqyemtzTrgaCsfBdYDtP0vAx4+8YdW1e6qmq6q6ampqUUNQpL0dHOGe1V9uKrWVdUG4Ergm1X1buBbwDtas23Aba28t23T9n+zqmqsvZYkPaPFnOf+98AHkhxksKZ+Q6u/ATir1X8A2Lm4LkqS5uv0uZs8paq+DXy7le8DLjpJm98A7xxD3yRJC+QVqpLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOzesiJj07bNj51ZHaHdp12RL3RNKzlTN3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDs15V8gkLwS+A7ygtb+lqj6a5BzgJuAsYD/wnqr6XZIXADcCrwEeBv6yqg4tUf+lJTXqHTjBu3Dq2WWUmftvgYur6nxgE/DWJJuBjwHXVtWrgEeA7a39duCRVn9taydJWkZzhnsN/KptPr89CrgYuKXV7wEub+WtbZu2f0uSjKvDkqS5jbTmnuS0JPcAx4A7gJ8Bj1bVk63JEWBtK68FDgO0/Y8xWLo58WfuSDKTZGZ2dnZRg5AkPd1I4V5Vv6+qTcA64CLg3MW+cFXtrqrpqpqemppa7I+TJA2Z19kyVfUo8C3gtcCqJMc/kF0HHG3lo8B6gLb/ZQw+WJUkLZM5wz3JVJJVrfwi4M3AAQYh/47WbBtwWyvvbdu0/d+sqhpjnyVJcxjlC7LXAHuSnMbgj8HNVfWVJD8GbkryT8D3gBta+xuAf0tyEPgFcOUS9FuS9AzmDPeq+j5wwUnq72Ow/n5i/W+Ad46ld5KkBfEKVUnqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHRrlCVRNqPl80ocUb9b+3X+qh5WC4P4sYxpLGxWUZSeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ3PeWybJeuBGYDVQwO6qui7JmcAXgA3AIeCKqnokSYDrgEuBJ4D3VtXdS9N9aWG8j496N8rM/Ungg1V1HrAZuCbJecBOYF9VbQT2tW2AS4CN7bEDuH7svZYkPaM5w72qHjg+866qXwIHgLXAVmBPa7YHuLyVtwI31sCdwKoka8bdcUnSqc1rzT3JBuAC4C5gdVU90HY9yGDZBgbBf3josCOt7sSftSPJTJKZ2dnZ+fZbkvQMRg73JC8BvgS8v6oeH95XVcVgPX5kVbW7qqaranpqamo+h0qS5jBSuCd5PoNg/2xVfblVP3R8uaU9H2v1R4H1Q4eva3WSpGUyZ7i3s19uAA5U1ceHdu0FtrXyNuC2ofqrMrAZeGxo+UaStAxG+Zq91wHvAX6Q5J5W9xFgF3Bzku3A/cAVbd/tDE6DPMjgVMirx9lhSdLc5gz3qvoPIKfYveUk7Qu4ZpH9kiQtgleoSlKHDHdJ6pDhLkkdGuUDVUljNOp9bQ7tumyJe6KeOXOXpA4Z7pLUIcNdkjrkmvsy8N7hkpabM3dJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtSh7yISXqW8gZjWgxn7pLUIcNdkjpkuEtShwx3SerQnOGe5NNJjiX54VDdmUnuSPLT9vzyVp8kn0xyMMn3k1y4lJ2XJJ3cKDP3fwXeekLdTmBfVW0E9rVtgEuAje2xA7h+PN2UJM3HnOFeVd8BfnFC9VZgTyvvAS4fqr+xBu4EViVZM6a+SpJGtNA199VV9UArPwisbuW1wOGhdkda3f+TZEeSmSQzs7OzC+yGJOlkFv2BalUVUAs4bndVTVfV9NTU1GK7IUkastBwf+j4ckt7PtbqjwLrh9qta3WSpGW00HDfC2xr5W3AbUP1V7WzZjYDjw0t30iSlsmc95ZJ8nngjcDZSY4AHwV2ATcn2Q7cD1zRmt8OXAocBJ4Arl6CPkuS5jBnuFfVu06xa8tJ2hZwzWI7JUlaHO8KKU047x6pk/H2A5LUIWfuizDqjEmSlpszd0nqkOEuSR0y3CWpQ4a7JHXID1Sl5whPmXxuceYuSR0y3CWpQ4a7JHXIcJekDvmB6kl45amkSefMXZI6ZLhLUocMd0nqkGvukp7Gi5364MxdkjpkuEtSh1yWkbTkXOpZfs7cJalDz5mZuxcmSXouWZJwT/JW4DrgNOBTVbVrKV5H0spZigmTyzfjM/ZwT3Ia8C/Am4EjwHeT7K2qH4/7tSQ9N/lHYG5LMXO/CDhYVfcBJLkJ2AosSbi73CLpVFbyj8BK/wFainBfCxwe2j4C/MWJjZLsAHa0zV8luXcJ+rKSzgZ+vtKdWCK9js1xTZ6xjC0fG0NPxvvao47rz061Y8U+UK2q3cDulXr9pZZkpqqmV7ofS6HXsTmuydPr2MYxrqU4FfIosH5oe12rkyQtk6UI9+8CG5Ock+QM4Epg7xK8jiTpFMa+LFNVTyb5G+DrDE6F/HRV/WjcrzMBul1yot+xOa7J0+vYFj2uVNU4OiJJehbx9gOS1CHDXZI6ZLiPQZL1Sb6V5MdJfpTkfa3+zCR3JPlpe375Svd1PpK8MMl/JfnvNq5/aPXnJLkrycEkX2gfnE+cJKcl+V6Sr7TtXsZ1KMkPktyTZKbVTfR7ESDJqiS3JPlJkgNJXtvJuF7dflfHH48nef9ix2a4j8eTwAer6jxgM3BNkvOAncC+qtoI7Gvbk+S3wMVVdT6wCXhrks3Ax4Brq+pVwCPA9pXr4qK8DzgwtN3LuADeVFWbhs6VnvT3IgzuV/W1qjoXOJ/B727ix1VV97bf1SbgNcATwK0sdmxV5WPMD+A2BvfWuRdY0+rWAPeudN8WMaY/Ae5mcLXxz4HTW/1rga+vdP8WMJ517X+Yi4GvAOlhXK3vh4CzT6ib6Pci8DLgf2gngfQyrpOM8y3Af45jbM7cxyzJBuAC4C5gdVU90HY9CKxeqX4tVFu6uAc4BtwB/Ax4tKqebE2OMLjlxKT5BPAh4A9t+yz6GBdAAd9Isr/d5gMm/714DjALfKYtpX0qyYuZ/HGd6Erg8628qLEZ7mOU5CXAl4D3V9Xjw/tq8Od34s47rarf1+Cfi+sY3BTu3JXt0eIleRtwrKr2r3Rflsjrq+pC4BIGS4RvGN45oe/F04ELgeur6gLg15ywTDGh4/qj9hnP24EvnrhvIWMz3MckyfMZBPtnq+rLrfqhJGva/jUMZr8TqaoeBb7FYLliVZLjF8BN4u0lXge8Pckh4CYGSzPXMfnjAqCqjrbnYwzWbi9i8t+LR4AjVXVX276FQdhP+riGXQLcXVUPte1Fjc1wH4MkAW4ADlTVx4d27QW2tfI2BmvxEyPJVJJVrfwiBp8jHGAQ8u9ozSZuXFX14apaV1UbGPwz+JtV9W4mfFwASV6c5KXHywzWcH/IhL8Xq+pB4HCSV7eqLQxuIz7R4zrBu3hqSQYWOTavUB2DJK8H/h34AU+t4X6Ewbr7zcArgPuBK6rqFyvSyQVI8ufAHga3kXgecHNV/WOSVzKY8Z4JfA/4q6r67cr1dOGSvBH4u6p6Ww/jamO4tW2eDnyuqv45yVlM8HsRIMkm4FPAGcB9wNW09yUTPC744x/i/wVeWVWPtbpF/c4Md0nqkMsyktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR16P8AH+QJ1BnwUZwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "posterior = simulation(50000)\n",
    "\n",
    "print(mean(posterior))\n",
    "\n",
    "plt.hist(posterior, bins=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}