MP_expenses.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MP expenses\n",
    "\n",
    "This dataset is from [http://www.theipsa.org.uk/mp-costs/your-mp/](http://www.theipsa.org.uk/mp-costs/your-mp/). Code for retrieving the dataset is given at the bottom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = 'https://www.cl.cam.ac.uk/teaching/2021/DataSci/data/mpexpenses_2018_office.csv'\n",
    "mpexpenses = pandas.read_csv(url)['Amount Claimed'].values\n",
    "# We'll be trying to model log(amount claimed). There are some claims <= 0,\n",
    "# which certainly merit understanding, but that's not the point of this exercise.\n",
    "mpexpenses = mpexpenses[mpexpenses>0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Initial exploration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQ/UlEQVR4nO3dbaxlVX3H8e+vFNRUW7DcknFm7BA7tsEmjuYGMPQFxYoDNkWb1mBSmRiS8QUkmJi04Bt8KAkmKq2pJUEhYmulpGqcKCkdEUNMCsxFERlGwi1CmMnIjIJPMaEZ/PfFXVOOw525T+eec+9Z309ycvf5733OXYuH31l37bX3SVUhSerDb4y7AZKk0TH0Jakjhr4kdcTQl6SOGPqS1JHfHHcDTuT000+vLVu2jLsZkrSuPPDAAz+qqqn59q3p0N+yZQszMzPjboYkrStJnjzePqd3JKkjC4Z+kpcmuT/Jd5PsTfKhVj8zyX1JZpP8e5JTWv0l7fls279l4L2uafVHk7x11XolSZrXYkb6zwEXVNXrgW3A9iTnAh8FbqiqPwCeBS5vx18OPNvqN7TjSHIWcCnwOmA78M9JThpiXyRJC1gw9GvOL9rTk9ujgAuA/2j1W4G3t+1L2nPa/jcnSavfVlXPVdUPgFng7GF0QpK0OIua009yUpIHgUPAbuB/gJ9U1ZF2yH5gY9veCDwF0Pb/FPjdwfo8rxn8XTuTzCSZOXz48JI7JEk6vkWFflU9X1XbgE3Mjc7/aLUaVFU3VdV0VU1PTc274kiStExLWr1TVT8B7gbeBJya5OiSz03AgbZ9ANgM0Pb/DvDjwfo8r5EkjcBiVu9MJTm1bb8MeAuwj7nw/6t22A7gK217V3tO2/+Nmrt/8y7g0ra650xgK3D/kPohSVqExVyctQG4ta20+Q3g9qr6apJHgNuS/D3wHeDmdvzNwL8kmQWeYW7FDlW1N8ntwCPAEeCKqnp+uN2RJJ1I1vKXqExPT5dX5GrSbLn6a7/2/Inr3zamlmhSJXmgqqbn2+cVuZLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHFvPF6JLGyO/U1TA50pekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdWTD0k2xOcneSR5LsTXJVq38wyYEkD7bHxQOvuSbJbJJHk7x1oL691WaTXL06XZIkHc9i1ukfAd5fVd9O8grggSS7274bqupjgwcnOQu4FHgd8Crg60le23Z/CngLsB/Yk2RXVT0yjI5Ikha2YOhX1UHgYNv+eZJ9wMYTvOQS4Laqeg74QZJZ4Oy2b7aqHgdIcls71tCXpBFZ0hW5SbYAbwDuA84DrkxyGTDD3F8DzzL3gXDvwMv288KHxFPH1M9ZXrOlyXXsFbjSMC36RG6SlwNfBN5XVT8DbgReA2xj7i+Bjw+jQUl2JplJMnP48OFhvKUkqVlU6Cc5mbnA/3xVfQmgqp6uquer6lfAp3lhCucAsHng5Zta7Xj1X1NVN1XVdFVNT01NLbU/kqQTWMzqnQA3A/uq6hMD9Q0Dh70DeLht7wIuTfKSJGcCW4H7gT3A1iRnJjmFuZO9u4bTDUnSYixmTv884N3A95I82GofAN6VZBtQwBPAewGqam+S25k7QXsEuKKqngdIciVwJ3AScEtV7R1aTyRJC1rM6p1vAZln1x0neM11wHXz1O840eskSavLK3IlqSOGviR1xG/OktY5v1lLS+FIX5I64khfGjOvwNUoOdKXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0JakjrtOX1pmF1vV7ha5OxJG+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEe+9I3XI+/P0y5G+JHXE0JekjiwY+kk2J7k7ySNJ9ia5qtVfmWR3ksfaz9NaPUk+mWQ2yUNJ3jjwXjva8Y8l2bF63ZIkzWcxI/0jwPur6izgXOCKJGcBVwN3VdVW4K72HOAiYGt77ARuhLkPCeBa4BzgbODaox8UkqTRWPBEblUdBA627Z8n2QdsBC4Bzm+H3Qp8E/i7Vv9cVRVwb5JTk2xox+6uqmcAkuwGtgNfGGJ/pDVnoS89kUZpSXP6SbYAbwDuA85oHwgAPwTOaNsbgacGXra/1Y5XP/Z37Ewyk2Tm8OHDS2meJGkBiw79JC8Hvgi8r6p+NrivjeprGA2qqpuqarqqpqempobxlpKkZlHr9JOczFzgf76qvtTKTyfZUFUH2/TNoVY/AGweePmmVjvAC9NBR+vfXH7TpbXJ6RytZYtZvRPgZmBfVX1iYNcu4OgKnB3AVwbql7VVPOcCP23TQHcCFyY5rZ3AvbDVJEkjspiR/nnAu4HvJXmw1T4AXA/cnuRy4EngnW3fHcDFwCzwS+A9AFX1TJKPAHvacR8+elJXkjQai1m98y0gx9n95nmOL+CK47zXLcAtS2mgJGl4vCJXkjpi6EtSRwx9SeqIoS9JHTH0JakjfomKtEJejKX1xJG+JHXE0Jekjhj6ktQRQ1+SOuKJXGnCeaJZgxzpS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1JEFQz/JLUkOJXl4oPbBJAeSPNgeFw/suybJbJJHk7x1oL691WaTXD38rkiSFrKYL1H5LPBPwOeOqd9QVR8bLCQ5C7gUeB3wKuDrSV7bdn8KeAuwH9iTZFdVPbKCtksakmO/aOWJ6982ppZotS0Y+lV1T5Iti3y/S4Dbquo54AdJZoGz277ZqnocIMlt7VhDX5JGaCVz+lcmeahN/5zWahuBpwaO2d9qx6u/SJKdSWaSzBw+fHgFzZMkHWu5oX8j8BpgG3AQ+PiwGlRVN1XVdFVNT01NDettJUks84vRq+rpo9tJPg18tT09AGweOHRTq3GCuiRpRJY10k+yYeDpO4CjK3t2AZcmeUmSM4GtwP3AHmBrkjOTnMLcyd5dy2+2JGk5FhzpJ/kCcD5wepL9wLXA+Um2AQU8AbwXoKr2JrmduRO0R4Arqur59j5XAncCJwG3VNXeYXdGknRii1m98655yjef4PjrgOvmqd8B3LGk1kmShsorciWpI8s6kSv1xAuXNEkc6UtSRxzpS0t07MhfL+ZfR2uXI31J6ogjfUkv4kh9cjnSl6SOGPqS1BFDX5I64py+dAxX52iSOdKXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHXHJpqQl8zYN65cjfUnqiCN9dc+LsVbOf4brhyN9SeqII31JI+c5gfFxpC9JHXGkL2lBztlPDkf6ktQRQ1+SOmLoS1JHFgz9JLckOZTk4YHaK5PsTvJY+3laqyfJJ5PMJnkoyRsHXrOjHf9Ykh2r0x1J0oksZqT/WWD7MbWrgbuqaitwV3sOcBGwtT12AjfC3IcEcC1wDnA2cO3RDwpJ0ugsGPpVdQ/wzDHlS4Bb2/atwNsH6p+rOfcCpybZALwV2F1Vz1TVs8BuXvxBIklaZcud0z+jqg627R8CZ7TtjcBTA8ftb7Xj1V8kyc4kM0lmDh8+vMzmSZLms+ITuVVVQA2hLUff76aqmq6q6ampqWG9rSSJ5V+c9XSSDVV1sE3fHGr1A8DmgeM2tdoB4Pxj6t9c5u+WtM54cdfasdyR/i7g6AqcHcBXBuqXtVU85wI/bdNAdwIXJjmtncC9sNUkSSO04Eg/yReYG6WfnmQ/c6twrgduT3I58CTwznb4HcDFwCzwS+A9AFX1TJKPAHvacR+uqmNPDkuSVlnmpuTXpunp6ZqZmRl3MzThnHpYe7zr5sokeaCqpufb5w3X1B1DXj0z9CWtOd5vf/V47x1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjriOn2ta67nlpbGkb4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQR772jiecXoUsvMPQ1UQx46cSc3pGkjhj6ktQRQ1+SOrKi0E/yRJLvJXkwyUyrvTLJ7iSPtZ+ntXqSfDLJbJKHkrxxGB2QJC3eMEb6f1pV26pquj2/GrirqrYCd7XnABcBW9tjJ3DjEH63JGkJVmN65xLg1rZ9K/D2gfrnas69wKlJNqzC75ckHcdKQ7+A/0ryQJKdrXZGVR1s2z8EzmjbG4GnBl67v9UkSSOy0nX6f1JVB5L8HrA7yfcHd1ZVJamlvGH78NgJ8OpXv3qFzZMkDVrRSL+qDrSfh4AvA2cDTx+dtmk/D7XDDwCbB16+qdWOfc+bqmq6qqanpqZW0jxJ0jGWHfpJfivJK45uAxcCDwO7gB3tsB3AV9r2LuCytornXOCnA9NAkqQRWMn0zhnAl5McfZ9/q6r/TLIHuD3J5cCTwDvb8XcAFwOzwC+B96zgd6tT3mZBWpllh35VPQ68fp76j4E3z1Mv4Irl/j5J/Tr2w/6J6982ppasf16RK0kdMfQlqSPeWllrin/GS6vLkb4kdcTQl6SOGPqS1BHn9CWtO577WT5H+pLUEUNfkjpi6EtSR5zT10g5F6u1YL57OPXy36KhrzXNG6xJw2XoS5o4/kV5fIa+pHXPvwgXz9DXWPk/qzRart6RpI4Y+pLUEad3JE08pxFf4EhfkjriSF8rstAIyqVy0triSF+SOuJIX6vKuVRpbTH0tSSGuCZVL1fxGvqStAiT8qHgnL4kdcSRvk7I6RxpsjjSl6SOjHykn2Q78I/AScBnqur6UbehZ66rlxZnUv/KHWnoJzkJ+BTwFmA/sCfJrqp6ZJTtkKSVWq8ndkc90j8bmK2qxwGS3AZcAhj6q2Spo5VJHd1Iq229/BU96tDfCDw18Hw/cM7gAUl2Ajvb018keXREbRuF04EfjbsRQzJJfYHJ6s8k9QUmpD/5KDC6vvz+8XasudU7VXUTcNO427EaksxU1fS42zEMk9QXmKz+TFJfYLL6sxb6MurVOweAzQPPN7WaJGkERh36e4CtSc5McgpwKbBrxG2QpG6NdHqnqo4kuRK4k7klm7dU1d5RtmHMJmnaapL6ApPVn0nqC0xWf8bel1TVuNsgSRoRr8iVpI4Y+pLUEUN/hJL8dZK9SX6VZN0uQUuyPcmjSWaTXD3u9qxEkluSHEry8LjbslJJNie5O8kj7b+zq8bdpuVK8tIk9yf5buvLh8bdpmFIclKS7yT56rjaYOiP1sPAXwL3jLshyzVwK42LgLOAdyU5a7ytWpHPAtvH3YghOQK8v6rOAs4FrljH/26eAy6oqtcD24DtSc4db5OG4ipg3zgbYOiPUFXtq6r1foXx/99Ko6r+Fzh6K411qaruAZ4ZdzuGoaoOVtW32/bPmQuXjeNt1fLUnF+0pye3x7pedZJkE/A24DPjbIehr6Wa71Ya6zJYJlmSLcAbgPvG3JRla1MhDwKHgN1VtW770vwD8LfAr8bZCEN/yJJ8PcnD8zzW7WhY60uSlwNfBN5XVT8bd3uWq6qer6ptzF25f3aSPx5zk5YtyZ8Dh6rqgXG3Zc3de2e9q6o/G3cbVpm30ljDkpzMXOB/vqq+NO72DENV/STJ3cyde1mvJ9zPA/4iycXAS4HfTvKvVfU3o26II30tlbfSWKOSBLgZ2FdVnxh3e1YiyVSSU9v2y5j7Do7vj7VRK1BV11TVpqrawtz/M98YR+CDoT9SSd6RZD/wJuBrSe4cd5uWqqqOAEdvpbEPuH0930ojyReA/wb+MMn+JJePu00rcB7wbuCCJA+2x8XjbtQybQDuTvIQcwON3VU1tmWOk8TbMEhSRxzpS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUkf8DMxKi4VuP5dMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots()\n",
    "ax.hist(np.log10(mpexpenses), bins=75)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A fitted probability model\n",
    "\n",
    "def rx(p=0.966, μ1=1.797, μ2=2.861, σ1=0.679, σ2=0.116):\n",
    "    # let k='a' with probability p, 'b' with probability 1-p\n",
    "    k = np.random.choice(['a','b'], p=[p,1-p])\n",
    "    if k == 'a':\n",
    "        y = np.random.normal(loc=μ1, scale=σ1)\n",
    "    else:\n",
    "        y = np.random.normal(loc=μ2, scale=σ2)\n",
    "    return np.power(10, y)\n",
    "\n",
    "xsample = [rx() for _ in range(len(mpexpenses))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A vectorized version of the same probability model\n",
    "# (Exercise: modify the code to accept μ,σ as arrays)\n",
    "\n",
    "def rx(size=1, p=0.966, μ1=1.797, μ2=2.861, σ1=0.679, σ2=0.116):\n",
    "    k = np.random.choice([1,2], p=[p,1-p], size=size)\n",
    "    loc = np.where(k==1, μ1, μ2)\n",
    "    scale = np.where(k==1, σ1, σ2)\n",
    "    return np.power(10, np.random.normal(loc=loc, scale=scale))\n",
    "\n",
    "xsample = rx(len(mpexpenses))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEMCAYAAAAoB2Y1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlKUlEQVR4nO3de3gV1b3/8fcXBESDXAQC5WIoooKg9jSCVEG8VKilx6rVU7UWREQtnMKvUkXrpYBVsNSneEALVoQq9fKIrVQ4Ug4tWi8QgQoq2IoWBIwkFEUIKALf3x8ziTubnWQH9i2Zz+t59uOeNWuvvSbB76x8Z80ac3dERCQaGmS7AyIikjkK+iIiEaKgLyISIQr6IiIRoqAvIhIhCvoiIhFSY9A3syPNrMjMVpvZ22Y2PizvYmbLzWy9mT1lZo3D8ibh9vpwf0FMW7eG5f8ws4FpOyoREUnIapqnb2YGHO3uu8ysEfAyMBr4CfCsuz9pZr8BVrv7Q2b2I+AUd7/BzL4PXOzu/2VmPYAngN7AV4D/A05w9/1VfXfr1q29oKAgBYcpIhIdK1eu3ObubRLtO6KmD3twVtgVbjYKXw6cC1wZls8Bfg48BFwUvgd4BpgWnjguAp5098+Bf5nZeoITwGtVfXdBQQErVqyoqYsiIhLDzDZWtS+pnL6ZNTSzN4ASYDHwHvCJu+8Lq2wGOoTvOwCbAML9O4BjY8sTfEZERDIgqaDv7vvd/TSgI8Ho/KR0dcjMRpjZCjNbUVpamq6vERGJpFrN3nH3T4C/An2BFmZWnh7qCGwJ328BOgGE+5sD/44tT/CZ2O+Y6e6F7l7Ypk3ClJSIiByiGnP6ZtYG+MLdPzGzpsA3gckEwf97wJPAEOC58CPzw+3Xwv1/cXc3s/nA783sfoILud2AohQfj0itTV2whq2f7KlUlt+iKaO/fUqWeiSSPjUGfaA9MMfMGhL8ZfC0uz9vZmuBJ83sbuDvwCNh/UeAx8ILtduB7wO4+9tm9jSwFtgHjKxu5o5Ipmz9ZA/5LY6KK9udpd6IpFcys3fWAF9LUP4+QX4/vvwz4LIq2voF8Ivad1NERFIh0nfkDh48mKFDh2a7GyIiGRPpoF8bS5cuxczYtm1btrsiInLIFPRFRCIkMkF/9+7dDB06lLy8PPLz87nnnnsq7X/88cc5/fTTadasGW3btuWyyy5jy5ZgRumGDRs455xzAGjTpg1mVpEWeuGFF+jXrx8tW7akVatWDBw4kHXr1mX02EREkhWZoD927FgWL17MvHnzWLJkCX//+9956aWXKvbv3buX8ePHs3r1ap5//nm2bdvGFVdcAUCnTp2YN28eAG+//TbFxcVMnToVgLKyMsaMGUNRURFLly6lefPmfOc732Hv3r2ZP0gRkRokM2Wzztu1axePPPIIs2bNYuDAYHHPRx99lI4dO1bUGTZsWMX7r371qzz00EN0796dzZs307FjR1q1agVA27Ztad26dUXdSy+9tNJ3PfrooxxzzDEUFRVx1llnpfOwRERqLRJB/7333mPv3r307du3oiwvL49evXpVbK9atYrx48fzxhtvsH37dspXH/3ggw8qnRwStX3HHXewfPlySktLOXDgAAcOHOCDDz5I3wFJzoi/sUs3dUmui0TQr0lZWRkDBw7k/PPP57HHHqNt27Zs27aNfv361ZimGTx4MB07dmTGjBl06NCBI444gh49eii9ExHxN3bppi7JdZHI6Xft2pVGjRqxbNmyirKysjLeeustAN555x22bdvGPffcQ//+/TnppJMoKSmp1Ebjxo0B2L//y5uI//3vf/POO+9w2223cf7559O9e3d27tzJvn37EBHJRZEI+nl5eVx77bXccsstLF68mLfffpthw4ZVBPDOnTvTpEkTpk2bxvvvv8+CBQu44447KrVx3HHHYWYsWLCA0tJSdu3aRcuWLWndujUPP/ww69ev58UXX+SGG27giCP0B5SI5KZIBH2AKVOmcM4553DxxRdzzjnn0LNnT/r37w8E0zDnzJnDH//4R3r06MH48eO5//77K32+Q4cOjB8/np/97Gfk5+czatQoGjRowFNPPcWaNWvo2bMnI0eOZOLEiTRp0iQbhygiUqMaH5eYTYWFha4nZ0m63TZ3ecIF1+65qk+lskSrcW4s3UXvbm2r/ZxIppnZSncvTLRPeQiRJCVajfPd4h1Z6o3IoYlMekdERDTSF0k7PaRFcomCvkgCG0t3cdvc5QeVxad3kqGHtEguUdAXSWDvvv3K30u9pJy+iEiEKOiLiESI0jsiKZTstYBE9XRxVzJBQV8khZK9FpConi7uSiYovSMiEiEK+iIiEVJj0DezTmb2VzNba2Zvm9nosLyVmS02s3fD/7YMy83MHjCz9Wa2xsz+I6atIWH9d81sSPoOS0REEklmpL8PuMndewBnACPNrAcwDlji7t2AJeE2wLeAbuFrBPAQBCcJ4C6gD9AbuKv8RCEiIplR44Vcdy8GisP3O81sHdABuAgYEFabAywFbgnLf+fB8p3LzKyFmbUP6y529+0AZrYYGAQ8kcLjEalX9DhGSbVazd4xswLga8ByID88IQB8BOSH7zsAm2I+tjksq6pcRKqgxzFKqiUd9M0sD5gHjHH3T82sYp+7u5mlZGF+MxtBkBaic+fOqWhSpJL40fOhrqkjUhclFfTNrBFBwJ/r7s+GxVvNrL27F4fpm/KHym4BOsV8vGNYtoUv00Hl5Uvjv8vdZwIzIXiIStJHIpJAMg8+0Zo6EiXJzN4x4BFgnbvHPkNwPlA+A2cI8FxM+Q/DWTxnADvCNNAi4AIzaxlewL0gLBNJm/L0SOxr7779NX9QpJ5KZqR/JnA18KaZvRGW3QZMAp42s2uBjcDl4b6FwIXAemA3cA2Au283s4nA62G9CeUXdUVEJDOSmb3zMmBV7D4vQX0HRlbR1ixgVm06KCIiqaO1d6TeqCp/r4u0Il9S0Jd6Qw8uF6mZgr5IjkjlIxpFqqKgL5Ij9IhGyQStsikiEiEK+iIiEaKgLyISIQr6IiIRoqAvIhIhCvoiIhGiKZsidUiiufx6sIrUhoK+SB2SaC6/HqwitaH0johIhCjoi4hEiNI7kvMSrZ6pPLbIoVHQl5yXaPXMondLtDhZGugEW/8p6EudpMXJvpTKGT2JTrC6UFy/KOiL1HGa0SO1oaAvOSc+xaC0jUjqKOhLzolPMUQ1bSOSDpqyKSISIRrpi0SEZuYIKOhLliUKRMrhp0cyU1/1s6//agz6ZjYLGAyUuHvPsKwV8BRQAGwALnf3j83MgKnAhcBuYKi7rwo/MwS4PWz2bnefk9pDkbooUSBSDj9z4mf+JPuz118NdVcyI/3ZwDTgdzFl44Al7j7JzMaF27cA3wK6ha8+wENAn/AkcRdQCDiw0szmu/vHqToQEflSorn7qRzFaz5/3VXjhVx3fwnYHld8EVA+Up8DfDem/HceWAa0MLP2wEBgsbtvDwP9YmBQCvovIgmUj+BjX3v37c92tyQHHOrsnXx3Lw7ffwTkh+87AJti6m0Oy6oqFxGRDDrsKZvu7gQpm5QwsxFmtsLMVpSWlqaqWRER4dBn72w1s/buXhymb0rC8i1Ap5h6HcOyLcCAuPKliRp295nATIDCwsKUnUxE5NCk+/qAZNahBv35wBBgUvjf52LKR5nZkwQXcneEJ4ZFwD1m1jKsdwFw66F3W0QyJdnF7fQox7ohmSmbTxCM0lub2WaCWTiTgKfN7FpgI3B5WH0hwXTN9QRTNq8BcPftZjYReD2sN8Hd4y8OSz2nOfn1W6KTQ/x9ADoJZF+NQd/dr6hi13kJ6jowsop2ZgGzatU7qVc0Jz964k8EmtaZfVp7R0QkQhT0RUQiRGvvSK3F5+aVpxWpOxT0pdbic/PK04rUHUrviIhEiIK+iEiEKOiLiESIcvqSFroRSyQ3KejLYatqbZbe3dpWKtONWCLZp6Avhy3ZtVlEJPuU0xcRiRAFfRGRCFHQFxGJEOX0pVqahSOppDX3s09BX6ql5ZAllRJd9NcyHpmloB9RiUbwGnGJ1H8K+hGVaASvEZdkg1I+maWgLyJZpZRPZino10Na716iSmnLmino10OHut59VcspaKaO5IJkArrSljVT0I+AZIO5llOQXJHsek4K6LWnoB8BCuZS1yT7bzb+5KC/TGumoJ+jEv0pW7JjD22bN61UpnylRFn8yUGDmZplPOib2SBgKtAQ+K27T8p0H+qCqm6K6nXcsZXKit4tUR5eRJKW0aBvZg2B6cA3gc3A62Y2393XZrIfueZwljpQ6kakesneBxCVmT+ZHun3Bta7+/sAZvYkcBEQqaAf/49LDxwRSZ9k7wOIysyfTAf9DsCmmO3NQJ8M9yElDifnHv+PSwFeJLOiPD3Z3D1zX2b2PWCQuw8Pt68G+rj7qJg6I4AR4WZP4K2MdTA3tAa2ZbsTGaZjjgYdc+Yc5+5tEu3I9Eh/C9ApZrtjWFbB3WcCMwHMbIW7F2aue9mnY44GHXM05OIxZ/ohKq8D3cysi5k1Br4PzM9wH0REIiujI31332dmo4BFBFM2Z7n725nsg4hIlGV8nr67LwQWJll9Zjr7kqN0zNGgY46GnDvmjF7IFRGR7NKD0UVEIkRBX0QkQhT0RUQiREFfRCRCFPRFRCJEQV9EJEIU9EVEIkRBX0QkQhT0RUQiREFfRCRCFPRFRCJEQV9EJEIU9EVEIkRBX0QkQjK+nn5ttG7d2gsKCrLdDRGROmXlypXbcuUZubVSUFDAihUrst0NEZE6xcw2VrVP6R0RkQhR0BcRiRAFfRGRCMnpnL5Itk1dsIatn+ypVJbfoimjv31KlnokcngU9EWqsfWTPeS3OCqubHdSn9UJQ3KRgr5ImhzOCUMkXZTTT7FRo0YxYMCAWn3GzHjmmWfS0yERkRga6UfEhg0b6NKlC6+//jqFhYUpb3/AgAH07NmTadOmpbztXLOxdBe3zV1eqaxkxx7aNm96UL34kb5Itinoi8SIz8MnCtx79+0/qOzd4h30Ou7Yg8riJTphKM8vmRSJ9M6AAQO48cYbuemmm2jVqhVt2rRh6tSpfP7554wcOZIWLVrQuXNnHnvssUqfe/PNNzn//PNp2rQprVq1YujQoezY8eX/yPv372fs2LG0bNmSli1bMmbMGPbv31+pDXfnvvvuo2vXrjRt2pRevXrx+OOP16r/Bw4cYOLEiXTq1IkmTZrQq1cvnnvuuYr9GzZswMwOuns5Nm3UpUsXAE4//XTMrCIFNXToUAYPHszdd99Nfn4+eXl5XHPNNezZ82XgGzBgAKNGjarUdvnnyt+/+OKLTJ8+HTPDzNiwYUOtjjFXlOfhy1979+2v+UO1UH7CiH3FX+wVSadIBH2AuXPn0qxZM5YvX864ceMYM2YM3/3udznhhBNYsWIFQ4YMYfjw4RQXFwNQVlbGwIEDycvLo6ioiD/84Q+8+uqrDBs2rKLNX/3qVzz88MPMmDGD1157jf379zN37txK33v77bfzyCOPMH36dNauXcutt97K9ddfz4IFC5Lu+9SpU/nlL3/J5MmTefPNN7n44ou55JJLeOONN5Juo6ioCIAXXniB4uJinn322Yp9L774IqtXr2bJkiXMmzePP//5z9xyyy216l/fvn255pprKC4upri4mE6dOiX9eRHJnMgE/ZNPPpmf//zndOvWjZ/85Ce0bt2aRo0aMXr0aI4//njuvPNO3J1XXnkFgN///veUlZXx2GOP0atXL84++2xmzpzJs88+y/r16wH49a9/zc0338zll1/OSSedxNSpU2nXrl3Fd5aVlXH//ffz29/+lkGDBtGlSxeuvPJKrrvuOqZPn55036dMmcLYsWO58sorOeGEE5gwYQL9+vVjypQpSbfRpk2w9tKxxx5Lu3btaNWqVcW+hg0b8uijj9KzZ08GDhzI5MmTmTFjBmVlZUm13bx5cxo3bsxRRx1Fu3btaNeuHQ0bNky6byKSOZHJ6Z9yypc5UzOjbdu29OrVq6KsUaNGtGzZkpKSEgDWrVvHKaecQrNmzSrqfOMb36BBgwasXbuWNm3aUFxcTN++fSv2N2jQgD59+rBp0yYA1q5dy2effcagQYMws4p6X3zxBcmuHvrpp5/y4YcfcuaZZ1YqP+uss1i4cGHyP4BqnHLKKeTl5VVs9+3bl7179/Lee+9V+rmJSN0XmaDfqFGjSttmlrDswIEDNbYVG8CrU97Wn/70Jzp37lxtfw5FeT8aNAj+YHP3in1ffPHFYbdfrkGDBpXaTnX7maAbpUQCkQn6tdW9e3dmzZrFzp07K0b7r776KgcOHKB79+40b96c9u3bs2zZMs4991wgCLpFRUW0b98egB49etCkSRM2btxYUae2jjnmGL7yla/wyiuvcN5551WUv/zyy/To0QP4MnVTfj0COCjf37hxY4CDLjRDcMG6rKyMo48+GoBly5bRuHFjunbtWtF+bNsAq1evrvTXSuPGjRO2nSsS3ShV9G7JQTNpsjHNUjN6JJNqDPpm1gn4HZAPODDT3aeaWSvgKaAA2ABc7u4fWzD8nApcCOwGhrr7qrCtIcDtYdN3u/uc1B5O6lx11VXcdddd/PCHP2TChAl8/PHHXH/99VxyySUcf/zxAIwePZp7772XE044gV69evHggw9SXFxcEfSbNWvG2LFjGTt2LO5O//792bVrF8uWLaNBgwaMGDEiqb789Kc/5c4776Rbt258/etf5/HHH+dvf/sbq1atAqBp06acccYZTJ48ma5du7Jjxw5uvfXWSm20bduWpk2bsmjRIgoKCjjyyCNp3rw5APv27WPYsGHceeedfPjhh4wbN47rrruu4iRw7rnnMmbMGObPn8+JJ57IjBkz2LRpU6WgX1BQQFFRERs2bCAvL49WrVpV/AWSq6qaepkL/dCdu5IuyYz09wE3ufsqM2sGrDSzxcBQYIm7TzKzccA44BbgW0C38NUHeAjoE54k7gIKCU4eK81svrt/nOqDSoWjjjqKRYsWMWbMGHr37s2RRx7JRRddxNSpUyvq3HTTTXz00UcMHz4cgKuvvpqrrrqKdevWVdSZOHEi+fn5TJkyhRtvvJFjjjmG0047jZtvvjnpvvz4xz9m586d3HzzzWzdupUTTzyRefPmceqpp1bUmTVrFsOHD+f000+na9euPPjgg/Tv379i/xFHHMEDDzzAhAkTGD9+PP369WPp0qUAnH322Zx88smcc8457N69m0svvZT77ruv4rPDhg1jzZo1FTOXRo4cycUXX8y2bdsq6owdO5YhQ4bQo0cP9uzZw7/+9a+kr1vIwTT6l3Sx+FxtjR8wew6YFr4GuHuxmbUHlrr7iWY2I3z/RFj/H8CA8pe7Xx+WV6qXSGFhoevJWek1dOhQtm3bxvPPP5/trqTVbXOXHzSafnldMWd1b19tWTJ1MlW29ZPd3HNVn6oOUaSCma1094S33tfq728zKwC+BiwH8t29PNH7EUH6B6ADsCnmY5vDsqrKRUQkQ5IO+maWB8wDxrj7p7H7PPhzoXZ/MlT9PSPMbIWZrSgtLU1FkyIiEkpq9o6ZNSII+HPdvfxWzq1m1j4mvVMSlm8BYm/H7BiWbSFI8cSWL43/LnefCcyEIL2T9JHIIZk9e3a2uyAiGZTM7B0DHgHWufv9MbvmA0OASeF/n4spH2VmTxJcyN0RnhgWAfeYWcuw3gVA5SkmIimQaE6+VrwUCSQz0j8TuBp408zeCMtuIwj2T5vZtcBG4PJw30KC6ZrrCaZsXgPg7tvNbCLwelhvgrtvT8VBiMRKNCc/G1MxRXJRjUHf3V8GqroF9bz4gjC/P7KKtmYBs2rTQRERSR3dkStSR2juvqSCgr5IHaE7dyUVcvs+eRERSSkFfRGRCFHQFxGJEOX0pU7TnHyR2lHQlzpNc/JFakfpHRGRCNFIX6QO09x9qS0FfZE6THP3pbaU3hERiRCN9KXO0EwdkcOnoC91hmbqiBw+pXdERCJEI32RekYzeqQ6Cvoi9Yxm9Eh1lN4REYkQjfQlJ2mmjkh6KOhLTtJMHZH0UHpHRCRCNNIXiQDN6JFyCvoiEaAZPVJO6R0RkQhR0BcRiRAFfRGRCFFOXySidHE3mhT0Jet0I1Z2HM7F3US/M50w6gYFfck63YiV26o6Kffu1rZSmWYD1Q0K+iJSLZ2U6xddyBURiZAaR/pmNgsYDJS4e8+wrBXwFFAAbAAud/ePzcyAqcCFwG5gqLuvCj8zBLg9bPZud5+T2kORukD5+9yW6OKufj/1SzLpndnANOB3MWXjgCXuPsnMxoXbtwDfArqFrz7AQ0Cf8CRxF1AIOLDSzOa7+8epOhCpG5QqyG2JLu4m+/vRbKC6ocb0jru/BGyPK74IKB+pzwG+G1P+Ow8sA1qYWXtgILDY3beHgX4xMCgF/ReRHFF+woh9xf9VJ9l3qDn9fHcvDt9/BOSH7zsAm2LqbQ7LqioXEZEMOuwLue7uBCmblDCzEWa2wsxWlJaWpqpZERHh0IP+1jBtQ/jfkrB8C9Appl7HsKyq8oO4+0x3L3T3wjZt2hxi90REJJFDDfrzgSHh+yHAczHlP7TAGcCOMA20CLjAzFqaWUvggrBMREQyKJkpm08AA4DWZraZYBbOJOBpM7sW2AhcHlZfSDBdcz3BlM1rANx9u5lNBF4P601w9/iLw1LPaHqmSO6pMei7+xVV7DovQV0HRlbRzixgVq16J3WapmeK5B7dkSsiEiEK+iIiEaKgLyISIVplU1JCF20lES3NkHsU9CUldNFWEjmcB7VIeii9IyISIQr6IiIRoqAvIhIhyulLremirUjdpaAvtaaLtiJ1l9I7IiIRoqAvIhIhSu9ItZS/F6lfFPSlWsrfi9QvCvoiklFamiG7FPRFJKO0NEN26UKuiEiEaKQvFXTRVqT+U9CXCrpoK1L/KehHlEb1ItGkoB9RGtWLRJOCvohkXfw0Tk3hTB8FfRHJuvhpnJrCmT6asikiEiEa6YtIztFdu+mjoB8BmqkjdY3u2k0fBf16pqoA37tb20plmqkjEk0K+vWMpmJKfaWUT2oo6ItInZAo5VP0bolOBLWkoF9HJErblOzYQ9vmTSuVKVcvUaLcf+1lPOib2SBgKtAQ+K27T8p0H+qiqtI2vY479qAykShTGqh6GQ36ZtYQmA58E9gMvG5m8919bSb7kWviR/EawYscumTTQIn+P4vCySHTI/3ewHp3fx/AzJ4ELgLqZdCvTUomdnaNRvAiqZXoRJDo/7MonBwyHfQ7AJtitjcDfTLch1pJNnAnE8xBAV0kl6X65JCLJwxz98x9mdn3gEHuPjzcvhro4+6jYuqMAEaEmz2BtzLWwdzQGtiW7U5kmI45GnTMmXOcu7dJtCPTI/0tQKeY7Y5hWQV3nwnMBDCzFe5emLnuZZ+OORp0zNGQi8ec6QXXXge6mVkXM2sMfB+Yn+E+iIhEVkZH+u6+z8xGAYsIpmzOcve3M9kHEZEoy/g8fXdfCCxMsvrMdPYlR+mYo0HHHA05d8wZvZArIiLZpYeoiIhESM4HfTO7zMzeNrMDZpZTV8FTzcwGmdk/zGy9mY3Ldn/SzcxmmVmJmUVmWq6ZdTKzv5rZ2vDf9ehs9ymdzOxIMysys9Xh8Y7Pdp8yxcwamtnfzez5bPclVs4HfYJ5+pcAL2W7I+kUs0TFt4AewBVm1iO7vUq72cCgbHciw/YBN7l7D+AMYGQ9/z1/Dpzr7qcCpwGDzOyM7HYpY0YD67LdiXg5H/TdfZ27/yPb/ciAiiUq3H0vUL5ERb3l7i8B27Pdj0xy92J3XxW+30kQFDpkt1fp44Fd4Waj8FXvLySaWUfg28Bvs92XeDkf9CMk0RIV9TYYCJhZAfA1YHkNVeu0MM3xBlACLHb3en28oV8DNwMHstyPg+RE0Dez/zOztxK86vVIV6LLzPKAecAYd/802/1JJ3ff7+6nEdyB39vMema5S2llZoOBEndfme2+JJITD1Fx9/Oz3YccUOMSFVI/mFkjgoA/192fzXZ/MsXdPzGzvxJcx6nPF+/PBP7TzC4EjgSOMbPH3f0HWe4XkCMjfQG0REUkmJkBjwDr3P3+bPcn3cysjZm1CN83JXiWxjtZ7VSaufut7t7R3QsI/j/+S64EfKgDQd/MLjazzUBfYIGZLcp2n9LB3fcB5UtUrAOeru9LVJjZE8BrwIlmttnMrs12nzLgTOBq4FwzeyN8XZjtTqVRe+CvZraGYGCz2N1zagpj1OiOXBGRCMn5kb6IiKSOgr6ISIQo6IuIRIiCvohIhCjoi4hEiIK+VGJms3NtVcBYZtbSzLaaWdds96UuMDM3s+/Vov5QM9tVc80a2ykys0sPtx1JPQV9Sbtwed3ZZrbGzL4ws6VV1DvbzFaa2Wdm9r6Z3ZCg2m3AQnd/L62djq6ngK+moJ2JwCQzU4zJMfqFSCY0BD4DpgELElUwsy4Ej9F8lWARsnuB/4kdLZrZUcBwgjtaJQ3cfY+7l6SgqYVAM4KlwiWHKOhLtcysiZn9OkypfGZmy8zsrLg63w4f/vKZmb1kZt8P0woFAO5e5u43uPtMgtVDE7kB+NDd/ztcTvthYA4wNqbOhQTL8r4S9/09zGyBme0MH8ryhJm1C/edaGZlZjYkpv4gM9trZn3D7dlm9ryZ3R4e5y4zezRcNqD8M2ZmN5vZe2a2x8zeNLMfxOwvCI/5UjNbbGa7wwelfDOmTiMze8DMPjSzz81sk5lNitnf2Mwmh3cn7zaz181sYA2/HzOzm8zs3bDNzWZ2bzX1J4W/qz1mtsHM7jOzI2P2V0rvmNnPLVj8cEhYvyz82TQ2sx+Fx/BvM7s/dlTv7vsJAv8V1fVfMk9BX2pyH/BfwDCCEfibwAtm1h7AzDoDzxKM4E8FHgg/U1t9gT/HlS0CCi1YoAygH7DSY24jD/vxEsECXr2B84E84DkzaxA+i+H/EfzV8FUza0Pw8JZfuPtrMd91dtj/84BLgQuAyTH77wauBUYSPOTmXmCGmX07rs+/IPgZnEqw7MCTFqyoCfBj4GKC9Vi6EfxcY58V8WjYjyuBngQnvT+Z2alV/MwA7gHuCPtzMnAZlZfojldG8LvsDvwo7MvPqqkPUEDwbIfBBA80uoxgXajTCX5Ow4H/Do8tVlF4PJJL3F0vvSpeBAHx+fD90cBe4Icx+xsC7wF3h9v3EiweFtvGbQQj8oIE7U8DliYo/ydwZ1xZ/7Cd9uH2H4E5cXUmAEviylqGn+sdU/YHYBnByekVoGHcMX8C5MWU/YDgqU9Hh689QL+47/k1wfUFCAKjA9fH7O8Qlp0Vbj8ALCFc/iSura4Ea693jiv/I/BgFb+rPIK02Q3V/D4d+F41+28geHhP+fZQYFfM9s/DY28eU/YMUAo0jilbCkyLa/s/w2M6Itv/rvX68pUTSytLzupK8KSjinSKu+83s9cIRrsAJxGMaGOl6yEZTYGtcWVfB/pXMeOkK8FoE4LR6D8IRsOneJB+iLXGv3zCEwQLwTUO22hCsETuC2YWu1hVI2BDfDsx7z8M/9s2/O9sYDHwTzP7M0H643/d/QDwH4ABa80str0mwF8SHBsEv4MmBCeSpIQzecYAxxOcNBqGr+p84O47Yra3Av/04AlvsWVtK3+MPQTHdCRw2DOCJDUU9OVQpXqlvo+A/LiyfIJnym4Lt7cRjOJjNSAYvY/lYLEniJ5A8/B9B+BftehbeRr0O8AHcfu+qGrb3T0M4A3C7VXhdY6BBGmkOcDqMO/fgOBnenqCNvfUoq9VsuDZtE8C4wlSXp8QjMan1PDR+P54FWXxJ49WwGdxJ1PJMgV9qc57BOmdM8P35Q9w7wv8PqzzDgc/y7f3IXzXaxycE/4msMLdywPM3wnSD7FWAZcDG2PqVWLBeu6PEQS3psBjZnaqV35iVS8zO9rdy8LtMwiO/T2CgPw5cJy7VzXqTooHz8V9BnjGzGYTpJyOD4/NgHbu/tckm1sX9us84N0k6p8JbHH3ieUFZnZc8r2vtZ4Evx/JIbqQK1UKA+BDwGQzu9DMuofb+cCDYbXfAF3NbEo4U+YS4PryJsrbCmfYnAa0BvLM7LRwm5h2OoQzhbqb2XCCAB87Cl0EdDezY2PKphOM4J8ysz7hxdrzzWymmTWLabsUuBO4BdgZfi7WEcAsMzs5HHlPAh72YObRzrAfU8xsmJkdH/b/BjMbkezP08x+YmZXhMd3PMEF20+Bze7+T2AuMNvMvhceR6GZjQ1/pgcJ+zUVuNfMrjGzrmbW28xurKIL/yT4GV8Vtn8j6Z1d0w94IY3ty6HI9kUFvXLrRcyF3HC7CcEFy60Eo8plhBcmY+oMJggonwF/A64hCPj5MXU2hGWVXnHtnE0wMvycIP1y0AVKgr8IRsaVdSMYPX9MkAr5B/A/BDn5q4HdwEkx9XuG9b4fe8wEJ4USgvzzHOComM8YwQyVtWH/Sgny898M9xeEx1QY17eKC6nAdeHx7SQI9i8C34ip24jgwun7BH9lfEQwS+br1fy+GgDjYj6ziWBm0kHfH27fG/Z9F8Gsqxtjfw8kvpD7Vtx3HnQxniBt9EzMdoewPx2z/W9ar8ovPURFUs7MRhPMqmnhKf4HZmaDCEa3Pfzgi7GH2uZsoLW7D05FewJm9kuCGT9J/yUkmaGcvhw2MxtJMIOnlCAXfgcwO9UBH8DdXzCz6QQPjt+Y6vYlZUqo+QKxZIGCvqTC8QRz848luOP2NwQj/bRw9wfS1bakhrv/Mtt9kMSU3hERiRDN3hERiRAFfRGRCFHQFxGJEAV9EZEIUdAXEYkQBX0RkQj5/7+vUsvtRgVCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Histogram of raw data, and of a synthetic sample from the probability model\n",
    "\n",
    "fig,(ax1,ax2) = plt.subplots(2,1, figsize=(6,4), sharex='col')\n",
    "ax1.hist(np.log10(mpexpenses), bins=75, alpha=.6, fc='steelblue', ec='steelblue')\n",
    "ax2.hist(np.log10(xsample), bins=75, alpha=.6, fc='steelblue', ec='steelblue')\n",
    "ax1.set_title('data', x=.05, y=.75, fontsize=14, ha='left')\n",
    "ax2.set_title('model output', x=.05, y=.75, fontsize=14, ha='left')\n",
    "ax1.set_xlim([-1,4.5])\n",
    "ax2.set_xlabel('log10(expense claim)', fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.6 Fitting a generative model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 53774.555345\n",
      "         Iterations: 573\n",
      "         Function evaluations: 900\n"
     ]
    }
   ],
   "source": [
    "φ = scipy.stats.norm.pdf\n",
    "\n",
    "def logPr(y, θ):\n",
    "    q,μ1,μ2,τ1,τ2 = θ\n",
    "    p = np.exp(q) / (1 + np.exp(q))\n",
    "    σ1,σ2 = np.exp([τ1,τ2])\n",
    "    lik = p*φ(y, loc=μ1, scale=σ1) + (1-p)*φ(y, loc=μ2, scale=σ2)\n",
    "    return np.log(lik)\n",
    "\n",
    "y = np.log10(mpexpenses)\n",
    "initial_guess = [1, 1.8,2.7, np.log(1.2),np.log(0.5)] # from eyeballing the histogram\n",
    "θhat = scipy.optimize.fmin(lambda θ: -np.sum(logPr(y,θ)), initial_guess, maxiter=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAADUCAYAAACMCNgJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAr+klEQVR4nO3de5xN9frA8c9jbkQkt9yJQi6hCSVFpVxKp5SjVJR7SC6pTuXIoSJRiTqKhESXUz+dFCrVKdSg3EMuFXItl2HMmJnn98das9tzYTb2zNqz9/N+vfbLrO/67rWfNcNjzXd91/cRVcUYY0zBV8jrAIwxxgSHJXRjjAkTltCNMSZMWEI3xpgwYQndGGPChCV0Y4wJE9FefXDp0qW1WrVqXn28McYUSCtWrNivqmVy2udZQq9WrRrLly/36uONMaZAEpFfTrbPhlyMMSZMWEI3xpgwYQndGGPChCV0Y4wJE5bQjTEmTHg2y8WY/DAvYTsHEpMztZUqFkeHy6t5E5AxecgSuglrBxKTKVeiSKa2PYeSPIrGmLxlQy7GGBMmAkroItJGRDaKyM8i8mgO+yeIyI/ua5OIHAx6pMYYY04p1yEXEYkCJgGtgR1AgojMU9X1GX1UdZBf/wFAozyI1RhjzCkEcoXeBPhZVbeqagowB7jlFP3vBN4ORnDGGGMCF0hCrwj85re9w23LRkSqAtWBL84+NGOMMacj2DdFOwPvqWpaTjtFpJeILBeR5fv27QvyRxtjTGQLJKHvBCr7bVdy23LSmVMMt6jqFFWNV9X4MmVyXP3RGGPMGQokoScAF4lIdRGJxUna87J2EpHaQElgaXBDNMYYE4hcE7qqpgL9gQXABuAdVV0nIiNFpINf187AHFXVvAnVGGPMqQT0pKiqzgfmZ2kbnmV7RPDCMsYYc7rsSVFjjAkTltCNMSZMWEI3xpgwYQndGGPChC2fa0wObB11UxBZQjcmB7aOuimIbMjFGGPChCV0Y4wJE5bQjTEmTFhCN8aYMGEJ3RhjwoTNcjGG7NMUN+86mG2WizGhLihFot0+nURkvYisE5HZwQ3TmLyVMU0x45WUkup1SMactqAUiRaRi4DHgOaq+qeIlM2rgI0xxuQsWEWiewKTVPVPAFXdG9wwjTHG5CZYRaIvBi4WkW9FZJmItAlWgMYYYwITrJui0cBFQEucmqNfi0h9VT3o30lEegG9AKpUqRKkjzbGGAOBJfRAikTvAL5T1RPANhHZhJPgE/w7qeoUYApAfHy8laozBZ4t4mVCSbCKRH+Ic3WOiJTGGYLZGrwwjQlNWWfHlCtRJFuCNya/5HqFrqqpIpJRJDoKmJZRJBpYrqrz3H03iMh6IA14WFUP5GXgxpypzbsO8sbijdnabN65KeiCUiRaVRUY7L6MCWlJKanZkvfq7fs9isaY4LEnRY1xJR07yqqEJWxet4oVq9bwWXoSaWmpxMYVpky5CpQsX43ahW+lSZMmREVFeR2uMdlYQjfG9c1n83nhqaGn7PPulOcoW7YsnTt35vHHH8+nyIwJjC3OZYyrSYtrEZFc++3du5epU6cSExOTD1EZEzi7QjcRJ/HP/UwY8SLdH3qc4ueV9LWXKFmK6zvcQbFixTlW5AKuuewSYmJjOXbsKLt3/Mo3S5bx67rvOfTHPprdcCsf/rjXbqaakGIJ3USUpYsXMOufQzl+9DBRUVE8+OSYTPsfGv4cAItW/call1bOtC+uViuurVeB1QlLqFi1OmVLFMl0M3XFkq+oUKU6W47EZJtFY3PTTX6whG4iQlpaGjMnj+Pd6ZN9bQv/by63d+tLhcrVAj5OVFQUjZq1yNa+89dtPPPoAxSSQtzQawRX1emYab8VmDb5wRK6CXsnUpIZ+/iDLPniU19b6XLlGTJywmkl85NRVcY98RBJRxMB+L8XH+bC84Rr29921sc25nTYTVET1pKPJzFycM9Mybxa/Wa8/PanNIi/IiifISL0GfYUpcuVByA9LY3nhw/i3emv4DyiYUz+sIRuwtbhw4eZ8GhPVi79ytd2y53302HgGM4tcV5QP6tWvYY8/8YHVK1Ry9c2feKzzHxlXFA/x5hTsYRuwlJSUhJt27Zl0+rlvrYuvQfRc8hwChXKm4eCSpcrz9ip71KxVkNf29ypL/P+jH/nyecZk5UldBOWBgwYwJIlS3zb3R96nLt6PRTQPPOzUezcEtw6eDzxzVv52qa9+DRff/xunn6uMWAJ3YSpxx57jBo1agDQY/CT3HZPr3z77OiYWB4b8wp1GzXxtb05fjj/+c9/8i0GE5ksoZuwVKNGDZYsWUKXAU9wa5ce+f75hYsU4Z8TplKjdj3AmQkzZ86cfI/DRJaApi26JeVexFk+93VVfTbL/m7Ac/xV+OJlVX09iHEac9rKli3Ldbfe7dnnFz23OP96eQZD77uNhle15u2Zr/j2WWEMkxdyTegiEgVMAlrjVCZKEJF5qro+S9e5qto/D2I0Jldbt27l6NGj1K9f3+tQMilRshQvvvUxR1KjMq3QmFEYw589fGTOViBDLk2An1V1q6qmAHOAW/I2LGMCd+zYMW677TaaNWvG22+/7XU42ZxTtJjXIZgIEciQS0XgN7/tHUDTHPp1FJGrgU3AIFX9LWsHKxJt8kK/fv1YtWoVAPd27cpOuYBS5SoAoVmJKC0tjY9mvUqD+g1ocvV1Xodjwkiwbop+BFRT1QbAIuDNnDqp6hRVjVfV+DJlygTpo00ke/vtt5k+fbpvu8/DT3HJxTV89T2TUlK9Cy4Hu3fvpnXr1nww7QUmjBjCgX17vA7JhJFAEvpOwH/ZuUr8dfMTAFU9oKoZd3heBy4LTnjG5Gxewnaee/sLuvfs7Wurd1Vb2tx2l4dR5S46OpqNG52VGA8f+pNxTwwkLS3N46hMuAgkoScAF4lIdRGJBToD8/w7iEh5v80OwIbghWhMdvsOHWPGuH+QdPQIABdUrMLVdw7K8weHzlbp0qV56623fHGuXr6U9958JZd3GROYXBO6qqYC/YEFOIn6HVVdJyIjRaSD2+1BEVknIquAB4FueRWwMQCfzJ3K2pXfAVAoKoqho14gtsg5HkcVmJYtW9K+Sx/f9ux/v8C2zXYNZM5eQGPoqjpfVS9W1RqqOtptG66q89yvH1PVuqp6qaq2UtWf8jJoE9lWr17Nh9Mn+rY7dx9AnQYFa5Tvlq79qFWvEQCpqSeY8M8hpKae8DgqU9DZk6KmQElNTeX+++8nzU1+teo1onP3AR5HdfqioqIZNGIcMbFxAGzZuI75s6d4HJUp6CyhmwJlwYIFrFixAoCY2DgGjRhHVHTBrNNSuXpN7uk7xLf90cxXfNMvjTkTltBNgdK+fXs++eQTzi9bnrt6DqRy9Zpeh3RW/talB7UbNAYgLS2V++67j9TU0JpqaQqOgnlpYyJOprVP4qrTbfQsrqlf1duggiAqKopB/xzHgLvakpKczIEDB/jll198K0UaczosoZsCIevaJ6uj4oiOifEwouCpVK0GXXoP5vc9e3lv6osULVrU65BMAWUJ3YS8gwcPkpx0DELsEf5gur1rH/YcSrJkbs6KjaGbkDdkyBAev689SxcvsKLLxpyCJXQT0r799lumTZvGH3t/Z9TQXmxc+4PXIZ2RzbsO8sbijb7X5l0Hc33Pp59+ytGjR/M+OBM2LKGbkJWamkrfvn1921e0vJHa9Rt7GNGZS0pJ9S0YltuiYTt27KBjx460bduWkSNH5mOUpqCzhG5C1sSJE1mzZg0AsYWL0GvocI8jyh+fffaZr/7o+PHjWbt2rccRmYLCEroJSTt37mT48L8SeId7HqBs+UoeRpR/unbtytVXXw04v6U8+OCDdu/ABMQSuglJQ4cOJTExEYA6depww+1dPY4o/4gIkydP9pWsW7x4MR988IHHUZmCIKCELiJtRGSjiPwsIo+eol9HEVERiQ9eiCbSfPPNN8yZM8e3PXnyZKJjYj2MKP/VrVuXfv36+baHDBlCUpLVHDWnlmtC9ysS3Ra4BLhTRC7Jod+5wEDgu2AHaSJHeno6AwcO9G136tSJli1beheQh0aMGEGpUqUA2L59O+PHj/c4IhPqglkk+l/AGOB4EOMzEWb69OmsXLkSgMKFCzN27FiPI/JOyZIlGTVqlG/76aefZseOHR5GZEJdIAk9pyLRFf07iEhjoLKqfhzE2EwESktLo3jx4gA8/PDDVK1a8NdrORs9e/akQYMGABw7doxHHnnE44hMKDvrm6IiUggYDwwJoG8vEVkuIsv37dt3th9twlDPnj3ZtGkTgwcPtuSFs3jXSy+95NuePXs23377rYcRmVAWjCLR5wL1gC9FZDvQDJiX041RVZ2iqvGqGl+mTJkzj9qEtXLlyvH888/buiaua665hjvuuAOAjh07UqlSZEzfNKcvkMW5fEWicRJ5Z8BXWl1VDwGlM7ZF5EtgqKouD26oxkSu5557jr59+9KqVSuvQzEhLNeErqqpIpJRJDoKmJZRJBpYnlFX1JgzlZCQQOHChalfv77XoYSsqlWrRvz9BJO7gJbPVdX5wPwsbTk+h62qLc8+LBMpMmqErl+/nt69ezNq1CjOP/98r8PyRMYCXv5KFYujw+XVvAnIFDj2pKjx1GuvvcbatWtJT09nxowZEf3wTNYFvMqVKPJXlaYsEhMTefzxx9m4cWOO+01ksgIXxjMHDx7kySef9G0/9thjVKxY8RTvMAALFy6kW7du/P7776xatYr//ve/XodkQoRdoRvPjBo1igMHDgDOGPHgwYM9jqhgKF26NLt37wbg448/5pNPPvE4IhMqLKEbT/z888+Z5lePHTuWIkXCt8RcMDVu3Jju3bv7tgcNGkRKSoqHEZlQYQndeGLYsGGcOHECgCuvvNI3z9oEZvTo0b4najdu3MikSZM8jsiEAkvoJt99+eWXmZaDnTBhAiLiYUQFT9myZfnnP//p2x4xYgR79+71MCITCiyhm3yVlpaWaaz8mja3sO5oiUz1NuclbPcuwAKkf//+1KpVC4DDhw/zxBNPeByR8ZrNcjH5asaMGfzwg1PouUiRItzc7SHKlcg8dv7Nht+zTdfbvOtgtn6RLjY2lgkTJtCuXTsAXn/9dfr27UujRo1y7D8vYXu276vNcw8vdoVu8lX9+vW58sorAacq0flly2frk9N87FMVVY5kbdu29SV0VT1luboDickBz3M3BZNdoZt8FR8fzzfffMP7779PmzZteDdhZ+5vimCBPD06fvx4Fi5cSGpqKt988w0LFy7kxhtvzOdITSiwK3ST70SE22+/nWLFinkdSsgL5OnRWrVq8eCDD1K8eHGmTJnCDTfc4FG0xmt2hW7yTU5juDY2HhzDhw9n0KBBtrRuhAtKkWgR6SMia0TkRxH5JqeaoyZyJSQkMGnSJPYeOmpj43mkRIkSlsxN0IpEz1bV+qraEBiLU8HIGFSVgQMH0r9/f4b3uIXN61d7HVLESElJ4ejRo16HYfJRUIpEq+phv82iQM632U3EmTt3LkuXLgVg365fKXZuCY8jigwrVqzg8ssvp3///l6HYvJRIGPoORWJbpq1k4j0AwYDscC1QYnOFGhJSUmZ6oJef+s9lK9sRRry2po1a2jatClpaWmsXr2a5s2b06NHD6/DMvkgaDdFVXUSMElE7gKeALpm7SMivYBeAFWqVAnWR5sQNWHCBH799VfAWSHwprv7ehxRZKhXrx533XUXM2fOBKBfv35Ur14dCp16jP3o0aPMmTOHbdu2ERcXR/v27WncuHF+hGyCJJCEnluR6KzmAK/ktENVpwBTAOLj421YJozt2rWLZ555xrc9cuRIChc718OIwlvWGUQ393yMVatWsXr1alJSUujQoQMPjJxEuVYn/+V5+PDhjB8/PtN2ly5dmDRpEiVK2FBZQRDIGLqvSLSIxOIUic5UR1RELvLbbA9sDl6IpiB65JFHSExMBKBu3br07NnT44jCW9anQH/49TB3PzqekmUuAODYsWM8/3B3pk8cw4F9e0hPTyc9PT3TMfr370+hQplTwltvvcUdd9xBaqrNRioIck3oqpoKZBSJ3gC8k1EkWkQ6uN36i8g6EfkRZxw923CLiRzffvsts2bN8m23vW8oM/+3hc27DnoXVIRJSkmlbq2LGPPvtylZqgwAqum8O30y97Zpws2XV2fDqhWZ3lO9enWGDBnC0KFDadOmja990aJFPPpottnKJgQFNA9dVeer6sWqWkNVR7ttw1V1nvv1QFWtq6oNVbWVqq7Ly6BN6EpLS2PAgAG+7ebXtaVVq2ttzrlHKla9kBdm/Zf6lzXLtm/l0q+ytY0dO5YWnfrRadgL3HzPA772559/npUrV+ZprObs2aP/JqhmzpzpW00xJjaOHoNsSVevlS57AaNfmU2bnsOpf1kzoqNjAPhp7Q859s8Yvun14MPEN2/pa7fleUOfPfpvgqpTp05s3bqVsWPH0qZzT8qWt6cXgy2nBbtyW0IhKiqK2lfcwIA+Tum61BMniIqOZu/h4yd9T6FChej+0OOsXPo16enpfPLJJ/zvf/+jRYsWwTkRE3R2hW6C6pxzzmHkyJGsX7+etp1t7nNeCMbywtExMQFViapy4cW0anurb3vChAmnHa/JP5bQTZ648MILiY0r7HUYJgj+3mMAVWrWYdKkSUybNs3rcMwp2JCLOWsZBRWsLmh4qlilOiOmfMB9rWp5HYrJhV2hm7M2bdo02rVrx9atW70OxZiIZgndnJV9+/YxbNgwPv30U+rWrcuXX37pdUjGRCwbcjFnZejQofzxxx8AlC9fniZNmngckTkdpztjRlVZuXIl27dvp2PHjvkRojkNltDNGVu8eDEzZszwbU+ePJlzzjnHw4jM6cqYMeNv9fb9OfbdsWMHLVu2ZMuWLZQqVYoOHToQExOTH2GaANmQizkjycnJ9OnTx7fdqVOnTI+Lm/BToUIFjh935q0fOHCAL774wuOITFZ2hW7OyJgxY9i0aRMAxYsXp03Xwaf9sIspODbvOsibX22m7hXXs/O9NwGneMmNN97ocWTGnyV0c9p++uknnn76ad/26NGjoch5Af/qbgqejKGZG2+6lYVuQn/3/Q+48q4hREU5aaRUsTg6XF7NwyiNDbmY05KWlsZ9991HcrKz9nZ8fDx9+1rhikhxcd1LKVXWWZI38fBB9m1Z43ta1X89duONgBK6iLQRkY0i8rOIZFtHU0QGi8h6EVktIp+LiNUZC1MTJkxg2bJlAMTExDB16lSioqI8jsrkl0KFCnFFq7+GWZYs/tTDaExWuSZ0EYkCJgFtgUuAO0XkkizdfgDiVbUB8B4wNtiBmtCwc+dfxaqefPJJGjRo4GE0xgtXtvrr5vfSxQuyFcow3gnkCr0J8LOqblXVFJwSc7f4d1DVxap6zN1chlOmzoShCRMm8Nlnn/G3v/3Nih5EqHqNmlC4mFOS7sDe3Wxat8rjiEyGQBJ6ReA3v+0dbtvJdAc+yWmHiPQSkeUisnzfvn2BR2lCynXXXccHH3xgc5AjVFR0NDUaXuXbtmGX0BHUWS4icjcQD1yT034rEh3ashYaBpu5YHJWo/HVrPvmYwCWfvEp9w2w39ZCQSAJfSdQ2W+7ktuWiYhcDzwOXKOqdru7AMqoVONv54HDDBs2jIcffpgyZcp4FJkJNVXqxlO7fiMaNbuaK6+1B8pCRSAJPQG4SESq4yTyzsBd/h1EpBHwb6CNqu4NepTGM/83fSIfz57CzJkzmTFjBq1bt/Y6JBMComPieH76h16HYbLIdQxdVVOB/sACYAPwjqquE5GRItLB7fYcUAx4V0R+FJF5eRaxyTerEpYw/+3XANi9ezerVtnNL2NCWUBj6Ko6H5ifpW2439fXBzku47HDB//k+eGDfMUrrr/+egYPHuxxVMaYU7EnRU02qspL/3qEA3t3A1C6dGlmzJhBoUL218XkbP/e3fyyeb3XYUQ8+xdqsvnPzCks/XKBb3vatGmUL1/ew4hMqNr+80aGdPsbXds2ZdaLI70OJ+JZQjeZrF35HdNfHuPbvv62e7n55ps9jMiEspKlyrBp/WoAtqz/kV27dnkcUWSzhG58Dh7Yy7OP9iM9LQ2A2g0a06n3UI+jMqGsRMnzqd+4qW977ty5HkZjbPlcA0BKSgqvjhzMnwecJ3iLn3c+jz07mbSY2Gx9c3oAydY+j1zXtLmFVQlLAJgxYwaDBg3yOKLIZVfoxueCKtUBEBGGPT2R0uXK+2pO+r+WbdrjWzI145WUkupx9MYrV13Xjti4OAB+/PFHVq9e7XFEkcuu0A0AsbGxdB08kjp16nLiRDKNmjprdZxOzUkTmYqeW5xmLW/k6wXO4yczZsxg3LhxHkcVmewK3fiICB3uvI+O9/bJvbMxfq5rf5vv6zfffNNXe9TkL0voEWz79u2kuTdAjTkbjZq24PyyztTW/fv3884773gcUWSyhB6htmzZQtOmTfn73/9uV1PmrEVFR9OqQ2ff9sSJE31PGZv8Ywk9Au3Zs4cbb7yRvXv38v7773PnnXd6HZIJA1e370Sce3N0xYoVbNq0yeOIIo8l9Ajz559/0rZtW7Zs2QJA4cKFGTrU5pqbs3duiZJ06dKFO+64gzVr1lCrVi2vQ4o4Ac1yEZE2wItAFPC6qj6bZf/VwAtAA6Czqr4X5DhNEBw8eJDWrVvzww8/AE7B37lz59K8eXOPIzPhYsqUKVY03EO5JnS/ItGtccrPJYjIPFX1X4nnV6AbYJd6IerQoUPccMMNrFixwtfW7/FnOXBuLd5YvBGwh4PM2dm86yAzvv45U5tVvMpfgVyh+4pEA4hIRpFoX0JX1e3uPiv/HYIOHDhAu3btSEhI8LW99tprRNVokSmB2/xyczZyemZhz6EkUlJSiI3N/sSxCb68KBJ9UlYkOv/t2LGDFi1a8P333/vaXn31VXr06OFhVCYSJB8/zofTX6J27docOnTI63AiQr4+KWpFovOH/1orLzzWmw0bNgDOg0OTJ0+md+/eXoZnIoCq8mS/u1n3o/Nb4QMPPMCsWbMQEY8jC2+BXKEHVCTahI6MYs/lShRhyIixlC1fiajoGGbPnk2fPvYUqMl7GU8dZ5g9ezbTpk3zMKLIEEhC9xWJFpFYnCLRVjO0gChd9gJGTZ7FQ0+/SufOnXN/gzFBctX17WnR7nbfdt++ffniiy88jCj8BaVItIhcLiI7gDuAf4vIurwM2uQsOTmZhQsXZmuvWKU6deNtaqLJf0079qdi9YsAOHHiBDd1uIVx0z/0NqgwFtCDRao6X1UvVtUaqjrabRuuqvPcrxNUtZKqFlXVUqpaNy+DNtlt3ryZK6+8kjZt2vDj0sVeh2MMAKkSzaiJb1KqTDkAko4m8sQDXfjqq688jiw82ZOiBZyq8tZbb9G4cWNWrlyJqjJl9FD2/r7D69CMAaBs+Yo8NfFNipcoCUBy0jFat27N1KlTPY4s/Nh66AVI1kpBBw/sZe7Ef/Hd14t8bbGxsdx2/0OUueCMZpYakyeqX1SHZ197hyce6MIf+/dy4sQJ39TZ7t27B3SMnCpl2YNLmVlCL0AyZq+kp6fz2Ufv8vqEURw9cti3v2bNmsydO5dVh4pmmx6WUXkoa5s9GWryS9UaF/P89A8YPrA7v235ierVq9OpUycge7LOKVFn/P33t+dQUp7HXZBYQi9gflrzA1PGPcXGtT9kau/duzfPPfcc5557LquyJG6wykMmNJQtX4lOj73CsndfokXbjry3fBfgXFxcVac8qoqIWKI+Q5bQC5Cli+bx2jPDMrWdV7YiPR8ZTZ1GzTL947ArbxOq0grF8I+nX8jUtnr7flSV0UN7U/OS+lzZ/i5vgivgLKEXIPWbtKBY8RIkHj5EdEwst3bpQdmmt9KyycWZ+tmVtymIln21iKVfLmDplwv4cPY0jvzjUR544AGKFi3qdWgFhs1yCVEJCQmsW5d5On+xEiW5u89grmh1I6++9zndBjxCTJxdiZvw8PlHf626feTgHwwbNowLKlamw739mPDe/9i866B3wRUQltBDSGJiIm+88QbNmzenSZMmPPnkk9n63NSpK0+Mm0L5SlU8iNCYvPPYmMn0f/yZTDO0Eg/9yUczJ/Nw52v5z8tPsiphCenptqjrydiQiwf87+innkhhww/f8eWn81i/7HOSjx/z9fvwww/ZunUrF154oa/NFjcy4SoqOpq2t93F9TffzsTJU1iz6G3f8xRpaalsXLaIfyxbRNnylXjkmZepXb+RxxGHHkvoHtj22y52bljOyqVfs+zrRZmmHmaIjo6hdrPreW/ZNsr8cgKwm50mMsTExFK/5S0M7NeHJYsX8NGcN3yrNgL8sW8PFSpXy/SetLQ00tPTiYmJyedoQ4sl9DyWdX6tqvJE91tIPJjzjcvK1WvSukMnWrW7lRU7k6lX+6+FLu1mp4kkUdHRtGjdnhat2/PzhjVMmzaNrSu+oP5lzSh+nvPUacbzFetWLOHlJ/tTs25DGsY3o9ddt9CkSRMKFy7s8VnkL0voQZSWlsYvv/zChg0bWLduHQkJCRSt0oBOXbpm6lepdmN+WvbXIlply1eiUoOr6HJnZ2rVa/jXsMrO3zDGQM069Wl192CeGv0sRw7/VSwj4/mKD5Z/TfLxY6xbsYR1K5bw1r/HExsbS7169WjYsCFRJStTqvJFVLrwYs4pVjxsnzANVpHoOGAGcBlwAPh7Rlm6cLRr1y6+++47fv31V99r8+bNbPjpJ1KSMz+aXOvyA9kSerUGzYhKPkijplcR37wVNevU57PVO6hdvzLGmJOLiY3j/NJls7Vv27whW1tKSgorV65k5cqVmdq79hvGNR3vz9S2YsUK4uLiqFChAiVLliyw96qCVSS6O/CnqtYUkc7AGODveRHw6UpLS6NQoUKZfkDJycls3ryZI0eOcOTIERITEzP9eeTIEf744w/2799PXFwcs2bNynTMzz//nHvvvTegz9+1JftKwrWb3cCA3oGtX2GMyd2Y195h7+87WPvD9yQsW8Ken1ezcWP2J6YBSl9QPltbt27dWLt2LYAvsZcpU4aSJUty3nnnUbJkSd/XXbp0oWLF0FwrKShFot3tEe7X7wEvi4io6knLzG3ZsoXrrruO9PT0k75UlfT0dN59911q1Kjhe+/x48dp0KCBb3/GKy0tjRMnTpCcnExKSgrJycmkp6dz5MgRihUr5nv/tm3bqF+/fkDfoHOKFeeNxRvZdSCRCqWcY2w8cPL+55UqQ+VqNahcvSYX1qrLHzEX+B5nNsbknbLlK3Ft+UrEVr+CiyqcR+Lhg/y25Sd+27KRdWtWkbh7Gzt/3UqFytWzvXfnzr+KsCUnJ7Nt2za2bduW4+e0atWqQCf0nIpENz1ZH1VNFZFDQCngpHfxjhw5EnD1kqSkzOs6iAibN28O6L3g/ID8E3pcXFzA7z2WeJhS50SzevtxGl1YBoCoWhdz+VXXkhp3HhdWr8r5ZStQ+oKKJMWWovXltTO9f9Gq3yyZG5OPMsbVy5UoQo3K5aFlKxat+o3Wl1b2zWFfsnGPb7G69LQ0ylW9mKgiuzny5z6Sk46d6vCULFkyz8/hTMkpLqKdDiK3A21UtYe7fQ/QVFX7+/VZ6/bZ4W5vcfvsz3KsXkAvd7MesDZYJ1JAlOYU/8mFKTvnyBBp5+zl+VZV1TI57QjkCj2QItEZfXaISDRQAufmaCaqOgWYAiAiy1U1PoDPDxt2zpHBzjn8her5BqtI9DwgYyrH7cAXpxo/N8YYE3y5XqG7Y+IZRaKjgGkZRaKB5W5d0anATBH5GfgDJ+kbY4zJRwHNQ1fV+cD8LG3D/b4+Dtxxmp895TT7hwM758hg5xz+QvJ8c70paowxpmCw5XONMSZMeJrQReQOEVknIukiEnJ3jINFRNqIyEYR+VlEHvU6nvwgItNEZK87pTXsiUhlEVksIuvdv9MDvY4pr4lIYRH5XkRWuef8lNcx5RcRiRKRH0Tkv17H4s/rK/S1wG3A1x7HkWf8lk5oC1wC3Ckil3gbVb6YDrTxOoh8lAoMUdVLgGZAvwj4OScD16rqpUBDoI2INPM2pHwzEMi+gIzHPE3oqrpBVXNecCF8+JZOUNUUIGPphLCmql/jzHiKCKr6u6qudL8+gvOPPTSfDw8SdSS6mzHuK+xvyolIJaA98LrXsWTl9RV6JMhp6YSw/oce6USkGtAI+M7jUPKcO/TwI7AXWKSqYX/OwAvAMCDkauHleUIXkc9EZG0Or7C/SjWRR0SKAe8DD6lq9lJUYUZV01S1Ic4T5E1EpJ7HIeUpEbkJ2KuqK7yOJSd5XuBCVa/P688IcYEsnWDCgIjE4CTzt1T1P17Hk59U9aCILMa5bxLON8KbAx1EpB1QGCguIrNU9W6P4wJsyCU/BLJ0gingxFlScyqwQVXHex1PfhCRMiJynvt1EZyaCT95GlQeU9XHVLWSqlbD+bf8Ragkc/B+2uKtIrIDuAL4WEQWeBlPXlDVVCBj6YQNwDuqmr3qRZgRkbeBpUAtEdkhIuFe0aM5cA9wrYj86L7aeR1UHisPLBaR1TgXLotUNaSm8UUae1LUGGPChA25GGNMmLCEbowxYcISujHGhAlL6MYYEyYsoRtjTJiwhB5BRGR6qK0O509ESorIHhGp4XUsBYGIqFvEPdD+3UQkMfeeuR7nexHpeLbHMcFnCd2cFXcJ1ekislpETojIlyfpd42IrBCR4yKyVUT65NDtH8B8Vd2Sp0FHrrnAhUE4zr+AZ0XE8keIsR+IOVtRwHHgZeDjnDqISHWcEoZLcBategaY6H+VJyLnAD1wnrY0eUBVk1R1bxAONR84F2dJaBNCLKFHMBGJE5EX3GGO4yKyTESuytKnvVuc47iIfC0ind1f9asBqOpRVe2jqlNwVpLMSR9gl6oOcJdMfg14Exjq16cdztKr32b5/EtE5GMROeIWzHhbRC5w99USkaMi0tWvfxsRSRGRK9zt6SLyXxF5wj3PRBF5w31UPeM9IiLDRGSLiCSJyBoRudtvfzX3nDuKyCIROeYWsmjt1ydGRF4SkV0ikiwiv4nIs377Y0VkjPvU7DERSRCRG3P5+YiIDBGRze4xd4jIM6fo/6z7s0oSke0iMlZECvvtzzTkIiIjxFkor6vb/6j7vYkVkQfcczggIuP9r8ZVNQ0nqd95qvhN/rOEHtnGAn8H7se5cl4DfCoi5QFEpArwH5wr70uBl9z3nK4rgIVZ2hYA8eIsaAXQAlihfo8uu3F8jbPYUxPgeqAY8H8iUshdS38QztX+hSJSBqewxmhVXer3Wde48V8HdARuAMb47R8FdAf64RQheQb4t4i0zxLzaJzvwaU4j7rPEWd1RYAHgVtx1ve4COf76r/W/xtuHHcB9XD+Q/tIRC49yfcM4GngSTeeujiF2H87Rf+jOD/LOsADbiyPn6I/QDWc9flvwik2cwfOWkOX43yfegAD3HPz9717PiaUqKq9IuSFk+z+635dFEgB7vXbHwVsAUa528/gLDblf4x/4FxJV8vh+C8DX+bQvgkYnqXtavc45d3tD4E3s/QZCXyepa2k+74mfm0fAMtw/uP5FojKcs4HgWJ+bXfjVNsp6r6SgBZZPucFnPF8cJKeAr399ld0265yt18CPsddTiPLsWrgrJ1dJUv7h8Dkk/ysiuEMZfU5xc9TgdtPsb8PTnGVjO1uQKLf9gj33Ev4tb0H7ANi/dq+BF7OcuwO7jlFe/332l5/vfJ8+VwTsmrgVJjxDXGoapqILMW5SgWojXMl6i+vChgUAfZkabsMuPokMzNq4FwlgnMVuRHnKraBOkMC/lbrX5V1wFk0LNY9RhzOMqifioj/wkYxwPasx/H7epf7Z1n3z+nAImCTiCzEGZL4RFXTgcaAAOtFxP94ccAXOZwbOD+DOJz/JALiznh5CKiJ8x9ClPs6lV9V9ZDf9h5gkzrVtfzbymZ+G0k451QYOOuZMyY4LKGbnAR7xbbdQLksbeVw6nDud7f341x9+yuEc9U9lOz8k389oIT7dUVg22nEljHseDPwa5Z9J062rarqJudC7vZK977CjThDO28Cq9xx9kI439PLczhm0mnEelLi1PKcAzyFMwx1EOcqelwub80aj56kLet/DOcDx7P8R2k8Zgk9cm3BGXJp7n6dUdD6CmC22+cnstc/bXIGn7WU7GOwrYHlqpqRPH7AGRLwtxLoBPzi1y8TcdbjnomTuIoAM0XkUs1cLai+iBRV1aPudjOcc9+Ck2yTgaqqerKr5YCoU0v0PeA9EZmOMwxU0z03AS5Q1cUBHm6DG9d1wOYA+jcHdqrqvzIaRKRq4NGftno4Px8TQuymaIRyk9srwBgRaSciddztcsBkt9urQA0RGefOKLkN6J1xiIxjuTNRGgKlgWIi0tDdxu84Fd0ZNXVEpAdO8va/elwA1BGRUn5tk3CuvOeKSFP3xuf1IjJFRM71O/Y+YDjwCHDEfZ+/aGCaiNR1r5ifBV5TZ4bOETeOcSJyv4jUdOPvIyK9Av1+ishgEbnTPb+aODc/DwM7VHUT8BYwXURud88jXkSGut/TbNy4XgSeEZH7RKSGiDQRkb4nCWETzve4i3v8vuTtLJQWwKd5eHxzJrwexLdX/r3wuynqbsfh3Pzbg3M1uAz3Jp9fn5twksVx4H/AfTjJvJxfn+1uW6ZXluNcg3NFl4wzJJLtZh/OlXy/LG0X4Vz1/okzPLERmIgzBn4PcAyo7de/ntuvs/854yT8vTjjvW8C5/i9R3Bmcqx349uHMx7e2t1fzT2n+Cyx+W5KAj3d8zuCk8i/Aq706xuDcxNyK85vB7txZpNcdoqfVyHgUb/3/IYzgyfb57vbz7ixJ+LMTurr/3Mg55uia7N8ZrYb2zhDOe/5bVd046nk9d9pe2V+WYELc1pEZCDO7JPzNMh/eUSkDc5V6SWa/cbmmR5zOlBaVW8KxvEMiMhzODNjAv4NxuQPG0M3pyQi/XBmuuzDGXt+Epge7GQOoKqfisgknELavwT7+CZo9pL7zVbjAUvoJjc1ceael8J5EvRVnCv0PKGqL+XVsU1wqOpzXsdgcmZDLsYYEyZslosxxoQJS+jGGBMmLKEbY0yYsIRujDFhwhK6McaECUvoxhgTJv4fj8voYDGtNL8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Density histogram of the data,\n",
    "# and fitted probability density function\n",
    "\n",
    "y = np.linspace(-1,4,200)\n",
    "lik = np.exp(logPr(y,θhat))\n",
    "\n",
    "fig,ax = plt.subplots(figsize=(6,3))\n",
    "ax.hist(np.log10(mpexpenses), bins=75, density=True, alpha=.4, fc='steelblue', ec='steelblue')\n",
    "ax.plot(y, lik, c='black', linestyle='dashed', linewidth=3)\n",
    "ax.set_xlim([-1,4.5])\n",
    "ax.set_xlabel('log10(expense claim)', fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Appendix: scraping the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "import requests_cache\n",
    "import lxml.html\n",
    "import urllib.parse\n",
    "import re\n",
    "import io\n",
    "import IPython.display\n",
    "import pandas\n",
    "requests_cache.install_cache('mp', backend='memory')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "URL = 'http://www.theipsa.org.uk/mp-costs/your-mp/'\n",
    "resp = requests.get(URL)\n",
    "doc = lxml.html.fromstring(resp.content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mps = doc.xpath('//*[@class=\"mp-letter\"]/ul/li/a')\n",
    "mps = pandas.DataFrame.from_records([(mp.text, mp.get('href')) for mp in mps], columns=['name','url'])\n",
    "print(\"MPs with no url\\n\", mps.loc[mps['url'] == ''])\n",
    "mps = mps.loc[mps['url'] != '']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = []\n",
    "for i, (mp_name, url) in enumerate(zip(mps['name'], mps['url'])):\n",
    "    IPython.display.clear_output(wait=True)\n",
    "    print(f\"{i+1}/{len(mps)}: {mp_name}\")\n",
    "\n",
    "    resp = requests.get(urllib.parse.urljoin(URL, url))\n",
    "    doc = lxml.html.fromstring(resp.content)\n",
    "\n",
    "    s = doc.xpath(\"//script[contains(text(),'downloadclaimscsvformp')]\")[0].text\n",
    "    download = re.search(r'window.open\\(\"(?P<url>/download/downloadclaimscsvformp/(?P<id>.+))\"\\)', s)\n",
    "    resp = requests.get(urllib.parse.urljoin(URL, download['url']))\n",
    "\n",
    "    res.append((download['id'], pandas.read_csv(io.BytesIO(resp.content))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = []\n",
    "for mp_id, mp_df in res:\n",
    "    mp_cleandf = mp_df.loc[~pandas.isnull(mp_df['Claim No.'])].copy(deep=False)\n",
    "    mp_cleandf['mp_id'] = mp_id\n",
    "    df.append(mp_cleandf)\n",
    "df = pandas.concat(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv('mpexpenses.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "df2 = df.loc[(df.Year=='18_19') & (df.Category=='Office Costs'),\n",
    "             ['Date','Claim No.',\"MP's Name\", \"MP's Constituency\", \"Category\", \"Expense Type\", \"Short Description\",\n",
    "              \"Details\",\"Amount Claimed\", \"Amount Paid\",\"Amount Not Paid\",\"Amount Repaid\",\"Status\",'mp_id',\"Reason If Not Paid\"]]\n",
    "\n",
    "df2 = pandas.concat([df2.loc[~pandas.isna(df2['Reason If Not Paid'])], df2.loc[pandas.isna(df2['Reason If Not Paid'])]])\n",
    "df2.to_csv('mpexpenses_2018_office.csv', index=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}