Skip to content

Instantly share code, notes, and snippets.

@benburrill

benburrill/life_expectancy_3.ipynb

Created July 20, 2016 19:54
Show Gist options
  • Save benburrill/c3c03c47006d1ba2107d297560b35be5 to your computer and use it in GitHub Desktop.
Save benburrill/c3c03c47006d1ba2107d297560b35be5 to your computer and use it in GitHub Desktop.
Download ZIP
life_expectancy_3.ipynb
Display the source blob
Display the rendered blob
Raw
life_expectancy_3.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# this better be the last life_expectancy(_\\d+)?!\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plot\n",
"import math\n",
"import itertools\n",
"import functools\n",
"import operator\n",
"import numpy\n",
"import mpl_toolkits.mplot3d as mp3d"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# \"acceleration of death\" is not a great name to describe the value, but it sounds pretty cool\n",
"# (it's really the beta from the Gompertz–Makeham law of mortality)\n",
"\n",
"def cdf_of_doom(age, inf_mort, acc_death):\n",
" return 1 - numpy.exp((1 - numpy.exp(acc_death * age)) * inf_mort / acc_death)\n",
"\n",
"def pdf_of_doom(age, inf_mort, acc_death):\n",
" return inf_mort * numpy.exp(acc_death * age) * numpy.exp((1 - numpy.exp(acc_death * age)) * inf_mort / acc_death)\n",
"\n",
"def life_table(specimens, inf_mort, acc_death):\n",
" for year in itertools.takewhile(lambda _: specimens, itertools.count()):\n",
" cdf_t = cdf_of_doom(year, inf_mort, acc_death)\n",
" cdf_tp1 = cdf_of_doom(year + 1, inf_mort, acc_death)\n",
" \n",
" died_this_year = round(specimens * (cdf_tp1 - cdf_t) / (1 - cdf_t))\n",
" specimens -= died_this_year\n",
" \n",
" yield died_this_year\n",
"\n",
"# old:\n",
"# died_this_year = (math.ceil(orig * cdf_of_doom(year + 1, inf_mort, acc_death)) - \n",
"# math.ceil(orig * cdf_of_doom(year , inf_mort, acc_death)))\n",
"\n",
"def life_expectancy(specimens, inf_mort, acc_death):\n",
" return sum(\n",
" itertools.starmap(operator.mul, enumerate(life_table(specimens, inf_mort, acc_death)))\n",
" ) / specimens\n",
"\n",
"le1000 = functools.partial(life_expectancy, 1000)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"[158.0, 314.0, 380.0, 143.0, 5.0]"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"list(itertools.islice(life_table(1000, 0.1, 1), 6))\n",
"# [158, 315, 379, 144, 4] (or similar)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"69.882999999999996"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"life_expectancy(1000, 0.0061, 0.0146228655953028) # from c=-(e^(a/b) Ei(-a/b))/b where c was set to 70 and a to 0.0061\n",
"# 69.396 (or as close to 70 as possible)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1000.0"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sum(list(life_table(1000, 0.0061, 0.0146228655953028))) # verifies we killed all the initial specimens\n",
"# 1000"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"life expectancy: 9.5004\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VNXWwOHfCr2G0KW30FEQAQGBIBB6VRCVK02vXEDl\nolIENIBeiqigAoogoh8IiErvYiJIlyqBEEC6BJBOKCn7++MMGDCZEjIlyXqfZx4nZ9bMrESysmef\nfdYWYwxKKaXSLj9vJ6CUUsq9tNArpVQap4VeKaXSOC30SimVxmmhV0qpNE4LvVJKpXFOFXoReU1E\n9tpur9qOBYjIahGJEJFVIuKfIH6oiESKyH4RCXZX8koppRxzWOhFpArQG3gMqA60EZGywBBgrTGm\nArAOGGqLrwx0ASoBLYEpIiLuSV8ppZQjzozoKwFbjDG3jDFxwC9AJ6AdMMsWMwvoYLvfDphrjIk1\nxhwFIoHaKZq1UkoppzlT6H8HGtimarIDrYDiQCFjTBSAMeYMUNAWXxQ4keD5p2zHlFJKeUFGRwHG\nmAMiMg5YA1wDdgJxiYWmcG5KKaVSgMNCD2CMmQnMBBCR97BG7FEiUsgYEyUihYGztvBTWCP+O4rZ\njt1DRPQPg1JKJYMxxqXzns6uuilg+28JoCMwB1gM9LCFdAcW2e4vBrqKSGYRKQ2UA7YmkazP3955\n5x2v56B5pq88IXXkmZZyTE15JodTI3rgexHJC8QAfY0xV2zTOfNFpBdwDGulDcaYcBGZD4QniNfR\nu1JKeYmzUzcNEzl2AWiaRPwYYMyDpaaUUiol6JWxDgQFBXk7BadonilL80w5qSFHSD15Jod4a1ZF\nRHRGR6lEiID+aqikiAjGxZOxzs7RK6VUiihVqhTHjh3zdho+r2TJkhw9ejRFXktH9Er5mLQ+oreN\nSL2dhs9L6ueUnBG9ztErpVQap4VeKaXSOC30SimVxmmhV0opJ4SFhVG8eHHHgT5IC71SSjnJ2a01\nevbsydtvv33PscmTJ1OrVi2yZs1Kr1693JFeknR5pVJKeUDRokUZMWIEq1at4saNGx59bx3RK6VU\nAqVLl2bs2LFUqVKFfPny0bt3b27fvv2PuAMHDtC4cWMCAgKoVq0aS5YsAeCLL75g9uzZjB8/nty5\nc9O+fXsAOnToQLt27cibN69Hvx/QQq+UUv8wZ84c1qxZw+HDh4mIiODdd9+95/HY2Fjatm1LixYt\nOHfuHB9//DHPP/88kZGRvPTSSzz//PMMGjSIK1eusGjRoiTexXO00CulfI5IytyS65VXXqFIkSLk\nyZOHYcOGMWfOnHse37RpE9evX2fw4MFkzJiRxo0b06ZNG7799tsH/M7dQ+folVL3iImLYdvpbfxy\n7Bdi42PJnik7OTLloFKBSjxR4gn8xP3jQ29fOFusWLG790uWLMmff/55z+N//vnnP1bglCxZklOn\n/rHHkk/QQq+UAuCnIz8xacskwo6FUSagDEElg8iROQcXblzg+u3rTN42mau3r/Kvh/9F90e6UzZv\nWW+n7DYnTvy97fWxY8coUqTIPY8XKVLknhiA48ePU6FCBcD51TmeooVeqXRu95ndDF47mEMXDjG8\n4XBmtJtBgRwF/hFnjGF31G5m7ZrF4zMeZ0CdAQx5YggZ/DJ4IWv3mjx5Mq1btyZbtmz873//o2vX\nrgB3e8/UqVOH7NmzM378eAYOHMiGDRtYunQpISEhABQqVIgjR47c85pxcXHExMQQFxdHbGwst27d\nImPGjGTI4IGfn5NbV/0X+B3YA8wGMgMBwGogAlgF+CeIHwpEAvuB4CRe0yil/slTvxq3Ym+Z/sv6\nm0LvFzKfbPnE3Iq95fRzT1w+YRp/1dg0nNnQHLt0zKX39fXf/VKlSpmxY8eaypUrm4CAANOzZ09z\n48YNExoaaooXL343Ljw83DRq1Mj4+/ubKlWqmEWLFt19LDIy0lSvXt0EBASYjh07GmOMCQkJMSJi\n/Pz87t5GjhyZZB5J/Zxsx13aftBh90oRKQJsACoaY26LyDxgOVAZ+MsYM15EBgMBxpghIlLZ9seg\nFtbG4GuBQHPfG2n3SqUS54nulX9e/ZOnv3uaAtkLMKPtV4TvyMOiRbBzJyxYAAEBjl8jLj6OCRsn\n8MGmD5j39Dwal27s1Hv7evfK0qVLM2PGDJ588kmv5uGN7pUZgBwikhHIBpwC2gOzbI/PAjrY7rcD\n5hpjYo0xR7FG9rVdSUop5T6bTmzisWm1KHG7OXlW/UDFUnl49VXIkQMKFIB33nHudTL4ZWDwE4OZ\n9/Q8nlnwDNtObXNv4irZHM7RG2NOi8gHwHEgGlhtjFkrIoWMMVG2mDMiUtD2lKLApgQvccp2TCnl\nRVFRMOa7VUw9040MS2dyLn8b2reHkSFQsqQV89dfUKkSvPgiPPywc6/buHRjZrSbQdtv27Ku+zoq\nF6jstu/BE3ztRGpKcFjoRSQP1ui9JHAZ+E5Engfu/0zh8mexOycuwNqvMS3v2aiUNxw4AIsWWbc9\nV9cR26EbQ0ov5L+h9cmT55/x+fLBqFHQvz+EhTm/Fr1thbZMuDWBFv/Xgl96/kKpPKVS9PvwpPtP\nonpbaGgooaGhD/QazszRPw00N8a8ZPv6X8DjwJNAkDEmSkQKAz8bYyqJyBCskwXjbPErgXeMMVvu\ne12do1cqEQ8yRx8XB5s2WYV98WK4fh3at4cyjdcz9kgnFnRZQKNSjRy+Ru3aMHAgPP+8a+//yZZP\nmLZjGtte2kbWjFkTjfH1OXpf4ek5+uPA4yKSVazPNE2AcGAx0MMW0x24c53vYqCriGQWkdJAOWCr\nK0kppZwXHW0V9l694KGHoF8/yJYN5syBEyegx1vbGPfHU3z79LcOizxAhgzw6acwaBBcueJaLv1r\n96di/ooMXzc8md+Ncgen9owVkXeArkAMsBN4EcgFzAeKA8eALsaYS7b4oUBvW/xrxpjVibymjuiV\nSoSzI/r582H2bPj5Z6hVC9q1s26lS/8dc/rqaWp/UZvJrSbTvmJ7l/Lo2RPy54f333ct/7+i/+Lh\nzx7m/zr+X6IrcXRE75yUHNHr5uBK+RhnCv2qVfDSSzBmDLRqlfhyyJuxN2n0VSPalW/HsIbDXM4j\nKgqqVoVffrFO0LpiReQK+izrw54+e/DP6n/PY1ronaOFXqk0zFGhv3oVqlWDadMgODjxGGMM3Rd2\n51bcLeY+NTfZK0kmTYKlS2H1atebhPVZ2ocbsTeY1WHWPce10DvHG+volVI+YsgQaNIk6SIP8NHm\nj9h7di9ftvvygZYL9usHZ87ADz+4/twJwRNYf2w96/5Yl+z39yW6laBSyiN++cU68frBB0nHbD21\nlXG/jmPhMwvJkTnHA71fxozwySfWCpzoaNeemzNzTsY2Hcsbq98g3sQ/UB6+IrlbCd6+fZsXX3yR\nUqVK4e/vz6OPPsrKlSvdleY/aKFXKpWIjobevWHKFBJdAw9w9dZVnvv+OSa3mkzJPCVT5H2DgqBe\nPet8gKs6V+5M5gyZmb1ndorkklrFxsZSokQJ1q9fz+XLlxk9ejRdunTh+PHjnknA1eY4KXXDxxsb\nKeUtSf1qvPGGMV272n9u9x+7m96Leqd4TidPGpMvnzGRka4/d8OxDab4h8VN9O1oY0zqaGo2ZswY\nU7lyZZM3b17Tq1cvc+vWrX80Ndu/f78JCgoyefLkMVWrVjWLFy82xhgzbdo0kylTJpMlSxaTK1cu\n065du0Tf5+GHHzY//PBDknkk9XMiGU3NtNAr5WMS+9XYssWYQoWMOXs26ed9u/dbU/6T8ubqratu\nyWvcOGNat07eczvN62TGrB9jjEkdhb5atWrm1KlT5uLFi6Z+/fpmxIgR9xT6mJgYU65cOTN27FgT\nExNj1q1bZ3LlymUOHjxojDGmR48eZsSIEUm+x5kzZ0y2bNlMREREkjEpWei1H71SPu7WLetiqIkT\nraZjiTlx+QSvrniVFc+vIGfmnG7JY8AAmDHDWoXTpo1rzx3bZCx1Z9Sld43eTsXLyJTpN2PeSd7q\nnjtbCQIMGzaMV155hSZNmtx9POFWgsA9WwkmnJtPTGxsLN26daNHjx6UL18+Wfm5Sgu9Uj7uvfeg\nbFl45pnEHzfG8PLSl3ml9ivULFLTbXlkzmydmP3Pf6BpU8iaeIeDRAXmC+TZqs8y7tdxTsUnt0Cn\nFHdtJWiMoVu3bmTJkoVPPvkk5RJ2QE/GKuXDdu+Gzz6DqVOTXsc+e+9sTl09xZAnhrg9n+BgeOQR\nmDDB9ecOqj+ImbtmpnxSbpDcrQSLFrUa9Sa1Oqd3796cP3+eH374wTM7S9looVfKR8XGWlM2Y8fC\nfXXmrqhrUby++nW+bPclmTJk8kheH34IH30Ex4659rzi/sVpU97FOR8vmTx5MqdOneLChQsOtxKM\njY0lNDSUpUuX8uyzzwKJbyXYp08fDhw4wOLFi8mcObNnvyFXJ/VT6oaPn5BRylvu/GqMGWNMs2bG\nxMcnHdvluy5m0OpBnkksgZEjjXnqKdeft+fMnlRxMjaltxI8duyYERGTLVs2kzNnTpMzZ06TK1cu\nM2fOnCTzSOrnhDu2EnQXbYGgVOJEYP9+eOIJ2L4dSpVKPG7RgUW8ueZNdvfZTbZM2Tya482bUKWK\nNa3UrJlrz/X1FgjpeStBpZQH9e4NISFJF/nrt6/z6spXmdZ2mseLPFgnYidOhFdegdu3Pf72ykVa\n6JXyQX5+0Ldv0o+P/mU0DUo0IKhUkMdyul+bNtZqoEmTvJaCW6TFrQR16kYpH3LkiFU8IyIgqSXW\n+8/tp+FXDdn7n70UzlnYswneJzIS6ta1VgcVdXJnaF+fuvEVOnWjVBpkjNVjHpIu8sYY+i3vx/AG\nw71e5AECA+Hll63dqJTvcljoRaS8iOwUkR22/14WkVdFJEBEVotIhIisEhH/BM8ZKiKRIrJfROw0\nU1VK3TF9uuOt++b+PpcLNy7Qr3Y/zyTlhLfegvXrrc6ayje5NHUjIn7ASaAO0B/4yxgzXkQGAwHG\nmCEiUhmYDdQCigFrgcD752l06kapv508CTVqwLp18PDDiW88cu32NSp8WoH5T8+nfon6nk/Sju++\ng9GjYccOq7WxPTp14xxvTt00BQ4bY04A7YE7W8fMAjrY7rcD5hpjYo0xR4FIoLaL76NUumEM9OkD\n/ftbO0clZeyGsTQu1djnijzA009bfXimTvV2Jioxrva6eQaYY7tfyBgTBWCMOSMiBW3HiwKbEjzn\nlO2YUioRc+bA8eP2d3E6eukoU7dPZXef3Z5LzAUiVh+cRo2snjwFCzp+jvIcpwu9iGTCGq0Pth26\n/zOFy5/FQkJC7t4PCgoiKCjI1ZdQKlWLirJ2b1q2zGoalpRBawbxWp3XKJa7WNJBXla5MnTvDkOH\nWl0u05qwsDC6dev2jx437hYaGkpoaOiDvYizl9BiFfmVCb7ejzWqBygM7LfdHwIMThC3EqiTyOsl\neemvUulF587GDLqvg8H9vxq/HP3FFP+wuLl++7rnEkumy5eNKVLEmE2bko5Jrb/797dAsCexfvTd\nunUzhQsXNv7+/qZChQpm+vTpdl8jqZ8TyWiB4Moc/bPAtwm+Xgz0sN3vDixKcLyriGQWkdJAOWCr\nC++jVLrwww/W+vMEH2z/Id7EM2DVAMY1HUf2TNk9llty5c4N48ZZ5xvi4rydjW8ZOnQof/zxB5cu\nXWLx4sUMHz6cnTt3euS9nSr0IpId60RswlnEcUAzEYkAmgBjAYwx4cB8IBxYDvS1/RVSStlcuGAV\nwxkzIJudDgZf7/6azBky07VqV88l94Cef976nlLr9E3p0qUZO3YsVapUIV++fPTu3ZvbifR5OHDg\nAI0bNyYgIIBq1aqxZMkSAL744gtmz57N+PHjyZ07N+3btwegcuXKZLU18TfGICIcPnzYM9+Uqx8B\nUupGKv34plRK6N7dmFdeSfyxO78a129fN8U+LGZ+Pf6rx/JKKTt3GlOwoDHnz//zMV//3XfnVoJ9\n+/Y12bNnNyJiatasaa5fT3o6LqmfE26eulFKpYAVKyAsDP73P/txEzdPpE7ROtQrXs8ziaWg6tWh\nc2cYMSKZLyCSMrdkurOVYJ48eRg2bBhz5sy55/GEWwlmzJjxnq0E7Zk8eTLXrl1jw4YNdOrUiSxZ\nsiQ7R1dooVfKg65csVoGTJsGOe1s7Xr2+lk+3PQhY5qM8VxyKWz0aOs8xI4dyXiy9cHmwW/J5K6t\nBMG64KlevXqcOHGCqR668EALvVIeNHiwtR2fox7uo8JG8Vy15wjMF+iZxNwgIMDa7/a117ydievc\ntZVgQrGxsR6bo9dCr5SHhIbCkiVO7Lea7yBzf5/L243e9kRabtWzp9Xh8r5d9XxeSm8leO7cOebN\nm8f169eJj49n1apVzJ07l6ZNm3rk+9FCr5QHREfDiy/ClCmQJ4+D4CZDeaPeG+TPnt8jubmTnx+0\nbm39gUtNnnvuOYKDgylXrhyBgYEMGzYM+HuknilTJpYsWcLy5cvJnz8//fv355tvviEw0PoE1rt3\nb/bt20fevHnp1KkTfn5+TJ06leLFi5M3b14GDRrEpEmTaN26tUe+H+1Hr5QHvP46/Pmn1e7Ans0n\nN1P3w85EjzvolZ2j3GHhQvj0U1i71vra15uapcWtBF3tdaOUctHBgzBrlrUPrD3GGIasHQJh76SZ\nIg/W+YgXXoDLl8Hf33G8Snk6daOUmw0fbvWzKVDAftzKQyuJuh4Fu3p4JC9PyZHD2uh81SpvZ+Ic\n3UowJd9Yp25UOrB9O7RrZ52QzJEj6bh4E0+Nz2sQ0iiETpU7PsjKQJ80dSps3AjffOP7Uze+QrcS\nVCqVGDrUumjIXpEHmLN3DtkzZadDxQ72A1Optm2tC8ViY72dSfqkc/RKucnatXD0qLXaxp7bcbcZ\n8fMIZnWYlSanDQCKFYMSJWDTJsexKuXpiF4pN4iPhyFD4N13IVMm+7HTfptGpfyVaFiyoWeS85K2\nbWHxYm9nkT7piF4pN/j+e+sK/M6d7cddv32d99a/x/LnlnsmMS9q2xa6dbNaBaTVTy4pqWTJkin2\nWlrolUphMTEwbJi1dtzPwWfmSVsm0ahkI2o8VMMzyXnRo49avX7Cwo4S6ERnhwErB5A3W940cYWw\nt+nUjVIp7MsvoXhxx/1sLty4wEebP2J049GeSczL/PygTRvnr5J9qtJTfL//e/cmlU5ooVcqBUVH\nw6hRMHas4y65438dT8eKHVN14zJXtW3rfKGvV7weUdeiOHThkHuTSgec3WHKX0S+E5H9IrJPROqI\nSICIrBaRCBFZJSL+CeKHikikLT7Yfekr5Vs+/hjq1oVatezH/Xn1T6b9Ni3dTUs0aQK//QYXLzqO\nzeCXgY4VO/LD/h8cByu7nB3RTwKWG2MqAY8AB7A2AV9rjKkArAOGAohIZaALUAloCUwRPfOi0oEL\nF+CDD6zWvI68+8u79Kzek2K5izkOTkOyZ4egIGtNvTM6Veqk0zcpwGGhF5HcQANjzEwAY0ysMeYy\n0B6YZQubBdy50qMdMNcWdxSIBGqndOJK+ZqxY6FjR6hQwX7ckYtHmLtvLkMbDPVMYj7GlemboFJB\nRP4Vyakrjjf0UElzZkRfGjgvIjNFZIeITLNtFl7IGBMFYIw5AxS0xRcFEnbkP2U7plSadfIkTJ8O\n77zjOHZk2Eheqf1KmmhDnBxt2lh9b2JiHMdmypCJFuVasCxymfsTS8OcWV6ZEXgU6GeM2S4iH2FN\n29zfhMHl5hUhISF37wcFBREUFOTqSyjlE0aOhJdegqIOhjT7zu5jReQKDr2afk8wPvQQlC0LGzZA\n48aO49uWb8uc3+fw75r/dn9yPig0NJTQ0NAHeg2HTc1EpBCwyRhTxvb1E1iFviwQZIyJEpHCwM/G\nmEoiMgRrl/JxtviVwDvGmC33va42NVNpwoED0KCB1Y44IMB+7FPzn6Jusbq8Ue+NJGNEHmi701Rh\n1Ci4dAk+/NBx7MUbFyk5sSRn3jhD9kzZ3Z+cj3NLUzPb9MwJESlvO9QE2AcsBnrYjnUHFtnuLwa6\nikhmESkNlAO2upKUUqnJ8OHwxhuOi/z209vZcnIL/Wr180xiPuzOPL0zf9ACsgVQs0hN1h5Z6/7E\n0ihnr4x9FZgtIpmAI0BPIAMwX0R6AcewVtpgjAkXkflAOBAD9NWhu0qrtm61GnV9/bXj2OHrhjOs\nwbA0talIclWvDjdvQkQEVKzoOL5t+bYsiVhCuwrt3J9cGqT96JVKJmOsdeFdu8K/HUwfhx0No9fi\nXuzvt5/MGTLbjU0PUzcAfftCqVIwaJDj2Mi/Imn4VUNODTyFn6Tv6zy1H71SHrRmDZw6Bb162Y8z\nxjBs3TBGBo10WOTTE1eWWQbmC8Q/iz+/nf7NvUmlUVrolUqGO22I33sPMjqYAF1xaAUXb17k2arP\neia5VKJxY9izB/76y7n4tuXbsuSgk38Z1D200CuVDPPnQ4YM8NRT9uPiTTzD1g3j3cbvksEvg2eS\nSyWyZoUnn4TlTnZobltBC31yaaFXykUxMdZKG2caly0IX0BGv4xpdovAB+Vqk7Pjl49z4vIJx8Hq\nHlrolXLR9OlQpox1Itae2PhYRvw8gv89+T/daCMJrVvD6tVw+7bj2Ix+GWlZrqVeJZsMWuiVcsH1\n6zB6NIwZ4zh21q5ZPJTzIZqWaer+xFKpQoWgUiUIC3MuvlVgK1YccrIjmrpLC71SLpg40boKtmZN\n+3E3Y28yMmwkY5qM0dG8A65M3wSXDSb0aCi3Ym+5N6k0Rgu9Uk766y/46CNrw29Hpm6bSvXC1alb\nvK77E0vlXLlKNn/2/FTKX4n1x9e7P7E0RAu9Uk4aM8ba7NvRfqdXb11l7K9jee9JJxrTK6pWtYr8\nvn3Oxbcs15IVkTp94wot9Eo54fhxay/YESMcx3646UOalWlGtULV3J9YGiDi2vRNy8CWOk/vIi30\nSjlh5Ejo0weKFLEfdz76PB9v/ZiRQSM9k1ga4Uqhf6zIY5yLPsexS8fcm1QaooVeKQfCw60i5ExP\nljHrx/BMlWcom7es+xNLQxo1sn7OZ886jvUTP5qXba6jehdooVfKgWHDrCKfJ4/9uOOXj/PV7q8Y\n0dCJ+R11jyxZoFkzWObkEvmW5XT6xhVa6JWyY/Nm2L4d+jnRQj4kNIQ+NfvwUK6H3J9YGuTK9E3z\ncs11maULtNArlQRjrMZlISGQzUEL+fBz4Sw9uJRB9Z2Y31GJatUKfvrJ6lPvyJ1llhuOb3B/YmmA\nU4VeRI6KyG4R2SkiW23HAkRktYhEiMgqEfFPED9URCJFZL+IBLsreaXcaeVKiIqC7t0dxw5bN4xB\n9Qfhn9XfcbBKVP78UK0aOLs9astyLVke6WRHtHTO2RF9PNb+sDWMMbVtx4YAa40xFYB1wFAAEamM\ntdtUJaAlMEX00kCVysTHw9ChzrUh3nRiE7+d/o3+tft7Jrk0zNVllisPr3RvQmmEs4VeEoltD8yy\n3Z8F3GnP1w6Ya4yJNcYcBSKB2iiVisyda7XR7djRfpwxhqE/DSUkKISsGbN6Jrk0zJWrZGs+VJMz\n185oN0snOFvoDbBGRLaJyIu2Y4VsG4djjDkDFLQdLwok/Mmfsh1TKlW4fdu6MMqZNsTLI5dz9vpZ\nXnjkBc8kl8ZVqgSZMsHu3Y5jM/hloGmZpqw5ssb9iaVyzhb6+saYR4FWQD8RaYBV/BNKB7tcqvRg\n2jQoXx6CguzHxcXHMeSnIYxtOpaMfg7md5RTRKBdOxdW35RtzqrDq9ybVBrg1L9OY8yftv+eE5GF\nWFMxUSJSyBgTJSKFgTuXOpwCiid4ejHbsX8ICQm5ez8oKIggR79ZSrnZxYtW07KVTkz9fr37a/Jk\nzUPb8m3dn1g60rattdrJmXYTwWWDeXPNm8TFx6XZHbxCQ0MJdfYMdRLEOJgME5HsgJ8x5pqI5ABW\nAyOBJsAFY8w4ERkMBBhjhthOxs4G6mBN2awBAs19byQi9x9SyuteeQViY2HqVPtxN2JuUP7T8sx/\nen6Kd6gUcW6OOq2KiYGCBa0rZR9y4pKEqlOqMqPdDOoUq+P+5HyAiGCMcWmBizMj+kLAjyJibPGz\njTGrRWQ7MF9EegHHsFbaYIwJF5H5QDgQA/TViq5Sg927rb1gw8Mdx3685WNqF62tbYjdIFMmaN7c\nukr2xRcdx9+ZvkkvhT45HI7o3fbGOqJXPsQYaNgQunWDl1+2H/tX9F9UnFyRDT03UCF/hRTPJb2P\n6AFmz7b+6C5a5Dh29eHVjAobxYZe6ePiqeSM6PXKWKWAOXMgOtq5EeS7v7zL05WedkuRV5aWLeHn\nn+HGDcexDUo0YHfUbi7fvOz+xFIpLfQq3btyxWpa9umnkMHB+bxDFw7xzZ5vGNlY2xC7U968UKOG\n1RLBkWyZslGveD1++sOJ4HRKC71K90aNsuaE6zox3T70p6EMrDuQgjkKOg5WD8TlZZaHdJllUnTx\nr0rXwsNh1iz4/XfHsRtPbGTLyS3M6jDLcbB6YG3bQuPG1vkKRxeuNS/bnI+3fIwxRjdjT4SO6FW6\nZQy8+qq1XrtQIUexhtdXv867T75L9kzZPZNgOle+POTMCTt2OI6tXKAysfGxHPzroPsTS4W00Kt0\n6/vvre6Uffs6jv0u/Dtuxd6i28Pd3J+YusvZJmciQnDZYG2HkAQt9Cpdun4dBg60TsA66k55M/Ym\nQ9YOYULwBPxEf2U8yZVulsFlg1l9eLV7E0ql9F+tSpfGjIEnnrD2KnVk4uaJPFL4EZ4s/aT7E1P3\nqF8fjh6FkycdxzYt05SwY2Hcjrvt9rxSGy30Kt05dAg++wzef99x7JlrZ5iwcQLvN3MiWKW4jBmh\nRQtYutRxbP7s+QnMG8jmk5vdn1gqo4VepTsDBljr5os60Tx7+Lrh9Kzek3J5y7k/MZUoV5ZZBpcN\nZs1hnae/nxZ6la4sWWKN6AcMcBy788+dLItcxvCGw92fmEpSixawfr11XsWR4LLBrD6i8/T300Kv\n0o2bN+E2w3KoAAAgAElEQVS11+DjjyFzZvuxxhgGrBrAyKCRug+sl/n7Q61asHat49i6xeqy/9x+\nLty44P7EUhEt9CrdeP9967L6YCe2q18QvoBLNy/Ru0Zv9yemHHJ29U2WjFloULIBPx3RdggJaaFX\n6cLRozBxInz4oePY67ev88aaN/i4xcdpdjOL1KZtW+uEbHy849jgMrrM8n5a6FW68Prr1rx8yZKO\nY8duGEu94vVoVMqJtZfKI8qWtRqdbdvmOPbOPL22Qf+bFnqV5q1eDbt2wZtvOo49cvEIU7ZP0eWU\nPsjZ6ZuK+SsSb+K1HUICThd6EfETkR0istj2dYCIrBaRCBFZJSL+CWKHikikiOwXESdmRJVyj9u3\nrX42EydC1qyO4weuGsjrdV+nWO5i7k9OucTZZZYiotM393FlRP8a1vaAdwwB1hpjKgDrgKEAtj1j\nuwCVgJbAFNF2cspLJk2yPva3aeM4dtWhVfx+9ncG1h3o/sSUyx5/HE6fhmPHHMfqMst7OVXoRaQY\n0AqYnuBwe+BOv9ZZQAfb/XbAXGNMrDHmKBAJ1E6RbJVywalTMG6cVewdDTVuxt6k/4r+TGwxkawZ\nnRj6K4/LkAFatXLuKtkmZZoQdlTbIdzh7Ij+I+BNIOHZjULGmCgAY8wZ4M5ODEWBEwniTtmOKeVR\nb74JffpAOScuan3/1/epWrAqbco7MfRXXuPsPH3+7PmpkL8Cm05scn9SqYDDQi8irYEoY8wuwN64\nSE9xK58RFga//gpDhzqOPXzhMJO2TGJi84nuT0w9kOBg2LgRrl51Ilbn6e9yZoep+kA7EWkFZANy\nicg3wBkRKWSMiRKRwsBZW/wpoHiC5xezHfuHkJCQu/eDgoIICgpy+RtQ6n6xsdC/P3zwAeTIYT/W\nGMOrK1/lzXpvUjKPE2svlVflzm3N1a9eDU89ZT82uGwwb6x5g/eavOeZ5NwkNDSU0NDQB3oNcWWt\nqYg0Al43xrQTkfHAX8aYcSIyGAgwxgyxnYydDdTBmrJZAwSa+95IRO4/pFSK+PhjWLwY1qxxPDf/\n4/4fGbZuGLv67CJzBgd9ETxExNr9SiXuk0/gt9/gq6/sx92Ou02B9wtw+NXD5M+e3yO5eYKIYIxx\naYHLg6yjHws0E5EIoInta4wx4cB8rBU6y4G+WtGVp0RFwejRVrF3VOSv3rrKaytfY3KryT5T5JVj\n7drBsmUQE2M/LnOGzDQq2Yi1R5xokpPGuTSiT9E31hG9coNevawrKCdMcBw7YOUALt+6zMz2M92f\nmAt0RO9Y7drw3nvQrJn9uE+3fsqOP3fwZfsvPZOYByRnRO/MHL1SqcLmzbBqFezf7zh226ltzP19\nLvv67nN/YirFde4M333nuNAHlw1m7IaxGGNIz5fzaAsElSbExVknYMeNs07Y2RMTF8NLS15iQvAE\n8mXP55kEVYp6+mn48UfrxLs9gXkDyeiXkf3nnfjrn4ZpoVdpwowZkC0bPP+849iJmydSMEdBnq/m\nRLDySaVLQ6lS4GgxiojopuFooVdpwF9/wYgR8Omnjk/AHrl4hHG/jmNq66np+qN8WnBn+saR5mWb\np/tCrydjVar3n/9Ym0h/8on9OGMMTb5uQotyLRhUf5BnkksGPRnrnD/+gDp1rP43Ge2cbbx44yIl\nJ5bk7Jtn00R7C08vr1TK63bssOZqR41yHPvFji+4dvuaNi1LI0qXhhIlrKug7QnIFkC1QtVYf2y9\nZxLzQVroVap16xa89JK1zC4gwH7s8cvHGbZuGDPbzySjny42Syucnb5pUbYFKw+tdH9CPkoLvUq1\n3njDOiHXq5f9OGMMLy99mdfqvEaVglU8kpvyjM6dnVt906JcC1Ye1kKvVKry3XewfLm12sbROdVZ\nu2dx5toZBtcf7JnklMeUKQPFisEvv9iPq1mkJmevn+X45eOeSczHaKFXqc6hQ9C3L8ybB3ny2I89\ncfkEb655k5ntZ5IpQybPJKg8ypnpGz/xI7hsMKsOrfJMUj5GC71KVW7etH6x33kHHnvMfmy8iafn\nop789/H/Ur1wdc8kqDyuc2f44Qfrojl7WpRNv9M3WuhVqjJwIAQGQr9+jmMnb53M9ZjrPr2UUj24\nsmWhaFHH0zfBZYP56chPxMQ56IaWBmmhV6nGvHlWH/IvvnA8Lx9xPoKRYSOZ1WGWrrJJB5yZvimU\nsxBl85Zl88nNnknKh2ihV6nCwYNWL5vvvgN/f/uxsfGxdF/YnZCgEMrnK++ZBJVXuTR9kw6XWWqh\nVz7vxg3rF3n0aKhRw3H8qLBR+Gf1p2+tvu5PTvmEcuXgoYdgvYNrotLrMkst9MrnvfYaVK4ML7/s\nOHb9sfV8seMLZnWYhZ/oP+/0xJnpm8eLPc7hC4eJuhblmaR8hDObg2cRkS0islNE9orIO7bjASKy\nWkQiRGSViPgneM5QEYkUkf0iEuzOb0ClbbNnWx0KP//c8bz8xRsX6fZjN75o+wWFcxb2SH7Kd3Tu\nDN9/b3/6JlOGTDQp0yTdTd84LPTGmFtAY2NMDaA60FJEagNDgLXGmArAOmAogG3P2C5AJaAlMEW0\nTaBKhgMHYMAAa5TmqMe8MYY+y/rQrnw72pRv45kElU8JDITChWHDBvtxrQNbsyxymWeS8hFOfbY1\nxkTb7mbB2pXKAO2BWbbjs4AOtvvtgLnGmFhjzFEgEqidUgmr9CE62hqh/e9/8MgjjuNn7ppJ+Llw\nxjcb7/7klM9yZvqmVWAr1hxZk66WWTpV6EXET0R2AmeANcaYbUAhY0wUgDHmDFDQFl4UOJHg6ads\nx5Ry2iuvWAX+xRcdx/5+9ncGrx3MvKfnkS1TNvcnp3yWM9M3hXMWplzecvx64lfPJeZlTi0wNsbE\nAzVEJDfwo4hUwRrV3xPm6puHhITcvR8UFERQUJCrL6HSoK+/ho0bYds2x/Py125fo/N3nZnQbAKV\nC1T2TILKZ5UvDwULwq+/QsOGSce1DmzNsoPLCCoV5LHckis0NJRQR1tpOeDyxiMiMgKIBl4Egowx\nUSJSGPjZGFNJRIYAxhgzzha/EnjHGLPlvtfRjUfUP4SHQ6NG8PPPULWq/VhjDC8sfIFMfpn4sv2X\nnknQA3TjkQfz3ntw5oz9jWi2ndrGCwtfYH+/1LeXrFs2HhGR/HdW1IhINqAZsB9YDPSwhXUHFtnu\nLwa6ikhmESkNlAO2upKUSp+uX7c+eo8b57jIA3y580t2/rmTT1t96v7kVKpxZ/omPj7pmJpFanLx\nxkWOXDziucS8yJk5+oeAn0VkF7AFWGWMWQ6MA5qJSATQBBgLYIwJB+YD4cByoK8O3ZUjxlgdKR97\nDHr2dBy/488dDPlpCPM7zyd7puzuT1ClGuXLQ4EC1vRNUvzEj5aBLVl2MH2svtE9Y5VP+PJL+OAD\n2LoVcuSwH3s++jyPTXuMCcETeLry055J0IN06ubBvfsunD0LH3+cdMyC8AVM3zGdld1S15r65Ezd\naKFXXrd3Lzz5pLX3Z2UH51Pj4uNoMbsFjxZ+lHHNxnkmQQ/TQv/gIiKsf1MnToBfEvMWl29epthH\nxTjz+hlyZHYwuvAhujm4SnWuXbPmVD/4wHGRBxi+bjjGGN5r8p77k1OpVoUKkC+ftXorKf5Z/ald\ntDY//fGT5xLzEi30ymuMgT59oH59eOEFx/Hz983n29+/5dunvtXWw8ohZy6eah3YmqUHl3omIS/S\nQq+8Zvp02L3b/jK4O7af3k6/5f1Y2HUhBXIUcH9yKtXr3BkWLLC/+qZt+bYsObiEeGMnKA3QQq+8\nYvdueOsta8SV3cGimdNXT9NxXkc+b/O5bgmonFaxIuTNC5s2JR0TmC+QfNnyseXklqSD0gAt9Mrj\nrlyxRlsTJ1q/jPbciLlBh7kd6FOzD50qdfJMgirNcGb6pkPFDiw8sNAzCXmJrrpRHmUMPPus1Y1y\n2jT7sfEmnue+fw4RYU6nOaSXJqi66ibl7N8PzZrB8eNJr77Zfno7z33/HBH9I1LFvzFddaN83mef\nWe2HJ01yHDt07VBOXjnJzPYzU8UvoPI9lSpBnjyw2c42sTUfqkl0TDQHzh/wXGIepoVeecyOHfD2\n2zB/PmRz0GRy6rapLIxYyKKui8iaMatnElRpkqPpGxFJ89M3WuiVR+zdC23bwpQp1iXq9iw9uJRR\nv4xi+XPLyZc9n2cSVGmWM6tvOlbsyMIILfRKJdvWrdC0qXVRVOfO9mM3ndhEz0U9WfjMQsrmLeuZ\nBFWaVrmydU7I3vRNw5INOXThEKeunPJcYh6khV65VVgYtGljrZnv2tV+7O9nf6fDvA583eFr6hSr\n45kEVbrgaPomU4ZMtApsxeKIxZ5LyoO00Cu3WbkSnn4avv3Wmrax54+Lf9BydksmNp9Iy8CWnklQ\npRvOTN90qNAhzU7faKFXbvH991Zbg0WLoEkT+7Fnrp0h+P+CGfrEUJ6t9qxnElTpSpUqkCsXbLFz\nXVTzcs3ZdGITF29c9FxiHqKFXqW4WbOgf39YtQrq1bMfe+76OZp83YQXHn6BvrX6eiZBlS45mr7J\nmTknTcs05ccDP3ouKQ/RQq9S1JQpMHw4rFsHNWrYj71w4wLNvmlGx4odGdFohGcSVOmWM9M3z1R5\nhnn75nkuKQ9xZivBYiKyTkT2icheEXnVdjxARFaLSISIrLqz3aDtsaEiEiki+0Uk2J3fgPId48bB\nhAnWCdhKlezHXrp5ieBvgmlWphmjG4/2TIIqXatSxdrUZqudjU3blG/D5pObOXf9nOcS8wBnRvSx\nwEBjTBWgLtBPRCoCQ4C1xpgKwDpgKICIVAa6AJWAlsAU0csa0zRjYNgw+OorWL8eypSxH3/xxkWa\n/19z6hevz/hm4/WqV+URIo6nb3JkzkGLci34Yf8PnkvMAxwWemPMGWPMLtv9a1gbgxcD2gOzbGGz\ngA62++2AucaYWGPMUSASqJ3CeSsfER8PAwbA8uXwyy9QtKj9+PPR52nydRPqFavHxBYTtcgrj7oz\nfWOvl1DXKl3T3PSNS3P0IlIKqA5sBgoZY6LA+mMAFLSFFQVOJHjaKdsxlcbExcGLL8K2bfDzz9aG\nzPZEXYviyVlP0rxscz5s/qEWeeVxVata7TfsTd+0DGzJzjM7OXPtjOcSczOnt+kRkZzAAuA1Y8w1\nEbn/b6LL/fZCQkLu3g8KCiIoKMjVl1Becvs2dOsGFy7A6tWQM6f9+JNXTtLsm2Z0rdKVtxu9rUVe\neUXC6Zs6SVyTlzVjVtqUb8OC8AX0r93fswkmIjQ0lNDQ0Ad6DafaFItIRmApsMIYM8l2bD8QZIyJ\nEpHCwM/GmEoiMgQwxphxtriVwDvGmC33vaa2KU6lbtywflkyZIB58yCrg55jEecjaP5/zelfuz9v\n1HvDM0mmYtqm2L327IF27eCPP6yfdWKWHlzK2A1j2dBrg2eTc4I72xR/CYTfKfI2i4EetvvdgUUJ\njncVkcwiUhooB9j5oKRSk6tXoXVrq3fIggWOi/y2U9sImhVESFCIFnnlE6pVgyxZrCnHpASXDWb/\n+f2cuHwi6aBUxJnllfWB54EnRWSniOwQkRbAOKCZiEQATYCxAMaYcGA+EA4sB/rq0D1tuHjR2sSh\nXDn45hvIlMl+/OrDq2k9pzWft/mcHtV7eCRHpRwRgS5d7K++yZwhMx0qdEgzJ2V1hynllKgoCA62\nulBOmJD0R947pu+YzvB1w1nQZQFPlHjCM0mmETp1437OTN+EHQ2j3/J+7P3PXp86p6Q7TCm3OHEC\nGjWCjh0dF/l4E89bP73FuF/H8UvPX7TIK590Z/pm+/akYxqWbMjN2JtsP20nKJXQQq/sOnwYGjaE\nl16CkBD7RT46Jprnvn+O0KOhbOy1kfL5HOwwopSXOHPxlIjQo3oPZu6a6bnE3EQLvUrSvn3WSH7o\nUHj9dfuxxy8fp8HMBmT0y8hPL/xEgRwOFtUr5WV3Cr29abLuj3Rn3r553Iy96bnE3EALvUrUli1W\ne+Fx4+Df/7Yfu+H4Bh6f/jjPVn2Wbzp+Q7ZMDjaEVcoHPPywtaDgt9+SjinuX5yaD9VM9fvJaqFX\n97h1y+o+2bYtTJsGzz+fdKwxhslbJ9NpXie+bP8lb9R7w6dOWilljzPTNwA9q/dM9dM3uupG3bV9\nO/ToAWXLwmefwUMPJR177fY1XlryEvvP7WdBlwWUy1vOY3mmdbrqxnN27YJOnaxzUUmNUW7E3KDY\nR8XY9fIuivsX92yCidBVNypZbt2Ct96yLoR66y1YuNB+kd93dh+1vqhFjkw52NR7kxZ5lWo98oh1\nhfeOHUnHZMuUjS6Vu/D17q89l1gK00Kfzm3dCo8+CgcOwO7d8NxzSY9sjDF8vv1zgmYFMajeIKa3\nm67z8SpVc3b6pleNXszYOYO4+DjPJJbCtNCnUzdvwpAh1kUjb79t7fFauHDS8X9F/8VT85/is98+\nY33P9fSs0dNzySrlRs6svqlVtBYFchRgeeRyzyWWgrTQp0ObN1vb/B06ZI3in3nG/vr4VYdWUf3z\n6pTOU5rNvTdTMX9FzyWrlJtVr279+9+5037cq7VfZdKWSfaDfJQW+nTkxg14803o0AFGjbKakhUq\nlHT81VtX6bO0D/9e+m9mtp/JB80/IEvGLJ5LWCkPcHb6pnOVzoSfC2ff2X2eSSwFaaFPJzZtskbx\nx4/D3r3WP2x71v2xjkc+e4TbcbfZ02cPTcs09UyiSnlBt24wcyacPJl0TOYMmenzWB8+2fqJ5xJL\nIbq8Mo2LjoYRI2DOHPjkE3j6afvxf0X/xZtr3mTtkbVMaT2FNuXbeCZRdZcur/SO//0PVq6Edesg\nYxJbMkVdi6Li5IocfvUwebPl9WyCNrq8Ut1jwwZr/vH0aWsUb6/IG2OYs3cOVadWJWfmnOzru0+L\nvEpXhgyx9lcYNSrpmEI5C9G2fFtm7JjhucRSgI7o06DoaBg2zNr96dNPrQtC7Nl3dh/9V/Tn4o2L\nfNbmMx4v9rhnElWJ0hG990RFWcuNv/7aagGSmN9O/8ZT85/i0KuHyOjn9G6sKUZH9Ir1662LQM6e\ntUbx9or85ZuXeX3V6wTNCuKpSk+x/d/btcirdK1QIZg1C154wSr6ialZpCbF/Yszf998zyb3AJzZ\nYWqGiESJyJ4ExwJEZLWIRIjIKhHxT/DYUBGJFJH9IhLsrsTVva5fh9deg65drZ7xs2dDvnyJx8bE\nxTB562QqfFqBy7cus6/vPvrX7u+V0YlSvqZpU+jZ0yr28fGJx7zd8G1GhY1KNRdQOTOinwk0v+/Y\nEGCtMaYCsA4YCiAilYEuQCWgJTBFtMuVW8XHw5IlVie+CxesUXz79onHGmNYdGAR1aZWY2HEQlZ1\nW8X0dtMpmKOgZ5NWyseFhFhToOPHJ/540zJNyZ89P3N/n+vRvJLNGOPwBpQE9iT4+gBQyHa/MHDA\ndn8IMDhB3AqgThKvaVTyRUcb8/nnxlSoYMyjjxqzZIn9+HVH1pnHpz9uqk2pZpZGLDXx8fGeSVS5\nTH81fMOJE8YUKmTMhg2JP7728FpT/pPyJiYuxqN52WqnU7X7zi25c/QFjTFRtmp9BrgzJCwKJNw2\n/ZTtmEoh587ByJFQqpQ1kv/sM6vrZJskFshsOL6BZt8046UlL/FK7VfY1WcXrcu31nbCSjlQrBh8\n8YXV/+nChX8+/mTpJymUo1CqGNWn1KRsstYIhISE3L0fFBREUFBQCqWT9kREwIcfwvz51sVOoaFQ\nqVLiscYYwo6FMSpsFEcvHeWtBm/R/ZHuZMqQyaM5K5XatW0LP/8MvXrBjz/e2ypERAgJCqHP0j50\nrdrVbee4QkNDCQ0NfaDXcGp5pYiUBJYYYx62fb0fCDLGRIlIYeBnY0wlERmC9bFinC1uJfCOMWZL\nIq9pnHnv9MwY+OUX6+Tqli3wn/9Av35QMIkp9XgTz+KIxYz/dTznos8xrMEwnq/2vBb4VEaXV/qW\n27ehfn3417/g1VfvfcwYQ9CsIHpV70X36t09kk9yllc6W+hLYRX6aravxwEXjDHjRGQwEGCMGWI7\nGTsbqIM1ZbMGCEysomuhT1pMjNWH5oMP4OpVGDjQWgGQLYmOwNEx0fzfnv/jw00fkitLLgbXH0zH\nih3J4JfBs4mrFKGF3vccPgx168KKFVCz5r2PbTyxkWcWPEN433ByZcnl9lzcUuhFZA4QBOQDooB3\ngIXAd0Bx4BjQxRhzyRY/FOgNxACvGWNWJ/G6Wujvc+UKTJ8OkyZZc/BvvGFtBuKXxJmUE5dPMHnb\nZGbsnEHdYnX57+P/JahUkM6/p3Ja6H3T/PnWxjw7dkDu3Pc+1mNhDwrmKMj4Zkks00lBbhvRu4MW\n+r+dOGEV95kzITgYXn8dHnss8di4+DhWH17NZ799xobjG3jh4RfoV7uf7vKUhmih9119+sDly1bv\nqITjqahrUVSdWpWwHmFULlDZrTlooU9lduywpmdWrLD2an3tNShZMvHYY5eOMWv3LL7c+SUFchSg\nT03rBFCOzDk8mrNyPy30vuvGDahTx5qrf/HFex/7dOun/LD/B3564Se3fqrWQp8KxMdbhf2DDyAy\n0iruL70E/v7/jL12+xoLDyzkq11fsevMLrpW7UqvGr149KFHPZ+48hgt9L5t/35o2NBajVO16t/H\nY+NjqfVFLQbXH0zXql3d9v5a6H3UtWuwZg0sW2bdihSxTrB26QKZ7lsQcyv2FqsPr2bO73NYEbmC\nJ0o8wQuPvEC7Cu3ImjGrd74B5VFa6H3fV1/B++/Dtm2QPfvfxzee2Ejn7zqzu89u8mfP75b31kLv\nQ44csYr60qWwcSM8/rh1UVPr1lDuvun0m7E3WXN4Dd+Ff8fSg0upUrAKz1V9js5VOrvtH4vyXVro\nfZ8x1kq4LFmsBRQJDVoziPBz4Sx5dolbpnC00HtRbKxV0O8U9/PnraLepg00awa57lt1deHGBZYd\nXMaiiEWsPbKWRwo/QufKnelUqRNFchXxzjehfIIW+tTh6lVr0cQ771hXz94RExdDg5kN6FKlCwPr\nDkzx99VC72EXLlg70ixdCqtWWSdS27Sxbo89du+yyHgTz+4zu1lxaAXLI5ezJ2oPTco0oX2F9rQO\nbE2BHAW8940on6KFPvXYtcsayG3cCIGBfx8/eukotb+ozbLnllGraK0UfU8t9G5mDISH/z1q37UL\nGje2Ru6tW0PR+7r6nLxykrVH1rLmyBrWHlmLfxZ/WpZrSavAVjQq1Ujn3FWitNCnLlOmWNM3mzZZ\nUzl3/LD/B95c8ya//fs38mTNk2Lvp4XeDW7ehLAwq7AvXWqtmrkzag8Kuvdq1WOXjrH++HrCjoYR\neiyUSzcv8WTpJ2lauinNyjajVJ5S3vo2VCqihT51McbqP1WkCHz88b2P/Xflf9l5Zicrnl9BtkxJ\nXNruIi30DygmBg4cgJ07rdH6zp3WWveqVf8u7lWrWr+It+NusydqD5tPbmbjiY1sOL6BW3G3eKLE\nEwSVDKJx6cZULlAZP9FNvJRrtNCnPpcuQY0aVuPBjh3/Ph5v4vnXj//i8s3L/PjMjynSd0oLvQuu\nXYM9e+4t6uHhULy49T+senXrv489BgF54zn410G2n97O9tPb2XZ6G7vP7KZMQBnqFK1D3eJ1aVCi\nAeXyltP2A+qBaaFPnbZsgXbtYOvWey98jImLodP8TuTKnItvOn7zwD2otNAnISrq3oK+a5fVdqBK\nlb8LevXq1i5NZL7G72d/Z2/UXnad2cXOMzvZe3YvBbIXoFbRWjz20GM8VuQxahapSe4suR2+t1Ku\n0kKfer3/vtXOOCzs3mtkbsTcoOXslpQJKMPnbT5/oJF9ui/08fHW+vX7i/rNm/cW9Bo14KFSVzhy\nOYID5w8Qfi6cfef2se/cPv68+ieVClSiWsFqVC9c/e4tJU+mKGWPFvrU6845vIcfhrFj733s6q2r\nPPv9s1yPuc53nb9L9jUyab7Qx8dbo/Njx+D4ceuW8P7hwxAQ8HdBr/jwNfKUPkJ01sMcvnCIyAuR\nRF6IJOJ8BJdvXaZ8vvJUzF+RKgWqULlAZSoXqEy5vOV0k2zlVVroU7dz56waNGMGNL9vt+24+DiG\nrRvG/H3zWdR1EdUKVXP59VN9ob9x4++ifX8RP3YMTp6EPHms+a8SJaBQiSvkKnacjPmOE5frODcy\nH+XMzaMcvXSUPy79wdVbVykdUJoyAWUIzBto3fIFUj5feYrlLqYnSpVP0kKf+oWFQdeu8Ntv1mqc\n+83eM5sBqwYwrMEw+tfu79LgMtUV+gEDzD0F/coV62Ro0VLR5C8VhX/RM2TN/yd+/n8Sm+001+Q0\nUdGnOHnlJCevnCTexFMyT0mK5y5OCf8SlMpT6p5b4ZyFtZirVEcLfdowejRMm2ZdUNWggXUrW/bv\n9sYR5yPot7wf56LPMaXVFOqXqO/U6/pUoReRFsBEwA+YcWd7wQSPmyffGwo5zhGb+RzRcpaLt88S\ndT2KmLgYCucsfPdWJFcRiuQqwkM5H6JY7mJ3b7mz5NZVLirN0UKfNhgDv/8O69fDhg3Wf2Nj4Ykn\nrKL/xBPw8MOG7w/MZ+DqgXc3D6pXvJ7duuYzhV5E/ICDQBPgNLAN6GqMOZAgxowOG02B7AUokKMA\nBXMUpFCOQhTMUdCnCnhoaGiq2LRc80xZ3szTlUKfGn6eqSFHcH+exlizFwkL/6lT1haFtepf40qZ\nmSw7/zF5svnTv3Z/2pZvS77s+f7xOskp9O6a16gNRBpjjhljYoC5QPv7g4Y3HM7Lj71Mp0qdeKLE\nEwTmC8Q/q7/PFHnggXdf9xTNM2VpniknNeQI7s9TxNoi9F//gs8/t67bOXwY/vMfuHE5J5s/foXT\nb0VwfelI3vt+ISU+KEO9Lxoy/tfxbDyxkeiY6GS/t7uWlxQFTiT4+iRW8VdKKWWTPz+0b2/dAK5f\n973jq9kAAATrSURBVGPLltasX9+asA032HzmZyIeXsaEkvO5nDmcUv5lkvU+uo5QKaV8RI4c8OST\n1g2yERPTil27WrFhA4RtuE3Y/t+Bmi6/rrvm6B8HQowxLWxfDwFMwhOyIqKnm5RSKhl85WRsBiAC\n62Tsn8BW4FljzP4UfzOllFJ2uWXqxhgTJyL9gdX8vbxSi7xSSnmB1y6YUkop5RleuWxURFqIyAER\nOSgig72RgyMiUkxE1onIPhHZKyKvejunpIiIn4jsEJHF3s4lKSLiLyLfich+28+0jrdzSoyI/FdE\nfheRPSIyW0QyezsnABGZISJRIrInwbEAEVktIhEiskpE/L2Zoy2nxPIcb/v/vktEvhcRr7d9TSzP\nBI+9LiLxIpLXG7ndl0uieYrIK7af6V4RGZvU8+/weKG3XUz1KdAcqAI8KyIVPZ2HE2KBgcaYKkBd\noJ+P5gnwGhDu7SQcmAQsN8ZUAh4BfG4qT0SKAK8AjxpjHsaa2uzq3azumon1O5PQEGCtMaYCsA4Y\n6vGs/imxPFcDVYwx1YFIfDdPRKQY0Aw45vGMEvePPEUkCGgLVDPGVAMmOHoRb4zonbqYytuMMWeM\nMbts969hFaai9p/lebZ/mK2A6d7OJSm2EVwDY8xMAGNMrDHmipfTSkoGIIeIZASyY13Z7XXGmA3A\nxfsOtwdm2e7PAjp4NKlEJJanMWatMSbe9uVmoJjHE7tPEj9PgI+ANz2cTpKSyPM/wFhjTKwt5ryj\n1/FGoU/sYiqfK6AJiUgpoDqwxbuZJOrOP0xfPtlSGjgvIjNtU0zTRCRlNtBMQcaY08AHwHHgFHDJ\nGLPWu1nZVdAYEwXWwAQo6OV8nNELWOHtJBIjIu2AE8aYvd7OxYHyQEMR2SwiP4vIY46eoK0dHRCR\nnMAC4DXbyN5niEhrIMr2yUNsN1+UEXgUmGyMeRSIxpp28CkikgdrlFwSKALkFJHnvJuVS3z5jz0i\nMgyIMcbM8XYu97MNPN4C3kl42EvpOJIRCDDGPA4MAuY7eoI3Cv0poESCr4vZjvkc28f3BcA3xphF\n3s4nEfWBdiJyBPgWaCwiX3s5p8ScxBopbbd9vQCr8PuapsARY8wFY0wc8ANQz8s52RMlIoUARKQw\ncNbL+SRJRHpgTTH66h/OskApYLeI/IFVl34TEV/8lHQC698mxphtQLyI/LP7WQLeKPTbgHIiUtK2\noqEr4KurRb4Ewo0xk7ydSGKMMW8ZY0oYY8rw/+3cP0rEQBzF8e8rBAsPYKuFvZ1gtyJYeQHBPxew\n8wpiL1hY2LlYiq1gaWGzinqBtfYAVs8is6K7rFaaOLxPFVKER8L8MpnJTHMfb2xvt51rXBleeJG0\nVE6t0c3J4yGwImlWzc56a3Rr0nj8q+0K2C3HO0BXOiNfcpYtyw+ATdtvraWa9JHT9pPteduLthdo\nOifLtrvw8hx/7pdAD6C0qRnbr99d4M8LfekpjRZTPQMXXVxMJWkV2AJ6kgZlbHmj7Vz/2D5wLume\n5q+bw5bzTLB9R/O1MQAeaBrXaauhCkl94BZYkjSUtAccAeuSRqvQf/zN7rdNyXkMzAHXpR2dtBqS\nqTk/Mx0YupmS8wxYlPQI9IEfO3dZMBURUblMxkZEVC6FPiKicin0ERGVS6GPiKhcCn1EROVS6CMi\nKpdCHxFRuRT6iIjKvQMHtDQXg9El0QAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x6b90a50>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"[(0, 19.0), (1, 29.0), (2, 45.0), (3, 70.0), (4, 108.0), (5, 164.0), (6, 244.0), (7, 356.0), (8, 498.0), (9, 656.0), (10, 785.0), (11, 811.0), (12, 669.0), (13, 389.0), (14, 135.0), (15, 21.0), (16, 1.0)]\n",
"0.1622\n"
]
}
],
"source": [
"spec, imort, rat = 5000, 0.003, 0.446142 # porcupines\n",
"lex = life_expectancy(spec, imort, rat)\n",
"plot.axvline(lex)\n",
"print(\"life expectancy:\", lex)\n",
"dat = list(life_table(spec, imort, rat))\n",
"plot1 = plot.plot(dat)\n",
"\n",
"x = numpy.linspace(0, len(dat)-1, 100)\n",
"plot2 = plot.plot(x, spec * pdf_of_doom(x, imort, rat))\n",
"\n",
"plot3 = plot.plot(x, 0 * x)\n",
"\n",
"plot.legend(plot1 + plot2 + plot3, [\"plot1\", \"plot2\", \"plot3\"])\n",
"\n",
"plot.show()\n",
"print(list(enumerate(dat)))\n",
"print(max(dat)/spec)\n",
"\n",
"# life expectancy: 9.4988\n",
"# [(0, 19), (1, 30), (2, 45), (3, 70), (4, 108), (5, 164), (6, 244), (7, 355), (8, 498), (9, 656), (10, 785), (11, 812), (12, 668), (13, 390), (14, 134), (15, 21), (16, 1)]\n",
"# 0.1624"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import sympy # hopefully in the actual program, once we switch to logistic, we won't need sympy\n",
"import sympy.abc as var\n",
"\n",
"def acceleration_of_death(life_exp, inf_mort):\n",
" return sympy.nsolve(-(sympy.exp(inf_mort / var.x) * sympy.Ei(-inf_mort / var.x)) / var.x - life_exp, 1)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.014622865595302821"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"float(acceleration_of_death(70, 0.0061))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Hooray, we got the same thing we got when we used WolframAlpha, which when plugged into `life_expectancy` gives us something really close to 70. This means that `sympy.nsolve` works as expected."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2.2940523520307313e-54"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"float(acceleration_of_death(15, 7 / 100)) # kittens!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Um what? That doesn't seem normal"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.0014571640933469553"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"float(acceleration_of_death(14, 7 / 100))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That's even stranger. If we reduce the life expectancy just by 1, the acceleration of death goes up by a factor of like 6e50. \n",
"\n",
"Graphing it in WolframAlpha helps to explain this a bit. It shows that the only solutions to the graph when the life expectancy is 15 are **negative**. I have avoided negative solutions for the acceleration of death in the past because they have always been wrong and the CDF seems to stop making sense.\n",
"\n",
"So, I would conclude that there are no _good_ solutions to that problem, but I don't understand why."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"[0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0,\n",
" 0.0]"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"list(itertools.islice(life_table(1000, 7 / 100, 2.2940523520307313e-54), 500))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And yeah, as expected, utterly terrible results from using the acceleration of death we got for a life expectancy of 15. 500 years later and **0** cats have died. I suppose some would consider that a good thing. Anyway, _maybe_ without rounding we might get a slightly better result, but not by much I don't think."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we look at WolframAlpha's graphs for these things, we find that the range of life expectancy for any given infant mortality rate appears to have a finite upper bound. I find this somewhat surprising as it seems from what I can remember from the graph of the CDF that you can have an acceleration of death that makes the graph take forever to get anywhere, making the life expectancy really high."
]
}
],
"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.1"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment