{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 9: Conditional Expectation\n",
    " \n",
    "This Jupyter notebook is the Python equivalent of the R code in section 9.8 R, pp. 408 - 410, [Introduction to Probability, 1st Edition](https://www.crcpress.com/Introduction-to-Probability/Blitzstein-Hwang/p/book/9781466575578), 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": [
    "## Mystery prize simulation\n",
    "\n",
    "We can use simulation to show that in Example 9.1.7, the example of bidding on a mystery prize with unknown value, any bid will lead to a negative payout on average. First choose a bid `b` (we chose 0.6); then simulate a large number of hypothetical mystery prizes using `rvs` in [`scipy.stats.uniform`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.uniform.html), and store them in `v`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# seed the random number generator\n",
    "np.random.seed(5702887)\n",
    "\n",
    "from scipy.stats import uniform\n",
    "\n",
    "b = 0.6\n",
    "nsim = 10**5\n",
    "v = uniform.rvs(size=nsim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bid is accepted if `b > (2/3)*v`. To get the average profit conditional on an accepted bid, we use [`numpy.where`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.where.html) to find the indices in `v` where the values satisfy the condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "average profit conditional on an accepted bid: -0.1499047321616897\n"
     ]
    }
   ],
   "source": [
    "bid_accepted = v[np.where(b > (2/3)*v)]\n",
    "\n",
    "avg_profit = np.mean(bid_accepted) - b\n",
    "\n",
    "print('average profit conditional on an accepted bid: {}'.format(avg_profit))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This value is negative regardless of `b`, as you can check by experimenting with different values of`b`.  Try changing the value for `b` in the code block above, and then hit SHIFT+ENTER to re-run the code blocks. You can also try changing the seed value for random number generation as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time until HH vs. HT\n",
    "\n",
    "To verify the results of Example 9.1.9, we can start by generating a long sequence of fair coin tosses. This is done with the [`numpy.random.choice`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.choice.html) function. `numpy.random.choice` will chooose with replacement items from the specified list of `H` and `T`, returning a sequence of length specified bye the `size` parameter. We use [`str.join`](https://docs.python.org/3.7/library/string.html), passing in an empty `str` argument `''` to be interspersed between the elements in the given sequence. This results in a single string of `H`'s and `T`'s:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'HHHHTHHTTHTHTHTHTTHHHHTTHTTHHHHTHHHHHTHHHHHTTHHHHHHTTHHTHHHTTHTHHHHTTHHHTHTHTTHHTHTTHHHTTTHHTHTTTHHH'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "''.join(np.random.choice(['H', 'T'], size=100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A sequence of length 100 is enough to virtually guarantee that both `HH` and `HT` will have appeared at least once.\n",
    "\n",
    "To determine how many tosses are required on average to see `HH` and `HT`, we need to generate many sequences of coin tosses. For this, we use our familiar friend [Python list comprehensions](https://docs.python.org/3/tutorial/datastructures.html#list-comprehensions) to execute a large number of iterations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(9227465)\n",
    "\n",
    "r = [''.join(np.random.choice(['H', 'T'], size=100)) for _ in range(10**3)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now `r` contains a thousand sequences of coin tosses, each of length 100. To find the first appearance of `HH` in each of these sequences, it is easiest to use [regular expressions](https://docs.python.org/3.7/howto/regex.html) in Python's [`re`](https://docs.python.org/3.7/library/re.html) library.  We load the `re` library with the `import` keyword, and use `re.search` to obtain the start and end indices of the first appearance of `HH` in each sequence `ht_seq` in `r`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix t has shape: (1000, 2)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[1, 3],\n",
       "       [1, 3],\n",
       "       [4, 6],\n",
       "       [2, 4],\n",
       "       [3, 5],\n",
       "       [0, 2],\n",
       "       [4, 6],\n",
       "       [0, 2],\n",
       "       [5, 7],\n",
       "       [1, 3]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import re\n",
    "\n",
    "# to learn more about re, un-comment ouf the following line\n",
    "#print(re.__doc__)\n",
    "\n",
    "target = 'HH'\n",
    "\n",
    "t = [re.search(target, ht_seq).span() for ht_seq in r]\n",
    "t = np.array(t)\n",
    "\n",
    "print('matrix t has shape: {}'.format(t.shape))\n",
    "\n",
    "t[0:10, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code above creates a two-column matrix `t`, whose columns contain the starting and ending positions of the first appearance of `HH` in each sequence of coin tosses. (Use `t[0:10, :]` to display the first 10 rows of the matrix and get an idea of what your results look like.) What we want are the ending positions, given by the second column. In particular, we want the average value of the second column, which is an approximation of the average waiting time for `HH`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "average waiting time for 'HH': 6.16\n"
     ]
    }
   ],
   "source": [
    "mean = t[:, 1].mean()\n",
    "\n",
    "print('average waiting time for \\'{}\\': {}'.format(target, mean))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is your answer around 6? Trying again with `HT` instead of `HH`, is your answer around 4? You can change the value of the `target` variable in the earlier code block, and then hit SHIFT+ENTER to execute the code blocks again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Example 9.3.10, we derived formulas for the slope and intercept of a linear regression model, which can be used to predict a response variable using an explanatory variable. Let's try to apply these formulas to a simulated dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(14930352)\n",
    "\n",
    "from scipy.stats import norm\n",
    "\n",
    "x = norm.rvs(size=100)\n",
    "y = 3 + 5*x + norm.rvs(size=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The array `x` contains 100 realizations of the random variable $X \\sim N(0, 1)$ and the array `y` contains 100 realizations of the random variable $Y = a + bX + \\epsilon$, where $\\epsilon \\sim N(0,1)$. As we can see, the true values of $a$ and $b$ for this dataset are 3 and 5, respectively. We can visualize the data as a scatterplot with [`matplotlib.pyplot.scatter(x,y)`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.scatter.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b = 4.869568946916432\n",
      "a = 2.9953152340248232\n"
     ]
    }
   ],
   "source": [
    "# numpy.cov(x, y) returns a 2 x 2 covariance matrix\n",
    "#   cov(x,x)    cov(x,y)\n",
    "#   cov(y,x)    cov(y,y)\n",
    "cov_xy = np.cov(x, y, ddof=1)[0][1]\n",
    "var_x = np.var(x, ddof=1)\n",
    "\n",
    "b = cov_xy / var_x\n",
    "a = np.mean(y) - b*np.mean(x)\n",
    "\n",
    "print('b = {}'.format(b))\n",
    "print('a = {}'.format(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here `numpy.cov(x, y, ddof=1)[0][1]`, `numpy.var(x, ddof=1)`, and `numpy.mean(x)` provide the sample covariance, sample variance, and sample mean, estimating the quantities $Cov(X, Y), Var(X)$, and $\\mathbb{E}(X)$, respectively. (We have discussed sample mean and sample variance in detail in earlier chapters. Sample covariance is defined analogously, and is a natural way to estimate the true covariance.)\n",
    "\n",
    "You should find that `b` is close to 5 and `a` is close to 3. These estimated values define the line of best fit. We use yet another list comprehension to calculate `y` values corresponding to `x` on the best-fit line using `a` and `b`, and then `matplotlit.pyplot.plot` lets use render the _line of best fit_ on top of our scatterplot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XucFOWZL/DfU93DDDIDiCCoSACRwUvMSTLxRNwkRrwQE1dxI144u9HkBNmsJyefT7IhOcTclF2zl+xn/ZgT9SRmTRYV8BYSLygqMQaSiH50vQyDExwNIgOIMDPKDHTXc/7orp7q7qrq6u6qruru3/cf6Le7a15bpp963+d9n1dUFURERG6MqDtARETxxkBBRESeGCiIiMgTAwUREXlioCAiIk8MFERE5ImBgoiIPDFQEBGRJwYKIiLylIy6A0GYPHmyzpw5M+puEBHVlWeffXavqk4p9bqGCBQzZ87Eli1bou4GEVFdEZHX/byOU09EROSJgYKIiDwxUBARkScGCiIi8sRAQUREnhgoiIjIEwMFERF5YqAgIqozur8f6VUrYG78BfTwSOg/ryE23BERNQNVE+ad140+fnMrZOQ9oKU11J/LQEFEVAfMZ34F3fb7vDY54cOQ9iND/9kMFEREMaZD78D85b8UtRsX/C/IkdNq0gcGCiKimEqvWlHcOGYsEpd+q6b9YKAgIooZfeMlmL+9q6jdWLQc2waBTZv7MDiSQkdrEvNnT0Ln1PGh9oeBgogoRrxGET39A3i8Zw9SpgIABkdSeLxnDwCEGiwYKIiIYsD89b9DD+wuajcu/x4kkfmq3rR9Xy5IWFKmYuO2PQwURESNSg8dhLn2hqJ2mTcfxoc/ndc2OJJyvMZIWtHTPxBasGCgICKKiH2aafhwGkMjaZiqeOQD12L+9EnoLHh9R2vSNVhs2r6PgYKIqFFo/2swN/wk93j4cBqDwyn84ZhP4vUJnYBL7mH+7ElY3108PQW4jzaCwEBBRFRDTsnqoZE01sxblteWMrVolNA5dTx+8+peDKfMomt0tIb3dc5AQURNr6d/AJu27wt1yam55UFoz6aiduOir2H1lrcd3+M0SvjEiZPzVj4BQNIQzJ89KbjOFmCgIKKmFvaS08L6THaJJSsBAB2tBxyDgtMowepT2IHNjoGCiJqa25LTIJLDjnsiABhXXg+R0eLd82dPch0luI12wt5kZ8dAQURNzS0JXE1yWAffhrnuh0XtMm0OjAVXF7W7jRIARLLBrhADBRE1Nbclp5Umh91GEdY0kxunUcLPNveFNtopBwMFETU1r2mfcpgvPAZ9aWNRu/EXl0Pe9/6K+hbGaKcSDBRE1NSCSA5XOoooJejRTqUYKIio6VWaHHZNVi/6OuSICdV2K7DRTrUiDRQicjuAzwDYraqnZtsmAVgNYCaAPgCLVfWdqPpIRFRIDx+CueZ7js/1nrMcm17Yh8GRPVUvXY1iKayTqEcU/wHgZgA/t7V9A8DjqnqjiHwj+3h5BH0jIirivuT1BmzbPRj4KqVaL4V1YpR+SXhU9SkA+wqaLwJwR/bvdwC4uKadIiJyoG/1OgYJOW4eEktWQkQ892TUs6hHFE6mqupbAKCqb4nI0U4vEpGlAJYCwIwZM2rYPSJqNn6T1XFZpRS0OAYKX1T1NgC3AUBXV5eWeDkRUdnM3/wndEd3Ubvxyasgx55Y1B6XVUpBi3TqyUW/iBwDANk/nWvqEhGFRFWRXrXCMUgklqx0DBJAZpVS0pC8tihWKQUtjmFuHYDPAbgx++cvo+0OETUT12T14m9DWlo93xuXVUpBi3p57F0AzgIwWUR2APgOMgFijYh8AcAbAC6NrodE1Czc6jMB5W2ci8MqpaBFGihU9QqXpxbUtCNE1NTC2lltqcV5F2GK49QTEVFNmM8/Cn35N7nH1rnV2ztOxNZZ52F+/0DVX+hhn3dRCwwURNSUCkcR1rnVuSNJA/pCD/O8i1phoCCiulbutI7bNNOj0y/Cn5NT8tqC+EJvhL0VDBREFEt+AkA50zpqpvHez1dgaCQNUxWGCNpbE2hrSSCxZCX+vLHXsR/VfqE3wt6K+ukpETUNvwHA77ROetWK3NSS9WpTFT+f9QWcPW8qOhHeF3pcKsBWI44b7oioyfmtmVRqWkf7X8tNNQ2NpGG/4tp5y3BYJXfNsDbLdU4djwWdU3IBp6M1iQWdU+omPwFwREFEMeR3Xt9rFFCYizA1EybWWsnqgmuGuVmu3vdWMFAQUez4nQZymtY5v28tpuMA0JLIe+2fpp2O5yZ+yPOa9f6FHhZOPRFR7PidBiqc1rms51ZMxwG0FQSJxJKVmHLGwoasw1QLHFEQUeyUMw3UOXU85mz4QeZB+5i854yLvgZpP9LXNet993SYGCiIKJb8TAPp8Lsw7/0Hx+ecym+4XbMRdk+HiYGCiOpSkPWZGmH3dJgYKIgoVEFP6Zjdv4M+95DDM4LEkhsqumYj7J4OEwMFEYUm6CmdsKq8NsLu6TDxUyCi0AQ1peN6mNBHL4Gc8OGq+gg0xu7pMDFQEDWxsFf6VDulo6ow7/yW43NBnRUBNO7JdEFhoCBqUrVY6VPNlI77kaTXQVraqu5bIW62c8dAQdSkqp0W8jMaqWRKR/fthPnwjxyfC3IUQf4xUBA1qWqmhfyORsqd0gn7SFKqDAMFUZOqZlqonNGInykdc+MvoG9uLWqXiVNhfPrLJftD4WKgIGpCPf0DOJw2i9r9rvQpZzRSaooq6FEES3EEj4GCqMkUThtZWhOCs+b6OyfB72jEa4pqzoYfYPhwuujEuSMWLoVMmx3IfxtLcQSD1WOJmozTtBEAjEkmfH+Z+q3u6vSzzNRhTHzg+tyJc9Y5EaYq7pj9RWyTyeX855T8eU4HHlF5OKIgajJBlKvwm6QuvOalW28BAJgABoZHn1vbuRQQA/Cx6spraomlOMLBQEHUZIIqV+EnSW39rDnvvIQP9j/t+Bq3E+eclJpaYimOcPDTI2oytSxXMX/2JEy4//sonugqDhAWry/1UqutWIojHAwURE2mVuUq0qtWYA6A4bZkLmENAG+Nm4Gnj7/A8T2lvtRLTS2xFEc4GCiImlDY5SrsS17bWhJoa0lg79AhrO68pui1AkABX1/qfqaWWIojeAwURA0oqr0ErvWZFv4t9qfGI+kwLbSg09+SXIBVXqPCQEHUYMrdS2APKuXc3dvp/n6YD97k+Jy1ca4z+7iaAMappWgwUBA1mHLKaxQGFetd5WxUK2dndRDTQpxaqj0GCqIGU85eArfNd0DpSrLp+34AHBxwfI5F/BoLAwVRgylnL0GpjWhuz7PKa3NhoCBqMDOPGosXdw7mtbklfN2Civ15O7cAIbM/DOOMSyroLdUDBgqiBtLTP4DuXUNF7SdNa3fMTzhVkLXYg4uqCfPO6xxfx1FE42OgIKpTTktg3XIOfW8fLHqvUwVZi301keuS17/6JqStvaJ+MhldX2IbKESkD8AggDSAlKp2RdsjovhwWwLr9sVfOL3kFlA6WpO4+oyZAADt3YL0hvsdr+d3FMGy340htoEi65OqujfqThDFjdsSWGsfRKHCXEOplVFBJaurPZeb4iHugYKIHLh90SsyuYXCL+fBkRRu2tgLAXDqsR2uSezLem5F+s0xjtfuPWc5Nm3uK2sKiWW/G0OcDy5SAI+KyLMisjTqzhDFiXi0L+ickhtBJAteqABe3DmICWMTRQcPLd56C9pbE0XXTCxZid5zluPxnj25L3hrCqmn33kfhcWtEizLfteXOP/fOlNVd4rI0QAeE5GtqvqU9WQ2eCwFgBkzZkTVR6JIOGciMu3WzuWe/gGs797t+Lo394/gvJOOxqbt+7DwhZszx5C2JdHWMhoojI9cCJn7UQCVTyGxNlNjiG2gUNWd2T93i8j9AE4H8JTt+dsA3AYAXV1dbr83RA3Ja/+DdZdvJY2dKIBne3Zg4Su3586qtgeJwlyE3ykkpxVOCzqncNVTnYtloBCRcQAMVR3M/v08AN+PuFtEsTF/9iTX0YJ1PrTbCihg9EhSIHNW9cBwCgPDKTxy2pcw/4TJuQJ+Fj+7vd1WOC3onJJbSUX1Ka45iqkAnhaRFwD8EcCDqvpIxH0iig2vO/LBkZTrCGBB3715QcJu7bxlGDxkOuYe5s+eVJTTKJxC8pqeovoWyxGFqm4H8IGo+0EUZ6XKbxTyChB2TrkHP+W9ucKpccUyUBA1uiB2Kzslip1YAUIAdLQlMTicyiXD3c6tdvpyL1Xeu5xihFRf+H+QqMaC2q1ceJfvxAoSVsIayCSy3QKEpaM1WXYw4wqnxsVAQVRjQe5Wtu7yf7a5Ly9Y2APE5PbRDXS/PvoCvN4y1fOaSUMw86ixZQcznj7XuBgoiGosyLl8+10/AEwbegMf2/EQgMxUk30DXWLJSry+sdfzel7FBf0EM54+15gYKIhqrC1pYDhVXN673Ln8wikse7LavjfCvifCK49gX8LqtvSWienmxEBBVEM9/QMYcQgShqDsuXzrrr9wNZN9uqlw45zfPAIT02TH/+tEAfNKAm/avs+x/EaLIWVP2QyOpByXvJqqrlVe/eYRmJgmOwYKogCVWtHkNnUzki6vCk161QpcNnQIhWOTd5Pt+O0pV+Fqj/f6ySMwMU12DBREVbKPIJzOg7Angd2mdATATRt7S34hazoF8+7vAMgkqgv3RCQNwYKA7vqZmCYLAwVRFQpHEG7jAus8iLak4RhM1PY6t2WohYcJWUX81s3+H9hjtvGun0LDQEFUBbcjRd0Mp0wYAowxBCNp5xPp7COQnv4BDD95J47e35s7g0Ixuqpp3FU34oqA/luI3DBQEFWhkuWipgKH0qVHIE9u68e8J/4Z1vjA/trVnddkppn6B3IBhfkECgsDBVEVvHIOXuOMUmOQxVtvgW51fs4qv2GvzBpESRAiNwwURFVwW0a6oHOKY2kNPxb7rPIKZIKC2y7qR7Ob5hgsqFoMFERVKLWM1G+FVyCzs9o6knRgeDS4/GniyXhu2scd3+NValwBjiwoEAwURFUqXEba0z+QG0l0tCZx0rR29L190HX5bPuhA/j0a3djckdrrs1a9upV5dXaAOdVPbbSYoNEdgwURAFy2nDXvWsoNxXlVJ9JALS3jf4qDh9O4555S2E6HECZECCtKBq5eI1aWJ+JqsVAQRSgUlVXrS/2iQ9cD1M1r3ifZU3nMpgOX+5tSQNL/2J2Ubt1zUe7dzsmyVmfiarFf0FEAfJTQnzOhh8AtjMiLFZ9pkGXUuBOFWctbiML1meiIDBQEAXIq+pq4c5qO7+lwL2wPhOFhYGCKEAzjxqLF3cOFrUveuX/Am0tRe1OVV6rqdzK+kwUBgYKogD1vX0w77FVBnxEBB22duOTV0GOPdHxGhwZUNwwUBAFyJoyOuGdl/Gh/t/m2k0dHR24nRVhx5EBxQkDBVGAOlqTWPjCzUXthoivAEEURyUDhYhcC2CVqr5Tg/4QRSKIonrpVSuw+HAag8jfVCcA9l98HaYG2WGiGvIzopgG4BkReQ7A7QDWq6r/uspEMVfqVDo/rBVN1n6IoZE0TFU88oFryw46rARLcVMyUKjqt0TkOgDnAbgawM0isgbAT1X1T2F3kChspTbJeXFa8trWkkDbxKOQuPjvPY8kdRJE0CIKmq8chaqqiOwCsAtACsCRAO4RkcdU9ethdpAobKU2yTnd4c89aizM1d91fF81uYhqghZRWPzkKL4M4HMA9gL4CYC/V9XDImIAeBUAAwXVNa8Nbk53+OPv/z76AbS1GBhJaa4Ux8AFX8WJx0+rqi9+dnYT1ZqfEcVkAJeo6uv2RlU1ReQz4XSLqHa8NrjZ7/DP3PEwjh3K/BoogIOHR0tqrO68BsnX3oU5ZqCqO/9Kd2UThclPjuLbHs91B9sdIv+CSvp6bXBbnz3851IfhwkFMUVUza5sorDwNoXqUtBJX7cNbpf13Jq3Wc7O7cS5anBXNsURAwXVpVokfdOrVqC9NZE7RMjidZiQfYqo0hEPd2VT3DBQUF0KM+lrX/Jq7YsYHE7hxcmno3vyh1zfZ58i4jJXaiQMFBRbXnfk1SZ9HZe8th6E+eBNRa9ta0lg3FU34kD/AHZk39OaEJgADqczgaA1IThr7pS8qSMuc6VGwUBBsVTqjryapK/TtSfc/32815bMO2kOAIwrb4CI5H6u03GmQOZ4Ujsuc6VGEttAISILAfw7gASAn6jqjRF3iWrI75GileQA7Ne2VjMpMmU37IEisWSl48jDz2iBy1ypkcTyX62IJAD8CMC5AHYgU2tqnaq+Em3PqFb83JE7JX39JJCtaxQueTVVsXfoEPZffJ3jyMEa1RQGCae+cZkrNZJYBgoApwPoVdXtACAidwO4CAADRZOo5I7cbwLZa8nr6s5rkMy+x23kIMivDuvWt4Rk6t0AQFvSwCdOnMz8BNUlI+oOuDgOwJ9tj3dk26hJzJ89CUlD8tpK3ZF7TQlZrCWvUvDetfOW5Za9Wu9xG9Voti9ufbMC1ogtceE2CiGqB3EdURT+HgMFN3EishTAUgCYMWNGLfpEISucNjppWjv63j7oOwfhNV2165av5f4BCUbrND06YxHeGXu043u8RjVWrsKpb1zxRI0mroFiB4DjbY+nA9hpf4Gq3gbgNgDo6uri7Vqdc5o26t41hAWdU/K+XN1yED39A45TQu870IPT33oyr92q0zSw6NtIbd8HeAQDtzyD16Y4rniiRhPXQPEMgBNFZBaANwFcDuDKaLtEYfJzF+6Wg9h54CC6dw0VBYlS9Zk6tu8rGQysvpWzsoornqjRxPJfrqqmskewrkdmeeztqvpyxN2iEPm5C3cLJi/tHMwLEm4BAsgvvzE4kioZDCopp8EVT9RoYhkoAEBVHwLwUNT9oNrwcxfulVy2+KnyWnjtoGsrsbAfNZrYBgpqLn7uwt2CiQD4rEuAeLN9FjZNP9/xPWHe4bOwHzUSBgqKBT934U7BpBUpfPbVn+KQwzWtUcT0ia3YO3QYw6nMQUOFdZmIyBsDBcWG1124tdrJvuHtsp5b0d6awJCZ/9r7534eKWMMBMB5Jx1dcUAI6mAkonrHQEGxV7ja6bT+36HznRfRni3iNzA8Oh1lz0UoKi/pzTLhRKMYKCgWvO7eSxXxM0SwuvOaomtWsxyVm+aIRjFQUORK3b0PjqQcVzNZ9Zr2X3wdkgEvR+WmOaJRDBQUuVJ374u33uJYhG/dKX+LMckEBrt3oy1pICHASFoDySdw0xzRKP6rp8i53aUvfOFmvNuTLAoSv51+AXa1z4BhKkay7x1OmUgagvOrSF7bcdMc0SgGCopc4d37EYcG8Ontd8IQwdBIOu+19mR1YUHWIHMI3DRHNIqBgiJnv3u3chECoL11dEWT085qJ0HmELhpjiiDgYIi1zl1PCY9/wDMvhdhYrTGvFeQ8Ht4EBFVj79VFLn0qhWYBADtYzB8OI3B4RQU7qOIpCE4aVo7uncNMYdAVAMMFBSZ9KoVRW1DI2m80XECfn/cuXnt1gjCnis4dsJY5hCIaoCBgmpOVWHe+S3H55w2zgGZIPHls+bktTGHQFQbDBRUU06jCAAwLlkOGTseHZv7uH+BKGaMqDtAzUF39+UFieHDaewdOoTdgyP4+ewvYttApn3+7ElIGvlHpjP3QBQt3qZR4ArrNi3uuQVtLYnc81bCeo2VrHYouMfcA1F8MFBQoOx1m87bvgYTDu3DYPY5K1gMjaRHg0SWfbMccw9E8cJAQYGy6jbZi/jZK70mlqzE6o29ju9lwT2ieGKgoEAtfOFmx/anjz0Xf3XheQBYcI+o3jCZTYHQdArpVStgiBQ9t3beMgxMnpt7zIQ1UX3hLRxVzb6aqb01kdtZfU/nUqgYRUGACWui+sJAQRXTN7fC3PiLvDYrYb2mcxnUIwgwYU1UPxgoqCJOG+eGD6expnNZZpQABHY2BBFFi4GCfLH2Rpz5yi8wITWA9tZE3t6IPdO7cJ+chlQ2SV14nCkR1S8GCirJ2hux6JUfAwBMAIPZEuDWkteHN/flgoQlZSrWd+/G+u7dzEMQ1TEGCipp4gPXY5Hmn/6gAH4160pc/onTAJTeA8ERBlH94vJYcqWHhpFetQKmFh8RtHbeMuzWI3KP/eyBsHZfE1F94YiCHNmT1YZILljYDxOyBwf7caZeuPuaqP4wUBCA0WR1+95efHzn+rxkdXtrAq+0zsLmYxbkvWfmUWNzf7emkx7t3u14RKmFu6+J6g9/a6lksnrcVTdiaFs/sHMw733du4Zw7ISxuSBh/ek2suDua6L6xEBBGPrNWiza153XpgA2HHchLjr3TABA39sHi95nr/hqKdx17XSEKRHVFwaKJmYdSXrs4EjRc2sLyoC75Rac2gt3XVvTWuu7d2PT9n0MGER1hoGiSaXv+Qdg5F0A+cnqe+f+T5hG5p+FPZ9QacVX+/kUAJfJEtUjBoomoweHYN73j3lt7a0JvKkdeGTWZbk2K59gP62ukJ+cg3U+hZ3TlBURxRcDRRNxqs8EAOOuuhHaP4COgmqugHti2m/OoZwpKyKKp9gFChH5LoAvAtiTbfo/qvpQdD2qf05VXgHA+NiVkBmnAHCu5vqzzX2uQeLqM2bmHheekW0PIDykiKj+xfW39d9U9V+i7kQjcBtFJJasLPleP6OBUjkIp414XCZLVF/iGiioSubv1kL7ni9qNy5ZDhnrLzfgZzRQKgfBQ4qI6l9cA8W1IvI3ALYA+KqqvhN1h+qFmibMu64rfsJIInHF98q6lp/RgJ9RBw8pIqpvkQQKEdkAYJrDUysA/BjA9cjs07oewL8C+LzDNZYCWAoAM2bMCK2v9cRtmsm48gaIw1nWpfgZDTAHQdT4RB0qg8aFiMwE8GtVPdXrdV1dXbply5aa9CmOdGAvzF/9W1G7nPxxGB88P9SfXZijADKjjgWdUziKIIo5EXlWVbtKvS52t30icoyqvpV9uAjAS1H2J+6qSVYHgTkIosYXu0AB4J9E5L8hM/XUB+CaaLsTT/rqMzD/+EBRu3HeNZAp/qbivJa1loM5CKLGFrtAoap/HXUf4i6IUQRLaxCRX7ELFOTO/OM66Kt/KGo3Fn8b0tJa1rVYWoOI/GKgqAOaOgxz9XeL2mXScTA+9aWKrsnSGkTkFwNFzKXv+g5gFn95V5us5rJWIvKL3woxpfv7YT54U1H7jvcvwhPvTsbgxt6qEtAsrUFEfjFQxJBbsrr3nOXZL/fMSKCaBDSXtRKRXwwUIStnCar2PgPzDw5LXi/+OmTcBGzaXFzNtZoENJe1EpEfDBQh8rsE1TqStMj4KUhc+JXcQyagiSgKDBQh8rME1XziDuhb24rea1xxPcQw8tqYgCaiKPAbJkReIwAdeQ/mPcUrl+QD58I49SzH9zEBTURRYKAIkdsI4LKeW2G+OaaovdSSVyagiSgKDBQhKhwBTH7vLZz9xi/R3pb/sZdTn4kJaCKqNQaKENlHAAtfuBmGCNrbkmhrSeReU6sqr0RElWKgCNncMe9izvb/B3Tk12IyPrsC0npERL0iIvKPgSIkqib08Z9B+7fntcuM98P42OUR9YqIqHwMFCHQndtgPnlHUXulR5ISEUWJgSJAmjoE894bgdRIXrtx0dcg7UdG1CsiouowUATE3LoJ+uyDeW1eeyKIiOoFA0WV9L0DMO//p6J2Y/F1kJa2CHpERBQsBooqmJvvg25/Nq/N+NiVkBmnFL02qPOpiYhqjYGiAvr2DpiP/DivTSYdB1m4DCJG0et5PjUR1TMGijKoaUIf/hF0/668duNTX4JMOs71fTyfmojqGQOFT/r6izCfvjuvTeZ8BMZ/v7jke1kenIjqGQNFCXpoGOba64vajUuWQ8b6Gw2wPDgR1TN+U3kw/+sJ6IuP57VJ14UwOj9a1nVYHpyI6hkDhQMdfBvmuh/mN445Asair0OSLWVfj+XBiaieMVDYqCr0qVXQHd157cbZV0OOmVPVtVkenIjqFQNFlva/BnPDT/LaZNocyNlXsT4TETW1pg8Umk5lppneO5DXbnzmK5AJUyLqFRFRfDR1oFBVmHd/J69NTv44jA+eH1GPiIjip6kDBQBg7Hjg4AAAHiZEROSkqQOFiMC44FpA0773RBARNZumDhQAIG3jou4CEVGsNXWgYEVXIqLSmjZQsKIrEZE/xTWxm4RXRVciIhrVtIGCFV2JiPxp2kDhVrmVFV2JiPJFEihE5FIReVlETBHpKnjumyLSKyI9IhLazrf5sychaeSX5mBFVyKiYlHdPr8E4BIAt9obReRkAJcDOAXAsQA2iMhcVU0H3QFWdCUi8ieSQKGq3QCciu1dBOBuVR0B8JqI9AI4HcDmMPrBiq5ERKXFLUdxHIA/2x7vyLYVEZGlIrJFRLbs2bOnJp0jImpGoY0oRGQDgGkOT61Q1V+6vc2hTR3aoKq3AbgNALq6uhxfQ0RE1QstUKjqORW8bQeA422PpwPYGUyPiIioEnGbeloH4HIRaRWRWQBOBPDHiPtERNTUoloeu0hEdgA4A8CDIrIeAFT1ZQBrALwC4BEAfxfGiiciIvJPVOt/el9E9gB4Pep+lGkygL1RdyJm+Jk44+fijJ9LsXI/k/epasmjPBsiUNQjEdmiql2lX9k8+Jk44+fijJ9LsbA+k7jlKIiIKGYYKIiIyBMDRXRui7oDMcTPxBk/F2f8XIqF8pkwR0FERJ44oiAiIk8MFBESkX8Wka0i8l8icr+ITIy6T1HzKkHfbERkYbbcfq+IfCPq/sSBiNwuIrtF5KWo+xInInK8iDwpIt3Z35//HeT1GSii9RiAU1X1NADbAHwz4v7EgVWC/qmoOxIlEUkA+BGATwE4GcAV2TL8ze4/ACyMuhMxlALwVVU9CcBHAfxdkP9eGCgipKqPqqp19urvkalt1dRUtVtVe6LuRwycDqBXVber6iEAdyNThr+pqepTAHiwfQFVfUtVn8v+fRBAN1wqb1eCgSI+Pg/g4ag7QbHhu+Q+kZ2IzATwQQB/COqaPCA6ZH7KrYvICmSGjqtq2beoVFiCvtnXmYKfAAABnklEQVT4LrlPZBGRdgD3AviKqg4EdV0GipCVKrcuIp8D8BkAC7RJ1ipXWIK+2bDkPpVFRFqQCRKrVPW+IK/NqacIichCAMsB/KWqvhd1fyhWngFwoojMEpExyJwlvy7iPlFMSeZc6Z8C6FbVHwZ9fQaKaN0MoAPAYyLyvIjcEnWHouZWgr7ZZBc5XAtgPTKJyTXZMvxNTUTuArAZQKeI7BCRL0Tdp5g4E8BfAzg7+13yvIhcENTFuTObiIg8cURBRESeGCiIiMgTAwUREXlioCAiIk8MFERE5ImBgoiIPDFQEBGRJwYKohCIyEey54y0ici47BkBp0bdL6JKcMMdUUhE5AYAbQDGAtihqv8YcZeIKsJAQRSSbI2mZwAMA5ivqumIu0RUEU49EYVnEoB2ZOp5tUXcF6KKcURBFBIRWYfMyXSzAByjqtdG3CWiivA8CqIQiMjfAEip6p3Z8683icjZqvpE1H0jKhdHFERE5Ik5CiIi8sRAQUREnhgoiIjIEwMFERF5YqAgIiJPDBREROSJgYKIiDwxUBARkaf/D68XlMU0zivmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x, y, color='#91bfdb')\n",
    "\n",
    "abline_values = [b * val + a for val in x]\n",
    "\n",
    "plt.plot(x, abline_values, lw=2.2, alpha=0.8, color='#fc8d59')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "\n",
    "&copy; Blitzstein, Joseph K.; Hwang, Jessica. Introduction to Probability (Chapman & Hall/CRC Texts in Statistical Science)."
   ]
  }
 ],
 "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": 2
}