01Simulating risk data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulating risk data\n",
    "\n",
    "*Author: Miao Cai* [miao.cai@slu.edu](miao.cai@slu.edu)\n",
    "\n",
    "## Data generating process\n",
    "\n",
    "The data generating process of real-life transportation risk is very complex. As an illustrating example, here we assume that the risk was generated from a Poisson distribution with the following data generating process:\n",
    "\n",
    "\\begin{align}\n",
    "Y_i & \\sim \\text{Poisson}(d_i\\cdot\\lambda_i)\\\\\n",
    "\\log(\\lambda_i) & = \\beta_0 + \\beta_1x_{1i} + \\beta_2x_{2i}  + \\epsilon_i \\\\\n",
    "\\epsilon_i & \\sim \\text{Normal}(0, 2^2),\n",
    "\\end{align}\n",
    "\n",
    "where $d_i$ is the distance traveled in the $i$-th trip, $x_1$ is precipitation, and $x_2$ is road traffic.\n",
    "\n",
    "We assume the sample size $N = 10,000$ and the parameters and data has following values or distributions:\n",
    "\n",
    "- $\\beta_0 = -3$,\n",
    "- $\\beta_1 = 0.5$,\n",
    "- $\\beta_2 = 0.9$,\n",
    "- $d \\sim \\text{Poisson}(1000)$\n",
    "- $x_1 \\sim \\text{Bernoulli}(0.15)$,\n",
    "- $x_2 \\sim \\text{Beta}(2, 2)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating data\n",
    "### Import packages and print package version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python version: 3.7.5 (default, Oct 31 2019, 15:18:51) [MSC v.1916 64 bit (AMD64)]\n",
      "pandas version: 0.25.3\n",
      "numpy version: 1.17.4\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import beta\n",
    "import sys\n",
    "\n",
    "print(\"Python version: \" + sys.version)\n",
    "print(\"pandas version: \" + pd.__version__)\n",
    "print(\"numpy version: \" + np.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### set parameters and generate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_distance(N_size):\n",
    "    return(np.random.poisson(lam = 1000, size = N_size))\n",
    "def simulate_precipitation(N_size):\n",
    "    return(np.random.binomial(n = 1, p = 0.15, size = N_size))\n",
    "def simulate_traffic(N_size):\n",
    "    return(np.random.beta(a = 2, b = 2, size = N_size))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(123) # set random seed\n",
    "N = 10**4\n",
    "b0, b1, b2 = -10, 0.5, 0.9\n",
    "\n",
    "d = simulate_distance(N)\n",
    "x1 = simulate_precipitation(N)\n",
    "x2 = simulate_traffic(N)\n",
    "epsilon = np.random.normal(loc = 0, scale = 1, size = N)\n",
    "\n",
    "lambda_i = np.exp(b0 + b1*x1 + b2*x2 + epsilon)\n",
    "y = np.random.poisson(d*lambda_i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot distributions of simulated data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum potential crash in a trip is 5\n",
      "The distribution of y is: \n",
      "Counter({0: 8910, 1: 937, 2: 117, 3: 24, 4: 6, 5: 6})\n"
     ]
    }
   ],
   "source": [
    "# The distribution of y\n",
    "from collections import Counter\n",
    "print(\"The maximum potential crash in a trip is \" + str(max(y)))\n",
    "print(\"The distribution of y is: \\n\" + str(Counter(y)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEaCAYAAADZvco2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAeqklEQVR4nO3dfZxdVX3v8c/XhAQREiDEVvJgYgnUgCIwjSiUq1IlgBBEkOADtOZlXq2kyPVybaiKlCttUBGKIDVCBKIQaHxglNgU5cHSQsgEFAgYGUOQISiRxBAQiIm/+8deA4eT87BnsvecOZnv+/U6r9ln7bX2+S3OML/stfZeWxGBmZlZEV7V6gDMzGzH4aRiZmaFcVIxM7PCOKmYmVlhnFTMzKwwTipmZlYYJxUzMyuMk4qZmRXGScXaiqSrJX2+1XHUI2mlpHe08PP/RdJZrfr8ijjukbR/q+OwgeekYlZF0hpJz0t6VtJvJH1D0q552kbE/hFxe8kh1iRpLHAa8LWKsjmSuiS9KOnqqvp7SvqupOckPSbpg3n3N2sLfAk4v/BO2qDnpGJW23ERsStwMPAXwGdaHE8efw0siYjnK8rWAp8HFtSofzmwGfgT4EPAFVVnF432N2vbCbxT0uu2t1PWXpxUbFCTdJCkeyVtknQDsPNAfn5EPAH8EDigIqY3Srpd0u/ScNfxFfvWSPqrtP0Pkp5Isa+SdGRFvZr7chz7bEn3S9oo6QZJlf89jgbuqIr/OxHxPeDpynJJrwHeD3w2Ip6NiDvJEsFHmu1v1jZ97gvACuA9ffxPbm3OScUGLUkjgO8BC4E9gX8n+2M2kDFMAI4B7kvvdwK+D/wn8Frg74FvSdqvqt1+wBzgLyJiN+AoYE2jfTmP/QFgOjAZeDPZ2UmvNwGrcnZtX2BrRPyiouxnwP459jdr2+th4MCc8dgOwknFBrNDgZ2ASyLiDxGxGFjerJGkcZK+IGmJpEsk/aWk4eks4NM5P/t7kn4H3En2r/9/rohpV2BeRGyOiFuBHwCnVrXfCowEpkraKSLWRMQvm+zLc+xLI2JtRKwnS0Bvqdi3O7ApZ/92BTZWlW0Edsuxv1nbXptSTDaEOKnYYLY38ES88vkMj+Vo9/+AXwEXAb8GLib7o7cIWJnzs0+IiN0j4vUR8fGKeYq9gccj4o9VMY2rbBwR3cBZwHnAU5IWSdq7yb48x/51xfbvyf7A99rAtn/Y63kWGFVVNoqXk1Kj/c3a9toN+F3OeGwH4aRig9mTwDhJqiibmKPdxyLisoj4cUTMi4iOiHhNRByY5he2x1pggqTK/3cmAk9UV4yI6yLicOD1QAAXNtmX+9h13E82NJXHL4DhkqZUlB3Iy0m30f5mbXu9kWxYzIYQJxUbzO4CtgBnpuGrE4FpzRpFxNYSY1oGPAd8StJO6Z6U48jOgl4iaT9J75I0EngBeJ5s2KvRvlzHbmAJ8L+q4hieJvOHAcMk7SxpeEQ8B3wHOF/SayQdBswgm7+i0f5mbdPnjgQOAW7JGbvtIJxUbNCKiM3AiWST0RuAU8j+mLU6puPJrrT6LfBV4LSI+HlV1ZHAvFTn12QT7//YaF8fjl3PtcAxkl5dUfYZsqQ1F/hw2u69PPrjwKuBp4Drgb+LiMqzjUb7m7U9Hrg9ItbmjN12EPLjhM12HJL+GXgqIi5pcRzLgFkR8WAr47CB56RiZmaF8fCXmZkVxknFzMwK46RiZmaFGd7qAAbCXnvtFZMmTWp1GGZmbWXFihW/jYixfWkzJJLKpEmT6OrqanUYZmZtRVKeFSxewcNfZmZWGCcVMzMrjJOKmZkVxknFzMwK46RiZmaFcVIxM7PCOKmYmVlhnFTMzKwwTipmZlaYIXFHvVlRJs29OVe9NfOOLTkSs8HJZypmZlYYn6mYkf8MxMwaK/VMRdJ0SaskdUuaW2P/SEk3pP3LJE1K5WMk3SbpWUmXVdTfRdLNkn4uaaWkeWXGb2ZmfVNaUpE0DLgcOBqYCpwqaWpVtVnAhojYB7gYuDCVvwB8Fji7xqG/FBF/DhwEHCbp6DLiNzOzvivzTGUa0B0RqyNiM7AImFFVZwZwTdpeDBwpSRHxXETcSZZcXhIRv4+I29L2ZuBeYHyJfTAzsz4oM6mMAx6veN+TymrWiYgtwEZgTJ6DS9odOA74cZ39syV1Sepat25dH0M3M7P+KDOpqEZZ9KPOtgeWhgPXA5dGxOpadSJifkR0RETH2LF9enCZmZn1U5lJpQeYUPF+PLC2Xp2UKEYD63Mcez7wSERcUkCcZmZWkDKTynJgiqTJkkYAM4HOqjqdwOlp+yTg1ohoeKYi6fNkyeesguM1M7PtVNp9KhGxRdIcYCkwDFgQESslnQ90RUQncBWwUFI32RnKzN72ktYAo4ARkk4A3gM8A3wa+DlwrySAyyLiyrL6YdYfvvPehqpSb36MiCXAkqqycyu2XwBOrtN2Up3D1pqHMTOzQcDLtJiZWWGcVMzMrDBOKmZmVhgnFTMzK4yTipmZFcZJxczMCuOkYmZmhXFSMTOzwvjJj2Yt5DvvbUfjMxUzMyuMk4qZmRXGScXMzArjpGJmZoVxUjEzs8I4qZiZWWGcVMzMrDBOKmZmVhgnFTMzK4yTipmZFcZJxczMCuOkYmZmhXFSMTOzwjipmJlZYZxUzMysMKUmFUnTJa2S1C1pbo39IyXdkPYvkzQplY+RdJukZyVdVtXmEEkPpDaXSlKZfTAzs/xKSyqShgGXA0cDU4FTJU2tqjYL2BAR+wAXAxem8heAzwJn1zj0FcBsYEp6TS8+ejMz648yz1SmAd0RsToiNgOLgBlVdWYA16TtxcCRkhQRz0XEnWTJ5SWSXgeMioi7IiKAa4ETSuyDmZn1QZlJZRzweMX7nlRWs05EbAE2AmOaHLOnyTEBkDRbUpekrnXr1vUxdDMz648yk0qtuY7oR51+1Y+I+RHREREdY8eObXBIMzMryvASj90DTKh4Px5YW6dOj6ThwGhgfZNjjm9yTLOXTJp7c6tDKETefqyZd2zJkZg1VuaZynJgiqTJkkYAM4HOqjqdwOlp+yTg1jRXUlNEPAlsknRouurrNOCm4kM3M7P+KO1MJSK2SJoDLAWGAQsiYqWk84GuiOgErgIWSuomO0OZ2dte0hpgFDBC0gnAeyLiIeDvgKuBVwM/TC8zMxsEyhz+IiKWAEuqys6t2H4BOLlO20l1yruAA4qL0szMiuI76s3MrDBOKmZmVhgnFTMzK4yTipmZFcZJxczMCuOkYmZmhXFSMTOzwpR6n4pZGXaUpVfMdkQ+UzEzs8I0TSpp+fgzJO0xEAGZmVn7ynOmMhPYG1guaZGko/wIXzMzq6VpUomI7oj4NLAvcB2wAPiVpH+StGfZAZqZWfvINaci6c3ARcAXgW+TLVP/DHBreaGZmVm7aXr1l6QVwO/IlqmfGxEvpl3LJB1WZnBmZtZe8lxSfHJErK4skDQ5Ih6NiBNLisvMzNpQnuGvxTnLzMxsiKt7piLpz4H9gdGSKs9IRgE7lx2YmZm1n0bDX/sB7wV2B46rKN8EfKzMoMzMrD3VTSoRcRNwk6S3RcRdAxiTmZm1qUbDX5+KiC8AH5R0avX+iDiz1MjMzKztNBr+ejj97BqIQMzMrP01Gv76fvp5TW+ZpFcBu0bEMwMQm5mZtZk8C0peJ2mUpNcADwGrJP3f8kMzM7N2k+c+lanpzOQEYAkwEfhIqVGZmVlbypNUdpK0E1lSuSki/gBEuWGZmVk7ypNUvgasAV4D/ETS68kWk2xK0nRJqyR1S5pbY/9ISTek/cskTarYd04qXyXpqIry/y1ppaQHJV0vyTdimpkNEnmWvr80IsZFxDGReQx4Z7N2koYBlwNHA1OBUyVNrao2C9gQEfsAFwMXprZTyZ7jsj8wHfiqpGGSxgFnAh0RcQAwLNUzM7NBIM8qxSOB9wOTquqf36TpNKC7dzFKSYuAGWST/b1mAOel7cXAZekBYDOARWlF5Ecldafj/SrF8GpJfwB2AdY264OZmQ2MPMNfN5H9kd8CPFfxamYc8HjF+55UVrNORGwBNgJj6rWNiCeAL5EllyeBjRHxn7U+XNLs9CjkrnXr1uUI18zMtleepe/HR8T0fhy71iOHqyf469WpWS5pD7IEN5nsGS//LunDEfHNbSpHzAfmA3R0dPjCAjOzAZDnTOV/JL2pH8fuASZUvB/PtkNVL9WRNBwYDaxv0PavgEcjYl26Cu07wNv7EZuZmZUgT1I5HFiRrsK6X9IDku7P0W45MEXSZEkjyCbUO6vqdAKnp+2TgFsjIlL5zHR12GRgCnAP2bDXoZJ2SXMvR/LycjJmZtZieYa/ju7PgSNii6Q5wFKyq7QWRMRKSecDXRHRSfaI4oVpIn496UquVO9Gskn9LcAZEbGV7BHGi4F7U/l9pCEuMzNrPWUnBk0qSYcDUyLiG5LGkq3/9Wjp0RWko6Mjurq8LuaOYtLcm1sdwqC1Zt6xrQ7BdiCSVkRER1/a5Fn763PAPwDnpKKdgG0mxs3MzPLMqbwPOJ50GXFErAV2KzMoMzNrT3mSyuY0eR4AabViMzOzbeRJKjdK+hqwu6SPAT8Cvl5uWGZm1o6aXv0VEV+S9G6yRST3A86NiFtKj8zMzNpOnkuKSUnEicTMzBqqm1QkbaLBc1MiYlQpEZmZWdtq9Iz63QDSzYq/BhaSrcn1IXz1l5mZ1ZBnov6oiPhqRGyKiGci4gqypfDNzMxeIc+cylZJHwIWkQ2HnQpsLTUqG5J8p7xZ+8tzpvJB4APAb9Lr5FRmZmb2CnkuKV5D9gwTMzOzhnJdUmxm7aEvQ4hefNLKkGf4y8zMLJc8qxQPG4hAzMys/eU5U+mW9EVJU0uPxszM2lqepPJm4BfAlZLuljRbku+mNzOzbTRNKummx69HxNuBTwGfA56UdI2kfUqP0MzM2kauORVJx0v6LvCvwEXAG4DvA0tKjs/MzNpInkuKHwFuA74YEf9TUb5Y0hHlhGVmZu0oT1I5LSLurCyQdFhE/HdEnFlSXGZm1obyTNRfWqPsK0UHYmZm7a/R81TeBrwdGCvpkxW7RgG+d8XMzLbRaPhrBLBrqlP5/JRngJPKDMrMzNpTo4d03QHcIenqiHisPweXNJ3sirFhwJURMa9q/0jgWuAQ4GnglLSAJZLOAWaRLbN/ZkQsTeW7A1cCB5Atxf/RiLirP/GZmVmxGg1/XRIRZwGXSdrmscIRcXyjA6flXS4H3g30AMsldUbEQxXVZgEbImIfSTOBC4FT0t37M4H9gb2BH0naNyK2kiWp/4iIkySNAHbpS4fNzKw8jYa/FqafX+rnsacB3RGxGkDSIrIl9CuTygzgvLS9mCyBKZUviogXgUcldQPTJK0EjgD+GiAiNgOb+xmfmZkVrNHw14r0845+Hnsc8HjF+x7grfXqRMQWSRuBMan87qq244DngXXANyQdCKwAPhERz1V/uKTZwGyAiRMn9rMLZmbWF3UvKZb0gKT7671yHFs1yqqH0erVqVc+HDgYuCIiDgKeA+bW+vCImB8RHRHRMXbs2BzhmpnZ9mo0/PXe7Tx2DzCh4v14YG2dOj2ShgOjgfUN2vYAPRGxLJUvpk5SMTOzgddo+KtfV3xVWA5MkTQZeIJs4r362fadwOnAXWSXKd8aESGpE7hO0pfJJuqnAPdExFZJj0vaLyJWAUfyyjkaMzNroUZXf90ZEYdL2sTLQ1Iv/YyIhsvfpzmSOcBSskuKF0TESknnA10R0QlcBSxME/HryRIPqd6NZAljC3BGuvIL4O+Bb6Urv1YDf9PfzpuZWbEanakcnn7uVq9OMxGxhKqVjCPi3IrtF4CT67S9ALigRvlPgY7+xmRmZuXJs6Akkg4GDic7U7kzIu4rNSozM2tLeZ6nci5wDdmlvnsBV0v6TNmBmZlZ+8lzpnIqcFAaqkLSPOBe4PNlBmZmZu0nz9L3a4CdK96PBH5ZSjRmZtbWGl399RWyOZQXgZWSbknv3w3cWa+dmZkNXY2Gv7rSzxXAdyvKby8tGjMza2uNLim+ZiADMTOz9td0ol7SFOBfgKlUzK1ExBtKjMvMzNpQnon6bwBXkN3Z/k6yh2otbNjCzMyGpDxJ5dUR8WNAEfFYRJwHvKvcsMzMrB3luU/lBUmvAh5Ja3k9Aby23LDMzKwd5TlTOYvskb1nkj1L/iNkKwubmZm9QtMzlYhYDpDOVs6MiE2lR2VmZm0pz9VfHWST9bul9xuBj/Y+btjM2tOkuTfnqrdm3rElR2I7kjxzKguAj0fEfwFIOpwsyby5zMDMzKz95JlT2dSbUAAi4k7AQ2BmZraNRmt/HZw275H0NeB6srW/TsFLtZiZWQ2Nhr8uqnr/uYrtKCEWMzNrc43W/nrnQAZiZmbtL8+TH0dL+rKkrvS6SNLogQjOzMzaS96rvx4EPpDef4Ts6q8TywrKdix5L101s/aXJ6n8WUS8v+L9P0n6aVkBmZlZ+8pzSfHz6d4UACQdBjxfXkhmZtau8pyp/C1wbcU8yga89peZmdXQMKmk9b72i4gDJY0CiIhnBiQyMzNrOw2HvyLij8CctP1MXxOKpOmSVknqljS3xv6Rkm5I+5dJmlSx75xUvkrSUVXthkm6T9IP+hKPmZmVK8+cyi2SzpY0QdKeva9mjSQNAy4HjiZ7FPGpkqZWVZsFbIiIfYCLgQtT26nATGB/YDrw1XS8Xp8AHs4Ru5mZDaA8SeWjwBnAT4AV6dWVo900oDsiVkfEZmARMKOqzgzgmrS9GDhSklL5ooh4MSIeBbrT8ZA0HjgWuDJHDGZmNoDyPE9lcj+PPQ54vOJ9D/DWenUiYktaVn9MKr+7qu24tH0J8CnSUvz1SJoNzAaYOHFi/3pgZmZ9kueO+p0lfVLSdyR9W9JZknbOcWzVKKteM6xenZrlkt4LPJXnWS4RMT8iOiKiY+zYsc2jNTOz7ZZn+OtasrmNrwCXkc2PLMzRrgeYUPF+PLC2Xh1Jw4HRwPoGbQ8Djpe0hmw47V2SvpkjFjMzGwB57lPZLyIOrHh/m6Sf5Wi3HJgiaTLwBNnE+wer6nSS3fNyF3AScGtEhKRO4DpJXwb2BqYA90TEXcA5AJLeAZwdER/OEYuZmQ2APEnlPkmHRsTdAJLeCvx3s0ZpjmQOsBQYBiyIiJWSzge6IqITuApYKKmb7AxlZmq7UtKNwEPAFuCMiNjaj/6ZmdkAypNU3gqcJulX6f1E4GFJDwAREXUfKxwRS4AlVWXnVmy/AJxcp+0FwAUNjn07fliYmdmgkiepTC89CjMbtPKuMr1m3rElR2LtIM8lxY8NRCBmZtb+8lz9ZWZmlouTipmZFcZJxczMCpNnot6sJj8m2Myq+UzFzMwK46RiZmaFcVIxM7PCOKmYmVlhnFTMzKwwTipmZlYYJxUzMyuMk4qZmRXGScXMzArjpGJmZoVxUjEzs8I4qZiZWWGcVMzMrDBOKmZmVhgnFTMzK4yTipmZFcZJxczMCuOkYmZmhSk1qUiaLmmVpG5Jc2vsHynphrR/maRJFfvOSeWrJB2VyiZIuk3Sw5JWSvpEmfGbmVnflJZUJA0DLgeOBqYCp0qaWlVtFrAhIvYBLgYuTG2nAjOB/YHpwFfT8bYA/yci3ggcCpxR45hmZtYiZZ6pTAO6I2J1RGwGFgEzqurMAK5J24uBIyUplS+KiBcj4lGgG5gWEU9GxL0AEbEJeBgYV2IfzMysD8pMKuOAxyve97BtAnipTkRsATYCY/K0TUNlBwHLan24pNmSuiR1rVu3rt+dMDOz/IaXeGzVKIucdRq2lbQr8G3grIh4ptaHR8R8YD5AR0dH9eeaWcEmzb05V701844tORJrpTKTSg8woeL9eGBtnTo9koYDo4H1jdpK2oksoXwrIr5TTuhDW94/DmZm1coc/loOTJE0WdIIson3zqo6ncDpafsk4NaIiFQ+M10dNhmYAtyT5luuAh6OiC+XGLuZmfVDaWcqEbFF0hxgKTAMWBARKyWdD3RFRCdZglgoqZvsDGVmartS0o3AQ2RXfJ0REVslHQ58BHhA0k/TR/1jRCwpqx9mZpZfmcNfpD/2S6rKzq3YfgE4uU7bC4ALqsrupPZ8i5mZDQK+o97MzArjpGJmZoUpdfjLzKyaLz3esflMxczMCuOkYmZmhXFSMTOzwjipmJlZYZxUzMysME4qZmZWGCcVMzMrjJOKmZkVxknFzMwK46RiZmaFcVIxM7PCeO2vIcRPdDSzsjmpmNmg5IUn25OHv8zMrDBOKmZmVhgnFTMzK4yTipmZFcZJxczMCuOrv3YAvlTYzAYLJxUza2u+9Hhw8fCXmZkVxmcqg5SHtMysHZWaVCRNB/4VGAZcGRHzqvaPBK4FDgGeBk6JiDVp3znALGArcGZELM1zTDOzWvryDzUPlfVfaUlF0jDgcuDdQA+wXFJnRDxUUW0WsCEi9pE0E7gQOEXSVGAmsD+wN/AjSfumNs2OOaj5DMRs8PM8Tf+VeaYyDeiOiNUAkhYBM4DKBDADOC9tLwYuk6RUvigiXgQeldSdjkeOYxbKScDM6nHy2VaZSWUc8HjF+x7grfXqRMQWSRuBMan87qq249J2s2MCIGk2MDu9fVbSqn70oSx7Ab9tdRADZKj0daj0E4ZOXwvrpy4s4iilqtfX1/f1QGUmFdUoi5x16pXXulqt+phZYcR8YH6jAFtFUldEdLQ6joEwVPo6VPoJQ6evQ6WfUGxfy7ykuAeYUPF+PLC2Xh1Jw4HRwPoGbfMc08zMWqTMpLIcmCJpsqQRZBPvnVV1OoHT0/ZJwK0REal8pqSRkiYDU4B7ch7TzMxapLThrzRHMgdYSnb574KIWCnpfKArIjqBq4CFaSJ+PVmSINW7kWwCfgtwRkRsBah1zLL6UKJBOSxXkqHS16HSTxg6fR0q/YQC+6rsxMDMzGz7eZkWMzMrjJOKmZkVxkmlBJI+IelBSSslnZXK9pR0i6RH0s89UrkkXSqpW9L9kg5ubfR9U6ev50l6QtJP0+uYivrnpL6uknRU6yJvTtICSU9JerCirM/fo6TTU/1HJJ1e67NaqY/9fIekjRXf7bkVbaan77Vb0txW9KWZOn09Of3+/lFSR1X9mr+vg72vfemnpEmSnq/4Tv+tYt8hkh5I/bw03ZzeWET4VeALOAB4ENiF7EKIH5FdvfYFYG6qMxe4MG0fA/yQ7N6cQ4Flre5DAX09Dzi7Rv2pwM+AkcBk4JfAsFb3o0H/jgAOBh6sKOvT9wjsCaxOP/dI23u0um/b0c93AD+ocYxh6ft8AzAifc9TW923nH19I7AfcDvQ0ez3tR362sd+TqqsV3Wce4C3pd/rHwJHN/tsn6kU743A3RHx+4jYAtwBvI9sOZlrUp1rgBPS9gzg2sjcDewu6XUDHXQ/1etrPS8tvxMRjwKVy+8MOhHxE7KrEiv19Xs8CrglItZHxAbgFmB6+dHn18d+1vPSskwRsRnoXUJpUKnV14h4OCJqrbhR7/d10Pe1j/2sKf3+joqIuyLLMNfS/PfASaUEDwJHSBojaReyf8FOAP4kIp4ESD9fm+rXWs5mHO2hXl8B5qRhoAW9Qye0d1979fV7bNc+1+snwNsk/UzSDyXtn8ratZ+N7GjfaSOTJd0n6Q5Jf5nKxpH1rVeufjqpFCwiHiZbbfkW4D/ITo23NGiSZzmbQalBX68A/gx4C/AkcFFq0rZ9zaGvSw61q3uB10fEgcBXgO+l8h2tnzB0vtMngYkRcRDwSeA6SaPoZz+dVEoQEVdFxMERcQTZKegjwG96h7XSz6dS9bZeeqZWXyPiNxGxNSL+CHydl4e42rqvSV+/x3btc81+RsQzEfFs2l4C7CRpL9q3n43saN9pTWl47+m0vYJsvmhfsn6Or6iaq59OKiWQ9Nr0cyJwInA9r1yS5nTgprTdCZyWrh46FNjYO+zQDmr1tWpO6H1kw2RQf/mddtLX73Ep8B5Je6RhwPekssGuZj8l/WnvFUCSppH9DXmaHXMJpSGxXJSkscqef4WkN5D1c3X6/d0k6dD0nZ/Gy7/v9bX6KoUd8QX8F9kSMz8DjkxlY4Afk521/BjYM5WL7MFjvwQeoOKqjHZ41enrwtSX+8n+Z3tdRf1Pp76uIseVJC3u2/VkQwN/IPtX26z+fI/AR8kmebuBv2l1v7azn3OAlen7vht4e8VxjgF+kf4bfLrV/epDX9+Xtl8EfgMsbfb7Otj72pd+Au+v+E7vBY6rOE4H2T8KfwlcRlqFpdHLy7SYmVlhPPxlZmaFcVIxM7PCOKmYmVlhnFTMzKwwTipmZlYYJxUzMyuMk4qZmRXGScVsEJB0laRjWx2H2fZyUjEbHN5CdkezWVsb3uoAzIYiSfsCC4DR6eefRkRP41Zmg5/PVMwGmKThwDeBT0bEm8gW8Pt5a6MyK4bPVMwG3onAwxHRu0LzSuD5FsZjVhifqZgNvDcDKyreH4LnU2wH4aRiNvCeBg4AkHQIcCrw05ZGZFYQL31vNsDSkxKXACPIntNxBNnjXP/Q0sDMCuCkYmZmhfHwl5mZFcZJxczMCuOkYmZmhXFSMTOzwjipmJlZYZxUzMysME4qZmZWmP8Pj+ylzivS+NkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The distribution of d\n",
    "from scipy.special import factorial\n",
    "count, bins, ignored = plt.hist(d, 30, density=True)\n",
    "plt.title(\"d $\\sim$ Poisson$(1000)$\")\n",
    "plt.xlabel('$d$')\n",
    "plt.ylabel('probability density')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({0: 8508, 1: 1492})"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The distribution of x1\n",
    "Counter(x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEaCAYAAAAcz1CnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAaWUlEQVR4nO3df5hdVX3v8fcHCKASwJrYRxPiYBusKaDYEVFoDUVLCBqqRUtEqErNtYo8VqvilQuK3iuCVC6KQNAYQPnlTyKNpV5FECSaiShCuNEYApmiZUCE+BOBz/1j73CPk5kzeyazz8mc/Xk9T545e++19/4uEuZ71lp7ryXbREREc+3Q7QAiIqK7kggiIhouiSAiouGSCCIiGi6JICKi4ZIIIiIaLokgIqLhkggixiBpo6SXdPB+H5L0tprv8V1Jf17nPWLqSCKIqNF4k4ikmcDxwIXl9i6SPiXpLkmbJd0i6YgxrlHlnI8Ap4+3PtGbkggiti+vA1ba/k25vROwCXgxsAfwP4CrJPW1uUaVc1YAh0p62uSFHlNVEkFslySdKelLLdtnSfq6pGldCun5ktZKekDSpyXt2hLb0yV9QdKQpDslnVTuvxSYA3xF0i8lvavcf7Kkn5Tf1tdKekXLfY4Art+yYftXtt9ne6Ptx2xfA9wJ/MVogVY5x/ZvgTXA30zGf5yY2pIIYnv1YYpvrM+V9CZgAfBK27/vUjzHAocDfwLsA5wCIGkH4CvAD4BZwGHA2yQdbvs44G7g5bZ3s31mea2fAH9J8W39/cBnWr6Z7wesGy0ISX9c3v/2qoG3OecO4DlVrxO9K4kgtku27wfOAS4B3gMstP2gpBdKulnS9ZIur9JCkDSrbGGslHSOpL+UtJOkZ0t6b8WQPm57k+2fA/8TWFzufz4w0/bpth+2vQG4CDimTd0+Z/ue8tv6lcCPgQPLw3sCm0epxzTgs8DFtv9vlaDHOGdzeb9ouCSC2J7dQvEN+T22N5X77gL+2vaLgQ3AURWu8wGKb+ZnAz8DPgo8CFxB9W/Wm1o+3wU8vfz8DODpkn6x5Q/w34E/Hu1Cko6X9P2W8vsCM8rDDwDTRzhnB+BS4GHgxCoBVzhnOvCLKteK3rZTtwOIGImk/YDzgYuBNwCXAdi+p6XYI8BjFS73RtuPlp+/DpwxgZD2avk8B9gSxybgTttzRznvD+Z5l/QMihbDYcDNth+V9H1AZZFbKbpxVrecI+BTFMllYZXusYrnPBv4zFjXit6XFkFsdyTNouh3fxPwZmA/SfOHldmbYmD1mrGu15IEtsVbJM2W9EcU3/ivLPd/F3hI0rslPUHSjpL2lfT88vh/Ac9suc6TKJLDUFmP11O0CLZYSfG0T6vzKX5pv7zlaSLK85dLWj5CvKOeU563C8Xg8dfGqHc0QBJBbFck7U7xy/Bfba+w/WvgLIp++dYyFwPH2X64Q6FdBvwHRXfUBuCD8HiSeTnwXIonc+4DPkkxEAzwIeCUshvoX2yvpeiiupkiSewH3NRyn0uAhZKeAI+3IP5bef2flU8f/VLSsWX5vYadX+UcgEXAN4e1sKKhlBXKYiqRtBNwNXC27W90O546SPpfwL22zxmj3M4UTyvtP96nqSR9BzjB9m0TjzR6RRJBTCmSjqMY7N3yC+z88smbiJigJIKIiIbLGEFERMMlEURENFwSQUREwyURREQ0XG1vFktaBryM4jG4fUcpM59iPplpwH3ltAFtzZgxw319fZMYaURE71uzZs19tmeOdKzOKSaWAx+neEFmK5L2BD4BLLB9t6SnVrloX18fAwMDkxZkREQTSLprtGO1dQ3ZvgH4eZsirwG+aPvusvy9dcUSERGj6+YYwT7AkyV9U9IaScd3MZaIiMbq5uyjO1FMenUY8ATgZkmrbP9oeEFJS4AlAHPmzOlokBERva6bLYJB4N/LZfXuA25glNWSbC+13W+7f+bMEcc6IiJigrqZCK4GtqwU9UTgBRRL50VERAfV+fjo5cB8YIakQeA0isdEsX2B7Tsk/TvFQhyPAZ/MTIgREZ1XWyKwvbhCmbMo5pqPiIguyZvFERENlzWLI8bQd/K/VS678Ywja4wkoh5pEURENFwSQUREwyURREQ0XMYIIiZR1fGEjCXE9iQtgoiIhkuLIKIL0nKI7UlaBBERDZdEEBHRcEkEERENl0QQEdFwSQQREQ2XRBAR0XBJBBERDZdEEBHRcEkEERENl0QQEdFwtSUCScsk3Sup7TrEkp4v6VFJR9cVS0REjK7OFsFyYEG7ApJ2BD4MXFtjHBER0Uadi9ffIKlvjGJvBb4APL+uOKJ3ZKK2iHp0bYxA0izgFcAFFcoukTQgaWBoaKj+4CIiGqSbg8XnAO+2/ehYBW0vtd1vu3/mzJkdCC0iojm6uR5BP3CFJIAZwEJJj9j+chdjih7QxC6kJtY5Jk/XEoHtvbd8lrQcuCZJIDqp6i/PiF5XWyKQdDkwH5ghaRA4DZgGYHvMcYGIiOiMOp8aWjyOsq+rK46IqSytluiEvFkcEdFwWbw+ui7fejsng8oxkrQIIiIaLokgIqLhkggiIhouiSAiouEyWBwREzaegf4MQG+/kgiiNnkaKGJqSNdQRETDJRFERDRcEkFERMMlEURENFwGiyNiKxnob5a0CCIiGi6JICKi4ZIIIiIaLokgIqLhkggiIhqutkQgaZmkeyXdNsrxYyXdWv75tqTn1BVLRESMrs4WwXJgQZvjdwIvtr0/8AFgaY2xRETEKOpcvP4GSX1tjn+7ZXMVMLuuWCIiYnTbyxjBCcBXRzsoaYmkAUkDQ0NDHQwrIqL3df3NYkmHUiSCQ0YrY3spZddRf3+/OxRajCJvnUb0lq4mAkn7A58EjrB9fzdjiYhoqq51DUmaA3wROM72j7oVR0RE09XWIpB0OTAfmCFpEDgNmAZg+wLgVOApwCckATxiu7+ueCIiYmR1PjW0eIzj/wj8Y133j4iIaraXp4YiIqJLxkwE5WObb5H05E4EFBERnVWlRXAM8HRgtaQrJB2uslM/IiKmvjETge31tt8L7ANcBiwD7pb0fkl/VHeAERFRr0qDxeXz/q8HFgJfAD5L8QLYN4Dn1hZdRPSMqi8ibjzjyJojieHGTASS1gC/AD4FnGz7d+Wh70g6uM7gIiKiflVaBK+yvaF1h6S9bd9p+5U1xRURER1SZbD48xX3RUTEFDRqi0DSnwF/DuwhqfWb/+7ArnUHFpMj/bIRMZZ2XUPPAl4G7Am8vGX/ZuCNdQYVnZcZRSOaa9REYPtq4GpJL7R9cwdjiogGSyu289p1Db3L9pnAayRtNW+Q7ZNqjSwiIjqiXdfQHeXPgU4EEhER3dGua+gr5c+Lt+yTtAOwm+2HOhBbRER0QJVJ5y6TtLukJwFrgXWS3ll/aBER0QlV3iOYV7YA/hZYCcwBjqs1qoiI6JgqiWCapGkUieBq278HsoB8RESPqJIILgQ2Ak8CbpD0DGDMMQJJyyTdK+m2UY5L0rmS1ku6VdLzxhN4RERMjirTUJ9re5bthS7cBRxa4drLgQVtjh8BzC3/LAHOr3DNiIiYZFVmH90F+Dugb1j509udZ/sGSX1tihwFXGLbwCpJe0p6mu2fjhVTRERMniqzj14NPAisAX43RtnxmAVsatkeLPclEUREdFCVRDDbdrsunokaabnLEQehJS2h6D5izpw5NYQSEdFcVQaLvy1pvxruPQjs1bI9G7hnpIK2l9rut90/c+bMGkKJiGiuKongEGCNpHXl0z0/lHTrJNx7BXB8+fTQQcCDGR+IiOi8Kl1DR0zkwpIuB+YDMyQNAqcB0wBsX0DxctpCYD3wa4o1kSMiosPGTAS275J0CDDX9qclzQR2q3DeVjOWDjtu4C2VI42IiFpUmWvoNODdwHvKXdOAz9QZVEREdE6VrqFXAAcA3wOwfY+k6bVGFRExhixgM3mqDBY/XHbjGKCchTQiInpElURwlaQLgT0lvRH4P8BF9YYVERGdUmWw+COSXkox0dyzgFNtf632yKKtLDYfEZOlyhgB5S/+/PKPiOhB7Rav30ybdQds715LRBER0VHt1iyeDiDpdOBnwKUU8wMdC+SpoYiIHlFlsPhw25+wvdn2Q7bPp5iWOiIiekCVRPCopGMl7ShpB0nHAo/WHVhERHRGlUTwGuDVwH+Vf15V7ouIiB5Q5fHRjRSriUVERA+q0iKIiIgelkQQEdFwVWYf3bETgURERHdUaRGsl3SWpHm1RxMRER1XJRHsD/wI+KSkVZKWSMpbxRERPWLMRFC+SHaR7RcB76JYcvKnki6W9Ke1RxgREbWqNEYgaZGkLwH/GzgbeCbwFYp1h9udu6Bc9H69pJNHOD5H0nWSbpF0q6SFE6xHRERMUJXZR38MXAecZfvbLfs/L+mvRjupHGQ+D3gpMAislrTC9tqWYqcAV9k+vxyDWAn0jbMOPSXTS0dEp1VJBMfbvrF1h6SDbd9k+6Q25x0IrLe9oTznCooX01oTgYEt4w17APdUjjwiIiZFlcHic0fY97EK580CNrVsD5b7Wr0PeK2kQYrWwFtHulA5QD0gaWBoaKjCrSMioqp26xG8EHgRMFPS21sO7Q5UebdAI+wbvr7BYmC57bPL+10qaV/bj/3BSfZSYClAf3//qGskRETE+LXrGtoZ2K0s07r+wEPA0RWuPQjs1bI9m627fk4AFgDYvlnSrsAM4N4K14+IiEnQbmGa64HrJS23fdcErr0amCtpb+A/gWPYetbSu4HDgOWSng3sCqTvJyKig9p1DZ1j+23AxyVt1R1je1G7C9t+RNKJwLUUXUnLbN9erng2YHsF8A7gIkn/TNFt9Drb6fqJiElT9Um8jWccWXMk2692XUOXlj8/MtGL217JsHcNbJ/a8nktcPBErx8REduuXdfQmvLn9Z0LJyIiOq1d19AP2fopn8fZ3r+WiCIioqPadQ29rGNRRERE17TrGprIk0IRETHFjPpmsaQby5+bJT00/GfnQoyIiDq1axEcUv6cPlqZiIiY+qpMOoek5wGHUAwe32j7llqjioiIjqmyHsGpwMXAUyimf1gu6ZS6A4uIiM6o0iJYDBxg+7cAks4Avgd8sM7AIiKiM6pMQ72RYg6gLXYBflJLNBER0XHtXij7GMWYwO+A2yV9rdx+KXDjaOdFRMTU0q5raKD8uQb4Usv+b9YWTUREdFy7x0cv7mQgERHRHWMOFkuaC3wImEfLWIHtZ9YYV0REdEiVweJPA+cDjwCHApfw/6eojoiIKa5KIniC7a8Dsn2X7fcBf11vWBER0SlV3iP4raQdgB+XK479J/DUesPqPVVXSYqI6LQqLYK3AU8ETgL+AjgO+Ic6g4qIiM4ZMxHYXm37l8BDwEm2X2l7VZWLS1ogaZ2k9ZJOHqXMqyWtlXS7pMvGF35ERGyrKk8N9VMMGE8vtx8E3rBlKcs25+0InEfxAtogsFrSinKd4i1l5gLvAQ62/YCkdDlFRFc0eZH7Kl1Dy4A32+6z3Qe8hSIxjOVAYL3tDbYfBq4AjhpW5o3AebYfALB9b+XIIyJiUlRJBJttf2vLhu0bgc0VzpsFbGrZHiz3tdoH2EfSTZJWSVow0oUkLZE0IGlgaGiowq0jIqKqdnMNPa/8+F1JFwKXU8w19PdUm2ZCI+zzCPefC8wHZgPfkrSv7V/8wUn2UmApQH9///BrRETENmg3RnD2sO3TWj5X+WU8COzVsj0buGeEMqts/x64U9I6isSwusL1IyJiErSba+jQbbz2amCupL0p3j04BnjNsDJfpljvYLmkGRRdRRu28b4RETEOVVYo20PSv27po5d0tqQ9xjrP9iPAicC1wB3AVbZvl3S6pEVlsWuB+yWtBa4D3mn7/olXJyIixqvKm8XLgNuAV5fbx1E8NfTKsU60vRJYOWzfqS2fDby9/BMREV1QJRH8ie2/a9l+v6Tv1xVQRER0VpXHR38j6ZAtG5IOBn5TX0gREdFJVVoEbwIuaRkXeIDMNRQR0TPaJoJy1tFn2X6OpN0BbD/UkcgiIqIj2nYN2X6M4skfbD+UJBAR0XuqdA19TdK/AFcCv9qy0/bPa4sqImI71YuT01VJBG8of76lZZ+BrFkcEdEDxkwEtvfuRCAREdEdVdYj2BV4M3AIRUvgW8AFtn9bc2wREdEBVbqGLqGYdvpj5fZi4FLgVXUFFRERnVMlETzL9nNatq+T9IO6Appqsih9REx1Vd4svkXSQVs2JL0AuKm+kCIiopOqtAheABwv6e5yew5wh6QfUswbt39t0UVERO2qJIIRl4+MiIjeUOXx0bs6EUhERHRHlTGCiIjoYUkEERENV2sikLRA0jpJ6yWd3Kbc0ZIsqb/OeCIiYmu1JQJJOwLnAUcA84DFkuaNUG46cBLwnbpiiYiI0dXZIjgQWG97g+2HgSuAo0Yo9wHgTCBTVkREdEGdiWAWsKlle7Dc9zhJBwB72b6mxjgiIqKNOhOBRtjnxw8Wq599FHjHmBeSlkgakDQwNDQ0iSFGRESdiWAQ2KtlezZwT8v2dGBf4JuSNgIHAStGGjC2vdR2v+3+mTNn1hhyRETz1JkIVgNzJe0taWfgGGDFloO2H7Q9w3af7T5gFbDI9kCNMUVExDC1JQLbj1Csd3wtcAdwle3bJZ0uaVFd942IiPGpMtfQhNleCawctu/UUcrOrzOWiIhOGs8U9d1e3zhvFkdENFwSQUREwyURREQ0XK1jBFNZlqCMiKZIiyAiouGSCCIiGi6JICKi4ZIIIiIaLokgIqLhkggiIhouiSAiouGSCCIiGi6JICKi4ZIIIiIaLokgIqLhkggiIhouiSAiouGSCCIiGq7WRCBpgaR1ktZLOnmE42+XtFbSrZK+LukZdcYTERFbqy0RSNoROA84ApgHLJY0b1ixW4B+2/sDnwfOrCueiIgYWZ0tggOB9bY32H4YuAI4qrWA7ets/7rcXAXMrjGeiIgYQZ2JYBawqWV7sNw3mhOAr450QNISSQOSBoaGhiYxxIiIqDMRaIR9HrGg9FqgHzhrpOO2l9rut90/c+bMSQwxIiLqXLN4ENirZXs2cM/wQpJeArwXeLHt39UYT0REjKDOFsFqYK6kvSXtDBwDrGgtIOkA4EJgke17a4wlIiJGUVsisP0IcCJwLXAHcJXt2yWdLmlRWewsYDfgc5K+L2nFKJeLiIia1Nk1hO2VwMph+05t+fySOu8fERFjy5vFERENl0QQEdFwtXYNRUTE2PpO/rdK5TaecWQt90+LICKi4ZIIIiIaLokgIqLhkggiIhouiSAiouEa9dRQ1ZH5iIgmSYsgIqLhkggiIhouiSAiouGSCCIiGi6JICKi4ZIIIiIaLokgIqLhkggiIhouiSAiouFqTQSSFkhaJ2m9pJNHOL6LpCvL49+R1FdnPBERsbXaEoGkHYHzgCOAecBiSfOGFTsBeMD2nwIfBT5cVzwRETGyOlsEBwLrbW+w/TBwBXDUsDJHAReXnz8PHCZJNcYUERHD1Dnp3CxgU8v2IPCC0crYfkTSg8BTgPtaC0laAiwpN38pad0EY5ox/NoN0cR6N7HO0Mx6N6bO+sM+k/HW+xmjHagzEYz0zd4TKIPtpcDSbQ5IGrDdv63XmWqaWO8m1hmaWe8m1hkmt951dg0NAnu1bM8G7hmtjKSdgD2An9cYU0REDFNnIlgNzJW0t6SdgWOAFcPKrAD+ofx8NPAN21u1CCIioj61dQ2Vff4nAtcCOwLLbN8u6XRgwPYK4FPApZLWU7QEjqkrntI2dy9NUU2sdxPrDM2sdxPrDJNYb+ULeEREs+XN4oiIhksiiIhouJ5MBE2d2qJCvd8uaa2kWyV9XdKozxVPFWPVuaXc0ZIsqSceM6xSb0mvLv++b5d0WadjnGwV/n3PkXSdpFvKf+MLuxHnZJK0TNK9km4b5bgknVv+N7lV0vMmdCPbPfWHYmD6J8AzgZ2BHwDzhpV5M3BB+fkY4Mpux92heh8KPLH8/E9Tvd5V6lyWmw7cAKwC+rsdd4f+rucCtwBPLref2u24O1DnpcA/lZ/nARu7Hfck1PuvgOcBt41yfCHwVYp3sg4CvjOR+/Rii6CpU1uMWW/b19n+dbm5iuLdjqmsyt81wAeAM4HfdjK4GlWp9xuB82w/AGD73g7HONmq1NnA7uXnPdj6vaUpx/YNtH+36ijgEhdWAXtKetp479OLiWCkqS1mjVbG9iPAlqktprIq9W51AsU3ialszDpLOgDYy/Y1nQysZlX+rvcB9pF0k6RVkhZ0LLp6VKnz+4DXShoEVgJv7UxoXTXe/+9HVOcUE90yaVNbTDGV6yTptUA/8OJaI6pf2zpL2oFiVtvXdSqgDqnyd70TRffQfIqW37ck7Wv7FzXHVpcqdV4MLLd9tqQXUryjtK/tx+oPr2sm5XdZL7YImjq1RZV6I+klwHuBRbZ/16HY6jJWnacD+wLflLSRog91RQ8MGFf9N3617d/bvhNYR5EYpqoqdT4BuArA9s3ArhQTs/WySv/fj6UXE0FTp7YYs95lN8mFFElgqvcZwxh1tv2g7Rm2+2z3UYyLLLI90J1wJ02Vf+Nfpng4AEkzKLqKNnQ0yslVpc53A4cBSHo2RSIY6miUnbcCOL58eugg4EHbPx3vRXqua8jb59QWtatY77OA3YDPlWPjd9te1LWgt1HFOvecivW+FvgbSWuBR4F32r6/e1Fvm4p1fgdwkaR/pugeed1U/4In6XKK7r0Z5djHacA0ANsXUIyFLATWA78GXj+h+0zx/04REbGNerFrKCIixiGJICKi4ZIIIiIaLokgIqLhkggiIhouiSAiouGSCCK2QTnt8UvLzx+UdG63Y4oYr557oSyiw04DTpf0VOAAYMq+oBfNlRfKIraRpOsp3tieb3uzpL8FjgSeSjEV9H90NcCIMSQRRGwDSfsBXwDus/2iYceeDHzE9gldCS6ioowRRExQuQDIZykWB/mVpMOHFTkFOK/jgUWMUxJBxARIeiLwReAdtu+gWAXtfeUxSfow8FXb3+telBHVpGsoYpJJOolimvPVwPfLWSIjtltJBBERDZeuoYiIhksiiIhouCSCiIiGSyKIiGi4JIKIiIZLIoiIaLgkgoiIhksiiIhouCSCiIiG+383XRiUK/DGNQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The distribution of x2\n",
    "count, bins, ignored = plt.hist(x2, 30, density=True)\n",
    "plt.title(\"$x_2$ $\\sim$ beta$(2,2)$\")\n",
    "plt.xlabel('$x_2$')\n",
    "plt.ylabel('probability density')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEaCAYAAAAL7cBuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAeDElEQVR4nO3de5xdVX338c/XcPPhJkrqo7mYaKOPURF0GrR4FzQUTbRiDd6wUlMrKfqyVsOjBY2PbQS13lChkop4iVRrnWp8IVZRaQUzEQQDRkIMZowKCgIiggnf54+9Rzcne87sGWbPOZl836/XeZ2z915rn98JzPmdtdbea8k2ERERne7T6wAiIqI/JUFEREStJIiIiKiVBBEREbWSICIiolYSRERE1EqCiIiIWkkQEbELSU+U9G1J35D0aUl79zqmmHpJEBFR53rgGbafCmwBlvY4nuiBJIiI3YCkrZKOrmz/k6TXtfV+trfbvqPc3AHcXXnv70h6VFvvHf0jCSL2eOWX788l7V/Z91eSLu5hWKOSNBN4OXB2Zd/9JX1e0u2Srpf04jHOsa+kc8uyt0m6XNKxNeXmA8cCX6zsfhewanI+TfSzJIiIwl7AayfjRJL2mozzdPEKYF3lFz7AWcBdwAOBlwAfHuNX/l7ANuCpwMHAPwAXSJo3UkDSQcB5wMts31WpOwg8XdKD7vUnib6WBBF9R9LB5cDoL8pft1dIavv/1TOBN0i63ygxPVLSxZJ+JWmjpCUdx7dKepOkK4HbJe1V7vt7SVeWv+zPlfRASV8uP9dXJR1S1l8p6bpy/9WSnt8l1mOBb1Tee3/gBcA/2P617UsovsRfNtoJbN9u+622t9q+2/YXgR8Bjy/PuRfwaeCttjd11P0tsAF4VpcYYxpIgoh+9E5gJ/AQil+3L7V9d/cq99oQcDHwhs4D5RU8/wl8Bfgj4G+BT0p6REfRE4DjgPvZ3lHuewFwDPBw4LnAl4H/CxxK8fd3SlnuOuDJFJ/3bcAnuvxCfwxQ/dJ+OLDT9g8r+74HNB4nkPTA8jwbK5/lSOC0MjG+qKPKNcBjm54/dk9JENGPfgdcC/ym/HX7/W6FJc2SdIakdZLeK+nJ5S/4R0p68zje9zTgb8s+/qonAAcAq23fZftrFH3yJ3SUe7/tbR1dPx+w/XPbPwG+BVxm+3LbdwKfB44AsP1v5cDw3bY/U37+RaPEeT/gtsr2AcAtHWVuAQ5s8qHLBPhJ4DzbPyjjOd/2obafVj4+01HttjKOmMaSIKIf/YDil/ztkv66Qfm3Az8G3g38DPhnii/ItfzhF/GYykT0RWBlx6EHA9s6WjHXA7M6ym2rOe3PK6/vqNk+AEDSy8uutF9J+hXwaIpWRp2bueeX/6+BgzrKHMQ9k0itsuvufIrxixVjla84EPjVOMrHbqjtwbSIcSn79lcAj7N9bcNqr7K9s3z9X8DqexHC6cB3KZLNiO3AHEn3qSSJucAPO+pOaPUtSQ8B/gV4JvBt2zslXQFolCpXUnQHrS+3fwjsJWlB5d/ssYyRHCUJOJdiYPvPbP9uHGE/EvjEOMrHbigtiOg3jwZ+QvlLW9LckYHc0VSSw71mezPwGf4wNgBwGXA78EZJe0t6GsV4wtpJetv9KZLLjQCS/pLi32E06yiuPhqJ+Xbg34FVkvaXdBTFjW3nl+f7mKSP1ZznwxRf9M/t6BbrStK+FIPZFzWtE7unJIjoN/9K0d2xTdItFF98U93SXUXxpQ1AeYnnEoqrh34BfAh4+Uh//b1l+2qKFsu3KRLjY4D/7lLl48CfSbpvZd9rgPsCN1BcffQ3tkdaEHM6z1e2Wv4aOBz4maRfl4+XNAh5CXCx7e0NysZuTFmTOmL3I+kfgRtsv3eMcvtQXNF02Di7kLqd8zLgpLEuHojdXxJERETUShdTRETUSoKIiIhaSRAREVErCSIiImpNmxvlDj30UM+bN6/XYURE7FY2bNjwC9ud08sA0yhBzJs3j6GhoV6HERGxW5F0/WjH0sUUERG1kiAiIqJWEkRERNRKgoiIiFqtJghJiyVtkrRZUucc+9Vyx0uypIHKvlPLepskPbvNOCMiYletXcUkaQbFQurHAMPAekmD5cyV1XIHUkytfFll30JgGcWSiQ8Gvirp4ZM5rXNERHTXZgtiEbDZ9pZyuuS1FHPUd3o7cAbw28q+pcBa23fa/hGwmdGXX4yIiBa0mSBmcc8lGIfpWKJR0hHAHNtfHG/dsv5ySUOShm688cbJiToiIoB2E0Tdcom/n1u8XAv3n4G/G2/d3++wz7E9YHtg5szaGwEjImKC2ryTephiJasRsynW9h1xIMWyihcXS+Pyv4HBck3isepG9IV5K7/UqNzW1ce1HEnE5GuzBbEeWCBpfrmq1TJgcOSg7VtsH2p7nu15wKXAEttDZbllkvaVNB9YAHynxVgjIqJDay0I2zskrQAuBGYAa2xvlLQKGLI92KXuRkkXAFcDO4CTcwVTRMTUanWyPtvrgHUd+04bpezTOrbfAbyjteAiIqKraTOba8Rkajq2EDGdZaqNiIiolQQRERG1kiAiIqJWEkRERNRKgoiIiFpJEBERUSsJIiIiaiVBRERErSSIiIiolQQRERG1kiAiIqJWEkRERNRKgoiIiFpJEBERUSsJIiIiaiVBRERErVYThKTFkjZJ2ixpZc3xV0u6StIVki6RtLDcP0/SHeX+KyR9pM04IyJiV62tKCdpBnAWcAwwDKyXNGj76kqxT9n+SFl+CfAeYHF57Drbh7cVX+x5erlKXNP33rr6uJYjiWiuzRbEImCz7S227wLWAkurBWzfWtncH3CL8URExDi0mSBmAdsq28PlvnuQdLKk64AzgFMqh+ZLulzSNyQ9ue4NJC2XNCRp6MYbb5zM2CMi9nhtJgjV7NulhWD7LNsPA94EvKXc/VNgru0jgNcDn5J0UE3dc2wP2B6YOXPmJIYeERFtJohhYE5lezawvUv5tcDzAGzfafuX5esNwHXAw1uKMyIiarSZINYDCyTNl7QPsAwYrBaQtKCyeRxwbbl/ZjnIjaSHAguALS3GGhERHVq7isn2DkkrgAuBGcAa2xslrQKGbA8CKyQdDfwOuBk4saz+FGCVpB3ATuDVtm9qK9bYvfXy6qSI6ay1BAFgex2wrmPfaZXXrx2l3ueAz7UZW0REdJc7qSMiolYSRERE1EqCiIiIWkkQERFRKwkiIiJqJUFEREStJIiIiKiVBBEREbWSICIiolYSRERE1EqCiIiIWkkQERFRKwkiIiJqJUFEREStJIiIiKiVBBEREbWSICIiolarCULSYkmbJG2WtLLm+KslXSXpCkmXSFpYOXZqWW+TpGe3GWdEROyqtQQhaQZwFnAssBA4oZoASp+y/RjbhwNnAO8p6y4ElgGPAhYDHyrPFxERU6TNFsQiYLPtLbbvAtYCS6sFbN9a2dwfcPl6KbDW9p22fwRsLs8XERFTZK8Wzz0L2FbZHgaO7Cwk6WTg9cA+wDMqdS/tqDurpu5yYDnA3LlzJyXoiIgotNmCUM0+77LDPsv2w4A3AW8ZZ91zbA/YHpg5c+a9CjYiIu6pzQQxDMypbM8GtncpvxZ43gTrRkTEJBszQUgaknSypEPGee71wAJJ8yXtQzHoPNhx7gWVzeOAa8vXg8AySftKmg8sAL4zzvePiIh7ockYxDLgL4H1koaAfwW+YnuXLp8q2zskrQAuBGYAa2xvlLQKGLI9CKyQdDTwO+Bm4MSy7kZJFwBXAzuAk23vnNhHjIiIidAY3/N/KCjdB3gO8GHgbmAN8D7bN7UXXnMDAwMeGhrqdRjRA/NWfqnXIUyarauP63UIsYeRtMH2QN2xRmMQkg4D3g2cCXwOOB64FfjaZAUZERH9ZcwuJkkbgF8B5wIrbd9ZHrpM0lFtBhd7tunUMojYHTUZg3ih7S3VHZLm2/6R7T9vKa6IPVLTpJiuqJgKTbqYPttwX0RETCOjtiAk/R+KuZAOllRtKRwE7Nd2YBER0VvdupgeQXHV0v2A51b23wa8qs2gIiKi90ZNELa/AHxB0hNtf3sKY4qIiD7QrYvpjbbPAF4s6YTO47ZPaTWyiIjoqW5dTNeUz7n7LCJiD9Sti+k/y+fzRvaVd1Mf0LGOQ0RETENNJuv7lKSDJO1PMTfSJkl/335oERHRS03ug1hYthieB6wD5gIvazWqiIjouSYJYm9Je1MkiC/Y/h01i/dERMT00iRBnA1spVgz+puSHkIxUV9ERExjY87FZPv9wPsru66X9PT2QoqIiH7QZDbXfYEXAPM6yq9qKaaIGEMm9Yup0GQ21y8AtwAbgDvHKBsREdNEkwQx2/biiZxc0mLgfRRLjn7U9uqO468H/opiWdEbgVfavr48thO4qiz6Y9tLJhJDRERMTJNB6v+R9JjxnljSDOAs4FhgIXCCpIUdxS4HBmwfRjGF+BmVY3fYPrx8JDlEREyxJgniScAGSZskXSnpKklXNqi3CNhse4vtu4C1wNJqAdtft/2bcvNSYPZ4go+IiPY06WI6doLnngVsq2wPA0d2KX8S8OXK9n6Shii6n1bb/o/OCpKWA8sB5s6dO8EwIyKizpgtiHJMYA7wjPL1b5rUA1R3utqC0kuBAeDMyu65tgeAFwPvlfSwmtjOsT1ge2DmzJkNQoqIiKaazMV0OvAm4NRy197AJxqce5gisYyYDWyvOf/RwJuBJbZ/f5WU7e3l8xbgYuCIBu8ZERGTpElL4PnAEuB2+P0X94EN6q0HFkiaL2kfYBkwWC0g6QiKO7WX2L6hsv+Q8v4LJB0KHEUxUWBEREyRJmMQd9m2JAOUs7qOyfYOSSuACykuc11je6OkVcCQ7UGKLqUDgH+TBH+4nPWRwNmS7qZIYqttJ0FEREyhJgniAklnA/eT9CrglcC/NDm57XUUM8BW951WeX30KPX+Bxj3pbWxe2h6F3BE9FaTuZjeJekYign6HgGcZvui1iOLiIieatKCoEwISQoREXuQUROEpNvosu6D7YNaiSgiIvpCtzWpDwQoB5V/BpxPcW/DS2h2FVNEROzGmlzm+mzbH7J9m+1bbX+YYvrviIiYxpokiJ2SXiJphqT7SHoJsLPtwCIioreaJIgXA38B/Lx8vLDcFxER01iTy1y30jELa0RETH9NWhAREbEHSoKIiIhaTWZznTEVgURERH9p0oLYLOnMmuVCIyJiGmuSIA4Dfgh8VNKlkpZLyl3UERHTXJMV5W6z/S+2/xR4I3A68FNJ50n649YjjIiInmg0BiFpiaTPA+8D3g08FPhPOqbyjoiI6aPJbK7XAl8HzizXaRjxWUlPaSesiIjotSYJ4uW2L6nukHSU7f+2fUpLcUVERI81GaR+f82+DzQ5uaTFkjZJ2ixpZc3x10u6WtKVkv5L0kMqx06UdG35OLHJ+0VExOTpth7EE4E/BWZKen3l0EEUa0x3Vd4/cRZwDDAMrJc02LG29OXAgO3fSPob4AzgRZLuTzEYPkCxJsWGsu7N4/t4ERExUd1aEPsAB1AkkQMrj1uB4xucexGw2fYW23cBa+mY08n2123/pty8FJhdvn42cJHtm8qkcBGwuNlHioiIydBtwaBvAN+Q9DHb10/g3LOAbZXtYeDILuVPAr7cpe6szgqSlgPLAebOnTuBECMiYjTdupjea/t1wAcl7bL0qO0lY5xbNftqlzCV9FKK7qSnjqeu7XOAcwAGBgZGXR41IiLGr9tVTOeXz++a4LmHgTmV7dnA9s5Cko4G3gw81fadlbpP66h78QTjiIiICejWxbShfP7GBM+9HlggaT7wE2AZHQsNSToCOBtYbPuGyqELgX+UdEi5/Szg1AnGERERE9Cti+kqRukSArB9WLcT294haQXFl/0MYI3tjZJWAUO2B4EzKQbC/00SwI9tL7F9k6S3UyQZgFW2bxrPB4uIiHunWxfTc+7tyW2vo2M6DtunVV4f3aXuGmDNvY0hIiImplsX00SuXIqIiGli1PsgJF1SPt8m6dbO56kLMSIieqFbC+JJ5fOBUxdORET0iyaT9SHpccCTKAatL7F9eatRRcSkmLfyS43KbV19XMuRxO6oyXoQpwHnAQ8ADgU+JuktbQcWERG91aQFcQJwhO3fAkhaDXwX+H9tBhYREb3VZLrvrcB+le19getaiSYiIvpGtxvlPkAx5nAnsFHSReX2McAlo9WLiIjpoVsX01D5vAH4fGX/xa1FExERfaPbZa7nTWUgERHRX8YcpJa0APgnYCGVsQjbD20xroiI6LEmg9T/CnwY2AE8Hfg4f5gKPCIipqkmCeK+tv8LkO3rbb8VeEa7YUVERK81uQ/it5LuA1xbTt/9E+CP2g0rIiJ6rUkL4nXA/wJOAR4PvAw4sc2gIiKi98ZsQdheD1C2Ik6xfVvrUUVERM81mYtpoFxd7krgKknfk/T49kOLiIheatLFtAZ4je15tucBJ1Nc2TQmSYslbZK0WdLKmuNPkfRdSTskHd9xbKekK8rHYJP3i4iIydNkkPo2298a2bB9iaQxu5kkzQDOopiaYxhYL2nQ9tWVYj8GXgG8oeYUd9g+vEF80SeaTi0dEbuHbnMxPa58+R1JZwOfppiL6UU0m25jEbDZ9pbyfGuBpcDvE4TtreWxuycQe0REtKhbC+LdHdunV167wblnAdsq28PAkQ3jAthP0hDFDXqrbf9HZwFJy4HlAHPnzh3HqSMiYizd5mJ6+r08t+pOO476c21vl/RQ4GuSrrJ9j2nGbZ8DnAMwMDAwnnNHRMQYmlzFdLCk90gaKh/vlnRwg3MPA3Mq27OB7U0Ds729fN5C0aV1RNO6ERFx7zW9iuk24C/Kx600u4ppPbBA0nxJ+wDLgEZXI0k6RNK+5etDgaOojF1ERET7mlzF9DDbL6hsv03SFWNVsr2jnJrjQmAGsMb2RkmrgCHbg5L+hGKtiUOA50p6m+1HAY8Ezi4Hr+9DMQaRBBERMYWaJIg7JD3J9iUAko4C7mhyctvrgHUd+06rvF5P0fXUWe9/gMc0eY+IiGhHkwTxauDjlXGHm8lcTBER017XBFHOv/QI24+VdBCA7VunJLKIiOiproPUtu8GVpSvb01yiIjYczTpYrpI0huAzwC3j+y0fVNrUUXElBrPNClbVx/XYiTRT5okiFeWzydX9hnImtQREdNYk/Ug5k9FIBER0V/GTBCS9gNeAzyJouXwLeAjtn/bcmwREdFDTbqYPk5xJ/UHyu0TgPOBF7YVVERE9F6TBPEI24+tbH9d0vfaCigiIvpDk7mYLpf0hJENSUcC/91eSBER0Q+atCCOBF4u6cfl9lzgmnKdats+rLXoIiKiZ5okiMWtRxEREX2nyWWu109FIBER0V+ajEFERMQeKAkiIiJqJUFEREStJIiIiKjVaoKQtFjSJkmbJa2sOf4USd+VtEPS8R3HTpR0bfnIAkUREVOstQQhaQZwFnAssBA4QdLCjmI/Bl4BfKqj7v2B0ynuwVgEnC7pkLZijYiIXTW5D2KiFgGbbW8BkLQWWApcPVLA9tby2N0ddZ8NXDSy5oSkiyjux/h0i/HGKMazVkBETB9tdjHNArZVtofLfW3XjYiISdBmglDNPk9mXUnLJQ1JGrrxxhvHFVxERHTXZoIYBuZUtmcD2yezru1zbA/YHpg5c+aEA42IiF21mSDWAwskzZe0D7AMGGxY90LgWZIOKQenn1Xui4iIKdJagrC9A1hB8cV+DXCB7Y2SVklaAiDpTyQNUyw+dLakjWXdm4C3UySZ9cCqkQHriIiYGm1exYTtdcC6jn2nVV6vp+g+qqu7BljTZnwRETG63EkdERG1kiAiIqJWEkRERNRqdQwiIqafpnfWb119XMuRRNvSgoiIiFpJEBERUSsJIiIiaiVBRERErQxS78EyjXdEdJMWRERE1EqCiIiIWkkQERFRKwkiIiJqJUFEREStJIiIiKiVBBEREbWSICIiolarCULSYkmbJG2WtLLm+L6SPlMev0zSvHL/PEl3SLqifHykzTgjImJXrd1JLWkGcBZwDDAMrJc0aPvqSrGTgJtt/7GkZcA7gReVx66zfXhb8UVERHdttiAWAZttb7F9F7AWWNpRZilwXvn6s8AzJanFmCIioqE252KaBWyrbA8DR45WxvYOSbcADyiPzZd0OXAr8Bbb32ox1oiYZFlYaPfXZoKoawm4YZmfAnNt/1LS44H/kPQo27feo7K0HFgOMHfu3EkIOSIiRrTZxTQMzKlszwa2j1ZG0l7AwcBNtu+0/UsA2xuA64CHd76B7XNsD9gemDlzZgsfISJiz9VmglgPLJA0X9I+wDJgsKPMIHBi+fp44Gu2LWlmOciNpIcCC4AtLcYaEREdWutiKscUVgAXAjOANbY3SloFDNkeBM4Fzpe0GbiJIokAPAVYJWkHsBN4te2b2oo1IiJ21eqCQbbXAes69p1Wef1b4IU19T4HfK7N2CIiorvcSR0REbWy5Og0k2VEI2KypAURERG1kiAiIqJWEkRERNRKgoiIiFoZpI6InsqcTf0rLYiIiKiVFsRuIpevRsRUSwsiIiJqpQUREbuFjFVMvbQgIiKiVhJERETUSoKIiIhaSRAREVErCSIiImrlKqYey/0NEdGvkiAiYlrJ5bCTp9UEIWkx8D6KNak/ant1x/F9gY8Djwd+CbzI9tby2KnASRRrUp9i+8I2Y51saRlExO6utQQhaQZwFnAMMAyslzRo++pKsZOAm23/saRlwDuBF0laCCwDHgU8GPiqpIfb3tlWvPnVEbFnyd/82NpsQSwCNtveAiBpLbAUqCaIpcBby9efBT4oSeX+tbbvBH4kaXN5vm+3GG8jaRlE7FnG8zc/3ZJJmwliFrCtsj0MHDlaGds7JN0CPKDcf2lH3VmdbyBpObC83Py1pE2TE/qkOhT4Ra+DGEVim5jENjHTPja9cxIi2VXb/24PGe1AmwlCNfvcsEyTutg+Bzhn/KFNHUlDtgd6HUedxDYxiW1iEtvE9DK2Nu+DGAbmVLZnA9tHKyNpL+Bg4KaGdSMiokVtJoj1wAJJ8yXtQzHoPNhRZhA4sXx9PPA12y73L5O0r6T5wALgOy3GGhERHVrrYirHFFYAF1Jc5rrG9kZJq4Ah24PAucD55SD0TRRJhLLcBRQD2juAk9u8gqll/dwFltgmJrFNTGKbmJ7FpuIHe0RExD1lLqaIiKiVBBEREbWSIKaApLdLulLSFZK+IunBvY5phKQzJf2gjO/zku7X65hGSHqhpI2S7pbUF5cgSlosaZOkzZJW9jqeEZLWSLpB0vd7HUuVpDmSvi7pmvK/5Wt7HdMISftJ+o6k75Wxva3XMXWSNEPS5ZK+2Iv3T4KYGmfaPsz24cAXgdN6HVDFRcCjbR8G/BA4tcfxVH0f+HPgm70OBO4xfcyxwELghHJamH7wMWBxr4OosQP4O9uPBJ4AnNxH/2Z3As+w/VjgcGCxpCf0OKZOrwWu6dWbJ0FMAdu3Vjb3p+amv16x/RXbO8rNSynuOekLtq+x3U93x/9++hjbdwEj08f0nO1vUlwJ2Fds/9T2d8vXt1F82e0yK0IvuPDrcnPv8tE3f5uSZgPHAR/tVQxJEFNE0jskbQNeQn+1IKpeCXy510H0sbrpY/riy253IGkecARwWW8j+YOyC+cK4AbgItt9ExvwXuCNwN29CiAJYpJI+qqk79c8lgLYfrPtOcAngRX9FFtZ5s0U3QGf7LfY+kijKWBiV5IOAD4HvK6jRd1TtneWXb+zgUWSHt3rmAAkPQe4wfaGXsaRBYMmie2jGxb9FPAl4PQWw7mHsWKTdCLwHOCZnuIbY8bx79YPMgXMBEjamyI5fNL2v/c6njq2fyXpYopxnH4Y6D8KWCLpz4D9gIMkfcL2S6cyiLQgpoCkBZXNJcAPehVLp3JRpzcBS2z/ptfx9Lkm08dERTl9/7nANbbf0+t4qiTNHLlqT9J9gaPpk79N26fanm17HsX/Z1+b6uQASRBTZXXZbXIl8CyKKxP6xQeBA4GLystwP9LrgEZIer6kYeCJwJck9XRVwXIwf2T6mGuAC2xv7GVMIyR9mmK9lEdIGpZ0Uq9jKh0FvAx4Rvn/1xXlr+J+8CDg6+Xf5XqKMYieXE7arzLVRkRE1EoLIiIiaiVBRERErSSIiIiolQQRERG1kiAiIqJWEkRERNRKgoiIiFpJEBEtKtePGLlB7DJJ+ZuL3UZulItokaRrgSfb/lmvY4kYr/yaiWjXOuAqSe/tdSAR45XZXCNaIulPKaYIf1BlUaaI3UZaEBHteSHwQ9s7VDio1wFFjEfGICJaImkRxVTXBu4AXtPrBWAixiMJIiIiaqWLKSIiaiVBRERErSSIiIiolQQRERG1kiAiIqJWEkRERNRKgoiIiFr/HwmffsKms+JsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The distribution of $\\epsilon$\n",
    "count, bins, ignored = plt.hist(epsilon, 30, density=True)\n",
    "plt.title(\"$\\epsilon$ $\\sim$ Normal$(0,2^2)$\")\n",
    "plt.xlabel('$\\epsilon$')\n",
    "plt.ylabel('probability density')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEaCAYAAADg2nttAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAab0lEQVR4nO3de7wkZX3n8c+Xi2DkJoKucnFQ0Q0KKo6igTWKq1wk4KoYlFVMSIgRX4Z1swrRFSW6wTtqDIqCgDdwvSyIZA1eUIkKzAioSJARBiGwgstdBAV++aNqSHPmnK6emVN9+pzzeb9e/equp56n6vdMn+lf1/NUV6WqkCRpmPXmOgBJ0uQzWUiSOpksJEmdTBaSpE4mC0lSJ5OFJKmTyUKS1MlkIUnqZLLQgpXk75IcMaXsmiS7zkEsc7LfUU16fMMkuSDJE+c6joXOZKGhklyW5Nr59p8xydbAq4CPDZRtATwK+Jdp6m+Z5MtJfp3k6iSvmMVYZtzvJJj0+EbwXuCYuQ5ioTNZqMuTgJ8BL5nrQNbQq4Gzq+o3A2U7A7+oqjunqf8R4LfAI4CDgeO7EmSStyV52wixDNtvr0aMcc7igzX6d5zJmcBzkzxylkLSNEwWGqqq7gXOA57cx/aTPCrJF5PcmOSqJK9vyx+b5KZVQyNtvV8leU67vDLJUUl+muTmJJ9MsvHApvcBvj1ldzsDlyd5f9vmiiR7JHkITTL8n1V1R1WdR/MB9MpZ6ua0+x3WIMm7k3x5YPk9Sb6RZMNZimmd4ht3jEk2T/K59m/g9iQXJ1kPoKruApYDL5jt/erfmSw0VJIHAwcBu3TUOyvJLTM8zpqhzXrAV4BLgG2A5wFHJNmrqn4OvAn4TJLfAz4JnFxV5w5s4mBgL+CxwOOBtwys2xm4fMoudwF2A74DPBz4NPDxtu29VfWzgbqXALM19DbTfod5F8235ackeQ2wN/DiqvrdLMW0rvGNO8Z3AfcCjwY2B/5rVd03sP4yevpCo4bJQl3eCfwr8Ngkm8xUqar2q6otZnjsN0OzpwNbV9UxVfXbqrqS5kPqoHabHweuAM4HHgm8eUr7v6+qa6rqpjbOlw+s2wK4fUr9nYEPVNX/aT/QPgE8AdgEuHVK3VuBTWfq7xqadr/tt+ULktyR5EmDDarq/wPHAacCRwH7VtWtw9r0EN8GSZ6V5PtJvt1+s7//qGFIjDO2WQe/o/lbuLOq7quqn0xZfzvNe66emCw0oyTPAl5GM0RzK838xWx6NPCowaMQ4G9o5g1W+Xi73w9X1d1T2l8z8PpqmknaVW5m9Q/7JwFfGFjeiqZfdwCbTam7GasnmwccQQFHAkd2HUEN2e+dwAunrBt0Ec0H+VFVtaqvXW3WJsZp46uqe2j+Xfesqj8ErgQOGCHGrjZrE+O/AH8N/DrJX0yzflPglun+PTQ7NpjrADSZ2vH/k4DXVNVNSS6hOcz/wQz1/xH4TzNs7rtVtc805dcAV1XVjjNscxOab64nAm9L8sX2KGKV7QZebw9cN7D8I5rhpQvbbT2aJgHcOFDnxcBZNBP4GyTZsaquaNc9Gbh0akyDR0mrJmWr6m3Txd+13/Zb/I1Jpmu3M3A8cArwp8Bn233N2GZtYuz4d6GqBv9N7wHuG2g7U4wztlnLGPcHXgfsOvD+TPX7NMNn6olHFprJMcD3q2rVt7yLGTJvUVX7VNUmMzymSxQAFwC3JXlTkgcnWT/Jk5I8vV3/QWB5Vf0Z8FXgo1PaH55k2yRb0hyRnD6w7mzgDweWd6b54HpFkvWS7Av8BXBMVf0a+BJwTJKHJNmd5tvwp4b+C41mxv3O1CDJNjRzOa8BXgvsnHZivwcjxZdkB5qTBs4aNcapbdbBk2iGQn/Zbnf7JA8d2M9GwNOAc9ZxPxrCZKHVJHkGcCDw3waKL2aWJxDbM63+CHgKcBXwK5rx8s2THEAzYfqatvobgF2THDywic8C/0Qz1HEl8I6BdacC+7YT9NB8KJ4K7E4zRPU24ICBb6qvBR4M3AB8DvjLqlrtyGItdO33AZJsRpPo3l9VZ7ans76HZk6mD53xtTGdAryyqn47SoxT26xjjJ+kOa35miS30iT2wVGR/YFzpxzRaJbF26pqPkqyEvizqvr6kDr/C7ihqo4bW2BrIcnJwHunmbSd1TZrI8kGwBnA+6rqm321WRdJzgcO7fvfYrEzWWheGiVZzAdJzqY5sroa+FhVndxHm3WI75XAB4BVH8THV9XpQ5qsVRtNPpOF5qWFkiyk+cJkIUnq5AS3JKnTgvydxVZbbVVLliyZ6zAkaV5Zvnz5r6pq6+nWLchksWTJEpYtWzbXYUjSvJLk6pnWOQwlSepkspAkdTJZSJI6mSwkSZ1MFpKkTiYLSVInk4UkqZPJQpLUyWQhSeq0IH/Bva6WHPnVkeqtPPaFPUciSZPBIwtJUieThSSpk8lCktTJZCFJ6mSykCR1MllIkjqZLCRJnUwWkqROJgtJUieThSSpk8lCktSp92SRZP0kFyU5q13eIcn5Sa5IcnqSB7XlG7XLK9r1Swa2cVRbfnmSvfqOWZL0QOM4svgr4LKB5XcBH6iqHYGbgUPb8kOBm6vqccAH2nok2Qk4CHgisDfwD0nWH0PckqRWr8kiybbAC4FPtMsB9gS+0FY5BXhR+/qAdpl2/fPa+gcAp1XV3VV1FbACeEafcUuSHqjvI4vjgDcC97XLDwNuqap72uVrgW3a19sA1wC0629t699fPk0bSdIY9JYskuwH3FBVyweLp6laHeuGtRnc32FJliVZduONN65xvJKkmfV5ZLE7sH+SlcBpNMNPxwFbJFl106Vtgeva19cC2wG06zcHbhosn6bN/arqhKpaWlVLt95669nvjSQtYr0li6o6qqq2raolNBPU36yqg4FvAS9tqx0CnNG+PrNdpl3/zaqqtvyg9mypHYAdgQv6iluStLq5uK3qm4DTkrwDuAg4sS0/EfhUkhU0RxQHAVTVpUk+D/wUuAc4vKruHX/YkrR4jSVZVNW5wLnt6yuZ5mymqroLOHCG9u8E3tlfhJKkYfwFtySpk8lCktTJZCFJ6mSykCR1MllIkjqZLCRJnUwWkqROJgtJUieThSSpk8lCktTJZCFJ6mSykCR1MllIkjqZLCRJnUwWkqROJgtJUieThSSpk8lCktTJZCFJ6mSykCR1MllIkjqZLCRJnUwWkqROJgtJUieThSSpk8lCktTJZCFJ6mSykCR1MllIkjqZLCRJnUwWkqROJgtJUieThSSpU2eySLIsyeFJHjqOgCRJk2eUI4uDgEcBFyY5LcleSdJzXJKkCdKZLKpqRVW9GXg88FngJOAXSd6eZMu+A5Qkzb2R5iyS7AK8D3gP8EXgpcBtwDf7C02SNCk26KqQZDlwC3AicGRV3d2uOj/J7n0GJ0maDJ3JAjiwqq4cLEiyQ1VdVVUv7ikuSdIEGWUY6gsjlj1Ako2TXJDkkiSXJnl7W75DkvOTXJHk9CQPass3apdXtOuXDGzrqLb88iR7jdY1SdJsmfHIIsl/BJ4IbJ5k8AhiM2DjEbZ9N7BnVd2RZEPgvCT/CLwB+EBVnZbko8ChwPHt881V9bgkBwHvAv44yU40Z2Q9keasrK8neXxV3bvGvZUkrZVhRxZPAPYDtgD+aOCxK/DnXRuuxh3t4obto4A9+fcjk1OAF7WvD2iXadc/rz1F9wDgtKq6u6quAlYAzxipd5KkWTHjkUVVnQGckeRZVfX9tdl4kvWB5cDjgI8APwduqap72irXAtu0r7cBrmn3fU+SW4GHteU/GNjsYJvBfR0GHAaw/fbbr024kqQZDBuGemNVvRt4RZKXT11fVa/v2ng7VPSUJFsAXwZ+f7pqq3Y5w7qZyqfu6wTgBIClS5eutl6StPaGnQ11Wfu8bF13UlW3JDkXeCawRZIN2qOLbYHr2mrXAtsB1ybZANgcuGmgfJXBNpKkMRg2DPWV9nnVPAJJ1gM2qarbujacZGvgd22ieDDwn2kmrb9F86O+04BDgDPaJme2y99v13+zqirJmcBnk7yfZoJ7R+CCNe2oJGntjXIhwc8m2SzJQ4CfApcn+R8jbPuRwLeS/Ai4EDinqs4C3gS8IckKmjmJE9v6JwIPa8vfABwJUFWXAp9v9/1/gcM9E0qSxmuUH+XtVFW3JTkYOJvmw345zaU/ZlRVPwKeOk35lUxzNlNV3QUcOMO23gm8c4RYJUk9GOVHeRu2v5N4EXBGVf2OaSaYJUkL1yjJ4mPASuAhwHeSPJrmIoKSpEWicxiqqj4EfGig6Ookz+0vJEnSpBnlqrMbAS8Blkypf0xPMUmSJswoE9xnALfSTGrf3VFXkrQAjZIstq2qvXuPRJI0sUaZ4P5ekp17j0SSNLFGObLYA3h1kqtohqFCc1HZXXqNTJI0MUZJFvv0HoUkaaJ1DkNV1dU0F/Lbs3195yjtJEkLxyjXhjqa5hIfR7VFGwKf7jMoSdJkGeUI4b8A+wO/Bqiq64BN+wxKkjRZRkkWv62qor0eVHv1WUnSIjJKsvh8ko/R3LToz4GvAx/vNyxJ0iQZ5dpQ703yfJqLBz4BeGtVndN7ZJKkiTHKqbO0ycEEIUmL1IzJIsntDLlvRVVt1ktEkqSJM+we3JsCJDkG+H/Ap2h+vX0wng0lSYvKKBPce1XVP1TV7VV1W1UdT3PJcknSIjFKsrg3ycFJ1k+yXnsv7nv7DkySNDlGSRavAF4G/LJ9HNiWSZIWiVFOnV0JHNB/KJKkSeUFASVJnUwWkqROo1x1dv1xBCJJmlyjHFmsSPKeJDv1Ho0kaSKNkix2AX4GfCLJD5IclsRfb0vSIjLKnfJur6qPV9UfAG8EjgauT3JKksf1HqEkac6NNGeRZP8kXwY+CLwPeAzwFeDsnuOTJE2AUa46ewXwLeA9VfW9gfIvJHl2P2FJkibJKMniVVV13mBBkt2r6p+r6vU9xSVJmiCjTHB/aJqyD892IJKkyTXsfhbPAv4A2DrJGwZWbQb42wtJWkSGDUM9CNikrTN4/4rbgJf2GZQkabIMu/nRt4FvJzm5qq4eY0ySpAkzbBjquKo6Avj7JKvdXrWq9u81MknSxBg2DPWp9vm94whEkjS5hg1DLW+fvz2+cCRJk2jYMNSPgdWGn1apql16iUiSNHGGDUPtty4bTrIdcCrwH4D7gBOq6oNJtgROB5YAK4GXVdXNSUJzOZF9gTuBV1fVD9ttHQK8pd30O6rqlHWJTZK0ZoYNQ63rGVD3AP+9qn6YZFNgeZJzgFcD36iqY5McCRwJvAnYB9ixfewGHA/s1iaXo4GlNEc6y5OcWVU3r2N8kqQRzfgL7iTntc+3J7lt6nPXhqvq+lVHBlV1O3AZsA3N/bxXHRmcAryofX0AcGo1fgBskeSRwF7AOVV1U5sgzgH2XqveSpLWyrAjiz3a501nqjOqJEuApwLnA4+oquvbbV+f5OFttW2AawaaXduWzVQ+dR+HAYcBbL/99usasiRpwCgXEiTJrsAeNMNA51XVRaPuIMkmwBeBI6rqtmZqYvqq05TVkPIHFlSdAJwAsHTp0hkn5iVJa26U+1m8lWa46GHAVsDJSd4yvNX9bTekSRSfqaovtcW/bIeXaJ9vaMuvBbYbaL4tcN2QcknSmIxy1dmXA0+vqqOr6mjgmcDBXY3as5tOBC6rqvcPrDoTOKR9fQhwxkD5q9J4JnBrO1z1NeAFSR6a5KHAC9oySdKYjDIMtRLYGLirXd4I+PkI7XYHXgn8OMnFbdnfAMcCn09yKPAL4MB23dk0p82uoDl19k8AquqmJH8LXNjWO6aqbhph/5KkWTLsR3kfppkbuBu4tD3ttYDnA+fN1G6V9oZJM01QPG+a+gUcPsO2TgJO6tqnJKkfw44slrXPy4EvD5Sf21s0kqSJNOzUWX8lLUkCRpizSLIj8HfATjRzFwBU1WN6jEuSNEFGORvqkzSX3rgHeC7N9Z4+NbSFJGlBGSVZPLiqvgGkqq6uqrcBe/YbliRpkoxy6uxdSdYDrkjyOuBfgYd3tJEkLSCjHFkcAfwe8HrgaTS/nThkaAtJ0oLSeWRRVRcCtEcXr2+vICtJWkRGuTbU0vaueT+i+TX2JUme1n9okqRJMcqcxUnAa6vquwBJ9qA5Q8rbqkrSIjHKnMXtqxIF3H8ZD4eiJGkRGXZtqF3blxck+RjwOZprQ/0xXvJDkhaVYcNQ75uyfPTAa28uJEmLyLBrQz13nIFIkibXKGdDbZ7k/UmWtY/3Jdl8HMFJkibDKBPcJ9FMaL+sfdxGczaUJGmRGOXU2cdW1UsGlt8+cOc7SdIiMMqRxW/a31YAkGR34Df9hSRJmjSjHFm8Bjh1YJ7iZrw2lCQtKkOTRXs9qCdU1ZOTbAZQVbeNJTJJ0sQYOgxVVfcBr2tf32aikKTFaZQ5i3OS/HWS7ZJsuerRe2SSpIkxypzFn7bPhw+UFeA9uCVpkRjlfhY7jCMQSdLk6kwWSTYGXgvsQXNE8V3go1V1V8+xSZImxCjDUKfS/IL7w+3yy4FPAQf2FZQkabKMkiyeUFVPHlj+VpJL+gpIkjR5Rjkb6qIkz1y1kGQ34J/7C0mSNGlGObLYDXhVkl+0y9sDl7X35a6q8vaqkrTAjZIs9u49CknSRBvl1NmrxxGIJGlyjTJnIUla5EwWkqROJgtJUieThSSpk8lCktTJZCFJ6mSykCR1MllIkjr1liySnJTkhiQ/GSjbMsk5Sa5onx/alifJh5KsSPKjJLsOtDmkrX9FkkP6ileSNLM+jyxOZvVLhRwJfKOqdgS+0S4D7APs2D4OA46HJrkAR9Ncn+oZwNGrEowkaXx6SxZV9R3gpinFBwCntK9PAV40UH5qNX4AbJHkkcBewDlVdVNV3Qycg9eqkqSxG/ecxSOq6nqA9vnhbfk2wDUD9a5ty2YqX02Sw5IsS7LsxhtvnPXAJWkxm5QJ7kxTVkPKVy+sOqGqllbV0q233npWg5OkxW7cyeKX7fAS7fMNbfm1wHYD9bYFrhtSLkkao3EnizOBVWc0HQKcMVD+qvasqGcCt7bDVF8DXpDkoe3E9gvaMknSGI1y86O1kuRzwHOArZJcS3NW07HA55McCvwCOLCtfjawL7ACuBP4E4CquinJ3wIXtvWOqaqpk+aSpJ71liyq6uUzrHreNHULOHyG7ZwEnDSLoUmS1tCkTHBLkiaYyUKS1MlkIUnqZLKQJHUyWUiSOpksJEmdTBaSpE4mC0lSp95+lLcYLDnyqyPVW3nsC3uORJL65ZGFJKmTyUKS1MlkIUnqZLKQJHUyWUiSOpksJEmdTBaSpE4mC0lSJ5OFJKmTyUKS1MlkIUnqZLKQJHUyWUiSOpksJEmdTBaSpE4mC0lSJ5OFJKmTyUKS1MlkIUnqZLKQJHUyWUiSOpksJEmdTBaSpE4bzHUAi8GSI786Ur2Vx76w50gkae14ZCFJ6mSykCR1MllIkjqZLCRJnUwWkqROng01QTxrStKkmjdHFkn2TnJ5khVJjpzreCRpMZkXySLJ+sBHgH2AnYCXJ9lpbqOSpMVjvgxDPQNYUVVXAiQ5DTgA+OmcRjVHRh2umksOlUkLy3xJFtsA1wwsXwvsNlghyWHAYe3iHUkuX8t9bQX8ai3bTrqx9S3vGsdeHmChvm8LtV9g3ybRo2daMV+SRaYpqwcsVJ0AnLDOO0qWVdXSdd3OJLJv889C7RfYt/lmXsxZ0BxJbDewvC1w3RzFIkmLznxJFhcCOybZIcmDgIOAM+c4JklaNObFMFRV3ZPkdcDXgPWBk6rq0p52t85DWRPMvs0/C7VfYN/mlVRVdy1J0qI2X4ahJElzyGQhSeq04JNF12VCkmyU5PR2/flJlgysO6otvzzJXqNuc1x66ttJSW5I8pPx9GJ1s92vJNsl+VaSy5JcmuSvxteb1WKf7b5tnOSCJJe0fXv7+HqzWuyz/vfYrls/yUVJzuq/F6vr6f/ZyiQ/TnJxkmXj6ck6qqoF+6CZDP858BjgQcAlwE5T6rwW+Gj7+iDg9Pb1Tm39jYAd2u2sP8o252vf2nXPBnYFfrKA3rNHAru2dTYFfrZQ3jOa3yBt0tbZEDgfeOZC6NtAuzcAnwXOWij9AlYCW427P+vyWOhHFvdfJqSqfgusukzIoAOAU9rXXwCelyRt+WlVdXdVXQWsaLc3yjbHoY++UVXfAW4aRwdmMOv9qqrrq+qHAFV1O3AZzVUBxq2PvlVV3dHW37B9zMVZK738PSbZFngh8Ikx9GE6vfRrPlroyWK6y4RM/ZC4v05V3QPcCjxsSNtRtjkOffRtEvTar3aI4Kk038DHrZe+tcM0FwM3AOdU1YLpG3Ac8EbgvtkPeSR99auAf0qyPM2liibeQk8WnZcJGVJnTcvHrY++TYLe+pVkE+CLwBFVddtaR7j2eulbVd1bVU+hubLBM5I8aZ2iXDuz3rck+wE3VNXydQ1uHfT197h7Ve1KcyXtw5M8e+1DHI+FnixGuUzI/XWSbABsTjMMM1PbSbn0SB99mwS99CvJhjSJ4jNV9aVeIu/W63tWVbcA5wJ7z2bQI+qjb7sD+ydZSTP8s2eST/cR/BC9vGdVter5BuDLzIfhqbmeNOnzQfML9StpJpdWTU49cUqdw3ng5NTn29dP5IGTU1fSTHZ1bnO+9m2g3RLmboK7j/cswKnAcQvw73FrYIu2zoOB7wL7LYS+TWn7HOZmgruP9+whwKZtnYcA3wP2nsu/zZH+LeY6gDG82fvSnP3yc+DNbdkxwP7t642B/00z+XQB8JiBtm9u210O7DNsmwuob58Drgd+R/PN6ND53i9gD5rD/x8BF7ePfRfCewbsAlzU9u0nwFsX0t/jwPrnMAfJoqf37DE0SeQS4NK5/AxZk4eX+5AkdVrocxaSpFlgspAkdTJZSJI6mSwkSZ1MFpKkTiYLSVInk4UkqZPJQhqTJDsnuTrJX851LNKaMllIY1JVP6a5HMSr5joWaU2ZLKTxuoHmmkHSvGKykMbrWGCjJI+e60CkNWGykMYkyd40Vxn9Kh5daJ4xWUhjkGRj4N0092v+MTAXNyiS1prJQhqPtwCnVtVKTBaah0wWUs+SPAF4Ps39pMFkoXnI+1lIkjp5ZCFJ6mSykCR1MllIkjqZLCRJnUwWkqROJgtJUieThSSp078BuqP+ZtPJEEsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The distribution of $\\lambda$\n",
    "count, bins, ignored = plt.hist(lambda_i, 30, density=True)\n",
    "plt.title(\"$\\lambda$ = exp($b0 + b_1x_1 + b_2x_2 + \\epsilon$)\")\n",
    "plt.xlabel('$\\lambda$')\n",
    "plt.ylabel('probability density')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create and save dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "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>y</th>\n",
       "      <th>Distance</th>\n",
       "      <th>Precipitation</th>\n",
       "      <th>Traffic</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1018</td>\n",
       "      <td>0</td>\n",
       "      <td>0.299886</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>973</td>\n",
       "      <td>0</td>\n",
       "      <td>0.565617</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>1021</td>\n",
       "      <td>0</td>\n",
       "      <td>0.414564</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>998</td>\n",
       "      <td>0</td>\n",
       "      <td>0.559767</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>985</td>\n",
       "      <td>0</td>\n",
       "      <td>0.777217</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0</td>\n",
       "      <td>994</td>\n",
       "      <td>0</td>\n",
       "      <td>0.716722</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0</td>\n",
       "      <td>991</td>\n",
       "      <td>0</td>\n",
       "      <td>0.782984</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1</td>\n",
       "      <td>967</td>\n",
       "      <td>0</td>\n",
       "      <td>0.587694</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0</td>\n",
       "      <td>1003</td>\n",
       "      <td>0</td>\n",
       "      <td>0.595906</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0</td>\n",
       "      <td>1012</td>\n",
       "      <td>0</td>\n",
       "      <td>0.742293</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   y  Distance  Precipitation   Traffic\n",
       "0  0      1018              0  0.299886\n",
       "1  0       973              0  0.565617\n",
       "2  0      1021              0  0.414564\n",
       "3  0       998              0  0.559767\n",
       "4  0       985              0  0.777217\n",
       "5  0       994              0  0.716722\n",
       "6  0       991              0  0.782984\n",
       "7  1       967              0  0.587694\n",
       "8  0      1003              0  0.595906\n",
       "9  0      1012              0  0.742293"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame({\n",
    "    'y': y, \n",
    "    'Distance': d, \n",
    "    'Precipitation': x1, \n",
    "    'Traffic': x2\n",
    "})\n",
    "df.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv(\"data/simulated_data.csv\", sep=',', encoding='utf-8')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate precipitation and traffic variables for links data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "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>Node A</th>\n",
       "      <th>Node Z</th>\n",
       "      <th>Distance</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>node 1</td>\n",
       "      <td>node 7</td>\n",
       "      <td>800</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>node 1</td>\n",
       "      <td>node 11</td>\n",
       "      <td>800</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>node 1</td>\n",
       "      <td>node 12</td>\n",
       "      <td>2400</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>node 2</td>\n",
       "      <td>node 6</td>\n",
       "      <td>1200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>node 2</td>\n",
       "      <td>node 10</td>\n",
       "      <td>900</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>node 3</td>\n",
       "      <td>node 6</td>\n",
       "      <td>1100</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>node 3</td>\n",
       "      <td>node 8</td>\n",
       "      <td>800</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>node 3</td>\n",
       "      <td>node 12</td>\n",
       "      <td>600</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>node 4</td>\n",
       "      <td>node 8</td>\n",
       "      <td>700</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>node 4</td>\n",
       "      <td>node 10</td>\n",
       "      <td>700</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>node 4</td>\n",
       "      <td>node 14</td>\n",
       "      <td>2800</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>node 5</td>\n",
       "      <td>node 6</td>\n",
       "      <td>2000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>node 5</td>\n",
       "      <td>node 7</td>\n",
       "      <td>500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>node 5</td>\n",
       "      <td>node 11</td>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>node 6</td>\n",
       "      <td>node 13</td>\n",
       "      <td>2000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>node 7</td>\n",
       "      <td>node 10</td>\n",
       "      <td>500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>node 9</td>\n",
       "      <td>node 12</td>\n",
       "      <td>1000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>node 9</td>\n",
       "      <td>node 13</td>\n",
       "      <td>600</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>node 9</td>\n",
       "      <td>node 14</td>\n",
       "      <td>1100</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>node 10</td>\n",
       "      <td>node 11</td>\n",
       "      <td>500</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     Node A   Node Z  Distance\n",
       "0    node 1   node 7       800\n",
       "1    node 1  node 11       800\n",
       "2    node 1  node 12      2400\n",
       "3    node 2   node 6      1200\n",
       "4    node 2  node 10       900\n",
       "5    node 3   node 6      1100\n",
       "6    node 3   node 8       800\n",
       "7    node 3  node 12       600\n",
       "8    node 4   node 8       700\n",
       "9    node 4  node 10       700\n",
       "10   node 4  node 14      2800\n",
       "11   node 5   node 6      2000\n",
       "12   node 5   node 7       500\n",
       "13   node 5  node 11       300\n",
       "14   node 6  node 13      2000\n",
       "15   node 7  node 10       500\n",
       "16   node 9  node 12      1000\n",
       "17   node 9  node 13       600\n",
       "18   node 9  node 14      1100\n",
       "19  node 10  node 11       500"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "links = pd.read_csv('data/links.csv')\n",
    "links.head(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "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>Node A</th>\n",
       "      <th>Node Z</th>\n",
       "      <th>Distance</th>\n",
       "      <th>Precipitation</th>\n",
       "      <th>Traffic</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>node 1</td>\n",
       "      <td>node 7</td>\n",
       "      <td>800</td>\n",
       "      <td>0</td>\n",
       "      <td>0.254345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>node 1</td>\n",
       "      <td>node 11</td>\n",
       "      <td>800</td>\n",
       "      <td>0</td>\n",
       "      <td>0.243435</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>node 1</td>\n",
       "      <td>node 12</td>\n",
       "      <td>2400</td>\n",
       "      <td>0</td>\n",
       "      <td>0.254188</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>node 2</td>\n",
       "      <td>node 6</td>\n",
       "      <td>1200</td>\n",
       "      <td>0</td>\n",
       "      <td>0.424037</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>node 2</td>\n",
       "      <td>node 10</td>\n",
       "      <td>900</td>\n",
       "      <td>0</td>\n",
       "      <td>0.573477</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>node 3</td>\n",
       "      <td>node 6</td>\n",
       "      <td>1100</td>\n",
       "      <td>0</td>\n",
       "      <td>0.834188</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>node 3</td>\n",
       "      <td>node 8</td>\n",
       "      <td>800</td>\n",
       "      <td>0</td>\n",
       "      <td>0.190679</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>node 3</td>\n",
       "      <td>node 12</td>\n",
       "      <td>600</td>\n",
       "      <td>1</td>\n",
       "      <td>0.689211</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>node 4</td>\n",
       "      <td>node 8</td>\n",
       "      <td>700</td>\n",
       "      <td>1</td>\n",
       "      <td>0.507221</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>node 4</td>\n",
       "      <td>node 10</td>\n",
       "      <td>700</td>\n",
       "      <td>0</td>\n",
       "      <td>0.129693</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>node 4</td>\n",
       "      <td>node 14</td>\n",
       "      <td>2800</td>\n",
       "      <td>0</td>\n",
       "      <td>0.192684</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>node 5</td>\n",
       "      <td>node 6</td>\n",
       "      <td>2000</td>\n",
       "      <td>0</td>\n",
       "      <td>0.334005</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>node 5</td>\n",
       "      <td>node 7</td>\n",
       "      <td>500</td>\n",
       "      <td>0</td>\n",
       "      <td>0.568231</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>node 5</td>\n",
       "      <td>node 11</td>\n",
       "      <td>300</td>\n",
       "      <td>1</td>\n",
       "      <td>0.499119</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>node 6</td>\n",
       "      <td>node 13</td>\n",
       "      <td>2000</td>\n",
       "      <td>0</td>\n",
       "      <td>0.441698</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>node 7</td>\n",
       "      <td>node 10</td>\n",
       "      <td>500</td>\n",
       "      <td>0</td>\n",
       "      <td>0.167268</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>node 9</td>\n",
       "      <td>node 12</td>\n",
       "      <td>1000</td>\n",
       "      <td>0</td>\n",
       "      <td>0.776060</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>node 9</td>\n",
       "      <td>node 13</td>\n",
       "      <td>600</td>\n",
       "      <td>0</td>\n",
       "      <td>0.514337</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>node 9</td>\n",
       "      <td>node 14</td>\n",
       "      <td>1100</td>\n",
       "      <td>0</td>\n",
       "      <td>0.122787</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>node 10</td>\n",
       "      <td>node 11</td>\n",
       "      <td>500</td>\n",
       "      <td>1</td>\n",
       "      <td>0.741051</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     Node A   Node Z  Distance  Precipitation   Traffic\n",
       "0    node 1   node 7       800              0  0.254345\n",
       "1    node 1  node 11       800              0  0.243435\n",
       "2    node 1  node 12      2400              0  0.254188\n",
       "3    node 2   node 6      1200              0  0.424037\n",
       "4    node 2  node 10       900              0  0.573477\n",
       "5    node 3   node 6      1100              0  0.834188\n",
       "6    node 3   node 8       800              0  0.190679\n",
       "7    node 3  node 12       600              1  0.689211\n",
       "8    node 4   node 8       700              1  0.507221\n",
       "9    node 4  node 10       700              0  0.129693\n",
       "10   node 4  node 14      2800              0  0.192684\n",
       "11   node 5   node 6      2000              0  0.334005\n",
       "12   node 5   node 7       500              0  0.568231\n",
       "13   node 5  node 11       300              1  0.499119\n",
       "14   node 6  node 13      2000              0  0.441698\n",
       "15   node 7  node 10       500              0  0.167268\n",
       "16   node 9  node 12      1000              0  0.776060\n",
       "17   node 9  node 13       600              0  0.514337\n",
       "18   node 9  node 14      1100              0  0.122787\n",
       "19  node 10  node 11       500              1  0.741051"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(0)\n",
    "links['Precipitation'] = simulate_precipitation(links.shape[0])\n",
    "links['Traffic'] = simulate_traffic(links.shape[0])\n",
    "links.head(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "links.to_csv(\"data/links_traffic_precipitation.csv\", sep=',', encoding='utf-8')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py37",
   "language": "python",
   "name": "py37"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}