{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 8: Transformations\n",
    " \n",
    "This Jupyter notebook is the Python equivalent of the R code in section 8.8 R, pp. 373 - 375, [Introduction to Probability, 2nd Edition](https://www.crcpress.com/Introduction-to-Probability-Second-Edition/Blitzstein-Hwang/p/book/9781138369917), Blitzstein & Hwang.\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Beta and Gamma distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import beta, gamma\n",
    "\n",
    "# to learn more about scipy.stats.beta, un-comment ouf the following line\n",
    "#print(beta.__doc__)\n",
    "\n",
    "# to learn more about scipy.stats.gamma, un-comment ouf the following line\n",
    "#rint(gamma.__doc__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Beta and Gamma distributions are implemented in [`scipy.stats.beta`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.beta.html) and [`scipy.stats.gamma`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gamma.html), respectively.\n",
    "\n",
    "* To evaluate the $Beta(a, b)$ PDF or CDF at $x$, we use `beta.pdf(x, a, b)` and `beta.cdf(x, a, b)`. To generate $n$ realizations from the $Beta(a, b)$ distribution, we use `beta.rvs(a, b, size=n)`. \n",
    "* To evaluate the $Gamma(a, \\lambda)$ PDF or CDF at $x$, we use`gamma.pdf(x, a, scale=1/lambd)` or `gamma.cdf(x, a, scale=1/lambd)`. To generate $n$ realizations from the $Gamma(a, \\lambda)$ distribution, we use `gamma.rvs(a, scale=1/lambd, size=n)`. \n",
    "   * The $\\lambda$ parameter in $Gamma(a, \\lambda)$ corresponds in `gamma` to using `scale` = $\\frac{1}{\\lambda}$ and default value of `loc` = 0.\n",
    "\n",
    "For example, we can check that the $Gamma(3, 2)$ distribution has mean $\\frac{3}{2}$ and variance $\\frac{3}{4}$. To do this, we generate a large number of $Gamma(3, 2)$ random variables using `gamma.rvs`, then compute their mean and var using their corresponding methods in `numpy.array` (you could of course also use `numpy.mean` and `numpy.var`, passing in the array of r.v.):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean of Gamma(3, 2) = 1.5012880719327166\n",
      "variance of Gamma(3, 2) = 0.7512319256875035\n"
     ]
    }
   ],
   "source": [
    "# seed the random number generator\n",
    "np.random.seed(317811)\n",
    "\n",
    "alpha = 3.0\n",
    "lambd = 2.0\n",
    "\n",
    "y = gamma.rvs(alpha, scale=1/lambd, size=10**5)\n",
    "\n",
    "mean = y.mean()\n",
    "#mean = np.mean(y)\n",
    "print('mean of Gamma(3, 2) = {}'.format(mean))\n",
    "\n",
    "var = y.var()\n",
    "#var = np.var(y)\n",
    "print('variance of Gamma(3, 2) = {}'.format(var))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOTE**: [`lambda`](https://docs.python.org/3/tutorial/controlflow.html#lambda-expressions) in Python is a reserved keyword for declaring small, anonymous functions, and so we cannot used the name `lambda` for variable or functions.\n",
    "\n",
    "Try changing the `numpy.random.seed` input value in the code block above, and hit SHIFT+ENTER to re-execute the code block. Did you get values that are close to 1.5 and 0.75, respectively?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convolution of Uniforms\n",
    "\n",
    "Using [`scipy.stats.uniform`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.uniform.html) and [`matplotlib.pyplot.hist`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.hist.html), we can quickly verify that for $X, Y \\stackrel{i.i.d.}{\\sim} Unif(0, 1)$, the distribution of $T = X + Y$ is triangular in shape:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEbCAYAAAAS4RmTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAG+xJREFUeJzt3X+8ZXVd7/HXWxAUUAFBowEcsAkD\nk6RJSe2KYgKSDFbe8GE5eimysLIfCmo3uP4ovI8S45oWCQVkAhIq/spGkOwX4PBDfooMMMEAwRgI\nggmBn/vH+h7cczjnzD7nrH32HOb1fDz246z1Xd+11mev2XM++/vjrJWqQpKkPjxh3AFIkh4/TCqS\npN6YVCRJvTGpSJJ6Y1KRJPXGpCJJ6o1JRZLUG5OKJKk3JhUtCkmuSXLAuOMYlyR7Jbk8ybeT/Oa4\n45GmY1LR2CVZm+QVk8remOSfJ9arap+qunC2x3kceTtwYVU9papOGtyQZGWS+9vru0keGVj/VpKt\n53rSJM9O8kCSXQbKXp/k9iS7zeO4O7QE+dJJ5WckOTdJ5npsjZdJRRpCki3HHMKzgGum2lBVp1XV\ndlW1HfCHwGcn1qtq+6p6cK4nraobgc8CbwVI8pPAh4DDq+rW6fZLcnyS42c47j3AycBvD+zzv4G9\ngV8s7x+1aJlUtCgMtkKSHJPktvZN9/okByY5A9gd+Ez7hv72VvdHklzYvrFfk+SwgWPuN9Cl9Ikk\nZyV576RzHpPkSuCBJFsmOTbJjW2fa5O8ZlL9tyW5sn27PyXJM5N8odX/UpIdZniPU8aa5ALgZcCH\n2nv74Rku1Y8BX5vTRZ7e+4FfTfJc4FzgzVV1SQ/H/QBwUGsNvRY4Cnh1VX2nh2NrTEwqWlSS7AW8\nBfiJqnoKcBCwtqp+CbiF7pfSdlX1f5M8EfgM8A/AM4DfAD7Wxie2Aj4J/DWwI/Bx4DWPOSG8DjgU\n2L6qHgZuBH4KeBrwf4C/GewaAn4O+Gngh4FXA18A3gnsRPf/bcrxkJliraqXA/8EvKW9t2/McIl+\nDLhiuo1JPtuS1lSvz061T1VdBlwCXAx8pKrOmuH8Q6uq2+iu+58BHwZWVNXtfRxb4zPuJr004VNJ\nHh5Y3wq4bIp6jwBbA3snWV9Va2c45v7AdsAJVfU94IL2i/N1wAV0n/+TWlfLuUmm+vZ90mA3T1V9\nYmDbWUneAbwA+HQr+39VdSdAkn8C7qqqy9v6J4ED5xDr8TO8x0cleSqwlBmSSlX9zDDHmnTcJ9Bd\n9+/RtVr69AHgKuAXWvLSImdLRZuKw1v///ZVtT3w61NVqqo1dP37xwN3JTkzyQ9Oc8wfBG5tv6Qn\n/DuwpG27bVLf/VRjBBuUJXlDkismvt0Dz6VrhUy4c2D5v6ZY324OsQ5rX+DbwM2z2GcYfwJsD9wA\nvH66SoOtIOBY4NiNtYLovjw8SNetpscBk4oWnar626p6Cd3gdfH9b8+TB3dvB3Zr37Qn7A7cBtwB\nLJk0y2iq2UyPHjPJs4C/pOt+e3pLflcDfcxUminWYf0YcOVMg9xtfOf+aV5fmKL+r9J1Cx5Od53f\nNt3MrKr6mYEvBSfQtbomvihM10LaF7i6dS3qccCkokWljYe8vE2T/S7dt/9H2uY7gT0Hql8MPAC8\nPckT0/2dy6uBM4F/a/u9pQ3Ar6DrxprJtnRJZn2L5U10LZU+zBTrsGYcTwGoqkMGZoZNfh0yWLdN\njPhDunGqO4Fz6FoWK2YR07xj1uJiUtFiszXdt+BvAv9BN6j9zrbtj4Dfb90tv1dVDwGHAYe0+h8G\n3lBVX2/bfhY4EvgW8It0U2ennX5bVdfSdQX9G10C+1HgX/p4UzPFOovD7EtPv6CTPIcuof1SVV3V\nYnyEbgzkmD7O0fQWszYNcTq41ElyMfDnVfVX445FWqxsqWizleSlSX6gdX+tBJ4H/P2445IWM6cU\na3O2F3A23YysG4Gfr6o7xhuStLiNrKWS5NQkdyW5eqBsxySrktzQfu7QypPkpCRr2l8j7zewz8pW\n/4b2bXKi/MeTXNX2OWm6GSnSdKrq5Kp6ZlVtW1XPq6rPjTsmabEbZffXXwMHTyo7Fji/qpYB57d1\n6AYnl7XXUcBHoEtCwHHAC+lm5hw3cJuLj7S6E/tNPpckaYGNrPurqr6SZOmk4hXAAW35NOBCupkk\nK4DT2/z6i5Js3259cQCwqqruBkiyCjg4yYXAU6vq31r56XTz6B8zz36ynXbaqZYunRyWJGk6l156\n6Teraudh6i70mMozJ/qsq+qOJM9o5UvY8C+X17WymcrXTVE+pSRH0bVq2H333Vm9evU834YkbT6S\n/PuwdTeV2V9TjYfUHMqn1PrOl1fV8p13HirZSpLmYKGTyp0Td3RtP+9q5evY8BYZu9LdtmKm8l2n\nKJckjdFCJ5XzgIkZXCv5/p1dzwPe0GaB7Q/c27rJvgi8Mt1T4nYAXgl8sW37dpL926yvNwwcS5I0\nJiMbU0nycbqB9p2SrKObxXUCcHaSI+meffHaVv3zwKuANcB3gDcBVNXdSd4DfLXVe/fEoD3wa3Qz\nzJ5MN0C/0UF6SdJobXa3aVm+fHk5UC9Jw0tyaVUtH6bupjJQL0l6HDCpSJJ6Y1KRJPXGpCJJ6o13\nKZZ6tPTYud+Tcu0Jh/YYiTQetlQkSb2xpSJtImzl6PHAlookqTcmFUlSb0wqkqTemFQkSb0xqUiS\nemNSkST1xqQiSeqNSUWS1BuTiiSpNyYVSVJvTCqSpN6YVCRJvTGpSJJ6Y1KRJPXGpCJJ6o1JRZLU\nG5OKJKk3JhVJUm9MKpKk3phUJEm92XLcAUibmqXHfm7cIUiLli0VSVJvTCqSpN7Y/SU9Dsyny27t\nCYf2GIk2d7ZUJEm9MalIknpjUpEk9WYsSSXJbye5JsnVST6e5ElJ9khycZIbkpyVZKtWd+u2vqZt\nXzpwnHe08uuTHDSO9yJJ+r4FTypJlgC/CSyvqucCWwBHAO8HTqyqZcA9wJFtlyOBe6rqh4ATWz2S\n7N322wc4GPhwki0W8r1IkjY0ru6vLYEnJ9kS2Aa4A3g5cE7bfhpweFte0dZp2w9MklZ+ZlU9WFU3\nA2uAFyxQ/JKkKSx4Uqmq24A/Bm6hSyb3ApcC36qqh1u1dcCStrwEuLXt+3Cr//TB8in22UCSo5Ks\nTrJ6/fr1/b4hSdKjxtH9tQNdK2MP4AeBbYFDpqhaE7tMs2268scWVp1cVcuravnOO+88+6AlSUMZ\nR/fXK4Cbq2p9Vf03cC7wImD71h0GsCtwe1teB+wG0LY/Dbh7sHyKfSRJYzCOpHILsH+SbdrYyIHA\ntcCXgZ9vdVYCn27L57V12vYLqqpa+RFtdtgewDLgkgV6D5KkKSz4bVqq6uIk5wCXAQ8DlwMnA58D\nzkzy3lZ2StvlFOCMJGvoWihHtONck+RsuoT0MHB0VT2yoG9GkrSBsdz7q6qOA46bVHwTU8zeqqrv\nAq+d5jjvA97Xe4CSpDnxL+olSb0xqUiSemNSkST1xqQiSeqNSUWS1BuTiiSpNz5OWI8783m0rqT5\nMalIm7n5JmGfca9Bdn9JknpjUpEk9cakIknqjUlFktQbk4okqTcmFUlSb0wqkqTemFQkSb0xqUiS\nemNSkST1xqQiSeqNSUWS1BuTiiSpNyYVSVJvTCqSpN6YVCRJvTGpSJJ6Y1KRJPXGxwlrk+Rz5qXF\nyZaKJKk3JhVJUm/s/pI0L/Ppqlx7wqE9RqJNgS0VSVJvTCqSpN6YVCRJvTGpSJJ6Y1KRJPVmqKSS\n5Ll9njTJ9knOSfL1JNcl+ckkOyZZleSG9nOHVjdJTkqyJsmVSfYbOM7KVv+GJCv7jFGSNHvDtlT+\nPMklSX49yfY9nPdPgb+vqucA+wLXAccC51fVMuD8tg5wCLCsvY4CPgKQZEfgOOCFwAuA4yYSkSRp\nPIZKKlX1EuD1wG7A6iR/m+Sn53LCJE8F/gdwSjv2Q1X1LWAFcFqrdhpweFteAZxenYuA7ZPsAhwE\nrKqqu6vqHmAVcPBcYpIk9WPoMZWqugH4feAY4KXASa376mdnec49gfXAXyW5PMlHk2wLPLOq7mjn\nugN4Rqu/BLh1YP91rWy68sdIclSS1UlWr1+/fpbhSpKGNeyYyvOSnEjXTfVy4NVV9SNt+cRZnnNL\nYD/gI1X1fOABvt/VNeXppyirGcofW1h1clUtr6rlO++88yzDlSQNa9iWyoeAy4B9q+roqroMoKpu\np2u9zMY6YF1VXdzWz6FLMne2bi3az7sG6u82sP+uwO0zlEuSxmTYpPIq4G+r6r8AkjwhyTYAVXXG\nbE5YVf8B3Jpkr1Z0IHAtcB4wMYNrJfDptnwe8IY2C2x/4N7WPfZF4JVJdmgD9K9sZZKkMRn2hpJf\nAl4B3N/WtwH+AXjRHM/7G8DHkmwF3AS8iS7BnZ3kSOAW4LWt7ufpktoa4DutLlV1d5L3AF9t9d5d\nVXfPMR5JUg+GTSpPqqqJhEJV3T/RUpmLqroCWD7FpgOnqFvA0dMc51Tg1LnGIUnq17DdXw9M+qPD\nHwf+azQhSZIWq2FbKm8FPpFkYiB8F+AXRhOSJGmxGiqpVNVXkzwH2ItuKu/Xq+q/RxqZFj2fMy9t\nfmbz5MefAJa2fZ6fhKo6fSRRSZIWpaGSSpIzgGcDVwCPtOICTCqSpEcN21JZDuzdZmJJkjSlYZPK\n1cAPAHeMMBZJm5n5jLutPeHQHiNRX4ZNKjsB1ya5BHhworCqDhtJVJKkRWnYpHL8KIOQJD0+DDul\n+B+TPAtYVlVfan9Nv8VoQ5MkLTbD3vr+V+juJvwXrWgJ8KlRBSVJWpyGvU3L0cCLgfvg0Qd2PWPG\nPSRJm51hk8qDVfXQxEqSLZnmgViSpM3XsEnlH5O8E3hyezb9J4DPjC4sSdJiNGxSOZbuufJXAb9K\n94yT2T7xUZL0ODfs7K/vAX/ZXpIkTWnYe3/dzBRjKFW1Z+8RSZIWrdnc+2vCk+ge9btj/+FIkhaz\nocZUquo/B163VdUHgZePODZJ0iIzbPfXfgOrT6BruTxlJBFJkhatYbu//mRg+WFgLfA/e49GkrSo\nDTv762WjDkSStPgN2/31OzNtr6oP9BOOJGkxm83sr58Azmvrrwa+Atw6iqAkSYvTbB7StV9VfRsg\nyfHAJ6rql0cVmCRp8Rk2qewOPDSw/hCwtPdotMmZz+NeJW1+hk0qZwCXJPkk3V/WvwY4fWRRSdJG\n+Hz7TdOws7/el+QLwE+1ojdV1eWjC0uStBgNe5digG2A+6rqT4F1SfYYUUySpEVq2McJHwccA7yj\nFT0R+JtRBSVJWpyGbam8BjgMeACgqm7H27RIkiYZNqk8VFVFu/19km1HF5IkabEaNqmcneQvgO2T\n/ArwJXxglyRpkmFnf/1xezb9fcBewB9U1aqRRiZJWnQ22lJJskWSL1XVqqp6W1X9Xh8JpR338iSf\nbet7JLk4yQ1JzkqyVSvfuq2vaduXDhzjHa38+iQHzTcmSdL8bDSpVNUjwHeSPK3nc/8WcN3A+vuB\nE6tqGXAPcGQrPxK4p6p+CDix1SPJ3sARwD7AwcCHk2zRc4ySpFkYdkzlu8BVSU5JctLEa64nTbIr\ncCjw0bYeuidJntOqnAYc3pZXtHXa9gNb/RXAmVX1YFXdDKwBXjDXmCRJ8zfsbVo+1159+SDwdr4/\nLfnpwLeq6uG2vg5Y0paX0O6GXFUPJ7m31V8CXDRwzMF9NpDkKOAogN13372/dyFJ2sCMSSXJ7lV1\nS1WdNlO92UjyM8BdVXVpkgMmiqeoWhvZNtM+GxZWnQycDLB8+fIp60iS5m9j3V+fmlhI8nc9nfPF\nwGFJ1gJn0nV7fZBuuvJEktsVuL0trwN2azFsCTwNuHuwfIp9JEljsLGkMtga2LOPE1bVO6pq16pa\nSjfQfkFVvR74MvDzrdpK4NNt+by2Ttt+QftDzPOAI9rssD2AZcAlfcQoSZqbjY2p1DTLo3AMcGaS\n9wKXA6e08lOAM5KsoWuhHAFQVdckORu4FngYOLrNVJMkjcnGksq+Se6ja7E8uS3T1quqnjqfk1fV\nhcCFbfkmppi9VVXfBV47zf7vA943nxgkSf2ZMalUlX/3IUka2myepyJJ0oxMKpKk3gz7x4+S9Lgx\nn+fbg8+4n4ktFUlSb2ypPM7N9xuZJM2GLRVJUm9MKpKk3phUJEm9MalIknpjUpEk9cakIknqjUlF\nktQbk4okqTcmFUlSb0wqkqTemFQkSb0xqUiSemNSkST1xqQiSeqNSUWS1BuTiiSpNyYVSVJvTCqS\npN74OOFFwEcCS1osbKlIknpjUpEk9cbuL0mapfl0Sa894dAeI9n02FKRJPXGpCJJ6o1JRZLUG5OK\nJKk3JhVJUm9MKpKk3ix4UkmyW5IvJ7kuyTVJfquV75hkVZIb2s8dWnmSnJRkTZIrk+w3cKyVrf4N\nSVYu9HuRJG1oHC2Vh4HfraofAfYHjk6yN3AscH5VLQPOb+sAhwDL2uso4CPQJSHgOOCFwAuA4yYS\nkSRpPBY8qVTVHVV1WVv+NnAdsARYAZzWqp0GHN6WVwCnV+ciYPskuwAHAauq6u6qugdYBRy8gG9F\nkjTJWMdUkiwFng9cDDyzqu6ALvEAz2jVlgC3Duy2rpVNVz7VeY5KsjrJ6vXr1/f5FiRJA8aWVJJs\nB/wd8Naqum+mqlOU1Qzljy2sOrmqllfV8p133nn2wUqShjKWpJLkiXQJ5WNVdW4rvrN1a9F+3tXK\n1wG7Dey+K3D7DOWSpDEZx+yvAKcA11XVBwY2nQdMzOBaCXx6oPwNbRbY/sC9rXvsi8Ark+zQBuhf\n2cokSWMyjrsUvxj4JeCqJFe0sncCJwBnJzkSuAV4bdv2eeBVwBrgO8CbAKrq7iTvAb7a6r27qu5e\nmLcgSZrKgieVqvpnph4PAThwivoFHD3NsU4FTu0vOknSfPgX9ZKk3phUJEm9MalIknpjUpEk9cak\nIknqzTimFEvSZmvpsZ+b875rTzi0x0hGw6SyQObzQZKkxcLuL0lSb0wqkqTemFQkSb0xqUiSemNS\nkST1xqQiSeqNSUWS1BuTiiSpNyYVSVJvTCqSpN6YVCRJvTGpSJJ6Y1KRJPXGpCJJ6o1JRZLUG5OK\nJKk3JhVJUm9MKpKk3phUJEm98Rn1krRILD32c3Ped+0Jh/YYyfRMKrMwn39QSdoc2P0lSeqNSUWS\n1BuTiiSpNyYVSVJvTCqSpN6YVCRJvTGpSJJ6s+iTSpKDk1yfZE2SY8cdjyRtzhZ1UkmyBfBnwCHA\n3sDrkuw93qgkafO1qJMK8AJgTVXdVFUPAWcCK8YckyRtthb7bVqWALcOrK8DXji5UpKjgKPa6v1J\nrp/j+XYCvjnHfUfJuGbHuGbHuGZnk4wr7wfmHtuzhq242JNKpiirxxRUnQycPO+TJauravl8j9M3\n45od45od45qdTTUuWJjYFnv31zpgt4H1XYHbxxSLJG32FntS+SqwLMkeSbYCjgDOG3NMkrTZWtTd\nX1X1cJK3AF8EtgBOraprRnjKeXehjYhxzY5xzY5xzc6mGhcsQGypeswQhCRJc7LYu78kSZsQk4ok\nqTcmFTZ+q5ckWyc5q22/OMnSgW3vaOXXJzlogeP6nSTXJrkyyflJnjWw7ZEkV7RXr5MXhojrjUnW\nD5z/lwe2rUxyQ3utXOC4ThyI6RtJvjWwbZTX69QkdyW5eprtSXJSi/vKJPsNbBvl9dpYXK9v8VyZ\n5F+T7DuwbW2Sq9r1Wr3AcR2Q5N6Bf68/GNg2sts2DRHX2wZiurp9pnZs20Z5vXZL8uUk1yW5Jslv\nTVFn4T5jVbVZv+gG+G8E9gS2Ar4G7D2pzq8Df96WjwDOast7t/pbA3u042yxgHG9DNimLf/aRFxt\n/f4xXq83Ah+aYt8dgZvazx3a8g4LFdek+r9BN7FjpNerHft/APsBV0+z/VXAF+j+7mp/4OJRX68h\n43rRxPnoboV08cC2tcBOY7peBwCfne9noO+4JtV9NXDBAl2vXYD92vJTgG9M8X9ywT5jtlSGu9XL\nCuC0tnwOcGCStPIzq+rBqroZWNOOtyBxVdWXq+o7bfUiur/TGbX53BrnIGBVVd1dVfcAq4CDxxTX\n64CP93TuGVXVV4C7Z6iyAji9OhcB2yfZhdFer43GVVX/2s4LC/f5GuZ6TWekt22aZVwL+fm6o6ou\na8vfBq6ju9vIoAX7jJlUpr7Vy+R/kEfrVNXDwL3A04fcd5RxDTqS7pvIhCclWZ3koiSH9xTTbOL6\nudbMPifJxB+obhLXq3UT7gFcMFA8qus1jOliH+X1mq3Jn68C/iHJpelug7TQfjLJ15J8Ick+rWyT\nuF5JtqH7xfx3A8ULcr3Sdc0/H7h40qYF+4wt6r9T6ckwt3qZrs5Qt4mZo6GPneQXgeXASweKd6+q\n25PsCVyQ5KqqunGB4voM8PGqejDJm+laeS8fct9RxjXhCOCcqnpkoGxU12sY4/h8DS3Jy+iSyksG\nil/crtczgFVJvt6+yS+Ey4BnVdX9SV4FfApYxiZyvei6vv6lqgZbNSO/Xkm2o0tkb62q+yZvnmKX\nkXzGbKkMd6uXR+sk2RJ4Gl0zeJS3iRnq2EleAbwLOKyqHpwor6rb28+bgAvpvr0sSFxV9Z8Dsfwl\n8OPD7jvKuAYcwaSuiRFer2FMF/vYb0OU5HnAR4EVVfWfE+UD1+su4JP01+27UVV1X1Xd35Y/Dzwx\nyU5sAtermenzNZLrleSJdAnlY1V17hRVFu4zNoqBo8X0omut3UTXHTIxuLfPpDpHs+FA/dlteR82\nHKi/if4G6oeJ6/l0A5PLJpXvAGzdlncCbqCnAcsh49plYPk1wEVteUfg5hbfDm15x4WKq9Xbi27Q\nNAtxvQbOsZTpB54PZcNB1EtGfb2GjGt3unHCF00q3xZ4ysDyvwIHL2BcPzDx70f3y/mWdu2G+gyM\nKq62feIL57YLdb3aez8d+OAMdRbsM9bbxV7ML7qZEd+g+wX9rlb2brpv/wBPAj7R/oNdAuw5sO+7\n2n7XA4cscFxfAu4Ermiv81r5i4Cr2n+qq4AjFziuPwKuaef/MvCcgX3/V7uOa4A3LWRcbf144IRJ\n+436en0cuAP4b7pvhkcCbwbe3LaH7mFzN7bzL1+g67WxuD4K3DPw+Vrdyvds1+pr7d/5XQsc11sG\nPl8XMZD0pvoMLFRcrc4b6SbvDO436uv1ErouqysH/q1eNa7PmLdpkST1xjEVSVJvTCqSpN6YVCRJ\nvTGpSJJ6Y1KRJPXGpCJJ6o1JRZLUG5OKtAlozwj563HHIc2XSUUaofbwpJ9uy+9NctIcjvGjSf5l\nYH2/JBfMtI80Lt6lWBqt44B3t7vTPh84bA7HuAZ4dpItqruz8p8Av9tjjFJvTCrSCFXVV9oD3X4H\nOKA2vN0+SS6muyHpdsCOSa5om46pqi+2Y3wvyTXAPkmWAbdUeyiTtKkxqUgjlORH6R73+s3qnsq3\ngap6Yat3APDGqnrjNIe6CHgx3aOte3v6o9Q3x1SkEWmPa/0Y3aNcH0hy0DwOdxHwXuCTVXVbH/FJ\no2BSkUagPVL2XOB3q+o64D10t92fq68DDwLvn3900uh463tpEUjyIeCrVXXauGORZmJLRdqEJXl2\nkq8DTzahaDGwpSJJ6o0tFUlSb0wqkqTemFQkSb0xqUiSemNSkST1xqQiSeqNSUWS1Jv/D8LjHQIC\nvipzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(514229)\n",
    "\n",
    "from scipy.stats import uniform\n",
    "\n",
    "x = uniform.rvs(size=10**5)\n",
    "y = uniform.rvs(size=10**5)\n",
    "\n",
    "t = x + y\n",
    "\n",
    "plt.hist(t, bins=20)\n",
    "\n",
    "plt.title(r'Histogram of $T = X + Y$')\n",
    "plt.xlabel(r'$x + y$')\n",
    "plt.ylabel('Frequency')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histogram looks like an ascending and then descending staircase, a discrete approximation to a triangle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayes' billiards\n",
    "\n",
    "In the Bayes' billiards story, we have $n$ white balls and 1 gray ball, throw them onto the unit interval completely at random, and count the number of white balls to the left of the gray ball. Letting $p$ be the position of the gray ball and $X$ be the number of white balls to the left of the gray ball, we have\n",
    "\n",
    "\\begin{align}\n",
    "  p &\\sim Unif(0, 1) \\\\\n",
    "  X|p &\\sim Bin(n, p)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By performing this experiment a large number of times, we can verify the results we derived in this chapter about the marginal PMF of $X$ and the posterior PDF of $p \\text{ given } X = x$. We'll let the number of simulations be called `nsims`, to avoid a name conflict with the number of white balls, `n`, which we set equal to 10:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsims = 10**5\n",
    "n = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We simulate 10<sup>5</sup> values of $p$, then simulate 10<sup>5</sup> values from the conditional distribution of $X$ given $p$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(832040)\n",
    "\n",
    "p = uniform.rvs(size=nsims)\n",
    "\n",
    "from scipy.stats import binom\n",
    "\n",
    "x = binom.rvs(n, p, size=nsims)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we feed the entire array `p` into `binom.rvs`. This means that the first element of `x` is generated using the first element of `p`, the second element of `x` is generated using the second element of `p`, and so forth. Thus, conditional on a particular element of `p`, the corresponding element of `x` is Binomial, but the elements of `p` are themselves Uniform, exactly as the model specifies.\n",
    "\n",
    "According to the Bayes' billiards argument, the marginal distribution of $X$ should be Discrete Uniform on the integers 0 through $n$. Is this in fact the case? We can make a histogram of `x` to check! Because the distribution of $X$ is discrete, we can specify `bins=np.arange(0, n+1, 0.5), align='left'` in the call to `matplotlib.pyplot.hist` so that each bar is centered at an integer value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEYCAYAAABLOxEiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAFztJREFUeJzt3XuUZWV95vHvIw1yUxuhIdCADUkP\nETUJPR0kwXEM7QU0AjoywTXR1mDILDFqkhkFV2ZwVNZoYsSQjCYouJCoCIiCl2haETPOGsHmErnp\n0ArSDUi34aagYONv/jhvmUOnuurstvapqq7vZ62zau/3vHu/v9OX89R+9z77pKqQJGlUj5vtAiRJ\n84vBIUnqxOCQJHVicEiSOjE4JEmdGBySpE4MDklSJwaHJKkTg0MLQpIbkzxntuuYLUkOSXJtkh8k\nef1s16P5zeDQvJfktiTP3aLtVUm+OrFeVU+rqiu67mc78ibgiqp6QlWdNdvFaH4zOKQxSLJolkt4\nCnDjLNeg7YTBoQVh+GgiyZuT3NGmbb6VZFWS84EDgU8n+WGSN7W+T01yRZL72nTXsUP7XDE0/XNR\nko8neccWY745yTeAB5MsSnJqkm+3bW5K8pIt+v/XJN9I8mCSc5Lsk+TvW/8vJtljitc4aa1JLgd+\nC/jr9tr+zSTb/lmSTw6t/3mSLyXZ8ef4Y9f2qqp8+JjXD+A24LlbtL0K+OqWfYBDgPXAfq19GfCL\nk+0H2BFYB7wF2Ak4CvhB28dOwHeBN7R+LwUeAd6xxZjXAQcAu7S2E4D9GPzS9jvAg8C+Q/2/BuwD\nLAU2AtcAhwGPBy4HTt/Kn8FWa23PXwG8Zoo/wz2B+4BfA/4zcD3wpNn+u/UxNx+zffgszZRPJdk8\ntL4TgzfdLT3K4E340CSbquq2KfZ5BLA78M6q+ilweZLPAC9n8Ca+CDirqgq4JMlVk+zjrKpaP7FS\nVRcNPffxJKcBhwOXtra/qqq7AZL8b2BjVV3b1j8JrNqGWt86xWucqOufk7wX+DDwJOBZVXX/dNtp\nYXKqStuL46tq8cQDeO1knapqHfBGBm+mG5NckGS/rexzP2B9eyOe8F0GRwP7AXe00Jiwnn/tMW1J\nXpnkujaddB/wdGCvoS53Dy3/aJL13beh1lFdCzwDOG047KQtGRxacKrqo1X1LAYnjAt418RTW3S9\nEzggyfD/kwOBO4C7gKVJMvTcAZMNN7GQ5CnAB4DXAXu2gLsByCTbdTVVrdNK8gzg/cB5wO/NQD3a\njhkcWlDa5xmOSvJ44McMfot/tD19N3DwUPcrGZyDeFOSHdvnQF4MXAD837bd69pJ7+MYTDlNZTcG\nQbKp1fJqBkccM2GqWqeUZCnwaQbnNl4LPGMhf+ZF0zM4tNA8Hngn8H3ge8DeDE4oA/xP4E/bNNJ/\nqapHgGOBY1r/9wGvrKpvtudeCpzE4KTy7wKfAR7e2sBVdRPwFwxC524G00L/ZyZe1FS1TrVdkicC\nnwPeU1WXVdVDwJ8DZ8xEXdo+5bFTtJK2VZIrgb+pqg/Ndi1SnzzikLZRkn+f5BfaVNVq4FeAz892\nXVLfvBxX2naHABcyuNLp28DLququ2S1J6p9TVZKkTpyqkiR1sl1OVe211161bNmy2S5DkuaVq6++\n+vtVtWS6fttlcCxbtoy1a9fOdhmSNK8k+e4o/ZyqkiR1YnBIkjoxOCRJnRgckqRODA5JUicGhySp\nE4NDktSJwSFJ6sTgkCR1sl1+cvzntezUz27Tdre980UzXImkYdv6fxP8/zmTDA6NlaEszX8GhzTD\n/K1YP6+5/m/I4NB2ba7/B5wpC+V1am4wOOYIp3A03yyUsFoor7MLg2MBM6wkbQsvx5UkdWJwSJI6\nMTgkSZ0YHJKkTgwOSVInBockqRODQ5LUicEhSerE4JAkdWJwSJI6MTgkSZ0YHJKkTgwOSVInBock\nqRODQ5LUicEhSerE4JAkdWJwSJI66TU4kvxRkhuT3JDkY0l2TnJQkiuT3JLk40l2an0f39bXteeX\nDe3ntNb+rSQv6LNmSdLUeguOJEuB1wMrq+rpwA7AicC7gDOrajlwL3BS2+Qk4N6q+iXgzNaPJIe2\n7Z4GHA28L8kOfdUtSZpa31NVi4BdkiwCdgXuAo4CLm7Pnwcc35aPa+u051clSWu/oKoerqpbgXXA\n4T3XLUnait6Co6ruAN4N3M4gMO4Hrgbuq6rNrdsGYGlbXgqsb9tubv33HG6fZJufSXJykrVJ1m7a\ntGnmX5AkCeh3qmoPBkcLBwH7AbsBx0zStSY22cpzW2t/bEPV2VW1sqpWLlmyZNuKliRNq8+pqucC\nt1bVpqr6CXAJ8JvA4jZ1BbA/cGdb3gAcANCefxJwz3D7JNtIksasz+C4HTgiya7tXMUq4Cbgy8DL\nWp/VwKVt+bK2Tnv+8qqq1n5iu+rqIGA5cFWPdUuSprBo+i7bpqquTHIxcA2wGbgWOBv4LHBBkne0\ntnPaJucA5ydZx+BI48S2nxuTXMggdDYDp1TVo33VLUmaWm/BAVBVpwOnb9H8HSa5KqqqfgycsJX9\nnAGcMeMFSpI685PjkqRODA5JUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIkjoxOCRJnRgckqRODA5J\nUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIkjoxOCRJnRgckqRODA5JUicGhySpE4NDktSJwSFJ6sTg\nkCR1YnBIkjoxOCRJnRgckqRODA5JUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIkjoxOCRJnRgckqRO\nDA5JUicGhySpE4NDktRJr8GRZHGSi5N8M8nNSX4jyZOTrElyS/u5R+ubJGclWZfkG0lWDO1ndet/\nS5LVfdYsSZpa30ccfwl8vqp+GfhV4GbgVOBLVbUc+FJbBzgGWN4eJwPvB0jyZOB04JnA4cDpE2Ej\nSRq/3oIjyROBZwPnAFTVI1V1H3AccF7rdh5wfFs+DvhwDXwNWJxkX+AFwJqquqeq7gXWAEf3Vbck\naWp9HnEcDGwCPpTk2iQfTLIbsE9V3QXQfu7d+i8F1g9tv6G1ba39MZKcnGRtkrWbNm2a+VcjSQL6\nDY5FwArg/VV1GPAg/zItNZlM0lZTtD+2oersqlpZVSuXLFmyLfVKkkYwUnAkefo27HsDsKGqrmzr\nFzMIkrvbFBTt58ah/gcMbb8/cOcU7ZKkWTDqEcffJLkqyWuTLB5lg6r6HrA+ySGtaRVwE3AZMHFl\n1Grg0rZ8GfDKdnXVEcD9bSrrC8Dzk+zRToo/v7VJkmbBolE6VdWzkiwHfg9Ym+Qq4ENVtWaaTf8Q\n+EiSnYDvAK9mEFYXJjkJuB04ofX9HPBCYB3wUOtLVd2T5O3A11u/t1XVPaO+QEnSzBopOACq6pYk\nfwqsBc4CDksS4C1VdclWtrkOWDnJU6sm6VvAKVvZz7nAuaPWKknqz6jnOH4lyZkMPodxFPDiqnpq\nWz6zx/okSXPMqEccfw18gMHRxY8mGqvqznYUIklaIEYNjhcCP6qqRwGSPA7Yuaoeqqrze6tOkjTn\njHpV1ReBXYbWd21tkqQFZtTg2Lmqfjix0pZ37ackSdJcNmpwPLjF3Wr/LfCjKfpLkrZTo57jeCNw\nUZKJT2zvC/xOPyVJkuayUT8A+PUkvwwcwuDeUd+sqp/0WpkkaU4a+QOAwK8Dy9o2hyWhqj7cS1WS\npDlrpOBIcj7wi8B1wKOtuQCDQ5IWmFGPOFYCh7bbgkiSFrBRr6q6AfiFPguRJM0Pox5x7AXc1O6K\n+/BEY1Ud20tVkqQ5a9TgeGufRUiS5o9RL8f9SpKnAMur6otJdgV26Lc0SdJcNOpt1X+fwVe//m1r\nWgp8qq+iJElz16gnx08BjgQegMGXOgF791WUJGnuGjU4Hq6qRyZWkixi8DkOSdICM2pwfCXJW4Bd\nkjwPuAj4dH9lSZLmqlGD41RgE3A98AfA5wC/+U+SFqBRr6r6KYOvjv1Av+VIkua6Ue9VdSuTnNOo\nqoNnvCJJ0pzW5V5VE3YGTgCePPPlSJLmupHOcVTVPw897qiq9wJH9VybJGkOGnWqasXQ6uMYHIE8\noZeKJElz2qhTVX8xtLwZuA34jzNejSRpzhv1qqrf6rsQSdL8MOpU1R9P9XxVvWdmypEkzXVdrqr6\ndeCytv5i4B+B9X0UJUmau7p8kdOKqvoBQJK3AhdV1Wv6KkySNDeNesuRA4FHhtYfAZbNeDWSpDlv\n1COO84GrknySwSfIXwJ8uLeqJElz1qhXVZ2R5O+Bf9eaXl1V1/ZXliRprhp1qgpgV+CBqvpLYEOS\ng3qqSZI0h4361bGnA28GTmtNOwJ/11dRkqS5a9QjjpcAxwIPAlTVnXjLEUlakEYNjkeqqmi3Vk+y\n26gDJNkhybVJPtPWD0pyZZJbknw8yU6t/fFtfV17ftnQPk5r7d9K8oJRx5YkzbxRg+PCJH8LLE7y\n+8AXGf1Lnd4A3Dy0/i7gzKpaDtwLnNTaTwLurapfAs5s/UhyKHAi8DTgaOB9SXYYcWxJ0gwb9bbq\n7wYuBj4BHAL896r6q+m2S7I/8CLgg209DG7HfnHrch5wfFs+rq3Tnl/V+h8HXFBVD1fVrcA64PBR\n6pYkzbxpL8dtv91/oaqeC6zpuP/3Am/iX86H7AncV1Wb2/oGYGlbXkq7hUlVbU5yf+u/FPja0D6H\ntxmu82TgZIADDzywY5mSpFFNe8RRVY8CDyV5UpcdJ/ltYGNVXT3cPNkQ0zw31TbDdZ5dVSurauWS\nJUu6lCpJ6mDUT47/GLg+yRralVUAVfX6KbY5Ejg2yQsZfN3sExkcgSxOsqgddewP3Nn6bwAOYPAZ\nkUXAk4B7htonDG8jSRqzUU+Ofxb4bwzuiHv10GOrquq0qtq/qpYxOLl9eVX9J+DLwMtat9XApW35\nsrZOe/7ydiXXZcCJ7aqrg4DlwFUj1i1JmmFTHnEkObCqbq+q86bq19GbgQuSvAO4FjintZ8DnJ9k\nHYMjjRMBqurGJBcCNzH49sFT2vSZJGkWTDdV9SlgBUCST1TVf9iWQarqCuCKtvwdJrkqqqp+DJyw\nle3PAM7YlrElSTNruqmq4RPTB/dZiCRpfpguOGory5KkBWq6qapfTfIAgyOPXdoybb2q6om9VidJ\nmnOmDI6q8tYekqTH6PJ9HJIkGRySpG4MDklSJwaHJKkTg0OS1InBIUnqxOCQJHVicEiSOjE4JEmd\nGBySpE4MDklSJwaHJKkTg0OS1InBIUnqxOCQJHVicEiSOjE4JEmdGBySpE4MDklSJwaHJKkTg0OS\n1InBIUnqxOCQJHVicEiSOjE4JEmdGBySpE4MDklSJwaHJKkTg0OS1InBIUnqxOCQJHVicEiSOukt\nOJIckOTLSW5OcmOSN7T2JydZk+SW9nOP1p4kZyVZl+QbSVYM7Wt1639LktV91SxJml6fRxybgT+p\nqqcCRwCnJDkUOBX4UlUtB77U1gGOAZa3x8nA+2EQNMDpwDOBw4HTJ8JGkjR+vQVHVd1VVde05R8A\nNwNLgeOA81q384Dj2/JxwIdr4GvA4iT7Ai8A1lTVPVV1L7AGOLqvuiVJUxvLOY4ky4DDgCuBfarq\nLhiEC7B367YUWD+02YbWtrX2Lcc4OcnaJGs3bdo00y9BktT0HhxJdgc+Abyxqh6YquskbTVF+2Mb\nqs6uqpVVtXLJkiXbVqwkaVq9BkeSHRmExkeq6pLWfHebgqL93NjaNwAHDG2+P3DnFO2SpFnQ51VV\nAc4Bbq6q9ww9dRkwcWXUauDSofZXtqurjgDub1NZXwCen2SPdlL8+a1NkjQLFvW47yOBVwDXJ7mu\ntb0FeCdwYZKTgNuBE9pznwNeCKwDHgJeDVBV9yR5O/D11u9tVXVPj3VLkqbQW3BU1VeZ/PwEwKpJ\n+hdwylb2dS5w7sxVJ0naVn5yXJLUicEhSerE4JAkdWJwSJI6MTgkSZ0YHJKkTgwOSVInBockqROD\nQ5LUicEhSerE4JAkdWJwSJI6MTgkSZ0YHJKkTgwOSVInBockqRODQ5LUicEhSerE4JAkdWJwSJI6\nMTgkSZ0YHJKkTgwOSVInBockqRODQ5LUicEhSerE4JAkdWJwSJI6MTgkSZ0YHJKkTgwOSVInBock\nqRODQ5LUicEhSerE4JAkdWJwSJI6mTfBkeToJN9Ksi7JqbNdjyQtVPMiOJLsAPwv4BjgUODlSQ6d\n3aokaWGaF8EBHA6sq6rvVNUjwAXAcbNckyQtSKmq2a5hWkleBhxdVa9p668AnllVrxvqczJwcls9\nBPhWT+XsBXy/p307pmM6pmPO5phPqaol03Va1MPAfcgkbY9JvKo6Gzi790KStVW1su9xHNMxHdMx\n59KYw+bLVNUG4ICh9f2BO2epFkla0OZLcHwdWJ7koCQ7AScCl81yTZK0IM2Lqaqq2pzkdcAXgB2A\nc6vqxlkqp/fpMMd0TMd0zDk45s/Mi5PjkqS5Y75MVUmS5giDQ5LUicExotm45UmSc5NsTHLDmMY7\nIMmXk9yc5MYkbxjDmDsnuSrJP7Ux/0ffYw6NvUOSa5N8Zkzj3Zbk+iTXJVk7pjEXJ7k4yTfb3+tv\n9DzeIe31TTweSPLGPsds4/5R+/dzQ5KPJdl5DGO+oY13Y1+vcbL3gCRPTrImyS3t5x59jD2lqvIx\nzYPBCflvAwcDOwH/BBw6hnGfDawAbhjT69wXWNGWnwD8v75fJ4PP6OzelncErgSOGNPr/WPgo8Bn\nxjTebcBe4xhraMzzgNe05Z2AxWMcewfgeww+VNbnOEuBW4Fd2vqFwKt6HvPpwA3ArgwuMvoisLyH\ncf7VewDwZ8CpbflU4F3j/DdVVR5xjGhWbnlSVf8I3NP3OEPj3VVV17TlHwA3M/hP2eeYVVU/bKs7\ntkfvV2wk2R94EfDBvseaLUmeyOCN5xyAqnqkqu4bYwmrgG9X1XfHMNYiYJckixi8mff9Oa+nAl+r\nqoeqajPwFeAlMz3IVt4DjmPwCwHt5/EzPe50DI7RLAXWD61voOc31NmWZBlwGIMjgL7H2iHJdcBG\nYE1V9T4m8F7gTcBPxzDWhAL+IcnV7RY5fTsY2AR8qE3JfTDJbmMYd8KJwMf6HqSq7gDeDdwO3AXc\nX1X/0POwNwDPTrJnkl2BF/LYDyn3aZ+qugsGv+wBe49p3J8xOEYz7S1PtidJdgc+Abyxqh7oe7yq\nerSqfo3BHQEOT/L0PsdL8tvAxqq6us9xJnFkVa1gcJfnU5I8u+fxFjGY5nh/VR0GPMhgaqN37YO6\nxwIXjWGsPRj8Fn4QsB+wW5Lf7XPMqroZeBewBvg8g+nrzX2OOZcYHKNZMLc8SbIjg9D4SFVdMs6x\n2zTKFcDRPQ91JHBsktsYTDseleTveh6Tqrqz/dwIfJLBFGifNgAbho7gLmYQJONwDHBNVd09hrGe\nC9xaVZuq6ifAJcBv9j1oVZ1TVSuq6tkMppNu6XvM5u4k+wK0nxvHNO7PGByjWRC3PEkSBvPhN1fV\ne8Y05pIki9vyLgzeBL7Z55hVdVpV7V9Vyxj8XV5eVb3+hppktyRPmFgGns9guqM3VfU9YH2SQ1rT\nKuCmPscc8nLGME3V3A4ckWTX9m94FYPzc71Ksnf7eSDwUsb3ei8DVrfl1cClYxr3Z+bFLUdmW83S\nLU+SfAx4DrBXkg3A6VV1To9DHgm8Ari+nXMAeEtVfa7HMfcFzmtf1vU44MKqGsvlsWO2D/DJwfsa\ni4CPVtXnxzDuHwIfab/wfAd4dd8Dtjn/5wF/0PdYAFV1ZZKLgWsYTBddy3huyfGJJHsCPwFOqap7\nZ3qAyd4DgHcCFyY5iUFonjDT405bV7ukS5KkkThVJUnqxOCQJHVicEiSOjE4JEmdGBySpE4MDklS\nJwaHJKkTg0Mak/ZdJ89ry+9IctZs1yRtCz85Lo3P6cDb2q0qDmNwE0Bp3vGT49IYJfkKsDvwnPad\nJ9K841SVNCZJnsHg3lwPGxqazwwOaQza7a8/wuB7Ix5M8oJZLknaZgaH1LN2t9hLgD9pXwD0duCt\ns1qU9HPwHIckqROPOCRJnRgckqRODA5JUicGhySpE4NDktSJwSFJ6sTgkCR18v8BkmSlt64WfKUA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(x, bins=np.arange(0, n+1, 0.5), align='left')\n",
    "\n",
    "plt.title(r'Histogram of $x$')\n",
    "plt.xlabel(r'$x$')\n",
    "plt.ylabel('Frequency')\n",
    "# also set the x-axis ticks to show integers 0, 1, 2, ..., 10\n",
    "plt.xticks(range(0, 10+1))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, all the histogram bars are approximately equal in height, consistent with a Discrete Uniform distribution. \n",
    "\n",
    "Now for the posterior distribution of $p$ given $X = x$. Conditioning is very simple since we are using `numpy.array`. To consider only the simulated values of $p$ where the value of $X$ was 3, we use [`numpy.where(x==3)`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.where.html) to find the indices where `x` = 3, and then use that within the square brackets of `p` to index only those corresponding values. In particular, we can create a histogram of these values using `matplotlib.pyplot.hist` to see what the posterior distribution of $p$ given $X = 3$ looks like.\n",
    "\n",
    "According to the Beta-Binomial conjugacy result, the true posterior distribution is $p|X = 3 \\sim Beta(4, 8)$. We can plot the histogram of `p[numpy.where(x==3)]` next to a histogram of simulated values from the $Beta(4, 8)$ distribution to confirm that they look similar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA00AAAF5CAYAAABDd8G+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XmYZHV97/H3R0ZRXFhHMw6QwTga\njdd1XBKDIaJGTQLGYERNHIQ4WVxjcgWXezUmuUGTG6LXXHUMI4NRXIiJ6EUNLijmEXBUgoIaRiQw\nMsIoiwsiIt/7xzkNNUX1merpqq7q7vfreerpOr/zq1Pfqq6qb33P73dOpaqQJEmSJA12h0kHIEmS\nJEnTzKJJkiRJkjpYNEmSJElSB4smSZIkSepg0SRJkiRJHSyaJEmSJKmDRZMkSZIkdbBokiRJ0tRL\n8qwkOyYdh5YniyZJkiRNtSR3AI4Crph0LFqeLJokSdJuSXJRksMmHcekJLl/ki8l+X6SF086HoAk\nf53kpQt0X+cn+YWFuC/g2cDpwC27c+NxvlaTXJbkCePY9oD7OiXJX454mwsW/2Jm0bRMmejml+iS\nHJbkliQ/SPLkccQ4DSb5QZrkgCRfS3LnSdz/MJL8c5KnDGjfs31t/GTUyU1aKIPe/0mOSfLZmeWq\n+oWqOnuu21lCXg6cXVV3r6o39a9Msm+Saj8PfpDk+iSnJ9lr2DtI8q0kDx2y70rgucDb+trXJrkx\nyT8Ne7/t7dYkOTPJtUm+neTNSVb0dPlb4HVz2N4Lk2xJ8uMkpwxYv1+Sf0nywyT/leTZbfsewO8A\n751L/L2Gea2O2xJ/Lyx5Fk1LkIluKJ2JbkhXVtXdquqjowxMt3o5cEpV3TipAJL8U/tF4fok5yX5\nxb4uJwL/q/92VfXjqrob8K4FCVRapvq+wE/CzwIXdax/KPDdNlfcDfh54LHA7w2z8SQHAPcEvjpk\nPMcAZ1bVj/ra/wH4/JDb6PV/gauBVTSP5VeAP+5Zfwbwq0lWDbm9K4G/BDbNsv4fgJuAewHPAd7S\njmT9LvC+qtqtUSZpFCyaNBGLINFpgpLckSb5z2mv6Bj8DXBIVe0N/A/gg+0eTwCq6vPA3ZI8clIB\nSpPUu3MtyfHtqMj3k3w9yeFJ3gkcDHyoHWl5edv3AUnOTnJdO/PhiJ5tPrxnJsD7k7y3d8S2vc/j\nk1wI/DDJiiQnJPlGe5uLk/xWX///nuTCdgTj5CT3SvKRtv/Hk+zb8RgHxprkk8CvAm9uH9v9Btz8\nofTkmqraTnNMzh17tv/8Nubr25ju2bbft+17B+C7Sb6b5J5JzkhyVZLvJflQknv03N9TgE/3xX80\ncB3widkeY4dDaIqVG6vq28BHgVun47U7tb4APGmYjVXVB6rqX4Hv9q9Lclfgt4H/UVU/qKrP0hRl\nvwc8EHhuko8Ca5PMurNz0Ouwbe99rc7pNZFmtPC+PcuzTpGb7bXY8V64d5pZCzuSfDN9s1+SPCzJ\nF9vtvRcYOPuivd/T+9reOPNcdb1HBmxr1sc7RLwDn/+lwKJpmTLRdSe69r5fn+TKJJcnGXav4GVJ\n/jzJ1iTXJPnD2Z73drmS3DfNSOCWJNuSvCvN9IRvJ/lvXfEkeXTbb4+ebZ6UjoTSF8+2JA/oa0v7\nnK1umw7teDz7pRmNuar98Hz+gMf7wiRfaP+HZwxzO+DRwPeqalvf9vZPckWSX2qXV7av3UcN83jn\nqqr+o6p+lCQ0X3JWAgf0dTsbeOo47l9aLJLcH3gh8Miqujvwa8BlVfV7wOXAb7ajLW9Is1PkQ8C/\n0YyivAh4V5pp03cC/gU4BdgPOA0Y9OXuWcCvA/tU1c3AN4BDgb2BPwf+KTuPfvw28ETgfsBvAh8B\nXknzfr4DMHCadlesVfV44Bzghe1j+88Bm3gYbdHU5rxnAgcB72vbXgn8IXAEzefLt2hGYqiqrcCf\nAae3298fuAfwf2jy85o2/j/oub//Bny9J/570Eyf+9NBj28IbwSOTrJXmxOeQlM49foq8JDd3H6v\n+wE/7Xse/wP4hao6vqqeVFVPBi6pqtn+XwNfh7Pc3269JoYw8LU4y3vhDjSvr/8AVgOHAy9N8mvt\n47kT8K/AO2neD+9v4x7kNOCp7f+8d0rju7vimssDGyLeuTz/i45F0zJnops10QH8Ms1Uit8G3prk\noFn69Xs08GDgOOANGX5U7SbgQe39vR04GXhmVzxVdR7wA+AwaAoemrMLnTbkfZ4HrOtrWwt8v6q+\nNcTjeSfwE5rk/QTgtUn6t/cHNHPs9+a2qWy7ut2D6Un8M6rqu20cp6TZK/k24G1VdT7cuiftugGX\nD/duZ9h+bd//C9wIfBh4f1Vd1ddlVF8YpGn0r73vEZrpWoP8FNgTeGCSO1bVZVX1jVn6Pga4G3Bi\nVd1UVZ+keX89q123AnhTVf2kqj4AnD9gG2+qqitmpqFV1fur6sqquqWq3gtcAvTuTPk/VXVV+7l2\nDnBeVX2pqn5Mk7sethuxDuOhwDHtc3cDzWf7c6rq6jQjSq8GnlVVW6vqJprP/d6R64cAF8wstP3O\naqcAXwOcBfTuPNwH+H7P8l8AJ1fV7p5x7tM0I0vfA7YBW2i+xPf6fnu/83U34Pq+tuuBu/c2VFV/\njuk1l9fh7r4mOg3xWuz1SGBlVb2ufX1dSvMaObpd/xiaHXZ/374fTmeWaZZV9V/AF4GntU2PB26o\nqnN3I67Z7CreuTz/i45F09JloptfogPYWFXfa6dgfQ4Y9oQPm6rqBuD/0XzY/8yQt/tGVV0HfIem\nYPgmzbzuXcXzbm77wPpF4Oaq+tyQ93ku8Igkd0xydTuy9Yi2vfPxtMXtU4CXVdWP2tfN6dy+WN5Y\nVRdV1c1Vde6Qt9uXnRP/rarq32iK53No5tn/Vc+6E6tqnwGX3+jbxlD92r5/3D7uZ9PuHe4zqi8M\n0jR6Wu97hJ2PZ7lVNaMiLwVeC1yd5D1J7j3LNu8NXFE7H5/yXzR7ru8NfKuqqmfdoC/8O7UleW6S\nC3py3oPYeVS4d2fHjwYs3203Yu2UZE/gAcDj2ufuLjTHQW5uuxwO3Ak4vyfuj7Jz4fBQmr36M9t8\nRpJ/bz+vrwNOAHp3/F1LW2SkOXnEE4CTdhXrLPHfAfgY8AHgrjTP577A6/u63p1m+t98/YBmJK3X\nPZglFwwyx9fh7r4mOg3xWuz1s8C9+76vvZLbcv+g98N/ddz9u7nte86zuW2Uaa5xzaYz3jk+/4uO\nRdPSZaLbzUTX4+q+6/earWOfawDavYYwy/zjAX7a/r25vX4zTSG6q3jeDTy9HV37HYYfZYKmOFpH\nU2ReSTOqt46di6bZHs/MyNs3e/5/x3H75+mSvuVhbndr4p/F22gK5k1V9dOOfiPRFt6nAa9L33RG\nmqQ+ii8M0qJWVe+uql+m+WJV3Pbluvq6Xgkc1H4pn3EwzdS07cDqdtR8xqBR/lu3meRnafZ2vxDY\nv815XwEy4HZz1RXrrjyI5vTYFwK0n1WfBu7Vfl7vB/xL346bvavqsPZx3aHdxgXt8uNpntOX0uS4\nA2hywQU993khzcwMaGYgrAEuT/Jtmql+v53ki0M+9v1onvs3tyNb3wXewe2nIz+AnsJuHv4TWJFk\nbU/bQ5jj8ccdr8P5uAHoPePhwJ2hQ7wW+98LVwDf7HsN3L2qZp7jQe+HgzvifD9wWJIDaXZEvnvI\nuIZ9vLuKd1zP/1SwaJKJbna9X+Lvyc5Fy+66kbYQys4H786m97kYGE9VfY2mIHwSzdS9uRRNW2iS\n8hOBVwGP4/ZF02yuoHk8+/d8eN6tqn6/r9/Nu3G7C4H7D7rTdp7224E3AX+e5Gd61r0yt53at/fy\nkb5tDNVvgD3oOQi69YA2XmnZSjNN+/Ht6MqNNDu1ZnZoXAXcp6f7ecAPgZe3o9yH0Uy/fg/NKPpP\ngRemOQboSHY9heiuNLllRxvL82g+10ahK9ZdeRjw5ZkdTu0o+6uAM6rqJzRTqX41ycPb9fdIcmRP\nHr1Le5nJYw+h+fz8Gs2IzyaaXHBxz32eSXOGO4CNwM/RjFY9FHgrzYyBX5vpnOYA/1MGBV9V36GZ\n8fBH7f9iH2A9O4987UkzO+GsXW2vXb8izc9I7AHskeTOaad8V9UPaUa1XpfkrkkeCxxJM517KLt4\nHc7HBcCzk+yR5mdGfmWWfrt6Lfa/F84HvpfmeO+7tNt/UG47udDnaHLoi9vn7ul0vB+qagfNcbbv\noCluZs66ONf3yGyPtzPeMT7/U8GiaZkz0XV6fpvEHkkz7e1jI4ird/rhUXO8bVc87wLeAFxfVV/u\nv2GS1WlO5vD03vZ2muQ3aKb3/RvN3OkH0ZwNqVM1Z4H6NHBim+DumOSXknQe3zPk7c4H9m73lvV7\nJXBdVb2E5kvAO2a+ZFTV/2oLsP7LTr+lNEy/JD+T5Lgke7ev6T+g2Wvb/9wcRvNFRVrO9qSZevYd\n4Ns0X+Zf2a77a+DVaUaW/6wtIo6gmab7HZrp48+tqq+1655OM/p8Hc2ppj8M/Hi2O66qi4H/TZOH\nrqI5GcK/j+JBdcU6xM0fCjy03SFzfRvTRcDz2m1/juYkDf+c5Ac0xc+TZ2ZstEXEW4GLk2yj+Zy/\nI83z+2GafHJxzywAgFNpTgZwl6q6oaq+PXOhmf52Y/vFesZBdD9XT6eZCr4D2ErzBf5PetYfQfPz\nHVcOub1X03zPOIHmf/ujtm3GH9MUilfT7AD8o6qay0hT1+twPl5C8x3iOppTofcf1wUM9Vrsfy/8\ntN3uQ2kK1O8A/0hzDDA974djaGZgPJOmsOzybpppmbdOzduN98jAx7ureBnf8z8dqsrLErvQnKnk\nCX1txwCf7e9Dc8D9+TRzhq+h+SC+d9vnSJqTQVwH/Fnb9gs0X3ivp/mA/62eba6j2TvxA5oh4g/Q\nnDq0K66/au/3O8Dftdv+/UH9aU4//dqe5d8HPt7xPHTFevbM/XQ8h6+nGbG6Ajimb/1hwLZdPfc0\nReF9e5YfSZPoPkPzoVTAfdv/zyk921gzoK0rntU0ReurZ3k8a9r7OmbAujcDH2qvvwr4/Bwez/40\n8/O3t6+Tc4BHdP3Ph7ld2+cNwAl9bY+i+cBf1S6voDko9kVjeB+tBD7Zvja/T/PaPmJAPF/q2MYp\nwF+O433uxctyudDsBHvepONYLBeaE+68dIh+d6I5kc0d5/m/edCotufFyzRfUtU/A0sajSTnAW+t\nqndMOpa5SnIZTVH18VnWP45mpOfHwDOrahSjUPOJ5840ewIfXlX9xxAtSml+1PGzwMPq9j/UOBWS\n/DPNmanO7Gvfk6a4uyPwhqr680nEJy1GSX6F5mQ436HZy/1W4D7VjFJL0kRM+gdGtYQMSHQP5va/\n57AkVNVnaKYQTIvnAxculYIJbp1P//OTjqNLVQ38vYxqzu7oGfWk3XN/mjNV3o1m+vBRFkySJs2i\nSaNkopuAJJfS/ObRcyYdiyTNV1VtpDmJgSRNDafnSZIkSVIHz54nSZIkSR0smiRJy1KSTUmuTvKV\nnra/SfK1JBcm+Zf2t2lm1r2iPXX/15P0/tbNk9u2rUlOWOjHIUkavyU5Pe+AAw6oNWvWTDoMSVrW\nvvCFL3ynqlZOOo7ZtGfB/AFwalU9qG17EvDJqro5yesBqur4JA+k+d2YRwH3Bj4O3K/d1H/S/ED0\nNppT4D+rmt9FmZV5SpKmw7C5akmeCGLNmjVs2bJl0mFI0rKW5L8mHUOXqvpMkjV9bf/Ws3gut/0I\n9ZHAe9ozI34zyVZu+6HqrVV1KUCS97R9O4sm85QkTYdhc5XT8yRJGuxY4CPt9dU0Pyw9Y1vbNlu7\nJGkJsWiSJKlPklcBNwPvmmka0K062gdtc0OSLUm27NixYzSBSpIWhEWTJEk9kqwHfgN4Tt124O82\n4KCebgcCV3a0305VbayqdVW1buXKqT3US5I0gEWTJEmtJE8GjgeOqKobeladARydZM8khwBrgfNp\nTvywNskhSe4EHN32lSQtIUvyRBCSJO1KktOAw4ADkmwDXgO8AtgTOCsJwLlV9YdVdVGS99Gc4OFm\n4AVV9dN2Oy8EPgbsAWyqqosW/MFIksbKokmStCxV1bMGNJ/c0f+vgL8a0H4mcOYIQ5MkTRmn50mS\nJElSB4smSZIkSepg0SRJkiRJHSyaJEmSJKmDRZMkSZIkdbBokiRJkqQOFk2SJEmS1MHfadKStP3Y\n+e8PWLXplhFEIknSYOsP3Tiv228+Z8OIIpG0K440SZIkSVIHiyZJkiRJ6mDRJEmSJEkdLJokSZIk\nqYNFkyRJkiR1sGiSJEmSpA4WTZIkSZLUwaJJkiRJkjpYNEmSJElSB4smSZIkSeqwYtIBSJIkLTbr\nD9046RAkLSCLJkmSpEVovoXb5nM2jCgSaelzep4kSZIkdRhb0ZRkU5Krk3ylp22/JGcluaT9u2/b\nniRvSrI1yYVJHt5zm/Vt/0uSrB9XvJIkSZI0yDhHmk4BntzXdgLwiapaC3yiXQZ4CrC2vWwA3gJN\nkQW8Bng08CjgNTOFliRJkiQthLEVTVX1GeCavuYjgc3t9c3A03raT63GucA+SVYBvwacVVXXVNW1\nwFncvhCTJEmSpLFZ6BNB3KuqtgNU1fYk92zbVwNX9PTb1rbN1q4lbvuxHm4nSZKk6TAt30wzoK06\n2m+/gWRDki1JtuzYsWOkwUmSJElavha6aLqqnXZH+/fqtn0bcFBPvwOBKzvab6eqNlbVuqpat3Ll\nypEHLkmSJGl5Wuii6Qxg5gx464EP9rQ/tz2L3mOA69tpfB8DnpRk3/YEEE9q2yRJkiRpQYztmKYk\npwGHAQck2UZzFrwTgfclOQ64HHhG2/1M4KnAVuAG4HkAVXVNkr8APt/2e11V9Z9cQpIkSZLGZmxF\nU1U9a5ZVhw/oW8ALZtnOJmDTCEOTJEmSpKFNy4kgJEmSJGkqWTRJkiRJUgeLJkmSJEnqYNEkSZIk\nSR0smiRJkiSpg0WTJEmSJHWwaJIkSZKkDhZNkiRJktTBokmSJEmSOqyYdABamrYfaz0uSZpe6w/d\nOOkQJC0ifrOVJEmSpA4WTZIkSZLUwaJJkiRJkjpYNEmSJElSB4smSZIkSepg0SRJWpaSbEpydZKv\n9LTtl+SsJJe0f/dt25PkTUm2JrkwycN7brO+7X9JkvWTeCySpPGyaJIkLVenAE/uazsB+ERVrQU+\n0S4DPAVY2142AG+BpsgCXgM8GngU8JqZQkuStHT4O03SLOb7W1OrNt0yokgkjUNVfSbJmr7mI4HD\n2uubgbOB49v2U6uqgHOT7JNkVdv3rKq6BiDJWTSF2GljDl+StIAcaZIk6Tb3qqrtAO3fe7btq4Er\nevpta9tma5ckLSEWTZIk7VoGtFVH++03kGxIsiXJlh07dow0OEnSeFk0SZJ0m6vaaXe0f69u27cB\nB/X0OxC4sqP9dqpqY1Wtq6p1K1euHHngkqTxsWiSJOk2ZwAzZ8BbD3ywp/257Vn0HgNc307f+xjw\npCT7tieAeFLbJklaQjwRhCRpWUpyGs2JHA5Iso3mLHgnAu9LchxwOfCMtvuZwFOBrcANwPMAquqa\nJH8BfL7t97qZk0JIkpYOiyZJ0rJUVc+aZdXhA/oW8IJZtrMJ2DTC0CRJU8bpeZIkSZLUwaJJkiRJ\nkjpYNEmSJElSB4smSZIkSepg0SRJkiRJHTx7niRJ0jK0/tCN87r95nM2jCgSafo50iRJkiRJHSya\nJEmSJKmDRZMkSZIkdbBokiRJkqQOFk2SJEmS1MGiSZIkSZI6WDRJkiRJUgeLJkmSJEnqYNEkSZIk\nSR0smiRJkiSpg0WTJEmSJHWwaJIkSZKkDhZNkiRJktTBokmSJEmSOlg0SZIkSVKHFZMOQNNp+7HW\n0/M13+dw1aZbRhSJJEmS5sNvxpIkSZLUwaJJkiRJkjpYNEmSJElSB4smSZIkSeowkaIpyZ8kuSjJ\nV5KcluTOSQ5Jcl6SS5K8N8md2r57tstb2/VrJhGzJEmSpOVpwYumJKuBFwPrqupBwB7A0cDrgZOq\nai1wLXBce5PjgGur6r7ASW0/SZIkSVoQk5qetwK4S5IVwF7AduDxwOnt+s3A09rrR7bLtOsPT5IF\njFWSJEnSMrbgRVNVfQv4W+BymmLpeuALwHVVdXPbbRuwur2+Griive3Nbf/9+7ebZEOSLUm27Nix\nY7wPQpIkSdKyMYnpefvSjB4dAtwbuCvwlAFda+YmHetua6jaWFXrqmrdypUrRxWuJEmSpGVuEtPz\nngB8s6p2VNVPgA8AvwTs007XAzgQuLK9vg04CKBdvzdwzcKGLEmSJGm5mkTRdDnwmCR7tccmHQ5c\nDHwKOKrtsx74YHv9jHaZdv0nq+p2I02SJEmSNA6TOKbpPJoTOnwR+HIbw0bgeOBlSbbSHLN0cnuT\nk4H92/aXAScsdMySJEmSlq8Vu+4yelX1GuA1fc2XAo8a0PdG4BkLEZckSZIk9ZvUKcclSZIkaVGY\nyEiTJEnSfKw/dOOkQ5C0jDjSJEmSJEkdLJokSZIkqYNFkyRJkiR18JgmSZIkzdl8jyvbfM6GEUUi\njZ8jTZIkSZLUwaJJkiRJkjpYNEmSJElSB4smSZIkSepg0SRJUp8kf5LkoiRfSXJakjsnOSTJeUku\nSfLeJHdq++7ZLm9t16+ZbPSSpFGzaJIkqUeS1cCLgXVV9SBgD+Bo4PXASVW1FrgWOK69yXHAtVV1\nX+Cktp8kaQnxlOPSlNp+7Pz2aazadMuIIpGWpRXAXZL8BNgL2A48Hnh2u34z8FrgLcCR7XWA04E3\nJ0lV1UIGLEkaH0eaJEnqUVXfAv4WuJymWLoe+AJwXVXd3HbbBqxur68Grmhve3Pbf//+7SbZkGRL\nki07duwY74OQJI2URZMkST2S7EszenQIcG/grsBTBnSdGUlKx7rbGqo2VtW6qlq3cuXKUYUrSVoA\nFk2SJO3sCcA3q2pHVf0E+ADwS8A+SWamtR8IXNle3wYcBNCu3xu4ZmFDliSNk0WTJEk7uxx4TJK9\nkgQ4HLgY+BRwVNtnPfDB9voZ7TLt+k96PJMkLS0WTZIk9aiq82hO6PBF4Ms0uXIjcDzwsiRbaY5Z\nOrm9ycnA/m37y4ATFjxoSdJYefY8SZL6VNVrgNf0NV8KPGpA3xuBZyxEXJKkyXCkSZIkSZI6WDRJ\nkiRJUgeLJkmSJEnqYNEkSZIkSR0smiRJkiSpg0WTJEmSJHWwaJIkSZKkDhZNkiRJktTBokmSJEmS\nOlg0SZIkSVIHiyZJkiRJ6mDRJEmSJEkdLJokSZIkqYNFkyRJkiR1sGiSJEmSpA4WTZIkSZLUwaJJ\nkiRJkjpYNEmSJElSB4smSZIkSepg0SRJkiRJHSyaJEmSJKmDRZMkSZIkdRiqaEryoHEHIknS7jJP\nSZLGadiRprcmOT/JHyfZZ6wRSZI0d+YpSdLYDFU0VdUvA88BDgK2JHl3kieONTJJkoZknpIkjdPQ\nxzRV1SXAq4HjgV8B3pTka0mePq7gJEkalnlKkjQuwx7T9OAkJwFfBR4P/GZVPaC9ftIY45MkaZfM\nU5KkcVoxZL83A28HXllVP5pprKork7x6LJFJkjQ885QkaWyGLZqeCvyoqn4KkOQOwJ2r6oaqeufY\nopMkaTjmKUnS2Ax7TNPHgbv0LO/VtkmSNA3MU5KksRm2aLpzVf1gZqG9vtfu3mmSfZKc3h6g+9Uk\nv5hkvyRnJbmk/btv2zdJ3pRka5ILkzx8d+9XkrRkjTRPSZLUa9ii6Ye9xUqSRwA/6ui/K28EPlpV\nPw88hObA3ROAT1TVWuAT7TLAU4C17WUD8JZ53K8kaWkadZ6SJOlWwx7T9FLg/UmubJdXAc/cnTtM\ncg/gccAxAFV1E3BTkiOBw9pum4GzaU4beyRwalUVcG47SrWqqrbvzv1LkpakkeUpSZL6DVU0VdXn\nk/w8cH8gwNeq6ie7eZ/3AXYA70jyEOALwEuAe80UQlW1Pck92/6rgSt6br+tbbNokiQBI89TkiTt\nZNiRJoBHAmva2zwsCVV16m7e58OBF1XVeUneyG1T8QbJgLa6XadkA830PQ4++ODdCEuStMiNKk9J\nkrSToYqmJO8Efg64APhp21zA7iSjbcC2qjqvXT6dpmi6ambaXZJVwNU9/Q/quf2BwJX0qaqNwEaA\ndevW3a6okiQtXSPOU5Ik7WTYkaZ1wAPb44rmpaq+neSKJPevqq8DhwMXt5f1wInt3w+2NzkDeGGS\n9wCPBq73eCZJUp+R5SlJkvoNWzR9BfgZRncc0YuAdyW5E3Ap8DyaM/m9L8lxwOXAM9q+Z9L8aOFW\n4Ia2ryRJvUadpyRJutWwRdMBwMVJzgd+PNNYVUfszp1W1QU0ewX7HT6gbwEv2J37kSQtGyPNUxq/\n9YdunHQIkjS0YYum144zCEmjt/3YYX+GbbBVm24ZUSTSgnjtpAOQJC1dw55y/NNJfhZYW1UfT7IX\nsMd4Q5MkaTjmKUnSOA21KzrJ82nOcve2tmk18K/jCkqSpLkwT0mSxmnY+TsvAB4LfA+gqi4B7tl5\nC0mSFo55SpI0NsMe0/TjqropaX5nNskKBvzArKbHfI9nkaRFxjwlSRqbYb9ZfzrJK4G7JHki8H7g\nQ+MLS5KkOTFPSZLGZtii6QRgB/Bl4A9ofjvp1eMKSpKkORppnkqyT5LTk3wtyVeT/GKS/ZKcleSS\n9u++bd8keVOSrUkuTPLwkTwiSdLUGPbsebcAb28vkiRNlTHkqTcCH62qo9ofYt8LeCXwiao6MckJ\nNIXa8cBTgLXt5dHAW9q/kqQlYqiiKck3GTA3vKruM/KIJEmao1HmqST3AB4HHNNu4ybgpiRHAoe1\n3TYDZ9MUTUcCp7Y/xn5uO0q1qqq2z/2RSJKm0bAngljXc/3OwDOA/UYfjiRJu2WUeeo+NFP93pHk\nIcAXgJcA95ophKpqe5KZs/OtBq7ouf22tm2noinJBmADwMEHH7yboUlLx/pDN87r9pvP2TCiSKRd\nG+qYpqr6bs/lW1X198DjxxybJElDGXGeWgE8HHhLVT0M+CHNVLzZZFBIA2LcWFXrqmrdypUrdzM0\nSdIkDDs9r/eg1jvQ7NG7+1htfYu6AAASO0lEQVQikiRpjkacp7YB26rqvHb5dJqi6aqZaXdJVgFX\n9/Q/qOf2BwJX7uZ9S5Km0LDT8/53z/WbgcuA3xl5NJIk7Z6R5amq+naSK5Lcv6q+DhwOXNxe1gMn\ntn8/2N7kDOCFSd5DcwKI6z2eSZKWlmHPnver4w5EkqTdNYY89SLgXe2Z8y4FnkczgvW+JMcBl9Mc\nNwXN6c2fCmwFbmj7SpKWkGGn572sa31V/d1owpEkae5Gnaeq6gJ2PrnEjMMH9C3gBXPZviRpcZnL\n2fMeSTMFAeA3gc+w89mCJEmaFPOUJGlshi2aDgAeXlXfB0jyWuD9VfX74wpMkqQ5ME9JksZmqFOO\nAwcDN/Us3wSsGXk0kiTtHvOUJGlshh1peidwfpJ/ofntid8CTh1bVJIkzY15SpI0NsOePe+vknwE\nOLRtel5VfWl8YUmSNDzzlCRpnIadngewF/C9qnojsC3JIWOKSZKk3WGekiSNxVBFU5LXAMcDr2ib\n7gj807iCkiRpLsxTkqRxGnak6beAI4AfAlTVlcDdxxWUJElzZJ6SJI3NsEXTTe2P9xVAkruOLyRJ\nkubMPCVJGpthi6b3JXkbsE+S5wMfB94+vrAkSZoT85QkaWyGPXve3yZ5IvA94P7A/6yqs8YamSRJ\nQzJPSZLGaZdFU5I9gI9V1RMAE5AkaaqYpyRJ47bL6XlV9VPghiR7L0A8kiTNiXlKkjRuQ03PA24E\nvpzkLNozEwFU1YvHEpUkSXNjnpIkjc2wRdP/ay+SJE0j85QkaWw6i6YkB1fV5VW1eaECkiRpWOYp\nSdJC2NUxTf86cyXJP485FkmS5so8JUkau10VTem5fp9xBiJJ0m4wT0mSxm5XRVPNcl2SpGlgnpIk\njd2uTgTxkCTfo9mTd5f2Ou1yVdU9xhqdJEndzFOSpLHrLJqqao+FCkSSpLkyT0mSFsIuf9xWkiRJ\nkpYziyZJkiRJ6mDRJEmSJEkdLJokSZIkqYNFkyRJkiR1sGiSJEmSpA4WTZIkSZLUwaJJkiRJkjpY\nNEmSJElSB4smSZIkSepg0SRJkiRJHSyaJEmSJKmDRZMkSZIkdbBokiRJkqQOEyuakuyR5EtJPtwu\nH5LkvCSXJHlvkju17Xu2y1vb9WsmFbMkSZKk5WeSI00vAb7as/x64KSqWgtcCxzXth8HXFtV9wVO\navtJkiRJ0oKYSNGU5EDg14F/bJcDPB44ve2yGXhae/3Idpl2/eFtf0mSJEkau0mNNP098HLglnZ5\nf+C6qrq5Xd4GrG6vrwauAGjXX9/2lyRJkqSxW/CiKclvAFdX1Rd6mwd0rSHW9W53Q5ItSbbs2LFj\nBJFKkiRJ0mRGmh4LHJHkMuA9NNPy/h7YJ8mKts+BwJXt9W3AQQDt+r2Ba/o3WlUbq2pdVa1buXLl\neB+BJEmSpGVjxa67jFZVvQJ4BUCSw4A/q6rnJHk/cBRNIbUe+GB7kzPa5c+16z9ZVbcbaZI0WtuP\nnd8+lVWbbtl1J0mSpEVgmn6n6XjgZUm20hyzdHLbfjKwf9v+MuCECcUnSZIkaRla8JGmXlV1NnB2\ne/1S4FED+twIPGNBA5MkSZKk1jSNNEmSJEnS1JnoSJMkSdMqyR7AFuBbVfUbSQ6hOe52P+CLwO9V\n1U1J9gROBR4BfBd4ZlVdNqGwpWVj/aEb53X7zedsGFEkWg4caZIkabCXAF/tWX49cFJVrQWuBY5r\n248Drq2q+wIntf0kSUuIRZMkSX2SHAj8OvCP7XJofiLj9LbLZuBp7fUj22Xa9Ye3/SVJS4RFkyRJ\nt/f3wMuBmXPn7w9cV1U3t8vbgNXt9dXAFQDt+uvb/jvxR9glafHymCZJknok+Q3g6qr6Qvt7ggCD\nRo5qiHW3NVRtBDYCrFu3btH/3uB8jyeRpMXEokmSpJ09FjgiyVOBOwP3oBl52ifJinY06UDgyrb/\nNuAgYFuSFcDewDULH7YkaVycnidJUo+qekVVHVhVa4CjgU9W1XOATwFHtd3WAx9sr5/RLtOu/2RV\nLfqRJEnSbSyaJEkazvHAy5JspTlm6eS2/WRg/7b9ZcAJE4pPkjQmTs+TJGkWVXU2cHZ7/VLgUQP6\n3Ag8Y0EDkyQtKEeaJEmSJKmDRZMkSZIkdbBokiRJkqQOFk2SJEmS1METQUyp7cdaz0qSJEnTwG/m\nkiRJktTBokmSJEmSOlg0SZIkSVIHiyZJkiRJ6mDRJEmSJEkdLJokSZIkqYNFkyRJkiR1sGiSJEmS\npA4WTZIkSZLUwaJJkiRJkjpYNEmSJElSB4smSZIkSeqwYtIBSFqath87v30yqzbdMqJIJEmS5seR\nJkmSJEnqYNEkSZIkSR2cnidJkqRlZ/2hG+e9jc3nbBhBJFoMHGmSJEmSpA4WTZIkSZLUwaJJkiRJ\nkjpYNEmSJElSB4smSZIkSepg0SRJkiRJHSyaJEmSJKmDRZMkSZIkdbBokiRJkqQOFk2SJEmS1MGi\nSZIkSZI6WDRJkiRJUgeLJkmSJEnqYNEkSZIkSR0smiRJkiSpw4pJByBJg2w/dv77dFZtumUEkUiS\npOXOkSZJkiRJ6mDRJEmSJEkdLJokSZIkqcOCF01JDkryqSRfTXJRkpe07fslOSvJJe3ffdv2JHlT\nkq1JLkzy8IWOWZIkSdLyNYmRppuBP62qBwCPAV6Q5IHACcAnqmot8Il2GeApwNr2sgF4y8KHLEmS\nJGm5WvCiqaq2V9UX2+vfB74KrAaOBDa33TYDT2uvHwmcWo1zgX2SrFrgsCVJkiQtUxM9pinJGuBh\nwHnAvapqOzSFFXDPtttq4Iqem21r2/q3tSHJliRbduzYMc6wJUmSJC0jE/udpiR3A/4ZeGlVfS/J\nrF0HtNXtGqo2AhsB1q1bd7v1kiQNI8lBwKnAzwC3ABur6o1J9gPeC6wBLgN+p6quTZPA3gg8FbgB\nOGZmRoWkpW39oRvndfvN52wYUSQat4mMNCW5I03B9K6q+kDbfNXMtLv279Vt+zbgoJ6bHwhcuVCx\nSpKWHY+9lSTtZBJnzwtwMvDVqvq7nlVnAOvb6+uBD/a0P7c9i95jgOtnpvFJkjRqHnsrSeo3iel5\njwV+D/hykgvatlcCJwLvS3IccDnwjHbdmTRTHrbSTHt43sKGK0larrqOvU2yq2Nv3cEnSUvEghdN\nVfVZBh+nBHD4gP4FvGCsQUmS1GfUx94m2UAzfY+DDz54VGFKkhbARM+eJ0nSNBrHsbdVtbGq1lXV\nupUrV44veEnSyFk0SZLUw2NvJUn9JnbKcUmSppTH3kqSdmLRJElSD4+9lST1c3qeJEmSJHWwaJIk\nSZKkDhZNkiRJktTBY5rGZPux1qOSJEnSUuA3e0mSJEnqYNEkSZIkSR0smiRJkiSpg8c0SZK0DK0/\ndOOkQ5CkRcORJkmSJEnqYNEkSZIkSR0smiRJkiSpg0WTJEmSJHWwaJIkSZKkDhZNkiRJktTBokmS\nJEmSOlg0SZIkSVIHf9xW0pK1/dj57RdatemWEUUiSZIWM0eaJEmSJKmDRZMkSZIkdbBokiRJkqQO\nFk2SJEmS1MGiSZIkSZI6WDRJkiRJUgeLJkmSJEnqYNEkSZIkSR38cVtJkiRpAtYfunFet998zoYR\nRaJdcaRJkiRJkjpYNEmSJElSB4smSZIkSepg0SRJkiRJHSyaJEmSJKmDZ8+TpFlsP3Z++5VWbbpl\nRJFIkqRJsmiSJEmSFiFPWb5wnJ4nSZIkSR0smiRJkiSpg0WTJEmSJHXwmCZJGhNPJCFJ0tLgSJMk\nSZIkdXCkaYD57h2WJEmStHRYHUiSJElSB4smSZIkSerg9DxJkhah+f6opSRpeBZNkiRJ0jI0350v\nm8/ZMKJIpp/T8yRJkiSpg0WTJEmSJHVYNNPzkjwZeCOwB/CPVXXihEOSJOlW5ilJy81ymt63KIqm\nJHsA/wA8EdgGfD7JGVV18WQjk6Txme9vxq3adMuIItGumKckaWlbFEUT8Chga1VdCpDkPcCRgMlI\nkmZh0bWgzFOSNEeLaaRqsRzTtBq4omd5W9smSdI0ME9J0hK2WEaaMqCtduqQbABmys0fJPl6z+oD\ngO+MKbZRMs7RMs7RMs7Rmv443xGYX5w/O7pgpt5881S/6X99TH+M0x4fTH+M0x4fTH+M0x4fTH+M\ns8Z3av5gFNsfKlctlqJpG3BQz/KBwJW9HapqIzBwjC/JlqpaN77wRsM4R8s4R8s4R8s4l5x55al+\ni+F5n/YYpz0+mP4Ypz0+mP4Ypz0+mP4YpyW+xTI97/PA2iSHJLkTcDRwxoRjkiRphnlKkpawRTHS\nVFU3J3kh8DGaU7luqqqLJhyWJEmAeUqSlrpFUTQBVNWZwJm7efP5nZpj4RjnaBnnaBnnaBnnEjPP\nPNVvMTzv0x7jtMcH0x/jtMcH0x/jtMcH0x/jVMSXqtp1L0mSJElaphbLMU2SJEmSNBFLqmhK8uQk\nX0+yNckJA9bvmeS97frzkqxZ+CiHivNxSb6Y5OYkR00ixjaOXcX5siQXJ7kwySeSTOT0wkPE+YdJ\nvpzkgiSfTfLAaYyzp99RSSrJRM4UM8TzeUySHe3zeUGS35/GONs+v9O+Ri9K8u5pizHJST3P438m\nuW6hYxwyzoOTfCrJl9r3+1MnEedStBjy1rTnrMWQq6Y9T017floMeWnac9JiyEdTn4uqaklcaA68\n/QZwH+BOwH8AD+zr88fAW9vrRwPvndI41wAPBk4Fjpri5/NXgb3a6380xc/nPXquHwF8dBrjbPvd\nHfgMcC6wbhrjBI4B3rzQse1GnGuBLwH7tsv3nLYY+/q/iObkAdP4XG4E/qi9/kDgskn+/5fKZTHk\nrWnPWYshV017npr2/LQY8tK056TFkI8WQy5aSiNNjwK2VtWlVXUT8B7gyL4+RwKb2+unA4cnGfSD\nhOO0yzir6rKquhC4ZYFj6zVMnJ+qqhvaxXNpfpdkoQ0T5/d6Fu9K3w9OLpBhXp8AfwG8AbhxIYPr\nMWyckzZMnM8H/qGqrgWoqqunMMZezwJOW5DIdjZMnAXco72+N32/P6Tdthjy1rTnrMWQq6Y9T017\nfloMeWnac9JiyEdTn4uWUtG0GriiZ3lb2zawT1XdDFwP7L8g0Q2IoTUozmkw1ziPAz4y1ogGGyrO\nJC9I8g2aD/wXL1BsvXYZZ5KHAQdV1YcXMrA+w/7ff7sdGj89yUED1o/bMHHeD7hfkn9Pcm6SJy9Y\ndI2h30PtdKFDgE8uQFz9honztcDvJtlGc3a4Fy1MaEveYshb056zFkOumvY8Ne35aTHkpWnPSYsh\nH019LlpKRdOgPW/9e2qG6TNu0xDDMIaOM8nvAuuAvxlrRIMNFWdV/UNV/RxwPPDqsUd1e51xJrkD\ncBLwpwsW0WDDPJ8fAtZU1YOBj3PbXvCFNEycK2imQxxGs9fsH5PsM+a4es3lvX40cHpV/XSM8cxm\nmDifBZxSVQcCTwXe2b5mNT+LIW9N+v53ZTHkqmnPU9OenxZDXpr2nLQY8tHU56KllPS2Ab17Fg7k\n9sN2t/ZJsoJmaO+aBYluQAytQXFOg6HiTPIE4FXAEVX14wWKrddcn8/3AE8ba0SD7SrOuwMPAs5O\nchnwGOCMhT7YliGez6r6bs//+u3AIxYotl7Dvt8/WFU/qapvAl+nSVgLZS6vzaOZzNQ8GC7O44D3\nAVTV54A7AwcsSHRL22LIW9OesxZDrpr2PDXt+Wkx5KVpz0mLIR9Nfy5ayAOoxnmhqeAvpRlSnDmA\n7Bf6+ryAnQ+ofd80xtnT9xQmdyKIYZ7Ph9EctLd2yv/va3uu/yawZRrj7Ot/NpM5EcQwz+eqnuu/\nBZw7pXE+GdjcXj+AZth//2mKse13f+Ay2t/Nm9Ln8iPAMe31B9AksonEu5QuiyFvTXvOWgy5atrz\n1LTnp8WQl6Y9Jy2GfLQYctGCPiEL8IQ/FfjP9sPxVW3b62j2LEFTkb4f2AqcD9xnSuN8JE3F/UPg\nu8BFUxrnx4GrgAvayxlTGucbgYvaGD/VlQwmGWdf3wVNSnN8Pv+6fT7/o30+f35K4wzwd8DFwJeB\no6ctxnb5tcCJk3gO5/BcPhD49/Z/fgHwpEnGu5QuiyFvTXvOWgy5atrz1LTnp8WQl6Y9Jy2GfDTt\nuShtEJIkSZKkAZbSMU2SJEmSNHIWTZIkSZLUwaJJkiRJkjpYNEmSJElSB4smSZIkSepg0SRJkiRJ\nHSyaJEmSJKmDRZMkSZIkdfj/oIUcqN7x43kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(1346269)\n",
    "\n",
    "fig = plt.figure(figsize=(14, 6))\n",
    "\n",
    "p_given_x_3 = p[np.where(x==3)]\n",
    "\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax1.hist(p_given_x_3, bins=20, color='#e66101')\n",
    "ax1.set_title(r'Histogram of $\\tt{p[numpy.where(x==3)]}$')\n",
    "ax1.set_ylabel('Frequency')\n",
    "\n",
    "brv = beta.rvs(4.0, 8.0, size=10**4)\n",
    "\n",
    "ax2 = fig.add_subplot(122)\n",
    "ax2.hist(brv, bins=20, color='#5e3c99')\n",
    "ax2.set_title(r'Histogram of $Beta(4, \\, 8)$, $10^4$ simulated values')\n",
    "ax2.set_ylabel('Frequency')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the side-by-side graph above, we obtain an instance of [`matplotlib.pyplot.figure`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.figure.html) from `pyplot`, specifying the figure's width and height. Next, we use [`figure.add_subplot(121)`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.subplot.html) to create graph axes in a figure with 1 row and 2 columns, in the 1st column position from the left, for the histogram of `p[numpy.where(x==3)]`. Likewise, `figure.add_subplot(122)` creates the graph axes for the histogram for $Beta(4, 8)$ in the 2nd column position of the 1-row, 2-column figure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating order statistics\n",
    "\n",
    "Simulating order statistics in Python/NumPy/SciPy is easy: we simply simulate i.i.d. r.v.s and sort them in order. For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-2.05800632 -1.14300612 -0.09108752 -0.08188422  0.45810851  0.55556877\n",
      "  0.59398037  0.9312065   0.98367818  1.12202827]\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(2178309)\n",
    "\n",
    "from scipy.stats import norm\n",
    "\n",
    "rv = np.sort(norm.rvs(size=10))\n",
    "print(rv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "produces one realization of $X_{(1)}, \\ldots, X_{(10)}$, where $X_{(1)}, \\ldots, X_{(10)}$ are i.i.d. $N(0, 1)$. If we want to plot a histogram of realizations of, say, $X_{(9)}$, we'll need to use iteration. Let's try using [Python's list comprehensions](https://docs.python.org/3/tutorial/datastructures.html#list-comprehensions):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "order_stats has shape: (10000, 10)\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(3524578)\n",
    "\n",
    "rv_arr = [np.sort(norm.rvs(size=10)) for _ in range(10**4)]\n",
    "order_stats = np.array(rv_arr)\n",
    "\n",
    "print('order_stats has shape: {}'.format(order_stats.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 10<sup>4</sup> iterations are done with\n",
    "\n",
    "    rv_arr = [np.sort(norm.rvs(size=10)) for _ in range(10**4)]\n",
    "\n",
    "This Python list comprehension does the following:\n",
    "* loops 10<sup>4</sup> times with `for _ in range(10**4)`, where we ignore the iteration counter by using `_`\n",
    "   * as we loop, we do not need to know which particular iteration we are on\n",
    "* for each iteration, calculate `np.sort(norm.rvs(size=10))`\n",
    "* gather up the results of all the iterations into a Python list (array), assigning this resulting list to `rv_arr`\n",
    "* this one-liner list comprehension is equivalent to all this code:\n",
    "<pre>\n",
    "    rv_arr = []\n",
    "    for _ in range(10**4):\n",
    "        rv = np.sort(norm.rvs(size=10))\n",
    "        rv_arr.append(rv)\n",
    "</pre>\n",
    "\n",
    "This creates a matrix, `order_stats`, with 10<sup>4</sup> rows and 10 columns. The _i_<sup>th</sup> column of the matrix contains 10<sup>4</sup> realizations of $X_{(i+1)}$ (don't forget that Python is zero-indexed, so the 8<sup>th</sup> column is $X_{(9)}$!). Now we can create a histogram of $X_{(9)}$, simply by selecting column 8 of the matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAGFlJREFUeJzt3Xu0XWV97vHvY7h5D0hQTIJBG1HK\n0Uoj4tG2CqhclKCVFo5HKaWNHYKXY1uJl3GwF8+J51ip1IqNhQqKIqIWFCxGLqU9Qy4REIFgiZjC\nJgixXLwgpIHf+WPNrSvJzs6aO2vttffO9zPGGnvNd75rzd9eg/Ds933nnCtVhSRJvXrcsAuQJE0v\nBockqRWDQ5LUisEhSWrF4JAktWJwSJJaMTgkSa0YHJpRktyc5BXDrmNYkuyT5PokP0nyjmHXo5nJ\n4NC0kWRNkkM2afu9JP86ul1Vv1pVV7R9nxnkPcAVVfXkqjpt051JnpPkZ0n27Gp7U5K1SeZPaqWa\ntgwOqY+S7DDkEp4F3LylnVX1feBrwLsAkrwU+DhwVFXdOSkVatozODSjdI8mkpyc5K5m2uZ7SQ5O\n8hlgL+CrSX6a5D1N3+cnuSLJA81015Fd77l/1/TPF5N8IclfbnLMk5PcCPwsyQ5Jlib5fvOaW5K8\nfpP+f5rkxuav/zOSPD3J15v+30yy6zi/45i1JrkMeCXw8eZ3e+4W3uLDwFuT7Ad8Gfijqrpmgh+5\ntkMGh2akJPsAJwEvrqonA68B1lTVm4E7gNdV1ZOq6v8k2RH4KvANYA/g7cA5zXrBTsBXgE8DuwGf\nB16/2QHhWOAIYHZVbQC+D/wG8FTgz4DPdk8PAb8NvAp4LvA64OvA+4Dd6fy7HHN9Yrxaq+og4F+A\nk5rf7d/Geo+qug64BrgaOL2qvrDlT1LanMGh6eYfm7+0H0jyAPCJLfR7FNgZ2DfJjlW1ppmmGcuB\nwJOAZVW1vqouozOdc2yzbwfgtKr6z6r6Mp3/6W7qtKq6s6p+DlBVX6yqtVX1WPM/5tuAA7r6/01V\n3VNVd9H5n/3VVXV9VT1CJ6heNIFae5LkcXQ+n8fojD669x2R5PeTzEpyTpLLk5w5OgWXZH6SU3s9\nlmYmg0PTzVFVNXv0AbxtrE5VtZrOPP4HgXuTnJvkmVt4z2cCd1bVY11t/w7MbfbdVRvfRnqstYCN\n2pK8JckNXQG3H53RxKh7up7/fIztJ02g1l79FTCbTpi9aZN9S4DP0RlV3V5VrwRuBd4A0KyDPCPJ\n7BbH0wxjcGjGqqrPVdXL6SwYF7/863rT7xJYC8xv/hIftRdwF3A3MDdJuvaNdfbRL94zybOAT9GZ\nKntaE3A3ARnjdW2NV+tWJXkrnVA4is7n8aejv1sTBo+rqoeB5wA3NC+7js6026h/oTP1p+2UwaEZ\nqVmfOCjJzsDDdP6Kf7TZfQ/w7K7uVwM/A96TZMfmOpDXAecC32ped1Kz6L2YjaecxvJEOkGyrqnl\neDojjn4Yr9ZxNScN/C866zv3AOcDOwGLmy7PpTN6AbgFOKh5fgjQvVh/O7Dvtv0ams4MDs1UOwPL\ngB8BP6SzkPy+Zt//Bj7QTCP9SVWtB44EDmv6fwJ4S1Xd2ux7A3AC8ADw3+msKTyypQNX1S10poO+\nRSek/gvw//rxS41X63ivS/I8OuHy5qr6bvNejwIfBU7uPkTz82vAw82ZWk9k46m0foycNI3FbwCU\n2klyNfDJqvqHYdfST81U1Wer6rWbtH8QuKyqrmy23wbcV1VbHeVoZnLEIW1Fkt9K8oxmquo44AXA\nPw27rn6rqgeAx5Ls0vy+VyS5FFg/GhqN3wQuGU6VmgqGfZWrNB3sA5xH50yn7wNvrKq7h1vSwPwd\n8KaqOgN4xaY7m9uS/LCq7p/swjR1OFUlSWrFqSpJUisGhySplRm5xrH77rvXggULhl2GJE0r3/72\nt39UVXO21m9GBseCBQtYuXLlsMuQpGklyb9vvZdTVZKklgwOSVIrBockqRWDQ5LUisEhSWrF4JAk\ntTKw4Gi+bvLeJDeNse9PklSS3ZvtJDktyeokNybZv6vvcUluax7HDapeSVJvBjni+DRw6KaNzU3S\nXgXc0dV8GLCweSwBTm/67gacAryEzpfnnJJkVyRJQzOw4Ghuw3zfGLtOBd7Dxl/fuRg4uzquAmYn\n2ZPO11OuqKr7mrtxrmCMMJIkTZ5JvXI8yZHAXVX1nY2/wpm5wJ1d2yNN25bax3rvJXRGK+y11159\nrFrbkwVLL5q0Y61ZdsSkHUvqp0kLjiRPAN4PvHqs3WO01TjtmzdWLQeWAyxatMh7xWvKm2hIGTga\ntsk8q+o5wN7Ad5KsAeYB1yV5Bp2RxPyuvvOAteO0S5KGZNKCo6q+W1V7VNWCqlpAJxT2r6ofAhcC\nb2nOrjoQeLD5hrVLgFcn2bVZFH81fmWlJA3VIE/H/TzwLWCfJCNJThin+8XA7cBq4FPA2wCq6j7g\nL4Brm8efN22SpCEZ2BpHVR27lf0Lup4XcOIW+p0JnNnX4iRJE+aV45KkVgwOSVIrBockqRWDQ5LU\nisEhSWrF4JAktWJwSJJaMTgkSa0YHJKkVgwOSVIrBockqRWDQ5LUisEhSWrF4JAktWJwSJJaMTgk\nSa0YHJKkVgwOSVIrBockqRWDQ5LUysCCI8mZSe5NclNX2/9NcmuSG5N8Jcnsrn3vTbI6yfeSvKar\n/dCmbXWSpYOqV5LUm0GOOD4NHLpJ2wpgv6p6AfBvwHsBkuwLHAP8avOaTySZlWQW8LfAYcC+wLFN\nX0nSkAwsOKrqSuC+Tdq+UVUbms2rgHnN88XAuVX1SFX9AFgNHNA8VlfV7VW1Hji36StJGpJhrnH8\nPvD15vlc4M6ufSNN25baN5NkSZKVSVauW7duAOVKkmBIwZHk/cAG4JzRpjG61TjtmzdWLa+qRVW1\naM6cOf0pVJK0mR0m+4BJjgNeCxxcVaMhMALM7+o2D1jbPN9SuyRpCCZ1xJHkUOBk4Miqeqhr14XA\nMUl2TrI3sBC4BrgWWJhk7yQ70VlAv3Aya5YkbWxgI44knwdeAeyeZAQ4hc5ZVDsDK5IAXFVVf1RV\nNyc5D7iFzhTWiVX1aPM+JwGXALOAM6vq5kHVLEnauoEFR1UdO0bzGeP0/xDwoTHaLwYu7mNpkqRt\n4JXjkqRWDA5JUisGhySpFYNDktSKwSFJasXgkCS1MulXjkvaNguWXjSh161ZdkSfK9H2yhGHJKkV\ng0OS1IrBIUlqxeCQJLVicEiSWjE4JEmtGBySpFYMDklSKwaHJKkVrxzXjDTRq6slbZ0jDklSKwaH\nJKkVg0OS1IrBIUlqZWDBkeTMJPcmuamrbbckK5Lc1vzctWlPktOSrE5yY5L9u15zXNP/tiTHDape\nSVJvBjni+DRw6CZtS4FLq2ohcGmzDXAYsLB5LAFOh07QAKcALwEOAE4ZDRtJ0nAMLDiq6krgvk2a\nFwNnNc/PAo7qaj+7Oq4CZifZE3gNsKKq7quq+4EVbB5GkqRJNNlrHE+vqrsBmp97NO1zgTu7+o00\nbVtq30ySJUlWJlm5bt26vhcuSeqYKovjGaOtxmnfvLFqeVUtqqpFc+bM6WtxkqRfmuzguKeZgqL5\neW/TPgLM7+o3D1g7TrskaUgmOzguBEbPjDoOuKCr/S3N2VUHAg82U1mXAK9OsmuzKP7qpk2SNCQD\nu1dVks8DrwB2TzJC5+yoZcB5SU4A7gCObrpfDBwOrAYeAo4HqKr7kvwFcG3T78+ratMFd0nSJBpY\ncFTVsVvYdfAYfQs4cQvvcyZwZh9LkyRtg6myOC5JmiYMDklSKwaHJKkVg0OS1IrBIUlqxeCQJLVi\ncEiSWjE4JEmtGBySpFZ6Co4k+w26EEnS9NDriOOTSa5J8rYkswdakSRpSuspOKrq5cCb6NzifGWS\nzyV51UArkyRNST2vcVTVbcAHgJOB3wJOS3JrkjcMqjhJ0tTT6xrHC5KcCqwCDgJeV1XPb56fOsD6\nJElTTK+3Vf848CngfVX189HGqlqb5AMDqUySNCX1GhyHAz+vqkcBkjwO2KWqHqqqzwysOknSlNNr\ncHwTOAT4abP9BOAbwH8dRFGS+m/B0osm9Lo1y47ocyWa7npdHN+lqkZDg+b5EwZTkiRpKus1OH6W\nZP/RjSS/Dvx8nP6SpBmq16mqdwFfTLK22d4T+N3BlCRJmsp6Co6qujbJ84B9gAC3VtV/TvSgSf4H\n8AdAAd8FjqcTRucCuwHXAW+uqvVJdgbOBn4d+A/gd6tqzUSPLUnaNm1ucvhi4AXAi4Bjk7xlIgdM\nMhd4B7CoqvYDZgHHAB8GTq2qhcD9wAnNS04A7q+qX6FzzciHJ3JcSVJ/9HoB4GeAjwAvpxMgLwYW\nbcNxdwAen2QHOovsd9O5mPD8Zv9ZwFHN88XNNs3+g5NkG44tSdoGva5xLAL2rara1gNW1V1JPgLc\nQWeB/RvAt4EHqmpD020EmNs8nwvc2bx2Q5IHgacBP+p+3yRLgCUAe+2117aWKUnagl6nqm4CntGP\nAybZlc4oYm/gmcATgcPG6DoaUmONLjYLsKpaXlWLqmrRnDlz+lGqJGkMvY44dgduSXIN8MhoY1Ud\nOYFjHgL8oKrWAST5Mp0LCWcn2aEZdcwDRs/gGqFzV96RZmrrqcB9EziuJKkPeg2OD/bxmHcAByZ5\nAp2pqoOBlcDlwBvpnFl1HHBB0//CZvtbzf7L+jFlJkmamF5Px/3nJM8CFlbVN5v/6c+ayAGr6uok\n59M55XYDcD2wHLgIODfJXzZtZzQvOQP4TJLVdEYax0zkuJKk/ugpOJL8IZ2F592A59BZsP4kndFC\na1V1CnDKJs23AweM0fdh4OiJHEeS1H+9Lo6fCLwM+DH84kud9hhUUZKkqavX4HikqtaPbjSL1K4z\nSNJ2qNfF8X9O8j46F+29Cngb8NXBlSV1TPRW4JIGp9cRx1JgHZ37Sr0VuJjO949LkrYzvZ5V9Rid\nr4791GDLkSRNdb2eVfUDxr5a+9l9r0iSNKW1uVfVqF3onB67W//LkSRNdT2tcVTVf3Q97qqqv6Zz\nN1tJ0nam16mq/bs2H0dnBPLkgVQkSZrSep2q+quu5xuANcDv9L0aSdKU1+tZVa8cdCGSpOmh16mq\nd4+3v6o+2p9yJElTXZuzql5M5xbnAK8DrqT5Zj5J0vajzRc57V9VPwFI8kHgi1X1B4MqTJI0NfV6\ny5G9gPVd2+uBBX2vRpI05fU64vgMcE2Sr9C5gvz1wNkDq0qSNGX1elbVh5J8HfiNpun4qrp+cGVJ\nkqaqXqeqAJ4A/LiqPgaMJNl7QDVJkqawnoIjySnAycB7m6Ydgc8OqihJ0tTV64jj9cCRwM8Aqmot\n3nJEkrZLvQbH+qoqmlurJ3ni4EqSJE1lvQbHeUn+Dpid5A+Bb7INX+qUZHaS85PcmmRVkpcm2S3J\niiS3NT93bfomyWlJVie5cZMbLkqSJlmvt1X/CHA+8CVgH+B/VtXfbMNxPwb8U1U9D3ghsIrO19Ne\nWlULgUubbYDDgIXNYwlw+jYcV5K0jbZ6Om6SWcAlVXUIsGJbD5jkKcBvAr8HUFXrgfVJFgOvaLqd\nBVxBZ0F+MXB2M1V2VTNa2bOq7t7WWiRJ7W11xFFVjwIPJXlqn475bGAd8A9Jrk/y982aydNHw6D5\nuUfTfy4b3xNrpGnbSJIlSVYmWblu3bo+lSpJ2lSvV44/DHw3yQqaM6sAquodEzzm/sDbq+rqJB/j\nl9NSY8kYbWN9//lyYDnAokWLNtsvSeqPXoPjoubRDyPASFVd3WyfTyc47hmdgkqyJ3BvV//5Xa+f\nB6ztUy2SpJbGDY4ke1XVHVV1Vr8OWFU/THJnkn2q6nvAwcAtzeM4YFnz84LmJRcCJyU5F3gJ8KDr\nG5I0PFsbcfwjnWklknypqn67T8d9O3BOkp2A24Hj6ay3nJfkBOAO4Oim78XA4cBq4KGmryRpSLYW\nHN3rC8/u10Gr6gY6Xw61qYPH6FvAif06tiRp22wtOGoLzyVtJxYsbb+8uWbZEQOoRFPF1oLjhUl+\nTGfk8fjmOc12VdVTBlqdJGnKGTc4qmrWZBUiSZoe2nwfhyRJBockqR2DQ5LUisEhSWrF4JAktWJw\nSJJaMTgkSa0YHJKkVgwOSVIrBockqRWDQ5LUisEhSWrF4JAktWJwSJJaMTgkSa0YHJKkVgwOSVIr\nBockqZWhBUeSWUmuT/K1ZnvvJFcnuS3JF5Ls1LTv3GyvbvYvGFbNkqThjjjeCazq2v4wcGpVLQTu\nB05o2k8A7q+qXwFObfpJkoZkh2EcNMk84AjgQ8C7kwQ4CPhvTZezgA8CpwOLm+cA5wMfT5Kqqsms\nWdtuwdKLhl2CpD4Y1ojjr4H3AI81208DHqiqDc32CDC3eT4XuBOg2f9g038jSZYkWZlk5bp16wZZ\nuyRt1yY9OJK8Fri3qr7d3TxG1+ph3y8bqpZX1aKqWjRnzpw+VCpJGsswpqpeBhyZ5HBgF+ApdEYg\ns5Ps0Iwq5gFrm/4jwHxgJMkOwFOB+ya/bEkSDGHEUVXvrap5VbUAOAa4rKreBFwOvLHpdhxwQfP8\nwmabZv9lrm9I0vBMpes4TqazUL6azhrGGU37GcDTmvZ3A0uHVJ8kiSGdVTWqqq4Armie3w4cMEaf\nh4GjJ7UwSdIWTaURhyRpGjA4JEmtDHWqStLMNNGLPdcsO6LPlWgQHHFIkloxOCRJrRgckqRWDA5J\nUisGhySpFYNDktSKwSFJasXgkCS1YnBIkloxOCRJrRgckqRWDA5JUisGhySpFYNDktSKt1VXaxO9\nZbakmcERhySpFYNDktTKpAdHkvlJLk+yKsnNSd7ZtO+WZEWS25qfuzbtSXJaktVJbkyy/2TXLEn6\npWGMODYAf1xVzwcOBE5Msi+wFLi0qhYClzbbAIcBC5vHEuD0yS9ZkjRq0oOjqu6uquua5z8BVgFz\ngcXAWU23s4CjmueLgbOr4ypgdpI9J7lsSVJjqGscSRYALwKuBp5eVXdDJ1yAPZpuc4E7u1420rRt\n+l5LkqxMsnLdunWDLFuStmtDC44kTwK+BLyrqn48Xtcx2mqzhqrlVbWoqhbNmTOnX2VKkjYxlOBI\nsiOd0Dinqr7cNN8zOgXV/Ly3aR8B5ne9fB6wdrJqlSRtbNIvAEwS4AxgVVV9tGvXhcBxwLLm5wVd\n7SclORd4CfDg6JSWpJlloheXrll2RJ8r0XiGceX4y4A3A99NckPT9j46gXFekhOAO4Cjm30XA4cD\nq4GHgOMnt1xJUrdJD46q+lfGXrcAOHiM/gWcONCiJEk988pxSVIrBockqRWDQ5LUisEhSWrF4JAk\ntWJwSJJaMTgkSa0YHJKkVvzO8e2Y3x0uaSIMDknTnve4mlxOVUmSWjE4JEmtGBySpFYMDklSKwaH\nJKkVg0OS1IrBIUlqxes4JG23JnL9h9d+OOKQJLVkcEiSWnGqagbwnlOSJtO0GXEkOTTJ95KsTrJ0\n2PVI0vZqWow4kswC/hZ4FTACXJvkwqq6ZbiV9Z+jB2lqm+x/o1NxMX5aBAdwALC6qm4HSHIusBiY\nccEhSd2m4p1/p0twzAXu7NoeAV7S3SHJEmBJs/lIkpsmqbbtwe7Aj4ZdxAzhZ9k/fpbjyIdbdR/9\nLJ/VS+fpEhwZo6022qhaDiwHSLKyqhZNRmHbAz/P/vGz7B8/y/5p+1lOl8XxEWB+1/Y8YO2QapGk\n7dp0CY5rgYVJ9k6yE3AMcOGQa5Kk7dK0mKqqqg1JTgIuAWYBZ1bVzeO8ZPnkVLbd8PPsHz/L/vGz\n7J9Wn2Wqauu9JElqTJepKknSFGFwSJJambHBkeToJDcneSyJp+xNgLd56Z8kZya51+uLtk2S+Uku\nT7Kq+ff9zmHXNJ0l2SXJNUm+03yef9bL62ZscAA3AW8Arhx2IdNR121eDgP2BY5Nsu9wq5rWPg0c\nOuwiZoANwB9X1fOBA4ET/e9ymzwCHFRVLwR+DTg0yYFbe9GMDY6qWlVV3xt2HdPYL27zUlXrgdHb\nvGgCqupK4L5h1zHdVdXdVXVd8/wnwCo6d5bQBFTHT5vNHZvHVs+YmrHBoW021m1e/AeqKSPJAuBF\nwNXDrWR6SzIryQ3AvcCKqtrq5zktruPYkiTfBJ4xxq73V9UFk13PDLPV27xIw5LkScCXgHdV1Y+H\nXc90VlWPAr+WZDbwlST7VdW4a3HTOjiq6pBh1zCDeZsXTUlJdqQTGudU1ZeHXc9MUVUPJLmCzlrc\nuMHhVJW2xNu8aMpJEuAMYFVVfXTY9Ux3SeY0Iw2SPB44BLh1a6+bscGR5PVJRoCXAhcluWTYNU0n\nVbUBGL3NyyrgvK3c5kXjSPJ54FvAPklGkpww7JqmqZcBbwYOSnJD8zh82EVNY3sClye5kc4fiyuq\n6mtbe5G3HJEktTJjRxySpMEwOCRJrRgckqRWDA5JUisGhySpFYNDktSKwSFJauX/A83eHogaeYn9\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x9 = order_stats[:, 8]\n",
    "\n",
    "plt.hist(x9, bins=20)\n",
    "\n",
    "plt.xticks(np.arange(-1, 4, 1))\n",
    "plt.title(r'Histogram of $X_{(9)}$')\n",
    "plt.ylabel('Frequency')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compute summaries like `numpy.mean(x9)` and `numpy.var(x9)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean of x9 = 1.0092010380841254\n",
      "variance of x9 = 0.2108867888683199\n"
     ]
    }
   ],
   "source": [
    "mean_x9 = np.mean(x9)\n",
    "print('mean of x9 = {}'.format(mean_x9))\n",
    "\n",
    "var_x9 = np.var(x9)\n",
    "print('variance of x9 = {}'.format(var_x9))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "\n",
    "Joseph K. Blitzstein and Jessica Hwang, Harvard University and Stanford University, &copy; 2019 by Taylor and Francis Group, LLC"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  },
  "notebook_info": {
   "author": "Joseph K. Blitzstein, Jessica Hwang",
   "chapter": "8. Transformations",
   "creator": "buruzaemon",
   "section": "8.8",
   "title": "Introduction to Probability, 1st Edition"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}