{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "As a scientist, watching the Brexit vote was a little bit painful. Though probably not for the reason you're thinking. No, it wasn't the politics that bothered me, but the method for making such an incredibly important decision. Let me explain...\n", "\n", "Scientists are a bit obsessed with the concept of error. In the context of collecting data and anaylzing it, this takes the form of our \"confidence\" in the results. If all the data say the same thing, then we are usually pretty confident in the overall message. If the data is more complicated than this (and it always is), then we need to define *how confident* we are in our conclusions.\n", "\n", "Which brings me to this gigantic nation-wide referendum vote. I couldn't help but notice that the cutoff for winning / losing the vote was set at 50%. To me, this sounds crazy. If I simply flipped a coin at 50% and tallied the results each time, I'd get some difference between # heads and # tails that would vary around 50%. In the context of voting, it means that a yes/no split that's really close to 50% might actually be too close to call.\n", "\n", "In science, saying that a number is *different* from some other number requires that the difference falls outside of a certain region of uncertainty. It's a way of saying \"yeah, I know that random fluctuations cause strange looking data sometimes, but my difference is *so far* from those fluctuations that I think there's something real going on.\"\n", "\n", "But this is all a little abstract, so let's try it out on some voting data...\n", "\n", "# Simulating a national referendum vote\n", "For a referendum vote to go through, it seems reasonable to say \"the people need to vote in numbers that are significantly different from random chance. To ask what \"random chance\" looks like, we can use computer simulations.\n", "\n", "We'll take on the task of assessing what national votes might look like if they happened completely randomly. Then, we can compare the actual results to our simulation in order to decide if we've got a \"real\" result or not." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# First, import a bunch of stuff that we'll use later\n", "import numpy as np\n", "from scipy import stats\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "from tqdm import tqdm\n", "import seaborn as sns\n", "sns.set_style('white')\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Initializing the simulation\n", "# We'll simulate ten thousand votes.\n", "# On each iteration, generate random votes and calculate the results \n", "n_votes = 10000\n", "\n", "# This is the actual difference in percentage points that we had between the sides\n", "actual_diff = 51.9 - 48.1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we'll create a completely random vote. Each person randomly chooses between the two options: *yes* and *no*. Then, we compare the difference in percentage points between the two." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Create 10,000 citizens and assign each a random vote\n", "total_population = int(1e4)\n", "diff = np.zeros(n_votes)\n", "for ii in tqdm(range(n_votes)):\n", " votes = np.random.rand(total_population)\n", " yes = np.sum(votes < .5) / float(total_population)\n", " no = np.sum(votes > .5) / float(total_population)\n", " \n", " # This is the difference in percentage points\n", " diff[ii] = (yes - no) * 100" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`diff` is a list of numbers representing the lead that \"yes\" has over \"no\". Remember, we've randomly chosen these values, so they are the results you'd get if every single person in the country voted completely randomly.\n", "\n", "How can we summarize the \"limits of uncertainty\" that `diff` defines? We can use percentiles to get an idea for the variability of this number. We'll take the 1st and the 99th percentile of our simulated differences as a proxy for the limits of what we'd expect if there were *no* true opinion in the population" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Here we calculate 98% confidence interval on the difference\n", "clo, chi = np.percentile(diff, [1, 99])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we'll make a plot with 3 things:\n", "\n", "1. The distribution of all our simulated differences\n", "2. A vertical black line for each limit of the confidence interval.\n", "3. A vertical red line representing the actual difference between yes/no that was reported" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAFJCAYAAACl2EgGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xuc1VW9+P/XMCg6gROSYIo1ytGlx1t5KU3QUiw1jEAz\nzLyhlGaePAfTxMq72E80QI+pJxIzBPMCpWSpmAqEX8kyry1TnBhFES1uggxz+f3x+exhs5lhLjCf\nPcO8no/HPGb22ms+6/3Zl9nvWWt91iqpr69HkiRJ2ehW7AAkSZK6EpMvSZKkDJl8SZIkZcjkS5Ik\nKUMmX5IkSRky+ZIkScpQ92IHIOULIVwGXFZQXA+sBhYBfwRujDHGgt/7JPAGMCPGOLwN7R4E9I4x\nPtqKGL8aY/ztprbdgvaOBv4VY3w2vX0EyeMwPsb4P5u7vc0thFAKXAd8E/goEGOMn9pI/fXOtw3t\n9QDOizHe2MbfnwycBnwqxvh8ez+/xRZCmAF8BaiIMS4sdjztqfC5TcvqgOdijAfk1dsd2D/GeF9R\nAtUWz54vdUT1wAzg8vTrKuB2oAo4G/hLCOG4gt9Zmtad1trG0mM9DezVwl95Im3r761tq7VCCOcC\nfwB2ziuuTNv/fXu3v5mcDYwG/g38FJjcVMUmzre1ngJ+uAm/X59+dRVd6XwbO9fLgVtzN0II+wEv\nAJ/LLix1NfZ8qaOaEWP8ZWFhCOEYksRsWgjhUzHGBQAxxmXAlW1sawegpKWVY4xPAk+2sa3W6kvB\nh0WM8Z+0/VyL4dMk53BejPGPzdTd4HzboO8m/r66kBhj4XupN7B1MWJR12HPlzqVGOPvgR8BPdPv\nm0MJrUi+MtZR42qNbdLv77eg7pZwvurcfA2q3dnzpc7oZuAK4IQQwlkxxrrG5uWkc41+CAwHBgBr\ngPnA/xdjfDytcwdwOklvy/gQwk+BXdOvPwLfAY4AhpIMbZ4AfJG8OV/5gYUQhqWx7Q4sBO4AxsUY\na/LqbDDHJC0/Pa1/QYxxYgjhj2nb9cCMEEJ9jLG0qTlf6TyVy4DBJP+9VwH3A9fEGJfn1ZtMMu9l\ne2As8FWSuVgvAdfGGB9oyZOQzs26CPgMsBXwCvBz4NYYY33ec0J6Ds+FEOqBL8QYn2rkeI2eb3pf\nCXAOMArYE6gGniF5Lh9L6+TaqwdK0sd5coxxZHr/3sAP0jb6AR+SDC/d2NJzTo/zTeCXwNUxxh8X\n3Lct8C7J8zuoid/PxXkVyfM0kmRO4zkxxvtDCH3SOL8MfDL9tTeAKen51qbHyb0OziD5W34Byevu\nPZLh9x/HGFfntduNZPj3LOATwD9IhtyaOs+NPr959SqBV4H/AcYBh5E8tvcD/52e40+BL5G8Bx8h\neY1vNBlvyfs3rTeZ5PXcD7gRGALUkQw/XxpjfLmZdhrej3nzOeuBC0II3yN9vYYQBpC8Xz4D7Ai8\nDfwOuDLGuHhjbUiF7PlSp5N+oPwF+AjQ5MRtkiTtMpIel5uAe0j+cP4hhHB4Wmc6yTAmJHOoLidJ\nsnIuAw4EJgJ/TtuFxofGPgf8GngNuAWoBa4l+aBuqfzj3sG64c1pbPyD8rPAX4GvA38iOd/FwPeB\neSGEjxa0UQ88SvKBeA/wK+A/gV+HEAY3F2QI4XySuVkHAg8Ak4DtgP8lSRJg3Ty8v6W3byVJTCub\nOGyj55smXvekx+6VtjUdOIjkuTynoL3lJB/+l5E+tyGEz5B8cB9H8jyPS78fDNzbyBzCjXkAWAmc\n3Mh9w4Ay4M4WHOdbwIkkr5V5wNMhhO1Iksr/IkmGx5M8njsC15B8+Bc6Pz3GC8AEkkRuNMk8yXx3\nAj8B1pI8F1XAfcAhhQds4fObUw/sBsxNb99CkpiMStucC+wC3EaS8J3SSGyNacn7N9d+PfAw8Pk0\n1kdJkrA5IYR9W9BWzhMkcxJLSOaBXg5UhhA+BjwOHEuS8N4AvAicC/wxTRSlFrPnS53VW+n3jzd2\nZwihF8kf/ydjjEfmlU8i+XA7D3gqvVqxN0nvz+9jjBPTerlf6QnsF2NckneMpmLaAfivGOP/pvXG\nADOBr4cQJsUYZ7XgvBqGPGKMvwwh7AocDkwr7GXLi6cbcBdJ78Rx+VdshhDGAhcD16ePR347NcB/\nxhg/TOs+TvLBOhJ4rKkA05huIEmivpDOQcv1+jyYnu/MGOMU4Mq0/n4kPSbPN3XcjZzvN0mSlIeB\nE3O9OSGECpIP9gkhhN/HGCvT9s4EymOMV+Ud/gqgFDg0xvhq3rmcSJIwf4OkF6NZMcZVIYQHgFND\nCAfHGOfn3X0KSa/cvS041A4kV9S9lBfPxUAFcHaM8Y688itIkvpvkPRG5dsfGBhjfCatey1JknNS\nCOHbabxfSGN7mKTHdm1a91yShCq/J6s1z2/OriQ9saPzYniLpKf41zHGk9PybiQXqnw1hLBN7rVX\nqKXv37xfKSHpyd0vxvjvtO4wkt63CcCRtECM8cn0/X0G8HTuNRRCOA/oD5yZPxc1hHATSe/4F0ke\nW6lF7PlSZ7Um/b5dE/d3I/mDvEsIoV+uMF2+YADJh1hLzM1PvJrxOsl//bm21gCXpnGc0sJjtMXn\ngP8A7m5kqYzLSD4ETwkhbJVXXg/cVPDhl0s+Kppp75skicwVuQ9maOiR/C+S8z2rtSexEWeQxPud\n/GG0NNm6hiTpPK2ZY9wInJKfeKVyH+CtnaR/J8l5NryO0t6Ro4EH0wtAmvNafuKV+j3J8Op6vaUx\nxreABU3E+WQu8UrrLifp/exOkjBA0ktXD/wwl3ildX/GhlfttvX5HZ9XdxmQG+77aV55HZBbQuST\nNK2179964Kpc4pXWnQ7MAY4IITT6T1or5OI5KE0gc8YAH48xmnipVez5UmfVK/2+srE7Y4zLQgj3\nkAzDLQwhzCX5z/ShGOMrrWjnjearNHg6fy5M6lmS+Sf7t+I4rfUpkg+f2YV3xBirQwjzSeas7Uky\nNJXzj4K6y9L/+ns0017uXBpr7+UQwlI27/nuD7yVnwjkmVMQU6NySWn6Qb4/yQf4nsDAtEqrho1i\njH8MIVSR9C79T/q8j0iPc1cLD7PBayvG+DfgbyGEj4QQPk2SVO9BMjy6exNxFiaUALnkL/dc7kcy\nDP63Rur+Ccjvzm3L87s2xlhVUPZB+r3wPHMJf5Ovsza+fzeYR0jSS3ZYGu/bTbXXAvcBPwa+C4wI\nIfwhjed3McZ3N+G46qLs+VJnVZF+X7CROqcCFwKRZJL1dcBLIYRnQggtTQ5WN1+lwQaTbtOJ9h+S\nDF+2l1zvX1O9LYvS72UF5WsKK6aau9qrJe0VtrUptmumLZprL4SwS7qY6CKSD82JJBcm/Dmt0pYr\n3H5FMhfrC+ntb5LMT2ppL8gGr60QQo8Qwo0kr6UnSeYvfTO93VQPbGPPY+6fgNx59QZWpz1Phf5V\ncLstz++qJuo2FV9LtPb9+1YjZe+k38vbGAMAMca3SeYYTmJdj+evgHdCCLcW9CpLzTL5UqeTztHa\nm2SCdZNXMsUYa2OMP40x7kcyxHE2ySTig4AH22GS7EcLC9IJ1GVs+AHX2HuvrQnLCpIPhKYWJu2d\nfm/JUg8tbY9m2ttcbeXa29Rz+x3J1YNXk/Qi9Ywx7s2mLVeSG3o8KYTwCZLJ4NPyr2xtgxtJrlp8\niGTyeJ8Y4ydijN+k6WSoJf4NlDXxmi/8xyDr57dRbXj/btvIYXLvyfc2Qzz/jDGOIhn6PYRkSH8R\nydy0zrTunjoAhx3VGX2b5LV7TyPDfEDDZOxRwJ9ijDNjjG+SXE13RwjhMZLeil1JJjFvrtW9D26k\nLLdK9p/zyqpJrtQs9B+NxNKS2J5Lvw8k6dFpkF4pOJBkeLaxYbu2eI7kqr6BwHoT6EMI/0FyEUSz\n2zQ1obHzfQ74fAjhPxtZNuCI9Hv+3Kn1jhGSFcv3Bu6NMRZuXfWf6fdW93zFGF8NITxDclXd82m7\nLR1ybMrJwOIY44j8whDCNmx8jlRzniV5fR7CuqsScwpft+35/LZIK9+/OQeTLGOR73MkF5a0Zquq\nDV6DIYTjgWOAi2OMK0munJ0fkqVqFgKNLisiNcWeL3UqIYQjSXorltP4Zfc5q0mu8rsyhNCwWnX6\n804kQyG5IYncBORNXdV63/TquVxbvUh6WupYf+mBvwO7hhD2yqv7SZJhlkItiW0OyYfQ8BDCsQX3\nXUlymf89+ROtN9GvSD7QxqRXxgEQQihj3ZVzLVlqoTGNne9kkuRoQtpGrr1dSebhVLP+tlJrSSbh\n5+TmGK03WT2EsD3JVaAU1G+NX5K8nr5PMoF+fjP1m/MhsE3+0iDpBO+JpD07IYS2/NOcez6uCyE0\n9HSFEEaQLCeRrz2f35ZqzfsXktfHFel7Llf3RJLewxkxxvzlY5rT2GtwT5JlJc4pqJt7fCpbcXzJ\nni91SCXAsLw//N1I5qEcQPIf5ipgRCMTfBvEGBeHZMHU/wFeDCHMJEmCjiGZXHxl+h8srJsr8p2Q\nLHA5oYUxFnodmJJe4r4EOJ5kbtrYGGN+z9f/kaxb9GQI4W6SFeBPIullOJz1vZW29aMQwgE0stZX\nTBY0PZ3kSrkHQwgPprF8jqSn4yU2XJ6gzWKMb4QQRpNc3faXdC7VSpI1kHYFpsYY727j4Tc43xjj\nXSGEr5Astvl8COFhkqGyoSQXXnw3xvhGwTH+I4RwF0lPyK9IJl4fHkJ4iqTn52Mky4v0IJkY3qeN\n8U4lGSr8BEkiuKl+RbJG15/Tx7U7yVpse5As3roDSaytWtQzxvhMCGEcyRyq50IID6UxDyVJ3Afk\n1W3P57el8bbm/ZsTgL+m57ZLem5VJI9na+T+Hnw9hPABSfL/fyTrsv0kXbbjeZJk/iSSYdrrWtmG\nujh7vtQR1QNfIfkw+zHJKtdnk6zjMxHYNybbDDX2e/lDBheR/Le6jGQV+1EkPWanxxivyFWKyWrr\nN5PMZTmPdUNRG9twuLHhwQfTOA8kGRr9gGS9pvU2eU7XAfseybyZb5OsQXQ1yWrghW3ek37tlp5L\nbuhpvXoxxnkkwy7TgENJ1h7qTdLz9dlW/Offok2WY4w3kXwY/5lkiOp0knk1Z6fzk9qq0fONMX6N\nZJmD5STrkA0hSaKOijHeVnCMi0kSzhOBb6ZD018h+RCtIFmUdBDJGmwHkgyh7ZHfy9OIRh+XdGmD\n3PpthYuPbkxTj/MYkrlEtSTn/1WSi0q+RLJgLyQLxTZ3nA3EGC8meX2uJHkv7E2yZMTMRuq29vlt\n6fukNVr0/s1r5xskw4tnkvzjcQdwSDpkuTGF76WFJEvE1JH8PfhM+v45HPgZyVWn3yOZQ/hQ2saL\nbT9NdUUl9fWba7qLJHUt6Zy6fwILYoyfL3I4XVI67+o04NNxI4v4Sh1Ji4YdQ7J1yXUxxi+EED5F\n0vtQQzLuflqMcUkIYRRJt+xakr3kZqaTRH9F0j2b+4+l3a+SkaSMjCJZyPSSYgciqfNodtgxhPB9\nkvHu3IJ444Hz0i0fpgMXpwsXnk8y3HEMMDZd9+Rc4PkY4+EkVwFtymXdktQhhBCmhRD+QjIM9TLr\nT/iXpI1qyZyv10jG/HO+HmPMrZLdneTqnM8Ac2KMNenWFv8g3W+MZBIwJAsPNrthryR1Au+STIJ/\nmmSvxNoixyOpE2l22DHGOD29DD53ezFACOFzJJMRDyfp7cpfAHAlyYrCvfLKV9D0PnwNQgg9SCYO\nv00y6VSSOpob0y+gYV0qFccV6ZfPg4qllGT9u/npnr7NatNSEyGEr5PMcTguxvh+CGE56ydWvUhW\nVF7Ouj34epGsSN6cg2lkTzFJkqQObBDr9pvdqFYnXyGEb5JMrP983uXrzwBXpwvgbUuyIN2LJBu2\nHkdyufJxtCypehtgypQp7Ljjjq0NT1KGamtrWbx4Md27u2Tg5lJTU0O/fv0oLd3cu181YWC6t/ic\nFn1mSCrwzjvvcMopp0ArNm9v1V/MdKXlCSSXVk8PIdQDT8YYrwghTCTJ+EqAMTHG6hDCz4A7Qwiz\nSa6M/EYLmqkF2HHHHenfv39rwpOUsZqaGkpKSky+NqOamhp22mmn7B9T/95Km6rFU6Va9O6OMf6T\ndXvUNboSdIxxEsmO7/llq0lWAJYkdUSVlcWOQOpyXOFekiQpQyZfkiRJGTL5kiRJypDJlyRJUoZM\nviRJkjJk8iVJXVlFRfIlKTMmX5IkSRky+ZIkqQOoqalh4MCBjBo1qkX1zzrrLJYubcmufY27+eab\nufrqqzconz59OgcddBDDhg1j2LBhHH/88Zx++uk8//zzDXW+/e1v8/rrrwNw2WWXMXjwYMaPH8/c\nuXM58sgj+drXvkZ1dXWbY9vSuSy1JEkdwKOPPsqee+7JSy+9xIIFC9htt902Wn/u3LntFstBBx3E\nrbfe2nB73rx5fPvb3+aBBx7g4x//OLfddlvDfb/+9a954okn6NevH2PGjOGkk07inHPOabfYtgT2\nfEmS1AHcfffdHH300Rx33HFMnjy5ofy+++5jyJAhDB06lDPOOIN33nmHSy65BIDTTjuNd955hyOP\nPJKXXnqp4Xfyb99666187WtfY+jQoXzxi1/ksccea3Vshx56KEcffTRTp05d7/jpnoaMGjWKW265\nhVmzZjF16lSuv/76hraHDx/OsGHD+O53v8uSJUsAOPXUUzn//PMZMmQIU6ZMYeXKlVxyySWccMIJ\nDB06lOuuu466ujoA9ttvP26++WZOPvlkBg8ezJ133tkQ12233caxxx7L8ccfz/nnn8/KlSsbHrPh\nw4czfPhwRo4cyYIFC1p9zu3Jni9JUtfT1EUGTW231Nr6rfTaa6/x/PPPc8stt1BZWclpp53G6NGj\nefvtt7nhhhuYMWMG/fr145e//CW33norY8eOZfr06dx1112Ul5c3edxFixbx9NNPM2XKFLbeemt+\n97vfMXHiRAYPHtzqGEMIzJ49e72yKVOmsOeeezbEUVVVxR577MGZZ57JjBkzePXVV7nvvvvo1q0b\nv/71r7n00ku5/fbbASgvL+ehhx4CYMyYMeyzzz6MHTuWuro6fvCDH3DHHXdw1llnUV1dzfbbb8/U\nqVN56aWXOPnkkzn55JOZPXs2M2bM4N5776Vnz5785Cc/YcqUKRxwwAHMmDGDqVOn0qNHD+bOncv5\n55/PzJkzW33O7cXkS5K6Mvd27BCmTZvGEUccQa9evdh3333ZeeedmTZtGj169GDQoEH069cPSHq6\n8tXX12/0uDvttBPXXXcdv/nNb1i4cCHPPfccq1atalOMJSUlbLPNNo3e11gcTzzxBC+88ALDhw8H\noK6ujjVr1jTcf9BBB21Q99577wVgzZo1dOu2bnDuqKOOAmDvvfdm7dq1rF69mnnz5nHMMcfQs2dP\nAC6++GIArr/+ehYuXMiIESMa4lq+fDnLly9nu+22a9O5b24mX5Kkrqe1SWc7JqmrV69mxowZbLPN\nNhx11FHU19fzwQcfcPfdd3P22WevV3fNmjW89dZbDfPBSkpKGr7nJ0Br164F4OWXX+Y73/kOZ5xx\nBgMHDuTggw/miiuuaFOcL7zwAnvssUeL69fV1TFq1ChGjBjRENPy5csb7i8rK1uv7oQJExrOa8WK\nFQ3nBtCjR4/1jl1fX0/37t3Xq7NixQqWL19OXV0dQ4cOZfTo0Q33LV68uMMkXuCcL0mSiuq3v/0t\n22+/PXPmzGHWrFk8/vjjPPbYY6xatYply5Yxb9483nvvPQCmTp3KuHHjACgtLW1Isvr06cOLL74I\nwHPPPddQf/78+ey7776cccYZHHzwwTz22GMNc6la48knn+Spp55qSKRaYuDAgdx7770N87DGjx/P\nRRdd1GTd3Dy36upqzj33XKZMmdJo3VySeeihh/Loo4/ywQcfAHDTTTcxefJkBg4cyMyZMxvml02Z\nMoUzzjijxXFnwZ4vSZKKaNq0aZx55pnrlfXq1YtTTz2VJ598kosuuoizzjqLkpISdthhB8aOHQvA\n0UcfzTe+8Q1uueUWRo8ezeWXX84999zD3nvvzd577w3AkCFDeOSRR/jyl7/M1ltvzSGHHMLSpUub\nHXp89tlnGTZsGJD0qvXt25dJkyax/fbbN5Tl5P+c72tf+xrvvvsuX//61+nWrRsf//jHue666xr9\nnUsvvZRrr72W448/npqaGg477LCGXr/CurnbRxxxBAsWLGDEiBGUlJSw++67c9VVV1FWVsbZZ5/N\nyJEj6datGz179uTmm2/e6PlmraS58eKshRAqgDdmzZpF//79ix2OpI2oqalh0aJFdO/u/3GbS01N\nDTvttJOPqdRJvPnmm7k5abvGGCtb8jsOO0qSJGXI5EuSujL3dpQyZ/IlSZKUIZMvSZKkDJl8SZIk\nZcjkS5IkKUNeyyxJ6nJqamra9fguFaKN8dUhSV1ZF9zbsaamhqqqKkpLS9vl+LW1teyyyy4mYGqS\nw46SpC6ntLSU7t27t8tXa5O6uro67rjjDk444QSGDRvGkCFDGDduHNXV1W0+v7q6Os4991yOOeYY\npkyZwrBhwxq2+cn3i1/8gksuuaTN7WwORx55JC+99NJG66xcuZLTTz89k3guueQS7rjjjnZtw7Rc\nkqQiuuyyy1ixYgV33nknPXv25MMPP2T06NH86Ec/4ic/+UmbjvnOO+/wpz/9ieeee46SkhJOOeWU\nzRx1tpYuXcoLL7xQ7DA2G5MvSZKK5M033+Shhx5i7ty5lJWVAbDNNttw5ZVX8te//hVIen2uuOIK\n/v73v1NSUsKgQYMYPXo03bp1Y7/99uNb3/oWc+fOZcmSJZx22mmccMIJjBo1ipqaGoYPH87EiRM5\n+uijefrpp+nZsydXXXUV8+bNo0+fPvTp04devXo1tHPNNdfw6quvUlNTw6GHHspFF13UaDunnnpq\nQ0/UbbfdxowZM+jevTsVFRWMHTuWnj17ct9993H33XcD8NGPfpQf/vCH7Lbbbht9PBo7n9NOO40x\nY8bw4YcfMmzYMB544AEWLFjAtddey9KlS6mrq+PUU09l+PDhPPPMM1xzzTVsu+22fPjhhwwYMIC9\n996bkSNHAsnG5PPnz+eGG27gmmuu4YUXXuCDDz6gvr6eq6++mk9/+tPt8jwXcthRkqQiefnll9l9\n990bEq+cPn36MHjwYACuvvpqevfuzYMPPsj999/P3//+dyZNmgRAdXU122+/PVOnTmXChAmMGzeO\nrbbaittvv50ePXowffp0dtlll4bNqKdMmcLChQt5+OGH+cUvfsGiRYsa2rz22mvZZ599uP/++5k+\nfTr/+te/GobfCtu54YYbqK6uZtasWcyYMYN7772XBx98kP79+zNlyhTmz5/PjBkzmDp1Kg888ABn\nnXUW559/frOPR2PnU11dzdixY9lmm22YPn06dXV1fO973+PCCy/k/vvv56677mLSpEk8//zzALz2\n2muMHz+eGTNmcNJJJzF9+vSG40+fPp2TTjqJv/3tb7z33nvcc889PPTQQwwdOpTbb799E57J1rHn\nS5KkIunWrRt1dXUbrfPUU08xbdo0ALbaaitOPvlk7rzzTkaNGgWQ29SZvffem7Vr17J69eomjzVv\n3jyGDBlCaWkp2267LV/5yleIMQLwxBNP8MILL3DvvfcCsGbNGrp1W9dH01g78+bN45hjjqFnz54A\nXHzxxQBcf/31LFy4kBEjRlBfXw/A8uXLWb58Odttt91Gz7e586msrGThwoWMGTOm4dhr1qzh5Zdf\nZrfddmPHHXdkxx13BOCzn/0s1dXVvPTSS2yzzTb8+9//5pBDDgHge9/7HlOnTmXhwoU888wzDeeQ\nBZMvSerKcvs6dsGrHjuCfffdl9dff51Vq1at1/u1ePFifvzjHzNx4sQNkrO6urr1lsro0aNHw8/1\n9fUNCUljSkpK1rs//+KA2tpaJkyY0DA0uGLFioYes8J2cm117959vTorVqxg+fLl1NXVMXToUEaP\nHr3eOTWXeLXkfGpra9luu+3W69F6//336dWrF88999wGvYgnnngi06dPZ+utt+bEE08EkkTz2muv\nZeTIkQwePJjddtuNBx98sNnYNheHHSVJXU5tbS01NTXt8lVbW9viOPr168fxxx/PmDFjGq5GzM3x\n2n777enRoweDBg1iypQpQDIsd88993DYYYc1e+z8pCX386BBg/jNb35DdXU1a9as4Xe/+11DnYED\nBzJ58uSGds4999yGdps69qGHHsqjjz7KBx98AMBNN93E5MmTGThwIDNnzmTJkiVAMtx5xhlntPhx\nKdS9e/eGJHTXXXelR48e/Pa3vwXg7bffZsiQIU1eMTls2DAef/xx/vCHPzB8+HAA/vSnP3HkkUcy\nYsQI9tlnH2bNmtVsD+TmZM+XJKlL6d69O7vssku7t9FSl19+Of/7v//LySefTPfu3amurmbw4MEN\nc6QuvfRSrrrqKo4//njWrl3LoEGDOOeccwDW63UqvN3YzyNGjGDhwoUMGTKE3r1788lPfrKhzg9/\n+EOuvfZajj/+eGpqajjssMM4++yzN9rOEUccwYIFCxgxYgQlJSXsvvvuXHXVVZSVlXH22WczcuRI\nunXrRs+ePbn55psbPf+mYs6/vcMOO7DXXntx3HHHMXXqVG655Rauvvpqfv7zn1NbW8t///d/8+lP\nf5pnnnlmg+N/7GMfY5999qG2tpYddtih4XG48MILGTp0KKWlpRx00EE88sgjjcbXHko21j1ZDCGE\nCuCNWbNm0b9//2KHI2kjampqWLRokYtJbkY1NTXstNNO2T2mDjtKm+TNN9/MzVPbNcZY2ZLfcdhR\nkiQpQyZfkiRJGXKsQJK6MocbpczZ8yVJkpQhky9JkqQMtWjYMYTwWeC6GOMXQggDgMlAHfBijPG8\ntM4o4FvAWuCaGOPMEMI2wK+AvsBy4PQY4/ub/zQkdXW1tbVUNjKEVlFRsd5CkpJUbM0mXyGE7wOn\nAivTohuBMTHG2SGEn4UQhgJPA+cDBwBlwJwQwiPAucDzMcYrQwhfB34EXNAO5yGpi6usrOTC8bMo\nK+/bUPbYGJ+BAAAUFklEQVTB0ne44MQ911vTyWRMUrG1pOfrNWAYcFd6+8AY4+z054eBL5L0gs2J\nMdYAy0MI/wD2BwYCP8mr+6PNFbgkFSor70vP3js33F61bDET7n+VsvKl6e13GXfBUQwYMKBYIUpS\n83O+YozTgZq8ovzlZ1cA2wG9gGV55SuB8oLyXF1JykwuIevZe+f1esWUqqhYt9CqpEy0ZcJ9/uZH\nvYClJPO5tiso/3da3qugriRJUpfVluTrLyGEw9OfjwVmA/OBgSGErUMI5cCewIvAn4Dj0rrHpXUl\nSZK6rLYkXxcCV4YQ5gJbAffFGBcDE4E5wGMkE/KrgZ8B+4QQZgNnA1dsnrAlSZI6pxYtNRFj/Cfw\nufTnfwCfb6TOJGBSQdlq4KRNjlKSJGkL4SKrkiRJGXJvR0nqytzbUcqcPV+SJEkZMvmSJEnKkMmX\nJElShky+JEmSMmTyJUmSlCGTL0nqytzbUcqcyZckSVKGTL4kSZIyZPIlSZKUIVe4l9Rl1NfVUVVV\ntUF5RUUFpaWlRYhIUldk8iWpy1i9YgkT7n+PsvKlDWWrlr3LuAuOYsCAAUWMTFJXYvIlqUspK+9L\nz947N9zu8r1h7u0oZc7kS1KXZm+YpKyZfEnq8gp7wySpPXm1oyRJUoZMviRJkjJk8iVJkpQhky9J\n6src21HKnMmXJElShky+JEmSMmTyJUmSlCGTL0mSpAyZfEmSJGXIFe4lqStzb0cpc/Z8SZIkZcjk\nS5IkKUMmX5IkSRlyzpekTqe2tpbKgrlKVVVVxQlGklrJ5EtSp1NZWcmF42dRVt63oez9N1+hT/+9\nihiVJLWMyZekTqmsvC89e+/ccHvVssVFjKYTy+3r6FWPUmac8yVJkpQhky9JkqQMmXxJkiRlyORL\nkiQpQyZfkiRJGfJqR0nqyrzKUcqcPV+SJEkZalPPVwihO3AnUAHUAKOAWmAyUAe8GGM8L607CvgW\nsBa4JsY4c5OjliRJ6qTa2vN1HFAaYzwMuAq4FrgRGBNjPALoFkIYGkLoB5wPHAocA4wNIWy1GeKW\nJEnqlNqafL0KdA8hlADlJL1aB8QYZ6f3PwwcDXwGmBNjrIkxLgf+Aey3iTFLkiR1Wm2dcL8S2BX4\nO9AHOB4YlHf/CmA7oBewrOD3ytvYpiRJUqfX1p6v/wZ+H2MMwP7AL4Gt8+7vBSwFlpMkYYXlkqSO\noKJi3f6OkjLR1uTrX6zr0VpK0oP21xDCEWnZscBsYD4wMISwdQihHNgTeHET4pUkSerU2jrsOB74\nRQjhKWAr4AfAs8DP0wn1rwD3xRjrQwgTgTlACcmE/OrNELckSVKn1KbkK8b4AfD1Ru76fCN1JwGT\n2tKOJEnSlsZFViVJkjJk8iVJkpQh93aUpK7MvR2lzNnzJUmSlCGTL0mSpAyZfEmSJGXI5EuSJClD\nJl+SJEkZMvmSpK7MvR2lzJl8SZIkZcjkS5IkKUMusipJBerr6qiqqlqvrKKigtLS0iJFJGlLYvIl\nSQVWr1jChPvfo6x8KQCrlr3LuAuOYsCAAUWOTNKWwORLkhpRVt6Xnr13LnYYkrZAJl+SOrTa2loq\nC/YfLBwS1CZwb0cpcyZfkjq0yspKLhw/i7Lyvg1l77/5Cn3671XEqCSp7Uy+JHV4hUOAq5YtLmI0\nkrRpXGpCkiQpQyZfkiRJGTL5kiRJypDJlyR1Ze7tKGXO5EuSJClDJl+SJEkZMvmSJEnKkMmXJElS\nhky+JEmSMuQK95LUlbm3o5Q5e74kSZIyZM+XJDWjvq6OqqqqDcorKiooLS0tQkSSOjOTL0lqxuoV\nS5hw/3uUlS9tKFu17F3GXXAUAwYMKGJkkjojky9JaoGy8r707L1zscOQtAVwzpckSVKGTL4kqStz\nb0cpcyZfkiRJGTL5kiRJypDJlyRJUoZMviRJkjJk8iVJkpShNq/zFUL4AfAVYCvgFuApYDJQB7wY\nYzwvrTcK+BawFrgmxjhzE2OWJG0u7u0oZa5NPV8hhCOAQ2OMnwM+D3wCuBEYE2M8AugWQhgaQugH\nnA8cChwDjA0hbLVZIpckSeqE2jrs+CXgxRDCDOC3wEPAATHG2en9DwNHA58B5sQYa2KMy4F/APtt\nYsySJEmdVluHHT9G0ts1BNiNJAHLT+RWANsBvYBleeUrgfI2tilJktTptTX5eh94JcZYA7waQvgQ\n6J93fy9gKbCcJAkrLJekRtXW1lKZNw+pqqqqeMFIUjtoa/I1B/gv4KchhJ2AjwCzQghHxBifBI4F\nHgfmA9eEELYGtgX2BF7c9LAlbakqKyu5cPwsysr7AvD+m6/Qp/9eRY5KkjafNiVfMcaZIYRBIYRn\ngBLgXKAS+Hk6of4V4L4YY30IYSJJslZCMiG/evOELmlLVVbel569dwZg1bLFRY5mC5fb19GrHqXM\ntHmpiRjjDxop/nwj9SYBk9rajiRJ0pbERVYlSZIyZPIlSZKUIZMvSZKkDJl8SZIkZajNE+4lSVsA\nr3KUMmfPlyRJUoZMviRJkjJk8iVJkpQhky9JkqQMmXxJkiRlyORLkrqyiop1+ztKyoTJlyRJUoZM\nviRJkjJk8iVJkpQhky9JkqQMmXxJkiRlyL0dJakrc29HKXP2fEmSJGXI5EuSJClDJl+SJEkZMvmS\nJEnKkMmXJElShky+JKkrc29HKXMmX5IkSRky+ZIkScqQyZckSVKGTL4kSZIyZPIlSZKUIfd2lKSu\nzL0dpcyZfEkqmtraWioLPvyrqqqKE4wkZcTkS1LRVFZWcuH4WZSV920oe//NV+jTf68iRtUy9XV1\njSaKFRUVlJaWFiEiSZ2FyZekoior70vP3js33F61bHERo2m51SuWMOH+9ygrX9pQtmrZu4y74CgG\nDBhQxMgkdXQmX5LURoWJoyS1hFc7SpIkZcjkS5K6Mvd2lDJn8iVJkpQhky9JkqQMmXxJkiRlyORL\nkiQpQ5u01EQIoS/wZ2AwUAtMBuqAF2OM56V1RgHfAtYC18QYZ25Km5IkSZ1Zm3u+QgjdgVuBVWnR\njcCYGOMRQLcQwtAQQj/gfOBQ4BhgbAhhq02MWZK0uVRWur+jlLFNGXYcB/wMWASUAAfEGGen9z0M\nHA18BpgTY6yJMS4H/gHstwltSpIkdWptSr5CCGcA78YYHyVJvAqPtQLYDugFLMsrXwmUt6VNSZKk\nLUFb53ydCdSFEI4G9gd+CeyQd38vYCmwnCQJKyyXJEnqktqUfKXzugAIITwOnANcH0I4PMb4FHAs\n8DgwH7gmhLA1sC2wJ/DiJkctSZLUSW3OjbUvBP4vnVD/CnBfjLE+hDARmEMyPDkmxli9GduUJEnq\nVDY5+YoxHpl38/ON3D8JmLSp7UiS2kFuX0eveJQy4yKrkiRJGTL5kiRJytDmnPMlSV1afV0dVVVV\n65VVVFRQWlpapIgkdUQmX5K0maxesYQJ979HWXmyos6qZe8y7oKjGDBgQJEjk9SRmHxJ0mZUVt6X\nnr13LnYYkjowky9J6sq8ylHKnBPuJUmSMmTyJUmSlCGTL0mSpAyZfEmSJGXI5EuSJClDJl+S1JVV\nVKzb31FSJky+JEmSMmTyJUmSlCGTL0mSpAyZfEmSJGXI5EuSJClD7u0oSV2ZeztKmbPnS5IkKUMm\nX5IkSRly2FFSJmpra6ksGOKqqqoqTjCSVEQmX5IyUVlZyYXjZ1FW3reh7P03X6FP/72KGJUkZc/k\nS1Jmysr70rP3zg23Vy1bXMRoJKk4nPMlSV2ZeztKmTP5kiRJypDJlyRJUoZMviRJkjJk8iVJkpQh\nky9JkqQMudSEJHVl7u0oZc6eL0mSpAyZfEmSJGXI5EuSJClDJl+SJEkZMvmSJEnKkFc7SmoXtbW1\nVOZdSVdVVVW8YNS03L6OXvUoZcbkS1K7qKys5MLxsygr7wvA+2++Qp/+exU5KkkqPpMvSe2mrLwv\nPXvvDMCqZYuLHE326uvqGu3xq6iooLS0tAgRSeoI2pR8hRC6A78AKoCtgWuAl4HJQB3wYozxvLTu\nKOBbwFrgmhjjzE2OWpI6gdUrljDh/vcoK1/aULZq2buMu+AoBgwYUMTIJBVTW3u+vgm8F2M8LYTw\nUeBvwHPAmBjj7BDCz0IIQ4GngfOBA4AyYE4I4ZEY49rNEbwkdXT5vX+SBG1Pvn4N3Jv+XArUAAfE\nGGenZQ8DXyTpBZsTY6wBlocQ/gHsBzzb9pAlSZI6rzYlXzHGVQAhhF4kSdilwLi8KiuA7YBewLK8\n8pVAeZsilSRtfl7lKGWuzet8hRB2AR4H7owxTiPp5crpBSwFlpMkYYXlkiRJXVKbkq8QQj/gD8BF\nMcY70+K/hhAOT38+FpgNzAcGhhC2DiGUA3sCL25izJIkSZ1WW+d8XQJ8FPhRCOHHQD3wPeCmEMJW\nwCvAfTHG+hDCRGAOUEIyIb96M8QtSZLUKbV1ztcFwAWN3PX5RupOAia1pR1JkqQtjXs7SpIkZcjk\nS5K6soqKdfs7SsqEyZckSVKGTL4kSZIyZPIlSZKUIZMvSZKkDJl8SZIkZaiti6xKkrYE7u0oZc6e\nL0mSpAyZfEmSJGXIYUdJylB9XR1VVVUblFdUVFBaWlqEiCRlzeRLkjK0esUSJtz/HmXlSxvKVi17\nl3EXHMWAAQOKGJmkrJh8SVLGysr70rP3zsUOQ1KRmHxJ2mS1tbVUFlw119jQmjqg3L6OXvUoZcbk\nS9Imq6ys5MLxsygr79tQ9v6br9Cn/15FjEqSOiaTL0mbReFQ2qpli4sYjSR1XC41IUmSlCGTL0mS\npAyZfEmSJGXIOV+S1JV5laOUOXu+JEmSMmTPlyQVWf6WQ7W1taxcuZIQgtsNSVsoky9JKrLCLYdW\nLXuCu8Z+gz322KPIkUlqDyZfktQBuOWQ1HU450uSJClDJl+S1JVVVKzb31FSJhx2lNQqtbW1vP76\n6wDU1NTw7rvvsmjRoiJHJUmdh8mXpFZ5/fXXOfWSu91EW5LayORLUqu5ibYktZ1zviRJkjJk8iVJ\nkpQhhx0lqStzb0cpcyZfktTB1NfV8cYbb2xQPmDAALcckrYAJl+S1MGsXrGEH9/+HmXlrzeUrVr2\nrlsOSVsIky9J6oDcbkjacjnhXpIkKUP2fEnaqPwV7YFG5yJJklrO5EvSRhWuaO9q9sXRbpPwc/s6\netWjlBmTL0kNCnu5IOnpyp9/5Gr2xeEkfGnL0e7JVwihBLgF2B/4EDg7xrigvduV1Hru29ixFU7C\nb6w3zOUopI4vi56vrwI9YoyfCyF8FrgxLZOUocJerdraWoD1PqgLe7nAnq6OrLA37IOl73DVtw9j\n1113Xa+eCZnUsWSRfA0Efg8QY/x/IYSDMmhTUoHG5m5t26uPvVydXOGQ8I9vn7fe0KQJmdTxZJF8\nbQcsy7tdE0LoFmOsa6J+KcA777zT7oFJW6q/PPc8N0y4jZJuJQ1l/T62PbXV21Pz4TYA1NesprZ6\nZcPtXNmKJQuo+XB5Q9mqf79FbfUHDWWFt1talvXvdYQYivF7PT7y0fWe09VLF3HhddPoUfbRhrI1\nq5Zyyaij2WWXXdgl7QHt/uabSGq9vHylxf/NlNTX17dPNKkQwg3AvBjjfenthTHGT2yk/kBgdrsG\nJUmStHkNijHOaUnFLHq+5gJDgPtCCIcALzRTfz4wCHgbqG3n2CRJkjZFKfBxkvylRbLo+cpd7bhf\nWnRmjPHVdm1UkiSpg2r35EuSJEnruLejJElShky+JEmSMmTyJUmSlKEOt7djCGEYcGKM8ZT09leB\nccDCtMplMUaXotiMGnnMPwtMANYCj8YYryxmfFuqEMKbQO7ik3kxxkuLGc+Wyi3OiiOE8Czr1nh8\nI8Z4VjHj2ZKlf7OvizF+IYQwAJgM1AEvxhjPK2pwW6iCx/xTwEOs+3v+sxjjvRv7/Q6VfIUQxgNf\nBJ7LKz4Q+H6McXpxotqyNfGY3woMizFWhhBmhhD2jzH+rTgRbpnSP5DPxhiHFjuWLsAtzjIWQugB\nEGM8stixbOlCCN8HTgVWpkU3AmNijLNDCD8LIQyNMf6meBFueRp5zA8Ebogx/rSlx+how45zgXML\nyg4ERoYQngohjAshdLSYO7v1HvMQQi9g6xhjZVr0B2BwEeLa0h0I9A8hPB5CeCiEsEexA9qCrbfF\nGeAWZ+1vf+AjIYQ/hBAeS5NetY/XgGF5tw/MGx16GP9+t4cNHnPgyyGEJ0MIPw8hfKS5AxQlkQkh\njAwhvBBCeD7v+4FNdNM9ApwfYzwc6Amck220W4ZWPObbAcvzbq8AyrOLdMvT2GNPsojwtWnPwFjg\nV8WNcovW6BZnxQqmi1gFXB9j/BLJP3dTfMzbRzoqVJNXVJL3s3+/20Ejj/n/IxmhOwJYAFze3DGK\nMuwYY/wF8IsWVr8jxpj7w/kbYHj7RLVla8VjvpzkwyqnF7C0XYLqIhp77EMI25K+eWOMc0MIHy9G\nbF3EcpLXcc7G9pbV5vEqSe8AMcZ/hBDeJ1kB/K2iRtU15L+2/fudjRl5ecp0YGJzv9AZ/hN5PoSw\nU/rzUcCzxQxmSxdjXAGsCSHsmk5U/hLutdkeLgMuAAgh7A9UFTecLdpc4DiAFm5xpk03ErgBIP37\n3Yukt1ft7y8hhMPTn4/Fv99Z+EMIITedoUV5SoeacN+Es4DpIYRVwMvA/xU5nq7gHOBukuT8kRhj\ni/erUotdB/wqhPBlkqtKzyhuOFu06cDRIYS56e0zixlMFzEJuCOEMJukJ2akvY2ZuRD4vxDCVsAr\nwH1FjqcrOBe4KYRQDbwDfKu5X3B7IUmSpAx1hmFHSZKkLYbJlyRJUoZMviRJkjJk8iVJkpQhky9J\nkqQMmXxJkiRlyORLkiQpQyZfkiRJGfr/AXzIDDvV1Q7BAAAAAElFTkSuQmCC\n" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Let's look at the distribution of differences when voting is *totally* random\n", "f, ax = plt.subplots(figsize=(10, 5))\n", "_ = ax.hist(diff, bins=30)\n", "_ = ax.axvline(actual_diff, color='r', ls='--')\n", "axfill = ax.fill_between([clo, chi], *ax.get_ylim(), alpha=.1, color='k')\n", "ax.set_title('Distribution of totally random splits', fontsize=20)\n", "ax.legend([ax.lines[0], axfill], ['Actual Difference', 'Confidence Interval'], fontsize=12)\n", "_ = plt.setp(ax, xlim=[-15, 15])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For a vote to be \"different\" than 50%, it'd need to be outside our margin of error described by the grey rectangle. In this case, it seems that a totally random vote yields about 2% points of spread around 0, and that the recorded vote difference (~4%) is outside of the margin of error for 50%. So maybe we can conclude that the Brexit vote was significantly different from a random 50/50 vote.\n", "\n", "BUT - we also know that people don't vote completely randomly. They are influenced by external factors, they talk to one another, they tend to vote similarly to those around them. This is why everybody could predict which districts would vote \"yes\" and which would vote \"no\" well before the election.\n", "\n", "So, let's build that in to our simulation...\n", "\n", "# Simulating a not-completely-random population\n", "So how exactly do we simulate the fact that people don't vote *totally* randomly? There are a lot of ways to do this, but I'll take the semi-arbitrary decision to say that we could expect the same pattern of voting to occur *within a district*. That is - we can simulate random *district* votes instead of random *individual* votes. Moreover, we'll then weight that district's percentage split by an amount proportional to that district's size. Intuitively it doesn't seem like this should make much difference in our simulation (we're still totally randomly choosing the yes/no split), but let's see what happens...\n", "\n", "First, I grabbed a list of each UK voting district, along with its size..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [ { "data": {}, "execution_count": 19, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEcCAYAAADgJkIVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecXHW9//HXzGzJ1vRKekI+BISQBAQCCS2oiAh4VVSw\n4A8Rrh3lKiiKXlGUC5aLcqWIoFIEBQSUIkjvPbQPCUlIIT1ks7spW39/fM8kJ5PZ3Ukysy3v5+Ox\nj90zp33m7Mz5nG8535NobW1FRERkVyW7OgAREekdlFBERCQvlFBERCQvlFBERCQvlFBERCQvlFBE\nRCQvdsuEYmYLzawl9rPJzN4ys4vNrDK23OFm1mxmI3Lc7qfNbFA787fZnpktMLPzdvG9HGxmM2LT\nLWb2qV3Z5i7GU2RmN5pZvZktaWe5Pc3sZTMrKlAcKTP7+i5uY4/oeM6Kpr9oZv+bnwh3OJZ/m9kV\nXbTvMdFxmNHx0l3HzD5rZg2x6clm9sHYdLvfNzP7gZm9uQv7LzOzszK2N3dnt9fGPtrcZua8tj4z\nZnaJmTWZ2Sn5jA1204QCtAI/BYZFP5OB7wCfAP4ZO8k9Bgx393c62mD0ZbsWKG9nsZy3twMeBibG\npocBt+Rx+zvqGODjwH8AB7Wz3JXABe7eVKA4TgYuycN24jdqXQkcZWaH5mG7O+ok4Owu2C/AIsLn\n6qku2n+ubgT2iE3fDhywA+tfDBy8C/v/BnBOxmuFuNGvvW22uz8z+ynwFeBUd/9zXqMCCnJ12EPU\nu/vK2PQCM5sHPAt8HrgiOtmtzLr29pJ08M/cwe3lKpGxj3xvf0cNAFrd/e62FjCzDwBj3P1vBYwj\nXxdLW46vu7eY2a+BnwCH52n7OXH3dZ25v4x9t5L/z23euftmYFXspURby7ax/gZgwy6E0K0v0M3s\nR8A3gVPc/eZC7GN3TijbcfcXzOxRQknlCjM7AngAGOnu75jZccB/A3sB6wglgW8BwwklhVZCYvoh\n8BBwN/BjwpXli8CPgH+ntxftdqSZ3QvMJFwJXuDuN0AowhKuJPZMx2hmFwCfcvdJZraA8CH+g5l9\nzt2PMrOWaJ3ro+U/T7hymgAsBX7p7r+J5n2WUDK7FDgveh9PA19wd892jMysDPgBoQQwDHgJOM/d\nH4ji/UG0XDPwQ3f/UZbNfJ1YKcrMksDPo+M+CHDgv939lmh+ivBFOB0YBcyN5t8cO06HA2uA9wFz\ngBmxOE5z9+vMbCZwETAVeAf4SxTj5mjZ0cDlwCzCCfQnbH+R8DfgMjOb5u7PZzk+mbH82t3PN7MT\ngQsAAxYAVwOXRidrzGxP4DdR3CuB7wPXAEe7+8Nm9m9grrufES1/GOGzNY1wErwJ+I67bzSzMdE+\nPgp8F9gHWAh8291vz/L/wMwmAf9LuEJvJXxOv+7ub8e2dxhQHM1rJZyw47/HuvtiM3sP8D+Ez/Ra\n4K5o3zXRvrJ9j85J/x9iMfUjJIjj0xcoZnYd8BGgr7s3R5+N1cDngH7AVe5eHB2vCcAF0XdjfLTZ\nkWZ2OzAbqAH+191/Gvvfnerue+7oMYy+Sz+K/m4GjoxmJczsu8CXgL7Av4DT3X1VtOxI4JeEkv3G\n6Nie7e7Lsv2fdlYUw7eBT7r7X/O57bhunVG7yBxg3+jv1ugHMxsI/BX4LTAJ+BShaue/CInghGid\nAwlfJoBS4Ijota/Gthn3BeAf0T6vBf5kZtNj8zOXj08fCLRE2/5I5hsxs7MJJ4lLo+3/HLjYzL4R\nW2x89F5OIlRRDYjWactNhC/ZF4ApwJPA3WZ2IKHK4MtRjMPYehziMVUQvmx3xV7+EuH4nUQ4tjcD\n10dfaoBfEBLKt6P3cQNwo5mdFNvGEYREMxU4LSOOm8xsf0KCv4Vwcjgd+BDh/0lUzXkP0Ac4hFBK\n/U5m/NGJ4Fngw+0co3gsV0X1+H+K3sfehM/MV4HvRfsuJ5xoNhD+p18gnJyyfj/N7CDgfkIV1AHA\nZ6Pjd2PGoj+P3sPehAuaP0QXBNlcTzhh7k9IHAMJSS8t/bl7jHBMh0e/9yFcqFwXJZM9gAej/U0h\nVH1OJiTi9r5HmVVF6VLZ44STf9pRhP/Re6PpQ4ES4N6MOD8SvZ//Ydtqr88Dd0bH5FfAhVFyznyf\nabkewxuBnwGLo+PyRPT6hOg4HAm8P4o7ncDKCceqjpDI30dI2Pfns23RzM4hJPDfFjKZgEoo2bwL\nVGd5fRThn73U3ZcAS6Kqmzp3bzWztdFyq919g5ml1/uZu8+H0CifZbt/cfdfRn//xMxmA18DPtNR\noO6+OtrP+jaqRM4hXAVfE02/ZWYTCCe0X0SvFQFfdPc3oxivAC7Mtj8zm0w4CR/j7vdHL389OsF9\ny91PNrOaKLZV2bZBuKIuAl6NvTaBcDJd5O4rgB+b2VPAWjOrAs4EznL3W6Plf2pmUwhf9PRrLWxb\n2nhvPA4z+yZwp7un3/cCMzsTeDRqqJ0G7AnMdvel0TpfJZx8Mr1C+3XtmbH8CfiNu18b23c1oU3m\nvwmlvb6Eq+M64HUz+wrw9za2/03gGXf/djT9ZvRe/hH9j9LVNj939/uiGH4CfIxwYnwuyzYnEhLq\noujK/1TCiTEtAdtW20alg+sJF1RnRMudCbzl7luSsYVOIoujz8lm2vgetfFe7wROjbazF1BJKP0f\nTjhpfwB4ICqZbVnJ3d+NSgp17r42tr2/uPuV0d8/M7PvEBLOo23sP6dj6O6bzawOaI595oje72fc\nfRPgZnYT4YIDQjItJ5Sg0xeupxBKZf9BuHjbVccRLg4eBP6fmf3G3eflYbtZKaFsr5pQDN+Gu79o\nZn8B7rLQe+le4DZ3z3bCSWslFJvb80TG9DOE4u8uMbPBwNAs238YOCean44x/gGrIVzxZfOeaPnM\nbT5C+ODmYmj0e3Xstd8SSidLzexZQkniz+5eG5V8Um28j+Nj08syq0wyTAUmmllt7LUE4eQ/mXCl\nvTqdTCJPkr0efhXtN/ZmxjIVOMDM/jP2WhIoNbOx0fzXo2SS9mgb+yaK9a6M1x6Jfr+HUG0JoZSU\nVhNtr63/7fcIJdkvmdkDhBN5Zokn028IpYzp7t4YvTYVmJpxnCF8bia7+x928Ht0J3CRmQ0BjiYk\nk+cIJ+WLgA9GceQqs4fUOqCtUlvm8h0dw2zeiZJJ2rux/e0PDAHWx5NhNH9yG9trpO2apWQ0P24g\noQR4P6F6+gYzO6RQnWFU5bW9acAL2Wa4+ycJVye/JJRYbjWzq7MtG7Oxg/nNGdNJwlVNW3K9CGhr\nv6nod/qD1+LuLRnLtHUia2+bmR/ktqSrFNJxEJWOxhNODuk2rFfN7Mhon9niydxnR8e5gVCluB+h\nCmJK9PckQtVRui0gc51sUoRE1JbMWBoI7TFTYj/7RvteCjSxY9/FbO81vX78mGT7HGX937r7ZYQe\nUt8ANhGSyzNmVpxteTP7GqHkcEJGabSBkCTix3kKofT312hfOX+P3P11YD7hIutoQpvmv4EZUZXo\nvsAd2dZtQ+b3DdpvvM/5GO7E/hoIpd3MYzUJ+HUb23uXUJrNpj+hzSruenf/u7vXE2o9phI+iwWh\nhBITVaPMINR3Z86bbmaXuPsb7n6pu7+f0JB9crTIznYPnJoxfSjhQwbhA1eVMX9SxnTW/UZXu0ui\n7cXNBJbvZK+h12Ixxh0am9eRdGNjuoSEhb77H3X3e939W4SrMycU++cRjkO299HePjOPy6uEK+QF\n7j4/qoYcSuhaXEWoHx8UVQmmHZhlO+nYd6Tr96vAnun9RvuewtYv9svAXlH1Xlq6cTyb14g6HcTM\njJZ/fQfiAsDMBljovVbi7te4+6cIJ/DJUZyZy3+A0LZwmrtnXny9Gq23KPZeWwntFaNy+B5l8w9C\n+8NMwpX2k4Rz1wXAs+6+vI31OvvZHDu6v1eBccDa2LFaRaiO3reNdZ4DBkSdKDIdSqjhiNtSEnH3\nxwltSmeb2dE7GGtOducqr0ozS1e/lBMapC8i1DXG+2enryZqCNUBmwiNldWEKpcno/npIv40M1uX\nsW5c5mufNrOXgPsIjcnTCY2sEKp5fhxdDd5GuIL/ALAitn4tsLeZDc7SbvFj4FIzmx+9r6MIjdXf\nyxJXh9x9flQHfHmUBBYBXySU6r7a7spbvURIEFPZmlwGEXrj1BE6RUwHxhLanzaZ2aWE47A2Wv+j\nhCqy9k5CtRAuBIA3CA2mz5nZJcAVhPaBK4HF7r4y6hX0HKFTxJcJHSp+1ca2pxJdbefox8CdZvZq\ntJ4B/0do02k0sxuAHwLXmdn3CNUg6Y4R2U5SPwOeN7OLo/cwDrgMuMvdPdaZIVfvAscC46L2pI2E\nz+K7hMQ+IL2gme1NqAr7BfBA7DsE4TtyGaGTxbVmdhGhAf0ywlX1m4T/a3vfo2zuJLQn1br7K1Ec\njwCfJupV2IZaYJKZDc93r6l29tc/Otm/ncPyfyYk05vN7FxCaegiQnXqq9lWcPenzexh4BYz+xbh\nsz0M+E9CKX+7zjkZzif8r68zs/3cfU0OceZsdy6hfJtwlfkOoYrrXMKX/Lh0A1mkFSBqyDqB0OPk\nJcKV0iJCwxqEq8a/EXogXRBfN0Pmtv+HcGJ8iZAwjnP3udE+HyR8Yb5N+IAdRehOGncR4cOUvu9j\ny/bd/QrCB/Y7hFLP1wldQS9t45jk4vRoX38knIAPJDTSP93uWltjqidUWRwZe/lCwsnlN4QT2M+A\n77t7uqT4feB3hJPYy4SG0ZO9/ftY/k1oZ3mM0A36FUI7zwzC//vGaJmPRHG1EI7/IkK1yk2Eap9t\nmNkAwtXjbbm832jb9xBOfp8kJMzLgT8QGrDT908cS+j2+gxwVfR+IUu1m7u/SugcMYvwubmakKg+\nHluso89efHuthPcOW3toTQbe5+7pC6X0uh8jlOjOITTOvxP7+XjUqWI2ofT3JPBPQm+rY9y9KYfv\nUTYPEo7Dv2Ov3U+4OGur4wKE/98HgZfMLN21OVNrG6/TxuvtlUL+Skgk6e9yu6K2lWOAesL7eYRw\nTj7S3Ve3s+pxhAvQywkJ5e+Ezgoz3D2eyLaLNWrr+jThIuH3HcW4oxJ6YqN0NjM7lnDSHJWl/aZb\nszCcywnufmSHC+e+zdHARHd/IPbawYRkODqjo4BIt1XQKq/oquC3hHrYTYQbeubH5h9PKII1Ate4\n+1VR/+vfE4rGJcCF7n6HhfsI7iQUmwEu9wLd7SmF5e7/jKrhPkHodtojRN1kzySU0vKpHLg36ir8\nD7a27TykZCI9SaGrvE4ESt19BqFKaUsVQpQ4LiUUfY8Azoi6sp5K6L45i1ANcFm0ynTgEnc/KvpR\nMunZTgfOy+cNXJ3gDOB+d2/rnoWd4u5vEKrDvkioOr2D0Lj+0XzuR6TQCv1lPoyobt/dnzKzeN/9\nyYShJNYDWBjyZBZhOIx0soj3q55OaGA7kdA3/GtRfbz0QO7uhHsmegx3v7yA276ZrZ97kR6p0Aml\nmtDzI63JzJJRvXnmvFrC+DwbAKIulDcTxtGBcK/AlR7G2zqP0PC93XAN0bqlhMbiZWTvBy4iIttL\nEYbVeaaDG4WzKnRCWc+291EkY42w69l2iJMqojvUzWwU0SB87p4efuA2jwaXIwy30daNPxCSySPt\nzBcRkbbNpO3haNpU6ITyGKF74y1Rr5U5sXmvE4bC6EcYe2gWYeDCoYQxhb7k7vFugveY2Zfd/VnC\nHbPZxiNKWwbw5z//mWHDhrWzmIiIpC1fvpxTTjkFtt4jtkMKnVBuBY4xs8ei6dPM7JNARdSj62zC\nMA0JwrDTy8zsl4T++Oeb2fcJfamPJfSuuczCE9mWs3UwumyaAYYNG8bIkSML8sZERHqxnWoq6JX3\noUQD7i24//77lVBERHK0ZMkSjj76aIBx7r5wR9ffne+UFxGRPFJCERGRvFBCERGRvFBCERGRvFBC\nERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGR\nvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBC\nERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGRvFBCERGR\nvFBCERGRvFBCERGRvFBCERGRvCgq5MbNLAH8FpgCbAJOd/f5sfnHA+cDjcA17n6VmRUBvwfGAiXA\nhe5+h5lNAP4AtACvuPuXOtp/U3NLft+QiIi0qdAllBOBUnefAZwLXJqeESWOS4HZwBHAGWY2GDgV\nWO3us4BjgcuiVS4FznP3w4GkmZ3Q0c6/dPEDzF38bh7fjoiItKXQCeUw4G4Ad38KOCA2bzIw193X\nu3sj8CgwC/gLodSSjq8x+nu6uz8S/f1PQiJqV3NzK8tXb9jlNyEiIh0raJUXUA3UxKabzCzp7i1Z\n5tUCfd19A4CZVQE3A9+N5icyl80lgOYWVXuJiHSGQpdQ1gNV8f1FySQ9rzo2rwpYB2Bmo4AHgGvd\n/aZofnO2ZTvS1Ny6E2GLiMiOKnRCeQz4IICZHQzMic17HZhoZv3MrIRQ3fWEmQ0F7gH+y92vjS3/\ngpnNiv4+FniEHDS3KKGIiHSGQld53QocY2aPRdOnmdkngYqoR9fZwL2E6qyr3H2Zmf0S6Aecb2bf\nB1oJCeRbwJVmVkxIRrfkEoCqvEREOkdBE4q7twJnZbz8Zmz+XcBdGet8Hfh6ls3NJfQG2yHNqvIS\nEekUvf7GRpVQREQ6R+9PKCqhiIh0il6fUJpUQhER6RS9PqG0qIQiItIpen1CaVK3YRGRTtHrE0qz\nBogUEekUvT+hqIQiItIplFBERCQven9CUZWXiEin6P0JRSUUEZFO0eHQK2Y2BriK8ATFWcCfgc+7\n+8KCRpYnurFRRKRz5FJC+R1wMeEZJMuBG4DrChlUPunGRhGRzpFLQhnk7vcCCXdvdfcr2fY5Jt2a\nbmwUEekcuSSUjWY2kjCMPGZ2GLC5oFHlkUooIiKdI5fh688G7gQmmNmLwADg4wWNKo/UhiIi0jk6\nTCju/oyZHQhMAlLAG+7eUPDI8kS9vEREOkcuvbx+n/FSq5ltJDw18crunlx0H4qISOfIpQ2lGegL\n3Bb9lAFDCCWW/ytcaPmhEoqISOfIpQ1lqrsfkJ4wszuAp9z942b2UuFCyw8lFBGRzpFLCaXCzIbF\npocQSilQ4GfS76pUKqEqLxGRTpJLQvgB8JyZPU5olD8A+JqZXQDcV8DYdlkykdTzUEREOkkuvbz+\nYmYPADMJ7SlnuPtqM3vI3dcWPMJdkEoldGOjiEgnyaWX1xDgFKASSADTzWycu3+m0MHtqlQyoRsb\nRUQ6SS5tKH8D9gdOBSqADwM94iydSiR0Y6OISCfJdSyvzwJ3EJLLEcA+hQwqX5KpBM0qoYiIdIpc\nEsq70W8Hprh7DVBcuJDyJ5VMqNuwiEgnyaWX1wNmdjPwLeBeM5sGbCpsWPmRSiZV5SUi0kk6LKG4\n+3eB77j728CnCCWVkwodWD4kk6ryEhHpLB0mFDP7q7u/BeDuz7n7L4A/FTyyPChKqVFeRKSztFnl\nZWa3AlOAPcxsfsY6iwsdWD4kkgmaG1VCERHpDO21oXyW8OyTXwFfjb3eBKwoZFD5UpRUCUVEpLO0\nmVDcfT2wHjjBzPYhJJdENHsC8HDhw9s1qaSGXhER6Sy53Cl/GeFmxvlEjwGOfh9VwLjyIplK0NLS\nSmtrK4lEouMVRERkp+XSbfj9gLn7xkIHk2+pKIk0t7RSlFJCEREppFxubJzP1qquHiWZ3JpQRESk\nsHIpoawFXouGr99yQ6O7f75gUeVJKpUAWsMzUYpTXR2OiEivlktCuTv66XHiVV4iIlJYuTwP5Voz\nG0sYEPIeYJS7L8hl42aWAH5LuJ9lE3C6u8+PzT8eOB9oBK5x96ti8w4CLnL3I6Pp/YE7gTejRS53\n95vb238ylQSa1XVYRKQT5HKn/MmEkYZ/Reg6/ISZnZrj9k8ESt19BnAucGlsu0XR9GzCCMZnmNng\naN45wJVAaWxb04FL3P2o6KfdZAJhcEhAw6+IiHSCXBrlvw3MAGrdfSUwlZAccnEYUXWZuz9FeHxw\n2mRgrruvd/dG4FFgVjRvHtuPFzYdOM7MHjKzq8ysoqOdp5Lh7amEIiJSeLkklGZ3r01PuPsycn/A\nVjVQE5tuMrNkG/Nqgb7RPm4l3JEf9xRwjrsfTuh5dkFHO09Fe9JTG0VECi+XRvlXzezLQHHUjvGf\nwIs5bn89UBWbTrp7S2xedWxeFbCunW3dFj2LBeBW4Ncd7TyVUglFRKSz5FJC+RKwB7ARuJpQqvjP\nHLf/GPBBADM7GJgTm/c6MNHM+plZCaG664mM9eP3v9xjZukqs6OB5zraue5DERHpPLmUUDYBT7j7\nuWY2iDAMS12O278VOMbMHoumTzOzTwIV7n6VmZ0N3EtIHFdF1Wlx8UxwJnCZmTUAy4EzOtr5lm7D\nzaryEhEptFwSylWEkszfo+kjgYOAL3a0oru3AmdlvPxmbP5dwF1trPs2oTNAevpFQiN/zrZUeamE\nIiJScLkklAPcfV8Ad18NfNrMXi5sWPmxpcpLbSgiIgWXSxtK0syGpyfMbAi59/LqUuleXroPRUSk\n8HIpoVwIvGBmjxLaOt4LfK2gUeVJUvehiIh0mlwSyivANOAQwhApX87SeN4tpYes130oIiKFl0tC\nucndJwN/LXQw+Za+U76pSQlFRKTQckkor5nZ9wl3qm95yJa7d/tHABcXhYTSqG7DIiIFl0tCGUDo\nKnxk7LUe8QjgoqhVvqFRCUVEpNByGb7+yI6W6a62lFBU5SUiUnAdJhQzG0O4uXEsMBO4Hvi8uy8s\naGR5kC6hNDU1d3EkIiK9Xy73ofwOuJgw3MoK4AbgukIGlS/FRaGXV4NKKCIiBZdLQhnk7vdCGErF\n3a9k21GCu63iVHiOvKq8REQKL5eEstHMRhIN1GhmhwGbCxpVnhSpDUVEpNPk0svrG4RnuU8ws5eA\n/sDHChpVnmxNKGpDEREptFx6eT1rZgcCkwglGnf3hoJHlgfF0Z3yKqGIiBReh1VeZjYauAV4EngI\n+L2ZDS50YPlQXKQ2FBGRzpJLG8qfgfuAEcA4wpMSry1kUPmSrvJqUJWXiEjB5dKGUu3ul8Wmf2Fm\nnytQPHmlKi8Rkc6TSwnlOTM7NT1hZscBLxQupPwpUpWXiEinyaWE8iHgc2b2O8KDtSoAzOwzQKu7\npwoY3y5J3ymvhCIiUni59PIa0hmBFEKxug2LiHSaXKq8eqxkIvxotGERkcLr1QklkUhQVJTS81BE\nRDpBr04oEKq99MRGEZHCy2X4+n7AKYQHbSXSr7v7jwoYV96UFCVpaFQbiohIoeXSy+tmoAZ4hWiA\nyJ6kuCipKi8RkU6QS0IZ5u7HFDySAikuSlK/qamrwxAR6fVyaUN5wcz2K3gkBVJclKJRVV4iIgWX\nSwnlPYSksgLYRGhHaXX38QWNLE+Ki5K6sVFEpBPkklBOKngUBZRuQ2ltbSWRSHS8goiI7JRcEsoi\n4Ezg6Gj5B4DL2l2jGykpStHaCk3NrVueMS8iIvmXS0L5ObAn8HtCdddphGHsv1HAuPIm/tTG9FAs\nIiKSf7kklPcBU929BcDM7gLm0EMSSkmxBogUEekMuVyyF7Ft4ikCeky3qeKUhrAXEekMuZRQ/gw8\naGY3RNOfBG5oZ/luZeuIw0ooIiKFlMvw9T8xsxeAowglmgvd/a6CR5YnxcV6DLCISGdos8rLzKZF\nv2cB9cAdwO1AbfRaj6ASiohI52ivhHIW8AXgh1nmtRJKLN1ecfqpjXomiohIQbWZUNz9C9GfX3H3\nV+LzzOzgXDZuZgngt8AUwl32p7v7/Nj844HzgUbgGne/KjbvIOAidz8ymp4A/IHwGOJX3P1LucQw\nqF8ZAAuX1TB53IBcVhERkZ3QXpXXoVHV1t/MbKaZzYp+jgKuy3H7JwKl7j4DOBe4NLb9omh6NnAE\ncIaZDY7mnQNcCZTGtnUpcJ67Hw4kzeyEXAJ47z7DAHj85WU5hiwiIjujvSqvY4DDgeFA/NknTcDv\nctz+YcDdAO7+lJkdEJs3GZjr7usBzOxRYBbwV2AeYciXP8aWn+7uj0R//zOK7/aOAhjSvxwb3Z+X\n31pNTd1m+laWdrSKiIjshPaqvC4AMLNPu/sf21quA9WEZ6mkNZlZMrpJMnNeLdA32vetZjamne1u\nWTYXM/Ybji96l6deXc77DmpvsyIisrNyuQ/laTP7FVBJGHolBYxz91x6eq0HqmLT6WSSnlcdm1cF\nrGtnW/FW9Y6W3caM/UZwzZ2v8fjL7yihiIgUSC53yt9EOHlPBV4EhhCe3piLx4APwpaG/Dmxea8D\nE82sn5mVEKq7nshYPz6a4wux7srHAo+Qo2EDK5gwsi8vzV3FmpqNua4mIiI7IJeEknT3HxDaQp4n\nNLQflOP2bwU2m9ljwCXAN8zsk2Z2urs3AWcD9xISz1XuntlyHn/k8LeAH0XbKgZuyTEGAI49ZCxN\nza1c/KfnWL6mfkdWFRGRHORS5bXBzEqBNwkN44+aWZ9cNu7urYT7WeLejM2/C8h61727vw3MiE3P\nJfQG2ynvO2gMz76+gidfWc45v36EK86bTVlpLm9fRERykUsJ5U+Eu+TvAr5iZv8ElhY0qgJIJBJ8\n+zMHcvjUkayr28yi5eu7OiQRkV6lw4Ti7pcB/+HuqwglhCsI1V49TlEqyb4TBwKwaHltF0cjItK7\ntFnnY2bfz5iOT+7Ltvem9Bijh4aOZYtWKKGIiORTe40IvfJ5uaOGhV7MKqGIiORXezc2ZhsUsser\nLCtmYN8+akMREcmzDrs5mVkL23bfBXjH3UcVJqTCG79HX555bQWvLVjD3uMGdnU4IiK9Qi6N8kl3\nT7l7CugDfAK4ueCRFdDHjpoEwNV/z/X+TBER6Ugu3Ya3cPdGd7+ZHvIslLZMHjeA/ScN5s1F66ip\n29zV4Yg4AwfdAAAaNElEQVSI9Aq5VHl9JjaZAPYBGgoWUSex0f158c1VvLW0hmk2pKvDERHp8XK5\nVfzI2N+twGrg5MKE03nG7xEGK56vhCIikhcdJhR3Py16GNZ+hGehzImGVOnR0gllwdKaDpYUEZFc\ndNiGYmazgUWEO+SvBeab2YGFDqzQhg4op6JPEfOWrKO1tcfnRxGRLpdLo/wvgWPd/QB3nwp8DLi8\nsGEVXiKRYN+Jg3hndT0vz13d1eGIiPR4uSSUze7+UnrC3Z+ll9xFf/LsMJzMjf/yLo5ERKTny6VR\n/ikzuwq4ktCG8glgYfphV+7+cAHjK6iJo/oxdng1cxeHaq9EolfkSRGRLpFLQpkc/b4o4/UfEnp9\n9eh7Ugb1K2PhsvXUb2qisqy4q8MREemxcunldSSAmVUBKXfP+VnuPcHAvuFZYWtqNiqhiIjsglx6\neY03s6eBhYQeXi+Y2Z4Fj6yTDOxbBsCamk1dHImISM+WS6P874Cfu/tAdx8A/JTQntIrDIpKKGtr\nNnZxJCIiPVsuCWWQu9+SnnD3vwADChdS51IJRUQkP3LqNmxm09ITZjYd2FC4kDrX1jYUJRQRkV2R\nSy+vrwN/NbO1hPtPBtALxvJKU0IREcmPXJ6H8iQwCfgM8Flgkrs/VejAOktFWTElxSlemb+an//x\nWV6au6qrQxIR6ZHaLKGY2QjgMmBP4FHg3N7WZRjCECyzDxzF/c8u5pEXl7JkZS2//uaRHa8oIiLb\naK+Ecg3wBnAO4UmNv+iUiLrAWf8xhb9ceBwTR/Vj8YpaGptaujokEZEep702lD3c/f0AZnY/8GLn\nhNQ1kskE40f0Zd7idSxdVcfY4dVdHZKISI/SXglly1MZ3b2RXvCUxo6MGxGSyMJ39IwUEZEdtSPP\nlO/1Dw1Jl0oWvLO+iyMREel52qvy2sfM5sem94imE0Cru48vbGidb+yIviSTCe57+m1aWlupKCtm\n1tQ9GDGosqtDExHp9tpLKJM6LYpuorKsmC9/dAq//evL3PbQWwDc9uA8jj5wNCfMmsCQAeVdHKGI\nSPfVZkJx97c7M5Du4piDxnDg3sNY+e4G5i1Zx9V/f5W/PzKfZ15fwSeOmcR79x5GZXlJV4cpItLt\n5HKn/G6nX1Up/apKmTS6P0dOH8UN9zq3PjiPX9zwAn1KUvSv7kNxUZKq8hIOnzaSaTaEwf3KSCb1\ngC4R2X0poXSgrLSI0z60N9P3GsIbb6/lweeWsGFTE3UbGliyopZX568BoLgoybCBFYwbUc3M/fdg\n+l5DKS7akT4PIiI9mxJKDhKJBFP2HMyUPQdveQ49wOp1G3n0paW8uWgdy9bUs2xVHYtX1PLwC0sZ\nMaiCC886lEH9yrowchGRzqOEsgsG9SvjxMMnbplubW3lraU1/OOxBdz39CLO+tn9nDBrAqceO7md\nrYiI9A6qk8mjRCLBxJH9+MrH9+f0E95DeZ9ibvrXmzz43GLW1GyktbXX38ojIruxgpZQzCwB/BaY\nAmwCTnf3+bH5xwPnA43ANe5+VVvrmNn+wJ3Am9Hql7v7zYWMf2clEglOmDWB94wfyNm/fIhLrn8e\ngMH9yzj2kLEcd+g4yvvo+fUi0rsUusrrRKDU3WeY2UHApdFrmFlRND0d2Ag8Zma3A4e1sc504BJ3\n7zGDVE4Y2Y8fnH4Iz72xglXrNvLim6u47h+vc+uDb3HSERM4YPJQxgyrVu8wEekVCp1QDgPuBnD3\np8zsgNi8ycBcd18PYGaPAIcDh2SsMz1afjowycxOBOYCX3P3+gLHv8um7TWEaXsNAaB+YyN3PDqf\n2x56i+v+8TrX/eN1Jo3ux6ePncyY4dVU9AnPZhER6YkKnVCqgfhIi01mlnT3lizz6oC+QFXG681m\nlgSeAq509xfM7DzgAsLQ+j1GRVkxnzjG+NBh4/n3s4t5ed4qnnxlOef/7gkgjHg8fa8hHDFtJMMG\nVtCvspT+1aUUFynJiEj3V+iEsp6QINLSySQ9Lz5GfBXwblvrmNlt7p5ONLcCvy5QzAVXWVbM8TPH\nc/zM8fjba3ng2cXU1DewbFU9z7y2gmdeW7Fl2eKiJNNsCIfsO5z37jOMKt2lLyLdVKETymPAh4Bb\nzOxgYE5s3uvARDPrB2wAZgIXR/OyrXOPmX3Z3Z8FjgaeK3DsncLGDMDGDNgyveCdGl7wlbxbu5l1\ndZt5a0kNT726nKdeXU4ymWDPkf0Yv0dfjpg+kr3HDezCyEVEtlXohHIrcIyZPRZNn2ZmnwQqoh5d\nZwP3EkYwvtrdl5nZdutEv88ELjOzBmA5cEaBY+8S40b0ZdyIvtu8tmRlLU/MWcYTc5Yxd8k6fNG7\n3PPkQmZNHcm4EX2ZNXUP3UApIl0u0RvvjTCzscCC+++/n5EjR3Z1OHnV1NzCnHmr+fVNL7C6ZhMA\n1RUlfHjmeKorSpg+eShD+mtUZBHZcUuWLOHoo48GGOfuC3d0fd0p38MUpZJMtSFccd5sVqzdwNOv\nLufaf7zOn+5+AwhtLkdMG8nRB45m73EDSCTUJVlEOocSSg9VXJRi5JAqRg6pYsZ+I1i6qo7lq+v5\n20Nvcd/Ti7jv6UUMH1TBkdNHsdeY/uy352BSut9FRApICaUXGDawgmEDKwA4dsY45sxbzb+eXcTj\nLy/j+ntCyWWvMf357mkH0a+qtCtDFZFeTAmll0kmE0yZNJgpkwZz1kcaecFX8fCLS3j85WWc/7vH\n+fjsSUweO0CN+CKSd0oovVh5n2IOnTKCQ/Ydzo+veYpnXlvBz//4LMlkguMPG8/Hjt6TvpUqsYhI\nfiih7AaSyQTnfvZAnn5tBSvWbODuJxZy+8NvcddjC9hv4iBGD6ti1NAqRg6pZI/BlVRXlKgxX0R2\nmBLKbqK4KMWh+40A4LjDxnHPkwu598m3ed5X8ryv3GbZ4QMr2G/PQYwYVMF+EwczYWRfJRgR6ZAS\nym6otDjFh2dO4MMzJ1C3oYHFK+pYvLKWpSvrWLKyjpfnreKeJ9/esvzY4dWMHVHN7ANGM2XS4C6M\nXES6MyWU3VxleQmTxw1g8ritw780NjXzzup6Fi2v5ZEXl/L0q8tZuGw9Dz63hP5VpUzbawj7jBvI\noVNG6LkuIrKFEopsp7goxZhh1YwZVs3M/fdg0+YmFryznhvvcxYuq+H+ZxZz/zOLufL2OQwfVMmR\n00dyzHvHUFGm5CKyO1NCkQ71KS1i8rgB/PCMQ2hpaWXeknW84Cv593OLWbyilqv//ip//OcbTNij\nL/2qStlzVD8G9i1jYHUfRg2rYkB1n65+CyLSCZRQZIckkwkmje7PpNH9OfkYY319A3c/sZAHn1/C\n6wvXAvDEnGVbl0/APuMHMXJoJZVlxQwdUMHEkX0ZPaya4qJkF70LESkEJRTZJdUVJXx89iQ+PnsS\nLS2trK7ZyIKlNayt3czamk087yuY89Zq5ry1epv1ilIJqitKef/BY/jIERPpU6qPokhPp2+x5E0y\nmWBI//JtRjs+5QN7UbexkTXrNlK7oYGlq+qYt6SGBUtreGd1HTfc69xwr1NVXsy4EX3Zc1Q/Rg+r\nZvTQKkYMrqCstEhdlkV6CCUUKbjKsmIqowb790wYxPuj12s3NHDjfc7i5bWsWLuBl+et5uV525dk\nqspLGDqgnGEDKxg+KPzYmP6MGFTZye9ERNqjhCJdpqq8hC+csO+W6fqNjcxfWsOiFbUsWVHLsjX1\n1G9sZF3dZuYuXscbb7+7ZdlkAiaPG0hpcYpUKsGA6j6MHlrF6GFVjBlWTb+qUpVsRDqZEop0GxVl\nxew7cRD7Thy03bzm5hZWrdvI8jX1LF1Zx91Pvs2r89e0ua3KsmIqy4spLkpSWpxi5NAqykqLKC1O\nsdfYAfSvKqUiKjmV9ymmT0lKCUhkFymhSI+QSiW3DNO//6QhHHfYeJpbWmlubqGxqYWV725g0fJa\nFq2oZdHy9SxZWcfGzU3hZ1MT85bUbN3YQ29tt/1kMkFpcYrS4hR9K0voW1lKZXkx5aXF9K0sYczw\nagZU9aG6soTB/cqoKCtWAhLJoIQiPVYqmSCVTFFSnGJcWV/GjeibdbnmllZWrt1AQ1Mz6+samLv4\nXeo2NlK3sZH66GfDpiY2NzSzqaGJ1TWbeHt5bbv77lOSoqqihKqyEkYPq2L8Hn0ZM6yastIiKsqK\nGDG4kqKUukXL7kUJRXq9VDLB8EEVW6azVallamxqYcOmRuo3NfLu+s3MX1pD3cZG1tVuYvW6Tayu\n2UjdhgaWrq5j/js1PPj8km3WL0olGTawnLLSIvqUFFFakqJPSYo+JUXhd2kR5X2KqK4oobK8hOry\nEqoqShg2oFxdqKXH0idXJIvioiR9K0vpW1nKiEGV7DN+YNblWlpaWb6mnreW1LBkVR2NTc2sq93M\nwmXrWb5mQ1QyatmhfQ8ZUM7IIZWURYloxKCKkHjKSqgoL6a8T2gLKi5K0qekKCSt0iI94lm6nBKK\nyC5IJhOMGFzJiMFtd2Fubmllc0O6Si1Uq23a3Ez9pkZqNzRQW9/A+g0NrK8L9+ksXlHL82+sbHN7\nbSkrTVFWWkRZaUg6ZVEpqLxP6KBQXV5Cn9JQQiqNSkrlfYoYUN2HAdV9dM+P7DIlFJECSyUTlPcp\n3qGRmTdubqKhsZn6jY0sW1NP7YZG6jc0UBe19zQ0NdPY1LK148HmJjZsCh0Q6jc2surdHS8ZhdEL\nSqiuKKWqvCT6O1TFVW/3U0p1RYl6x8k2lFBEuqFQ0igKVW7tlH7a09QcEk79xkbqNoTS0KaGpqiU\n1MzmhibqNjaytmYT79Zu3lJSWrVuIwuXrc9pH8VFSarKS7Z0we5XVUq/ytItJaR0dVxZ7Cc+r6y0\niOKipJJSL6GEItJLFaXCyb6qvASyNwG1qbm5hdoNjayv38z6+gZqNzSwvn7bn5q6zVv+Xle7iaUr\na2lp3Zk4E1RE9wOVFCVJpZIUpRIUpZKUFKcYNbSKflE37r6VIWH1rdyaxIqLUju+UykIJRQR2U4q\nlQyljarSnNdpaWmldkNINPGquPS9QBsyX4u9nu6+XVvfQnNLC03N4R6jllZ48c1V7e63pChJRVnx\nNj9l6bai4tQ2paHK8uItSba0JHQ5T9/8mv6dUnfvnaaEIiJ5kUwmtvSMy5eNm5tYvKKWug2NrI+S\nVfhp2JKE6jaF3+vrG1i2up7mnSkmxaSSCUqKk1GySVFanGTPUf058fAJDB+kAUvbo4QiIt1WWWkR\nk0b3z3n51tbWLT3ptvSqS3da2NwUElP9Zuo3NrK5oZmGphYaGptpaGqmoTH6uzH+egvr6xt58Pkl\nW+41KkolGT6ogoHVfSgtSVFakqKirJj+laW87+AxDOxbVqjD0e0poYhIr5FIJLZUb+VLS0srj738\nDi/NXcWamk3U1G1mycpaFq/YfjSFOx5dwEH7DKN/dWjrGTKgnOEDKxgxuHK3eKCcEoqISDuSyQQz\n99+Dmfvvsc3rjU0tbG4MJaD6jY288OYqrr/ndf71zKLttlGUSjC4fzkVfYoYMqAcGz2AstJwP9Be\nY3vPoxiUUEREdkJxUZLioiSVZcUM6lfGmOHVfOiwcSxbXU9N3WbW1W1m5doNLF1Vz9vL17Ni7QbW\n1Gxi3pIaHn952TbbOXn2JPYeP5CRQyrpX9WnC9/VrlFCERHJk6JUklFDqxg1tKrNZZasrGXJyjo2\nNzSzrm4zN93n/OnuN7bMHzU0JJWiVJK9xvTnkP1GMHZ4dWeEv8uUUEREOtHIIVWMHLI14Rx9wCge\nn7OM5Wvqmbd4Ha8tXMviFXUAPO8ruf5eZ9peQ/j6yVPpX929Sy9KKCIiXaiyvIT3HTRmy3Rrayst\nraHL9AtvrOTuJxfy/Bsr+dkfn+XHZ87o1o9F6L6RiYjshhKJBKlkgsqyYmZO3YMfnzmDQ/cbwavz\n1/DNXz5M3YaGrg6xTUooIiLdWCKR4GufmMqhU0Yw/50aHn3pna4OqU1KKCIi3VxZaREnz54EgL/9\nbhdH07aCtqGYWQL4LTAF2ASc7u7zY/OPB84HGoFr3P2qttYxswnAH4AW4BV3/1IhYxcR6U5GR4+Y\nfuPttV0dSpsKXUI5ESh19xnAucCl6RlmVhRNzwaOAM4ws8HtrHMpcJ67Hw4kzeyEAscuItJtpJIJ\nJo3ux5KVddR203aUQieUw4C7Adz9KeCA2LzJwFx3X+/ujcAjwOFZ1pkeLT/d3R+J/v4nIRGJiOw2\n9hk/CIArbptDa+uuDYJZCIXuNlwN1MSmm8ws6e4tWebVAX2BqozXm80sBcSH96yNlm1LCmD58uW7\nELqISPdy4IRSHu3fzH2PvMz79u+b9/tSYufMnXrITKETynpCgkhLJ5P0vPjtn1XAu22s02xmLRnL\nrmtnv8MBTjnllJ2NW0SkW/voAxcVcvPDgbd2dKVCJ5THgA8Bt5jZwcCc2LzXgYlm1g/YAMwELo7m\nZVvneTOb5e4PA8cCD7Sz32ei7S0DmvP1ZkREerkUIZk8szMrJwpZDxfrsbVf9NJphDaRiqhH13HA\nDwjVWVe7+/9lW8fd3zSzPYErgWJCMvqCu3e/SkQRkd1UQROKiIjsPnRjo4iI5IUSioiI5IUSioiI\n5EWvG76+o+Feuisze46t998sAH5CDxhqxswOAi5y9yPbGh7HzL4AnEEYYudCd7+rq+LNlBH//sCd\nwJvR7Mvd/ebuGH800sTvgbFACXAh8Bo94Pi3Efties6xTxI6CBnhWJ8JbKYHHHtoM/4S8nD8e12j\nvJmdBBzv7p+PThbnuvuJXR1Xe8ysFHjc3afHXrsd+B93f8TMLgfudvfbuyzILMzsHODTQJ27z8gW\nM/AkcB8wDSgHHiWMetDYVXGnZYn//wHV7v6L2DJD6Ybxm9nngP3c/eyo6/1LwIv0gOOfEXt/Qtw/\nBPr2kGN/AuEcc7qZHQ58g9BTtdsfe2gz/jvIw2e/15VQyBi6xcwO6GD57mAKUGFm9xD6gX8XmJYx\n1MwxQLdKKMA84CTgj9F05vA47yNcAT3q7k3AejObS+gS/lxnB5vFdvEDk8zsRMKV2jeA99I94/8L\ncHP0dwpoYvvPTHc9/vHYk4Sr3+nAXj3h2Lv77WZ2RzQ5hnBD9uwecuwz4x9LiH86YLt6/HtjG0rW\n4V66KpgcbQAudvf3A2cBf2bHhprpEu5+K+FElpYZczXbD6WTHmKny2WJ/yngnGgA0vmEe6TaGiKo\nS7n7BnevN7Mqwsn5u/SQ458l9u8BTwPf6gnHHsDdW8zsD8CvgevpIcc+LRb/rwjnm6fIw/Hv7ifa\nndHecC/d1ZuEfyruPhdYAwyNze9oqJnuItvwONmG2Omu7+U2d38h/TewP+EL1S3jN7NRhBEjrnX3\nG+lBxz9L7D3q2AO4++eAScBVQFlsVrc+9mkZ8d+bj+PfGxPKY8AHAbIM99JdfR64BMDMRhD+ifdG\n9ZsQhpp5pI11u5PnzWxW9Hc65meAw8ysxMz6AnsBr3RVgB24J1ZFejShaN8t44/qt+8B/svdr41e\nfqEnHP82Yu9Jx/5UM/tONLmJMLzTs1m+rz0l/hbgb2Z2YPTaTh//3tiGcitwjJk9Fk2f1pXB5Ohq\n4Boze4Twz/0coZRylZmlh5q5pevCy9m3gCvjMbt7q5n9mtCglyA806Z7PswhVDf+r5k1AMuBM9y9\nrpvGfy7QDzjfzL4PtAJfI8Tf3Y9/tti/Afyyhxz7vxG+rw8RzqFfBd4g4/vaTY89bB//1wi97C7b\n1ePf63p5iYhI1+iNVV4iItIFlFBERCQvlFBERCQvlFBERCQvlFBERCQvlFBERCQveuN9KCLbMLMx\nhNEIXiX0p+8DvAx8xd1Xmtl04IvufkYb648Fvufup2eZ90Wg1d2vaGPdA4H/cPfvZJufb2a2ADjc\n3Rd1xv5E4pRQZHex1N2npSfM7CeEm0VnuftzhCG62zIWGJ9thrv/roP97g0M2bFQd4luLJMuoxsb\npdeLSij/dvfxsdeKgRXALGAgcEH0TJSzgc8QhtN42t3PMrOXgHHAtYQkdDGhpPMKsBDA3X9oZp8i\nDNLYQhi24r8IQ1hUAJe4+09j+/8s8AFgACFZ3ePuX46G77jA3Y+MlrsG+DfwEGGMpfnAvsCzwIOE\nURX6ASe5u0cllAcJI1hvBM509zlmNgT4HTAyiu9cd3/AzH4AHAyMAi5z9//btaMtuzO1ochuKXqm\nw1zC+EQArWaWAr5DGMr7AKDFzIYThtZ41t2/Ei07ETjS3U+LrTsCuJQwjPm+hCHlZwDnA3+PJ5OY\nQwjD5+8HfNjM9klvr42w9wN+6O6TgAOBMe4+A7iRbUtYHpXGfkxIghBGlb3a3Q8ETgCuMLOKaF6p\nu79HyUR2lRKK7M5aCVfxALh7M2Fw0WcJw3f/xt2XZVnP3b0u47VDCM+OWBYt8Fl3/zvbDmue6fFo\nKPeNhJLHgA7iXebuL0d/LwHuj/5+G+gfW+7qKIZ/AqPNrBqYDfzIzF4gPK8jBUyIln+qg/2K5EQJ\nRXZLZlZCeATqa/HX3f0kwiNRIYyAOzPL6huzvNZILHmY2SAzG9RBGJtif7dG66d/pxXH/s4cmK+J\n7DJfbyQkkKPcfaq7TyWUntIjx2Z7PyI7TAlFdhfxk32C8MjZJ9x9Qez1QWb2OjDH3S8A7iVUMzWx\n7Yk9m2eA90ZtFQC/AD6c47pxq4Hx0ZDhA4B4QmuvtBN3Cmx5HPYbUQnofiD9nPO9Cb3cytrcgshO\nUEKR3cVwM3s+qvJ5ERgOfCq+gLuvBv6P8GyLZwiN3X8gDEfe18yupQ1RVdfXCM+xeZnwFM5rCE8i\nPCjqVdae1mg7rwH/IHRxvgl4OHOZLH9nbmdS9D6/Dnw2ev2rwMFRB4MbgFPcvb6DmER2iHp5iYhI\nXqiEIiIieaGEIiIieaGEIiIieaGEIiIieaGEIiIieaGEIiIieaGEIiIieaGEIiIiefH/ARNrQ8BH\nnItmAAAAAElFTkSuQmCC\n" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# UK Population data pulled from\n", "# https://en.wikipedia.org/wiki/List_of_English_districts_by_population\n", "populations = pd.read_csv('../data/uk_population.csv', header=None, index_col=0)\n", "header = ['District', 'Population', 'Type', 'Ceremonial_County', 'Historic_County', 'English_Region']\n", "populations.columns = header\n", "\n", "# Convert population to numbers\n", "populations['Population'] = pd.to_numeric(populations['Population'].str.replace(',', ''))\n", "\n", "# Now, turn these populations into percentages\n", "n_areas = populations.shape[0]\n", "total_population = populations['Population'].sum()\n", "populations['percent'] = populations['Population'].astype(float) / total_population\n", "plt.plot(populations['percent'])\n", "plt.title('Distribution of (sorted) region sizes within the UK', fontsize=15)\n", "plt.xlabel('District number')\n", "plt.ylabel('Population percentage')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we'll run the simulation. On each iteration, all the districts vote totally randomly. Then, the difference between \"yes\" and \"no\" is calculated for each. Finally, these differences are combined in a weighted average, where the weights are proportional to the district sizes. This means that bigger districts have a larger influence on the outcome, mimicking the way that the UK tallies votes." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "diff = np.zeros(n_votes)\n", "for ii in tqdm(range(n_votes)):\n", " # Define a random split yes vs. no for each area\n", " area_splits = np.random.rand(n_areas)\n", " area_splits = np.vstack([area_splits, 1 - area_splits])\n", " yes, no = area_splits\n", " \n", " # Now, calculate the difference and average these together, weighted by the area size\n", " diffs = yes - no\n", " diffs = np.average(diffs, weights=populations['percent'])\n", " diff[ii] = diffs * 100" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Here we calculate 99% confidence interval on the difference\n", "clo, chi = np.percentile(diff, [.5, 99.5])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAFJCAYAAACl2EgGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//FXCILmG4gIoiKtwXz1o1+ttooLitUqWmxR\nBBVR6wISlyq/0srXBW1dWfoVLaC1LqC4RKqAIIo7VkXESm2pqPSgYCTUsiokECRk+f1x7oQhTLaZ\nyZ0k834+HnlMcu+Zc89dMvcz55x7TkZVVRUiIiIiEo42qS6AiIiISDpR8CUiIiISIgVfIiIiIiFS\n8CUiIiISIgVfIiIiIiFS8CUiIiISorapLoA0D2Z2G3BbjcVVwFbga+AvwH3OOVfjfQcAXwKznXMD\n49huT6CTc+6NRpTxHOfcnES33YDtnQ5845z7KPj7ZPxxmOCc+02yt5dsZpYJjAN+AewJOOfcD0Mu\nQ4PPb2tmZkcC/wCmOueGpqgMjwOXAT90zn2cgu3PBs4Gcp1zK+tINxW4lBSVs6mZWSWw2Dl3VKrL\nIqmj4EuiVQEvAIuDv9sAHYEjgWHAL8zsfOfcy1Hv2QjcDvyrsRszs58Bc4DfAA25Ob8dlLHR22os\nM7sG+CNwDvBRsLgQv68fNPX2k2QYcD3+eD0OrA1z43GcX2las/BfVlanaPtVwU+y0rVUt5O6cyDN\nhIIvqWm2c+7JmgvNrC8wG/izmf3QObcCwDm3Cbgzzm3tDWQ0NLFz7h3gnTi31VhdqXEDcM59Rfz7\nmgo/wu/Dtc65v6Rg+406v9K0nHNz8MGwpJBzriV9hkgTUZ8vaRDn3KvAb4Hs4DUZMmi+N+fmWq7G\n2D143ZCi7Tfn8ysikjKq+ZLGeAC4AzjXzK5wzlXG6ncV9DW6FRgI5AHbgEXA/znn3grSRPqfVAET\nzOwPQI/g5y/AL4GTgf74ps1zgTOI6vMVXTAzGxCU7SBgJb6ZbbxzrjwqTcy+FmZ2WZB+hHNukpn9\nJdh2FTDbzKqcc5m19fkys4OCcvUBOgFFwExgtHOuOCrdVHxflr2AsfgmzT2BT4ExzrnnG3ISgr5o\nNwDHArsBS4HJwEPOuaqoc0KwD4vNrAr4iXPu3Rj5zQH6Aeac+7zGusHAM8ANzrnxDd3f2s5vpK+P\nmZ0K3Awcg/8c+hi41zk3s8b284JjdSywL/Af4GXgTufcmnqO01T88T4WeBJ/bX3knOsdrO8HXAsc\njT8PG4EFwO3OuX9G5fM28H3gJOAe/HW4B/A34HdBjWz0dn8A3A30xgefs4CnayljR3b8r3wP+BZ4\nE7gj+lxEXaOnACcAVwbH41/ATc65181sKDASOABYHuQxMyqPyPH4oXPu46jrvDZvO+dOjXr/+cCv\ngR8Alfj/6budc2/X2Kc2+ObuK4Lj9jm+qa2xOgfX0Tn4ioJ3gN9Gzo2Z9QbeBZ52zl1a881mthzI\ndM7l1raB4DPhCWAZ/n8K/Pmf2Jh9DtJehb+e8vDNig8Ba4CpwCmR/71Yn0NxXAd9gKPw18H3gVXA\nY8A451xlbfsrzYdqvqTBnHNbgb8D/wXU1XH7AfzNeQNwP/As/gb4mpn9OEgzC9+MCfAq/sN5Y1Qe\nt+FvipPwN7m/B8tj9QU5AXgO+AJ4EKgAxuBvuA0Vne/j7Gje/DN13DjM7Dh8R+oLgPfx+7sG+F9g\noZntWWMbVfj+Tz/FH5engf8BnjOzPvUV0syGA6/hj83zwBR8v7w/AgVBskg/vEgA8RA+MC2sJdun\ngtdBMdYNxt90Chq5v7WeXzMbhj8Gh+OP70P4JsrpZnZT1L52Ad4CzsQHvfcCnwDXAH8Jgvy6RI73\ni/ib65+C/DCz6/BNcHn44PIP+CC4P/COme1TI59sYD7+Jjw12L8TgVfN7NCoMv8QH8D9FB8kFgCn\nB687XbtmthfwIT5QWRMcy/fxx3aRmR0TY58m4IOBF4FpQXnmmNmEYN17+JtwD3wXgSNjHI+Ix/Hn\npebP8iDdgqiy3om/XvcJ3jcVf92+aWYX1SjjE8Dvge34c1sEzACOj7E/tcnAH7PT8Nf4XPwxXWBm\nPwJwzr2H/5LR38x2j36zmZ0QHIOYQW8NffGB11T8tfpBkEeD9zn4cvEnfG3zI/jzeDfwO+rpvxbn\ndfD7IO938Z+3ewTbu6MB+yvNgGq+pLH+HbzuF2ulmXUA8oF3anxrnoL/gLkWeDd4WrET/lvtq865\nSUG6yFuygSOcc+ui8qitTHsD/88598cg3Sj8h/UFZjbFOTevAftV3TzmnHvSzHoAPwb+XLOWLao8\nbfCBy27Az6Kf6DOzscCN+JqS/BrbKQf+xzn3XZD2LfyNZij+225MQZnuxQdRPwn6oGFme+BvxheY\n2VznXAFwZ5D+CHyNWF1Pjc0BivEf9qOjttcRf8N72zn3n8bsbx3nd3/8zeUz4CTnXCQguwWYB9xl\nZnOcc58F5ekODInuh2hm9+NrRs8AXqljv8Af7/nOuerA0sza4W9U/wKOipyHYN0fgauBs/C1iRFd\n8MHXIOdcRZD20yCfS4BRQbqJ+Bvw6ZEaMTO7HR/MRwd0BMfqIOAu59ztUWXoi79+nzKzQ51z0Tfv\nHsDhzrl/B2lX42sQfxnsyyfB8kX4YGEwO4LwndTSt/N8fED6Jv7mTnDzvwUfuP7cObctar/+Cjxs\nZq855zaY2U+Ai/Hn5Rzn3PYgbeQBlsZ0pC8Bjg36lWJmkYB2Er4WEvz1+Fv8+Zoe9d6Lg201JPjq\nCpwV/SBRI/f5aOD/4QOmM5xzpUHaZ/Dnsb59juc6yAOOdM59GaS9H/8F4wqS1y1EmpBqvqSxtgWv\nHWtZ3wZ/w/tedO1BMFxDHlDzW3JtFkQHXvVYjq/ximxrG/6DMwP/IdxUTgD+G3gmxlAKt+ED1YvN\nbLeo5VXA/dE3fPwNBSC3nu39AsjEN0V8FVkY1Ej+P/z+XtHYnQiO10zgsOhaHGAA0J4dNWrx7G9N\nlwDtgNsigVdUGW4L9u+yYHHkWuoZBH4Ro4D9nHP1BV7gj3fN5txM/JOg+TXOA/gnajPwN+Sa+dwX\nCbwCLwdpcwHMrBs+KHgluinSObcB/6BGdYAfHKPBwFfRN9wg/av483EQO4KMiOcjgVcgUjv1RiTw\nCvw1eM2lgYJau6n42qTBUTf7yDV1QyQICcr5Lb4GJosdtaYX4o/VrZHAK0j7Jxr3lHIVPhjZFJXH\na8DrwAlm9v1g8ZP441r9uWJmbYHzgb875xqyza3sGsQ3Zp8j1+stkcArSPsK9Tzlm8B1MCMSeAVp\nv8J/odkn+HIhzZxqvqSxOgSvm2OtdM5tMrNn8bUWK81sAf6D7SXn3NJGbOfL+pNU+6DGt0Lww0NU\n4ofJaCo/xN8k5tdc4ZwrC2of+gOHAEuiVn9eI+2moFavfT3bi+xLrO19ZmYbiX9/nwaG4M/b7cGy\nwcB3+BsAxL+/0SL9XPoEfaOiRa6tSJP2DHzty3XAYDN7DX8tveyca8ywGTtdS0GwOgOq+6/9D/6L\nweH4Zq4qfIBW07Iaf0cCg8h5OyJ4/YhdvV/jb8M3Fb1XS5nfw/dzPBLftBTxRY10W4LXwhrLI0Fl\nfdeUL4xv4n0Bv+8DnXPfRK2OnLPzzOysGm/tjg9+IufsCHyzf6zatvfx+91QNY8Z+NrzM/DHZaVz\nbkXwGdPXzHKCYK0vvqbyrgZupyjG50dj9rln8LooRt4L8P2zahPvdfB5jLTR12NZHduUZkDBlzRW\nbvC6oo40l+A/iIbgO/SeDPzezP6Gr22I2QxSw9ZGlGmXjtfOuXIz+w7ffNlUIrV/m2pZ/3XwmlVj\n+baaCQP1PRnYkO3l1ZNHbd7Gd9q9ALg96IdyGv5BishDA/Hub7Q98ft5VS3rq/Cd+AmaOnviOyKf\ng6/duBgoCzqPD4+uXanDLtdS0PfwD+wYjuM7fMDwN3bcXGuqed4iN+xI2k7Ba0mM935T4++GHMsM\ndj2WW2KkjVW2Bgtqip7H7/clMf4/I/34bqwli+pzFrxuraXTd81jUJ9YD1REjm30//WT+P535+L7\nu/0C39/szw3cTqzPmsbscxdgS3StV5SvYyyLFu91EOt817wepRlT8CUNFvThOQzfcfqz2tIFTTN/\nAP5gZt3xHY4H4fsPvWhmPWo03yRqz5oLgv5KWez6gR+rqb2uYKEuJfgPuv1rWR/5cE7WUA+RG8/+\nteTZKd5tOf+U5DRgZFAj1Qtf+1MQlSwZ+7sZf5M4MLrptI5yfQXkm9mV+BqGvvigPh//RNjN9eVR\nU9Bk9QpQim9+XAAsC47BIHxzazy+DV5zYqyr+SUg+lzGkuxrpy5/xD+ZOcE590yM9ZvxtVm7N+BJ\num+BA80sM8b/eGO/CO3JrsFmt+A1+v/6OXxfu0FmVoB/cve1RnRbiKUx+1wM5Nayz7V1z4hoTteB\nhEh9vqQxrsIH7M/GqKYHwMxyzWy0mf0cwDm3yjn3uHPuTHzn1f3xnYYheaNYx3oa6ITg9W9Ry8rw\nT2rW9N8xytKQskVmAuhdc4WZZQTLNwP1BhkNtBgf/MTa3n/jH4L4NIH8nw7y7w+chw+yo2czaOz+\nxjqGkY7/u5wzM/tvM7sncu2Y2Vlm9kczy3bOVTnnFjnn7sI/CJHBrv1gGuocfKf43zrnHnNepKz/\nE7zGU3vwD/w+nxhjXc39dfjatmNq6SMXGeokkfNZLzP7JT6Q/Qt+mIpYPsYH4rtMh2Nmx5nZWDOL\n7PNH+PtKrCcbY/2f1qW2/+sq/LEGqgd6noM/Zufgv0w9FeO9jdHYfc7EP4FcU31PeDaL60DCp+BL\nGsT8uEy/xX/LG1tH0q34qvo7ozt+Br93w1eXR6bWiDQZJdpB9Admdl7Utjrgn0KLjOET8S+gh+08\nNMAB+GbSmhpStvfwfXAGmtmZNdbdiR+v59kGNo01xNP4JyVHBU8yAmBmWex4kuyJWt5bL+fcEnxf\nrcH48aSm1yh7Y/c31jF8Gn9eRkc/kGF+2IgH8FMRdQ4WH4IfVuLqGtuK7HthI3Yv2nf44Grf6IVm\ndgT+wYUq/BOdjeL8uGOvAqeaWfVco0Et7G1EBaPOuTL8UBH7U2N4gOApt0HA5865hY0tR0OZ2Sn4\n4SlWABfUUcMzFX+8/hD8b0Xe3wE/lMQN7OgjF7n+xplZdlTawcQOTmqTAdxsUUNImNnF+CFr5sbo\n8/ckvq/TOPxnVKIj+U+l4fv8eJB2tPknjyNpf4IPBmvVHK4DSQ01O0q0DGBA1I09MrfjUfhahlL8\nU1BFtWXgnFtjfsyb3wCfmNlc/M22L75z6Z3OuUhn/chTW780s874poOGlLGm5UCB+YFW1+EfO88F\nxjrnomu+HsUPc/BO8Bj47vgPt4/xtSnR/h1s67dmdhQxxvoKmqkuw99wXzSzF4OynID/xvspOwZu\nTJhz7kszux5/w/y7+YmKN+PHweoBTKul2agxnsY/zVWFH/8qevuN3d9dzq9z7gszuwEYD3xqZi/g\nm6rOxAdbL7JjeIBH8YNI/j64kX2MfwpxEL65Zlyc+/gSvlZvVBCIL8c/UdYvWJ7BjgCwsa7DN2M+\nF5yfVfjrsYJdr90b8MfuxiAQeh84ED/59CZ836UmEQSE0/EBxGvA0CDQ2amMzrk7nHNvm9kkYDj+\nnM3Ff4kagO8n9icXDCDqnPvQzMbja9EWm9lL+EFA++MD98b0SdwzyGMOvna6P74P1PAYaV/D9xH7\nPvBY9BOK8WjkPn9gZg/hWwYWm9kr+GFFzsVf213w5782KbsOJHVU8yXRqvD/8L8Lfm7F94nZCz+2\nzg+Cx59jvS+6iekGfI3FJvxj2Pn4b6OXOeeqv90FH14P4Ps1XMuOJp+6JtaN1Tz4YlDOo/EfgFuA\nYc65W6MTOj8O2K/w/SeuAk7F15D9OsY2nw1+Dgz25YBYZQu+kR6D79zbCz/eUid8TdBx0cMp1KOu\nfY7eh/vxgcrf8DeCy4D1wf4m40P6GfyNosjFGA2/Mftb2/l1zv0B+Dm+6WggPsAqwwfs50dqYIK8\nfowfvPIg/Ln7OT54Or7G0AoN5pz7Gv8wwVv4a+CaIP8J+ABwA75/YrS6rsfo6+FLfCA6Df+FZQi+\nWersGGk3BGnvxd+sr8X3a3sc6Fnji0ODy9CA5eDPyV7B71fja7NvY8f//u+IGi/KOTcCX0O8Eh8M\nXIafbWCIc+666Iydczfi/x834//3D8MP3TC3lrLEUokPWv+J/189GX9tHu+CWRJqbLOCHUOKNKbJ\nsdZj1Jh9xp+7/w3KfRX+f2Qk/lyC/+Iac5tJvA7qWyfNSEZVlc6ViIi0bMGQE92ccz3qTZzc7e4D\nlAVjgNVc9wQ+cNvHObc+zHJJ89agZkfzU4qMc879xPxca1PxEf4nzrlrgzT5+G+w2/FzvM0NqrGf\nxjcVRGo+9NSGiIgkjZmdga+JTcXo7r8A7jGzy93OMzHk4ft8farAS2qqt+bLzP4XX/W62Tl3QtBH\nY7xzbr6Z/Ykdc2G9ge8blIXvmHs0vv9DB+fcnWZ2AdArqMoVERFJSNC/tDd+ENJvgEMa0dSfrDLs\nj++PmIXv6L8c/zDHQPzDJn1jNeFLemtIn68v2Hncm6Odc5ERrl/Bj+F0LPCec648GJDxc/w/Q298\ncBZJW+/EwSIiIg30Nf5BnqXA2WEHXgDOT/d0DL6VpycwAv+A0avACQq8JJZ6mx2dc7OCx/Ejop+G\nKcE/DdeBnUfo3YwfaDB6eSRtncysPf5C/g91PyEiIiLpbTpRE2qbWW6KylGJn5R+dM0VKSyThCcT\nP87iooY+aRvPUBPRY8F0wD+aXczOgVUH/CO2xeyYry2Stj7HEGPuOBEREZFm7CRqn6dzJ/EEX383\nsx8HVamRUcsX4QeYa4efJPQQ4BP8eCU/wz8W/zMaFlT9B6CgoIB99923vrSSgIqKCtasWUPbthru\nTVq3r776itGPLWCPDl0SymdryXpuGXoiBxxwQP2JW4jO5/hxQDfMnp3ikkhjlJeXs88++5CZGWsO\neAnT6tWrufjiiyGIXxoinrvuSODRYCqEpcCMYPDFSfiILwMY5ZwrCzrkP2Fm8/ED1F3UgPwrAPbd\nd1+6d+8eR/GkocrLy8nIyFDwJa3eli1b6LD3gWR3qm0KvYZpu/u/6dKlS6v6Ytg1eOhqt1a0T+mg\nvLycbt266fO7eWlwV6kGnbVgctsTgt8/x089UjPNFGBKjWVb8aNRi4hIM7T2r39NdRFE0o5GuBcR\nEREJkeorRaTZqqiooLCwMOF8iopqnY5URCR0Cr5EpNkqLCxk5IR5ZOV0TSifDauW0rn7oUkqlYhI\nYhR8iUizlpXTNeGO8qWb1iSpNCIiiVOfLxEREZEQKfgSEUljXY87jq7HHZfqYoikFQVfIiIiIiFS\n8CUiItIMVFRUcN5553HTTTc1KH1+fj4bN8Y/l/gDDzzA3XffvcvyWbNm0bNnTwYMGMCAAQM466yz\nuOyyy/j444+r01x11VUsX74cgNtuu40+ffowYcIEFixYwKmnnsr5559PWVlZ3GVr7dThXkREpBmY\nP38+eXl5LFu2jJUrV/L973+/zvQLFy5ssrL07NmThx56aKdtXXXVVTz//PPst99+PPzww9Xrnnvu\nOd5++2322WcfRo0axaBBg7j66qubrGytgWq+REREmoEXXniB3r17c8oppzBjxozq5S+//DJDhgxh\n2LBhXH/99axbt47x48cDcOmll7J69WpOPfVUPv300+r3RP/90EMPcf7559O/f3/OOOMM3nzzzUaX\nrVevXpx++ulMmzZtp/yDOQ3Jz8/nwQcfZN68eUybNo177rmnetsDBw5kwIABXHfddaxbtw6ASy65\nhOHDh9OvXz8KCgrYvHkzN998M+eeey79+/dn3LhxVFZWAnDEEUfwwAMPcOGFF9KnTx+eeOKJ6nI9\n/PDDnHnmmZx11lkMHz6czZs3AzBjxgwGDhzIwIEDGTp0KCtWrGj0Pjcl1XyJiEjaqe0hg9qmW2ps\n+sYqLCxk6dKl3HXXXRx88MGMGDGC/Px81qxZw6OPPsqjjz5Kly5dmDlzJk8//TQjR47k9ddf56mn\nniInJ6fWfL/++ms++OADCgoKaNeuHS+//DKTJk2iT58+jS6jmTF//vydlhUUFHDIIYdUl6OoqIiD\nDz6YIUOGMHv2bJYtW8aMGTNo06YNzz33HLfccguPPPIIADk5Obz00ksAjBo1isMPP5yxY8dSWVnJ\nTTfdxOOPP84VV1xBWVkZe+21F9OmTePTTz/lwgsv5MILL2T+/PnMnj2b6dOnk52dze9//3sKCgo4\n6qijmD17NtOmTaN9+/YsWLCA4cOHM3fu3Ebvc1NR8CUiksY0t2PzMGfOHI4//niys7MxM/bdd1/m\nzJlDu3btOOaYY+jSpQsA5557LuAn1gaoCiZGr023bt0YN24cL7zwAitXrmTx4sWUlpbGVcaMjAx2\n3333mOtilePtt99myZIlDBw4EIDKykq2bdtWvb5nz567pJ0+fToA27Zto02bHY1zp512GgCHHXYY\n27dvZ+vWrSxcuJC+ffuSnZ0NwI033gjAPffcw8qVKxk8eHB1uYqLiykuLqZjx45x7XuyKfgSEZG0\n09igsymD1O+++47XX3+d9u3bc9FFF1FVVUVpaSkvvPACgwcP3iltWVkZq1evplu3boAPiCKv0QHQ\n9u3bAfjss8/45S9/yeWXX07v3r055phjuOOOO+Iq55IlSzj44IMbnL6yspL8/Pzqfdi+fTvFxcXV\n67OysnZKO3HiRA488EAASkpKqvcNoH379jvlXVVVRdu2bXdKU1JSQnFxMZWVlfTv35/rr7++et2a\nNWuaTeAF6vMlIiKSUm+88QZ77rknM2fO5JlnnmHatGkUFBSwdetWSkpK+Mc//sE333wD+BqySLNd\nZmZmdZDVuXNnPvnkEwAWL17M+vXrAVi0aBE/+MEPuPzyyznmmGN48803q/tSNcY777zDu+++u0sw\nWJfevXszffr06n5YEyZM4IYbbqg17dSpUwEfYF5zzTUUFBTETBsJMnv16sUbb7zBli1bALj//vuZ\nOnUqvXv3Zu7cudX9ywoKCrj88ssbXO4wqOZLREQkhV588UUGDRq007Ls7GwGDhzIBx98wFVXXcUN\nN9xARkYGnTt3rg5g+vTpw0UXXcSDDz7I9ddfz+23386zzz7LYYcdxmGHHQZAv379eP311/n5z39O\nu3btOP7449m4cWO9TY8fffQRAwYMAHytWteuXZkyZQp77bVX9bKI6N+jnX/++axdu5YLLriANm3a\nsN9++zFu3LiY77nlllsYM2YMZ511FuXl5Zx44okMGzYsZtrI3yeffDIrVqxg8ODBZGRkcNBBB3HX\nXXeRlZXFsGHDGDp0KG3atCE7O5sHHnigzv0NW0Z97cVhM7Nc4Mt58+bRvXv3VBenVSsvL+frr7+m\nbVvF4NI8LV++nN89viThuR3XFv6drJx9Es5n87f/5s4hPyAvLy+hfEQSVV5eTrdu3fT53QysWrUq\n0ieth3OusCHv0VkTkaSrqKigsLAw4XyKiooSL4yISDOj4EtEkq6wsJCRE+aRldM1oXw2rFpK5+6H\nJqlUEktkCAU99SgSHgVfItIksnK6JtzMV7ppTZJKIyLSfOhpRxEREZEQqeZLRKSBqiork9YPLTc3\nl8zMzKTkJSIti4IvEZEG2lqyjokz15OVszGhfEo3rWX8iNP01KRImlLwJSLSCMnoyyapF5mep6lo\nCAipi64OEZE0lo5POUbGOIyeOzCZKisrNQaX1Ekd7kVEJO20adOGtm3bNslPY4O6yspKnnvuOa6+\n+mquvPJKhgwZwiOPPFI9dVA8Kisrueaaa+jbty8FBQUMGDCgepqfaI899hg333xz3NtJhlNPPZVP\nP/20zjSbN2/msssuC6U8N998M48//niTbkNhuYiISArdd999bNmyhfvuu4+srCy2bdvG3Xffzfjx\n4+MOjFavXs3777/P4sWLycjI4OKLL05yqcO1ceNGlixZkupiJI2CLxERkRRZvXo1b731FjNnzmSP\nPfYAoH379vzmN7+prg3asmULEydO5IsvviAjI4Njjz2WIUOGAHDEEUdw5ZVXsmDBAtatW8ell17K\nueeeS35+PuXl5QwcOJBJkyZx+umn88EHH5Cdnc1dd93FwoUL6dy5M507d6ZDhw6Ar10aPXo0y5Yt\no7y8nF69enHDDTfQpk2bXbZzySWXVNdEPfzww8yePZu2bduSm5vL2LFjyc7OZsaMGTzzzDMA7Lnn\nntx6660ceOCBdR6PWPtz6aWXMmrUKL777jsGDBjA888/z4oVKxgzZgwbN26ksrKSSy65hIEDB/Lh\nhx8yevRo9thjD7777jvy8vI47LDDGDp0KADTpk1j0aJF3HvvvYwePZolS5awZcsWqqqquPvuu/nR\nj36U/JMcg5odRUREUmTZsmXk5uZWB14RnTp1onfv3gDcf//95OTk8Nhjj/Hwww+zfPlypk+fDkBZ\nWRl77bUX06ZNY+LEiYwfP57ddtuNRx55hPbt2zNr1iy+973vVU9GXVBQwMqVK3nllVd47LHH+Prr\nr6u3OWbMGA4//HBmzpzJrFmz+Oabb6qb32pu595776WsrIx58+Yxe/Zspk+fzosvvkj37t0pKChg\n0aJFzJ49m2nTpvH8889zxRVXMHz48HqPR6z9KSsrY+zYsey+++7MmjWLyspKfvWrXzFy5EhmzpzJ\nU089xZQpU/j4448B+OKLL5gwYQKzZ89m0KBBzJo1qzr/WbNmMWjQIP75z3+yfv16nn32WV566SX6\n9+/PI488ksCZbBzVfImIiKRImzZtqKqqqjPNhx9+yP333w/4pyjPPvtsZsyYUb0+mNSZww47jO3b\nt7N169Za81q4cCH9+vUjMzOTPfbYg7PPPhvnHABvv/02S5YsqQ7stm3btlP/tVjbWbhwIX379iU7\nOxuAG2+8EYB77rmHlStXMnjw4Or9Ky4upri4mI4dO9a5v/XtT2FhIStXrmTUqFHVeW/bto3PPvuM\nAw88kH3LLiqqAAAXf0lEQVT33Zd9990XgOOOO46ysjI+/fRTdt99d7799luOP/54AH71q18xbdo0\nVq5cyYcffli9D2FQ8CUiksY0t2NqHXLIIXz11Vds3bp1p9qvdevWcd9993HHHXdQWVm503sqKyt3\nGiqjffv21b9XVVXVGcxlZGTstD56oN+KigomTpxY3TRYUlJSXWNWczuRbbVt23anNCUlJRQXF1NZ\nWUn//v25/vrrq9etWbOm3sCrIftTUVFBx44dd6rR2rBhAx06dGDx4sVkZWXtlP68885j1qxZtGvX\njvPOOw/wgeaYMWMYOnQoffr04cADD+TFF1+st2zJomZHERFJO5EApil+agZLdenSpQt9+vTh//7v\n/ygtLQV29PHac889adeuHcceeyyzZ88GfLPcSy+9xNFHH11v3tFBS+T3k046iRdeeIGysjK2bdvG\nyy+/XJ2md+/eTJ06tXo711xzDQUFBXXm3atXL9544w22bNkC+CbSqVOn0rt3b+bOncu6desA39x5\n+eWXN/i41NS2bdvq49qjRw/at2/PnDlzAPjPf/5Dv379an1icsCAAbz11lu89tprDBw4EID333+f\nU089lcGDB3P44Yczb968Rp23RKnmS0RE0krbtm3p1q1bk2+joUaMGMGTTz7JddddR9u2bSkrK+Ok\nk06qDlauu+46Jk2axNChQykvL+fYY4/loosuAtip1qnm37F+Hzx4MCtXrqRfv3506tSJAw44oDrN\nrbfeypgxYzjrrLMoLy/nxBNPZNiwYXVu5+STT2bFihUMHjyYjIwMDjroIO666y6ysrIYNmwYQ4cO\npU2bNmRnZ/PAAw/E3P/ayhz99957782hhx7Kz372M6ZNm8aDDz7I3XffzeTJk6moqODXv/41P/rR\nj/jwww93yb9Lly4cfvjhVFRUsPfee1cfh5EjR9K/f38yMzPp2bMnr7/+eszyNYWM+tqaw2ZmucCX\n8+bNo3v37qkuTqsWGWhQAwFKsi1fvpzfPb4k4ZHg1xb+naycfVpdPpu//Td3DvlBs5heSM2OLVN5\nebkGcm0mVq1aFemn1sM5V9iQ96jZUURERCRECr5EREREQqT6ShGRNKbmRpHwqeZLREREJEQKvkRE\nRERCpGZHEalWUVFBYWFhwvkUFRUlXhgRkVZKwZeIVCssLGTkhHlk5XRNKJ8Nq5bSufuhSSqViEjr\nouBLRHaSldM14XGsSjetSVJpRERaH/X5EhFJY12PO656oFURCYeCLxEREZEQKfgSERERCZGCLxER\nEZEQKfgSERERCZGCLxEREZEQxTXUhJm1BZ4AcoFyIB+oAKYClcAnzrlrg7T5wJXAdmC0c25uwqUW\nEZGk0NyOIuGLt+brZ0Cmc+5E4C5gDHAfMMo5dzLQxsz6m9k+wHCgF9AXGGtmuyWh3CIiIiItUrzB\n1zKgrZllADn4Wq2jnHPzg/WvAKcDxwLvOefKnXPFwOfAEQmWWURERKTFineE+81AD+BfQGfgLOCk\nqPUlQEegA7Cpxvty4tymiIiISIsXb83Xr4FXnXMGHAk8CbSLWt8B2AgU44OwmstFRERE0lK8NV/f\n4JsawQdTbYF/mNnJzrl3gDOBt4BFwGgzawfsARwCfJJYkUVEWraqykqKioqSkldubi6ZmZlJyUtE\nwhFv8DUBeMzM3gV2A24CPgImBx3qlwIznHNVZjYJeA/IwHfIL0tCuUVEWqytJeuYOHM9WTmJNQSU\nblrL+BGnkZeXF3cekXkd9dSjSHjiCr6cc1uAC2KsOiVG2inAlHi2IyLSWmXldCW70/6pLoaIpIAG\nWRUREREJkYIvERERkRAp+BIREREJkYIvERERkRDF+7SjiIi0AnrKUSR8qvkSERERCZGCLxEREZEQ\nKfgSERERCZGCLxEREZEQKfgSERERCZGCLxGRNNb1uOOq53cUkXAo+BIREREJkYIvERERkRAp+BIR\nEREJkYIvERERkRAp+BIREREJkeZ2FBFJY5rbUSR8qvkSERERCZGCLxEREZEQKfgSERERCZGCLxER\nEZEQKfgSERERCZGCLxGRNKa5HUXCp+BLREREJEQKvkRERERCpEFWRVqBiooKCgsLE86nqKgo8cKI\niEidFHyJtAKFhYWMnDCPrJyuCeWzYdVSOnc/NEmlEhGRWBR8ibQSWTldye60f0J5lG5ak6TSiIhI\nbRR8iYikMc3tKBI+dbgXERERCZGCLxEREZEQKfgSERERCZGCLxEREZEQKfgSERERCZGCLxGRNKa5\nHUXCp+BLREREJEQKvkRERERCpOBLREREJEQKvkRERERCpOBLREREJESa21FEJI1pbkeR8KnmS0RE\nRCRECr5EREREQqTgS0RERCRE6vMlItJCVVVWUlRUlJS8cnNzyczMTEpeIlI3BV8iIi3U1pJ1TJy5\nnqycjQnlU7ppLeNHnEZeXl6SSiYidVHwJSLSgmXldCW70/5xv3/y5HyqKiv4fMRpSSyViNRFfb5E\nREREQhR3zZeZ3QScDewGPAi8C0wFKoFPnHPXBunygSuB7cBo59zcBMssIiIi0mLFVfNlZicDvZxz\nJwCnAN8H7gNGOedOBtqYWX8z2wcYDvQC+gJjzWy3pJRcREREpAWKt9nxp8AnZjYbmAO8BBzlnJsf\nrH8FOB04FnjPOVfunCsGPgeOSLDMIiIiIi1WvM2OXfC1Xf2AA/EBWHQgVwJ0BDoAm6KWbwZy4tym\niIiISIsXb/C1AVjqnCsHlpnZd0D3qPUdgI1AMT4Iq7lcRESagWHDHmXzt//mzlQXRCSNxNvs+B6+\nDxdm1g34L2Be0BcM4ExgPrAI6G1m7cwsBzgE+CSxIouIiIi0XHHVfDnn5prZSWb2IZABXAMUApOD\nDvVLgRnOuSozm4QP1jLwHfLLklN0ERERkZYn7qEmnHM3xVh8Sox0U4Ap8W5HREREpDXRIKsiIiIi\nIVLwJSIiIhIize0oIpLGqud2HDI71UURSRuq+RIREREJkYIvERERkRAp+BIREREJkYIvERERkRCp\nw72ISLqrgqKioqRklZubS2ZmZlLyEmmtFHyJiKSxYcMeZW3h32HmMrJyEpt6t3TTWsaPOI28vLwk\nlU6kdVLwJSIiZOV0JbvT/qkuhkhaUJ8vERERkRAp+BIREREJkYIvERERkRCpz5dIClVUVFBYWJhw\nPsl6Uk1ERJqegi+RFCosLGTkhHlk5XRNKJ8Nq5bSufuhSSqVpJPJk/OpKC/j4gvGpLooImlDwZdI\niiXjKbPSTWuSVBoREWlq6vMlIiIiEiIFXyIiIiIhUvAlIiIiEiIFXyIiIiIhUod7EZE0FpnbMSvV\nBRFJI6r5EhEREQmRgi8RERGRECn4EhEREQmRgi8RERGRECn4EhEREQmRgi8RkTQ2eXI+s16dmOpi\niKQVBV8iIiIiIVLwJSIiIhIiBV8iIiIiIVLwJSIiIhIiBV8iIiIiIdLcjiIiaUxzO4qETzVfIiIi\nIiFS8CUiIiISIgVfIiIiIiFS8CUiIiISIgVfIiIiIiFS8CUiksY0t6NI+BR8iYiIiIRIwZeIiIhI\niBR8iYiIiIRIwZeIiIhIiBR8iYiIiIRIczuKiKQxze0oEj7VfImIiIiESMGXiIiISIgUfImIiIiE\nKKE+X2bWFfgb0AeoAKYClcAnzrlrgzT5wJXAdmC0c25uItsUERERacnirvkys7bAQ0BpsOg+YJRz\n7mSgjZn1N7N9gOFAL6AvMNbMdkuwzCIiIiItViI1X+OBPwE3AxnAUc65+cG6V4Az8LVg7znnyoFi\nM/scOAL4KIHtiohIkkyenE9FeRkXXzAm4byqKispKipKQqkgNzeXzMzMpOQl0tzEFXyZ2eXAWufc\nG2Y2KlgcXYtWAnQEOgCbopZvBnLi2aaIiDRvW0vWMXHmerJyNiaUT+mmtYwfcRp5eXlJKplI8xJv\nzdcQoNLMTgeOBJ4E9o5a3wHYCBTjg7Cay0VEpBXKyulKdqf9U10MkWYtruAr6NcFgJm9BVwN3GNm\nP3bOvQucCbwFLAJGm1k7YA/gEOCThEstkmIVFRUUFhYmnE+ymmhERKTlSOYI9yOBR4MO9UuBGc65\nKjObBLyH7xc2yjlXlsRtiqREYWEhIyfMIyuna0L5bFi1lM7dD01SqUREpCVIOPhyzp0a9ecpMdZP\nAaYkuh2R5iYZzSulm9YkqTQiItJSaG5HEZE0prkdRcKnEe5FREREQqTgS0RERCRECr5EREREQqTg\nS0RERCRECr5EREREQqTgS0QkjU2enM+sVyemuhgiaUXBl4iIiEiIFHyJiIiIhEjBl4iIiEiIFHyJ\niIiIhEjBl4iIiEiINLejiEga09yOIuFTzZeIiIhIiBR8iYiIiIRIwZeIiIhIiBR8iYiIiIRIwZeI\niIhIiBR8iYikMc3tKBI+BV8iIiIiIVLwJSIiIhIiBV8iIiIiIVLwJSIiIhIiBV8iIiIiIdLcjiIi\naUxzO4qETzVfIiIiIiFS8CUiIiISIgVfIiIiIiFS8CUiIiISIgVfIiIiIiFS8CUiksY0t6NI+DTU\nhIiINCtVlZUUFRUlJa/c3FwyMzOTkpdIsij4EhGRZmVryTomzlxPVs7GhPIp3bSW8SNOIy8vL0kl\nE0kOBV8iItLsZOV0JbvT/qkuhkiTUJ8vERERkRAp+BIREREJkZodJa1UVFRQWFiYcD7J6gwskmqa\n21EkfAq+JK0UFhYycsI8snK6JpTPhlVL6dz90CSVSkRE0omCL0k7yejIW7ppTZJKIyIi6UZ9vkRE\nRERCpOBLREREJEQKvkRERERCpOBLRCSNaW5HkfAp+BIREREJkYIvERERkRAp+BIREREJkYIvERER\nkRAp+BIREREJUVwj3JtZW+AxIBdoB4wGPgOmApXAJ865a4O0+cCVwHZgtHNubsKlFhGRpNDcjiLh\ni7fm6xfAeufcj4G+wAPAfcAo59zJQBsz629m+wDDgV5BurFmtlsSyi0iIiLSIsU7t+NzwPTg90yg\nHDjKOTc/WPYKcAa+Fuw951w5UGxmnwNHAB/FX2QRERGRliuu4Ms5VwpgZh3wQdgtwPioJCVAR6AD\nsClq+WYgJ66SioiIiLQCcXe4N7PvAW8BTzjn/oyv5YroAGwEivFBWM3lIiIiImkpruAr6Mv1GnCD\nc+6JYPE/zOzHwe9nAvOBRUBvM2tnZjnAIcAnCZZZREREpMWKt8/XzcCewG/N7HdAFfAr4P6gQ/1S\nYIZzrsrMJgHvARn4DvllSSi3iIgkweTJ+VSUl3HxBWNSXRSRtBFvn68RwIgYq06JkXYKMCWe7YiI\niIi0NvHWfImIiDRrVZWVFBUVJSWv3NxcMjMzk5KXiIIvERFplbaWrGPizPVk5ST2nFfpprWMH3Ea\neXl5SSqZpDsFXyIi0mpl5XQlu9P+qS6GyE40t6OIiIhIiFTzJSKSxjS3o0j4VPMlIiIiEiLVfEmL\nUFFRQWFhYcL5JOvJJxERkXgp+JIWobCwkJET5pGV0zWhfDasWkrn7ocmqVQiIiKNp+BLWoxkPLVU\numlNkkojIiISH/X5EhEREQmRgi8RkTQ2eXI+s16dmOpiiKQVBV8iIiIiIVLwJSIiIhIiBV8iIiIi\nIVLwJSIiIhIiBV8iIiIiIdI4XyIiaUxzO4qETzVfIiIiIiFSzZeIiEgdqiorkzYvbG5uLpmZmUnJ\nS1ouBV8iIiJ12Fqyjokz15OVszGhfEo3rWX8iNPIy8tLUsmkpVLwJSIiUo9kzC0rEqE+XyIiIiIh\nUvAlIpLGNLejSPjU7ChNqqKigsLCwoTzSVZnVxERkVRT8CVNqrCwkJET5pGV0zWhfDasWkrn7ocm\nqVQiIiKpo+BLmlwyOqqWblqTpNKIiIiklvp8iYiIiIRIwZeIiIhIiNTsKCKSxjS3o0j4VPMlIiIi\nEiIFXyIiIiIhUvAlIiIiEiL1+ZKYNDiqiEhyVVVWJu0zsXv37knJR1JDwZfEpMFRRUSSa2vJOibO\nXE9WzsaE8indtJbfDz+F733ve0kqmYRNwZfUSoOjirR+kyfnU1FexsUXjEl1UdJCMj5XpeVTny8R\nERGRECn4EhEREQmRgi8RERGRECn4EhEREQmRgi8RERGREOlpx1amoqKC5cuXNyhteXk5a9euJTMz\nc5d1Gp9LJD1obkeR8Cn4amWWL1/OJTc/o/G5REREmikFX62QxucSEWm9IiPlL1u2jLZtE7uN5+Xl\nxWz9kKal4EtERKQF2VqyjgdmryfrLyUJ5VO6aS1Pjb2Igw8+OEklk4ZS8CUiItLCaKT8lk1PO4qI\niIiESMGXiEgamzw5n1mvTkx1MUTSipodm4nGDBFRly+//DIJpRERkdauqrIyafcMddxvnCYPvsws\nA3gQOBL4DhjmnFvR1NsNSzKDpt89slBDRIiISCi2lqzjd4+sJysnsXvYlo2rueuqE+nRo0fCZUqX\nIC6Mmq9zgPbOuRPM7DjgvmBZXJIV7EByTnKyx9XSEBEiIhKWZA1N5CsPErs3p9PTl2EEX72BVwGc\nc381s56JZJasYCdZkfqXX36pcbVERCStJeM+mE7NoGEEXx2BTVF/l5tZG+dcZS3pMwHOuyifzLbt\ndlm5detm2nT5EeXf7Z5QobZu/JqR4/5M+6w9E8qnZP1KOnTtQfl3xQnlU/rtv6ko26J8lI/yUT6h\n5vN1RhWVGRmUrFvRLMqjfNI3n42rv2DkuI8Svi9vK93Ig3cOTUozaEOsXr068muDo72MqqqqpilN\nwMzuBRY652YEf690zn2/jvS9gflNWigRERGR5DrJOfdeQxKGUfO1AOgHzDCz44El9aRfBJwE/Aeo\naOKyiYiIiCQiE9gPH780SBg1X5GnHY8IFg1xzi1r0o2KiIiINFNNHnyJiIiIyA4a4V5EREQkRAq+\nREREREKk4EtEREQkRM1ubkczGwCc55y7OPj7HGA8sDJIcptzTkNRJFGMY34cMBHYDrzhnLszleVr\nrcxsFRB5+GShc+6WVJantWrtU5w1V2b2ETvGePzSOXdFKsvTmgWf2eOccz8xszxgKlAJfOKcuzal\nhWulahzzHwIvsePz/E/Ouel1vb9ZBV9mNgE4A1gctfho4H+dc7NSU6rWrZZj/hAwwDlXaGZzzexI\n59w/U1PC1in4gPzIOdc/1WVJA0md4kzqZ2btAZxzp6a6LK2dmf0vcAmwOVh0HzDKOTffzP5kZv2d\ncy+kroStT4xjfjRwr3PuDw3No7k1Oy4Arqmx7GhgqJm9a2bjzay5lbml2+mYm1kHoJ1zrjBY9BrQ\nJwXlau2OBrqb2Vtm9pKZtf7JzFJnpynOgISmOJMGORL4LzN7zczeDIJeaRpfAAOi/j46qnXoFfT5\n3RR2OebAz83sHTObbGb/VV8GKQlkzGyomS0xs4+jXo+upZrudWC4c+7HQDZwdbilbR0accw7AtFz\nRJQAOeGVtPWJdezxgwiPCWoGxgJPp7aUrVrMKc5SVZg0UQrc45z7Kf7LXYGOedMIWoXKoxZlRP2u\nz+8mEOOY/xXfQncysAK4vb48UtLs6Jx7DHisgckfd85FPjhfAAY2Talat0Yc82L8zSqiA7CxSQqV\nJmIdezPbg+Cf1zm3wMz2S0XZ0kQx/jqOqGtuWUmOZfjaAZxzn5vZBvwI4P9OaanSQ/S1rc/vcMyO\nilNmAZPqe0NL+CbysZl1C34/DfgolYVp7ZxzJcA2M+sRdFT+KZprsyncBowAMLMjgaLUFqdVWwD8\nDKCBU5xJ4oYC9wIEn98d8LW90vT+bmY/Dn4/E31+h+E1M4t0Z2hQnNKsOtzX4gpglpmVAp8Bj6a4\nPOngauAZfHD+unOuwfNVSYONA542s5/jnyq9PLXFadVmAaeb2YLg7yGpLEyamAI8bmbz8TUxQ1Xb\nGJqRwKNmthuwFJiR4vKkg2uA+82sDFgNXFnfGzS9kIiIiEiIWkKzo4iIiEiroeBLREREJEQKvkRE\nRERCpOBLREREJEQKvkRERERCpOBLREREJEQKvkRERERCpOBLREREJET/Hzr6vyZIW8dbAAAAAElF\nTkSuQmCC\n" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Let's look at the distribution of differences\n", "# This time, voting is randomized for each region\n", "f, ax = plt.subplots(figsize=(10, 5))\n", "_ = ax.hist(diff, bins=30)\n", "_ = ax.axvline(actual_diff, color='r', ls='--')\n", "axfill = ax.fill_between([clo, chi], *ax.get_ylim(), alpha=.1, color='k')\n", "ax.set_title('Distribution of votes randomized by region', fontsize=20)\n", "ax.legend([ax.lines[0], axfill], ['Actual Difference', 'Confidence Interval'], fontsize=12)\n", "_ = plt.setp(ax, xlim=[-15, 15])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we see a different sort of picture. Randomizing votes by district instead of by individual greatly increased the variability in the outcome. So much so that the \"true\" results from the Brexit now fall well within our confidence interval.\n", "\n", "# Concluding thoughts?\n", "So what can we conclude from something like this? The point of this article isn't to say that this particular simulation proves anything about the Brexit vote, but it does bring up an important point: we need to account for randomness whenever we aggregate data that we've collected.\n", "\n", "When deciding whether to make a gigantic decision that will affect millions of people, we should be reasonably certain that the people's opinion is clear. Choosing a 50/50 split as a cutoff means that we could potentially make such a decision because of random chance. Doesn't sound like a great way to conduct national policy to me.\n", "\n", "What could we do instead? There's the hard part. But the short answer is that we could include some idea of random variability in our voting rules. For example, we could require that this kind of \"should we deviate from the norm\" decision exceeds the results expected from a totally random vote. Settling on this uncertainty limit is not a simple task, but then again I wouldn't want to bet 2 trillion dollars worth of global economy on a coin flip." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Extra: Adding a short term swing\n", "As I mentioned above, deciding how to simulate the votes involves making assumptions about how things will go. I tried to keep the simulation as simple as possible in order to make a point, but you could include extra components as well.\n", "\n", "For example, what if we chose a random subset of districts in the UK and swung their vote several percentage points in one random direction? This might happen if an eye-catching event caused sentiment to momentarily swing in one direction or another. In the long run each district's sentiment would probably ease back into it's \"natural\" split, but since votes happen on one day, these short-term factors can play an important role.\n", "\n", "We can build this into our simulation..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Here we define a percentage of districts that undergo a sudden swing in voter opinion\n", "perc_swing_amt = .1\n", "perc_swing_districts = .2\n", "n_perc_swing_districts = int(n_areas * perc_swing_districts)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Now, re-run the simulation including the random swing.\n", "diff = np.zeros(n_votes)\n", "for ii in tqdm(range(n_votes)):\n", " # Define a random split yes vs. no for each area\n", " area_splits = np.random.rand(n_areas)\n", " \n", " # Define a random swing across a subset of random districts\n", " swing = perc_swing_amt * np.random.choice([-1, 1])\n", " ixs_swing = np.random.choice(range(n_areas), n_perc_swing_districts, replace=False)\n", " area_splits[ixs_swing] = np.clip(swing + area_splits[ixs_swing], 0, 100)\n", " \n", " # Now calculate the opposing side amount and average\n", " area_splits = np.vstack([area_splits, 1 - area_splits])\n", " yes, no = area_splits\n", " \n", " # Now, calculate the difference and average these together, weighted by the area size\n", " diffs = yes - no\n", " diffs = np.average(diffs, weights=populations['percent'])\n", " diff[ii] = diffs * 100" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Here we calculate 99% confidence interval on the difference\n", "clo, chi = np.percentile(diff, [.5, 99.5])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAFJCAYAAACoxjhBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//FXCEKJYEAp8hVqo7F8bGu1BUHrT6sV69Kq\niCsubRUFtYpSpS5YtwqVryIiaKkLiha0FiigKIJf0Aq0iqK2LvhRkQhKRUXZwZBMfn+cO2GIk2QS\n5mZ9Px+PPJK598y55y6Z+5lzz5JTVlaGiIiIiGRXi/ougIiIiEhTpCBLREREJAYKskRERERioCBL\nREREJAYKskRERERioCBLREREJAYt67sAUrfM7EbgxgqLy4DNwErgOWCUu3uF930bWAZMd/eTa7Hd\nA4EO7v5sDcp4krs/saPbzmB7PwO+cPfF0evDCcdhtLtfke3tZZuZ5QIjgHOA9oC7+w/ruAwZn9+m\nzMwOAF4DJrh7/3oqw0PAr4Efuvt/6mH704ETgQJ3X15FugnAr6incsbNzBLA6+7evb7LIvVHQVbz\nVAbMAF6PXrcAdgEOAC4AzjGz09z96ZT3rAFuAt6p6cbM7OfAE8AVQCY34eejMtZ4WzVlZhcD9wAn\nAYujxUWEfX0x7u1nyQXAlYTj9RDwaV1uvBbnV+I1jfCl5JN62n5Z9JOtdI3VTdTfOZAGQkFW8zXd\n3R+puNDMjgWmA381sx+6+wcA7r4W+EMtt/VNICfTxO7+D+AftdxWTXWiwge9u39I7fe1PvyIsA+X\nuPtz9bD9Gp1fiZe7P0EIeqUeuXtj+gyRmKhNlmzH3Z8BrgfaRr+zIYeGexNuqOWqiW9Ev1fX0/Yb\n8vkVEak3qsmSdO4GbgZOMbPz3T2Rrl1U1Bbo98DJQCHwFfAycJu7z4vSJNuHlAGjzexOYK/o5zng\nN8DhQB/CI8lTgKNJaZOVWjAz6xuV7TvAcsLjsZHuXpKSJm1bCDP7dZR+sLuPMbPnom2XAdPNrMzd\ncytrk2Vm34nKdRTQAVgBTAWGu/u6lHQTCG1NdgVuJTyKbA+8BfzR3f+eyUmI2opdBfQCdgKWAA8A\nf3b3spRzQrQPr5tZGfBTd38hTX5PAMcD5u7vVVjXD3gUuMrdR2a6v5Wd32RbHDM7ErgW6En4vPkP\ncIe7T62w/cLoWPUCOgP/BZ4G/uDuq6o5ThMIx7sX8Ajh2lrs7odG648HLgF6EM7DGmAhcJO7/zsl\nn+eBPYHDgNsJ12Eb4BXghqiGNXW7PwCGAYcSgsxpwMRKyrgL2/5XvgV8CfwfcHPquUi5Ro8ADgEG\nRsfjHeAad59jZv2BIcC3gaVRHlNT8kgejx+6+39SrvPKPO/uR6a8/zTgt8APgAThf3qYuz9fYZ9a\nEB5Tnx8dt/cIj8hqarfoOjqJ8MX/H8D1yXNjZocCLwAT3f1XFd9sZkuBXHcvqGwD0WfCw8C7hP8p\nCOf/rprsc5T2QsL1VEh4HPhnYBUwATgi+b+X7nOoFtfBUUB3wnWwJ/AR8CAwwt0Tle2vNByqyZKv\ncffNwKvAzkBVDajvJtyEVwNjgccJN7rZZvaTKM00wuNHgGcIH8JrUvK4kXDzG0O4mb0aLU/XVuMQ\n4G/A+8CfgFLgj4Qba6ZS832IbY8l/0oVNwgzO4jQoPkM4J+E/V0F/A74l5m1r7CNMkL7pGMIx2Ui\n8D3gb2Z2VHWFNLNBwGzCsfk7MJ7Qbu4eYFKULNlOLhko/JkQgBZVku1fot+np1nXj3BzmVTD/a30\n/JrZBYRjsB/h+P6Z8Ghxspldk7KvHYF5wHGE4PYO4E3gYuC5KJivSvJ4P0m4iY6L8sPMLiU8Oisk\nBJF3EoLdPsA/zGz3Cvm0BeYTbrYTov37f8AzZvbdlDL/kBCoHUMIBicBP4t+b3ftmtmuwCJCQLIq\nOpb/JBzbl82sZ5p9Gk246T8JPBaV5wkzGx2tW0C42e5FeLR/QJrjkfQQ4bxU/FkapVuYUtY/EK7X\n3aP3TSBct/9nZmdVKOPDwP8CWwnndgUwBTg4zf5UJodwzHoTrvGnCMd0oZn9CMDdFxC+TPQxs2+k\nvtnMDomOQdrgtoJjCQHWBMK1+mKUR8b7HH2JGEeoPb6PcB6HATdQTfuyWl4H/xvl/QLh87ZNtL2b\nM9hfaQBUkyWV+Tj6/T/pVppZO2AA8I8K34LHEz5ILgFeiHoHdiB8S33G3cdE6ZJvaQvs7+6fpeRR\nWZm+CVzm7vdE6YYSPpTPMLPx7j43g/0qf6zl7o+Y2V7AT4C/Vqw1SylPC0KAshPw89QedGZ2K3A1\noeZjQIXtlADfc/ctUdp5hBtKf8K317SiMt1BCJZ+GrURw8zaEG66Z5jZU+4+CfhDlH5/Qg1XVb20\nngDWET7Uh6dsbxfCje15d/9vTfa3ivPbhXATeRs4zN2Tgdd1wFzgFjN7wt3fjsrTFTgvtZ2gmY0l\n1HQeDcyqYr8gHO/57l4eQJpZK8IN6R2ge/I8ROvuAS4CTiDUDiZ1JARZp7t7aZT2rSifXwJDo3R3\nEW60P0vWcJnZTYSgPTVwIzpW3wFucfebUspwLOH6/YuZfdfdU2/SewH7ufvHUdpPCDWCv4n25c1o\n+cuEoKAf24Lt7VTS9vI0QuD5f4SbONFN/jpCgPoLd/8qZb9eAu41s9nuvtrMfgqcTTgvJ7n71iht\nsiNJTRq0rwd6Re0+MbNk4DqGUKsI4Xq8nnC+Jqe89+xoW5kEWZ2AE1I79NRwn3sAlxECo6PdfVOU\n9lHCeaxun2tzHRQCB7j7sijtWMIXifPJXnMOiZFqsqQyX0W/d6lkfQvCje1bqbUB0TAIhUDFb72V\nWZgaYFVjKaEGK7mtrwgfkDmED9u4HALsAzyaZoiCGwkB6dlmtlPK8jJgbOqNnXDjACioZnvnALmE\nRwgfJhdGNYyXEfb3/JruRHS8pgLfT62VAfoCrdlWQ1ab/a3ol0Ar4MZkgJVShhuj/ft1tDh5LR0Y\nBXhJQ4H/cffqAiwIx7viY9hcQs/LARXOA4QerDmEG2/FfEYlA6zI01HaAgAz24Nw85+V+gjR3VcT\nOkyUB/LRMeoHfJh6Y43SP0M4H99hWzCR9PdkgBVJ1jY9mwywIi9FvwvIUFQLN4FQO9Qv5aaevKau\nSgYbUTm/JNSo5LGtFvRMwrH6fTLAitKOo2a9gssIQcfalDxmA3OAQ8xsz2jxI4TjWv65YmYtgdOA\nV909k21u5uvBek32OXm9XpcMsKK0s6imV+0OXAdTkgFWlPZDwheX3aMvEdLAqSZLKtMu+r0h3Up3\nX2tmjxNqIZab2ULCB9hMd19Sg+0sqz5JuRcrfMuDMOxCgjD8RFx+SLgZzK+4wt2Lo9qEPsC+wBsp\nq9+rkHZtVEvXuprtJfcl3fbeNrM11H5/JwLnEc7bTdGyfsAWwgc91H5/UyXboRwVtV1Klby2ko+i\npxBqUy4F+pnZbMK19LS712Q4iu2upSgonQLl7cu+R/gCsB/h8VQZIRCr6N0Kr5MBQPK87R/9XszX\n/bPCayM84llQSZkXENohHkB4JJT0foV0G6PfRRWWJ4PH6q6pUJjwaHYGYd9PdvcvUlYnz9mpZnZC\nhbd2JQQ5yXO2P+Fxfbras38S9jtTFY8ZhNrwownHZbm7fxB9xhxrZvlRUHYsoebxlgy3syLN50dN\n9vnA6PfLafJeSGg/VZnaXgfvpUmbej0WV7FNaQAUZEllCqLfH1SR5peED5zzCA1rDwf+18xeIdQe\npH18UcHmGpTpaw2g3b3EzLYQHjvGJVmbt7aS9Suj33kVln9VMWGkup54mWyvsJo8KvM8ofHsGcBN\nUTuR3oQODcnG+7Xd31TtCft5YSXrywiN6YkeUR5IaBB8EqG24mygOGrEPSi1tqQKX7uWoraBd7Jt\nmIsthMDgFbbdRCuqeN6SN+Zk2g7R7/Vp3vtFhdeZHMscvn4sN6ZJm65sGYtqfv5O2O9fpvn/TLaz\nu7qSLMrPWfR7cyWNryseg+qk69iQPLap/9ePENrHnUJoj3YOoT3YXzPcTrrPmprsc0dgY2otVoqV\naZalqu11kO58V7wepQFTkCVfE7Wx+T6hAfPblaWLHqncCdxpZl0JDX9PJ7TvedLM9qrw2GVHta+4\nIGpPlMfXP9jTPQqvKiioynrCB1qXStYnP4SzNYRC8gbTpZI8O9R2Wx56JT4GDIlqmH5MqM2ZlJIs\nG/u7gXAz2Dv1kWcV5foQGGBmAwk1BscSgvcBhB5Y11aXR0XRo6ZZwCbCY8OFwLvRMTid8Ji0Nr6M\nfuenWVcx2E89l+lk+9qpyj2EnpCj3f3RNOs3EGqnvpFBz7Uvgb3NLDfN/3hNv/C05+tB5R7R79T/\n678R2sKdbmaTCD1lZ9eguUE6NdnndUBBJftcWbOKpIZ0HUgdUpssSedCQgD+eJrqdQDMrMDMhpvZ\nLwDc/SN3f8jdjyM0Iu1CaLwL2RvVOV3vm0Oi36+kLCsm9IysaJ80ZcmkbMmR8Q+tuMLMcqLlG4Bq\ng4kMvU4IctJtbx9CZ4S3diD/iVH+fYBTCcF06uj+Nd3fdMcw2QD/a+fMzPYxs9uT146ZnWBm95hZ\nW3cvc/eX3f0WQoeEHL7eTiVTJxEap1/v7g96kCzr96LftakNeI2wz/8vzbqK++uE2rOelbRhSw4h\nsiPns1pm9htCwPocYfiHdP5DCLi/Ng2MmR1kZreaWXKfFxPuH+l6Eqb7P61KZf/XZYRjDZQPiPwE\n4ZidRPjS9Jc0762Jmu5zLqHHb0XV9ahsENeB1D0FWbIdC+MaXU/41nZrFUk3E6rY/5DaADP6ew9C\nNXdySonko54dbaj5AzM7NWVb7Qi9vpJj4CS9A+xl23e5/zbh8WZFmZRtAaGNzMlmdlyFdX8gjHfz\neIaPtDIxkdAzcWjUcxAAM8tjW8+thyt5b7Xc/Q1CW6p+hPGYJlcoe033N90xnEg4L8NTO0ZYGI7h\nbsIUPLtFi/clDNdwUYVtJfe9qAa7l2oLIYjqnLrQzPYndCAoI/SgrBEP43Y9AxxpZuVzaUa1qjeS\nEnS6ezFhCIYuVOh2H/UqOx14z93/VdNyZMrMjiAM+/ABcEYVNTYTCMfrzuh/K/n+doQhGq5iWxu2\n5PU3wszapqTtR/ogpDI5wLWWMjSDmZ1NGArmqTRt8h4htEUaQfiM2tGR7SeQ+T4/FKUdbqGnbzLt\nTwlBX6UawnUg9UOPC5unHKBvyg08OXdhd0KtwSZCr6MVlWXg7qssjBlzBfCmmT1FuKkeS2jk+Qd3\nTzaaT/aS+o2Z7Uao8s+kjBUtBSZZGJD0M0J37gLgVndPrcm6nzB8wD+i7tXfIHyI/YdQO5Lq42hb\n15tZd9KMlRU9Xvo14cb6pJk9GZXlEMI32LfYNsDhDnP3ZWZ2JeHG+KqFCXc3EMaR2gt4rJLHPTUx\nkdB7qowwflTq9mu6v187v+7+vpldBYwE3jKzGYRHTMcRgqon2dbt/n7CYIv/G92w/kPo9Xc64THL\niFru40xCLd3QKOBeSujBdXy0PIdtgV5NXUp4/Pi36Px8RLgeS/n6tXsV4dhdHQU8/wT2JkyivJbQ\ntigWUeA3mRAozAb6RwHNdmV095vd/XkzGwMMIpyzpwhflvoS2nGN82igTXdfZGYjCbVir5vZTMJg\nmX0IAXpN2gy2j/J4glDb3IfQRmlQmrSzCW249gQeTO0RWBs13OcXzezPhJr+181sFmG4jlMI13ZH\nwvmvTL1dB1J/VJPVPJUR/rFviH5+T2izsithbJofRN2K070v9dHQVYQaiLWE7s0DCN8uf+3u5d/W\nog+puwntDi5h26OaqiaITfdY78monD0IH3QbgQvc/fepCT2Mo3U5oX3DhcCRhBqv36bZ5uPRz97R\nvnw7Xdmib5g9CY1sf0wYr6gDoWbnoNRhCqpR1T6n7sNYQkDyCuED/9fA59H+ZuPD+FHCDWGFpxkd\nvib7W9n5dfc7gV8QHvmcTAikigmB+WnJGpUor58QBnn8DuHc/YIQJB1cYciCjLn7SkKj/nmEa+Di\nKP/RhEBvNaH9YKqqrsfU62EZIeB8jPDF5DzC46QT06RdHaW9g3BTvoTQ7uwh4MAKXxAyLkMGyyGc\nk12jvy8i1E7fyLb//RtIGW/J3QcTanyXE276vyaMvn+eu1+amrG7X034f9xA+N//PmFIhKcqKUs6\nCUJw+m/C/+rhhGvzYI9mDaiwzVK2DdVRk0eFlR6jmuwz4dz9Lir3hYT/kSGEcwnhC2rabWbxOqhu\nnTQgOWVlOlciItI4REM57OHue1WbOLvb3R0ojsbQqrjuYUKAtru7f16X5ZKGrdrHhVG334cJj2VK\nCN9YSgnPshPAm+5+SZR2AOHb6lbC/GY1+UYjIiJSKTM7mlCzWh+jnZ8D3G5m5/r2MxMUEtpkvaUA\nSyqqtibLzE4EznL3fhbmXLuI0Fh0pLvPN7NxbJsH6llCu548QuPZHllsDCwiIs1Q1P7zUMJgnV8A\n+9bgEX22ytCF0F4wj9DgfimhU8XJhE4fx6Z79C7NWyZtst4FWkZdt/MJtVTd3T05GvQswvhIvYAF\n7l4SDWr4HttGRhYREamtlYQONUuAE+s6wALwMM1RT0KHjQOBwYSOPs8AhyjAknQy6V24gdCj6R1C\nT5wT2H7cmvWEnmnt2H402w2kH6wPADNrTbhg/0vVPTJERKR5m0zKxNBmVlBP5UgQJlcfXnFFPZZJ\n6k4uYZzClzPt2ZpJkPVb4Bl3vy6qLn2e7cfDaUfoDr2O7Ue9TS6vTE/SzI0mIiIi0oAdRuXzUG4n\nkyDrC7YNNrgmes9rZna4hxnokyN8v0wYpK0VYSLMfYGqul7/F2DSpEl07ty5imT1q7S0lFWrVtGy\npYYUE5HGbbeTwpiZq6dPr+eSiNRcSUkJu+++O7m56eZ1j98nn3zC2WefDVH8kolMIofRwINm9gKh\nwfs1hPFgHoimB1gCTIkGMBxDiO5ygKHRKLeVKQXo3LkzXbt2zbS8da6kpIScnBwFWSLS6HWKOjrt\n1IC/2IpUpqSkhD322KMh3I8zbuJUbUndfSNwRppVR6RJOx4Yn+nGRUSk7nz60kv1XQSRZkUjvouI\niIjEQEGWiIiISAwUZImIiIjEQEGWiIiISAwUZImIiIjEQEGWiEgz0emgg+h00EH1XQyRZkNBloiI\niEgMFGSJiIjUodLSUk499VSuueaajNL/7ne/Y926dbXe3sMPP8yYMWO+tvyZZ57h+OOPZ+DAgQwc\nOJDzzz+fK664gnfeeac8zbXXXsvy5csBuPPOOzn77LN58MEHeeWVV+jXrx8XX3wxxcVVjTvevNX7\nsKkiIiLNyfz58yksLOTdd99l+fLl7LnnnlWmX7x4cWxlOeCAAxg+fNt814sXL+baa6/l3nvvpVOn\nTtx6663l62bOnMnjjz9Ox44due222zj++OM555xzYitbU6AgS0REpA7NmDGDI488ki5dujBlyhSu\nuOIKAJ5++mkmT55Mbm4u+fn5XHPNNTz44IMA/Pa3v2XEiBFcdtll3HzzzXTr1g2AM888s/z1xIkT\nWbhwIVu3bmXLli1cdNFFHHrooTUqW48ePTj00EOZMWMGAwYMKM//nnvuAeDqq6/miCOOYOHChbRu\n3ZqNGzdy4YUXMnHiRObPn09ZWRmdO3dm8ODB7Lrrrvz2t7+lXbt2rFixghNPPJGjjz6au+++m2XL\nllFSUkL37t256KKLaNGiBccccwxnnXUWr7zyCl988QV9+/bl1FNPBcI8x7Nnz6ZNmzYUFBRw6623\n0rZtW6ZMmcKjjz4KQPv27fn973/P3nvvnZXzlA0KskQka0pLSykqKoot75ycHFq0iKeVQ0FBQb1N\nPCvxqayhf2VTDNU0fU0VFRWxZMkSbrnlFrp168bgwYMZMGAAq1at4v777+f++++nY8eOTJ06lYkT\nJ3L11Vcze/ZsRo8eTbt27SrNd9WqVbz22mvcddddtGrVinnz5vHQQw/VOMgCKCwsZNGiRdstu+uu\nuzjyyCPLy7Fy5Ur22msvTj/9dObMmcOyZcsYN24cLVq0YObMmdx2222MGDECgF122YWHHnoIgNtu\nu41u3bpx9dVXk0gkGDFiBJMnT+aMM85g69attG/fnrFjx/Luu+8yaNAg+vTpw6JFi5gzZw5jx47l\nO9/5DnfccQeTJk2ie/fuTJ8+nccee4zWrVuzcOFCBg0axFNPPVXjfY6LgiwRyZqioiKGjJ5LXn6n\nrOe9+qMltGm3Wyx5b1r7KSMH96awsDDreTckmruw/j3xxBMcfPDBtG3bFjOjc+fOPPHEE7Rq1Yqe\nPXvSsWNHAE455ZTt3lcWTe5dmd13351rrrmGZ599lpUrV/L222+zZcuWWpUxJyeH1q1bp12Xrhz/\n+te/cHcuvPBCABKJxHbttH7wgx+U//3iiy/i7jz99NMAFBcXb/fF6ZBDDgGgW7dulJSUsGXLFl59\n9VUOP/xwdt55ZyDUpgHcfvvtLF++nH79+pWXa926daxbt45ddtmlVvuebQqyRJqZOGubVqxYQV5+\nJ9p26JL1vDetXRVb3tJ01TSwjDMQ3bJlC3PmzKF169acddZZlJWVsWnTJmbMmEG/fv22S1tcXMwn\nn3xS3l4rJyen/HdqoFNSUgLAe++9x+9//3tOO+00evbsyQEHHMDo0aNrVc533nmnRo/cEokE/fr1\n48QTTywv0/r168vXt2nTZru0N954Y/l+bdiwYbsgq2JwV1ZWRm5ubvn+A6xfv55169aRSCTo06cP\nV155Zfm6VatWNZgACxRkiTQ7cdc27db1u1nPV6QpePbZZ2nfvj0TJ04sX7ZhwwbOPPNM1q9fz2uv\nvcYXX3zBrrvuyhNPPMHrr7/OsGHDaNGiBVu3bgVCuyN3x8x4++23+eKLLwD497//jZlx6qmnkkgk\nuPPOO0kkEtWWqWLN1IsvvshLL71U3hYsEz179mTmzJkcddRR5OXlMX78eN5//31uv/32tGknT57M\nlVdeSXFxMddddx29evXi7LPPrrRsPXr04N577y1vnzV27FjKyso44ogjuP766/nVr37FN7/5TSZN\nmsTEiROZNWtWxmWPm4IskWYoztomEUnvySef5PTTT99uWdu2bTn55JN58cUXufDCC7nqqqvIyclh\nt91246qrrgLgsMMO47LLLmPYsGEMHDiQO++8k5kzZ9KtW7fyBvC9e/dm/vz5nHvuubRq1Yru3buz\nbt06Nm/eXGWZ3njjDQYOHFj+OtlzsH379gDb1SCl/p3qF7/4BatXr+Y3v/kNLVq0oFOnTlx77bVp\n33PppZdyzz330L9/f0pLS+nRo0d5LV7FtMnXBx10EB9++CGXX345O+20E926deOWW24hLy+PCy64\ngP79+9OiRQvatm3L3XffXeX+1rWc6p7zxsXMCoBlc+fOpWvXrvVShkyUlJSwcuVKWrZUPCpNw9Kl\nS7nhoTdiCbI+LXqVvPzdG13eG778mD+c94Mm3yZLpDErKSlhjz32qLf78UcffUTv3r0B9nL3okze\no8FIRURERGKgIEtEpJnQ3IUidUtBloiIiEgM1NBIRJq9skSCFStWxJK3BjkVab4UZIk0UHGNZxVX\nMNGYbV7/GXdN/Zy8/DVZzbe5DHIqIukpyBJpoOIaz0pjWaWngU5FJNsUZIk0YHHc+DWWlTQnyRHR\n46LhfaQqujpERJqJ5jZ3YXKcw7gmFU8kEvU6bpM0fLoyRERiEmeDelCj+ky0aNEitiCoprVkiUSC\nKVOmMG/ePBKJBFu3buXHP/4x5513HjvttFOtypBIJLj++utZsWIFffv2ZdasWdx5553lkykn/e1v\nf2PZsmXlkyvXhzPPPJObb765fJT6dDZu3Mj111/PqFGjYi/PtddeS7du3TjvvPNi24aCLBGRmMTV\noB7UqL4xGjVqFBs3bmTUqFHk5eXx1VdfMWzYMEaOHFk+DU1NffbZZyxevJhZs2aRk5ND3759s1zq\nurV+/Xrcvb6LkTUKskREYqQG9QLwySefMG/ePKZOnUqbNm0AaN26NVdccQVvvfUWEGpx7rrrLt5/\n/31ycnLo1asXAwYMoEWLFhxzzDGcddZZvPLKK3zxxRecfPLJ/PznP+fqq6+mpKSECy+8kJtuuolz\nzjmH6dOns/POOzNmzBgWL15Mhw4d6NChQ3nt1saNG7n77rtZtmwZJSUldO/enYsuuijtdvr27Vs+\nMfOkSZOYM2cOLVu2pEuXLlxzzTXk5eXx9NNPM2PGDMrKysjPz2fQoEHsueeeVR6PdPtzyimncNtt\nt7FlyxYGDhzIvffey/Lly7n77rtZv349paWlnHvuuZx22mksWrSI4cOH06ZNG7Zs2UJhYSHf//73\n6d+/PwCPPfYYL7/8MnfccQfDhw/njTfeYOPGjZSVlTFs2DB+9KMfxXWqt1NtkGVmvwbOBcqANsAB\nwGHAaCABvOnul0RpBwADga3AcHd/Kp5ii4iINB7vvvsuBQUF5QFWUocOHTj00EMBGDt2LPn5+Tz4\n4IOUlJQwdOhQHn/8cc4880y2bt1K+/btGTt2LO+++y6DBg3ixBNPZMSIEZx//vncd999wLZJladN\nm8bHH3/Mww8/zNatW7n88svZe++9Abjnnnvo1q0bV199NYlEghEjRjB58mTOOOOMtNvp06cPixYt\nYs6cOfzpT39i5513Zty4cUybNo399tuP2bNnM3bsWFq1asUrr7zCDTfcwIQJE6o8HpXtz1VXXVW+\nP6Wlpdx0001cd9117LPPPqxdu5Yrr7wSMwPg/fffZ+7cuXTu3JmXXnqJYcOGlQdZ06ZN44orruDf\n//43n3/+OY8//jgA9913H/fddx/jxo3LzomtRrVBlrs/DDwMYGZ3A+OBG4Ch7j7fzMaZWR/gRWAQ\n0B3IAxaY2Rx33xpb6UVERBqBFi1aUFZWVmWaRYsWMXbsWCD0WjzxxBOZOnUqZ555JgCHHHIIAN26\ndaOkpIQtW7ZUmtdrr71G7969yc3NJTc3l6OOOooPPvgAgBdffBF35+mnnwaguLh4u84B6bbz6quv\ncvjhh5emz9fFAAAZoklEQVTXhl188cUA3HvvvaxcuZJLL720fP82btzIhg0baNu2bZX7W93+fPTR\nR6xcuZLbbruNsrIyysrKKC4u5u2332bvvfemc+fOdO7cGYCDDjqI4uJi3nrrLb7xjW/w5ZdfcvDB\nBwNw+eWX89hjj7F8+XIWLVpUbbmyKePHhWZ2IPA9d7/UzG5y9/nRqlnA0YRarQXuXgKsM7P3gP2B\nxdkutIiI1Fxy3sLm1suwIdh333358MMP2bx583a1WZ999hmjRo3i5ptvJpFIbPeeRCKxXeP61q1b\nl/+dDDqqkro+tYNEaWkpN954Y/kjvQ0bNmwXZKVuJ5lPbm5ueS1Z8j0bNmwgkUhw9NFHM2DAgO32\nKZNAprr9SSQStGvXrryWrqSkhNatW9OhQwdef/118vLytkt/6qmnMm3aNFq1alX+iPP555/nj3/8\nI/379+eoo45i77335sknn6y2bNlSk36t1wI3pVm+HtgFaAesTVm+AcivdclERER2UDJQieOnYlBU\nlY4dO3LUUUdx2223sWnTJmBbG6z27dvTqlUrevXqxfTp04FQuzRz5kwOPPDAavNODU6Sf/fq1Ys5\nc+ZQXFxMcXExzz33XHmanj17Mnny5PLtXHfddUybNq3KvHv06MH8+fPZvHkzAA8//DBTpkyhZ8+e\nzJ07ly+++AKA6dOnM2TIkIyPS0W5ubnlx/Vb3/oWrVq14tlnnwXg008/5aSTTipvw1ZR3759mTdv\nHrNnz+bkk08G4J///CdHHnkk/fr1Y7/99mPu3Lk1Om87KqOaLDPLB7q5+wvRotQStgPWAOsIwVbF\n5SJNVlxT34Cmv5Gq1WZ4iF2jWpGlS5dWm7YpDA/RsmVL9thjj9i3kanBgwfzyCOPcOmll9KyZUuK\ni4s57LDDOPfccwG49NJLGTNmDP3796ekpIRevXpxzjnnAGxXi1Txdbq/TzjhBD7++GP69+9Pfn4+\nXbps63wxaNAg7rnnHvr3709paSk9evSgX79+VW7noIMO4sMPP+TSSy8FwvUxZMgQ2rRpw5lnnsmQ\nIUNo0aIFO++8M7fcckva/a+szKmvd9ttN/bZZx/OPfdcxo4dy7Bhwxg7dix//etfKS0t5bLLLuNH\nP/oRixYt+lr+HTt2ZL/99qO0tJRvfvObAPTr148hQ4bQp08fcnNzOfDAA5kzZ07a8sUhp7rqRgAz\nOwHo7e6Do9czgDvc/QUzGwfMA14A5gA9CQ3k/wX80N2LK8mzAFg2d+5cunbtmo19iUVyMDsNNifp\nLF26NJapb2Db9DfZ7pn2adGr5OXvHkuPN+VdN/km84acGl17j02+HoAzT0t/E0zS8BDSEJWUlNTr\n4K8fffQRvXv3BtjL3YsyeU+mJTXgg5TXQ4D7zWwnYAkwxd3LzGwMsADIITSMTxtgiTQlcXXR1/Q3\nUp2aXns5LULNlIaUEKkbGQVZ7j6ywuv3gCPSpBtP6H0oIiIi0qzpGZiISDNxwQX313cRRJqVeGbN\nFBEREWnmFGSJiIiIxEBBloiIiEgMFGSJiIiIxEBBloiIiEgMFGSJiDQTDzwwgAceGFB9QhHJCgVZ\nIiIiIjFQkCUiIiISAwVZIiIiIjFQkCUiIiISA02rIyIi2ylLJFixYkUseRcUFJCbmxtL3iINjYIs\nEZFmItO5Czev/4y7pn5OXv6arG5/09pPGTm4N4WFhVnNV6ShUpAlIiJfk5ffibYdutR3MUQaNbXJ\nEhEREYmBgiwRERGRGCjIEhEREYmBgiwRERGRGCjIEhFpJjR3oUjdUpAlIiIiEgMFWSIiIiIxUJAl\nIiIiEgMFWSIiIiIxUJAlIiIiEgNNqyMi0kxkOnehiGSHarJEREREYpBRTZaZXQOcCOwE/Al4AZgA\nJIA33f2SKN0AYCCwFRju7k/FUGYRERGRBq/amiwzOxz4sbsfAhwB7AmMAoa6++FACzPrY2a7A4OA\nHwPHArea2U6xlVxERESkAcvkceExwJtmNh14ApgJdHf3+dH6WcDPgF7AAncvcfd1wHvA/jGUWURE\nRKTBy+RxYUdC7dXxwN6EQCs1OFsP7AK0A9amLN8A5GenmCIiIiKNSyZB1mpgibuXAO+a2Raga8r6\ndsAaYB0h2Kq4XEREGoDkvIXqZShSNzJ5XLiA0MYKM9sD2BmYG7XVAjgOmA+8DBxqZq3MLB/YF3gz\n+0UWERERafiqrcly96fM7DAzWwTkABcDRcADUcP2JcAUdy8zszGEoCyH0DC+OL6ii4iIiDRcGQ3h\n4O7XpFl8RJp044HxO1gmERERkUZPg5GKiIiIxEBBloiIiEgMNHehiEgzoV6FInVLNVkiIiIiMVCQ\nJSIiIhIDBVkiIiIiMVCQJSIiIhIDBVkiIiIiMVCQJSLSTDzwwIDy+QtFJH4KskRERERioCBLRERE\nJAYajFREROpEWSLBihUrYsu/oKCA3Nzc2PIXqSkFWSIiUic2r/+Mu6Z+Tl7+mqznvWntp4wc3JvC\nwsKs5y1SWwqyRESkzuTld6Jthy71XQyROqEgS0SkmdDchSJ1S0GWNHmlpaUUFRXFknec7UtERKRx\nU5AlTV5RURFDRs8lL79T1vNe/dESduv63aznKyIijZ+CLGkW4moHsmntqqznKSIiTYPGyRIRERGJ\ngYIsERERkRgoyBIRaSY0d6FI3VKQJSIiIhIDBVkiIiIiMVCQJSIiIhIDBVkiIiIiMVCQJSIiIhID\nDUYqItJMaO5CkbqVUZBlZouBtdHLZcAfgQlAAnjT3S+J0g0ABgJbgeHu/lS2CywiIiLSGFQbZJlZ\nawB3PzJl2QxgqLvPN7NxZtYHeBEYBHQH8oAFZjbH3bfGU3QRERGRhiuTmqwDgJ3NbDaQC1wHdHf3\n+dH6WcDRhFqtBe5eAqwzs/eA/YHF2S+2iIiISMOWScP3TcDt7n4McDEwCchJWb8e2AVox7ZHigAb\ngPwslVNERESkUckkyHqXEFjh7u8Bq4HdU9a3A9YA6wjBVsXlIiIiIs1OJkFWf+AOADPbgxBIzTGz\nw6P1xwHzgZeBQ82slZnlA/sCb2a/yCIiUhuau1CkbmXSJms88JCZzSe0uzqXUJv1gJntBCwBprh7\nmZmNARYQHicOdffieIotIiIi0rBVG2RFvQPPSbPqiDRpxxOCMhEREZFmTSO+i4iIiMRAQZaIiIhI\nDBRkiYiIiMRAcxeKiDQTmrtQpG6pJktEREQkBgqyRERERGKgIEtEREQkBgqyRERERGKgIEtEREQk\nBgqyRESaCc1dKFK3NISDiIg0emWJBCtWrIgt/4KCAnJzc2PLX5omBVkiItLobV7/GXdN/Zy8/DVZ\nz3vT2k8ZObg3hYWFWc9bmjYFWSIi0iTk5XeibYcu9V0MkXJqkyUiIiISAwVZIiIiIjHQ40JpEEpL\nSykqKool7zgbw4o0Jpq7UKRuKciSBqGoqIgho+eSl98p63mv/mgJu3X9btbzFRERqYqCLGkw4mq0\numntqqznKSIiUh21yRIRERGJgYIsERERkRgoyBIRERGJgYIsEZFmQnMXitQtBVkiIiIiMVCQJSIi\nIhIDBVkiIiIiMVCQJSIiIhIDBVkiIiIiMchoxHcz6wS8AhwFlAITgATwprtfEqUZAAwEtgLD3f2p\nOAosIiK1o7kLRepWtTVZZtYS+DOwKVo0Chjq7ocDLcysj5ntDgwCfgwcC9xqZjvFVGYRERGRBi+T\nx4UjgXHASiAH6O7u86N1s4CfAb2ABe5e4u7rgPeA/WMor4iIiEijUGWQZWbnAp+6+7OEAKvie9YD\nuwDtgLUpyzcA+dkrpoiIiEjjUl2brPOAhJn9DDgAeAT4Zsr6dsAaYB0h2Kq4XERERKRZqjLIitpd\nAWBm84CLgNvN7Cfu/gJwHDAPeBkYbmatgDbAvsCbsZVaREREpIHLqHdhBUOA+6OG7UuAKe5eZmZj\ngAWEx4pD3b04i+UUEZEdlJy3UL0MRepGxkGWux+Z8vKINOvHA+OzUCYRERGRRk+DkYqIiIjEQEGW\niIiISAwUZImIiIjEQEGWiIiISAxq07tQREQaIfUqFKlbqskSERERiYGCLBEREZEYKMgSERERiYGC\nLBEREZEYKMgSERERiYGCLBGRZuKBBwaUz18oIvFTkCUiIiISAwVZIiIiIjFQkCUiIiISAwVZIiIi\nIjFQkCUiIiISA81dKCLSTGjuQpG6pZosERERkRgoyBIRERGJgYIsERERkRgoyBIRERGJgYIsERER\nkRgoyBIRaSY0d6FI3dIQDlIjpaWlFBUVZT3fFStWZD1PERGR+qQgS2qkqKiIIaPnkpffKav5rv5o\nCbt1/W5W8xQREalPCrKkxvLyO9G2Q5es5rlp7aqs5iciIlLf1CZLREREJAbV1mSZWQvgfsCABHAR\n8BUwIXr9prtfEqUdAAwEtgLD3f2peIotIiIi0rBl8rjwBKDM3Q81s8OBPwI5wFB3n29m48ysD/Ai\nMAjoDuQBC8xsjrtvjavwIiKSOc1dWDtliUSsnXMKCgrIzc2NLX+pP9UGWe4+w8yejF5+G/gSOMrd\n50fLZgFHE2q1Frh7CbDOzN4D9gcWZ7/YIiIidWPz+s+4a+rn5OWvyXrem9Z+ysjBvSksLMx63lL/\nMmr47u4JM5sAnAScBvwsZfV6YBegHbA2ZfkGID87xRQREak/cXT4kaYv44bv7n4u0A14AGiTsqod\nsAZYRwi2Ki4XERERaXaqDbLM7BwzuyZ6uQUoBV6J2mcBHAfMB14GDjWzVmaWD+wLvBlDmUVEREQa\nvEweF/4deMjM/hGlvwx4B3jAzHYClgBT3L3MzMYAC9jWML44pnKLiIiINGiZNHzfBJyRZtURadKO\nB8bveLFERCTbkvMWqpehSN3QYKQiIiIiMVCQJSIiIhIDBVkiIiIiMVCQJSIiIhIDBVkiIiIiMcho\nxHcREWn81KtQpG6pJktEREQkBgqyRERERGKgIEtEREQkBgqyRERERGKgIEtEREQkBgqyRESaiQce\nGFA+f6GIxE9BloiIiEgMFGSJiIiIxEBBloiIiEgMFGSJiIiIxEBBloiIiEgMNHehiEgzobkLReqW\narJEREREYqAgS0RERCQGCrJEREREYqAgS0RERCQGCrJEREREYqAgS0SkmdDchSJ1S0M4iIiI1JOy\nRIIVK1bEkndBQQG5ubmx5C2ZUZAlIiJSTzav/4y7pn5OXv6arOa7ae2njBzcm8LCwqzmKzWjIEtE\nRKQe5eV3om2HLvVdDIlBlUGWmbUEHgQKgFbAcOBtYAKQAN5090uitAOAgcBWYLi7PxVbqUVEREQa\nuOoavp8DfO7uPwGOBe4GRgFD3f1woIWZ9TGz3YFBwI+jdLea2U4xlltERESkQavuceHfgMnR37lA\nCdDd3edHy2YBRxNqtRa4ewmwzszeA/YHFme/yCIiUhuau1CkblUZZLn7JgAza0cItq4DRqYkWQ/s\nArQD1qYs3wDkZ7WkIiIiIo1IteNkmdm3gHnAw+7+V0KtVVI7YA2wjhBsVVwuIiIi0ixV1/B9d2A2\ncIm7Pxctfs3MfuLuLwDHEQKwl4HhZtYKaAPsC7wZX7GlKqWlpRQVFcWSd1zjuYiIiDQ11bXJuhZo\nD1xvZjcAZcDlwNioYfsSYIq7l5nZGGABkENoGF8cY7mlCkVFRQwZPZe8/E5Zz3v1R0vYret3s56v\niIhIU1Ndm6zBwOA0q45Ik3Y8MD47xZIdFde4K5vWrsp6niIiIk2R5i4UEWkmNHehSN1SkCUiIiIS\nAwVZIiIiIjFQkCUiIiISAwVZIiIiIjFQkCUiIiISg+rGyRIRkSZCcxeK1C3VZImIiIjEQEGWiIiI\nSAwUZImIiIjEQEGWiIiISAwUZImIiIjEQEGWiEgzobkLReqWgiwRERGRGCjIEhEREYmBgiwRERGR\nGCjIEhEREYmBgiwRERGRGGjuQhGRZkJzF4rULdVkiYiIiMRAQZaIiIhIDPS4UEREpIkpSyRYsWJF\nbPkXFBSQm5sbW/5NhYIsERGRJmbz+s+4a+rn5OWvyXrem9Z+ysjBvSksLMx63k2NgiwREZEmKC+/\nE207dKnvYjRrapMlItJMaO5CkbqlIEtEREQkBgqyRERERGKQUZssMzsIGOHuPzWzQmACkADedPdL\nojQDgIHAVmC4uz8VT5FFREREGr5qgywz+x3wS2BDtGgUMNTd55vZODPrA7wIDAK6A3nAAjOb4+5b\nYyp3o1daWkpRUVEsecfZbVdEREQyk0lN1vtAX+Av0ese7j4/+nsWcDShVmuBu5cA68zsPWB/YHGW\ny9tkFBUVMWT0XPLyO2U979UfLWG3rt/Ner4iIiKSuWqDLHefZmbfTlmUk/L3emAXoB2wNmX5BiA/\nKyVswuLqXrtp7aqs5ykijZ/mLhSpW7Vp+J5I+bsdsAZYRwi2Ki4XERERaZZqE2S9amY/if4+DpgP\nvAwcamatzCwf2Bd4M0tlFBEREWl0ajPi+xDgfjPbCVgCTHH3MjMbAywgPE4c6u7FWSyniIiISKOS\nUZDl7h8Ch0R/vwcckSbNeGB8NgsnIiIi0lhpMFIRERGRGCjIEhFpJjR3oUjdUpAlIiIiEgMFWSIi\nIiIxUJAlIiIiEoPaDOEgIiIizVRZIhHrHLkFBQXk5ubGln9dUpAlIiIiGdu8/jPumvo5efnZn9hl\n09pPGTm4N4WFhVnPuz4oyKpGaWkpH3zwQdaj6ji/BYiIpKO5CyVb4pp7t6lRkFWNpUuXcvXY58nL\n75TVfFd/tITdun43q3mKiIhIw6EgKwNxROyb1q7Kan4iIiLSsKh3oYiIiEgMFGSJiIiIxEBBloiI\niEgMFGSJiDQTmrtQpG4pyBIRERGJgXoXioiISINQ1WjypaWlbNiwgZYtax+6FBYW1ulo8gqyRERE\npEGofjT5t2ud96a1n/KXW8+iW7dutc6jphRkiYiISIPRlEaTV5ssERERkRg0iZqs0tJSli5dGkve\nRUVFseQrIlLXNHehSN1qEkHW0qVL+eW1j2Z9fkHQHIMiIiJSO00iyIL4nuFqjkERERGpDbXJEhER\nEYmBgiwRERGRGDSZx4UiIiIilSlLJFi2bFmt379qVc2bD9V7kLVs2TI2bdq0w3mIiEjVkvMWqpeh\nNEeb13/GDfd9Tl5+7UYjKNmyrsbvqfcga+i4hbT8xi47lId6AIqIiEh1dqST3NZNbWr8nqwGWWaW\nA/wJOADYAlzg7h9U9Z6d83dnp7xdd2i76gEoIiIiDU22G76fBLR290OAa4FRWc5fREREpFHIdpB1\nKPAMgLu/BByY5fxFREREGoVst8naBVib8rrEzFq4eyJN2lyAdZ++T27rdju00U1ffkxp8cZaNUqr\nr7wbY5njzLsxlrmx5t0Yy9xY825oZV6ZUwbAmv8uyXremWhox6Mp590Yy9zQ8y79an3yz9xM35Pt\nIGsdkBoxVRZgAfwPwOo3p2Rlw8XA+mpTNay8G2OZ48y7MZa5sebdGMvcWPNuSGXut3vU/nXxQ1nP\nO1MN6Xg09bwbY5kbSd7/A2TURTHbQdZC4HhgipkdDLxRRdqXgcOA/wKlWS6HiIiISDblEgKslzN9\nQ05ZWVnWtp7Su3D/aNF57v5u1jYgIiIi0khkNcgSERERkUBzF4qIiIjEQEGWiIiISAwUZImIiIjE\noN7mLjSzvsCp7n529PokYCSwPEpyo7vPr6/yNUVpjvlBwF3AVuBZd/9DfZavKTOzj4BkJ5B/uft1\n9Vmepqo2U3vJjjOzxWwbI3GZu59fn+VpyqLP7RHu/lMzKwQmAAngTXe/pF4L10RVOOY/BGay7fN8\nnLtPruy99RJkmdlo4Gjg9ZTFPYDfufu0+ihTU1fJMf8z0Nfdi8zsKTM7wN3/XT8lbLqiD8LF7t6n\nvsvSDJRP7RV9MI6KlklMzKw1gLsfWd9laerM7HfAL4EN0aJRwFB3n29m48ysj7vPqL8SNj1pjnkP\n4A53vzOT99fX48KFwMUVlvUA+pvZC2Y20sz0KDO7tjvmZtYOaOXuRdGi2cBR9VCu5qAH0NXM5pnZ\nTDPrVt8FasI0tVfdOwDY2cxmm9n/RcGtxON9oG/K6x4pT3xmoc/wOHztmAO/MLN/mNkDZrZzVW+O\nNZAxs/5m9oaZ/Sfld49KqtbmAIPc/SdAW+CiOMvWVNXgmO9CGKE/aT2QX3clbZrSHX/CgLt/jL7p\n3wpMrN9SNmlpp/aqr8I0E5uA2939GMIXuUk65vGInvSUpCzKSflbn+ExSHPMXyI8dTsc+AC4qar3\nx/q40N0fBB7MMPlD7p78cJwBnBxPqZq2GhzzdYQbUlI7YE0shWpG0h1/M2tD9E/q7gvN7H/qo2zN\nRE2m9pLseJfwbR93f8/MVhNGxf64XkvVPKRe2/oMrxvTU2KVacCYqhI3pG8b/zGzPaK/ewOL67Mw\nTZ27rwe+MrO9osbCxwDqaBCPG4HBAGZ2ALCifovTpC0Efg6QwdRekh39gTsAos/wdoTaW4nfq2b2\nk+jv49BneF2YbWbJZgjVxir11rswjfOBaWa2CXgbuL+ey9McXAQ8Sgi257h7xvMxSY2MACaa2S8I\nPTnPrd/iNGnTgJ+Z2cLo9Xn1WZhmYjzwkJnNJ9Ss9FftYZ0ZAtxvZjsBS4Ap9Vye5uBiYKyZFQOf\nAAOrSqxpdURERERi0JAeF4qIiIg0GQqyRERERGKgIEtEREQkBgqyRERERGKgIEtEREQkBgqyRERE\nRGKgIEtEREQkBgqyRERERGLw/wFhdZiyaFWrOAAAAABJRU5ErkJggg==\n" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Let's look at the distribution of differences\n", "f, ax = plt.subplots(figsize=(10, 5))\n", "_ = ax.hist(diff, bins=30)\n", "_ = ax.axvline(actual_diff, color='r', ls='--')\n", "axfill = ax.fill_between([clo, chi], *ax.get_ylim(), alpha=.1, color='k')\n", "ax.set_title('Distribution of votes randomized by region', fontsize=20)\n", "ax.legend([ax.lines[0], axfill], ['Actual Difference', 'Confidence Interval'], fontsize=12)\n", "_ = plt.setp(ax, xlim=[-15, 15])" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "It looks like now our confidence intervals are even wider than before. This is becase basically any change to our voting system that deviates away from a completely random 50/50 split will increase the variability in the outcome." ] } ], "metadata": { "category": "til", "date": "2016-07-08", "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "2.7.11" }, "redirect": "brexit-voting-randomness", "tags": "python, voting, statistics, computation, fft", "title": "Could Brexit have happened by chance?" }, "nbformat": 4, "nbformat_minor": 0 }