notes_03_random numbers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Numbers and Applications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Deterministic Randomness\n",
    "Some problems are physically unpredictable. Examples:\n",
    "  * thermal motion\n",
    "  * games of chance\n",
    "  * radioactive decay\n",
    "\n",
    "How do we deal with randomness numerically?\n",
    "  * Computers are deterministic; no chance involved.\n",
    "  * Always same output for same input; unless hardware error.\n",
    "  * Generate pseudo-random numbers.\n",
    "  \n",
    "Monte Carlo calculations: simulate random events.\n",
    "  * Solve equations statistically\n",
    "  * Solve intractable problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Drawing Random Numbers\n",
    "  * Python has a random module for drawing random numbers\n",
    "  * `random.random()` draws random numbers in `[0, 1)` :  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import random "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['BPF',\n",
       " 'LOG4',\n",
       " 'NV_MAGICCONST',\n",
       " 'RECIP_BPF',\n",
       " 'Random',\n",
       " 'SG_MAGICCONST',\n",
       " 'SystemRandom',\n",
       " 'TWOPI',\n",
       " '_BuiltinMethodType',\n",
       " '_MethodType',\n",
       " '_Sequence',\n",
       " '_Set',\n",
       " '__all__',\n",
       " '__builtins__',\n",
       " '__cached__',\n",
       " '__doc__',\n",
       " '__file__',\n",
       " '__loader__',\n",
       " '__name__',\n",
       " '__package__',\n",
       " '__spec__',\n",
       " '_acos',\n",
       " '_ceil',\n",
       " '_cos',\n",
       " '_e',\n",
       " '_exp',\n",
       " '_inst',\n",
       " '_log',\n",
       " '_pi',\n",
       " '_random',\n",
       " '_sha512',\n",
       " '_sin',\n",
       " '_sqrt',\n",
       " '_test',\n",
       " '_test_generator',\n",
       " '_urandom',\n",
       " '_warn',\n",
       " 'betavariate',\n",
       " 'choice',\n",
       " 'expovariate',\n",
       " 'gammavariate',\n",
       " 'gauss',\n",
       " 'getrandbits',\n",
       " 'getstate',\n",
       " 'lognormvariate',\n",
       " 'normalvariate',\n",
       " 'paretovariate',\n",
       " 'randint',\n",
       " 'random',\n",
       " 'randrange',\n",
       " 'sample',\n",
       " 'seed',\n",
       " 'setstate',\n",
       " 'shuffle',\n",
       " 'triangular',\n",
       " 'uniform',\n",
       " 'vonmisesvariate',\n",
       " 'weibullvariate']"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dir(random)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.37314316199120035"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random.random()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7225804702953376"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random.random()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9863691538857358"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random.random()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.451502960567976,\n",
       " 0.9339710674906362,\n",
       " 0.924059436730153,\n",
       " 0.04034346018637547,\n",
       " 0.4917288063196501,\n",
       " 0.5928789817332918,\n",
       " 0.09606744546454482,\n",
       " 0.5564895172881258,\n",
       " 0.9714023594754505,\n",
       " 0.732155107408455]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[random.random() for i in range(10)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  * The sequence of random numbers is produced by a deterministic algorithm - the numbers just appear random."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distributions of Random Numbers\n",
    "  * `random.random()` generates random numbers that are uniformly\n",
    "distributed in the interval $[0, 1)$.\n",
    "  * `random.uniform(a, b)` generates random numbers uniformly\n",
    "distributed in $[a, b)$.\n",
    "  * \"Uniformly distributed\" means that if we generate a large set\n",
    "of numbers, no part of $[a, b)$ gets more numbers than others.\n",
    "  * See `random_distribution` demo."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vectorized Drawing of Random Numbers\n",
    "  * `random.random()` generates one number at a time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.513089412173477"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random.random() # one number  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-46.075876850491014"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random.uniform(-100, 10) # one number    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  * `numpy` has a random module that efficiently generates a (large)\n",
    "number of random numbers at a time:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "r = np.random.random(size=1000) # array\n",
    "r = np.random.uniform(-1, 10, size=1000) # array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  98.,   85.,  110.,   97.,  124.,   88.,  104.,  101.,  100.,   93.]),\n",
       " array([-0.99302548,  0.10434511,  1.2017157 ,  2.29908629,  3.39645688,\n",
       "         4.49382747,  5.59119807,  6.68856866,  7.78593925,  8.88330984,\n",
       "         9.98068043]),\n",
       " <a list of 1 Patch objects>)"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEXpJREFUeJzt3WuspVV9x/HvT0ZUahGpYYBhGoaGiUDUKgGptfFU0RBq\nwV7CJdWOSpoYvGEaK6NRp2+IrfHWNr7wAp2mQouIBBovnCInxdQ6WhBGhumAZZTBzsF6i05rhPDv\ni/1gD4eZM7MvZz9n1nw/yWT2Ws9l/XfOmd+svfbz7J2qQpLUlif1XYAkafIMd0lqkOEuSQ0y3CWp\nQYa7JDXIcJekBi0Z7kmuTDKfZOtetv1pkkeTHL2gb2OSe5NsT/KK5ShYkrR/+5u5XwWcs7gzyVrg\n5cC3F/SdClwInNod89EkvjKQpB4sGb5VdRvww71s+iDwZ4v6zgeuqaqHq2oncB9w5iSKlCQNZ+iZ\ndZLzgV1VddeiTccDuxa0dwFrxqhNkjSiVcPsnOQI4J0MlmR+0b3EIX62gST1YKhwB34NOBG4MwnA\nCcC/J3kh8CCwdsG+J3R9j5PEwJekEVTVUpPpxxkq3KtqK7D6sXaS+4HTq+oHSW4Erk7yQQbLMScD\nW8Yt8GCTZFNVbeq7juXi8zu4tfz8Wn5uMPzEeH+XQl4D/CuwPskDSV63aJdfDFZV24BrgW3A54FL\ny4+clKReLDlzr6qL97P9pEXtK4ArJlCXJGkMXoc+eXN9F7DM5vouYJnN9V3AMpvru4BlNNd3AStJ\npr1ykqRaXnOXpOUwbHY6c5ekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ\n7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KAlvyBbOlgkPB84vccS/rGK\nn/Q4vvQ4hrta8UrgD4Cv9zD2RcAtYLhr5TDc1ZKbqnj3tAdNeNm0x5T2Z8k19yRXJplPsnVB3/uT\n3JPkziTXJ3nGgm0bk9ybZHuSVyxn4ZKkfdvfG6pXAecs6rsZOK2qngfsADYCJDkVuBA4tTvmo0l8\nw1aSerBk+FbVbcAPF/XNVtWjXfOrwAnd4/OBa6rq4araCdwHnDnZciVJB2LcmfXrgc91j48Hdi3Y\ntgtYM+b5JUkjGPkN1STvAn5eVVcvsVvt49hNC5pzVTU3ah2S1KIkM8DMqMePFO5JXgucC4+7SuBB\nYO2C9gld3xNU1aZRxpWkQ0U36Z17rJ3kvcMcP/SyTJJzgLcD51fVzxZsuhG4KMnhSdYBJwNbhj2/\nJGl8S87ck1wDvAR4VpIHgPcyuDrmcGA2CcBXqurSqtqW5FpgG/AIcGlV7XVZRpK0vDLt/E1SVZWp\nDqrmJbwbOLynm5juB15axf3THluHjmGz0+vQJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lq\nkOEuSQ0y3CWpQYa7JDXI71BtUMKR9PezrarHf8GLpOkz3Nv0aeBFwMNTHvdJwB78khapd4Z7u/6w\nii9Oc8CENfgxz9KK4Jq7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkN8jp36SCVsA54bY8l\nXF/FnT2OryU4c5cOXicCG4BHe/jze8Bzlv0ZamTO3KWD2/1V/Pm0B004edpjajhLztyTXJlkPsnW\nBX1HJ5lNsiPJzUmOWrBtY5J7k2xP8orlLFyStG/7W5a5CjhnUd/lwGxVrQdu6dokORW4EDi1O+aj\nSVz2kaQeLLksU1W3JTlxUfd5wEu6x5uBOQYBfz5wTVU9DOxMch9wJvBvE6xX0srxjoTX9DDuI1X8\nTg/jHlRGWXNfXVXz3eN5YHX3+HgeH+S78KNfpVb9JXBcD+M+Gfhsd6VQH35SxX/3NPZQxnpDtaoq\nSS21y946k2xa0Jyrqrlx6pA0XVXcBdw17XETDmcwcfzStMcGjmTwXQlvmMZgSWaAmVGPHyXc55Mc\nW1W7kxwHPNT1PwisXbDfCV3fE1TVphHGlXSIq+Ln0M+sPeENwK9Pa7xu0jv3/+PnvcMcP8obnjcy\nuLaW7u8bFvRflOTwJOuAk/GLGySpF0vO3JNcw+DN02cleQB4D/A+4NoklwA7gQsAqmpbkmuBbcAj\nwKVVtdSSzVQk/BdwRE/Df7aq1zsIJR2i9ne1zMX72HT2Pva/Arhi3KIm7EgGryL2THnc3+eJl5Ee\nCo5I2NfvzXJ6LrC9h3GlFelQuUP1x1XTDfeE/5nmeCvE/wJfYHC57LQ9Amzd717SIeJQCXdNQRU/\ngF5m7ZIW8Q5SSWqQ4S5JDTLcJalBhrskNcg3VCXpwP1JcnDcu2K4S9KB+QTwtz2N/WXg9GEOMNwl\n6QBU8QiD+ymmLuHRYY9xzV2SGmS4S1KDDHdJapDhLkkNMtwlqUFeLSNNxn8mvYw718uoWvEMd2l8\nJ/VdgLSY4S6NqWrvXwQv9ck1d0lqkDP35fUrCS/sYdyjehhT0gpiuC+fHwC/DPxVT+P/uKdxJa0A\nhvsyqWIWmO27DkmHJtfcJalBhrskNWjkcE/ytiTfTLI1ydVJnpLk6CSzSXYkuTmJb+xJUg9GCvck\na4A3A6dX1XOAw4CLgMuB2apaD9zStSVJUzbOsswq4Igkq4AjgO8C5wGbu+2bgVeNV54kaRQjhXtV\nPQh8APgOg1D/UVXNAqurar7bbR5YPZEqJUlDGelSyCTPZDBLP5HB9dSfTvLqhftUVSXZ623ZSTYt\naM5V1dwodUhSq5LMADOD1mVrhj1+1Ovczwbur6rvd0VcD/wGsDvJsVW1O8lxwEN7O7iqNo04riQd\nErpJ7xxAwrnw4eOHOX7UNfdvA2cleVqSMAj7bcBNwIZunw3ADSOeX5I0hpFm7lW1Jcl1wO0Mvg38\nduBjDG63vzbJJcBO4IIJ1SlJGkKqpvtppUmqqqb2tQYJe4BjqtgzrTElaZIStkDOGCY7vUNVkhpk\nuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7\nJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAaNHO5JjkpyXZJ7kmxL\n8sIkRyeZTbIjyc1JjppksZKkAzPOzP0jwOeq6hTgucB24HJgtqrWA7d0bUnSlKWqhj8oeQZwR1Wd\ntKh/O/CSqppPciwwV1XPXrRPVVXGKXq4WtkDHFPFnmmNKUmTlLAFcsYw2TnqzH0d8L0kVyW5PcnH\nk/wSsLqq5rt95oHVI55fkjSGVWMc9wLgTVX1tSQfZtESTFVVkr2+LEiyaUFzrqrmRqxDkpqUZAaY\nGbQuWzP08SMuyxwLfKWq1nXtFwMbgZOA366q3UmOA251WUaSxjO1ZZmq2g08kGR913U2cDdwE7Ch\n69sA3DDK+SVJ4xl1WQbgzcCnkhwOfAt4HXAYcG2SS4CdwAVjVyhJGtrI4V5VdwJn7GXT2aOXI0ma\nBO9QlaQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KD\nDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBY4V7ksOS\n3JHkpq59dJLZJDuS3JzkqMmUKUkaxrgz97cC24Dq2pcDs1W1Hrila0uSpmzkcE9yAnAu8AkgXfd5\nwObu8WbgVWNVJ0kayTgz9w8BbwceXdC3uqrmu8fzwOoxzi9JGtGqUQ5K8krgoaq6I8nM3vapqkpS\ne9uWZNOC5lxVzY1ShyS1qsvWmUHrsjVDH1+11/zd36BXAK8BHgGeChwJXA+cAcxU1e4kxwG3VtWz\nFx1bVZXF51wuCXuAY6rYM60xJWmSErZAzhgmO0dalqmqd1bV2qpaB1wEfKmqXgPcCGzodtsA3DDK\n+SVJ45nUde6PTf/fB7w8yQ7gpV1bkjRlIy3LjDWgyzKSNJSpLctIklY2w12SGmS4S1KDDHdJapDh\nLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S\n1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDVopHBPsjbJrUnuTvLNJG/p+o9OMptkR5Kbkxw1\n2XIlSQdi1Jn7w8Dbquo04CzgjUlOAS4HZqtqPXBL15YkTdlI4V5Vu6vqG93jnwL3AGuA84DN3W6b\ngVdNokhJ0nDGXnNPciLwfOCrwOqqmu82zQOrxz2/JGl4q8Y5OMnTgc8Ab62qnyT5xbaqqiS1j+M2\nLWjOVdXcOHVIUmuSzAAzg9Zla4Y+vmqv+XsgAz8Z+Cfg81X14a5vOzBTVbuTHAfcWlXPXnRcVVWe\neMblkbAHOKaKPdMaU5ImKWEL5IxhsnPUq2UCfBLY9liwd24ENnSPNwA3jHJ+SdJ4Rpq5J3kx8C/A\nXcBjJ9gIbAGuBX4V2AlcUFU/WnSsM3dJGsIoM/eR1tyr6svse9Z/9ijnlCRNjneoSlKDDHdJapDh\nLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S\n1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBk083JOck2R7knuTvGPS55ck7d9Ewz3J\nYcDfAOcApwIXJzllkmOsdElm+q5hOfn8Dm4tP7+Wn9soJj1zPxO4r6p2VtXDwD8A5094jJVupu8C\nltlM3wUss5m+C1hmM30XsIxm+i5gJZl0uK8BHljQ3tX1SZKmaNLhXhM+3yR8A3i07yIkaQzbhz0g\nVZPL4yRnAZuq6pyuvRF4tKr+YsE+K/E/AEla8aoqB7rvpMN9FfAfwMuA7wJbgIur6p6JDSJJ2q9V\nkzxZVT2S5E3AF4HDgE8a7JI0fROduUuSVoZe7lBN8v4k9yS5M8n1SZ7RRx2T1vINXEnWJrk1yd1J\nvpnkLX3XNGlJDktyR5Kb+q5l0pIcleS67t/dtu79sWYkeVv3e7k1ydVJntJ3TeNIcmWS+SRbF/Qd\nnWQ2yY4kNyc5aqlz9PXxAzcDp1XV84AdwMae6piYQ+AGroeBt1XVacBZwBsbe34AbwW2sTKv+hrX\nR4DPVdUpwHOBZpZLk6wB3gycXlXPYbAkfFG/VY3tKgZZstDlwGxVrQdu6dr71Eu4V9VsVT12eeJX\ngRP6qGPCmr6Bq6p2V9U3usc/ZRAOx/db1eQkOQE4F/gEcMBXJBwMulfGv1VVV8LgvbGq+nHPZU3a\nKuCI7qKOI4AHe65nLFV1G/DDRd3nAZu7x5uBVy11jpXwwWGvBz7XdxETcMjcwJXkROD5DP5jbsWH\ngLfT5j0R64DvJbkqye1JPp7kiL6LmpSqehD4APAdBlfp/aiq/rnfqpbF6qqa7x7PA6uX2nnZwr1b\nG9q6lz+/u2CfdwE/r6qrl6uOKWrxpfwTJHk6cB3w1m4Gf9BL8krgoaq6g8Zm7Z1VwAuAj1bVC4A9\n7Ocl/cEkyTMZzGpPZPBq8ulJ/qjXopZZDa6EWTJzJnop5KLBX77U9iSvZfAy+GXLVcOUPQisXdBe\ny2D23owkTwY+A/x9Vd3Qdz0T9CLgvCTnAk8Fjkzyd1X1xz3XNSm7gF1V9bWufR0NhTtwNnB/VX0f\nIMn1DH6mn+q1qsmbT3JsVe1Ochzw0FI793W1zDkMXgKfX1U/66OGZfB14OQkJyY5HLgQuLHnmiYm\nSYBPAtuq6sN91zNJVfXOqlpbVesYvBH3pYaCnaraDTyQZH3XdTZwd48lTdq3gbOSPK37PT2bwRvj\nrbkR2NA93gAsOcFatpn7fvw1cDgwO/hZ8JWqurSnWibiELiB6zeBVwN3Jbmj69tYVV/osabl0uIS\n25uBT3UTj28Br+u5nompqi1JrgNuBx7p/v5Yv1WNJ8k1wEuAZyV5AHgP8D7g2iSXADuBC5Y8hzcx\nSVJ7VsLVMpKkCTPcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lq0P8B4pyz3o5V4ZAAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108dae470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(r, histtype='step')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.99302548, -0.98761579, -0.9602745 , -0.9504779 , -0.94450209,\n",
       "       -0.93288549, -0.93129153, -0.93044691, -0.92091803, -0.91877624,\n",
       "       -0.90329292, -0.89441737, -0.88126331, -0.87417934, -0.87031556,\n",
       "       -0.8622009 , -0.86037184, -0.83505392, -0.82326597, -0.81464077,\n",
       "       -0.80345092, -0.78893974, -0.76613161, -0.75548421, -0.75131649,\n",
       "       -0.72850547, -0.64359718, -0.63908161, -0.63461537, -0.61539378,\n",
       "       -0.60177248, -0.59582897, -0.59176409, -0.58469372, -0.56681861,\n",
       "       -0.53118404, -0.52476354, -0.5193742 , -0.50592642, -0.49907395,\n",
       "       -0.49630181, -0.49432204, -0.49024554, -0.47831178, -0.47488964,\n",
       "       -0.42993944, -0.42489313, -0.41651328, -0.41219473, -0.40360417,\n",
       "       -0.40329467, -0.40243511, -0.39000204, -0.38664694, -0.36382394,\n",
       "       -0.33409915, -0.31499019, -0.31497072, -0.29898558, -0.24197389,\n",
       "       -0.23036803, -0.22842294, -0.2185465 , -0.21137609, -0.19478051,\n",
       "       -0.18556998, -0.18215362, -0.17657992, -0.17101128, -0.16853141,\n",
       "       -0.16626663, -0.14268024, -0.13080945, -0.12793159, -0.12247863,\n",
       "       -0.11618439, -0.10837897, -0.08831257, -0.08722219, -0.08215761,\n",
       "       -0.06580178, -0.05476295, -0.04859993, -0.04055641, -0.03114233,\n",
       "       -0.01040139,  0.0271867 ,  0.0346155 ,  0.04434243,  0.05863683,\n",
       "        0.05934531,  0.06010559,  0.06356934,  0.06575736,  0.07825514,\n",
       "        0.08326836,  0.08865348,  0.10069934,  0.1106506 ,  0.12850612,\n",
       "        0.13796046,  0.14642767,  0.15703395,  0.16798237,  0.18466432,\n",
       "        0.19482645,  0.20685475,  0.2091683 ,  0.22097414,  0.22434922,\n",
       "        0.22483589,  0.23070757,  0.2328387 ,  0.2328963 ,  0.24219039,\n",
       "        0.25435424,  0.25666572,  0.2741592 ,  0.29635564,  0.31366531,\n",
       "        0.3207991 ,  0.34294595,  0.35007365,  0.36613707,  0.38654263,\n",
       "        0.42036536,  0.42076139,  0.42659311,  0.43340417,  0.44594389,\n",
       "        0.46111945,  0.46794417,  0.47976155,  0.49380501,  0.4989108 ,\n",
       "        0.52177268,  0.54542165,  0.54971106,  0.57635236,  0.59819126,\n",
       "        0.60775473,  0.63903606,  0.64666217,  0.69817215,  0.73113155,\n",
       "        0.73482239,  0.7649351 ,  0.76878568,  0.78105862,  0.78647773,\n",
       "        0.78769031,  0.80093109,  0.81392291,  0.83286755,  0.84009354,\n",
       "        0.84327425,  0.85247339,  0.86321697,  0.88033215,  0.88265485,\n",
       "        0.88837148,  0.89718892,  0.91181223,  0.91857385,  0.91946216,\n",
       "        0.93578302,  0.94172204,  0.94924698,  0.98224359,  0.99657498,\n",
       "        1.00254728,  1.00583193,  1.01099652,  1.04558437,  1.07100156,\n",
       "        1.07411782,  1.07534828,  1.07640162,  1.08121212,  1.08127683,\n",
       "        1.0860841 ,  1.11689448,  1.12053658,  1.2105081 ,  1.21769802,\n",
       "        1.23093629,  1.2352844 ,  1.26092281,  1.26544522,  1.27332789,\n",
       "        1.27460466,  1.28780781,  1.28852284,  1.29811244,  1.33450294,\n",
       "        1.33548936,  1.33738745,  1.36243033,  1.36382718,  1.36816244,\n",
       "        1.37864178,  1.40691631,  1.4103135 ,  1.41317465,  1.41816859,\n",
       "        1.44910411,  1.44974081,  1.45482785,  1.45487091,  1.46236382,\n",
       "        1.46828587,  1.48195397,  1.4830075 ,  1.49493724,  1.49561246,\n",
       "        1.50658984,  1.50803798,  1.51898021,  1.52031372,  1.53877616,\n",
       "        1.53957621,  1.56169379,  1.56176626,  1.56713315,  1.56725488,\n",
       "        1.57031373,  1.57795916,  1.58519914,  1.61848347,  1.63993703,\n",
       "        1.64258133,  1.6441278 ,  1.65979307,  1.66714936,  1.67012093,\n",
       "        1.68122139,  1.68900176,  1.72354172,  1.72500105,  1.7277881 ,\n",
       "        1.74660063,  1.7478925 ,  1.77347745,  1.79684394,  1.82529071,\n",
       "        1.87544253,  1.8874441 ,  1.88898948,  1.90136955,  1.90137727,\n",
       "        1.92423732,  1.93023455,  1.93254453,  1.93381868,  1.93698218,\n",
       "        1.93857302,  1.94434109,  1.96196478,  1.96251642,  1.96663638,\n",
       "        1.96781486,  1.96871063,  1.97165793,  1.97314507,  1.99366611,\n",
       "        2.00285865,  2.01855542,  2.03506669,  2.03578982,  2.04648342,\n",
       "        2.06097929,  2.06992234,  2.07003147,  2.07498643,  2.07820517,\n",
       "        2.08358242,  2.11097901,  2.11886097,  2.11951184,  2.13762861,\n",
       "        2.14937   ,  2.1775655 ,  2.18476807,  2.18511659,  2.20481689,\n",
       "        2.21405699,  2.21480143,  2.22498873,  2.24093971,  2.2501798 ,\n",
       "        2.25496344,  2.27960297,  2.28509968,  2.30080143,  2.31722164,\n",
       "        2.3265428 ,  2.34209157,  2.38079488,  2.40931925,  2.41186993,\n",
       "        2.42135217,  2.43002329,  2.44924833,  2.4513168 ,  2.46192834,\n",
       "        2.47177654,  2.48064659,  2.48072282,  2.49466217,  2.50275824,\n",
       "        2.50614319,  2.51629533,  2.52646307,  2.53455854,  2.53648542,\n",
       "        2.53669508,  2.54745606,  2.56712292,  2.59977575,  2.60661116,\n",
       "        2.61882795,  2.63305092,  2.64358774,  2.66234795,  2.7442657 ,\n",
       "        2.75013747,  2.79942497,  2.83210515,  2.83739957,  2.85544831,\n",
       "        2.86075468,  2.86121855,  2.86289869,  2.86322309,  2.89229484,\n",
       "        2.90590677,  2.91163534,  2.92326474,  2.92511494,  2.92804672,\n",
       "        2.93997324,  2.94116476,  2.94720065,  2.96624002,  2.9939288 ,\n",
       "        2.99481643,  2.99712887,  3.00696545,  3.01824508,  3.02124157,\n",
       "        3.03335715,  3.07688402,  3.07977449,  3.09007557,  3.11182355,\n",
       "        3.11898895,  3.12523334,  3.13076983,  3.13964602,  3.14104221,\n",
       "        3.14789281,  3.15210504,  3.16372793,  3.16452898,  3.16563252,\n",
       "        3.17966105,  3.18704286,  3.19973241,  3.20211107,  3.21443839,\n",
       "        3.23377837,  3.24241696,  3.24669988,  3.24996887,  3.25460236,\n",
       "        3.2692675 ,  3.27407961,  3.2755814 ,  3.2894182 ,  3.30929982,\n",
       "        3.33026954,  3.33569983,  3.35993685,  3.36517984,  3.3664525 ,\n",
       "        3.37820456,  3.37946803,  3.379594  ,  3.38293924,  3.38344803,\n",
       "        3.39860119,  3.4033102 ,  3.40644611,  3.40780676,  3.40854883,\n",
       "        3.40988497,  3.41595605,  3.41616183,  3.41942086,  3.43282171,\n",
       "        3.43304325,  3.43612179,  3.43771893,  3.44037982,  3.478903  ,\n",
       "        3.4808821 ,  3.48134873,  3.4950097 ,  3.50969797,  3.50977555,\n",
       "        3.51689419,  3.55808219,  3.56054313,  3.56637207,  3.57036934,\n",
       "        3.59985866,  3.61395524,  3.63494484,  3.64104055,  3.67781521,\n",
       "        3.70765751,  3.73820079,  3.73890685,  3.74369413,  3.76485697,\n",
       "        3.76786032,  3.77714216,  3.77965629,  3.78226656,  3.78374364,\n",
       "        3.79583361,  3.80172158,  3.80964068,  3.81818473,  3.82456057,\n",
       "        3.83167607,  3.8405459 ,  3.85062973,  3.85109141,  3.85722256,\n",
       "        3.86268967,  3.86707031,  3.88093453,  3.88498591,  3.89676347,\n",
       "        3.89695161,  3.90062641,  3.90759241,  3.91732907,  3.91809858,\n",
       "        3.91843295,  3.93334532,  3.96534684,  3.97532705,  3.97727972,\n",
       "        3.98680749,  3.99518802,  3.99568931,  3.9967673 ,  3.99819069,\n",
       "        4.01126262,  4.01842991,  4.02871705,  4.0295221 ,  4.03006651,\n",
       "        4.03678536,  4.04437063,  4.06500976,  4.10418514,  4.10455721,\n",
       "        4.10660005,  4.11110087,  4.11504492,  4.1390256 ,  4.14060583,\n",
       "        4.15807223,  4.16357077,  4.18703546,  4.18790091,  4.19079016,\n",
       "        4.19782654,  4.2037529 ,  4.20618829,  4.20682949,  4.22196623,\n",
       "        4.24520397,  4.24899168,  4.26116805,  4.26819544,  4.27236684,\n",
       "        4.2777076 ,  4.27849923,  4.29294776,  4.29333736,  4.29446267,\n",
       "        4.30010381,  4.31433847,  4.31653174,  4.37493882,  4.37907691,\n",
       "        4.38284434,  4.38942165,  4.39982957,  4.40089832,  4.40196995,\n",
       "        4.40782448,  4.41918245,  4.41921246,  4.42060354,  4.44569947,\n",
       "        4.4705891 ,  4.478532  ,  4.48487554,  4.49335988,  4.49986905,\n",
       "        4.50116726,  4.50188148,  4.50715275,  4.51553438,  4.54090919,\n",
       "        4.54176312,  4.56130746,  4.57658898,  4.58300488,  4.58433619,\n",
       "        4.61212893,  4.6484168 ,  4.65298157,  4.65556208,  4.71122108,\n",
       "        4.71435134,  4.71628742,  4.72446679,  4.72619434,  4.72827452,\n",
       "        4.75685436,  4.75983168,  4.77451367,  4.7877414 ,  4.7932592 ,\n",
       "        4.79330405,  4.79466082,  4.80225924,  4.82066104,  4.82737205,\n",
       "        4.84064285,  4.87405855,  4.87682458,  4.87862292,  4.89801202,\n",
       "        4.90874204,  4.92971624,  4.93943032,  4.95746367,  4.96042303,\n",
       "        4.96418509,  4.97949317,  5.00429427,  5.01790888,  5.0263083 ,\n",
       "        5.02948494,  5.03142348,  5.03905564,  5.05974746,  5.0809539 ,\n",
       "        5.08140786,  5.08590667,  5.10820382,  5.11638198,  5.11934471,\n",
       "        5.1309116 ,  5.15151982,  5.15183362,  5.16071022,  5.16979866,\n",
       "        5.17950765,  5.18745475,  5.21370989,  5.21943443,  5.22142602,\n",
       "        5.22330833,  5.24417853,  5.24652016,  5.26421405,  5.28300046,\n",
       "        5.30442897,  5.36590348,  5.38890007,  5.39049744,  5.40669979,\n",
       "        5.42947444,  5.44599887,  5.47413222,  5.47628409,  5.47964832,\n",
       "        5.52048454,  5.55117336,  5.55140186,  5.55461539,  5.56761217,\n",
       "        5.56772336,  5.57510764,  5.59222983,  5.60255261,  5.62265021,\n",
       "        5.6251262 ,  5.66327222,  5.67165663,  5.67551702,  5.68020178,\n",
       "        5.68777784,  5.70611061,  5.71123242,  5.71406475,  5.71606435,\n",
       "        5.7167724 ,  5.71714628,  5.74133742,  5.75552163,  5.75900489,\n",
       "        5.75901702,  5.75987269,  5.76165913,  5.7621241 ,  5.76640733,\n",
       "        5.76650666,  5.78414961,  5.78880976,  5.7897779 ,  5.79433917,\n",
       "        5.79776106,  5.82487799,  5.82840532,  5.85695345,  5.85965459,\n",
       "        5.86976701,  5.87236359,  5.87499143,  5.92174253,  5.93582815,\n",
       "        5.95922626,  5.96409919,  5.97831934,  5.99270585,  6.00387387,\n",
       "        6.00617428,  6.02735325,  6.03407275,  6.03491933,  6.04073053,\n",
       "        6.04958675,  6.05131793,  6.05177833,  6.05352525,  6.08682484,\n",
       "        6.11258122,  6.11778407,  6.11952545,  6.12708142,  6.13225292,\n",
       "        6.14358294,  6.16403484,  6.16645624,  6.17940485,  6.18001339,\n",
       "        6.18023818,  6.18855754,  6.19491536,  6.19509139,  6.20895976,\n",
       "        6.22897334,  6.23423949,  6.23956267,  6.29365634,  6.30229683,\n",
       "        6.30242582,  6.30789518,  6.36277554,  6.36647159,  6.37113296,\n",
       "        6.38316288,  6.38806258,  6.39972202,  6.42798239,  6.43286258,\n",
       "        6.46409395,  6.47206933,  6.50454567,  6.50724964,  6.5080522 ,\n",
       "        6.53220893,  6.53837566,  6.5410563 ,  6.56644477,  6.59129313,\n",
       "        6.5962785 ,  6.60309159,  6.60641512,  6.61816523,  6.6209375 ,\n",
       "        6.64209317,  6.64246855,  6.65465582,  6.6555963 ,  6.66177035,\n",
       "        6.68782212,  6.69464116,  6.6976801 ,  6.71421595,  6.71768814,\n",
       "        6.72148452,  6.72425742,  6.72543708,  6.72820771,  6.73552621,\n",
       "        6.74468886,  6.77311125,  6.77963944,  6.78173483,  6.79516906,\n",
       "        6.79564464,  6.7978205 ,  6.81300365,  6.82313698,  6.86614477,\n",
       "        6.88126558,  6.88914172,  6.90296344,  6.90952665,  6.91583399,\n",
       "        6.92079094,  6.9413799 ,  6.99428537,  7.01193997,  7.01687803,\n",
       "        7.01969454,  7.02491808,  7.02888456,  7.04402929,  7.06018947,\n",
       "        7.06263313,  7.07741312,  7.08684413,  7.09583931,  7.10134811,\n",
       "        7.10698641,  7.11773066,  7.12766694,  7.1309009 ,  7.13378253,\n",
       "        7.15052066,  7.17169111,  7.19982758,  7.20348711,  7.20860298,\n",
       "        7.21087488,  7.22059979,  7.23440789,  7.23685396,  7.24298955,\n",
       "        7.24691467,  7.25885514,  7.28945265,  7.2986341 ,  7.30105988,\n",
       "        7.30650719,  7.32309936,  7.32784936,  7.33409428,  7.34700568,\n",
       "        7.34715856,  7.38564883,  7.3924266 ,  7.43564399,  7.44075708,\n",
       "        7.45153035,  7.45469833,  7.48833843,  7.4985198 ,  7.50697706,\n",
       "        7.51297938,  7.51814808,  7.52931985,  7.53008304,  7.53950438,\n",
       "        7.54264706,  7.54502271,  7.54906863,  7.58634178,  7.59633545,\n",
       "        7.59768725,  7.60238401,  7.60893084,  7.61420897,  7.63195233,\n",
       "        7.63272094,  7.64911965,  7.65380909,  7.66662472,  7.67095506,\n",
       "        7.68175668,  7.71142242,  7.71634537,  7.71886376,  7.72901062,\n",
       "        7.77398603,  7.78430462,  7.79074283,  7.79520159,  7.81025714,\n",
       "        7.81067754,  7.81687116,  7.82338439,  7.84132461,  7.84212571,\n",
       "        7.84676358,  7.85727838,  7.87562444,  7.8838897 ,  7.89027781,\n",
       "        7.89274851,  7.8997841 ,  7.92924733,  7.93564714,  7.93569241,\n",
       "        7.93917859,  7.96300319,  7.96520327,  8.01384587,  8.02045858,\n",
       "        8.02122143,  8.05382343,  8.06237501,  8.06418277,  8.06622413,\n",
       "        8.07714655,  8.09952065,  8.10831258,  8.11175841,  8.11359285,\n",
       "        8.12767111,  8.12810255,  8.13244375,  8.13729263,  8.14442299,\n",
       "        8.14626581,  8.15967977,  8.18638371,  8.19207641,  8.23960399,\n",
       "        8.25417199,  8.2631335 ,  8.26977139,  8.27448342,  8.28667759,\n",
       "        8.29268367,  8.30278935,  8.30618804,  8.3254671 ,  8.33927012,\n",
       "        8.35594983,  8.36076167,  8.36601487,  8.37432284,  8.38662493,\n",
       "        8.39785867,  8.41376365,  8.4236327 ,  8.42642705,  8.43836477,\n",
       "        8.44953965,  8.44984592,  8.46459355,  8.47807097,  8.49910768,\n",
       "        8.504441  ,  8.50769712,  8.52392476,  8.53798986,  8.54157672,\n",
       "        8.58011848,  8.58838577,  8.60467633,  8.63685036,  8.64451769,\n",
       "        8.6493155 ,  8.65271944,  8.65309001,  8.65955551,  8.66967407,\n",
       "        8.67590353,  8.69397108,  8.6991573 ,  8.70130835,  8.71372039,\n",
       "        8.7166161 ,  8.727886  ,  8.73062748,  8.73194627,  8.75173617,\n",
       "        8.75874406,  8.77232918,  8.80793995,  8.83772145,  8.85071009,\n",
       "        8.86250731,  8.87087518,  8.88410271,  8.90605   ,  8.91802326,\n",
       "        8.91894298,  8.92090063,  8.95159919,  8.97598753,  9.00018232,\n",
       "        9.03299192,  9.03583258,  9.0636271 ,  9.075008  ,  9.10561889,\n",
       "        9.12233661,  9.12965742,  9.13941225,  9.1410814 ,  9.14586842,\n",
       "        9.15295467,  9.15337524,  9.17052126,  9.17756084,  9.18858413,\n",
       "        9.19304371,  9.2057327 ,  9.20650682,  9.21298403,  9.21723901,\n",
       "        9.22810619,  9.24033569,  9.25808615,  9.28542748,  9.29304549,\n",
       "        9.29808327,  9.2995773 ,  9.3107325 ,  9.31586744,  9.33273332,\n",
       "        9.34376584,  9.3584737 ,  9.362098  ,  9.36507839,  9.37623514,\n",
       "        9.37837956,  9.39201238,  9.39248528,  9.39605497,  9.39857888,\n",
       "        9.41251534,  9.4300362 ,  9.45232435,  9.45798445,  9.46686942,\n",
       "        9.50437243,  9.50881302,  9.51142108,  9.52421792,  9.53077019,\n",
       "        9.54555577,  9.56853476,  9.58270109,  9.60222784,  9.63693751,\n",
       "        9.64546635,  9.64785008,  9.65145651,  9.66724008,  9.67762271,\n",
       "        9.67937723,  9.68004931,  9.69611067,  9.69961905,  9.70007482,\n",
       "        9.70820123,  9.71527232,  9.76537662,  9.80015862,  9.81000926,\n",
       "        9.83239421,  9.83771072,  9.83891625,  9.85605404,  9.87408516,\n",
       "        9.88537438,  9.89572611,  9.89884315,  9.90431059,  9.92018414,\n",
       "        9.92758251,  9.94480244,  9.96204112,  9.97385309,  9.98068043])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r[np.argsort(r)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1090b4208>]"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEACAYAAABPiSrXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHKxJREFUeJzt3X28VWWd9/HPV4FEJdDRAQQczcCkGlEKfAx8CIkatPIp\nR3OS7FZybHJqAmfupOl+dftwTz7Urc0rn9BRJrREvCUVzVNNpeSEhSICk6Qc5eBYPpYJ+bv/uK7j\n2R4PsM/Ze5+1H77v12u/ztrX3mvt31nK/p7rWutaSxGBmZlZObYrugAzM2scDg0zMyubQ8PMzMrm\n0DAzs7I5NMzMrGwODTMzK1vFoSFpnaRfSVouaVlu21XSUkmrJd0jaVjJ++dKWiNplaRpJe0TJa3I\nr11eaV1mZlZ91ehpBDA1Ig6IiEm5bQ6wNCLGAffl50gaD5wEjAemA1dKUl7nKmBWRIwFxkqaXoXa\nzMysiqo1PKVuz2cC8/PyfOC4vHwssCAiNkXEOmAtMFnSSGBIRCzL77uhZB0zM6sT1epp3CvpIUln\n5rbhEdGRlzuA4Xl5D2B9ybrrgVE9tLfndjMzqyMDqrCNQyPiGUm7A0slrSp9MSJCkq9VYmbWBCoO\njYh4Jv98VtJtwCSgQ9KIiNiQh5425re3A2NKVh9N6mG05+XS9vbun+XwMTPrm4jofhihTyoKDUk7\nAttHxEuSdgKmAV8BFgOnAxfln4vyKouBmyV9nTT8NBZYlnsjL0qaDCwDTgOu6Okzq/WLNzpJ8yJi\nXtF11APviy7eF128L7pU8w/uSnsaw4Hb8glQA4CbIuIeSQ8BCyXNAtYBJwJExEpJC4GVwGZgdnRd\nZnc2cD0wGFgSEXdVWJuZmVVZRaEREU8AE3po/y1w9BbW+RrwtR7a/xN4byX1mJlZbXlGeONqK7qA\nOtJWdAF1pK3oAupIW9EFNCM10k2YJIWPaZiZ9U41vzvd0zAzs7I5NMzMrGwODTMzK5tDw8zMyubQ\nMDOzsjk0zMysbA4NMzMrm0PDzMzK5tAwM7OyOTTMzKxsDg0zMyubQ8PMrIlJVPV6fQ4NM7MmJbET\nUNV7Ezk0zMyakMQg4Fbg6Wpu16FhZtZkJLYDrgNeA86s5rYrvd2rmZnVkXwM41JgT2BaBJtVxaMa\nDg0zs+YyFzgC+EAEf6j2xh0aZmZNQuJMYBZwWATP1+IzHBpmZk1A4mPAV0g9jGdq9TkODTOzBidx\nBPAt4JgI1tbys3z2lJlZA5M4EPgOcGIEy2v9eQ4NM7MGJbEv8P+AsyJo64/PdGiYmTUgib2ApcD5\nEXyvvz7XoWFm1mAkRgL3AhdHcH1/frZDw8ysgUjsRgqMayL4Zn9/vkPDzKxBSAwlXYBwcQT/u4ga\nqhIakraXtFzSHfn5rpKWSlot6R5Jw0reO1fSGkmrJE0raZ8oaUV+7fJq1GVm1iwkdiQd9H4QOL+o\nOqrV0/gcsBKI/HwOsDQixgH35edIGg+cBIwHpgNXSm9cFeUqYFZEjAXGSppepdrMzBqaxNuA24Bf\nA38b8cZ3bb+rODQkjQZmAFfDGzf7mAnMz8vzgePy8rHAgojYFBHrgLXAZEkjgSERsSy/74aSdczM\nWpbEAGAB8DIwK4LXi6ynGj2NS4Evwpt+keER0ZGXO4DheXkPYH3J+9YDo3pob8/tZmYtK1/i/Fpg\nR+CUCDYXXFJloSHpI8DGiFgOPd9SMCICiutKmZk1onyJ828CewEfi+CPxVaUVHrtqUOAmZJmADsA\nb5d0I9AhaUREbMhDTxvz+9uBMSXrjyb1MNrzcml7e08fKGleydO2iGir8HcwM6srOTAuAiYBR0Xw\n+96tr6nA1OpXBkodgSpsSJoCfCEi/krSxcBzEXGRpDnAsIiYkw+E30zaEaNI5xq/MyJC0oPAucAy\n4E7gioi4q9tnRERU9SbpZmb1RuKrpGPAR0TwXOXbq953Z7WvctuZQBcCCyXNAtYBJwJExEpJC0ln\nWm0GZkdXas0GrgcGA0u6B4aZWSuQ+J/Ax6hSYFRb1Xoa/cE9DTNrZhJzgL8BpkawoXrbrd+ehpmZ\n9YHEecCngSnVDIxqc2iYmRVM4hzgHFIPo8eTgOqFQ8PMrEAS/4M0121qBE8WXc+2ODTMzAoicQbw\nT6SD3k8UXU85HBpmZgWQOBX4KnBkre/rXU0ODTOzfiZxEnAJaeLe40XX0xsODTOzfiRxPHA5MC2C\nlUXX01u+CZOZWT/JPYxvAtMj+FXR9fSFQ8PMrB9InAJcRuphPFx0PX3l0DAzqzGJTwL/Bzi6UXsY\nnRwaZmY1JPEp4Gukg96PFl1PpXwg3MysRiTOBL5MA54ltSUODTOzGpA4CzifNHGvYeZhbItDw8ys\nyvK1pL5AujTIr4uup5ocGmZmVSTxd8DnSIGxruByqs6hYWZWBfkWrf8InE66vHndX3ywLxwaZmYV\nyoFxIfBh4AMRPFNwSTXj0DAzq4DEdsA3gMmkHkbd3aK1mhwaZmZ9JDEAuBrYh3Ra7QsFl1RzDg0z\nsz6QGATcBAwlXUvqlYJL6heeEW5m1ksSg4FFwEBgZqsEBjg0zMx6RWIIsAT4HXBCBK8WXFK/cmiY\nmZVJYhdgKbAW+GQEmwouqd85NMzMyiAxEmgDfgZ8JoI/FVtRMRwaZmbbILEP8B/ALcB5EUTBJRXG\nZ0+ZmW2FxATSMYx/juBbRddTNIeGmdkWSEwh9S4+G8EtRddTDxwaZmY9kDgW+DZwSgT3Fl1PvfAx\nDTOzbvLd9r4FzHBgvFlFoSFpB0kPSnpY0iOS5uX2XSUtlbRa0j2ShpWsM1fSGkmrJE0raZ8oaUV+\n7fJK6jIz6yuJLwIXkC5t/lDR9dSbikIjIl4FjoiICcAEYLqkycAcYGlEjAPuy8+RNB44CRgPTAeu\nlKS8uauAWRExFhgraXoltZmZ9YaEJC4B/gY4tFluz1ptFQ9PRcTv8+Ig0pT6AGYC83P7fOC4vHws\nsCAiNkXEOtIEmcmSRgJDImJZft8NJeuYmdWUxEDgGuAw4PAI2gsuqW5VHBqStpP0MNAB3JO/+IdH\nREd+SwcwPC/vAawvWX09MKqH9vbcbmZWUxI7A7cDI4CjI/htwSXVtYrPnoqI14EJkoYCt0l6T7fX\nQ1LVJsJ0HjfJ2iKirVrbNrPWIjECuBP4BXB2BJsLLqkqJE0FptZi21U75TYiXpB0P3AM0CFpRERs\nyENPG/Pb2oExJauNJvUw2vNyaXuP3cOImFetms2sdUnsC3wfuA74X800yzv/Md3W+VzSBdXadqVn\nT+3WeWaUpMHAB4HHgMWk++SSfy7Ky4uBkyUNkrQ3MBZYFhEbgBclTc4Hxk8rWcfMrKokDgV+SJrl\n/dVmCoxaq7SnMRKYL2l7UgB9JyKWSHoAWChpFrAOOBEgIlZKWgisBDYDsyOi8z/WbOB6YDCwJCLu\nqrA2M7O3kPg46WzN0yK4u+h6Go26vrPrn6SICG37nWZmbyVxLvAPwF9FsLzoevpLNb87fRkRM2t6\nEtsBFwMzSHMwflNwSQ3LoWFmTU1iB9J8sZHAYT6ltjK+9pSZNS2J3Uh32hMwzYFROYeGmTUliXcB\nDwI/Bk5utXt514qHp8ys6UgcBSwAvhTBdUXX00zc0zCzpiJxJnAzcIIDo/rc0zCzppDPkLqQdLHT\nwyNYXXBJTcmhYWYNT2In4N+AXYCDI3iu4JKaloenzKyhSYwCfgQ8TzpDyoFRQw4NM2tYEgcAPwNu\nAc6I4LWCS2p6Hp4ys4Yk8VHgX4HZEdxadD2twqFhZg0lH/D+J+DTwAzfx7t/OTTMrGHku+xdT7rb\n56QINhRbUevxMQ0zawgSewE/AV4EjnBgFMOhYWZ1T2IK6YD3tcCsCP5YcEkty8NTZlbXJM4G5gF/\nHcG9BZfT8hwaZlaXJAYBVwCHk+6BsbbgkgyHhpnVIYnhpLkXz5NmeL9YcEmW+ZiGmdUViYOBh4A2\n4DgHRn1xT8PM6oKEgLOAr5AOdt9RcEnWA4eGmRVOYjBwJfA+0vGLNQWXZFvg4SkzK1Sef/EfwA7A\nQQ6M+ubQMLPCSEwDHgBuBE6J4JWCS7Jt8PCUmfW7fP2oOcA5wEkR/LDgkqxMDg0z61cSuwHzSTdM\nmhTB+oJLsl7w8JSZ9RuJQ4FfAI8CUxwYjcc9DTOruTwc9ffAF4BP+3TaxuXQMLOaktiVNBy1O2k4\n6jcFl2QVqGh4StIYSfdLelTSI5LOze27SloqabWkeyQNK1lnrqQ1klZJmlbSPlHSivza5ZXUZWb1\nQeIgYDmwGviAA6PxVXpMYxPw+Yh4N3AQ8FlJ+5HOilgaEeOA+/JzJI0HTgLGA9OBKyUpb+sqYFZE\njAXGSppeYW1mVhAJSZwHLAbOjeDvff/u5lBRaETEhoh4OC+/DDwGjAJmkrqj5J/H5eVjgQURsSki\n1gFrgcmSRgJDImJZft8NJeuYWQOR2AW4DTgZmBzB7QWXZFVUtbOnJO0FHAA8CAyPiI78UgcwPC/v\nAW86W2I9KWS6t7fndjNrIPnsqOXAOuCwCJ4otiKrtqocCJe0M/Bd4HMR8VLXiBNEREiKanxO/qx5\nJU/bIqKtWts2s76RGAD8I3A2cKbPjiqWpKnA1Fpsu+LQkDSQFBg3RsSi3NwhaUREbMhDTxtzezsw\npmT10aQeRnteLm1v7+nzImJepTWbWfVI7AncBLwGHBjB0wWX1PLyH9Ntnc8lXVCtbVd69pSAa4CV\nEXFZyUuLgdPz8unAopL2kyUNkrQ3MBZYFhEbgBclTc7bPK1kHTOrUxLHk+59cScwzYHR/BTR95Ej\nSYcBPwJ+BXRuaC6wDFgI7Eka2zwxIp7P65wPnAFsJg1n3Z3bJwLXA4OBJRFxbg+fFxGh7u1m1r8k\ndgQuA44CPhHBsm2sYgWq5ndnRaHR3xwaZsWT2B9YQLocyGzfWa/+VfO709eeMrOySGwvMQe4F7gw\nglMdGK3HlxExs22S2Ic0f+o14H2e2d263NMwsy3KM7vPIs2/ugU4yoHR2tzTMLMeSexBOjtyd9J1\no1YWXJLVAfc0zOwtJE4izex+EDjYgWGd3NMwszfky5j/X2AC8JEIfl5wSVZn3NMwMwAkPgY8Amwg\nzex2YNhbuKdh1uIk/hz4Bql3cUIEPym4JKtj7mmYtah8ZtTJpCs6rAMmODBsW9zTMGtBEiNJNz57\nJzDTlwGxcrmnYdZCcu/idOCXwApgogPDesM9DbMWIbEXqXcxEjgmguXFVmSNyD0NsyYnMSDfr/sh\n0lWp3+/AsL5yT8OsiUkcCHwbeJ40SW9NwSVZg3NPw6wJSewkcQnwfdLptEc7MKwaHBpmTUbiGNIk\nvZHAeyO4PoLGuXGO1TUPT5k1iTxJ7+vAocDZEdxVcEnWhNzTMGtw+eZIs4FHSZcAeY8Dw2rFPQ2z\nBiYxCbgS+D1wZAQrCi7Jmpx7GmYNSOLPJP4VuB24ApjiwLD+4NAwayB5KOozwEpgE7BfBDf4QLf1\nFw9PmTUIiSOBS4EXgA9F8IuCS7IW5NAwq3MSY4FLgL8Evgh8zz0LK4qHp8zqlMQwiX8BfpYf4yP4\nrgPDiuTQMKsz+VpRZwOPA0OAd0dwUQSvFlyamYenzOqFhIAZwEXAs6Qr0T5cbFVmb+bQMKsDEgeR\nwmJ3YC6w2MNQVo88PGVWIIn9JG4DbgHmA38Zwe0ODKtXFYeGpGsldUhaUdK2q6SlklZLukfSsJLX\n5kpaI2mVpGkl7RMlrcivXV5pXWb1TGK0xNWk+1v8FBgXwbURbC64NLOtqkZP4zpgere2OcDSiBgH\n3JefI2k8cBIwPq9zpSTlda4CZkXEWGCspO7bNGt4ErtIXES63ep/k8Likgj+UHBpZmWpODQi4sfA\n77o1zyR1tck/j8vLxwILImJTRKwD1gKTJY0EhkRE572KbyhZx6zhSQyRmAOsBnYhDUPNiXjLvx2z\nularYxrDI6IjL3cAw/PyHsD6kvetB0b10N6e280aWklY/BewP3B4BJ+JoL3g0sz6pOZnT0VESKra\nQT1J80qetkVEW7W2bVYtEjsDnwXOA34ATI1gZbFVWauQNBWYWott1yo0OiSNiIgNeehpY25vB8aU\nvG80qYfRnpdL23v8Sywi5lW/XLPq6CEsjnBYWH/Lf0y3dT6XdEG1tl2r4anFwOl5+XRgUUn7yZIG\nSdobGAssi4gNwIuSJucD46eVrGNW9yR2lvgSaRhqAiksPuHAsGZTcU9D0gJgCrCbpKeALwMXAgsl\nzQLWAScCRMRKSQtJl3XeDMyOiM6hq9nA9cBgYElE+M5jVvckdgHOyQ/3LKzpqes7u/5JiojQtt9p\nVlsSewCfB84g3Qjp4ghWFVuVWc+q+d3pGeFmvSCxT75j3iPAQGBCBGc4MKxVODTMyiCxv8QC4AHS\naeTjIvi7CJ4quDSzfuXQMNsCCUkcKXEnsAT4T+AdEXw5gv8uuDyzQvgqt2bdSAwCTiadNjuIdIvV\nj/t+FmYODbM3SPwZcBZpnsUK0jXT7vYVZ826eHjKWp7EvhJXka6F9g7SzY+OieAuB4bZm7mnYS0p\n3yXvCNIQ1PuBbwHviqBjqyuatTiHhrWUfJmPU0mT8QRcBpzgS5OblcehYS1BYhzpqgOnAT8E/hZo\n8/CTWe84NKxpSWwPfIjUqzgQuBo4IIInCy3MrIE5NKzpSIwkXSjzTNINwr4BHOdTZs0q59CwpiAx\nEPgwMAs4DLgVOAVY5iEos+pxaFhDk3gX6aKBnySdMnsN8IkIXi60MLMm5dCwhpPPgDqR1Kt4B+me\n8lMieLzQwsxagC+Nbg0hz6s4BPgU8HHgR8C1wJIINhVZm1m9q+Z3p3saVtck9gX+Oj9eI92oa3wE\nzxRZl1mrcmhY3ZEYTrpg4Kmk+8UvAE4AlvugtlmxPDxldUFiJ+A4UlAcTLob3k3ADyLYXGRtZo3O\nw1PWFCR2BKYDxwMzgJ+QDmofH8ErRdZmZj1zT8P6Ve5RzCAFxTHAQ6Q5Fd+LYGORtZk1q2p+dzo0\nrObyKbIfJh2X+CDplqm3AosieLbI2sxagUPD6p7EnqQexQxgCmno6Vbg9gieK7I2s1bj0LC6ky/j\ncQipRzEDGA58n3Rv7bsj+F2B5Zm1NIeGFS5PttsHOJI05HQ06TIeS/LjoQj+VFyFZtbJoWGFkBhF\nCokjgaOA7YH78uMu3/XOrD45NKzmSnoShwCHko5L7AbcD/wgP1Z7sp1Z/XNoWNVJ7ABMJIVE5+OP\nwE/z48fALyN4vbAizaxPHBpWkdyLeAfwPuD9pIDYH3iMFBA/AX4awVOFFWlmVdO0oSFpOnAZaaz8\n6oi4qNvrDo1eygHxF8ABpIB4X368AvycNLnuZ6SbFXkWtlkTasrQkLQ98DjpLJx20hfaJyLisZL3\nODS2QmJ34D3Ae0t+vht4GVhOV0g8FMGGouo0s/7VrNeemgSsjYh1AJL+HTiWNGRimcR2wChgXH7s\nSwqI9wBvAx4BVgAPA/8GPBLBb4up1syaTT2Fxih40xj6emByQbUUKh+UHg2MIQ0tjaUrJN4JvACs\nJvXMVgN3kYLiaZ/NZGa1VE+hUdaXncT9wDP5saHb8gbg+XqdVCYxCPhz0mzp0scoUkB0PoaShuie\nyo/VwHfzzzURvNTvxZuZUV+h0U76wuw0htTb6Ob9T8LQIbDze+HY0fCpPwEjgJH559slfk/6a/wF\n4MWS5ReAPwCvkk4nfbXbcufppFHyKH2+HTBoC4+BwA7AEODtPfwcCuwEPAtsBDryYyPwX0AbKSCe\nBDb61FYz6ytJU4GpNdl2HR0IH0AabjkKeBpYRh8OhOcx/51JX9LdH28HBpPG/nfIj9JldXvQbTlI\nAfNayWNTyfIfgZdIQdX5s3P5BVIvyGFgZv2qKQ+ER8RmSecAd5NOub2mNDDK3w6v0/Vl7XkGZmZV\nVDc9jXL4lFszs96r5nfndtXYiJmZtQaHhpmZlc2hYWZmZXNomJlZ2RwaZmZWNoeGmZmVzaFhZmZl\nc2iYmVnZHBpmZlY2h4aZmZXNoWFmZmVzaJiZWdkcGmZmVjaHhpmZlc2hYWZmZXNomJlZ2RwaZmZW\nNoeGmZmVzaFhZmZlc2iYmVnZHBpmZlY2h4aZmZXNoWFmZmVzaJiZWdkcGmZmVjaHhpmZlc2hYWZm\nZXNomJlZ2focGpJOkPSopD9JOrDba3MlrZG0StK0kvaJklbk1y4vaX+bpO/k9gck/UVf6zIzs9qp\npKexAvgo8KPSRknjgZOA8cB04EpJyi9fBcyKiLHAWEnTc/ss4LncfilwUQV1tQRJU4uuoV54X3Tx\nvujifVEbfQ6NiFgVEat7eOlYYEFEbIqIdcBaYLKkkcCQiFiW33cDcFxengnMz8vfBY7qa10tZGrR\nBdSRqUUXUEemFl1AHZladAHNqBbHNPYA1pc8Xw+M6qG9PbeTfz4FEBGbgRck7VqD2szMrAIDtvai\npKXAiB5eOj8i7qhNSWZmVq+2GhoR8cE+bLMdGFPyfDSph9Gel7u3d66zJ/C0pAHA0Ij4bU8blxR9\nqKkpSbqg6BrqhfdFF++LLt4X1bfV0OgFlSwvBm6W9HXSsNNYYFlEhKQXJU0GlgGnAVeUrHM68ABw\nPHBfTx8SEeqp3czM+kefQ0PSR0lf+rsBd0paHhEfioiVkhYCK4HNwOyI6OwdzAauBwYDSyLirtx+\nDXCjpDXAc8DJfa3LzMxqR13f52ZmZlvXMDPCJU3PkwXXSPpS0fXUkqQxku7PkycfkXRubt9V0lJJ\nqyXdI2lYyTo9TqhsFpK2l7Rc0h35eUvuC0nDJN0q6TFJKyVNbuF98fn872OFpJvzJOGW2BeSrpXU\nIWlFSVuvf/ctTbjeqoio+wewPWm+x17AQOBhYL+i66rh7zsCmJCXdwYeB/YDLgb+Ibd/CbgwL4/P\n+2Rg3kdrge2K/j2qvE/OA24CFufnLbkvSPOZzsjLA4ChrbgvSMdLfw28LT//Dum4aEvsC+Bw4ABg\nRUlbb373zlGmZcCkvLwEmL6tz26UnsYkYG1ErIuITcC/kyYRNqWI2BARD+fll4HHSP9ISidBzqdr\ncmRPEyon9WvRNSRpNDADuJquky5abl9IGgocHhHXQprTFBEv0IL7IhsA7JjPuNwReJoW2RcR8WPg\nd92ae/O7b2vC9RY1Smi8Mfkv65ww2PQk7UX6i+JBYHhEdOSXOoDheXlLEyqbxaXAF4HXS9pacV/s\nDTwr6TpJv5D0bUk70YL7IiLagX8BniSFxfMRsZQW3Bclevu7b23C9RY1Smi05NF6STuTLqvyuYh4\nqfS1SP3Jre2Xpthnkj4CbIyI5bz51O43tMq+IP1lfSBwZUQcCLwCzCl9Q6vsC0m7kP6y3ov05bez\npFNL39Mq+6InZfzufdYoodF9wuAY3pyQTUfSQFJg3BgRi3Jzh6QR+fWRwMbc3tOEyvb+qrXGDgFm\nSnoCWAAcKelGWnNfrAfWR8TP8/NbSSGyoQX3xdHAExHxXKRLD30POJjW3BedevNvYksTrre5Txol\nNB4iXRV3L0mDSFfRXVxwTTWTrwp8DbAyIi4realzEiT556KS9pMlDZK0N3lCZX/VW0sRcX5EjImI\nvUnzd34QEafRmvtiA/CUpHG56WjgUeAOWmxfAL8BDpI0OP97OZo0N6wV90WnXv2byP8/vZjPwBNp\nwvWi7ht9i6LPAujF2QIfIp1FtBaYW3Q9Nf5dDyON3z8MLM+P6cCuwL3AauAeYFjJOufnfbMKOKbo\n36FG+2UKXWdPteS+APYHfg78kvTX9dAW3hfzSCeJrCAd+B3YKvuC1Ot+GniNdLz3U3353YGJef+t\nBa4o57M9uc/MzMrWKMNTZmZWBxwaZmZWNoeGmZmVzaFhZmZlc2iYmVnZHBpmZlY2h4aZmZXNoWFm\nZmX7/5EW/0aVw+2XAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108fc55c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cs = np.cumsum(r[np.argsort(r)])\n",
    "plt.plot(cs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  * Vectorized calculations are important for speeding up programs!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Drawing Integers\n",
    "  * Quite often we want to draw an integer from $[a, b]$ and not a\n",
    "real number.\n",
    "  * Python's `random` module and `numpy.random` have functions for\n",
    "drawing uniformly distributed integers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "(a, b) = (-1, 10)\n",
    "import random  \n",
    "r = random.randint(a, b) # a, a+1, ..., b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n"
     ]
    }
   ],
   "source": [
    "print(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100\n",
    "import numpy as np  \n",
    "r = np.random.randint(a, b+1, N)  # b+1 is not included  \n",
    "r = np.random.random_integers(a, b, N)  # b is included  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5030"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(np.random.random_integers(0, 1, 10000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Throwing a Die; Vectorized Version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5 6 2 3 4 3 5 6 2 2]\n",
      "[False  True False False False False False  True False False]\n",
      "[False  True False False False False False  True False False]\n",
      "Got six 2 times out of 10\n"
     ]
    }
   ],
   "source": [
    " import numpy as np  \n",
    " N = 10   \n",
    " eyes = np.random.randint(1, 6+1, N)  \n",
    "print(eyes)\n",
    "print(eyes == 6)\n",
    "success = (eyes == 6)\n",
    "print(success)\n",
    " \n",
    " # True/False array   \n",
    " success = eyes == 6   \n",
    " # treats True as 1, False as 0   \n",
    " num_six = np.sum(success)   \n",
    " print('Got six %d times out of %d' % (num_six, N)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  * Important: use `numpy.sum` and not the default `sum` function! (The latter is slow.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fixing the Seed Fixes the Random Sequence\n",
    "  * Debugging programs with random numbers is difficult because\n",
    "the numbers produced vary each time we run the program.\n",
    "  * For debugging it is important that a new run reproduces the\n",
    "sequence of random numbers in the last run.\n",
    "  * This is possible by fixing the seed of the random module:  \n",
    "`random.seed(121) # int argument`\n",
    "  * Or the numpy module:\n",
    "`np.random.seed(121) # int argument`\n",
    "  * The value of the seed determines the random sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "None\n",
      "[ 0.4359949   0.02592623  0.54966248]\n",
      "None\n",
      "[ 0.4359949   0.02592623  0.54966248]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np   \n",
    "print(np.random.seed(2))\n",
    "print(np.random.random(3))\n",
    "\n",
    "print(np.random.seed(2))\n",
    "print(np.random.random(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distributions of Random Numbers\n",
    "  * Sometimes we want uniformly distributed random numbers,\n",
    "and sometimes not.\n",
    "  * Example: it is more likely have normal (mean) blood pressure\n",
    "than large deviations from the mean.\n",
    "  * We can use the Gaussian (or \"normal\") distribution to get random\n",
    "numbers clustered around a mean value `m` with a std deviation `s`:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "43.20908995946115\n"
     ]
    }
   ],
   "source": [
    "m, s = 50, 15\n",
    "import random   \n",
    "r = random.normalvariate(m, s)\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  * Vectorized drawing of $N$ Gaussian/normal numbers:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  7.,  13.,  24.,  11.,  28.,   9.,   4.,   2.,   1.,   1.]),\n",
       " array([  18.76022096,   28.0472378 ,   37.33425464,   46.62127148,\n",
       "          55.90828832,   65.19530516,   74.482322  ,   83.76933884,\n",
       "          93.05635568,  102.34337252,  111.63038936]),\n",
       " <a list of 1 Patch objects>)"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD4dJREFUeJzt3XusZWV5x/HvTwaqiGUktsMUSCBBDGor2KhYajy1qGPT\nAKaplWo7VdL4h1VCmlawF07bNNFaqU0a/KMCmWChGC+TIaF0RstJIQ1e2gGHm6MGWlDmDPUKtY2D\nPP1jr9Hj5Fz33mvvc975fpKdWbe93uchO79ZvHutPakqJElteMa0C5AkjY+hLkkNMdQlqSGGuiQ1\nxFCXpIYY6pLUkGVDPckzk3w2yd1J7k0y220/KcmeJPuT7E6yeSLVSpKWlZXuU09yfFV9L8km4E7g\nMuDXgP+uqr9K8h7guVV1Rf/lSpKWs+L0S1V9r1s8DjgWKOBCYEe3fQdwcS/VSZLWZMVQT/KMJHcD\n88DuqvocsKWq5rtD5oEtPdYoSVql1VypP11V5wCnAq9I8uIj9heDq3dJ0pRtWu2BVfWdJLcDrwfm\nk5xcVQeSbAUOHnl8EoNekoZQVRn2vcuGepLnAU9V1beTPAt4LfA+YBewHXh/9+fOcRe23iWZrarZ\nadfRF/vb2Frur+XeYPQL4pWu1LcCO5Icw2Cq5uaqujXJXcDHklwKPAy8aZQiJEnjsWyoV9U+4KWL\nbP8mcEFfRUmShuMTpcObm3YBPZubdgE9m5t2AT2bm3YBPZqbdgHr2YoPHw194qRanlOXpD6Mmp1e\nqUtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1JBV/6CXtF4lhAl+lqs4NKmxpLUy1NWCE4Dv\nAk/1PE4Y/NbRmT2PIw3NUFcrnqziOX0OkHAmcFufY0ijck5dkhpiqEtSQwx1SWqIoS5JDTHUJakh\nhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktSQZUM9yWlJ\nbk9yX5J7k7y72z6b5NEke7vXtsmUK0lazkr/8tEh4PKqujvJCcC/J9kDFHB1VV3de4WSpFVbNtSr\n6gBwoFt+MskDwCnd7vRcmyRpjVY9p57kdOBc4K5u07uS3JPk2iSbe6hNkrRGqwr1burl48BlVfUk\n8GHgDOAc4DHgg71VKElatZXm1ElyLPAJ4KNVtROgqg4u2P8R4JYl3ju7YHWuquZGKVaSWpNkBpgZ\n2/mqarnBAuwAvlFVly/YvrWqHuuWLwdeVlW/ecR7q6qcd1fvEp4DfL2K5/Q8zpnAbVWc2ec4OrqN\nmp0rXamfD7wV+GKSvd229wKXJDmHwV0wDwHvGLYASdL4rHT3y50sPu/+T/2UI0kahU+USlJDDHVJ\nasiKd7+oTQnnA8+c0HBzVfxgQmNJRzVD/ej1UQbPGPxvz+P8EnAC8L2ex5GEoX60e0sVD/U5QGKY\nS5PknLokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12S\nGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakh\ny4Z6ktOS3J7kviT3Jnl3t/2kJHuS7E+yO8nmyZQrSVrOSlfqh4DLq+pFwHnAO5OcDVwB7Kmqs4DP\ndOuSpClbNtSr6kBV3d0tPwk8AJwCXAjs6A7bAVzcZ5GSpNVZ9Zx6ktOBc4HPAluqar7bNQ9sGXtl\nkqQ127Sag5KcAHwCuKyqnkjyw31VVUlqiffNLlidq6q54UuVpPYkmQFmxnW+FUM9ybEMAv2GqtrZ\nbZ5PcnJVHUiyFTi42HuranZchUpSi7qL3bnD60muGuV8K939EuBa4P6q+tCCXbuA7d3ydmDnke+V\nJE3eSlfq5wNvBb6YZG+37UrgfcDHklwKPAy8qbcKJUmrtmyoV9WdLH01f8H4y5EkjcInSiWpIYa6\nJDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtS\nQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXE\nUJekhhjqktSQFUM9yXVJ5pPsW7BtNsmjSfZ2r239lilJWo3VXKlfDxwZ2gVcXVXndq/bxl+aJGmt\nVgz1qroD+NYiuzL+ciRJoxhlTv1dSe5Jcm2SzWOrSJI0tE1Dvu/DwJ93y38BfBC49MiDkswuWJ2r\nqrkhx5OkJiWZAWbGdb6hQr2qDh5eTvIR4JYljpsdrixJOjp0F7tzh9eTXDXK+YaafkmydcHqG4F9\nSx0rSZqcFa/Uk9wEvBp4XpJHgKuAmSTnMLgL5iHgHb1WKUlalRVDvaouWWTzdT3UIkkakU+USlJD\nDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQ4b9lUb1JOHXgS0TGOonJzCGpAkz\n1Nef3wceBQ70PM5NwBM9jyFpwgz19emvq7hr2kVI2nicU5ekhhjqktQQQ12SGmKoS1JDDHVJaoih\nLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQ/xBL03CXyYc6vH8x/V4bmlDMdTVtz9mMp+z\nP53AGNK6Z6irV1VcPe0apKOJc+qS1JAVQz3JdUnmk+xbsO2kJHuS7E+yO8nmfsuUJK3Gaq7Urwe2\nHbHtCmBPVZ0FfKZblyRN2YqhXlV3AN86YvOFwI5ueQdw8ZjrkiQNYdg59S1VNd8tzwNbxlSPJGkE\nI9/9UlWVpBbbl2R2wepcVc2NOp4ktSTJDDAzrvMNG+rzSU6uqgNJtgIHFzuoqmaHrkySjgLdxe7c\n4fUkV41yvmGnX3YB27vl7cDOUYqQJI3Ham5pvAn4N+AFSR5J8jbgfcBrk+wHXtOtS5KmbMXpl6q6\nZIldF4y5FknSiHyiVJIaYqhLUkP8QS9pbU5OuHkC4+yv4k8mMI4aY6hLqzcPXDqBcV4AvB4Mda2d\noS6tUhVPQP9X6QnnMwh1ac2cU5ekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCX\npIYY6pLUEENdkhpiqEtSQ/xBr1VK+ADwqgkM9eIJjCGpUYb66p0FfBL41wmMdd8ExpDUIEN9bb5U\nxV3TLkKSluKcuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhIz18lORh4LvAD4BD\nVfXycRQlSRrOqE+UFjBTVd8cRzGSpNGMY/olYziHJGkMRg31Aj6d5AtJfnccBUmShjfq9Mv5VfVY\nkp8C9iR5sKruGEdhkqS1GynUq+qx7s/Hk3wKeDnww1BPMrvg8LmqmhtlPElqTZIZYGZc5xs61JMc\nDxxTVU8keTbwOuDPFh5TVbOjlSdJbesuducOrye5apTzjXKlvgX4VJLD5/mHqto9SjGSpNEMHepV\n9RBwzhhrkSSNyCdKJakhhrokNcR/o1Ran05N+MMJjPPlKj41gXE0IYa6tP58DbgZeF7P45wFvBIM\n9ZYY6tI6U8XD0P9VesIbgd/uexxNlnPqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCX\npIYY6pLUEENdkhpiqEtSQwx1SWqIP+glHd1ekvB3Exjnzir+cQLjHPUMdeno9UXg6gmM8ypgBgz1\nSTDUpaNUFV+F/q/SE57Cf894YpxTl6SGGOqS1BBDXZIaYqhLUkP8olTSJJyYcOa0ixijx6v4zrSL\nWIyhLqlv3wVeBtw27ULG5KeBK4Brpl3IYgx1Sb2q4kbgxmnXMS7J+gzzw5xTl6SGDB3qSbYleTDJ\nl5O8Z5xFSZKGM1SoJzmGwZNo24AXApckOXucha13SWamXUOf7G9ja7m/lnsbh2Gv1F8OfKWqHq6q\nQwx+0+Gi8ZW1IcxMu4CezUy7gJ7NTLuAns1Mu4AezUy7gPVs2FA/BXhkwfqj3TZJ0hQNG+o11io2\nhq/C+rwvVdJEPQI8Pu0ilpKqtedzkvOA2ara1q1fCTxdVe9fcMzRGPySNLKqyrDvHTbUNwFfAn4Z\n+DrwOeCSqnpg2EIkSaMb6uGjqnoqye8B/wwcA1xroEvS9A11pS5JWp96eaK0tQeTkpyW5PYk9yW5\nN8m7u+0nJdmTZH+S3Uk2T7vWYSU5JsneJLd06y31tjnJx5M8kOT+JK9orL/Lu8/lviQ3JvmJjdxf\nkuuSzCfZt2Dbkv0kubLLmgeTvG46Va/eEv19oPt83pPkk0lOXLBvTf2NPdQbfTDpEHB5Vb0IOA94\nZ9fTFcCeqjoL+Ey3vlFdBtzPj+5saqm3vwVuraqzgZ8DHqSR/pKcArwL+Pmq+lkG06FvZmP3dz2D\n/Fho0X6SvBD4DQZZsw24Jsl6//mTxfrbDbyoql4C7AeuhOH666P55h5MqqoDVXV3t/wk8ACD+/Iv\nBHZ0h+0ALp5OhaNJcirwK8BHgMPfurfS24nAq6rqOhh8H1RV36GR/jqbgOO7GxiOZ3Dzwobtr6ru\nAL51xOal+rkIuKmqDlXVw8BXGGTQurVYf1W1p6qe7lY/C5zaLa+5vz5CvekHk5KcDpzL4D/8lqqa\n73bNA1umVNao/gb4A+DpBdta6e0M4PEk1yf5jyR/n+TZNNJfVX0N+CDwXwzC/NtVtYdG+ltgqX5+\nhkHGHNZC3rwduLVbXnN/fYR6s9+8JjkB+ARwWVU9sXBfDb5x3nC9J/lV4GBV7eVHV+k/ZqP21tkE\nvBS4pqpeCvwPR0xFbOT+kjyXwVXs6QwC4IQkb114zEbubzGr6GfD9prkj4DvV9VyP1W8bH99hPrX\ngNMWrJ/Gj/9NsyElOZZBoN9QVTu7zfNJTu72bwUOTqu+EfwCcGGSh4CbgNckuYE2eoPBZ+/Rqvp8\nt/5xBiF/oJH+LgAeqqpvVNVTwCeBV9JOf4ct9Xk8Mm9O7bZtOEl+h8E06FsWbF5zf32E+heA5yc5\nPclxDCb5d/UwzsQkCXAtcH9VfWjBrl3A9m55O7DzyPeud1X13qo6rarOYPAF279U1W/RQG8w+D4E\neCTJWd2mC4D7gFtooD/gP4Hzkjyr+5xewOAL71b6O2ypz+Mu4M1JjktyBvB8Bg9DbihJtjGYAr2o\nqv5vwa6191dVY38Bb2DwxOlXgCv7GGOSL+AXGcw33w3s7V7bgJOATzP4tno3sHnatY7Y56uBXd1y\nM70BLwE+D9zD4Er2xMb6m2Xw5f0+Bl8iHruR+2Pwf4xfB77P4Pu5ty3XD/DeLmseBF4/7fqH6O/t\nwJcZ/AV9OF+uGbY/Hz6SpIas9/s5JUlrYKhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktSQ\n/wdH86sGWHef1QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x109873240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "(m, s, N) = (50, 15, 100)\n",
    "import numpy as np   \n",
    "samples = np.random.normal(m, s, N) \n",
    "# print(samples)\n",
    "plt.hist(samples, histtype='step')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probability via Monte Carlo Simulation\n",
    "  * What is the probability that a certain event A happens?\n",
    "    1. Simulate $N$ events\n",
    "    2. Count how many times $M$ the event $A$ happens\n",
    "    3. The probability of the event $A$ is then $M/N$ (as $N \\rightarrow \\infty$).\n",
    "  * Example: what is the probability of getting 6 on two or more\n",
    "dice if we throw 4 dice?  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    " N = 1000 # number of experiments   \n",
    " M = 0 # number of successful events   \n",
    "\n",
    "def count_six():\n",
    "     six = 0\n",
    "     dice_1 = random.randint(1, 6)   \n",
    "     dice_2 = random.randint(1, 6)   \n",
    "     dice_3 = random.randint(1, 6)   \n",
    "     dice_4 = random.randint(1, 6)   \n",
    "\n",
    "     if dice_1 == 6: six += 1   \n",
    "     if dice_2 == 6: six += 1   \n",
    "     if dice_3 == 6: six += 1   \n",
    "     if dice_4 == 6: six += 1   \n",
    "    \n",
    "     return six"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 149,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count_six()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "probability of getting 6 on two or more dice out of 4 dice is:\n",
      "0.138000 based on 1000 trials\n"
     ]
    }
   ],
   "source": [
    "    \n",
    "for i in range(N):   \n",
    "     # successful event?   \n",
    "     if count_six() >= 2:   \n",
    "         M += 1 \n",
    "        \n",
    "print('probability of getting 6 on two or more dice out of 4 dice is:')\n",
    "print('%f based on %d trials' % (float(M)/N, N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Walks\n",
    "\n",
    "### Random Walk in One Spatial Dimension\n",
    "  * A particle moves to the left and right with equal probability\n",
    "  * $n$ particles start at $x = 0$ at time $t = 0$ - how do the particles\n",
    "get distributed over time?\n",
    "  * This is called a random walk and constitutes a simple model for\n",
    "molecular motion in many situations:\n",
    "    - heat transport\n",
    "    - quantum mechanics\n",
    "    - polymer chains\n",
    "    - population genetics\n",
    "    - pricing of financial instruments\n",
    "    - ...\n",
    "  * We'll make a program for simulating random walk.\n",
    "  \n",
    "### Random Walk as a Difference Equation\n",
    "  * Let $x_n$ be the position of one particle at time $n$.\n",
    "  * Updating rule: $x_{n+1} = x_n + s$  \n",
    "where $s = 1$ or $s = -1$, both with probability $1/2$.\n",
    "  * For $n_p$ particles, we need $n_p$ such difference equations.\n",
    "  * Would like to calculate the statistics (mean position, \"width\" of the cluster of particles, how particles are distributed throughout space).\n",
    "  * Random walk in 2-D:\n",
    "    - $x_{n+1} = x_n + r_x$\n",
    "    - $y_{n+1} = y_n + r_y$\n",
    "    - $r_x = r \\cos\\theta$, $r_y = r \\sin\\theta$, with random $\\theta$\n",
    "    \n",
    "### How Far After $N$ Steps (in 2-D)?\n",
    "  * Start at $(x, y) = (0, 0)$.\n",
    "  * After $N$ steps:  \n",
    "$R^2 = (\\Delta x_1 + \\cdots + \\Delta x_N)^2 + (\\Delta y_1 + \\cdots + \\Delta y_N)^2$\n",
    "  * Expand:  \n",
    "$R^2 = \\Delta x_1^2 + \\cdots + \\Delta x_N^2 + 2 \\Delta x_1 \\Delta x_2 + \\cdots  \n",
    "+ \\Delta y_1^2 + \\cdots + \\Delta y_N^2 + 2 \\Delta y_1 \\Delta y_2 + \\cdots$\n",
    "  * Cross terms average to zero:  \n",
    "$R^2 \\simeq \\Delta x_1^2 + \\cdots + \\Delta x_N^2 + \\Delta y_1^2 + \\cdots + \\Delta y_N^2$\n",
    "  * $\\Rightarrow R^2 \\simeq N \\langle r^2 \\rangle$\n",
    "  * $\\Rightarrow R \\simeq \\sqrt{N} r_{RMS}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spontaneous (e.g. Radioactive) Decay\n",
    "  * \"Spontaneous\" process $\\rightarrow$ no external stimulus.\n",
    "  * Time of decay is random, independent of:\n",
    "    - How long it existed.\n",
    "    - How many other are around.\n",
    "  * Describes nuclear decay, electronic relaxation, etc.\n",
    "  * Theory: probability of decay per $\\Delta t$ per particle is:   \n",
    "$P(t)=\\lambda$ (constant)    \n",
    "$\\Rightarrow$ $N(t)$, $dN/dt$ decrease with time.  \n",
    "\n",
    "### Model: Continuous Decay\n",
    "\n",
    "  * If $N\\rightarrow \\infty$, $\\Delta t \\rightarrow 0$:  \n",
    "$\\frac{\\Delta N(t)}{\\Delta t} \\rightarrow \\frac{dN(t)}{dt} = -\\lambda N(t)$\n",
    "  * Can integrate differential equation:  \n",
    "$N(t) = N(0) e^{-\\lambda t} = N(0) e^{-t/\\tau}$  \n",
    "with $\\lambda = 1/\\tau$\n",
    "  * $\\frac{dN}{dt}(t) = -\\lambda N(0) e^{-\\lambda t} = \\frac{dN}{dt}(0) e^{-t/\\tau}$\n",
    "  * But this is an approximation!  \n",
    "Nature can have small $N$, stochastic behavior.\n",
    "\n",
    "### Method: Decay Simulation\n",
    "  * Algorithm:  \n",
    "Loop through remaining nuclei   \n",
    "$r_i < \\lambda$? $\\Rightarrow$ decays   \n",
    "$t = t + \\Delta t$   \n",
    "Repeat loop  \n",
    "  * $r_i$ is random number in $[0, 1]$.\n",
    "  * Higher $\\lambda\\Rightarrow$ more frequent decays."
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}