monty_{hall,carlo}

PCMLAI_M3_Monty.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "3cd1829e",
      "metadata": {
        "id": "3cd1829e"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import random\n",
        "from random import randrange\n",
        "from matplotlib import pyplot as plt"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "349dffeb",
      "metadata": {
        "id": "349dffeb"
      },
      "source": [
        "## Monty Hall Problem and Monte Carlo Simulation\n",
        "### Let's play some Monty Hall problem with a host!"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "739c01a8",
      "metadata": {
        "id": "739c01a8"
      },
      "source": [
        "#### Step 1: Construct three 'imaginary' doors, one with a car behind. We don't see which one!"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "27f976e5",
      "metadata": {
        "id": "27f976e5",
        "outputId": "7d59e2fb-a8df-48f4-86a3-2b0c50bbd084"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "You picked door 1.\n",
            "The remaining doors are therefore: [0, 2]\n"
          ]
        }
      ],
      "source": [
        "random.seed(23) # Set the random number generator \n",
        "doors = [0,1,2]\n",
        "winning_door = randrange(3) # it can be 0, 1, or 2! We don't know the truth yet.\n",
        "picked_door = 1 # user input, we pick door #1\n",
        "print(\"You picked door\", str(picked_door) + \".\")\n",
        "doors.remove(picked_door)\n",
        "print(\"The remaining doors are therefore:\", doors)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "c7b12ba1",
      "metadata": {
        "id": "c7b12ba1"
      },
      "outputs": [],
      "source": [
        "def reveal_empty_door(doors, winning_door):\n",
        "    \"\"\" \n",
        "    The host knows which door is winning, so opens one of the empty doors!\n",
        "    \"\"\"\n",
        "    sample_from = [] #Host cannot open the door that is the winning door. So the host has 1 or 2 candidate doors.\n",
        "    if len(doors) != 2:\n",
        "        print(\"There is an issue!\")\n",
        "        return None\n",
        "    else:\n",
        "        if doors[0]!= winning_door:\n",
        "            sample_from.append(doors[0])\n",
        "        if doors[1]!= winning_door:\n",
        "            sample_from.append(doors[1])\n",
        "    return random.choice(sample_from)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "82ae721a",
      "metadata": {
        "id": "82ae721a",
        "outputId": "2c0ee029-19f8-4660-fc3e-b5659224b32a"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "The host reveals door 0\n"
          ]
        }
      ],
      "source": [
        "open_door = reveal_empty_door(doors, winning_door)\n",
        "print(\"The host reveals door\", open_door)\n",
        "doors.remove(open_door)\n",
        "remaining_door = doors[0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "75df1912",
      "metadata": {
        "id": "75df1912",
        "outputId": "881d8c99-e3f1-498c-eb8f-4029180ecaa0"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "You initially picked door 1 but you can switch to door 2\n"
          ]
        }
      ],
      "source": [
        "print(\"You initially picked door\", picked_door, \"but you can switch to door\", remaining_door)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "cbc659dd",
      "metadata": {
        "id": "cbc659dd",
        "outputId": "606a380e-6295-4489-af1f-91eb03ca2bb3"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Type 'yes' if you would like to change your door.\n",
            "You keep your door 1\n"
          ]
        }
      ],
      "source": [
        "yes_no = input(\"Type 'yes' if you would like to change your door.\")\n",
        "if yes_no == \"yes\":\n",
        "    print(\"You switched to door\", remaining_door)\n",
        "    picked_door = remaining_door\n",
        "else:\n",
        "    print(\"You keep your door\", picked_door)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "38911454",
      "metadata": {
        "id": "38911454",
        "outputId": "16b7d473-9c2c-41a5-8f3d-bcec1feaa113"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Let's see if we won: (True means you are the winner!)\n",
            "True\n"
          ]
        }
      ],
      "source": [
        "print(\"Let's see if we won: (True means you are the winner!)\")\n",
        "print(picked_door == winning_door)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "bda15511",
      "metadata": {
        "id": "bda15511"
      },
      "source": [
        "#### Step 2: Now let's write a function which returns $1$ if we win, otherwise $0$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "a44ddb0d",
      "metadata": {
        "id": "a44ddb0d"
      },
      "outputs": [],
      "source": [
        "def simulate_round(picked_door, yes_no):\n",
        "    \"\"\"\n",
        "    Inputs: picked_door (0 1 or 2), decision of switching the door (\"yes\" or \"no\")\n",
        "    Returns: whether or not we win the prize (1: win, 0: lose)\n",
        "    \"\"\"\n",
        "    doors = [0,1,2]\n",
        "    winning_door = randrange(3) # it can be 0, 1, or 2! We don't know the truth yet.\n",
        "    doors.remove(picked_door)\n",
        "    open_door = reveal_empty_door(doors, winning_door)\n",
        "    doors.remove(open_door)\n",
        "    remaining_door = doors[0]\n",
        "    if yes_no == \"yes\":\n",
        "        picked_door = remaining_door\n",
        "    return int(picked_door == winning_door)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "88bb1f85",
      "metadata": {
        "id": "88bb1f85"
      },
      "source": [
        "#### Step 3: Let us simulate our decisions!"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "686ed1a9",
      "metadata": {
        "id": "686ed1a9"
      },
      "outputs": [],
      "source": [
        "picked_door = 1 #fix a door (0,1, or 2)\n",
        "yes_no = \"no\" #let's start with the strategy of keeping our door!\n",
        "simulate = 30000 #number of times to simulate\n",
        "results = np.zeros(simulate) #a lot of zeros\n",
        "for i in range(simulate):\n",
        "    results[i] = simulate_round(picked_door, yes_no)\n",
        "ratios = np.cumsum(results) / (np.arange(1,simulate+ 1)) #running ratios. ratios[10] would give us the percentage of times we won by round 10"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "aeba4efe",
      "metadata": {
        "id": "aeba4efe",
        "outputId": "2ee95388-0dc8-4fae-9041-db83bd3b4faa"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAA78klEQVR4nO3dd3hUVfrA8e+bkF4ooRdpUqSJEKoIWHAR3bWuqyL2gnVdy6q7v7Xturq7umtfLMvaxa64YpcqIE2qtFA0AQUCJCQhIe39/XFvwmQySe5AJsmE9/M8eXLn3nPPnDOTzDun3HNFVTHGGGO8iqjvAhhjjAkvFjiMMcYExQKHMcaYoFjgMMYYExQLHMYYY4JigcMYY0xQLHAYROQEEVlfj89/lIjkikhkfZXhUInIGyJylse074nI+GqOV3gfRKSXiHwnIjkicnMtFNeYWmGBoxESkbtFZIbfvo1V7LtAVeeqaq86LN9WETml7LGq/qiqiapaUldlqA0iMgA4FvjQ4ykPAw9WdTDA+/B7YJaqJqnqE4de0oZFRGaJyNhayOdFEbnssAtkgmaBo3GaAxxf9g1eRNoCUcAgv31Hu2lrjYg0qc38GoJq6nQt8Jp6vIpWVRcBySKS6vGpOwNrPKY1ps5Y4GicFuMEioHu49HATGC9375NqrpdRMaKSEbZyW6L4HYRWSki2SLypojEBnoiEblMRL4RkX+JyB7gPhHpLiJfi8huEckUkddEpJmb/hXgKOAjt3vq9yLSRUS07ANaRNqLyHQR2SMiaSJydVUVFZGmIvKyiOwSkR9E5P9EJEJEYkQkS0T6+aRtJSL5ItLafXyGiCx30813WxC+r8GdIrISyKsieJwGzPZ7LeaJyCMisldEtojIaX7nzAJOr6Iu5e+DiHwNnAg85b5OPf3Snigiq3wefykii3wezyvrQnNfz3fd12hLdd1eIjJBRL53u8e2icjt7v7mIvI/N4+97nZH99ivRWSpXz63icgHVT2PT7r7ROQt9z3MEZE1voFVRI5xWyhZ7rFf1ZSnqQOqaj+N8AcnUPzO3X4KuAKnm8R331R3eyyQ4XPuVmAR0B5oAawFJlfxPJcBxcBNQBMgDqclMw6IAVrhtGoe88v/FJ/HXQAFmriPZwPPALE4gW4XcHIVz/8yTldRkpvPBuBK99hU4EGftDcAn7rbg4CdwDAgErjULVeMTxmXA52AuADPm+CWuZXfa1EEXO3meR2wHRCfNLcC71VRF//3YRZwVRVpY4F8oKX7uv/sPleS+x7kAyk4Xw6XAvcA0UA3YDPwiyry/Qk4wd1uDgxyt1OAc4F49zneBj5wj8UAe4BjfPL5DjjXw9/pfUABMMF9zR4CFrrHooA04A9u2U8CcoBe9f3/daT/WIuj8ZqN06oAOAGY6/747psd4LwyT6jqdlXdA3zEwZZKINtV9UlVLVbVfFVNU9UvVPWAqu4C/gmM8VJoEekEjALuVNUCVV0OvABMCpA2EvgNcLeq5qjqVuBRn7SvAxf6nHKRuw+cD/dnVfVbVS1R1ZeAA8Bwv9cgXVXzAxS1mfs7x2//D6r6vDrjNS8B7YA2PsdzfM49ZKpaACzBeT9TgZXAPOB4tw4bVXU3MAQnuD2gqoWquhl4HrigiqyLgD4ikqyqe1V1mft8u1X1XVXdr6o5OF9CxrjHDgBvAhcDiEhfnCD+P4/VmaeqM9zX7BWccSPceiQCD7tl/9rN88Iq8jF1xAJH4zUHGCUizXE+ODYC84GR7r5+VD++8bPP9n6cf+CqpPs+EJHWIjLN7erYB7yK883Yi/bAHvfDqcwPQIcAaVvifBP9oYq0XwNxIjJMRDrjBL/33WOdgdvcLpAsEcnCaV20r6pefrLc30l++8tfN1Xd7276vnZJPucertk4rZTR7vYsnA/zMRz8UtAZaO9Xzz9QMZj5Ohfn2/8PIjJbREYAiEi8iDzrdgfuw/nbaSYHZ8K9BFwkIoITuN9yA4oX/n9rsW7XYHsgXVVLfY5X9bdg6pAFjsZrAdAUuAb4BkBV9+F0Z1yD00rYUkvP5T84/JC7b4CqJuN8E5Vq0vvaDrQQEd8P5KOAbQHSZuJ8Q+4cKK37gfMWzjfUi4D/+QSkdJxurGY+P/Gq+oaXcqpqHrAJ6FlVmiocA6wI8pyq+AeO2VQOHOnAFr96JqnqhEAZqupiVT0TaA18gPP6AdwG9AKGue9pWctV3PMWAoU4LdmLcFoOh2s70ElEfD+nqvpbMHXIAkcj5XavLMHpU5/rc2ieu69WZ1P5SQJygSwR6QDc4Xd8B05feyWqmo7TMnpIRGLdAesrgdcCpC3B+WB7UESS3FbFrTgtnDKv43RnTeRgNxU43TWT3daIiEiCiJzuF7BqMgOPXXA+xgCfBHlOVebjfJgPBRap6hqcIDqMg+/vImCfO9AfJyKRItJPRIb4ZyYi0SIyUUSaqmoRsA8omyKdhDNukiUiLYB7A5TnZZyxs2JVnVcL9fsWyAN+LyJR4kzh/SUwrRbyNofBAkfjNhvnm6PvP/Fcd18oA8f9OIPP2cDHwHt+xx8C/s/tOrk9wPkX4vSRb8fpWrpXVb+o4rluwvlw2YxTz9dxBsUBUNWyD5/2+Hxgq+oSnHGOp4C9OIOwlwVRR4DngIlu90yN3A/rPHWm5R42t9WzDFijqoXu7gU44yw73TQlOB+2A4EtOK20F3Bao4FMAra63VGTccctgMdwBt0zgYXApwHOfQWnC7Q2Whu4dfoVzuy1TJwJE5eo6rrayN8cOlG1GzkZc6hE5HWc/vwPPKR9F/iPqs6oKW04EpE4nJlqg9wxNdNIWeAwxtQKEbkVOENVT6rvspjQCmlXlYiMF5H14lzEdVcVacaKcxHWGhGpbnqoMaaBEpGtwG9xBtFNIxeyFoc7TW8DzoVgGThXM1+oqt/7pGmGM8A3XlV/FJHWZX2zxhhjGqZQtjiGAmmqutkd5JoGnOmX5iKcq2h/BLCgYYwxDV8oF6TrQMULqDJwpgn66glEicgsnOl+j6vqy/4Zicg1ONcekJCQMLh3794hKbAxxjRWS5cuzVTVVrWRVygDR6Apiv79Yk2AwcDJOFP9FojIQlXdUOEk1edwpj6SmpqqS5YsCUFxjTGm8RKRH2pO5U0oA0cGzhIOZTrizMv3T5PpzkfPE5E5OOvUbMAYY0yDFMoxjsVADxHpKiLROIuqTfdL8yFwgog0EZF4nK6stSEskzHGmMMUshaHqhaLyI3AZzjLJU9V1TUiMtk9PkVV14rIpzgre5YCL6jq6lCVyRhjzOELuwsAbYzDHMmKiorIyMigoKCgvotiGqjY2Fg6duxIVFRUhf0islRVvd59slqN7jafxjRmGRkZJCUl0aVLFzwukWWOIKrK7t27ycjIoGvXriF7Hlvk0JgwUlBQQEpKigUNE5CIkJKSEvIWqQUOY8KMBQ1Tnbr4+7DAYYwxJigWOIwxQYmMjGTgwIH069ePX//61+zfv7/mk1wvvvgiN954Y1DPl5gY+K7F99xzD19++SUAY8eOpWzSzIQJE8jKyiIrK4tnnnkmqOeqzh133EHfvn254w7/+5JVb8mSJdx88801phs5cuShFq3O2eC4MSYocXFxLF++HICJEycyZcoUbr311vLjJSUlREZGVnF27XnggQcC7p8xw7ndydatW3nmmWe4/vrra+X5nn32WXbt2kVMTExQ56WmppKaWvNkpvnz5x9q0eqctTiMMYfshBNOIC0tjVmzZnHiiSdy0UUX0b9/fwoKCrj88svp378/xx13HDNnziw/Jz09nfHjx9OrVy/uv//+8v1nnXUWgwcPpm/fvjz33HMVnue2225j0KBBnHzyyezatQuAyy67jHfeeadSmbp06UJmZiZ33XUXmzZtYuDAgdxxxx1MmjSJDz/8sDzdxIkTmT694jXJqsodd9xBv3796N+/P2+++SYAv/rVr8jLy2PYsGHl+8r079+frKwsVJWUlBReftlZbm/SpEl8+eWXzJo1izPOOAOA++67jyuuuIKxY8fSrVs3nnjiifJ8ylpWs2bNYuzYsZx33nn07t2biRMnUnbZxIwZM+jduzejRo3i5ptvLs+3rlmLw5gw9uijj1baN3jwYMaOHUthYSFPPvlkpeMjRoxg5MiR5Obm8uyzz1Y4dttt3m+nUVxczCeffML48eMBWLRoEatXr6Zr167l5Vq1ahXr1q3j1FNPZcOGDRXSxcfHM2TIEE4//XRSU1OZOnUqLVq0ID8/nyFDhnDuueeSkpJCXl4egwYN4tFHH+WBBx7g/vvv56mnnqqxfA8//DCrV68ubx3Nnj2bf/3rX5x55plkZ2czf/58XnrppQrnvPfeeyxfvpwVK1aQmZnJkCFDGD16NNOnTycxMbE8L1/HH38833zzDZ07d6Zbt27MnTuXSy65hIULF/Lvf/8b/+vO1q1bx8yZM8nJyaFXr15cd911la65+O6771izZg3t27cvzz81NZVrr72WOXPm0LVrVy688EJP71MoWIvDGBOU/Px8Bg4cSGpqKkcddRRXXnklAEOHDi2/dmDevHlMmjQJgN69e9O5c+fywDFu3DhSUlKIi4vjnHPOYd68eQA88cQTHHvssQwfPpz09HQ2bnTuPhsREcFvfvMbAC6++OLy9MEaM2YMaWlp7Ny5kzfeeINzzz2XJk0qfneeN28eF154IZGRkbRp04YxY8awePHiavM94YQTmDNnDnPmzOG6665j1apVbNu2jRYtWgQcnzn99NOJiYmhZcuWtG7dmh07dlRKM3ToUDp27EhERAQDBw5k69atrFu3jm7dupW/xvUZOKzFYUwYq66FEB0dXe3xxMTEoFoYZXzHOHwlJCSUb1e3IoX/dFERYdasWXz55ZcsWLCA+Ph4xo4dW+W1CIcz3XTSpEm89tprTJs2jalTp1Y6figraYwePZqnn36aH3/8kQcffJD333+fd955hxNOOCFget8xksjISIqLiz2laUirfFiLwxhT60aPHs1rr70GwIYNG/jxxx/p1asXAF988QV79uwhPz+fDz74gOOPP57s7GyaN29OfHw869atY+HCheV5lZaWlo9lvP7664waNcpTGZKSksjJyamw77LLLuOxxx4DoG/fvgHL/eabb1JSUsKuXbuYM2cOQ4cOrfZ5OnXqRGZmJhs3bqRbt26MGjWKRx55pMrAcah69+7N5s2b2bp1K0ClsZa6ZC0OY0ytu/7665k8eTL9+/enSZMmvPjii+XfokeNGsWkSZNIS0vjoosuIjU1lf79+zNlyhQGDBhAr169GD58eHleCQkJrFmzhsGDB9O0aVPPH5gpKSkcf/zx9OvXj9NOO41//OMftGnThmOOOYazzjor4Dlnn302CxYs4Nhjj0VE+Pvf/07btm1rfK5hw4ZRUlICOF1Xd999t+cA51VcXBzPPPMM48ePp2XLljUGtFCyRQ6NCSNr167lmGOOqe9ihK39+/fTv39/li1bRtOmTeu7OEHLzc0lMTERVeWGG26gR48e/O53v6uULtDfSW0ucmhdVcaYI8KXX35J7969uemmm8IyaAA8//zzDBw4kL59+5Kdnc21115bL+WwFocxYcRaHMYLa3EYYyoIty97pm7Vxd+HBQ5jwkhsbCy7d++24GECKrsfR2xsbEifx2ZVGRNGOnbsSEZGRvmyG8b4K7sDYChZ4DAmjERFRYX0zm7GeGFdVcYYY4JigcMYY0xQLHAYY4wJigUOY4wxQbHAYYwxJigWOIwxxgTFAocxxpigWOAwxhgTFAscxhhjgmKBwxhjTFAscBhjjAmKBQ5jjDFBscBhjDEmKBY4jDHGBMUChzHGmKBY4DDGGBOUkAYOERkvIutFJE1E7gpwfKyIZIvIcvfnnlCWxxhjzOEL2R0ARSQSeBoYB2QAi0Vkuqp+75d0rqqeEapyGGOMqV01tjhEpI2I/EdEPnEf9xGRKz3kPRRIU9XNqloITAPOPLziQsbefO75cPXhZmOMMeYQeemqehH4DGjvPt4A3OLhvA5Aus/jDHefvxEiskJEPhGRvoEyEpFrRGSJiCzZu7+Qlxf84OHpjTHGhIKXwNFSVd8CSgFUtRgo8XCeBNinfo+XAZ1V9VjgSeCDQBmp6nOqmqqqqR6e1xhjTAh5CRx5IpKC+6EvIsOBbA/nZQCdfB53BLb7JlDVfaqa627PAKJEpKWXghtjjKkfXgbHbwWmA91F5BugFXCeh/MWAz1EpCuwDbgAuMg3gYi0BXaoqorIUJxAtjuI8htjjKljNQYOVV0mImOAXjjdT+tVtcjDecUiciPO+EgkMFVV14jIZPf4FJwAdJ2IFAP5wAWq6t+dZYwxpgGpMXCIyDl+u3qKSDawSlV3Vneu2/00w2/fFJ/tp4CnvBfXGGNMffPSVXUlMAKY6T4eCyzECSAPqOorISqbMcaYBshL4CgFjlHVHeBc1wH8GxgGzAEscBhjzBHEy6yqLmVBw7UT6Kmqe4AaxzqMMcY0Ll5aHHNF5H/A2+7jc4E5IpIAZIWqYMYYYxomL4HjBpxgcTzOrKqXgXfd2U8nhrBsxhhjGiAv03EVeMf9McYYc4TzssjhcBFZLCK5IlIoIiUisq8uCmeMMabh8TI4/hRwIbARiAOuwllXyhhjzBHI0/04VDVNRCJVtQT4r4jMD3G5jDHGNFBeAsd+EYkGlovI34GfgITQFssYY0xD5aWrapKb7kYgD2fFW/9lSIwxxhwhvASOs1S1wF0C/X5VvRWwW70aY8wRykvguDTAvstquRzGGGPCRJVjHCJyIc79M7qKyHSfQ0nYPTOMMeaIVd3g+HycgfCWwKM++3OAlaEslDHGmIarysChqj8AP+AsqW6MMcYA3q4cP0dENopItojsE5Ecu3LcGGOOXF6u4/g78EtVXRvqwhhjjGn4vMyq2mFBwxhjTBkvLY4lIvIm8AFwoGynqr4XqkIZY4xpuLwEjmRgP3Cqzz4FLHAYY8wRyMv9OC6vi4IYY4wJD15mVfUUka9EZLX7eICI/F/oi2aMMaYh8jI4/jxwN1AEoKorgQtCWShjjDENl5fAEa+qi/z2FYeiMMYYYxo+L4EjU0S64wyIIyLn4SxFYowx5gjkZVbVDcBzQG8R2QZsAS4OaamMMcY0WF5mVW0GThGRBCBCVXNCXyxjjDENlZdZVX8VkWaqmqeqOSLSXET+UheFM8YY0/B4GeM4TVWzyh6o6l5gQshKZIwxpkHzEjgiRSSm7IGIxAEx1aQ3xhjTiHkZHH8V+EpE/oszs+oK4KWQlsoYY0yDVW3gEBEB3sC5498pgAB/VtXP6qBsxhhjGqBqA4eqqoh8oKqDgU/rqEzGGGMaMC9jHAtFZMihZC4i40VkvYikichd1aQbIiIl7sWFxhhjGjAvYxwnApNFZCuQh9Ndpao6oLqTRCQSeBoYB2QAi0Vkuqp+HyDd3wDr/jLGmDDgJXCcdoh5DwXS3AsIEZFpwJnA937pbgLeBQ6pVWOMMaZu1dhVpao/AJ2Ak9zt/V7OAzoA6T6PM9x95USkA3A2MKW6jETkGhFZIiJLyvZl5h6o7hRjjDEh4uXK8XuBO3GWVgeIwpmiW+OpAfap3+PHgDtVtaS6jFT1OVVNVdXUsn35hdWeYowxJkS8dFWdDRwHLANQ1e0ikuThvAyclkqZjsB2vzSpwDRn1i8tgQkiUqyqH3jI3xhjTD3wEjgK3Wm5ZcuqJ3jMezHQQ0S6Attwbv50kW8CVe1ati0iLwL/8xo01L/tYowxpk54Gat4S0SeBZqJyNXAlzh3BayWqhYDN+LMlloLvKWqa0RksohMPpxCA2ilXi9jjDF1wcuy6o+IyDhgH9ATuEdVv/CSuarOAGb47Qs4EK6ql3nJ82D6YFIbY4ypLV66qgBWAXE4g9urQlcc7yxuGGNM/fAyq+oqYBFwDnAezpXkV4S6YDVRa3IYY0y98NLiuAM4TlV3A4hICjAfmBrKgtXEwoYxxtQPL4PjGYDv7WJzqHhhnzHGmCOIlxbHNuBbEfkQ54v+mcAiEbkVQFX/GcLyVcl6qowxpn54CRyb3J8yH7q/vVwEGEIWOYwxpj54mY57f10UJFjW4jDGmPrhZYyjQbK4YYwx9SN8A4dFDmOMqRfhGziszWGMMfWixjEOEWkFXA108U2vqvV+EaAxxpi652VW1YfAXJzFDRvMTTCsq8oYY+qHl8ARr6p3hrwkQfISOEpLlZyCYprGR4W+QMYYc4TwMsbxPxGZEPKSBMnLGMfDn67j2Ac+Z+e+gjookTHGHBm8BI7f4gSPfBHZJyI5IrIv1AWriZcWx3NzNgOQvjc/xKUxxpgjh5cLAOv5CvFDU1p6MLLY/cmNMab2VBk4RKS3qq4TkUGBjqvqstAVq2Y1tTgOFJeWb+cUFFU4dvvbK0jfs5+Lhh3FmQM7hKJ4xhjTaFXX4rgVuAZ4NMAxBU4KSYk8qmmMo6DoYCvjpQVbOa1/u/LH7yzNAODbLXv4w3urOLpNEh/ecHxoCmqMMY1MlYFDVa9xf59Yd8XxrqYWR75P4Fi4eU+V6fIKS1iRnkVpqRIRIbVVPGOMabS83AFwrog8KCLjRSRsxjt8Wxy+9vl1W5X515cbQlkcY4xpNLzMqroUWA+cC8wXkSUi8q/QFqtmNU2qyvcLHGW3mt2XHzhwvLHI7k1ljDFeeJlVtVlE8oFC9+dE4JhQF6wmpTX0VRUUOYPjLRNjyMw9wL5850LAzNzCgOkzcw9QUFRCbFRkrZfVGGMaEy9rVW0CMoHXgf8AN6lqafVnhd6TL79Hv7gsAAYPHszYsWMpLCzkySefBCC9MB44irbFP5NJc3blHiCy9ABXTfkKiGVQ3G42FyaSVRJTnmfvP33KTa3WVXqucePG0aZLT1Zt3saKrz6odHzChAkcc8wxpKen89Zbb1U6ftZZZ9G9e3c2bdrEBx9UPv/888+nU6dOrF27lhkzZlQ6fvHFF9OmTRtWrlzJF198Uen4FVdcQfPmzVmyZAmzZ8+udPzaa68lMTGR+fPns2DBgkrHb7rpJqKjo5k1axZLly6tdPy2224D4PPPP2fVqlUVjkVFRXHzzTcD8PHHH7NuXcXXLyEhgcmTJwPw/vvvs3nz5grHmzdvzhVXOMuevfXWW6SnV2z5tWnThosvvhiAV199lR07dlQ43qlTJ84//3wApk6dyt69eysc79atG2effTYAU6ZMIS8vr8Lx3r17c/rppwPwxBNPUFRUsUXav39/Tj31VAAefbTyPJFAf3u+RowYwciRI8nNzeXZZ5+tdHzMmDGkpqayd+9epk6dWun4uHHjGDBgADt27ODVV1+tdNz+9uxv71D/9g6HlyVHngBGARcCxwGzRWSOqm6q/rTQmpnbtjxwBFKszkB3UqTzZuzOPUDb+BjiI4qhBLrH5NAlOpf3sjtXOE8VxG+M/OYZP7M1y/mjujolgtiIeo+bDU72/sBdgLVBVVmRkc2GnCYkB3h/jDF1S9TjaoEikghcDtwOdFTVeunTiWnXQ9td+hgAWx8+vcp0H6/8iRteX8bjFwzkt9OW88zEQUzo347Uv3xBZm4hWx6agLifQD9l5zPioa8B+PYPJ9MmObZCXl3u+rjC4w1/OY3oJhGUlCqRR9hMLFWlpFTJKywhITqSyAih690Vv6me3Ls1T08cxMOfrCMuOpLbxvWkSaS3Ffw37colbWcuv39nJf84bwD5RSX8dtrySuluOLE7k8d0JynW1iEzxgsRWaqqqbWRl5euqkdxWhyJwALgHpzVchu0ssHxjs3jAbj+tWUsv2dc+RiH+Hxtbdc0jnOO68B7321j087cCoEjUGC96Y1lnJ/aiStfWsJXt42he6vEUFal3qgqaTtzKVFle1Y+3Vomcv1ry/j+p+pXnPlq3U56/+nT8sdpO3N5/pKKf6/5hSVERQp/+nANAzs15c53V/lnwzWvVO66KPP0zE08PXMTvdok8e+LB3FUi3hEhAiB/YUlfLl2Byf0aMV/5m3m6ZmbmHLxYMb3a1t+fmmpUlBcwpC/fElegJUFEqIjGdWjJYkxUZw7qAMjj25ZbZ2NOZJ46apaCPxdVXfUmLIBKClVbntrOfExTtXaNzsYBHIKiqs8787TevPed9u46IVvWfzHU2iV5Ix9fLV2JwD/d/oxbM7M4/Vvf2T1tn18u2UFAC/M3cxD5wwIVXXq1aT/LGJeWqantCvvO5XP1+xgzoZdbNyZy1o3uMRGRfDF9zvoctfH3H5qT0SEpT/s5et1O8vPfWNR5fwuGNKJaYud7sFrR3fjrtN6IyJk7y9i6Y97+N2bK8jOL2L9jhxOerRy37q/ya9WHYQCySss4bM1zp/8u8syyvf/ccIxXHZ8F6I8tqCMaYy8zKp6uy4KUlt+3LOfD5ZvL3+cEHOwinmFVQeO1kkHB8m/3bKbMwa0Z/3POVz18hLneHIsV53QjSVb97BhR2552jcWpXP/r/oR3SS8P0hUlZ/3FRAdGcHX63bSpWVClUHjmYmDSEmIplurRF6Yu5mxvVqTHBvFeYM7ct7gjsDBFoWIMOYfM8nYm88jn1d9rcwlIzrzwJn9Kux7+NwBlboDm8ZHcVLvNqy491RUlVcW/sA9H66ptm5vTx7B3e+tIm1nbqVjzeOjePDs/uQUFNElJYGhXVtQUqp8/9M+/vy/79mZc4CUhGiW/ZgFwIMz1vLgjLUAnDmwPX8/bwClpRAXbbPxzJHDS4sjrGzPqrgSbpzP9NrMHKeb6uUrhlY6z7fr6qcsZxn2a15ZUr6vVaITWNokx1YIHADjH5/D17eNPbyC17Nz/z2//MPR10tXDGVMz1YALNi0G0UZ2f1gt83dEwLPzPb9IJ1350l8uHwbv522nObxUSTHRTGwUzN+d0pPOjSPq/bbe3VjSCLCJSO6cMmILhQUlRAhQnSTiPIFLiMiBFVFRPj8ltHsLyqhSYQQGSHVPmeTSGFAx2a8PXlk+b5dOQdIim3CI5+t54V5WwD4cPl2PnS/pFw49ChO7t2a9s3iuO3tFRSVlFYIVP+9fAgn9mpd5XOaxq/sb7EqM9fv5JuNmazclk33VonER0eyJTOPpnFRnN6/HaldmrM5M4+tmXnszDlA66QYerdNJr+ohPQ9+5mzcRerMrIpKC6hXdM4erdNonfbZHbsK2Duxl21WhfPg+MNRU2D4795dgHfbjm4xMjWh0/n7SXp3PHOSpJimpBzoJh3rxvJ4M7NA+af+pcvaJkYw4juKfz3m63l+7++bQzdWiWSd6CYvvd+BsAtp/TgsS83VlmW2vb1uh1c8eISpl0znOHdUgKmWb0tm335RZ765EtLlUv/u4joyAi+8uk68lUX9Qo3JaVK7oFiZq3fGXDgvjrHH53CkC4t6NUmieJSZWyvVlz8n0WsSM8CoGlcFFeN6srVo7sRGxWJqvLjnv3kFBRzTLtk1v+cw4vzt3BMu2TOT+1UoUXtRXUfXiWlSk5BEc3io4PKs66VvSbtmsbVeks/p6CIz9bsYH9hMZm5hXRrmcD8TZmcMaA9w7ullD/fz9kFPDtnEz3bJHFsx2a0To4ha38hS7buJWNvPisysmgeH02vtkl892MWM9fvpKRUiYoUmsZF06FZLBl78ykqKUWpvhs9GL3bJpEQ04StmXns2V9YYWmmH/52Rt0NjgOIyCigh6r+170HeaKqbqmNAtQ236BRprU72J1zwHlzkmKrrnab5FjWbN/Hup9zyvf959JUurkD4L7/qO19BtW3ZeXToVlcrdQhkM27crniRacFdMFzC+nTLpkZvz2hUrpfPTWPUoV/nn8s5wzqiKry+FcbmbV+F5cf34UxPVuRENOEt5dk8If3Kw5IP3h2P47t2IyOzeP44/uruW5s95DVJ5xFRghN46I4c2AHxvRsRXJsFO99t4273l1JcakyumcruqTEc/nxXenYPI7s/CLu/+h7PlqxnW/SdvNN2u4q887OL+LRLzbw6Bc1L4Fz/0ffA5Ac24Q/ndGHU/u05fm5mzk/tRMdm8dVWnvtkc/W8+aSdNo3iyMnv4jNmXmBsgVgXJ829G6bxPh+beneKpH0PfvJyi8iKjKC6MgIPl3zM91aJnBKnzYs/WEvI7unVNmKU1X2FRTTNC74GXDFJaU88dVGTunThk9W/8yrC3+o8CEbHx3J8G4pzN6wi2FdW5CZe4D9hSXs3HeAkUenMHlMd1olxVSYwFJQVMLbS9LZlVvIO0vSKSwppVVSLEe1iGN+2u7yzwlfby1xxrkiI4SS0iC+bK84uBkXFUnz+Ci2ZxeQmXsAcL4olK1mMbBTM0Z2T2FC/3aUlCpzN+5iUOfm9GidxMcrtzNrwy46No9jTM/WdEmJJzO3kI07c2gSEUG7ZrEM75pSoaWfe6CYtJ25JMc2oVurRORv3otdkxpbHCJyL5AK9FLVniLSHnhbVetlOdmaWhz+U2e3Pnw6qzKy+eVT88r3fXPXSVV+yPuf/9a1IxjatUWFfcvTs/jrjLVMu3o4G3fm8ovH5vD3cwfQp30y01dsJ7ZJBLec0rP8H3fnvgI+W/Mzvxly1CF/Q/IvF8DsO8bSOSWh/HH2/iKOfeDz8sfz7jyR615dxqpt2Z6ew3eKsqldqkrG3nzeW7aN2Rt2Eh/dhH0FRazMyKZPu2Q+umkUAszeuIu3l6Tz9bqd5asfAHRtmcAW94P+4uFHsW1vPjPXB9/9kBjThNwAH4xlRLzdJK0mbZNjySkoqjBjbWjXFvzp9D7079gUoHxh0bSduWzYkUOrpBiObpXIX2esZUfOAVrER1UYr/TVvmksIsK2LO83aaup7mV+PbgjXVslEB0ZwdbdeQzs1Jw9eQf4cPl21mx3Jn10aBbHLaf0IL+ohG/SMtlfWELW/iI6NIsjPjqS0T1bkdqlOQs27SYuOpIJ/dpVCOSFxaUUlpSSGGSL8XDU5nRcL4FjOc6Ff8tU9Th330pVrZepRIcSOHbsK2DYX78q37fi3lOr/Pbzzy828MRXG8sfL/vTOFokVN10V1X63vsZ+wNM6Ty5d2uaJ0Tz3Y972bTL+af3nxZaE1Vl6jdb+fP/vi+vz+wNu7h0qjMVac39vyAhpgmvf/tjeQvipN6tK8xaqkqHZnG8dMVQjm7dOKcThzPfMZuaLP1hD5+s+pkX5m2hb/tkmkRGlHd9+UpJiGbmHWOZvnw7+wuLGdm9JS0TY2gWH0VMk4jyLw2qyjdpu/n+p2ymLUrnxz37KS5VjmmXjLpTs8f2as3uvAMs3LyH0T1asvanHH6u5hbNiTFNKHGnQLdKjCG/qMRz90zH5nEcd1RzVJVxfdpwxoD2Fca+0nbm0jQuipJSJSICYiIjSYptQm5hMf+Zu4Xl6Vms3pbN7ryDyw2N69OGbi0TuGRkFzo0c1qFO/cVEBkh5b0LjU2dXscBFKqqioi6T55Q0wn1Zbfb/PPXJjmWP5/Vjz99sBpw5uhX5ZaTe3DOcR0Y+8gswJl1Ux0RCRg0gIDjBpNfXcox7ZL5JEA3UyBPfJVWvnLvzSf3AGBMz1YM7dqCRVv28OTXaUyZXfEi/v9cmspZz8wv//C45ZQe3HJKTw4Ul/D12p3c/f4q2ibH8uktoz2VwdS9YNZMG9y5BYM7t+D/zuhTvq+0VMkpcNZn83fx8M6V9vkSEUb1aMmoHi25ZvShdVfmHSimqKSUqMiI8u7dfQVFPD0zjWdnV1z6Iy4qkj7tk+neKoE5GzI5tW8bTuzVmm/SMvlFv7YM6dIi0FOUq+qLT3JsFL8b17PCvqou2m0aF3VIXWlHKi8tjtuBHsA44CHgCuB1Va28OEodqK7FsTUzr/wDH+DEXq347+UHZ1CVtUa8DPiuyshmf2Exw6oYhPb1TVomE1/41kPpD/roxlFsz85nWNcW1Q5G+rag1v9lPDFNnA8UVeWSqYuYu7HilNnXrxrGyKNboqoc/cdPKClV1v15fIUPoppmdxgTSgeKS9iTV0i7pqEbEzSV1WmLQ1UfEZFxwD6gF3CPqlZe7SwAERkPPA5EAi+o6sN+x88E/gyUAsXALao6r1JGHvl/Ft7+i16HmlV5P6wXx/vMYNr44GlERUawY18By9OzeHXhD8zdmMnrVw9jeXoWz87eTHZ+UYUxl2bxUfzvplHlV7m/NH8rrZNi2LTr4HRO36ABzrfCG088ujxwPHHhcfzq2PYVjm/664SA5bWgYepTTJNICxphzsuSI7/DGQz3FCx8zosEnsZpqWQAi0Vkuqp+75PsK2C62xU2AHgL6O31Ofzv2ldUUnHxQf+BpzeuHs7mzMoXgdWGd68bwZKte8tnlrRJjuUXfdtyQo+WLP8xi5HdWzKye0uuH3t0pXGYrP1FjPrbTDb9dQLL07O4d3rFC9r+PXFQhaBRZli3FL66bQwdmsXZcvDGmDrjZYwjGfhMRPYA04B3PC4/MhRIU9XNACIyDTgTKA8cqur7KZ5AzfdnqqBElQgOBo4DxU7g+OvZ/QEqzDgCGNE9hRHda+56OhRl/cz+4qObVLqmYsW9pzJ34y76tEvm5H/OLp/F8tyczfzt08rLuvveL91fY10nyxjTcNU4ZUNV71fVvsANQHucZdW/9JB3B8B3gfsMd18FInK2iKwDPsYZP6lERK5x7zy4xHe//3zqohLncdumMVw07CgPRawfTeOiOGNAe7q1SmTLQ6ez8cHTaJkYEzBofP/AL+qhhMYYU7VgJhHvBH4GdgNe1k4I1JFeqUWhqu8D74vIaJzxjlMCpHkOeA6cwfGy/f53ASx0WxzRkeHVbRMVGcHQrs2ZsepnADa7YxP+F3AZY0xDUGOLQ0SuE5FZOOMRLYGrPV7DkQF08nncEQh8NQ+gqnOA7iLief1q/xbHnA3OBVFRkeH3gfvQ2QPo1jKBj28eRUSEWNAwxjRYXlocnXFmOy0PMu/FQA8R6QpsAy4ALvJNICJHA5vcwfFBQDROi8aTUr8b8T01Mw0gLFeqbRofxde3j63vYhhjTI2qDBwikqyq+4C/u48rjPyqauVFoSoeLxaRG4HPcKbjTlXVNSIy2T0+BTgXuEREioB84DcaxKqLJVUktXslGGNM6FTX4ngdOANYijM24dt3okC3mjJX1RnADL99U3y2/wYc8tJbVS02FhOGLQ5jjAkXVQYOVT3D/d217ooTnKoCh7U4jDEmdLwMjn/lZV998O2q8g0i1d38xxhjzOGpbowjFogHWopIcw52VSXjXM9R70p9goX/VePGGGNCo7oxjmuBW3CCxFIOBo59OEuJ1DvfVsYOd0nnPu2S6dQivr6KZIwxjV51YxyPA4+LyE31tRJuTXy7qs55Zj4Av07tWF/FMcaYI4KX1XGfFJF+QB8g1mf/y6EsmBe+XVVlN2kJ6raOxhhjguZlddx7gbE4gWMGcBowD6j3wBHoOo6yhQ6NMcaEhpd5q+cBJwM/q+rlwLFATEhL5VFxSeXAUVAU+G58xhhjaoeXwJGvqqVAsYgk4yx2WOPFf3XBf5FDsMBhjDGh5mWtqiUi0gx4Hmd2VS6wKJSF8irQeEZBkXVVGWNMKHkZHL/e3ZwiIp8Cyaq6MrTF8sa3xdG1ZQJbMvM4b7DNqjLGmFCq7gLAQdUdU9VloSmSd77X/G3JzAPg2E7N6qcwxhhzhKiuxfFoNccUOKmWyxK0sq6qVRnZ9VwSY4w5clR3AeCJdVmQQ1HWVfXLp+bVc0mMMebI4eU6jksC7W8IFwAW28V+xhhT57zMqhrisx2Lc03HMhrABYClFjiMMabOeZlVdZPvYxFpCrwSshLVwHfBdP8Vca8b271uC2OMMUegQ7nj0X6gR20XxCvfe234d1W1bxrrn9wYY0wt8zLG8RHOLCpwAk0f4K1QFsor/xaHdVwZY0zoeRnjeMRnuxj4QVUzQlSeoBQWl1YY5zhgV40bY0zIeRnjmA3grlPVxN1uoap7Qly2gEQqdlWt2b6v/LGtU2WMMaHnpavqGuDPQD5QijM+rTSAhQ6LSkorXMNhS6obY0zoeemqugPoq6qZoS5MsIr8llW3FocxxoSel1lVm3BmUjUI1U3Htcs6jDEm9Ly0OO4G5ovIt8CBsp2qenPISlUdn8ixZVde+XZ0ZAS/PaXeZgkbY8wRw0vgeBb4GliFM8ZRr3xbHG8uSS/f/tMv+9A0LqruC2SMMUcYL4GjWFVvDXlJDlOkSM2JjDHGHDYvYxwzReQaEWknIi3KfkJesiAFuo2sMcaY2uelxXGR+/tun331Oh23Q7M4tmXlV9hnU3GNMaZueLkAsGtdFCQYAzs1qxQ4bCquMcbUjbC+H4cva3EYY0zdCOv7cfg6UGwtDmOMqQthdz+OqtgCh8YYUzfC7n4cVfnlse3quwjGGHNECOn9OERkPPA4EAm8oKoP+x2fCNzpPswFrlPVFd6KftDFw49icOcGN0PYGGMapZDdj0NEIoGngXFABrBYRKar6vc+ybYAY1R1r4icBjwHDPNcepdd/GeMMXWnysAhIkcDbcrux+Gz/wQRiVHVTTXkPRRIU9XN7nnTgDOB8sChqvN90i8EOnoqtV+ciIiwwGGMMXWlujGOx4CcAPvz3WM16QCk+zzOcPdV5Urgk0AH3CvXl4jIkpKSyrOnrMVhjDF1p7rA0UVVV/rvVNUlQBcPeQf6NA+4LoiInIgTOO4MdFxVn1PVVFVNjYyMrHQ80locxhhTZ6ob44it5lich7wzgE4+jzsC2/0TicgA4AXgNFXd7SHfSqyryhhj6k51LY7FInK1/04RuRJY6iHvxUAPEekqItHABcB0v7yOAt4DJqnqBu/Frsi6qowxpu5U1+K4BXjfnTJbFihSgWjg7JoyVtViEbkR+AxnOu5UVV0jIpPd41OAe4AU4BlxPvyLVTXVS8FbJESzJ68QgBJbGdcYY+pMlYFDVXcAI93xh37u7o9V9WuvmavqDGCG374pPttXAVcFVWLXrDvGMuC+zwFb4NAYY+qSlyVHZgIz66AsQUmKOVj0AltuxBhj6syhLDlS7wQQn3GNA9biMMaYOhOWgcNfga2Ma4wxdaZRBA5bGdcYY+pOowgc1uIwxpi60zgCh7U4jDGmzjSKwNEmOaa+i2CMMUeMRhE4Hj53QH0XwRhjjhhhGTjEb4mR5NioeiqJMcYcebzcyKnB+vjmUWTmFtZ3MYwx5ogS1oGjb/um9V0EY4w54oRlV5Uxxpj6Y4HDGGNMUCxwGGOMCYoFDmOMMUEJu8Bht2wyxpj6FXaBA5xl1Y0xxtSPsAwcxhhj6o8FDmOMMUGxwGGMMSYoFjiMMcYExQKHMcaYoIRf4LD5uMYYU6/CL3AAYvNxjTGm3oRl4DDGGFN/LHAYY4wJigUOY4wxQbHAYYwxJigWOIwxxgTFAocxxpigWOAwxhgTlLAMHHYZhzHG1J+wDBzGGGPqjwUOY4wxQbHAYYwxJighDRwiMl5E1otImojcFeB4bxFZICIHROT2UJbFGGNM7WgSqoxFJBJ4GhgHZACLRWS6qn7vk2wPcDNwVqjKYYwxpnaFssUxFEhT1c2qWghMA870TaCqO1V1MVAUwnIYY4ypRSFrcQAdgHSfxxnAsEPJSESuAa5xH+Y+fuGg3Y9fSOZhlq+hagmNtm5g9Qt3jbl+jbluAL1qK6NQBo5Al1sc0m2YVPU54LnyjEWWqGrqoRasIWvMdQOrX7hrzPVrzHUDp361lVcou6oygE4+jzsC20P4fMYYY+pAKAPHYqCHiHQVkWjgAmB6CJ/PGGNMHQhZV5WqFovIjcBnQCQwVVXXiMhk9/gUEWkLLAGSgVIRuQXoo6r7asj+uRqOh7PGXDew+oW7xly/xlw3qMX6ieohDTsYY4w5QtmV48YYY4JigcMYY0xQwipw1LSESbgQka0iskpElpdNkRORFiLyhYhsdH8390l/t1vn9SLyi/oreWAiMlVEdorIap99QddHRAa7r0uaiDwhIvW+gn4VdbtPRLa5799yEZngcyxs6gYgIp1EZKaIrBWRNSLyW3d/2L9/1dStUbx/IhIrIotEZIVbv/vd/aF/71Q1LH5wBtg3Ad2AaGAFzkB6vZftEOqyFWjpt+/vwF3u9l3A39ztPm5dY4Cu7msQWd918Cv7aGAQsPpw6gMsAkbgXAP0CXBaA63bfcDtAdKGVd3ccrUDBrnbScAGtx5h//5VU7dG8f65ZUl0t6OAb4HhdfHehVOLo8YlTMLcmcBL7vZLHFy/60xgmqoeUNUtQBrOa9FgqOocnHXHfAVVHxFpBySr6gJ1/pJfpgGsYVZF3aoSVnUDUNWfVHWZu50DrMVZ9SHs379q6laVsKkbgDpy3YdR7o9SB+9dOAWOQEuYVPdH0JAp8LmILBVnORWANqr6Ezh/8EBrd3+41jvY+nRwt/33N1Q3ishKtyurrCsgrOsmIl2A43C+uTaq98+vbtBI3j8RiRSR5cBO4AtVrZP3LpwCR60tYdIAHK+qg4DTgBtEZHQ1aRtTvaHq+oRTPf8NdAcGAj8Bj7r7w7ZuIpIIvAvcotVfRxV2dQxQt0bz/qlqiaoOxFmZY6iI9Ksmea3VL5wCR6NZwkRVt7u/dwLv43Q97XCbjLi/d7rJw7XewdYnw93239/gqOoO9x+2FHieg12HYVk3EYnC+WB9TVXfc3c3ivcvUN0a2/sHoKpZwCxgPHXw3oVT4GgUS5iISIKIJJVtA6cCq3Hqcqmb7FLgQ3d7OnCBiMSISFegB85AVkMXVH3cJnWOiAx3Z3Rc4nNOg1L2T+k6G+f9gzCsm1ue/wBrVfWfPofC/v2rqm6N5f0TkVYi0szdjgNOAdZRF+9dfc8MCHIWwQScmRGbgD/Wd3kOsQ7dcGY2rADWlNUDSAG+Aja6v1v4nPNHt87raQCzOQLU6Q2cJn8RzreXKw+lPkAqzj/xJuAp3JUNGmDdXgFWASvdf8Z24Vg3t1yjcLolVgLL3Z8JjeH9q6ZujeL9AwYA37n1WA3c4+4P+XtnS44YY4wJSjh1VRljjGkALHAYY4wJigUOY4wxQbHAYYwxJigWOIwxxgTFAodplERkloik1sHz3OyuvvpaqJ+riufvIj4r9x5GPi+KyHm1UaYq8q+VcpqGIWS3jjUmXIlIE1Ut9pj8epz58FtqKT9jGjxrcRhP3G+Ma0XkeXft/8/dq1VDkqdvi0FEWorIVnf7MhH5QEQ+EpEtInKjiNwqIt+JyEIRaeHzFBeLyHwRWS0iQ93zE9yF7Ra755zpk+/bIvIR8HmAst7q5rNaRG5x903BuaBzuoj8zi99hfzEuUfCB+7CegtFZICb7j4Rud3nvNXu61LdazNYnHswLABu8Dm3rzj3Z1juPk+PAPXIFZFHRWSZiHwlIq0CpLnHfX1Wi8hz4uguIst80vQQkaU+5ZktzqKdn8nB5S4CltOEPwscJhg9gKdVtS+QBZzrn0BEJsrBG+T4/rxzqHkG0A+4CGeNoQeB/ap6HLAAZ7mEMgmqOhKnVTDV3fdH4GtVHQKcCPxDnKVfwLkfwaWqepJfnQYDlwPDcO53cLWIHKeqk3HW9DlRVf8VoJy++d0PfKeqA4A/4CxdXZOqXpv/Ajer6gi/9JOBx9VZ9C6ViiuelkkAlqmzyOZs4N4AaZ5S1SGq2g+IA85Q1U1AtogMdNNcDrwozlpQTwLnqepgnNf5wRrKacKcdVWZYGxR1eXu9lKgi38CVX0NCKa/v8Y8A5ipzv0VckQkG/jI3b8KZxmGMm+4ZZojIsnirOtzKvArn2/5scBR7vYXqhro3hujgPdVNQ9ARN4DTsBZ7qE6vvmNwv3gV9WvRSRFRJrWcH6l18Y9p5mqznb3v4KzyjI4gfOPItIReE9VNwbIsxR4091+FXgvQJoTReT3QDzQAmdpnI+AF4DLReRW4Dc4gbsXTiD/QpybxkUCP9VQThPmLHCYYBzw2S7B+TZagYhMBO4IcG6aqgYafK0qz2IOtohjqzmn1OdxKRX/pv3X0ylbQvpcVV3vV+5hQF6A8kHgZae98M2vqqWrfesJFesa6LURqljyWlVfF5FvgdOBz0TkKlX9uoYyVshLRGKBZ4BUVU0Xkft8yvQuTgvla2Cpqu4WkfbAGv9WhRukbT2jRsq6qkytUtXXVHVggJ9gZ+xsBQa724c62+c3ACIyCshW1WzgM+Amcb8ei8hxHvKZA5wlIvFut9bZwNwgyzIHmOg+51ggU517Q2zFuTUtIjII55aeVVJn+exst06U5eme3w3YrKpP4CzeN6ByDkRw8PW8CJjnd7wsSGSKcx+L8tdeVQtwXr9/43RDgbNYXisRGeGWIUpE+lZXThP+rMVhGqpHgLdEZBLON9xDsVdE5gPJwBXuvj8DjwEr3eCxFTijukxUdZmIvMjB5exfUNWauqn83Qf8V0RWAvs5uOz1u8Al4tzFbTHO6s81uRyYKiL7cT7Iy/wGZ0JAEfAz8ECAc/OAvu7AdrZ7TjlVzRKR53G6/ba6ZfL1GnAO7gQCVS0UZxrvE273VBOc13dNNeU0Yc5WxzXmCCIiuaqaeBjn3w40VdU/1WKxTJixFocxxhMReR/nlqsn1ZTWNG7W4jDGGBMUGxw3xhgTFAscxhhjgmKBwxhjTFAscBhjjAmKBQ5jjDFB+X+0SU5rBLE8cAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "plot_limit = 3000 #plot how the ratio changes oby increasing $n$\n",
        "x = np.arange(1, plot_limit + 1) \n",
        "y = ratios[:plot_limit]\n",
        "plt.title(r\"Win ratio over (n) if we say 'no'\") \n",
        "plt.xlabel(\"n = number of rounds played\") \n",
        "plt.ylabel(r\"Cumulative win percentage\") \n",
        "plt.axhline(y = 1/3, color = 'k', linestyle='--', alpha = 0.6, label = r\"Probability of wining\")\n",
        "plt.plot(x,y) \n",
        "plt.xlim(-5,plot_limit)\n",
        "plt.ylim(0.1,0.6)\n",
        "plt.legend()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "276cca7b",
      "metadata": {
        "id": "276cca7b"
      },
      "outputs": [],
      "source": [
        "picked_door = 1 #fix a door (0,1, or 2)\n",
        "yes_no = \"yes\" #let's start with the strategy of not keeping our door!\n",
        "simulate = 30000 #number of times to simulate\n",
        "results = np.zeros(simulate) #a lot of zeros\n",
        "for i in range(simulate):\n",
        "    results[i] = simulate_round(picked_door, yes_no)\n",
        "ratios = np.cumsum(results) / (np.arange(1,simulate+ 1)) #running ratios. ratios[10] would give us the percentage of times we won by round 10"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "59ecc212",
      "metadata": {
        "id": "59ecc212",
        "outputId": "6d02a8c1-e2de-425a-8c74-cd3ad6719013"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6wklEQVR4nO3dd5xU1f3/8dfbFaSqCIhRaRosARRhQY2ooNEQu2IXjCUiGmNiS/Sbb2yJPxPLN/YChthbrJhYQBER0UgRKQqCiLCiKNKRtvD5/XHusrOzM7N3lzvb+Dwfj3nMLeece+7M7v3MOffec2VmOOecc0nYqqYr4Jxzrv7woOKccy4xHlScc84lxoOKc865xHhQcc45lxgPKs455xLjQcWVI+lgSTNrcPvtJK2UVFBTdagqSU9JOiFm2hck9cuxvsz3IGlPSR9JWiHp0gSq61ziPKhsASRdI+nVtGWzsiw73czeNbM9q7F+cyX9rGTezOaZWTMz21BddUiCpH2AfYGXY2b5K3BTtpUZvoffA6PNrLmZ3VX1mtYukkZL6pPH8h+WdE6+yndleVDZMowBDir55S9pJ6AB0D1t2Y+jtImRtHWS5dUGOfbpQuAJi3lHsZl9CGwrqTDmptsD02Omda5GeFDZMownBJFu0fwhwNvAzLRln5vZAkl9JBWVZI5aEldKmiJpmaRnJDXKtCFJ50h6T9LfJS0Grpe0u6RRkr6XtEjSE5K2j9I/BrQDXom6vH4vqYMkKzl4S9pZ0nBJiyXNlnRBth2VtJ2kRyV9J+lLSf8raStJ20haKqlLStrWklZL2jGaP0bS5CjduKjlkfoZ/EHSFGBVlsDyC+CdtM9irKTbJC2R9IWkX6TlGQ0cnWVfNn0PkkYBfYF7os9pj7S0fSVNTZl/U9KHKfNjS7rlos/z+egz+iJXV5qkoyR9EnW5fSXpymh5C0n/jspYEk3vGq07RdLEtHKukPRStu1EaXaS9IOklinLekTbaBDNnyfp02ibb0hqHy1X9Df3bfQ3OiX1u3bVyMz8tQW8CEHksmj6HuA8QtdL6rJh0XQfoCgl71zgQ2BnYAfgU2Bwlu2cAxQDvwG2BhoTWkBHANsArQmtoTvSyv9ZynwHwICto/l3gPuARoQg+B1weJbtP0rofmoelfMZcH60bhhwU0raXwOvR9PdgW+B/YEC4JdRvbZJqeNkoC3QOMN2m0Z1bp32WawHLojKvAhYACglzeXAC1n2Jf17GA38KkvaRsBqoFX0uX8Tbat59B2sBloSfkhOBK4FGgK7AXOAn2cp92vg4Gi6BdA9mm4J9AeaRNv4F/BStG4bYDGwd0o5HwH9Y/ydvgpclDL/d+DuaPoEYDawd7SP/wuMi9b9PNqv7QFFaX5U0/93W+Krxivgr2r6ouF64MVo+mOgE9Avbdkvo+n0g9lcYEDK/C3AA1m2cw4wr4K6nAB8lFZ+xqASHcQ3AM1T1t8MPJyh3AJgLfCTlGUXEs5DAPwMmJOy7j3g7Gj6fuDPaeXNBA5NqeN5OfZpl6jOjdI+i9kp802iNDulLLsAGJWlzPTvYTRZgkq0/l3gJOAAYATwbPQd9wWmRGn2T/9+gGuAf2Ypc170GW5bwXfaDViSMn8/UQAHOgNLiAJ0BeWcBryX8n1+A/SK5l8j+oEQzW8F/EDoFjyM8APiAGCr6vq/8lf5l3d/bTnGAL0ltSD8mp4FjAN+Gi3rQu7zKd+kTP8ANMuRdn7qjKQdJT0ddZ8sBx4n/KKOY2dgsZmtSFn2JeEgnq4V4df3l1nSjgIaS9o/6jbpBrwYrWsPXBF1fS2VtJQQ0HbOtl9plkbvzdOWb/rczOyHaDL1s2uekndzvUMIRIdE06OBQ6NXSbdce2DntP38H6BNljL7A0cBX0p6R9KBAJKaSHow6mJcTvjb2V6lV+w9ApwpScBA4FkzWxtjH14GfiJpN0LrdpmFc08ldb8zpd6LCa2SXcxsFKG1fS+wUNIQSdvG2J5LmAeVLcf7wHbAIMIvdMxsOaGLZBCwwMy+SGhb6Seqb46W7WNm2wIDCAeDbOlTLQB2kJR6sG4HfJUh7SJCd1P7TGnNbCPh1/sZwJnAv1OC1XzCL+vtU15NzOypOPU0s1XA58Ae2dJksTehlZiE9KDyDuWDynzgi7T9bG5mR2Uq0MzGm9nxwI7AS4TPD+AKYE9g/+g7PSRarijfB8A64GDCZ/1YnB0wszXRNs4iBKPUfPOBC9Pq3tjMxkV57zKzHoSW0R7AVXG26ZLlQWULYWargQmEPvx3U1aNjZYletVXmubASmCppF0o/8++kNC3X46ZzSe0qG6W1Cg6eX4+8ESGtBsIB6SbJDWPWiOXE1pGJZ4kdLGcFU2XGAoMjloxktRU0tFpwawirxIO4JVxKKFbJwnjCAf6XsCHZjadEGD3p/T7/RBYHl100FhSgaQuknqmFyapoaSzJG1nZuuB5YSuSAjf6WrCd7oDcF2G+jxKaD0Um9nYSuzHo4Suw+Mo+909AFwjqXNUv+0knRJN94y+uwbAKmBNSl1dNfKgsmV5h/CLM/Uf/N1oWT6Dyg2EE+HLgP8AL6Stvxn436hb48oM+c8gnGdZQOiuus7MRmbZ1m8IB5U5hP18knCCHgAz+2+0fmdSDuZmNoFwfuMeQv//bMKBrTKGAGdFXT4Vig7kq1K6dzZL1FqaBEw3s3XR4veBL83s2yjNBuBYQtffF4TW3UOEVmwmA4G5URfXYEIrE+AOwgUAi4APgNcz5H2M0K0aq5WSsh/vARuBSWY2N2X5i8DfgKej+kwjXHEHsC3hh8ESQpfn98BtldmuS4bM/CFdziVF0pOE8wcvxUj7PPAPM3u1orR1kaTGhCvqukfn8CqTdxTwpJk9lJfKubzxoOKcywtJlwPHmNlhlczXExgJtE27QMPVAXnt/pLUT9JMhRvWrs6wfjtJr0j6WNJ0SefGzeucq70kzQV+SzihX5l8jwBvAr/zgFI35a2lEl1a+BnhssAiwl3dZ5jZJylp/gfYzsz+IKk14b6AnQgn2HLmdc45V/vks6XSi3Dj15zopOHTwPFpaQxoHp3YbEa47rw4Zl7nnHO1TD4H+9uFsjeLFREubUx1DzCc0uEkTjOzjdFlpxXlBUDSIMJ9FjRt2rTHXnvtlUztnXNuCzBx4sRFZtY6qfLyGVQyXVaZ3tf2c8J4SocBuwMjJb0bM29YaDaEcCknhYWFNmHChKrW1znntjiSvqw4VXz57P4qIgxzUWJXQosk1bmEwfTMzGYTrpvfK2Ze55xztUw+g8p4oJOkjpIaAqcTurpSzQMOB5DUhnA38JyYeZ1zztUyeev+MrNiSZcAbxBGGx1mZtMlDY7WPwD8GXhY4TkQAv5gZosAMuXNV12dc84lo17d/OjnVFx9t379eoqKilizZk1NV8XVMY0aNWLXXXelQYMGZZZLmmhmcZ8+WqF696hX5+qzoqIimjdvTocOHYg5xJhzmBnff/89RUVFdOzYMa/b8gElnatD1qxZQ8uWLT2guEqRRMuWLaulhetBxbk6xgOKq4rq+rvxoOKccy4xHlScc5VSUFBAt27d6NKlC6eccgo//PBDxZkiDz/8MJdcckmlttesWeYnV1977bW8+eabAPTp04eSi3SOOuooli5dytKlS7nvvvsqta1crrrqKjp37sxVV1XugZITJkzg0ksvrfJ2U/ezLvAT9c65SmncuDGTJ08G4KyzzuKBBx7g8ssv37R+w4YNFBQUZMmdnBtvvDHj8ldfDY+nmTt3Lvfddx8XX3xxItt78MEH+e6779hmm20qla+wsJDCwqpfXJVtP2srb6k456rs4IMPZvbs2YwePZq+ffty5pln0rVrV9asWcO5555L165d2W+//Xj77bc35Zk/fz79+vVjzz335IYbbti0/IQTTqBHjx507tyZIUOGlNnOFVdcQffu3Tn88MP57rvvADjnnHN47rnnytWpQ4cOLFq0iKuvvprPP/+cbt26cdVVVzFw4EBefvnlTenOOusshg8ve0+1mXHVVVfRpUsXunbtyjPPPAPAcccdx6pVq9h///03LSvRtWtXli5dipnRsmVLHn30UQAGDhzIm2++yejRoznmmGMAuP766znvvPPo06cPu+22G3fddRcQAuDee+/NBRdcQOfOnTnyyCNZvXp1uf3s0KED1113Hd27d6dr167MmDEDgO+++44jjjiC7t27c+GFF9K+fXsWLVoU6ztMmrdUnKvDbr/99nLLevToQZ8+fVi3bh133313ufUHHnggP/3pT1m5ciUPPvhgmXVXXBH/8SfFxcW89tpr9OvXD4APP/yQadOm0bFjx031mjp1KjNmzODII4/ks88+K5OuSZMm9OzZk6OPPprCwkKGDRvGDjvswOrVq+nZsyf9+/enZcuWrFq1iu7du3P77bdz4403csMNN3DPPfdUWL+//vWvTJs2bVOr6p133uHvf/87xx9/PMuWLWPcuHE88sgjZfK88MILTJ48mY8//phFixbRs2dPDjnkEIYPH06zZs02lZXqoIMO4r333qN9+/bstttuvPvuu5x99tl88MEH3H///aTfOzdjxgzefvttVqxYwZ577slFF10EwKxZs3jqqacYOnQop556Ks8//zwDBgwot71WrVoxadIk7rvvPm677TYeeughbrjhBg477DCuueYaXn/99XJBuTp5S8U5VymrV6+mW7duFBYW0q5dO84//3wAevXqtekeiLFjxzJw4EAA9tprL9q3b78pqBxxxBG0bNmSxo0bc9JJJzF27FgA7rrrLvbdd18OOOAA5s+fz6xZ4QnEW221FaeddhoAAwYM2JS+sg499FBmz57Nt99+y1NPPUX//v3Zeuuyv6vHjh3LGWecQUFBAW3atOHQQw9l/PjxOcs9+OCDGTNmDGPGjOGiiy5i6tSpfPXVV+ywww4ZzwcdffTRbLPNNrRq1Yodd9yRhQsXAtCxY0e6desGhB8Gc+fOzbi9k046qVyasWPHcvrppwPQr18/WrRoEfdjSZy3VJyrw3K1LBo2bJhzfbNmzSrVMimRek4lVdOmTTdN5xqpI/3SVkmMHj2aN998k/fff58mTZrQp0+frPdUbM6lsQMHDuSJJ57g6aefZtiwYeXWV2WEkUMOOYR7772XefPmcdNNN/Hiiy/y3HPPcfDBB2dMn3pOpqCggOLi4ozLS7q/suVPzVubRkbxlopzLnGHHHIITzzxBACfffYZ8+bNY8899wRg5MiRLF68mNWrV/PSSy9x0EEHsWzZMlq0aEGTJk2YMWMGH3zwwaayNm7cuOmcwpNPPknv3r1j1aF58+asWFH2icTnnHMOd9xxBwCdO3fOWO9nnnmGDRs28N133zFmzBh69eqVcztt27Zl0aJFzJo1i912243evXtz2223ZQ0q+dC7d2+effZZAEaMGMGSJUuqbdvpPKg45xJ38cUXs2HDBrp27cppp53Gww8/vOkXdu/evRk4cCDdunWjf//+FBYW0q9fP4qLi9lnn33405/+xAEHHLCprKZNmzJ9+nR69OjBqFGjuPbaa2PVoWXLlhx00EF06dJl02XAbdq0Ye+99+bcc8/NmOfEE09kn332Yd999+Wwww7jlltuYaeddqpwW/vvvz977LEHELrDvvrqq9jBLwnXXXcdI0aMoHv37rz22mv86Ec/onnz5tW2/VQ+oKRzdcinn37K3nvvXdPVqLN++OEHunbtyqRJk9huu+1qujqJWbt2LQUFBWy99da8//77XHTRRRm7KDP9/fiAks45VwVvvvkm5513Hpdffnm9CigA8+bN49RTT2Xjxo00bNiQoUOH1lhdPKg457YIP/vZz5g3b15NVyMvOnXqxEcffVTT1QD8nIpzdU596rJ21ae6/m48qDhXhzRq1Ijvv//eA4urlJLnqTRq1Cjv2/LuL+fqkF133ZWioqJNQ5U4F1fJkx/zLa9BRVI/4E7Cc+YfMrO/pq2/CjgrpS57A63NbLGkucAKYANQnOTVCc7VVQ0aNMj7k/uc2xx5CyqSCoB7gSOAImC8pOFm9klJGjO7Fbg1Sn8scJmZLU4ppq+Z1cyoaM455yotn+dUegGzzWyOma0DngaOz5H+DOCpPNbHOedcnuUzqOwCzE+ZL4qWlSOpCdAPeD5lsQEjJE2UNChvtXTOOZeYfJ5TyTTqW7ZLVo4F3kvr+jrIzBZI2hEYKWmGmY0pt5EQcAYBtGvXbnPr7JxzbjPks6VSBLRNmd8VWJAl7emkdX2Z2YLo/VvgRUJ3WjlmNsTMCs2ssHXr1ptdaeecc1WXz6AyHugkqaOkhoTAMTw9kaTtgEOBl1OWNZXUvGQaOBKYlse6OuecS0Deur/MrFjSJcAbhEuKh5nZdEmDo/UPRElPBEaY2aqU7G2AF6PnJmwNPGlmr+errs4555LhoxQ759wWLOlRin2YFuecc4nxoOKccy4xHlScc84lxoOKc865xHhQcc45lxgPKs455xLjQcU551xiPKg455xLjAcV55xzifGg4pxzLjEVBhVJbST9Q9Jr0fxPJJ2f/6o555yra+K0VB4mDAq5czT/GfC7PNXHOedcHRYnqLQys2eBjRBGHwY25LVWzjnn6qQ4QWWVpJZET22UdACwLK+1cs45VyfFeZ7K5YSHa+0u6T2gNXByXmvlnHOuTqowqJjZJEmHAnsSnjs/08zW571mzjnn6pwKg4qkk9IW7SFpGTA1en68c845B8Tr/jofOBB4O5rvA3xACC43mtljeaqbc865OiZOUNkI7G1mCyHctwLcD+wPjAE8qDjnnAPiXf3VoSSgRL4F9jCzxUDOcyuS+kmaKWm2pKszrL9K0uToNU3SBkk7xMnrnHOu9onTUnlX0r+Bf0Xz/YExkpoCS7NlklQA3AscARQB4yUNN7NPStKY2a3ArVH6Y4HLzGxxnLzOOedqnzhB5deEQHIQ4eqvR4HnzcyAvjny9QJmm9kcAElPA8cD2QLDGcBTVczrnHOuFohzSbEBz0WvytgFmJ8yX0Q4D1OOpCZAP+CSKuQdBAwCaNeuXSWr6JxzLklxBpQ8QNJ4SSslrYvOeyyPUbYyLLMsaY8F3ovO01Qqr5kNMbNCMyts3bp1jGo555zLlzgn6u8hdE3NAhoDvwLujpGvCGibMr8rsCBL2tMp7fqqbF7nnHO1RKznqZjZbKDAzDaY2T/JfS6lxHigk6SOkhoSAsfw9ESStgMOBV6ubF7nnHO1S5wT9T9EB/bJkm4BvgaaVpTJzIolXUIYNr8AGGZm0yUNjtY/ECU9ERhhZqsqyluZHXPOOVf9FM7D50ggtQcWAg2By4DtgHvN7PP8V69yCgsLbcKECWUXzp8PbdtmzuCcc1s4SRPNrDCp8uJ0f51gZmvMbLmZ3WBmlwPHJFWBvHruOWjXDt58s6Zr4pxzW4Q4QeWXGZadk3A98uO998L71Kk1Ww/nnNtCZD2nIukM4Eygo6TUk+TNge/zXbFEFBeH94KCmq2Hc85tIXKdqB9HOCnfCrg9ZfkKYEo+K5WYkqCydZzrEZxzzm2urEdbM/sS+JIw7H3dtGFDePeWinPOVYs4d9SfJGmWpGWSlktaEfOO+prnLRXnnKtWcY62twDHmtmn+a5M4kpaKh5UnHOuWsS5+mthnQwo4CfqnXOumsX5CT9B0jPAS8DakoVm9kK+KpWYkpaKMo1P6ZxzLmlxgsq2wA/AkSnLDKj9QaWkpVISXJxzzuVVnOepnFsdFcmLkqBS8u6ccy6vKgwqkvYA7gfamFkXSfsAx5nZX/Jeu0qySZP414UXctyjj/LIlVdy+GefsTtAcTHr1q3j7rvLj9h/4IEH8tOf/pSVK1fy4IMPllt/6KGHUlhYyJIlSxg2bFi59UcccQT77LMPCxcu5PHHHy+3/qijjmLvvfdm/vz5PPvss+XWn3DCCey+++58/vnnvPTSS+XWn3rqqbRt25ZPP/2UV199tdz6AQMG0KZNG6ZMmcLIkSPLrT/vvPNo0aIFEyZM4J133im3/sILL6RZs2aMGzeO999/v9z63/zmNzRs2JDRo0czceLEcuuvuOIKAEaMGMHUtJELGjRowKWXXgrAf/7zH2bMmFFmfdOmTRk8eDAAL774InPmzCmzvkWLFpx33nkAPPvss8yfP7/M+jZt2jBgwAAAHn/8cRYuXFhmfdu2bTn11FMBGDZsGEuWLCmzfrfdduPEE08E4IEHHmDVqlVl1u+1114cffTRANx1112sX7++zPquXbty5JGhAX/77beTrkePHvTp08f/9vxvr1b/7SUtzon6ocA1wHoAM5tCGIq+1pEZ/YcOZZs1a+g0ZQoqGSzTWyrOOVct4oxSPN7Mekr6yMz2i5ZNNrNu1VHByiiUbIIEZnDHHfDaa/DGG3DnnRD9anHOubyZORPefx8+/DCMjr7ffvDzn4djUsmxaautwvunn4bpdu1g40Zo0gQefRSuuw5Wr4bOneGvf4WOHeFf/4LBg/NyJWvSoxTHOVG/SNLuRI/zlXQyYfiW2qkkSJZ8geAtFee2dKtWhQN98+Zw111w7rlw330waBAcemg4uK9ZE9ItXAh77w3ffQcjR8LatfDWW7ByJbz7LqxYEcps2xY6dYJRo/JT59Gj4YADSucvuaR0etttQ7A6+2w4+ujS493TT8Mrr8D06bBoEWy3HfzoR7DzzrD77vDf/8Lrr5eWs99+iVc7TlD5NTAE2EvSV8AXwIDEa5IPGzeGdw8qztUNK1ZA06bhIL85ioth4kQ46ihYvLj8+sceC+/PPVf1bcyfH16p+vSBffYJrZDdd4cXXoBvv82cv3VraNUqpIXQCvnHP+C008LjOh55JNSvb9+wDx9/XJp3+fLQevnXv3LXccGC0vIz+eijCnezsuJc/TUH+JmkpsBWZrYi8Vrkg+RBxbmaUNLVEzftN9/A5Mlw773wn/+UrnvySdhtN3j7bXjoIZg7F554Ikz37x+6gyC0ML76KrQivvgidH0vzzKS1LnnwpFHws03h1ZLyeMxSrRoASUn1bffHrp2hSOOCF1U7dvDwQeH1sz69aFenTuH4LHVVplH7rj//nifQ7pjjgmvdMXF4bP98MPwObz1Fnz5ZVi37bZw3HGh9dW1a2hhzZwJ48eH9YsXh/3/8Y/D/Nq14bEgPXtWrY7ZmFnOF/D/gO1T5lsAf6koX028eoQ/0fC6+26zPn3C9A03mHMuxcaNZg89ZHbooaX/M+PGhXUrV5qtXm324Ydm69ebjRplJpn171+aNvW1yy5mTZqUX3788WaHHx6mu3c3u/12s0WLzA47LHM5lX01bpx7/UknmS1ZUoMfct0ATLAEj8NxgspHGZZNilU49ANmArOBq7Ok6QNMBqYD76QsnwtMjdbF2ukyQeXee8169w7Tf/pTYl+AmZlNn262fHmyZdY1zz9v9sgjNV2LLdeSJWZ/+EPp3/vNN+dOv3Gj2axZIXCcfLLZ2Wcnc2BP4lVQYNaqVfibMjNbt85s7FizFi3MttrKbMAAs6efNvvLX8xatzb74x8zl3PAAWG/li7N+8dfnyQdVOKcUymQtI2ZrQWQ1BjYpqJMkgqAe4EjgCJgvKThZvZJSprtgfuAfmY2T9KOacX0NbNFMeqYqQKwbl2YXrwYiopg112rVFQZixaFJi+EP+WbbgrzJ5yw+WXXBV9/Db/7HZTc99CkCZx8co1WKXGffBL+VrbdtnL5JkwI3SCtWoXumPbtw5WH11wTulKSsGEDXHYZpN97cM01MGtWuFqodevwt7l0Kfz2t6XnD9L16weNG8Oxx4YTuv37Z063115he8OGwZlnQo8eodvl0ENDd9N//wsHHhi2u8ceIc+YMaGrap99YJddYMoUuP32cM5k7Fho1gwaNYI99yy/vQYN4KCDMp8L+eMfw/vAgeF8wPHHh/Q+aGztUVHUAX4PjAXOB86Lpn8fI9+BwBsp89cA16SluZgsXWmElkqrykTIMi2V++8369at7C+ZJOy+e2l5//M/pdN33lma5p//NJsxI5nt1RYrV5p17pz5F+I559R07ZIxeXL5ffvLX8w++CB0AX37bfk8GzaElnDcX+WPP27273/H75Z5441Qr9WrzYqLzc48s7Ss444Ly6dMib/9Ro3C+403htaL2+JRnd1fgIC2hG6s2whPgPx5rILhZOChlPmBwD1pae4gtGZGAxOBs1PWfQFMipYPyrGdQcAEYEKZoPLAA2Y/+UnyQSXXP+zGjaFbDMx23jmZ7dW09evNvvuu/L6ef37Z+euuq746/fvfZbc9e3bFeSZNMhs/vuyya6+NfzDOx2vkyPDZ3nyz2Wefla3bunVmhxySOV+/fmarVpVN/9VXmdPefXcIZOnlOxep1qAStsfEKhUMp2QIKnenpbkH+ABoSnhs8Sxgj2jdztH7jsDHwCEVbbNcUPnxj5MNKqtX5z5ILF8e+qzBbOut6/Z5l40bwyt9H3v3Dr/OzcJ7YWHpunnzzIYONfv9780GDy5Nl15u3O2bhVbf/Plhes2a7J/9N99kLufmm8un3bjR7NhjM5dz881h/bRpZsccEy84bLNNaE1ks2GDWVFRaNHlKmfgQLOGDc1++cvSZVLZNFdc4S0Ml6iaCCr3Aj0rXXC87q+rgetT5v8BnJKhrOuBKyvaZpmg8uCDZu3aJRtURo3KfVAYP778svXrM5e1dm04qfjll5tfr2eeMZs5MwS9lSvN3nvP7NZbq17etGmZ92/16vJpN24023ff7J/JwoVmp51mduKJ5VuOUHqyf8MGswsvDPWeOLHiA/kdd4TWR8l8Sctw9WqzuXPNXnvN7LLLKi7nhRdCt9LGjaF1kElJcFy92mzFCrNPPzX77W9Dt1i27zebDRvC6+qrS+vQtm3mut10U9nPec2aym3LuRhqIqh8AmwAPgemEK7ImhIj39bAHKAj0DBqbXROS7M38FaUtgkwDegStVyaR2maAuMIJ/PjB5WhQ8122qnyQWXuXLMRIzKve/75UE6nTmZ//nNpua+9lv2gtf32mX9ZPvpoMsGuQYOy22vevHT68ssztxZyWb++/D4sWWK2eHH2PJlaNPl8pR5c166tOP3FF5ftmgSzHXYIwaS2uOeecFXWgAHhPE5lg5VzVVQTQaV9pleswuEo4LMoIP0xWjYYGJyS5qoocE0Dfhct2y0KQh8TLjX+Y5ztlQkqDz0UDhypy+IcYEtOZGbyf/8X1n3zTThItWkTlqX+Ys70ev31suXMn192/YoVFdcrky+/rPiAOnRo/O6SOXPK5584MV7ee+8169LF7IcfwvyyZZnr88Yb4b1zZ7O+fbPXu3PncCL6uefMttvOrGnTcODN9Gs914nqWbPKpl20yOyllyofbJ2rp6o9qIRt0hs4N5puDXRMshJJvcoElX/8w6xZs7IHmIq6D4qLcwegSy8NZaYfpNNPkr7yStk++8ceC+W9/nrIe9ddZdNPnRrK+eYbs1NOCQe+XNauLXsuo6LXLbfkLu+//zX7/vuyeZL4FV9cbPazn5k9+2y4d6AkeKZ+tvPmhSvnRo4M8xs3Vu2Av3FjCC5+vsG5SqmJlsp1wCvAZ9H8zsB7SVYiqVeZoDJsWDjpmbqsohPnCxfmTnvssWZdu2bOm3p58RdfhGWpVyndckt4v/pqs1NPDdOnnx7ejz46HAwHDw7zAwfmrmd60DAL5X/yiW3q2tmwoWya1IPtypVmb70Vpt9+u3x5VW05OefqnJoIKpMJlxZ/lLKswnMqNfEqF1TSD5Z33ZX70506tTTt0KGhNZDat921awgsmaSeVyj5pZ3rXEO7dmUv1e3SJVwqmhooMikJHCWvIUOyp73++rLlm4WrkHK1ar7+Ovdn5JyrV5IOKnGGAl0XbdgAooEla79Mz6XP9kyVDRvCHfhdu5Yuu+CCcId0gwZhALrFi8Pga23bZi4jdQC9khFWcw2qN29euPO6xLRppUNSd+tWujwE8VLR0+i44oqw7oILsm/juuvCIHwl5Q8ZEu7yTte2LTz/fHiGw047ZS/POecqECeoPCvpQWB7SRcAbxKeBlm7pT1+M6e33869fuRIaNkyTH/2WfZ0M2eGV6oMj2EFwlAgUDoiaqrJk8MDxi66KAQoKWxXgg8+CGluvTV3nUu0bx8e8gNw4YWlwbZv39I0X3wBJ50Uhs1wzrnNEGfo+9skHQEsB/YArjWz8g+krm0qE1Sefz5+2kMOyb6uZNyjVPvuWzr997+H5yTceWcYLhvCmFCrV4eHAs2dG1pFI0aE50CkSh0j6eKL4w8tDvD552WfTzFuXBiryTnnElbh44QBJO0E9CJ0gY03s2/yXbGqKJRsQsnMLbfA739fPlH6/n79dXgqWlzFxZV/pOedd0KvXvEP5LkCxnHHwcsvV277EAYVPPvsEBTfeafy+Z1z9VLSjxOusPtL0q+AD4GTCON5fSDpvKQqkDclI5ymH6DHjCmdvvLK8BCguHr1qtozon/726q3DL77DoZGvY2vvlq1gAJhJNpTTgnPwHbOuTypsKUiaSbwUzP7PppvCYwzswxjVtesMi2VEgUFZU/a/+1vpS2YbC2CnXcOj+FM9/Ofl32+c75MmxYuGhg1quy5D+ecS1i1t1QIz0JJfYTwCmB+lrS1T/qzrnM9d+GKK2D06PDsh+uvD8+C+O1vS9c3a5aPGpbXpUvopvOA4pyrY+K0VB4FugIvE86pHE/oDvsMwMz+L891jC1jS6Vhw9KHdQH84helV2SltlQ6doQ5czIX/NBD4dLdc86Bf/4zwRo751zNSrqlEudxaZ9HrxIlnfrNk6pEXqUHzddeC+9r15ZdPnly9jKaR7ua6eou55xzm8S5pPiG6qhI3mRriS1bVnY+16NjTzklXPV12mnJ1cs55+qhOOdU6raSoJLeyhg/Pn4ZW20FZ53lz8F2zrkK1P+gsnFjeE+9Z6V3bzjmmJqpj3PO1WP1P6iUtFQaNy5d9t57ZdOkrnPOOVdlFfbnSGoNXAB0SE1vZrX/BshUzXNcV5B62bBzzrkqi3OS4GXgXcJAkhmG/q0jcgWV7t2rrx7OOVePxQkqTczsD3mvSb786U9hHLDU4eRTffxxuMnROefcZotzTuXfko6qOFl5kvpJmilptqSrs6TpI2mypOmS3qlM3lguuwzWrAmjAWfiAcU55xIT5476FUBTYC2wnvAUSDOzHDd2gKQCwl33RxCGehkPnGFmn6Sk2R4YB/Qzs3mSdjSzb+PkzSTjHfWp+9e6NSxalH29c85tYap97C8za25mW5lZYzPbNprPGVAivYDZZjbHzNYBTxOGeEl1JvCCmc2LtvVtJfJW3t13b3YRzjnnsst6TkXSXmY2Q1LGs9hmNqmCsneh7MCTRcD+aWn2ABpIGk0Y9uVOM3s0Zt6Seg4CBgH0qKBCHHts2fknnqgoh3POuUrIdaL+csLB+vYM6ww4rIKyM40rn97XtDUhFhwONAbel/RBzLxhodkQYAiE7q8yK6+7rmzipk1h+PDwoCuAM8/MvQfOOecqJWtQMbNB0XtVx18vAtqmzO8KpD+kpAhYZGargFWSxgD7xsxbsUzDqqS3VpxzziUmzpMf35V0U3Q1VmVGJh4PdJLUUVJD4HRgeFqal4GDJW0tqQmhi+vTmHkr5mN1OedctYpz1P0l0BvoD9wqaS3wrpldliuTmRVLugR4AygAhpnZdEmDo/UPmNmnkl4HpgAbgYfMbBpApryV3zsPKs45V53iDH0/R9JqYF306gvsHadwM3sVeDVt2QNp87cCt8bJW2keVJxzrlrF6f76HHgJaAP8A+hiZv3yXK9keFBxzrlqFeeoexeh++sMYD/gHUljzOzz3NlqgWxB5c03cz+UyznnXJXE6f66E7hTUjPgXOB6wtVYBfmtWgKyBZXDD6/eejjn3BYiztD3txNaKs2A94FrCaMW137e/eWcc9UqzlH3A+AWM1uY78okzoOKc85VqzjdX/+qjorkhQcV55yrVvX7ccIeVJxzrlp5UHHOOZeYWEFFUm9J50bTrSV1zG+1EuJBxTnnqlWcmx+vA/4AXBMtagA8ns9KJWar+t0Qc8652ibOUfdE4DhgFYCZLSA8+6T286DinHPVKs5Rd52FZw4bgKSm+a1SgpTpsSzOOefyJU5QeVbSg8D2ki4A3gSG5rdaCfGWinPOVas496ncJukIYDmwJ3CtmY3Me82S4C0V55yrVnGGabkM+FedCSSpvKXinHPVKs5Rd1vgjegJkL+W1CbflUqMt1Scc65aVRhUzOwGM+sM/BrYmTD0/Zt5r1kSvKXinHPVqjJH3W+Bb4DvgR3zU52EeUvFOeeqVZybHy+SNBp4C2gFXGBm+8QpXFI/STMlzZZ0dYb1fSQtkzQ5el2bsm6upKnR8gnxdymFt1Scc65axRnHpD3wOzObXJmCJRUA9wJHAEXAeEnDzeyTtKTvmtkxWYrpa2aLKrPdtEpUOatzzrnKyxpUJG1rZsuBW6L5HVLXm9niCsruBcw2szlR/qeB44H0oJI/nTpV26acc87l7v56MnqfCEyI3iemzFdkF2B+ynxRtCzdgZI+lvSapM4pyw0YIWmipEHZNiJpkKQJZbrIBgwAM9ixbpz6cc65+iJrS6WkS8rMqjoicaa+J0ubnwS0N7OVko4CXgJKmhcHmdkCSTsCIyXNMLMxGeo5BBgCUCiF8r3byznnakScE/VvxVmWQRHQNmV+V2BBagIzW25mK6PpV4EGklpF8wui92+BFwndafH4CXrnnKsRWY++khpF51FaSWohaYfo1YFwv0pFxgOdJHWU1BA4HRieto2dpNCskNQrqs/3kppKah4tbwocCUyLv1ceVJxzribkuvrrQuB3hAAykdLurOWEq7pyMrNiSZcAbwAFwDAzmy5pcLT+AeBk4CJJxcBq4HQzs+iu/RejeLM18KSZvR57rzyoOOdcjVAY1T5HAuk3ZnZ3NdVnsxRKNgHgV7+CoXVjIGXnnKtJkiaaWWFi5VUUVKKNdgF+AjQqWWZmjyZViaRsCioQrv5yzjmXU9JBJc4oxdcBfQhB5VXgF8BYoNYFFeecczUrzsmHk4HDgW/M7FxgX2CbvNbKOedcnRQnqKw2s41AsaRtCQNL7pbfajnnnKuL4oz9NUHS9oRHCE8EVgIf5rNSzjnn6qZYJ+o3JQ73qGxrZlPyVqPN4CfqnXOucqrtRL2k7rnWmdmkpCrhnHOufsjV/XV7jnUGHJZwXZxzztVxuQaU7FudFXHOOVf3xblP5exMy2vjzY/OOedqVpyrv3qmTDci3LMyCb/50TnnXJoKg4qZ/SZ1XtJ2wGN5q5Fzzrk6qyrD+f5A6YO0nHPOuU3inFN5hdInNm5FGAPs2XxWyjnnXN0U55zKbSnTxcCXZlaUp/o455yrw+KcU3kHIBr3a+toegczW5znujnnnKtj4nR/DQL+THgy40bCEyANH1TSOedcmjjdX1cBnc1sUb4r45xzrm6Lc/XX54QrvipNUj9JMyXNlnR1hvV9JC2TNDl6XRs3r3POudonTkvlGmCcpP8Ca0sWmtmluTJJKgDuBY4AioDxkoab2SdpSd81s2OqmNc551wtEieoPAiMAqYSzqnE1QuYbWZzACQ9DRwPxAkMm5PXOedcDYkTVIrN7PIqlL0LMD9lvgjYP0O6AyV9DCwArjSz6ZXIW3IhwSCAHptKPLAK1XXOObe54pxTeVvSIEk/krRDyStGPmVYlv7krElAezPbF7gbeKkSecNCsyFmVljmITO33BKjes4555IWp6VyZvR+TcqyOJcUFwFtU+Z3JbRGSgsxW54y/aqk+yS1ipM3p4KC2Emdc84lJ87Njx2rWPZ4oJOkjsBXwOmUBigAJO0ELDQzk9SL0HL6HlhaUV7nnHO1T96ep2JmxZIuAd4ACoBhZjZd0uBo/QPAycBFkooJN1eebmYGZMxbif1yzjlXAxSO4TkSSHenzG56noqZnZzPilVFoWQTAMaN85P1zjkXg6SJZc5Jb6b6+TwVZTrP75xzLt/8eSrOOecSUz+fp1JBl55zzrn88OepOOecS0zWoCLpx0CbkueppCw/WNI2ZvZ53mtXVX5OxTnnakSucyp3ACsyLF8drXPOOefKyBVUOpjZlPSFZjYB6JC3GjnnnKuzcgWVRjnWNU66Is455+q+XEFlvKQL0hdKOh+YmL8qOeecq6tyXf31O+BFSWdRGkQKgYbAiXmul3POuTooa1Axs4XATyX1BbpEi/9jZqOqpWbOOefqnDjDtLwNvF0NdXHOOVfHVWWYltrP71NxzrkaUT+DinPOuRpRP4OKj/3lnHM1on4GFeecczWifgYVP6finHM1on4GFeecczUir0FFUj9JMyXNlnR1jnQ9JW2QdHLKsrmSpkqaLGlCPuvpnHMuGXGep1IlkgqAe4EjgCLCsC/DzeyTDOn+BryRoZi+Zrao0hvv0qXiNM455xKXz5ZKL2C2mc0xs3XA08DxGdL9Bnge+DaxLTdtmlhRzjnn4stnUNkFmJ8yXxQt20TSLoRxxB7IkN+AEZImShqUt1o655xLTN66v4BMl2Cl30ByB/AHM9ug8ldsHWRmCyTtCIyUNMPMxpTbSAg4gwB6bH6dnXPObYZ8BpUioG3K/K7AgrQ0hcDTUUBpBRwlqdjMXjKzBQBm9q2kFwndaeWCipkNAYYAFEp+16NzztWgfHZ/jQc6SeooqSFwOjA8NYGZdTSzDmbWAXgOuNjMXpLUVFJzAElNgSOBaXmsq3POuQTkraViZsWSLiFc1VUADDOz6ZIGR+sznUcp0YbwLJeSOj5pZq/nq67OOeeSIatH42QVSjYBfOwv55yLSdJEMytMqjy/o94551xiPKg455xLjAcV55xzifGg4pxzLjEeVJxzziXGg4pzzrnEeFBxzjmXGA8qzjnnEuNBxTnnXGI8qDjnnEuMBxXnnHOJ8aDinHMuMR5UnHPOJcaDinPOucR4UHHOOZcYDyrOOecS40HFOedcYjyoOOecS4wHFeecc4nJa1CR1E/STEmzJV2dI11PSRsknVzZvM4552qPvAUVSQXAvcAvgJ8AZ0j6SZZ0fwPeqGxe55xztUs+Wyq9gNlmNsfM1gFPA8dnSPcb4Hng2yrkdc45V4tsnceydwHmp8wXAfunJpC0C3AicBjQszJ5U8oYBAyKZtcKpiFtXs1rr1bAopquRB7V5/2rz/sGvn912Z5JFpbPoJLpyG5p83cAfzCzDSobCOLkDQvNhgBDACRNMLPCyle1bvD9q7vq876B719dJmlCkuXlM6gUAW1T5ncFFqSlKQSejgJKK+AoScUx8zrnnKtl8hlUxgOdJHUEvgJOB85MTWBmHUumJT0M/NvMXpK0dUV5nXPO1T55CypmVizpEsJVXQXAMDObLmlwtP6ByuaNsdkhCVS9NvP9q7vq876B719dlui+ySzjqQrnnHOu0vyOeuecc4nxoOKccy4x9Sao1IdhXSTNlTRV0uSSy/wk7SBppKRZ0XuLlPTXRPs7U9LPa67mmUkaJulbSdNSllV6fyT1iD6X2ZLukmrHjUhZ9u96SV9F3+FkSUelrKsz+yepraS3JX0qabqk30bL68X3l2P/6vz3J6mRpA8lfRzt2w3R8ur57syszr8IJ/M/B3YDGgIfAz+p6XpVYT/mAq3Slt0CXB1NXw38LZr+SbSf2wAdo/0vqOl9SKv7IUB3YNrm7A/wIXAg4f6l14Bf1PS+5di/64ErM6StU/sH/AjoHk03Bz6L9qFefH859q/Of39RPZpF0w2A/wIHVNd3V19aKvV5WJfjgUei6UeAE1KWP21ma83sC2A24XOoNcxsDLA4bXGl9kfSj4Btzex9C3/lj6bkqVFZ9i+bOrV/Zva1mU2KplcAnxJGuqgX31+O/cumzuyfBSuj2QbRy6im766+BJVMw7rk+gOprQwYIWmiwvAzAG3M7GsI/wjAjtHyurrPld2fXaLp9OW12SWSpkTdYyVdDHV2/yR1APYj/OKtd99f2v5BPfj+JBVImkwYU3GkmVXbd1dfgkrsYV1quYPMrDthdOZfSzokR9r6ss8lsu1PXdvP+4HdgW7A18Dt0fI6uX+SmhEGfP2dmS3PlTTDsrq4f/Xi+zOzDWbWjTAaSS9JXXIkT3Tf6ktQqRfDupjZguj9W+BFQnfWwqgZSvReMppzXd3nyu5PUTSdvrxWMrOF0T/0RmAopV2SdW7/JDUgHHCfMLMXosX15vvLtH/16fsDMLOlwGigH9X03dWXoLJpSBhJDQnDugyv4TpViqSmkpqXTANHAtMI+/HLKNkvgZej6eHA6ZK2URjOphPhpFptV6n9iZrpKyQdEF15cnZKnlqn5J82ciLhO4Q6tn9RXf4BfGpm/5eyql58f9n2rz58f5JaS9o+mm4M/AyYQXV9dzV5lUKSL+AowhUcnwN/rOn6VKH+uxGuwPgYmF6yD0BL4C1gVvS+Q0qeP0b7O5NacEVNhn16itCFsJ7wq+f8quwPYeDRadG6e4hGgqjpV5b9ewyYCkyJ/ll/VBf3D+hN6OqYAkyOXkfVl+8vx/7V+e8P2Af4KNqHacC10fJq+e58mBbnnHOJqS/dX84552oBDyrOOecS40HFOedcYjyoOOecS4wHFeecc4nxoOK2OJJGSyqshu1cGo2C+0S+t5Vl+x2UMoLyZpTzsKSTk6hTlvITqaerHfL5jHrn6h1JW5tZcczkFxOu+f8iofKcq/W8peI2W/RL81NJQ6PnN4yI7uTNS5mpLQ1JrSTNjabPkfSSpFckfSHpEkmXS/pI0geSdkjZxABJ4yRNk9Qryt80GkRwfJTn+JRy/yXpFWBEhrpeHpUzTdLvomUPEG5oHS7psrT0ZcpTeM7FS9Eghh9I2idKd72kK1PyTYs+l1yfTQ+F52i8D/w6JW9nhWdsTI620ynDfqyUdLukSZLektQ6Q5pro89nmqQhCnaXNCklTSdJE1Pq847CIKlvqHSYkIz1dHWfBxWXlE7AvWbWGVgK9E9PIOkslT78KPX1XFXLzKALcCZhzKabgB/MbD/gfcIwEyWamtlPCa2JYdGyPwKjzKwn0Be4VWHIHAjPlPilmR2Wtk89gHOB/QnPrLhA0n5mNpgwTlJfM/t7hnqmlncD8JGZ7QP8D2GI8Ypk+2z+CVxqZgempR8M3GlhkMFCyo4+W6IpMMnCoKbvANdlSHOPmfU0sy5AY+AYM/scWCapW5TmXOBhhbG17gZONrMehM/5pgrq6eo47/5ySfnCzCZH0xOBDukJzOwJoDLnFyosM4O3LTwfY4WkZcAr0fKphOErSjwV1WmMpG0Vxko6EjgupXXQCGgXTY80s0zPTukNvGhmqwAkvQAcTBgmI5fU8noTBQUzGyWppaTtKshf7rOJ8mxvZu9Eyx8jjHgNIaj+UdKuwAtmNitDmRuBZ6Lpx4EXMqTpK+n3QBNgB8KQQq8ADwHnSrocOI0Q1PckBPmRCg8MLAC+rqCero7zoOKSsjZlegPhV2wZks4CrsqQd7aZZToRnK3MYkpb2Y1y5NmYMr+Rsn/v6eMTlQz13d/MZqbVe39gVYb6QebhweNILS/bEOOp+wll9zXTZyOyDE1uZk9K+i9wNPCGpF+Z2agK6limLEmNgPuAQjObL+n6lDo9T2jZjAImmtn3knYGpqe3RqIA7uND1VPe/eWqjZk9YWbdMrwqe2XRXKBHNF3Vq5JOA5DUG1hmZsuAN4DfKPpZLWm/GOWMAU6Q1CTqKjsReLeSdRkDnBVtsw+wyMKzPeYSHleMpO6ER71mZWGY82XRPlFSZpR/N2COmd1FGChxn/IlsBWln+eZwNi09SUBZJHCc0g2ffZmtobw+d1P6NqCMDhha0kHRnVoIKlzrnq6us9bKq4uug14VtJAwi/jqlgiaRywLXBetOzPwB3AlCiwzAWOyVWImU2S9DCljx14yMwq6vpKdz3wT0lTgB8oHZ78eeBshSf4jSeMwl2Rc4Fhkn4gHORLnEa4OGE98A1wY4a8q4DO0Un2ZVGeTcxsqaShhK7EuVGdUj0BnER0MYOZrVO4FPmuqMtra8LnOz1HPV0d56MUO+eAcPWXmTXbjPxXAtuZ2Z8SrJarY7yl4pzbbJJeJDyG97CK0rr6zVsqzjnnEuMn6p1zziXGg4pzzrnEeFBxzjmXGA8qzjnnEuNBxTnnXGL+PzMCiXFiQHCyAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "plot_limit = 3000 #plot how the ratio changes oby increasing $n$\n",
        "x = np.arange(1, plot_limit + 1) \n",
        "y = ratios[:plot_limit]\n",
        "plt.title(r\"Win ratio over (n) if we say 'yes'\") \n",
        "plt.xlabel(\"n = number of rounds played\") \n",
        "plt.ylabel(r\"Cumulative win percentage\") \n",
        "plt.axhline(y = 2/3, color = 'k', linestyle='--', alpha = 0.6, label = r\"Probability of winning\")\n",
        "plt.plot(x,y, color='r') \n",
        "plt.xlim(-5,plot_limit)\n",
        "plt.ylim(0.4,0.8)\n",
        "plt.legend()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1caf8548",
      "metadata": {
        "id": "1caf8548"
      },
      "source": [
        "### Modification of the Game\n",
        "#### Every time the host asks if we want to keep our door, let's throw a loaded coin where probability of \"H\" is $\\theta$. If \"H\" shows up, then we say \"yes\"; otherwise we say \"no\" and keep our door."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "85876bf8",
      "metadata": {
        "id": "85876bf8"
      },
      "source": [
        "Question: Derive the probability of winning with this strategy as a function of $\\theta$. "
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2e3c547d",
      "metadata": {
        "id": "2e3c547d"
      },
      "source": [
        "Answer: Let us simulate for, e.g., $\\theta = 0.6$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "df290924",
      "metadata": {
        "id": "df290924"
      },
      "outputs": [],
      "source": [
        "def simulate_with_coins(picked_door, theta):\n",
        "    \"\"\"\n",
        "    Inputs: picked_door (0 1 or 2), decision of switching the door (\"yes\" or \"no\")\n",
        "    Returns: whether or not we win the prize (1: win, 0: lose)\n",
        "    \"\"\"\n",
        "    doors = [0,1,2]\n",
        "    winning_door = randrange(3) # it can be 0, 1, or 2! We don't know the truth yet.\n",
        "    doors.remove(picked_door)\n",
        "    open_door = reveal_empty_door(doors, winning_door)\n",
        "    doors.remove(open_door)\n",
        "    remaining_door = doors[0]\n",
        "    flip = np.random.uniform(0,1) \n",
        "    if flip <= theta:\n",
        "        picked_door = remaining_door\n",
        "    return int(picked_door == winning_door)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "a711587e",
      "metadata": {
        "id": "a711587e",
        "outputId": "4e64a27c-eb3c-4a71-883b-ab686e9938d2"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "We win with a 0.5305666666666666 fraction of the time!\n"
          ]
        }
      ],
      "source": [
        "picked_door = 1 #fix a door (0,1, or 2)\n",
        "theta = 0.6\n",
        "simulate = 30000 #number of times to simulate\n",
        "results = np.zeros(simulate) #a lot of zeros\n",
        "for i in range(simulate):\n",
        "    results[i] = simulate_with_coins(picked_door, theta)\n",
        "ratios = np.cumsum(results) / (np.arange(1,simulate+ 1))\n",
        "print(\"We win with a\", ratios[-1], \"fraction of the time!\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e0a2c42a",
      "metadata": {
        "id": "e0a2c42a"
      },
      "source": [
        "This is not surprising! We previously observed that $\\mathbb{P}[\\text{Win} | \\text{\"yes\"}] = 2/3$ and $\\mathbb{P}[\\text{Win} | \\text{\"no\"}] = 1/3$. Hence, in this new strategy, by the law of total probability, we can derive:\n",
        "\\begin{align}\n",
        "\\mathbb{P}[\\text{Win} | \\text{\"yes\" if \"H\"}] &=  \\mathbb{P}[\\text{Win} | \\text{\"yes\"}] \\times \\theta +  \\mathbb{P}[\\text{Win} | \\text{\"no\"}] \\times (1- \\theta) \\\\\n",
        "& = 2/3 \\times 0.6 + 1/3 \\times 0.4 \\\\\n",
        "& = 0.534.\n",
        "\\end{align}"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3bef6b17",
      "metadata": {
        "id": "3bef6b17"
      },
      "source": [
        "### Optional exercise questions\n",
        "1. In the above paragraph, why did we say \"this is not surprising\"? What is the relationship between the ratio of rounds that we win versus the probability of winning a single game?\n",
        "2. The codes in this notebook are not optimized, and there are several inefficiencies. This is due to teaching purposes, as it is clearer to read what happens in a single iteration by going in a for-loop. After you understand the notebook well, please improve the efficiencies and share your attempts with us/students.\n",
        "3. We always fix picking the door \\#1 in the beginning. Simulate the case where we also randomly pick a door in the beginning of each round. Does it make any difference?\n",
        "4. You play a variant of this game where you pay £1 to play a round. If the door you pick wins, then you make a profit of £1. Otherwise, you lose your £1. You have £5, and you play with the above mentioned strategy of saying \"yes\" (i.e., switch the door) if a loaded coin comes as \"H\". If $\\theta = 0.4$, what is the probability that you will survive at least until the $10$-th round? Simulate and then mathematically derive this probability."
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8.8"
    },
    "colab": {
      "name": "PCMLAI-M3-Monty.ipynb",
      "provenance": [],
      "collapsed_sections": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}