orbeckst/pandemic_season_2.ipynb

## pandemic_season_2.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pandemic Season 2 (contains SPOILERS) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## December\n",
    "\n",
    "We want to know how many supply cubes we should put on _Utopia_ at the beginning of the December game. Because there are 8 Utopia cards in the infection deck, there might be a high chance to start with a plague in Utopia right away during step 2 of game setup when the grid is infected. Because we do not have many supply cubes in the reserve, we do not want to over-provide Utopia (also, they don't really deserve it...), just enough so that we can handle the situation.\n",
    "\n",
    "We want to assess the chances to draw Utopia in setup so that we can make a good choice for how to provide.\n",
    "\n",
    "The parameters for the calculation are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Ninfect = 40  # count all cards in the infection deck (including Utopia)\n",
    "NUtopia = 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(The number of cards in the infection deck will depend on how you played the game. We had 40 in December.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Python packages\n",
    "We use Python and numpy for calculations and matplotlib for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.special\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.matplotlib.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probability to draw exactly $k$ Utopia cards during _Infect the Grid_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw $n$ cards from $N$. In how many ways will this draw contain exactly $k$ copies out of the $K$ total Utopia cards: $m = \\binom{K}{k} \\binom{N-K}{n-k}$. Thus the probability for the draw is $m$ divided by the total number of ways to draw $n$ from $N$:\n",
    "\\begin{gather}\n",
    "p_k = \\frac{\\binom{K}{k} \\binom{N-K}{n-k}}{\\binom{N}{n}}\n",
    "\\end{gather}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def draw_copies(k, n=6, N=Ninfect, K=NUtopia):\n",
    "    \"\"\"Calculate the number of ways to draw k copies of a card\n",
    "    in n draws, out N total cards with K copies of the card in the deck.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    k : int\n",
    "        In how many ways will you get exactly k copies\n",
    "        of the card of interest (of which there are K\n",
    "        copies in the deck). (k <= K)\n",
    "    n : int\n",
    "        Number of cards drawn; n >= k and n <= N.\n",
    "    N : int\n",
    "        Total number of cards in the deck.\n",
    "    K : int\n",
    "        Total number of identical cards of interest.\n",
    "        k <= K <= N\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        Total number of ways to draw the k cards.\n",
    "    \"\"\"\n",
    "    # assert all(k <= n)\n",
    "    # assert all(k <= K)\n",
    "    # assert n <= N\n",
    "    g = n - k\n",
    "    G = N - K\n",
    "    c_g = scipy.special.comb(G, g)\n",
    "    c_k = scipy.special.comb(K, k)  # draw k from K\n",
    "    return c_g * c_k\n",
    "\n",
    "def p_draw_copies(k, n=6, N=Ninfect, K=NUtopia):\n",
    "    \"\"\"Probability to draw exactly k copies of a card\n",
    "    in n draws, out N total cards with K copies of the card in the deck.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    k : int or array of type int\n",
    "        In how many ways will you get exactly k copies\n",
    "        of the card of interest (of which there are K\n",
    "        copies in the deck). (k <= K)\n",
    "    n : int\n",
    "        Number of cards drawn; n >= k and n <= N.\n",
    "    N : int\n",
    "        Total number of cards in the deck.\n",
    "    K : int\n",
    "        Total number of identical cards of interest.\n",
    "        k <= K <= N\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        probability to draw the k cards.\n",
    "    \"\"\"\n",
    "    N_draws = scipy.special.comb(N, n)\n",
    "    return draw_copies(k) / N_draws"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, the probability for drawing exactly 0  Utopia cards is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.23608709924499399"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_draw_copies(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "so about 24%.\n",
    "\n",
    "What is the probability to get exactly 4 of those cards:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.009045482729693256"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_draw_copies(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is only 0.9%. Given that the probability for any single Utopia card is about 1 in 5, you might have guessed that the probability of drawing 4 would have been higher and you would have put 4 cubes on the city – almost certainly a waste that you won't be able to afford in the December game."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis of the probabilities to draw exactly $k$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now calculate the probabilities for all values of $k$, the exact number of Utopia cards to be drawn:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 6\n",
    "k_values = np.arange(n+1)\n",
    "p = p_draw_copies(k_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The array `p` contains the probabilities $p_k$ for drawing exactly $k$ Utopia cards during setup:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  2.36087099e-01,   4.19710399e-01,   2.62318999e-01,\n",
       "         7.23638618e-02,   9.04548273e-03,   4.66863625e-04,\n",
       "         7.29474414e-06])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write out the probabilities as a table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.2361\n",
      "1 0.4197\n",
      "2 0.2623\n",
      "3 0.0724\n",
      "4 0.0090\n",
      "5 0.0005\n",
      "6 0.0000\n"
     ]
    }
   ],
   "source": [
    "print(\"\\n\".join([\"{} {:.4f}\".format(kk, pp) for \n",
    "      kk,pp in zip(k_values, p)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot of the probability $p_k$ that exactly $k$ supply cubes are removed from Utopia in step 2 of game set-up:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAENCAYAAAAYIIIKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcFPX+P/DXXkBURGFXWBGsh+s9DxFthXg5IvstT3aM\nLNN6qEcwtRCvleYltQy/dLxVClmEZEfr4CUsq6MdNLMkCzW8YQZpJoLg7qKAgrDM/P7w53xdWXFH\nl2XXXs/Hw8djZ/YzM69Zd/fNzOx8PgpRFEUQERE5SNncAYiIyLOwcBARkSwsHEREJAsLBxERycLC\nQUREsrBwEBGRLCwcREQkCwsHERHJwsJBRESysHAQEZEs6uYO0FSKi4ubOwK0Wi1MJlNzx3CYp+UF\nmNlVPC2zp+UF3CNzcHCwQ+14xEFERLKwcBARkSwsHEREJAsLBxERycLCQUREsrBwEBGRLCwcREQk\nCwsHERHJwsJBRESy3LF3jv+Z1I8f6pT1lDplLYAq7XMnrYmI3BGPOIiISBYWDiIikoWFg4iIZGHh\nICIiWVg4iIhIFpcVjry8PEydOhWTJ0/Gli1bbthu7969ePrpp/Hbb79J87KysjB58mRMnToVeXl5\nrohLREQ34JLCIQgC0tPTMWfOHKxYsQJ79uxBUVFRg3bV1dX4z3/+g65du0rzioqKkJOTg+XLl2Pu\n3LlIT0+HIAiuiE1ERHa4pHAUFhZCp9MhKCgIarUaUVFRyM3NbdAuMzMTQ4cOhZeXlzQvNzcXUVFR\n8PLyQmBgIHQ6HQoLC10Rm4iI7HDJDYAWiwUajUaa1mg0KCgosGlz8uRJmEwm3H///di6davNstce\ngQQEBMBisTTYRnZ2NrKzswEAycnJ0Gq1zt4N2dRqtUtyOOvGPWdx5WvvqtfYmZi56XlaXsCzMruk\ncIii2GCeQqGQHguCgLVr1yIhIcGhZe0xGo0wGo3SdHOP3Qu4xxjCzcGV++yJrzEzNz1Pywu4R2ZH\nxxx3SeHQaDQwm83StNlshr+/vzRdU1OD06dP47XXXgMAnD9/Hv/85z8xc+bMBstaLBYEBAS4IjYR\nEdnhksKh1+tRUlKCsrIyBAQEICcnB1OmTJGeb9WqFdLT06XphQsXYvTo0dDr9fD29sY777yDxx57\nDOXl5SgpKUGXLl1cEZuIiOxwSeFQqVSIj49HUlISBEFAdHQ0QkNDkZmZCb1eD4PBcMNlQ0ND0adP\nH8yYMQNKpRLjxo2DUsnbT4iImotCdPQigocpLi5u7gguO2fprN5xncWVveO6w3lhuZi56XlaXsA9\nMjt6jYN/uhMRkSwsHEREJAsLBxERycLCQUREsrBwEBGRLCwcREQkCwsHERHJwsJBRESysHAQEZEs\nLBxERCQLCwcREcnCwkFERLKwcBARkSwsHEREJAsLBxERyeKSgZwAIC8vDxkZGRAEATExMYiNjbV5\n/uuvv8b27duhVCrh4+ODiRMnIiQkBGVlZZg+fbrUT3zXrl0xYcIEV8UmIqLruKRwCIKA9PR0zJs3\nDxqNBrNnz4bBYEBISIjUpl+/fnj44YcBAPv27cPatWsxd+5cAIBOp8OSJUtcEZWIiG7CJaeqCgsL\nodPpEBQUBLVajaioKOTm5tq0adWqlfS4pqYGCoXCFdGIiEgmlxxxWCwWaDQaaVqj0aCgoKBBu23b\ntuHLL7+E1WrF/PnzpfllZWWYOXMmWrZsiZEjR6Jnz56uiE1ERHa4pHDYG9bc3hHF4MGDMXjwYHz/\n/ffYvHkzEhMT4e/vj9TUVLRp0wYnTpzAkiVLsGzZMpsjFADIzs5GdnY2ACA5ORlarbZpdkYGtVrt\nkhylTb4FeVz52rvqNXYmZm56npYX8KzMLikcGo0GZrNZmjabzfD3979h+6ioKKSlpQEAvLy84OXl\nBQDo3LkzgoKCUFJSAr1eb7OM0WiE0WiUppt70HfAPQafbw6u3GdPfI2Zuel5Wl7APTJf/RHSzbjk\nGoder0dJSQnKyspgtVqRk5MDg8Fg06akpER6fODAAXTo0AEAUFFRAUEQAAClpaUoKSlBUFCQK2IT\nEZEdLjniUKlUiI+PR1JSEgRBQHR0NEJDQ5GZmQm9Xg+DwYBt27bh8OHDUKlU8PX1xaRJkwAA+fn5\n2LBhA1QqFZRKJcaPHw9fX19XxCYiIjtcdh9HREQEIiIibOaNGDFCehwXF2d3ucjISERGRjZpNiIi\nchzvHCciIllYOIiISBYWDiIikoWFg4iIZGHhICIiWVg4iIhIFhYOIiKShYWDiIhkYeEgIiJZWDiI\niEgWFg4iIpKFhYOIiGRh4SAiIllYOIiISBYWDiIikoWFg4iIZHHZQE55eXnIyMiAIAiIiYlBbGys\nzfNff/01tm/fDqVSCR8fH0ycOBEhISEAgKysLOzcuRNKpRJxcXEIDw93VWwiIrqOSwqHIAhIT0/H\nvHnzoNFoMHv2bBgMBqkwAEC/fv3w8MMPAwD27duHtWvXYu7cuSgqKkJOTg6WL1+O8vJyLFq0CG+/\n/TaUSh4sERE1B5d8+xYWFkKn0yEoKAhqtRpRUVHIzc21adOqVSvpcU1NDRQKBQAgNzcXUVFR8PLy\nQmBgIHQ6HQoLC10Rm4iI7HDJEYfFYoFGo5GmNRoNCgoKGrTbtm0bvvzyS1itVsyfP19atmvXrlKb\ngIAAWCyWpg9NRER2uaRwiKLYYN7VI4prDR48GIMHD8b333+PzZs3IzEx0e6y9mRnZyM7OxsAkJyc\nDK1We3uhnUCtVrskR2mTb0EeV772rnqNnYmZm56n5QU8K7NLCodGo4HZbJamzWYz/P39b9g+KioK\naWlpdpe1WCwICAhosIzRaITRaJSmTSaTM6LfFq1W6xY5XM2V++yJrzEzNz1Pywu4R+bg4GCH2rnk\nGoder0dJSQnKyspgtVqRk5MDg8Fg06akpER6fODAAXTo0AEAYDAYkJOTg7q6OpSVlaGkpARdunRx\nRWwiIrLDJUccKpUK8fHxSEpKgiAIiI6ORmhoKDIzM6HX62EwGLBt2zYcPnwYKpUKvr6+mDRpEgAg\nNDQUffr0wYwZM6BUKjFu3Dj+ooqIqBkpREcvIniY4uLi5o7gskPP+vFDm3wbcqjSPnfZttzh8F4u\nZm56npYXcI/MbnWqioiI7hwsHEREJAsLBxERycLCQUREsrBwEBGRLA4XjsrKyqbMQUREHsLh+zhe\neOEFhIWFYcCAATAYDFCrXdYjOxERuRGHjzhSU1PRu3dvfPbZZxg/fjzee+89/PLLL02ZjYiI3JDD\nhw1+fn549NFH8eijj6K4uBi7d+/GypUroVAo0L9/fwwaNAjt27dvyqxEROQGbuni+Pnz53H+/HlU\nV1cjKCgIFosFM2fOxJYtW5ydj4iI3IzDRxynT5/Gd999h++++w4+Pj7461//iqVLl0o91T755JN4\n+eWXGwwJS0REdxaHC8eCBQvQt29fvPjii3Z7pw0MDMSjjz7q1HBEROR+HC4cL730Enr16tVgfmFh\noVRIRowY4bxkRETklhwuHG+++SbWrl3bYH5SUhIyMjKcGorufM7q0ddZox+6skdfIk9308IhCAKA\nK8O/Xv13VWlpKVQqVdOlIyIit3PTwvHMM89Ij0eOHGnznFKpxBNPPOH8VERE5LZuWjhWrVoFURSx\ncOFCvPbaa9J8hUIBPz8/eHt7O7ShvLw8ZGRkQBAExMTENPj11RdffIEdO3ZApVLBz88PL7zwgnRf\nyIgRI9CpUycAVwY7mTVrlsM7SEREznXTwnH1yzs1NfWWNyIIAtLT0zFv3jxoNBrMnj0bBoMBISEh\nUpu7774bycnJaNGiBb7++musW7cO06dPBwB4e3tjyZIlt7x9IiJynkYLx3vvvYeJEycCuHLkcSOJ\niYmNbqSwsBA6nQ5BQUEAgKioKOTm5toUjt69e0uPu3btiu++++7m6YmIyOUaLRyBgYHS46tf+rfC\nYrFAo9FI0xqNBgUFBTdsv3PnToSHh0vTdXV1eOWVV6BSqfD444/jwQcfbLBMdnY2srOzAQDJycnQ\narW3nNdZ1Gq1S3I465dFzuLIPntiZmdx1fvCmTwts6flBTwrc6OF49oL38OHD7/ljVz7S6yrFAqF\n3ba7d+/GiRMnsHDhQmleamoqAgICUFpaitdffx2dOnWCTqezWc5oNMJoNErTzT3oO+Aeg883B0/c\nZ1dm9sT3hadl9rS8gHtkDg4Odqhdo4XjyJEjDq3k2tNM9mg0GpjNZmnabDbD39+/QbtDhw4hKysL\nCxcuhJeXlzT/arcmQUFB6NWrF37//fcGhYOIiFyj0cLx7rvv3nQFCoWi0esfAKDX61FSUoKysjIE\nBAQgJycHU6ZMsWlz8uRJpKWlYc6cOWjbtq00v6qqCi1atICXlxcqKipw/PhxPP744zfNRURETaPR\nwpGSkuKUjahUKsTHxyMpKQmCICA6OhqhoaHIzMyEXq+HwWDAunXrUFNTg+XLlwP4v5/dnjlzBu+/\n/z6USiUEQUBsbKzNRXUiInItlw3jFxERgYiICJt51/Zt9eqrr9pdrnv37li2bFmTZiMiIsc1Wjim\nT5+OFStWALgydOyNOHJKi4iI7gyNFo6r93AAwOTJk5s8jLtwpw742PkeEbmbRgtHjx49pMf2ulQn\nIqI/H4evcVitVmzevBl79uxBeXk5/P39ERUVhWHDhjncXxUREXk+hwtHWloaiouLERcXh/bt2+Pc\nuXPYsmULPvjgAyQkJDRlRiIiciMOF47c3FysXLkSrVu3BgCEhISga9euf6prH0REBCgdbdiuXTtc\nvnzZZl5tba3dO8CJiOjO5XCXIwMGDMDixYsxePBgqQuR7du3Y8CAAU0ekoiI3IfsLkeysrJsprOz\nsxsMykRERHcul3Q5QkREdw6Hr3EQEREBMn5VdenSJWzcuBH5+fmorKy0GWODXY4QEf15OHzE8cEH\nH+DkyZN46qmnUFVVhfj4eGi1WgwZMqQp8xERkZtxuHAcOnQIL774Ih544AEolUo88MADmD59OscG\nJyL6k3G4cIiiiFatWgEAfHx8cPHiRbRr1w5nz55tsnBEROR+HL7GcddddyE/Px9/+ctf0KNHD6Sn\np8PHxwcdOnRwaPm8vDxkZGRAEATExMQ0+AnvF198gR07dkClUsHPzw8vvPAC2rdvDwDYtWsXPv30\nUwDAsGHDMHDgQEdjExGRkzl8xDFx4kTpizw+Ph5eXl64ePEiEhMTb7qsIAhIT0/HnDlzsGLFCuzZ\nswdFRUU2be6++24kJydj6dKliIyMxLp16wBcGTp206ZNWLx4MRYvXoxNmzahqqpKzj4SEZETOXzE\nERQUJD2+ekTgqMLCQuh0OmkdUVFRyM3NtRkCtnfv3tLjrl27StdO8vLyEBYWBl9fXwBAWFgY8vLy\n0K9fP4e3T0REziNr6NidO3fadKvet29fREdHQ6FQNLqcxWKBRqORpjUaDQoKChrdTnh4uN1lAwIC\nYLFY5MQmIiIncrhwrFu3Drm5uRgyZAi0Wi1MJhO2bt2K4uJijBo1qtFlr73n46obFZvdu3fjxIkT\nWLhw4Q3XZ2/Z7OxsZGdnAwCSk5Oh1WobzdQYZ4zc5yyO7Ic75QXu3MzOolarXbo9Z/C0zJ6WF/Cs\nzA4Xjl27duHNN9+0+es/IiICs2bNumnhuNop4lVms9lur7qHDh1CVlYWFi5cCC8vLwBXjjDy8/Ol\nNhaLxe5ohEajEUajUZo2mUyO7ppb88T9YObGXf3Dy5N4WmZPywu4R+bg4GCH2jl8cbxly5Zo2bJl\ng3lXf6LbGL1ej5KSEpSVlcFqtSInJwcGg8GmzcmTJ5GWloaZM2eibdu20vzw8HAcPHgQVVVVqKqq\nwsGDB6XTWERE5HqNHnGUlv7fCYVHH30US5cuRWxsLAICAmA2m/H55587dOe4SqVCfHw8kpKSIAgC\noqOjERoaiszMTOj1ehgMBqxbtw41NTVYvnw5gCvVd9asWfD19cWTTz6J2bNnAwCeeuop6UI5ERG5\nXqOFY8qUKQ3mHT161Gb6yJEjGDx48E03FBERgYiICJt5I0aMkB6/+uqrN1x20KBBGDRo0E23QURE\nTa/RwpGZmemqHERE5CFk/RwXuHIR0WKxICAgwGN+AUBERM7jcOEoLy/HW2+9hV9//RVt2rRBZWUl\nunXrhqlTpyIgIKApMxIRkRtx+FdVaWlpuOuuu5CRkYH3338fGRkZuPvuu5GWltaU+YiIyM04XDiO\nHz+OMWPGwMfHB8CVHnJHjRqFX3/9tcnCERGR+3G4cLRu3bpBx4TFxcUO3cdBRER3DoevcQwdOhSL\nFi3CoEGD0L59e5w7dw67du2y+UktERHd+RwuHEajETqdDt9//z3++OMP+Pv7Y+rUqTa92hIR0Z3P\nocIhCAJSU1MxceJEFgoioj85h65xKJVKHDp06KbdpxMR0Z3P4YvjQ4YMwYYNG2C1WpsyDxERuTmH\nr3Fs27YN58+fx5dffgk/Pz+b5959912nByMiIvfkcOGYPHlyU+YgIiIP4XDh6NatGzZv3mwzdGxU\nVBSGDRvWlPmIiMjNOFw40tLSUFxcjLi4OOk+ji1btuCDDz5AQkJCU2YkIiI34nDhyM3NxcqVK9G6\ndWsAQEhICLp27cpTWEREfzIOF4527drh8uXLUuEAgNraWrtjh9uTl5eHjIwMCIKAmJgYxMbG2jyf\nn5+PtWvX4tSpU5g2bRoiIyOl50aMGIFOnToB+L+RAYmIqHk4XDgGDBiAxYsXY/DgwdBoNDCbzdi+\nfTsGDBiAI0eOSO3s3SAoCALS09Mxb948aDQazJ49GwaDASEhIVIbrVaLhIQEbN26tcHy3t7eWLJk\nidx9IyKiJuBw4fjvf/8LAMjKymow/+pzCoUCq1atarBsYWEhdDodgoKCAABRUVHIzc21KRyBgYHS\nOoiIyH05XDhSUlJueSMWiwUajUaa1mg0KCgocHj5uro6vPLKK1CpVHj88cfx4IMP3nIWIiK6PbKH\njr0Voig2mCfnyCI1NRUBAQEoLS3F66+/jk6dOkGn09m0yc7ORnZ2NgAgOTn5toa1Lb3lJZ3Pkf1w\np7zAnZvZWdRqtccNu+xpmT0tL+BZmV1SOK5eE7nKbDY7fFEdgDQ0bVBQEHr16oXff/+9QeEwGo0w\nGo3StMlkus3U7sET94OZG6fVaj3uNfK0zJ6WF3CPzMHBwQ61c7ivqtuh1+tRUlKCsrIyWK1W5OTk\nwGAwOLRsVVUV6urqAAAVFRU4fvy4zbURIiJyLZcccahUKsTHxyMpKQmCICA6OhqhoaHIzMyEXq+H\nwWBAYWEhli5diosXL2L//v3YsGEDli9fjjNnzuD999+HUqmEIAiIjY1l4SAiakYuKRwAEBERgYiI\nCJt5144e2KVLF6xevbrBct27d8eyZcuaPB8RETnGJaeqiIjozsHCQUREsrBwEBGRLCwcREQkCwsH\nERHJwsJBRESysHAQEZEsLBxERCQLCwcREcnCwkFERLKwcBARkSwsHEREJAsLBxERycLCQUREsrBw\nEBGRLCwcREQki8sGcsrLy0NGRgYEQUBMTAxiY2Ntns/Pz8fatWtx6tQpTJs2DZGRkdJzu3btwqef\nfgoAGDZsGAYOHOiq2EREdB2XHHEIgoD09HTMmTMHK1aswJ49e1BUVGTTRqvVIiEhAf369bOZX1VV\nhU2bNmHx4sVYvHgxNm3ahKqqKlfEJiIiO1xSOAoLC6HT6RAUFAS1Wo2oqCjk5ubatAkMDMRdd90F\nhUJhMz8vLw9hYWHw9fWFr68vwsLCkJeX54rYRERkh0tOVVksFmg0Gmlao9GgoKDglpYNCAiAxWJp\n0C47OxvZ2dkAgOTkZGi12lvOW3rLSzqfI/vhTnmBOzezs6jVapduzxk8LbOn5QU8K7NLCocoig3m\nXX9kIYe9ZY1GI4xGozRtMpluef3uxBP3g5kbp9VqPe418rTMnpYXcI/MwcHBDrVzyakqjUYDs9ks\nTZvNZvj7+zu0bEBAgM2yFovF4WWJiMj5XFI49Ho9SkpKUFZWBqvVipycHBgMBoeWDQ8Px8GDB1FV\nVYWqqiocPHgQ4eHhTZyYiIhuxCWnqlQqFeLj45GUlARBEBAdHY3Q0FBkZmZCr9fDYDCgsLAQS5cu\nxcWLF7F//35s2LABy5cvh6+vL5588knMnj0bAPDUU0/B19fXFbGJiMgOl93HERERgYiICJt5I0aM\nkB536dIFq1evtrvsoEGDMGjQoCbNR0REjuGd40REJAsLBxERycLCQUREsrBwEBGRLCwcREQkCwsH\nERHJwsJBRESysHAQEZEsLrsBkMjT1Y8f6pT1OKtnYFXa505aE5E8POIgIiJZWDiIiEgWFg4iIpKF\nhYOIiGRh4SAiIllYOIiISBaX/Rw3Ly8PGRkZEAQBMTExiI2NtXm+rq4Oq1atwokTJ9CmTRtMmzYN\ngYGBKCsrw/Tp06WxcLt27YoJEya4KjYREV3HJYVDEASkp6dj3rx50Gg0mD17NgwGA0JCQqQ2O3fu\nROvWrbFy5Urs2bMH69evx/Tp0wEAOp0OS5YscUVUIiK6CZecqiosLIROp0NQUBDUajWioqKQm5tr\n02bfvn0YOHAgACAyMhJHjhyBKIquiEdERDK45IjDYrFAo9FI0xqNBgUFBTdso1Kp0KpVK1RWVgIA\nysrKMHPmTLRs2RIjR45Ez549XRGbiIjscEnhsHfkoFAoHGrj7++P1NRUtGnTBidOnMCSJUuwbNky\ntGrVyqZtdnY2srOzAQDJycnQarW3nNdZXUI4gyP74U55AWZ2ldt5j8ulVqtdur3b5Wl5Ac/K7JLC\nodFoYDabpWmz2Qx/f3+7bTQaDerr63Hp0iX4+vpCoVDAy8sLANC5c2cEBQWhpKQEer3eZnmj0Qij\n0ShNm0ymJtwj1/HE/WBm13BlZq1W61GvkaflBdwj89UfId2MS65x6PV6lJSUoKysDFarFTk5OTAY\nDDZt7r//fuzatQsAsHfvXtxzzz1QKBSoqKiAIAgAgNLSUpSUlCAoKMgVsYmIyA6XHHGoVCrEx8cj\nKSkJgiAgOjoaoaGhyMzMhF6vh8FgwKBBg7Bq1SpMnjwZvr6+mDZtGgAgPz8fGzZsgEqlglKpxPjx\n4+Hr6+uK2EREZIfL7uOIiIhARESEzbwRI0ZIj729vTFjxowGy0VGRiIyMrLJ8xERkWN45zgREcnC\nwkFERLKwcBARkSwsHEREJAsLBxERycLCQUREsrBwEBGRLCwcREQkCwsHERHJwsJBRESysHAQEZEs\nLBxERCQLCwcREcnCwkFERLKwcBARkSwuG4+DiFyvfvxQp6zHWeOtq9I+d9KaqDm5rHDk5eUhIyMD\ngiAgJiYGsbGxNs/X1dVh1apVOHHiBNq0aYNp06YhMDAQAJCVlYWdO3dCqVQiLi4O4eHhropNRETX\nccmpKkEQkJ6ejjlz5mDFihXYs2cPioqKbNrs3LkTrVu3xsqVKzFkyBCsX78eAFBUVIScnBwsX74c\nc+fORXp6ujQGORERuZ5LCkdhYSF0Oh2CgoKgVqsRFRWF3Nxcmzb79u3DwIEDAVwZLvbIkSMQRRG5\nubmIioqCl5cXAgMDodPpUFhY6IrYRERkh0sKh8VigUajkaY1Gg0sFssN26hUKrRq1QqVlZUNlg0I\nCGiwLBERuY5LrnGIothgnkKhcKiNvfn2ZGdnIzs7GwCQnJyM4ODgW0j6/32579aXbQ6elhdgZlfx\nxMxOclvfAc3EUzK75IhDo9HAbDZL02azGf7+/jdsU19fj0uXLsHX17fBshaLBQEBAQ22YTQakZyc\njOTk5CbaC/leeeWV5o4gi6flBZjZVTwts6flBTwrs0sKh16vR0lJCcrKymC1WpGTkwODwWDT5v77\n78euXbsAAHv37sU999wDhUIBg8GAnJwc1NXVoaysDCUlJejSpYsrYhMRkR0uOVWlUqkQHx+PpKQk\nCIKA6OhohIaGIjMzE3q9HgaDAYMGDcKqVaswefJk+Pr6Ytq0aQCA0NBQ9OnTBzNmzIBSqcS4ceOg\nVPK+RSKi5uKy+zgiIiIQERFhM2/EiBHSY29vb8yYMcPussOGDcOwYcOaNF9TMBqNzR1BFk/LCzCz\nq3haZk/LC3hWZoXo6NVnIiIisK8qIiKSiX1VNYGbda/iblJTU3HgwAG0bdsWy5Yta+44DjGZTEhJ\nScH58+ehUChgNBrx6KOPNnesG6qtrcWCBQtgtVpRX1+PyMhIPP30080dyyGCIOCVV15BQECAR/zy\nZ9KkSfDx8YFSqYRKpXKrX1reyMWLF7F69WqcPn0aCoUCL7zwArp169bcsW6IhcPJrnavMm/ePGg0\nGsyePRsGgwEhISHNHe2GBg4ciMGDByMlJaW5ozhMpVJh9OjR6Ny5M6qrq/HKK68gLCzMbV9nLy8v\nLFiwAD4+PrBarZg/fz7Cw8Pd+svhqq+++godO3ZEdXV1c0dx2IIFC+Dn59fcMRyWkZGB8PBwvPji\ni7Barbh8+XJzR2oUT1U5mSPdq7ibXr16wdfXt7ljyOLv74/OnTsDAFq2bImOHTu6dY8CCoUCPj4+\nAK7cp1RfX9/gJlh3ZDabceDAAcTExDR3lDvWpUuXcOzYMQwaNAgAoFar0bp162ZO1TgecTiZve5V\nCgoKmjHRna+srAwnT550+/t7BEHArFmzcPbsWTzyyCPo2rVrc0e6qQ8//BCjRo3yqKMNAEhKSgIA\n/M///I/b/1qprKwMfn5+SE1NxalTp9C5c2eMHTtW+kPDHfGIw8kc6V6FnKempgbLli3D2LFj0apV\nq+aO0yilUoklS5Zg9erV+O233/DHH380d6RG7d+/H23btpWO7DzFokWL8Oabb2LOnDnYvn078vPz\nmztSo+rr63Hy5Ek8/PDD+Oc//4kWLVpgy5YtzR2rUSwcTuZI9yrkHFarFcuWLUP//v3x0EMPNXcc\nh7Vu3Rq9evVCXl5ec0dp1PHjx7Fv3z5MmjQJb731Fo4cOYJ33nmnuWPd1NUuidq2bYsHHnjA7XvT\n1mg00GguKX0wAAAQgElEQVQ00hFoZGQkTp482cypGsfC4WSOdK9Ct08URaxevRodO3bEY4891txx\nbqqiogIXL14EcOUXVocPH0bHjh2bOVXjnn32WaxevRopKSmYNm0aevfujSlTpjR3rEbV1NRIp9Vq\nampw6NAhdOrUqZlTNa5du3bQaDQoLi4GABw+fNhtf+RxFa9xONmNuldxZ2+99Rby8/NRWVmJ559/\nHk8//bR0oc5dHT9+HLt370anTp3w8ssvAwCeeeaZBr0TuIvy8nKkpKRAEASIoog+ffrg/vvvb+5Y\nd5wLFy5g6dKlAK6cAurXr59HjBgaHx+Pd955B1arFYGBgUhISGjuSI3ineNERCQLT1UREZEsLBxE\nRCQLCwcREcnCwkFERLKwcBARkSwsHOQ0CxcuxI4dO5pl27W1tUhOTsY//vEPLF++vFkyXOvYsWOY\nOnVqc8do4OjRo3j++edvez27du3Cq6++6oREzmMymTB69GgIgtDcUe54LBwezB0/vM1l7969uHDh\nAtasWWN3JMmUlBT8+9//tplXVlaGp59+GvX19QCudMd96NAhp+Tp2bMn3n77baesixyj1Wrxr3/9\ny+lDSzvzfXG9X3/9FYsWLUJcXBzGjRuH5cuXo7y8vEm25UwsHOR2RFGU/VfjuXPn0KFDB6hUqiZK\n5XmuFsQ/y3Y90cWLF2E0GpGSkoLU1FT4+PggNTW1uWPdFO8cdyKLxYI1a9bg2LFj8PHxwZAhQ6TB\nhf73f/8XHTt2xJgxYwAAK1asQIsWLZCQkICzZ8/ivffew6lTp6BQKHDvvfdi3LhxUtfKJpMJH374\nIY4dOwZRFNG3b1888sgjSEtLg9VqxejRo6FSqfDhhx82yLRw4UL06NEDR48exalTp9CtWzdMmTIF\nfn5+OHr0KFauXInVq1dL7SdNmoSJEyciLCwMGzZsQFFREdRqNfbt24f27dvjxRdfxI8//ogvv/wS\nXl5eeP7553HvvfdKy5eWlmL27NkoLi7GPffcg4SEBKnL9l9//RUfffQRioqK0L59e4wdOxb33HOP\nlLN79+7Iz8/HiRMnsGzZMuh0Opt9KSoqwgcffIDff/8dAQEBePbZZ2EwGLBhwwZkZWUBAHJzcxEX\nFyf7zveVK1fCZDLhzTffhFKpxFNPPYXHH38c+/btw8cffwyLxYK7774bzz33nNQdxKRJk2A0GrF7\n926cP38eDzzwAJ577jl4e3s3eG23bNmCHTt24MKFC9BoNHjmmWfw4IMP2s0iCAK2bNmCb775Bhcu\nXECHDh3w8ssvQ6vVIiMjAz/99BMuXboEnU6HsWPHomfPngCADRs24PTp0/Dy8sL+/fsxZswY9O/f\nH2lpadi3bx/atWuH6Ohom21t2bIF//nPf1BdXQ1/f38899xz+Mtf/tIgU2VlJVJTU5Gfn4/g4GCb\n/3MAePrppxEfH4+vvvoK9fX1SElJuWHW2tpaxMXF4d1334Wfnx82b96MjRs3Ys2aNWjVqhX+/e9/\no6amBmPHjkVKSgpatGiBc+fO4dixYwgJCcGUKVMavDeAK0eQiYmJ+OSTT6BSqRp971+voqICqamp\n+OWXX6BQKBAaGoqFCxciJSXF7vviZu/lbt264fDhw3Y/B9e67777bKYHDx6MhQsX2n1fuBWRnKK+\nvl6cOXOmuHHjRrGurk48e/asOGnSJPHnn38WRVEUy8vLxXHjxomHDx8Wd+/eLU6aNEm8dOmSKIqi\nWFJSIh48eFCsra0VL1y4IM6fP1/MyMiQ1vvSSy+JGRkZYnV1tXj58mXx2LFjoiiK4jfffCPOmzev\n0VwLFiwQExMTxTNnzoiXL18WFyxYIK5bt04URVE8cuSIOHHiRJv2CQkJ4sGDB0VRFMXMzEzx2Wef\nFX/++WfRarWKK1euFBMSEsTNmzeLdXV14n//+18xISHBZlsTJkwQT506JVZXV4tLliwR3377bVEU\nRdFsNotxcXHi/v37xfr6evHgwYNiXFyceOHCBWnZ559/Xvzjjz9Eq9Uq1tXV2eSqq6sTExMTpW0f\nPnxYHD16tHjmzBkp69Vt2bNq1Srxk08+sZlXWloqDh8+XLRarQ32XRRF8cyZM+KoUaPEgwcPinV1\ndeKWLVvExMREKVtCQoI4Y8YM8dy5c2JlZaU4b948aRvXv7Y5OTmi2WwW6+vrxT179oijRo0SLRaL\n3ayfffaZOGPGDPHMmTOiIAjiyZMnxYqKClEURfHbb78VKyoqRKvVKn7++efic889J16+fFl6DUaO\nHCn++OOPYn19vXj58mVx3bp14quvvipWVlaK586dE2fMmCHlOnPmjPj888+LZrNZej1KSkrsZlqx\nYoW4bNkysbq6Wjx16pQ4YcIEm/fe8OHDxddff12srKyU8jSWdf78+eIPP/wgiqIoLlq0SExMTBQP\nHDggPffjjz9K/29jx44VCwoKRKvVKr799tviihUr7Ga8/v+zsff+9davXy++9957Yl1dnVhXVyfm\n5+eLgiBI/8/Xvi8ceS/f6HNwM1988YU4Z84ch9o2J56qcpLffvsNFRUVeOqpp6BWqxEUFISYmBjk\n5OQAuNKR2fjx45GSkoIPP/wQiYmJaNmyJQBAp9MhLCwMXl5e8PPzw5AhQ6SuoAsLC2GxWDB69Gj4\n+PjA29sbPXr0kJVt4MCBCA4Ohre3N/r06YPff//d4WV79OiB8PBwqFQqREZGoqKiArGxsVCr1ejb\nty/OnTsndd4HAAMGDECnTp3g4+ODkSNH4ocffoAgCNi9ezfuu+8+REREQKlUIiwsDHq9HgcOHLDJ\nGRoaCpVKBbXa9mC4oKAANTU10rZ79+6NiIgIfP/997JeCzlycnJw3333ISwsDGq1Gn//+99RW1uL\n48ePS20eeeQRaLVa+Pr64oknnsCePXvsrqtPnz4ICAiAUqlEVFQUdDrdDXtt3bFjB0aOHIng4GAo\nFArcfffdaNOmDYArr2+bNm2gUqnw97//HVarVeocDwC6deuGBx98EEqlEt7e3vjhhx8wbNgw+Pr6\nQqvV4m9/+5vUVqlUoq6uDkVFRVIfSfb+khcEAT/++CNGjBgBHx8fdOrUCX/9618btHviiSfg6+sL\nb2/vm2bt2bMn8vPzUV9fj1OnTuFvf/sb8vPzUVtbi8LCQpv3+EMPPYQuXbpApVKhX79+st6/jr73\nVSoVzp8/D5PJBLVajZ49e95wOARH3ss3+hw05tSpU9i0aRNGjRrl8P41F56qcpJz586hvLwcY8eO\nleYJgiCdRgCA+++/H2vWrEFwcLDNB+PChQvIyMjAsWPHUFNTA0EQpMNak8mE9u3b39a5+3bt2kmP\nW7RogZqaGoeXbdu2rfTY29sbfn5+0sXHq18QNTU10mm1awex0mq1qK+vR0VFBUwmE/bu3Yv9+/dL\nz9fX10uH99cve73y8nJotVqbC5/t27d3eNQ/lUoFq9VqM+/qKHw3+oIoLy9H+/btpWmlUgmtVmuz\nTa1W61Ceb7/9Fl988QXOnTsH4MprVllZabet2WxGUFCQ3ee2bt2KnTt3wmKxQKFQoLq62mY917+G\n5eXlDf5Prrp6+mjjxo0oKirCvffeizFjxkjdkl9VUVGB+vp6m/W0b98ex44ds2l3/bYby9qrVy98\n9NFHOHnyJDp16oSwsDC8++67KCgogE6nszmddDvvX0eXHTp0KDZu3Ig33ngDAGA0GhEbG2u3rdz3\n8rWfg2vzXOvs2bNYvHgx4uLibL4z3BULh5NotVoEBgY2Ol7BJ598go4dO6KsrAzff/89+vXrBwD4\n+OOPAQBLly5FmzZt8NNPP2HNmjXSek0mE+rr651+4bdFixY2YxsLgoCKiorbWue1Y5GYTCaoVCr4\n+flBo9Ggf//+jf4UtLEBr/z9/WEymSAIglQ8TCYTOnTo4FAurVaL06dP28wrKytrUIyu3+a1gy2J\nogiTyWTzxWoymWweX/+lC1z5o+K9997D/Pnz0a1bNyiVSrz88st2B/0CrnzplJaWNugO/NixY/js\ns88wf/58hISEQKlUIi4u7obrAa58cZrNZqmH5mvzAkC/fv3Qr18/XLp0Ce+//z7Wr1+PyZMn27Tx\n8/ODSqWC2WyWuoK/fj2A7f/fzbJ2794dxcXF+Omnn9CrVy+EhITAZDLhwIED6NWr1w33p6m0bNkS\nY8aMwZgxY3D69Gm89tpr0Ov1dq/3OPJevtHnwJ5z585h0aJFePLJJzFgwIDb3xkX4KkqJ+nSpQta\ntmyJLVu2oLa2FoIg4I8//pBOR+Tn52PXrl1ITExEYmIiMjIypL9Oq6ur4ePjg9atW8NisWDr1q02\n6/X398f69etRU1OD2tpa/PLLLwCufClYLJYGf0k7Kjg4GHV1dThw4ACsVis2b96Murq623odvvvu\nOxQVFeHy5cvYsGEDIiMjoVQq0b9/f+zfvx95eXkQBAG1tbU4evSozQesMV27doWPjw8+//xzWK1W\nHD16FPv370ffvn0dWv6hhx7CgQMHcPDgQQiCAIvFgk8//RRRUVFSm3bt2qGsrEyajoqKws8//4zD\nhw/DarVi69at8PLyQvfu3aU227dvh9lsRlVVFbKystCnT58G2758+TIUCoX0xfHNN980KGLXiomJ\nQWZmJkpKSiCKIk6dOoXKykpUV1dLX0CCIGDTpk24dOlSo/vdp08fZGVloaqqCmazGdu2bZOeKy4u\nxpEjR1BXVwdvb294e3vbLaJKpRIPPvggNm7ciMuXL6OoqAjffvtto9u9WdYWLVqgc+fO2L59u1Qo\nunfvjuzs7GYpHPv378fZs2chiiJatmwJpVIpvRbXvy8ceS/f6HNwPYvFgtdffx2PPPIIHn744abf\nUSfhEYeTKJVKzJo1Cx999BEmTZoEq9WK4OBgjBgxApcuXUJKSgri4+MREBCAgIAAREdHIzU1FXPn\nzsXw4cOxatUq/OMf/4BOp8OAAQPw5Zdf2qx3zZo1SEhIgEKhQN++fdGjRw/07t0bISEhGD9+PJRK\nJdLT02VlbtWqFZ577jmsXr0agiBg6NChjZ4ucsSAAQOQkpKC4uJi9OzZUxpXQKvVYubMmVi3bh3e\nfvttKJVKdOnSBePHj3dovWq1GjNnzsQHH3yArKwsBAQEIDEx0eHBkEJDQzF16lR8/PHHOHv2LFq3\nbo0+ffpg+PDhUpvY2FisWbMG69atw7BhwzB06FBMnjwZa9askX5VNWvWLJvrL/369cMbb7yB8vJy\nGAwGPPnkkw22HRISgsceewxz586FUqnEgAEDbIrP9R577DHU1dXhjTfeQGVlJTp27IiXXnoJ4eHh\nCA8Px9SpU9GiRQsMGTLE5tSTPcOHD0daWhoSExPh7++P6OhofPXVVwCAuro6rF+/HmfOnIFKpUL3\n7t0xYcIEu+sZN24cUlNTMWHCBAQHB2PgwIE4evToDbfrSNZevXrZjBXfq1cv7N27t1lO1ZSUlGDN\nmjWoqKhA69at8fDDD0unnuy9L272Xr7R5+B6O3bsQGlpKTZu3IiNGzdK8//1r3817Q7fJo7HQXSL\nrv3pMtFVCxcuRP/+/RETE9PcUZoMT1UREZEsLBxERCQLT1UREZEsPOIgIiJZWDiIiEgWFg4iIpKF\nhYOIiGRh4SAiIllYOIiISJb/B9NnVFnPweCuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10cf33278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(k_values, p)\n",
    "plt.xlabel(\"exact number of Utopia cards drawn in step 2\")\n",
    "plt.ylabel(\"probability\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Average number of cubes removed\n",
    "Expectation value: average number of supply cubes to be removed during step 2\n",
    "$$\n",
    "\\langle k \\rangle = \\sum_{k=0}^n p_k \\, k\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expectation = np.sum(p * k_values)\n",
    "expectation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Averaged over many games, only 1.2 cubes will be removed, so we will probably just put 2 cubes down.\n",
    "\n",
    "However, we can analyze the outcomes in more detail:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cumulative distributions: How bad could it be?\n",
    "\n",
    "We want to know what the probabilities are for certain events, for instance, that _at least 2 cubes are removed_ or that _no more than two cubes are removed_. \n",
    "\n",
    "We obtain these probabilities by summing the probabilities over the individual events with their probabilities $p_k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Probability to have _at a maximum_ $\\kappa$ cubes removed\n",
    "\n",
    "One might also want to know the probability that no more than $\\kappa$ cubes will be removed:\n",
    "\\begin{gather}\n",
    "P(k \\leq \\kappa) = \\sum_{k=0}^\\kappa p_k\n",
    "\\end{gather}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.2361\n",
      "1 0.6558\n",
      "2 0.9181\n",
      "3 0.9905\n",
      "4 0.9995\n",
      "5 1.0000\n",
      "6 1.0000\n"
     ]
    }
   ],
   "source": [
    "g = p.cumsum()\n",
    "print(\"\\n\".join([\"{} {:.4f}\".format(kk, pp) for \n",
    "      kk,pp in zip(k_values, g)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAENCAYAAADzFzkJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlUU3eiB/BvFhAQQUgEimtBqIOOpYgd0WpFmWq17bM+\nV6ZaRVsXcBtfcTla7etzBmutjojVUWSsdtzGp3Xp1Fe0PrVM+1ChitoK1g1BIAEVZQ339/7wmENk\n+4GQQPv9nOM5Jrm595twk2/uvcnvqoQQAkRERBLUtg5AREQtB0uDiIiksTSIiEgaS4OIiKSxNIiI\nSBpLg4iIpLE0iIhIGkuDiIiksTSIiEgaS4OIiKRpbR3gaWVlZdl0+Xq9HgaDwaYZ6ouZm15Lywsw\ns7U0h8ze3t4Nvi+3NIiISBpLg4iIpLE0iIhIGkuDiIiksTSIiEiaVb49tWHDBpw7dw6urq5YvXp1\nlduFEEhISEBKSgpatWqFmTNnwsfHxxrRiIioHqyypTFw4EAsXry4xttTUlJw584drFu3Du+++y62\nbNlijVhERFRPVimNgIAAODs713j7mTNnMGDAAKhUKvj7++Phw4coKCiwRjQiIqqHZnFMIz8/H3q9\n3nxZp9MhPz/fhomIiKg6zeIX4UKIKtepVKpqp01MTERiYiIAICYmxqJsbEGr1do8Q30xc9OzZt6c\nN/s2znwaZS6A5/6kupfFzE9NJnNTaBalodPpLH5WbzQa4ebmVu20YWFhCAsLM1+29c/xm8OQAPXF\nzE2vpeVtTC3xcf/aMj/NMCLNojSCg4Px1VdfoV+/fkhPT4eTk1ONpUHU1CreeeOp59FYnyY1mw82\n0pyIGodVSmPt2rW4dOkSCgsLMX36dIwZMwYmkwkA8Morr+CFF17AuXPnMHv2bNjb22PmzJnWiEVE\nRPVkldKYO3durberVCpMnTrVGlGIiOgpNItvTxERUcvA0iAiImksDSIiksbSICIiaSwNIiKSxtIg\nIiJpLA0iIpLG0iAiImksDSIiksbSICIiaSwNIiKSxtIgIiJpzWJodPrlaoxhxgEONU7UXHBLg4iI\npLE0iIhIGkuDiIiksTSIiEgaS4OIiKSxNIiISBpLg4iIpLE0iIhIGkuDiIiksTSIiEgaS4OIiKSx\nNIiISBpLg4iIpLE0iIhIGkuDiIiksTSIiEgaS4OIiKSxNIiISBpLg4iIpFntHOGpqalISEiAoigY\nPHgwRowYYXG7wWBAXFwcHj58CEVREB4ejqCgIGvFIyIiCVYpDUVREB8fjyVLlkCn02HRokUIDg5G\nhw4dzNPs27cPISEheOWVV5CZmYk///nPLA0iombGKrunMjIy4OXlBU9PT2i1WvTt2xfJyckW06hU\nKhQVFQEAioqK4ObmZo1oRERUD1bZ0sjPz4dOpzNf1ul0SE9Pt5hm9OjR+K//+i989dVXKC0txdKl\nS60RjYiI6sEqpSGEqHKdSqWyuPztt99i4MCBeP3113HlyhXExsZi9erVUKstN4YSExORmJgIAIiJ\niYFer2+64BK0Wq3NM9SXNTPnWGUp8mQed3PK3NLyAsxsLbZ637FKaeh0OhiNRvNlo9FYZffT8ePH\nsXjxYgCAv78/ysvLUVhYCFdXV4vpwsLCEBYWZr5sMBiaMHnd9Hq9zTPUV0vM3Fha2uNuaXkBZraW\np8ns7e3d4Pta5ZiGr68vsrOzkZubC5PJhKSkJAQHB1tMo9frkZaWBgDIzMxEeXk5XFxcrBGPiIgk\nWWVLQ6PRICIiAitWrICiKAgNDUXHjh2xe/du+Pr6Ijg4GBMnTsSmTZtw5MgRAMDMmTOr7MIiIiLb\nstrvNIKCgqp8hXbs2LHm/3fo0AEffvihteIQEVED8BfhREQkjaVBRETSWBpERCSNpUFERNJYGkRE\nJI2lQURE0lgaREQkjaVBRETSWBpERCSNpUFERNJYGkREJI2lQURE0lgaREQkjaVBRETSWBpERCSN\npUFERNJYGkREJI2lQURE0lgaREQkjaVBRETSWBpERCSNpUFERNJYGkREJE26NAoLC5syBxERtQBa\n2QlnzJiBnj17YsCAAQgODoZWK31XIiL6hZB+59+wYQNOnz6NL774Aps2bUKfPn3w8ssvo1u3bk2Z\njyqpeOeNRplPTqPMBdBsPthIcyKilkK6NFxcXDBs2DAMGzYMWVlZOHnyJGJjY6FSqdC/f38MGjQI\n7dq1a8qsRERkYw06EH737l3cvXsXxcXF8PT0RH5+PqKjo3HgwIHGzkdERM2I9JbGrVu3cOrUKZw6\ndQoODg54+eWX8fHHH8Pd3R0A8O///u947733MGLEiCYLS0REtiVdGsuWLUO/fv0wf/58dO3atcrt\nHh4eGDZsWKOGIyKi5kW6NP7jP/4DAQEBVa7PyMgwl8jYsWMbLxkRETU70sc0Vq5cWe31K1asaLQw\nRETUvNW5paEoCgBACGH+91hOTg40Gk3TpSMiomalztIYP368+f/jxo2zuE2tVuPNN9+UWlBqaioS\nEhKgKAoGDx5c7QHzpKQk7N27FyqVCp07d8acOXOk5k1ERNZRZ2msX78eQggsX74cH3zwgfl6lUoF\nFxcX2Nvb17kQRVEQHx+PJUuWQKfTYdGiRQgODkaHDh3M02RnZ+PAgQP48MMP4ezsjHv37jXwIRER\nUVOpszQe/2Bvw4YNDV5IRkYGvLy84OnpCQDo27cvkpOTLUrj2LFjGDJkCJydnQEArq6uDV4eERE1\njVpLY9OmTZg2bRqAR1scNYmKiqp1Ifn5+dDpdObLOp0O6enpFtNkZWUBAJYuXQpFUTB69GgEBgbW\nnp6IiKyq1tLw8PAw///xVkJDVD54/phKpbK4rCgKsrOzsWzZMuTn5+P999/H6tWr0bp1a4vpEhMT\nkZiYCACIiYmBXq9vcK7GoNVqrZahscaMaiwyj5uZn05Lywsws7XY6r2v1tKofJB79OjRDV6ITqeD\n0Wg0XzYajXBzc7OYxt3dHf7+/tBqtfDw8IC3tzeys7Or/JAwLCwMYWFh5ssGg6HBuRqDXq+3eQZb\naYmPu6Vlbml5AWa2lqfJ7O3t3eD71loaaWlpUjPp0aNHrbf7+voiOzsbubm5cHd3R1JSEmbPnm0x\nzYsvvojTp09j4MCBuH//PrKzs59q64aIiBpfraXx6aef1jkDlUpV6/EOANBoNIiIiMCKFSugKApC\nQ0PRsWNH7N69G76+vggODsbzzz+PH374AfPmzYNarcZbb72FNm3a1O/REBFRk6q1NOLi4hptQUFB\nQQgKCrK4rvKwIyqVCm+//TbefvvtRlsmERE1Lp4jnIiIpNW6pTFv3jysWbMGwKPTvdZEZjcWERG1\nfLWWxuPfaADArFmzmjwMERE1b7WWRuXzf1c3LDoREf26SJ9Pw2QyYd++ffj2229RUFAANzc39O3b\nFyNHjpQaf4qIiFo+6dLYvHkzsrKyMHnyZLRr1w55eXk4cOAAtmzZgpkzZzZlRiIiaiakSyM5ORmx\nsbHmYT06dOgAPz8/HusgIvoVkf7Kbdu2bVFaWmpxXVlZWZXhQIiI6JdLehiRAQMG4E9/+hOGDh1q\nHkvq6NGjGDBgQJOHJCKi5qHew4js37/f4nJiYmK1Z+EjIqJfHqsNI0JERC0fhxEhIiJp0t+eKioq\nwt69e3Hp0iUUFhZanFiJw4gQEf06SG9pbNmyBdeuXcOoUaPw4MEDREREQK/XY/jw4U2Zj4iImhHp\n0jh//jzmz5+P3r17Q61Wo3fv3pg3bx5OnTrVlPmIiKgZkS4NIQScnJwAAA4ODnj48CHatm2LO3fu\nNFk4IiJqXqSPaXTu3BmXLl3Cb3/7W3Tr1g3x8fFwcHDAM88805T5iIioGZHe0pg2bRratWsHAIiI\niICdnR0ePnyIqKioJgtHRETNi/SWhqenp/n/Li4utZ6UiYiIfpmkSwMAjh8/bjE0er9+/RAaGgqV\nStVU+YiIqBmRLo0dO3YgOTkZw4cPh16vh8FgwKFDh5CVlYW33nqrKTMSEVEzIV0aJ06cwMqVK6HT\n6czXBQUFYcGCBSwNIqJfCekD4Y6OjnB0dKxy3eOv4RIR0S9frVsaOTk55v8PGzYMH3/8MUaMGAF3\nd3cYjUYcPHiQvwgnIvoVqbU0Zs+eXeW6ixcvWlxOS0vD0KFDGzcVERE1S7WWxu7du62Vg4iIWoB6\nfeUWAAwGA/Lz8+Hu7g69Xt8UmYiIqJmSLo2CggKsXbsWV65cQZs2bVBYWAh/f3/MmTMH7u7uTZmR\niIiaCelvT23evBmdO3dGQkIC/vrXvyIhIQFdunTB5s2bmzIfERE1I9Kl8dNPP2HixIlwcHAA8Gik\n27feegtXrlxpsnBERNS8SJdG69atkZmZaXFdVlYWf6dBRPQrIn1M44033sCHH36IQYMGoV27dsjL\ny8OJEycwduzYpsxHRETNiPSWRlhYGObNm4fCwkKcPXsWhYWFmDNnDsLCwqTun5qaijlz5mDWrFk4\ncOBAjdN99913GDNmDK5evSobjYiIrERqS0NRFGzYsAHTpk1Djx496r0QRVEQHx+PJUuWQKfTYdGi\nRQgODkaHDh0spisuLsY///lP+Pn51XsZRETU9KS2NNRqNc6fP9/gIdAzMjLg5eUFT09PaLVa9O3b\nF8nJyVWm2717N9544w3Y2dk1aDlERNS0pHdPDR8+HHv27IHJZKr3QvLz8y1Gx9XpdMjPz7eY5tq1\nazAYDOjVq1e9509ERNYhfSD8q6++wt27d3HkyBG4uLhY3Pbpp5/Wel8hRJXrKm+1KIqCbdu2YebM\nmXXmSExMRGJiIgAgJibG5r9K12q1VsuQU/ckViXzuJn56bS0vAAzW4ut3vukS2PWrFkNXohOp4PR\naDRfNhqNcHNzM18uKSnBrVu38MEHHwAA7t69i48++gjR0dHw9fW1mFdYWJjFwXeDwdDgXI3h8Qmp\nfo1a4uNuaZlbWl6Ama3laTJ7e3s3+L7SpeHv7499+/ZZnO61b9++GDlyZJ339fX1RXZ2NnJzc+Hu\n7o6kpCSLEXSdnJwQHx9vvrx8+XJMmDChSmEQEZFtSZfG5s2bkZWVhcmTJ5t/p3HgwAFs2bKlzt1K\nGo0GERERWLFiBRRFQWhoKDp27Ijdu3fD19cXwcHBT/1AiIio6UmXRnJyMmJjY9G6dWsAQIcOHeDn\n5ye92yooKAhBQUEW19X0w8Dly5fLxiIiIiuS/vZU27ZtUVpaanFdWVmZxbEJIiL6ZZPe0hgwYAD+\n9Kc/YejQoeYD20ePHsWAAQOQlpZmnq4hP/4jIqKWQbo0vv76awDA/v37q1z/+DaVSoX169c3Yjwi\nImpOpEsjLi6uKXMQEVELIH1Mg4iIiKVBRETSWBpERCSNpUFERNJYGkREJI2lQURE0lgaREQkjaVB\nRETSWBpERCSNpUFERNJYGkREJI2lQURE0lgaREQkjaVBRETSWBpERCSNpUFERNJYGkREJI2lQURE\n0lgaREQkjaVBRETSWBpERCSNpUFERNJYGkREJI2lQURE0lgaREQkTWvrALZS8c4bjTKfnEaZC6DZ\nfLCR5kRE1HS4pUFERNJYGkREJI2lQURE0qx2TCM1NRUJCQlQFAWDBw/GiBEjLG4/fPgwjh07Bo1G\nAxcXF8yYMQPt2rWzVjwiIpJglS0NRVEQHx+PxYsXY82aNfj222+RmZlpMU2XLl0QExODjz/+GH36\n9MGOHTusEY2IiOrBKqWRkZEBLy8veHp6QqvVom/fvkhOTraYpkePHmjVqhUAwM/PD/n5+daIRkRE\n9WCV3VP5+fnQ6XTmyzqdDunp6TVOf/z4cQQGBlZ7W2JiIhITEwEAMTEx0Ov1DcrUWF+VbSwyj4OZ\nn15Ly9zS8gLMbC0Nfe97WlYpDSFEletUKlW10548eRI///wzli9fXu3tYWFhCAsLM182GAyNktHW\nWuLjYOam19LyAsxsLU+T2dvbu8H3tcruKZ1OB6PRaL5sNBrh5uZWZbrz589j//79iI6Ohp2dnTWi\nERFRPVilNHx9fZGdnY3c3FyYTCYkJSUhODjYYppr165h8+bNiI6OhqurqzViERFRPVll95RGo0FE\nRARWrFgBRVEQGhqKjh07Yvfu3fD19UVwcDB27NiBkpISfPLJJwAe7a9bsGCBNeIREZEkq/1OIygo\nCEFBQRbXjR071vz/pUuXWisKERE1EH8RTkRE0lgaREQkjaVBRETSWBpERCSNpUFERNJYGkREJI2l\nQURE0lgaREQkjaVBRETSWBpERCSNpUFERNJYGkREJI2lQURE0lgaREQkjaVBRETSWBpERCSNpUFE\nRNJYGkREJI2lQURE0lgaREQkjaVBRETSWBpERCSNpUFERNJYGkREJI2lQURE0lgaREQkjaVBRETS\nWBpERCSNpUFERNJYGkREJI2lQURE0lgaREQkTWutBaWmpiIhIQGKomDw4MEYMWKExe3l5eVYv349\nfv75Z7Rp0wZz586Fh4eHteIREZEEq2xpKIqC+Ph4LF68GGvWrMG3336LzMxMi2mOHz+O1q1bIzY2\nFsOHD8fnn39ujWhERFQPVimNjIwMeHl5wdPTE1qtFn379kVycrLFNGfOnMHAgQMBAH369EFaWhqE\nENaIR0REkqxSGvn5+dDpdObLOp0O+fn5NU6j0Wjg5OSEwsJCa8QjIiJJVjmmUd0Wg0qlqvc0AJCY\nmIjExEQAQExMDLy9vRsW6siZht3PlpjZOlpa5paWF2DmFswqWxo6nQ5Go9F82Wg0ws3NrcZpKioq\nUFRUBGdn5yrzCgsLQ0xMDGJiYpo2tKSFCxfaOkK9MXPTa2l5AWa2lpaYuTKrlIavry+ys7ORm5sL\nk8mEpKQkBAcHW0zTq1cvnDhxAgDw3XffoXv37tVuaRARke1YZfeURqNBREQEVqxYAUVREBoaio4d\nO2L37t3w9fVFcHAwBg0ahPXr12PWrFlwdnbG3LlzrRGNiIjqwWq/0wgKCkJQUJDFdWPHjjX/397e\nHn/84x+tFafRhIWF2TpCvTFz02tpeQFmtpaWmLkyleD3WomISBKHESEiImlW2z31S1TX0CjNzYYN\nG3Du3Dm4urpi9erVto5TJ4PBgLi4ONy9excqlQphYWEYNmyYrWPVqqysDMuWLYPJZEJFRQX69OmD\nMWPG2DqWFEVRsHDhQri7u7eIb/hERkbCwcEBarUaGo2m2XyjsiYPHz7Exo0bcevWLahUKsyYMQP+\n/v62jlVvLI0Gejw0ypIlS6DT6bBo0SIEBwejQ4cOto5Wo4EDB2Lo0KGIi4uzdRQpGo0GEyZMgI+P\nD4qLi7Fw4UL07NmzWT/HdnZ2WLZsGRwcHGAymfD+++8jMDCwRbw5fPnll2jfvj2Ki4ttHUXasmXL\n4OLiYusYUhISEhAYGIj58+fDZDKhtLTU1pEahLunGkhmaJTmJiAgoNrfvjRXbm5u8PHxAQA4Ojqi\nffv2VUYSaG5UKhUcHBwAPPq9UUVFRYv46rjRaMS5c+cwePBgW0f5RSoqKsLly5cxaNAgAIBWq0Xr\n1q1tnKphuKXRQNUNjZKenm7DRL9subm5uHbtGrp27WrrKHVSFAULFizAnTt3MGTIEPj5+dk6Up3+\n9re/4a233mpRWxkAsGLFCgDA73//+2b9raTc3Fy4uLhgw4YNuHHjBnx8fDBp0iTzB4yWhFsaDSQ7\n7Ak9vZKSEqxevRqTJk2Ck5OTrePUSa1WY9WqVdi4cSOuXr2Kmzdv2jpSrc6ePQtXV1fzVl1L8eGH\nH2LlypVYvHgxjh49ikuXLtk6Uo0qKipw7do1vPLKK/joo4/QqlUrHDhwwNaxGoSl0UAyQ6PQ0zOZ\nTFi9ejX69++P3/3ud7aOUy+tW7dGQEAAUlNTbR2lVj/99BPOnDmDyMhIrF27FmlpaVi3bp2tY9XJ\n3d0dAODq6orevXsjIyPDxolqptPpoNPpzFudffr0wbVr12ycqmFYGg0kMzQKPR0hBDZu3Ij27dvj\ntddes3UcKffv38fDhw8BPPom1YULF9C+fXsbp6pdeHg4Nm7ciLi4OMydOxc9evTA7NmzbR2rViUl\nJeZdaSUlJTh//jw6depk41Q1a9u2LXQ6HbKysgAAFy5caNZf6KgNj2k0UE1DozRna9euxaVLl1BY\nWIjp06djzJgx5gNzzdFPP/2EkydPolOnTnjvvfcAAOPHj68yskBzUlBQgLi4OCiKAiEEQkJC0KtX\nL1vH+sW5d+8ePv74YwCPdv289NJLCAwMtHGq2kVERGDdunUwmUzw8PDAzJkzbR2pQfiLcCIiksbd\nU0REJI2lQURE0lgaREQkjaVBRETSWBpERCSNpdFM5ObmYsyYMaioqLDJ8pcvX45jx45Ve5vBYMCE\nCROgKIqVU8nbs2ePTX+QtmvXLkyZMgXvvPOOzTJUNmHCBOTk5Ng6RhVjxozBnTt3nmoetn6t1OSP\nf/wjLl68aOsYTY6lYSORkZE4f/68TZZd3zdYvV6P7du3Q63m6lIdg8GAQ4cO4ZNPPsHmzZur3H7i\nxAksXbq0yvWV14G4uDjs2rWr0TJt374dnp6ejTY/qtsnn3yC7t27N+o8G3u9qOzevXtYu3Ytpk2b\nhrfffhtLly6VGj+P7wKNoLl94qGnU9+/p8FgQJs2beDq6tpEiVoeW70mhBDNeou4OSkpKUHXrl0R\nExODhIQEvPzyy4iJiUFJSUmt92v2vwiPjIzEkCFDcPLkSeTl5SEwMBCRkZGwt7cHACQmJuKLL77A\ngwcP0K1bN7zzzjvmMWkqy83NRVRUFGbMmIE9e/agpKQE48ePh4+PDzZu3AiDwYD+/ftjypQpAB6N\nVLp//34cO3YMZWVlCAwMREREBJycnMzzmj59Ovbu3QsPDw988MEHuHLlCj777DNkZmaiXbt2mDRp\nUrWfPGJjY2EwGLBy5Uqo1WqMGjUKISEhAIBTp05h9+7dKCsrw/DhwzFy5EgAj4ZiT0hIwO3bt2Fv\nb4/f/e53ePvtt6HVPvoTjhkzBlOnTsXhw4dRWFiIfv36YcqUKVUGUUxNTcX+/fsBAMnJyfDy8sKq\nVasAAHl5eVi6dClu3LgBf39/zJ49Gy4uLubHu3PnTmg0Gpw4cQL/+Mc/cP/+fbRp0wbjxo1D//79\nqzzOPXv2IDMzE/b29vi///s/6PV6REZGwtfX15x53bp18PLyAvDoU5VOp8O4ceNw8eJFxMbG4tVX\nX8WhQ4egVqsxdepUaLVabNu2Dffv38frr79ufn4AoLy8HGvWrEFKSgqeeeYZzJgxA126dAHwaFTi\nrVu34vLly3BwcMDw4cPNJ3Tas2cPbt26BTs7O5w9exYTJ06sMkR4UVERtm7dipSUFLRq1QqDBw/G\nm2++ibS0NKxcuRImkwkTJkxAnz59EBkZWf3KXIPExEScPn0aAHDkyBF0794dCxcuRGZmJrZs2YLr\n16/D3d0d4eHh5qFq4uLiYGdnh5ycHKSnp+PZZ59FVFQU2rVrV+W5PXfuHHbt2oWcnBw4OTkhNDS0\n1hNDJScnY8+ePeaRWadMmYLAwEB88803OHjwIIxGI1xcXPBv//Zv+P3vfw8A5r/X0KFDceTIEfTs\n2ROzZs3CwYMHcfjwYahUKowdO9ZiOefOncP27dthNBrh6OiI4cOH44033qiSR1EU7NixA//7v/8L\nR0fHKkPKLF++HM899xwuXbqEn3/+GatXr8bly5drzLps2TK8+uqr6NOnD3788Ue8//77WLhwIYKC\ngnD+/Hls374dq1atwokTJ3Ds2DH4+fnhm2++gZOTE6ZOnYoXXnih2uctMjIS06ZNQ8+ePetc9ysT\nQmDbtm04ffo0ysvL0a5dO8yePRtXrlypdr2QWZfVanW1r4PKPD09LZ7LsLAwbN++HVlZWbUPXima\nuZkzZ4qFCxcKo9EoCgsLxdy5c8XRo0eFEEJcuHBBREREiKtXr4qysjIRHx8v3n///Wrnk5OTI0aP\nHi02bdokSktLRWpqqggPDxcrV64Ud+/eFUajUUyZMkVcvHhRCCHEsWPHRFRUlLhz544oLi4Wq1at\nEuvWrbOYV2xsrCguLhalpaXCaDSKyZMni7Nnz4qKigrxww8/iMmTJ4t79+7V+Lh++OGHKvk+/fRT\nUVpaKq5duybGjx8vbt26JYQQ4urVq+Knn34SJpNJ5OTkiLlz54rDhw+b7z969Gjx5z//WTx48EDk\n5eWJiIgIkZKSUu2yd+/eLf7yl79YXLds2TIRFRUlbt++LUpLS8WyZcvEjh07LLKZTCZRXFwsJk6c\nKG7fvi2EECI/P1/cvHmzxuWEh4ebn5PPP/9cLF682CJzdna2+fL69evFzp07hRBCpKWlibFjx4q9\ne/eK8vJy8fXXX4uIiAixdu1aUVRUJG7evCnCw8PFnTt3zMsaN26c+Ne//iXKy8vFF198IWbOnCnK\ny8tFRUWFiI6ONs/rzp07IjIy0vz8PL7v999/LyoqKkRpaWmVxxIbGytWrlwpioqKRE5Ojpg9e7Y4\nduyYOeu0adOqfQ6EEOKbb74RS5YsqXJ95XWg8mMXQojy8nIRFRUl9u3bJ8rLy8WFCxfEhAkTzM/7\n+vXrxYQJE8TFixdFWVmZ2Lp1q8UyKj+3aWlp4saNG6KiokJcv35dTJ06VXz//ffVZk1PTxcTJ04U\nP/zwg6ioqBBGo1FkZmYKIYQ4e/asyM7OFoqiiIsXL4o//OEP4urVqxZ/r+3bt4uysjJRWloqUlJS\nxNSpU8WNGzdEcXGxWLt2rUWud955R1y6dEkIIURhYaF5Xk86evSomDNnjsjLyxOFhYVi+fLl5vVR\niEfr7vTp08XNmzeFyWQS5eXltWbdtWuXiI+PF0IIsW/fPhEVFSW2b99uvm3r1q3mv9u4cePE119/\nLSoqKsTRo0fFu+++KxRFqTZn5b9nXet+ZSkpKSI6Olo8ePBAKIoibt26JfLz881/58rrhey6XN3r\noC7Xrl0T4eHh4uHDh7VO1yJ2T7366qtwd3eHs7MzevXqhevXrwN49Kk8NDQUPj4+sLOzQ3h4OK5c\nuYLc3NyxWV18AAAJ+klEQVQa5zVq1CjY29vj+eefR6tWrfDSSy/B1dUV7u7u6Natm3nkydOnT+O1\n116Dp6cnHBwcEB4ejqSkJIvN7tGjR8PBwQH29vY4efIkXnjhBQQFBUGtVqNnz57w9fXFuXPn6vVY\nR48eDXt7e3Tp0gWdO3fGjRs3AAA+Pj7w9/eHRqOBh4cHwsLCqgwFPWLECLRu3Rp6vR7du3c3P0+y\nBg4cCG9vb9jb2yMkJKTG+6tUKty8eRNlZWVwc3Ordcytbt26mZ+TAQMG1CuTRqPByJEjodVq0a9f\nPxQWFmLYsGFwdHREx44d0aFDB/PzAzx6jvr06QOtVovXXnsN5eXlSE9Px9WrV3H//n2MGjUKWq0W\nnp6eGDx4MJKSksz39ff3x4svvgi1Wm3ein1MURQkJSUhPDwcjo6O8PDwwGuvvYaTJ09KP5b6Sk9P\nR0lJCUaMGAGtVosePXogKCjI/MkTAIKCghAQEAA7OzuMHz8eV65cgcFgqDKv7t27o1OnTlCr1ejc\nuTP69etX4zDix48fR2hoKHr27Am1Wg13d3fzgItBQUHw8vKCSqVCQEAAevbsiR9//NF8X5VKhTFj\nxsDOzg729vZISkrCwIED0alTJzg4OGD06NEWy9JoNMjMzERRURGcnZ1r/HT7r3/9C8OGDYNer4ez\ns3O1p1UeOHAgOnbsCI1GA61WW2vWgIAA8+O/fPkyRowYgcuXLwMALl26hICAAPN89Xo9wsLCoFar\n8fLLL6OgoAD37t2rNueTZNd9rVaLkpIS3L59G0IIdOjQocYRs2XW5ZpeB7UpKipCbGwsRo0aVefp\nB5r97ing0QiRj9nb25vP3lZQUIBnn33WfJuDgwOcnZ2Rn58PDw+PaudVeb+zvb19lcuP9+cVFBSY\nN/WBRytPRUWFxQpT+SRMBoMB3333Hc6ePWu+rqKiot4Hxio/1latWpnzZGVl4bPPPsPVq1dRVlaG\nioqKKi+ymu77tMuuzMHBAXPnzsWhQ4ewceNGPPfcc5g4cWKNI7k++fyWl5ejoqICGo2mzjxt2rQx\nH3x//EZe098LsPx7qNVq6HQ6FBQUAHj095w0aZL5dkVR8Jvf/Kba+z7p/v37MJlM0Ov15uvatWsn\nfRZBtVpd7T5+k8lk3r34pIKCAuj1eosvHzy5zMqZH6/7j+9XWXp6Ov7+97/j5s2bMJlMMJlM6NOn\nT7XLNRqNNe5+SUlJwT/+8Q9kZWVBCIHS0lKLkWVdXFwsCregoMBiHa38egKA+fPn47//+7/x97//\nHZ06dcIf/vCHak+L++RjenI+Tz4XdWX19/dHdnY27t69i+vXr2PBggXYs2cP7t+/j4yMDIv14snX\nBADp15Xsut+jRw8MGTIE8fHxMBgMePHFFzFhwoRq37zz8vLqtS4/+TqoTllZGVauXAk/Pz+8+eab\ndT6uFlEaNXFzc7P4ZFVSUoIHDx5Ue0yjIfPOy8szXzYYDNBoNHB1dTWfR6Py8QKdTof+/ftj+vTp\nT73s6mzZsgVdunTBnDlz4OjoiCNHjuC7775r0Lye9mRRgYGBCAwMRFlZGXbt2oVNmzbhP//zP+s9\nn1atWlmcJ/nu3bu1vnnXpfL5TRRFMZ/j5PHWWUO/kuvi4gKNRgODwWAeztpgMEivZ3q9HgaDAUII\n83NfWlqK+/fvm98Mn/ybPF63FUUxF4fBYMAzzzxT7eN9vO5X9wl13bp1GDJkCBYtWgR7e3v87W9/\nw/3796vNqtPpqv1KbHl5OVavXo2oqCgEBwdDq9Xio48+spimusdQOeOTW0Fdu3ZFdHQ0TCYTvvrq\nK6xZswaffvpplWU/+Tqvbmuq8rLrytqqVSv4+Pjgyy+/RKdOnaDVavHcc8/h8OHD8PLyssk5x4cN\nG4Zhw4bh3r17WLNmDQ4ePIhx48ZVeU71en2d63JNr4PqlJeXY9WqVXB3d8e7774rlbVF7J6qyUsv\nvYRvvvkG169fR3l5OXbu3ImuXbvWuJVRH/369cORI0eQm5uLkpIS7Ny5EyEhITV+Qu7fvz/Onj2L\n1NRUKIqCsrIyXLx40eIPWFnbtm1r3Y32pOLiYjg5OcHBwQG3b9/G//zP/zTocQGPPgHl5eU16Fsm\nd+/exZkzZ1BSUgKtVgsHB4cGfxW3S5cuOH36NBRFQWpq6lOfee3nn3/G999/j4qKCnz55Zews7OD\nn58funbtCkdHRxw4cABlZWVQFAU3b96UPmmPWq1GSEgIdu7cieLiYuTl5eHw4cPVHvyvjp+fH+zs\n7MzLLykpweeffw4fHx/zp2ZXV1eL31X4+fnBwcEBBw8ehMlkwsWLF3H27Fn069fPPE1KSgp+/PFH\nmEwm7Nq1C35+flW2MoBH646zszPs7e2RkZFhsYvrSYMGDcKJEydw4cIFKIqC/Px83L59GyaTCeXl\n5eYCTUlJqfMr4yEhIThx4gQyMzNRWlqKvXv3mm8zmUw4deoUioqKoNVq4eTkVON6FBISgn/+858w\nGo148OBBnWe8k8n6m9/8BkePHjV/Qg8ICLC4bE0ZGRlIT0+HyWRCq1atYGdnZ34unlwvZNblml4H\nT3p8gjM7OztERUVJv45b9JbGb3/7W4wdOxarV6/GgwcP8Nxzz2Hu3LmNMu/Q0FAUFBRg2bJlKCsr\nw/PPP4+IiIgap9fr9YiOjsaOHTvwl7/8BWq1Gl27dq3xx14jRozA1q1bsWPHDowcObLG3QWPTZgw\nAX/961/xxRdf4Nlnn0Xfvn2RlpbWoMcWEhKCU6dOYcqUKfDw8MDKlSul7yuEwKFDhxAbGwuVSoUu\nXbpg6tSpDcoxadIkxMXF4ejRo+jduzd69+7doPk8FhwcjKSkJMTFxcHLywvz58837/5ZsGABPvvs\nM0RGRsJkMsHb27vKt3lqExERga1btyIqKgr29vYYPHgwQkNDpe5rZ2eHRYsWYdu2bTh8+DDUajW6\ndeuGefPmmT9JDho0CJ988gkmTZqEgIAAREdHIzo6Glu2bMH+/fvh7u6OqKgoi92A/fr1w969e3Hl\nyhX4+PjUeOKkqVOn4rPPPsPWrVsREBCAkJAQ84mintS1a1fMnDkT27ZtQ25uLlxdXTFlyhS0b98e\nkydPxpo1a1BeXo5evXrVedKxF154AcOHD8cHH3wAtVqNsWPHWhTWyZMnsXXrViiKAm9vb8yaNava\n+QwePBhZWVl477334OjoiNdff73Wdd/R0bHOrAEBAThw4ID5+EVAQACKi4stjmdYS3FxMbZt24ac\nnBzz8dbH3yKrbr2oa12u7XVQ2ZUrV3Du3DnY29tb7O5avHhxreXJ82kQtUCVv55M9NiePXtw586d\nJj3zYovePUVERNbF0iAiImncPUVERNK4pUFERNJYGkREJI2lQURE0lgaREQkjaVBRETSWBpERCTt\n/wFisHpGMZHtdgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10cf33320>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(k_values, g)\n",
    "plt.xlabel(\"no more than this number of Utopia cards drawn in step 2\")\n",
    "plt.ylabel(\"probability\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The chance of Utopia not getting hit is only 24%. The probability of getting hit not more often than twice is 92% – this could mean it gets hit not at all, once, or twice. So this analysis suggests that a prudent course of action is to put 2 supply cubes on Utopia at the beginning of the game so that even when it is hit twice, no incident occurs. \n",
    "\n",
    "This analysis confirmed our hunch from looking at the expectation value where we already thought that 2 cubes might be a good strategy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Probability to have _at least_ $\\kappa$ supply cubes removed\n",
    "You can also look at the previous analysis for the question of \"how bad could it get?\", i.e., the probability that _at least a certain number of cubes are removed_.\n",
    "\n",
    "The probability that _at least_ $\\kappa$ supply cubes are removed is\n",
    "\n",
    "$$\n",
    "P(k \\ge \\kappa) = \\sum_{k=\\kappa}^6 p_k = 1 - \\sum_{k=0}^{\\kappa-1} p_k\n",
    "$$\n",
    "and probability to have 0 supply cubes removed\n",
    "$$\n",
    "P(k=0) = p_0\n",
    "$$\n",
    "which I added to the table because one would like to know \n",
    "1. What is the probability to have no cubes removed?\n",
    "2. What is the probability to have at least $k$ cubes removed?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.2361\n",
      "1 0.7639\n",
      "2 0.3442\n",
      "3 0.0819\n",
      "4 0.0095\n",
      "5 0.0005\n",
      "6 0.0000\n"
     ]
    }
   ],
   "source": [
    "f = 1 - p.cumsum()\n",
    "f = np.concatenate([p[:1], f[:-1]])  # prepend k = 0\n",
    "print(\"\\n\".join([\"{} {:.4f}\".format(kk, pp) for \n",
    "      kk,pp in zip(k_values, f)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAENCAYAAAAykHOlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8DHe/B/DPXnIjEskuiUi0R0Q1POqJcIhWRVJVtE0d\nl9aDElRLUPW41q2nJ4TWpSKtloi4PI+4nCraoqGopJe4pBrRSooSCbEbBLnIZub84WWOldtKMrO7\nfN6vl9fLzP525rOT2f3u/HbmNypRFEUQERHVM7W1AxAR0aOJBYaIiGTBAkNERLJggSEiIlmwwBAR\nkSxYYIiISBYsMEREJAsWGCIikgULDBERyYIFhoiIZKG1doC6ys3Nter69Xo9DAaDVTM8LGaWn73l\nBZhZKbaQ2cfHR5H18AiGiIhkwQJDRESyYIEhIiJZsMAQEZEsWGCIiEgWip1Flp6ejoSEBAiCgLCw\nMERERJg9bjAYEBcXh9u3b0MQBAwZMgRBQUFKxSMionqmSIERBAHx8fGYPXs2dDodZs6cieDgYPj6\n+kpttm/fjq5du6JXr17IycnBwoULWWCIiOyYIl1k2dnZ8Pb2hpeXF7RaLUJCQpCWlmbWRqVSoaio\nCABQVFQEDw8PJaIREZFMFDmCKSgogE6nk6Z1Oh2ysrLM2gwcOBD/8z//gz179qC0tBRz5sxRIhoR\nEclEkQIjimKFeSqVymw6JSUFPXr0wMsvv4wzZ84gNjYWS5YsgVptfpCVnJyM5ORkAEBMTAz0er18\nwS2g1WplyfDjjz8CALp27Vrvy5Yrs5zsLbO95QWYWSn2mLm2FCkwOp0ORqNRmjYajRW6wA4cOIBZ\ns2YBAFq3bo2ysjLcvHkT7u7uZu3Cw8MRHh4uTVt7yAW5hn0ICAgAIM/rs4WhKh6WvWW2t7wAMyvF\nFjI/UkPF+Pv7Iy8vD/n5+TCZTEhNTUVwcLBZG71ej4yMDABATk4OysrK4ObmpkQ8m5SWllbhdyoi\nInuiyBGMRqNBZGQkoqOjIQgCQkND4efnh6SkJPj7+yM4OBjDhw/H559/jq+//hoAMG7cuArdaI+T\nRYsWAQC2bdtm5SRERLWj2HUwQUFBFU47Hjx4sPR/X19ffPjhh0rFISIimfFKfiIikgULDBERyYIF\nhoiIZGH3d7R8VM2fP9/aEYiI6oQFxka1a9euwrzyMa/Uy7Kv1MtSAM3qnfW0JCJ6FLGLzEYdPnwY\nhw8ftnYMIqJa4xGMjVqxYgUAoHv37lZOQkRUOzyCISIiWbDAEBGRLFhgiIhIFiwwREQkC/7Ib6Ni\nYmKsHYGIqE5YYGxUq1atrB2BiKhO2EVmo/bt24d9+/ZZOwYRUa3xCMZGffHFFwCAXr16WTkJEVHt\n8AiGiIhkwQJDRESyUKyLLD09HQkJCRAEAWFhYYiIiDB7fN26dTh16hQA4M6dO7hx4wbWrVunVDwi\nIqpnihQYQRAQHx+P2bNnQ6fTYebMmQgODoavr6/UZsSIEdL/v/32W5w7d06JaEREJBNFCkx2dja8\nvb3h5eUFAAgJCUFaWppZgblfSkoKBg0apEQ0m/XJJ59YOwIRUZ0oUmAKCgqg0+mkaZ1Oh6ysrErb\nXr16Ffn5+ZXeD+Vx0rx5c2tHICKqE0UKjCiKFeapVKpK26akpKBLly5Qqys//yA5ORnJyckA7l7t\nrtfr6y9oLWi1WlkybN26FQAwcOBAaV593Sisvii57eXaznKxt7wAMyvFHjPXliIFRqfTwWg0StNG\noxEeHh6Vtk1NTcWoUaOqXFZ4eDjCw8OlaYPBUH9Ba0Gv18uSIS4uDgAQGhpa78uuL0pue7m2s1zs\nLS/AzEqxhcw+Pj6KrEeR05T9/f2Rl5eH/Px8mEwmpKamIjg4uEK73Nxc3L59G61bt1YiFhERyUiR\nIxiNRoPIyEhER0dDEASEhobCz88PSUlJ8Pf3l4rNkSNHEBISUmX3GRER2Q/FroMJCgpCUFCQ2bzB\ngwebTT/uZ44RET1KeCU/ERHJgoNd2qh7g10SEdkrFhgb5enpae0IRER1wi4yG5WUlISkpCRrxyAi\nqjUWGBu1detW6WJLIiJ7xAJDRESyYIEhIiJZsMAQEZEsWGCIiEgWPE3ZRm3YsMHaEYiI6oQFxka5\nuLhYOwIRUZ2wi8xGrVu3DuvWrbN2DCKiWmOBsVG7d+/G7t27rR2DiKjWWGCIiEgWLDBERCQLFhgi\nIpIFCwwREclCsdOU09PTkZCQAEEQEBYWhoiIiAptUlNTsXXrVqhUKjzxxBOYNGmSUvFszrZt26wd\ngYioThQpMIIgID4+HrNnz4ZOp8PMmTMRHBwMX19fqU1eXh527NiBDz/8EK6urrhx44YS0YiISCaK\ndJFlZ2fD29sbXl5e0Gq1CAkJQVpamlmb/fv348UXX4SrqysAwN3dXYloNmvVqlVYtWqVtWMQEdWa\nIkcwBQUF0Ol00rROp0NWVpZZm9zcXADAnDlzIAgCBg4ciA4dOlRYVnJyMpKTkwEAMTEx0Ov1Miav\nmVarlSXDoUOHAACzZ8+W5l2p97XUjZLbXq7tLBd7ywsws1LsMXNtKVJgRFGsME+lUplNC4KAvLw8\nzJs3DwUFBZg7dy6WLFmChg0bmrULDw9HeHi4NG0wGOQJbSG9Xi9LhrKyMgDWf33VUTKbXNtZLvaW\nF2BmpdhCZh8fH0XWo0gXmU6ng9FolKaNRiM8PDzM2nh6eqJTp07QarVo2rQpfHx8kJeXp0Q8IiKS\ngSIFxt/fH3l5ecjPz4fJZEJqaiqCg4PN2nTu3BkZGRkAgMLCQuTl5cHLy0uJeEREJANFusg0Gg0i\nIyMRHR0NQRAQGhoKPz8/JCUlwd/fH8HBwXjmmWfw66+/YvLkyVCr1Rg6dCgaNWqkRDyb5OzsbO0I\nRER1ohIr+4HEjtw7OcBalOxPLR/ziiLrsZRm9U7F1mUL/dYPw97yAsysFFvI/Ej9BkNERI8fFhgb\ntWzZMixbtszaMYiIao0FxkalpKQgJSXF2jGIiGqNBYaIiGTBAkNERLJggSEiIlkoNlw/PZwHRzog\nIrI3LDA2avXq1daOQERUJ+wiIyIiWbDA2KiFCxdi4cKF1o5BRFRr7CKzUceOHbN2BCKiOuERDBER\nycLiAnPz5k05cxAR0SPG4i6yd955B+3bt0f37t0RHBwMrZa9a0REVDWLj2A+/fRTtGvXDl999RXG\njBmDzz//HL///ruc2R5rzZo1Q7Nmzawdg4io1iw+DHFzc0OfPn3Qp08f5Obm4vDhw4iNjYVKpcJz\nzz2Hnj17okmTJnJmfazExsZaOwIRUZ3U6kf+69ev4/r16yguLoaXlxcKCgowbdo07Nixo77zERGR\nnbL4CObixYv44Ycf8MMPP8DZ2RnPP/88Pv74Y3h6egIA/uu//gtTp05FREREpc9PT09HQkICBEFA\nWFhYhXYHDx7Ehg0bpOX17t0bYWFhtX1ddm/u3LkAgP/+7/+2chIiotqxuMDMmzcP3bp1w5QpU9Cq\nVasKjzdt2hR9+vSp9LmCICA+Ph6zZ8+GTqfDzJkzERwcDF9fX7N2ISEhGDVq1EO+hEdTZmamtSMQ\nEdWJxQXmn//8JwIDAyvMz87OlgrO4MGDK31udnY2vL294eXlBeBuIUlLS6tQYIiI6NFhcYFZtGgR\nEhMTK8yPjo5GQkJCtc8tKCiATqeTpnU6HbKysiq0+/nnn3H69Gk0a9YMb775JvR6fYU2ycnJSE5O\nBgDExMRU2kZJWq1WlgwODg4AYLbsK/W+lrpRctvLtZ3lYm95AWZWij1mrq0aC4wgCAAAURSlf/dc\nuXIFGo2mxpXc/5x7VCqV2XTHjh3RrVs3ODg4YN++fYiLi8O8efMqPC88PBzh4eHStMFgqHH9ctLr\n9bJkKCsrA2D911cdJbPJtZ3lYm95AWZWii1k9vHxUWQ9NRaYN954Q/r/66+/bvaYWq3Ga6+9VuNK\ndDodjEajNG00Givc76RRo0bS/8PDw7Fp06Yal/soa9mypbUjEBHVSY0FZuXKlRBFEfPnz8cHH3wg\nzVepVHBzc4Ojo2ONK/H390deXh7y8/Ph6emJ1NRUTJw40azNtWvXpKJz9OjRx/73mcWLF1s7AhFR\nndRYYO5dPPnpp5/WeiUajQaRkZGIjo6GIAgIDQ2Fn58fkpKS4O/vj+DgYHz77bc4evQoNBoNXF1d\nMW7cuFqvj4iIrK/aAvP5559j7NixAO4eyVQlKiqqxhUFBQUhKCjIbN79Z50NGTIEQ4YMqXE5j4tp\n06YB4JEMEdmvagtM06ZNpf/fO8WYlHH27FlrRyAiqpNqC8z9P+APHDhQ9jBERPToqLbAZGRkWLSQ\ndu3a1UsYIiJ6dFRbYD777LMaF6BSqar9fYaIiB5P1RaYuLg4pXLQAyobloeIyJ7wtpQ2iqMoE5G9\nq7bATJ48GcuWLQNw95bJVbGkK42IiB4v1RaYe9fAAMCECRNkD0P/79725p0ticheVVtg2rRpI/2f\nvwkoKy8vz9oRiIjqxOLfYEwmE7Zv346UlBRp3LCQkBD079/fovHIiIjo8WJxgVm9ejVyc3MxcuRI\nNGnSBFevXsWOHTuwZs0ajhtGREQVWFxg0tLSEBsbi4YNGwIAfH19ERAQwN9miIioUhYXmMaNG6O0\ntFQqMABw586dCvd1ofrRsWNHa0cgIqoTi4eK6d69OxYsWIDevXtLNxDbu3cvunfvLnvIx9HMmTOt\nHYGIqE4eeqiYL7/80mw6OTkZERER9ZuKiIjsHoeKsVFjxowBcPfkCiIie8ShYmzUtWvXrB2BiKhO\nLC4wRUVF2Lp1KzIzM3Hz5k2Ioig9ZslQMenp6UhISIAgCAgLC6uyW+2nn37C0qVLsXDhQvj7+1sa\nj4iIbIza0oZr1qzBuXPnMGDAANy6dQuRkZHQ6/Xo27dvjc8VBAHx8fGYNWsWli1bhpSUFOTk5FRo\nV1xcjG+//RYBAQEP9yqIiMjmWFxgTp48iSlTpqBTp05Qq9Xo1KkTJk+ejB9++KHG52ZnZ8Pb2xte\nXl7QarUICQlBWlpahXZJSUl45ZVX4ODg8HCvgoiIbI7FXWSiKKJBgwYAAGdnZ9y+fRuNGzfG5cuX\na3xuQUEBdDqdNK3T6ZCVlWXW5ty5czAYDOjYsSN27dpV5bKSk5ORnJwMAIiJiYFer7f0JchCq9XK\nkqFXr14AYLbsK/W+lrpRctvLtZ3lYm95AWZWij1mri2LC8wTTzyBzMxM/O1vf0ObNm0QHx8PZ2dn\nNGvWrMbn3v97zT0qlUr6vyAISExMtGjImfDwcISHh0vTBoPBwlcgD71eL0uGt956C4D1X191lMwm\n13aWi73lBZhZKbaQ2cfHR5H1WNxFNnbsWDRp0gQAEBkZCQcHB9y+fRtRUVE1PvfehZn3GI1GsxEA\nSkpKcPHiRXzwwQcYP348srKysHjxYvz5558P81qIiMiGWHwE4+XlJf3fzc2t2huQPcjf3x95eXnI\nz8+Hp6cnUlNTMXHiROnxBg0aID4+XpqeP38+hg0b9lifRTZ06FAAwMaNG62chIiodh7qOpgDBw6Y\nDdffrVs3hIaGmnV3VUaj0SAyMhLR0dEQBAGhoaHw8/NDUlIS/P39ERwcXKcX8SgqKSmxdgQiojqx\nuMBs3LgRaWlp6Nu3r9SHuGvXLuTm5krftqsTFBSEoKAgs3mDBw+utO38+fMtjUVERDbK4gJz8OBB\nLFq0yOxssKCgIEyfPt2iAkNERI8Xi3/kd3FxgYuLS4V5905dJiIiul+1RzBXrvz/lRd9+vTBxx9/\njIiICHh6esJoNGLnzp0WXclPD+/+U7GJiOxRtQXm/jO97jl16pTZdEZGBnr37l2/qQhvv/22tSMQ\nEdVJtQUmKSlJqRxERPSIeejh+g0GAwoKCuDp6fnYDHdgDQMGDAAAbNu2zcpJiIhqx+ICc+3aNSxf\nvhxnzpxBo0aNcPPmTbRu3RqTJk2Cp6ennBmJiMgOWXwW2erVq/HEE08gISEBX3zxBRISEvDkk0/y\njotERFQpiwvMH3/8geHDh8PZ2RnA3RGVhw4dijNnzsgWjoiI7JfFXWQNGzZETk4OnnzySWlebm4u\nr4OhapWPeaVellNftyrQrN5ZT0sioppYXGBeeeUVfPjhh+jZsyeaNGmCq1ev4uDBg1UO90J1069f\nP2tHICKqE4sLTHh4OLy9vXHkyBFcuHABHh4emDRpEtq1aydnvsfWiBEjrB2BiKhOLCowgiDg008/\nxdixY1lQFFJcXAwAFYbnISKyFxYVGLVajZMnT9Y4LD/Vn2HDhgHgdTBEZL8sPousb9++2LJlC0wm\nk5x5iIjoEWHxbzB79uzB9evX8fXXX8PNzc3ssc8++6zegxERkX2zuMBMmDChTitKT09HQkICBEFA\nWFgYIiIizB7ft28f9u7dC7VaDWdnZ4wdOxa+vr51WicREVmPxQWmdevW2L59u9ktk0NCQtC/f/8a\nnysIAuLj4zF79mzodDrMnDkTwcHBZgXk2WefRa9evQAAR48eRWJiIt5///1avCQiIrIFFheY1atX\nIzc3FyNHjpSug9mxYwfWrFmDcePGVfvc7OxseHt7w8vLCwAQEhKCtLQ0swJz/wWbJSUlj/0JBQMH\nDrR2BCKiOrG4wKSlpSE2NhYNGzYEAPj6+iIgIMCirrOCggKzWy3rdDpkZWVVaLdnzx58/fXXMJlM\nmDt3rqXRHkm8gJWI7J3FBaZx48YoLS2VCgwA3LlzBx4eHjU+VxTFCvMqO0Lp3bs3evfujSNHjmD7\n9u2Iioqq0CY5ORnJyckAgJiYGKvfMkCr1cqSwWAwAIDZsutruJT6YsnrtsfM9UGu/UJOzKwMe8xc\nWxYXmO7du2PBggXo3bs3dDodjEYj9u7di+7duyMjI0NqV9mFmPfa32M0GqstTCEhIVWO0hweHm52\nO+F7H8TWotfrZclgD/eDsfa2rw2lMsu1X8iJmZVhC5l9fHwUWY/FBea7774DAHz55ZcV5t97TKVS\nYeXKlRWe6+/vj7y8POTn58PT0xOpqakVbsecl5eHZs2aAQCOHz8u/Z+IiOyTxQUmLi6u1ivRaDSI\njIxEdHQ0BEFAaGgo/Pz8kJSUBH9/fwQHB2PPnj347bffoNFo4OrqivHjx9d6fUREZH0Pfcvk2goK\nCkJQUJDZvPt/yB45cqRSUYiISAEWDxVDRET0MBQ7gqGHc2+wSyIie8UCY6NeffVVa0cgIqoTdpHZ\nqEuXLuHSpUvWjkFEVGs8grFRkyZNAmDb18EQEVWHRzBERCQLFhgiIpIFCwwREcmCBYaIiGTBH/lt\n1FtvvWXtCEREdcICY6Pu3d2TiMhesYvMRmVnZyM7O9vaMYiIau2xPYIpH/NKvSynvm6opVm902x6\nxowZAHgdDBHZLx7BEBGRLFhgiIhIFiwwREQkCxYYIiKShWI/8qenpyMhIQGCICAsLAwRERFmj+/e\nvRv79++HRqOBm5sb3nnnHTRp0kSpeDZn4sSJ1o5ARFQnihQYQRAQHx+P2bNnQ6fTYebMmQgODoav\nr6/U5sknn0RMTAycnJywb98+bNy4EZMnT1Yink3q3r27tSMQEdWJIl1k2dnZ8Pb2hpeXF7RaLUJC\nQpCWlmbWpl27dnBycgIABAQEoKCgQIloNisjIwMZGRnWjkFEVGuKHMEUFBRAp9NJ0zqdDllZWVW2\nP3DgADp06KBENJs1f/58ALwOhojslyIFRhTFCvNUKlWlbQ8fPoyzZ89KH7APSk5ORnJyMgAgJiYG\ner2+Vpnq6wLJ+vLg63BwcKgw39YzV8YeM9cHrVar2LrqCzMrwx4z15YiBUan08FoNErTRqMRHh4e\nFdqdPHkSX375JebPny99wD4oPDwc4eHh0rTBYKj/wFbw4OsoKyurdL4tseVsVVEqs16vt7vtw8zK\nsIXMPj4+iqxHkd9g/P39kZeXh/z8fJhMJqSmpiI4ONiszblz57B69WpMmzYN7u7uSsQiIiIZKXIE\no9FoEBkZiejoaAiCgNDQUPj5+SEpKQn+/v4IDg7Gxo0bUVJSgqVLlwK4W+WnT5+uRDwiIpKBYtfB\nBAUFISgoyGze4MGDpf/PmTNHqSh2gcWViOzdYzuasq3r1KmTtSMQEdUJh4qxUWlpaRWuFSIisic8\ngrFRixYtAsDrYIjIfvEIhoiIZMECQ0REsmCBISIiWbDAEBGRLPgjv42qaiw2IiJ7wQJjo9q1a2ft\nCEREdcIuMht1+PBhHD582NoxiIhqjUcwNmrFihUAeGdLIrJfPIIhIiJZsMAQEZEsWGCIiEgWLDBE\nRCQL/shvo2JiYqwdgYioTlhgbFSrVq2sHYGIqE4UKzDp6elISEiAIAgICwtDRESE2eOZmZlITEzE\nX3/9hXfffRddunRRKppN2rdvHwCgV69eVk5CRFQ7ivwGIwgC4uPjMWvWLCxbtgwpKSnIyckxa6PX\n6zFu3Dg8++yzSkSyeV988QW++OILa8cgIqo1RY5gsrOz4e3tDS8vLwBASEgI0tLS4OvrK7Vp2rQp\nAEClUikRiYiIZKbIEUxBQQF0Op00rdPpUFBQoMSqiYjIShQ5ghFFscK82h6pJCcnIzk5GcDdM630\nen2tlnOlVs+Sz4Ovw8HBocJ8W89cGXvMXB+0Wq1i66ovzKwMe8xcW4oUGJ1OB6PRKE0bjUZ4eHjU\nalnh4eEIDw+Xpg0GQ53z2YIHX0dZWVml822JLWerilKZ9Xq93W0fZlaGLWT28fFRZD2KFBh/f3/k\n5eUhPz8fnp6eSE1NxcSJE5VYtd365JNPrB2BiKhOFCkwGo0GkZGRiI6OhiAICA0NhZ+fH5KSkuDv\n74/g4GBkZ2fj448/xu3bt3Hs2DFs2bIFS5cuVSKeTWrevLm1IxAR1Yli18EEBQUhKCjIbN7gwYOl\n/7dq1QqrVq1SKo7N++qrrwAAr776qpWTEBHVDq/kt1EbNmwAwAJjDeVjXqnzMurr5AbN6p31tCQi\n5XGwSyIikgULDBERyYIFhoiIZMECQ0REsuCP/DaKA10Skb1jgbFRnp6e1o5ARFQn7CKzUUlJSUhK\nSrJ2DCKiWmOBsVFbt27F1q1brR2DiKjWWGCIiEgWLDBERCQLFhgiIpIFCwwREcmCpynbqHuDXRIR\n2SsWGBvl4uJi7QhERHXCLjIbtW7dOqxbt87aMYiIao0Fxkbt3r0bu3fvtnYMIqJaU6yLLD09HQkJ\nCRAEAWFhYYiIiDB7vKysDCtXrsTZs2fRqFEjvPvuu2jatKlS8YiIqJ4pUmAEQUB8fDxmz54NnU6H\nmTNnIjg4GL6+vlKbAwcOoGHDhoiNjUVKSgo2bdqEyZMnKxGPyK7Vxx04Ad6Fk+qfIl1k2dnZ8Pb2\nhpeXF7RaLUJCQpCWlmbW5ujRo+jRowcAoEuXLsjIyIAoikrEIyIiGShSYAoKCqDT6aRpnU6HgoKC\nKttoNBo0aNAAN2/eVCIeERHJQJEussqORFQq1UO3AYDk5GQkJycDAGJiYuDj41O7UF8frd3zFJKa\nmlpxpo1nrhQzy8/e8tazWn8GWJE9Zq4NRY5gdDodjEajNG00GuHh4VFlm/LychQVFcHV1bXCssLD\nwxETE4OYmBh5Q1toxowZ1o7w0JhZfvaWF2Bmpdhj5tpSpMD4+/sjLy8P+fn5MJlMSE1NRXBwsFmb\njh074uDBgwCAn376CW3btq30CIaIiOyDIl1kGo0GkZGRiI6OhiAICA0NhZ+fH5KSkuDv74/g4GD0\n7NkTK1euxIQJE+Dq6op3331XiWhERCQTxa6DCQoKQlBQkNm8wYMHS/93dHTEe++9p1ScehMeHm7t\nCA+NmeVnb3kBZlaKPWauLZXIc4GJiEgGHCqGiIhkwdGU66Cm4W9szaefforjx4/D3d0dS5YssXac\nGhkMBsTFxeH69etQqVQIDw9Hnz59rB2rWnfu3MG8efNgMplQXl6OLl26YNCgQdaOZRFBEDBjxgx4\nenraxZlO48ePh7OzM9RqNTQajc2cWVqV27dvY9WqVbh48SJUKhXeeecdtG7d2tqxZMUCU0uWDH9j\na3r06IHevXsjLi7O2lEsotFoMGzYMLRs2RLFxcWYMWMG2rdvb9Pb2MHBAfPmzYOzszNMJhPmzp2L\nDh062MUHyTfffIPmzZujuLjY2lEsNm/ePLi5uVk7hkUSEhLQoUMHTJkyBSaTCaWlpdaOJDt2kdWS\nJcPf2JrAwMBKry2yVR4eHmjZsiWAu/fHad68eYURIGyNSqWCs7MzgLvXc5WXl9vF6fZGoxHHjx9H\nWFiYtaM8koqKinD69Gn07NkTAKDVatGwYUMrp5Ifj2BqqbLhb7KysqyY6NGWn5+Pc+fOoVWrVtaO\nUiNBEDB9+nRcvnwZL774IgICAqwdqUbr1q3D0KFD7eroBQCio6MBAC+88IJNn52Vn58PNzc3fPrp\np/jrr7/QsmVLjBgxQvoy8qjiEUwtWTq0DdVdSUkJlixZghEjRqBBgwbWjlMjtVqNjz76CKtWrcKf\nf/6JCxcuWDtStY4dOwZ3d3fpaNFefPjhh1i0aBFmzZqFvXv3IjMz09qRqlReXo5z586hV69eWLx4\nMZycnLBjxw5rx5IdC0wtWTL8DdWdyWTCkiVL8Nxzz+E///M/rR3noTRs2BCBgYFIT0+3dpRq/fHH\nHzh69CjGjx+P5cuXIyMjAytWrLB2rBp5enoCANzd3dGpUydkZ2dbOVHVdDoddDqddDTbpUsXnDt3\nzsqp5McCU0uWDH9DdSOKIlatWoXmzZujX79+1o5jkcLCQty+fRvA3TPKfvvtNzRv3tzKqao3ZMgQ\nrFq1CnFxcXj33XfRrl07TJw40dqxqlVSUiJ155WUlODkyZNo0aKFlVNVrXHjxtDpdMjNzQUA/Pbb\nbzZ9skp94W8wtVTV8De2bPny5cjMzMTNmzfx9ttvY9CgQdKPjrbojz/+wOHDh9GiRQtMnToVAPDG\nG29UGBHClly7dg1xcXEQBAGiKKJr167o2LGjtWM9cm7cuIGPP/4YwN3up2effRYdOnSwcqrqRUZG\nYsWKFTBs0TDPAAAUrElEQVSZTGjatCnGjRtn7Uiy45X8REQkC3aRERGRLFhgiIhIFiwwREQkCxYY\nIiKSBQsMERHJwi4LzPz587F//35rx1DMli1brHrh2+bNmzFq1CiMGTPGahnuN2zYMFy5csXaMSoY\nNGgQLl++XKdl5OfnY9CgQSgvL6+nVPXjvffew6lTp6wdQ1Y5OTl2MYr0w5DrsyMxMRH79u2rsZ3N\nFZiDBw9izpw5Vs1gq29yazAYDNi1axeWLl2K1atXV3i8qr/X+PHjcfLkSQBAXFwcNm/eXG+ZNmzY\nAC8vr3pbHtVs6dKlaNu2bb0us773i/vduHEDy5cvx9ixY/Hmm29izpw5NY4VuHnzZrz88suy5FHC\nqVOn8Pbbb9fb8s6fP4/p06dj6NChmD59Os6fPy899sorr+DLL7+EyWSqdhk2V2BIXg9bNA0GAxo1\nagR3d3eZEtkfa33xEEURgiBYZd32pqSkBK1atUJMTAwSEhLw/PPPIyYmBiUlJZW2v3btGk6dOoVO\nnTrVan2P2t/FZDJh8eLFeO6556Ttt3jxYqmgeHh4wMfHB0ePHq12OVa5kn/Hjh3Yv38/bty4AZ1O\nhzfeeAOdO3dGTk4OVq9eDZPJhGHDhkGj0WDdunU1Lu/AgQPYtWsXrl+/jlatWuGtt95CkyZNANy9\nB8Mvv/yCoqIieHt7Y8SIEXj66acB3B1yf82aNcjLy4OjoyOeffZZvPnmm5g3bx4AYMSIEQCAOXPm\nVLifx5YtW5CTkwNHR0f88ssv0Ov1GD9+PPz9/QHc7S5ZsWIFvL29Adz9tqbT6fD666/j1KlTiI2N\nxUsvvYRdu3ZBrVZj9OjR0Gq1SExMRGFhIV5++WX0799fWl9ZWRmWLVuGEydOoFmzZnjnnXfw5JNP\nArg7svPatWtx+vRpODs7o2/fvtKNubZs2YKLFy/CwcEBx44dw/DhwysMyV5UVIS1a9fixIkTcHJy\nQlhYGF577TVkZGRg0aJF0t+jS5cuGD9+/EP8pYHk5GQcOXIEAPD111+jbdu2mDFjBnJycrBmzRqc\nP38enp6eGDJkiDTUTlxcHBwcHHDlyhVkZWXhP/7jPxAVFSX9Te/ftsePH8fmzZtx5coVNGjQAKGh\nodXe4CstLQ1btmyRRrcdNWoUOnTogO+//x47d+6E0WiEm5sbXn31VbzwwgsAIP29evfuja+//hrt\n27fHhAkTsHPnTuzevRsqlQqDBw82W8/x48exYcMGGI1GuLi4oG/fvnjllVcq5BEEARs3bsShQ4fg\n4uJSYUic+fPn46mnnkJmZibOnj2LJUuW4PTp01VmnTdvHl566SV06dIFv//+O+bOnYsZM2YgKCgI\nJ0+exIYNG/DRRx/h4MGD2L9/PwICAvD999+jQYMGGD16NP7+979Xut3Gjx+PsWPHon379jXu+/cT\nRRGJiYk4cuQIysrK0KRJE0ycOBFnzpypdL+wZF9Wq9WVvg/u5+XlZbYtw8PDsWHDBuTm5lY6qOfJ\nkyfRsmVLODo6SvNq2kcdHR1hMBiQmZmJqVOnon379mbLLCoqQmJiIk6cOAGVSiXtm2q1GqtXr0Zh\nYSGmTJkCANi4cSPOnj2LOXPm4Pbt21i5ciWysrIgCAKeeuopjBkzRhq9/datW1i/fj1+/fVX3Llz\nB08//TQmTpyIBQsWSO9VAPjkk0/M8ixcuBAdOnTASy+9JM375z//iUGDBqFz585mbU+dOoXy8nL0\n7dsXKpUKffr0wa5du5CRkSGNmNC2bVscP34cXbp0qbA9JaIVpKamikajUSwvLxdTUlLEoUOHigUF\nBaIoiuL3338vzp49u9rnz5s3T0xOThZFURR//vlnMSoqSrx48aJoMpnEbdu2ie+//77U9tChQ2Jh\nYaFoMpnEnTt3iqNHjxZLS0tFURTFWbNmiYcOHRJFURSLi4vFP/74QxRFUbxy5Yo4cOBA0WQyVZkh\nKSlJHDJkiHjs2DGxvLxc3LRpkzhr1izp8YEDB4p5eXnS9MqVK8V///vfoiiKYkZGhjh48GBx69at\nYllZmfjdd9+JkZGR4vLly8WioiLxwoUL4pAhQ8TLly9L63r99dfFH3/8USwrKxO/+uorcdy4cWJZ\nWZlYXl4uTps2TVrW5cuXxfHjx4snTpwwe+7PP/8slpeXS6/9frGxseKiRYvEoqIi8cqVK+LEiRPF\n/fv3S1nHjh1b5Xao6u81btw48ddff63w2kVRFMvKysSoqChx+/btYllZmfjbb7+Jw4YNEy9duiS1\nHzZsmHjq1Cnxzp074tq1a83Wcf+2zcjIEP/66y+xvLxcPH/+vDh69Gjx559/rjRrVlaWOHz4cPHX\nX38Vy8vLRaPRKObk5IiiKIrHjh0T8/LyREEQxFOnTon/+Mc/xD///NPs77Vhwwbxzp07YmlpqXji\nxAlx9OjR4l9//SUWFxeLy5cvN8s1ZswYMTMzUxRFUbx586a0rAft3btXnDRpknj16lXx5s2b4vz5\n8832vXnz5olvv/22eOHCBdFkMollZWXVZt28ebMYHx8viqIobt++XYyKihI3bNggPbZ27Vrp7/b6\n66+L3333nVheXi7u3btXfOutt0RBECrNef/fs6Z9/34nTpwQp02bJt66dUsUBEG8ePGi9F5/cL+w\ndF+u7H1Qk3PnzolDhgwRb9++Xenj69evF1evXi1NW7KPDh8+XDx9+nSV76tFixaJn3/+uVhcXCxe\nv35dnDFjhrhv3z5RFEWxpKREnDhxovj999+LmZmZYmRkpGgwGERRFMXCwkLxxx9/FEtKSsSioiJx\nyZIl4qJFi6TlLliwQFy6dKl48+ZNsaysTDx16pQoipW/V5OSksRPPvlEFEVRTElJEWfOnGm2TUaO\nHFnp9tu1a5cYHR1tNm/hwoXizp07pemffvpJnDZtWqXb8x6rdJF17doVnp6eUKvVCAkJgbe3d61H\nQk1OTsZrr70GX19faDQavPbaazh//jyuXr0KAOjevTsaNWoEjUaDl19+GSaTSRpwTqvV4vLlyygs\nLISzs/ND33WwTZs2CAoKglqtRvfu3c36KGui0WjQv39/aLVadOvWDTdv3kSfPn3g4uICPz8/+Pr6\n4q+//pLat2zZEl26dIFWq0W/fv1QVlaGrKws/PnnnygsLMSAAQOg1Wrh5eWFsLAwpKamSs9t3bo1\nOnfuDLVabfYNDbj7DTo1NRVDhgyBi4sLmjZtin79+uHw4cMPtS0eRlZWFkpKShAREQGtVot27doh\nKChI+kYLAEFBQQgMDISDgwPeeOMNnDlzBgaDocKy2rZtixYtWkCtVuOJJ55At27dqhy2/cCBAwgN\nDUX79u2hVqvh6ekpDUQZFBQEb29vqFQqBAYGon379vj999+l56pUKgwaNAgODg5wdHREamoqevTo\ngRYtWsDZ2RkDBw40W5dGo0FOTg6Kiorg6upa5VD4P/74I/r06QO9Xg9XV9dKb7vdo0cP+Pn5QaPR\nQKvVVps1MDBQev2nT59GREQETp8+DQDIzMxEYGCgtFy9Xo/w8HCo1Wo8//zzuHbtGm7cuFFpzgdZ\nuu9rtVqUlJTg0qVLEEURvr6+VY46bsm+XNX7oDpFRUWIjY3FgAEDqrzdw+3bt+Hi4iJNW7KPdurU\nCW3atKn0fXX9+nWkp6dL93xxd3dH3759pdfi5OSECRMmIDExEbGxsRg5cqR0hNKoUSN06dIFTk5O\ncHFxQf/+/aW/4bVr15Ceno4xY8bA1dUVWq3W7G9anU6dOiEvLw95eXkAgMOHDyMkJARabcWOrJKS\nkgrbqkGDBmb3C3JxcZEGdq2KVbrIDh06hN27d0tFoKSkBDdv3qzVsq5evYqEhASsX79emieKIgoK\nCtCkSRPs2rULBw4cQEFBAVQqFYqLi6V1vf3220hKSsLkyZPRtGlTDBgw4KEGJrz/dwlHR0eUlZWh\nvLwcGo2mxuc2atQIarVaem5ly7u/v/j+m5up1WrodDpcu3YNwN2d7l53HnC3aNzrBnzwuQ8qLCyE\nyWSCXq+X5jVp0sTiO0eq1epKf5MwmUyV7rj38ur1eun1V7bO+zM7OzvD1dVVet79srKy8K9//QsX\nLlyAyWSCyWSq8pDdaDRW2QV04sQJbNu2Dbm5uRBFEaWlpWaj87q5uZl9iFy7ds2saNzrvrtnypQp\n+N///V/861//QosWLfCPf/yj0i8wD76mB5fz4LaoKWvr1q2Rl5eH69evSz/SbtmyBYWFhcjOzjbb\nLxo3biz938nJCQCq/I3iQZbu++3atcOLL76I+Ph4GAwGdO7cGcOGDav0g/7q1asPtS8/+D6ozJ07\nd7Bo0SIEBATgtddeq7Kdq6ur2Yfnw+6jDzIYDCgvL8dbb70lzRNF0ew5rVq1gpeXF27cuIGQkBBp\nfmlpKRITE5Geni59gBcXF0MQBBiNRri6utbqzrQODg7o2rUrfvjhBwwYMAApKSlSF92DnJ2dK9x8\nrqioyKwIFxcX13hXTsULzNWrV/H5559j7ty5aN26NdRqNaZOnVrpDbwsodfr0b9/fzz33HMVHjt9\n+jS++uorzJ07F76+vlCr1Rg5cqS0rmbNmuHdd9+FIAj45ZdfsHTpUsTHx9fLjcOcnJzM7rl9/fr1\nanfImtx/75l7O5qHhwc0Gg2aNm1a61MR3dzcoNFoYDAYpOHDDQaDdK+Nmuj1ehgMBoiiKG230tJS\nFBYWSh+cD25PDw8PGAwGCIIgvYENBgOaNWtW6estKSnBrVu3Kv3mu2LFCrz44ouYOXMmHB0dsW7d\nOhQWFlaaVafTVXoacVlZGZYsWYKoqCgEBwdDq9Vi8eLFZm0qew33Z3zw6KpVq1aYNm0aTCYT9uzZ\ng2XLluGzzz6rsO5726Kq5Ty47pqyOjk5oWXLlvjmm2/QokULaLVaPPXUU9i9eze8vb2tcv/6Pn36\noE+fPrhx4waWLVuGnTt34vXXX6+wTfV6fY37clXvg8qUlZXho48+gqenp9kHfWVatGiBQ4cOSdOW\n7KPVfU7odDpotVrEx8dX+YVzz549KCsrg6enJ7766iupAO7atQu5ublYsGABGjdujPPnz2PatGlS\ngbp16xZu375d4cPdks+tHj16IDY2Fm3atIGTk1OVvTZ+fn7YvXu32fv6woUL6N27t9Tm0qVLeOKJ\nJ6pdn+JdZKWlpVCpVNKO/v333+PixYvS440bN0ZBQUGNp7/d88ILL2DHjh3SMoqKivDjjz8CuFth\nNRoN3NzcIAgCtm3bhqKiIum5hw8fRmFhIdRqtfSNSq1Ww83NDSqVqk7XWjz55JM4cuQIBEFAenp6\nne+2d/bsWfz8888oLy/HN998AwcHBwQEBKBVq1ZwcXHBjh07cOfOHQiCgAsXLljc5ahWq9G1a1f8\n+9//RnFxMa5evYrdu3dXWrArExAQAAcHB2n9JSUl2LRpE1q2bCl9G3d3dzfblgEBAXB2dsbOnTth\nMplw6tQpHDt2DN26dZPanDhxAr///jtMJhM2b96MgICACkcvwN2/saurKxwdHZGdnW3WhfGgnj17\n4uDBg/jtt98gCAIKCgpw6dIlmEwmlJWVScX2xIkT0inWVenatSsOHjyInJwclJaWYuvWrdJjJpMJ\nP/zwA4qKiqDVatGgQQOzb8IPLufbb7+F0WjErVu3arzLoSVZn376aezdu1f65h8YGGg2raTs7Gxk\nZWXBZDLByckJDg4O0rZ4cL+wZF+u6n3woHs3qnNwcEBUVFSV2/+e9u3b49y5c7hz5w4Ay/bR6nh4\neOCZZ57B+vXrUVRUBEEQcPnyZelzIDc3F0lJSZgwYQKioqKwc+dOqZuxpKQEjo6OaNCgAW7dumW2\nb3l4eKBDhw5Ys2YNbt26BZPJJC3T3d0dN2/eNPuMe9C9L/Xr169H9+7dq2zXtm1bqNVqfPvttygr\nK8OePXsA3D0ivSczM7PGWyQofgTj6+uLfv364f3335f6b5966inp8Xbt2sHX1xdjxoyBWq1GfHx8\ntcvr3LkzSkpKsHz5chgMBjRo0AB/+9vf0LVrV3To0AEdOnTApEmT4OTkhL59+5p9SKWnp2P9+vUo\nLS1FkyZNMGnSJKkbpH///pgzZw7Ky8sxa9ash/59ZsSIEYiLi8PevXvRqVOnWp/+eE9wcDBSU1MR\nFxcHb29vTJkyReqCmj59OtavX4/x48fDZDLBx8enwllN1YmMjMTatWsRFRUFR0dHhIWFITQ01KLn\nOjg4YObMmUhMTMTu3buhVqvRpk0bTJ48Wfrm07NnTyxduhQjRoxAYGAgpk2bhmnTpmHNmjX48ssv\n4enpiaioKLMbc3Xr1g1bt27FmTNn0LJlyypvgDV69GisX78ea9euRWBgILp27Vplv3CrVq0wbtw4\nJCYmIj8/H+7u7hg1ahSaN2+OkSNHYtmyZSgrK0PHjh1rvHnc3//+d/Tt2xcffPAB1Go1Bg8ebFbc\nDh8+jLVr10IQBPj4+GDChAmVLicsLAy5ubmYOnUqXFxc8PLLLyMjI6PK9bq4uNSYNTAwEDt27JD6\n5gMDA1FcXGxxX319Ki4uRmJiIq5cuQJHR0c888wz0tl0le0XNe3L1b0P7nfmzBkcP34cjo6OZl1u\ns2bNqrTQNm7cGO3atcPRo0el3yVq2kdrEhUVhU2bNuG9995DcXExvLy88Oqrr6K8vByxsbF49dVX\npTPg3njjDcTGxiImJgZ9+vTBihUrMGrUKHh6eqJfv35IS0uTljthwgSsW7cOkydPhslkQtu2bREY\nGIjmzZujW7duiIqKgiAIWLp0aaW5unfvjqSkJOkeS5XRarWYOnUqVq1ahU2bNsHX1xdTp06VtvW1\na9eQk5NT4eyzB/F+MGRz7j+lm+ieLVu24PLly7LdbTMnJwdxcXFYsGBBvXST26pDhw4hOTkZH374\nYa2XsX79enh5eeHFF1+sth3vaElEhLu9KwsXLrR2DFmVlpZi37596NWrV52WM3z4cIva8Up+IqLH\nQHp6OkaPHg13d3c8++yziqyTXWRERCQLHsEQEZEsWGCIiEgWLDBERCQLFhgiIpIFCwwREcmCBYaI\niGTxf7dVpWKb++5AAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10cf331d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(k_values, f)\n",
    "plt.xlabel(\"at least number of Utopia cards drawn in step 2 (or exactly 0)\")\n",
    "plt.ylabel(\"probability\")\n",
    "plt.vlines(0.5, ymin=0, ymax=0.8, linestyles=\"--\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot shows that with 76% probability at least 1 cube will get removed and with 34% at least 2 cubes will be removed. To say it in another way: in about 1/3 of all games, Utopia will get hit _at least twice_ during setup.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}