Mean Misanthrope Density.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Puzzle of the Misanthropic Neighbors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider [this puzzle](http://fivethirtyeight.com/features/can-you-solve-the-puzzle-of-your-misanthropic-neighbors/)  from [The Riddler](http://fivethirtyeight.com/tag/the-riddler/):\n",
    "    \n",
    "> *The misanthropes are coming. Suppose there is a row of some number, N, of houses in a new, initially empty development. Misanthropes are moving into the development one at a time and selecting a house at random from those that have nobody in them and nobody living next door. They keep on coming until no acceptable houses remain. At most, one out of two houses will be occupied; at least one out of three houses will be. But what’s the expected fraction of occupied houses as the development gets larger, that is, as N goes to infinity?*\n",
    "\n",
    "# Consider *N*=4 Houses\n",
    "\n",
    "![N houses](http://www.publicdomainpictures.net/pictures/180000/nahled/street-doodle.jpg)\n",
    "\n",
    "To make sure we understand the problem, let's try a simple example, with *N*=4 houses. We will represent the  originally empty row of four houses by four dots:\n",
    "\n",
    "     ....\n",
    "     \n",
    "Now the first person chooses one of the four houses (which are all acceptable). We'll indicate an occupied house by a `1`, so the four equiprobable choices are:\n",
    "\n",
    "     1...\n",
    "     .1..\n",
    "     ..1.\n",
    "     ...1\n",
    "     \n",
    "When a house is occupied, any adjacent houses become unacceptable. We'll indicate that with a `0`:\n",
    "\n",
    "     10..\n",
    "     010.\n",
    "     .010\n",
    "     ..01\n",
    "     \n",
    "In all four cases, a second occupant has a place to move in, but then there is no place for a third occupant. \n",
    "The occupancy is 2 in all cases, and thus the *expected occpancy* is 2. The occupancy fraction, or *density*, is 2/*N* = 2/4 = 1/2.\n",
    "     \n",
    "\n",
    "# Consider *N*=7 Houses\n",
    "\n",
    "With *N*=7 houses, there are 7 equiprobable choices for the first occupant:\n",
    "\n",
    "     10.....\n",
    "     010....\n",
    "     .010...\n",
    "     ..010..\n",
    "     ...010.\n",
    "     ....010\n",
    "     .....01\n",
    "     \n",
    "Now we'll add something new: the lengths of the **runs** of consecutive acceptable houses to the left and right of the chosen house:\n",
    "\n",
    "     10..... runs = (0, 5)\n",
    "     010.... runs = (0, 4)\n",
    "     .010... runs = (1, 3)\n",
    "     ..010.. runs = (2, 2)\n",
    "     ...010. runs = (3, 1)\n",
    "     ....010 runs = (4, 0)\n",
    "     .....01 runs = (5, 0)\n",
    "     \n",
    "# Defining `occ(n)`\n",
    "\n",
    "This gives me a key hint as to how to analyze the problem. I'll define `occ(n)` to be the expected number of occupied houses in a row of `n` houses. So:\n",
    "\n",
    "     10..... runs = (0, 5); occupancy = occ(0) + 1 + occ(5)\n",
    "     010.... runs = (0, 4); occupancy = occ(0) + 1 + occ(4)\n",
    "     .010... runs = (1, 3); occupancy = occ(1) + 1 + occ(3)\n",
    "     ..010.. runs = (2, 2); occupancy = occ(2) + 1 + occ(2)\n",
    "     ...010. runs = (3, 1); occupancy = occ(3) + 1 + occ(1)\n",
    "     ....010 runs = (4, 0); occupancy = occ(4) + 1 + occ(0)\n",
    "     .....01 runs = (5, 0); occupancy = occ(5) + 1 + occ(0)     \n",
    "     \n",
    "So we can say that `occ(n)` is:\n",
    "\n",
    "- 0 when `n` is 0 (because no houses means no occupants),\n",
    "- 1 when `n` is 1 (because one isolated acceptable house has one occupant),\n",
    "- else the mean, over the `n` choices for the first occupied house,\n",
    "of the sum of that house plus the occupancy of the runs to the left and right.\n",
    "\n",
    "We can implement that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from statistics import mean\n",
    "\n",
    "def occ(n):\n",
    "    \"The expected occupancy for a row of n houses (under misanthrope rules).\"\n",
    "    return (0 if n == 0 else\n",
    "            1 if n == 1 else\n",
    "            mean(occ(L) + 1 + occ(R)\n",
    "                 for (L, R) in runs(n)))\n",
    "\n",
    "def runs(n):\n",
    "    \"\"\"A list [(L, R), ...] where the i-th tuple contains the lengths of the runs\n",
    "    of acceptable houses to the left and right of house i.\"\"\"\n",
    "    return [(max(0, i - 1), max(0, n - i - 2))\n",
    "            for i in range(n)]\n",
    "\n",
    "def density(n): return occ(n) / n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check that `occ(4)` is 2, as we computed it should be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "occ(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that `runs(7)` is what we described above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, 5), (0, 4), (1, 3), (2, 2), (3, 1), (4, 0), (5, 0)]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "runs(7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check on `n = 7`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.3238095238095235"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "occ(7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.47482993197278905"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "density(7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That seems reasonable, but I don't know for sure that it is correct. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dynamic Programming Version of `occ`\n",
    "\n",
    "The computation of `occ(n)` makes multiple calls to `occ(n-1), occ(n-2),` etc. To avoid re-computing the same calls over and over, we will modify `occ` to save previous results using [dynamic programming](https://en.wikipedia.org/wiki/Dynamic_programming). I could implement that in one line with the decorator [`@functools.lru_cache`](https://docs.python.org/3/library/functools.html#functools.lru_cache), but instead I will explicitly manage a list, `cache`, such that `cache[n]` holds `occ(n)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def occ(n, cache=[0, 1]):\n",
    "    \"The expected occupancy for a row of n houses (under misanthrope rules).\"\n",
    "    # Store occ(i) in cache[i] for all as-yet-uncomputed values of i up to n:\n",
    "    for i in range(len(cache), n+1):\n",
    "        cache.append(mean(cache[L] + 1 + cache[R]\n",
    "                          for (L, R) in runs(i)))\n",
    "    return cache[n]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make sure this new version gets the same results as the old version:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "occ(4) == 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.47482993197278905"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "density(7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make sure the caching makes computation faster the second time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.23 s, sys: 115 ms, total: 6.34 s\n",
      "Wall time: 6.97 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "864.9617138385324"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time occ(2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 5 µs, sys: 0 ns, total: 5 µs\n",
      "Wall time: 9.06 µs\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "864.9617138385324"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time occ(2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting `density(n)`\n",
    "\n",
    "To get a feel for the limit of `density(n)`, start by drawing a plot over some small values of `n`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def plot_density(ns):\n",
    "    \"Plot density(n) for each n in the list of numbers ns.\"\n",
    "    plt.xlabel('n houses'); plt.ylabel('density(n)')\n",
    "    plt.plot(ns, [density(n) for n in ns], 's-')\n",
    "    return density(ns[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4353323288377046"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEPCAYAAAC3NDh4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGpRJREFUeJzt3X20XXV95/H3J7mXhgeJohBHIshD88TicY0RhzpcxUro\nWMNgR4nVFmtLlojDVMchuDrl0M6sgrZaW2inqRQfFppSQQlFIFPghFFRwkMIkEeF8BBJ0FoimgTy\n8J0/9r65Oyf33LPvzf2dc/Y9n9daZ+XsffbZ93c2l/O5399v799WRGBmZjacSZ1ugJmZdS+HhJmZ\nNeWQMDOzphwSZmbWlEPCzMyackiYmVlTSUNC0vWStkhaNcI2fyVpg6SVkk5L2R4zMxud1JXEDcC5\nzV6UdB5wQkT8KrAQ+D+J22NmZqOQNCQi4jvAv42wyXzgK/m2PwCmSpqWsk1mZlZep8ckjgaeLSxv\nyteZmVkX6HRImJlZF+vr8M/fBLyxsDw9X7cfSZ5kysxsDCJCY31vOyoJ5Y/hLAV+B0DSmcCLEbGl\n2Y7OPvtKIPZ7nH32lUREzzyuvLK3Pq+PhY+Fj8XYHwcqaSUh6WvAAPBaSc8AVwIHARERiyPi25J+\nQ9IPgV8CH07ZHjMzG52kIRERHyixzaVl9zdjxhSgxv33w+mnw5QpxfVmZjbeOj0mMSqLFy8C4Pjj\n4cYb4YQTOtygDhkYGOh0E7qGj8UQH4shPhbjR+PRZ9UOkmKwrTNnwtKl2b9mZtacJKLLB67HXV8f\n7NzZ6VaYmU18lQyJ/n7YtavTrTAzm/gqGRKuJMzM2qOSIdHf75AwM2uHyoaEu5vMzNKrZEi4u8nM\nrD0qGRKuJMzM2qOSIeFKwsysPSoZEh64NjNrj8qGhLubzMzSq2RIuLvJzKw9KhkSriTMzNqjkiHh\nSsLMrD0qGRIeuDYza4/KhoS7m8zM0qtkSLi7ycysPSoZEu5uMjNrj0qGRF+fu5vMzNqhkiHhSsLM\nrD0qGxKuJMzM0qtkSHjg2sysPSoZEu5uMjNrj0qGhAeuzczao5Ih4UrCzKw9KhsSriTMzNKrZEh4\n4NrMrD2Sh4SkeZLWSlov6fJhXn+1pFskPSrp+5LmtNqnu5vMzNojaUhImgRcC5wLnAQskDSrYbNP\nA49ExKnA7wJ/1Wq/7m4yM2uP1JXEXGBDRDwdETuBJcD8hm3mAPcARMQ64E2Sjhxpp+5uMjNrj9Qh\ncTTwbGH5uXxd0aPABQCS5gLHANNH2qkrCTOz9uiGgeurgddIehj4GPAIsHukN7iSMDNrj77E+99E\nVhkMmp6v2ysiXgJ+b3BZ0lPAk8PtrFarAfCjH8GWLQPAwHi21cys8ur1OvV6fdz2p4gYt53tt3Np\nMrAOOAd4HngAWBARawrbTAW2RcROSX8AnBURFw2zrxhs6733wp/8SfavmZk1J4mI0Fjfn7SSiIjd\nki4FlpF1bV0fEWskLcxejsXAbODLkvYATwAfadlodzeZmbVF6u4mIuJOYGbDur8rPP9+4+ut+DoJ\nM7P26IaB61HzBH9mZu1RyZBwJWFm1h6VDQlXEmZm6VUyJDxwbWbWHpUMCXc3mZm1RyVDwgPXZmbt\nUcmQcCVhZtYelQ0JVxJmZulVMiQ8cG1m1h6VDAl3N5mZtUclQ8ID12Zm7VHJkBgck0g4ga2ZmVHR\nkJBg8mTYPeKticzM7EBVMiTAg9dmZu1Q2ZDw4LWZWXqVDQkPXpuZpVfZkHAlYWaWnkPCzMyaqmxI\nuLvJzCy9yoaEKwkzs/QqGxKuJMzM0qtsSLiSMDNLzyFhZmZNVTYk3N1kZpZeZUPClYSZWXqVDglX\nEmZmaVU2JDzBn5lZepUNCXc3mZmlV9mQ8MC1mVl6yUNC0jxJayWtl3T5MK8fLmmppJWSHpN0UZn9\nupIwM0svaUhImgRcC5wLnAQskDSrYbOPAU9ExGnA24G/kNTXat8euDYzSy91JTEX2BART0fETmAJ\nML9hmwBelT9/FfCvEdHy698D12Zm6aUOiaOBZwvLz+Xriq4F5kj6MfAocFmZHbu7ycwsvZbdOm1w\nLvBIRLxD0gnA/5V0SkT8onHDWq229/lPfjLArl0DbWukmVkV1Ot16vX6uO1PETFuO9tv59KZQC0i\n5uXLi4CIiGsK2/wz8GcR8d18+W7g8oh4sGFfUWzrRz8KJ58Ml1ySrPlmZpUniYjQWN+furtpBXCi\npGMlHQRcCCxt2OZp4J0AkqYBM4AnW+3Y3U1mZukl7W6KiN2SLgWWkQXS9RGxRtLC7OVYDPwv4EuS\nVuVv+x8R8bNW+/Z1EmZm6SUfk4iIO4GZDev+rvD8ebJxiVFxJWFmlp6vuDYzs6YqGxKuJMzM0nNI\nmJlZU5UNCXc3mZmlV9mQcCVhZpZeZUPClYSZWXqVDQlXEmZm6TkkzMysqcqGhLubzMzSG1VISDpU\n0uRUjRkNVxJmZumNGBKSJkn6gKTbJb0ArAWel7Ra0mclndieZu7PlYSZWXqtKol7gROAK4DXR8Qb\nI+Io4NeA7wPXSPpg4jYOy5WEmVl6rSb4e2d+29F95LO03gzcLKk/SctacEiYmaU3YkgUAyIfi5hW\nfE9EPDNciLSDu5vMzNIrNVW4pI8DVwJbgD356gBOSdSullxJmJmlV/Z+EpcBMyPiX1M2ZjQcEmZm\n6ZU9BfZZYGvKhoyWu5vMzNIrW0k8CdQl3Q68PLgyIj6XpFUluJIwM0uvbEg8kz8Oyh8d50rCzCy9\nUiEREVelbshouZIwM0uv1RXXfy/p5CavHSrp9yT9dpqmjcwhYWaWXqtK4jrgf+ZB8TjwE2AK8KvA\n4cA/ADcmbWET7m4yM0uv1cV0K4H3SToM+PfAvwO2A2siYl0b2teUKwkzs/TKDly/Hbg9Iva03LJN\nXEmYmaVX9jqJ9wMbJH1G0qyUDSrLlYSZWXqlQiIiPgicDvwI+JKk+yVdLOlVSVs3AoeEmVl6pW86\nFBE/B74BLCEbm/jPwMP5vE5t5+4mM7P0SoWEpPmSvgnUgX5gbkScB5wKfDJd85pzJWFmll7ZSuIC\n4PMRcXJEfDYiXgCIiG3AR0Z6o6R5ktZKWi/p8mFe/++SHpH0sKTHJO2S9OpWDXIlYWaWXtmQ2BwR\n9xVXSLoGICLubvYmSZOAa4FzgZOABY0D3xHx5xFxekScQXYHvHpEvNiqQZMnw5492cPMzNIoGxK/\nPsy680q8by6wISKezm9OtASYP8L2C4Cvl2mQ5C4nM7PUWk3L8VFJjwGzJK0qPJ4CVpXY/9Fk04wP\nei5fN9zPOhiYR3Zb1FLc5WRmllari+m+BtwB/BmwqLD+pfw+1+PpN4HvjNTVVKvV9j4fGBigv3/A\nlYSZWUG9Xqder4/b/hQRzV+UDo+In0s6YrjXWwWFpDOBWkTMy5cXZW+La4bZ9hbgpohY0mRf0djW\n174W1q6FI48cqRVmZr1LEhGhsb6/TCXxbuAhsntaF39QAMe3eP8K4ERJxwLPAxeSjTvsQ9JU4Gxg\nVDPK9ve7u8nMLKVWE/y9O//3uLHsPCJ2S7oUWEY2/nF9RKyRtDB7ORbnm54P3BUR20ezfw9cm5ml\nNWJ3096NpLOAlRHxS0kfBM4A/jIinkndwEIb9utuOu44uPtuOL5VPWNm1qMOtLup7CmwfwtskzR4\nhfWPgK+O9YeOF1cSZmZplQ2JXfmf8fOBayPiOqBjk/sN6utzSJiZpVT2fhIvSboC+CDwH/MrqfvT\nNascD1ybmaU1mvtJvAx8JCI2A9OBzyZrVUnubjIzS6tUJZEHw+cKy88AX0nVqLJ8xbWZWVplpwq/\nQNIGSVsl/VzSS5J+nrpxrbiSMDNLq+yYxGeA34yINSkbM1oeuDYzS6vsmMSWbgsI8MC1mVlqZSuJ\nByX9I/AtsgFsACLiliStKsndTWZmaZUNicOBbcC7CusC6GhIeODazCytsmc3fTh1Q8bClYSZWVpl\nz26aIeluSY/ny6dI+qO0TWvNIWFmllbZgeu/J7v/9E6AiFhFNu13R7m7ycwsrbIhcUhEPNCwruNf\nz64kzMzSKhsSP5V0AtlgNZJ+i+wmQh3l6yTMzNIqe3bTx4DFwCxJm4CnGOVd5FLwdRJmZmmNGBKS\nPlFY/DZwL1n18UvgvRTmc+oEdzeZmaXVqpIYvGfETODNwK1k97n+ENA4RtF2Hrg2M0ur1T2urwKQ\ndB9wRkS8lC/XgNuTt64FVxJmZmmVHbieBrxSWH4lX9dRHrg2M0ur7MD1V4AHJH0zXz4f+FKSFo2C\nKwkzs7TKTsvxvyXdAbwtX/XhiHgkXbPK6e+Hbds63Qozs4mrbCVBRDwMPJywLaPmgWszs7TKjkl0\nJXc3mZmlVemQ8MC1mVlalQ4JX3FtZpZW5UPClYSZWTqVDgkPXJuZpZU8JCTNk7RW0npJlzfZZkDS\nI5Iel3Rv2X27kjAzS6v0KbBjIWkScC1wDvBjYIWkWyNibWGbqcB1wLsiYpOk15XdvweuzczSSl1J\nzAU2RMTTEbETWALMb9jmA8DNEbEJICJ+WnbnHrg2M0srdUgcDTxbWH4uX1c0AzhC0r2SVkj6UNmd\nu7vJzCytpN1NJfUBZwDvAA4F7pd0f0T8sHHDWq229/nAwAB9fQMOCTOzgnq9Tr1eH7f9KSLGbWf7\n7Vw6E6hFxLx8eREQEXFNYZvLgSmFacm/CNwRETc37Csa23rPPfCnfwr3lh7qNjPrLZKICI31/am7\nm1YAJ0o6VtJBwIXA0oZtbgV+TdJkSYcAbwHWlNm5B67NzNJK2t0UEbslXQosIwuk6yNijaSF2cux\nOCLWSroLWAXsBhZHxOoy+/fAtZlZWkm7m8bTcN1NDz4ICxfCQw91qFFmZl2u27ubknJ3k5lZWpUO\nCXc3mZmlVemQcCVhZpZWpUPClYSZWVqVDwlXEmZm6VQ6JNzdZGaWVqVDwt1NZmZpVT4kXEmYmaVT\n6ZDwnenMzNLqhllgx2ywkrj44qtZv37Hfq/PmDGFxYsXdaBlZmYTQ6VDYrCSWLduB/fdVxtmi+HW\nmZlZWZXubpo0KXuYmVkalf+K7e+HPXs63Qozs4mp8iHR1wcVmcjWzKxyKh8S/f0OCTOzVCo9cA1Z\nJXHccVPo76/xwANZYBx/PBx5ZHZ2k5mZjV2lbzoE8IY3wIoVWShMnQqXXJL9+8d/3IFGmpl1mZ6+\n6RAMXSuxYQMccwycdhqsKXWHbDMza6XyITE4yd+aNTBnDsyeDatL3SHbzMxaqXxIDE7yt3p1FhCz\nZmVVxe7dnW6ZmVn1VT4kGiuJww7Lxic2bux0y8zMqq/yIdFYSYC7nMzMxsuECImXX866mGbNytbN\nmePBazOz8VD5kOjrg/Xr4aij4NBDs3WuJMzMxkflQ6K/H1atGupqAlcSZmbjpfIh0deXhcScOUPr\nZs/OQqIi1wmamXWtyofEcJXEEUfAIYfApk2da5eZ2UQwIULihRf2rSTA4xJmZuMheUhImidpraT1\nki4f5vWzJb0o6eH88Uej2X9fPkVhsZIAj0uYmY2HpLPASpoEXAucA/wYWCHp1ohY27DpfRHxnrH8\njP5+mDYNXvOafdfPng2PPTaWPZqZ2aDUU4XPBTZExNMAkpYA84HGkBj1DIUXX3w169fvYPVqeOUV\nGBjI1g9OD/6DH+xg40ZYt27oPTNmTGHx4kWj/hBmZr0qdUgcDTxbWH6OLDgavVXSSmAT8KmIaDma\nsH79DpYvr+1dXr588Fm2btWqWsP6odfMzKycbrjp0EPAMRGxTdJ5wLeAGeO3+6uBHQCsXLmRgYEa\n4KrCzKyM1CGxCTimsDw9X7dXRPyi8PwOSX8j6YiI+Fnjzmq12t7nL764sWQTdjBYQWzdun/FYWY2\nkdTrder1+rjtL3VIrABOlHQs8DxwIbCguIGkaRGxJX8+l+xuefsFBOwbEvV6bbhNzMx62sDAAAOD\ng7TAVVdddUD7SxoSEbFb0qXAMrLTba+PiDWSFmYvx2LgtyR9FNgJbAfen7JNZmZWXmXvcT14dlOj\nwbObBl9buXIjW7d+KX91aHxi6tSNnHbam/a+x+MTZjYRHeg9rrth4HpMyn6pDwzUCuMQHp8wMxuN\nyoZEWVllUQMGq4qONsfMrFIq2900FllVUSusybqfil1P4O4nM5s4era7aXxk3U/7dj2Bu5/MzDI9\nHhJFvujOzKxRT4VEcXwCGscoioPaV7N8+WBgrGX9+tre9zswzKyX9FRINH7B73vmU5HPgjIzgx4L\nibHLuqKK3VDgysLMJr6eDonyp8cWB7iLYxfuijKzia2nQ6L4pd6866lR87GL2277L2zffigHHzyJ\nmTOH5jV0gJhZVfV0SBSN7aK7fccutm6tMVhxbN48VHF873t1bropu89SMUAcHmbW7RwSueKXdTYv\nVG3v8tiu1B4KkJ07a3mAsE+AFKsPcICYWfdxSAyj/FlQYzU0xjFYfUDWfbV58/DdV/CSw8TM2s4h\nUUL75n9q3n2VrR96rVk18sorWzjooGkA+4SLg8bMxsIhUUKzrqjOThg4fDUCNbZvz54Xw2X/oClX\ntZQNnZG28yC+WXU5JEZppLGLdeueYPv2izj44Els376ny2ecLVe1lA2dkbYbaRB/LKEzlu3Gax+u\nvqzXOCQOwEhfEsUAGQwPoCIBMt6aD+KPJXTGst347KN91Vc3bNeNberksejVPxAcEomMJUC2bdvB\nzp1taJyNUfuqr27Yrhvb1LljUf6kkm4IuGIX74FySHRAswAZqfsq+wW4COjVasSsk8qfVNINAbd5\nc7HtV43ys+7LIdFFRlPCNqtGsr8usufFcGkMGlctZlaGQ6KiDrRPtGzVUjZ0RtrOVY9ZdTkkelQ7\nB95GGsQfS+iMZbvx2IerL+tFPXWPa7MDkYXdjr3LWeB1zxk4Prsp3XbZHwhLCr8NNYpjA0PPR3qt\nndsVHdg9rh0SZmYtlP0DoVsCrnh20/LlVzkkzMxseNKBVRKTxrMxZmY2sTgkzMysKYeEmZk1lTwk\nJM2TtFbSekmXj7DdmyXtlHRB6jaZmVk5SUNC0iTgWuBc4CRggaRZTba7GrgrZXsminq93ukmdA0f\niyE+FkN8LMZP6kpiLrAhIp6OiJ3AEmD+MNt9HPgG8ELi9kwI/h9giI/FEB+LIT4W4yd1SBwNPFtY\nfi5ft5ekNwDnR8TfAmM+TcvMzMZfNwxc/yVQHKtwUJiZdYmkF9NJOhOoRcS8fHkREBFxTWGbJwef\nAq8DfglcHBFLG/blK+nMzMaga6+4ljQZWAecAzwPPAAsiIg1Tba/AbgtIm5J1igzMyst6SywEbFb\n0qXAMrKuresjYo2khdnLsbjxLSnbY2Zmo1OZuZvMzKz9umHguqWyF+RNRJKmS7pH0hOSHpP0X/P1\nr5G0TNI6SXdJmtrptraDpEmSHpa0NF/u1eMwVdI/SVqT/268pYePxR9KelzSKkk3Sjqol46FpOsl\nbZG0qrCu6eeXdIWkDfnvzrta7b/rQ6LsBXkT2C7gExFxEvBW4GP5518E/EtEzATuAa7oYBvb6TJg\ndWG5V4/DF4BvR8Rs4FRgLT14LPJT6D8OnBERp5B1oS+gt47FDWTfj0XDfn5Jc4D3AbOB84C/kTTi\noHbXhwTlL8ibkCJic0SszJ//AlgDTCc7Bl/ON/sycH5nWtg+kqYDvwF8sbC6F4/D4cDbIuIGgIjY\nFRFb6cFjkZsMHCqpDzgY2EQPHYuI+A7wbw2rm33+9wBL8t+ZjcAGsu/YpqoQEi0vyOsVkt4EnAZ8\nH5gWEVsgCxLgqM61rG0+D3yKfU9w6MXjcBzwU0k35F1viyUdQg8ei4j4MfAXwDNk4bA1Iv6FHjwW\nDY5q8vkbv0830eL7tAohYYCkw8imLrksrygazziY0GcgSPpPwJa8qhqpPJ7QxyHXB5wBXBcRZ5Bd\nW7SIHvudAJD0arK/mo8F3kBWUfw2PXgsWhjz569CSGwCjiksT8/X9Yy8jP4G8NWIuDVfvUXStPz1\n1zPx5706C3hPfvHl14F3SPoqsLnHjgNk1fSzEfFgvnwzWWj02u8EwDuBJyPiZxGxG/gm8B/ozWNR\n1OzzbwLeWNiu5fdpFUJiBXCipGMlHQRcCCxt8Z6J5h+A1RHxhcK6pcBF+fPfBW5tfNNEEhGfjohj\nIuJ4st+BeyLiQ8Bt9NBxAMi7EZ6VNCNfdQ7wBD32O5F7BjhT0pR8APYcshMbeu1YiH0r7Gaffylw\nYX4G2HHAiWQXOTffcRWuk5A0j+xsjsEL8q7ucJPaRtJZwH3AY2QlYwCfJvsPexPZXwVPA++LiBc7\n1c52knQ28MmIeI+kI+jB4yDpVLIB/H7gSeDDZAO4vXgsriT7w2En8Ajw+8Cr6JFjIelrwADwWmAL\ncCXwLeCfGObzS7oC+AjZ8bosIpaNuP8qhISZmXVGFbqbzMysQxwSZmbWlEPCzMyackiYmVlTDgkz\nM2vKIWFmZk05JMwa5HMiXdDpdph1A4eEmZk15ZCwnpFP7bI6nzX1cUl3SvqVJpufLem7kn5YrCok\nfTa/+dOjkt6Xrztb0m2Fbf5a0u/kz6/Of9ZKSZ/J171O0jck/SB/vLWwn0fymV0fknRosoNhVlLS\ne1ybdaETgfdHxMWS/hF4L/C1YbZ7fUScJWk22Xw3t0h6L3BKRJws6ShghaTl+fb7TV2QTxlyfkTM\nypcPz1/6AvC5iPiepDcCdwFzgE8Cl0TE/fnU3zvG7VObjZFDwnrNUxHxWP78IeBNTbb7FkBErMkD\nAbKZaL+er39BUh14M/BSk31sBbZL+iJwO/DP+fp3ArMLdwQ7LA+F7wKfl3QjcEtE9NRsx9ad3N1k\nveblwvPdNP9Dqbhds/tXDK7fRTa53qApAPnU1XPJpnl/N3Bn4X1viYjT88cxEbEtIq4hm3jtYOC7\nhVlezTrGIWG9ZsT7+bZ4z/8D3i9pkqQjgbeRzcb7NFll0J/fBOccgLw6eHVE3Al8Ajgl388ysnt1\nk293av7v8RHxRER8hmyK/F66l7t1KXc3Wa8pM+3xsHc1i4hvSjoTeBTYA3wqIl4AkHQT8DjwFPBw\n/r7DgVslTcmX/zD/9zLgOkmPklUg9wGXAP9N0tvJKpwngDtG//HMxpenCjczs6bc3WRmZk05JMzM\nrCmHhJmZNeWQMDOzphwSZmbWlEPCzMyackiYmVlTDgkzM2vq/wNokhP4XD3+tQAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1043b8eb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_density(range(1, 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is something funny going on with the first few values of `n`. Let's separately look at the first few:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.46203174603174607"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEPCAYAAABGP2P1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHbxJREFUeJzt3Xu0XHV99/H3JwkhJJAgQSNyFTEQKCIXI4KGoWAJPBYo\nbbmoD0VYhaJYnnpLcNXFca22EPXx0gdsTaVUqBgFuVVQoMrBCyCBcM8JQUMTwi3IPYSEEL7PH3uf\nZOZkzpnrPnvPzOe11qyZ2fObPd9z1jnnc/b+7v3bigjMzMwGjcm7ADMzKxYHg5mZVXAwmJlZBQeD\nmZlVcDCYmVkFB4OZmVXINBgkXSLpGUkPjDDmnyU9Kuk+Se/Nsh4zM6st6y2GS4GjhntR0tHAuyLi\n3cBZwL9mXI+ZmdWQaTBExK+BF0YYchxwWTr2t8AUSdOyrMnMzEaWd49hR+DxsudPpMvMzCwneQeD\nmZkVzLicP/8JYOey5zulyzYjyZM6mZk1ISLUyPjR2GJQeqvmeuBUAEkHAy9GxDPDryqA4LDDzici\ncr+df34x6uiEulyTa+qFuopYUzMy3WKQdAVQAqZKWgGcD4wHIiLmR8SNko6R9DvgVeATWdZjZma1\nZRoMEfHROsacU+/6dt21j912g+nTJ7RUl5mZDS/vHkND3vvePq69Nu8qNimVSnmXUFUR63JN9XFN\n9StiXUWsqRlqdh/UaJMUO+4YrFyZdyVmZp1DElHA5nPbrFkDz4zQmjYzs9Z1VDAccADcc0/eVZiZ\ndbeOCoaDDnIwmJllraOC4cADHQxmZlnruGC4++68qzAz624dFQzvfKcb0GZmWeuoYJDcgDYzy1pH\nBQO4AW1mlrWOCwY3oM3MstWRweAGtJlZdjouGNyANjPLVscFgxvQZmbZ6rhgADegzcyy1JHB4Aa0\nmVl2OjYY3IA2M8tGRwaDG9BmZtnpyGBwA9rMLDsdGQzgBrSZWVY6NhjcgDYzy0ZHB4Mb0GZm7dex\nweAGtJlZNjo2GNyANjPLRscGA7gBbWaWhY4OBjegzczar+ODwQ1oM7P26uhgcAPazKz9Mg8GSbMl\nLZG0VNKcKq9vK+lqSfdLulPS3vWv2w1oM7N2yzQYJI0BLgKOAvYBTpG015BhXwTujYj9gL8C/rmR\nz3AD2sysvbLeYpgJPBoRyyNiPbAAOG7ImL2BXwBExCPAbpLeWu8HuAFtZtZeWQfDjsDjZc9XpsvK\n3Q+cACBpJrALsFO9H+AGtJlZexWh+Xwh8BZJi4BPAfcCG+p9sxvQZmbtNS7j9T9BsgUwaKd02UYR\n8Qpw+uBzSY8By6qtrK+vb+PjUqlEqVSqaEAfc0wbKzcz60D9/f309/e3tA5FRHuqqbZyaSzwCHAE\n8BRwF3BKRAyUjZkCrImI9ZL+Gjg0Ik6rsq4Yrta5c2HSJPjSlzL4IszMOpgkIkKNvCfTXUkRsQE4\nB7gZeBhYEBEDks6SdGY6bAbwkKQBkqOXzm30c9yANjNrn0y3GNpppC2GZctg1ixYuXKUizIzK7jC\nbTGMFjegzczapyuCwWdAm5m1T1cEA/gMaDOzdumaYHAD2sysPboqGHwGtJlZ67omGNyANjNrj64J\nBjegzczao2uCAdyANjNrh64KBjegzcxa13XB4Aa0mVlruioY3IA2M2tdVwWDG9BmZq3rqmAAN6DN\nzFrVdcHgBrSZWWu6MhjcgDYza17XBYMb0GZmrem6YHAD2sysNV0XDOAGtJlZK7oyGNyANjNrXtcG\ngxvQZmbN6cpgcAPazKx5XRkMbkCbmTWvK4MB3IA2M2tW1waD+wxmZs3p6mDwFoOZWeO6NhjcgDYz\na07XBoMb0GZmzenaYAA3oM3MmtHVweAGtJlZ4zIPBkmzJS2RtFTSnCqvT5Z0vaT7JD0o6bR2fbYb\n0GZmjVNEZLdyaQywFDgCeBJYCJwcEUvKxpwHTI6I8yRtDzwCTIuIN4asKxqtNQKmToWBAZg2rcUv\nxsysA0kiItTIe7LeYpgJPBoRyyNiPbAAOG7ImAC2SR9vAzw3NBSa5Qa0mVnjsg6GHYHHy56vTJeV\nuwjYW9KTwP3Aue0swA1oM7PGjMu7AOAo4N6I+GNJ7wJukfSeiFg9dGBfX9/Gx6VSiVKpVHPlBx4I\n//mf7SvWzKzI+vv76e/vb2kdWfcYDgb6ImJ2+nwuEBExr2zMT4ALIuI36fOfA3Mi4u4h62q4xwCw\nbBnMmgUrV7bwhZiZdagi9hgWAntI2lXSeOBk4PohY5YDRwJImgZMB5a1qwCfAW1m1phMgyEiNgDn\nADcDDwMLImJA0lmSzkyH/QNwiKQHgFuAL0TE8+2qwQ1oM7PGZLorqZ2a3ZUEMHcuTJoEX/pSm4sy\nMyu4Iu5KKgSfAW1mVr+eCQbvSjIzq09PBIMb0GZm9euJYHAD2sysfj0RDOAzoM3M6tUzweAGtJlZ\nfXoqGLzFYGZWW88EgxvQZmb16ZlgcAPazKw+DQWDpEmSxmZVTNbcgDYzq23EYJA0RtJHJd0gaRWw\nBHhK0mJJX5W0x+iU2R5uQJuZ1VZri+FW4F3AecDbI2LniHgb8EHgTmCepI9nXGPbuAFtZlbbiJPo\nSdoivSRnS2PaoZVJ9Ab5GtBm1mvaPole+R98SWMlvUPSLoO3oWOKzg1oM7Pa6mo+S/o08AzJ9RJu\nSG8/ybCuzLgBbWY2snqv+XwusGdEPJdlMaPB14A2MxtZvYerPg68lGUho8UNaDOzkdW7xbAM6Jd0\nA7BucGFEfD2TqjJUfga0G9BmZpurd4thBUl/YTywTdmt47gBbWY2srq2GCLiy1kXMpoGG9DHHJN3\nJWZmxVPrzOd/k7TvMK9NknS6pI9lU1p2fAa0mdnwap3g9l7gi8C+wEPAs8AE4N3AZODfgX+NiHXD\nrqRdhbbhBLdBy5bBrFmwcmVbVmdmVljNnOA2YjCUrXhr4CBgB+A1YCAiHmmqyia1Mxh8BrSZ9Ypm\ngqHeo5IOB26IiDcbL6t4yhvQ7jOYmVWq96ikk4BHJX1F0l5ZFjRafAa0mVl1dQVDRHwc2B/4PfAf\nku6QdKakjjxkFdyANjMbTt0X6omIl4GrgAUkvYY/Axal8yh1HJ8BbWZWXb2T6B0n6RqgH9gCmBkR\nRwP7AZ/Nrrzs+BrQZmbV1bvFcALwjYjYNyK+GhGrACJiDXDGSG+UNFvSEklLJc2p8vrnJN0raZGk\nByW9IWnbhr+SBvkMaDOz6uoNhqcj4pflCyTNA4iInw/3JkljgIuAo4B9gFOGNq8j4msRsX9EHEBy\npbj+iHixga+haW5Am5ltrt5g+HCVZUfX8b6ZwKMRsTy9oM8C4LgRxp8C/KDOmlrmBrSZ2eZqTYlx\ntqQHgb0kPVB2ewx4oI7170gyZfeglemyap+1FTAb+HF9pbfODWgzs83VOsHtCuCnwAXA3LLlr0TE\n822u5U+BX4+0G6mvr2/j41KpRKlUaukDPQW3mXWb/v5++vv7W1pHrbmSJkfEy5K2q/Z6rXCQdDDQ\nFxGz0+dzk7fFvCpjrwZ+FBELhllX26bEKHfkkfCZz/gMaDPrTs1MiVGrx3BFen8PcHd6f0/Z81oW\nAntI2lXSeOBk4PqhgyRNAQ4Drquz7rZxA9rMrNKIu5Ii4iPp/TubWXlEbJB0DnAzSQhdEhEDks5K\nXo756dDjgZsi4rVmPqcVvga0mVmlemdXPRS4LyJelfRx4ADgmxGxIusCy2rIZFeSp+A2s26Wxa6k\nQf8CrJE0eKbz74HLG6yvkHwGtJlZpXqD4Y303/XjgIsi4mI69JrPQ/kMaDOzSvUGwyuSzgM+DtyQ\nntG8RXZljS43oM3MNmnkegzrgDMi4mlgJ+CrmVU1ynwGtJnZJnU1n4sgq+YzuAFtZt0rs+azpBMk\nPSrpJUkvS3pF0svNlVk8bkCbmW1S766krwDHRsSUiJgcEdtExOQsCxtNbkCbmW1SbzA8ExEDmVaS\nMzegzcwStSbRG3S3pB8C15I0oQGIiKszqSoHPgPazCxR75nPl1ZZHBFxevtLGraGzJrP4Aa0mXWn\nZprPPiopFQFTp8LAgKfgNrPukeVRSdMl/VzSQ+nz90j6+2aKLCo3oM3MEvU2n/+N5HrM6wEi4gGS\nKbS7ykEH+UQ3M7N6g2FiRNw1ZNkb7S4mb77Up5lZ/cHwB0nvAgJA0l8AT2VWVU4cDGZm9R+VtDsw\nHzgEeAF4DPhYRCzPtryKGjJtPoMb0GbWfdp+VJKkzwxZtBXJVsarABHx9UaLbNZoBAP4GtBm1l2y\nOCppm/R2EHA28BZgW+BvSK7i1nXcgDazXlfrms9fBpD0S+CAiHglfd4H3JB5dTnwGdBm1uvqbT5P\nA14ve/56uqzruAFtZr2u3rmSLgPuknRN+vx44D8yqShn5VNwuwFtZr2ori2GiPhH4BMkRyS9AHwi\nIi7IsrC8+AxoM+t19W4xEBGLgEUZ1lIYgw1oH5lkZr2o3h5DT3Gfwcx6mYOhCgeDmfUyB0MVvga0\nmfUyB0MVbkCbWS9zMAzDZ0CbWa/KPBgkzZa0RNJSSXOGGVOSdK+khyTdmnVN9XCfwcx6VaaX9pQ0\nBlgKHAE8CSwETo6IJWVjpgC3A38SEU9I2j4i/lBlXaMyid4gXwPazLpBZpf2bMFM4NGIWB4R64EF\nwHFDxnwU+HFEPAFQLRTy4Aa0mfWqrINhR+Dxsucr02XlpgPbSbpV0kJJ/zvjmuriBrSZ9aq6z3zO\n0DiSKbz/GJgE3CHpjoj43dCBfX19Gx+XSiVKpVKmhfkMaDPrNP39/fT397e0jqx7DAcDfRExO30+\nF4iImFc2Zg4woWyK7+8CP42IHw9Z16j2GACuvDKZgvu660b1Y83M2qaIPYaFwB6SdpU0HjgZuH7I\nmOuAD0oaK2ki8H5gIOO66uIjk8ysF2W6KykiNkg6B7iZJIQuiYgBSWclL8f8iFgi6SbgAWADMD8i\nFmdZV708BbeZ9aJMdyW1Ux67ksDXgDazzlbEXUkdz2dAm1mvcTDU4D6DmfUaB0MNDgYz6zUOhhp8\nBrSZ9RoHQw0+A9rMeo2DoQ5uQJtZL3Ew1MF9BjPrJQ6GOjgYzKyXOBjq4Aa0mfWSIsyuWnjlDeii\nngF95pkXsnTp2s2WT58+gfnz5+ZQkZl1KgdDnYo+BffSpWu57ba+Kq9UW2ZmNjzvSqqT+wxm1isc\nDHUqcjC89hosX553FWbWLRwMdSpiA/rNN+Hyy2HPPWH16upjli9PgsPMrF4OhjoV7QzoW29N+h7f\n/jYsWAD77FN93OrVSXBcfnkSJGZmtbj53IAiNKAHBuALX4CHH4YLL4S//MsktKZPn0C1RvP06RM4\n7TT47Gfhm9+Er30NDj98tKs2s07iC/U0IM9rQK9aBX19SQ3nnQef+hRsuWX9749I3jt3LvzRH8G8\neTBjRmblmllB+EI9GcujAf3aa/BP/wR7750EwZIlyRXlGgkFSLYqTjwx2eI47DCYNQs++ckkcMzM\nyjkYGjCaDejyxvK998Kdd8I3vgFTp7a23i23THYrLVmSPN57b7jgAjeozWwTB0MDRqsBPbSxfOWV\nsMce7f2MqVOToLnzTli0yA1qM9vEwdCgLKfgHhiAY4+FM85IegG33w6HHJLNZw3aY48keBYsSILo\nfe9LgsnMepeDoUFZ9BlWrUr298+aBaVSEhAnnphsoYyWQw5JgmjOnCSYjj022d1kZr3HwdCgdgbD\na68l+/dbbSy3y9AG9Yc+lBz95Aa1WW9xMDSoHQ3o8sbyokXtayy3S3mDevx4N6jNeo2DoUGtNqBv\nvTXZj59lY7ld3KA2600OhiY004AubyzPmTM6jeV2qdag7u/Puyozy4qDoQmN9BmK0Fhul/IG9emn\nu0Ft1q0cDE2oJxiK1lhuFzeozbpf5sEgabakJZKWSppT5fXDJL0oaVF6+/usa2rVSA3oojeW28UN\narPulWkwSBoDXAQcBewDnCJprypDfxkRB6S3f8iypnYYrgHd3985jeV2cYParPtkOruqpIOB8yPi\n6PT5XCAiYl7ZmMOAz0XEn9ZYV+6zqwKceeaFLF26lmXLYMwY2G23wa2HCYwdO7diKuxedPvtyZbE\n66/DDjtcyOrVazcbM336BObPn5tDdWa9p5nZVbO+HsOOwONlz1cCM6uM+4Ck+4AngM9HxOKM62ra\n0qVrue22vo3PBy+pufvufSxe3Pk9hFYNNqivvBJOPXUt69b1VRlVbZmZFUURLtRzD7BLRKyRdDRw\nLTA955oatvPODoVBgw3qiy6CX/1q89cXL4azz4a3vGXk2+TJ7d/yGtziG8pbMWabZB0MTwC7lD3f\nKV22UUSsLnv8U0nflrRdRDw/dGV9fX0bH5dKJUqlUrvrtTYaM0wHa+pU2HdfeOGFpIG/ZEnyeOht\nzRqYMmVTUGy3Xe0wqRUqQ7f4Nqm2zKzz9Pf309/iiUZZB8NCYA9JuwJPAScDp5QPkDQtIp5JH88k\n6XtsFgpQGQzWuaZNS87tqGX9enjxxeqhMVKoPP98cnRUeagMBsvSpdU/y81y6xZD/2n+8pe/3PA6\nMg2GiNgg6RzgZpIjoC6JiAFJZyUvx3zgLySdDawHXgNOyrIm6xxbbAFvfWtya9RwoTLc+Se/+hXs\ntBPsskty23XXzR9PmdK7BxVYb8m8xxARPwP2HLLsO2WPLwYuzrqOdpk+fQLVdjsky61cnt+r4ULl\nO9+B3/1u8/GzZiWH2a5YkdyWL4eHHoIbb0weL1+ehEK1wBh8vMMOMK6J3yj3PaxoitB87ij+Ra1f\nJ32vpE1/5KuJSLZAyoNjxQq4775Nj599NgmHwcCoFiJbb735ut33sKJxMFhPaXYrRtrUq9hvv+pj\nXn8dVq6sDI+774arr94UHhMmbB4Ynk7EisbBYD0ly62Y8eNh992TWzUR8Ic/bL7V8eyz1cffcw98\n+MP1H5E1efLwR4I1yru3epuDwWyUSJv6HgceuGn5okVw222bj3/3u+Fzn6tsnj/7bHJkVbWjtF59\nNQmHeg/pLb9NmVIZKt691dscDGYFNXkyHHVU/ePfeANeemn4w3ufey5pvFd7bfVq2GabTUGxbFn1\nz3j6abjmmmTs5MmVt0mTsj9qy1syo8PBYJazdh29NW5ccvJgMzP5bthQGSqnn54024d64QW47DJ4\n5RV4+eVNt1deSU5I3HrryrAYGiDVAqXasvHjq9dZxC2ZbgwrB4NZzorwx2Ps2KSPsd12yfPB+6Fm\nzEi2GKrZsCHZ8igPi/LwGFy2alWy5TLSOKl6gCweZha1p5+GH/0o2WqZNAkmTtz8fuLE9vVgyhUx\nrFrlYDCzthg7NulVTJnS+rrWrds8UF5+GX7/++rN+hdfTCZuXLMm6bVUu1+zJpnPbLjgaPZ+3brW\nv952G24rpl4OBjPbTN4ncm65ZfUTFL/2teqXk91rryQYRhIBa9cmQTFceAy9f+kleOqpkd+zYkX1\nz7vjDpg+PTlEeautKu+rLWtlzPjxlf2dyq2Ygk2JYWadqQi7t9pNSv6IbrUVbL99+9ZbKlU/qmz/\n/eF730vCaO3aZP6uavflj196qfaYaq+tX18ZFM8919rX5GAws46R95ZMIyZMSK5oOBo2bEh2aQ2G\nxQknwF13Nb8+B4OZdYwibskUIazGjt3UYIdkq6EVDgYzsxYUMaxa5WAwM+sy5Vsx1foftSgi2lpQ\nViRFp9RqZlYUkoiIhs5Jz+B0DzMz62QOBjMzq+BgMDOzCg4GMzOr4GAwM7MKDgYzM6vgYDAzswoO\nBjMzq+BgMDOzCg4GMzOr4GAwM7MKDgYzM6vgYDAzswqZB4Ok2ZKWSFoqac4I494nab2kE7KuyczM\nhpdpMEgaA1wEHAXsA5wiaa9hxl0I3JRlPe3W39+fdwlVFbEu11Qf11S/ItZVxJqakfUWw0zg0YhY\nHhHrgQXAcVXGfRq4CliVcT1tVdQfgiLW5Zrq45rqV8S6ilhTM7IOhh2Bx8uer0yXbSTpHcDxEfEv\nQEMXkzAzs/YrQvP5m0B578HhYGaWo0wv7SnpYKAvImanz+cCERHzysYsG3wIbA+8CpwZEdcPWZev\n62lm1oRGL+2ZdTCMBR4BjgCeAu4CTomIgWHGXwr8V0RcnVlRZmY2onFZrjwiNkg6B7iZZLfVJREx\nIOms5OWYP/QtWdZjZma1ZbrFYGZmnacIzecRSbpE0jOSHsi7lkGSdpL0C0kPS3pQ0t8WoKYtJf1W\n0r1pTefnXdMgSWMkLZJ0fe3Ro0PS/0i6P/1+3ZV3PQCSpki6UtJA+rP1/pzrmZ5+fxal9y8V5Gf9\n7yQ9JOkBSd+XNL4ANZ2b/t7l+veg2t9LSW+RdLOkRyTdJGlKrfUUPhiAS0lOkCuSN4DPRMQ+wAeA\nT1U7cW80RcQ64PCI2B94L3C0pJl51lTmXGBx3kUM8SZQioj9I6Io36dvATdGxAxgP6BqL260RMTS\n9PtzAHAgyYEh1+RZU3p4+6eBAyLiPSS7w0/OuaZ9gDOAg0h+9z4iafecyqn293Iu8N8RsSfwC+C8\nWispfDBExK+BF/Kuo1xEPB0R96WPV5P8Au848ruyFxFr0odbkvzC5L6fUNJOwDHAd/OuZQhRoJ9/\nSZOBD0XEpQAR8UZEvJxzWeWOBH4fEY/XHJm9scAkSeOAicCTOdczA/htRKyLiA3AL4FcpvYZ5u/l\nccD30sffA46vtZ7C/GJ0Kkm7kfyX8Nt8K9m4y+Ze4GnglohYmHdNwDeAz1OAkBoigFskLZT013kX\nA7wT+IOkS9NdN/MlbZV3UWVOAn6QdxER8STwf4EVwBPAixHx3/lWxUPAh9JdNhNJ/hHaOeeayr0t\nIp6B5J9a4G213uBgaIGkrUmm8jg33XLIVUS8me5K2gl4v6S986xH0v8Cnkm3rkSxTl48NN1FcgzJ\nrsAP5lzPOOAA4OK0rjUkuwByJ2kL4FjgygLUsi3Jf8C7Au8Atpb00TxrioglwDzgFuBG4F5gQ541\n1VDznzQHQ5PSzdirgMsj4rq86ymX7oK4FZidcymHAsemJzH+ADhc0mU51wRARDyV3j9Lst887z7D\nSuDxiLg7fX4VSVAUwdHAPen3Km9HAssi4vl0t83VwCE510REXBoRB0VECXgRWJpzSeWekTQNQNLb\nqWNOuk4JhqL9twnw78DiiPhW3oUASNp+8GiDdBfEh4EledYUEV+MiF0iYneSBuEvIuLUPGsCkDQx\n3dpD0iTgT0h2B+Qm3dR/XNL0dNERFKdhfwoF2I2UWgEcLGmCJJF8n3Jt0gNIemt6vwvwZ8AVeZZD\n5d/L64HT0sd/BdT8RzbTE9zaQdIVQAmYKmkFcP5ggy7Hmg4FPgY8mO7TD+CLEfGzHMvaAfheOoX5\nGOCHEXFjjvUU2TTgmnSalXHA9yPi5pxrAvhb4PvprptlwCdyrod0n/mRwJl51wIQEXdJuopkd836\n9H7oibJ5+LGk7Uhq+mReBw5U+3tJckmDKyWdDiwHTqy5Hp/gZmZm5TplV5KZmY0SB4OZmVVwMJiZ\nWQUHg5mZVXAwmJlZBQeDmZlVcDCYDZHOV5TLJGhmReBgMDOzCg4G6xmSdpW0OJ259CFJP5O05TDD\nD5P0G0m/K996kPTV9GIs90s6MV12mKT/Khvz/ySdmj6+MP2s+yR9JV22vaSr0gsr/VbSB8rWM3hh\nnHvS6TrMRl3hp8Qwa7M9gJMi4kxJPwT+nOrz2rw9Ig6VNINkrpmrJf058J6I2FfS24CFkm5Lx282\nhUA6RcLxEbFX+nxy+tK3gK9HxO2SdgZuAvYGPksyncId6VQUa9v2VZs1wMFgveaxiHgwfXwPsNsw\n464FiIiBNAQgmS32B+nyVZL6gfcBrwyzjpeA1yR9F7gB+Em6/EhgRjoJHCRTR08EfgN8Q9L3gasj\n4okmvj6zlnlXkvWadWWPNzD8P0fl44ab2Xdw+RskVxUbNAEgnRZ6JskU2h8Bflb2vvenl83cP52B\ndk1EzCO5RORWwG/KZlo1G1UOBus1zUzfPvieXwEnpVfKeyvwIeAukhkrZ0jaIr2QzBGwcWbSbdNZ\ndz8DvCddz80k18EmHbdfer97RDwcEV8BFgK5Xkfcepd3JVmvqWc64aFjAiAirpF0MHA/8Cbw+YhY\nBSDpRyTXdHgMWJS+bzJwnaQJ6fO/S+/PBS6WdD/JlsYvgU8C/0fS4SRbMg8DP238yzNrnafdNjOz\nCt6VZGZmFRwMZmZWwcFgZmYVHAxmZlbBwWBmZhUcDGZmVsHBYGZmFRwMZmZW4f8DwIc9pERldEwA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c5ac940>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_density(range(1, 11))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And at a wider range:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.43240754751464183"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEPCAYAAAAEfBBiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+UVeV97/H3Z7AgQRFdJP4i/iAEf6RNlAhJamqGGg32\nR+CSrqj39uYmthdaozdts7xqflwGb26LdtWYRk2Csdb0JkGXtYr1tzeOtlaBBFQM4OAviihqEwWV\nAKN87x/7OePmcGbmzHD22WdmPq+1Zrn3s5+993dvdb7zPPvZz1ZEYGZm1mxtZQdgZmYjkxOQmZmV\nwgnIzMxK4QRkZmalcAIyM7NSOAGZmVkpCk9AkmZJWiepS9KFfdSbLqlb0ty0PkbSMkmrJK2WtCBX\nd4Gk5yWtTD+zUvmRkrblyq8u+vrMzGxw9iny4JLagCuBU4EXgBWSbo2IdTXqLQLurpRFxA5JMyNi\nm6RRwEOS7oyI5anK5RFxeY3TPhUR0wq5IDMza5iiW0AzgPURsSEiuoElwOwa9c4HbgJezhdGxLa0\nOIYsWebfmlUv5+yt3MzMWkjRCehwYGNu/flU1kPSYcCciPgOVclDUpukVcBm4N6IWJHbfJ6kRyV9\nX9KEXPlRqfvtfkkfb+jVmJlZw7TCIIQrgPyzoZ4kFBG7IuJEYBLwEUnHp01XA5Mj4gSy5PQ3qfxF\n4IjUBfdl4EeS9iv6AszMbOAKfQYEbAKOyK1PSmV5JwFLJAmYCJwhqTsillYqRMRWSfcDs4A1EfFK\nbv9rgNtSvZ3AzrS8UtLTwFRgZf6EkjwBnpnZIEREwx5zFN0CWgFMSaPTRgNnAUvzFSJicvo5muw5\n0LkRsVTSREkHAEgaC5wGrEvrh+QOMRd4IpVPTAMakDQZmAI8UyuwiGj5nwULFpQeg+N0nEM5zqEQ\n41CKs9EKbQFFxNuSzgPuIUt210bEWknzs82xuHqX3PKhwPUpobQBN0TEHWnbZZJOAHYBzwHzU/kp\nwCWSdqZt8yPitSKuzczM9k7RXXBExF3AMVVl3+ul7jm55dVAzeHUEfG5XspvBm4edLBmZtY0rTAI\nwXrR3t5edgh1cZyN5TgbZyjECEMnzkZTEf16rU5SjMTrNjPbG5KIITQIwczMrCYnIDMzK4UTkJmZ\nlaLwUXBDybx5i+jq2r5H+dSp+7J48UUlRGRmNnw5AeV0dW3ngQc6amypVWZmZnvDXXBmZlYKJyAz\nMyuFE5CZmZXCCcjMzErhQQg5U6fuC3TQ1QXjxsHhh+fLzcyskTwVTw1f+UqWgL761SYGZWbW4jwV\nTxMceCC8+mrZUZiZDW9OQDU4AZmZFc8JqAYnIDOz4jkB1eAEZGZWPCegGpyAzMyK5wRUgxOQmVnx\nnIBqcAIyMyueE1AN++8P27bBW2+VHYmZ2fDlBFRDWxsccAC89lrZkZiZDV9OQL2YMMHdcGZmRSo8\nAUmaJWmdpC5JF/ZRb7qkbklz0/oYScskrZK0WtKCXN0Fkp6XtDL9zMptu1jSeklrJZ0+2Lj9HMjM\nrFiFTkYqqQ24EjgVeAFYIenWiFhXo94i4O5KWUTskDQzIrZJGgU8JOnOiFieqlweEZdXHec44LPA\nccAk4D5J7+9z4rdeOAGZmRWr6BbQDGB9RGyIiG5gCTC7Rr3zgZuAl/OFEbEtLY4hS5b5RFJrQrzZ\nwJKIeCsingPWpxgGzAnIzKxYRSegw4GNufXnU1kPSYcBcyLiO1QlFUltklYBm4F7I2JFbvN5kh6V\n9H1JB/Ryvk3V56uXE5CZWbFa4XtAVwD5Z0M9SSgidgEnShoP3CLp+IhYA1wNXBIRIekbwN8AfzyQ\nk3Z0dPQst7e3097evtt2JyAzG+k6Ozvp7Ows7PhFJ6BNwBG59UmpLO8kYIkkAROBMyR1R8TSSoWI\n2CrpfmAWsCYiXsntfw1wW+587+3nfMDuCaiWAw+E//iPPquYmQ1r1X+cL1y4sKHHL7oLbgUwRdKR\nkkYDZwFL8xUiYnL6OZrsOdC5EbFU0sRK15qkscBpwLq0fkjuEHOBJ9LyUuAsSaMlHQ1MAZYzCG4B\nmZkVq9AWUES8Lek84B6yZHdtRKyVND/bHIurd8ktHwpcn0bItQE3RMQdadtlkk4AdgHPAfPT+dZI\nuhFYA3STJbNBffLVCcjMrFj+JHcv7rsP/vIv4Sc/aVJQZmYtzp/kbpIDD/RUPGZmRXIC6oW74MzM\niuUE1AsnIDOzYvkZUC927YLRo2HHDhg1qkmBmZm1MD8DapK2tuy7QFu2lB2Jmdnw5ATUB3fDmZkV\nxwmoD05AZmbFcQLqgxOQmVlxnID64ARkZlYcJ6A+OAGZmRXHCagPTkBmZsVxAuqDE5CZWXGcgPrg\nBGRmVhwnoD44AZmZFccJqA9OQGZmxXEC6oMTkJlZcZyA+uAEZGZWHCegPjgBmZkVx59j6MPbb8OY\nMbBzZzY7tpnZSObPMTTRqFGw337+JIOZWRGcgPrhbjgzs2I4AfXDCcjMrBhOQP2YMMEJyMysCIUn\nIEmzJK2T1CXpwj7qTZfULWluWh8jaZmkVZJWS1pQY58vS9ol6aC0fqSkbZJWpp+r9zZ+t4DMzIqx\nT5EHl9QGXAmcCrwArJB0a0Ssq1FvEXB3pSwidkiaGRHbJI0CHpJ0Z0QsT/tMAk4DNlSd9qmImNao\na3ACMjMrRtEtoBnA+ojYEBHdwBJgdo165wM3AS/nCyNiW1ocQ5Ys82OnvwlcUONYDRsiCE5AZmZF\nKToBHQ5szK0/n8p6SDoMmBMR36EqeUhqk7QK2AzcGxErUvmngY0RsbrGOY9K3W/3S/r43l6AE5CZ\nWTEK7YKr0xVA/tlQTxKKiF3AiZLGA7dIOh54FvgKWfdb9T4vAEdExKuSplX2iYg3qk/a0dHRs9ze\n3k57e3vN4A48EDZUd/KZmY0AnZ2ddHZ2Fnb8QmdCkPRRoCMiZqX1i4CIiEtzdZ6pLAITgTeBeRGx\ntOpYX0/b7gHuA7alfSYBm4AZEfFy1T73A1+OiJVV5XXNhACwZAncfDPceGN912xmNlw1eiaEoltA\nK4Apko4EXgTOAs7OV4iIyZVlSdcBt0XEUkkTge6I2CJpLFmLZ1FEPAEcktvnWWBaavVMBH4ZEbsk\nTQamAM+wF9wFZ2ZWjEITUES8Lek8slZLG3BtRKyVND/bHIurd8ktHwpcn0bItQE3RMQdtU7DO11w\npwCXSNoJ7ALmR8Rre3MNTkBmZsXwZKT9WL8eZs2Cp58uOCgzsxbnyUibzC0gM7NiuAXUj7feyj7J\n0N3tTzKY2cjmFlCT7bMPjBsHW7eWHYmZ2fDiBFQHd8OZmTWeu+D6MG/eIrq6tvPTn8Ixx8D++2fl\nU6fuy+LFFxUcpZlZaxlq7wENaV1d23nggQ4AVu72KmtHCdGYmQ0v7oIzM7NSOAGZmVkpnIDMzKwU\nTkBmZlYKD0Low9Sp+wIdbN0KXV1w0kn5cjMz2xsehl2HX/wCpkzxu0BmNrJ5JoQSHHQQ7NrlBGRm\n1khOQHWQYPJkeGavvixkZmZ5TkB1cgIyM2ssJ6A6OQGZmTWWE1CdnIDMzBrLCahORx/tBGRm1khO\nQHVyC8jMrLH8HlCdduyA8ePhzTezj9SZmY00pb4HJGmcpFGNOvlQMmYMHHwwbNxYdiRmZsNDnwlI\nUpuk/yzpdkkvA+uAFyWtkfTXkqY0J8zW4G44M7PG6a8FdD/wPuBi4JCIeG9EvAf4OPAIcKmkPyw4\nxpbhBGRm1jj9JaBPRsT/jojHI2JXpTAifhkR/xgRnwFu6OsAkmZJWiepS9KFfdSbLqlb0ty0PkbS\nMkmrJK2WtKDGPl+WtEvSQbmyiyWtl7RW0un9XN+AOAGZmTVOn4/TI6K7spye/Ryc3yci/j1fp5qk\nNuBK4FTgBWCFpFsjYl2NeouAu3PH3iFpZkRsS+d+SNKdEbE87TMJOA3YkDvOccBngeOAScB9kt4/\n4BEHvZg8GW69tRFHMjOzugYhSDofeAm4F7g9/fxzHbvOANZHxIaUqJYAs2vUOx+4CXg5XxgR29Li\nGLLEl08k3wQuqDrObGBJRLwVEc8B61MMDeEWkJlZ49Q7oPhLwDER8YsBHv9wID9u7HmqEoKkw4A5\nETFTUvW2NuBnZM+hroqIFan808DGiFgt7TYi8HDg4dz6plTWEJMnw7PPNupoZmYjW70JaCOwpaAY\nrgDyz4Z6Mkp67nSipPHALZKOB54FvkLW/TZoHR0dPcvt7e20t7f3u8+73w3bt8OWLXDAAXtzdjOz\n1tfZ2UlnZ2dhx6/rRVRJ1wLHkHW97aiUR8Tl/ez3UaAjImal9Yuy3eLSXJ1Kp5aAicCbwLyIWFp1\nrK+nbfcA9wHb0j6TyFo6M4BzUlyL0j53AQsiYlnVsQb9WOiDH4Qf/ABOOGFQu5uZDVllvYj672TP\nf0YD++d++rMCmCLpSEmjgbOA3RJLRExOP0eTPQc6NyKWSpoo6QAASWPJWjzrIuKJiDgkt8/zwIkR\n8XI69pmSRks6GpgCLK/zGuvi50BmZo1RVxdcRCwczMEj4m1J55G1WtqAayNiraT52eZYXL1LbvlQ\n4Pr0HKgNuCEi7qh1GlK3XUSskXQjsAboJktmDZ1ryJOSmpk1Rp9dcJKuAf42IlbX2DYOOBPYERE/\nLC7ExtubLrhvfxvWroWrr25wUGZmLa7RXXD9tYCuAr4u6TeAJ4BXgH2B9wPjgb8DhlTy2VuTJ8Pt\nt5cdhZnZ0Nffi6iPAp+VtB9wElm32K+AtRHxZBPiazl+BmRm1hj1joL7feD2/HQ8Q9nedMH96lcw\nYQJs2wajRuS84GY2UpU1Cu5MYL2kyyQd26iTD0Vjx8LEibBpU9mRmJkNbXV/kC69DHo28AWykWfX\nAT+OiNeLC68Yg20BzZu3iK6u7axalY2GmzAhK586dV8WL76owVGambWWZg9C6BERWyXdBIwF/gz4\nT8AFkv42Ir7dqIBaWVfXdh54oAOAxx7Lb+koIRozs6Gt3slIZ0v6J6AT+DVgRkScAXwI+HJx4ZmZ\n2XBVbwtoLvDNiHgwX5g+lfBHjQ/LzMyGu3oHIWyuTj6SLgWIiP/X8KjMzGzYqzcB1Zp5+oxGBmJm\nZiNLn11wkv4UOBd4n6THc5v2Bx4qMrBWNHXqvlQGHKxZAwcdBIccUik3M7OB6G8uuAOAA4G/AvLj\njF+PiF8WHFth9uZF1IorroB16+C7321QUGZmLa7ZL6JG+rT1F4HXcz9IOqhRQQxFH/sYPPJI2VGY\nmQ1d/bWA/jkifk/Ss+Q+e5BEREwuOsAiNKIFtGNH1gW3eTPsX8+XkczMhrimvogaEb+X/nl0o044\nXIwZAx/6EPz0pzBzZtnRmJkNPfW+iHpy+v4Pkv5Q0uWSjig2tNb3sY/Bww+XHYWZ2dBU7zDs7wDb\nJFVmPnga+IfCohoiPvpRJyAzs8GqNwG9lR6azAaujIiryIZij2iVgQiN/ei3mdnIUG8Cel3SxcAf\nArdLaiObE25EmzQpexb09NNlR2JmNvQM5HtAO4A/iojNwCTgrwuLagjxcGwzs8Gp+3tAw0kjhmFX\nXH551gK66qqGHM7MrGWV8kVUSXMlrZe0RdJWSa9L2tqoIIYyD0QwMxuculpAkp4Cfj8i1hYfUvEa\n2QLavj17IfWVV2DcuIYc0sysJZX1RdSXBpt8JM0CriBrbV0bEZf2Um868G/AmRFxs6QxwIPA6BTn\nTRGxMNW9hGxE3i7gJeDzEbFZ0pHAWmBdOuwjEXHuYOKuR+UT3aNGwckn+xPdZmYDUW8C+qmkG4Bb\nyAYjABARN/e1UxotdyVwKvACsELSrRGxrka9RcDduWPvkDQzffRuFPCQpDsjYjlwWUT8r7Tv+cAC\n4E/Trk9FxLQ6r2uv+BPdZmaDV28CGg9sA07PlQXQZwICZgDrI2IDgKQlZC2XdVX1zgduAqbnCyNi\nW1ock2KNVP5Grto4spZQRcOah2ZmVpy6ElBEfGGQxz8c2Jhbf54sKfWQdBgwJyJmSqre1gb8DHgf\ncFVErMht+wbwOeA1ID8b21GSVgJbgK9HxL8OMnYzMytQXQlI0lSy6XgOjohfl/RB4NMR8Y0GxHAF\ncGH+dJWFiNgFnChpPHCLpOMjYk3a9jXga5IuJGtBdQAvAkdExKuSpuX2ybeYAOjo6OhZbm9vp729\nvQGXYmY2fHR2dtLZ2VnY8evtgrsGuAD4HkBEPC7pR0B/CWgTkJ+0dFIqyzsJWCJJwETgDEndEbG0\nUiEitkq6H5gFrKna/0fAHUBHROwEdqZ9Vkp6GpgKrKwOLJ+AzMxsT9V/nC9cuLChx683Ab0rIpZn\nOaLHW3XstwKYkkanvQicBZydr5D/ppCk64DbImKppIlAd0RskTQWOI1soAKSpkTEU2m3OWQj30j7\n/DIidkmaDEwBnqnzGgcs/4nunTth2TL4zd/0J7rNzOpRbwL6D0nvIw0CkPQHZAmlTxHxtqTzgHt4\nZxj2Wknzs82xuHqX3PKhwPXpOVAbcENE3JG2LUrdgruADcCfpPJTgEsk7Uzb5kfEa3Ve44BVD7U+\n9VT44hdh7tyizmhmNnzU+yLqZGAx8JvAq8CzwH+pjG4bahr5Imred78LDzwAP/5xww9tZla6Rr+I\n2t8nuf+iqmgsWWvkTYCIuLxRgTRTUQno5Zdh6lR48UUYO7bhhzczK1Wz54LbP/2cRPai54HABLIu\nr6a87DmUvOc9MG0a3HVX2ZGYmbW+ervgHgR+NyJeT+v7A7dHxCkFx1eIolpAAFdfDQ89BD/8YSGH\nNzMrTVlzwR1MGt6c7ExlljNv3iKeeGI7y5fDxo3QltqXnhvOzGxP9SagHwDLJf1TWp8D/H0hEQ1h\nXV3befjhDgD+5V/yWzpKiMbMrLXVOxXP/5F0J/BbqegLEbGquLDMzGy4q7cFRESspMaMAmZmZoNR\n1xdRzczMGs0JyMzMSlF3F5z1Lz83HMDjj8O73+254czMaqnrPaDhpsj3gPI6O2H+fFizBkaNKvx0\nZmaFavZMCLYXPvEJmDABbr217EjMzFqPW0AF++QnF7Fs2XamTYP81yz8cqqZDTVlzYRgg9TdvZ03\n3ujgwQert3SUEI2ZWetwF1zB1LC/FczMhhcnIDMzK4UTkJmZlcIJyMzMSuFBCAXLv5y6bRusWgUf\n/rBfTjUz8zDsJvvwhxexYcN2fv3Xdy/3sGwza3Uehj3EjRu3nV/8ooMHHqje0lFCNGZm5fEzoCZr\n8x03MwOcgMzMrCSFJyBJsyStk9Ql6cI+6k2X1C1pblofI2mZpFWSVktakKt7iaTH0ra7JB2S23ax\npPWS1ko6vdirMzOzwSo0AUlqA64EPgV8ADhb0rG91FsE3F0pi4gdwMyIOBE4AThD0oy0+bKI+FDa\ndjuwIB3neOCzwHHAGcDVkuciMDNrRUUPQpgBrI+IDQCSlgCzgXVV9c4HbgKm5wsjYltaHEMWa6Ty\nN3LVxgG70vKngSUR8RbwnKT1KYZljbqgvbXnN4N+zquvjmPNmjba2zt2q+dRcWY2nBWdgA4HNubW\nnydLCD0kHQbMiYiZuRZOZVsb8DPgfcBVEbEit+0bwOeA14CZufM9nDvEplTWMqqTSnt7Bw880MEr\nr1A1Mq6jmWGZmTVdKwzDvgLIPxvq6TKLiF3AiZLGA7dIOj4i1qRtXwO+lp4rnc8Af2N3dLxTvb29\nnfb29kGGb2Y2PHV2dtLZ2VnY8YtOQJuAI3Lrk1JZ3knAkvSsZiLZs57uiFhaqRARWyXdD8wC1lTt\n/yOy50Ad6djv7ed8wO4JyMzM9lT9x/nChQsbevyiE9AKYIqkI4EXgbOAs/MVImJyZVnSdcBtEbFU\n0kSgOyK2SBoLnEY2UAFJUyLiqbTbHN55prQU+KGkb5J1vU0Blhd2dYVYBGxn1arn/EzIzIa1QhNQ\nRLwt6TzgHrIRd9dGxFpJ87PNsbh6l9zyocD16TlQG3BDRNyRti2SNJVs8MEG4E/S+dZIupGsldQN\nnFvanDuDth3oYOtWPxMys+Gt8GdAEXEXcExV2fd6qXtObnk1MK2Xen/Qx/n+CvirQQVbgupRcY8+\n+hxbtpQWjplZ03gy0hZTGRX3jqxL7oADnuOEE47qKXWXnJk1mycjHXGyLrktW9wlZ2bDi+eCMzOz\nUrgF1GLqeya0iEcf3X2UXGVfd8uZ2VDhBNRias+UUF1rO1u2/L2/KWRmQ5oT0LCQDVSobhW5RWRm\nrcwJqMVVd8lBrW45D1Qws6HHCajF1WrB1O6WMzMbWpyAhiV3yZlZ63MCGoL6HylX3SVXSUjr6Orq\n2O04TkhmVhYnoCGovpFyeX5GZGatxwloGBj4fHLuojOz8nkuuGFoz/nkOti9tVO97vnmzKx/ngvO\nCuAuOjNrPiegYagxn3jwdD9mVix3wY0A8+Ytoqtre896lpD+Plejgz1bO9Vl7qYzG+ncBWcDNvBR\nc7V4aLeZNZYT0AjUmC66vhPSk0/+nF/9ahxjx7ZxzDFH7HZuJygzA3fBGXt20UE93XQDXQdYxAEH\nrNutCw+clMyGCnfBWcM1b7656s9IuNVkNpI5AVlNjemm6091N15Hz/rmzeDnTGbDmxOQ1VT9Cz7r\npuvoWS8mIVXrf+DDk0/+HNh/txYTOEmZDQVOQFaX1kxIUGk1ZS0m6K9bz918Zq2j8AQkaRZwBdAG\nXBsRl/ZSbzrwb8CZEXGzpDHAg8DoFOdNEbEw1b0M+H1gB/A08IWI2CrpSGAtsC4d9pGIOLe4qxu5\n+ktI2S/6z+/2i76cJNXB7t161ev9P4dyK8usGIUmIEltwJXAqcALwApJt0bEuhr1FgF3V8oiYoek\nmRGxTdIo4CFJd0bEcuAe4KKI2CVpEXBx+gF4KiKmFXldtqd6fhGX02rqT38J650yt7LMGqvoFtAM\nYH1EbACQtASYzTstlIrzgZuA6fnCiNiWFseQxRqp/L5ctUeAz+TWGzZE0BproK2m1khQtexdKwv2\nbFU5idlIVHQCOhzYmFt/niwp9ZB0GDAnImZKqt7WBvwMeB9wVUSsqHGOc4AlufWjJK0EtgBfj4h/\n3fvLsCL094u0NVtMg1HPs6sO9qarsJ6uQyc1azWtMAjhCuDC3HpPCyYidgEnShoP3CLp+IhY01NR\n+irQHRE/SkUvAEdExKuSpuX2eaP6pB0dHT3L7e3ttLe3N/CSrBH6azFB5ZfsOU1+zlSGgba63inr\nPck1Pqm5JTe8dHZ20tnZWdjxi05Am4D8k9tJqSzvJGCJJAETgTMkdUfE0kqFNMDgfmAWsAZA0ueB\n3wF+O1evG3g1La+U9DQwFVhZHVg+AdnQMJjnTNXdekOnm68Zikhq/a3vfXfkYJKeW4ODU/3H+cKF\nCxt6/KIT0ApgShqd9iJwFnB2vkJETK4sS7oOuC0ilkqaSNa62SJpLHAa2X+9lZF1FwCnRMSO3P4T\ngV+mwQmTgSnAM4VeobWUgf7yqGf03shpZTVDI7ojB7pezzma3xos4hwVQyWRFpqAIuJtSeeRjVqr\nDMNeK2l+tjkWV++SWz4UuD49B2oDboiIO9K2b5MNz743azj1DLc+BbhE0k5gFzA/Il4r6vps6BvM\n/6RuZQ1HZbQGizjHwF8rGEhSa7TCnwFFxF3AMVVl3+ul7jm55dVAzeHUEfH+XspvBm4edLBmddjb\nVhbs2apyErPGGMxrBdV1au1TKW+sVhiEYDasNaIrZKCtrkqdvpKck5qVzQnIbAgooj+/iKTmlpwN\nhBOQ2QhVxkPqRnRHDnS9nnM4MZbDCcjMmqZVR2aV0Ros4hxDLZH6i6hmZsNE9deNGz0K7pprLm7o\nF1GdgMzMrC6N/iR3W6MOZGZmNhBOQGZmVgonIDMzK4UTkJmZlcIJyMzMSuEEZGZmpXACMjOzUjgB\nmZlZKZyAzMysFE5AZmZWCicgMzMrhROQmZmVwgnIzMxK4QRkZmalcAIyM7NSOAGZmVkpCk9AkmZJ\nWiepS9KFfdSbLqlb0ty0PkbSMkmrJK2WtCBX9zJJayU9KukfJY3PbbtY0vq0/fRir87MzAar0AQk\nqQ24EvgU8AHgbEnH9lJvEXB3pSwidgAzI+JE4ATgDEkz0uZ7gA9ExAnAeuDidJzjgc8CxwFnAFdL\natjX+5qts7Oz7BDq4jgby3E2zlCIEYZOnI1WdAtoBrA+IjZERDewBJhdo975wE3Ay/nCiNiWFscA\n+wCRyu+LiF1p2yPApLT8aWBJRLwVEc+RJacZDFFD5T9Kx9lYjrNxhkKMMHTibLSiE9DhwMbc+vOp\nrIekw4A5EfEdQFXb2iStAjYD90bEihrnOAe4o5fzbao+n5mZtYZWGIRwBZB/NtSThCJiV+qCmwR8\nJHWxvVNR+irQHRE/bkqkZmbWMIqI4g4ufRToiIhZaf0iICLi0lydZyqLwETgTWBeRCytOtbXgTcj\n4vK0/nngvwO/nZ4X7XF8SXcBCyJiWdWxirtoM7NhLCIa9ly96AQ0CngSOBV4EVgOnB0Ra3upfx1w\nW0TcLGkiWetmi6SxZAMUFkXEHZJmAX8DnBIRv8jtfzzwQ+AjZF1v9wLvjyIv0szMBmWfIg8eEW9L\nOo9s1FobcG1ErJU0P9sci6t3yS0fClyfRsi1ATdEROVZz7eB0cC9aZDbIxFxbkSskXQjsAboBs51\n8jEza02FtoDMzMx60wqDEJqq3hdjmxTLc5IeSy/bLk9lB0q6R9KTku6WdECuflNespV0raSXJD2e\nKxtwXJKmSXo83esrmhTnAknPS1qZfma1QJyTJP1E0s/TS9X/I5W31D2tEef5qbxl7ql6eUG9Be9l\nb3G2zL2sirctxbM0rTfnfkbEiPkhS7hPAUcCvwY8ChxbYjzPAAdWlV0K/M+0fCHZcy+A44FVZN2m\nR6XrUEFxfZzs5d/H9yYuYBkwPS3fAXyqCXEuAP6iRt3jSozzEOCEtLwf2XPRY1vtnvYRZ0vdU+Bd\n6Z+jyN4DnNFq97KPOFvqXubO/+fA/wWWpvWm3M+R1gKq98XYZhF7tkJnA9en5euBOWm5aS/ZRsS/\nAq/uTVzSz699AAAE7ElEQVSSDgH2j3fe3fpBbp8i44Sq98mS2SXGuTkiHk3LbwBryV4taKl72kuc\nlffoWuaeRu0X1FvqXvYRJ7TQvYSs5Qv8DvD9qngKv58jLQH1+2JskwXZQIoVkv44lR0cES9B9gsB\neE8qL/sl2/cMMK7Dye5vRTPv9XnK5gn8fq7roCXilHQUWavtEQb+77ppsebirLzC0DL3VLVfUG+5\ne9lLnNBC9zL5JnABuw8Ca8r9HGkJqNWcHBHTyP76+KKk32L3/wiosd4qWjWuq4HJkc0TuJlsuH5L\nkLQf2ZRTX0otjJb8d10jzpa6p7H7C+ozJH2AFryXNeI8nha7l5J+F3gptXz7er+nkPs50hLQJuCI\n3PqkVFaKiHgx/fMV4BayLrWXJB0MkJq1lfnxNgHvze3e7NgHGlcp8UbEK5E6oYFreKebstQ4Je1D\n9kv9HyLi1lTccve0Vpytek8jYivQCcyiBe9lrThb8F6eDHxa2YQAPwZ+W9I/AJubcT9HWgJaAUyR\ndKSk0cBZwNJ+9imEpHelvzSRNA44HVid4vl8qvbfgMovq6XAWZJGSzoamEL2Ym9hIbL7X0QDiis1\n27dImiFJwOdy+xQWZ/qfpWIu8ESLxPl3wJqI+FaurBXv6R5xttI9lTSx0m2l7AX108ieVbXUvewl\nznWtdC8BIuIrEXFEREwm+334k4j4r8BtNON+Nno0Rav/kP219CTZw7OLSozjaLJReKvIEs9Fqfwg\n4L4U4z3AhNw+F5ONOlkLnF5gbD8CXgB2AP8OfAE4cKBxAR9O17Ye+FaT4vwB8Hi6t7eQ9WWXHefJ\nwNu5f98r03+HA/53XWSsfcTZMvcU+I0U16Mppq8O9v+bgu9lb3G2zL2sEfMneGcUXFPup19ENTOz\nUoy0LjgzM2sRTkBmZlYKJyAzMyuFE5CZmZXCCcjMzErhBGRmZqVwAjIrkKTrJM0tOw6zVuQEZGZm\npXACMhuENJ3TGkmLJT0h6S5JY3qp/glJD0l6Kt8akvTXyj5W9pikz6ayT0i6LVfn25I+l5YXpXM9\nKumyVDZR0k3KPn62TNLHcsdZpewjYz9L0z2ZtZR9yg7AbAibApwZEfMk3QB8hmx6oGqHRMTJko4j\nm0vrZkmfAT4YEb8h6T3ACkkPpPp7TE8i6SBgTkQcm9bHp03fAi6PiH+T9F7gbrKPhn0ZODciHpb0\nLmB7w67arEGcgMwG79mIWJ2Wf0b2hchabgGIiLUp2UA279qPU/nLkjqB6cDrvRxjC/ArSd8Hbgf+\nOZV/EjguTQAJsF9KOA8B35T0Q+DmiCht1nez3rgLzmzwduSW36b3P+jy9Xr75kql/C2yTzhX7AsQ\nEW+TTd1/E/B7wF25/T4SESemnyMiYltEXAr8ETAWeEjS1DqvyaxpnIDMBq+vD3j1t8+/AGemr2a+\nG/gtss9rbCBr0fyapAnAqZB9voNsRuK7gL8APpiOcw/wpZ6DSx9K/5wcET+PiMvIPkNy7CBiNSuU\nu+DMBq+eqeRrfqkzIv5J0keBx4BdwAUR8TKApBvJvhPzLNmU/gDjgVsl7ZvW/zz980vAVZIeI2s5\nPQicC/yZpJlkLbOfA3cO/PLMiuXPMZiZWSncBWdmZqVwAjIzs1I4AZmZWSmcgMzMrBROQGZmVgon\nIDMzK4UTkJmZlcIJyMzMSvH/AUDe4DLY0Sj0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ca875c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_density(range(100, 4000, 50))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The density is going down, and the curve is almost but not quite flat. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "# lim <sub>n &rarr; &infin;</sub> density(n)\n",
    "\n",
    "[Thes puzzle](http://fivethirtyeight.com/features/can-you-solve-the-puzzle-of-your-misanthropic-neighbors/)\n",
    "is to figure out the limit of `density(n)` as `n` goes to infinity. The plot above makes it look like  0.432+something, but we can't answer the question just by plotting;  we'll need to switch modes from *computational* thinking to *mathematical* thinking.\n",
    "\n",
    "At this point I started playing around with `density(n)`, looking at various values, differences of values, ratios of values, and ratios of differences of values, hoping to achieve some mathematical insight. Mostly I got dead ends.\n",
    "But eventually I hit on something promising. I looked at the difference between density values (using the function `diff`), and particularly the difference as you double `n`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0014849853757253895"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def diff(n, m): return density(n) - density(m)\n",
    "\n",
    "diff(100, 200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And compared that to the difference when you double `n` again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0007424926878626947"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(200, 400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmm&mdash;I noticed that the first difference is just about twice as much as the second. Let's check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(100, 200) / diff(200, 400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow&mdash;not only is it *close* to twice as much, it is *exactly* twice as much (to the precision of floating point numbers). Let's try other starting values for `n`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 500; diff(n, 2*n) / diff(2*n, 4*n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9999999999992524"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 1000; diff(n, 2*n) / diff(2*n, 4*n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, I'm convinced this is real! \n",
    "\n",
    "Now, what mathematical function behaves like this? I figured out that *f*(*n*) = (1 / *n*) does. The ratio of the differences would be:\n",
    "\n",
    "$$\\frac{f(n) - f(2n)}{f(2n) - f(4n)} = \\frac{1/n - 1/(2n)}{1/(2n) - 1 / (4n)}\\;\\;$$\n",
    "\n",
    "Multiplying top and bottom by *n* you get:\n",
    "\n",
    "\n",
    "$$\\frac{1 - 1/2}{1/2 - 1 /4} = \\frac{1/2}{1/4} = 2\\;\\;$$\n",
    "\n",
    "If the function (1 / *n*) fits the pattern, then so does any *affine transformation* of (1 / *n*), because we are taking the ratio of differences.  So that means a density function of the form\n",
    "\n",
    "    density(n) = A + B / n \n",
    "    \n",
    "would fit the patterm.  I can try a `curve_fit` to estimate the parameters *A* and *B*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  7.64810902e-34,  -2.42654280e-31],\n",
       "       [ -2.42654280e-31,   4.67778766e-28]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "Ns = list(range(100, 10001, 100))\n",
    "\n",
    "def f(n, A, B): return A + B / n\n",
    "\n",
    "((A, B), covariance) = curve_fit(f=f, xdata=Ns, ydata=[density(n) for n in Ns])\n",
    "\n",
    "covariance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `curve_fit` function returns a sequence of parameter values, and a covariance matrix. The fact that all the numbers in the covariance matrix are really small indicates that the parameters are a really good fit.  Here are the parameters, `A` and `B`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.43233235838169343, 0.29699707514520357)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A, B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plug them into a function that estimates the density:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def estimated_density(n): return A + B / n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can test how close this function is to the true `density` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.8857805861880479e-16"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max(abs(density(n) - estimated_density(n))\n",
    "    for n in range(200, 4000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That says that, for all values of `n` from 200 to 4,000, `density(n)` and `estimated_density(n)` agree at least through the first 15 decimal places!\n",
    "\n",
    "We now have a plausible answer to the puzzle:\n",
    "\n",
    "> lim<sub style=\"font-size:large\"> <tt>n</tt>&rarr;&infin;</sub> <tt>density(n) &cong; A = 0.43233235838169343</tt>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Why?\n",
    "\n",
    "Theis answer is empirically strong (15 decimal places of accuracy) but theoretically weak: we don't have a **proof**, and we don't have an explanation for **why** the density function has this form. We need some more mathematical thinking. \n",
    "\n",
    "I didn't have any ideas, so I looked to see if anyone else had written something about the number 0.43233235838169343. I tried several searches and found two interesting formulas:\n",
    "\n",
    "- Search: [`[0.4323323]`](https://www.google.com/search?q=0.4323323) Formula: `sinh(1) / exp(1)` [(Page)](http://arxiv.org/pdf/1310.4360.pdf)\n",
    "- Search: [`[0.432332358]`](https://www.google.com/search?q=0.432332358) Formula: `0.5(1-e^(-2))` [(Page)](http://www.actuarialoutpost.com/actuarial_discussion_forum/archive/index.php/t-52095.html)\n",
    "\n",
    "I can verify that the two formulas are equivalent, and that they are indeed equal to `A` to at least 15 decimal places:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.43233235838169365, 0.43233235838169365, 0.43233235838169343)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from math import sinh, exp, e\n",
    "\n",
    "S = sinh(1) / exp(1)\n",
    "E = 0.5 * (1 - e ** (-2))\n",
    "\n",
    "assert S == E\n",
    "\n",
    "S, E, A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So I now have a suspicion that\n",
    "\n",
    "> lim <sub>n &rarr; &infin;</sub> density(n) = (1 - e<sup>-2</sup>) / 2\n",
    "\n",
    "but I still have no proof, nor any intuition as to why this is so."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# John Lamping to the Rescue\n",
    "\n",
    "I reported my results to Anne Paulson and John Lamping, who had originally related the problem to me, and the next day John wrote back with the following:\n",
    "\n",
    "<blockquote>\n",
    "I got a derivation of the formula!\n",
    "\n",
    "<p>Suppose that each house has a different \"attractiveness\", and that when it is a misanthrope's turn to pick a house, they consider the houses in order, from most attractive to least, and pick the first house not next to any other house. If the attractiveness are chosen independently, this process gives the same probabilities as each misanthrope picking from the available houses randomly.\n",
    "\n",
    "<p>To be more concrete, let the attractivenesses range from 0 to 1 with uniform probability, with lower numbers being being considered earlier, and hence more attractive. (This makes the math come out easier.)\n",
    "\n",
    "<p>Given the attractiveness, a, of a house, we can compute the chance that it will end up getting picked. It will get picked if and only if neither house on either side gets considered earlier and gets picked. Lets consider one side, and assume the houses, and their attractivenesses, are labeled a, b, c, d, e, f, g, ... with the house we are considering being a. There are an infinite number of cases where b won't get considered before a and picked. Here are the first few:\n",
    "\n",
    "<p>   a < b  (a is considered before b)\n",
    "<p>   a > b > c < d  (Someone considered c first among these 4 and picked it. A later person considered b before a, but rejected it because c had already been chosen.)\n",
    "<p>   a > b > c > d > e < f (Someone considered e first among these 6, and picked it. A later person considered d, but rejected it because e was chosen, then considered c, which they picked. Still later, someone considered b, which they rejected because c had been chosen.)\n",
    "\n",
    "<p>We can write down the probabilities of these cases as a function of the attractiveness of a:\n",
    "\n",
    "<p>   a < b: The chance that b is greater than a: 1 - a.\n",
    "<p>   a > b > c < d: Let y = a - c, so y is between 0 and a. The probability is the integral over the possible y of the chance that b is between a and c, times the chance that d is greater than c, or integral from 0 to a of y * (1 - a + y) dy.\n",
    "<p>   a > b > c > d > e < f: Let y = a - e. The probability is the integral of (the chance that b, c, and d are all between a and e, and ordered right, which is y^3 / 3!, times the chance that f is greater than e, which is (1 - a + y))\n",
    "\n",
    "<p>If you work out the definite integrals, you get\n",
    "\n",
    "<p>   a < b: 1 - a\n",
    "<p>   a > b > c < d: a^2 / 2 - a^3 / 3!\n",
    "<p>   a > b > c > d > e < f: a^4 / 4! - a^5 / 5!\n",
    "\n",
    "<p>Add them all up, and you have 1 - a + a^2 / 2 - a^3 / 3! + a ^4 / 4! ... the Taylor expansion for e^-a. \n",
    "\n",
    "<p>Now, there will be a house at a if both adjacent houses are not picked earlier, so the chance is the square of the chance for one side: e^(-2a). Integrate that from 0 to 1, and you get 1/2 (1 - e^-2).\n",
    "</blockquote>\n",
    "\n",
    "You can compare John Lamping's solution to the solutions by [Jim Ferry](http://fivethirtyeight.com/features/the-perplexing-puzzle-of-the-proud-partygoers/) and [Andrew Mascioli](http://andrewmascioli.com/index.php/2016/04/24/the-misanthropic-neighborhood/). \n",
    "\n",
    "# Styles of Thinking\n",
    "\n",
    "It is clear that different write-ups of this problem display different styles of thinking. I'll attempt to name and describe them:\n",
    "\n",
    "- **Mathematical Publication Style**: This style uses sophisticated mathematics (e.g. generating functions, differentiation, asymptotic analysis, and manipulation of summations). It defines new abstractions without\n",
    "necessarily trying to motivate them first, it is terse and formal, and it gets us to the conclusion in a way that is\n",
    "clearly correct, but does not describe all the steps of how the author came up with the ideas. (Ferry and Mascioli)\n",
    "- **Mathematical Exploration Style**: Like the Mathematical Publication Style, but with more explanatory prose. (Lamping)\n",
    "- **Computational Thinking**: A more concrete style; tends to use programs with specific values of `n` rather than \n",
    "creating a proof for all values of `n`; produces tables or plots to help guide intuition; verifies results\n",
    "with tests or alternative implementations. (Norvig)\n",
    "\n",
    "In this specific puzzle, my steps were:\n",
    "- Analyze the problem and implement code to solve it for small values of `n`.\n",
    "- Plot results and examine them for insight.\n",
    "- Play with the results and get a guess at a partial solution <tt>density(n) &cong; A + B / n</tt>.\n",
    "- Solve numerically for `A` and `B`.\n",
    "- Do a search with the numeric value of `A` to find a formula with that value.\n",
    "- Given that formula, let Lamping figure out how it corresponds to the problem.\n",
    "\n",
    "This is mostly  computational thinking, with a little mathematical thrown in. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Validation by Anticipating Objections\n",
    "\n",
    "Is our implementation of `occ(n)` correct? I think it is, but I can anticipate some objections and answer them:\n",
    "\n",
    "- In `occ(n)`, is it ok to start from all empty houses, rather than considering layouts of partially-occupied houses? Yes, I think it is ok, because the problem states that initially all houses are empty, and each choice of a house breaks the street up into runs of acceptable houses, flanked by unacceptable houses. If we get the computation right for a run of `n` acceptable houses, then we can get the whole answer right. A key point is that the chosen first house breaks the row of houses into 2 runs of *acceptable* houses, not 2 runs of *unoccupied* houses. If it were unoccupied houses, then we would have to also keep track of whether there were occupied houses to the right and/or left of the runs. By considering runs of acceptable houses, eveything is clean and simple.\n",
    "\n",
    "- In `occ(7)`, if the first house chosen is 2, that breaks the street up into runs of 1 and 3 acceptable houses. There is only one way to occupy the 1 house, but there are several ways to occupy the 3 houses. Shouldn't the average give more weight to the 3 houses, since there are more possibilities there? No, I don't think so. We are caclulating occupancy, and there is a specific number (5/3) which is the expected occupancy of 3 houses; it doesn't matter if there is one combination or a million combinations that contribute to that expected value, all that matters is what the expected value is.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Validation by Simulation\n",
    "\n",
    "\n",
    "A simulation can add credence to our implementation of `density(n)`, for two reasons:\n",
    "- The code for the simulation can  have a more direct corresponance to the problem statement.\n",
    "- When two very different implementations get the same result, that is evidence supporting both of them.\n",
    "\n",
    "\n",
    "The simulation will start with an empty set of occupied houses. Following Lamping's suggestion, we sample <i>n</i> houses (to get a random ordering) and go through them, occupying just the ones that have no neighbor\" "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "def simulate(n):\n",
    "    \"Simulate moving in to houses, and return a sorted tuple of occupied houses.\"\n",
    "    occupied = set()\n",
    "    for house in random.sample(range(n), n):\n",
    "        if (house - 1) not in occupied and (house + 1) not in occupied:\n",
    "            occupied.add(house)\n",
    "    return sorted(occupied)\n",
    "\n",
    "def simulated_density(n, repeat=10000):\n",
    "    \"Estimate density by simulation, repeated `repeat` times.\"\n",
    "    return mean(len(simulate(n)) / n \n",
    "                for _ in range(repeat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see if the simulation returns results that match the actual `density` function and the `estimated_density` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  n   simul    density  estimated\n",
      " 25   0.444    0.444    0.444\n",
      " 50   0.438    0.438    0.438\n",
      "100   0.435    0.435    0.435\n",
      "200   0.434    0.434    0.434\n",
      "400   0.433    0.433    0.433\n"
     ]
    }
   ],
   "source": [
    "print('  n   simul    density  estimated')\n",
    "for n in (25, 50, 100, 200, 400):\n",
    "    print('{:3}   {:.3}    {:.3}    {:.3}'\n",
    "          .format(n, simulated_density(n), density(n), estimated_density(n)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We got perfect agreement (at least to 3 decimal places), suggesting that either our three implementations are all correct, or we've made  mistakes in all three. \n",
    "\n",
    "The `simulate` function can also give us insights when we look at the results it produces:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 2, 4, 6]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simulate(7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's repeat that multiple times, and store the results in a `Counter`, which tracks how many times it has seen each result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({(0, 2, 4, 6): 3303,\n",
       "         (0, 2, 5): 1255,\n",
       "         (0, 3, 5): 959,\n",
       "         (0, 3, 6): 783,\n",
       "         (1, 3, 5): 1426,\n",
       "         (1, 3, 6): 974,\n",
       "         (1, 4, 6): 1300})"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from collections import Counter\n",
    "\n",
    "Counter(tuple(simulate(7)) for _ in range(10000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That says that about 1/3 of the time, things work out so that the 4 even-numbered houses are occupied. But if anybody ever chooses an odd-numbered house, then we are destined to have 3 houses occupied (in one of 6 different ways, of which (1, 3 5) is the most common, probably because it is the only one that has three chances of getting started with an odd-numbered house)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Verification by Test\n",
    "\n",
    "Another way to gain more confidence in the code is to run a test suite:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ok'"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def test():\n",
    "    assert occ(0) == 0\n",
    "    assert occ(1) == occ(2) == 1\n",
    "    assert occ(3) == 5/3\n",
    "    assert density(3) == occ(3) / 3\n",
    "    assert density(100) == occ(100) / 100\n",
    "    assert runs(3) == [(0, 1), (0, 0), (1, 0)]\n",
    "    assert runs(7) == [(0, 5), (0, 4), (1, 3), (2, 2), (3, 1), (4, 0), (5, 0)]\n",
    "    for n in (3, 7, 10, 20, 100, 101, 200, 201):\n",
    "        for repeat in range(500):\n",
    "            assert_valid(simulate(n), n)         \n",
    "    return 'ok'\n",
    "\n",
    "def assert_valid(occupied, n):\n",
    "    \"\"\"Assert that, in this collection of occupied houses, no house is adjacent to an\n",
    "    occupied house, and every unoccupied position is adjacent to an occupied house.\"\"\"\n",
    "    occupied = set(occupied) # coerce to set\n",
    "    for house in range(n):\n",
    "        if house in occupied:\n",
    "            assert (house - 1) not in occupied and (house + 1) not in occupied\n",
    "        else:\n",
    "            assert (house - 1) in occupied or (house + 1) in occupied\n",
    "\n",
    "test()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion\n",
    "\n",
    "I'm happy with my result:\n",
    "\n",
    "> <tt>lim <sub>n &rarr; &infin;</sub> density(n) = (1 - e<sup>-2</sup>) / 2 &cong; 0.43233235838169365 </tt>\n",
    "\n",
    "I got the fun of working on the puzzle, the satisfaction of seeing John Lamping work out a proof, and the awe of seeing the mathematically sophisticated solutions of Jim Ferry and Andrew Mascioli, solutions that I know I could never come up with, but that I could get near to by coming at the problem with a different style of thinking."
   ]
  }
 ],
 "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
}