Galton-Watson.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Member 0 had this number of children: 3\n",
      "Population is now: 3\n",
      "Member 0 had this number of children: 3\n",
      "Member 1 had this number of children: 2\n",
      "Member 2 had this number of children: 5\n",
      "Population is now: 10\n",
      "Member 0 had this number of children: 3\n",
      "Member 1 had this number of children: 2\n",
      "Member 2 had this number of children: 2\n",
      "Member 3 had this number of children: 2\n",
      "Member 4 had this number of children: 2\n",
      "Member 5 had this number of children: 0\n",
      "Member 6 had this number of children: 0\n",
      "Member 7 had this number of children: 2\n",
      "Member 8 had this number of children: 0\n",
      "Member 9 had this number of children: 3\n",
      "Population is now: 16\n",
      "Member 0 had this number of children: 0\n",
      "Member 1 had this number of children: 4\n",
      "Member 2 had this number of children: 0\n",
      "Member 3 had this number of children: 0\n",
      "Member 4 had this number of children: 0\n",
      "Member 5 had this number of children: 3\n",
      "Member 6 had this number of children: 1\n",
      "Member 7 had this number of children: 0\n",
      "Member 8 had this number of children: 1\n",
      "Member 9 had this number of children: 1\n",
      "Member 10 had this number of children: 2\n",
      "Member 11 had this number of children: 0\n",
      "Member 12 had this number of children: 1\n",
      "Member 13 had this number of children: 0\n",
      "Member 14 had this number of children: 5\n",
      "Member 15 had this number of children: 2\n",
      "Population is now: 20\n",
      "Member 0 had this number of children: 2\n",
      "Member 1 had this number of children: 1\n",
      "Member 2 had this number of children: 2\n",
      "Member 3 had this number of children: 2\n",
      "Member 4 had this number of children: 3\n",
      "Member 5 had this number of children: 0\n",
      "Member 6 had this number of children: 1\n",
      "Member 7 had this number of children: 2\n",
      "Member 8 had this number of children: 4\n",
      "Member 9 had this number of children: 0\n",
      "Member 10 had this number of children: 1\n",
      "Member 11 had this number of children: 3\n",
      "Member 12 had this number of children: 4\n",
      "Member 13 had this number of children: 3\n",
      "Member 14 had this number of children: 4\n",
      "Member 15 had this number of children: 0\n",
      "Member 16 had this number of children: 2\n",
      "Member 17 had this number of children: 1\n",
      "Member 18 had this number of children: 1\n",
      "Member 19 had this number of children: 6\n",
      "Population is now: 42\n",
      "Now for plot\n",
      "Population is now: 2\n",
      "Population is now: 2\n",
      "Population is now: 4\n",
      "Population is now: 14\n",
      "Population is now: 27\n",
      "Population is now: 46\n",
      "Population is now: 132\n",
      "Population is now: 266\n",
      "Population is now: 519\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl0VeW9//H3NxOBAAmQMGTAQEAUkRCNiNhqneoMtld7kVZxqLSrUu3P1lZvb1fv7e/21g6rtV78VbnOVaqtQ4nUeaK2KAomYVYIIGQAEiABAiHDeX5/nB1IICQBzsk+5+TzWivr7P2cnZyvR/jk4bv32Y855xARkdgV53cBIiISXgp6EZEYp6AXEYlxCnoRkRinoBcRiXEKehGRGKegFxGJcQp6EZEYp6AXEYlxCX4XAJCenu5yc3P9LkNEJKosW7asxjmX0dVxERH0ubm5LF261O8yRESiipl93p3j1LoREYlxCnoRkRinoBcRiXEKehGRGKegFxGJcQp6EZEe9tCiMhaX1bQbW1xWw0OLysLyegp6EZEeNjE7lTnziw+G/eKyGubML2ZidmpYXi8irqMXEelNpual8z8zCrjl8Y+5piCLN1ZvY+7MAqbmpYfl9TSjFxHxQZ/EOBqaAzz78Ra+cfbIsIU8KOhFRHzx8KINAHz7/NE8vWTzET37UFLQi4j0sPfXVfPWmm2cPXow91x+KnNnFrTr2Yeagl5EpIe9XFqJA245dxQQ7NnPnVnA8vK6sLxet07GmtkmYA/QAjQ75wrNbDDwHJALbAK+5pzbZWYG/B64AtgH3OSc+yT0pYuIRKfmFseA5AS+NO7QjSen5qVHxMnYC5xzk5xzhd7+PcDbzrmxwNvePsDlwFjvazbwh1AVKyIS7fY3tvD6qq1cMWEEfRLie+Q1T6R1Mx140tt+ErimzfhTLuhDIM3MRpzA64iIxIx31m6nvrGFaZMye+w1uxv0DnjDzJaZ2WxvbJhzrgrAexzqjWcBW9p8b7k31o6ZzTazpWa2tLq6+viqFxGJMgtKKsgY0Icpo4f02Gt29wNT5zrnKs1sKPCmma3t5FjrYMwdMeDcPGAeQGFh4RHPi4jEmrr9Tbz3aTXfmHIS8XEdRWV4dGtG75yr9B63Ay8Bk4FtrS0Z73G7d3g5kNPm27OBylAVLCISrV5fuZXGlkCPtm2gG0FvZilmNqB1G/gysBIoAmZ5h80CFnjbRcCNFjQFqGtt8YiI9GYLSis4aUg/8sN0T5uj6U7rZhjwUvCqSRKA+c6518zsY+DPZnYrsBm4zjv+FYKXVq4neHnlzSGvWkQkymzf3cDish1894IxeHnaY7oMeufcBiC/g/EdwEUdjDvg9pBUJyISIxYur8I5erxtA/pkrIhIj1hQWsn4EQMZM3RAj7+2gl5EJMw21dRTuqWW6T7M5kFBLyISdi+XBi88vDpfQS8iEnOcc/y1pILJuYPJTOvrSw0KehGRMFpdtZuy6npfTsK2UtCLiIRRUWklCXHGFaf7d8svBb2ISJgEAo6XSyr54th0Bqck+VaHgl5EJEyWfr6LyroGpk864r6OPUpBLyISJkWlFSQnxnHJ+GG+1qGgFxEJg6aWAH9bXsXFpw4jpU93bxQcHgp6EZEw+Me6Gnbta/K9bQMKehGRsCgqrSS1byLnn5zR9cFhpqAXEQmx1nVhL58wnKQE/2PW/wpERGLMW2u2sa+H14XtjIJeRCTEikorGTawD2eP6rl1YTujoBcRCaG6fU289+l2rpqY2aPrwnZGQS8iEkKvrqyiqcX5dkvijijoRURCqKi0klHpKZye1bPrwnZGQS8iEiLbdjfwwYYdXJ2f2ePrwnZGQS8iEiIvl1YG14X1aYGRo1HQi4iEyMullUzIGsiYof39LqUdBb2ISAhsrKmntLwu4mbzoKAXEQmJopJKzPxbF7YzCnoRkRPknGNBaXBd2BGp/qwL2xkFvYjICVpVuZsNPq8L2xkFvYjICTq4LuwE/9aF7YyCXkTkBAQCjpdLKzn/5AwG+bgubGcU9CIiJ+DjTTupqmuI2LYNKOhFRE7IgtJK+ibGc/Gp/q4L25luB72ZxZtZsZkt9PZHmdkSM1tnZs+ZWZI33sfbX+89nxue0kVE/NXYHOCVFVVcMt7/dWE7cywz+juBNW32fwn8zjk3FtgF3OqN3wrscs6NAX7nHSciEnP+sb6a2n1NEfkhqba6FfRmlg1cCTzi7RtwIfC8d8iTwDXe9nRvH+/5iyyS7u4jIhIiC0qC68KeFwHrwnamuzP6+4EfAgFvfwhQ65xr9vbLgdalzrOALQDe83Xe8SIiMWNfYzNvrNrGFaePiIh1YTvTZXVmdhWw3Tm3rO1wB4e6bjzX9ufONrOlZra0urq6W8WKiESKt9ZsZ39TS8S3baB7M/pzgWlmtgl4lmDL5n4gzcxazz5kA5XedjmQA+A9nwrsPPyHOufmOecKnXOFGRmR/c8eEZHDFZVUMHxgMpNHDfa7lC51GfTOuXudc9nOuVxgBvCOc+7rwLvAtd5hs4AF3naRt4/3/DvOuSNm9CIi0ap2XyOLPqvm6vwREbMubGdOpLH0I+AuM1tPsAf/qDf+KDDEG78LuOfEShQRiSyvrtxKU4tjWn5W1wdHgGO68NM59x7wnre9AZjcwTENwHUhqE1EJCItKKlgdHoKE7IG+l1Kt0T2qWIRkQizta6BJRt3Mm1SZK0L2xkFvYjIMVi4PDLXhe2Mgl5E5BgsKKnk9KxURmdE1rqwnVHQi4h004bqvayoqGN6BN+psiMKehGRbioqDa4Le9VEBb2ISMxxzlFUUsnZowYzPDXZ73KOiYJeRKQbVlbsZkNNPdMnRce1820p6EVEuqGotILEeOPyCcP9LuWYKehFRLrQEnAUeevCpvWLzHVhO6OgFxHpwkcbd7Jt9wGmRWHbBhT0IiJdKjq4LuxQv0s5Lgp6EZFOtK4L++XThtEvKXLXhe2Mgl5EpBN//6yauv1NUfchqbYU9CIinSgqrSStXyJfGBO9CyQp6EVEjqL+QDNvro6OdWE7E72Vi4iE2VtrtrG/qYXpUXSnyo4o6EVEjqKopJIRqcmclRv568J2RkEvItKBXfWt68JmEhcF68J2RkEvItKBV1ZW0RxwUbXAyNEo6EVEOlBUUkleRgqnZUbHurCdUdCLiBymsnY/H23aybT8rKhZF7YzCnoRkcMcXBc2ij8k1ZaCXkTkMEWlleRnpzIqPcXvUkJCQS8i0sb67XtZWbGbq2PgJGwrBb2ISBut68Iq6EVEYpBzjpdLKzln9BCGDYyudWE7o6AXEfGsqKhjY019TFw735aCXkTEs6Ck0lsXdoTfpYSUgl5EhOC6sAuXV/KlcUNJ7Zfodzkh1WXQm1mymX1kZqVmtsrM/tMbH2VmS8xsnZk9Z2ZJ3ngfb3+993xueP8TRERO3JKNO4LrwsZY2wa6N6M/AFzonMsHJgGXmdkU4JfA75xzY4FdwK3e8bcCu5xzY4DfeceJiES0opJK+iXFc/Gpw/wuJeS6DHoXtNfbTfS+HHAh8Lw3/iRwjbc93dvHe/4ii4XPEItIzDrQ3MKrK7dy6WnD6ZsU73c5IdetHr2ZxZtZCbAdeBMoA2qdc83eIeVAlredBWwB8J6vA4aEsmgRkVD6+2c11O1vism2DXQz6J1zLc65SUA2MBk4taPDvMeOZu/u8AEzm21mS81saXV1dXfrFREJuQUlFQzql8gXxqb7XUpYHNNVN865WuA9YAqQZmYJ3lPZQKW3XQ7kAHjPpwI7O/hZ85xzhc65woyM6F10V0SiW/2BZt5as40rJ44gMT42L0TszlU3GWaW5m33BS4G1gDvAtd6h80CFnjbRd4+3vPvOOeOmNGLiESCN1dvo6EpwLT8rK4PjlIJXR/CCOBJM4sn+Ivhz865hWa2GnjWzP4LKAYe9Y5/FPijma0nOJOfEYa6RURCYkFJBZmpyRSeNMjvUsKmy6B3zi0HCjoY30CwX3/4eANwXUiqExEJo531jby/roZbvzgq6teF7UxsNqRERLrhlRWxsy5sZxT0ItJrFZVUMmZof8aPiP51YTujoBeRXqnCWxd2en5mTKwL2xkFvYj0SgtLg1eEx9ICI0ejoBeRXmlBSSX5OWnkxsi6sJ1R0ItIr/HQojIWl9WwfvseVlftZnp+JovLanhoUZnfpYWVgl5Eeo2J2anMmV/Mg++uJ85g+MA+zJlfzMTsVL9LCysFvYj0GlPz0vnt1/JZUFJJZlpf/n3BKubOLGBqXmze46aVgl5EepUPNuwg4KB8136+cfbImA95UNCLSC9SsqWWeYs20CchjjsuHMPTSzazuKzG77LCTkEvIr3CgeYWbn9mGWbw4NfP4K4vj2PuzALmzC+O+bBX0ItIr/DA2+uoqG3gB5eOO7hc4NS8dObOLGB5eZ3P1YVXd+5eKSIS1VaU1/HQog1ce2Y23/nSmHbPTc1Lj/k+vWb0IhLTGpsD3P18KUNSkvjJleP9LscXmtGLSEyb++561m7dwyM3FpLaL9HvcnyhGb2IxKxVlXX8v3fX85WCLC4eP8zvcnyjoBeRmNTUEuDuvywnrV8SP726d7ZsWql1IyIx6Q/vlbG6ajcPfeNM0vol+V2OrzSjF5GYs3brbv7nnXVcnZ/JZROG+12O7xT0IhJTmr2WzcDkRP5z2ml+lxMR1LoRkZgy7/0NrKio48GZZzA4pXe3bFppRi8iMWPdtj3c/+Y6Lp8wnCsnjvC7nIihoBeRmNAScNz9/HJS+sTzs+kT/C4noqh1IyIx4dF/bKBkSy2/nzGJjAF9/C4nomhGLyJRr6x6L7954zMuGT+Mab1gse9jpaAXkajWEnD88Pnl9E2M5+fXTMDM/C4p4ijoRSSqPbF4E8s+38VPrx7P0IHJfpcTkRT0IhK1NtXU8+vX13LhKUP5SkGW3+VELAW9iESlQMDxwxeWkxgfx39/5XS1bDrRZdCbWY6ZvWtma8xslZnd6Y0PNrM3zWyd9zjIGzcze8DM1pvZcjM7I9z/ESLS+/zxw8/5aONOfnLVeIanqmXTme7M6JuB7zvnTgWmALeb2XjgHuBt59xY4G1vH+ByYKz3NRv4Q8irFpFebfOOffzytbWcd3IG152Z7Xc5Ea/LoHfOVTnnPvG29wBrgCxgOvCkd9iTwDXe9nTgKRf0IZBmZvqImoiERCDg+NELy4kz476vqmXTHcfUozezXKAAWAIMc85VQfCXATDUOywL2NLm28q9MRGREzb/o818sGEHP77yVDLT+vpdTlTodtCbWX/gBeB7zrndnR3awZjr4OfNNrOlZra0urq6u2WISC9Wvmsfv3hlDV8Yk86Ms3L8LidqdCvozSyRYMg/45x70Rve1tqS8R63e+PlQNv/A9lA5eE/0zk3zzlX6JwrzMjION76RaSXcM5x74srAPiFWjbHpDtX3RjwKLDGOffbNk8VAbO87VnAgjbjN3pX30wB6lpbPCIix+u5j7fw/roa7rniVHIG9/O7nKjSnZuanQvcAKwwsxJv7N+A+4A/m9mtwGbgOu+5V4ArgPXAPuDmkFYsIr1OVd1+fv63NUwZPZivTx7pdzlRp8ugd879g4777gAXdXC8A24/wbpERIBDLZvmgONX/5JPXJxaNsdKn4wVkYj2/LJy3vu0mh9dNo6RQ9SyOR4KehGJWNt2N/B/F65mcu5gbjwn1+9yopaCXkQiknOOH7+0ggPNAX557US1bE6Agl5EItKCkkreWrOduy8dx6j0FL/LiWoKehGJONv3NPDTolWcMTKNm88d5Xc5UU9BLyIRxTnHT/66kv1NLfzq2nzi1bI5YQp6EYkoC5dX8fqqbXz/kpMZM7S/3+XEBAW9iESMmr0H+GnRKvJz0vjmF0f7XU7MUNCLSMT46YJV7G1o5jfXTlTLJoQU9CISEV5dUcXfVlRx58VjGTtsgN/lxBQFvYj4bmd9Iz9ZsJLTs1L51nlq2YRad25qJiISVv9RtIq6/U08/c2zSYjX/DPU9I6KiK/eWLWVotJK5lwwllOGD/S7nJikoBcR39Tua+THf13J+BED+c4FeX6XE7PUuhER3/xs4Wp21TfyxM1nkaiWTdjonRURX7yzdhsvflLBd76Ux2mZqX6XE9MU9CLS4+r2N3Hviys4ZfgA5lw41u9yYp5aNyLS437+t9XU7G3kkRvPIilB881w0zssImH30KIyFpfVAPDep9v589Jyrjp9BP/0xiS8FPQiEnYTs1OZM7+Yt9ds494XV5CVlszf11UzMVu9+Z6g1o2IhN3UvHQemFHALU98TGNLgIHJCTx0w5lMzUv3u7ReQTN6EQm75pYAL3xSTmNLAICbpuYq5HuQgl5EwqqpJcCdz5bwUnEFfRPjuePCMTy9ZPPBnr2En1o3IhI2B5pbmDO/mDdXb6NfUjyPzCpkal46U/KGMGd+MXNnFmhm3wM0oxeRsGhoauFbf1zGm6u3ccn4YQdDHoI9+7kzC1heXudzlb2DZvQiEnL7GpuZ/dQy/llWwy++ejrXTx55xDFT89I1m+8hCnoRCam9B5q55YmPWbppJ7+5Np9/OTPb75J6PQW9iIRM3f4mbnr8I5aX1/H7GQVcnZ/pd0mCgl5EQqR2XyM3PPoRa7fu5sGZZ3DZhOF+lySeLk/GmtljZrbdzFa2GRtsZm+a2TrvcZA3bmb2gJmtN7PlZnZGOIsXkciwY+8Brv/fJXy6bQ8P33CmQj7CdOeqmyeAyw4buwd42zk3Fnjb2we4HBjrfc0G/hCaMkUkUm3f08CMeR+ysWYvj9xYyIWnDPO7JDlMl0HvnPs7sPOw4enAk972k8A1bcafckEfAmlmNiJUxYpIZKmq28+Mhz+konY/j980mfNOzvC7JOnA8V5HP8w5VwXgPQ71xrOALW2OK/fGRCTGbNm5j689/AHVew7wx1snc07eEL9LkqMI9QemrIMx1+GBZrPNbKmZLa2urg5xGSISTp/vqGfGvA+p29fE0988mzNPGux3SdKJ4w36ba0tGe9xuzdeDuS0OS4bqOzoBzjn5jnnCp1zhRkZ+ueeSLQoq97L1x7+gH2Nzfxp9hTyc9L8Lkm6cLxBXwTM8rZnAQvajN/oXX0zBahrbfGISPT7dOse/vXhD2kJOJ6dfY7Weo0SXV5Hb2Z/Ar4EpJtZOfBT4D7gz2Z2K7AZuM47/BXgCmA9sA+4OQw1i4gPVlbUccOjS0hKiOOZb57DmKH9/S5JuqnLoHfOXX+Upy7q4FgH3H6iRYlIZCndUssNjy5hQHIi8287m5OGpPhdkhwDfTJWRDq1dNNObnr8YwanJDH/trPJHtTP75LkGOk2xSJyVB+U7eDGxz5i6IA+PPetKQr5KKUZvYh06P111dz21FJyBvXjmdvOZuiAZL9LkuOkoBeRI7yzdhvffvoT8jL68/StkxnSv4/fJckJUNCLSDuvrdzKd//0CacMH8gfb51MWr8kv0uSE6SgF5GDXi6t5HvPlZCfncoTt0xmYHKi3yVJCCjoRQSAF5aVc/fzpRTmDuaxm86ifx/FQ6zQ/0kR4dmPNnPvSys4Ny+d/72xkL5J8X6XJCGkyytFermnPtjEPS+u4PyTM3hklkI+FmlGL9KLPfL+Bv7rb2u4ZPww5s4soE+CQj4WKehFeqkH313Pr1//lCtPH8H9MyaRGK9/4McqBb1IL+Oc43dvfsYD76znKwVZ/PraiSQo5GOagl6kF3HOcd9ra3l40Qb+tTCH//7q6cTHdbRekMQS/RoXiWEPLSpjcVkNEAz5ny1czcOLNnDGyDR+oZDvNRT0IjFsYnYqc+YX8891Nfz7X1fy+D83kZwQxw++PI44hXyvodaNSAza39jCyso6VlXsZuzQ/tzw2BICDpIT43hs1llMHZPud4nSgxT0IlEuEHBsqNlL8eZaSrYEv9Zu3UNLwAGQldaX0Rn9Wb99L7d9cbRCvhdS0ItEmeo9B7xA30XJllqWb6ljz4FmAAb0SWBiTirfPn80k3IGMSknjXXb9zBnfjF3XDiGp5ds5py8IUzNU9j3Jgp6kQjW0NTCyoo6SrbUUryllpLNtVTU7gcgPs44ZfgApk3KZFJOGgUj0xid3r9d731xWQ1z5hczd2YBU/PSmZI3pN2+9A4KepEIcXgLprS8lrVVe2hu04KZlJPGTVNzmTQyjQmZqV3ermB5eV27UJ+al87cmQUsL69T0PciFlzP21+FhYVu6dKlfpchEjIPLSpjYnZquzBdXFbD8vI6vn1+HgA1ew9Q0qavXlpey56G9i2YSTlpTMoZRH5OqlZ4kiOY2TLnXGFXx2lGLxIGrZc1ts6m3/t0O3f8qZhp+ZnMmf8JJVtqKd91qAUzbtgArs73WjA5aeRl9NfljxIyCnqRENp7oJlNNfXU7G3kwlOGctPjH5PaN4HqPY0APL1k88EWzKxzut+CETkRCnqRY9TQ1MLmnfvYUF3Pph31bKyuZ+OOejbW1FO950C7Y/v3iad6TyOFJw1i9nmjmTQyTS0Y6XEKepEONLUEKN+1n0019WyoqWdTTTDIN9bUU1m3n7anttL7JzEqPYULxmWQm57C6PQUctNT2FrbwF1/KeWOc0fx9JLN9E9OUMiLLxT0EnO6cyIUgle5VO1uOCLMN9XUs3nnvoNXuwAMSE5gdHoKZ+UOIjc9m1HpKYzyAr2jdVUXl9Vw119KdVmjRARddSMxp+214+eMHsJrK7fywxeWc/1ZI4mLs0OBvqOeA82Bg9+XnBhH7pAURmekkDsk5WCYj0pPYXBKEmbdPzna3V82Iieiu1fdKOglqjW1BNi2u4GtdQ1U1jVQVbufqroGVlbUUbyllniDxpZDf8YT442cwf0Y3WZG3hrmwwYk60oXiSq6vFKiXkvAUbP3AJVeeLc+VtXtp7I2+Fi95wCBw+YqA/okMDw1may0vmzeuY/zxqZzyxdGMSo9hay0vlpkQ3odBb2EzLG0K5xz7KhvpKq2gcq6/Qdn4m1n5dt2N7TrkwP0TYxnRFoymal9+eLYDDJTkxmR1pcRqclkeo8DkhMPtm9a7++SlBDHSUNSeuR9EIk0YQl6M7sM+D0QDzzinLsvlD8/UvqfqqO9gx8Sur6A8ZkDeXXFVn7+yhquK8zmV6+tbTcr31rXQGNLoN33JyXEMSI1mRGpyUweNTi4ndY3GOapfclMSya1b2KXvXLd30WkvZD36M0sHvgMuAQoBz4GrnfOrT7a9xxrj/7wv8iH73dXIOBocY6WgCPQ+hgguO1c++cDHHFs8eZd3PfqWu6+dBwTslJZUVHHb17/lB9cOo7Ts1K7XceJan3d7186jtNGDKR4yy4eeHs93z5vNHlD+9PY4mhqDtAcCBzcbmoJ0BxwNHrbwS/XbruxJUBzm/HGDr/P0dzi/dyWAA1NLe1OcLZKiDOGDUwmMy0Y2q2BHgzyvoxIS2bIMZ7wPJpI+cUnEm6+nYw1s3OA/3DOXert3wvgnPvF0b7neE7GLi6r4bYnl5KUEEfd/ibS+/chKSGuTTjTJrwPBbZzhwJbghLjjcT4OO+r8+2khOB+QpyRmBBHkvdcQvyh7eLNtSz9fBeXTxjO7PNGk5nWl/T+fbRsnUiI+XkyNgvY0ma/HDj78IPMbDYwG2DkyJHH/CJT89K54JShLFxexdhh/TltxEDizIiLM+JbH+M4tG1GfFxwO87ajwePbbNteMcdNh5HuzHztotKK3i5tIpp+Zl85Yys43zbTtxLn1RQVFrJVwuy+PqUkW1Cun1oJ8XHkZhgJMQFx0Mxi261uKyGFz6pONgbv+Gckxg2UB8SEvFTOIK+o9Q4YvrsnJsHzIPgjP5YX2RxWQ2Ly3YcDJSvnZXjS/91cVkN/1x/qI4Zk/2r4x/raw7WcW1hNmeeNLjHa1BvXCTyhOM6s3Igp81+NlAZyhdoGyh3fXkcc2cWMGd+8cHV7nuK6mivs3ufi4h/wtGjTyB4MvYioILgydiZzrlVR/ueY+3RR8rJNtUhIn7y9ZOxZnYFcD/Byysfc879vLPj9clYEZFj5+snY51zrwCvhONni4jIsdFnwUVEYpyCXkQkxinoRURinIJeRCTGRcT96M2sGvj8OL89HejZC8Yjm96P9vR+HKL3or1YeD9Ocs5ldHVQRAT9iTCzpd25vKi30PvRnt6PQ/RetNeb3g+1bkREYpyCXkQkxsVC0M/zu4AIo/ejPb0fh+i9aK/XvB9R36MXEZHOxcKMXkREOhHVQW9ml5nZp2a23szu8bsev5hZjpm9a2ZrzGyVmd3pd02RwMzizazYzBb6XYvfzCzNzJ43s7Xen5Nz/K7JL2b2f7y/JyvN7E9mFvMr40Rt0Htr0z4IXA6MB643s/H+VuWbZuD7zrlTgSnA7b34vWjrTmCN30VEiN8DrznnTgHy6aXvi5llAXcAhc65CQTvsDvD36rCL2qDHpgMrHfObXDONQLPAtN9rskXzrkq59wn3vYegn+J/VvTMAKYWTZwJfCI37X4zcwGAucBjwI45xqdc7X+VuWrBKCvt3ZGP0K8MFIkiuag72ht2l4dbgBmlgsUAEv8rcR39wM/BAJ+FxIBRgPVwONeK+sRM0vxuyg/OOcqgN8Am4EqoM4594a/VYVfNAd9t9am7U3MrD/wAvA959xuv+vxi5ldBWx3zi3zu5YIkQCcAfzBOVcA1AO98pyWmQ0i+C//UUAmkGJm3/C3qvCL5qAP+9q00cTMEgmG/DPOuRf9rsdn5wLTzGwTwZbehWb2tL8l+aocKHfOtf4r73mCwd8bXQxsdM5VO+eagBeBqT7XFHbRHPQfA2PNbJSZJRE8oVLkc02+MDMj2H9d45z7rd/1+M05d69zLts5l0vwz8U7zrmYn7UdjXNuK7DFzMZ5QxcBq30syU+bgSlm1s/7e3MRveDEdFiWEuwJzrlmM5sDvM6htWmPugB5jDsXuAFYYWYl3ti/eUs6igB8F3jGmxRtAG72uR5fOOeWmNnzwCcEr1Yrphd8QlafjBURiXHR3LoREZFuUNCLiMQ4Bb2ISIxT0IuIxDgFvYhIjFPQi4jEOAW9iEiMU9CLiMS4/w8pqnaIAAAABElEQVRGLWzlfeEHCQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1e864d99c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Pass a distribution of class rv_discrete that says how many children a member of generation i has\n",
    "# also pass the initial population, usually this is 1\n",
    "class Galton_Watson():\n",
    "    def __init__(self,dist,initial_pop=1):\n",
    "        self.distr=dist\n",
    "        self.population=initial_pop\n",
    "        return\n",
    "    def advance_generation(self):\n",
    "        num_new=self.distr.rvs(self.population).sum()\n",
    "        print(\"Population is now: %i\"%num_new)\n",
    "        self.population=num_new\n",
    "        return self.population\n",
    "    def advance_generation_individually(self):\n",
    "        num_new=0\n",
    "        for i in range(self.population):\n",
    "            current=self.distr.rvs()\n",
    "            print(\"Member %i\"%(i)+\" had this number of children: %i\"%current)\n",
    "            num_new+=current\n",
    "        print(\"Population is now: %i\"%num_new)\n",
    "        self.population=num_new\n",
    "        return self.population\n",
    "    def plot_generations(self,num_gens):\n",
    "        ys=list(range(num_gens))\n",
    "        for x in range(num_gens-1):\n",
    "            ys[x]=self.population\n",
    "            self.advance_generation()\n",
    "        ys[num_gens-1]=self.population\n",
    "        plt.plot(range(num_gens),ys,marker='x')\n",
    "        plt.show()\n",
    "        return\n",
    "        \n",
    "\n",
    "rv=scipy.stats.poisson(2)\n",
    "my_gw=Galton_Watson(rv,1)\n",
    "for i in range(5):\n",
    "    my_gw.advance_generation_individually()\n",
    "print(\"Now for plot\")\n",
    "my_gw=Galton_Watson(rv,1)\n",
    "my_gw.plot_generations(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}