2021_02_10.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Hypothesis testing\n",
    "\n",
    "We have two biased coins, referred to as coins 1 and 2, with probabilities of heads equal to $p_1$ and $p_2$, respectively. We choose a coin at random\n",
    "(either coin is equally likely to be chosen) and we want to infer its identity, based\n",
    "on the outcome of a several tosses. Let $\\Theta=1$ and $\\Theta=2$ be the hypotheses that coin\n",
    "1 or 2, respectively, was chosen. Let X be equal to 1 or 0, depending on whether\n",
    "the outcome of the toss was a head or a tail, respectively. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "text": [
      "WARNING (theano.gof.compilelock): Overriding existing lock by dead process '16320' (I am process '21600')\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "BinaryGibbsMetropolis: [coin]\n",
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 20 seconds.\n"
     ],
     "output_type": "stream"
    },
    {
     "data": {
      "text/plain": "<IPython.core.display.HTML object>",
      "text/html": "\n    <div>\n        <style>\n            /* Turns off some styling */\n            progress {\n                /* gets rid of default border in Firefox and Opera. */\n                border: none;\n                /* Needs to be in here for Safari polyfill so background images work as expected. */\n                background-size: auto;\n            }\n            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n                background: #F44336;\n            }\n        </style>\n      <progress value='0' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      \n    </div>\n    "
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaPklEQVR4nO3dd5hU1f3H8fdhqYquIIKCGESsWAAVQVRsRCMajA1rIBqV2LGLBRJFgwpizE9BLNjBjl0jlqBgQ9FYAmrAgoIi0pvA/f0xwzq7LLBtdnbPvl/PMw87Z86993vu3bOfvXfvDCFJEiRJik2tXBcgSVI2GHCSpCgZcJKkKBlwkqQoGXCSpCgZcJKkKBlwUhUQQkhCCANyXYcUk9q5LkASAJ2Bb3NdhBST4Bu9JUkx8hKlVE4hhF1CCE+EEH4KISwOIUwOIVyWfi2EEPqm25aFEL4PIfwzhLBhkXUUukQZQhiQbts6hPBsCGFBCOGrEMJVIQTnrVQCThSpHEIIHYEJwFZAX6A7MATYPN1lYPr5v4DDgOuB3sCzJQyqJ4BXgMOBJ4G/Ar0qqn4pZv4NTiqfG4GfgE5JkixKt70CEEJoDJwP3JMkyVnp114MIfwI3AccCjy1jvUPTpLk7vTXL4cQ9geOA+5eyzKS8AxOKrMQwnpAF+CBjHDL1AmoB9xfpH0UsBzoWoLNPFvk+cfAFqUsVaqRDDip7BqRmkNruvuxcfrf7zMbkyRZTuqsr/FqS6xudpHnS4H6pahRqrEMOKnsfgZWAi3W8PqqcNo0szGEUBvYmFTIScoSA04qo/RlyTeAE0MIDYrp8hapM65ji7T3JPX379ezW6FUs3mTiVQ+F5IKqgkhhMGkLle2BtolSXJ2CGEIcFkIYSHwHLA9cA2pYCz69zVJFciAk8ohSZJ3QwhdgL8Bt5C6qeQrfr3L8XLgR6APcAapy5L3ApclSbKy8iuWag4/yUSSFCX/BidJipIBJ0mKkgEnSYqSASdJipIBJ0mKUq7eJuCtm5KkVUI2VuoZnCQpSgacJClKBpwkKUoGnCQpSgacJClKBpwkKUoGnCQpSgacJClKBpwkKUoGnCQpSv6P3lKOjJk0PdclSFVCj3YtsrJez+AkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CTVKO/++2UOb785M7/7JuvbGn37UK46/ViO22u7StumfmXASVKWvPTY/axYsZyddtsz16XUSLVzsdExk6bnYrOSVKlGPP8OtWrV4t1/v8w7r7+U63JqHM/gJOXUzVf15YLjD+G9cWM564j9OKZzG64++4/Mn/sz3389lStOPZqenbfmguMPYdqUTwuWW7p4MSMGXUnvA9tz9B5bceEJ3flgwuuF1p0kCQ8NG0yv/Xfh2C7bMvSKc1m0cP5qNSxbuoSRQ6/hlIN356iOrTnvmG68N25sucdWq5Y/YnPJvS8p52bNmM6Dt93ICWdexBlXDOK/H07k1qsv4cZLz2Tvg3pw8Q3DWbFiOTdeeiZJkgDwf1dfxNinHuaoU87m0sEjaLLpZlxzTi8+/eCdgvU+89CdPHz7UH575AlccsNw6tavz71DB662/esvOp1Xn3qEo04+m8tvvps2bXfh2r4n87/JnxT0WblyJSuWL1/7Y8WK7O8slVhOLlFKUqb58+bw93vGsFnLVgBM+/wznrxnGOf+bSj7HXYUkDobu+acXnw79QsAxr0whrMHDGH/3x8NQPs99+W8Y7rx8IibGXDrA6xYsYLH776Ng448kRPOvLigT/8+x/HTDzMKtv3h22/w3rixXDPiEXbcrXOqX+eufPfVVB694x9cfMNwAG4ZcAGvPv3IWsfRdtdODLzj0YrbMSoXA05SzjVtvnlBuAEFX+/U8debMzbbYksAZv8wg9k/ziRJErp0O7Tg9Vq1arFnt+48MfI2AGbN/I6fZ82k476/LbStTvv/jg/fHlfw/KO3x9GoSVO2b7c7K5YvL2jfeY8uvPLUr4F2bJ/z6d6z91rH0WD9hiUbsCqFAScp59ZvmF/oee06dVPtG+RntNUBYNmypfw8ayb111ufeg0aFFpuo8absHTJYn5ZtpQ5s34AIL9xk0J98htvXOj5vDmz+XnWDxy5e6vV6qqVl1fw9SabtqBJ083WPpAQ1v66KpUBJ6naadSkGUsWLWTp4sWFQm7O7B+pV78BderWY6MmTQGYO3tWoWXnzv6p0POG+RuxcdNNuWzInWvdppcoqx8DTlK106btLoQQGP/ys4X+Rjf+5WfZvn1HAJo0a06jJk1557WX6NBlv4Jl33rl+ULr2rnjXoy573bqr7c+m2/ZZo3b9BJl9WPASap2Wrbemr0P7sHtg65g0cL5bNayFS89/iDTp31Jn37XAZCXl8cfev2FkTddzQYbNWaHDh2ZMPY5vp36eaF1teu0D+07d6X/X47jiN5nsMVW27BowQKmTvmEX5Yu5aRzLgOgWfOWNGveslR1fvzeBOb9/BNffvYfAN5/41XyGzWmZettaLnVNhWwJ7Q2BpykaunMK2/gnpsH8vCIm1k4fx6/abMdV9w8kh3SZ3AAh53wZ+bPncOLj97H0w/eQceu3fjjeZdzU7+zC/qEELh08AgevfMWnn7wDmZ9/x0N8zdiy213oPuxJ5erxoeGDeaTiW8VPB9+XT8Aep7el+O2uqBc69a6hVXvKalMYyZNr/yNSpKqpB7tWmTl7hzf6C1JipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBpxrvrVdf4NxjDuSojq05rXtnxtx3+1r733FDfw5vvzl3D7m6UPu3U7/gopMO5fi9t+fGS85g8aKFhV7/ZOJbnPzbXVdrL87Ypx7m8PabF9v3oWGDOWm/nQqez/zuGw5vv3nBo2fnrTn9sC4M6XcWn7z/9mrL33xVXy44/pB11iBVdwacarTPJr3LoAtOZeu27bj85rs5oEdP7v3HtTz1wB3F9v/myymMHTOa9RpusNpr/+jfl81atuKiQbfxzdQpPHrnLQWvrVy5kjtvHMBJZ19Gg/XWz8pYeve9kkH3jOHKW+7lmFPPZf7cn7n8lCMZNXxIVrYnVXW1c12AlEujbx/K9u06clb/GwFo37krC+fN5eHbb+J3x/yROnXqFuo/4vqrOPS4U3jt2ccKtS9etJAp//mAy4eOJL/xxiycP48n7x3OSWdfCsDLTz5EXu3a7HvokVkbS4tWW7HtzrsCsCOdOeD3x/DgrTcwatgQ2u7aiZ122zNr25aqIs/gVKNNnfwJu+yxV6G2dp33YcG8uUz+cGKh9vH/eoZvp37OkSefudp6lv+yDIC69eun/21Q0LZowXwevPVG/nzRXwkhZGMYa9Tz9PNpvEkzXnzk/krdrlQVGHCq0X5ZtpTadeoUaqtTN3XW9u3ULwrali5ZzN1Druakcy6jfoP1VlvPBvmNaNZiC5596G7mz/2Zlx5/gK122BmAh0cMZZc99io4uyqNlStWsGL58kKPZOXKEi+fl5fHTh27MPk/75d621J15yVK1WibtmzF5598WKjt848nATB/3pyCtsfu+ieNNmnKvt3XfInxtEuv4YaL+3D/P/9O8y225PRLB/L911N5+clRDH34X2Wq74R9dii2fYONGpV4HU2absbc2T+WaftSdWbAqUY7+KgTGXZtP156/AH2PLA7n388qeAuyrxaqQscM6d/zZP3Defq4aPXeolx1732556xHzLrh+/ZdPPfkJeXx8Bz/8RhJ5xKk2bNeW70SB4feSsAR/zpTA45ptc66xt452PUq1e/UNuLjz/AW688X+IxJklS4r5STAw41WgH9DiWaVM+Zdi1/bj16kuoV78Bfzy3HyMGXUn+xpsAcO8/rqPDnvvRYss2LJg/F4AkWckvy5ayYP5c1m+4YUHw1WvQgBa/aQ3Ah2+NY9qUT7lw0K1MnfwpD956A38fOQaAS3r1YId2u9Nqm+LP0FZpvd2Oq911+e64l0s1xp9+nEF+401KtYwUAwNONVpeXh6nXTqQ48+4iFkzv6dZiy2Ynv7b27Y7dQBg+rQvmTbl09XOmp4bPZLnRo/kjhfeoUmz5oVeW7F8OXfeOIBe511OvfoN+HjieHbavQubb9kGgJ07duHjiW+tM+DKa8Xy5fznnTdp26FTVrcjVUUGnAQ03HAjGm64EQDPP3Iv2+2yW0EYnXXVDSxeXPgN14MvPZO2u3bi4KNPIr/Rxqut74VH72P9DfPZ66DfF7QtXbL4168XL4JKuHQ4+vabmP3jTA466sSsb0uqagw41WiTP5rIZ5PeZctt2rJo4XzGvTCGDya8znV3PV7Qp03bXVZbrk7dejRp1rzY95YtmDeH0cNvov+tDxS0te3QiXuGXsvLT44C4KN3x3PSOf0qdCzTp33Jhhs1YvkvvzBz+je88eIY3h//Gsf2OZ8dd+tcoduSqgMDTjVaXu06vPHi04waNoRQqxY7tO/IdXc/Qautty/zOh+6bTC7d+3GVtv/+nFarbfbkV7n9eOB/7segN59r2DLbSv28uTIm1IfHVa3Xj0aNWnGtjt3YOCdj9G2wx4Vuh2pugi5uMNqzKTp3tYlSQKgR7sWWfkEBN/oLUmKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKUu1cF6DK882XU7j9+iuZ/NFE1m+YT7c/HEvP088nLy8v16VJleb7r6fyxL3DmPLR+3z95WS2b9+RgXc8muuylAUGXA2xYN4crupzHC1bb02/m+5ixjdfcfeQv5EkCSeceXGuy5MqzddfTmHiG6+w7U4dWL78l1yXoywy4GqIFx65n2VLl3Dp4BGs13AD6ASLFs5n1PAh/KHXX1JtUg2we9du7LHfQQAMuvA05s2ZneOKlC3+Da6GeP/NV2nfuWuhINv7oB4sW7KEjye+lcPKpMpVq5Y/9moKj3QNMX3aF7TYcqtCbZts1oJ69RswfdoXOapKkrLHgKshFsyfy/ob5K/W3nDDfBbMm5uDiiQpuwy4GiSE1duSJCEU94IkVXMGXA3RcIN8Fs6ft1r7ogXzWX+DDXNQkSRllwFXQ7Ro1YbpU78s1PbjjO9YsngRLVq1yVFVkpQ9BlwN0aHLfnww4TUWL1xQ0PbmS09Rt359dty1Uw4rk6Ts8H1wNcTBR5/Is6Pu4u8XnMoRvc9gxvSvGDVsCD1OPM33wKlGWbp4MRPfGAvA7B9msGjhAsb/6xkAdt3rAOo1aJDL8lSBQpIklb7RMZOmV/5GlfqorkFXpD6qa4N8Djz8OI7t40d1qWaZ+d03nN69c7GvDX92As2at6zkitSjXYus3OlmwEmScipbAeff4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVHysyhrqIeGDWb08JtyXYaUUz1P78txfS7IdRnKEs/gJElRMuAkSVHyw5YlSTnlhy1LklQKBpwkKUo5uYty0pg7vgTm5GLbFaQ58F2uiygnx5B71b1+qP5jqO71QwRjmDSGL/r379+mwlecJEmlPwYMGJDkYrvW7xiq0qO61x/DGKp7/Y5h7Q8vUUqSopSrgPtrjrZbUap7/eAYqoLqXj9U/zFU9/rBMaxRTt4mIElStnmJUpIUJQNOkhQlA06SFKUKD7gQQr0Qwi0hhFkhhIUhhKdCCJuvY5neIYSkmEf9Iv3OCCFMDSEsCSFMDCHsXUXqPzWEMC6EMDuEMCeE8GoIYa8ifQYUM74ZFVRzqfZLCGGnEMLrIYTFIYTpIYSrQgihSJ+u6XUtCSH8L4TQpyJqrYgxhBD2DSGMCSF8H0JYFEL4KIRwcjF9ivue2q4K1N9qDbUdXKRfVT4GxX0/r3o0TfeptGMQQtgnPVenp7fRuwTLVKl5UNoxVNF5UNoxZHcuVPT7DoDbSL3psBvQAXgNmATkrWWZ3sBCYNPMR5E+PYFfgFOB7YFbgAXAFlWg/geAs4D2wLbAsPR4ts7oMwD4b5ExblIB9ZZqvwAbAjOAh4EdgSOB+cAFGX22TNd/S3qdp6a3cWQ23qtShjH0A64BugCtgb8Ay4HjM/rsCyTADkX2+RqPYyXW3ypd20FFaqtbjY5Bw6LzNT1XXs3RMTgEuBY4ClgE9F5H/6o4D0o7hio1D8o4hqzOhYoeXD6wDDgho60lsBI4aC3L9QYWrGPdbwMjirR9DlyX6/qLWU9IT56zM9oGAB9n4RuqVPslPQnmAQ0y2q4ApvPrXbWDgM+LLHcHMKGi66+oY5v+QfVYxvNVE7tJNmou5zFYNal3W8s6q9UxSM+TFWv44Zr1Y1CklgUl+MFa5eZBacewhuVyNg/KeByyOhcq+hLlrkAd4KVVDUmSfAN8Buy5jmUbhBC+CiF8G0J4JoTQftULIYS66XW/VGSZl0qw3tIoT/2Z6gL1gZ+LtLdOn7pPDSGMCiG0Lk+xZdwvnYFxSZIszmh7kdTH/bTK6FN0nS8Cu4UQ6pSn5qIq8NhuyOr7G+C99CWcsSGE/cpY5hqVs/7HQwg/hBDeDCEcVeS16nYMTiH18XuPFfNaVo9BGVWpeVCBcjIPKkBW5kJFB9ympH6Lm1WkfWb6tTWZDJwM9ACOA5YAb4YQtk6/3gTIS6+nNOstrbLWX9Q1pH57eSqj7W1SZ6q/I3WKvSkwPoSwcVmLpWz7ZdM19F/12tr61E5vsyKV+9iGEA4FDgBuz2j+ntRv6UcCR5D6HhsbQtinvAUXUZb6FwAXAseQuqQzFhgdQjgxo0+1OQYhhFqk5u+9SZIszXipso5BWVS1eVBuOZ4HZZXVuVCiD1sOIVwDXL6Obmv7rSCQOg0tVpIkE4AJGdsbT+rvXmcD52R2Lc16M9aX1fqLbOtc4HTgwCRJ5q1qT5Lk+SL93gL+B/QChpRk3WtR2v1SXP+i7SXpU5HKemy7AA8C5yRJ8k7BypJkMqnJvMqEEEIrUpPp3+WudnUlrj9JklnA4Iym90IITYCLgfvXsc7i2itKmY4BqR9MLUldNvp1ZZV/DEqrKs6DMqlC86BUsj0XSvq/CQwtsrHifA10IvWbYBPgx4zXmlKKnZkkyYoQwnvAqjO4WaTOrIr+NtmU1ZO9OEOphPrT4XYN8LvMb7LiJEmyIITwCb+OsSzKsl9mrKE/Gcusqc9y4KcyVbpmZT62IXWn6nPAVUmS3FaCbb0NHFuWIteivN+bq7wN/CnjebU4BmmnAuOTJPmkBH2zcQzKoqrNgzKrIvOgIlXYXCjRJcokSWYlSfLfdTwWARNJ3d3SbdWyIXWL/fbA+JJsK71MAHYmdXpNkiTL0uvuVqRrt5KstzLqDyGcDwwEuidJ8kYJxlgf2G7VGMuijPtlArB3KPwWjG6k7hydltHnwGLW+V6SJL+Utd7ilPXYpi+xPA/8NUmSoSXcXDvKsb+LU97vzQztKFxblT8GACGE5kB3YEQJN9eOCj4GZVSl5kFZVZV5UMHaUVFzIQt3ztxG6k6kA0ndNv8qRW6zJ3Wd9bqM5/1J3SbaOj24u0gFTceMPj1J3eH4Z1KBczOp67e/qQL1X5Su7RgK3+qan9HnRqArqVte9wCeIXUXV7nqX9d+Aa4Dxmb0zyf1G9EoUrdHH5Guo7jbo4em1/nn9DayeYt6acawb7q+G4rs700y+pwHHE7qDLlteh0JcEQVqL8XcHy677akLhctA/pWl2OQsdwVwFxgvWJeq8xj0JDUz452pG5Pvyr99RbVaB6UdgxVah6UcQxZnQvZGGB9Uu9X+Ck9wKeBlkX6TANGZjy/CfgKWAr8QOoOmc7FrPuM9LJLSf3GuU8VqX9a+pum6COzzyhSvx0uIxWgjwE7VFDNa9wvwEhgWpH+O5G65LqE1G9K/UnfGp3RpyvwfnqdU4E+2ZgQZRlD+nlx+zuzz8XAF8BiYDYwDjikitTfC/g0PWnnAe8BJxazzip7DNJtIV3XrWtYX6UdA369Hb7YOVgd5kFpx1BF50Fpx5DVueD/JiBJipKfRSlJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYqSASdJipIBJ0mKkgEnSYrS/wOvbkWK4CfL6AAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "text": [
      "       mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_mean   ess_sd  \\\ncoin  0.532  0.499     0.0      1.0      0.004    0.003   14075.0  14075.0   \n\n      ess_bulk  ess_tail  r_hat  \ncoin   14075.0    4000.0    1.0  \n"
     ],
     "output_type": "stream"
    }
   ],
   "source": [
    "import pymc3 as pm\n",
    "import theano\n",
    "import numpy as np\n",
    "import arviz as az\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# These are the initial probabilities\n",
    "p_1 = 0.46\n",
    "p_2 = 0.52\n",
    "\n",
    "# This is the data\n",
    "number_tosses = 1\n",
    "heads = 1\n",
    "\n",
    "with pm.Model():\n",
    "    # Here we indicate that we have 2 coins, and the probability\n",
    "    # of selecting each of the coins is 0.5\n",
    "    coin = pm.Categorical('coin', p=[0.5, 0.5])\n",
    "    \n",
    "    # This is just a way to tell that if we select the first coin,\n",
    "    # the probability of getting heads is p_1 and if we select\n",
    "    # the second coin, the probability of heads is p_2\n",
    "    p = theano.shared(np.array([p_1, p_2]))[coin]\n",
    "    \n",
    "    # Here the likelihood is binomial, because we don't have the result\n",
    "    # for each individual toss\n",
    "    y = pm.Binomial('y', p=p, n=number_tosses, observed=[heads])\n",
    "    \n",
    "    trace = pm.sample(1000, progressbar=True)\n",
    "    \n",
    "    # The mode o\n",
    "    az.plot_posterior(trace, point_estimate='mode')\n",
    "\n",
    "    plt.show()\n",
    "    \n",
    "    print(az.summary(trace))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "source": [
    "This indicates, that, according to the data, the selected coin was the second one.\n",
    "\n",
    "**NOTE:** For problems with only discrete random variables, don't use the find_MAP function, it will break.\n",
    "\n",
    "If now we have 50 tosses, 30 of them being tail"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "BinaryGibbsMetropolis: [coin]\n",
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 18 seconds.\n"
     ],
     "output_type": "stream"
    },
    {
     "data": {
      "text/plain": "<IPython.core.display.HTML object>",
      "text/html": "\n    <div>\n        <style>\n            /* Turns off some styling */\n            progress {\n                /* gets rid of default border in Firefox and Opera. */\n                border: none;\n                /* Needs to be in here for Safari polyfill so background images work as expected. */\n                background-size: auto;\n            }\n            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n                background: #F44336;\n            }\n        </style>\n      <progress value='0' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      \n    </div>\n    "
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbHklEQVR4nO3dd5hU1f3H8fdhKUuRFUVQEEMziljAgiAqNqIGE4waRdFATIxEYy+xQxQ0ShFjoiIW1KhYI7aoEU2Cgg1Foyag/sCCgiJSll7u748Z1m3AttnZPft+Pc88MmfOPfd75u7Zz8zdO2NIkgRJkmJTL9sFSJKUCQacJClKBpwkKUoGnCQpSgacJClKBpwkKUoGnFQDhBCSEMKwbNchxaR+tguQBEAv4ItsFyHFJPhBb0lSjDxFKVVSCGGPEMLfQgjfhhBWhBBmhhAuTT8WQgjnpdtWhxC+CiH8OYTQvNgYRU5RhhCGpdt2DCE8E0LIDyF8GkK4KoTgupXKwIUiVUIIoQcwDegEnAf0A8YA26e7jEjf/wfwE+AGYDDwTBmD6m/AS8DRwBPAH4BBVVW/FDP/BidVzijgW6BnkiTL020vAYQQtgLOB+5JkuR36ceeDyF8A9wHHAU8uZnxRydJcnf63y+GEA4BTgTu3sQ2kvAdnFRhIYQmQG/g/kLhVlhPoBHw12LtE4G1QJ8y7OaZYvffB3YoZ6lSnWTASRXXgtQa2tjVj1ul//tV4cYkSdaSete3VYktSlpY7P4qILccNUp1lgEnVdx3wHqg7UYe3xBO2xZuDCHUB7YmFXKSMsSAkyoofVryFeDkEELjUrq8Ruod14Bi7SeQ+vv3vzJboVS3eZGJVDkXkgqqaSGE0aROV3YEuiVJclYIYQxwaQhhGfAs0AUYTioYi/99TVIVMuCkSkiS5M0QQm/gauBmUheVfMr3VzleDnwDDAHOIHVa8l7g0iRJ1ld/xVLd4TeZSJKi5N/gJElRMuAkSVEy4CRJUTLgJElRMuAkSVHK1scEvHRTkrRByMSgvoOTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEUpK1/VNWnG3GzsVqpR+ndrm+0SpKj5Dk6SFCUDTpIUJQNOkhQlA06SFCUDTpIUJQNOkhQlA06SFCUDTpIUJQNOkhQlA06SFCUDTpIUJQNOkhQlA06SFCUDTpIUJQNOkhQlA06SFCUDTpIUJQNOkhQlA06SFCUDTpIUJQNOkhQlA06SFCUDTpIUJQNOkhQlA06SFCUDrg55898vcnT37Zn/5ecZ39ea1au4a/TVDDpkD07otSPXnPWLatmvJG1gwCkjxt9wFS8/9TCDz7uSi0eOY8mihQwbciKrV63MdmmS6ggDTlVuwfwvefGJiZx6wTAO/slx7LX/IVwyejzffDWXfz3zeLbLk1RHGHBZdNNV53HBST/mrSmT+d0xB3N8r85cc9YvWLr4O776bDZXnPZzTui1Ixec9GPmzPqwYLtVK1Yw/vorGXxYd36+bycuHNiPd6b9q8jYSZLw4G2jGXTIHgzovRNjrziH5cuWlqhh9aqVTBg7nF8dsQ/H9ejIucf35a0pkys1rxnT/g1Az0OPLGjbutV2dOm+D9NffblSY0tSWRlwWbZg3lweuHUUA8+8iDOuuJ7/vTudW675PaMuOZMDDu/PxSPHsW7dWkZdciZJkgDwl2suYvKTD3Pcr87iktHjabntdgw/exAfvvNGwbhPP3gnD98+lh8dO5DfjxxHw9xc7h07osT+b7jodF5+8hGOO/UsLr/pbjp33YNrzzuV/5v5QUGf9evXs27t2k3f1q0r6P/FnI/ZutV2NG7StMi+tu+wI3PnfFLVT6Eklap+tguo65YuWcQf75nEdu3aAzDno//yxD23cc7VYzn4J8cBqXdjw88exBezPwZgynOTOGvYGA756c8B6L7fQZx7fF8eHn8Tw265n3Xr1vH43bdy+LEnM/DMiwv6DB1yIt9+Pa9g3+++/gpvTZnM8PGPsOvevVL9evXhy09n8+gdf+LikeMAuHnYBbz81CObnEfXvXoy4o5HAVi2ZDFNt2heok+z5nnkL1lc0adKksrFgMuyVm22Lwg3oODfu/XY7/u2HToAsPDreSz8Zj5JktC771EFj9erV4/9+vbjbxNuBVJ/A/tuwXx6HPSjIvvqeciRvPv6lIL7770+hRYtW9Gl2z6sW7u2oH33fXvz0pPfB9qAIefT74TBm5xH46bNitwPIZTokyQJpTRLUkYYcFnWtFlekfv1GzRMtW+RV6itAQCrV6/iuwXzyW3SlEaNGxfZbsuttmHVyhWsWb2KRQu+BiBvq5ZF+uRttXWR+0sWLeS7BV9z7D7tS9RVLyen4N/bbNuWlq222/RECiVX0+Z5LFu6pESXZUuXFJmXJGWSAVfLtGjZmpXLl7FqxYoiIbdo4Tc0ym1Mg4aN2LJlKwAWL1xQZNvFC78tcr9Z3pZs3WpbLh1z5yb3Wd5TlNu378yC+V+ycsVychs3KejzxZyPadu+0+YnKUlVwICrZTp33YMQAlNffKbI3+imvvgMXbr3AKBl6za0aNmKN/75Anv2Prhg29de+nuRsXbvsT+T7rud3CZN2b5D543us7ynKLv1OrBgfwf1OxZInV7979tvcPqlJS90kaRMMOBqmXYdd+SAI/pz+/VXsHzZUrZr154XHn+AuXM+Ychl1wGQk5PDzwb9lgk3XsMWW27FLnv2YNrkZ/li9kdFxurW80C69+rD0N+eyDGDz2CHTj9keX4+s2d9wJpVqzjl7EsBaN2mHa3btCtzjS1bt+Gwowdw56hhJElCXoutmThuDNts15Y+/Y6puidDkjbBgKuFzrxyJPfcNIKHx9/EsqVL+EHnnbnipgnskn4HB/CTgb9m6eJFPP/ofTz1wB306NOXX5x7OTdedlZBnxACl4wez6N33sxTD9zBgq++pFnelnTYaRf6DTi1UjWedvHV5OY24e7RV7Nq5Qq67tWTC677Cw0b5VZqXEkqq7Dhs1XVadKMudW/U6mG6d+tbbZLkGqKjFxf7Qe9JUlRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuBU57328nOcc/xhHNejI7/p14tJ992+yf53jBzK0d235+4x1xRp/2L2x1x0ylGcdEAXRv3+DFYsX1bk8Q+mv8apP9qrRHtpJkyYQAiB/Pz8Eo8NGzaMli1bFtyfM2cOIYSCW9OmTenUqRMDBw5kypQpJbYfPHgwe++992ZrkGq7+tkuQMqm/854k+svOI1D+5/A4POuZNZ/3uHeP11LqFePnw78dYn+n38yi8mTHqJJsy1KPPanoeexXbv2nPTbC5kwdjiP3nkzp5x1CQDr16/nzlHDOOWsS2ncpCkAk2bM3Whdb3/2HQBPv/dlQf8N/jdvCavXrS/Yfv6X8wAYfN6VdOm2N2tWr2b+l5/zyvOTOPDAAxkw5HwGnH5+wfafLVzOouVrNrl/qTr179Y2I+MacKrTHrp9LF269eB3Q0cB0L1XH5YtWczDt9/Ikcf/ggYNGhbpP/6GqzjqxF/xz2ceK9K+YvkyZv3nHS4fO4G8rbZm2dIlPHHvuIKAe/GJB8mpX5+Djjo2Y3Np274TO+2+FwC70otDf3o8D9wykom3jaHrXj3Zbe/9MrZvqSbyFKXqtNkzP2CPffcv0tat14HkL1nMzHenF2mf+o+n+WL2Rxx76pklxlm7ZjUADXNz0/9tXNC2PH8pD9wyil9f9AdCCJmYxkadcPr5bLVNa55/5K/Vul+pJjDgVKetWb2K+g0aFGlr0DD1ru2L2R8XtK1auYK7x1zDKWdfSm7jJiXG2SKvBa3b7sAzD97N0sXf8cLj99Npl90BeHj8WPbYd/+Cd1flsX7dOtatXVvklqxfX+btc3Jy2K1Hb2b+5+1y71uq7TxFqTpt23bt+eiDd4u0ffT+DACWLllU0PbYXX+mxTatOKjfxk8x/uaS4Yy8eAh//fMfabNDB06/ZARffTabF5+YyNiH/1Gh+gYeuEup7Vts2aLMY7RstR2LF35Tof1LtZkBpzrtiONO5rZrL+OFx+9nv8P68dH7MwquosyplzrBMX/uZzxx3ziuGffQJk8x7rX/Idwz+V0WfP0V227/A3Jychhxzi/5ycDTaNm6Dc8+NIHHJ9wCwDG/PJMfHz9os/WNuPMxGjXKLdL2/OP389pLfy/zHJMkKXNfKSYGnOq0Q/sPYM6sD7nt2su45Zrf0yi3Mb845zLGX38leVtvA8C9f7qOPfc7mLYdOpO/dDEASbKeNatXkb90MU2bNS8IvkaNG9P2Bx0BePe1KcyZ9SEXXn8Ls2d+yAO3jOSPEyYB8PtB/dml2z60/2Hp79A26LjzriWuonxzyovlmuO338wjb6ttyrWNFAMDTnVaTk4Ov7lkBCedcREL5n9F67Y7MDf9t7eddtsTgLlzPmHOrA9LvGt69qEJPPvQBO547g1atm5T5LF1a9dy56hhDDr3chrlNub96VPZbZ/ebN+hMwC79+jN+9Nf22zAVda6tWv5zxuv0nXPnhndj1QTGXAS0Kz5ljRrviUAf3/kXnbeY++CMPrdVSNZsaLoh7NHX3ImXffqyRE/P4W8FluXGO+5R++jafM89j/8pwVtq1au+P7fK5ZDNZw6fOj2G1n4zXwOP+7kjO9LqmkMONVpM9+bzn9nvEmHH3Zl+bKlTHluEu9M+xfX3fV4QZ/OXfcosV2Dho1o2bpNqZ8ty1+yiIfG3cjQW+4vaOu6Z0/uGXstLz4xEYD33pzKKWdfVqVzmTvnE5pv2YK1a9Ywf27qg95vT/0nA4acz65796rSfUm1gQGnOi2nfgNeef4pJt42hlCvHrt078F1d/+N9jt2qfCYD946mn369KVTl90K2jruvCuDzr2M+/9yAwCDz7uCDjtV7enJCTemvjqsYaNGtGjZmp1235MRdz5G1z33rdL9SLVFyMYVVpNmzPWyLkkSAP27tc3INyD4QW9JUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpTqZ7sAVZ/PP5nF7Tdcycz3ptO0WR59fzaAE04/n5ycnGyXJlWbrz6bzd/uvY1Z773NZ5/MpEv3Hoy449Fsl6UMMODqiPwli7hqyIm067gjl914F/M+/5S7x1xNkiQMPPPibJcnVZvPPpnF9FdeYqfd9mTt2jXZLkcZZMDVEc898ldWr1rJJaPH06TZFtATli9bysRxY/jZoN+m2qQ6YJ8+fdn34MMBuP7C37Bk0cIsV6RM8W9wdcTbr75M9159igTZAYf3Z/XKlbw//bUsViZVr3r1/LVXV3ik64i5cz6mbYdORdq22a4tjXIbM3fOx1mqSpIyx4CrI/KXLqbpFnkl2ps1zyN/yeIsVCRJmWXA1SEhlGxLkoRQ2gOSVMsZcHVEsy3yWLZ0SYn25flLabpF8yxUJEmZZcDVEW3bd2bu7E+KtH0z70tWrlhO2/ads1SVJGWOAVdH7Nn7YN6Z9k9WLMsvaHv1hSdpmJvLrnv1zGJlkpQZfg6ujjji5yfzzMS7+OMFp3HM4DOYN/dTJt42hv4n/8bPwKlOWbViBdNfmQzAwq/nsXxZPlP/8TQAe+1/KI0aN85meapCIUmSat/ppBlzq3+nSn1V1/VXpL6qa4s8Djv6RAYM8au6VLfM//JzTu/Xq9THxj0zjdZt2lVzRerfrW1GrnQz4CRJWZWpgPNvcJKkKBlwkqQoGXCSpCgZcJKkKBlwkqQoGXCSpCgZcJKkKBlwkqQoGXCSpCj5XZR11IO3jeahcTdmuwwpq044/TxOHHJBtstQhvgOTpIUJQNOkhQlv2xZkpRVftmyJEnlYMBJkqKUlasoZ0y64xNgUTb2XUXaAF9mu4hKcg7ZV9vrh9o/h9peP0QwhxmT+Hjo0KGdq3zgJEmq/TZs2LAkG/u1fudQk261vf4Y5lDb63cOm755ilKSFKVsBdwfsrTfqlLb6wfnUBPU9vqh9s+httcPzmGjsvIxAUmSMs1TlJKkKBlwkqQoGXCSpChVecCFEBqFEG4OISwIISwLITwZQth+M9sMDiEkpdxyi/U7I4QwO4SwMoQwPYRwQA2p/7QQwpQQwsIQwqIQwsshhP2L9RlWyvzmVVHN5XpeQgi7hRD+FUJYEUKYG0K4KoQQivXpkx5rZQjh/0IIQ6qi1qqYQwjhoBDCpBDCVyGE5SGE90IIp5bSp7SfqZ1rQP3tN1LbEcX61eRjUNrP84Zbq3SfajsGIYQD02t1bnofg8uwTY1aB+WdQw1dB+WdQ2bXQlV/7gC4ldSHDvsCewL/BGYAOZvYZjCwDNi28K1YnxOANcBpQBfgZiAf2KEG1H8/8DugO7ATcFt6PjsW6jMM+F+xOW5TBfWW63kBmgPzgIeBXYFjgaXABYX6dEjXf3N6zNPS+zg2E59VqcAcLgOGA72BjsBvgbXASYX6HAQkwC7FnvONHsdqrL99urbDi9XWsBYdg2bF12t6rbycpWPwY+Ba4DhgOTB4M/1r4joo7xxq1Dqo4BwyuhaqenJ5wGpgYKG2dsB64PBNbDcYyN/M2K8D44u1fQRcl+36SxknpBfPWYXahgHvZ+AHqlzPS3oRLAEaF2q7ApjL91fVXg98VGy7O4BpVV1/VR3b9C+qxwrd37CwW2ai5koegw2Leu9NjFmrjkF6nazbyC/XjB+DYrXkl+EXa41bB+Wdw0a2y9o6qOBxyOhaqOpTlHsBDYAXNjQkSfI58F9gv81s2ziE8GkI4YsQwtMhhO4bHgghNEyP/UKxbV4ow7jlUZn6C2sI5ALfFWvvmH7rPjuEMDGE0LEyxVbweekFTEmSZEWhtudJfd1P+0J9io/5PLB3CKFBZWourgqPbXNKPt8Ab6VP4UwOIRxcwTI3qpL1Px5C+DqE8GoI4bhij9W2Y/ArUl+/91gpj2X0GFRQjVoHVSgr66AKZGQtVHXAbUvqVdyCYu3z049tzEzgVKA/cCKwEng1hLBj+vGWQE56nPKMW14Vrb+44aRevTxZqO11Uu9UjyT1FntbYGoIYeuKFkvFnpdtN9J/w2Ob6lM/vc+qVOljG0I4CjgUuL1Q81ekXqUfCxxD6mdscgjhwMoWXExF6s8HLgSOJ3VKZzLwUAjh5EJ9as0xCCHUI7V+702SZFWhh6rrGFRETVsHlZbldVBRGV0LZfqy5RDCcODyzXTb1KuCQOptaKmSJJkGTCu0v6mk/u51FnB24a7lGbfQeBmtv9i+zgFOBw5LkmTJhvYkSf5erN9rwP8Bg4AxZRl7E8r7vJTWv3h7WfpUpYoe297AA8DZSZK8UTBYkswktZg3mBZCaE9qMf270tWWVOb6kyRZAIwu1PRWCKElcDHw182MWVp7VanQMSD1i6kdqdNG3w9W/cegvGriOqiQGrQOyiXTa6Gs/zeBscV2VprPgJ6kXgm2BL4p9FgryvFkJkmyLoTwFrDhHdwCUu+sir+abEXJZC/NWKqh/nS4DQeOLPxDVpokSfJDCB/w/RwroiLPy7yN9KfQNhvrsxb4tkKVblyFj21IXan6LHBVkiS3lmFfrwMDKlLkJlT2Z3OD14FfFrpfK45B2mnA1CRJPihD30wcg4qoaeugwmrIOqhKVbYWynSKMkmSBUmS/G8zt+XAdFJXt/TdsG1IXWLfBZhaln2ltwnA7qTeXpMkyer02H2Lde1blnGro/4QwvnACKBfkiSvlGGOucDOG+ZYERV8XqYBB4SiH8HoS+rK0TmF+hxWyphvJUmypqL1lqaixzZ9iuXvwB+SJBlbxt11oxLPd2kq+7NZSDeK1lbjjwFACKEN0A8YX8bddaOKj0EF1ah1UFE1ZR1UsW5U1VrIwJUzt5K6EukwUpfNv0yxy+xJnWe9rtD9oaQuE+2YntxdpIKmR6E+J5C6wvHXpALnJlLnb39QA+q/KF3b8RS91DWvUJ9RQB9Sl7zuCzxN6iquStW/uecFuA6YXKh/HqlXRBNJXR59TLqO0i6PHpse89fpfWTyEvXyzOGgdH0jiz3f2xTqcy5wNKl3yF3TYyTAMTWg/kHASem+O5E6XbQaOK+2HINC210BLAaalPJYdR6DZqR+d3QjdXn6Vel/71CL1kF551Cj1kEF55DRtZCJCeaS+rzCt+kJPgW0K9ZnDjCh0P0bgU+BVcDXpK6Q6VXK2Gekt11F6hXngTWk/jnpH5rit8J9JpJ6dbiaVIA+BuxSRTVv9HkBJgBzivXfjdQp15WkXikNJX1pdKE+fYC302POBoZkYkFUZA7p+6U934X7XAx8DKwAFgJTgB/XkPoHAR+mF+0S4C3g5FLGrLHHIN0W0nXdspHxqu0Y8P3l8KWuwdqwDso7hxq6Dso7h4yuBf9vApKkKPldlJKkKBlwkqQoGXCSpCgZcJKkKBlwkqQoGXCSpCgZcJKkKBlwkqQoGXCSpCj9P0UgmUFr4QZfAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "text": [
      "       mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_mean  ess_sd  \\\ncoin  0.265  0.441     0.0      1.0      0.006    0.005    4707.0  4707.0   \n\n      ess_bulk  ess_tail  r_hat  \ncoin    4707.0    4000.0    1.0  \n"
     ],
     "output_type": "stream"
    }
   ],
   "source": [
    "import pymc3 as pm\n",
    "import theano\n",
    "import numpy as np\n",
    "import arviz as az\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# These are the initial probabilities\n",
    "p_1 = 0.46\n",
    "p_2 = 0.52\n",
    "\n",
    "# This is the data\n",
    "number_tosses = 50\n",
    "heads = 20\n",
    "\n",
    "with pm.Model():\n",
    "    # Here we indicate that we have 2 coins, and the probability\n",
    "    # of selecting each of the coins is 0.5\n",
    "    coin = pm.Categorical('coin', p=[0.5, 0.5])\n",
    "    \n",
    "    # This is just a way to tell that if we select the first coin,\n",
    "    # the probability of getting heads is p_1 and if we select\n",
    "    # the second coin, the probability of heads is p_2\n",
    "    p = theano.shared(np.array([p_1, p_2]))[coin]\n",
    "    \n",
    "    # Here the likelihood is binomial, because we don't have the result\n",
    "    # for each individual toss\n",
    "    y = pm.Binomial('y', p=p, n=number_tosses, observed=[heads])\n",
    "    \n",
    "    trace = pm.sample(1000, progressbar=True)\n",
    "    \n",
    "    # The mode o\n",
    "    az.plot_posterior(trace, point_estimate='mode')\n",
    "\n",
    "    plt.show()\n",
    "    \n",
    "    print(az.summary(trace))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Signal detection\n",
    "\n",
    "We have 3 sources (labeled 0, 1 and 2) and one detector; when the detector is \"paired\" with one source it only reads the messages from that source. Let's assume that our detector is paired but we don't know what source is paired with. When the detector reads the message it sees some white noise that is modeled as a standard normal distribution.\n",
    "\n",
    "$$X = \\mu + N$$\n",
    "\n",
    "In order to find out with what source the detector is paired, we send a different message (a number) from each source and we repeat that message 10 times. Then, we collect the message that the receptor received.\n",
    "\n",
    "We know that the probability that the detector is paired with the source number 2 is twice  the probability that it is paired with source 0 or 1, and we also know that the probability that it is paired with source 0 is the same as the probability that it is paired with source 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "CategoricalGibbsMetropolis: [source]\n",
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 18 seconds.\n"
     ],
     "output_type": "stream"
    },
    {
     "data": {
      "text/plain": "<IPython.core.display.HTML object>",
      "text/html": "\n    <div>\n        <style>\n            /* Turns off some styling */\n            progress {\n                /* gets rid of default border in Firefox and Opera. */\n                border: none;\n                /* Needs to be in here for Safari polyfill so background images work as expected. */\n                background-size: auto;\n            }\n            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n                background: #F44336;\n            }\n        </style>\n      <progress value='0' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      \n    </div>\n    "
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAakUlEQVR4nO3deXhVxcHH8e8YgbAGEEE2BQWt4gIuFFxxodqqxaoVqVKorXXXopUiLvAWl6IsUlsU0Io7uFVQ67607gsWrEtBKCgiIIhA2MJ23j9uiAkJIcRcbjL5fp7nPsmdM+ecmcnAL3PuuTchSRIkSYrNDplugCRJ6WDASZKiZMBJkqJkwEmSomTASZKiZMBJkqJkwEmSomTASZKiZMBJlVAIoVam2yBVdQacVIIQwp4hhL+HEL4OIawJIXwRQngkhLBj/va98rcvDSGsDiG8HUI4YbNjjA8hzCnh2K+GEF4t9LxbCCEJIZwaQhgXQlgELCy0/dwQwgf55/k2hPDPEMKhhbbXCSEMDSHMDiGszf96dQjBf9+q1nbMdAOkSuopYClwAbAYaAn8BNghhNACeB3IBS4GlgEXAU+HEE5KkuSZcp7zNuAZoDeQDRBCGAZcAdwFDAI2Al2AXYE38wP3OWAfYAjwn/zt1wKN8/eVqiUDTtpMCKEJ0B7okSTJ5EKbHszffjnQCOiaJMnM/LJ/AJ8AN5AKqfJ4N0mS3xRqRzugHzAySZLLC9V7utD3vYDDgaOSJPlXftlLIQSAQSGEoUmSfF3O9khVmpcwpOK+Af4H/Cn/8mD7zbYfCby9KdwAkiTZADwEdAwhNCjnef++2fPjSP0bHVvKPicAn5O/mtv0AJ4HapBazUnVkgEnbSZJ/YmN7sD7wE3AjBDC/0IIF+RXaQzML2HXBUAgtborj82PuVP+1y9L2acpsBuwbrPHu5sdQ6p2vEQplSBJkv8Bvwypa30HkHqtbXT+TSNLgF1K2G0XIMnfDrAGqFlCvZ1IrRKLnXaz54vzv7YEpm+hqd8As4EztrB9zhbKpei5gpNKkaRMBTa9BrYv8E+gSwihzaZ6IYQsoCfw7yRJcvOLPwea5b+mt6neHsBeZTz9i6RuKvltKXWeBVoDK5Ikeb+Ex+JS9pWi5gpO2kwIYX9gFDARmAlkAX2B9cDLpC4l9gVeCCEMApYDFwJ7AicWOtQjpO5sfCCEMAJoAlzFdyuzUiVJMiuEMBK4PIRQH5gMbAA6A/9NkmQi8ADwK1I3lgwHppFaNe4B/BQ4JUmSVeUaCKmKM+Ck4hYAX5BatbUidanxP8BJSZJMAQghHA4MBW4HagFTgROTJHl200GSJJkZQjgduB54ApiRf8yBZW1IkiS/DyHMJBWgfYCVwIekbiIhSZJ1IYTjgQGkVnpt8+vMInW35dryDIAUg5B6PV2SpLj4GpwkKUoGnCQpSgacJClKBpwkKUoGnCQpStvjbQLepilJ2hahIg7iCk6SFCUDTpIUJQNOkhQlA06SFCUDTpIUJQNOkhQlA06SFCUDTpIUJQNOkhQlA06SFKW0f1TXpKnz0n0KVTE9OrbMdBMkVQOu4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLgJElRMuAkSVEy4CRJUTLg0ui9f73IKZ1asfCruWk/17q1efxt+B/pc8wB9OzaniGX/HK7nFeSKisDLhLjbr6OV558mL79rqX/LWNYvnQJg8/vxdq8NZlumiRlhAEXgcULv+LFJyZwzhWDOfrk0zno8GMYMHwci+bP459PP57p5klSRuyY6QZsT6Ou68cXM6fT64IrGD/yer6eP5f9Dj6U310/ihXLlvLXIf357KOptGrbnksGD6PNnvsAkLd6Nff++UbeeOEpVuYuZ7d2P+Csi/vTqetRBcdOkoQJY0bw7MP3kpe3hi5Hn0Cnw7oVa8PavDU8ePswXnt2EsuWfEOrNntw9iUDOPiIY8vdr6lv/QuALsf+uKBsp6bN2bvTIUx54xW6n/qLch87HSZNnZfpJqgS6dGxZaaboEhVuxXc4gXzePD2YZx10ZVceM1Q/jttCqOH/IFhAy7iiON70P+WMWzYsJ5hAy4iSRIA/jrkSl6a/DCn//oSBgwfR5NdmnP9pX345N/vFhz3qYfu4uGxt/Kj087iD7eMoWZ2NvfeekOx89985Xm8MvkRTj/nEq4edTftOhzAjf3O4X/TPy6os3HjRjasX1/6Y8OGgvpfzpnJTk2bU7tO3SLnatW2PfPmzKroIZSkKqFareAAcpcv5U/3TKJ56zYAzPnsU5645w4u++OtHH3y6UBqNXb9pX34cvZMAF57dhKXDB7BMT/9OQCdDu3G787ozsPjRjF49ANs2LCBx+++neNPO5uzLupfUGfQ+b345usFBeee9s7rvP/aS1w/7hH2Pbhrql7Xo/jq89k8euef6X/LGABuG3wFrzz5SKn96HBQF26481EAVi5fRt36DYrVqdcghxXLl5V3qCSpSqt2Ade0RauCcAMKvt+v86Hfle3aFoAlXy9gyaKFJEnCYd1PKti+ww47cGj3E/n7+NuB1Gtg3y5eSOduPypyri7H/Jhp77xW8PzDd16jUZOm7N3xEDasX19Qvv8PD+Plyd8F2pnnX86JPfuW2o/adesVeR5CKFYnSRJKKJakaqHaBVzdejlFnu9Yo2aqvH5OobIaAKxdm8e3ixeSXacutWrXLrJfw8Y7k7dmNevW5rF08dcA5DRuUqROTuOdijxfvnQJ3y7+mtMOaVOsXTtkZRV8v/MuLWnStHnpHSmUXHUb5LAyd3mxKitzlxfplyRVJ9Uu4LZVoybNWLNqJXmrVxcJuaVLFlEruzY1ataiYZOmACxbsrjIvsuWfFPkeb2chuzUdBeuGnFXqefc1kuUrdq0Y/HCr1izehXZtesU1Plyzkxattlj652UpAgZcFvRrsMBhBB488Wni7xG9+aLT7N3p84ANGnWgkZNmvLuq89z4GFHF+z79svPFDnW/p0PZ9J9Y8muU5dWbdtt8ZzbeomyY9cjC87X7cTTgNTl1U8/eJfzrip+o4skVQcG3Fa03r09R5zQg7FDr2HVylyat27D848/yLw5szh/4E0AZGVl8bM+FzB+5BDqN2zMPgd25q2X/sGXsz8rcqyOXY6kU9ejGHRBL07teyG77rEnq1asYPaMj1mXl0fvS68CoFmL1jRr0brMbWzSrAXHnXImdw0bTJIk5DTaiQljRrBz85YcdeKpFTcYklSFGHBlcNG1t3DPqBt4eNyogvfBXTNqPPvkr+AATj7rN+QuW8pzj97Hkw/eSeejuvPL313NyIGXFNQJITBg+Dgeves2nnzwThbP/4p6OQ1pu9c+nHjmOd+rjef2/yPZ2XW4e/gfyVuzmg4HdeGKm/5KzVrZ3+u4klRVhU3v9UqXSVPnpfcEkqo03+itElTI/d/V7o3ekqTqwYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYBTtff2K89y2RnHcXrn3fntiV2ZdN/YUuvfecsgTunUirtHDClS/uXsmVzZ+yR+ccTeDPvDhaxetbLI9o+nvM05PzqoWHlJXpr8MKd0alVi3YfuGE7vo/creL7wq7mc0qlVwaNn1/acd/JhjBh4MR9/8E6x/Udd148rfvGTrbZBquoMOFVrn059j6FXnEv7Dh25etTdHNujJ/f++UYmP3BnifXnzprBS5MmUqde/WLb/jyoH81bt+HKobczd/YMHr3rtoJtGzdu5K5hg+l9yVXUrlM3LX3p2+9aht4ziWtvu5czzr2M3GXfcvWvT2PCmBFpOZ9U2e2Y6QZImTRx7K3s3bEzFw8aBkCnrkexcvkyHh47kh+f8Utq1KhZpP64m6/jpF6/5tWnHytSvnrVSmb8599cfet4chrvxMrc5Txx7xh6XzIAgBefeIisHXek20mnpa0vLdvswV77HwTAvnTl2J+ewYOjb2HCHSPocFAX9jv40LSdW6qMXMGpWps9/WMO+OHhRco6dj2SFcuXMX3alCLlb77wFF/O/ozTzrmo2HHWr1sLQM3s7PyvtQvKVq3I5cHRw/jNlf9HCCEd3diinuddTuOdm/HcI/dv1/NKlYEBp2pt3do8dqxRo0hZjZqpVduXs2cWlOWtWc3dI4bQ+9KryK5dp9hx6uc0olnLXXn6obvJXfYtzz/+AHvssz8AD4+7lQN+eHjB6mpbbNywgQ3r1xd5JBs3lnn/rKws9ut8GNP/88E2n1uq6rxEqWptl9Zt+OzjaUXKPvtoKgC5y5cWlD32t7/QaOemdDtxy5cYfzvgem7pfz73/+VPtNi1LecNuIH5X8zmxScmcOvDL5SrfWcduU+J5fUbNirzMZo0bc6yJYvKdX6pKjPgVK2dcPrZ3HHjQJ5//AEOPe5EPvtoasFdlFk7pC5wLJz3BU/cN4YhYyaWeonxoMOP4Z6XprH46/ns0mo3srKyuOGyX3HyWefSpFkL/jFxPI+PHw3Aqb+6iJ+c0Wer7bvhrseoVSu7SNlzjz/A2y8/U+Y+JklS5rpSTAw4VWvH9jiTOTM+4Y4bBzJ6yB+olV2bX142kHFDryVnp50BuPfPN3HgoUfTsm07VuQuAyBJNrJubR4rcpdRt16DguCrVbs2LXfbHYBpb7/GnBmf8Puho5k9/RMeHH0Lfxo/CYA/9OnBPh0Poc2eJa/QNtn9B/sWu+vyvdde3KY+frNoATmNd96mfaQYGHCq1rKysvjtgBv4xYVXsnjhfJq13JV5+a+97bXfgQDMmzOLOTM+KbZq+sfE8fxj4njufPZdmjRrUWTbhvXruWvYYPr87mpqZdfmoylvst8hh9GqbTsA9u98GB9NeXurAfd9bVi/nv+8+wYdDuyS1vNIlZEBJwH1GjSkXoOGADzzyL384ICDC8Lo4utuYfXqom+4Hj7gIjoc1IUTft6bnEY7FTves4/eR90GORx+/E8LyvLWrP7u+9WrYDtcOpw4diRLFi3k+NPPTvu5pMrGgFO1Nv3DKXw69T3a7tmBVStzee3ZSfz7rX9y098eL6jTrsMBxfarUbMWTZq1KPG9ZSuWL2XimJEMGv1AQVmHA7twz6038uITEwD48L036X3pwArty7w5s2jQsBHr161j4by5vP7cJD5481XOPP9y9j24a4WeS6oKDDhVa1k71uD1555kwh0jCDvswD6dOnPT3X+nTfu9y33Mh24fziFHdWePvb/7OK3df7AvfX43kAf+ejMAfftdQ9u9Kvby5PiRqY8Oq1mrFo2aNGOv/Q/khrseo8OBP6zQ80hVRUj3HVaTps7zFi5JW9SjY8tMN0GVT4V8IoJv9JYkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEXJgJMkRcmAkyRFyYCTJEVpx0w3QFL1NmnqvEw3QZVMj44tK+Q4BlyazZ01g7E3X8v0D6dQt14O3X92Jj3Pu5ysrKxMN01Shr3xwlO8+tRjzPr0Q1atyKVFmz04pfd5HPnjUzLdtCgYcGm0YvlSrju/F613b8/AkX9jwdzPuXvEH0mShLMu6p/p5knKsMn3jaVpy9ac8/vBNGjYmCmvv8yIgRezfOkSTup1TqabV+UZcGn07CP3szZvDQOGj6NOvfrQBVatzGXCmBH8rM8FqTJJ1dbVo8bToFHjguf7dz6MJYsWMPn+cQZcBfAmkzT64I1X6NT1qCJBdsTxPVi7Zg0fTXk7gy2TVBkUDrdNdv/BvixbsjgDrYmPAZdG8+bMpGXbPYqU7dy8JbWyazNvzswMtUpSZfbfae/Tevc9M92MKBhwabQidxl16+cUK6/XIIcVy5dloEWSKrNp77zOu68+z4979sl0U6JgwKVZCMXLkiQhlLRBUrW18Ku5jBh4MZ27/Yhjf3pGppsTBQMujerVz2Fl7vJi5atW5FK3foMMtEhSZZS77FuGXNybnXdpSb8bbst0c6LhXZRp1LJNO+bNnlWkbNGCr1izehUt27TLUKskVSZ5q1dz/aV9WbduHdffeQ/ZtetkuknRcAWXRgcedjT/futVVq9cUVD2xvOTqZmdzb4HdclgyyRVBhvWr+fm/ucx/4vZXPeX+2jYuEmmmxQVAy6NTvj52dSoWYs/XXEu095+jeceu58Jd4ygx9m/9T1wkrjjpoFMef1lzjj3MlYs+5bpH04peKxbm5fp5lV5IUmStJ5g0tR56T1BJTd31gzGDr0m9VFd9XM47pRenHm+H9UlCc79SRcWzf+yxG1jnn6LZi1ab+cWVQ49OraskLvwDDhJUqVSUQHnJUpJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKUDDhJUpQMOElSlAw4SVKU/GsC29FDdwxn4piRmW6GpEqq53n96HX+FZluRjRcwUmSomTASZKi5IctS5IqFT9sWZKkUhhwkqQopf0uyqmT7pwFLE33eaqIFsBXmW5EJeA4pDgOKY7DdxwLYOokZg4aNKjd9z5QkiRpfQwePDhJ9zmqysOxcBwcB8fBsdh+4+AlSklSlLZHwP3fdjhHVeFYpDgOKY5DiuPwHccipULGIe1vE5AkKRO8RClJipIBJ0mKkgEnSYpShQdcCKFWCOG2EMLiEMLKEMLkEEKrrezTN4SQlPDIruj2pUsI4cIQwuwQwpoQwpQQwhFbqb9fCOGfIYTVIYR5IYTrQggV8vE0mbYtYxFCaLOFn/0J27PNFS2EcGT+3J+X35++ZdgnujmxreMQ8Xy4KoTwXghheQhhUQjhyRDCvmXYL6o5UZ5x+D5zIh0ruFuB04BewBFAA+CpEELWVvZbBTQv/EiSZE0a2lfhQgg9gVHAjUAn4E3gmRDCrluo3wB4AVgIHAJcClwJXL5dGpxG2zoWhZxA0Z//y+ls53ZQD/gIuAxYvbXKEc+JbRqHQmKbD92A0cChwDHAeuDFEELjLe0Q6ZzoxjaOQyHbPicq8s15QA6wFjirUFlrYCNwfCn79QVWZPrNhd+j3+8A4zYr+wy4aQv1LwCWA7ULlV0DzCP/ztaq+ijHWLQBEuDgTLc9jWOyAui7lTrRzoltHIfo50N+P+sBG4CTq/mcKMs4lHtOVPQK7iCgBvD8poIkSeYCn5JK7NLUDiF8HkL4MoTwVAihUwW3LS1CCDVJ9fv5zTY9z5b73BV4LUmSwr/RPkfqY3raVHQbt5dyjsUmj4cQvg4hvBFCOD0tDazcopwT30Ps86E+qSto35ZSpzrMibKMwybbPCcqOuB2IZXGizcrX5i/bUumA+cAPUhd2lwDvBFCaF/B7UuHJkAWqT4WVlqfd9lC/U3bqqryjMUK4PfAGcBPgJeAiSGEs9PVyEoq1jmxrarLfBgFTAXeKqVOdZgTZRmHcs+JMn3YcgjheuDqrVQ7urRDkFpilihJkrco1MEQwpukOn0JqevOVcHm/Su1z1uoX1J5VVTmsUiSZDEwvFDR+yGEJkB/4P70NK/SinlOlEl1mA8hhBHA4cDhSZJs2Er1aOdEWcfh+8yJsq7gbgX23srjXWABqd/gm2y2f1OK/yayRfmdfR+oCiu4xaRWrZv/RlVanxdsoT6l7FMVlGcsSvIOVeNnX5FinRMVIZr5EEIYSeoq1TFJkvxvK9WjnRPbOA4lKdOcKFPAJUmyOEmS/27lsQqYAqwDuhfqSCtSAfhmWVuefxvs/sD8su6TKUmSrCXV7+6bberOlvv8FnDEZm+D6E7qz2TMqeg2bi/lHIuSdKQK/OwrWJRzooJ0JIL5EEIYBfyC1H/q/y3DLlHOiXKMQ0k6UpY5kYa7Ym4ndZfPcaRuE3+F1OXGrEJ1XqLQXXXAIOB4YPf8hv+NVFB2zvRdPmXsc09Sd4/+hlSYjyJ13Xi3/O03AS8Vqp9D6rezCcC+wKmk7pa6ItN9ycBY9CE12fcG9iJ1rX0t0C/Tffme41Avfy53JPUWmOvyv9+1Os2JcoxDrPPhr/k/z2NIrco2PeoVqhP9nCjnOJR7TqSjA9nAbcA3+RP6SaD1ZnXmAOMLPR8JfA7kAV+TulOoa6Z/GNvY7wvz+5VHahVzZKFt44E5m9XfD/gXqRtq5pMK+Vhu/S3zWORP3k+AlfkT/33g7Ez3oQLGoBup10k2f4yvTnNiW8ch4vlQ0hgkwOBCdaKfE+UZh+8zJ/xrApKkKPlZlJKkKBlwkqQoGXCSpCgZcJKkKBlwkqQoGXCSpCgZcJKkKBlwkqQoGXCSpCj9P3B+8hhWhE4jAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "text": [
      "         mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_mean  ess_sd  \\\nsource  0.564  0.611     0.0      2.0       0.01    0.007    3946.0  3946.0   \n\n        ess_bulk  ess_tail  r_hat  \nsource    3946.0    3983.0    1.0  \n"
     ],
     "output_type": "stream"
    }
   ],
   "source": [
    "# These are the values of the messages\n",
    "mu_0 = -0.3\n",
    "mu_1 = 0.2\n",
    "mu_2 = 0.7\n",
    "\n",
    "data = [-1.26318642,  0.5758068 , -0.16657075,  0.08515534, -0.89754746,\n",
    "    0.20262975,  0.19462712,  0.52539873, -0.34519777,  0.34857591]\n",
    "\n",
    "with pm.Model() as signal_detection_model:\n",
    "    # Here we indicate that we have 3 different sources, and the prior probabilities\n",
    "    source = pm.Categorical('source', p=[0.25, 0.25, 0.5])\n",
    "\n",
    "    # This is just a way to tell that if the data comes from source 0\n",
    "    # the mean should be mu_0 and the same with the other sources\n",
    "    mu = theano.shared(np.array([mu_0, mu_1, mu_2]))[source]\n",
    "\n",
    "    # The likelihood is normal\n",
    "    y = pm.Normal('y', mu=mu, sigma=1, observed=data)\n",
    "\n",
    "    trace = pm.sample(1000, progressbar=True)\n",
    "\n",
    "    # The mode is the MAP\n",
    "    az.plot_posterior(trace, point_estimate='mode')\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "    print(az.summary(trace))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data and the model indicates that, contrary to our prior beliefs, the detector is most likely paired with source 0.\n",
    "\n",
    "**Q:** What's the probability that the signal is from sources 0 or 1?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true,
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "text": [
      "C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\arviz\\data\\io_pymc3.py:89: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n  FutureWarning,\n"
     ],
     "output_type": "stream"
    },
    {
     "data": {
      "text/plain": "array([<AxesSubplot:title={'center':'source'}>], dtype=object)"
     },
     "metadata": {},
     "output_type": "execute_result",
     "execution_count": 4
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnZklEQVR4nO3dd3gVVeLG8e9JbyS0JJACCYReFWmCioVFRcS2iouuWNYCdsS1i3XVpfkDUVEQd1eRFQv2unZEFEGl14QkECCQkF7v/P644ZqQQsCECyfv53nuk3tnzsycOTnJmzNz7o1xHAcRERHb+Hi7AiIiIo1BASciIlZSwImIiJUUcCIiYiUFnIiIWEkBJyIiVlLAiYiIlRRwIiJiJQWcyFHIGBPo7TqIHOsUcCI1MMZ0Nsa8ZYzZZYwpMsZsM8a8bozxq1jfpWJ9tjGm0Biz1Bhz5gH7mG+MSa5h318aY76s9HqYMcYxxlxgjHnBGLMb2Flp/d+MMT9XHCfLGPOVMebESutDjDFPGmO2GmNKKr7ea4zRz7c0aX7eroDIUeo9IBu4AcgEYoGzAR9jTAzwLZAL3AjsAyYA7xtjznEc58PDPOZM4EPgciAIwBgzBZgIzAUeBFzAIKAdsKQicD8GugOPAL9VrL8faFmxrUiTpIATOYAxpjXQCRjtOM47lVa9WrH+dqAFMNhxnE0Vyz4A1gCP4Q6pw7HMcZxrKtUjCbgNmO44zu2Vyr1f6fmlwFDgFMdxvq5Y9rkxBuBBY8yTjuPsOsz6iBzTdAlDpLo9wBbgiYrLg50OWH8ysHR/uAE4jlMOLAD6GmPCD/O4bx3w+gzcP6Nz6tjmTCCFitHc/gfwCeCPezQn0iQp4EQO4Lj/xcZw4CfgH8AGY8wWY8wNFUVaAjtq2DQDMLhHd4fjwH22qviaVsc2UUB7oPSAx7ID9iHS5OgSpUgNHMfZAvzVuK/19cF9r212xaSRvUCbGjZrAzgV6wGKgIAayrXCPUqsdtgDXmdWfI0F1tdS1T3AVuDiWtYn17JcxHoawYnUwXFbCey/B9YT+AoYZIxJ2F/OGOMLXAKscBwnt2JxChBdcU9vf7mOQJd6Hv4z3JNKrq2jzEdAPJDnOM5PNTwy69hWxGoawYkcwBjTG3gaWAhsAnyBcUAZ8D/clxLHAZ8aYx4EcoDxQGdgZKVdvY57ZuMrxphpQGvgbn4fmdXJcZzNxpjpwO3GmGbAO0A5MABY5zjOQuAV4ErcE0umAr/gHjV2BM4FznMcp+CwGkLkGKeAE6kuA9iGe9QWh/tS42/AOY7jLAcwxgwFngSeBQKBlcBIx3E+2r8Tx3E2GWMuAh4F3gY2VOzznvpWxHGcO4wxm3AH6BVAPvAr7kkkOI5TaowZAdyFe6SXWFFmM+7ZliWH0wAiNjDu++kiIiJ20T04ERGxkgJORESspIATERErKeBERMRKCjgREbHSkXibgKZpiojIoTANsRON4ERExEoKOBERsZICTkRErKSAExERKyngRETESgo4ERGxkgJORESspIATERErKeBERMRKCjgREbFSo39U1+KV6Y19CDnGjO4b6+0qiEgToBGciIhYSQEnIiJWUsCJiIiVFHAiImIlBZyIiFhJASciIlZSwImIiJUUcCIiYiUFnIiIWEkBJyIiVlLAiYiIlRRwIiJiJQWciIhYSQEnIiJWUsCJiIiVFHAiImIlBZyIiFhJASciIlZSwImIiJUUcCIiYiUFnIiIWEkBJyIiVlLAiYiIlRRwIiJiJQWciIhYSQEnx7zi4mImTpxIVFQUoaGhjBw5kuTk5INuN2zYMIwx1R5FRUXVys6ZM4eePXsSFBREdHQ0l1xySSOciYg0JD9vV0Dkj7r55ptZtGgR06dPJzIyksmTJzN8+HB+++03goKC6tz21FNP5fHHH6+yLDAwsMrr++67j1mzZnHffffRv39/du7cyVdffdXg5yEiDcs4jtOoB1i8Mr1xDyDHnNF9YxtsX2lpaSQkJDBv3jz++te/ApCenk5iYiKzZ8/mmmuuqXXbYcOG0bp1axYtWlRrmdWrV9O7d28++ugjhg8f3mD1FpE6mYbYiUZwTcTTD9zGtk3rufSGicyf/ii7dqTS64QTufXRp8nbl80zj9zJxlUriUvsxE2Tp5DQuTsALpeLt+bP5tO3F5CZsYOotrFcdPXNnHbunz37/umbz3n3lRdJ3rCGkpJi4jt04tIb7uC4wad4yix4biofLJzPQ88uoMtfR5G8cS2x7TtyzZ0P0+P4gYd9Xp+9vRAHCEgaxOKV6Z7lXfv2Z+6CN4k84axat83MK6bEr7DKdgeaP2MWbeLaUxDZvc5ycvga8g8ekcp0D64JycxI59VnpzB2wiTG3/ck635ZzuxH/s6UuyZw0ojR3PnP5ykvL2PKXRPYP7J/4cn7ef3F/+NPF4zl/v97mYGnncmshyby49efefa7M30b/U8+g1sffZq/T5lD194n8MiNl7N25Y9Vjl9cVMjTD9zKiAsv4+//nIN/QABPTLyG4sJCTxmXy0V5WVndj/JyT/m05E20impLcEholWPFJXYiPXnzQdtk5dKvuXhwEhcPTmLy+LEkb1hTZf3G31bQLqkrC56byuWn9uKiAR144LoxpG7ZWP+GFxGv0AiuCcnNyeaJlxfTNj4BgOSNa3n75ee45eEZnDrqIgAcx+HRm68gbesm/Pz8+Oj1f3HT5GmeEVufQSeRtXsXC5+fTv+TzwBg5JgrPcdwuVz0OuFEtm3ZwGdvv0a3vv0960qKirj6jofoPWAIAC0io7l9zAhW/7yU44ecCsDMyRP54t3X6zyPHv0G8diL7suK+Tn7CG0WXq1MWHgEeTn76txPz36DOG3Un2kTn8DuHWksmjuTe66+kOkLPyE6Jh6ArD272bxuFWlbNzHh/qfw8fXj1Wee4uEJl/HM218REFj3PT4R8R4FXBMSFRPnCTfA87zXgBN/X9YuEYC9uzLISEvB+Pgw6LQzKS8r85TpPWAo33y8mPLycnx9fcncuZ1XZj3FLz98Q1bmLs/or3K4Afj5+dPzhMGe1/EdOgGwZ+cOz7Ix19/OyEvG1XkewaFhVV4bU/1yveM41LC4iktvuKPSq4H0GXgSE84fxruvvMg1kx5y78floriwgL9PmeOpb7uOnRl/3sl89cFbDD//0roPIiJeo4BrQkLDIqq89vMPcC9vFlFpmT8AJSXF5GTvxVVezl9O6lbj/rIyd9Iysg2P33oVhfl5XHrDHbSNTyAoOIRXn53Cvr17qpQPDgvDx+f3q+L+FccvKSn2LItsE0vrqLZ1n0il5AoNjyA/N6dakfzcnCrnVR8tWkfRre8JbFm7yrMsLKI5zQsLPOEG0CauPVExcbpMKXKUU8BJrcLCm+Pr58c/XnoLH1P9dm1Ey9bsSE1my7pVPDDr357LjAAlxdXfS1Yfh3qJMi4hicyd2ykqLCAoOMRTJi15E7EJHQ+rDpVHfnGJncjM2F6tjOOAj0+DTPQSkUaigJNa9R4wBFd5OQV5ufQddHKNZUoq3hTtHxDgWbZrexrrVv5E+041j/zqcqiXKPsOdtdr6f8+ZNjICwH35dW1Py/jursfO6RjZ+/ZzdqVP3HG6N/fxN3/pNP54t3XSd28gfiOnQHYkZrM7h1pnpmmInJ0UsBJrWITOjLiosuZetd4zr/iBjp270NpSTHbNq9ne8oWbnxwCnGJHWkV3ZaXpj3CX8ZPojA/jwXPTaVlVJvDOmZ0TLxngkd9tI6O4YzzxjB3ymQcxyGiRStee34akW1jOWXkBZ5yC5+fzsIXZvDmTykAJG9Yw79nPsGJw88hqm0su3ds542XZuHjYxg19vf3zg087Sw6dOvFE3f8jbHjJ+Hj68urz04hpl0Hhv7p3MM6RxE5MhRwUqfr7n6MmPYd+PTNV3n12amEhIYR36ETZ5znnlzhHxDIXVNfYM4/7uWpSdfSKqotF11zM6t++p5tm9YfkTr+7c6HCQoK4aWpD1NcVEiPfoOY+I9nqsxwdDkuXJXeXtCseQscx+E/M58gNzuL4NBQevQbzN3T5hLZ9vf3Zfn6+vLAzH8x95+TmfXwJFwuF30GDuWaSQ977leKyNFJn2QiIl6lN3pLDRrkBrfe6C0iIlZSwImIiJUUcCIiYiUFnIiIWEkBJyIiVlLAiYiIlRRwIiJiJQWciIhYSQEnIiJWUsCJiIiVFHAiIkeDlQvcD2kw+rBlEZGjwb5Ub9fAOhrBiYiIlRRwIiJiJQWciIhYSQEnIiJWUsCJiIiVFHAiImIlBZyIiFhJASciIlZSwImIiJUUcCIiYiUFnIiIWEkBJyIiVlLAiYiIlRRwIiJiJQWciIhYSQEnIiJWUsCJiIiVFHAiImIlBZyIiFhJASciIlZSwImIiJUUcCIiYiUFnIiIWEkBJyIiVvLzdgUEojK+ovPGF4jIXgvGkNW8F+u63cjeVv2qlOuz8kFa7l1JUOFO/MoKKPVvRnbzHmzsdA17IgfUeYz4lLc4fsW9ta4vCI7h0xGfAdBizwp6/fYEzXI3URASy5rut7Gz7ameskEFOzj985GkxY3il+Me+gNnXrvoHV8Ql/Y+zbNXEVS0G5dPALnNOrC54xXsiP1ToxwTIKgwgy7rZhO961sCizIp9Q9nX0RX1na/hewWverctuvamURnfEVIQTp+ZXmU+YWyL6IrWxP/UqXOR0P7ijQFGsF5Wfy2xQxaOp5We37Gr7wQv7ICIjN/YMi3V9J69/fVyobnbCSgNAcfp4zAkiyid33LiUuupsWeFX+oHmV+IQD4leQw8IcJBBZn8v2Jc3D5BtL/x1sJyd/mKdt9zTQc48fa7jf/oWPWJXHrAuLSPyAsfxt+5YUElO6j1d4VDPjxVpI2vNgoxwzL3cywLy4iIWURwYUZFW28l6jdSwjft/6g28ekf0TzfWsIKN2Hj1NOQGkOkZnLGPDjrcSkfQgcPe0r0hRoBOdFxlVGj1VPYXDID4nnuyHzMLgY8u2VhBRup8/Kh/n8jA/AGAA2dL6WjLankR/aDv/SXHr/8jBtM77AxyknNv1DslodV+uxUtufT2r786ssi874kkFLxwOQFjcSgJZZKwksyWZzx3PZ26ofqXGj6LXqCaJ2LSE5sR3Ns34jLu0D1vS4nZLAVg3WFs2zVuFbXsie1v0BcPkGsjHpKra1v5DC4GjiU9+hzy8PA9B5wxw2J43D8alf9w3L3UJofio725xSZ7njl99NYMleSvzD+aXvZHZFDcHHVUrLvSsoDmx90OOkJPyZvS2PI7dZR3xcpXRbO4OE5NcBiEt7j+1xZ3mtfUWaIgWcFzXL2UhgSRYA22OGUxga63metPllwvJTaJ79G9ktegOwoet4z7blfiGktL+IthlfANT7l31lCVsXAuAyfmxrfyEAPq5S9zIff/dX34Aqy3uuepKCkDi2dLj8kI93IP+SbOJS36N9yhtE5KxnXZfxnoBb3u8Jyv1CPWWTE8eQuHUB4Tkb8S/LI6Aki+KgyFr37VtWQGz6R7RLeZNWe39mW/x5dQZciz0raJG9CoA1PW5ne+yZnnUZbU+v1/lsThpX5fWWxLGegHMq2vNItq9IU6eA8yJfV/HvLypGaQeKyF7rCTgPxyGoMIP2KYsAKPMNJjV+9CEdO6hgB9E7vwFgR9vTKQ5yj1Cym/ei3CeQqJ3fsqHztURnfImDYU+rfsSkf0SrPT+zbMDTnl/Mh8xxaJ25lPYpb9J2+6f4ukpwMOxuPYDdUYM9xSqH234+5e72KvcJpCQgosbdN8/6lfbJbxCb/gH+ZfkAZDXvxY6YM+qsVus9P3qeN8vZzOmfjiC4cCf5oe3Y2Oka0tqde0inGViUSYctrwDgMr4kt78IOALtKyIeCjgvygtLxGX88HHKiEn/hK2JYzCOQ8z2Tz1lAkqyq2zTY9VTJG2a73ld6hfGsoFPkxPR5ZCOnZD8OgYX4B4d7VcUHMXK4x6m9y8PM/L9gZT7+LO6x0RywpPov+w2Mlv1Z0fMcACMq9QzMjmYoMIM2m17m3YpbxJakAZAblgiafGjSI0/l8KQmDq3j0t9l7CK+1Sp8aNwfH4PAP+SbOJT36V9yiLCczYCUBDchq2Jl5La7jzymnU4aP2CCzM8zztu+bfneXjuJvr9fBc+rhK2JVx00P0kbnmF3r8+5nld7uPPyuMeZXf0UKDx2ldEqlPAeVFpQARbOowlafPLhBak8qdPhlcrc7BLj/5leQz44Ra+G/oS+5r3qNdxjauM9tveANwhkxk5sMr6tPhRpMeeRXBhBkVBUbh8A0jaOJeQgnR+HDCdiOw19F35IBHZaynzCyU54c+s6XE7mNrnLA395nJCC9Ip8Q9na+IYUuNHk9WyT73qG7lrCX1XPABATrMkVvecVGV9/2W3Epm5rGIkey7b2o0ms/XAOutTU5vsl9W8J0sHP0tQ4S6GfnsF/mV5dF03q14BdyBfVyl9V9xHqX+45xJpY7SviFSnnxgvW91zEmu630p+SCzlPgHkhHciOeHPnvWFwW0PKH8ni0ev4uMRX7Cp4xWAO+S6rn2m3sdsk/EFQUW7AUhOuLjGMo6PHwWhcbh8AwgozqLz+jlsa38+OeFJDPjhJsJyt/Bzv3+wO2ownTbNo33FvabaOBW/nH1cpfiV5uJXmgtO+UHrGrnrOwYunYCvq5jcsA4sGTKXMv9mVQsdsG//0jyMU1bD3mpXGtDc8zwtfhQlga3Iad6N3ZGDAAgu2kVAcdZB97O1w1gWj17Nh2d9y+oeEwF3yHVfM61KuYZuXxGpTiM4bzM+bOx8LRs7X+tZ1H21+5ehgw97Wh1f4zZFwdGs7zqepM0vAxCan1LvQ+6fXFLmG0Rqu/MOWr7rulmAi7XdbiEsN5mQwh3saHMaafGjyG3WkZjtnxC1awkpiZfUuo+vT/kvcWnv0S7lDeLT3ic+7X0Kg6JIizuH1HajyQ3vVG0bd7jdiK+rmH3hnVkyZG6NMwt/GDiLmPSPab/tDdpmfEHbjC8o8Y8gPe6seo8UsyO6eZ471Hw/tLy+98WMoSSwJZs6XU2nDS8QUJpDaF7t35+GaF8RqU4B52WRu5bg8vFnX0Q3jFNGTPrHdKgIre0xZ1AU3AaANjs+JyQ/jV3RJ1EQEot/aS5Jm+Z59lMQGud53mXtLLqunw3AJ8M/9czOBAjJTyWy4v1122POpLSWyRr7heVuoX3y66zrdhPFQa099wRdFZdO988GdHx869xPaUAEWzuMZWuHsURkr6F9yhvEpr1Pp03z6LRpHtkR3Vjd8w4yI90TTSJ3fsvAH27C11VMdkR3lgx5scooq7JyvxDP2yBC85Jpn/IG8dsWk7j1NRK3vkZeWALrut5IetzZtdZvV/TJlPqF4V+WR3zqO6THnUVQ0W4idy8FILt5D8/El5rat8WeFUTuXlrxNo54fFwltNv2NgGlOUDV709jtK+IVKeA87JWmT/SZcPz1Zbnhbbjt96/f/JISH46vVY9CauerFa23CeQ9Z2vr9fxEpL/i8EBILkeI4Ieq56iMLgNmysuh+aFJZAXlkDk7qW03v0DsekfAJARPaxexwfY17w7vzbvzqqedxKz/RPapbxJ68xltMpc7gm4zhte8Mwybb5vDWd/cGKVfXw7ZH6Nn96SH5bAmh4TWdvtFqJ3fkW7lDeJ3vkNUTu/rTPgyvxDWdXz7xy38n5aZK/irA9P8qxzGT9W9ZhU67YAgcV76LZuJt3Wzay2zsGHtV1vqnG7xmhfEXFTwHlZVove7G3Rh7C8rfiWF1IUFM2Otqezoct1VUYsWS37sL3tcJpnryagZC8+rjKKgqLY0/oENiZdRW5E54Mey7hKaJfyFgD7wrsc9NJd693f02bn1yzrP90zbd3x8eOHgTPp/etjDPxhAiX+EazpdsshT6MH95u50+JHkRY/ipD81GozRv8Ix8ePjLank9H2dAKLdhOWl3zQbbYlXEipfzM6bZxLeM4GXD5+ZLXsy7quE8hq2bfObXPDk0iNG0mLrF8JKsrEx1VKcWArslr2YXPHy6t97Bo0fvuKNHXGcZxGPcDilemNewAROaaN7ht78EJNwVdPub+ecqd363F0qPlG+CHSLEoREbGSAk5ERKykgBMRESsp4ERExEoKOBERsZICTkRErKSAExERKyngRETESgo4ERGxkgJORESspIATERErKeBERMRKCjgREbGSAk5ERKykgBMRESsp4ERExEoKOBERsZICTkRErKSAExERKyngRETESgo4ERGxkgJORESspIATERErKeBERMRKCjgREbGSAk5ERKykgBMRESsp4ERExEoKOBERsZICTkRErKSAExERKyngRETESgo4ERGxkgJORESspIATERErKeBERMRKCjgREbGSAk5ERKykgBMRESsp4ERExEoKOBERsZICTkRErKSAExERKyngRETESgo4ERGxkgJORESspIATERErKeBERMRKCjgREbGSAk5ERKykgBMRESsp4ERExEoKOBERsZICTkRErKSAExERKyngRETESgo4ERGxkgJORESspIATERErKeBERMRKCrgm6puPFnP7pWcy5sTOXPWnfsy47xb27sqoUsZxHF6fO5Orz+zPxYM6cs9VF7Jl/eoqZdK2bmLS5efwl5O6MeXv4yksyK+yfvXypVz1p37Vltfk83f+y3nHxdVYdsFzU7n81F6e1zu3p3LecXGexyWDO3HdqCFMu+dGVv/8Q7Xtn37gNib+5eyD1kFE7KGAa4KWffkJU++eQNc+/bh7+jz+ess9rP75Bx69ZRwul8tT7o15z/D6CzO4YNx47p0xn6CQEB68fgxZmbs8Zf7vwdtoG5/ApCefJXXrBhbNnelZ53K5mDtlMpffdDfBIaGNci7jbrufJ19ezP0z/8XFf7uF3H1Z3Hv1hbz2/LRGOZ6IHDv8vF0BOfK+/vBtOnTrxbV3PeZZFhLajMdvu4r05M3Ed+hESXERb85/hguvvJGRY64EoEufflx79iA+WDifsRPupLAgnw2/reDeGfOJaNmK/Nwc3v7X81x+010AfPb2Anz9/Bh2zoWNdi6xCR3p0rsfAD0ZzOnnXsyrs//Ja89No0e/QfQ64cRGO7aIHN00gmuCyspKCQ1rVmVZaLNw9xPHAWDdL8spyMtlyJ9GecoEBYfQ/5ThLP/uC/d+SksACAgKqvga7FlWkJfLq7OncM2khzDGNOr5HOiS626nZWQ0H7/+nyN6XBE5uijgmqAzzhvDmhXL+OLdRRTk5ZKesoVXnvknvfqfSHzHzgCkJW/Cx9eXtu0Sq2wbl5hE+tZNADSLaEF0bDveX/ASufuy+OTNV+jYvTcA/31hBn0GDvWMrg6Fq7yc8rKyKg+n0qXTg/H19aXXgCGs/+3nQz62iNhDlyiboBNOOp2bH5rGrIcm8fQDtwLQtc8J3DN1nqdMfs4+goJD8fX1rbJtWLMIiosKKS0twd8/gGvvepR/3nk9/5n1BDHtErnursfYsW0rn739GjP+++lh1W/syd1rXN6seYt676N1VFv27d19WMcXETso4Jqg3378jmcfu5tz/nIVxw85lX17Mnnt+Wk8cfs1PPTca55Qq+nSooP7EqbBva7f0NN4+fNfyNy1gzZx7fH19eWxW65k1Ni/0To6hg8WzufN+bMBuODKCZx98RUHrd9jc98gMDCoyrKP33yFpf/7sN7n6FRcahWRpksB1wS9NO0RBpwynCtuudezLLFLDyacfwrLvvyYwaefTWh4BIUFeZSXl1cZxeXn5hAYFIyfv79nWWBwMLHtOwDwy9JvSN6whjuenM3W9Wt4dfY/eWL+YgD+fsVouvftT0Lnmkdo+3Xo2rParMsfv/nskM5xz+4MIlpGHtI2ImIX3YNrgtKSN5HYpUeVZbEJHQkICiIjLQWAuIQkXOXlZKQmV9126yZiE5Nq3G95WRlzp0zmilvvJTAomFXLl9Cr/xDiEpOIS0yi94AhrFq+tFHO6cB6/LbsO7oexv0/EbGHAq4Jimobx5a1q6osS92ykZKiIqJi4gHo2qcfIWHN+O7T9zxligsL+enrz+g35NQa9/vRon8TGh7B0BHn/r5NUWGl7Qs8szQb08I509m7eycjLrqs0Y8lIkcvXaJsgkZcdBnzpjxEi8ho+g05ley9u1k4ZwZRMfH0G3oaAAGBQVwwbgL/fXEGYeERxCYk8c5/5uByXJ73xVWWl5PNwuen8+DsVzzLehw/iJdnPM5nb78GwK8/LuHym+9p0HNJT95MePMWlJWWsjM9lW8/XszPS75kzPW30/OEwQ16LBE5tijgmqBzLr0aP78APlr0Lz5e9G9Cm4XT7bgBXH7TXQQFh3jKXXjVBBzHxRvzZpG7L4uO3fvw0LMLaN6q+r2tBc9Opf8pw+nY7feP0+rQtSdX3HoPrzzzFADjbruPxC513387VPOnPwJAQGAgLVpH06X38Tw29w16HD+wQY8jIsce09izzRavTNd0NhGp1ei+sd6uwtHhK/cfgpxyp3frcXRokE+H0D04ERGxkgJORESspIATERErKeBERMRKCjgREbGSAk5ERKykgBMRESsp4ERExEoKOBERsZICTkRErKSAExERKyngRETESgo4ERGxkgJORESspIATERErKeBERMRKCjgREbGSAk5ERKykgBMRESsp4ERExEoKOBERsZICTkRErKSAExERKyngRETESgo4ERGxkgJORESspIATERErKeBERMRKCjgREbGSAk5ERKykgBMRESsp4ERExEoKOBERsZICTkRErKSAExERKyngRETESgo4ERGxkgJORESspIATERErKeBERMRKCjgREbGSAk5ERKzk5+0KiEjTtnhlurercFTovCMHgA1qD0b3jW2Q/SjgGlnq5g3Meep+1v+6nNCwCIafP4ZLrrsdX19fb1dNvGDHtq289a/n2PDrz2zbvJ5uxw3gsRcXebtaIlZSwDWivJxsHrj+UuI7dOKe6fPISE3hpWkP4zgOYyfc6e3qiRds27yB5d/+jy69jqesrNTb1RGxmgKuEX30+n8oKS7irqkvEBLWDAZBQX4urz0/jfOvuMG9TJqU/qcMZ+CpIwB48o5rycne6+UaidhLk0wa0c/ffcFxg0+pEmQnjRhNSVERq5Yv9WLNxFt8fPQjJ3Kk6KetEaUnbyI2sWOVZZFtYwkMCiY9eZOXaiUi0jQo4BpRXu4+QptFVFseFh5BXs4+L9RIRKTpUMA1MmOqL3McB1PTChERaTAKuEYU1iyC/NycassL8nIJbRbuhRqJiDQdCrhGFJuQRPrWzVWW7c7YTlFhAbEJSV6qlYhI06CAa0THDzmVFd9/SWF+nmfZd5+8Q0BQED37DfJizURE7Kf3wTWiM/98Ge+/No8nJv6NC8aNJyM9hdeem8boy67Ve+CaqOLCQpZ/+zkAe3dlUJCfx5JP3wOg39DTCQwO9mb1RKxiHMdp1AMsXpneuAc4yqVu3sCcJ+9zf1RXswjOOO9Sxlyvj+pqqnZuT+W6kYNrXPf8+98THRN/hGskR4vO654FYEPXG7xcE+8b3Te2QWbhaQTXyOI7duaROf/1djXkKBEdE8/bK9K8XQ2RJkH34ERExEoKOBERsZICTkRErKR7cCIiR4HCkBhvV8E6CjgRkaNAarvR3q6CdXSJUkRErKSAExERKyngRETESgo4ERGxkiaZHEELnpvKwuene7sacpS45LrbuPT6id6uhoi1NIITERErKeBERMRK+m8CIiJyVGmo/yagEZyIiFhJASciIlZq9FmUKxe/uBnIbuzjHCNigO3ersRRQO3gpnZwUzv8Tm0BrFzMpgcffDDpD+/IcZxGfUyePNlp7GMcKw+1hdpB7aB2UFscuXbQJUoREbHSkQi4h47AMY4Vags3tYOb2sFN7fA7tYVbg7RDo79NQERExBt0iVJERKykgBMRESsp4ERExEoNHnDGmEBjzExjTKYxJt8Y844xJu4g24wzxjg1PIIaun6NxRgz3hiz1RhTZIxZbow56SDlexljvjLGFBpj0o0xDxhjGuTjabztUNrCGJNQy/f+zCNZ54ZmjDm5ou+nV5zPuHpsY12fONR2sLg/3G2M+dEYk2OM2W2MedcY07Me21nVJw6nHf5In2iMEdwM4ELgUuAkIBx4zxjje5DtCoC2lR+O4xQ1Qv0anDHmEuBp4HHgOGAJ8KExpl0t5cOBT4GdQH/gZmAScPsRqXAjOtS2qORMqn7//9eY9TwCwoBVwC1A4cEKW9wnDqkdKrGtPwwDZgMnAqcBZcBnxpiWtW1gaZ8YxiG2QyWH3ica8s15QARQAoyttCwecAEj6thuHJDn7TcX/oHz/gF44YBlG4F/1FL+BiAHCK607D4gnYqZrcfq4zDaIgFwgBO8XfdGbJM8YNxByljbJw6xHazvDxXnGQaUA6OaeJ+oTzscdp9o6BFcP8Af+GT/AsdxUoG1uBO7LsHGmBRjTJox5j1jzHENXLdGYYwJwH3enxyw6hNqP+fBwDeO41T+i/Zj3B/Tk9DQdTxSDrMt9nvTGLPLGPOdMeaiRqng0c3KPvEH2N4fmuG+gpZVR5mm0Cfq0w77HXKfaOiAa4M7jTMPWL6zYl1t1gNXAaNxX9osAr4zxnRq4Po1htaAL+5zrKyuc25TS/n9645Vh9MWecAdwMXA2cDnwEJjzGWNVcmjlK194lA1lf7wNLAS+L6OMk2hT9SnHQ67T9Trw5aNMY8C9x6k2Kl17QL3ELNGjuN8T6UTNMYswX3SN+G+7nwsOPD86jznWsrXtPxYVO+2cBwnE5haadFPxpjWwJ3Afxqnekctm/tEvTSF/mCMmQYMBYY6jlN+kOLW9on6tsMf6RP1HcHNALod5LEMyMD9F3zrA7aPovpfIrWqONmfgGNhBJeJe9R64F9UdZ1zRi3lqWObY8HhtEVNfuDY+N43JFv7REOwpj8YY6bjvkp1muM4Ww5S3No+cYjtUJN69Yl6BZzjOJmO46w7yKMAWA6UAsMrnUgc7gBcUt+aV0yD7Q3sqO823uI4Tgnu8x5+wKrh1H7O3wMnHfA2iOG4/01GckPX8Ug5zLaoSV+Oge99A7OyTzSQvljQH4wxTwN/wf1LfV09NrGyTxxGO9SkL/XpE40wK+ZZ3LN8zsA9TfwL3JcbfSuV+ZxKs+qAB4ERQIeKis/DHZQDvD3Lp57nfAnu2aPX4A7zp3FfN25fsf4fwOeVykfg/uvsNaAncAHu2VITvX0uXmiLK3B39m5AF9zX2kuA27x9Ln+wHcIq+nJf3G+BeaDiebum1CcOox1s7Q/PVHw/T8M9Ktv/CKtUxvo+cZjtcNh9ojFOIAiYCeyp6NDvAvEHlEkG5ld6PR1IAYqBXbhnCg329jfjEM97fMV5FeMexZxcad18IPmA8r2Ar3FPqNmBO+Rtmfpb77ao6LxrgPyKjv8TcJm3z6EB2mAY7vskBz7mN6U+cajtYHF/qKkNHGBypTLW94nDaYc/0if03wRERMRK+ixKERGxkgJORESspIATERErKeBERMRKCjgREbGSAk5ERKykgBMRESsp4ERExEoKOBERsdL/A2CJYAvPs4HlAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "az.plot_posterior(trace, hdi_prob=0.80, ref_val=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**R:** The probability that the signal is from sources 0 or 1 if 93.2%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Groups comparison\n",
    "\n",
    "One pretty common statistical analysis is group comparison. We may be interested in how\n",
    "well patients respond to a certain drug, the reduction of car accidents by the introduction of\n",
    "a new traffic regulation, student performance under different teaching approaches, and so\n",
    "on.\n",
    "\n",
    "Sometimes, this type of question is framed under the hypothesis testing scenario with the\n",
    "goal of declaring a result statistically significant. Relying only on statistical significance can\n",
    "be problematic for many reasons: on the one hand, statistical significance is not equivalent\n",
    "to practical significance; on the other hand, a really small effect can be declared significant\n",
    "just by collecting enough data. The idea of hypothesis testing is connected to the concept of\n",
    "p-values. This is not a fundamental connection but a cultural one; people are used to\n",
    "thinking that way mostly because that's what they learn in most introductory statistical\n",
    "courses. There is a long record of studies and essays showing that, more often than not, p-\n",
    "values are used and interpreted the wrong way, even by people who are using them on a\n",
    "daily basis.\n",
    "\n",
    "Instead of doing hypothesis testing, we are going to take a different route and we are going\n",
    "to focus on estimating the **effect size**, that is, quantifying the difference between two\n",
    "groups. One advantage of thinking in terms of effect size is that we move away from the\n",
    "yes-no questions like; Does it work?, Is there any effect? to the more nuance type of\n",
    "question like; How well does it work? How large/small is the effect?\n",
    "\n",
    "Sometimes, when comparing groups, people talk about a control group and a treatment\n",
    "group (or maybe more than one control and treatment groups). This makes sense, for\n",
    "example, when we want to test a new drug: because of the placebo effect and other reasons,\n",
    "we want to compare the new drug (the treatment) against a control group (a group not\n",
    "receiving the drug). In this case, we want to know how well one drug works compared to\n",
    "doing nothing (or, as is generally done, against the placebo effect). One interesting\n",
    "alternative question will be to ask how good a new drug is compared with the (already\n",
    "approved) most popular drug to treat that illness. In such a case, the control group cannot\n",
    "be a placebo; it should be the other drug. Bogus control groups are a splendid way to lie\n",
    "using statistics.\n",
    "\n",
    "For example, imagine you work for a dairy product company that wants to sell overly\n",
    "sugared yogurts to kids by telling their dads and moms that this particular yogurt boosts\n",
    "the immune system or help their kids grow stronger. One way to cheat and falsely back up\n",
    "your claims with data is by using milk or even water as a control group, instead of another\n",
    "cheaper, less sugary, less marketed yogurt. It may sound silly when I put this way, but\n",
    "there is a lot of actual research done this way. In fact, I am describing actual papers, not\n",
    "imaginary hypothetical scenarios. When someone says something is harder, better, faster,\n",
    "stronger, remember to ask what the baseline used for the comparison was.\n",
    "\n",
    "To compare groups, we must decide which feature (or features) we are going to use for the\n",
    "comparison. A very common feature is the mean of each group. Because we are Bayesian,\n",
    "we will work to obtain a posterior distribution of the differences of means between groups\n",
    "and not just a point-estimate of the differences. To help us see and interpret such a\n",
    "posterior, we are going to use three tools:\n",
    "\n",
    "* A posterior plot with a reference value\n",
    "* The Cohen's d\n",
    "* The probability of superiority\n",
    "\n",
    "### Cohen's d\n",
    "\n",
    "A common way to measure the effect size is Cohen's d, which is defined as follows:\n",
    "\n",
    "$$\\frac{{{\\mu _2} - {\\mu _1}}}{{\\sqrt {\\frac{{{\\sigma _1}^2 + {\\sigma _2}^2}}{2}} }}.$$\n",
    "\n",
    "According to this expression, the effect size is the difference of the means with respect to\n",
    "the pooled standard deviation of both groups. Because we can get a posterior distribution\n",
    "of means and standard deviations, we can compute a posterior distribution of Cohen's d\n",
    "values. Of course, if we just need or want one single value, we can compute the mean of\n",
    "that posterior distribution and get a single Cohen's d value. Generally, when computing a\n",
    "pooled standard deviation, we take into account the sample size of each group explicitly,\n",
    "but the previous formula is omitting the sample size of both groups. The reason for this is\n",
    "that we are getting the values of the standard deviation from the posterior and thus we are\n",
    "already accounting for the standard deviations' uncertainty.\n",
    "\n",
    "Cohen's d introduces the variability of each group by using their standard deviations. This\n",
    "is really important, as differences of one when you have a standard deviation of 0.1 seems\n",
    "large compared to the same difference when the standard deviation is 10. Also, a change of $x$\n",
    "units from one group to another could be explained by every individual data point\n",
    "changing exactly $x$ units or by half of them not changing and the other half changing\n",
    "units $2x$, and by many other combinations. Thus, including the intrinsic variations of groups is\n",
    "a way to put the differences in context. Re-scaling (standardizing) the differences helps us\n",
    "make sense of the importance of the different between groups, even when we are not very\n",
    "familiar with the scale used for the measurements.\n",
    "\n",
    "A Cohen's d can be interpreted as a Z-score (a standard score). A Z-\n",
    "score is the signed number of standard deviations by which a value differs\n",
    "from the mean value of what is being observed or measured. Thus, a\n",
    "Cohen's d of 0.5 could be interpreted as a difference of 0.5 standard\n",
    "deviation of one group with respect to the other.\n",
    "\n",
    "Even when the differences of means are standardized, we may still need to calibrate\n",
    "ourselves based on the context of a given problem to be able to say if a given value is big,\n",
    "small, medium, and so on. Fortunately, this calibration can be acquired with enough practice.\n",
    "Just as an example, if we are used to performing several analyses for more or less the same\n",
    "type of problems, we can get used to a Cohen's d of say 1, so when we get a Cohen's d of\n",
    "say 2, we know we have something important (or someone made a mistake somewhere!). If\n",
    "you do not have this practice yet, you can ask a domain expert for their valuable input. A\n",
    "very nice web page to explore what different values of Cohen's d look like is\n",
    "http://rpsychologist.com/d3/cohend. On that page, you will also find other ways to\n",
    "express an effect size; some of them could be more intuitive, such as the probability of\n",
    "superiority, which we will discuss next.\n",
    "\n",
    "### Probability of superiority\n",
    "\n",
    "This is another way to report the effect size, and this is defined as the probability that a data\n",
    "point taken at random from one group has a larger value than one also taken at random\n",
    "from the other group. If we assume that the data we are using is distributed as a normal, we\n",
    "can compute the probability of superiority from the Cohen's d using the following\n",
    "expression:\n",
    "\n",
    "$${p_s} = \\Phi \\left( {\\frac{\\delta }{{\\sqrt 2 }}} \\right)$$\n",
    "\n",
    "Here $\\Phi$, is the cumulative normal distribution and $\\delta$ is the Cohen's d. We can compute a\n",
    "point-estimate of the probability of superiority (what is usually reported), or we can\n",
    "compute the whole posterior distribution of values. If we are OK with the normality\n",
    "assumption, we can use this formula to get the probability of superiority from the Cohen's\n",
    "d. Otherwise, as we have samples from the posterior, we can directly compute it. This is a very nice advantage of using Markov chain Monte Carlo\n",
    "(MCMC) methods; once we get samples from the posterior, we can compute many\n",
    "quantities from it.\n",
    "\n",
    "## An example\n",
    "\n",
    "To explore the subject matter of this section, we are going to use the tips dataset. This data\n",
    "was reported for the first time by Bryant, P. G. and Smith, M (1995) in Practical Data\n",
    "Analysis: Case Studies in Business Statistics. You can download the data file from here https://teams.microsoft.com/l/file/AEB26E16-B2F1-41B9-99CC-627BD65E2374?tenantId=a37861db-695a-43bf-9e29-5179a055805b&fileType=csv&objectUrl=https%3A%2F%2Fparisnanterrefr.sharepoint.com%2Fsites%2FStatistiqueetanalyseBaysienne%2FDocuments%20partages%2FGeneral%2Ftips.csv&baseUrl=https%3A%2F%2Fparisnanterrefr.sharepoint.com%2Fsites%2FStatistiqueetanalyseBaysienne&serviceName=teams&threadId=19:e1e15013961d4502a255bfd237146fd0@thread.tacv2&groupId=81d7aa5f-4325-42a3-8142-b6919c20c742\n",
    "\n",
    "We want to study the effect of the day of the week on the amount of tips at a restaurant. For\n",
    "this example, the different groups are the days. Notice there is no control group or\n",
    "treatment group. If we wish, we can arbitrarily establish one day, for example, Thursday, as\n",
    "the reference or control. For now, let's start the analysis by loading the dataset as a pandas\n",
    "DataFrame using just one line of code. If you are not familiar with pandas, the tail\n",
    "command is used to show the last rows of a DataFrame (you can also try using head):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "     total_bill   tip     sex smoker   day    time  size\n239       29.03  5.92    Male     No   Sat  Dinner     3\n240       27.18  2.00  Female    Yes   Sat  Dinner     2\n241       22.67  2.00    Male    Yes   Sat  Dinner     2\n242       17.82  1.75    Male     No   Sat  Dinner     2\n243       18.78  3.00  Female     No  Thur  Dinner     2",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>total_bill</th>\n      <th>tip</th>\n      <th>sex</th>\n      <th>smoker</th>\n      <th>day</th>\n      <th>time</th>\n      <th>size</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>239</th>\n      <td>29.03</td>\n      <td>5.92</td>\n      <td>Male</td>\n      <td>No</td>\n      <td>Sat</td>\n      <td>Dinner</td>\n      <td>3</td>\n    </tr>\n    <tr>\n      <th>240</th>\n      <td>27.18</td>\n      <td>2.00</td>\n      <td>Female</td>\n      <td>Yes</td>\n      <td>Sat</td>\n      <td>Dinner</td>\n      <td>2</td>\n    </tr>\n    <tr>\n      <th>241</th>\n      <td>22.67</td>\n      <td>2.00</td>\n      <td>Male</td>\n      <td>Yes</td>\n      <td>Sat</td>\n      <td>Dinner</td>\n      <td>2</td>\n    </tr>\n    <tr>\n      <th>242</th>\n      <td>17.82</td>\n      <td>1.75</td>\n      <td>Male</td>\n      <td>No</td>\n      <td>Sat</td>\n      <td>Dinner</td>\n      <td>2</td>\n    </tr>\n    <tr>\n      <th>243</th>\n      <td>18.78</td>\n      <td>3.00</td>\n      <td>Female</td>\n      <td>No</td>\n      <td>Thur</td>\n      <td>Dinner</td>\n      <td>2</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "metadata": {},
     "output_type": "execute_result",
     "execution_count": 5
    }
   ],
   "source": [
    "import pandas as pd\n",
    "tips = pd.read_csv('tips.csv')\n",
    "tips.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this problem we assume that the amount of tips received on a given day is a normal random variable.\n",
    "\n",
    "Just to simplify things, we are going to create three variables: the `tip` variable, representing\n",
    "the tips, the `idx` variable, a categorical dummy variable to encode the days with numbers,\n",
    "that is, [0, 1, 2, 3] instead of [Thursday, Friday, Saturday, Sunday], and finally\n",
    "the `groups` variable, with the number of groups (4):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "text": [
      "[3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2\n 2 3 3 3 3 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n 0 0 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3\n 3 3 3 3 3 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1\n 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0]\n[ 1.01  1.66  3.5   3.31  3.61  4.71  2.    3.12  1.96  3.23  1.71  5.\n  1.57  3.    3.02  3.92  1.67  3.71  3.5   3.35  4.08  2.75  2.23  7.58\n  3.18  2.34  2.    2.    4.3   3.    1.45  2.5   3.    2.45  3.27  3.6\n  2.    3.07  2.31  5.    2.24  2.54  3.06  1.32  5.6   3.    5.    6.\n  2.05  3.    2.5   2.6   5.2   1.56  4.34  3.51  3.    1.5   1.76  6.73\n  3.21  2.    1.98  3.76  2.64  3.15  2.47  1.    2.01  2.09  1.97  3.\n  3.14  5.    2.2   1.25  3.08  4.    3.    2.71  3.    3.4   1.83  5.\n  2.03  5.17  2.    4.    5.85  3.    3.    3.5   1.    4.3   3.25  4.73\n  4.    1.5   3.    1.5   2.5   3.    2.5   3.48  4.08  1.64  4.06  4.29\n  3.76  4.    3.    1.    4.    2.55  4.    3.5   5.07  1.5   1.8   2.92\n  2.31  1.68  2.5   2.    2.52  4.2   1.48  2.    2.    2.18  1.5   2.83\n  1.5   2.    3.25  1.25  2.    2.    2.    2.75  3.5   6.7   5.    5.\n  2.3   1.5   1.36  1.63  1.73  2.    2.5   2.    2.74  2.    2.    5.14\n  5.    3.75  2.61  2.    3.5   2.5   2.    2.    3.    3.48  2.24  4.5\n  1.61  2.   10.    3.16  5.15  3.18  4.    3.11  2.    2.    4.    3.55\n  3.68  5.65  3.5   6.5   3.    5.    3.5   2.    3.5   4.    1.5   4.19\n  2.56  2.02  4.    1.44  2.    5.    2.    2.    4.    2.01  2.    2.5\n  4.    3.23  3.41  3.    2.03  2.23  2.    5.16  9.    2.5   6.5   1.1\n  3.    1.5   1.44  3.09  2.2   3.48  1.92  3.    1.58  2.5   2.    3.\n  2.72  2.88  2.    3.    3.39  1.47  3.    1.25  1.    1.17  4.67  5.92\n  2.    2.    1.75  3.  ]\n"
     ],
     "output_type": "stream"
    }
   ],
   "source": [
    "tip = tips['tip'].values\n",
    "idx = pd.Categorical(tips['day'], categories=['Thur', 'Fri', 'Sat', 'Sun']).codes\n",
    "groups = len(np.unique(idx))\n",
    "print(idx)\n",
    "print(tip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model for this problem is almost the same as the normal model we saw last class; the only difference is that now $\\mu$\n",
    "and $\\sigma$ are going to be vectors instead of scalar variables. PyMC3 syntax is extremely\n",
    "helpful for this situation: instead of writing `for loops`, we can write our model in a\n",
    "vectorized way. This means that for the priors, we pass a shape argument and for the\n",
    "likelihood, we properly index the means and sds variables using the idx variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [sigma, mu]\n"
     ],
     "output_type": "stream"
    },
    {
     "data": {
      "text/plain": "<IPython.core.display.HTML object>",
      "text/html": "\n    <div>\n        <style>\n            /* Turns off some styling */\n            progress {\n                /* gets rid of default border in Firefox and Opera. */\n                border: none;\n                /* Needs to be in here for Safari polyfill so background images work as expected. */\n                background-size: auto;\n            }\n            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n                background: #F44336;\n            }\n        </style>\n      <progress value='0' class='' max='24000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      \n    </div>\n    "
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mRemoteTraceback\u001b[0m                           Traceback (most recent call last)",
      "\u001b[1;31mRemoteTraceback\u001b[0m: \n\"\"\"\nTraceback (most recent call last):\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\pymc3\\parallel_sampling.py\", line 114, in _unpickle_step_method\n    self._step_method = pickle.loads(self._step_method)\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\theano\\compile\\function_module.py\", line 1082, in _constructor_Function\n    f = maker.create(input_storage, trustme=True)\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\theano\\compile\\function_module.py\", line 1715, in create\n    input_storage=input_storage_lists, storage_map=storage_map)\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\theano\\gof\\link.py\", line 699, in make_thunk\n    storage_map=storage_map)[:3]\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\theano\\gof\\vm.py\", line 1091, in make_all\n    impl=impl))\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\theano\\gof\\op.py\", line 955, in make_thunk\n    no_recycling)\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\theano\\gof\\op.py\", line 858, in make_c_thunk\n    output_storage=node_output_storage)\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\theano\\gof\\cc.py\", line 1217, in make_thunk\n    keep_lock=keep_lock)\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\theano\\gof\\cc.py\", line 1157, in __compile__\n    keep_lock=keep_lock)\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\theano\\gof\\cc.py\", line 1624, in cthunk_factory\n    key=key, lnk=self, keep_lock=keep_lock)\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\theano\\gof\\cmodule.py\", line 1189, in module_from_key\n    module = lnk.compile_cmodule(location)\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\theano\\gof\\cc.py\", line 1527, in compile_cmodule\n    preargs=preargs)\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\theano\\gof\\cmodule.py\", line 2399, in compile_str\n    (status, compile_stderr.replace('\\n', '. ')))\nException: ('The following error happened while compiling the node', AdvancedSubtensor1(mu, TensorConstant{[3 3 3 3 3..2 2 2 2 0]}), '\\n', 'Compilation failed (return status=3): ', '[AdvancedSubtensor1(mu, TensorConstant{[3 3 3 3 3..2 2 2 2 0]})]')\n\nDuring handling of the above exception, another exception occurred:\n\nTraceback (most recent call last):\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\pymc3\\parallel_sampling.py\", line 135, in run\n    self._unpickle_step_method()\n  File \"C:\\Users\\carlo\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\pymc3\\parallel_sampling.py\", line 116, in _unpickle_step_method\n    raise ValueError(unpickle_error)\nValueError: The model could not be unpickled. This is required for sampling with more than one core and multiprocessing context spawn or forkserver.\n\"\"\"",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;31mValueError\u001b[0m: The model could not be unpickled. This is required for sampling with more than one core and multiprocessing context spawn or forkserver.",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[1;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-7-689138b867b2>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      8\u001b[0m     \u001b[1;31m# We use the variable idx to index the mean and the standard deviation\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      9\u001b[0m     \u001b[0my\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mNormal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'y'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmu\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mmu\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0midx\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msd\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msigma\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0midx\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mobserved\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mtip\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 10\u001b[1;33m     \u001b[0mtrace_cg\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msample\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m5000\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     11\u001b[0m     \u001b[0maz\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot_trace\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtrace_cg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     12\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\pymc3\\sampling.py\u001b[0m in \u001b[0;36msample\u001b[1;34m(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs)\u001b[0m\n\u001b[0;32m    543\u001b[0m         \u001b[0m_print_step_hierarchy\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstep\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    544\u001b[0m         \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 545\u001b[1;33m             \u001b[0mtrace\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_mp_sample\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0msample_args\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mparallel_args\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    546\u001b[0m         \u001b[1;32mexcept\u001b[0m \u001b[0mpickle\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mPickleError\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    547\u001b[0m             \u001b[0m_log\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwarning\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Could not pickle model, sampling singlethreaded.\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\pymc3\\sampling.py\u001b[0m in \u001b[0;36m_mp_sample\u001b[1;34m(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs)\u001b[0m\n\u001b[0;32m   1484\u001b[0m         \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1485\u001b[0m             \u001b[1;32mwith\u001b[0m \u001b[0msampler\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1486\u001b[1;33m                 \u001b[1;32mfor\u001b[0m \u001b[0mdraw\u001b[0m \u001b[1;32min\u001b[0m \u001b[0msampler\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1487\u001b[0m                     \u001b[0mtrace\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtraces\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mdraw\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mchain\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mchain\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1488\u001b[0m                     \u001b[1;32mif\u001b[0m \u001b[0mtrace\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msupports_sampler_stats\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mdraw\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstats\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\pymc3\\parallel_sampling.py\u001b[0m in \u001b[0;36m__iter__\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    490\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    491\u001b[0m         \u001b[1;32mwhile\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_active\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 492\u001b[1;33m             \u001b[0mdraw\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mProcessAdapter\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrecv_draw\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_active\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    493\u001b[0m             \u001b[0mproc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mis_last\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdraw\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtuning\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstats\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mwarns\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdraw\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    494\u001b[0m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_total_draws\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\miniconda3\\envs\\bayesianStats\\lib\\site-packages\\pymc3\\parallel_sampling.py\u001b[0m in \u001b[0;36mrecv_draw\u001b[1;34m(processes, timeout)\u001b[0m\n\u001b[0;32m    363\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    364\u001b[0m                 \u001b[0merror\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Chain %s failed.\"\u001b[0m \u001b[1;33m%\u001b[0m \u001b[0mproc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mchain\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 365\u001b[1;33m             \u001b[1;32mraise\u001b[0m \u001b[0merror\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0mold_error\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    366\u001b[0m         \u001b[1;32melif\u001b[0m \u001b[0mmsg\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m\"writing_done\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    367\u001b[0m             \u001b[0mproc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_readable\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mRuntimeError\u001b[0m: Chain 3 failed."
     ],
     "ename": "RuntimeError",
     "evalue": "Chain 3 failed.",
     "output_type": "error"
    }
   ],
   "source": [
    "with pm.Model() as comparing_groups:\n",
    "    # mu and sigma are vectors because we are passing the parameter shape.\n",
    "    # With this, we are telling PyMC3 that we want to analyze 4 different mu's\n",
    "    # and 4 different sigma's\n",
    "    mu = pm.Normal('mu', mu=0, sd=10, shape=groups)\n",
    "    sigma = pm.HalfNormal('sigma', sd=10, shape=groups)\n",
    "    \n",
    "    # We use the variable idx to index the mean and the standard deviation\n",
    "    y = pm.Normal('y', mu=mu[idx], sd=sigma[idx], observed=tip)\n",
    "    trace_cg = pm.sample(5000)\n",
    "    az.plot_trace(trace_cg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code is just a way of plotting the difference without repeating the\n",
    "comparison. Instead of plotting the all-against-all matrix, we are just plotting the upper\n",
    "triangular portion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "dist = stats.norm()\n",
    "\n",
    "# This is just preparations for the plots\n",
    "_, ax = plt.subplots(3, 2, figsize=(14, 8), constrained_layout=True)\n",
    "comparisons = [(i, j) for i in range(4) for j in range(i+1, 4)]\n",
    "pos = [(k, l) for k in range(3) for l in (0, 1)]\n",
    "\n",
    "\n",
    "for (i, j), (k, l) in zip(comparisons, pos):\n",
    "    # Calculate the difference between the means\n",
    "    means_diff = trace_cg['mu'][:, i] - trace_cg['mu'][:, j]\n",
    "    # Obtain Cohen's d\n",
    "    d_cohen = (means_diff / np.sqrt((trace_cg['sigma'][:, i]**2 +trace_cg['sigma'][:, j]**2) / 2)).mean()\n",
    "    # Obtain the probability of superiority\n",
    "    ps = dist.cdf(d_cohen/(2**0.5))\n",
    "    # Plot the posterior distribution of the difference between the means\n",
    "    az.plot_posterior(means_diff, ref_val=0, ax=ax[k, l])\n",
    "    \n",
    "    # Draw the legend\n",
    "    ax[k, l].set_title(f'$\\mu_{i}-\\mu_{j}$')\n",
    "    ax[k, l].plot(0, label=f\"Cohen's d = {d_cohen:.2f}\\nProb sup = {ps:.2f}\",alpha=0)\n",
    "    ax[k, l].legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to interpret these results is by comparing the reference value with the HPD\n",
    "interval. According to the previous diagram, we have only one case when the 94% HPD\n",
    "excludes the reference value of zero, that is, the difference in tips between Thursday and\n",
    "Sunday ($\\mu_0-\\mu_3$). For all the other examples, we cannot rule out a difference of zero (according to the\n",
    "HPD-reference-value-overlap criteria). However, the plots suggest that tips on Sundays are a little bit better.\n",
    "But even for that case, is an average difference of\n",
    "≈0.5 dollars large enough? Is that difference enough to accept working on Sunday and\n",
    "missing the opportunity to spend time with family or friends? Is that difference enough to\n",
    "justify averaging the tips over the four days and giving every waitress and waiter the same\n",
    "amount of tip money? Those kinds of questions cannot be answered by statistics; they can\n",
    "only be informed by statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "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.7.9"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}