CS146 - Election modeling

7-2_pre-class-work.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling elections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as sts\n",
    "import matplotlib.pyplot as plt\n",
    "import pystan\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data\n",
    "\n",
    "The `electoral_votes` variable is a dictionary containing the number of Electoral College votes for each state. For example\n",
    "```\n",
    "  >>> electoral_votes['Indiana']\n",
    "  11\n",
    "```\n",
    "Data from [Wikipedia: United_States_Electoral_College](https://en.wikipedia.org/wiki/United_States_Electoral_College)\n",
    "\n",
    "The `survey_results` variable is a dictionary mapping from states to an array of survey results for each candidate. Each row in a survey results array represents one survey and each column represents one candidate. There are 4 columns, representing Trump (Republican), Biden (Democrat), Jorgensen (Libertarian), and Hawkins (Green) in that order. In the example below, Trump got 340 votes in the first survey, Biden got 258, Jorgensen got 27, and Hawkins got 13.\n",
    "```\n",
    "  >>> survey_results['Indiana']\n",
    "  array([[340, 258,  27,  13],\n",
    "         [240, 155,   5,   5],\n",
    "         [235, 155,  50,  20],\n",
    "         [308, 266,  49,  35],\n",
    "         [222, 161,  80,  30]])\n",
    "```\n",
    "Data from [Wikipedia: Statewide opinion polling for the 2020 United States presidential election](https://en.wikipedia.org/wiki/Statewide_opinion_polling_for_the_2020_United_States_presidential_election)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Modeling 5 states with 38 electoral college votes\n"
     ]
    }
   ],
   "source": [
    "electoral_votes = {\n",
    "    'Alabama': 9,\n",
    "    'Alaska': 3,\n",
    "    'Arizona': 11,\n",
    "    'Arkansas': 6,\n",
    "    'Colorado': 9,\n",
    "}\n",
    "\n",
    "survey_results = {\n",
    "    'Alabama': np.array([[611, 397, 0, 0], [799, 528, 0, 0], [793, 403, 0, 0], [288, 264, 0, 0], [353, 219, 0, 0], [997, 554, 0, 0], [312, 232, 0, 0], [409, 227, 0, 0], [319, 234, 0, 0]]),\n",
    "    'Alaska': np.array([[348, 320, 0, 0], [298, 253, 0, 0], [283, 277, 0, 0], [269, 198, 0, 0], [227, 177, 0, 0], [442, 389, 0, 0], [519, 486, 0, 0], [325, 318, 0, 0], [84, 74, 0, 0]]),\n",
    "    'Arizona': np.array([[522, 478, 22, 0], [313, 356, 7, 7], [291, 304, 0, 0], [270, 288, 0, 0], [236, 264, 16, 0], [180, 184, 0, 0], [133, 151, 0, 0], [269, 321, 20, 0], [230, 250, 5, 0], [3337, 3621, 0, 0], [360, 392, 0, 0], [235, 235, 0, 0], [364, 396, 8, 0], [383, 409, 9, 9], [221, 216, 0, 0], [113, 128, 0, 0], [284, 278, 0, 0], [168, 212, 0, 0], [258, 270, 0, 0], [260, 266, 0, 0], [359, 402, 9, 0], [185, 202, 17, 0], [261, 320, 26, 0], [519, 584, 0, 0], [328, 342, 0, 0], [487, 520, 0, 0], [252, 312, 0, 0], [752, 768, 0, 0], [414, 441, 0, 0], [212, 230, 0, 0], [357, 398, 0, 8], [309, 378, 23, 0], [3357, 3034, 0, 0], [396, 490, 0, 0], [162, 169, 0, 0], [325, 402, 9, 9], [445, 426, 0, 0], [311, 350, 0, 0], [188, 193, 0, 0], [466, 456, 30, 0], [271, 295, 0, 0], [204, 192, 0, 0], [522, 547, 24, 12], [2547, 2348, 0, 0], [164, 172, 0, 0], [381, 445, 0, 0], [393, 428, 0, 0], [326, 395, 17, 9], [372, 413, 0, 0], [432, 470, 0, 0], [315, 343, 0, 0], [155, 176, 0, 0], [500, 500, 0, 0], [264, 294, 0, 0], [1230, 1088, 0, 0], [270, 282, 0, 0], [137, 159, 0, 0], [258, 237, 0, 0], [337, 372, 17, 9], [266, 312, 0, 0], [616, 670, 0, 0], [88, 90, 0, 0], [421, 461, 0, 0], [148, 145, 0, 0], [368, 353, 0, 0], [180, 188, 0, 0], [388, 426, 0, 0], [258, 300, 0, 0], [230, 235, 0, 0], [258, 312, 0, 0]]),\n",
    "    'Arkansas': np.array([[478, 293, 0, 0], [462, 220, 0, 0], [493, 239, 0, 0], [209, 135, 0, 0], [408, 391, 0, 0]]),\n",
    "    'Colorado': np.array([[408, 510, 0, 0], [1114, 1549, 0, 0], [283, 322, 0, 0], [320, 400, 0, 0], [312, 400, 32, 8], [978, 1359, 0, 0], [262, 325, 0, 0], [252, 306, 0, 0], [246, 307, 0, 0], [246, 306, 0, 0], [240, 312, 0, 0], [935, 1355, 0, 0], [240, 320, 0, 0], [246, 306, 0, 0], [365, 481, 0, 0], [328, 470, 0, 0], [457, 620, 0, 0], [240, 286, 0, 0], [280, 371, 0, 0], [216, 330, 0, 0], [133, 201, 0, 0]]),\n",
    "}\n",
    "\n",
    "states = sorted(survey_results.keys())\n",
    "print('Modeling', len(states), 'states with', sum(electoral_votes[s] for s in states), 'electoral college votes')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generative model\n",
    "\n",
    "1. For each state we generate an $\\vec{\\alpha}$ vector, which defines a Dirichlet distribution over the proportion of votes that go to each of the 4 candidates whenever we do a survey – including the final survey, namely the election itself which we want to predict. The prior over each component of $\\vec{\\alpha}$ is taken as a Cauchy distribution with location 0 and scale 1. Since the components of $\\vec{\\alpha}$ are positive, we actually use the positive half-Cauchy distribution.\n",
    "\n",
    "2. For each survey in a state we generate a probability vector $\\vec{p_i} \\sim \\text{Dirichlet}(\\vec{\\alpha})$ for the probability that a voter selects each of the 4 candidates.\n",
    "\n",
    "3. For each survey, we then generate the number of votes going to each candidate as $\\vec{k_i} \\sim \\text{Multinomial}(\\vec{p_i})$.\n",
    "\n",
    "### Tasks\n",
    "\n",
    "* Use Stan to sample from the posterior distribution over $\\alpha$ and visualize your results. There are 10 states, so you will have 10 posteriors.\n",
    "* The posteriors over $\\alpha$ show a lot of variation between different states. Explain the results you get in terms of the model and the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_af1aafcd9e9c2bdfe0200572b3f4f7fd NOW.\n"
     ]
    }
   ],
   "source": [
    "stan_code = '''\n",
    "data {\n",
    "    int<lower=1> n_surveys;  // number of surveys \n",
    "    int<lower=1> n_candidates;   // number of categories (candidates)\n",
    "    int survey_counts[n_surveys,n_candidates];  // \n",
    "}\n",
    "\n",
    "parameters {\n",
    "    vector[n_candidates] alpha;\n",
    "    simplex[n_candidates] theta;\n",
    "}\n",
    "\n",
    "model {\n",
    "    alpha ~ cauchy(0,1);  // Sample alpha from a Cauchy distribution\n",
    "    for(i in 1:n_surveys) {\n",
    "        theta ~ dirichlet(alpha);  // Likelihood function\n",
    "        survey_counts[i] ~ multinomial(theta);  // Multinomial prior\n",
    "    }\n",
    "}\n",
    "'''\n",
    "\n",
    "stan_model = pystan.StanModel(model_code=stan_code)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "stan_data = {}\n",
    "for state in survey_results:\n",
    "    stan_data[state] = {'survey_counts' : survey_results[state],\n",
    "                       'n_surveys': len(survey_results[state]),\n",
    "                        'n_candidates' : 4}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:n_eff / iter below 0.001 indicates that the effective sample size has likely been overestimated\n",
      "WARNING:pystan:Rhat above 1.1 or below 0.9 indicates that the chains very likely have not mixed\n",
      "WARNING:pystan:22 of 4000 iterations ended with a divergence (0.55 %).\n",
      "WARNING:pystan:Try running with adapt_delta larger than 0.8 to remove the divergences.\n",
      "WARNING:pystan:3978 of 4000 iterations saturated the maximum tree depth of 10 (99.5 %)\n",
      "WARNING:pystan:Run again with max_treedepth larger than 10 to avoid saturation\n",
      "WARNING:pystan:Chain 3: E-BFMI = 0.00107\n",
      "WARNING:pystan:E-BFMI below 0.2 indicates you may need to reparameterize your model\n",
      "WARNING:pystan:n_eff / iter below 0.001 indicates that the effective sample size has likely been overestimated\n",
      "WARNING:pystan:Rhat above 1.1 or below 0.9 indicates that the chains very likely have not mixed\n",
      "WARNING:pystan:2143 of 4000 iterations ended with a divergence (53.6 %).\n",
      "WARNING:pystan:Try running with adapt_delta larger than 0.8 to remove the divergences.\n",
      "WARNING:pystan:1856 of 4000 iterations saturated the maximum tree depth of 10 (46.4 %)\n",
      "WARNING:pystan:Run again with max_treedepth larger than 10 to avoid saturation\n",
      "WARNING:pystan:n_eff / iter below 0.001 indicates that the effective sample size has likely been overestimated\n",
      "WARNING:pystan:Rhat above 1.1 or below 0.9 indicates that the chains very likely have not mixed\n",
      "WARNING:pystan:4000 of 4000 iterations saturated the maximum tree depth of 10 (100 %)\n",
      "WARNING:pystan:Run again with max_treedepth larger than 10 to avoid saturation\n",
      "WARNING:pystan:n_eff / iter below 0.001 indicates that the effective sample size has likely been overestimated\n",
      "WARNING:pystan:Rhat above 1.1 or below 0.9 indicates that the chains very likely have not mixed\n",
      "WARNING:pystan:1337 of 4000 iterations ended with a divergence (33.4 %).\n",
      "WARNING:pystan:Try running with adapt_delta larger than 0.8 to remove the divergences.\n",
      "WARNING:pystan:2657 of 4000 iterations saturated the maximum tree depth of 10 (66.4 %)\n",
      "WARNING:pystan:Run again with max_treedepth larger than 10 to avoid saturation\n",
      "WARNING:pystan:Chain 3: E-BFMI = 0.000257\n",
      "WARNING:pystan:E-BFMI below 0.2 indicates you may need to reparameterize your model\n",
      "WARNING:pystan:n_eff / iter below 0.001 indicates that the effective sample size has likely been overestimated\n",
      "WARNING:pystan:Rhat above 1.1 or below 0.9 indicates that the chains very likely have not mixed\n",
      "WARNING:pystan:36 of 4000 iterations ended with a divergence (0.9 %).\n",
      "WARNING:pystan:Try running with adapt_delta larger than 0.8 to remove the divergences.\n",
      "WARNING:pystan:3964 of 4000 iterations saturated the maximum tree depth of 10 (99.1 %)\n",
      "WARNING:pystan:Run again with max_treedepth larger than 10 to avoid saturation\n"
     ]
    }
   ],
   "source": [
    "## This code block takes several minutes to run ##\n",
    "alabama_results = stan_model.sampling(data = stan_data['Alabama'])\n",
    "alaska_results = stan_model.sampling(data = stan_data['Alaska'])\n",
    "arizona_results = stan_model.sampling(data = stan_data['Arizona'])\n",
    "arkansas_results = stan_model.sampling(data = stan_data['Arkansas'])\n",
    "colorado_results = stan_model.sampling(data = stan_data['Colorado'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation time\n",
    "\n",
    "Use the posterior samples to predict the outcome of the presidential elections.\n",
    "\n",
    "* Predict the probability that each candidate will win each state.\n",
    "   * Use the posterior $\\alpha$ samples to generate posterior predictive samples for $p$ – the proportion of votes each candidate would get in each state in an election.\n",
    "   * Use these $p$ samples to estimate the probability that each candidate will win each state.\n",
    "* Predict the probability that each candidate will win the presidential election.\n",
    "   * Use the posterior predictive probability that each candidate will win each state to generate samples over the total number Electoral College votes each candidate would get in an election.\n",
    "   * Use the total number of votes to generate samples over who would win the election."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAF1CAYAAADiNYyJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoBklEQVR4nO3de7xcVX338e+XE8I1hECOxdy5ViOolcPFApK2ooD4pL6kBUGubTEo1Vp8AFEp9VJpn7agLRoRKUVFtF4wIogURQUJ5MQiCAiNAUkIkRDuiGDg9/yx1kl2JjNn9klmMvGsz/v1Oq8ze6+19157zZq9v7NnzzmOCAEAAACl2azXDQAAAAB6gSAMAACAIhGEAQAAUCSCMAAAAIpEEAYAAECRCMIAAAAoEkF4hGzPtf2hDq1rmu2nbffl6Rts/2Un1p3Xd43tEzq1vtHK9qW2P9rrdpTA9qm2f5XH/Y416p9o+8aN0bYW2x9Rezcm27NsL13PZWfYDttjWpSfbfviZnU35nHF9kdtP2J7eZOyDdn/g2zfs+EtHNE217u9WJftc21/IT9e61w6XN2N0K6OZYQm6669H5zX6iMIV9i+3/aztp+y/bjtH9ueY3t1P0XEnIj4SM11vX64OhHxQERsGxEvdKDt67xAIuKwiPjPDV03ILUPTzWW31zSv0p6Qx73Kzu5/ibb26AQ3a69o1lE/ENENH1TXj2udPONiu2pkk6XNDMidurkuiPiRxHx+51cJ3qnw+fSDQqQdTNCjTassj1pQ9aDegjC63pzRIyTNF3SeZLOlPS5Tm+kUyf70Yr+GZV+T9KWku7sdUNqWu/2Otmg4yuvAU2XtDIiHu51Q4CNxfY2kt4q6QlJx/a4OUUgCLcQEU9ExDxJR0k6wfae0trvFm1PtH1Vvnr8qO0f2d7M9uclTZP0rfxxzRmVq11/YfsBSd9rcQVsV9u32n7C9jdt75C3tc7HakNXnW0fKulsSUfl7f00l6++1SK364O2f2n7YduX2R6fy4bacYLtB/JHkR9o1Te2x+flV+T1fTCvf4vcF3tW6vbnq+wvydNH2L6tcsX9lQ37c6bt2yU90ywI2H6Z7etyf99j+88rZW+y/T+2n7S9xPa5DcsemLf5eC4/sVI8wfa386cBt9jedZj9b7qeVv2Sy060fZPt8/Nyi23/YZ6/JD8nJ1S2canTR2zX5Tb9wPb0Svkf2l6Qx8kC239YKbvB9kfy9p6y/V3bEyvl+1fa/1Pbs2ou+8P8+/E8zl7bpG+2sH2B7WX554I8bw9J91SW/16Trm25ftv/bPsx2/fZPqwyf7ztz9l+yPaDTh+l99l+uaS5kl6b1/V4rj/sGKmst2l7a/T7x2zfJOnXknZpst77bb/f9l15f/7D9pa5bJbtpfk1sFzSf7Tqz4Z1nu30mr3f9rGV+XX29eS83odsn15ZtuVHsHk//7JZH9vex+lWkjGV+m+1fVuLdbU6lrxe0nWSJuV1X9ps+Tb7v0UeNw/kNs21vVW1rxuel/fZvj0/t18eel5y+Rm5j5blfQ/bu7Vozw75eV2Wn+MrG8pPd3q9P2T7pHZ9kct2czoGPJH39cuV5YY7Jl5q+0LXP7b9l+3leTs/tP2KStnhedw+5fRae98w6/kr23fnunfZfk2ef5btX1Tmv6WyzIm2b3Tr1/rOuQ+esn2dpOoxba1z6XB1h9tP26cohc8z8rj7Vp4/yfbX8nNzn+13D7Pv1Yww9Jpu+py38FZJj0v6sKRhb0Ea7vnKJrr1OeQTTseFJ20vtH1QpezcvO4v5GXvsL2H07Hr4bzcGyr1T6o834ttv6PNPm5aIoKf/CPpfkmvbzL/AUmn5seXSvpofvxxpRPB5vnnIEluti5JMySFpMskbSNpq8q8MbnODZIelLRnrvM1SV/IZbMkLW3VXknnDtWtlN8g6S/z45MlLVI6OW8r6euSPt/Qts/mdr1K0nOSXt6iny6T9E1J4/Ky90r6i1x2iaSPVeq+S9J38uPXSHpY0n6S+pRe5PdL2qKyP7dJmippqybb3UbSEkknSRqT1/eIpFdU+mgvpTd4r5T0K0l/msumSXpK0tvyc7WjpFdXntNHJe2b1/tFSVe02Pfh1jNcv5woaVVue5+kjyqNqwslbSHpDXm921ba9JSk1+XyT0i6MZftIOkxScfl9r4tT+9Yed5/IWmP/HzeIOm8XDZZ0kpJh+d+OiRP99dYdoYq47VF/3xY0nxJL5HUL+nHkj5SZ/lm5bnffivpr3K/nSppmda8zq6U9BmlsfESSbdKekdl2RsbtjFLLcZIu/bU7PcHJL0il2/e4hjzM6UxvoOkm7TmeDJLaYz8Y37Ot2rTn0P1/zXXP1jSM5J+v8brYWjfvpT7bi9JK9TkeNKkH27QmuNKsz6+S9JhlelvSDp9PY4ls9RwzGvyXA63/xdImpf7eZykb0n6eLN15+flVkmTcv27Jc3JZYdKWp6f160lfT73x24t2vVtSV+WNEHpGHFwQ3s/nOcfrvSGaUKNvviSpA/k53JLSQfWPCZeqprHtsp5Ylzuzwsk3VYpe0jSQfnxBEmvabGOP1M6j+0jyZJ2kzS9UjYp78dR+fl6ac3X+s2V5/p1SsfHVmO0Zd0a+3mp8msyT28maaGkcySNVTqHLpb0xhb7v3r5ds95i+Wvl/RPSp9Irar2sxrO8zX2o+k5JJe/Xen8NUbpFqTlkrasbOc3kt6Yyy+TdJ/SGNw8P0f3Vdb1Jkm75uf74LyPTcfHpvjT8wZsSj9qHYTnS/pAZXANDfIPKx241jkgNq6r8kLdpcm86gnmvEr5TEnPKx0UZmnDgvD1kt5ZKft9pYPOmEo7plTKb5V0dJP96lMKyTMr894h6Yb8+PWSFlfKbpJ0fH78aeWTeKX8Hq05Udwv6eRhnp+jJP2oYd5nJP1di/oXSDo/P36/pG+0qHeppIsr04dL+nmLuk3XU6NfTpT0v5WyvXKf/15l3kqtHc6vqJRtK+kFpQB1nKRbG7Z/s6QTK8/7Bytl79SaNyNnKr8BqpRfK+mEGssOjZPhgvAvJB1emX6jpPvrLN+sPPfbosr01rnOTkoniudUedOkFE6/X1n2xlZtbRwj7dpTs98/3GZ79ysHrMpY+0V+PEvp9b5lzf6cpXSi3KZS/hVJH6rxehjat5dVyv9J0ufy43O1/kH4TElfzI93UDopvnQ9XjOzVC8Ir7P/SifkZyTtWil7rfLJu3Hd+Xl5e0NfzM2PL1EO0Hl6N7UIwpJeKulFNQk6eZvPau3x/bCk/Wv0xWWSLlLlGJ3nD3tM1AiObU3au33ez/F5+oHcpu3aLHetpPfU3MZtkmZXxlKr1/q0Js/15c3GaLu6NfbzUq0dhPeT9EDDMu+X9B8t1rd6+eGe8xbLTsvjZ+g8cK2kT1TKzx3hfjQ9h7RY/jFJr6ps57pK2ZslPS2pL0+Py9vavsW6rqw7BjaFH26NqGey0rvqRv9P6Srrd/PHAWfVWNeSEZT/Uund18QWdUdiUl5fdd1jlMLEkOo3s3+t9MJpNFHpXXHjuibnx9+TtJXt/fLHMK9WuiIkpXv+Tnf6CPVxp4+rp+a2DRmuf6ZL2q9h+WOVDpTK2/x+/vjqCUlztKbvpiqFilbq7Ptw62nXL1K6IjfkWUmKiMZ51e2u7ouIeFppDE7Sus9ls2212p/pkv6soQ8PVDqBt1u2jmbjbEO/8LG6PRHx6/xwW6V92VzSQ5V9+YzS1dOm2oyRdur0e7vXd2Odxv5ZERG/GWabjfUfi4hnmpXX3Nfh2rK+viDpzba3lfTnSkHtoSb16rxm2mm1//1KQWphZWx8J89vpdW4n6S1+2m453iqpEcj4rEW5SsjYlWT7bTrizOUwv2ttu+0fXKeP+wxsc1+rcXplqLz8q0LTyq9OZDWjJm3KgXpX+aP2de5NSpreay1fbzX3Br3uNKnn9Ux2eq1PknNn+tmhq1bYz8bTVe6Rafax2dr7XPncFo9580cJ+nuiLgtT39R0jFOX9xdS839aHUOGbpF5+58W8XjksY3LNt4bnok1nwZ8dn8e9u8rsNsz3e6PedxpXHSidyyURCE27C9j9LBaJ1vRkfEUxFxekTsovSO6W9t/8lQcYtVtpo/ZGrl8TSlq7aPKF3d2LrSrj6tfVBvt95lSi/o6rpXae3BXscjuU2N63pQkiLiRaWrMm+TdIykqyLiqVxvidJtE9tXfraOiC/V3I8lkn7QsPy2EXFqLr9c6aPQqRExXum2FVeWbXlv3Ai0Ws+w/bKeVo+FHCp2UHoeG5/LkWxridIV4WofbhMR59VYtt0YU5O2Tcvz6qiz/qolSlfRJlb2ZbuIGLpPrtn6hhsj7dTp9zr70Pgar/ZP4/Lt+nOC05drmpXX2dfh2lLHOvsbEQ8qXSl/i9KJ/fMtlu3Ea6bV/j+idLJ+RWVsjI+IkbypG/KQpCmV6amtKiqNyR1sbz/CbbQ7ri6PiL+KiElKV2U/5XSPcrtj4kgcI2m20qd645Wuskp5zETEgoiYrfRG80ql43wzTY+R+cLIZyWdpnQ70fZKtwnVef09pObP9frUHXY/te6YXqL0SUK1j8dFxOE12j1Sx0vaJd/3u1zp9o6Jkg5rUrfdfkgtziH5fuAzld6oTsjPxROqfyxczek7C1+T9M9Kn3BuL+nq9VlXrxCEW7C9ne0jJF2h9FHEHU3qHOH0JQZLelLpY4ehd0y/UpMvy9TwdtszbW+tdOvFV/O7sHslben0BZjNJX1Q6b6fIb+SNMOtv6n+JUnvdfoSwbaS/kHSlxveqbaV2/IVSR+zPS4f3P5W6SrQkMuVPrI7Nj8e8llJc/KVKtveJu/PuJqbv0rSHraPs715/tnH6Us7Uvq45tGI+I3tfZUOFEO+KOn1tv/c9hjbO9p+9Uj2fbj11OyXkTrc6Yt5YyV9RNItEbFE6SCzh+1jchuOUrqN5qoa6xy6WvfGfEVhS6cvdExpu2S6h/RFDT+uvyTpg05fkpyodF9d3T6os/7V8lXG70r6l/x63cz2rrYPzlV+JWlK7r8hw42Rdjak36veZXuK0xdhz1a6n7SVOv3597bH5pPbEZL+K8+vs68fsr2105dsTmrTlmaa9bGUPso/Q+kWoG+ss5RqH0vqWGf/8xvyz0o632u+qDvZ9htHuG7lNp5k++X5uHxOq4p5TF6jFFQn5GPU69ptoF1f2P6zymv0MaWw9oLaHxNHYpzSG8uVShdd/mGoIPfvsbbHR8RvteZ818zFkt5ne+98nN8t7882ud0r8jpPUroi3FZE/FLSoNY81wcqXXxan7ot9zNrPHffKulJpy+xbpWPm3s6XSTrGKcr7Lsq3c/96vyzp9I59IQmi7TbD6n1OWSc0oWwFZLG2D5H0nbr2fSxSllkhaRVTl9wfMPwi2xaCMLr+pbtp5TeBX5A6R1Zq2957i7pv5XunblZ0qci4oZc9nGlE9jjHubbtU18XunenuVKX4p4t5T+ioXS/ZoXK10leEZS9a9IDJ38Vtr+SZP1XpLX/UOlm95/I+mvR9Cuqr/O21+sdKX88rx+5bbekssnKZ0UhuYPKt1k/+9KB/NFSveF1ZKvLL9B0tFKV32Wa80Xi6TUPx/Oz985qlyxiIgHlD6uOV3p46HblL4UOCJt1jNsv6yHyyX9Xd7O3sp/SifS37M9IrdhpVLgOCIiHqnR/iVKVxHOVjpwLZH0f1XjWJA/qvyYpJvyuN6/SbWPKp2Ebpd0h6Sf5Hlt1Vx/o+OVDsR3KY2pr2rNbR7fU/rTZ8ttD/VNyzFSo33r3e8NLlcK8Ivzz3D9064/lyvt9zKlN2lzIuLnuazOvv5A6XV4vaR/jojvjnBfmvWxlMLvdKX76Z9pumSyoa+Z4fb/TKV9m+/00fF/K303YkQi4hpJn5T0/by+m3PRcy0WOU7p6u7Ple4H/ZuamxquL/aRdIvtp5Wu8r8nIu6rcUwcicuUbiF4UOn1NL/Jft2f+3KO0pet1hER/6X0Or5c6ctaV0raISLukvQvSv33K6U3STeNoH3HKN2v+6jScfGy9azbbj8/J2lmPgZdmd+kvFkpmN6ndPX+YqWrsJ10gqRvRsQd+ROA5RGxXOlLbkfkN85V7fZDanEOUbr3+BqlC2y/VMoDdW7rWkceg+9WOr48ptT389ZnXb0y9G1MAJsQpz8XtTQiPtjrtqBzbN+v9EWz/+51W7rN9i+U/oLHqNrXfLX1Z0p/7WZEn6gB2PRwRRgA0FG236r0MXizvxf9O8f2W/LH7BOUrrh+ixAMjA4EYQBAx9i+QelPJb4r36s7GrxD6VaiXyjdG7s+X0YDsAni1ggAAAAUiSvCAAAAKBJBGAAAAEUa06sNT5w4MWbMmNGrzQMAAKAQCxcufCQi1vnvkj0LwjNmzNDg4GCvNg8AAIBC2G76b7m5NQIAAABFIggDAACgSARhAAAAFIkgDAAAgCIRhAEAAFAkgjAAAACKRBAGAABAkQjCAAAAKBJBGAAAAEUiCAMAAKBIBGEAAAAUiSC8EZx/3b29bgIAAAAajOl1A3rhqM/crLuWPaGTD9xF8xevlCR9+R2vlSQdcN71kqSbzvqT1XX332VHfXXhEh2591S995A91pr/3kP20PnX3av5i1dq/112lCR9deESSdKRe09da5tD5ZL03kP20AHnXa8pE7bW/rvsqPmLV65uAwAAALrPEdGTDQ8MDMTg4GBPtn3AedevFVLnL16ppY/9WkfuPVWfvmGRTp212+p5UyZsLUn6nwce0x9MmyBpTWg+/7p79ekbFql/3BZ68tnf6uQDd9ElNy7WzEnj19re0LqHwnI1NA8F8aHpoaANAACAzrC9MCIG1plfahAeCri33PeoJGkzp7Ixm1l/MG2CbrnvUW3mNN0/bgsdufdU/dv3/ldjNrOefyE0bos+zZw0Xrfc96gmb7+lHnz8N5KkcVv06ZnnX9CYzaxVL4a2Gdunp557QeO2SL+HjO3z6sdD6zv5wF0IwgAAAB1GEK6Ycda3e7LddiZvv+XqWzIAAADQGa2CMF+W24Q8+exve90EAACAYhT5Zbn7z3tTr5sAAACAHuOKMAAAAIpEEAYAAECRCMIAAAAoEkEYAAAARSIIAwAAoEgEYQAAABSJIAwAAIAi1QrCtg+1fY/tRbbPalFnlu3bbN9p+wedbSYAAADQWW3/oYbtPkkXSjpE0lJJC2zPi4i7KnW2l/QpSYdGxAO2X9Kl9gIAAAAdUeeK8L6SFkXE4oh4XtIVkmY31DlG0tcj4gFJioiHO9tMAAAAoLPqBOHJkpZUppfmeVV7SJpg+wbbC20f36kGAgAAAN3Q9tYISW4yL5qsZ29JfyJpK0k3254fEfeutSL7FEmnSNK0adNG3loAAACgQ+pcEV4qaWpleoqkZU3qfCcinomIRyT9UNKrGlcUERdFxEBEDPT3969vmwEAAIANVicIL5C0u+2dbY+VdLSkeQ11vinpINtjbG8taT9Jd3e2qQAAAEDntL01IiJW2T5N0rWS+iRdEhF32p6Ty+dGxN22vyPpdkkvSro4In7WzYYDAAAAG8IRjbf7bhwDAwMxODjYk20DAACgHLYXRsRA43z+sxwAAACKRBAGAABAkQjCAAAAKBJBGAAAAEUiCAMAAKBIBGEAAAAUiSAMAACAIhGEAQAAUCSCMAAAAIpEEAYAAECRCMIAAAAoEkEYAAAARSIIAwAAoEgEYQAAABSJIAwAAIAiEYQBAABQJIIwAAAAikQQBgAAQJEIwgAAACgSQRgAAABFIggDAACgSARhAAAAFIkgDAAAgCIRhAEAAFAkgjAAAACKRBAGAABAkQjCAAAAKBJBGAAAAEUiCAMAAKBIBGEAAAAUiSAMAACAIhGEAQAAUCSCMAAAAIpEEAYAAECRCMIAAAAoEkEYAAAARSIIAwAAoEgEYQAAABSpVhC2fajte2wvsn1Wk/JZtp+wfVv+OafzTQUAAAA6Z0y7Crb7JF0o6RBJSyUtsD0vIu5qqPqjiDiiC20EAAAAOq7OFeF9JS2KiMUR8bykKyTN7m6zAAAAgO6qE4QnS1pSmV6a5zV6re2f2r7G9is60joAAACgS9reGiHJTeZFw/RPJE2PiKdtHy7pSkm7r7Mi+xRJp0jStGnTRtZSAAAAoIPqXBFeKmlqZXqKpGXVChHxZEQ8nR9fLWlz2xMbVxQRF0XEQEQM9Pf3b0CzAQAAgA1TJwgvkLS77Z1tj5V0tKR51Qq2d7Lt/HjfvN6VnW4sAAAA0Cltb42IiFW2T5N0raQ+SZdExJ225+TyuZKOlHSq7VWSnpV0dEQ03j4BAAAAbDLcq7w6MDAQg4ODPdk2AAAAymF7YUQMNM7nP8sBAACgSARhAAAAFIkgDAAAgCIRhAEAAFAkgjAAAACKRBAGAABAkQjCAAAAKBJBGAAAAEUiCAMAAKBIBGEAAAAUiSAMAACAIhGEAQAAUCSCMAAAAIpEEAYAAECRCMIAAAAoEkEYAAAARSIIAwAAoEgEYQAAABSJIAwAAIAiEYQBAABQJIIwAAAAikQQBgAAQJEIwgAAACgSQRgAAABFIggDAACgSARhAAAAFIkgDAAAgCIRhAEAAFAkgjAAAACKRBAGAABAkQjCAAAAKBJBGAAAAEUiCAMAAKBIBGEAAAAUiSAMAACAIhGEAQAAUCSCMAAAAIpEEAYAAECRagVh24favsf2IttnDVNvH9sv2D6yc00EAAAAOq9tELbdJ+lCSYdJminpbbZntqj3j5Ku7XQjAQAAgE6rc0V4X0mLImJxRDwv6QpJs5vU+2tJX5P0cAfbBwAAAHRFnSA8WdKSyvTSPG8125MlvUXS3OFWZPsU24O2B1esWDHStgIAAAAdUycIu8m8aJi+QNKZEfHCcCuKiIsiYiAiBvr7+2s2EQAAAOi8MTXqLJU0tTI9RdKyhjoDkq6wLUkTJR1ue1VEXNmJRgIAAACdVicIL5C0u+2dJT0o6WhJx1QrRMTOQ49tXyrpKkIwAAAANmVtg3BErLJ9mtJfg+iTdElE3Gl7Ti4f9r5gAAAAYFNU54qwIuJqSVc3zGsagCPixA1vFgAAANBd/Gc5AAAAFIkgDAAAgCIRhAEAAFAkgjAAAACKRBAGAABAkWr91QgAo8C543vdgtbOfaLXLQAAFIgrwgAAACgSQRgoRd9YyZuteTz9gHXrbLFd+pHW1HXDYaJv7Lp1hqaHfjfTrKxvbPoBAKAHCMJAKabsI203WTr4rPR4xkEpDE8/QBo/NQXS/d+Z6o6fKr3ujLV/D9X70IpUb/xUaey2af5Oe6Vpae11Tj8g3fZw8FnSc0+m3+OnrgnFY7aUDvzb3vQHAKB4BGGgFI8/IL362PR4xkHSbV9ce/6UfdbU3X6adP+P0m9Jeu/PpJOuTtPn75nK3vsz6f1L0rpmHJTKdtor1RtaZsZBqb6UAvAfvT/V2/+dKSDv/840DwCAHnBE9GTDAwMDMTg42JNtA0X6/seHD53V8qHHzeYBAPA7xvbCiBhYZz5BGAAAAKNZqyDMrREAAAAoEkEYAAAARSIIAwAAoEgEYQAAABSJIAwAAIAiEYQBAABQJIIwAAAAikQQBgAAQJEIwgAAACgSQRgAAABFIggDAACgSARhAAAAFIkgDAAAgCIRhAEAAFAkgjAAAACKRBAGAABAkQjCAAAAKBJBGAAAAEUiCAMAAKBIBGEAAAAUiSAMAACAIhGEAQAAUCSCMAAAAIpEEAYAAECRCMIAAAAoEkEYAAAARaoVhG0favse24tsn9WkfLbt223fZnvQ9oGdbyoAAADQOWPaVbDdJ+lCSYdIWippge15EXFXpdr1kuZFRNh+paSvSHpZNxoMAAAAdEKdK8L7SloUEYsj4nlJV0iaXa0QEU9HROTJbSSFAAAAgE1YnSA8WdKSyvTSPG8ttt9i++eSvi3p5GYrsn1KvnVicMWKFevTXgAAAKAj6gRhN5m3zhXfiPhGRLxM0p9K+kizFUXERRExEBED/f39I2ooAAAA0El1gvBSSVMr01MkLWtVOSJ+KGlX2xM3sG0AAABA19QJwgsk7W57Z9tjJR0taV61gu3dbDs/fo2ksZJWdrqxAAAAQKe0/asREbHK9mmSrpXUJ+mSiLjT9pxcPlfSWyUdb/u3kp6VdFTly3MAAADAJse9yqsDAwMxODjYk20DAACgHLYXRsRA43z+sxwAAACKRBAGAABAkQjCAAAAKBJBGAAAAEUiCAMAAKBIBGEAAAAUiSAMAACAIhGEAQAAUCSCMAAAAIpEEAYAAECRCMIAAAAoEkEYAAAARSIIAwAAoEgEYQAAABSJIAwAAIAiEYQBAABQJIIwAAAAikQQBgAAQJEIwgAAACgSQRgAAABFIggDAACgSARhAAAAFIkgDAAAgCIRhAEAAFAkgjAAAACKRBAGAABAkQjCAAAAKBJBGAAAAEUiCAMAAKBIBGEAAAAUiSAMAACAIhGEAQAAUCSCMAAAAIpEEAYAAECRCMIAAAAoEkEYAAAARSIIAwAAoEi1grDtQ23fY3uR7bOalB9r+/b882Pbr+p8UwEAAIDOaRuEbfdJulDSYZJmSnqb7ZkN1e6TdHBEvFLSRyRd1OmGAgAAAJ1U54rwvpIWRcTiiHhe0hWSZlcrRMSPI+KxPDlf0pTONhMAAADorDpBeLKkJZXppXleK38h6ZoNaRQAAADQbWNq1HGTedG0ov1HSkH4wBblp0g6RZKmTZtWs4kAAABA59W5IrxU0tTK9BRJyxor2X6lpIslzY6Ilc1WFBEXRcRARAz09/evT3sBAACAjqgThBdI2t32zrbHSjpa0rxqBdvTJH1d0nERcW/nmwkAAAB0VttbIyJile3TJF0rqU/SJRFxp+05uXyupHMk7SjpU7YlaVVEDHSv2QAAAMCGcUTT2327bmBgIAYHB3uybQAAAJTD9sJmF2n5z3IAAAAoEkEYAAAARSIIAwAAoEgEYQAAABSJIAwAAIAiEYQBAABQJIIwAAAAikQQBgAAQJEIwgAAACgSQRgAAABFIggDAACgSARhAAAAFIkgDAAAgCIRhAEAAFAkgjAAAACKRBAGAABAkQjCAAAAKBJBGAAAAEUiCAMAAKBIBGEAAAAUiSAMAACAIhGEAQAAUCSCMAAAAIpEEAYAAECRCMIAAAAoEkEYAAAARSIIAwAAoEgEYQAAABSJIAwAAIAiEYQBAABQJIIwAAAAikQQBgAAQJEIwgAAACgSQRgAAABFIggDAACgSARhAAAAFIkgDAAAgCIRhAEAAFCkWkHY9qG277G9yPZZTcpfZvtm28/Zfl/nmwkAAAB01ph2FWz3SbpQ0iGSlkpaYHteRNxVqfaopHdL+tNuNBIAAADotDpXhPeVtCgiFkfE85KukDS7WiEiHo6IBZJ+24U2AgAAAB1XJwhPlrSkMr00zxsx26fYHrQ9uGLFivVZBQAAANARdYKwm8yL9dlYRFwUEQMRMdDf378+qwAAAAA6ok4QXippamV6iqRl3WkOAAAAsHHUCcILJO1ue2fbYyUdLWled5sFAAAAdFfbvxoREatsnybpWkl9ki6JiDttz8nlc23vJGlQ0naSXrT9N5JmRsST3Ws6AAAAsP7aBmFJioirJV3dMG9u5fFypVsmAAAAgN8J/Gc5AAAAFIkgDAAAgCIRhAEAAFAkgjAAAACKRBAGAABAkQjCAAAAKBJBGAAAAEUiCAMAAKBIBGEAAAAUiSAMAACAIhGEAQAAUCSCMAAAAIpEEAYAAECRCMIAAAAoEkEYAAAARSIIAwAAoEgEYQAAABSJIAwAAIAiEYQBAABQJIIwAAAAikQQBgAAQJEIwgAAACgSQRgAAABFIggDAACgSARhAAAAFIkgDAAAgCIRhAEAAFAkgjAAAACKRBAGAABAkQjCAAAAKBJBGAAAAEUiCAMAAKBIBGEAAAAUiSAMAACAIhGEAQAAUCSCMAAAAIpEEAYAAECRCMIAAAAoUq0gbPtQ2/fYXmT7rCbltv3JXH677dd0vqkAAABA57QNwrb7JF0o6TBJMyW9zfbMhmqHSdo9/5wi6dMdbicAAADQUWNq1NlX0qKIWCxJtq+QNFvSXZU6syVdFhEhab7t7W2/NCIe6niLAQCoaa//3KvXTWjpjhPu6HUTMErd/bKX97oJLb3853f3uglrqROEJ0taUpleKmm/GnUmS1orCNs+RemKsaZNmzbStgIAMCKETZRoUwubm7I69wi7ybxYjzqKiIsiYiAiBvr7++u0DwAAAOiKOkF4qaSplekpkpatRx0AAABgk1EnCC+QtLvtnW2PlXS0pHkNdeZJOj7/9Yj9JT3B/cEAAADYlLW9RzgiVtk+TdK1kvokXRIRd9qek8vnSrpa0uGSFkn6taSTutdkAAAAYMPV+bKcIuJqpbBbnTe38jgkvauzTQMAAAC6h/8sBwAAgCIRhAEAAFAkgjAAAACKRBAGAABAkQjCAAAAKBJBGAAAAEUiCAMAAKBIBGEAAAAUiSAMAACAIjn9U7gebNheIemXPdl4Z02U9EivG1Eo+r536Pveoe97h77vHfq+d0ZL30+PiP7GmT0LwqOF7cGIGOh1O0pE3/cOfd879H3v0Pe9Q9/3zmjve26NAAAAQJEIwgAAACgSQXjDXdTrBhSMvu8d+r536Pveoe97h77vnVHd99wjDAAAgCJxRRgAAABFIgjXYPsS2w/b/lmLctv+pO1Ftm+3/ZqN3cbRqkbfz7L9hO3b8s85G7uNo5Xtqba/b/tu23fafk+TOoz9LqjZ94z9LrC9pe1bbf809/3fN6nDuO+Cmn3PuO8i2322/8f2VU3KRuW4H9PrBvyOuFTSv0u6rEX5YZJ2zz/7Sfp0/o0Nd6mG73tJ+lFEHLFxmlOUVZJOj4if2B4naaHt6yLirkodxn531Ol7ibHfDc9J+uOIeNr25pJutH1NRMyv1GHcd0edvpcY9930Hkl3S9quSdmoHPdcEa4hIn4o6dFhqsyWdFkk8yVtb/ulG6d1o1uNvkeXRMRDEfGT/PgppYPj5IZqjP0uqNn36II8lp/Ok5vnn8Yv0zDuu6Bm36NLbE+R9CZJF7eoMirHPUG4MyZLWlKZXipOWhvTa/NHadfYfkWvGzMa2Z4h6Q8k3dJQxNjvsmH6XmLsd0X+ePg2SQ9Lui4iGPcbSY2+lxj33XKBpDMkvdiifFSOe4JwZ7jJPN7Fbhw/Ufq3ia+S9G+Sruxtc0Yf29tK+pqkv4mIJxuLmyzC2O+QNn3P2O+SiHghIl4taYqkfW3v2VCFcd8lNfqecd8Fto+Q9HBELByuWpN5v/PjniDcGUslTa1MT5G0rEdtKUpEPDn0UVpEXC1pc9sTe9ysUSPfp/c1SV+MiK83qcLY75J2fc/Y776IeFzSDZIObShi3HdZq75n3HfNAZL+j+37JV0h6Y9tf6Ghzqgc9wThzpgn6fj8jcr9JT0REQ/1ulElsL2TbefH+yqN6ZW9bdXokPv1c5Lujoh/bVGNsd8Fdfqesd8dtvttb58fbyXp9ZJ+3lCNcd8Fdfqecd8dEfH+iJgSETMkHS3pexHx9oZqo3Lc81cjarD9JUmzJE20vVTS3yndxK+ImCvpakmHS1ok6deSTupNS0efGn1/pKRTba+S9Kyko4P/EtMpB0g6TtId+Z49STpb0jSJsd9ldfqesd8dL5X0n7b7lELWVyLiKttzJMZ9l9Xpe8b9RlTCuOc/ywEAAKBI3BoBAACAIhGEAQAAUCSCMAAAAIpEEAYAAECRCMIAAAAoEkEYAAAARSIIAwAAoEgEYQAAABTp/wN1qk2jAXgxSgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "alabama_samples = alabama_results.extract()\n",
    "\n",
    "# Make a new array with same dimensions as alpha\n",
    "p_predicted_AL = np.empty(alabama_samples['alpha'].shape)\n",
    "\n",
    "# Generate one p sample for each alpha sample\n",
    "for i in range(alabama_samples['alpha'].shape[0]):\n",
    "    p_predicted_AL[i] = sts.dirichlet(alabama_samples['alpha'][i]).rvs()\n",
    "    \n",
    "plt.figure(figsize=(12, 6))\n",
    "for i in range(4):\n",
    "    plt.plot(sts.uniform.rvs(loc=i+1-0.1, scale=0.2, size=alabama_samples['theta'].shape[0]), alabama_samples['theta'][:,i], ',', alpha=0.5)\n",
    "plt.title('Distribution over each component of theta for probability of being chosen as candidate in Alabama')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAF1CAYAAADiNYyJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuxElEQVR4nO3de7xcVX338e83JyThkiukSm4EFKoRFSXcChQsUiHFh/qSFhSVS1vEWxX10eAVRQvt04q1XgAVKIqCxWpRo5giKCCBJBZFQGgMSA7hEgPkIibxJL/nj7X3OftMZs7ZJ5lzBrI+79frvM7M3muvvfaaPTPf2bP2HkeEAAAAgNyM6nQDAAAAgE4gCAMAACBLBGEAAABkiSAMAACALBGEAQAAkCWCMAAAALJEEG4j2xfb/nCb6pple73truL+Tbb/th11F/V93/Zp7apvR2X7Ctuf6HQ7cmD7LbYfK/b73WuUP932LSPRthbrH1J7R5Lto213b+Oys22H7dEt5n/A9pealR3J1xXbn7D9W9uPNpm3Pdt/pO37tr+FQ1rnNrcXW7N9nu2vFrf7vZcOVHYE2tW2jNCk7u3ejlz3Q4JwTbYftP172+tsP2X7p7bPtt3bhxFxdkScX7OuVw5UJiIeiojdImJzG9q+1RMkIo6PiH/f3roBafDwVGP5nSR9StKfF/v96nbW32R92xWiB2vvjiwi/iEimn4or76uDOcHFdszJb1H0pyIeG47646ImyPij9tZJzqnze+l23VgpG5GqNGGHtvTtqce9CEID82rI2K8pL0kXSjp/ZK+3O6VtOvNfkdF/+yQniNpnKS7O92Qmra5vU6267WX54D2krQ6Ih7vdEOAkWJ7V0mvlbRG0qkdbs4OgyC8DSJiTURcJ+lkSafZ3l/q/2nR9h62v1scPX7C9s22R9n+iqRZkr5TfF3zvsrRrr+x/ZCkH7U4AvY823fYXmP7v2xPKda11dcZ5VFn28dJ+oCkk4v1/byY3zvUomjXh2z/xvbjtq+0PbGYV7bjNNsPFV9FfrBV39ieWCy/qqjvQ0X9Y4u+2L9SdmpxlP2Pivsn2L6zcsT9JQ3b837bv5D0u2ZBwPYLbC8s+vs+239dmfcXtv/H9lrbK2yf17DsEcU6nyrmn16ZPdn294pvA263/bwBtr9pPa36pZh3uu1bbV9ULLfc9p8U01cUj8lplXVc4fQV28KiTT+2vVdl/p/YXlzsJ4tt/0ll3k22zy/Wt872D23vUZl/aKX9P7d9dM1lf1L8f6rYzw5r0jdjbX/a9sri79PFtP0k3VdZ/kdNurZl/bb/2faTth+wfXxl+kTbX7b9iO2Hnb5K77L9QkkXSzqsqOupovyA+0il3qbtrdHvn7R9q6SnJe3TpN4HbZ9r+55iey63Pa6Yd7Tt7uI58Kiky1v1Z0OdH3B6zj5o+9TK9DrbemZR7yO231NZtuVXsMV2/m2zPrZ9kNNQktGV8q+1fWeLulq9lrxS0kJJ04q6r2i2/CDbP7bYbx4q2nSx7Z2rfd3wuLzX9i+Kx/aa8nEp5r+v6KOVxbaH7ee3aM+U4nFdWTzG326Y/x6n5/sjts8YrC+Kec93eg1YU2zrNZXlBnpNvML251z/te0/bD9arOcntl9UmTev2G/XOT3X3jtAPX9n+96i7D22X15Mn2/715Xpr6ksc7rtW9z6ub530QfrbC+UVH1N6/deOlDZgbbT9llK4fN9xX73nWL6NNvfLB6bB2z//QDbXs0I5XO66WPewmslPSXp45IGHILUjsfL9t8X5WbYPtj2bU7P5Udsf9b2mEHa++wQEfzV+JP0oKRXNpn+kKS3FLevkPSJ4vYFSm8EOxV/R0pys7okzZYUkq6UtKuknSvTRhdlbpL0sKT9izLflPTVYt7RkrpbtVfSeWXZyvybJP1tcftMScuU3px3k/Sfkr7S0LYvFu16qaSNkl7Yop+ulPRfksYXy94v6W+KeZdJ+mSl7Nsk/aC4/XJJj0s6RFKX0pP8QUljK9tzp6SZknZust5dJa2QdIak0UV9v5X0okofvVjpw99LJD0m6S+LebMkrZP0uuKx2l3SAZXH9AlJBxf1XiXp6hbbPlA9A/XL6ZJ6irZ3SfqE0n71OUljJf15Ue9ulTatk/Snxfx/lXRLMW+KpCclvbFo7+uK+7tXHvdfS9qveDxvknRhMW+6pNWS5hX9dGxxf2qNZWersr+26J+PS1ok6Y8kTZX0U0nn11m+2fyi3/4g6e+KfnuLpJXqe559W9IlSvvGH0m6Q9KbK8ve0rCOo9ViHxmsPTX7/SFJLyrm79TiNeaXSvv4FEm3qu/15GilfeQfi8d850H6syz/qaL8UZJ+J+mPazwfym37etF3L5a0Sk1eT5r0w03qe11p1sf3SDq+cv9bkt6zDa8lR6vhNa/JYznQ9n9a0nVFP4+X9B1JFzSru3hc7pA0rSh/r6Szi3nHSXq0eFx3kfSVoj+e36Jd35N0jaTJSq8RRzW09+PF9HlKH5gm1+iLr0v6YPFYjpN0RM3XxCtU87Wt8j4xvujPT0u6szLvEUlHFrcnS3p5izr+Sul97CBJlvR8SXtV5k0rtuPk4vHas+Zz/bbKY/2nSq+PrfbRlmVrbOcVKp6Txf1RkpZK+oikMUrvocslvarF9vcuP9hj3mL5GyT9k9I3Uj3VflbD+/y2PF6q7PuSPizpZ+p7/T9Q0qHFvjJb6XnwrlZtfTb9dbwBz5Y/tQ7CiyR9sLhd3ck/rvTCtdULYmNdlSfqPk2mVd9gLqzMnyNpk9KLQu/O22wdjU+QSn3lG9YNkt5amffHSi865Q4fkmZU5t8h6ZQm29WlFJLnVKa9WdJNxe1XSlpemXerpDcVt7+g4k28Mv8+9b1RPCjpzAEen5Ml3dww7RJJH21R/tOSLipunyvpWy3KXSHpS5X78yT9qkXZpvXU6JfTJf1vZd6Liz5/TmXaavUP51dX5u0mabNSgHqjpDsa1n+bpNMrj/uHKvPeqr4PI+9X8QGoMv96SafVWLbcTwYKwr+WNK9y/1WSHqyzfLP5Rb8tq9zfpSjzXKU3io2qfGhSCqc3Vpa9pVVbG/eRwdpTs98/Psj6HlQRsCr72q+L20crPd/H1ezPo5XeKHetzP+GpA/XeD6U2/aCyvx/kvTl4vZ52vYg/H5JVxW3pyi98e+5Dc+Zo1UvCG+1/UoB7HeSnleZd5ikB5rVXTwub2joi4uL25epCNDF/eerRRCWtKekLWoSdIp1/l799+/HlYLHYH1xpaRLVXmNLqYP+JqoIby2NWnvpGI7Jxb3HyraNGGQ5a6X9M6a67hT0omVfanVc31Wk8f6a8320cHK1tjOK9Q/CB8i6aGGZc6VdHmL+nqXH+gxb7HsrGL/Kd8Hrpf0r5X55w1hO5o+XkWbHlb6oHBLWb5Fne9Si/fNZ9sfQyO233SlT9WN/p/SUdYfOn3VPb9GXSuGMP83Sp8i92hRdiimFfVV6x6tFCZK1TOzn1YKX432UPpU3FjX9OL2jyTtbPsQp6/yD1A6IiSlMX/vKb52ecrp6+qZRdtKA/XPXpIOaVj+VKUXShXrvLH4+mqNpLPV13czlUJFK3W2faB6BusXKR2RK/1ekiKicVp1vb19ERHrlfbBadr6sWy2rlbbs5ekv2rowyOU3sAHW7aOZvvZ9p7w0dueiHi6uLmb0rbsJOmRyrZconT0tKlB9pHB1On3wZ7fjWUa+2dVRGwYYJ2N5Z+MiN81m19zWwdqy7b6qqRX295N0l8rBbVHmpSr85wZTKvtn6oUpJZW9o0fFNNbabXfT1P/fhroMZ4p6YmIeLLF/NUR0dNkPYP1xfuUwv0dtu+2fWYxfcDXxEG2qx+nIUUXFkMX1ip9OJD69pnXKgXp3xTDDrYaGlVo+Vpr+03uGxr3lNK3n9V9stVzfZqaP9bNDFi2xnY22ktpiE61jz+g/u+dA2n1mDfzRkn3RsSdxf2rJL3e6cTdfrbz8Zok6SylD3hrKnXu5zTc89Gizn9Qe/JHxxGEt4Ptg5RejLY6Mzoi1kXEeyJiH0mvlvRu28eUs1tU2Wp6aWbl9iylo7a/VTq6sUulXV3q/6I+WL0rlZ7Q1bp71D+c1fHbok2NdT0sSRGxRemozOskvV7SdyNiXVFuhdKwiUmVv10i4us1t2OFpB83LL9bRLylmP81pa9CZ0bERKVhK64s23Js3BC0qmfAftlGvftCESqmKD2OjY/lUNa1QumIcLUPd42IC2ssO9g+piZtm1VMq6NO/VUrlI6i7VHZlgkRUY6Ta1bfQPvIYOr0e51taHyOV/uncfnB+nOy08k1zebX2daB2lLHVtsbEQ8rHSl/jdIb+1daLNuO50yr7f+t0gfLF1X2jYkRMZQPdaVHJM2o3J/ZqqDSPjnF9qQhrmOw19VHI+LvImKa0lG+zzuNUR7sNXEoXi/pRKVv9SYqHWWVin0mIhZHxIlKHzS/rfQ630zT18jiwMgXJb1daTjRJKVhQnWef4+o+WO9LWUH3E5tvU+vUPomodrH4yNiXo12D9WbJO1TBNFHlY7a7iHp+CZlt+fxelLSCUrnIRxemf4FSb+StG9ETFAK/HVfH5/RCMLbwPYE2ydIulrpq4i7mpQ5wekkBktaq/TVdXn5lsfU5GSZGt5ge47tXZSGXlwb6ZIw90sa53QCzE6SPqQ0Lqj0mKTZbn2m+tclneN0EsFuSp/0rmn4pDqooi3fkPRJ2+OLF7d3Kx0FKn1N6Su7U4vbpS9KOrs4UmXbuxbbM77m6r8raT/bb7S9U/F3kNNJO1IaK/VERGywfbDSC0XpKkmvtP3Xtkfb3t32AUPZ9oHqqdkvQzXP6cS8MZLOl3R7RKyQtECpH15ftOFkpWE0361RZ3m07lXFEYVxTid0zBh0yTSGdIsG3q+/LulDTidJ7qE0rq5uH9Spv1dxlPGHkv6leL6Osv0820cVRR6TNMP9T/YYaB8ZzPb0e9XbnE5MmaL0RnPNAGXr9OfHbI+xfaTSm9t/FNPrbOuHbe/idJLNGYO0pZlmfSylr/LfpzQE6FtbLaXaryV1bLX9xQfyL0q6yH0n6k63/aoh1q2ijWfYfmHxuvyRVgWLffL7SkF1cvEa9aeDrWCwvrD9V5Xn6JNKYW2zBn9NHIrxSh8sVysddPmHckbRv6fanhgRf1Df+10zX5L0XtsHFq/zzy+2Z9ei3auKOs9QOiI8qIj4jaQl6nusj1A6+LQtZVtuZ6HxvfsOSWudTmLduXjd3N/pIFnbOB2xfZ7SeO4Dir/9ld5DT2uyyHY9XhFxk9J79LdsH1Kpc62k9bZfoDROe4dAEB6a79hep/Qp8INKn8haneW5r6T/lrRe6QjI54udS0on0n3I6auUlmfXNvEVpTFGjyqdFPH3UrqKhdJ4zS8pHSX4naTqVSTKN7/Vtn/WpN7Lirp/IukBSRskvWMI7ap6R7H+5UpHyr9W1K+irbcX86cpvSmU05conQjxWaUX82VK48JqKY4s/7mkU5SO+jyqvhOLpNQ/Hy8ev4+o8gk4Ih5S+proPUpDDO5UOilwSAapZ8B+2QZfk/TRYj0HqriUTqTr2Z5QtGG1UuA4ISJ+W6P9K5SOInxA6Q1phaT/qxqvE8VXlZ+UdGuxXx/apNgnlN6EfiHpLqUTMWpdk7Nm/Y3epPSV8j1K+9S16hvm8SOlS589arvsm5b7SI32bXO/N/iaUoBfXvwN1D+D9eejStu9UulD2tkR8atiXp1t/bHS8/AGSf8cET8c4rY062Mphd+9lMYX/q7pksn2PmcG2v73K23bIqevef9b6dyIIYmI70v6jKQbi/puK2ZtbLHIG5WO7v5KaTzou2quaqC+OEjS7bbXKx3lf2dEPFDjNXEorlQaQvCw0vNpUZPterDoy7MlvaFZJRHxH0rP468pnaT2bUlTIuIeSf+i1H+PKX1IunUI7Xu90njdJ5ReF6/cxrKDbeeXJc0pXoO+XXxIebVSMH1A6ej9l5SOwrbTaZL+KyLuKr4BeDQiHlU6UfqE4oNz1XY/XhGxUCnfXGf7QEnvVeq7dUofJIf6wfgZqzzjEsCzhNPloroj4kOdbgvax/aDSiea/Xen2zLcbP9a6QoeO9S2Fkdbf6l0tZshfaMGoDM4IgwAGDG2X6v0NXiz60U/69h+TfF182SlI67fIQQDzx4EYQDAiLB9k9JJN28rxuruCN6sNJTo10pjLXeYsZNADhgaAQAAgCxxRBgAAABZIggDAAAgS6M7teI99tgjZs+e3anVAwAAIBNLly79bURs9QuSHQvCs2fP1pIlSzq1egAAAGTCdtOf3mZoBAAAALJEEAYAAECWCMIAAADIEkEYAAAAWSIIAwAAIEsEYQAAAGSJIAwAAIAsEYQBAACQJYIwAAAAskQQBgAAQJYIwgAAAMgSQXgAFy28v9NNAAAAwDDJPghXw25j8D3n2P22q46LFt7fNEyX0xr/AwAAYORkG4TL8FkNuwMF31ahtrrcRQvv36qOc47dTydfclvv/cMvvKG3zDnH7kcIBgAA6JAsg/DJl9zWL7yefMltvWH15Etu2yqcVstXleWqy5T1laH4ooX3q/vJp3uXOenAmb3LSNKi5av7tQUAAAAjI8sgXCrD6qH77K5D99ldknTNmw/TouWr+4XSct5gykB7zZsP6w3B5xy7n26df0y/9V208H5d8+bDesuWR5vrDsUAAADA9hvd6QZ0wjVvPkyHX3iDZkzeRRctvF+Llq/WofvsrpMvua03oF52y3Jdu3SFZkzeRVIKsZfdsry3jkXLV/ce6T10n917w3M1RH/hpmW9Ze9ZuUZnHrFP7zpe/NEfaGPPFt3/yXm96wcAAMDIyfKI8EUL79dJB85U95NP69qlK3qDbPeTT+vwC2/Q7Q88obs+dpweWbNBh+6zu7qffFqLlq/uXb4aiK9dukKX3bK8NxR3P/l0b4h+2azJvfXPmTZxq8D7lqOf3zt04t9+9L/9xhIDAABgeGUZhP/tR/+rRctX69b5x2jt7//QG4LLIQyH7D1FFy28XwfNnqJrl67Q2t//od/yZx6xj7qffLp3+pxpE3XSgTN17dIVOunAmZozbWLvkeQZk3fROcfup+4nn9ah++zeOzRizrSJklKQLpVHowEAADD8shwaIan36K8k3bNyTe/0GZN3KY7qpmlzpk1Ut9LR3nUbN+uyW5Zrws476aQDZ/Ye4b126QpduzQF42qwnTF5l94xwKX9PrhAU8eP1a3zj+ldfxmaAQAAMHIcER1Z8dy5c2PJkiUdWfeLP/oDrdu4WWO6rE2b0/aPsrQl1DttTJclqfd2z5bQlpCmTxqnVes2aur4sTrpwJn6wk3LNHb0KK3buLlfPaMs7TqmS5I0Yeed9PBTG3TI3lN0+wNPaPzYrt71v2zWZEnS4gef0Dv+bF9OmAMAAGgz20sjYu5W03MMwrPnf2+76yjDbCtlIB5qnXd97LjtbBkAAACqWgXhLMcIt8NAIVgaegiW0pFjAAAAjAyCcINRHtr0qnI4xUDLlfeb1cdYYQAAgJGTdRCuhtHpk8ZJkka3SLzl+OGBVMcbV5erKscN7zlx3Fbt4FrCAAAAIyfLIDymyxrT5X4h9eGnNkjqOzmuDKfjx3b1lilPbJP6h93G7FzWO8r9lx/lviEV5fqkvnBcvVYxAAAAhleWQfhlsyb3uzLEmC7rncfs23tUeNPm6L3yw7qNmzXK6fbtDzyhQ/aeoumTxmn0KGv82K5+Qbd6+53H7CtJ2tizpXdaeRS4XG81QI/pMtcRBgAAGEFZBmGpbyjEO4/Zt3dIw8NPbdD0SeN6j+TuOTHdLgPsKKfrD5904Ez1bInek9vKcmcesY/Gj+3SmC73/vrcy2ZN1juP2VdjuqyTDpyp8WO7NHX82N4wPcrpWsVjR2f7UAAAAHRElpdPO/zCG3qHJjReBm36pHF6ZM0GjR7Vdz3hTZuj93Jo1WsAl9cPLo/wVq8nLG19beLyiHHjNYfLOg/ZewpHhQEAANqMy6dVVMfn/m7T5t6juKUtIfVsSeG3PFpcjvstQ2zPlugXiDdtjq0uqVaeeFce7Z2w805NL7tWDr/ofvLp9m0kAAAABpRlEC7H/Ep9Abd6xFZKJ7CVQxdK5fzxY7t6l5s6fuxWR5TLdUwdP1ZSCtuStGrdxn7tqJ6st+fEcTrpwJnbvW0AAACoJ8sgvOfEcdpz4rjeYLuxZ4sO2XuKNvZs0cNPbdD4sV3a2LNFY0eP0kGzp/SOGa6e+CalsLv293/Q9EnjdMjeqdzDT23oPVFu1bqNGtPl3kA9dvQojelKJ9mVgbn8DwAAgJE1utMN6IQZk3fRPSvXaM60ib1jci9aeL+uefNhumjh/Trn2P16y1anL1q+WvesXKMzj9hHUrrcWePyJ19ymxYtX62DZk+RpN75h194gyTppANnatHy1f2uGXzt0hWS1G+9AAAAGF5ZniwnSSdfclvvFSDKAHryJbf1ht5SGVobw3FjaD38wht06/xjtppX1gkAAIDOaHWyXLZBeCDNgu72lAMAAEDnbNdVI2wfZ/s+28tsz28y/2jba2zfWfx9pB2N7pS64ZYQDAAA8Ow16Bhh212SPifpWEndkhbbvi4i7mkoenNEnDAMbQQAAADars4R4YMlLYuI5RGxSdLVkk4c3mYBAAAAw6tOEJ4uaUXlfncxrdFhtn9u+/u2X9SW1gEAAADDpM7l09xkWuMZdj+TtFdErLc9T9K3Je27VUX2WZLOkqRZs2YNraUAAABAG9U5ItwtqfqTZzMkrawWiIi1EbG+uL1A0k6292isKCIujYi5ETF36tSp29FsAAAAYPvUCcKLJe1re2/bYySdIum6agHbz7Xt4vbBRb2r291YAAAAoF0GDcIR0SPp7ZKul3SvpG9ExN22z7Z9dlHsJEm/tP1zSZ+RdEp06gLFAFq78YL0//J5fbdLl8/r/79avpWL9m9e/40XpL/G9TTeHqx+AACGUa3rCEfEgojYLyKeFxGfLKZdHBEXF7c/GxEvioiXRsShEfHT4Ww0gDZ4xbl9QfSCmdLsI9P92UemaZfPkx68Od2uBujq/XN+2T/8lv/vvKr/usp5d17VPwDfeVX/4A0AwAjil+WAXFy0v7T+MWnGQdKjd0kb10pHze8LrWtWSF1jpNHj0ryxE6SeDWnelh5p1Oi0bPfivjo3b+orN3pcmtazQTri3dItn0r3P7wqBe1N66Uxu0nnrugLv089JK19WProkyPTBwCALLX6Zbk6V40AsCPYsCb9716cAqwk/fjCFH7L+1VlCN7tOWnZTetTgN7Sk6bHFsmjUmguy60prrS46PN9Ifm8iWkdo0ansudNTMvFljR9QrOrMQIAMPxqDY0AsAMog2019HaNSUd5S5s39QXb8v76x/qW3bg2BdoyxI4a3VfP+sf6ltu0vq9815gUnpuF7WbTAAAYIQRhIBebN6W/rjHpvkel+7+5Nd32qL7pZfmxE/qWK0NvGV43b+o7OlyWkVL9saX/esv7Zd2jKl9GTeKa4gCAziAIA7k4ar40cWYawjB2gjTrsBRaJ85MQXXC9HR7wvQ0vWuMNG5iWm7shDTud6/D059Hpf9SWuao+X0Be0tPmlausyxT1j1xZiozdkK6/dRDI98XAACIk+WAfFw+L40PLo/ijhqdwu1P/indLsf4luN3J85Mwx2qwxeq8zasSUMmWs0vleOLy7HA5TjisRPS0Im9DpfOWDC82w4AyBonywG5K4dAlMMUtvSkk+XKIRJlQC3nr1mRwuqWnr5p1XlVZWgu622cX12mbMPGtek2IRgA0CEMjQByMXZC/7G7jeG2qhzLu3Ft3/xy6ENVOW3NihSmR4/rPx7YTV5iqvXFFq4jDADoGIIwkIueDf2HLAykWTge7AoP1UuplXVU62kM0mV9jBEGAHQIQRjIRfU6v+1SDcfNwnOrslXl9Y0BABhhBGEgF2sfHrxMs6EM7SjbqHp0uHoUGQCAEUQQBnIx67A0NKIaQsdO6F+mOn63VbnydvVawB619bCLsRP6/qrKk/PKIF1eYg0AgBFGEAZy8+FVfdfwPXeFdN6adHuvw/uuNTx6XN+1fvc6PF1PeK/DU9meDSkkH/Hu9P+o+Slkb6jUM3GmdOhb09/GtX0B/Kj5KUCXP6183hrpFed2ukcAAJkiCAO5eOoh6ZxfSudPTffLsbnlVRvOWJBC6QGnpuA7YXrfsIUNa6TZR0o3XpB+kvnDq9L0GQdJiz6f6n7ui9Oys4/sv16PSr8et9tz0v3dnpOWnzQrrfvGC4Z3uwEAaIEf1AByUQbeR+9KofWh21LYnTQr/dDGEe9O8++8Ko0nLn8A45xfprB651Xp/vrHUpjdsCbV0704BeJmV38455fSRfungHzLp9K0GQf1L8N1hAEAw6zVD2pwRBjIxRkLUmg99K3pqO2sw9L02UemEPyKc6UHb07hdcL0FF4POFW6oDL294BTUwg+4NQUgs9YkILtGQtSoC6vSjFpVvp/+bxU9sGb0/0j3p3WN/vIFMgbjx4DADCC+GU5IBc3XpCGJJRDEcojsTde0DdOt5xWHgWWUuCVUqAt/7/iXOlG9V+mWX2ly2/uG05x+bwUgA99a9s2DQCAbcERYSAXZTgtj/zeeEEatlAqg285hKIs33jU9hXnbl222Vjf6v1yfHF5uzxCzIlyAIAOYowwkJvGwFqG0fOn9h21bSxfht/G4FqdV9Z1+bz6436b1QkAQJu1GiPM0AggN9XgWQ3F5cly1XmvOLf/keRmyzXOawzBhF0AwDMUR4QBtDZYiCXkAgCeBbhqBIChGyzkEoIBAM9iBGEAAABkiSAMAACALBGEAQAAkCWCMAAAALJEEAYAAECWCMIAAADIEkEYAAAAWSIIAwAAIEsEYQAAAGSJIAwAAIAsEYQBAACQJYIwAAAAskQQBgAAQJYIwgAAAMgSQRgAAABZIggDAAAgSwRhAAAAZIkgDAAAgCwRhAEAAJAlgjAAAACyRBAGAABAlgjCAAAAyBJBGAAAAFmqFYRtH2f7PtvLbM8foNxBtjfbPql9TQQAAADab9AgbLtL0uckHS9pjqTX2Z7Totw/Srq+3Y0EAAAA2q3OEeGDJS2LiOURsUnS1ZJObFLuHZK+KenxNrYPAAAAGBZ1gvB0SSsq97uLab1sT5f0GkkXt69pAAAAwPCpE4TdZFo03P+0pPdHxOYBK7LPsr3E9pJVq1bVbCIAAADQfqNrlOmWNLNyf4aklQ1l5kq62rYk7SFpnu2eiPh2tVBEXCrpUkmaO3duY5gGAAAARkydILxY0r6295b0sKRTJL2+WiAi9i5v275C0ncbQzAAAADwTDJoEI6IHttvV7oaRJekyyLibttnF/MZFwwAAIBnnTpHhBURCyQtaJjWNABHxOnb3ywAAABgePHLcgAAAMgSQRgAAABZIggDAAAgSwRhAAAAZIkgDAAAgCwRhAEAAJAlgjAAAACyRBAGAABAlgjCAAAAyBJBGAAAAFkiCAMAACBLBGEAAABkiSAMAACALBGEAQAAkCWCMAAAALJEEAYAAECWCMIAAADIEkEYAAAAWSIIAwAAIEsEYQAAAGSJIAwAAIAsEYQBAACQJYIwAAAAskQQBgAAQJYIwgAAAMgSQRgAAABZIggDAAAgSwRhAAAAZIkgDAAAgCwRhAEAAJAlgjAAAACyRBAGAABAlgjCAAAAyBJBGAAAAFkiCAMAACBLBGEAAABkiSAMAACALBGEAQAAkCWCMAAAALJEEAYAAECWCMIAAADIEkEYAAAAWSIIAwAAIEsEYQAAAGSJIAwAAIAs1QrCto+zfZ/tZbbnN5l/ou1f2L7T9hLbR7S/qQAAAED7jB6sgO0uSZ+TdKykbkmLbV8XEfdUit0g6bqICNsvkfQNSS8YjgYDAAAA7VDniPDBkpZFxPKI2CTpakknVgtExPqIiOLurpJCAAAAwDNYnSA8XdKKyv3uYlo/tl9j+1eSvifpzPY0DwAAABgedYKwm0zb6ohvRHwrIl4g6S8lnd+0IvusYgzxklWrVg2poQAAAEA71QnC3ZJmVu7PkLSyVeGI+Imk59neo8m8SyNibkTMnTp16pAbCwAAALRLnSC8WNK+tve2PUbSKZKuqxaw/XzbLm6/XNIYSavb3VgAAACgXQa9akRE9Nh+u6TrJXVJuiwi7rZ9djH/YkmvlfQm23+Q9HtJJ1dOngMAAACecdypvDp37txYsmRJR9YNAACAfNheGhFzG6fzy3IAAADIEkEYAAAAWSIIAwAAIEsEYQAAAGSJIAwAAIAsEYQBAACQJYIwAAAAskQQBgAAQJYIwgAAAMgSQRgAAABZIggDAAAgSwRhAAAAZIkgDAAAgCwRhAEAAJAlgjAAAACyRBAGAABAlgjCAAAAyBJBGAAAAFkiCAMAACBLBGEAAABkiSAMAACALBGEAQAAkCWCMAAAALJEEAYAAECWCMIAAADIEkEYAAAAWSIIAwAAIEsEYQAAAGSJIAwAAIAsEYQBAACQJYIwAAAAskQQBgAAQJYIwgAAAMgSQRgAAABZIggDAAAgSwRhAAAAZIkgDAAAgCwRhAEAAJAlgjAAAACyRBAGAABAlgjCAAAAyBJBGAAAAFkiCAMAACBLBGEAAABkiSAMAACALNUKwraPs32f7WW25zeZf6rtXxR/P7X90vY3FQAAAGifQYOw7S5Jn5N0vKQ5kl5ne05DsQckHRURL5F0vqRL291QAAAAoJ3qHBE+WNKyiFgeEZskXS3pxGqBiPhpRDxZ3F0kaUZ7mwkAAAC0V50gPF3Sisr97mJaK38j6fvb0ygAAABguI2uUcZNpkXTgvYrlILwES3mnyXpLEmaNWtWzSYCAAAA7VfniHC3pJmV+zMkrWwsZPslkr4k6cSIWN2sooi4NCLmRsTcqVOnbkt7AQAAgLaoE4QXS9rX9t62x0g6RdJ11QK2Z0n6T0lvjIj7299MAAAAoL0GHRoRET223y7pekldki6LiLttn13Mv1jSRyTtLunztiWpJyLmDl+zAQAAgO3jiKbDfYfd3LlzY8mSJR1ZNwAAAPJhe2mzg7T8shwAAACyRBAGAABAlgjCAAAAyBJBGAAAAFkiCAMAACBLBGEAAABkiSAMAACALBGEAQAAkCWCMAAAALJEEAYAAECWCMIAAADIEkEYAAAAWSIIAwAAIEsEYQAAAGSJIAwAAIAsEYQBAACQJYIwAAAAskQQBgAAQJYIwgAAAMgSQRgAAABZIggDAAAgSwRhAAAAZIkgDAAAgCwRhAEAAJAlgjAAAACyRBAGAABAlgjCAAAAyBJBGAAAAFkiCAMAACBLBGEAAABkiSAMAACALBGEAQAAkCWCMAAAALJEEAYAAECWCMIAAADIEkEYAAAAWSIIAwAAIEsEYQAAAGSJIAwAAIAsEYQBAACQJYIwAAAAskQQBgAAQJYIwgAAAMgSQRgAAABZIggDAAAgS7WCsO3jbN9ne5nt+U3mv8D2bbY32n5v+5sJAAAAtNfowQrY7pL0OUnHSuqWtNj2dRFxT6XYE5L+XtJfDkcjAQAAgHarc0T4YEnLImJ5RGySdLWkE6sFIuLxiFgs6Q/D0EYAAACg7eoE4emSVlTudxfTAAAAgGetOkHYTabFtqzM9lm2l9hesmrVqm2pAgAAAGiLOkG4W9LMyv0ZklZuy8oi4tKImBsRc6dOnbotVQAAAABtUScIL5a0r+29bY+RdIqk64a3WQAAAMDwGvSqERHRY/vtkq6X1CXpsoi42/bZxfyLbT9X0hJJEyRtsf0uSXMiYu3wNR0AAADYdoMGYUmKiAWSFjRMu7hy+1GlIRMAAADAswK/LAcAAIAsEYQBAACQJYIwAAAAskQQBgAAQJYIwgAAAMgSQRgAAABZIggDAAAgSwRhAAAAZIkgDAAAgCwRhAEAAJAlgjAAAACyRBAGAABAlgjCAAAAyBJBGAAAAFkiCAMAACBLBGEAAABkiSAMAACALBGEAQAAkCWCMAAAALJEEAYAAECWCMIAAADIEkEYAAAAWSIIAwAAIEsEYQAAAGSJIAwAAIAsEYQBAACQJYIwAAAAskQQBgAAQJYIwgAAAMgSQRgAAABZIggDAAAgSwRhAAAAZIkgDAAAgCwRhAEAAJAlgjAAAACyRBAGAABAlgjCAAAAyBJBGAAAAFkiCAMAACBLBGEAAABkiSAMAACALBGEAQAAkCWCMAAAALJEEAYAAECWCMIAAADIEkEYAAAAWaoVhG0fZ/s+28tsz28y37Y/U8z/he2Xt7+pAAAAQPsMGoRtd0n6nKTjJc2R9DrbcxqKHS9p3+LvLElfaHM7AQAAgLYaXaPMwZKWRcRySbJ9taQTJd1TKXOipCsjIiQtsj3J9p4R8UjbWwwAQE0v/vcXd7oJLd112l2dbgJ2UPe+4IWdbkJLL/zVvZ1uQj91gvB0SSsq97slHVKjzHRJ/YKw7bOUjhhr1qxZQ20rAABDQthEjp5pYfOZrM4YYTeZFttQRhFxaUTMjYi5U6dOrdM+AAAAYFjUCcLdkmZW7s+QtHIbygAAAADPGHWC8GJJ+9re2/YYSadIuq6hzHWS3lRcPeJQSWsYHwwAAIBnskHHCEdEj+23S7peUpekyyLibttnF/MvlrRA0jxJyyQ9LemM4WsyAAAAsP3qnCyniFigFHar0y6u3A5Jb2tv0wAAAIDhwy/LAQAAIEsEYQAAAGSJIAwAAIAsEYQBAACQJYIwAAAAskQQBgAAQJYIwgAAAMgSQRgAAABZIggDAAAgS04/CteBFdurJP2mIytvrz0k/bbTjcgUfd859H3n0PedQ993Dn3fOTtK3+8VEVMbJ3YsCO8obC+JiLmdbkeO6PvOoe87h77vHPq+c+j7ztnR+56hEQAAAMgSQRgAAABZIghvv0s73YCM0fedQ993Dn3fOfR959D3nbND9z1jhAEAAJAljggDAAAgSwThGmxfZvtx279sMd+2P2N7me1f2H75SLdxR1Wj74+2vcb2ncXfR0a6jTsq2zNt32j7Xtt3235nkzLs+8OgZt+z7w8D2+Ns32H750Xff6xJGfb7YVCz79nvh5HtLtv/Y/u7TebtkPv96E434FniCkmflXRli/nHS9q3+DtE0heK/9h+V2jgvpekmyPihJFpTlZ6JL0nIn5me7ykpbYXRsQ9lTLs+8OjTt9L7PvDYaOkP4uI9bZ3knSL7e9HxKJKGfb74VGn7yX2++H0Tkn3SprQZN4Oud9zRLiGiPiJpCcGKHKipCsjWSRpku09R6Z1O7YafY9hEhGPRMTPitvrlF4cpzcUY98fBjX7HsOg2JfXF3d3Kv4aT6Zhvx8GNfsew8T2DEl/IelLLYrskPs9Qbg9pktaUbnfLd60RtJhxVdp37f9ok43Zkdke7akl0m6vWEW+/4wG6DvJfb9YVF8PXynpMclLYwI9vsRUqPvJfb74fJpSe+TtKXF/B1yvycIt4ebTONT7Mj4mdLPJr5U0r9J+nZnm7Pjsb2bpG9KeldErG2c3WQR9v02GaTv2feHSURsjogDJM2QdLDt/RuKsN8Pkxp9z34/DGyfIOnxiFg6ULEm0571+z1BuD26Jc2s3J8haWWH2pKViFhbfpUWEQsk7WR7jw43a4dRjNP7pqSrIuI/mxRh3x8mg/U9+/7wi4inJN0k6biGWez3w6xV37PfD5vDJf0f2w9KulrSn9n+akOZHXK/Jwi3x3WS3lScUXmopDUR8UinG5UD28+17eL2wUr79OrOtmrHUPTrlyXdGxGfalGMfX8Y1Ol79v3hYXuq7UnF7Z0lvVLSrxqKsd8Pgzp9z34/PCLi3IiYERGzJZ0i6UcR8YaGYjvkfs9VI2qw/XVJR0vaw3a3pI8qDeJXRFwsaYGkeZKWSXpa0hmdaemOp0bfnyTpLbZ7JP1e0inBr8S0y+GS3ijprmLMniR9QNIsiX1/mNXpe/b94bGnpH+33aUUsr4REd+1fbbEfj/M6vQ9+/0IymG/55flAAAAkCWGRgAAACBLBGEAAABkiSAMAACALBGEAQAAkCWCMAAAALJEEAYAAECWCMIAAADIEkEYAAAAWfr/LNUJQw/DGccAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "alaska_samples = alaska_results.extract()\n",
    "\n",
    "# Make a new array with same dimensions as alpha\n",
    "p_predicted_AK = np.empty(alaska_samples['alpha'].shape)\n",
    "\n",
    "# Generate one p sample for each alpha sample\n",
    "for i in range(alaska_samples['alpha'].shape[0]):\n",
    "    p_predicted_AK[i] = sts.dirichlet(alaska_samples['alpha'][i]).rvs()\n",
    "    \n",
    "plt.figure(figsize=(12, 6))\n",
    "for i in range(4):\n",
    "    plt.plot(sts.uniform.rvs(loc=i+1-0.1, scale=0.2, size=alaska_samples['theta'].shape[0]), alaska_samples['theta'][:,i], ',', alpha=0.5)\n",
    "plt.title('Distribution over each component of theta for probability of being chosen as candidate in Alaska')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAF1CAYAAADiNYyJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmhklEQVR4nO3debgkVX3/8feXARQRWccozAyL4AK4jwgRlRiJiBj00ShuCMYgRo3xhz9EYwxxiSS/xCVRg4iEoCIaF4KKQaOioqAMBhdAzYjojAMyLA7gAg58f3+cc2dqerrvrXun79xxzvv1PP3c6qpTVadOn6r6dHV138hMJEmSpNZsMdcVkCRJkuaCQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg/AGiohTI+Kvx7SsRRFxW0TMq88vjIgXj2PZdXmfjYgXjmt5m6uIODMi3jzX9WhBRLw0In5e+/3OPcofExEXbYy6jVj/tOq7MUXEIRGxfIbz7hERGRFbjpj+uog4fVjZjXlciYg3R8QNEXHdkGkbsv2PjYgfbHgNp7XOGddX64uIkyPig3V4nXPpZGU3Qr3GlhF6rm/Nvqp+hh70VETENcDvAauBO4ErgbOA0zLzLoDMPH4ay3pxZv73qDKZ+VPgnhtW6zXrOxnYOzOf31n+k8exbAlKIAJ+DGyVmatnMP9WwNuAAzPz2+Ne/pDlHUPZBw+e4fyT1ndzlpl/N8m0NceVDW3jyUTEQuAEYPfMvH6cy87MrwIPGOcyNXfGfC49E1iema+fYV16ZYQedXg+sCgzV0yxvpH7qobzivDUnpqZ2wG7A6cArwHeP+6VjLoSo8L22Sz9HnB34Iq5rkhPM65vFBt0vHUfYHfgxnGHYGlTFhHbAs8AVgHPm6Js68eImclMHyMewDXAEwfGHQDcBexfn58JvLkO7wJ8GvgFcBPwVcqbjQ/UeX4N3AacCOwBJPCnwE+Br3TGbVmXdyHwVuCblJ3gP4Gd6rRDKO9S16svcBhwB/Dbur5vd5b34jq8BfB64CfA9ZQr3dvXaRP1eGGt2w3AX03STtvX+VfW5b2+Lv9utS3275SdX9vh3vX5EcDltdzXgYcMbM9rgO8At0+0y8C6Hwh8vrb3D4BndaY9Bfgf4BZgGXDywLwH13X+ok4/pvOavhv4DHAr8A3gfpNs/6jlDG2XOu0Y4GvA2+t8VwO/X8cvq6/JCzvrOBM4tW7rrcCXKVfGJqb/PnAppZ9cCvx+Z9qFwJvq+m4FPgfs0pl+YKf+3wYO6TNv7RtJ6WO3AQcNaZu7Ae8AVtTHO+q4+wO/7Mz/xSHzrrf82j4XAf8I3Ey5Yvzkgb74fuBa4GfAm4F5wIOA31A+2bkN+EWfPtJZ7tD69mj3t9S2+zXlE5phx5jXUj5tuhn4N+Du3X2csg9cRzmODG3PgfKvo+yz1wDP67M/sHafP64u91rghM70k4EPDpTtHqdePKyNgUcBP6ez71JO6pdP81jyxNqGd9Vlnzlk3qm2/26UfvPTWqdTgW2GHU/rvK+mHHtWAR+ZeF3q9BNrG62o257DXt9adqf6uq6or/G5A/U9gbK/XwscO1Vb1Gl7U44Bq+q2fqTnMfFMpnds+w9K31tFOUft15l2OKXf3krZ1149yXL+DLiqlr0SeEQdfxLwo874p3fmOYbJ9/U9axvcWrf3XYzuoyPLTradlP3ht5Tz6W3Ap+r4XYGP19fmx8BfTLLtZ7I2I0z6mo+Y/2jK/vpK4HsD004GPgZ8kLJfv5h199V3sfb4eRvl0+2T67QHUfbdX1De3P9x334CvLPW6RbgMuCxk23Dpv6Y8wpsyg+GBOE6/qfASzsdZqKTv5VycN2qPh4LxLBldXbUs4BtgW2G7LwXUg4w+9cyH+908EMYEYTr8JqdoTP9QtYG4RcBS4G9KB8hfQL4wEDd3lfr9VBKEH3QiHY6ixLSt6vz/hD40zrtDOAtnbIvA/6rDj+CcjB4NCWsvLBuw90623M5sJB6whpY77Z1ZzyWcpvPIygnhf06bfRgyon0IZST39PqtEWUHfw59bXaGXhY5zW9ifKmZ0vgQ8A5I7Z9suVM1i7HUA5Kx9ZtfzOlX72bcsL+o7rce3bqdCvwuDr9ncBFddpOlBPFC2p9n1Of79x53X9ECXPb1Oen1Gm7ATdSTmpbAIfW5/N7zLsHnf46on3eCFwC3JvyJujrwJv6zD9sem2331JOrPOAl1ICxsR+di7wXkrfuDflTeRLOvNeNLCOQxjRR6aqT892/ymwX52+1YhjzPcofXwnSmjunjRXA39fX/NtpmjPifJvq+UfTwnvD+ixP0xs24dr2z2YcpJf73gypB0uZO1xZVgbX8m6AeaTdEL2NI4lhzBwzBvyWk62/e8AzqvtvB3wKeCtw5ZdX5dvUgLPTpQQd3yddhglNO0H3IPyBmWyIPwZSpDekXKMePxAfd9Yxx8O/ArYsUdbfBj4q/pa3h04uOcx8Ux6Hts654ntWPsG7PLOtGupAahu2yNGLONPKOexRwFBCfG7d6btWrfj2fX1um/Pff3izmv9OMrxcVQfHVm2x3aeSd0n6/MtKOHvDcDWlHPo1cCTRmz/mvmnes1HzP8F4B9Ye5vmIzrTTq5t9LRar20Ycu6vZR9G2acfXte9lPKmcWvgCbVNHtCp88h+QrlNY+c67QTK/nD3UduwqT/mvAKb8oPRQfgS6hXSgU7+RsqBa9SVn2FBeK8h47onmFM60/elvDOdx4YH4S8Af96Z9oC6Q23ZqceCzvRvAkcN2a55lJC8b2fcS4AL6/ATgas7074GHF2H/5V6Eu9M/wFrTxTXAC+a5PV5NvDVgXHvBf5mRPl3AG+vw68FPjmi3JnA6Z3nhwPfH1F26HJ6tMsxwP92pj24tvnvdcbdyLrhvHsguiflyttCShD75sD6L2btlekLgdd3pv05a9+MvIb6Bqgz/QLq1egp5p3oJ5MF4R8Bh3eePwm4ps/8w6bXdlvaeX6PWuY+lBPF7XTeNFHC6Zc68140qq6DfWSq+vRs9zdOsb5rqAGr09d+VIcPoezv3SuRk7XnIZQT5bad6R8F/rrH/jCxbQ/sTP8H4P11+GRmHoRfA3yoDu9EOfHfdwb7zCH0C8LrbT8lgP2Sda9qHQT8eNiy6+vy/IG2OLUOn0EN0PX53owIwsB9KVex1ws6dZ2/Zt3+fT3lE5qp2uIs4DQ6x+g6ftJjItM4tg2p7w51O7evz39a63SvKea7AHhlz3VcDhzZ6Uuj9vVFQ17rs4f10anK9tjOM1k3CD8a+OnAPK8F/m3E8tbMP9lrPmLeRbX/TJwHLgDe2Zl+MvCVgXlOHtw2ypvma6jncMpFuuuonzDUcR9m7dXiafUTygWAh/Z5jTfFh/cIz8xulHdLg/4f5V3W5yLi6og4qceylk1j+k8o7+R26VXLye1al9dd9paUMDGh+83sXzH8ywe7UN5RDi5rtzr8RWCbiHh0ROxOeVf6yTptd+CEiPjFxIMS7HbtLGuy9tkdePTA/M+jHCip6/xSRKyMiFXA8axtu4WUUDFKn22fbDlTtQuUK3ITfg2QmYPjuutd0xaZeRulD+7K+q/lsHWN2p7dgT8ZaMODKSfwqebtY1g/23VE2b7W1Cczf1UH70nZlq2Aazvb8l7K1dOhpugjU+nT7lPt34NlBttnZWb+ZpJ1Dpa/OTN/OWx6z22drC4z9UHgqRFxT+BZlKB27ZByffaZqYza/vmUIHVZp2/8Vx0/yqh+vyvrttNkr/FC4KbMvHnE9Btz3S+CTqxnqrY4kRLuvxkRV0TEi+r4SY+JU2zXOiJiXkScEhE/iohbKEEK1vaZZ1AC0k8i4ssRcdCIbRx5rI2IoyPi8k5d92fdPjlqX9+V4a/1MJOW7bGdg3YHdh1o49ex7rlzMqNe82FeAFyVmZfX5x8Cnlu/uDth0mNMLfsx4OzMPKeO3hVYlvVL/1XfcwYRcUJEXBURq+r2b894csmcMAhPU0Q8itJZ1vsJp8y8NTNPyMy9gKcC/yci/nBi8ohFjho/YWFneBHlqu0NlKsb9+jUax7rHtSnWu4Kyg7dXfZq1g1nfdxQ6zS4rJ8B1B3to5Qrc88FPp2Zt9Zyyyi3TezQedwjMz/cczuWAV8emP+emfnSOv1sykehCzNze8ptK9GZ937T3NZRdRi2nEnbZYbW9IUaKnZi7b2iuw+U7buuZZQrwt023DYzT+kx71R9jCF1W1TH9dFn+V3LKFfRdulsy70yc79JljdZH5lKn3bvsw2D+3i3fQbnn6o9d6xfrhk2vc+2TlaXPtbb3sz8GeVK+dMpJ/YPjJh3HPvMqO2/gfLGcr9O39g+M2fyywLXAgs6zxeOKkjpkztFxA7TXMdUx9XrMvPPMnNXylXZ90TE3kx9TJyO5wJHUj7V255ylRVqn8nMSzPzSMobzXMpx/lhhh4j64WR9wEvp9xOtAPlNqE++9+1DH+tZ1J20u1k/T69jPJJQreNt8vMw3vUe7qOBvaKiOvqTwa+jRI4u78ANdUx5l8otz10f/ViBbBw4Au8vfa1iHgs5VOeZ1E+6diBcm913+PmJscg3FNE3CsijgDOoXzs8N0hZY6IiL0jIig3kd9ZH1AC5l4zWPXzI2LfiLgH5daLj2XmnZT7xe4eEU+p7/heT7m/acLPgT0m+ab6h4FXRcSeNVT9HeULF9P6mapal48Cb4mI7erB7f9QrgJNOJvykd3z6vCE9wHH1ytVERHb1u3ZrufqPw3cPyJeEBFb1cejIuJBdfp2lKsxv4mIAygHvAkfAp4YEc+KiC0jYueIeNh0tn2y5fRsl+k6PCIOjoitKV9g+0ZmLgPOp7TDc2sdnk25jebTPZY5cbXuSfXKyN2j/L7pginnLPeb3cXk/frDwOsjYn5E7EK5r65vG/RZ/hr1KuPngH+q++sWEXG/iHh8LfJzYEFtvwmT9ZGpbEi7d70sIhZExE6UK0sfmaRsn/b824jYup6wjqB8EQj6betfR8Q9ImI/yn2mk9VlmGFtDOWj/BMptwB9cr256H0s6WO97a9vyN8HvD0i7g0QEbtFxJOmuWxqHY+NiAfV4/IbRhWsffKzlKC6Yz1GPW6qFUzVFhHxJ5199GZKGLqTqY+J07Ed5Y3ljZSLLmt+lqu27/MiYvvM/C1rz3fDnA68OiIeWY/ze9ft2bbWe2Vd5rGUK8JTysyfAEtY+1ofTLn4NJOyI7ezGjx3fxO4JSJeExHb1OPm/lEuko1NlCvs96Pcp/uw+tifcg59Yc9lvIRyr/xzB67+foNyMe3E2kcOobTJOestZH3bUS6arQS2jIg3APfqU59NlUF4ap+KiFsp7wL/ivKO7NgRZfcB/pvy7cyLgfdk5oV12lspJ7BfRMSrp7H+D1Du17mO8qWIvwDIzFWU+zVPp7yL+yXl26gTJk5+N0bEt4Ys94y67K9QvvX6G+AV06hX1yvq+q+mXCk/uy6fWteJnW5XyklhYvwSyhch3kU5mC+l3BfWS72y/EfAUZR3uNex9otFUNrnjfX1ewOdKxZZfmfycMqN/jdR7k17aN9191zOpO0yA2cDf1PX80jqT+lk5o2UE/4JlIP5icARmXlDj/ovo1wNeR3lwLYM+L/0ODbUjyrfAnyt9usDhxR7M+Uk9B3gu8C36rgp9Vz+oKMpHylP/ArDx1h7m8cXKd+Ovi4iJtpmZB/pUb8Zt/uAsykB/ur6mKx9pmrP6yjbvYLyJu34zPx+ndZnW79M2Q+/APxjZn5umtsyrI2hhN/dKffT/3LonMWG7jOTbf9rKNt2SZSPwP+bGfx2cGZ+Fvhn4Et1eRfXSbePmOUFlKu736fcD/qXPVc1WVs8CvhGRNxGucr/ysz8cY9j4nScRfm4/GeU/emSIdt1TW3L4ylfoFpPZv4HZT8+m3Jl8lzKrx9dCfwTpf1+TnmT9LVp1O+5lPt1b6IcF8+aYdmptvP9wL71GHRufZPyVEow/THl6v3plKvJ4/RC4D8z87v1E4DrMvM6yhelj6hvnKfyHEqIXxHlH4zcFhGvy8w7gD+mXFm+AXgP5bs7359kWRMuoJzHf0hpt9/Q7xawTdbEty8lbcJiA3/UXZum6PGPdjYXEfEjyi94bFbbWq+2fo/yazcb/I9fJG1cXhGWJM2qiHgG5WPwL851XcYhIp5eP2bfkXLF9VOGYOl3k0FYkjRrIuJCyk8lvmzgPsXfZS+h3Er0I8q9sTP5MpqkTYC3RkiSJKlJXhGWJElSkwzCkiRJatKWc7XiXXbZJffYY4+5Wr0kSZIacdlll92Qmev9N8k5C8J77LEHS5YsmavVS5IkqRERMfTfcHtrhCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQVjS+H3preUhSdImbMu5roCkjeTk7Tf+Or98Sr9yJ6+a3XpIkjSEV4SlVtztXuXv9gvLcNTdf97WZXj3xwyfL7ZYd3j7hesvc2J5E2W7y++Wm1hfd3mPP2lm2yNJ0gYyCEutuP2W8nfVsjKcd5Xnd95Rhn/yteHzTZSbGF61bP1lTixvomx3+d1yE+vruuQ9098WSZLGwCAsae7kXeuGZEmSNiKDsNSK2GLd2yMmdG9VmAvd2yYkSdqIen1ZLiIOA94JzANOz8xTBqYfAvwn8OM66hOZ+cbxVVPSBnvcieXvNV8tf1/1PXj7/rDDItjjsfAHry3j/+3w8nyi3MTwdd+F+zy4PO8up1t2wrHnl1+N+IPXlnU87Hnr1+ear5ZykiTNkcjMyQtEzAN+CBwKLAcuBZ6TmVd2yhwCvDozj+i74sWLF+eSJUtmUGVJkiSpv4i4LDMXD47vc2vEAcDSzLw6M+8AzgGOHHcFJUmSpI2pTxDeDeh8TZzlddyggyLi2xHx2YjYbyy1kyRJkmZJn3uEY8i4wfspvgXsnpm3RcThwLnAPustKOI44DiARYsWTa+mkiRJ0hj1uSK8HOh8xZwFwIpugcy8JTNvq8PnA1tFxC6DC8rM0zJzcWYunj9//gZUW5IkSdowfYLwpcA+EbFnRGwNHAWc1y0QEfeJiKjDB9Tl3jjuykqSJEnjMmUQzszVwMuBC4CrgI9m5hURcXxEHF+LPRP4XkR8G/hn4Kic6ucoNhFv//wPefvnf7je8LBy3b8Tw4855Qsjyw4ODzNsfkmSJM2+KX8+bbbM5c+n7fXaz3DXiM3eIuCuhK3nBXfcOXnbTJQd9Xy6rjnlKTOfWZIkSUNtyM+nbXYmC6sT06YKwcOWsyEhGGCPkz6zYQuQJElSb00G4U3V1vOG/UCHJEmSZkOvf7G8uZm4BeHZ7714zbiPvOQgnv3ei/nISw5ac1/vJVffyIF77QzAqw69P4855Qss2PEe68zzmFO+wNdO+sM19/o+85EL+dhly3jmIxfyqkPvz7PfezEH7rXzOsMTyz3joqu51zZbsfLW23npIXvzqkPvv7GaQJIkqXlN3iO8qd6C8Og9d+IjLzlorqshSZK0WRl1j3CTV4Rf+Yf7rLkaC3Dtqt+w7dbz2HfX7QHWXLUFuPSam3jFE/ZZc0UYYOWtt/PwRTty5YpVa5Z5++q7ePiiHQG4csUq9t11+zV/l9/8qzVXkg/ca2f+5Yv/y6P22GnN8zMuuprbV9/F8pt/tXEaQJIkSW1eEZYkSVI7/NUISZIkqcMgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUpF5BOCIOi4gfRMTSiDhpknKPiog7I+KZ46uiJEmSNH5TBuGImAe8G3gysC/wnIjYd0S5vwcuGHclJUmSpHHrc0X4AGBpZl6dmXcA5wBHDin3CuDjwPVjrJ8kSZI0K/oE4d2AZZ3ny+u4NSJiN+DpwKmTLSgijouIJRGxZOXKldOtqyRJkjQ2fYJwDBmXA8/fAbwmM++cbEGZeVpmLs7MxfPnz+9ZRUmSJGn8tuxRZjmwsPN8AbBioMxi4JyIANgFODwiVmfmueOopCRJkjRufYLwpcA+EbEn8DPgKOC53QKZuefEcEScCXzaECxJkqRN2ZRBODNXR8TLKb8GMQ84IzOviIjj6/RJ7wuWJEmSNkV9rgiTmecD5w+MGxqAM/OYDa+WJEmSNLv8z3KSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqUq8gHBGHRcQPImJpRJw0ZPqREfGdiLg8IpZExMHjr6okSZI0PltOVSAi5gHvBg4FlgOXRsR5mXllp9gXgPMyMyPiIcBHgQfORoUlSZKkcehzRfgAYGlmXp2ZdwDnAEd2C2TmbZmZ9em2QCJJkiRtwvoE4d2AZZ3ny+u4dUTE0yPi+8BngBcNW1BEHFdvnViycuXKmdRXkiRJGos+QTiGjFvvim9mfjIzHwg8DXjTsAVl5mmZuTgzF8+fP39aFZUkSZLGqU8QXg4s7DxfAKwYVTgzvwLcLyJ22cC6SZIkSbOmTxC+FNgnIvaMiK2Bo4DzugUiYu+IiDr8CGBr4MZxV1aSJEkalyl/NSIzV0fEy4ELgHnAGZl5RUQcX6efCjwDODoifgv8Gnh258tzkiRJ0iYn5iqvLl68OJcsWTIn65YkSVI7IuKyzFw8ON7/LCdJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDWpVxCOiMMi4gcRsTQiThoy/XkR8Z36+HpEPHT8VZUkSZLGZ8ogHBHzgHcDTwb2BZ4TEfsOFPsx8PjMfAjwJuC0cVdUkiRJGqc+V4QPAJZm5tWZeQdwDnBkt0Bmfj0zb65PLwEWjLeakiRJ0nj1CcK7Acs6z5fXcaP8KfDZDamUJEmSNNu27FEmhozLoQUj/oAShA8eMf044DiARYsW9ayiJEmSNH59rggvBxZ2ni8AVgwWioiHAKcDR2bmjcMWlJmnZebizFw8f/78mdRXkiRJGos+QfhSYJ+I2DMitgaOAs7rFoiIRcAngBdk5g/HX01JkiRpvKa8NSIzV0fEy4ELgHnAGZl5RUQcX6efCrwB2Bl4T0QArM7MxbNXbUmSJGnDRObQ231n3eLFi3PJkiVzsm5JkiS1IyIuG3aR1v8sJ0mSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCb1CsIRcVhE/CAilkbESUOmPzAiLo6I2yPi1eOvpiRJkjReW05VICLmAe8GDgWWA5dGxHmZeWWn2E3AXwBPm41KSpIkSePW54rwAcDSzLw6M+8AzgGO7BbIzOsz81Lgt7NQR0mSJGns+gTh3YBlnefL67hpi4jjImJJRCxZuXLlTBYhSZIkjUWfIBxDxuVMVpaZp2Xm4sxcPH/+/JksQpIkSRqLPkF4ObCw83wBsGJ2qiNJkiRtHH2C8KXAPhGxZ0RsDRwFnDe71ZIkSZJm15S/GpGZqyPi5cAFwDzgjMy8IiKOr9NPjYj7AEuAewF3RcRfAvtm5i2zV3VJkiRp5qYMwgCZeT5w/sC4UzvD11FumZAkSZJ+J/if5SRJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktSkLfsUiojDgHcC84DTM/OUgelRpx8O/Ao4JjO/Nea6SpI0LQ/+9wfPdRVG+u4LvzvXVZCaN2UQjoh5wLuBQ4HlwKURcV5mXtkp9mRgn/p4NPCv9a8kSXPGsClpMn2uCB8ALM3MqwEi4hzgSKAbhI8EzsrMBC6JiB0i4r6Zee3YayxJkqSRrnrgg+a6CiM96PtXzXUV1tEnCO8GLOs8X876V3uHldkNWCcIR8RxwHEAixYtmm5dJUmSNIVNLWxuyvp8WS6GjMsZlCEzT8vMxZm5eP78+X3qJ0mSJM2KPkF4ObCw83wBsGIGZSRJkqRNRp8gfCmwT0TsGRFbA0cB5w2UOQ84OooDgVXeHyxJkqRN2ZT3CGfm6oh4OXAB5efTzsjMKyLi+Dr9VOB8yk+nLaX8fNqxs1dlSZIkacP1+h3hzDyfEna7407tDCfwsvFWTZIkSZo9/mc5SZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1Kco/hZuDFUesBH4yJysfr12AG+a6Eo2y7eeObT93bPu5Y9vPHdt+7mwubb97Zs4fHDlnQXhzERFLMnPxXNejRbb93LHt545tP3ds+7lj28+dzb3tvTVCkiRJTTIIS5IkqUkG4Q132lxXoGG2/dyx7eeObT93bPu5Y9vPnc267b1HWJIkSU3yirAkSZKaZBDuISLOiIjrI+J7I6ZHRPxzRCyNiO9ExCM2dh03Vz3a/pCIWBURl9fHGzZ2HTdXEbEwIr4UEVdFxBUR8cohZez7s6Bn29v3Z0FE3D0ivhkR365t/7dDytjvZ0HPtrffz6KImBcR/xMRnx4ybbPs91vOdQV+R5wJvAs4a8T0JwP71MejgX+tf7XhzmTytgf4amYesXGq05TVwAmZ+a2I2A64LCI+n5lXdsrY92dHn7YH+/5suB14QmbeFhFbARdFxGcz85JOGfv97OjT9mC/n02vBK4C7jVk2mbZ770i3ENmfgW4aZIiRwJnZXEJsENE3Hfj1G7z1qPtNUsy89rM/FYdvpVycNxtoJh9fxb0bHvNgtqXb6tPt6qPwS/T2O9nQc+21yyJiAXAU4DTRxTZLPu9QXg8dgOWdZ4vx5PWxnRQ/SjtsxGx31xXZnMUEXsADwe+MTDJvj/LJml7sO/Pivrx8OXA9cDnM9N+v5H0aHuw38+WdwAnAneNmL5Z9nuD8HjEkHG+i904vkX5t4kPBf4FOHduq7P5iYh7Ah8H/jIzbxmcPGQW+/6YTNH29v1Zkpl3ZubDgAXAARGx/0AR+/0s6dH29vtZEBFHANdn5mWTFRsy7ne+3xuEx2M5sLDzfAGwYo7q0pTMvGXio7TMPB/YKiJ2meNqbTbqfXofBz6UmZ8YUsS+P0umanv7/uzLzF8AFwKHDUyy38+yUW1vv581jwH+OCKuAc4BnhARHxwos1n2e4PweJwHHF2/UXkgsCozr53rSrUgIu4TEVGHD6D06Rvntlabh9qu7weuysy3jShm358Ffdrevj87ImJ+ROxQh7cBngh8f6CY/X4W9Gl7+/3syMzXZuaCzNwDOAr4YmY+f6DYZtnv/dWIHiLiw8AhwC4RsRz4G8pN/GTmqcD5wOHAUuBXwLFzU9PNT4+2fybw0ohYDfwaOCr9LzHj8hjgBcB36z17AK8DFoF9f5b1aXv7/uy4L/DvETGPErI+mpmfjojjwX4/y/q0vf1+I2qh3/uf5SRJktQkb42QJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkpr0/wE1UVQw9P2zowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "arizona_samples = arizona_results.extract()\n",
    "\n",
    "# Make a new array with same dimensions as alpha\n",
    "p_predicted_AZ = np.empty(arizona_samples['alpha'].shape)\n",
    "\n",
    "# Generate one p sample for each alpha sample\n",
    "for i in range(arizona_samples['alpha'].shape[0]):\n",
    "    p_predicted_AL[i] = sts.dirichlet(arizona_samples['alpha'][i]).rvs()\n",
    "    \n",
    "plt.figure(figsize=(12, 6))\n",
    "for i in range(4):\n",
    "    plt.plot(sts.uniform.rvs(loc=i+1-0.1, scale=0.2, size=arizona_samples['theta'].shape[0]), arizona_samples['theta'][:,i], ',', alpha=0.5)\n",
    "plt.title('Distribution over each component of theta for probability of being chosen as candidate in Arizona')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAF1CAYAAADiNYyJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAo9klEQVR4nO3debgkVX3/8feXGRZBVhkXmGFTiCIo0csWMZBEIiAGfTQBQVGMAYwYY/Qn4EpcIskvEbNoABEJbmhcCCqKBkFFAbkYXADhNw7IjIAM+yICA9/fH6fuTE1P9711Z/pOT+a8X8/Tz+2qOlV16tSp6k9XV/eNzESSJEmqzTqjroAkSZI0CgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzC0xARp0bEu4a0rG0i4v6ImNUMXxwRrxvGspvlfT0iXj2s5a2tIuKsiHj/qOtRg4h4fUT8uun3T+hQ/jURccnqqNuA9U+rvqtTROwXEYtWct7tIiIjYvaA6W+PiDP6lV2d55WIeH9E3B4Rt/aZtirb//yIuG7Vazitda50fbWiiDgpIj7VPF/utXSysquhXkPLCB3Xt9q2bW1mEG5ExI0R8WBE3BcRd0fEDyLi2IhY2kaZeWxmvq/jsl4wWZnMvCkzH5+Zjw6h7iscDJl5YGb+x6ouW4Kpw1OH+dcFPgT8cdPv7xjm8vusb5VC9FT1XZtl5t9lZt835e3zyky+UYmIecBbgJ0z88nDXHZmfi8zf2eYy9ToDPm1dJUujHTNCB3qsCQitlqV5ag7g/DyXpyZGwPbAicDxwMfH/ZKhvViv7ayfdZKTwI2AK4edUU6Wun6RrFK51aPAbYF7sjM20ZdEWl1iYiNgJcB9wBHTFG29nPE8GSmj/Lf9W4EXtAzbg/gMWCXZvgs4P3N8y2BrwJ3A3cC36O8sfhkM8+DwP3A24DtgAT+HLgJ+G5r3OxmeRcDHwR+SDkI/gvYopm2H7CoX32BA4CHgUea9f24tbzXNc/XAd4J/BK4DTgb2LSZNlGPVzd1ux14xyTttGkz/+Jmee9slr9+0xa7tMrOadrhic3wwcBVTbkfAM/q2Z7jgZ8AD020S8+6nw58q2nv64A/a017EfA/wL3AQuCknnn3adZ5dzP9Na19+hHga8B9wOXAUyfZ/kHL6dsuzbTXAN8HTmnmWwD8XjN+YbNPXt1ax1nAqc223gd8B9i2Nf33gCso/eQK4Pda0y4G3tes7z7gm8CWrel7ter/Y2C/LvM2fSMpfex+YO8+bbM+8GHg5ubx4WbcTsADrfm/3WfeFZbftM8lwD8CdwE3AAf29MWPA7cAvwLeD8wCngH8Fni0WdbdXfpIa7l969uh3T/QtN2DwNMGnGNOBK5ptucTwAbtY5xyDNxKOY/0bc+e8m+nHLM3Akd0OR5Ydswf3Sz3FuAtreknAZ/qKds+T72uXxsDuwO/pnXsUl7Ur5rmueQFTRs+1iz7rD7zTrX961P6zU1NnU4FHtfvfNrM+1bKuece4HMT+6WZ/ramjW5utj377d+m7BbNfr252cfn9tT3LZTj/RbgqKnaopn2NMo54J5mWz/X8Zx4FtM7t/0npe/dQ3mNemZr2kGUfnsf5Vh76yTL+Qvg2qbsNcBzmvEnAL9ojX9pa57XMPmxvn3TBvc12/tvDO6jA8tOtp2U4+ERyuvp/cBXmvFbAV9s9s0NwF9Nsu1nsSwjTLrPB8x/JOV4fRPws55pJwFfAD5FOa5fx/LH6rrAZ5u6rgcc1doPC4Bj+hw/g/pj3/0NbE7JPYub/fRVYG7PflzQzHcDrWNyTX6MvAJryoM+QbgZfxPw+uZ5u5N/kHJyXbd5PB+IfstqHahnAxsBj+tz8F7cdLhdmjJfbHXw/RgQhJvnSw+G1vSLWRaEXwvMB3YAHg98CfhkT90+1tTr2ZQg+owB7XQ2JaRv3Mx7PfDnzbQzgQ+0yr4B+Ebz/DnNAbcnJay8utmG9VvbcxUwj+YFq2e9G1FOEEcBs5vl3c6yk9h+wK6UF9JnUV78XtJM26Y5MF/R7KsnALu19umdlDc9s4FPA+cM2PbJljNZu7wGWNLUfRYlsN1EeZFaH/jjZrmPb9XpPuD3m+n/DFzSTNuCcgJ6VVPfVzTDT2jt919QwtzjmuGTm2lbA3dQTnLrAPs3w3M6zLsdrf46oH3eC1wGPJHyJugHwPu6zN9vetNuj1BeWGcBr6cEjInj7FzgNErfeCLlTeQxrXkv6VnHfgzoI1PVp2O73wQ8s5m+7oBzzM8ofXwLSmhuv2guAf6+2eePm6I9J8p/qCm/LyW8/06H42Fi2z7btN2ulBe2Fc4nfdrhYpadV/q18TUsH2C+TCtkT+Ncsh8957w++3Ky7f8wcF7TzhsDXwE+2G/ZzX75ISXwbEEJD8c20w6ghKZnAhtS3qBMFoS/RgnSm1POEfv21Pe9zfiDgN8Am3doi88C72j25QbAPh3PiWfR8dzWep3YmGVvwK5qTbsFeH7zfHOacNtnGX9KeR3bHQhKiN+2NW2rZjsObfbXUzoe65e29vXvU86Pg/rowLIdtvMsmmOyGV4HuBJ4NyVc7kAJei8csP1L559qnw+Y/0LgHyifSC1ptzPluHwEeElTr8c14z7VPP9as/5ZTfkXAU9t9sO+zbqf06Vug/Y35TXvZZRjYWPKm4pzW/3xXpYdg0+h9WZqTX6MvAJryoPBQfgymiukPZ38vZQT16ArP/2C8A59xrVfYE5uTd+Z8s50FqsehC8E/rI17XeaA2p2qx7td3U/BA7rs12zKCF559a4Y4CLm+cvABa0pn0fOLJ5/u80L+Kt6dex7IXiRuC1k+yfQ4Hv9Yw7DXjPgPIfBk5pnp8IfHlAubOAM1rDBwE/H1C273I6tMtrgP/XmrZr0+ZPao27g+XD+TmtaY+nXHmbRwliP+xZ/6UsuzJ9MfDO1rS/ZNmbkeNp3gC1pl9AczV6inkn+slkQfgXwEGt4RcCN3aZv9/0pt3mt4Y3bMo8mfJC8RCtN02UcHpRa95LBtW1t49MVZ+O7f7eKdZ3I03AavW1XzTP96Mc7+0rkZO1536UF7KNWtM/D7yrw/EwsW1Pb03/B+DjzfOTWPkgfDzw6eb5FpQX16esxDGzH92C8ArbT3nhf4DW1U/KJww39Ft2s19e2dMWpzbPz6QJ0M3w0xgQhCkv/I/RJ+g063yQ5fv3bZRPaKZqi7OB02mdo5vxk54Tmca5rU99N2u2c9Nm+KamTptMMd8FwJs6ruMq4JBWXxp0rG/TZ19/pl8fnapsh+08i+WD8J7ATT3znAh8YsDyls4/2T4fMO82Tf+ZeB24APjn1vSTgO/2zHMS5Q3fd4B/oXnjMGD5507sm6nqNo39vRtwV/N8I8onQy+jz8WsNfnhPcJT25ryrrrX/6VcZf1mRCyIiBM6LGvhNKb/kvJObctOtZzcVs3y2sueTQkTE9rfzP4NJXz12pLyrrh3WVs3z78NPC4i9oyIbSkHyZebadsCb2m+iHh3RNxNCXbtLwRM1j7bAnv2zH8E5URJs86LImJxRNwDHMuytptHCRWDdNn2yZYzVbtAuSI34UGAzOwd117v0rbIzPspfXArVtyX/dY1aHu2Bf60pw33obyATzVvF/362ap+4WNpfTLzN83Tx1O2ZV3glta2nEa5etrXFH1kKl3afarju7dMb/sszszfTrLO3vJ3ZeYD/aZ33NbJ6rKyPgW8OCIeD/wZJajd0qdcl2NmKoO2fw4lSF3Z6hvfaMYPMqjfb8Xy7TTZPp4H3JmZdw2YfkdmLumznqna4m2UcP/DiLg6Il7bjJ/0nDjFdi0nImZFxMkR8YuIuJfy5gCW9ZmXUYL0LyPiOxGx94BtHHiujYgjI+KqVl13Yfk+OehY34r++7qfSct22M5e2wJb9bTx21n+tXMyg/Z5P68Crs3Mq5rhTwOHN1/cndCv/+1F+dTn5GwSKUBEHBgRl0XEnU29D2L57Zysbn33d0RsGBGnRcQvm/b7LrBZRMxq2vxQyrnmloj4WkQ8fcC2rlEMwpOIiN0pJ6MVvhmdmfdl5lsycwfgxcDfRMQfTUwesMhB4yfMaz3fhnLV9nbK1Y0NW/WaxfIn9amWezPlgG4vewnLh7Mubm/q1LusXwFk5mOUqzKvAA4HvpqZ9zXlFlJum9is9dgwMz/bcTsWAt/pmf/xmfn6ZvpnKO+M52XmppTbVqI171Onua2D6tBvOZO2y0pa2heaULEFy+4V3banbNd1LaRcEW634UaZeXKHeafqY/Sp2zbNuC66LL9tIeUq2patbdkkM585yfIm6yNT6dLuXbah9xhvt0/v/FO15+bNl2v6Te+yrZPVpYsVtjczf0W5Uv5Sygv7JwfMO4xjZtD23055Y/nMVt/YNDOn86Zuwi3A3NbwvEEFKX1yi4jYbJrrmOq8emtm/kVmbkW5SvfRiHgaU58Tp+Nw4BDKp3qbUq6yQtNnMvOKzDyE8kbzXMp5vp++58jmwsjHgOMotxNtRrlNqMvxdwv99/XKlJ10O1mxTy+kfJLQbuONM/OgDvWeriOBHSLi1ig/GfghSnA9sFWm3znmm5RbNS+MiCcBRMT6lNsr/5HyyeNmwPl0PN9Nsr/fQvlEec/M3IRy6wks6ycXZOb+lIsrP6fs8zWeQbiPiNgkIg4GzqF8pPLTPmUOjoinRURQ7ot5tHlACZg7rMSqXxkRO0fEhpRbL76Q5Sdhrgc2iIgXNe8O30m5v2nCr4HtJvmm+meBN0fE9k2o+jvKFy6WDCjfV1OXzwMfiIiNm5Pb31CuAk34DOVd4RHN8wkfA45trlRFRGzUbM/GHVf/VWCniHhVRKzbPHaPiGc00zemXI35bUTsQTnhTfg08IKI+LOImB0RT4iI3aaz7ZMtp2O7TNdBEbFPRKxH+QLb5Zm5kHIy2ykiDm/qcCjlNpqvdljmxNW6FzZXRjaI8vumc6ecs9xD+hiT9+vPAu+MiDkRsSXlvrqubdBl+Us1Vxm/CfxTc7yuExFPjYh9myK/BuY27Tdhsj4ylVVp97Y3RMTciNiCcmXpc5OU7dKefxsR60XE8ylfRv3PZnyXbX1Xc4XnmZT7TCerSz/92hjKR/lvo9wC9OUV5qLzuaSLFba/eUP+MeCUiHgiQERsHREvnOayaep4VEQ8ozkvv3tQwaZPfp0SVDdvzlG/P6h8a75J2yIi/rR1jN5FCUOPMvU5cTo2pryxvINy0eXvJiY07XtERGyamY+w7PWunzOAt0bEc5vz/NOa7dmoqffiZplHUa4ITykzfwmMs2xf70O5+LQyZQduZ6P3tfuHwL0RcXxEPK45b+4S5SLZ0ES54vpUyv3cuzWPXSivoa+eav7M/Iem7IXNuWI9SkZYDCyJiAMp30XpUpfJ9vfGlDeZdzfnsPe05ntSRPxJlDchD1G+cLjKP2m3OhiEl/eViLiP8i7wHZR3ZEcNKLsj8N+UnX0p8NHMvLiZ9kHKC9jdEfHWaaz/k5R7jG6lfCnirwAy8x7K/ZpnUK4SPED5xueEiRe/OyLiR32We2az7O9Svsn5W+CN06hX2xub9S+gXCn/TLN8mrpe3kzfivKiMDF+nPJFiH+jnMznU+4L66S5svzHwGGUqz63suyLRVDa573N/ns3rSsWmXkT5WOet1BuMbiK8qXAaZliOZO2y0r4DOUkcyfwXJqf0snye7YHN3W4gxI4Ds7M2zvUfyHlasjbKSfIhcD/ocN5oPmo8gPA95t+vVefYu+nvAj9BPgp8KNm3JQ6Lr/XkZQT/sSvMHyBZbd5fJvy02e3RsRE2wzsIx3qt9Lt3uMzlAC/oHlM1j5TteetlO2+mfIm7djM/Hkzrcu2fodyHF4I/GNmfnOa29KvjaGE320p99M/0HfOYlWPmcm2/3jKtl0W5SPc/6ZcyZqWzPw65d7Li5rlXdpMemjALK+iXN39OeWey7/uuKrJ2mJ34PKIuJ9ylf9NmXlDh3PidJxNuYXgV5Tj6bI+23Vj05bHAq/st5DM/E/KcfwZypfUzqX8+tE1wD9R2u/XlDdJ359G/Q6n3K97J+W8ePZKlp1qOz8O7Nycg85t3qS8mBJMb6BcvT+DcjV5mF4N/Fdm/rT5BODWzLyV8kXpg5vQOaksv198LqWvr0vJD5+nHCOHU/pOV4P294cpX8y7ndJ232jNsw7l/Hgzpe33pZyH1ngT38iUtIaIiLMoX+Z556jrouGJiBspXzT771HXZaZFxC8ov+CxVm1rc7X1Z5Rfu5nWJ2qS1kxeEZYkDU1EvIzyMfi3R12XYYiIlzYfF29OueL6FUOwtPYwCEuShiIiLqb8VOIbmnt11wbHUG4l+gXlnseV+TKapDWUt0ZIkiSpSl4RliRJUpUMwpIkSarS7FGteMstt8zttttuVKuXJElSJa688srbM3OF/zA5siC83XbbMT4+PqrVS5IkqRIR0fdfc3trhCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQbjHKd+6ftRVkCRJ0mpgEO7x5v13GnUVJEmStBoYhHuc8q3r+14V7h3Xb7g97tDTLu27/EHjJUmStHpVGYQnAutEeG0H2Dfvv9NyV4Wfd/KFS8f3MxFsJ+abWN5eOzyhb/nPHbM3h5526ZTBWpIkSTOryiD85v134nknX8hlC+5Y+pgIohPBdmL4+yf8EQC7vucbKyzn0NMuXRp4T/nW9ez0jvOXC8wTgbd9FfjQ0y7lc8fsvXSdu77nGyusU5IkSTMvMnMkKx4bG8vx8fGRrHvX93yDh5Y8xu9uszmfO2bvpeO3O+Fr7Ln9FgAsuus33PvgI0vLTYx7+XPnLQ3SL3/uPM68ZAGv3WcHvnDlQuZuviF77fAELltwB3vt8ISl4z53zN6c8q3rl45ve/P+O3HoaZdyzc338NO/PWD1NYIkSVIlIuLKzBxbYXytQfiBhx9l9+224PIb7uTGk1/EDid+jccSNl5/Fg8teYzrP3AQh552KZffcCd7br8F/3PTXQDM2Xh95m6+IZffcCcbrz+LBx5+FIDZ6wTrz16H1+6zA5ctuINrbr6H+x56lHUCHktYJ+CNf7gjX7hyIb+6+7dLl/nwo8k6ARutN4vX7rODX9aTJEkaMoNwy3YnfG0k6+3ixpNfNOoqSJIkrVUGBeEq7xFeU603K0ZdBUmSpGoYhNcgDz86mqvzkiRJNTIID9B7dXad1XCx9k1/tOPMr0SSJEkAzB51BUbB+3AlSZLkFWFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVeoUhCPigIi4LiLmR8QJA8rsFxFXRcTVEfGd4VZTkiRJGq4p/6FGRMwCPgLsDywCroiI8zLzmlaZzYCPAgdk5k0R8cQZqq8kSZI0FF2uCO8BzM/MBZn5MHAOcEhPmcOBL2XmTQCZedtwqylJkiQNV5cgvDWwsDW8qBnXthOweURcHBFXRsSR/RYUEUdHxHhEjC9evHjlaixJkiQNQZcgHH3GZc/wbOC5wIuAFwLvioidVpgp8/TMHMvMsTlz5ky7spIkSdKwTHmPMOUK8LzW8Fzg5j5lbs/MB4AHIuK7wLOB64dSS0mSJGnIulwRvgLYMSK2j4j1gMOA83rK/Bfw/IiYHREbAnsC1w63qpIkSdLwTHlFODOXRMRxwAXALODMzLw6Io5tpp+amddGxDeAnwCPAWdk5s9msuKSJEnSqojM3tt9V4+xsbEcHx8fybolSZJUj4i4MjPHesf7n+UkSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVXqFIQj4oCIuC4i5kfECX2m7xcR90TEVc3j3cOvqiRJkjQ8s6cqEBGzgI8A+wOLgCsi4rzMvKan6Pcy8+AZqKMkSZI0dF2uCO8BzM/MBZn5MHAOcMjMVkuSJEmaWV2C8NbAwtbwomZcr70j4scR8fWIeGa/BUXE0RExHhHjixcvXonqSpIkScPRJQhHn3HZM/wjYNvMfDbwr8C5/RaUmadn5lhmjs2ZM2daFZUkSZKGqUsQXgTMaw3PBW5uF8jMezPz/ub5+cC6EbHl0GopSZIkDVmXIHwFsGNEbB8R6wGHAee1C0TEkyMimud7NMu9Y9iVlSRJkoZlyl+NyMwlEXEccAEwCzgzM6+OiGOb6acCLwdeHxFLgAeBwzKz9/YJSZIkaY0Ro8qrY2NjOT4+PpJ1S5IkqR4RcWVmjvWO9z/LSZIkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKq1CkIR8QBEXFdRMyPiBMmKbd7RDwaES8fXhUlSZKk4ZsyCEfELOAjwIHAzsArImLnAeX+Hrhg2JWUJEmShq3LFeE9gPmZuSAzHwbOAQ7pU+6NwBeB24ZYP0mSJGlGdAnCWwMLW8OLmnFLRcTWwEuBUydbUEQcHRHjETG+ePHi6dZVkiRJGpouQTj6jMue4Q8Dx2fmo5MtKDNPz8yxzBybM2dOxypKkiRJwze7Q5lFwLzW8Fzg5p4yY8A5EQGwJXBQRCzJzHOHUUlJkiRp2LoE4SuAHSNie+BXwGHA4e0Cmbn9xPOIOAv4qiFYkiRJa7Ipg3BmLomI4yi/BjELODMzr46IY5vpk94XLEmSJK2JulwRJjPPB87vGdc3AGfma1a9WpIkSdLM8j/LSZIkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwpP4u+uBwy0mStIbp9KsRktYCJ206/Xm+c/Jwy/Uzaz14l/9yXZK0+nlFWNJoPfrwqGsgSaqUQViSJElVMghLkiSpSgZhqRb7ngAn3VPuyY11yvD6m5S/m84r43uf9xsX65Th9Tcp4zadt2x44u/EPBPj9j0Btn1eWf9E+W2fV/6edM+oW0aSVCmDsFSLPzgRTtkF9vkb2GZvuPF7cOJCuORDsNsRZfxlH132HGDu7svGzd6gzPOeu8q0vf6yfMlttyPKY4Pmy3hP3rX83e2I8ny3I8q6J+x2RPl7903w5p/BJw5aPdsvSVIPfzVCqsVFHyzBE4AmmJ6ySwm7V326NY1WcD2xzHfVp0uo3e75Jbhuts2ysld9ugxPBNyJ+S/6IBx1fvn7iYPK8wkTZSfKSJI0AgZhqRbtq7IXfbAMT1yRffPPlo1rmxi33LwsC7qw4hXf9nwTJsJu73h/g1iSNEKRmSNZ8djYWI6Pj49k3ZKmqV+A7Q2/K6P3SrEkSTMgIq7MzLHe8d4jLGl5/a7S9obeYYRgMARLkkbKICxpecMKuZIkreEMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqdQrCEXFARFwXEfMj4oQ+0w+JiJ9ExFURMR4R+wy/qpIkSdLwzJ6qQETMAj4C7A8sAq6IiPMy85pWsQuB8zIzI+JZwOeBp89EhSVJkqRh6HJFeA9gfmYuyMyHgXOAQ9oFMvP+zMxmcCMgkSRJktZgXYLw1sDC1vCiZtxyIuKlEfFz4GvAa/stKCKObm6dGF+8ePHK1FeSJEkaii5BOPqMW+GKb2Z+OTOfDrwEeF+/BWXm6Zk5lpljc+bMmVZFJUmSpGHqEoQXAfNaw3OBmwcVzszvAk+NiC1XsW6SJEnSjOkShK8AdoyI7SNiPeAw4Lx2gYh4WkRE8/w5wHrAHcOurCRJkjQsU/5qRGYuiYjjgAuAWcCZmXl1RBzbTD8VeBlwZEQ8AjwIHNr68pwkSZK0xolR5dWxsbEcHx8fybolSZJUj4i4MjPHesf7n+UkSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVXqFIQj4oCIuC4i5kfECX2mHxERP2keP4iIZw+/qpIkSdLwTBmEI2IW8BHgQGBn4BURsXNPsRuAfTPzWcD7gNOHXVFJkiRpmLpcEd4DmJ+ZCzLzYeAc4JB2gcz8QWbe1QxeBswdbjUlSZKk4eoShLcGFraGFzXjBvlz4Ov9JkTE0RExHhHjixcv7l5LSZIkaci6BOHoMy77Foz4A0oQPr7f9Mw8PTPHMnNszpw53WspSZIkDdnsDmUWAfNaw3OBm3sLRcSzgDOAAzPzjuFUT5IkSZoZXa4IXwHsGBHbR8R6wGHAee0CEbEN8CXgVZl5/fCrKUmSJA3XlFeEM3NJRBwHXADMAs7MzKsj4thm+qnAu4EnAB+NCIAlmTk2c9WWJEmSVk1k9r3dd8aNjY3l+Pj4SNYtSZKkekTElf0u0vqf5SRJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVeoUhCPigIi4LiLmR8QJfaY/PSIujYiHIuKtw6+mJEmSNFyzpyoQEbOAjwD7A4uAKyLivMy8plXsTuCvgJfMRCUlSZKkYetyRXgPYH5mLsjMh4FzgEPaBTLztsy8AnhkBuooSZIkDV2XILw1sLA1vKgZN20RcXREjEfE+OLFi1dmEZIkSdJQdAnC0WdcrszKMvP0zBzLzLE5c+aszCIkSZKkoegShBcB81rDc4GbZ6Y6kiRJ0urRJQhfAewYEdtHxHrAYcB5M1stSZIkaWZN+asRmbkkIo4DLgBmAWdm5tURcWwz/dSIeDIwDmwCPBYRfw3snJn3zlzVJUmSpJU3ZRAGyMzzgfN7xp3aen4r5ZYJSZIk6X8F/7OcJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKnYJwRBwQEddFxPyIOKHP9IiIf2mm/yQinjP8qkqSJEnDM2UQjohZwEeAA4GdgVdExM49xQ4EdmweRwP/PuR6SpIkSUM1u0OZPYD5mbkAICLOAQ4BrmmVOQQ4OzMTuCwiNouIp2TmLUOvsSRJHe36H7uOugoD/fTVPx11FbSWuvbpzxh1FQZ6xs+vHXUVltMlCG8NLGwNLwL27FBma2C5IBwRR1OuGLPNNttMt66SJE2LYVM1WtPC5pqsyz3C0WdcrkQZMvP0zBzLzLE5c+Z0qZ8kSZI0I7oE4UXAvNbwXODmlSgjSZIkrTG6BOErgB0jYvuIWA84DDivp8x5wJHNr0fsBdzj/cGSJElak015j3BmLomI44ALgFnAmZl5dUQc20w/FTgfOAiYD/wGOGrmqixJkiStui5fliMzz6eE3fa4U1vPE3jDcKsmSZIkzRz/s5wkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqoU5Z/CjWDFEYuBX45k5cO1JXD7qCtRKdt+dGz70bHtR8e2Hx3bfnTWlrbfNjPn9I4cWRBeW0TEeGaOjboeNbLtR8e2Hx3bfnRs+9Gx7UdnbW97b42QJElSlQzCkiRJqpJBeNWdPuoKVMy2Hx3bfnRs+9Gx7UfHth+dtbrtvUdYkiRJVfKKsCRJkqpkEO4gIs6MiNsi4mcDpkdE/EtEzI+In0TEc1Z3HddWHdp+v4i4JyKuah7vXt11XFtFxLyIuCgiro2IqyPiTX3K2PdnQMe2t+/PgIjYICJ+GBE/btr+b/uUsd/PgI5tb7+fQRExKyL+JyK+2mfaWtnvZ4+6Av9LnAX8G3D2gOkHAjs2jz2Bf2/+atWdxeRtD/C9zDx49VSnKkuAt2TmjyJiY+DKiPhWZl7TKmPfnxld2h7s+zPhIeAPM/P+iFgXuCQivp6Zl7XK2O9nRpe2B/v9THoTcC2wSZ9pa2W/94pwB5n5XeDOSYocApydxWXAZhHxlNVTu7Vbh7bXDMnMWzLzR83z+ygnx617itn3Z0DHttcMaPry/c3gus2j98s09vsZ0LHtNUMiYi7wIuCMAUXWyn5vEB6OrYGFreFF+KK1Ou3dfJT29Yh45qgrszaKiO2A3wUu75lk359hk7Q92PdnRPPx8FXAbcC3MtN+v5p0aHuw38+UDwNvAx4bMH2t7PcG4eGIPuN8F7t6/IjybxOfDfwrcO5oq7P2iYjHA18E/joz7+2d3GcW+/6QTNH29v0ZkpmPZuZuwFxgj4jYpaeI/X6GdGh7+/0MiIiDgdsy88rJivUZ97++3xuEh2MRMK81PBe4eUR1qUpm3jvxUVpmng+sGxFbjrhaa43mPr0vAp/OzC/1KWLfnyFTtb19f+Zl5t3AxcABPZPs9zNsUNvb72fM84A/iYgbgXOAP4yIT/WUWSv7vUF4OM4Djmy+UbkXcE9m3jLqStUgIp4cEdE834PSp+8Yba3WDk27fhy4NjM/NKCYfX8GdGl7+/7MiIg5EbFZ8/xxwAuAn/cUs9/PgC5tb7+fGZl5YmbOzcztgMOAb2fmK3uKrZX93l+N6CAiPgvsB2wZEYuA91Bu4iczTwXOBw4C5gO/AY4aTU3XPh3a/uXA6yNiCfAgcFj6X2KG5XnAq4CfNvfsAbwd2Abs+zOsS9vb92fGU4D/iIhZlJD1+cz8akQcC/b7Gdal7e33q1EN/d7/LCdJkqQqeWuEJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklSl/w+WCA/xogl94AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "arkansas_samples = arkansas_results.extract()\n",
    "\n",
    "# Make a new array with same dimensions as alpha\n",
    "p_predicted_AR = np.empty(arkansas_samples['alpha'].shape)\n",
    "\n",
    "# Generate one p sample for each alpha sample\n",
    "for i in range(arkansas_samples['alpha'].shape[0]):\n",
    "    p_predicted_AR[i] = sts.dirichlet(arkansas_samples['alpha'][i]).rvs()\n",
    "    \n",
    "plt.figure(figsize=(12, 6))\n",
    "for i in range(4):\n",
    "    plt.plot(sts.uniform.rvs(loc=i+1-0.1, scale=0.2, size=arkansas_samples['theta'].shape[0]), arkansas_samples['theta'][:,i], ',', alpha=0.5)\n",
    "plt.title('Distribution over each component of theta for probability of being chosen as candidate in Arkansas')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAF1CAYAAADiNYyJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAApNElEQVR4nO3de7gkVXnv8e/LHgZQhvsYZW5AxCjgJboFDBhJ1IioQY8mosQLxhCMxkv0KBijiZdgzkmCJ8cLoiLBGyYaDVEMMSoqCspgiAqIZxyRGQdkABnACzjwnj9W9UxNT/fetWe6p8dZ38/z7Gd3Va2qWrV6VdWvq6u7IzORJEmSarPTpCsgSZIkTYJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWD8BxFxJkR8RcjWtbSiLgjIqaa4Ysi4oWjWHazvM9ExPNGtbwdVUScExFvnnQ9ahARL4qIHzX9ft8O5Z8fERdvi7oNWf+c6rstRcQxEbF6C+c9ICIyIuYNmf7aiHjvoLLb8rgSEW+OiJsi4oYB07Zm+x8dEddsfQ3ntM4trq82FxF/GREfbB5vci6dqew2qNfIMsIc1zvjPj2G9Y00r0ySQbglIq6NiJ9FxO0RcWtEfDUiTomIDe2Umadk5ps6LutxM5XJzOsyc/fMvHsEdd9sR8/MJ2bmP27tsiXY+gNtROwM/D3wO02/v3mUyx+wvq0K0bPVd0eWmX+dmQNPcu3jyjhfqETEEuCVwCGZed9RLjszv5yZvzbKZWpyRnwu3aoLI10zwgzrf0JEfKnJIWsj4osR8btbujzNziC8uadk5gJgGfBW4DXA+0a9km31qu2Xle2zQ/oVYFfgyklXpKMtrm8UW3V8dR9gGXBzZt446YpI20JEPAP4Z+BcYDHlGPR64CljXm/dx5rM9K/5A64FHtc37nDgHuCwZvgc4M3N4/2ATwG3ArcAX6a8uPhAM8/PgDuAVwMHAAn8IXAd8KXWuHnN8i4CTge+DqwD/hXYp5l2DLB6UH2BY4G7gF806/vv1vJe2DzeCXgd8APgRsqOtmczrVeP5zV1uwn48xnaac9m/rXN8l7XLH+Xpi0Oa5Vd2LTDfZrhJwNXNOW+Cjykb3teA3wTuLPXLn3rfiDw2aa9rwF+vzXtScB/AbcBq4C/7Jv36GadtzbTn996Tt8BfBq4Hfga8KszbP+w5Qxsl2ba84GvAGc0860EfqMZv6p5Tp7XWsc5wJnNtt4OfBFY1pr+G8BllH5yGfAbrWkXAW9q1nc78B/Afq3pR7bq/9/AMV3mbfpGUvrYHcCjBrTNLsDbgDXN39uacQ8AftKa//MD5t1s+U37XAz8LfBj4PvAE/v64vuA64EfAm8GpoAHAT8H7m6WdWuXPtJa7sD6dmj3tzRt9zPg/kOOMacBVzXb835g1/Y+TtkHbqAcRwa2Z1/511L22WuBE7vsD2zc509ulns98MrW9L8EPthXtn2ceuGgNgYeCfyI1r4LPB24Yo7Hksc1bXhPs+xzBsw72/bvQuk31zV1OhPYbdDxtJn3VZRjzzrgo73npZn+6qaN1jTbnoOe36bsPs3zuqZ5jj/ZV99XUvb364GTZmuLZtr9KceAdc22frTjMfEc5nZs+2dK31tHOUcd2pp2HKXf3k7Z1141w3L+CLi6KXsV8PBm/KnA91rjn9aa5/nMvK8f2LTB7c32vp3hfXRo2Zm2k7I//IJyPr0D+Ldm/P7Ax5vn5vvAS2fY9nPYmBFmfM775gtKX/2fMyy7y3l8XqvO5zf9YgXwR33798eAD1KODy+kZJ1LKPvx9U2bzW/N83jgO02bvb1p31nzxS/D38QrsD39MSAIN+OvA17UPG538tMpB9edm79HAzFoWa1Oei5wb2C3AR33IsoB5rCmzMfZuKMfw5Ag3OrYH+ybflGro76g2RkOAnYH/gX4QF/d3tPU66GUIPqgIe10LiWkL2jm/S7wh820s4G3tMq+GPj35vHDm53kCEpYeV6zDbu0tucKYAnNCatvvfemnNBPAuY1y7uJjQexY4AHNzvlQygnv6c205ZSDorPap6rfYGHtZ7TWygHgnnAh4Dzhmz7TMuZqV2eD6xv6j5FCWzXUU5SuwC/0yx391adbgd+s5n+f4CLm2n7UE4Uz2nq+6xmeN/W8/49SpjbrRl+azNtEXAz5aS2E+XgdjOwsMO8B9Dqr0Pa543ApcB9KC+Cvgq8qcv8g6Y37fYLyol1CngRJWD09rNPAu+m9I37UF5E/nFr3ov71nEMQ/rIbPXp2O7XAYc203cecoz5NqWP70MJze2T5nrgb5rnfLdZ2rNX/u+b8o+hhPdf67A/9LbtI03bPZhykt/seDKgHS5i43FlUBtfxaYB5hO0QvYcjiXH0HfMG/BczrT9b6MEgX2a5f8bcPqgZTfPy9cp4WEfSog7pZl2LCU0HQrci/ICZaYg/GlKkN6bcox4TF9939iMPw74KbB3h7b4CPDnzXO5K3B0x2PiOXQ8trXOEwvY+ALsita064FHN4/3pgm3A5bxe5Tz2CMp4e7+NC/im2n7N9vxzOb5ul/Hff2S1nP9m5Tj47A+OrRsh+08h2afbIZ3Ai6nXJmdTzmHrgSeMGT7N8w/23PeN98Dm204cJbnZ7bzeK8Nvgi8s+kvD6Ps349t7d+/AJ7abN9uwCMoF0nmNcu6Gnh5U34/SmB+RrMdr2i2a9Z88cvwN/EKbE9/DA/Cl9JcIe3r5G+kHLiGXfkZFIQPGjCufYJ5a2v6IZRXplNsfRD+HPAnrWm/1uwIvU6fwOLW9K8DJwzYrilKSD6kNe6PgYuax48DVramfQV4bvP4XTQn8db0a9h4orgWeMEMz88zgS/3jXs38IYh5d8GnNE8Pg34xJBy5wDvbQ0fB3xnSNmBy+nQLs8H/l9r2oObNv+V1rib2TScn9eatjvlytsSShD7et/6L2HjlemLgNe1pv0JG1+MvIa+AxRwIc3V6Fnm7fWTmYLw94DjWsNPAK7tMv+g6U27rWgN36spc1/K24Z30nrRRAmnX2jNe/Gwuvb3kdnq07Hd3zjL+q6lCVitvva95vExlP29fSVypvY8hnIyundr+j8Bf9Fhf+ht2wNb0/8X8L7m8V+y5UH4NcCHmsf7UE7899uCfeYYugXhzbafEsB+QuvqJ+Udhu8PWnbzvPxBX1uc2Tw+myZAN8P3Z0gQBu5HuYo9KOgcQ7nK3e7fN1LCx2xtcS5wFq1jdDN+xmMiczi2DajvXs127tkMX9fUaY9Z5rsQeFnHdVwBHN/qS8P29aUDnusPD+qjs5XtsJ3nsGkQPgK4rm+e04D3D1nehvlnes4HzHdUU49dBy23KdPlPD6Pcp64G1jQKns6zTsrlP37S7M8Ny+nOdcBzwUubU0LypXuWfNFl34w6T/vEe5mEeVVdb//TXkV9B8RsTIiTu2wrFVzmP4Dyquv/TrVcmb7N8trL3seJUz0tD+Z/VNK+Oq3H+VVcf+yFjWPPw/sFhFHRMQyyivRTzTTlgGvbD6IeGtE3ErZYfdvLWum9lkGHNE3/4mUAyXNOr/QfMBgHXAKG9tuCSVUDNNl22dazmztAuWKXM/PADKzf1x7vRvaIjPvoPTB/dn8uRy0rmHbswz4vb42PJpyAp9t3i4G9bP9h5TtakN9MvOnzcPdKduyM3B9a1veTbl6OtAsfWQ2Xdp9tv27v0x/+6zNzJ/PsM7+8j/OzJ8Mmt5xW2eqy5b6IPCUiNgd+H1KULt+QLku+8xshm3/QkqQurzVN/69GT/MsH6/P5u200zP8RLglsz88ZDpN2fm+gHrma0tXk0JH1+PiCsj4gXN+BmPibNs1yYiYioi3hoR34uI2ygvDmBjn3k6JUj/oPkA16OGbOPQY21EPDcirmjV9TA27ZPD9vX9GfxcDzJj2Q7b2W8ZsH9fG7+WTc+dMxn2nG9Wrvl/vwHTerqcx3vlbsnM2/vKDj1WRcQDIuJTEXFD0y5/zcY22WQfyJJ22/N3rdd2ySA8i4h4JKXzbPbJ6My8PTNfmZkHUW5m/7OIeGxv8pBFDhvfs6T1eCnlVdVNlKsb92rVa4pND+qzLXcNZYduL3s9m4azLm5q6tS/rB8CZOY9lKsyzwKeDXyqtTOuotw2sVfr716Z+ZGO27EK+GLf/Ltn5oua6R+mvBW6JDP3pNy2Eq15f3WO2zqsDoOWM2O7bKENfaEJFfuw8V7RZX1lu65rFeWKcLsN752Zb+0w72x9jAF1W9qM66LL8ttWUa6i7dfalj0y89AZljdTH5lNl3bvsg39+3i7ffrnn609946Iew+Z3mVbZ6pLF5ttb2b+kHKl/GmUq+gfGDLvKPaZYdt/E+WF5aGtvrFnZs7lRV3P9ZQPLvUsGVaQ0if3iYi95riO2Y6rN2TmH2Xm/pSrsu+MiPsz+zFxLp4NHE95V29PyhVGaPpMZl6WmcdTXmh+knKcH2TgMbK5MPIe4CWU24n2otwm1GX/u57Bz/WWlJ1xO9m8T6+ivJPQbuMFmXlch3rPxTXNup4+Q5mu5/E1lH64oK/sTMeqd1HuAT44M/eghP1em1zPpuejYNP9YFT5YiIMwkNExB4R8WTgPMpbKt8aUObJEXH/plPcRnkrovf1LT+i3C8zV38QEYdExL0ot158LMtXwnwX2DUintR8rdPrKPc39fwIOGCGT6p/BHhFRBzYhKq/pnzgYv2Q8gM1dfkn4C0RsaA5uP0Z5SpQz4cpb9md2DzueQ9wSnOlKiLi3s32tHfWmXwKeEBEPCcidm7+HhkRD2qmL6C8Cv55RBxOOeD1fAh4XET8fkTMi4h9I+Jhc9n2mZbTsV3m6riIODoi5lM+wPa1zFwFXEBph2c3dXgm5TaaT3VYZu9q3ROaKyO7Rvl+08WzzlnuMbuHmfv1R4DXRcTCiNiPcl9d1zbosvwNmquM/wH8XbO/7hQRvxoRj2mK/AhY3LRfz0x9ZDZb0+5tL46IxRGxD+Vk89EZynZpz7+KiPkR8WjKh1H/uRnfZVv/IiLuFRGHUu4znakugwxqYyhv5b+acgvQJzabi87Hki422/7mBfl7gDMi4j4AEbEoIp4wx2XT1PGkiHhQc1x+/bCCTZ/8DCWo7t0co35zthXM1hYR8XutffTHlBBzN7MfE+diAeWF5c2Uiy5/3ZvQtO+JEbFnZv6Cjee7Qd4LvCoiHtEc5+/fbM+9m3qvbZZ5EuWK8Kwy8wfAcjY+10cz5JsUOpQdup2N/nP314HbIuI1EbFbc9w8LMpFspFprrL+GWWfPKl1TDs6Is5qinU6jzfnia8CpzfH+IdQPqj/oRmqsIDyvN4REQ+k3KPd82ng0Ij4H1G+YeKlbPquw0jyxaQYhDf3bxFxO+WV2Z9Tbrg/aUjZg4H/pHy69BLgnZl5UTPtdMoJ7NaIeNUc1v8Byj1GN1Bucn8pQGauo9yv+V7Kq7qfUO7R6emd/G6OiG8MWO7ZzbK/RPnU68+BP51Dvdr+tFn/SsqV8g83y6ep69ea6ftTTgq98cspH4R4O+VgvoJyX1gnzZXl3wFOoLwCvYGNHyyC0j5vbJ6/19O6YpGZ11He1nsl5RaDKygfCpyTWZYzY7tsgQ8Db2jW8wjKCwuyfJ/tk5s63EwJHE/OzJs61H8V5WrIayknpFXA/6TDsaB5q/ItwFeafn3kgGJvppyEvgl8C/hGM25WHZff77mUt5R738LwMTa+tfh5ylef3RARvbYZ2kc61G+L273PhykBfmXzN1P7zNaeN1C2ew3lJHdKZn6nmdZlW79I2Q8/B/xtZv7HHLdlUBtDCb/LKPcY/mTgnMXW7jMzbf9rKNt2aZS3ev+Tcu/inGTmZ4B/AL7QLO+SZtKdQ2Z5DuXq7nco94O+vOOqZmqLRwJfi4g7KFf5X5aZ3+9wTJyLcylvaf+Qsj9dOmC7rm3a8hTgDwYtJDP/mbIff5jyIbVPUr796Crg7yjt9yPKi6SvzKF+z6bcr3sL5bh47haWnW073wcc0hyDPtm8SHkK5Ta/71Ou3r+XcjV5pDLzY5SLSC+gPJ8/ouzv/9oUmct5/FmUq91rKPvjGzLzszOs/lWUdrud8iJyw4vi5hj3e5SvlL2Zkn3az90o88U21/s0pqTtSEScQ/kwz+smXReNTkRcS/mAyX9Oui7jFhHfo3yDxw61rc3V1m9Tvu3ml+KKl6ThvCIsSRqpiHg65W3wz0+6LqMQEU9r3mbfm3LF9d8MwdKOwSAsSRqZiLiI8sGbFzf36u4I/phyK9H3KPfGbsmH0SRth7w1QpIkSVXyirAkSZKqZBCWJElSleZNasX77bdfHnDAAZNavSRJkipx+eWX35SZm/265MSC8AEHHMDy5csntXpJkiRVIiIG/ix3p1sjIuLYiLgmIlZExKlDyhwT5TfEr4yIL25NZSVJkqRxm/WKcERMAe8AHk/5JbPLIuL85ldiemX2At4JHJuZ1/V+0lKSJEnaXnW5Inw4sCIzV2bmXcB5lJ9obXs28C/Nz8+SmTeOtpqSJEnSaHUJwouAVa3h1c24tgcAe0fERRFxeUQ8d1QVlCRJksahy4flYsC4/l/hmAc8AngssBtwSURcmpnf3WRBEScDJwMsXbp07rWVJEmSRqTLFeHVwJLW8GJgzYAy/56ZP8nMm4AvAQ/tX1BmnpWZ05k5vXDhZt9gIUmSJG0zXYLwZcDBEXFgRMwHTgDO7yvzr8CjI2JeRNwLOAK4erRVlSRJkkZn1lsjMnN9RLwEuBCYAs7OzCsj4pRm+pmZeXVE/DvwTeAe4L2Z+e1xVlySJEnaGpHZf7vvtjE9PZ3+oIYkSZLGLSIuz8zp/vGdflBDkiRJ2tEYhCVJklQlg7Ckwb5w+qRrIEnSWBmEpVp84fSNfz3vP27zMj2/ddqm8/TKDgrI7z9u02W3y7TH9T82bEuSJsgPy0m16A+9ALdeBz9fB7vuCa/4dgmm134ZTrqgPL7iQ6XcK769cf4DHl3+X/vlMv8rvg1nHAZ7LS3DDzuxTAO44Vtl2Q87sSyrV6Zd9ooPlWVIkjQmwz4sZxCWanHGYbBuFSw7Cq67BPKeMn7PJWU8wC57wPqfwz3rN04HiNabR73xU/Ph7rvK/D9fB3fetnF5P1+3cTk7zYN5u5ZpvTI9y44qwfq3Thv99kqS1BgWhLv8xLKkHUEv7P7gK4PHw+ZBtacdinvuvmvz+QcN333XxrL9fvCV8mcQliRNgPcIS5IkqUoGYUmSJFXJICxpsqbmT7oGkqRKGYQlbSoGHBbmElZ32WNu6+t9kE6SpG3MICzVYtlR5Rsd+kPtnkvK/974+buXx8uO2hhqd/+VjY932aNMg42huT3uztvK/P3r6g3vuWRj2WVHwWl9H66TJGkb8VsjpFoc8Ojy/b4POxEufScc+Sfl2xrOOKwE0vZ3AJ90wcb52t8/fNIFG4cfc+rG8b3vDe4F3F653ncTQ/m+4L9YWx5/4fRN1yFJ0gQYhKVa/NZpQPNrcaet2hhQH3Zi39eX9X2VWTuwDguwX2ivo9H74Y3eL9S1Q3GvXP+wJEnbkD+oIUmSpB3asB/U8B5hSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklSlTkE4Io6NiGsiYkVEnDpg+jERsS4irmj+Xj/6qkqSJEmjM2+2AhExBbwDeDywGrgsIs7PzKv6in45M588hjpKkiRJI9flivDhwIrMXJmZdwHnAcePt1qSJEnSeHUJwouAVa3h1c24fo+KiP+OiM9ExKEjqZ0kSZI0JrPeGgHEgHHZN/wNYFlm3hERxwGfBA7ebEERJwMnAyxdunRuNZUkSZJGqMsV4dXAktbwYmBNu0Bm3paZdzSPLwB2joj9+heUmWdl5nRmTi9cuHArqi1JkiRtnS5B+DLg4Ig4MCLmAycA57cLRMR9IyKax4c3y7151JWVJEmSRmXWWyMyc31EvAS4EJgCzs7MKyPilGb6mcAzgBdFxHrgZ8AJmdl/+4QkSZK03YhJ5dXp6elcvnz5RNYtSZKkekTE5Zk53T/eX5aTJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCMzjjs9/d8L/3WJIkSTuGKoNwf7DtEnTPvnglR731c0OX1/v/zHdfMnC5/cvvlRs2XZIkSeMVmTmRFU9PT+fy5csnsu4H/PkFLFywC1859bGbBFeAq9asY4/ddua2n/2CQ/bfk6vWrOMFRx8EwMcuX8Xa2+/k15fuzeof/xRgw/CRB+3Lxy5fxeK978WRB+27oXxvOat//FOe8YglvOLxD+CZ776Eq9as487192xYVm+aJEmSRisiLs/M6c3G1xiEDzrt09yTsFPAPc3mz58K7rp707aYPxUAm40HWLDLFLffefcm43YqxTcsc5CdAh55wD587fu3sFPAvJ02rnfRXrvylVMfu4VbJUmSpEGGBeEqb43oBdV2YB0Udmca3x+CF+wyxT05cwjurfNr379lw+Nhy5ckSdJ4VRmEu5pLSO0PxlviGY9YstXLkCRJUjcG4e3I2RevnHQVJEmSqmEQ3o4csv+ek66CJElSNQzC25Het01IkiRp/KoMwkccuA/XvvVJzJ8Krn3rk1iwyxRHHLgPRxy4Dy977MEs2GWKRXvtyk7BhvG9x/OnggW7TG0YXrTXrhxx4D5A+cBc72+nKMO98b31tpe9aK9dN/z1vnFCkiRJ28a8SVdgUs747Hf57luO44zPfpc719+zydXYQ/bfc8Pw2Rev5AVHH8RH//hRG+a7dOXNG4af+e5LOPKgfTeZ3v99wL15ess88qB9NyyjV77/+4wlSZI0XlV+j/Ag7QDb/xjgFY9/wMCQO2j+9ny9eWda72xlJEmStOX8QQ1JkiRVyR/UkCRJkloMwpIkSaqSQViSJElVMghLkiSpSp2CcEQcGxHXRMSKiDh1hnKPjIi7I+IZo6uiJEmSNHqzBuGImALeATwROAR4VkQcMqTc3wAXjrqSkiRJ0qh1uSJ8OLAiM1dm5l3AecDxA8r9KfBx4MYR1k+SJEkaiy5BeBGwqjW8uhm3QUQsAp4GnDm6qkmSJEnj0yUIx4Bx/b/C8TbgNZl594wLijg5IpZHxPK1a9d2rKIkSZI0evM6lFkNLGkNLwbW9JWZBs6LCID9gOMiYn1mfrJdKDPPAs6C8styW1hnSZIkaat1CcKXAQdHxIHAD4ETgGe3C2Tmgb3HEXEO8Kn+ECxJkiRtT2YNwpm5PiJeQvk2iCng7My8MiJOaaZ7X7AkSZJ+6XS5IkxmXgBc0DduYADOzOdvfbUkSZKk8fKX5SRJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqtQpCEfEsRFxTUSsiIhTB0w/PiK+GRFXRMTyiDh69FWVJEmSRmfebAUiYgp4B/B4YDVwWUScn5lXtYp9Djg/MzMiHgL8E/DAcVRYkiRJGoUuV4QPB1Zk5srMvAs4Dzi+XSAz78jMbAbvDSSSJEnSdqxLEF4ErGoNr27GbSIinhYR3wE+Dbxg0IIi4uTm1onla9eu3ZL6SpIkSSPRJQjHgHGbXfHNzE9k5gOBpwJvGrSgzDwrM6czc3rhwoVzqqgkSZI0Sl2C8GpgSWt4MbBmWOHM/BLwqxGx31bWTZIkSRqbLkH4MuDgiDgwIuYDJwDntwtExP0jIprHDwfmAzePurKSJEnSqMz6rRGZuT4iXgJcCEwBZ2fmlRFxSjP9TODpwHMj4hfAz4Bntj48J0mSJG13YlJ5dXp6OpcvXz6RdUuSJKkeEXF5Zk73j/eX5SRJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqtQpCEfEsRFxTUSsiIhTB0w/MSK+2fx9NSIeOvqqSpIkSaMzaxCOiCngHcATgUOAZ0XEIX3Fvg88JjMfArwJOGvUFZUkSZJGqcsV4cOBFZm5MjPvAs4Djm8XyMyvZuaPm8FLgcWjraYkSZI0Wl2C8CJgVWt4dTNumD8EPjNoQkScHBHLI2L52rVru9dSkiRJGrEuQTgGjMuBBSN+ixKEXzNoemaelZnTmTm9cOHC7rWUJEmSRmxehzKrgSWt4cXAmv5CEfEQ4L3AEzPz5tFUT5IkSRqPLleELwMOjogDI2I+cAJwfrtARCwF/gV4TmZ+d/TVlCRJkkZr1ivCmbk+Il4CXAhMAWdn5pURcUoz/Uzg9cC+wDsjAmB9Zk6Pr9qSJEnS1onMgbf7jt309HQuX758IuuWJElSPSLi8kEXaf1lOUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKnUKwhFxbERcExErIuLUAdMfGBGXRMSdEfGq0VdTkiRJGq15sxWIiCngHcDjgdXAZRFxfmZe1Sp2C/BS4KnjqKQkSZI0al2uCB8OrMjMlZl5F3AecHy7QGbemJmXAb8YQx0lSZKkkesShBcBq1rDq5txcxYRJ0fE8ohYvnbt2i1ZhCRJkjQSXYJwDBiXW7KyzDwrM6czc3rhwoVbsghJkiRpJLoE4dXAktbwYmDNeKojSZIkbRtdgvBlwMERcWBEzAdOAM4fb7UkSZKk8Zr1WyMyc31EvAS4EJgCzs7MKyPilGb6mRFxX2A5sAdwT0S8HDgkM28bX9UlSZKkLTdrEAbIzAuAC/rGndl6fAPllglJkiTpl4K/LCdJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVCWDsCRJkqpkEJYkSVKVDMKSJEmqkkFYkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVMghLkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVUyCEuSJKlKBmFJkiRVySAsSZKkKhmEJUmSVKVOQTgijo2IayJiRUScOmB6RMQ/NNO/GREPH31VJUmSpNGZNQhHxBTwDuCJwCHAsyLikL5iTwQObv5OBt414npKkiRJIzWvQ5nDgRWZuRIgIs4DjgeuapU5Hjg3MxO4NCL2ioj7Zeb1I6+xJEkdPfgfHzzpKgz1red9a9JV0A7q6gc+aNJVGOpB37l60lXYRJcgvAhY1RpeDRzRocwiYJMgHBEnU64Ys3Tp0rnWVZKkOTFsqkbbW9jcnnW5RzgGjMstKENmnpWZ05k5vXDhwi71kyRJksaiSxBeDSxpDS8G1mxBGUmSJGm70SUIXwYcHBEHRsR84ATg/L4y5wPPbb494khgnfcHS5IkaXs26z3Cmbk+Il4CXAhMAWdn5pURcUoz/UzgAuA4YAXwU+Ck8VVZkiRJ2npdPixHZl5ACbvtcWe2Hifw4tFWTZIkSRoff1lOkiRJVTIIS5IkqUoGYUmSJFXJICxJkqQqGYQlSZJUJYOwJEmSqmQQliRJUpUMwpIkSaqSQViSJElVivKjcBNYccRa4AcTWflo7QfcNOlKVMq2nxzbfnJs+8mx7SfHtp+cHaXtl2Xmwv6REwvCO4qIWJ6Z05OuR41s+8mx7SfHtp8c235ybPvJ2dHb3lsjJEmSVCWDsCRJkqpkEN56Z026AhWz7SfHtp8c235ybPvJse0nZ4due+8RliRJUpW8IixJkqQqGYQ7iIizI+LGiPj2kOkREf8QESsi4psR8fBtXccdVYe2PyYi1kXEFc3f67d1HXdUEbEkIr4QEVdHxJUR8bIBZez7Y9Cx7e37YxARu0bE1yPiv5u2/6sBZez3Y9Cx7e33YxQRUxHxXxHxqQHTdsh+P2/SFfglcQ7wduDcIdOfCBzc/B0BvKv5r613DjO3PcCXM/PJ26Y6VVkPvDIzvxERC4DLI+KzmXlVq4x9fzy6tD3Y98fhTuC3M/OOiNgZuDgiPpOZl7bK2O/Ho0vbg/1+nF4GXA3sMWDaDtnvvSLcQWZ+CbhlhiLHA+dmcSmwV0Tcb9vUbsfWoe01Jpl5fWZ+o3l8O+XguKivmH1/DDq2vcag6ct3NIM7N3/9H6ax349Bx7bXmETEYuBJwHuHFNkh+71BeDQWAataw6vxpLUtPap5K+0zEXHopCuzI4qIA4BfB77WN8m+P2YztD3Y98eieXv4CuBG4LOZab/fRjq0Pdjvx+VtwKuBe4ZM3yH7vUF4NGLAOF/FbhvfoPxs4kOB/wt8crLV2fFExO7Ax4GXZ+Zt/ZMHzGLfH5FZ2t6+PyaZeXdmPgxYDBweEYf1FbHfj0mHtrffj0FEPBm4MTMvn6nYgHG/9P3eIDwaq4ElreHFwJoJ1aUqmXlb7620zLwA2Dki9ptwtXYYzX16Hwc+lJn/MqCIfX9MZmt7+/74ZeatwEXAsX2T7PdjNqzt7fdjcxTwuxFxLXAe8NsR8cG+MjtkvzcIj8b5wHObT1QeCazLzOsnXakaRMR9IyKax4dT+vTNk63VjqFp1/cBV2fm3w8pZt8fgy5tb98fj4hYGBF7NY93Ax4HfKevmP1+DLq0vf1+PDLztMxcnJkHACcAn8/MP+grtkP2e781ooOI+AhwDLBfRKwG3kC5iZ/MPBO4ADgOWAH8FDhpMjXd8XRo+2cAL4qI9cDPgBPSX4kZlaOA5wDfau7ZA3gtsBTs+2PWpe3t++NxP+AfI2KKErL+KTM/FRGngP1+zLq0vf1+G6qh3/vLcpIkSaqSt0ZIkiSpSgZhSZIkVckgLEmSpCoZhCVJklQlg7AkSZKqZBCWJElSlQzCkiRJqpJBWJIkSVX6/wZ1cLqTGFHYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "colorado_samples = colorado_results.extract()\n",
    "\n",
    "# Make a new array with same dimensions as alpha\n",
    "p_predicted_CO = np.empty(colorado_samples['alpha'].shape)\n",
    "\n",
    "# Generate one p sample for each alpha sample\n",
    "for i in range(colorado_samples['alpha'].shape[0]):\n",
    "    p_predicted_CO[i] = sts.dirichlet(colorado_samples['alpha'][i]).rvs()\n",
    "    \n",
    "plt.figure(figsize=(12, 6))\n",
    "for i in range(4):\n",
    "    plt.plot(sts.uniform.rvs(loc=i+1-0.1, scale=0.2, size=colorado_samples['theta'].shape[0]), colorado_samples['theta'][:,i], ',', alpha=0.5)\n",
    "plt.title('Distribution over each component of theta for probability of being chosen as candidate in Colorado')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding up what we predict will be victories for each state, we obtain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The most likely scenario from our results implies Republicans with 18 electoral college votes and\n",
      "Democrats with 20 electoral college votes.\n"
     ]
    }
   ],
   "source": [
    "republican_votes = electoral_votes['Alabama'] + electoral_votes['Arkansas'] + electoral_votes['Alaska']\n",
    "democrat_votes = electoral_votes['Arizona'] + electoral_votes['Colorado']\n",
    "\n",
    "print(f'The most likely scenario from our results implies Republicans with {republican_votes} electoral college votes and\\nDemocrats with {democrat_votes} electoral college votes.')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}