{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import scipy" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEOCAYAAABB+oq7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xu8VXWd//HXW8QwEy/ACIE8YJQ0lCI7iY5Oo1aKZWIThToaMT7CSn/ZTNN4mcryMlG/ytE0kxLESi4/yiC1Icd0unnDJI+AjeQlDkEiIOgYKvT5/bG+Gze0zzlrn7PX3mef834+Hutx9v7s71rruxdwPnzX97IUEZiZmdXCbo2ugJmZ9R5OKmZmVjNOKmZmVjNOKmZmVjNOKmZmVjNOKmZmVjOFJxVJ/SQ9LOm29H60pPslrZI0X9IeKf6a9H5V+nxU2TEuTvHfSjqpLD4xxVZJuqjo72JmZh2rR0vlAmBl2fsvAVdFxMHAJuCcFD8H2JTiV6VySBoLnA4cBkwEvpESVT/gOuBkYCxwRiprZmYNUmhSkTQCeA/w7fRewAnAwlRkDnBaej0pvSd9/o5UfhIwLyJeiogngVXAkWlbFRFPRMTLwLxU1szMGqTolsp/AP8K/Dm9HwQ8FxHb0vs2YHh6PRxYDZA+35zK74jvsk97cTMza5DdizqwpFOAZyLiIUnHFXWenHWZDkwH2Guvvd566KGHNrI6ZmZNZfDgwSxZsmRJREzsrGxhSQU4BjhV0ruBAcBA4GpgX0m7p9bICGBNKr8GOBBok7Q7sA+woSxeUr5Pe/GdRMRMYCZAS0tLLF26tPvfzsysD5E0OE+5wm5/RcTFETEiIkaRdbT/NCL+AbgbmJyKTQUWpdeL03vS5z+NbLXLxcDpaXTYaGAM8ADwIDAmjSbbI51jcVHfx8zMOldkS6U9FwLzJF0BPAzcmOI3At+RtArYSJYkiIjlkhYAK4BtwHkRsR1A0vnAEqAfMCsiltf1m5iZ2U7U15a+9+0vM7PqSXooIlo6K9eIloqZWdN65ZVXaGtrY+vWrY2uSiEGDBjAiBEj6N+/f5f2d1IxM6tCW1sbe++9N6NGjSKbStd7RAQbNmygra2N0aNHd+kYXvvLzKwKW7duZdCgQb0uoQBIYtCgQd1qhTmpmJlVqTcmlJLufjcnFTOzJvTDH/4QSTz22GMAPPXUU+y5556MHz9+x/byyy9z0003MWTIEMaPH8+hhx7KVVddVWi9nFTMzJrQ3LlzOfbYY5k7d+6O2EEHHcSyZct2bHvssQcAU6ZMYdmyZfzyl7/kyiuvZPXq1e0dttucVMzMmswLL7zAL37xC2688UbmzZuXe79BgwZx8MEHs3bt2sLq5tFfZmZd9IUfLWfFH7bU9JhjXz+QS997WIdlFi1axMSJE3nDG97AoEGDeOihhxg0aBC/+93vGD9+PADHHHMM11133U77/f73v2fr1q286U1vqmmdyzmpmJk1mblz53LBBRcAcPrppzN37lzOP//8Hbe/djV//nx+9rOf8dhjj3HttdcyYMCAwurmpGJm1kWdtSiKsHHjRn7605/S2tqKJLZv344kzjvvvHb3mTJlCtdeey1Lly7lxBNP5NRTT2Xo0KGF1M99KmZmTWThwoWcffbZPP300zz11FOsXr2a0aNH5+p8b2lp4eyzz+bqq68urH5OKmZmTWTu3Lm8733v2yn2/ve/ny9+8Yu59r/wwguZPXs2zz//fBHV84KSZmbVWLlyJW984xsbXY1CVfqOeReUdEvFzMxqxknFzMxqxknFzMxqxknFzMxqxknFzMxqprCkImmApAck/UbScklfSPGbJD0paVnaxqe4JF0jaZWkRyQdUXasqZIeT9vUsvhbJbWmfa5Rb16P2sysCRTZUnkJOCEi3gyMByZKOip99umIGJ+20poCJwNj0jYduB5A0v7ApcAE4EjgUkn7pX2uBz5Stt/EAr+PmVmP0K9fv52WuJ8xYwYAo0aN4tlnn91R7p577uGUU04BqNsS+IUt0xLZBJgX0tv+aetoUswk4Oa0332S9pU0DDgOuDMiNgJIupMsQd0DDIyI+1L8ZuA04McFfB0zsx5jzz33rLjGV2dKy7Vs2LCBQw45hMmTJ3PggQfWtG6F9qlI6idpGfAMWWK4P310ZbrFdZWk16TYcKB8nYG2FOso3lYhbmZmHShyCfxCF5SMiO3AeEn7ArdKOhy4GFgH7AHMBC4ELiuyHpKmk91SY+TIkUWeysz6kh9fBOtaa3vMoePg5BkdFvnTn/60Y4l7gIsvvpgpU6bkPkWRS+DXZZXiiHhO0t3AxIj4Sgq/JGk28C/p/RqgvB02IsXWkN0CK4/fk+IjKpSvdP6ZZAmMlpaWvrUujZn1Ou3d/qo0Vqk8Vo8l8AtLKpKGAK+khLIn8C7gS5KGRcTaNFLrNODRtMti4HxJ88g65TenckuAfy/rnD8RuDgiNkrakjr/7wc+BHy9qO9jZvYXOmlR1NugQYPYtGkTgwcPBrJl8kuvoT5L4BfZpzIMuFvSI8CDZH0qtwHfk9QKtAKDgStS+TuAJ4BVwLeAjwOkDvrL0zEeBC4rddqnMt9O+/wOd9KbWR923HHH8Z3vfAeA7du3893vfpfjjz/+L8oVuQR+kaO/HgHeUiF+QjvlA6j4lJmImAXMqhBfChzevZqamTWXXftUJk6cyIwZM/jsZz/Lxz72Md785jcTEUycOJGzzjqr4jEuvPBCjjjiCC655BL23nvvmtXNS9+bmVXBS993zMu0mJlZzTipmJlZzTipmJlZzTipmJlVqTf3RXf3uzmpmJlVYcCAAWzYsKFXJpaIYMOGDd2aFFmXGfVmZr3FiBEjaGtrY/369Y2uSiEGDBjAiBEjOi/YDicVM7Mq9O/fn9GjRze6Gj2Wb3+ZmVnNOKmYmVnNOKmYmVnNOKmYmVnNOKmYmVnNOKmYmVnNeEixmfVpt9z/exYtq/jQ2E5NGj+cMyf4EeXl3FIxsz5t0bI1rFi7per9Vqzd0uVk1Ju5pWJmfd7YYQOZf+7RVe0z5YZ7C6pNc+u0pSLpAEk3Svpxej9W0jnFV83MzJpNnttfNwFLgNen9/8DfLKoCpmZWfPKk1QGR8QC4M8AEbEN2N7ZTpIGSHpA0m8kLZf0hRQfLel+SaskzZe0R4q/Jr1flT4fVXasi1P8t5JOKotPTLFVki6q6pubmVnN5elT+V9Jg4AAkHQUsDnHfi8BJ0TEC5L6A79It9D+GbgqIuZJ+iZwDnB9+rkpIg6WdDrwJWCKpLHA6cBhZK2l/5L0hnSO64B3AW3Ag5IWR8SKfF/dzKx7Vqzd0qW+ld48aixPUvkUsBg4SNIvgSHABzrbKbKHDbyQ3vZPWwAnAGem+Bzg82RJZVJ6DbAQuFaSUnxeRLwEPClpFXBkKrcqIp4AkDQvlXVSsZrqzpDTkt78S6SvmjR+eJf2K400661/HzpNKhHxkKS/Aw4BBPw2Il7Jc3BJ/YCHgIPJWhW/A55Lt9Aga2GU/mSGA6vTObdJ2gwMSvH7yg5bvs/qXeIT2qnHdGA6wMiRvfMP0opTGnI6dtjALu3f23+J9FVnThjZpT/TKTfc2+UWTneMff1ALn3vYYWfp9OkIul3wP+NiG+WxW6LiFM62zcitgPjJe0L3Aoc2p3KdlVEzARmArS0tPS+x7VZ4boy5LSku79E3MrpXbrawmkWeW5/vQIcL2kCcG5EvMyrLYVcIuI5SXcDRwP7Sto9tVZGAKX7CmuAA4E2SbsD+wAbyuIl5fu0FzfrMbrzS8StnN6nqy2cZpEnqbwYEVMk/Svwc0kfIHXad0TSEOCVlFD2JOtQ/xJwNzAZmAdMBRalXRan9/emz38aESFpMXCLpK+RddSPAR4guxU3RtJosmRyOq/21Zj1GN35JeIJdtZs8iQVAUTElyX9GvgJsH+O/YYBc1K/ym7Agoi4TdIKYJ6kK4CHgRtT+RuB76SO+I1kSYKIWC5pAVkH/DbgvHRbDUnnk82h6QfMiojleb60mZkVI09S+VzpRUT8V5onMrWznSLiEeAtFeJP8OrorfL4VtoZVRYRVwJXVojfAdzRWV3MrA9ZOhtaF+Yu/rkNaYbE7H1g3GRomVZQxfqGdic/Sip1qq+RdERpIxuRdVtdamdmVq3WhbCutfr91rVWlYysso5aKv9MNgz3qxU+K803MbNqdOd/0eD/Sec1dBxMuz1X0ctSv9X8Pa4oskZ9RrtJJSKmp5/H1686Zr1c6X/RQ8dVv2/pf99OKtaD5Zmn8gHgPyPieUmfAY4ALo+IhwuvnVlv1JX/RU87Gma/p8hamdVEngUlP5sSyrHAO8lGaX2zk33MzKwPypNUSisSvweYGRG3A3sUVyUzM2tWeZLKGkk3AFOAOyS9Jud+ZmbWx+RJDh8km2B4UkQ8Rzbx8dOF1srMzJpSnlWKXwR+UPZ+LbC2yEqZmVlzyjOj3szKVTnXZCddHU5s9bGutWuj7Dx/aAf3jZhVq6sztiFLKOMm17Y+VhvjJnd9/pBn4u/glopZV1Qx18SaRMu0rrU2PH9oJ522VCT9vaTHJW2WtEXS85K21KNyZmbWXPK0VL4MvDciVhZdGTMza255ksofnVDMrCe75f7fs2hZ9uDX0iKcl+V8wNmKtVsYO2xgYXXra/IklaWS5gM/BF4qBSPiB+3vYmZWP4uWrelychg7bGCvf258PeVJKgOBF4ETy2JB2dwVM7NGGztsIPPPPXrHYwLmTzu6wTXqm/JMfvTgazMzyyXP0vcjgK8Dx6TQz4ELIqKtk/0OBG4GDiBr2cyMiKslfR74CLA+Fb0kPRYYSRcD55AtYvmJiFiS4hOBq8meRf/tiJiR4qOBeWRPo3wIODsiXs731a0vKL/X3lW+526WX57Jj7OBxcDr0/ajFOvMNuBTETEWOAo4T9LY9NlVETE+baWEMhY4HTgMmAh8Q1I/Sf2A64CTgbHAGWXH+VI61sHAJrKEZLZD6V57d/ieu1l+efpUhkREeRK5SdInO9upfI2w9DyWlUBH/zInAfMi4iXgSUmrgCPTZ6si4gkASfOASel4JwBnpjJzgM8D1+f4TtaH7LjXbmaFy9NS2SDprFKrQdJZwIZqTiJpFPAW4P4UOl/SI5JmSdovxYYDq8t2a0ux9uKDgOciYtsu8Urnny5pqaSl69evr1TEzMxqIE9S+Uey5e/XkbU8JgO5O+8lvQ74PvDJiNhC1pI4CBifjvfVKutctYiYGREtEdEyZMiQok9nZtZn5Rn99TRwalcOLqk/WUL5XmleS0T8sezzbwG3pbdrgAPLdh+RYrQT3wDsK2n31FopL29mZg3QblKR9K8R8WVJXycbvbWTiPhERweWJLLn2a+MiK+VxYel/haA9wGPpteLgVskfY1sQMAY4AFAwJg00msNWWf+mRERku4maznNA6YCi3J8ZzMzK0hHLZXS0ixLu3jsY4CzgVZJy1LsErLRW+PJEtVTwLkAEbFc0gJgBdnIsfMiYjuApPPJnj7ZD5gVEcvT8S4E5km6AniYLImZmVmDtJtUIuJH6eWLEfH/yj+T9IHODhwRvyBrZezqjg72uRK4skL8jkr7pRFhR+4aNzOzxsjTUX9xzpiZmfVxHfWpnAy8Gxgu6ZqyjwaS3Z4yszpYsXYLU264t+rVdwEmjR/OmRNGFlU1s7/QUZ/KH8j6U04lWwKl5Hngn4qslJllujOTv7SSgJOK1VNHfSq/AX4j6ZaIeKWOdTKz5MwJI19NClWuvjulihaNddO61vo/VnjoODh5Rn3PmUOeZVpGSfoi2bpbA0rBiPjrwmplZtYsxk1udA16lDxJZTZwKXAVcDzZbPo8HfxmZr1fy7RsMyBfctgzIu4CFBFPR8TngTq388zMrBnkaam8JGk34PE0CXEN8Lpiq2VmZs0oT1K5AHgt8AngcrLl5qcWWSkza0cVHcKlIcilDn7GTfZtGitcngUlH0wvX6CK1YnNrMa60yG8rjX76aRiBeto8uOPqLCQZElEdGnlYjProio7hEuTJOdPO7r+w12tz+qopfKVutXCzMx6hY4mP/53PStiZmbNr9M+FUlPUvl5Kp78aGZmO8kz+qul7PUA4APA/sVUx8zMmlmnkx8jYkPZtiYi/gNPfjQzswry3P46ouztbmQtlzwtnJ7p2ce7PhLG4/zNzDqUJzl8tez1NrJHAH+wkNr0ZB7nb2bWqTyTH4/vyoElHQjcDBxA1tE/MyKulrQ/MB8YRUpQEbFJkoCryR4M9iLw4Yj4dTrWVOAz6dBXRMScFH8rcBOwJ9njhi+IiHbn1gAweAxMu736L+Rx/mZmneq0T0XSIEnXSPq1pIckXS1pUI5jbwM+FRFjgaOA8ySNBS4C7oqIMcBd6T3AycCYtE0Hrk/n359sleQJZM+jv1TSfmmf64GPlO03Mc+XNjOzYuRZpXgesB54PzA5vZ7f2U4RsbbU0oiI54GVwHBgEjAnFZsDnJZeTwJujsx9wL6ShgEnAXdGxMaI2ATcCUxMnw2MiPtS6+TmsmOZmVkD5EkqwyLi8oh4Mm1XkN3Syk3SKOAtwP3AARGxNn20ruxYw4HVZbu1pVhH8bYK8Urnny5pqaSl69evr6bqZmZWhTxJ5SeSTpe0W9o+CCzJewJJrwO+D3wyIraUf5ZaGB33gdRARMyMiJaIaBkyZEjRpzMz67PyJJWPALcAL6dtHnCupOclbeloR0n9yRLK9yLiByn8x3TrivTzmRRfAxxYtvuIFOsoPqJC3MzMGiTP5Me9I2K3iNg9bbul2N4RMbC9/dJorhuBlRHxtbKPFvPq81imAovK4h9S5ihgc7pNtgQ4UdJ+qYP+RGBJ+myLpKPSuT5UdiwzM2uAXJMYJZ0KvD29vScibsux2zHA2UCrpGUpdgkwA1gg6RzgaV6d83IH2XDiVWRDiqcBRMRGSZcDpee6XBYRG9Prj/PqkOIfp83MzBokz4z6GcDbgO+l0AWSjomIizvaLyJ+Aaidj99RoXwA57VzrFnArArxpcDhHdXD7C8snQ2tC7u+/7pWGDqudvUx60XytFTeDYyPiD8DSJoDPAx0mFSsZ7nl/t+zaFmxXU6Txg/nzAkjCz1HTbQu7F5iGDque09hNOvF8q7htS9QuuW0T0F1sQItWraGFWu3MHZYu91g3bJibTZmoymSCmSJoSsrK5hZh/IklS8CD0u6m+x21tt5dRa8NZGxwwYy/9yjCzn2lPToWjPr2/Ks/TVX0j1k/SoBXBgR64qumJmZNZ+8t7+OBo4lSyq7A7cWViMzM2taeUZ/fQM4GJibQudKemdEVBypZWbWbVWO0Pvchs3Zi9n7eHReg+VpqZwAvLG0pHwa/bW80FqZWd/WnRF6Hp3XUHmSyipgJNlERciWTFlVWI3MzKCqEXqXpYEi86cVMxDF8suTVPYGVkp6gKxP5UhgqaTFABFxaoH1MzOzJpInqXyu8FqYmVmvkGdI8X/XoyJmVnsr1m5hyg337ujIvqyK+URNs0KC9Sh5hxSbWZOZNL7iM+tyaboVEqzHcFIx66XOnDDy1aQwO1tdKW9HtldIsK7K85AuJO0p6ZCiK2NmZs0tz+TH9wJfAfYARksaT/ZME4/6MrOaKq2mXW0fUJGLpVp18rRUPk82jPg5gIhYBowusE5m1keVVtOu1thhA7vVh2S1k6dP5ZWI2Jw9sXeHKKg+ZlaUda0w+z25iu607AlkM9RbphVUsZ2NHTaQw/aorg/Ieo48SWW5pDOBfpLGAJ8AflVstcysprqzbMm61uxnnZKKNbc8SeX/AP8GvATcAiwBruhsJ0mzgFOAZyLi8BT7PPARYH0qdklE3JE+uxg4B9gOfCIilqT4ROBqoB/w7YiYkeKjgXnAIOAh4OyIeDnH9zHre1qmVZUUdlr2JGfrxgzy9akcGhH/FhFvS9tnImJrjv1uAiZWiF8VEePTVkooY4HTgcPSPt+Q1E9SP+A64GRgLHBGKgvwpXSsg4FNZAnJzMwaKE9S+aqklZIul3R43gNHxM949RHEnZkEzIuIlyLiSbIFK49M26qIeCK1QuYBk5R18JwAlNbGngOclrduZmZWjE6TSkQcDxxPdsvqBkmtkj7TjXOeL+kRSbMk7Zdiw4HVZWXaUqy9+CDguYjYtku8IknTJS2VtHT9+vXtFTMzs27KNfkxItZFxDXAR4FldH2RyeuBg4DxwFrgq108TlUiYmZEtEREy5AhQ+pxSjOzPinP5Mc3AlOA9wMbgPnAp7pysoj4Y9lxvwXclt6uIXtOS8mIFKOd+AZgX0m7p9ZKefkeoTSJq6fw5DAzq4c8o79mkSWSkyLiD905maRhEbE2vX0f8Gh6vRi4RdLXgNcDY4AHAAFj0kivNWSd+WdGREi6G5hM1s8yFVjUnbrVWmkSV0/5Rd5sk8NqlZR70p9Bn9TVxwLr934kcJPKs/R9l2YfSZoLHAcMltQGXAocl5Z5CeAp4Nx0juWSFgArgG3AeRGxPR3nfLJhzP2AWRFRepTxhcA8SVcADwM3dqWeRRo7bCDzz/Xkra6oVVJutmTak9Rk2fyuPhbYjwRuWu0mFUkLIuKDklrZeQa9gIiIN3V04Ig4o0K43V/8EXElcGWF+B3AHRXiT5CNDrNeykm5cWq6bL4fC9yndNRSuSD9PKUeFTGznsPL5ltXtTv6q6zv4+MR8XT5Bny8PtUzM7NmkmdI8bsqxE6udUXMzKz5ddSn8jGyFslBkh4p+2hvvKCkWd/S1RWOu9JJb02toz6VW4AfA18ELiqLPx8ReZdf6XGeWP+/Xbrn64cGWZ/VhVFYL768neVrNwMj+eVzR3CX/930Ge0mlYjYDGyWdDWwMSKeB5A0UNKEiLi/XpXsKUa98gSf2/DpfIX3gMEvvQZmDyiuQnV8xoX1YVWucPybbswx8hDw5pdn8uP1wBFl71+oEGsafz1kr64NU116DrQu5LDaV6lr/IwL66F2GjlmfU6epKKI2DFPJSL+LCnPfr1Llf9bK5yfcWFmPVCe0V9PSPqEpP5puwB4ouiKmZlZ88mTVD4K/A3Z2lttwARgepGVMjOz5pRn7a9nyBZyNDMz61CnLRVJb5B0l6RH0/s3dfMhXWZm1kvluf31LeBi4BWAiHgEt1zMzKyCPKO4XhsRD2SPhd9hW3uFre8qLZVeq2N5EpxZ88mTVJ6VdBBp+XtJk8keBWy2Q60nrHkSnFlzypNUzgNmAodKWgM8CfxDobWypuMJb2YGHS8oeUFEXA0Mi4h3StoL2K20XIuZmdmuOuqoL00f/zpARPyvE4qZmXWko6SyUtLjZLe9HinbWndZCr8iSbMkPVMaipxi+0u6U9Lj6ed+KS5J10halc5xRNk+U1P5xyVNLYu/NdVlVdpXmJlZQ3X05MczgL8FHgfeW7adkn525iZg4i6xi4C7ImIMcBevLql/MjAmbdPJFqxE0v7ApWSz+I8ELi0lolTmI2X77XouMzOrs3aTiqS7ImIdsGTXxwmnRwp3KCJ+Buz63JVJwJz0eg5wWln85sjcB+wraRhwEnBnRGyMiE3AncDE9NnAiLgvLXZ5c9mxzMysQToa/TVM0t8A75U0F9jp9lJE/LoL5zsgIkrDkdcBB6TXw4HVZeXaUqyjeFuFuJmZNVBHSeVzwGeBEcDXdvksgBO6c+KICEnRecnukzSdtAjmyJEe9mpmVpSOnvy4EFgo6bMRcXmNzvdHScMiYm26hfVMiq8BDiwrNyLF1gDH7RK/J8VHVChfUUTMJJtrQ0tLS10SmZlZX5Rn7a8rJZ0l6XMAkkZKOrKL51sMlEZwTQUWlcU/lEaBHQVsTrfJlgAnStovddCfSNbHsxbYIumoNOrrQ2XHMjOzBskzo/464M9kt7suA54Hvg+8raOdUj/MccBgSW1ko7hmAAsknQM8DXwwFb8DeDewCniRNEcmIjZKuhx4MJW7LCJKnf8fJxthtifw47SZmVkD5UkqEyLiCEkPA0TEJkl7dLZTGpJcyTsqlA2y5WAqHWcWMKtCfClweGf1MDOz+slz++sVSf14dUHJIWQtFzMzs53kSSrXALcCfyXpSuAXwL8XWiszM2tKeR4n/D1JD5HdthJwWkSsLLxmZmbWdPL0qRARjwGPFVwXMzNrcnluf5mZmeWSq6Vi1mMsnQ2tC7t3jHWtMHRcbepjZjtxS8WaS+vCLCl0x9BxMG5ybepjZjtxS8Waz9BxMO32RtfCzCpwS8XMzGrGScXMzGrGScXMzGrGScXMzGrGScXMzGrGScXMzGrGScXMzGrGScXMzGrGScXMzGrGM+qb2bpWmP2e4o4/bjK0TCvu+GbW6zSkpSLpKUmtkpZJWppi+0u6U9Lj6ed+KS5J10haJekRSUeUHWdqKv+4pKmN+C4NM25ysYsirmvt/sKNZtbnNLKlcnxEPFv2/iLgroiYIemi9P5C4GRgTNomANcDEyTtD1wKtJA96vghSYsjYlM9v0TDtEwrthVRZAvIzHqtntSnMgmYk17PAU4ri98cmfuAfSUNA04C7oyIjSmR3AlMrHelzczsVY1KKgH8RNJDkqan2AERsTa9XgcckF4PB1aX7duWYu3FzcysQRp1++vYiFgj6a+AOyXt9KjiiAhJUauTpcQ1HWDkyJG1OqyZme2iIS2ViFiTfj4D3AocCfwx3dYi/XwmFV8DHFi2+4gUay9e6XwzI6IlIlqGDBlSy69iZmZl6p5UJO0lae/Sa+BE4FFgMVAawTUVWJReLwY+lEaBHQVsTrfJlgAnStovjRQ7McXMzKxBGnH76wDgVkml898SEf8p6UFggaRzgKeBD6bydwDvBlYBLwLTACJio6TLgQdTucsiYmP9voaZme2q7kklIp4A3lwhvgF4R4V4AOe1c6xZwKxa19HMzLrGM+qtPpbOrs1kynWtxU76NLNu6UnzVKw3a12YJYTuGjouW03AzHokt1SsfoaOg2m3N7oWZlYgJxVrXy0XrPRtK7M+wUnFKqv1LSbftjLrE5xUrLKiF6w0s17JHfVmZlYzTipmZlYzTipmZlYzTipmZlYzTipmZlYzTipmZlYzTipz3dVrAAAGXUlEQVRmZlYzTipmZlYzTipmZlYzTipmZlYzTipmZlYzTipmZlYzTipmZlYzTZ9UJE2U9FtJqyRd1Oj6mJn1ZU2dVCT1A64DTgbGAmdIGtvYWpmZ9V1NnVSAI4FVEfFERLwMzAMmNbhOZmZ9VrM/pGs4sLrsfRswYddCkqYD09PbrZKWd3DMfYDNXfi8UjxPbDDwbAfnq7XOvl+t989TvqMy1VzvSvFK5ep5zbt7vas9RlHXu73P/He8vn/HG3W98x8vIpp2AyYD3y57fzZwbSf7zCzi80rxPDFgaZ2vWYffr9b75ynfUZlqrnc717fSn0Hdrnl3r3e1xyjqendwLf13vI5/x3vi9d51a/bbX2uAA8vej0ixjvyooM8rxfPG6qm75692/zzlOypTzfWuFG/2613tMYq63u195r/j9f073hOv906UMltTkrQ78D/AO8iSyYPAmRHR0e2tHkXS0ohoaXQ9+hJf8/ry9a6vRl/vpu5TiYhtks4HlgD9gFnNlFCSmY2uQB/ka15fvt711dDr3dQtFTMz61mavU/FzMx6ECcVMzOrGScVMzOrGSeVHkjSXpKWSjql0XXp7SQdJ+nnkr4p6bhG16cvkLSbpCslfV3S1EbXp7eT9Lfp7/e3Jf2q6PM5qdSBpFmSnpH06C7x9hbDvBBYUN9a9h5VXu8AXgAGkK3IYF1Q5TWfRDan7BV8zbukmusdET+PiI8CtwFzCq+bR38VT9LbyX5x3RwRh6dYP7I5Nu8i+4f1IHAG2dIzg8h+yT0bEbc1pNJNrMrr/VhE/FnSAcDXIuIfGlTtplblNT8V2BQRN0haGBGTG1TtplXN9Y6IFenzBcA5EfF8kXVr6nkqzSIifiZp1C7hHYthAkgqLYb5OmAvslWX/yTpjoj4cx2r2/Squd6lf3DAJuA1datkL1Pl3/HVwMupzPZ61bE3qfJ6r5A0EthcdEIBJ5VGqrgYZkScDyDpw2QtFSeU2qh4vSX9PXASsC9wbSMq1ou1t+Dr1cDXJf0t8LNGVKyX6miB3XOA2fWohJNKDxURNzW6Dn1BRPwA+EGj69GXRMSLZL/krE4i4tJ6ncsd9Y3TlcUwret8vevP17y+esT1dlJpnAeBMZJGS9oDOB1Y3OA69Wa+3vXna15fPeJ6O6nUgaS5wL3AIZLaJJ0TEduA0mKYK4EFTbgYZo/k611/vub11ZOvt4cUm5lZzbilYmZmNeOkYmZmNeOkYmZmNeOkYmZmNeOkYmZmNeOkYmZmNeOkYmZmNeOkYtZgaclys17BScUsJ0mXSfpk2fsrJV2QXn9a0oOSHpH0hbIyP5T0kKTlkqaXxV+Q9FVJvwGOljRD0oq0/1cqnPtISfdKeljSryQdkuKvlbQg7XurpPsltaTPzpDUKulRSV8q8NKY7eAZ9WY5pedX/CAijpC0G/A42TMs3gpMBs4FRLbe0pfTMy/2j4iNkvYkW5vp7yJig6QApkTEAkmDgF8Bh0ZESNo3Ip7b5dwDgRcjYpukdwIfi4j3S/oXYExEnCvpcGAZcBTwB+C+VLdNwE+AayLih4VeJOvzvPS9WU4R8ZSkDZLeAhwAPJwSxInAicDDqejrgDFkzwr5hKT3pfiBKb6B7OFU30/xzcBW4EZJt5E99nVX+wBzJI0hewRy/xQ/luz5JETEo5IeSfG3AfdExHoASd8D3g44qVihnFTMqvNt4MPAUGBWign4YkTcUF5Q0nHAO4GjI+JFSfeQPSYaYGtEbAdIrY8jgXeQtXjOB07Y5byXA3dHxPtSi+meWn4ps1pxn4pZdW4FJpK1BJak2BLgHyW9DkDScEl/Rda62JQSyqFkt6X+Qtpvn4i4A/gn4M0Viu3Dq8/G+HBZ/JfAB9NxxgLjUvwB4O8kDU4DAc4A/rv6r2tWHbdUzKoQES9Luht4rqyl8RNJbwTulQTwAnAW8J/ARyWtBH5L1sdRyd7AIkkDyFo9/1yhzJfJbn99Bri9LP6NFF8BPAYsJ3sW+VpJFwF3p2PeHhGLuvPdzfJwR71ZFVIH/a+BD0TE4z2gPv2A/hGxVdJBwH8Bh0TEyw2umvVRbqmY5ZRuL90G3NoTEkryWuBuSf3JWiQfd0KxRnJLxczMasYd9WZmVjNOKmZmVjNOKmZmVjNOKmZmVjNOKmZmVjNOKmZmVjP/HyEd8swMBhWOAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ " mu = 1.25e-8\n", " gen = 30\n", " afrDat = pd.read_csv(\"/home/training/share/MSMC-tutorial-files/results/AFR.msmc2.final.txt\", delim_whitespace=True)\n", " eurDat = pd.read_csv(\"/home/training/share/MSMC-tutorial-files/results/EUR.msmc2.final.txt\", delim_whitespace=True)\n", " plt.step(afrDat[\"left_time_boundary\"]/mu*gen, (1/afrDat[\"lambda\"])/(2*mu), label=\"AFR\")\n", " plt.step(eurDat[\"left_time_boundary\"]/mu*gen, (1/eurDat[\"lambda\"])/(2*mu), label=\"EUR\")\n", " plt.ylim(0,40000)\n", " plt.xlabel(\"years ago\");\n", " plt.ylabel(\"effective population size\");\n", " plt.gca().set_xscale('log')\n", " plt.legend()" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEOCAYAAABmVAtTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAGUBJREFUeJzt3Xu0XlV57/HvT0C5SEAlahqIcmwK5uClGG7KOdICCq2HDK0iSD1HiwRUrByrlR4tKNpRtXrOoILKRUqxctOiRIiAMLh4QUxQihCLIiIk4OAiEAHl5nP+eFdeXmOy99rJXvvdO/v7GWOP7DXXXHM9Oyt5nz3XXGvOVBWSJAE8ZdgBSJImD5OCJKnPpCBJ6jMpSJL6TAqSpD6TgiSpr7OkkOS0JHcluWEt+w9Jcn2SHyb5TpKXdBWLJKmdLnsKpwP7jbD/Z8Arq+pFwEeAkzuMRZLUwsZdNVxVVyV5/gj7vzOw+V1g265ikSS1M1nGFA4Fvj7sICRpuuusp9BWkj+hlxT2HKHOQmAhwBZbbPGyHXfccYKik6QNw7XXXntPVc0crd5Qk0KSFwOnAvtX1b1rq1dVJ9OMOcyfP7+WLl06QRFK0oYhyc/b1Bva7aMkc4DzgDdX1Y+HFYck6Umd9RSSnAXsBWyTZDlwLLAJQFV9DjgGeBbwmSQAj1fV/K7ikSSNrsunjw4eZf/bgLd1dX5J0thNlqePJEmTgElBktRnUpAk9ZkUJEl9JgVJUp9JQZLUZ1KQJPWZFCRJfSYFSVLf0GdJlSSN7MxrbuP861ZMyLnsKUjSJHf+dStYdufKCTmXPQVJmgLmzZrBOYfvsc7Hn3tEu3r2FCRJfSYFSVKfSUGS1GdSkCT1mRQkSX0mBUlSn0lBktRnUpAk9ZkUJEl9vtEsSeNsvOcqWnbnSubNmjFu7Y3EnoIkjbPxnqto3qwZLHjp7HFrbyT2FCSpA+s7V9Gw2FOQJPWZFCRJfd4+kjRlTORiM+tjIgeGx1tnPYUkpyW5K8kNa9mfJP+c5OYk1yfZuatYJG0YJnKxmfUxkQPD423UnkKSzYG/AeZU1WFJ5gI7VNUFoxx6OnACcMZa9u8PzG2+dgM+2/wpSWs1VQdwp4o2PYV/AR4BVl2FFcBHRzuoqq4CfjlClQXAGdXzXWDrJLNaxCNJ6kibpPCCqvoE8BhAVT0MZBzOPRu4fWB7eVMmSRqSNknh0SSbAQWQ5AX0eg4TJsnCJEuTLL377rsn8tSSNK20SQofAi4CtkvyReAy4P3jcO4VwHYD29s2Zb+nqk6uqvlVNX/mzJnjcGpJ0pqMOtBcVZckuRbYnd5to3dX1T3jcO5FwJFJzqY3wPxAVd05Du1KktZRm6ePLquqvYEL11A20nFnAXsB2yRZDhwLbAJQVZ8DFgN/BtwMPAy8dR1/BknSOFlrUkiyKbA5vQ/1Z/Dk4PIMWgwIV9XBo+wv4J3tQ5UkdW2knsLhwFHAHwDX8mRSWEnv/QNJ0gZmrUmhqo4Hjk/yrqr69ATGJEkakjYDzZ9OshMwD9h0oHxtbypLkqaoNgPNx9IbMJ5Hb3B4f+BbrH36CknSFNXmPYXXA3sDv6iqtwIvAbbqNCpJ0lC0SQq/rqrfAo8nmQHcxe++dCZJ2kC0WU9haZKtgVPoPYX0IHB1p1FJkoZixKSQJMA/VtX9wOeSXATMqKrrJyQ6SdKEGjEpVFUlWQy8qNm+dSKCkiQNR5sxhe8n2aXzSCRJQ9dmTGE34JAkPwceovdmc1XVizuNTJI04dokhVd3HoUkaVJo80bzzyciEEnS8LUZU5AkTRMmBUlSX6ukkOR5SfZpvt8syZbdhiVJGoZRk0KSw4AvAyc1RdsCX+0yKEnScLTpKbwTeAW9xXWoqp8Az+4yKEnScLRJCo9U1aOrNpJsDFR3IUmShqVNUrgyyf8BNkuyL/Al4GvdhiVJGoY2L68dDRwK/JDeus2LgVO7DErS1HXmNbdx/nUrOml72Z0rmTdrRidtq6dNUtgMOK2qTgFIslFT9nCXgUmams6/bkVnH97zZs1gwUtnj3u7elKbpHAZsA+9dRSglxAuAV7eVVCSprZ5s2ZwzuF7DDsMrYM2YwqbVtWqhEDz/ebdhSRJGpY2SeGhJDuv2kjyMuDX3YUkSRqWNrePjgK+lOQOetNmPxd4Y6dRSZKGos0sqUuS7Ajs0BTdVFWPdRuWJGkY2vQUAHYBnt/U3zkJVXVGZ1FJkoaizdxHXwA+CexJLznsAsxv03iS/ZLclOTmJEevYf+cJJcn+UGS65P82RjjlySNozY9hfnAvKoa09QWzfsMJwL7AsuBJUkWVdWygWofBM6tqs8mmUfvxbjnj+U8kqTx0+bpoxvoDS6P1a7AzVV1SzN30tnAgtXqFLDqDZetgDvW4TySpHHSpqewDbAsyfeAR1YVVtUBoxw3G7h9YHs5sNtqdT4EXJLkXcAW9F6S+z1JFgILAebMmdMiZEnSumiTFD7U4fkPBk6vqk8l2QP4QpKdquq3g5Wq6mTgZID58+c7Q6skdaTNI6lXJnkeMLeqLk2yObBRi7ZXANsNbG/blA06FNivOc/VSTal1zO5q03wkqTxtS4rr82m3cprS4C5SbZP8lTgIGDRanVuA/ZuzvNCYFPg7nahS5LGW2crr1XV48CRwMXAj+g9ZXRjkuOSrBqP+BvgsCT/AZwFvGWsTzlJksZPmzGFR6rq0STA2FZeq6rF9B4zHSw7ZuD7ZfQSjiRpEnDlNUlSX5ukcDS9+/yDK699sMugJEnD4cprkqS+Nj2Fy+glgVU2Ay7tJhxJ0jC58pokqc+V1yRJfa68Jknqc+U1SVLfqEkhyRuAi6rqhiQfpLfy2ker6vvdhydpPJx5zW2cf93qU491Y9mdK5k3a8boFTUptRlT+Puq+lWSPenNU/R54LPdhiVpPJ1/3QqW3blyQs41b9YMFrx09oScS+OvzZjCE82ffw6cUlUXJvlohzFJ6sC8WTM45/A9hh2GJrk2PYUVSU6iN7i8OMnTWh4nSZpi2ny4H0hvptNXV9X9wDOB93UalSRpKEZNClX1ML1Fb/Zsih4HftJlUJKk4WizyM6xwPuBv2uKNgH+rcugJEnD0eb20WuBA4CHAKrqDmDLLoOSJA1Hm6TwaLMaWgEk2aLbkCRJw9ImKZzbPH20dbNe86XAKd2GJUkahjbTXHyyWXFtJb2pLo6pqm90HpkkacK1eXmNJgmYCCRpA7fWpJDkVzTjCKvvAqqqnNxEkjYwa00KVeUTRpI0zbS6fQSQ5NnApqu2q+q2TiKSJA1Nm5fXDkjyE+BnwJXArcDXO45LkjQEbR5J/QiwO/Djqtqe3vTZ3+00KknSULRJCo9V1b3AU5I8paouB+Z3HJckaQjajCncn+TpwDeBLya5i2bKC0nShqVNT2EB8GvgKOAi4KfA/2jTeJL9ktyU5OYkR6+lzoFJliW5McmZbQOXJI2/Nm80P5TkOcAuwL3A15vbSSNKshFwIrAvsBxYkmRRVS0bqDOX3uyrr6iq+5onnCRJQ9Lm6aMDge8Bb6C34M41SV7fou1dgZur6paqehQ4m16vY9BhwIlVdR9AVd01luAlSeOrzZjCB4BdVn1gJ5lJb1K8L49y3Gzg9oHt5cBuq9X5o6bNbwMbAR+qqotWbyjJQmAhwJw5c1qELElaF23GFJ6y2m/w97Y8ro2NgbnAXsDBwClJtl69UlWdXFXzq2r+zJkzx+nUkqTVtekpXJTkYuCsZvuNtHt5bQWw3cD2tk3ZoOXANVX1GPCzJD+mlySWtGhfkjTO2qzR/D7gJODFzdfJVfW3LdpeAsxNsn2SpwIHAYtWq/NVer0EkmxD73bSLa2jlySNq1F7Ckm2BxZX1XnN9mZJnl9Vt450XFU9nuRI4GJ64wWnVdWNSY4DllbVombfq5IsA54A3tfmySZJUjfa3D76EvDyge0nmrJdRjuwqhYDi1crO2bg+wLe03xJkoaszYDxxs0jpQA03z+1u5AkScPSJincneSAVRtJFgD3dBeSJGlY2tw+OoLenEcnNNvLgTd3F5IkaVjaTHPxU2D3ZlI8qurBzqOSJA1F65XXTAaStOEbrzeTJUkbAJOCJKmvzSypb0iyZfP9B5Ocl2Tn7kOTJE20Nj2Fv6+qXyXZE9gH+Dzw2W7DkiQNQ5uk8ETz55/Tm/foQnx5TZI2SG2SwookJ9GbHXVxkqe1PE6SNMW0+XA/kN7Eda+uqvuBZwLv6zQqSdJQtHlPYRZwYVU9kmQvetNnn9FpVJKkoWjTU/h34IkkfwicTG/hnDM7jUqSNBRtksJvq+px4HXAp5tFd2Z1G5YkaRjaJIXHkhwM/E/ggqZsk+5CkiQNS5uk8FZgD+AfqupnzUpsX+g2LEnSMLRZo3kZ8F7gh0l2ApZX1cc7j0ySNOHarNG8F/CvwK1AgO2S/K+quqrb0CRJE63NI6mfAl5VVTcBJPkj4CzgZV0GJkmaeG3GFDZZlRAAqurHONAsSRukNj2FpUlOBf6t2T4EWNpdSJKkYWmTFN4OvBP462b7m8BnOotIkjQ0IyaFJBsBp1XVIcD/nZiQJEnDMuKYQlU9ATwviVNlS9I00Ob20S3At5MsAh5aVVhV9hwkaQPTJin8tPl6CrBlt+FIkoZp1KRQVR9e18aT7AccD2wEnFpVH1tLvb8AvgzsUlU+2SRJQzLqewpJvpFk64HtZyS5uMVxGwEnAvsD84CDk8xbQ70tgXcD14wlcEnS+Gvz8trMZsU1AKrqPuDZLY7bFbi5qm6pqkeBs4EFa6j3EeDjwG9atClJ6lCbpPBEkjmrNpI8D6gWx80Gbh/YXt6U9SXZGdiuqi5s0Z4kqWNtBpo/AHwryZX0JsT7b8DC9T1xkqfQe/fhLS3qLlx1zjlz5oxSW5K0rtoMNF/U/Ea/e1N0VFXd06LtFfSW7lxl26ZslS2BnYArkgA8F1iU5IDVB5ur6mR6S4Eyf/78Nr0USdI6aNNToEkCF4xa8XctAeY2i/KsAA4C3jTQ5gPANqu2k1wBvNenjyRpeNqMKayTZl3nI4GLgR8B51bVjUmOS3JAV+eVJK27Vj2FdVVVi4HFq5Uds5a6e3UZiyRpdK16Ckn2TPLW5vuZzS0hSdIGps3La8cC7wf+rinahCfXVpAkbUDa9BReCxxAMxleVd2BcyBJ0gapTVJ4tKqK5oW1JFt0G5IkaVjaJIVzk5wEbJ3kMOBS4JRuw5IkDUObl9c+mWRfYCWwA3BMVX2j88gkSRNu1KSQ5D3AOSYCSdrwtbl9tCVwSZJvJjkyyXO6DkqSNByjJoWq+nBV/VfgncAs4Mokl3YemSRpwo1lmou7gF8A99JuPQVJ0hTT5uW1dzST1V0GPAs4rKpe3HVgkqSJ12buo+3oTZd9XdfBSJKGa61JIcmMqloJ/FOz/czB/VX1y45jkyRNsJF6CmcCrwGupfc2cwb2FfBfOoxLkjQEa00KVfWa5k9nRJWkaaLNQPNlbcokSVPfSGMKmwKbA9skeQZP3j6aAcyegNgkSRNspDGFw4GjgD+gN66wKimsBE7oOC5J0hCMNKZwPHB8kndV1acnMCZJ0pC0mSX100l2AuYBmw6Un9FlYJKkiddmltRjgb3oJYXFwP7AtwCTgiRtYNrMffR6YG/gF1X1VuAlwFadRiVJGoo2SeHXVfVb4PEkM+hNjLddt2FJkoahzdxHS5NsTW8JzmuBB4GrO41KkjQUbQaa39F8+7kkFwEzqur6bsOSJA3DSC+v7TzSvqr6fjchSZKGZaSewqdG2FfAn45zLNIG68xrbuP861YM7fzL7lzJvFkzhnZ+TR0jvbz2J+vbeJL9gOOBjYBTq+pjq+1/D/A24HHgbuCvqurn63teabI5/7oVQ/1gnjdrBgte6uw0Gl2b9xQ2B94DzKmqhUnmAjtU1QWjHLcRcCKwL7AcWJJkUVUtG6j2A2B+VT2c5O3AJ4A3ruPPIk1q82bN4JzD9xh2GNKI2jyS+i/Ao8DLm+0VwEdbHLcrcHNV3VJVjwJnAwsGK1TV5VX1cLP5XWDbVlFLkjrRJim8oKo+ATwG0HyIZ+RDgN5MqrcPbC9n5NlVDwW+3qJdSVJH2ryn8GiSzegNLpPkBcAj4xlEkr8E5gOvXMv+hcBCgDlz5oznqSVJA9r0FI4FLgK2S/JF4DLgb1sct4LfffN526bsdyTZB/gAcEBVrTHZVNXJVTW/qubPnDmzxaklSetixJ5CkgD/CbwO2J3ebaN3V9U9LdpeAsxNsj29ZHAQ8KbV2v9j4CRgv6q6a+zhS5LG04hJoaoqyeKqehFw4VgarqrHkxwJXEzvkdTTqurGJMcBS6tqEfBPwNOBL/XyD7dV1QHr8oNIktZfmzGF7yfZpaqWjLXxqlpMb7rtwbJjBr7fZ6xtSpK60yYp7AYckuTnwEP0biFVVb2408gkSROuTVJ4dedRSJImhTazpDrthCRNE20eSZUkTRMmBUlSn0lBktRnUpAk9ZkUJEl9JgVJUp9JQZLUZ1KQJPWZFCRJfSYFSVKfSUGS1GdSkCT1mRQkSX0mBUlSn0lBktRnUpAk9ZkUJEl9JgVJUp9JQZLUZ1KQJPWZFCRJfSYFSVKfSUGS1GdSkCT1mRQkSX0bd9l4kv2A44GNgFOr6mOr7X8acAbwMuBe4I1VdWuXMWnqOvOa2zj/uhXDDmOdLLtzJfNmzRh2GNKoOuspJNkIOBHYH5gHHJxk3mrVDgXuq6o/BP4f8PGu4tHUd/51K1h258phh7FO5s2awYKXzh52GNKouuwp7ArcXFW3ACQ5G1gALBuoswD4UPP9l4ETkqSqqsO4JpUPf+1Glt0xNT/oJtqq37bPOXyPYYcibbC6TAqzgdsHtpcDu62tTlU9nuQB4FnAPYOVkiwEFjabjyS5oZOI189WwAOTrM2xHt+2/mj1Rto/1n3bMPDv4Qbg3CNaRNg9r3e7/et1vSeRDeF6z21Vq6o6+QJeT28cYdX2m4ETVqtzA7DtwPZPgW1GaXdpVzGv58978mRrc6zHt60/Wr2R9o91n9fb6z0ZvqbT9e7y6aMVwHYD29s2ZWusk2Rjepnv3g5j6tLXJmGbYz2+bf3R6o20f133TTZe73b7vd7dtdnJ9U6TQcZd8yH/Y2Bveh/+S4A3VdWNA3XeCbyoqo5IchDwuqo6cJR2l1bV/E6C1qTj9Z5evN7D19mYQvXGCI4ELqb3SOppVXVjkuPodREXAZ8HvpDkZuCXwEEtmj65q5g1KXm9pxev95B11lOQJE09vtEsSeozKUiS+kwKkqS+KZ8UkrwwyeeSfDnJ24cdj7qXZIskS5O8ZtixqFtJ9kryzeb/+F7Djmc6mJRJIclpSe5a/c3lJPsluSnJzUmOBqiqH1XVEcCBwCuGEa/Wz1iud+P9wLkTG6XGyxivdwEPApvSmxVBHZuUSQE4HdhvsGCkCfaSHABcCCye2DA1Tk6n5fVOsi+9+bPumuggNW5Op/3/729W1f70fhH48ATHOS1NyqRQVVfRe29hUH+Cvap6FFg1wR5Vtaj5h3PIxEaq8TDG670XsDvwJuCwJJPy37DWbizXu6p+2+y/D3jaBIY5bXW6nsI4W+MEe819xtfR+wdjT2HDscbrXVVHAiR5C3DPwIeGpra1/f9+HfBqYGvghGEENt1MpaSwRlV1BXDFkMPQBKuq04cdg7pXVecB5w07julkKnW920ywpw2H13t68XpPElMpKSwB5ibZPslT6c2TtGjIMak7Xu/pxes9SUzKpJDkLOBqYIcky5McWlWPA6sm2PsRcO7gjKuaurze04vXe3JzQjxJUt+k7ClIkobDpCBJ6jMpSJL6TAqSpD6TgiSpz6QgSeozKUiS+kwK0npqpn2WNggmBU0bSY5LctTA9j8keXfz/fuSLElyfZIPD9T5apJrk9yYZOFA+YNJPpXkP4A9knwsybLm+E+u4dy7Jrk6yQ+SfCfJDk355knObY79SpJrksxv9h2c5IdJbkjy8Q7/aqQ+32jWtJHk+cB5VbVzsw7DT+jN4/8y4PXA4UDozbnziaq6Kskzq+qXSTajNz/PK6vq3iQFvLGqzk3yLOA7wI5VVUm2rqr7Vzv3DODhqno8yT7A26vqL5K8F5hbVYcn2Qm4jt56EXcA321iuw+4BPjnqvpqp39Jmvam/NTZUltVdWuSe5P8MfAc4AfNB/yrgFcBP2iqPh2YC1wF/HWS1zbl2zXl9wJPAP/elD8A/Ab4fJILgAvWcPqtgH9NMpfeEpObNOV7Asc38d2Q5PqmfBfgiqq6GyDJF4H/DpgU1CmTgqabU4G3AM8FTmvKAvxjVZ00WLFZwGkfYI+qejjJFfTWCgb4TVU9AdD89r8rsDe9HseRwJ+udt6PAJdX1WubHssV4/lDSePFMQVNN1+htz7wLvRm5KT586+SPB0gyewkz6b32/19TULYkd5tnd/THLdVVS0G/jfwkjVU24on1wd4y0D5t4EDm3bmAS9qyr8HvDLJNs1A9sHAlWP/caWxsaegaaWqHk1yOXD/wG/6lyR5IXB1EoAHgb8ELgKOSPIj4CZ69/jXZEvg/CSb0ut1vGcNdT5B7/bRB4ELB8o/05QvA/4TuBF4oKruTHI0cHnT5oVVdf76/OxSGw40a1ppBpi/D7yhqn4yCeLZCNikqn6T5AXApcAOzeL10oSzp6Bpo7k9cwHwlcmQEBqbA5cn2YRej+AdJgQNkz0FSVKfA82SpD6TgiSpz6QgSeozKUiS+kwKkqQ+k4Ikqe//A1xSxYppGDHsAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ " mu = 1.25e-8\n", " gen = 30\n", " crossPopDat = pd.read_csv(\"/home/training/share/MSMC-tutorial-files/results/EUR_AFR.combined.msmc2.final.txt\", delim_whitespace=True)\n", " plt.step(crossPopDat[\"left_time_boundary\"]/mu*gen, 2 * crossPopDat[\"lambda_01\"] / (crossPopDat[\"lambda_00\"] + crossPopDat[\"lambda_11\"]))\n", " plt.xlim(1000,500000);\n", " plt.ylim(0, 1.2)\n", " plt.xlabel(\"years ago\");\n", " plt.ylabel(\"relative cross coalescence rate\");\n", " plt.gca().set_xscale('log')" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1, 2, 3]" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[1,2,3]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "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.5.5" } }, "nbformat": 4, "nbformat_minor": 2 }