bayesian-statistics/0_simple-bayesian-estimation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "latin-question",
   "metadata": {},
   "source": [
    "# Simple Example for Bayesian Parameter Estimation\n",
    "We are going to show how to implement MCMC for a simple example: determining whether a coin is fair"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "forty-finance",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import stats\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "applicable-employer",
   "metadata": {},
   "source": [
    "Configuration parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "historical-kidney",
   "metadata": {},
   "outputs": [],
   "source": [
    "actual_theta = 0.76\n",
    "n_samples = 4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "super-bible",
   "metadata": {},
   "source": [
    "## Problem Statement\n",
    "\n",
    "Our goal is to determine whether a coin is fair by using Bayesian statistics to evaluate some experimental evidence. \n",
    "\n",
    "Bayesian methods are built around a process of updating your prior beliefs about a system given evidence about how it performs to form a new \"posterior\" set of beliefs.\n",
    "\n",
    "We define these beliefs as models and probability distributions. For example, let's model our coin flipping example as model:\n",
    "\n",
    "$P(H) = \\theta$\n",
    "\n",
    "where we define the fairness of the coin as a variable, $\\theta$.\n",
    "\n",
    "## Prior distribution\n",
    "\n",
    "Starting out, we make no assuptions about the fairness of the coin and assume $\\theta$ can take on any value between 0 and 1.\n",
    "\n",
    "We represent this prior belief using a uniform probability distribution and a so-called diffuse prior: $P(\\theta) = U(0, 1)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "median-record",
   "metadata": {},
   "outputs": [],
   "source": [
    "prior = stats.uniform(0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "static-mambo",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Prior Distribution')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEYCAYAAACz2+rVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAU6klEQVR4nO3df7BfdX3n8eerCbhrAUETKoTEgBvRtBUGI6AjllarCdZJnWGHH65YWjfLCG53uzsLOlXacejqOG0tCk0zlLLurmB3YSTUKG53q1gFl7ANPwKCafgVQksAEQU6EHjvH9/Dh28uN9xvknvPvTc8HzN38j3nfL7nvD/35vt9fT/nnO85qSokSQL4mekuQJI0cxgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBc16SX6a5Iget3dCkjsncX1fT/Lh7vFvJPnbSVz3B5N8c7LWp72foaAZJ8k9SZ7q3uz/MclfJNlvZ+2rar+q2jxJ2/69JM8k+Un3c1eSLyY5ZGh736mqI0dc13+bqF1Vraiq/zIJtS9OUknmDq37v1fVe/Z03Xr5MBQ0U72/qvYDjgHeCvzu2AbDb3674yWe/5Wq2h94NfAB4LXATcPBMBky4GtQM4r/ITWjVdUDwNeBXwDoPgmfneSHwA+H5v2L7vGrknwpybYk9yb53effeLtdM99N8sdJHgV+b4JtP1NVG4FTgG3Af+jWc2KSLc+3S3Jukge6kcWdSd6VZDnwCeCUbsRzc9f2W0kuSPJd4EngiG7eR4Y2nSRfSPLjJD9I8q6hBfckeffQ9PBo5Lru38e6bb5t7O6oJG9PcmO37huTvH1o2beSfLr7Hf0kyTeTzJvgT6S9jKGgGS3JQuAk4O+GZv86cBywdJynfAF4FXAE8EvAGcCZQ8uPAzYDBwMXjFJDVT0LXA2cME59RwLnAG/tRhfvBe6pqm8Af8Bg1LFfVR019LQPAauA/YF7x9nk8zXOA84Hrkry6hFKfWf374HdNq8fU+urga8BFwKvAf4I+FqS1ww1O53B7+tgYF/gP46wXe1FDAXNVF9N8hjwt8C3GbzBPu8/V9WjVfXU8BOSzGHwqf7jVfWTqroH+EMGb8LP21pVX6iq7WOfP4GtDHYnjfUs8ApgaZJ9quqeqvr7CdZ1WVVt7Gp4ZpzlDwGf70YqXwHuBN63C7XuzPuAH1bVf+22fTnwA+D9Q23+oqru6n43fwkcPQnb1SxiKGim+vWqOrCqXldVHx3zBn7/Tp4zj8Gn2+FP3/cCC0Z47kQWAI+OnVlVm4B/x2BX1ENJrkhy6ATrmqiGB2rHK1XeC0y0zlEcyotHJmN/P/8w9PhJYKcH+LV3MhQ0G+3s0r4PA88Arxuatwh4YITn7lR3TOL9wHfGLabqy1X1jm67BXx2gm1NVMOCJBmaXsRgpALwBPDKoWWv3YX1bmXH383z635gnLZ6mTIUtNfo9v3/JXBBkv2TvA74HWDC00LHk2SfJG8CLmfw5vtH47Q5MsmvJHkF8E/AUwx2KQH8I7B4N84wOhj4t932/yXwJmBdt2wDcGq3bBlw8tDztgHPMTieMp51wBuSnJ5kbpJTGByX+atdrE97MUNBe5uPMfg0vZnB8YgvA5fu4jpOSfJT4DFgLfAI8Jaq2jpO21cAn2EwSvkHBm/on+iW/Y/u30eS/L9d2P73gSXdOi8ATq6qR7plnwReD/wI+H0G/QOgqp7s2n83yWNJjh9eabeOX2NwFtUjwH8Cfq2qHt6F2rSXizfZkSQ9z5GCJKkxFCRJjaEgSWoMBUlSs0cXFJtu8+bNq8WLF093GZI0q9x0000PV9X88ZbN6lBYvHgx69evn+4yJGlWSTLeNbcAdx9JkoYYCpKkxlCQJDWGgiSpMRQkSU0voZDk0iQPJbltJ8uT5MIkm5LckuSYPuqSJO2or5HCZcDyl1i+gsFVIZcwuE3hn/ZQkyRpjF6+p1BV1yVZ/BJNVgJf6u42dUOSA5McUlUP9lHfVPny9+/j6g3ev0R6uVl59AJOP27RdJexW2bKMYUF7HiLwi3seIvAJsmqJOuTrN+2bVsvxe2uqzc8wO0PPj7dZUjq0e0PPj6rPwzOlG80Z5x5497ooarWAGsAli1bNuNvBrH0kAP4yr9523SXIaknp/zZ9dNdwh6ZKSOFLcDCoenDeOGetJKknsyUUFgLnNGdhXQ88OPZfjxBkmajXnYfJbkcOBGYl2QLcD6wD0BVrWZwQ/GTgE3Ak8CZfdQlSdpRX2cfnTbB8gLO7qMWSdLOzZTdR5KkGcBQkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktQYCpKkxlCQJDWGgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSmt5CIcnyJHcm2ZTkvHGWvyrJNUluTrIxyZl91SZJGuglFJLMAS4CVgBLgdOSLB3T7Gzg9qo6CjgR+MMk+/ZRnyRpoK+RwrHApqraXFVPA1cAK8e0KWD/JAH2Ax4FtvdUnySJ/kJhAXD/0PSWbt6wLwJvArYCtwK/XVXPjV1RklVJ1idZv23btqmqV5JelvoKhYwzr8ZMvxfYABwKHA18MckBL3pS1ZqqWlZVy+bPnz/ZdUrSy1pfobAFWDg0fRiDEcGwM4GramATcDfwxp7qkyTRXyjcCCxJcnh38PhUYO2YNvcB7wJI8nPAkcDmnuqTJAFz+9hIVW1Pcg5wLTAHuLSqNiY5q1u+Gvg0cFmSWxnsbjq3qh7uoz5J0kAvoQBQVeuAdWPmrR56vBV4T1/1SJJezG80S5IaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktQYCpKkxlCQJDWGgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJTW+hkGR5kjuTbEpy3k7anJhkQ5KNSb7dV22SpIG5fWwkyRzgIuBXgS3AjUnWVtXtQ20OBC4GllfVfUkO7qM2SdIL+hopHAtsqqrNVfU0cAWwckyb04Grquo+gKp6qKfaJEmdvkJhAXD/0PSWbt6wNwAHJflWkpuSnDHeipKsSrI+yfpt27ZNUbmS9PLUVyhknHk1Znou8BbgfcB7gU8mecOLnlS1pqqWVdWy+fPnT36lkvQy1ssxBQYjg4VD04cBW8dp83BVPQE8keQ64Cjgrn5KlCT1NVK4EViS5PAk+wKnAmvHtLkaOCHJ3CSvBI4D7uipPkkSPY0Uqmp7knOAa4E5wKVVtTHJWd3y1VV1R5JvALcAzwGXVNVtfdQnSRroa/cRVbUOWDdm3uox058DPtdXTZKkHfmNZklSYyhIkhpDQZLUGAqSpMZQkCQ1uxwKSX62u8CdJGkvM2EoJPmZJKcn+VqSh4AfAA92l7f+XJIlU1+mJKkPo4wU/gZ4PfBx4LVVtbCqDgZOAG4APpPkX01hjZKknozy5bV3V9UzY2dW1aPAlcCVSfaZ9MokSb2bcKQwXiDsThtJ0sw34UghyWLgbAa7kB4FNgDXVNW9U1qZJKl3oxxTuJrBweXnb6d5FHBdkouSvGIqi5Mk9WuUUJhTVX9eVf8beLSq/jWDUcM9wJqpLE6S1K9RQuGvu8teQ3e3tKra3l3R9G1TVpkkqXejnH30O8DHk6wHDk2yCniSQSA8MpXFSZL6NcpIoarqAuCdwCrgtQzupXwbsAIgyXj3YJYkzTKjjBT+JsmVwNVVtZbuNprdbTXfkeTDDL7gdtmUVSlJ6sUoobAc+E3g8iSHA48B/4zBbTW/CfxxVW2YqgIlSf2ZMBSq6p+Ai4GLu28uzwOeqqrHprg2SVLPRrkg3oeTPJzkUeAS4KcGgiTtnUY50PxJBl9aeyNwH/AHU1qRJGnajHJM4fGq+rvu8SeTfH8qC5IkTZ9RQuGQ7rsJdzC43IVXRJWkvdQooXA+8Gbgg8AvAvslWQfcDNxSVZdPYX2SpB6NcvbRDtc3SnIYg5D4ReAkwFCQpL3EKCOFHVTVFmALsG7yy5EkTadRzj6SJL1MGAqSpMZQkCQ1hoIkqTEUJElNb6GQZHmSO5NsSnLeS7R7a5Jnk5zcV22SpIFeQiHJHOAiBjflWQqclmTpTtp9Fri2j7okSTvqa6RwLLCpqjZX1dPAFcDKcdp9DLgSeKinuiRJQ/oKhQXA/UPTW7p5TZIFwAeA1S+1oiSrkqxPsn7btm2TXqgkvZz1FQrj3cO5xkx/Hji3qp59qRVV1ZqqWlZVy+bPnz9Z9UmS2I3LXOymLcDCoenDgK1j2iwDrkgCg7u7nZRke1V9tZcKJUm9hcKNwJLuHs8PAKcCpw83qKrDn3+c5DLgrwwESepXL6FQVduTnMPgrKI5wKVVtTHJWd3ylzyOIEnqR18jBapqHWOurLqzMKiq3+ijJknSjvxGsySpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktQYCpKkxlCQJDWGgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1PQWCkmWJ7kzyaYk542z/INJbul+vpfkqL5qkyQN9BIKSeYAFwErgKXAaUmWjml2N/BLVfVm4NPAmj5qkyS9oK+RwrHApqraXFVPA1cAK4cbVNX3qupH3eQNwGE91SZJ6vQVCguA+4emt3Tzdua3gK+PtyDJqiTrk6zftm3bJJYoSeorFDLOvBq3YfLLDELh3PGWV9WaqlpWVcvmz58/iSVKkub2tJ0twMKh6cOArWMbJXkzcAmwoqoe6ak2SVKnr5HCjcCSJIcn2Rc4FVg73CDJIuAq4ENVdVdPdUmShvQyUqiq7UnOAa4F5gCXVtXGJGd1y1cDnwJeA1ycBGB7VS3roz5J0kBfu4+oqnXAujHzVg89/gjwkb7qkSS9mN9oliQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktQYCpKkxlCQJDWGgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktT0FgpJlie5M8mmJOeNszxJLuyW35LkmL5qkyQN9BIKSeYAFwErgKXAaUmWjmm2AljS/awC/rSP2iRJL5jb03aOBTZV1WaAJFcAK4Hbh9qsBL5UVQXckOTAJIdU1YOTXczvX7OR27c+PtmrfZHbH3ycpYccMOXbkTSz3P7g45zyZ9dP6TaWHnoA57//5yd9vX2FwgLg/qHpLcBxI7RZAOwQCklWMRhJsGjRokkvdDItPeQAVh69YLrLkNSj2f6a7ysUMs682o02VNUaYA3AsmXLXrR8FFORrpIEcPpxizj9uJn9gfWl9HWgeQuwcGj6MGDrbrSRJE2hvkLhRmBJksOT7AucCqwd02YtcEZ3FtLxwI+n4niCJGnnetl9VFXbk5wDXAvMAS6tqo1JzuqWrwbWAScBm4AngTP7qE2S9IK+jilQVesYvPEPz1s99LiAs/uqR5L0Yn6jWZLUGAqSpMZQkCQ1hoIkqcng+O7slGQbcO901zGBecDD013EJNlb+rK39APsy0w0G/rxuqqaP96CWR0Ks0GS9VW1bLrrmAx7S1/2ln6AfZmJZns/3H0kSWoMBUlSYyhMvTXTXcAk2lv6srf0A+zLTDSr++ExBUlS40hBktQYCpKkxlCYBEmWJ7kzyaYk542zPEku7JbfkuSY6ahzFCP05YNdH25J8r0kR01HnaOYqC9D7d6a5NkkJ/dZ364YpS9JTkyyIcnGJN/uu8ZRjPD/61VJrklyc9ePGXu15CSXJnkoyW07WT5rXvc7qCp/9uCHwaXA/x44AtgXuBlYOqbNScDXGdxd7njg+9Nd9x705e3AQd3jFbO5L0Pt/g+DK/iePN1178Hf5UAG9zxf1E0fPN1172Y/PgF8tns8H3gU2He6a99Jf94JHAPctpPls+J1P/bHkcKeOxbYVFWbq+pp4Apg5Zg2K4Ev1cANwIFJDum70BFM2Jeq+l5V/aibvIHBHfJmolH+LgAfA64EHuqzuF00Sl9OB66qqvsAqmom9meUfhSwf5IA+zEIhe39ljmaqrqOQX07M1te9zswFPbcAuD+oekt3bxdbTMT7Gqdv8Xgk9BMNGFfkiwAPgCsZmYb5e/yBuCgJN9KclOSM3qrbnSj9OOLwJsY3Ir3VuC3q+q5fsqbdLPldb+D3m6ysxfLOPPGnuc7SpuZYOQ6k/wyg1B4x5RWtPtG6cvngXOr6tnBB9MZa5S+zAXeArwL+OfA9UluqKq7prq4XTBKP94LbAB+BXg98L+SfKeqHp/i2qbCbHnd78BQ2HNbgIVD04cx+JSzq21mgpHqTPJm4BJgRVU90lNtu2qUviwDrugCYR5wUpLtVfXVXioc3aj/xx6uqieAJ5JcBxwFzKRQGKUfZwKfqcFO+U1J7gbeCPzffkqcVLPldb8Ddx/tuRuBJUkOT7IvcCqwdkybtcAZ3dkIxwM/rqoH+y50BBP2Jcki4CrgQzPsU+hYE/alqg6vqsVVtRj4n8BHZ2AgwGj/x64GTkgyN8krgeOAO3qucyKj9OM+BqMdkvwccCSwudcqJ89sed3vwJHCHqqq7UnOAa5lcHbFpVW1MclZ3fLVDM5sOQnYBDzJ4NPQjDNiXz4FvAa4uPuEvb1m4BUhR+zLrDBKX6rqjiTfAG4BngMuqapxT5WcLiP+TT4NXJbkVga7X86tqhl5GeoklwMnAvOSbAHOB/aB2fW6H8vLXEiSGncfSZIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBWmSJZmT5E+6+wHcmuSI6a5JGpWhIE2+jwObq+rngQuBj05zPdLIvMyFNImS/Czwgap6SzfrbuB901iStEsMBWlyvRtYmGRDN/1q4K+nrxxp17j7SJpcRwOfqqqjq+po4JsM7g8gzQqGgjS5DmJwRUySzAXeA1wzrRVJu8BQkCbXXQxu0g7w74GvVdXd01iPtEu8dLY0iZIcxOC+1fOA64FVVfXU9FYljc5QkCQ17j6SJDWGgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1Px/zSH3yGQz4wgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "x = np.linspace(-0.1, 1.1, 128)\n",
    "ax.step(x, prior.pdf(x), where='mid')\n",
    "\n",
    "ax.set_xlabel('$\\\\theta$')\n",
    "ax.set_ylabel('$P(\\\\theta)$')\n",
    "ax.set_title('Prior Distribution')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "subtle-creator",
   "metadata": {},
   "source": [
    "## Introducing Bayes's Rule\n",
    "Bayes's rule lets us update this prior distribution as we gather more data about the problem to form a \"posterior\" distribution that defines our belief about the parameters after we gather data.\n",
    "\n",
    "Expressed mathematically, we can write Bayes rule as\n",
    "\n",
    "$P(\\theta|D, M) = \\frac{P(D|\\theta,M)P(\\theta|M)}{P(D|M)}$\n",
    "\n",
    "where \n",
    "- $P(\\theta|D, M)$ is our posterior - a function that captures the probability of parameters $\\theta$ given our observations ($D$) and chosen model form ($M$).\n",
    "- $P(D|\\theta,M)$ is our likelihood function, the probability of observing data given the model form and choice of parmaeters.\n",
    "- $P(D|M)$ is evidence, our total probability of observing the data given the chosen model.\n",
    "- $P(\\theta|M)$ is our prior beliefs - what we assumed for the parameters before gathering any data.\n",
    "\n",
    "## Performing some experiments\n",
    "Now that we have set up our problem, let's flip a coin many times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "alternative-designation",
   "metadata": {},
   "outputs": [],
   "source": [
    "coin = stats.bernoulli(actual_theta)  # Make an engine that flips a heads (1), 60% of the time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "catholic-runner",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Flipped the coin 4 times\n"
     ]
    }
   ],
   "source": [
    "samples = coin.rvs(n_samples)\n",
    "print(f'Flipped the coin {len(samples)} times')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "extensive-december",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The coin flipped heads 75.0% of the time. So, it seems biased but we can be more sure\n"
     ]
    }
   ],
   "source": [
    "print(f'The coin flipped heads {samples.mean()*100:.1f}% of the time. So, it seems biased but we can be more sure')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fantastic-cache",
   "metadata": {},
   "source": [
    "## Sampling the Posterior\n",
    "Many analyses you want to perform involve computing properties of the posterior distribution (again, likely distribution of model parameters given your data). \n",
    "\n",
    "Solving that distribution (i.e., the function that gives probability as a function of parameters) analytically is rarely possible. But, there are techniques to be able to draw samples from this distribution without knowing its true shape, such as [Markov-Chain Monte Carlo (MCMC) methods](https://www.nature.com/articles/s43586-020-00001-2/tables/1).\n",
    "\n",
    "Here, we implement a simple one [\"Metropolis-Hasting Monte Carlo\"](https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm).\n",
    "\n",
    "This algorithm is as follows:\n",
    "1. Begin with a point, $\\theta$, chosen at random from the prior distribution\n",
    "2. Generate a new point, $\\theta_{new}$ by adding random values to $\\theta$.\n",
    "3. Compute the relative posterior probability of $\\theta_{new}$ and $\\theta$.\n",
    "   $r = \\frac{P(\\theta_{new}|D,M)}{P(\\theta|D,M)}=\\frac{P(D|\\theta_{new},M)P(\\theta_{new}|M)}{P(D|\\theta,M)P(\\theta|M)}$\n",
    "   \n",
    "   The second step in that derivation is by applying Bayes's Rule. Note that the evidence term, $P(D|M)$, cancels out\n",
    "4. Pick a random number ($\\beta$) from a uniform distribution between 0 and 1.\n",
    "5. If $\\beta > min(1, r)$, then let $\\theta = \\theta_{new}$\n",
    "6. Repeat from 2.\n",
    "\n",
    "Due to a principle known as \"detailed balance\" the distribution $\\theta$ will be equivalent to the posterior distribution.\n",
    "\n",
    "Here, we show how to implement this algorithm\n",
    "\n",
    "### Step 1: Generate a new point\n",
    "We define a function to generate a new parameter given the previous one. There are few constraints on this function, except that the function should be symmetric. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "intelligent-petroleum",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New random point given 0.5: 0.63\n"
     ]
    }
   ],
   "source": [
    "def perturb(th):\n",
    "    return th + np.random.uniform(-0.2, 0.2)\n",
    "print(f'New random point given 0.5: {perturb(0.5):.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "undefined-joseph",
   "metadata": {},
   "source": [
    "### Step 2: Compute Posterior Probabilities\n",
    "Our posterior probability has two terms: a likelihood for parameters given data ($P(D|\\theta,M)$) and given the prior distribution ($P(\\theta,M)$). \n",
    "\n",
    "The prior is easy. We defined all points to be equally likely as long as they are between -0.5 and 0.5.\n",
    "\n",
    "The likelihood is more difficult. Here, we must take the model form - a Bernoulii function - into account. \n",
    "\n",
    "See [MathWorld](https://mathworld.wolfram.com/BernoulliDistribution.html) for more details, but we can express the number of heads ($n$) given the number of flips ($N$) observed as a function of the parameter $\\theta$ as:\n",
    "\n",
    "$P(n|\\theta,N) = {N \\choose n}\\theta^n(1-\\theta)^{N-n}$\n",
    "\n",
    "So, we can write the ration between the probabilities as:\n",
    "\n",
    "$r = \\frac{P(D|\\theta_{new},M)P(\\theta_{new}|M)}{P(D|\\theta,M)P(\\theta|M)} = \\frac{\\theta_{new}^n(1-\\theta_{new})^{N-n}}{\\theta^n(1-\\theta)^{N-n}}$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "excited-bathroom",
   "metadata": {},
   "outputs": [],
   "source": [
    "def relative_prob(a, b, n, N):\n",
    "    \"\"\"\n",
    "    Args:\n",
    "        a - New point\n",
    "        b - Old point (must be within -0.5 and 0.5)\n",
    "    Returns:\n",
    "        Relative probability of a and b\n",
    "    \"\"\"\n",
    "    if a < 0 or a > 1:\n",
    "        return 0\n",
    "    return (a ** n) * (1 - a) ** (N - n) / (b ** n * (1 - b) ** (N - n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "choice-mayor",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8665364583333336"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "relative_prob(0.55, 0.6, samples.sum(), len(samples))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "legendary-abraham",
   "metadata": {},
   "source": [
    "### Step 3: Stringing it together\n",
    "Make a function that runs the MCMC loop a certain number of times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "plain-welsh",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_mcmc(n_steps):\n",
    "    \"\"\"Run MCMC\n",
    "    \n",
    "    Args:\n",
    "        n_steps - Number of steps to run\n",
    "    \"\"\"\n",
    "    \n",
    "    # Draw a point from the prior distribution\n",
    "    th = prior.rvs()\n",
    "    \n",
    "    # Loop to create new samples\n",
    "    output = np.zeros(n_steps)\n",
    "    for i in range(n_steps):\n",
    "        # Generate a new point\n",
    "        th_new = perturb(th)\n",
    "        \n",
    "        # Get the relative probability\n",
    "        r = relative_prob(th_new, th, samples.sum(), len(samples))\n",
    "        \n",
    "        # Determine if we select the new point\n",
    "        if np.random.random() < r:\n",
    "            th = th_new\n",
    "        \n",
    "        # Store the point\n",
    "        output[i] = th\n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "passive-canadian",
   "metadata": {},
   "source": [
    "## Perform the MCMC Sample\n",
    "We will run a large number of samples to generate statistics about the posterior distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "solar-width",
   "metadata": {},
   "outputs": [],
   "source": [
    "post_samples = run_mcmc(32768)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "thrown-horse",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$P(\\\\theta)$')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEICAYAAACwDehOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWFklEQVR4nO3df7Ad9Xnf8feDEOGn+SUZKEKWMRTjllhRVIQhoaTItcHqYDzEQwixcTvWhBimJnYcoA7+0SnFiVOMIYaRHYYybmFag21AFBLRYlEIYJCFBQgysgzxBWpAAgmQFFvw9I+zWi/iSjqS7vnu+fF+zdy5zzlndfXsXOl8zve7u9+NzESSJIBd2m5AktQ/DAVJUs1QkCTVDAVJUs1QkCTVdm27gZ0xZcqUnDFjRtttSNpZTz7Z+X7UUTxZ1UcddVSLDQ23hx9++MXMnDreawMdCjNmzOChhx5quw1JO+ukkzrf776bk6r67rvvbquboRcRT2/pNaePJEm1gR4pSBoSn/98o/z8VjZUrxkKkto3d26jnLuVDdVrTh9Jat/SpZ0vYOnSpSytapXnSEFS+z796c73u+/m01XtgeZ2OFKQJNUMBUlSzVCQJNUMBUlSzQPNktp36aWNslPPuHAhAE9d9qFWWhpVhoKk9h1/fKOs6lsWttTMaDMUJLXvvvs6348/nvs21WpFDPI9mmfPnp0uiCcNgXEWxHvquD95y2abppKcWto5EfFwZs4e7zUPNEuSaoaCJKnmMQVJrbt/5SoAjmu5DzlSkCQ1OFKQ1LovnzwfgNuBr33tawB8+MZn2mtohBkKkoradOYQ/OrsoccPOrx+bubMmZ3CUGiFoSCpNZsC4oSnllaPYX1V7zFjZjtNjThDQVLrzr/vRgDunTGTNVXdTSiMN+rQzvFAsySpZihIkmpOH0kaGU43bZsjBUlDa8aFC98UBNo2RwqSWnfxB86r6wMbtcozFCS1buWB0+p6cqNWeU4fSWrdySse4OQVDwCwbsUDrKtqledIQVLrPvngdwG464g5rK3qPY+Ys0M/y2MIO8eRgiSp5khB0g7pt7ufOUKYGI4UJEm1IiOFiDgMuB44GHgDWJCZV2y2TQBXAKcC64BzMnNJif4kTYxejx4cDfReqemjjcBnMnNJROwDPBwRf5uZjze2OQU4svqaA1xdfZc0wLoJigvmfaaupzRqlVckFDLzOeC5qn4lIpYDhwLNUDgNuD4zE7g/IvaLiEOqPytpwG3tU/5zb5ta17s2apVX/JhCRMwAfgPY/ETkQ4GfNR6PVc9JGnLzli9m3vLFALy2fDGvVbXKK3r2UUTsDdwEfDoz127+8jh/JMf5GfOB+QDTp0+f8B4lvVWvF5I7+0e3A3Db0SfySlXvdfSJE/73aNuKhUJETKYTCP8tM28eZ5Mx4LDG42nAs5tvlJkLgAUAs2fPfktoSNo+m0/r9MsppmpHkemj6syivwaWZ+Z/2cJmtwAfi47jgDUeT5CkskqNFE4A/gBYFhFLq+cuBqYDZOY1wO10TkddQeeU1E8U6k3STvA00eFS6uyj/8v4xwya2yTwqRL9SOotg2JwucyFpNad++GL6npqo+6lflumo18YCpJa99Ke+9b1pEat8lz7SFLrzli2iDOWLQLg1WWLeLWqVZ6hIKl1hkL/MBQkSTWPKUh6E88cGm2GgqShZ9B1z+kjSVLNkYKk1p3zu1+s67c36hJ6vdjfoDEUJLVuw+Td63qXRq3ynD6S1Lqzlyzk7CWdT+yvLFnIK0vaPQYw48KFbzkOMd5zw8hQkIbMIL55zXviHuY9cQ8Arz1xD69Vtcpz+kgaQYMWGirHUJC0XQyU4WYoSCPEN3Rti6EgSVswiiFqKEhq3ZlnXVbXBzdqlefZR5KkmiMFacgNwhTIJx+4GYBvzvkIa6p63zkfabOlkeVIQVLrTv7Jg5z8kwcBWP+TB1lf1SrPUJAk1QwFSVLNUJAk1TzQLKl1G3b9tbqORq3yDAVpwAzj+v/nfPRLdX1Qo1Z5Th9JkmqGgtTnBnEp7O11/r03cP69NwDw8r038HJVqzxDQVLrTnj6EU54+hEANjz9CBuqWuUZCpKkmgeaJaky7NN03XCkIEmqOVKQBtiwfLJ9aY+31fWkRt2PhvGU4CZDQVLrzj394rqe2qhVnqEgDYhhGRWov3lMQVLrPveD6/jcD64D4KUfXMdLVa3yHClIQ2qQRhaznnmirv+xUau8IiOFiLg2Ip6PiEe38PpJEbEmIpZWX5eU6EuS9GalRgrXAVcB129lm3syc16ZdiRJ4ykyUsjMxcDqEn+XJGnH9dMxhfdFxCPAs8BnM/OxthuSVMZz+0yp610btcrrl1BYArwjM1+NiFOB7wFHjrdhRMwH5gNMnz69WIPSzhj2C5521gX/5rN1PaVRq7y+OCU1M9dm5qtVfTswOSLG/biQmQsyc3Zmzp46dWrRPiVp2PVFKETEwRERVX0snb5WtduVpFIuWbSASxYtAGD1ogWsrmqVV2T6KCJuAE4CpkTEGPAFYDJAZl4DnAGcGxEbgfXAmZmZJXqT1L73PL+yrn/RqFVekVDIzN/bxutX0TllVRKDdeGZhktfTB9Jo2QUbq+pwWUoSJJq/XJKqqQRtvKAQ+t6cqNWeYaCpNZd/MHz6/rARq3ynD6SJNUMBUmtu/SOK7n0jisBWHXHlayqapXn9JGk1h2++pm6/mWjVnmGgtQy10VSPzEUpB7a3usRvH5BbfOYgiSp5khBaomjgl95/O2H1/VujVrlGQqSWvflufPr+oBGrfK2e/ooIvaKiEm9aEaS1K5thkJE7BIRZ0XEwoh4HngCeC4iHouIv4iIce+QJknduvzWr3L5rV8F4MVbv8qLVd3vhnFxw26mj/4PsAi4CHg0M98AiIgDgN8BLouI72bmt3vXpqRhdsgrL9b1xkY9KDYFwzCcUtxNKMzNzF9u/mRmrgZuAm6KiMkT3pkkqbhtTh+NFwg7so0kqf9tc6QQETOATwHvAlYDS4FbM/PpnnYmSSqum+mj7wNfB+4ArgUS+JOIuA3448z8xx72Jw2cYTvwWMKSQ99d17/WqFVeN6EwKTP/GiAiVmfmJyNiV+ACYAHw8V42KGn4/fm/PKeu92/Ug2xQDz53c53Coog4r6oTIDM3ZuZfAO/rWWeSpOK6GSn8MXBRRDwE/JOImA+soxMIq3rZnKTRcPV3LwXg3NMv5oWqnnr6xW22NLK6CYXMzP8UEZcDc4GZwP7Ao8B/AIiIyMzsWZfSAPBYwo7bf/3aun69Uau8ri5ei4ibgO9n5i3ALQARsRvwWxHxcToXuF3Xsy4lSUV0EwofBP4tcENEvBN4GdgdmAT8DXB5Zi7tVYOSpHK2GQqZuQH4BvCN6srlKcD6zHy5x71Jrer27BGnjTRMurl47ePAX9I5U+lW4LzMfKXXjUkaHfe+4711vXujHjTD8AGhm+mjPwPeDzwDnA9cWn2XpAlx5Qm/V9f7NWqV100orM3MH1X1n0XEA71sSJLUnm5C4ZDq2oTldO6l4IqokibUdf/jCwCc89Ev8fOqPuijX2qzpZHVTSh8Afh14PeBY4C9I+J24BHgx5l5Qw/7kzQCdt/4qyXUcqPLqbWpm7OPFjQfR8Q0OiFxDHAqYChI0pDoZqTwJpk5BowBt098O5KkNnWzIJ4kaURs90hB0nCcj95P7nrXsXW9R6NWeYaCpNZ9c85H6nrfRq3yDAVpOzhC0LArEgoRcS0wD3g+M//5OK8HcAWds5nWAedk5pISvUnw5jf7zdc6Mgh678b/fiEAZ551Gf+vqg8+67I2WxpZpQ40X0dntdUtOQU4svqaD1xdoCdJ0maKhEJmLgZWb2WT04Drs+N+YL+IOKREb5JUyowLF/b9yLNfTkk9FPhZ4/FY9dxbRMT8iHgoIh564YUXijQnSaOiX0Ihxnlu3Nt7ZuaCzJydmbOnTp3a47YkabT0y9lHY8BhjcfTgGdb6kVSYbe9+7freq9GrfL6JRRuAc6LiBuBOcCazHyu5Z4kFfLtWb8642ufWVu/0516q9QpqTcAJwFTImKMzsqrkwEy8xo66yidCqygc0rqJ0r0Jak/7P7LDQBsmLw7b1T1LpN3b7OlkVUkFDJzq7dSyswEPlWiF42ebu+1rPZc9z+/CHSuU3i+qr1OoR39Mn0k9Y1+P2VQ6iVDQZJ6aNA+ZPTLKamSpD7gSEEjZ2vrHEmjzlCQ1LrvHDO3rvdu1CrPUJDUulENhX4ctXpMQVLr9l+3hv3XrQHg9XVreL2qVZ4jBY20QTszZFhd/b3/DHSuU3ihqr1OoR2OFCRJNUNBklQzFCRJNUNBklTzQLOk1n37N06t630atcozFCS17rajT6zrvRq1ynP6SFLrDln7Aoes7dxzfePaF9i41vuvt8WRgqTWXX7bXwKd6xRerOphvk6hn6+PMRQ0lPr5P53Uz5w+kiTVHCloqDhCkHaOIwVJUs2RgqTWffPY0+v6bY1a5RkKklp31xFz6nrPRq3yDAWNDI839K/DV40BsPLAafyyqicfOK3NlkaWxxQkte7SO6/i0juvAmDVnVexqqpVnqEgSaoZCpKkmqEgSap5oFkDzwPI0sQxFCS17srjz6zrfRu1yjMUJLXu3hkz63qPRq3yDAX1tU1TQ09d9qEtvqbB956frwTg8YMO5xdVvdtBh7fZ0sgyFDRwDIPhc8ldC4DO/RRWV/Uw309hPFv7AFSSZx9JkmqGgiSp5vSR+o7TQ1J7HClIkmrFRgoR8UHgCmAS8K3MvGyz108Cvg/8tHrq5sz8cqn+JLXnz0/8eF3v16hVXpFQiIhJwF8B7wfGgB9GxC2Z+fhmm96TmfNK9KT+0S9nXag9S6YdXde7N2qVV2r66FhgRWauzMxfADcCpxX6uyX1uVljy5k1thyADWPL2VDVKq/U9NGhwM8aj8eA8W6v9L6IeAR4FvhsZj62+QYRMR+YDzB9+vQetKp+5MHn4fa5xf8V6Fyn8HJVj9p1Cv2i1EghxnkuN3u8BHhHZr4XuBL43ng/KDMXZObszJw9derUie1SkkZcqVAYAw5rPJ5GZzRQy8y1mflqVd8OTI6IKYX6kyRRbvroh8CREfFO4BngTOCs5gYRcTDw88zMiDiWTmCtKtSf+oBTRFL7ioRCZm6MiPOAO+mcknptZj4WEX9YvX4NcAZwbkRsBNYDZ2bm5lNMkqQeKnadQjUldPtmz13TqK8CvFu3NIK+fPL8uj6gUas8l7mQ1LrHG8tku2R2uwwFtcLjB2o64amlQOdmO+ur2pvttMNQkNS68++7EeiEwpqqNhTa4YJ4kqSaoaCemXHhQqeJpAHj9JGKMiSk/mYoSFIfaX5wamPlYENBUusu/sB5dX1go1Z5hoJ6zikjbcvKA6fV9eRGrfI80Kwd4kFkTaSTVzzAySseAGDdigdYV9Uqz5GCumYIqFc++eB3AbjriDmsreo9jxjvlivqNUcKkqSaoaCd4jSSNFwMBUnqU2186PKYgrbJkYA0OgwFTQiDQzvjgnmfqespjVrlGQqSWvfc26bW9a6NWuV5TEFS6+YtX8y85YsBeG35Yl6rapXnSEFS687+UedOvbcdfSKvVPVeR5/YZksjy5GCJKnmSGFEtL3yoqTBYChoXJ5NJI0mp48kqc+VvIjNkYKk1p374YvqemqjVnmGggCni9Sul/bct64nNWqVZyhIat0ZyxYB8J1j5vJqVe99zNw2W+pLJU4Y8ZjCCHJlU/WbM5YtqoPh1WWL6mBQeYbCkPKNX9KOMBQkSTWPKYwwRxKSNmcoDDnf+CVtD0NhQGx6c9/aGQcGgAbVOb/7xbp+e6NWeYbCEDAMNOg2TN69rndp1CrPA82SWnf2koWcvaTz4eaVJQt5ZYkfdNriSKGPjTcCcFSgYTTviXsA+PasD/FaVe8zy9V822Ao9CHf+CW1pdj0UUR8MCKejIgVEXHhOK9HRHy9ev3HETGrVG+SpI4iI4WImAT8FfB+YAz4YUTckpmPNzY7BTiy+poDXF19H3ieOSRpUJSaPjoWWJGZKwEi4kbgNKAZCqcB12dmAvdHxH4RcUhmPleox65t7Q180xu/xwMkDaJSoXAo8LPG4zHeOgoYb5tDgTeFQkTMB+ZXD1+NiCd3sKcpwIs7+Ge3KL4y0T9xQvVkn/uc+zwA3rep+Mq8+rmnG3UXBm6fd1Z8Zaf2+R1beqFUKMQ4z+UObENmLgAW7HRDEQ9l5uyd/TmDxH0eDe7zaOjVPpc60DwGHNZ4PA14dge2kST1UKlQ+CFwZES8MyJ2A84Ebtlsm1uAj1VnIR0HrOnH4wmSNMyKTB9l5saIOA+4E5gEXJuZj0XEH1avXwPcDpwKrADWAZ/ocVs7PQU1gNzn0eA+j4ae7HN0TvaRJMm1jyRJDYaCJKk29KEwistrdLHPv1/t648j4r6IeG8bfU6kbe1zY7t/ERGvR8QZJfvrhW72OSJOioilEfFYRPygdI8TrYt/2/tGxK0R8Ui1z70+NtlTEXFtRDwfEY9u4fWJf//KzKH9onNQ+yfA4cBuwCPAezbb5lTgf9G5TuI44IG2+y6wz8cD+1f1KaOwz43t/jedkxrOaLvvAr/n/eisGjC9evz2tvsusM8XA1+p6qnAamC3tnvfiX0+EZgFPLqF1yf8/WvYRwr18hqZ+Qtg0/IaTfXyGpl5P7BfRBxSutEJtM19zsz7MvOl6uH9dK4JGWTd/J4BzgduAp4v2VyPdLPPZwE3Z+Y/AGTmoO93N/ucwD4REcDedEJhY9k2J05mLqazD1sy4e9fwx4KW1o6Y3u3GSTbuz//js4njUG2zX2OiEOB04FrCvbVS938nv8psH9E3B0RD0fEx4p11xvd7PNVwNF0LnxdBvz7zHyjTHutmPD3r2G/n8KELa8xQLren4j4HTqh8Fs97aj3utnnrwF/mpmvdz5EDrxu9nlX4DeBk4E9gL+LiPsz8+973VyPdLPPHwCWAv8KeBfwtxFxT2au7XFvbZnw969hD4VRXF6jq/2JiF8HvgWckpmrCvXWK93s82zgxioQpgCnRsTGzPxekQ4nXrf/tl/MzNeA1yJiMfBeYFBDoZt9/gRwWXYm3FdExE+BdwMPlmmxuAl//xr26aNRXF5jm/scEdOBm4E/GOBPjU3b3OfMfGdmzsjMGcB3gD8a4ECA7v5tfx/47YjYNSL2pLMy8fLCfU6kbvb5H+iMjIiIg4CjgJVFuyxrwt+/hnqkkP25vEZPdbnPlwAHAt+oPjlvzAFeYbLLfR4q3exzZi6PiDuAHwNvAN/KzHFPbRwEXf6e/yNwXUQsozO18qeZObBLakfEDcBJwJSIGAO+AEyG3r1/ucyFJKk27NNHkqTtYChIkmqGgiSpZihIkmqGgiSpZihIkmqGgiSpZihIEywiJkXEFdV6/ssi4vC2e5K6ZShIE+8iYGVm/jPg68AftdyP1LWhXuZCKi0i9gJOz8zfrJ76KfChFluStouhIE2sucBhEbG0enwAsKi9dqTt4/SRNLFmApdk5szMnAn8DZ31/aWBYChIE2t/OqtVEhG7Av8auLXVjqTtYChIE+vv6dxAHeACYGFm/rTFfqTt4tLZ0gSKiP3p3PN6CvB3wPzMXN9uV1L3DAVJUs3pI0lSzVCQJNUMBUlSzVCQJNUMBUlSzVCQJNUMBUlS7f8DacC8ZVcAdHYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.hist(post_samples, bins=np.linspace(0, 1, 128), density=True)\n",
    "\n",
    "# Plot the actual and the mean\n",
    "ax.set_ylim(ax.get_ylim())\n",
    "ax.plot([actual_theta]*2, ax.get_ylim(), 'k--')\n",
    "ax.plot([samples.mean()]*2, ax.get_ylim(), 'r--')\n",
    "\n",
    "ax.set_xlabel('$\\\\theta$')\n",
    "ax.set_ylabel('$P(\\\\theta)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "enormous-assets",
   "metadata": {},
   "source": [
    "As shown in the above plot, the peak in the distribution is at the observed fraction of heads (red dashed) and the actual value (black dashed) is still within the probability distribution.\n",
    "\n",
    "Change a few parameters see whether: \n",
    "- The distributions moves with the value you choose for the true $\\theta$\n",
    "- More samples makes the distribution converage around the correct value\n",
    "\n",
    "Metropolis-Hastings Monte Carlo is simple to implement and understand, but has some issues. The samples are strongly-correlated (sampling based on the current position means that the points are not truely-random), it also takes a large number of samples to fully explore the distribution of parameters. These issues are addressed in more advanced MCMC methods, which we will use in the later notebooks. That said, the concept of defining a \"likelihood function\" to determine how likely a point is and sampling points from the posterior are constant through these methods and are things you can understand with our simpler method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abandoned-numbers",
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}