Game_of_chance_simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Game of chance\n",
    "\n",
    "### Dr. Tirthajyoti Sarkar, Fremont, CA, May 2019\n",
    "\n",
    "Two players play a turn-based game. The player who deals first wins the game with a 90% probability and gets to choose the dealer for the next game.\n",
    "\n",
    "Any rational person would choose herself as the dealer to maximize the chance of winning the next game too.\n",
    "\n",
    "Player 'A' gets the first chance of the first game. Obviously, player 'B' is at a disadvantage.\n",
    "\n",
    "But is she?\n",
    "\n",
    "**What happens if this game is played infinite times?**\n",
    "\n",
    "Can we simulate the game outcomes and show that this **game is symmetric in the limit of infinite repetition** i.e. no matter which players starts the game, for a large number of repeated plays, they both win the almost same number of games?\n",
    "\n",
    "<br>\n",
    "![game_of_chance](https://www.gamerbolt.com/wp-content/uploads/2017/11/dice-780x405.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from random import randint\n",
    "from matplotlib.pyplot import plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_games = 1000\n",
    "A_wins=0\n",
    "B_wins=0\n",
    "A_deals=True\n",
    "B_deals=False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate a game with the given probability and assuming that the winner of a game keeps herself as the dealer for the next game"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(n_games):\n",
    "    if A_deals:\n",
    "        rv = randint(1,10)\n",
    "        if rv < 10:\n",
    "            A_wins+=1\n",
    "            A_deals = True\n",
    "            B_deals=False\n",
    "        else:\n",
    "            B_wins+=1  \n",
    "            B_deals=True\n",
    "            A_deals=False\n",
    "    if B_deals:\n",
    "        rv = randint(1,10)\n",
    "        if rv < 10:\n",
    "            B_wins+=1\n",
    "            B_deals = True\n",
    "            A_deals=False\n",
    "        else:\n",
    "            A_wins+=1  \n",
    "            A_deals=True\n",
    "            B_deals=False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print the number of total wins. If the number of games is large (on the order of thousand) then the win numbers should be close to each other"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "477\n",
      "579\n"
     ]
    }
   ],
   "source": [
    "print(A_wins)\n",
    "print(B_wins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define a function to simulate the game for a given number of times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def play_game(n_games,A_first=True,B_first=False):\n",
    "    \"\"\"\n",
    "    Simulates a game play for a certain number of times and computes the number of wins for the players\n",
    "    \"\"\"\n",
    "    A_wins=0\n",
    "    B_wins=0\n",
    "    A_deals=A_first\n",
    "    B_deals=B_first\n",
    "    for i in range(int(n_games)):\n",
    "        if A_deals:\n",
    "            rv = randint(1,10)\n",
    "            if rv < 10:\n",
    "                A_wins+=1\n",
    "                A_deals = True\n",
    "                B_deals=False\n",
    "            else:\n",
    "                B_wins+=1  \n",
    "                B_deals=True\n",
    "                A_deals=False\n",
    "        if B_deals:\n",
    "            rv = randint(1,10)\n",
    "            if rv < 10:\n",
    "                B_wins+=1\n",
    "                B_deals = True\n",
    "                A_deals=False\n",
    "            else:\n",
    "                A_wins+=1  \n",
    "                A_deals=True\n",
    "                B_deals=False\n",
    "    \n",
    "    return (A_wins,B_wins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(572, 486)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "play_game(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the game using `play_game` function over a range of repetition and plot the A/B winning ratio. Initially the ratio will be random but it should tend to one as the number of repetitions grow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_b_ratio=[]\n",
    "for i in range(40):\n",
    "    a,b = play_game(10**(1+i/8))\n",
    "    if b!=0:\n",
    "        a_b_ratio.append(a/b)\n",
    "    else:\n",
    "        a_b_ratio.append(a/(b+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1fb4a34c400>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VOXZP/Dvk0lIkCQsJiwBWcJiiAllVwzVElQUqxYVqiiXtXXBYqtt1bpdatGqb923qvwqUjap+gpIVGT3dSlgQgAjQXaURQIKZCELmbl/f8zM8cxkkjlntnOGfD/XNVfmTM7MPHPEb+65n7MoEQEREcWPBKsHQERE5jC4iYjiDIObiCjOMLiJiOIMg5uIKM4wuImI4gyDm4gozjC4iYjiDIObiCjOJEbjRTMyMqR3797ReGkiolNSSUnJERHJNLJuVIK7d+/eKC4ujsZLExGdkpRSe42uy1YJEVGcYXATEcUZBjcRUZxhcBMRxRkGNxFRnGFwExHFGQZ3BIgItm7disbGRquHQkStAIM7Au666y4MHDgQBQUF4KXgiCjaGNwRsGDBAgDA+vXrsXev4X3oiYhCwuCOgOPHj2v3a2pqLBwJEbUGDO4wOZ1On7A+ceKEhaMhotaAwR2m6upqn2UGNxFFG4M7TJWVlT7LbJUQUbQxuMNUVVXls8yKm4iijcEdJv+Km8FNRNHG4A4TK24iijUGd5hYcRNRrDG4w+RfcXNykoiijcEdJlbcRBRrDO4wMbiJKNYY3GHi5CQRxRqDO0ysuIko1hjcYeLkJBHFGoM7TKy4iSjWEo2spJTaA6AKgBNAo4gMj+ag4gmDm4hizVBwe4wRkSNRG0mc4uQkEcUaWyVhYsVNRLFmNLgFwDKlVIlS6pZoDijecHKSiGLNaKukQEQOKKU6A1iulNoqIv+nX8ET6LcAQM+ePSM8TPtixU1EsWao4haRA56fFQAWAhgZYJ0ZIjJcRIZnZmZGdpQ21djYiNraWp/HGNxEFG1Bg1sp1U4plea9D+AiAGXRHlg88G+TAEBdXR1cLpcFoyGi1sJIq6QLgIVKKe/680VkaVRHFSf82yRetbW1aNeuXYxHQ0StRdDgFpFdAH4Wg7HEnUAVN+CeoGRwE1G0cHfAMDRXcbPPTUTRxOAOQ3MVN4ObiKKJwR0GVtxEZAUGdxiaC24ehENE0cTgDgNbJURkBQZ3GNgqISIrMLjDwIqbiKzA4A6DvuL2HKAEgMFNRNHF4A6DvuLOyMjQ7nNykoiiicEdBn3F3aVLF+0+K24iiiYGdxj0wd21a1ftPoObiKKJwR0GfauEwU1EscLgDgMrbiKyAoM7DM1V3JycJKJoYnCHgZOTRGQFBneIGhoaUF9fDwBwOBw4/fTTtd8xuIkomhjcIdK3SdLT030unMDgJqJoYnCHSN8mSUtLY3ATUcwwuEPkX3Gfdtpp2jInJ4komhjcIfKvuPXBzYqbiKKJwR2ilipuBjcRRRODO0SsuInIKgzuEOmDOz09HW3bttWWa2tr4XK5rBgWEbUCDO4Q+bdKEhISmoQ3EVE0MLhD5N8qAcB2CRHFBIM7RP4VN8DgJqLYYHCHiBU3EVnFcHArpRxKqVKlVFE0BxQvAlXcPHqSiGLBTMV9B4DyaA0k3gSruHn0JBFFi6HgVkr1AHApgH9Fdzjxw393QICtEiKKDaMV9/MA7gHAnZM9ODlJRFYJGtxKqV8CqBCRkiDr3aKUKlZKFR8+fDhiA7QrTk4SkVWMVNwFAC5XSu0BsABAoVJqrv9KIjJDRIaLyPDMzMwID9N+ODlJRFYJGtwicp+I9BCR3gCuAbBKRK6P+shsjpOTRGQV7scdgvr6epw8eRIAkJSUhOTkZABslRBRbCSaWVlE1gBYE5WRxBH/PUqUUgAY3EQUG6y4QxCoTQIwuIkoNhjcIQg0MQlwcpKIYoPBHQIjFTcnJ4koWhjcIWiu4marhIhigcEdAva4ichKDO4QBDpPCcDgJqLYYHCHgJOTRGQlBncIODlJRFZicIeAk5NEZCUGdwg4OUlEVmJwh4AVNxFZicEdAqMVt4jEdFxE1DowuEPQ3O6ACQkJSElJ0ZZra2tjOi4iah0Y3CForlUCsF1CRNHH4A5Bc60SgMFNRNHH4A4BK24ishKD2yQRabHi5tGTRBRtDG6T6urq4HQ6AQDJyclo06aNz+959CQRRRuD26Tm9ijxYquEiKKNwW1SS20SgMFNRNHH4DappYlJgMFNRNHH4DYpWMXNyUkiijYGt0lmKm5OThJRNDC4TWKPm4isxuA2iXuVEJHVGNwmcXKSiKzG4DaJk5NEZLWgwa2USlFKrVdKbVJKfa2U+lssBmZXnJwkIqslGlinHkChiFQrpZIAfKaU+khE1kZ5bLbEyUkislrQ4Bb3ZVyqPYtJnlurvbQLe9xEZDVDPW6llEMptRFABYDlIrIuusOyL1bcRGQ1Q8EtIk4RGQygB4CRSqk8/3WUUrcopYqVUsWHDx+O9DhtI9jugJycJKJoM7VXiYgcA7AGwMUBfjdDRIaLyPDMzMwIDc9+ODlJRFYzsldJplKqg+d+WwAXANga7YHZFVslRGQ1I3uVdAPwb6WUA+6gf1tEiqI7LPvi5CQRWc3IXiWbAQyJwVhsT0R8gpsVNxFZgUdOmnDixAm4XC4AQNu2bZGY2PTvnn9wu/emJCKKHAa3CcH2KAEAh8OB5ORkAO4Kva6uLiZjI6LWg8FtQrCJSS+2S4gomhjcJgSbmPRicBNRNDG4TWDFTUR2wOA2wWjFzaMniSiaGNwmhFJx8+hJIoo0BrcJRvYqAdgqIaLoYnCbwMlJIrIDBrcJnJwkIjtgcJsQyuQke9xEFGkMbhNYcRORHTC4TWCPm4jsgMFtAituIrIDBrcJ3B2QiOyAwW0CJyeJyA4Y3CawVUJEdsDgNoGTk0RkBwxug1wul09wp6amNrsug5uIoonBbVB1dbV2v127dnA4HM2uy+AmomhicBtktE0CcHKSiKKLwW2Q0YlJgBU3EUUXg9sgMxU3g5uIoonBbRArbiKyCwa3Qay4icguGNwGhVpx19TUQESiNi4ian0Y3AYZPU8JACQmJqJNmzYAABFBfX19VMdGRK1L0OBWSp2hlFqtlCpXSn2tlLojFgOzGzOtEoDtEiKKHiMVdyOAv4jIQADnAJimlMqN7rDsx0yrBGBwE1H0BA1uETkoIhs896sAlAPoHu2B2Q0rbiKyC1M9bqVUbwBDAKwL8LtblFLFSqniw4cPR2Z0NmK24ubRk0QULYaDWymVCuB/AdwpIpX+vxeRGSIyXESGZ2ZmRnKMtsCKm4jswlBwK6WS4A7teSLyXnSHZE/scRORXRjZq0QBeANAuYg8G/0h2ZOZ3QEBBjcRRY+RirsAwBQAhUqpjZ7b+CiPy3bitVWydu1aXHHFFZg5c6ZlYyCiyEoMtoKIfAZAxWAsthavk5O33norNm/ejCVLluC8885Dv379LBsLEUUGj5w0KB4r7qNHj2Lz5s0A3Edwvvdeq5yeIDrl2Cq4a2pqsGnTJjQ0NFg9FB9Op9OnatZX082JZHB/8cUXWL16telznnz55Zc+ywsXLgxrHERkD7YJ7tGjRyM1NRWDBw/Gtm3brB6OD321nZaWhoSE4JstUsG9atUqFBQUoLCw0HTFvH79ep/ltWvXYv/+/SGPhYjswTbBre8bb9++3cKRNGW2TQJELriLioq0+wsWLDD1XP/gBoBFixaFPBYisgfbBHf//v21+3aruM1OTAKRm5wsLy/X7n/++eeG2yUiEjC42ecmin+2DG5W3D/ZunWrdv/gwYPYvXu3oed99913OHToEACgbdu22uOffPIJfvjhh5DHQ0TWY3AbEErFHYngPnHiBPbu3evz2GeffWbouevW/XQ6mXPPPRejRo0C4J5oXbJkSUjjISJ7YHAbYFXFvW3btiatkc8//9zQc/VtkpEjR+LKK6/UltkuIYpvtgnuXr16ISkpCYC7JVBdXW3xiH5iVcWt7297Ga24/YN7woQJ2vKyZct8/hgRUXyxTXAnJiYiOztbW7ZT1W32PCVAZCYn9f1try1btuDHH39s8XmNjY0oLi7WlkeOHIm+ffti0KBBAID6+nosXbo0pDERkfVsE9yAfdslVrVKAlXcgPuAnGDP875n9+7dkZWVBQBslxCdIhjcBljVKtFX3Oeee652P1i7RD8xefbZZ2v39cFdVFSEurq6kMZFRNZicBtgRcXtdDp99me/6aabtPvBJij9+9teeXl52kmmqqursXLlStPjIiLr2Sq4BwwYoN2300E4VlTce/bsQX19PQCga9euGD/+pzPprl+/vsVqubngVkr5TFLy3CVE8clWwW3Xilt/Dc0OHToYeo4+uGtqakyfIErf3x44cCC6dOmibZ+GhgaUlJQEfF5NTQ3KysoAuIN62LBhPr/Xt0sWL16MxsZGU+MiIuvZKrh79OiBlJQUAMCRI0dw7Ngxw8/99NNPcc8992DXrl0RH9c333yj3df/cWlJUlKStnujy+UyfcZDfX87JycHAFBQUKA91ly7pLS0FE6nE4A78P1bOyNHjtQmK48cOWJ490Iisg9bBXdCQgL69u2rLRutuisrK3H55Zfjqaeewrnnnovvv/8+YmOqqanBt99+CwBwOBw+uywGE067xL/iBtxnUPRqLnD1bRL9xKRXQkKCT7uEe5cQxR9bBTfg2+c2GtwlJSVadX7o0CFMnjxZqzrDpR9D37590aZNG8PPDSe4g1XcX3zxBVwuV5Pn6fco0fe39fz73GbbOERkLdsFdyhnCSwtLfVZXr16NR577LGIjEcfoGeeeaap54Ya3CISsOI+88wzcfrppwMAfvjhB58WjldzE5N65513Hjp16gQA2Ldvn8/BOkRkf7YObqMV94YNG5o89re//Q2rVq0Kezz6cPRWvkaFevTk4cOHcfToUQBAamoqunfvDsA92dhSu6SiogJ79uwBACQnJyM/Pz/g6yclJeHyyy/XltkuIYovp0Rw6yvu3r17A3BXrdddd512atNQWVFx66vtnJwcKPXTtZpbmqDUX6ps6NCh2uRoIP5HUbJdQhQ/bB/cwQKlpqZGC1elFFasWIHOnTsDAL7//ntcf/31YfW7w6m4Qw3uQP1tr5YqbiNtEq8LL7xQ+0awbdu2Zg+vJyL7sV1wd+vWTQuUY8eOBT3p/+bNm7VJupycHPTt2xfz5s3TqtQVK1bg8ccfD2ksLpfLJ7itqLi9/W2voUOHartM7ty502cPmmB7lOilpKT4HNTDdglR/LBdcCulTE1Q6tskQ4cOBQBccMEFePDBB7XHH3nkEaxZs8b0WPbt26cF7umnn46MjAxTz49GxZ2cnIwRI0Zoy952if+lyoJV3IBvu2T69Om48sorsWjRItP7nFvlxIkT2LFjh9XDIIo52wU3YK7PrZ+YHDJkiHb/4Ycfxvnnnw/AXTlPnjwZFRUVpsYRTrUNhD452VLFDQRul+zcuVM73WunTp0M7W8+fvx4bYwnT57EwoULMWHCBHTv3h1//OMfUVxcbMvet4jgjTfeQLdu3dC/f39MnTo14K6RVtq9ezfWrFkTsd1SifROqeD2VtyA+2CZ+fPnIzMzE4D74gxTpkwx9T94S5WvEaFU3P4H/OgPSPIKNEHpX23rJzSbk56ejvfee69JW+XIkSN46aWXMGLECOTl5eHJJ59EeXm5LUJ87969GDduHG666SbtHDKvv/46/vznP1s+vrq6Orz11lsoLCxEdnY2xowZgzFjxpguGIiCievgbmho0M7LAQCDBw/2+X1WVhbmzp2rhdiyZcvw8ssvGx5HOBOTQGjBrX/Pfv36BTzgR3+K1w0bNqCmpsZ0m8Troosuwtq1a7F161bcf//9OOOMM3x+v2XLFtx3333Izc1FdnY2pk2bhqKioqDfIKqqqrBhwwa8++672Lx5s+HxNMflcuGf//wn8vLysHz58ia/f+GFF/DII4+E/T6h+Prrr3HnnXeie/fumDx5MlavXq397tNPP8WIESOaHGvQki+//BJXXXUVbr75Zhw4cCAaQ6Z4JyIt3gDMBFABoCzYut7bsGHDJByff/65ABAAMmTIkGbX27Bhg7Zenz59ml3v/vvv19Y7++yzDY9j7Nix2vMWL15s6jOIiDz88MPa8x966CFDz5k3b572nCuuuKLZ9fLy8rT1Vq5cKaNGjdKWP/jgA9Nj9XI6nbJy5Uq54YYbpF27dtpr+t+Sk5Nl3Lhx8vzzz8v7778vzzzzjNxyyy3yi1/8Qrp16+azrlJKZs+eHfKYtm/fLueff77PayYkJMhf/vIXufrqq30ef/rpp0N+H736+npZuXKlLFq0SBYvXixLliyRDz74QD766CNZunSpLFu2TP71r3/5bHf/8SmltOW2bdvKggULWnzPY8eOybRp03ye16FDB3nzzTfF5XJF5HORfQEoFoMZm2gg22cBeBnAbNN/FULkPzkpIgG/+geamAzkzjvv1PYs2bhxIxoaGgwdum5FxR2sv+1VUFCgfdtYs2aNT8tIP3lpVkJCAgoLC1FYWIiXX34ZCxcuxKJFi7B8+XKf85LX19fj448/xscffxz0NUUEN9xwA5RSuP766w2Pxel04sUXX8QDDzyA2tpa7fGBAwdi5syZOOecc9DQ0ICamhp89NFHAIC77roL6enpuPnmm0186p+cPHkSs2fPxqOPPoq9e/eafn7v3r3xu9/9DjfeeCM2bdqEa6+9FpWVlaitrcU111yDzZs349FHH0VCwk9fdkUEb7/9Nu68884m59k5duwYbrzxRrz99tt4/fXXm3wjCmTPnj0oKSlB9+7dkZ+f7zPXEm8aGhrwzTffoKysDF999RXKyspw9OhRpKeno3379tpP7y09PR0OhwP19fWoq6vTfnrv19fXIyUlBWlpaUhPT29ya9euHRoaGlBbW4va2lqcOHHC52ddXR0aGxu1m9PpbLI8ePBgU//OQ2Ik3QH0RgwrbpfLJe3bt9eqjgMHDgRcb9q0ado6jz32WIuv2atXL23dkpKSoGOoqqrS1k9MTJSGhgbTn+Pll1/WXuO2224z9JyrrrpKe86sWbOaXW/OnDnaevoKt6VvHuGor6+X1atXy9133y1nnXVWs5W495aUlCQ5OTnSs2dPnyp07ty5ht5v9+7dMnr0aJ/XdDgccv/990ttba3PujU1NfLzn//cp8J/6623TH2+kydPyptvvinZ2dlBP1ugzzpx4kRZtmyZOJ1On9ctLy+XAQMG+Kx/2WWXyfHjx0XE/W3ioosuavKaF154YZOxpKWlyYwZMwJW3xUVFfLKK69IQUFBk287AwYMkEmTJsnjjz8uH374oRw4cMBwBd/Y2CgVFRVSVlYmq1evlv/85z/y0ksvyXPPPSevvPKKzJgxQ2bNmiXz58+Xd999V95//31ZunSpfPbZZ1JWVibfffedVFVVBXw/l8slVVVVsmvXLlm3bp0UFRXJrFmzZPr06TJp0iTJzc2VxMRE0/89rL5de+21pv7tecFExR2x4AZwC4BiAMU9e/YMaeB6w4cP1zbEJ598EnAd/dfUDz/8sMXXmzhxorbua6+9FvT9S0pKtPVzcnJC+gwzZ87UXuOGG24w9Bx9KK5bt67Z9Xbv3h3wH82vf/3rkMZq1p49e+TVV1+VCRMmyNixY2Xq1Kny3HPPyQcffCDbt2+XkydPiojI4cOHJT8/3ye8582b1+Jrz58/X9LT030+16BBg6S4uLjZ5xw/flyGDRvm88e2qKgo6OdobGyUOXPmSL9+/Zpsy4yMDLnsssvkl7/8pYwfP14uueQSGTdunFx44YVywQUXyCWXXCJPP/20VFRUtPgeR48elYsvvtjntQcOHCgPPPCAJCcn+zyelZUl77zzjrhcLqmurpY77rjDp3XiDfU9e/ZIZWWlzJkzRy655BJxOBymwqV9+/bStWtX6dGjh/Tp00f69+8vubm5MmjQIBk2bJjk5eVJ586dJSEhISJh5nA4pFOnTpKdnS35+fnSo0cPSUlJsTxko3GbOHFi0H93gZgJbiOtEkNEZAaAGQDgCd2wDBgwQDv50fbt23Heeef5/N7pdGLTpk3ackutEsDdPnjnnXcAuCd/br311hbXD+dQdy+zrZLGxkaf/dZbas/06tULWVlZTSavzExMhqNXr16YOnUqpk6d2uJ6GRkZWLlyJcaOHYuvvvoKLpcLU6ZMAQBMnjzZZ93KykrcfvvtmDNnjvaYw+HAgw8+iPvvv7/F9lZ6ejqWLl2K888/H1u2bEFjYyOuuuoqfPTRRzj77LNRVVWl3aqrq1FVVYV9+/bh2WefbXKyrk6dOuHuu+/G7bffjtTUVLObpokOHTqgqKgI9913H5566ikA7pbY3//+d22dhIQETJs2DY899ph2DvV27drh+eefx9VXX43f/va32kT98uXLcdZZZ8Hlcvm0kPTbrKCgAEeOHMHWrVsD7kl1/PhxHD9+POzPZpTT6cSPP/6o7bJqRp8+fZCXl6fdsrKyUFVVhePHj6OyslL7LN6biCA5ORkpKSnazbvcpk0b1NfXo7KyMuCturoaycnJaNu2Ldq2bYvTTjvN535KSgqSkpLgcDiQmJio3fTLoeaFGREL7kgLtmfJtm3btDDs1q0bunTp0uLr6fu++nN6NCfc/jZgPrh3796NkydPAnDvEdPS9S29J5x6++23fR4PdsSkFTIzM7Fy5UoUFhairKxMC2+lFK699loAwNq1a3Hdddf5XAijT58+mD9/Ps455xxD75ORkYHly5fj5z//OXbt2oX6+noUFhYaHmeHDh1w11134Q9/+IPha4sa5XA48I9//AODBg3CTTfdpF2WDgCGDRuG1157DcOHDw/43NGjR2Pjxo146KGH8Oyzz0JEAu7VU1BQgMmTJ2PixInabrC1tbUoKyvDpk2bsHHjRmzcuBGbN2/2ma8IpmPHjujcuTM6d+6MzMxMZGZmIiUlBQ0NDQFv9fX1WrAeO3YMx48fD/gHBnAfwet9zczMTGRkZKBLly7Izc1FXl4ecnNzDV8usFUxUpYjxj1uEZG5c+dqXz0mTJjQ4u8vvfTSoK9XWVmpfeV0OBxSU1PT4vqTJk3SXn/mzJkhfYYVK1ZorzFmzJig6y9evFhbf+zYsUHXf+GFF5p8HQ32uax06NAhn1aQt+c9ffr0Jl/1p0yZovWBzdq1a5dkZWUZ/mqbnp4ujzzyiBw7dizCnziw9evXS25urnTt2lVefPFFaWxsNPzcL774QnJycrSx5+fnyxNPPCG7d+82/BpOp1MOHz4sBw4ckL1798qOHTtk69at8tVXX0lpaamsX79eSktLZf/+/VJfXx/CJ2yqvr5eKioqZPv27VJaWiq7d++W6upq7i2jg0i2SpRSbwH4BYAMpdQ+AA+LyBvh/LEwIljFrd+jRH/EZHPS0tKQk5OD8vJyOJ1OlJaW+hzI4i8SFbfZIyfNHvCjP4ISAPLz832qfLvp3LkzVq1ahcLCQnz99ddwuVxNZt/bt2+PV199VavEQ9GnTx8sX74cl19+OXbu3Ink5GSkpaX53FJTU5GWloYhQ4bgtttuQ8eOHcP9eIaNGDHC57qgZowaNQqlpaVYuXIlevXqhby8PNPvn5CQYPr0DeFq06aNVlVT+IIGt4iE/n9QGPTBvWPHDrhcLp9dqJo7YrIlI0aM0Ha3+/LLL5sNbpfL5dNrjlWP2+iugF6DBg1CamoqqqurAcSuvx0Ob3iPGTMGW7Zs8fnd6NGjMXfuXPTq1Svs98nNzcX27dvhdDqRmGi/jqDZwNZLSUnBpZdeGsHRULyx5ZGTgLuv5q0K6urqsH//fu13ImJ4H249fZ9bf6Shv++++07ryWVmZmpXizHLbHCbrbgTExMxatQobTkeghv4Kbxzc3MBuPu/06dPx+rVqyMS2l5KKVuGNlG4bP2vun///jhy5AgA92Sk9+CDPXv2aNeY7NixI3r27Gno9YxOUEZijxLAXHCLBL5cWTD33HMP1q5dizPOOAOTJk0KbaAW6NKlC9atW4eFCxdi+PDhhj8vEdm44gaa73P7t0mMfu382c9+pl0VZseOHdrlwfyFe3IpLzPBfejQIW33rLS0NHTr1s3Qe1xwwQU4fPgwysrK4m72PTU1FVOmTGFoE5kUl8EdSpsEcPcGBw0apC03d5HccE/n6uU/OSktnL3Ov9o20wNNTk4Oq2dKRPElLoO7uXNwG2GkXRKpijspKUnrsTqdTm0f7Wi+JxGd+mwd3AMGDNDut9QqMcNIcEdiV0Avo+2SUPrbRNQ62Tq4+/Xrp93fuXMnGhsbcfDgQe3K7e3atfOpyo0IFtyVlZXaYeRJSUnaVeNDZTS4WXETkVG2Du60tDR07doVgPt0m99++61PtT148GCffbuNGDhwoBam+/fvx8GDB31+r99/u3///mHvTsaKm4gizdbBDTTtc4c6MemVmJjo8zz/qjtSuwJ6GTl60nvCI+/4jFwvkohaL9sHt3+fO5yJSa+W2iWR7G8Dxipu/Xv2799f22WRiCgQ2we3f8UdzsSkl/4IQ/8jKCNdcRsJbva3iciMuArutWvXapeTatOmjXbItFn6iru4uNhn/2orKm72t4nIjLgKbn11nJ+fH3JLITs7Wzv/yI8//qidA9rpdEbk5FJ6+uBursfNipuIzLB9cPft2zfg46G2SQD3yYf0J6339rm//fZb7QT3Xbp0QYcOHUJ+Dy/95CQrbiKKBNsH92mnnRbwytahTkx6BZqgjHR/GwjeKjl58iR27NgR8fclolOX7YMbQMCDbMKpuIHgwR2plkWw4N6xY4d2KHyPHj3i7kRRRBR7cRncDofD52RRodAH94YNG+B0OiN2cim9loJbRHDPPfdoy6FczYSIWp+4DO6cnBy0bds2rNfMyspCVlYWAPekYXl5edQrbv/JyRdffBFFRUXa8h133BGR9ySiU1tcBne4bRIv/3ZJNCru5iYnS0pKcPfdd2vLf/rTn3DxxRdH5D2J6NQWF8GtP3oSiFxw6w/EWbFiBb7//nsA7n3Ewz03p9bSAAAG0ElEQVS5lFegVklVVRWuueYarbc9bNgwPPHEExF5PyI69cVFcGdnZ/ucTCrcPUq89BX3woULtfsDBgyAw+GIyHv4B7eI4LbbbtP2JElNTcWCBQuQnJwckfcjolNfXAR3mzZttJDt0KFDxCpu/b7c3osDA5HdJc8/uGfPno158+Zpj73++us+p68lIgomLoIbAGbPno2//vWvWLJkScR2mevYsWPA0Izk0Yv64C4vL8e0adO05d/85jeYPHlyxN6LiFoHW1/lXW/AgAF48sknI/66I0aM8DkABohsxa2fnNy5c6fPe7z00ksRex8iaj3ipuKOFn2f2ytaFbdXcnIyFixYgNTU1Ii9DxG1HgzuAMEdrR631zPPPIPBgwdH7D2IqHUxFNxKqYuVUt8opXYope6N9qBiaciQIT57rHTr1g3p6ekRe33/4P7Vr36F3//+9xF7fSJqfYIGt1LKAeAVAJcAyAVwrVIqtBNh21C7du18DjWP9EmeOnfujIyMDADAGWecgTfeeANKqYi+BxG1LkYq7pEAdojILhFpALAAwBXRHVZs6dslkT4fdlJSElatWoUnnngC//3vf7XzgBMRhcpIcHcH8J1ueZ/nsVPGxIkTtfuXXnppxF8/Pz8f9957L7p3P6U2GxFZxMjugIG+10uTlZS6BcAtANCzZ88whxVb48aNw2effQaXy4XRo0dbPRwiohYZCe59APRXMugB4ID/SiIyA8AMABg+fHiTYLe7goICq4dARGSIkVbJlwD6K6X6KKXaALgGwPvRHRYRETUnaMUtIo1KqdsBfAzAAWCmiHwd9ZEREVFAhg55F5EPAXwY5bEQEZEBrf7ISSKieMPgJiKKMwxuIqI4w+AmIoozSiTyu1wrpQ4D2Bvi0zMAHIngcCKJYwsNxxYaji008Tq2XiKSaeRFohLc4VBKFYvI8OBrxh7HFhqOLTQcW2haw9jYKiEiijMMbiKiOGPH4J5h9QBawLGFhmMLDccWmlN+bLbrcRMRUcvsWHETEVELbBPcdr6upVJqj1LqK6XURqVUsQ3GM1MpVaGUKtM91kkptVwptd3zs6ONxvaIUmq/Z/ttVEqNt2BcZyilViulypVSXyul7vA8bvl2a2FsdthuKUqp9UqpTZ6x/c3zeB+l1DrPdvuP58yhdhnbLKXUbt12s+zK3Eoph1KqVClV5FmOzHYTEctvcJ91cCeAbABtAGwCkGv1uHTj2wMgw+px6MZzHoChAMp0j/0DwL2e+/cC+B8bje0RAHdZvM26ARjquZ8GYBvc11C1fLu1MDY7bDcFINVzPwnAOgDnAHgbwDWex18DcJuNxjYLwNVWbjfdGP8MYD6AIs9yRLabXSruU/66lpEkIv8H4Ee/h68A8G/P/X8D+FVMB+XRzNgsJyIHRWSD534VgHK4L8Fn+XZrYWyWE7dqz2KS5yYACgG863ncqu3W3NhsQSnVA8ClAP7lWVaI0HazS3Db/bqWAmCZUqrEc4k2O+oiIgcBdxAA6GzxePzdrpTa7GmlWNLG8VJK9QYwBO4KzVbbzW9sgA22m+fr/kYAFQCWw/3t+JiINHpWsez/V/+xiYh3u/3ds92eU0olWzE2AM8DuAeAy7N8OiK03ewS3Iaua2mhAhEZCuASANOUUudZPaA48yqAvgAGAzgI4BmrBqKUSgXwvwDuFJFKq8YRSICx2WK7iYhTRAbDfdnCkQAGBlottqPyvKnf2JRSeQDuA5ADYASATgD+GutxKaV+CaBCREr0DwdYNaTtZpfgNnRdS6uIyAHPzwoAC+H+x2s3h5RS3QDA87PC4vFoROSQ538wF4D/B4u2n1IqCe5gnCci73ketsV2CzQ2u2w3LxE5BmAN3H3kDkop74VYLP//VTe2iz2tJxGRegBvwprtVgDgcqXUHrhbv4VwV+AR2W52CW7bXtdSKdVOKZXmvQ/gIgBlLT/LEu8DuMFz/wYAiy0ciw9vMHpMgAXbz9NffANAuYg8q/uV5dutubHZZLtlKqU6eO63BXAB3D341QCu9qxm1XYLNLatuj/ECu4ecsy3m4jcJyI9RKQ33Hm2SkSuQ6S2m9WzrrrZ1/Fwz6bvBPCA1ePRjSsb7r1cNgH42g5jA/AW3F+dT8L9beV3cPfPVgLY7vnZyUZjmwPgKwCb4Q7KbhaMazTcX0s3A9jouY23w3ZrYWx22G6DAJR6xlAG4CHP49kA1gPYAeAdAMk2Gtsqz3YrAzAXnj1PrLoB+AV+2qskItuNR04SEcUZu7RKiIjIIAY3EVGcYXATEcUZBjcRUZxhcBMRxRkGNxFRnGFwExHFGQY3EVGc+f/ztrBBCUkzcQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(a_b_ratio,color='k',linewidth=3)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}