{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparisons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "This notebook is part of [Bite Size Bayes](https://allendowney.github.io/BiteSizeBayes/), an introduction to probability and Bayesian statistics using Python.\n",
    "\n",
    "Copyright 2020 Allen B. Downey\n",
    "\n",
    "License: [Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The following cell downloads `utils.py`, which contains some utility function we'll need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "\n",
    "download('https://github.com/AllenDowney/BiteSizeBayes/raw/master/utils.py')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "If everything we need is installed, the following cell should run with no error messages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Review\n",
    "\n",
    "[The previous notebook](https://colab.research.google.com/github/AllenDowney/BiteSizeBayes/blob/master/10_joint.ipynb) I introduced cross tabulation, joint distributions, conditional disrtribution, and marginal distributions.\n",
    "\n",
    "In this notebook, we'll apply these ideas to Bayesian inference.\n",
    "\n",
    "But first I want to introduce a computational tool we will need, outer operations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outer operations\n",
    "\n",
    "Suppose you have two sequences, like `t1` and `t2`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = [1,3,5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "t2 = [2,4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many useful operations can be expressed in the form of an \"outer operation\" of these sequences.\n",
    "\n",
    "The most common outer operation is the outer product, which computes the product of every pair of values, one from each sequence.\n",
    "\n",
    "For example, here is the outer product of `t1` and `t2`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.multiply.outer(t1, t2)\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can understand this result more easily if we put it in a DataFrame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(a, index=t1, columns=t2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The values from `t1` appear along the rows of the result; the values from `t2` appear along the columns.  \n",
    "\n",
    "Each element of the result is the product of an element from `t1` and an element from `t2`.\n",
    "\n",
    "The outer sum is similar, except that each element is the *sum* of an element from `t1` and an element from `t2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "a = np.add.outer(t1, t2)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(a, index=t1, columns=t2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do the same thing with almost any operation.  For example, the \"outer greater than\" operation compares each element from `t1` to each element of `t2`; the result is an array of Boolean values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "a = np.greater.outer(t1, t2)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(a, index=t1, columns=t2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These outer operations work with Python lists and tuples, and NumPy arrays, but not Pandas Series.\n",
    "\n",
    "However, the following function works with Pandas Series, and puts the result into a DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def outer_product(s1, s2):\n",
    "    \"\"\"Compute the outer product of two Series.\n",
    "    \n",
    "    First Series goes down the rows;\n",
    "    second goes across the columns.\n",
    "    \n",
    "    s1: Series\n",
    "    s2: Series\n",
    "    \n",
    "    return: DataFrame\n",
    "    \"\"\"\n",
    "    a = np.multiply.outer(s1.to_numpy(), s2.to_numpy())\n",
    "    return pd.DataFrame(a, index=s1.index, columns=s2.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It might not be obvious yet why these operations are useful, but we'll see some examples soon.\n",
    "\n",
    "With that, we are ready to take on a new Bayesian problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How tall is A?\n",
    "\n",
    "Suppose I choose two people from the population of adult males in the U.S.; I'll call them A and B.  If we see that A taller than B, how tall is A?\n",
    "\n",
    "To answer this question:\n",
    "\n",
    "1. I'll use background information about the height of men in the U.S. to form a prior distribution of height,\n",
    "\n",
    "2. I'll construct a joint distribution of height for A and B, and update it with the information that A is taller, and \n",
    "\n",
    "3. I'll extract from the posterior joint distribution the posterior distribution of height for A."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the U.S. the average height of male adults in 178 cm and the standard deviation is 7.7 cm.  The distribution is not exactly normal, because nothing in the real world is, but the normal distribution is a pretty good model of the actual distribution, so we can use it as a prior distribution for A and B.\n",
    "\n",
    "Here's an array of equally-spaced values from roughly 3 standard deviations below the mean to 3 standard deviations above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = 178\n",
    "std = 7.7\n",
    "\n",
    "xs = np.arange(mean-24, mean+24, 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SciPy provides a function called `norm` that represents a normal distribution with a given mean and standard deviation, and provides `pdf`, which evaluates the probability distribution function (PDF), which we will use as the prior probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "ps = norm(mean, std).pdf(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll store the `xs` and `ps` in a Series that represents the prior PMF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior = pd.Series(ps, index=xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And normalize it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior /= prior.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior.plot()\n",
    "\n",
    "plt.xlabel('Height in cm')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Distribution of height for men in U.S.');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can think of this prior distribution as the marginal distribution for A and B, but what we want is the joint probability of their heights."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Joint distribution\n",
    "\n",
    "As we saw in the previous notebook, it is not *generally* possible to construct a joint distribution if we only have the marginals, because the marginals don't contain information about correlations between the variables.\n",
    "\n",
    "However, in the special case where there are no correlations, or they are small enough to ignore, it *is* possible to construct the joint distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see how, let's consider one element of the joint distribution,\n",
    "\n",
    "$P(A_y~\\mathrm{and}~B_x)$\n",
    "\n",
    "which is the probability that `A` is $y$ cm tall and `B` is $x$ cm tall.  We can rewrite this conjuction in terms of conditional probability:\n",
    "\n",
    "$P(A_y)~P(B_x~|~A_y)$\n",
    "\n",
    "We can compute $P(A_y)$ from the marginal distribution, but how should we compute the conditional probability, $P(B_x~|~A_y)$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the heights of `A` and `B` are \"independent\", which means that knowing the height of `A` provides no additional information about the height of `B`.\n",
    "\n",
    "And that means that the conditional probability, $P(B_x~|~A_y)$, is just the marginal probability $P(B_y)$.\n",
    "\n",
    "Which means that in this case, the joint probability is just the product of the marginal probabilities.\n",
    "\n",
    "$P(A_y~\\mathrm{and}~B_x) = P(A_y)~P(B_x)$\n",
    "\n",
    "Now, to compute the joint distribution, we have to compute this product for all values of $x$ and $y$.  And we can do that by computing the outer product of the marginal distributions, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "joint = outer_product(prior, prior)\n",
    "joint.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the prior is normalized, the joint prior should also be normalized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "joint.to_numpy().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function uses `pcolormesh` to plot the joint distribution.\n",
    "\n",
    "Recall that `outer_product` puts the values of `A` along the rows and the values of `B` across the columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_joint(joint):\n",
    "    \"\"\"Plot a joint distribution.\n",
    "    \n",
    "    joint: DataFrame representing a joint PMF\n",
    "    \"\"\"\n",
    "    plt.pcolormesh(joint.index, joint.index, joint)    \n",
    "    plt.ylabel('A height in cm')\n",
    "    plt.xlabel('B height in cm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "And here's what the result looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_joint(joint)\n",
    "plt.colorbar()\n",
    "plt.title('Joint prior distribution of height for A and B');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you might expect, the probability is highest near the mean and drops off away from the mean.\n",
    "\n",
    "Another way to visualize the joint distribution is a contour plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_contour(joint):\n",
    "    \"\"\"Plot a joint distribution.\n",
    "    \n",
    "    joint: DataFrame representing a joint PMF\n",
    "    \"\"\"\n",
    "    plt.contour(joint.index, joint.index, joint)\n",
    "    plt.ylabel('A height in cm')\n",
    "    plt.xlabel('B height in cm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_contour(joint)\n",
    "plt.title('Joint prior distribution of height for A and B');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each circle represents a level of equal probability.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Likelihood\n",
    "\n",
    "Now that we have a joint PMF that represents the prior distribution, we can update it with the data, which is that `A` is taller than `B`.\n",
    "\n",
    "Each element in the joint distribution represents a hypothesis about the heights of `A` and `B`; for example:\n",
    "\n",
    "* The element `(180, 170)` represents the hypothesis that `A` is 180 cm tall and `B` is 170 cm tall.  Under this hypothesis, the probability that `A` is taller than `B` is 1.\n",
    "\n",
    "* The element `(170, 180)` represents the hypothesis that `A` is 170 cm tall and `B` is 180 cm tall.  Under this hypothesis, the probability that `A` is taller than `B` is 0.\n",
    "\n",
    "To compute the likelihood of every pair of values, we can extract the values from the prior, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = prior.index.to_numpy()\n",
    "X = prior.index.to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then apply the `outer` version of `np.greater`, which compares every element of `Y` (height of `A`) to every element of `X` (height of `B`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.greater.outer(Y, X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is an array, which we can put in a DataFrame with the corresponding `index` and `columns`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "likelihood = pd.DataFrame(a, index=Y, columns=X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what it looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_joint(likelihood)\n",
    "plt.title('Likelihood of A>B');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The likelihood of the data is 1 where `Y>X` and 0 otherwise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The update\n",
    "\n",
    "We have a prior, we have a likelihood, and we are ready for the update.  As usual, the unnormalized posterior is the product of the prior and the likelihood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "unnorm_posterior = joint * likelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can get the normalized posterior by dividing through by the total."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "total = unnorm_posterior.to_numpy().sum()\n",
    "joint_posterior = unnorm_posterior / total\n",
    "total"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The total probability of the data is a little less than $1/2$.\n",
    "\n",
    "Here's what the normalized posterior looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_joint(joint_posterior)\n",
    "plt.colorbar()\n",
    "\n",
    "plt.title('Joint posterior distribution of height for A and B');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like a sunrise as seen from the deck of a [heeling sailboat](https://en.wikipedia.org/wiki/Sailing#Heeling)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The marginals\n",
    "\n",
    "From the posterior joint distribution we can extract the posterior marginal distribution of `A` and `B`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def marginal(joint, axis):\n",
    "    \"\"\"Compute a marginal distribution.\n",
    "    \n",
    "    axis=0 returns the marginal distribution of the second variable\n",
    "    axis=1 returns the marginal distribution of the first variable\n",
    "    \n",
    "    joint: DataFrame representing a joint PMF\n",
    "    axis: int axis to sum along\n",
    "    \n",
    "    returns: Series representing a marginal PMF\n",
    "    \"\"\"\n",
    "    return joint.sum(axis=axis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "marginal_A = marginal(joint_posterior, axis=1)\n",
    "marginal_B = marginal(joint_posterior, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what they look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior.plot(label='Prior')\n",
    "marginal_A.plot(label='Posterior for A')\n",
    "marginal_B.plot(label='Posterior for B')\n",
    "\n",
    "plt.xlabel('Height in cm')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Prior and posterior distributions for A and B')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you might expect, the posterior distribution for `A` is shifted to the right and the posterior distribution for `B` is shifted to the left.\n",
    "\n",
    "We can summarize the results by computing the posterior means:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pmf_mean(pmf):\n",
    "    \"\"\"Compute the mean of a PMF.\n",
    "    \n",
    "    pmf: Series representing a PMF\n",
    "    \n",
    "    return: float\n",
    "    \"\"\"\n",
    "    return np.sum(pmf.index * pmf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_mean(prior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_mean(marginal_A), pmf_mean(marginal_B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the observation that `A` is taller than `B`, we are inclined to believe that `A` is a little taller than average, and `B` is a little shorter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the posterior distribution are a little taller and narrower than the prior.  We can quantify that my computing their standard deviations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pmf_std(pmf):\n",
    "    \"\"\"Compute the standard deviation of a PMF.\n",
    "    \n",
    "    pmf: Series representing a PMF\n",
    "    \n",
    "    return: float\n",
    "    \"\"\"\n",
    "    deviation = pmf.index - pmf_mean(pmf)\n",
    "    var = np.sum(deviation**2 * pmf)\n",
    "    return np.sqrt(var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_std(prior), pmf_std(marginal_A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The standard deviation of the posterior distributions are a little smaller, which means we are a little more certain about the heights of `A` and `B` after we compare them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conditional posteriors\n",
    "\n",
    "Now suppose we measure `B` and find that he is 185 cm tall.  What does that tell us about `A`?\n",
    "\n",
    "We can answer that question by extracting the conditional posterior distribution for `A`, conditioned on `B=185`.\n",
    "\n",
    "Possible heights for `A` run down the rows of the joint PMF, so each row is an unnormalized posterior distribution conditioned on `A`.\n",
    "\n",
    "And possible heights for `B` run across the columns, so each column is an unnormalized posterior distribution conditioned on `B`.\n",
    "\n",
    "So we can condition on `B` by selecting a column and normalizing it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "cond_A = joint_posterior[185].copy()\n",
    "cond_A /= cond_A.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior.plot(label='Prior')\n",
    "marginal_A.plot(label='Posterior for A')\n",
    "cond_A.plot(label='Posterior for A given B=185', color='C4')\n",
    "\n",
    "plt.xlabel('Height in cm')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Prior, posterior and conditional distribution for A')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The posterior conditional distribution is cut off at 185 cm, because we have established that `A` is taller than `B` and `B` is 185 cm.\n",
    "\n",
    "And the posterior conditional is substantially different from the unconditional posterior; that is, for each value of $y$\n",
    "\n",
    "$P(A_y | B_x) \\ne P(A_y)$\n",
    "\n",
    "which means that in the posterior distribution, `A` and `B` are not independent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Elo rating\n",
    "\n",
    "[The Elo rating system](https://en.wikipedia.org/wiki/Elo_rating_system) is a way to quantify the skill level of players for games like chess.\n",
    "\n",
    "It is based on a model of the relationship between the ratings of players and the outcome of a game.  Specifically, if $R_A$ is the rating of player $A$ and $R_B$ is the rating of player $B$, the probability that $A$ beats $B$ is given by the [logistic function](https://en.wikipedia.org/wiki/Logistic_function):\n",
    "\n",
    "$P(A~\\mathrm{wins}) = 1 / (1 + 10^{(R_B-R_A)/400})$\n",
    "\n",
    "The parameters $10$ and $400$ are arbitrary choices that determine the range of the ratings.  In chess, values range from 100 to 2800.\n",
    "\n",
    "Notice that the probability of winning depends only on the difference in rankings.  As an example, if $R_A$ exceeds $R_B$ by 100 points, the probability that $A$ wins is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "1 / (1 + 10**(-100/400))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Suppose `A` has a current rating of 1600, but we are not sure it is accurate.  We could describe their true rating with a normal distribution with mean 1600 and standard deviation 100, to indicate our uncertainty.\n",
    "\n",
    "And suppose `B` has a current rating of 1800, with the same level of uncertaintly.\n",
    "\n",
    "Finally, `A` and `B` play and `A` wins.  How should we update their ratings?\n",
    "\n",
    "To answer this question:\n",
    "\n",
    "1. Construct prior distributions for `A` and `B`.\n",
    "\n",
    "2. Use them to construct a joint distribution, assuming that the prior distributions are independent.\n",
    "\n",
    "3. Use the logistic function above to compute the likelihood of the outcome under each joint hypothesis.  Hint: use `np.subtract.outer`.\n",
    "\n",
    "4. Use the joint prior and likelihood to compute the joint posterior. \n",
    "\n",
    "5. Extract and plot the marginal posteriors for `A` and `B`.\n",
    "\n",
    "6. Compute the posterior means for `A` and `B`.  How much should their ratings change based on this outcome?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this notebook I started with the \"outer\" operations, like outer product and outer sum; then we used them to construct a joint distribution.\n",
    "\n",
    "In general, you cannot construct a joint distrubution from two marginal distributions, but in the special case where the distributions are independent, you can.\n",
    "\n",
    "We extended the Bayesian update process we've seen in previous notebook and applied it to a joint distribution.  Then from the posterior joint distribution we extracted posterior marginal distributions and posterior conditional distributions.\n",
    "\n",
    "As an exercise, you had a chance to apply the same process to a slightly more difficult problem, updating Elo ratings based on the outcome of a chess game.\n",
    "\n",
    "[In the next notebook](https://colab.research.google.com/github/AllenDowney/BiteSizeBayes/blob/master/12_binomial.ipynb) we'll get back to a problem we left half-finished: the Euro problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}