csc421_tzanetakis_learning_probabilistic_models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CHAPTER 20 - Learning Probabilistic Models\n",
    "\n",
    "### George Tzanetakis, University of Victoria \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WORKPLAN \n",
    "\n",
    "The section number is based on the 4th edition of the AIMA textbook and is the suggested\n",
    "reading for this week. Each list entry provides just the additional sections. For example the Expected reading include the sections listed under Basic as well as the sections listed under Expected. Some additional readings are suggested for Advanced. \n",
    "\n",
    "1. Basic: Sections **20.1**, **20.2.1**, **20.2.2**, and **Summary**\n",
    "2. Expected: Same as Basic plus **20.3**, **20.3.1**, **20.3.3**\n",
    "3. Advanced: All the chapter including bibligraphical and historical notes \n",
    "\n",
    "\n",
    "We have covered a variety of probabilistic models that model uncertainty and allow us to do inference in different ways. In this notebook we describe some of the ways we can estimate probabistic modesl from data. \n",
    "These techniques provide the connection between statistics, probability, and machine learning. \n",
    "\n",
    "The ideas are based on Chapter 20 of the Artificial Intelligence: a Modern Approach textbook and specifically section 20.2 Learning from Complete Data. \n",
    "\n",
    "**Density estimation** refers to the task of learning the probabiity density function (for continuous models) or the probability distribution function (for discrete models) given some data that we assume was generated from that modal. **Complete data** means that we have data for all the **variables** in our model. \n",
    "\n",
    "The most common type of learning is **parameter learning** where we assume a particular structure for our model and characterize it by estimating a set of parameters. For example we might assume a normal or Gaussian multi-variate distribution and estimate the mean vector and the covariance matrix that characterize it. As another example, we might be given the structure of a Bayesian network (in terms of parent/child coonditional relationships) and learn the conditional probabilty tables. We will also briefly discuss the problem of learning structure as well as non-parametric density estimation in which we don't need to make any assumptions about the model. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning with Complete Data \n",
    "\n",
    "## A random variable class \n",
    "\n",
    "Define a helper random variable class based on the scipy discrete random variable functionality providing both numeric and symbolic RVs. You don't need to look at the implementation - the usage will be obvious through the examples below. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline \n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import numpy as np \n",
    "\n",
    "\n",
    "class Random_Variable: \n",
    "    \n",
    "    def __init__(self, name, values, probability_distribution): \n",
    "        self.name = name \n",
    "        self.values = values \n",
    "        self.probability_distribution = probability_distribution  \n",
    "        if all(issubclass(type(item), np.integer) for item in values): \n",
    "            self.type = 'numeric'\n",
    "            self.rv = stats.rv_discrete(name = name, values = (values, probability_distribution))\n",
    "        elif all(type(item) is str for item in values): \n",
    "            self.type = 'symbolic'\n",
    "            self.rv = stats.rv_discrete(name = name, values = (np.arange(len(values)), probability_distribution))\n",
    "            self.symbolic_values = values \n",
    "        else: \n",
    "            self.type = 'undefined'\n",
    "            \n",
    "    def sample(self,size): \n",
    "        if (self.type =='numeric'): \n",
    "            return self.rv.rvs(size=size)\n",
    "        elif (self.type == 'symbolic'): \n",
    "            numeric_samples = self.rv.rvs(size=size)\n",
    "            mapped_samples = [self.values[x] for x in numeric_samples]\n",
    "            return mapped_samples \n",
    "        \n",
    "    def probs(self): \n",
    "        return self.probability_distribution\n",
    "    \n",
    "    def vals(self): \n",
    "        print(self.type)\n",
    "        return self.values \n",
    "    \n",
    "    def prob_of_value(self, value): \n",
    "        indices = np.where(self.values == value)\n",
    "        return self.probability_distribution[indices[0][0]]\n",
    "\n",
    "\n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Likelihood of model given some data \n",
    "\n",
    "First let's review the concept of likelihood of a model given some data \n",
    "\n",
    "\n",
    "Lets start by creating a random variable corresponding to a 6-faced dice where there are two faces with the numbers 1,2 and 3 therefore each number appears with equal probability. We can generate random samples from this model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 1 3 3 2 2 3 2 2 2 3 1 2 3 3 2 1 3 3 2 2 1 1 1 2 1 1 2 2 1]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "values = np.array([1, 2, 3])\n",
    "probabilities = [2/6., 2/6., 2/6.]\n",
    "dice1 = Random_Variable('dice1', values, probabilities)\n",
    "samples = dice1.sample(30)\n",
    "print(samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets also create a random variable where three of the faces have the number 2, 2 faces have the number 1 and 1 face has the number 3. We can also generate random samples from this model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2 1 2 3 2 1 3 2 2 1 2 1 3 3 2 2 2 2 1 3 3 1 1 1 2 2 2 1 2 2]\n"
     ]
    }
   ],
   "source": [
    "values = np.int_([1, 2, 3])\n",
    "probabilities = [2./6, 3./6, 1./6]\n",
    "dice2 = Random_Variable('dice2', values, probabilities)\n",
    "samples = dice2.sample(30)\n",
    "print(samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The likelihood of a sequence of samples given a model can be obtained by taking the product of the corresponding \n",
    "probabilities. We can see that for this particular sequence of data the likelihood of the model for dice2 is higher. So if we have some data and some specific models we can select the model with the highest likelihood. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.3333333333333333\n",
      "0.3333333333333333\n",
      "0.16666666666666666\n",
      "Likelihood for dice1: 0.000017\n",
      "Likelihood for dice2: 0.000021\n"
     ]
    }
   ],
   "source": [
    "data = [1,2,2,1,1,3,1,2,3,2]\n",
    "print(dice1.prob_of_value(1))\n",
    "print(dice1.prob_of_value(3))\n",
    "print(dice2.prob_of_value(3))\n",
    "\n",
    "def likelihood(data, model):\n",
    "    likelihood = 1.0 \n",
    "    for d in data: \n",
    "        likelihood *= model.prob_of_value(d)\n",
    "    return likelihood \n",
    "    \n",
    "print(\"Likelihood for dice1: %f\" % likelihood(data,dice1))\n",
    "print(\"Likelihood for dice2: %f\" % likelihood(data,dice2))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that even with only 10 values the likelihood gets relatively small and we can expect it will get smaller as the sequences of data get smaller. We can also use log-likelihood to avoid this problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.3333333333333333\n",
      "0.3333333333333333\n",
      "0.16666666666666666\n",
      "Likelihood for dice1: -10.986123\n",
      "Likelihood for dice2: -10.750557\n"
     ]
    }
   ],
   "source": [
    "data = [1,2,2,1,1,3,1,2,3,2]\n",
    "print(dice1.prob_of_value(1))\n",
    "print(dice1.prob_of_value(3))\n",
    "print(dice2.prob_of_value(3))\n",
    "\n",
    "def log_likelihood(data, model):\n",
    "    likelihood = 0.0 \n",
    "    for d in data: \n",
    "        likelihood += np.log(model.prob_of_value(d))\n",
    "    return likelihood \n",
    "    \n",
    "print(\"Likelihood for dice1: %f\" % log_likelihood(data,dice1))\n",
    "print(\"Likelihood for dice2: %f\" % log_likelihood(data,dice2))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case above we examined two possible models. One could ask the question of all possible models for a particular problem can we find the one with the highest likelihood ? If we have a dice with six faces that can only have the number 1, 2, and 3 then there is a finite amount of models and we can calculate their likelihoods as we did above. However if we relax the requirement to have a dice and simply have the values 1,2 and 3 but with arbitrary associated probabilities then we have an infinte number of possible models. Without going into the math it turns out that at least for this particular case the model that will have the maximum likelihood can be simply obtained by counting the relative frequencies of the values in the data. This is called maximum likelihood estimation of model parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Counter({2: 7, 1: 4, 3: 2})\n",
      "[0.3076923076923077, 0.5384615384615384, 0.15384615384615385]\n"
     ]
    }
   ],
   "source": [
    "import collections \n",
    "\n",
    "data = [1,2,2,1,1,3,1,2,3,2,2,2,2]\n",
    "counts = collections.Counter(data)\n",
    "print(counts)\n",
    "est_probability_distribution = [counts[1]/float(len(data)), counts[2]/float(len(data)), counts[3]/float(len(data))]\n",
    "print(est_probability_distribution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2 1 3 2 1 1 2 1 2 2 2 2 2 2 2 1 3 2 2 2 2 1 3 2 2 1 2 2 3 2]\n"
     ]
    }
   ],
   "source": [
    "values = np.int_([1, 2, 3])\n",
    "probabilities = est_probability_distribution \n",
    "model = Random_Variable('model', values, probabilities)\n",
    "samples = model.sample(30)\n",
    "print(samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximum-likelihood parameter learning for Discrete Models  \n",
    "\n",
    "\n",
    "Lets start by creating a random variable corresponding to a bag of candy with two types lime and cherry similar to what we did in the previous notebook.  We can easily generate random samples from this model. For example in the code below we geenerate 100 samples. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['l', 'l', 'l', 'l', 'l', 'l', 'c', 'c', 'l', 'l', 'l', 'c', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'c', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'c', 'l', 'l', 'c', 'l', 'l', 'c', 'l', 'l', 'l', 'l', 'l', 'c', 'l', 'l', 'l', 'l', 'l', 'c', 'l', 'c', 'l', 'l', 'l', 'l', 'c', 'c', 'l', 'l', 'l', 'l', 'c', 'l', 'l', 'c', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'c', 'c', 'l', 'c', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'c', 'l', 'l', 'l', 'l', 'l', 'l', 'l', 'c', 'l']\n"
     ]
    }
   ],
   "source": [
    "values = ['c', 'l']\n",
    "probabilities = [0.2, 0.8]\n",
    "dice1 = Random_Variable('bag1', values, probabilities)\n",
    "samples = dice1.sample(100)\n",
    "print(samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now imagine that you are just given these samples and you are told that they were from a bag of candy but you don't know the percentage of each candy type in the bag and you need to estimate it. Let's call the probability a candy from the bag is cherry $\\theta$. Then our task of parameter learning is to estimate $\\theta$ from the provided samples. In the previous notebook without much explanation I stated that the \"best\" possible model in a maximum likelihood sense can be easily obtained by simply counting the percentage of each candy type in our bag. \n",
    "As you can see the estimated paremeter $\\theta$ is close but not the same as the original value which was $0.2$. If we had more samples this estimate becomes more accurate.  \n",
    "\n",
    "We can see that with this simple example we have the ability to \"learn\" a model. Once we have a \"learned\" model from the data we can use it to make predictions or inference in general as well as generate samples if needed. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Counter({'l': 81, 'c': 19})\n",
      "[0.19, 0.81]\n"
     ]
    }
   ],
   "source": [
    "import collections \n",
    "\n",
    "counts = collections.Counter(samples)\n",
    "print(counts)\n",
    "est_probability_distribution = [counts['c']/float(len(samples)), counts['l']/float(len(samples))]\n",
    "print(est_probability_distribution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the counts seems intuitive and I told you that for the case of discrete random variables this provides the maximum likelihood estimate but can we prove this assertion? \n",
    "\n",
    "Here is how we can do it. Each time we have a candy of a particular type we multiply the associated probability top get the likelihood of the sequence (assuming i.i.d. samples). If there are $c$ cherry candies and $l=N-c$ limes then we can write the likelihood as follows: \n",
    "\n",
    "$$ P({\\bf d} | h_{\\theta}) = \\prod_{j=1}^{N} P(d_j | h_{\\theta}) = \\theta^{c} * (1-\\theta)^{l}$$\n",
    "\n",
    "Note: check how the mathematical expression above is notated. It is using LaTeX notation which can be embedded in markdown cells. It is a useful thing to learn to produce nice looking equations in both notebook and papers. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The maximum-likelihood hypothesis is given by the value of $\\theta$ that maximizes the exression above. The same value can be obtained by maximizing the **log likelihood**. Note that we have used log-likelihood before to avoid small numerical likelihood values when computing over long sequences. Here we use it because it allows us to simplify our expression to prove our approach to maximum likelihood parameter estimation. By taking the log we convert the product to a sum which is easer to maximize. \n",
    "\n",
    "$$ L({\\bf d}| h_{\\theta}) = \\log{P({\\bf d} | h_{\\theta})} = \\sum_{j=1}^{N}\\log{P(d_j| h_{\\theta})} = c \\log{ \\theta} + l \\log{(1-\\theta)} $$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find the maximum-likelihood value of $\\theta$, we differentiate the $L$ with respect to $\\theta$ and se the resulting expression to zero: \n",
    "\n",
    "$$ \n",
    "\\frac{L({\\bf d} / h_{\\theta})}{d \\theta} = \\frac{c}{\\theta} - \\frac{l}{1-\\theta} = 0 \n",
    "$$\n",
    "\n",
    "Solving for $\\theta$ we get: \n",
    "$$ \n",
    "\\theta = \\frac{c}{c+l} = \\frac{c}{N}\n",
    "$$ \n",
    "\n",
    "This might seem like a lot of work to prove something obvious but now we actually know that of all the infinite possible models of bags we could have - the one we estimate by counting the proportion of candy is the \"best\" in a maximum likelihood sense. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The approach we followed can be used for a variety of probabilistic models. The steps are as follows: \n",
    "\n",
    "1. Write down an expression for the likelihood of the data as a function of the parameters and use log to simplify it for step 2 \n",
    "2. Write down the derivative of the log likelihood with respect to each parameter \n",
    "3. Find the parameter values such that the derivatives are zero. \n",
    "\n",
    "Note: If we are lucky we are able to perform steps 2 and 3 analytically and derive an exact ML parameter estimate. There are many cases especially when dealing with continuous models (which we cover below) in which \n",
    "maximimizing the likelihood function analytically is not possible and one needs to resort to numerical methods \n",
    "which do not provide an exact solution. \n",
    "\n",
    "As another example of analytical ML parameter estimation, the book has one more example in which there is an extra random variable wrapper and the model has three parameters $\\theta_1, \\theta_2, \\theta_3$. \n",
    "\n",
    "For example by basicallly filtering the data and counting we can do ML parameter estimation for Naive Bayes models as well as Baysian Networks with discrete random variables. You have already seen to some extent how this can be done during lectures as well as in the assignments. Notice that the structure of the Bayesian network allows us to simplify the problem of ML parameter estimation by factoring different groups of variables based on their conditional structural relationships. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximum-Likelihood parameter learning for continuous models \n",
    "\n",
    "\n",
    "Continuous probability models are heavily used in real-world applications. Like I mentioned in the previous video - in many cases we need to resort to numerical optimization methods to perform parameter estimatino. However in some cases we can get the exact answer analytically. Let's consider the simple example of learning the parameters of a Gaussian density function on a single variable. Similarly to what we did in the previous video for learning the parameters of a disrete random variable we will first generate some data and then estimate the parameters from the data. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data will be generated using a Gaussian density function on a single variable. The corresponding equation is: \n",
    "\n",
    "$$ P(x) = \\frac{1}{\\sqrt{2 \\pi \\sigma}} e^{-\\frac{x-\\mu}{2\\sigma^2}}$$\n",
    "\n",
    "The parameters of this model are the mean $\\mu$ and the standard deviation $\\sigma$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.90383961 3.11383124 2.74047631 2.81324639 2.85427437 2.83708547\n",
      " 3.27187123 2.75284102 3.29324794 3.17367798]\n"
     ]
    }
   ],
   "source": [
    "mu = 3.0 \n",
    "sigma = 0.2 \n",
    "s = np.random.normal(mu, sigma, 10)\n",
    "print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = np.random.normal(mu, sigma, 10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let the observed values by $x_1, \\dots, x_N$. Then the log-likelihood is: \n",
    "\n",
    "$$ \n",
    "L = \\sum_{j=1}^{N} \\log{\\frac{1}{\\sqrt{2 \\pi \\sigma}}} e^{-\\frac{x-\\mu}{2\\sigma^2}} = N(-\\log \\sqrt{2\\pi} - \\log \\sigma) - \\sum_{j=1}^{N} \\frac{x-\\mu}{2\\sigma^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setting the derivatives to zero we obtain: \n",
    "\n",
    "$$ \n",
    "\\frac{\\partial L}{\\partial \\mu} = - \\frac{1}{\\sigma^2}\\sum_{j=1}^{N} (x_j-\\mu) = 0 \n",
    "$$\n",
    "which implies: \n",
    "$$ \n",
    "\\mu = \\frac{\\sum_j x_k}{N}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the maximum likelihood value of the mean is the sample average. Similarly you can find that the maximum likelihood value of the standard deviation is the square root of the sample variance. You can check the textbook for the details of the standard deviation $\\sigma$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check how we can calculate these ML parameter estimates for the data that we have. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.9967265838293855\n",
      "2.9967265838293855\n",
      "0.19763576241858616\n"
     ]
    }
   ],
   "source": [
    "estimated_mean1 = np.sum(samples) / len(samples)\n",
    "print(estimated_mean1)\n",
    "estimated_mean2 = np.mean(samples)\n",
    "print(estimated_mean2)\n",
    "estimated_std = np.std(samples)\n",
    "print(estimated_std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So armed with simple filtering, counting and calculating sample mean and sample standard deviation we have everything we need to estimate the probabilities of a Naive Bayes model that contains a mixture of continuous and discrete variables. For the discrete variables we count and estimate directly the probabilities. For the continuous variables, we first estimate the ML parameters (sample mean and standard deviation) and then for a particular value of the feature we use the single varaible Gaussian density equation to derive a probability value for that value. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some notes for further reading - not needed for the \"final\" assignment for those interested in digging deeper. \n",
    "\n",
    "**Note1**: In a Baysian network with continuous variables you have the problem of having continuous parent and a continuous child variable. These can be addressed with linear Gaussian models. More details in the textbook \n",
    "\n",
    "**Note2**: Similarly to Baysian learning in discrete models one can follow a similar approach and use a hypothesis prior to guide the learning. The textbook shows an example that use **beta distributions** you can check out. \n",
    "\n",
    "**Note3**: If you remember when we covered Baysian network we looked at approximate inference using direct sampling and rejection sampling. You will notice that the approach we followed was similar to statistical learning in the sense that we generated samples and then used counting to estimate probabilities. So at a basic level inference and learning can be considered the same process. We start with a few things that we know and then using data we update what we know. \n",
    "\n",
    "**Note4**: It is also possible to learn the structure of a Bayesian network from data. The basic idea is to search over the space of possible model. To do so we will need some method to determine when a good structure has been found. More details can be found in the book. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Learning\n",
    "\n",
    "The two important terms we will cover are data and hypotheses or models. The hypotheses are different probabilistic theories. Let's consider the example described by the book. \n",
    "\n",
    "We have a candy manufacturer that produces bags of candy wrapped in the same opaque wrapper. The flavors are cherry and lime. There are 5 kinds of bags: \n",
    "\n",
    "h1: 100% cherry\n",
    "h2: 75% cherry and 25% lime \n",
    "h3: 50% cherry and 50% lime \n",
    "h4: 25% cherry and 75% lime \n",
    "h5: 100% lime \n",
    "\n",
    "Given a new bag of candy the random variable *H* takes one of these 5 values: h1, h2, h3, h4, h5. We don't know which type it is and we gradually unwrap candy D1, D2, D3,..., DN where each of those is a random varilable with value $cherry$ and $lime$ (the bags are really big so replacement does not make a difference). \n",
    "\n",
    "The task is given a sequence of observations D1, ... DN to predict the flavor of the next piece of candy. \n",
    "As an extreme example consider we observe a sequence of 100 lime candies then we have a high confidence the bag is of type h1 and therefore the next candy will be lime. \n",
    "\n",
    "One approach to solving this type of problem is the maximum likelihood described above. This consists of selecting the \"best\" hypothesis using the maximum likelihood and then doing the prediction using that hypothesis. If each candy bag is equally likely then this works. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 10\n",
      "3 7\n"
     ]
    }
   ],
   "source": [
    "data_book = ['l'] * 10\n",
    "data_other = ['l','c','l','l','l','l','l','c','l','c']\n",
    "\n",
    "counts_book = collections.Counter(data_book)\n",
    "print(counts_book['c'], counts_book['l'])\n",
    "counts_other = collections.Counter(data_other)\n",
    "print(counts_other['c'], counts_other['l'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Likelihood of h1 for data_book: 0.0\n",
      "Likelihood of h5 for data_book: 1.0\n",
      "Likelihood of h3 for data_book: 0.0009765625\n",
      "Likelihood of h3 for data_other: 0.0020856857299804688\n"
     ]
    }
   ],
   "source": [
    "h1 = [1.0, 0.0] \n",
    "h2 = [0.75, 0.25]\n",
    "h3 = [0.5, 0.5]\n",
    "h4 = [0.25, 0.75] \n",
    "h5 = [0.0, 1.0]\n",
    "\n",
    "prior = [0.1, 0.2, 0.4, 0.2, 0.1]\n",
    "#prior = [0.2, 0.2, 0.2, 0.2, 0.2]\n",
    "\n",
    "\n",
    "def likelihood(d, h): \n",
    "    counts = collections.Counter(d)\n",
    "    return np.power(h[0], counts['c']) * np.power(h[1], counts['l'])\n",
    "\n",
    "print('Likelihood of h1 for data_book:', likelihood(data_book, h1))\n",
    "print('Likelihood of h5 for data_book:', likelihood(data_book, h5))\n",
    "print('Likelihood of h3 for data_book:', likelihood(data_book, h3))\n",
    "print('Likelihood of h3 for data_other:', likelihood(data_other, h4))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now consider that we have a prior probability distribution for the hypotheses. For the candy bag scenario \n",
    "let's say we know from the manufacturer that 10% of candy bags are h1, 20% are h2, 40% h3, 20% h4 and 10% h5. \n",
    "\n",
    "An alternative more general approach, called $Baysian Learning$ is to calculate the probability of each hypothesis given the data and then use all the hypotheses, weighted by their probabilitities rather than \n",
    "just selecting the \"best\" by maximum likelihood to perform the prediction. Notice that with this approach we can take into account the prior probability over the hypotheses. \n",
    "\n",
    "Mathematically we can calculate the probability of each hypothesis by weighting the likilehood by the prior: \n",
    "\n",
    "$${\\bf P}(h_{i} | {\\bf d} ) = \\alpha {\\bf P}({\\bf d} | h_{i}){\\bf P}(h_i) $$\n",
    "\n",
    "Suppose we want to make a prediction about an unknown quantity X such as predicting what the next candy will be. Then we have: \n",
    "\n",
    "$$ {\\bf P}(X | {\\bf d}) = \\sum_{i} {\\bf P}(X| {\\bf d}, h_i){\\bf P}(h_i|{\\bf d}) = \\sum_{i} {\\bf P}(X|h_i) {\\bf P}(h_i | {\\bf d})$$\n",
    "\n",
    "where we assume that predictions are weighted averages over the predictions of the individual hypotheses. The key terms are the **hypothesis** prior $P(h_i)$ and the **likelihood** of the data under each hypothesis $P({\\bf d} | h_i)$. \n",
    "\n",
    "Assuming that the observatois are i.i.d we can calculate easily the likelihood like we did before and the priors for the hypotheses are given: \n",
    "\n",
    "$$ P({\\bf d} | h_i) = \\prod_{j} P(d_k, h_i) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0  1  2  3  4  5  6  7  8  9 10]\n",
      "[0.4        0.4        0.30769231 0.21052632 0.1322314  0.07804878\n",
      " 0.0440468  0.0240692  0.01285076 0.00674661 0.09423024]\n"
     ]
    }
   ],
   "source": [
    "data_book = ['l','l','l','l','c','c','c','c','c','c']\n",
    "data_book = ['l','l','l','l','l','l','l','l','l','c']\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "posterior = np.zeros(5)\n",
    "posteriors = np.zeros((11, 5))\n",
    "posteriors[0] = prior\n",
    "\n",
    "for n in range(1,11): \n",
    "    for (i,h) in enumerate([h1,h2,h3,h4,h5]):\n",
    "        posterior[i] = (prior[i] * likelihood(data_book[:n], h))\n",
    "    posterior /= np.sum(posterior)\n",
    "    posteriors[n] = posterior\n",
    "\n",
    "    \n",
    "x = np.arange(0,11)  \n",
    "y1 = posteriors[:,0]\n",
    "y2 = posteriors[:,1]\n",
    "y3 = posteriors[:,2]\n",
    "y4 = posteriors[:,3]\n",
    "y5 = posteriors[:,4]\n",
    "\n",
    "print(x)\n",
    "print(y3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to create a figure simular to 20.1 from the textbook that plots the posterior probabilities as a function of the number of observatons for each hypothesis. Each colored line corresponds to one of the five hypotheses (bag types of candy). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "        .bk-notebook-logo {\n",
       "            display: block;\n",
       "            width: 20px;\n",
       "            height: 20px;\n",
       "            background-image: url(data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABQAAAAUCAYAAACNiR0NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAABx0RVh0U29mdHdhcmUAQWRvYmUgRmlyZXdvcmtzIENTNui8sowAAAOkSURBVDiNjZRtaJVlGMd/1/08zzln5zjP1LWcU9N0NkN8m2CYjpgQYQXqSs0I84OLIC0hkEKoPtiH3gmKoiJDU7QpLgoLjLIQCpEsNJ1vqUOdO7ppbuec5+V+rj4ctwzd8IIbbi6u+8f1539dt3A78eXC7QizUF7gyV1fD1Yqg4JWz84yffhm0qkFqBogB9rM8tZdtwVsPUhWhGcFJngGeWrPzHm5oaMmkfEg1usvLFyc8jLRqDOMru7AyC8saQr7GG7f5fvDeH7Ej8CM66nIF+8yngt6HWaKh7k49Soy9nXurCi1o3qUbS3zWfrYeQDTB/Qj6kX6Ybhw4B+bOYoLKCC9H3Nu/leUTZ1JdRWkkn2ldcCamzrcf47KKXdAJllSlxAOkRgyHsGC/zRday5Qld9DyoM4/q/rUoy/CXh3jzOu3bHUVZeU+DEn8FInkPBFlu3+nW3Nw0mk6vCDiWg8CeJaxEwuHS3+z5RgY+YBR6V1Z1nxSOfoaPa4LASWxxdNp+VWTk7+4vzaou8v8PN+xo+KY2xsw6une2frhw05CTYOmQvsEhjhWjn0bmXPjpE1+kplmmkP3suftwTubK9Vq22qKmrBhpY4jvd5afdRA3wGjFAgcnTK2s4hY0/GPNIb0nErGMCRxWOOX64Z8RAC4oCXdklmEvcL8o0BfkNK4lUg9HTl+oPlQxdNo3Mg4Nv175e/1LDGzZen30MEjRUtmXSfiTVu1kK8W4txyV6BMKlbgk3lMwYCiusNy9fVfvvwMxv8Ynl6vxoByANLTWplvuj/nF9m2+PDtt1eiHPBr1oIfhCChQMBw6Aw0UulqTKZdfVvfG7VcfIqLG9bcldL/+pdWTLxLUy8Qq38heUIjh4XlzZxzQm19lLFlr8vdQ97rjZVOLf8nclzckbcD4wxXMidpX30sFd37Fv/GtwwhzhxGVAprjbg0gCAEeIgwCZyTV2Z1REEW8O4py0wsjeloKoMr6iCY6dP92H6Vw/oTyICIthibxjm/DfN9lVz8IqtqKYLUXfoKVMVQVVJOElGjrnnUt9T9wbgp8AyYKaGlqingHZU/uG2NTZSVqwHQTWkx9hxjkpWDaCg6Ckj5qebgBVbT3V3NNXMSiWSDdGV3hrtzla7J+duwPOToIg42ChPQOQjspnSlp1V+Gjdged7+8UN5CRAV7a5EdFNwCjEaBR27b3W890TE7g24NAP/mMDXRWrGoFPQI9ls/MWO2dWFAar/xcOIImbbpA3zgAAAABJRU5ErkJggg==);\n",
       "        }\n",
       "    </style>\n",
       "    <div>\n",
       "        <a href=\"https://bokeh.org\" target=\"_blank\" class=\"bk-notebook-logo\"></a>\n",
       "        <span id=\"b7a3933f-d58b-4a87-b20c-759f916bf692\">Loading BokehJS ...</span>\n",
       "    </div>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "(function(root) {\n",
       "  function now() {\n",
       "    return new Date();\n",
       "  }\n",
       "\n",
       "  const force = true;\n",
       "\n",
       "  if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n",
       "    root._bokeh_onload_callbacks = [];\n",
       "    root._bokeh_is_loading = undefined;\n",
       "  }\n",
       "\n",
       "const JS_MIME_TYPE = 'application/javascript';\n",
       "  const HTML_MIME_TYPE = 'text/html';\n",
       "  const EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n",
       "  const CLASS_NAME = 'output_bokeh rendered_html';\n",
       "\n",
       "  /**\n",
       "   * Render data to the DOM node\n",
       "   */\n",
       "  function render(props, node) {\n",
       "    const script = document.createElement(\"script\");\n",
       "    node.appendChild(script);\n",
       "  }\n",
       "\n",
       "  /**\n",
       "   * Handle when an output is cleared or removed\n",
       "   */\n",
       "  function handleClearOutput(event, handle) {\n",
       "    const cell = handle.cell;\n",
       "\n",
       "    const id = cell.output_area._bokeh_element_id;\n",
       "    const server_id = cell.output_area._bokeh_server_id;\n",
       "    // Clean up Bokeh references\n",
       "    if (id != null && id in Bokeh.index) {\n",
       "      Bokeh.index[id].model.document.clear();\n",
       "      delete Bokeh.index[id];\n",
       "    }\n",
       "\n",
       "    if (server_id !== undefined) {\n",
       "      // Clean up Bokeh references\n",
       "      const cmd_clean = \"from bokeh.io.state import curstate; print(curstate().uuid_to_server['\" + server_id + \"'].get_sessions()[0].document.roots[0]._id)\";\n",
       "      cell.notebook.kernel.execute(cmd_clean, {\n",
       "        iopub: {\n",
       "          output: function(msg) {\n",
       "            const id = msg.content.text.trim();\n",
       "            if (id in Bokeh.index) {\n",
       "              Bokeh.index[id].model.document.clear();\n",
       "              delete Bokeh.index[id];\n",
       "            }\n",
       "          }\n",
       "        }\n",
       "      });\n",
       "      // Destroy server and session\n",
       "      const cmd_destroy = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n",
       "      cell.notebook.kernel.execute(cmd_destroy);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  /**\n",
       "   * Handle when a new output is added\n",
       "   */\n",
       "  function handleAddOutput(event, handle) {\n",
       "    const output_area = handle.output_area;\n",
       "    const output = handle.output;\n",
       "\n",
       "    // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n",
       "    if ((output.output_type != \"display_data\") || (!Object.prototype.hasOwnProperty.call(output.data, EXEC_MIME_TYPE))) {\n",
       "      return\n",
       "    }\n",
       "\n",
       "    const toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n",
       "\n",
       "    if (output.metadata[EXEC_MIME_TYPE][\"id\"] !== undefined) {\n",
       "      toinsert[toinsert.length - 1].firstChild.textContent = output.data[JS_MIME_TYPE];\n",
       "      // store reference to embed id on output_area\n",
       "      output_area._bokeh_element_id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n",
       "    }\n",
       "    if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n",
       "      const bk_div = document.createElement(\"div\");\n",
       "      bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n",
       "      const script_attrs = bk_div.children[0].attributes;\n",
       "      for (let i = 0; i < script_attrs.length; i++) {\n",
       "        toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\n",
       "        toinsert[toinsert.length - 1].firstChild.textContent = bk_div.children[0].textContent\n",
       "      }\n",
       "      // store reference to server id on output_area\n",
       "      output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function register_renderer(events, OutputArea) {\n",
       "\n",
       "    function append_mime(data, metadata, element) {\n",
       "      // create a DOM node to render to\n",
       "      const toinsert = this.create_output_subarea(\n",
       "        metadata,\n",
       "        CLASS_NAME,\n",
       "        EXEC_MIME_TYPE\n",
       "      );\n",
       "      this.keyboard_manager.register_events(toinsert);\n",
       "      // Render to node\n",
       "      const props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n",
       "      render(props, toinsert[toinsert.length - 1]);\n",
       "      element.append(toinsert);\n",
       "      return toinsert\n",
       "    }\n",
       "\n",
       "    /* Handle when an output is cleared or removed */\n",
       "    events.on('clear_output.CodeCell', handleClearOutput);\n",
       "    events.on('delete.Cell', handleClearOutput);\n",
       "\n",
       "    /* Handle when a new output is added */\n",
       "    events.on('output_added.OutputArea', handleAddOutput);\n",
       "\n",
       "    /**\n",
       "     * Register the mime type and append_mime function with output_area\n",
       "     */\n",
       "    OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n",
       "      /* Is output safe? */\n",
       "      safe: true,\n",
       "      /* Index of renderer in `output_area.display_order` */\n",
       "      index: 0\n",
       "    });\n",
       "  }\n",
       "\n",
       "  // register the mime type if in Jupyter Notebook environment and previously unregistered\n",
       "  if (root.Jupyter !== undefined) {\n",
       "    const events = require('base/js/events');\n",
       "    const OutputArea = require('notebook/js/outputarea').OutputArea;\n",
       "\n",
       "    if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n",
       "      register_renderer(events, OutputArea);\n",
       "    }\n",
       "  }\n",
       "  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n",
       "    root._bokeh_timeout = Date.now() + 5000;\n",
       "    root._bokeh_failed_load = false;\n",
       "  }\n",
       "\n",
       "  const NB_LOAD_WARNING = {'data': {'text/html':\n",
       "     \"<div style='background-color: #fdd'>\\n\"+\n",
       "     \"<p>\\n\"+\n",
       "     \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n",
       "     \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n",
       "     \"</p>\\n\"+\n",
       "     \"<ul>\\n\"+\n",
       "     \"<li>re-rerun `output_notebook()` to attempt to load from CDN again, or</li>\\n\"+\n",
       "     \"<li>use INLINE resources instead, as so:</li>\\n\"+\n",
       "     \"</ul>\\n\"+\n",
       "     \"<code>\\n\"+\n",
       "     \"from bokeh.resources import INLINE\\n\"+\n",
       "     \"output_notebook(resources=INLINE)\\n\"+\n",
       "     \"</code>\\n\"+\n",
       "     \"</div>\"}};\n",
       "\n",
       "  function display_loaded() {\n",
       "    const el = document.getElementById(\"b7a3933f-d58b-4a87-b20c-759f916bf692\");\n",
       "    if (el != null) {\n",
       "      el.textContent = \"BokehJS is loading...\";\n",
       "    }\n",
       "    if (root.Bokeh !== undefined) {\n",
       "      if (el != null) {\n",
       "        el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n",
       "      }\n",
       "    } else if (Date.now() < root._bokeh_timeout) {\n",
       "      setTimeout(display_loaded, 100)\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function run_callbacks() {\n",
       "    try {\n",
       "      root._bokeh_onload_callbacks.forEach(function(callback) {\n",
       "        if (callback != null)\n",
       "          callback();\n",
       "      });\n",
       "    } finally {\n",
       "      delete root._bokeh_onload_callbacks\n",
       "    }\n",
       "    console.debug(\"Bokeh: all callbacks have finished\");\n",
       "  }\n",
       "\n",
       "  function load_libs(css_urls, js_urls, callback) {\n",
       "    if (css_urls == null) css_urls = [];\n",
       "    if (js_urls == null) js_urls = [];\n",
       "\n",
       "    root._bokeh_onload_callbacks.push(callback);\n",
       "    if (root._bokeh_is_loading > 0) {\n",
       "      console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n",
       "      return null;\n",
       "    }\n",
       "    if (js_urls == null || js_urls.length === 0) {\n",
       "      run_callbacks();\n",
       "      return null;\n",
       "    }\n",
       "    console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n",
       "    root._bokeh_is_loading = css_urls.length + js_urls.length;\n",
       "\n",
       "    function on_load() {\n",
       "      root._bokeh_is_loading--;\n",
       "      if (root._bokeh_is_loading === 0) {\n",
       "        console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n",
       "        run_callbacks()\n",
       "      }\n",
       "    }\n",
       "\n",
       "    function on_error(url) {\n",
       "      console.error(\"failed to load \" + url);\n",
       "    }\n",
       "\n",
       "    for (let i = 0; i < css_urls.length; i++) {\n",
       "      const url = css_urls[i];\n",
       "      const element = document.createElement(\"link\");\n",
       "      element.onload = on_load;\n",
       "      element.onerror = on_error.bind(null, url);\n",
       "      element.rel = \"stylesheet\";\n",
       "      element.type = \"text/css\";\n",
       "      element.href = url;\n",
       "      console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n",
       "      document.body.appendChild(element);\n",
       "    }\n",
       "\n",
       "    for (let i = 0; i < js_urls.length; i++) {\n",
       "      const url = js_urls[i];\n",
       "      const element = document.createElement('script');\n",
       "      element.onload = on_load;\n",
       "      element.onerror = on_error.bind(null, url);\n",
       "      element.async = false;\n",
       "      element.src = url;\n",
       "      console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n",
       "      document.head.appendChild(element);\n",
       "    }\n",
       "  };\n",
       "\n",
       "  function inject_raw_css(css) {\n",
       "    const element = document.createElement(\"style\");\n",
       "    element.appendChild(document.createTextNode(css));\n",
       "    document.body.appendChild(element);\n",
       "  }\n",
       "\n",
       "  const js_urls = [\"https://cdn.bokeh.org/bokeh/release/bokeh-3.1.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.1.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.1.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.1.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-mathjax-3.1.1.min.js\"];\n",
       "  const css_urls = [];\n",
       "\n",
       "  const inline_js = [    function(Bokeh) {\n",
       "      Bokeh.set_log_level(\"info\");\n",
       "    },\n",
       "function(Bokeh) {\n",
       "    }\n",
       "  ];\n",
       "\n",
       "  function run_inline_js() {\n",
       "    if (root.Bokeh !== undefined || force === true) {\n",
       "          for (let i = 0; i < inline_js.length; i++) {\n",
       "      inline_js[i].call(root, root.Bokeh);\n",
       "    }\n",
       "if (force === true) {\n",
       "        display_loaded();\n",
       "      }} else if (Date.now() < root._bokeh_timeout) {\n",
       "      setTimeout(run_inline_js, 100);\n",
       "    } else if (!root._bokeh_failed_load) {\n",
       "      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n",
       "      root._bokeh_failed_load = true;\n",
       "    } else if (force !== true) {\n",
       "      const cell = $(document.getElementById(\"b7a3933f-d58b-4a87-b20c-759f916bf692\")).parents('.cell').data().cell;\n",
       "      cell.output_area.append_execute_result(NB_LOAD_WARNING)\n",
       "    }\n",
       "  }\n",
       "\n",
       "  if (root._bokeh_is_loading === 0) {\n",
       "    console.debug(\"Bokeh: BokehJS loaded, going straight to plotting\");\n",
       "    run_inline_js();\n",
       "  } else {\n",
       "    load_libs(css_urls, js_urls, function() {\n",
       "      console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n",
       "      run_inline_js();\n",
       "    });\n",
       "  }\n",
       "}(window));"
      ],
      "application/vnd.bokehjs_load.v0+json": "(function(root) {\n  function now() {\n    return new Date();\n  }\n\n  const force = true;\n\n  if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n    root._bokeh_onload_callbacks = [];\n    root._bokeh_is_loading = undefined;\n  }\n\n\n  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n    root._bokeh_timeout = Date.now() + 5000;\n    root._bokeh_failed_load = false;\n  }\n\n  const NB_LOAD_WARNING = {'data': {'text/html':\n     \"<div style='background-color: #fdd'>\\n\"+\n     \"<p>\\n\"+\n     \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n     \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n     \"</p>\\n\"+\n     \"<ul>\\n\"+\n     \"<li>re-rerun `output_notebook()` to attempt to load from CDN again, or</li>\\n\"+\n     \"<li>use INLINE resources instead, as so:</li>\\n\"+\n     \"</ul>\\n\"+\n     \"<code>\\n\"+\n     \"from bokeh.resources import INLINE\\n\"+\n     \"output_notebook(resources=INLINE)\\n\"+\n     \"</code>\\n\"+\n     \"</div>\"}};\n\n  function display_loaded() {\n    const el = document.getElementById(\"b7a3933f-d58b-4a87-b20c-759f916bf692\");\n    if (el != null) {\n      el.textContent = \"BokehJS is loading...\";\n    }\n    if (root.Bokeh !== undefined) {\n      if (el != null) {\n        el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n      }\n    } else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(display_loaded, 100)\n    }\n  }\n\n  function run_callbacks() {\n    try {\n      root._bokeh_onload_callbacks.forEach(function(callback) {\n        if (callback != null)\n          callback();\n      });\n    } finally {\n      delete root._bokeh_onload_callbacks\n    }\n    console.debug(\"Bokeh: all callbacks have finished\");\n  }\n\n  function load_libs(css_urls, js_urls, callback) {\n    if (css_urls == null) css_urls = [];\n    if (js_urls == null) js_urls = [];\n\n    root._bokeh_onload_callbacks.push(callback);\n    if (root._bokeh_is_loading > 0) {\n      console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n      return null;\n    }\n    if (js_urls == null || js_urls.length === 0) {\n      run_callbacks();\n      return null;\n    }\n    console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n    root._bokeh_is_loading = css_urls.length + js_urls.length;\n\n    function on_load() {\n      root._bokeh_is_loading--;\n      if (root._bokeh_is_loading === 0) {\n        console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n        run_callbacks()\n      }\n    }\n\n    function on_error(url) {\n      console.error(\"failed to load \" + url);\n    }\n\n    for (let i = 0; i < css_urls.length; i++) {\n      const url = css_urls[i];\n      const element = document.createElement(\"link\");\n      element.onload = on_load;\n      element.onerror = on_error.bind(null, url);\n      element.rel = \"stylesheet\";\n      element.type = \"text/css\";\n      element.href = url;\n      console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n      document.body.appendChild(element);\n    }\n\n    for (let i = 0; i < js_urls.length; i++) {\n      const url = js_urls[i];\n      const element = document.createElement('script');\n      element.onload = on_load;\n      element.onerror = on_error.bind(null, url);\n      element.async = false;\n      element.src = url;\n      console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n      document.head.appendChild(element);\n    }\n  };\n\n  function inject_raw_css(css) {\n    const element = document.createElement(\"style\");\n    element.appendChild(document.createTextNode(css));\n    document.body.appendChild(element);\n  }\n\n  const js_urls = [\"https://cdn.bokeh.org/bokeh/release/bokeh-3.1.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.1.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.1.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.1.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-mathjax-3.1.1.min.js\"];\n  const css_urls = [];\n\n  const inline_js = [    function(Bokeh) {\n      Bokeh.set_log_level(\"info\");\n    },\nfunction(Bokeh) {\n    }\n  ];\n\n  function run_inline_js() {\n    if (root.Bokeh !== undefined || force === true) {\n          for (let i = 0; i < inline_js.length; i++) {\n      inline_js[i].call(root, root.Bokeh);\n    }\nif (force === true) {\n        display_loaded();\n      }} else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(run_inline_js, 100);\n    } else if (!root._bokeh_failed_load) {\n      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n      root._bokeh_failed_load = true;\n    } else if (force !== true) {\n      const cell = $(document.getElementById(\"b7a3933f-d58b-4a87-b20c-759f916bf692\")).parents('.cell').data().cell;\n      cell.output_area.append_execute_result(NB_LOAD_WARNING)\n    }\n  }\n\n  if (root._bokeh_is_loading === 0) {\n    console.debug(\"Bokeh: BokehJS loaded, going straight to plotting\");\n    run_inline_js();\n  } else {\n    load_libs(css_urls, js_urls, function() {\n      console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n      run_inline_js();\n    });\n  }\n}(window));"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "  <div id=\"e4a9b01e-acfb-4409-9598-6c2aeac37a08\" data-root-id=\"p1001\" style=\"display: contents;\"></div>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "(function(root) {\n",
       "  function embed_document(root) {\n",
       "  const docs_json = {\"94d2e61d-8368-4d7e-b498-0607a0971987\":{\"version\":\"3.1.1\",\"title\":\"Bokeh Application\",\"defs\":[],\"roots\":[{\"type\":\"object\",\"name\":\"Figure\",\"id\":\"p1001\",\"attributes\":{\"x_range\":{\"type\":\"object\",\"name\":\"DataRange1d\",\"id\":\"p1002\"},\"y_range\":{\"type\":\"object\",\"name\":\"DataRange1d\",\"id\":\"p1003\"},\"x_scale\":{\"type\":\"object\",\"name\":\"LinearScale\",\"id\":\"p1014\"},\"y_scale\":{\"type\":\"object\",\"name\":\"LinearScale\",\"id\":\"p1016\"},\"title\":{\"type\":\"object\",\"name\":\"Title\",\"id\":\"p1005\"},\"renderers\":[{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1053\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1047\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1049\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1048\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAEAAAACAAAAAwAAAAQAAAAFAAAABgAAAAcAAAAIAAAACQAAAAoAAAA=\"},\"shape\":[11],\"dtype\":\"int32\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"mpmZmZmZuT8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA==\"},\"shape\":[11],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1054\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1055\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1050\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"red\",\"line_width\":2}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1051\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"red\",\"line_alpha\":0.1,\"line_width\":2}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1052\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"red\",\"line_alpha\":0.2,\"line_width\":2}}}},{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1062\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1056\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1058\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1057\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAEAAAACAAAAAwAAAAQAAAAFAAAABgAAAAcAAAAIAAAACQAAAAoAAAA=\"},\"shape\":[11],\"dtype\":\"int32\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"mpmZmZmZuT8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA==\"},\"shape\":[11],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1063\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1064\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1059\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"red\"},\"fill_color\":{\"type\":\"value\",\"value\":\"red\"},\"hatch_color\":{\"type\":\"value\",\"value\":\"red\"}}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1060\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"red\"},\"line_alpha\":{\"type\":\"value\",\"value\":0.1},\"fill_color\":{\"type\":\"value\",\"value\":\"red\"},\"fill_alpha\":{\"type\":\"value\",\"value\":0.1},\"hatch_color\":{\"type\":\"value\",\"value\":\"red\"},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.1}}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1061\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"red\"},\"line_alpha\":{\"type\":\"value\",\"value\":0.2},\"fill_color\":{\"type\":\"value\",\"value\":\"red\"},\"fill_alpha\":{\"type\":\"value\",\"value\":0.2},\"hatch_color\":{\"type\":\"value\",\"value\":\"red\"},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.2}}}}},{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1071\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1065\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1067\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1066\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAEAAAACAAAAAwAAAAQAAAAFAAAABgAAAAcAAAAIAAAACQAAAAoAAAA=\"},\"shape\":[11],\"dtype\":\"int32\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"mpmZmZmZyT+amZmZmZm5PxQ7sRM7saM/KK+hvIbyij+QnOZr9exwP/sTsD8B+1M/8kWfWE2NNj/5zEWZmKUYPyU4v+B/Ufo+My8MglSi2z4YG/bLmhciPw==\"},\"shape\":[11],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1072\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1073\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1068\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"green\",\"line_width\":2}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1069\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"green\",\"line_alpha\":0.1,\"line_width\":2}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1070\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"green\",\"line_alpha\":0.2,\"line_width\":2}}}},{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1080\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1074\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1076\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1075\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAEAAAACAAAAAwAAAAQAAAAFAAAABgAAAAcAAAAIAAAACQAAAAoAAAA=\"},\"shape\":[11],\"dtype\":\"int32\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"mpmZmZmZyT+amZmZmZm5PxQ7sRM7saM/KK+hvIbyij+QnOZr9exwP/sTsD8B+1M/8kWfWE2NNj/5zEWZmKUYPyU4v+B/Ufo+My8MglSi2z4YG/bLmhciPw==\"},\"shape\":[11],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1081\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1082\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1077\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"green\"},\"fill_color\":{\"type\":\"value\",\"value\":\"green\"},\"hatch_color\":{\"type\":\"value\",\"value\":\"green\"}}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1078\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"green\"},\"line_alpha\":{\"type\":\"value\",\"value\":0.1},\"fill_color\":{\"type\":\"value\",\"value\":\"green\"},\"fill_alpha\":{\"type\":\"value\",\"value\":0.1},\"hatch_color\":{\"type\":\"value\",\"value\":\"green\"},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.1}}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1079\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"green\"},\"line_alpha\":{\"type\":\"value\",\"value\":0.2},\"fill_color\":{\"type\":\"value\",\"value\":\"green\"},\"fill_alpha\":{\"type\":\"value\",\"value\":0.2},\"hatch_color\":{\"type\":\"value\",\"value\":\"green\"},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.2}}}}},{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1089\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1083\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1085\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1084\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAEAAAACAAAAAwAAAAQAAAAFAAAABgAAAAcAAAAIAAAACQAAAAoAAAA=\"},\"shape\":[11],\"dtype\":\"int32\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"mpmZmZmZ2T+amZmZmZnZPxQ7sRM7sdM/KK+hvIbyyj+QnOZr9ezAP/sTsD8B+7M/8kWfWE2Npj/5zEWZmKWYPyU4v+B/UYo/My8MglSiez/KzvIPeR+4Pw==\"},\"shape\":[11],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1090\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1091\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1086\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"blue\",\"line_width\":2}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1087\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"blue\",\"line_alpha\":0.1,\"line_width\":2}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1088\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"blue\",\"line_alpha\":0.2,\"line_width\":2}}}},{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1098\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1092\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1094\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1093\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAEAAAACAAAAAwAAAAQAAAAFAAAABgAAAAcAAAAIAAAACQAAAAoAAAA=\"},\"shape\":[11],\"dtype\":\"int32\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"mpmZmZmZ2T+amZmZmZnZPxQ7sRM7sdM/KK+hvIbyyj+QnOZr9ezAP/sTsD8B+7M/8kWfWE2Npj/5zEWZmKWYPyU4v+B/UYo/My8MglSiez/KzvIPeR+4Pw==\"},\"shape\":[11],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1099\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1100\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1095\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"blue\"},\"fill_color\":{\"type\":\"value\",\"value\":\"blue\"},\"hatch_color\":{\"type\":\"value\",\"value\":\"blue\"}}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1096\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"blue\"},\"line_alpha\":{\"type\":\"value\",\"value\":0.1},\"fill_color\":{\"type\":\"value\",\"value\":\"blue\"},\"fill_alpha\":{\"type\":\"value\",\"value\":0.1},\"hatch_color\":{\"type\":\"value\",\"value\":\"blue\"},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.1}}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1097\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"blue\"},\"line_alpha\":{\"type\":\"value\",\"value\":0.2},\"fill_color\":{\"type\":\"value\",\"value\":\"blue\"},\"fill_alpha\":{\"type\":\"value\",\"value\":0.2},\"hatch_color\":{\"type\":\"value\",\"value\":\"blue\"},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.2}}}}},{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1107\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1101\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1103\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1102\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAEAAAACAAAAAwAAAAQAAAAFAAAABgAAAAcAAAAIAAAACQAAAAoAAAA=\"},\"shape\":[11],\"dtype\":\"int32\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"mpmZmZmZyT80MzMzMzPTP3ZiJ3ZiJ9Y/ymsor6G81j8l3o+c5mvVP/cidC9C99I/jGNXUBgO0D/C6azU1VHKP+8l+PIVFMU/ot0VD2SZwD/F5lRk7/rsPw==\"},\"shape\":[11],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1108\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1109\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1104\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"yellow\",\"line_width\":2}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1105\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"yellow\",\"line_alpha\":0.1,\"line_width\":2}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1106\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"yellow\",\"line_alpha\":0.2,\"line_width\":2}}}},{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1116\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1110\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1112\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1111\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAEAAAACAAAAAwAAAAQAAAAFAAAABgAAAAcAAAAIAAAACQAAAAoAAAA=\"},\"shape\":[11],\"dtype\":\"int32\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"mpmZmZmZyT80MzMzMzPTP3ZiJ3ZiJ9Y/ymsor6G81j8l3o+c5mvVP/cidC9C99I/jGNXUBgO0D/C6azU1VHKP+8l+PIVFMU/ot0VD2SZwD/F5lRk7/rsPw==\"},\"shape\":[11],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1117\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1118\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1113\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"yellow\"},\"fill_color\":{\"type\":\"value\",\"value\":\"yellow\"},\"hatch_color\":{\"type\":\"value\",\"value\":\"yellow\"}}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1114\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"yellow\"},\"line_alpha\":{\"type\":\"value\",\"value\":0.1},\"fill_color\":{\"type\":\"value\",\"value\":\"yellow\"},\"fill_alpha\":{\"type\":\"value\",\"value\":0.1},\"hatch_color\":{\"type\":\"value\",\"value\":\"yellow\"},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.1}}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1115\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"yellow\"},\"line_alpha\":{\"type\":\"value\",\"value\":0.2},\"fill_color\":{\"type\":\"value\",\"value\":\"yellow\"},\"fill_alpha\":{\"type\":\"value\",\"value\":0.2},\"hatch_color\":{\"type\":\"value\",\"value\":\"yellow\"},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.2}}}}},{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1125\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1119\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1121\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1120\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAEAAAACAAAAAwAAAAQAAAAFAAAABgAAAAcAAAAIAAAACQAAAAoAAAA=\"},\"shape\":[11],\"dtype\":\"int32\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"mpmZmZmZuT+amZmZmZnJPxQ7sRM7sdM/KK+hvIby2j+QnOZr9ezgP/sTsD8B++M/8kWfWE2N5j/5zEWZmKXoPyU4v+B/Ueo/My8MglSi6z8AAAAAAAAAAA==\"},\"shape\":[11],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1126\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1127\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1122\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"orange\",\"line_width\":2}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1123\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"orange\",\"line_alpha\":0.1,\"line_width\":2}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1124\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"orange\",\"line_alpha\":0.2,\"line_width\":2}}}},{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1134\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1128\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1130\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1129\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAEAAAACAAAAAwAAAAQAAAAFAAAABgAAAAcAAAAIAAAACQAAAAoAAAA=\"},\"shape\":[11],\"dtype\":\"int32\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"mpmZmZmZuT+amZmZmZnJPxQ7sRM7sdM/KK+hvIby2j+QnOZr9ezgP/sTsD8B++M/8kWfWE2N5j/5zEWZmKXoPyU4v+B/Ueo/My8MglSi6z8AAAAAAAAAAA==\"},\"shape\":[11],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1135\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1136\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1131\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"orange\"},\"fill_color\":{\"type\":\"value\",\"value\":\"orange\"},\"hatch_color\":{\"type\":\"value\",\"value\":\"orange\"}}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1132\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"orange\"},\"line_alpha\":{\"type\":\"value\",\"value\":0.1},\"fill_color\":{\"type\":\"value\",\"value\":\"orange\"},\"fill_alpha\":{\"type\":\"value\",\"value\":0.1},\"hatch_color\":{\"type\":\"value\",\"value\":\"orange\"},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.1}}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1133\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"orange\"},\"line_alpha\":{\"type\":\"value\",\"value\":0.2},\"fill_color\":{\"type\":\"value\",\"value\":\"orange\"},\"fill_alpha\":{\"type\":\"value\",\"value\":0.2},\"hatch_color\":{\"type\":\"value\",\"value\":\"orange\"},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.2}}}}}],\"toolbar\":{\"type\":\"object\",\"name\":\"Toolbar\",\"id\":\"p1008\",\"attributes\":{\"tools\":[{\"type\":\"object\",\"name\":\"PanTool\",\"id\":\"p1032\"},{\"type\":\"object\",\"name\":\"WheelZoomTool\",\"id\":\"p1033\"},{\"type\":\"object\",\"name\":\"BoxZoomTool\",\"id\":\"p1034\",\"attributes\":{\"overlay\":{\"type\":\"object\",\"name\":\"BoxAnnotation\",\"id\":\"p1035\",\"attributes\":{\"syncable\":false,\"level\":\"overlay\",\"visible\":false,\"left_units\":\"canvas\",\"right_units\":\"canvas\",\"bottom_units\":\"canvas\",\"top_units\":\"canvas\",\"line_color\":\"black\",\"line_alpha\":1.0,\"line_width\":2,\"line_dash\":[4,4],\"fill_color\":\"lightgrey\",\"fill_alpha\":0.5}}}},{\"type\":\"object\",\"name\":\"SaveTool\",\"id\":\"p1036\"},{\"type\":\"object\",\"name\":\"ResetTool\",\"id\":\"p1037\"},{\"type\":\"object\",\"name\":\"HelpTool\",\"id\":\"p1038\"}]}},\"left\":[{\"type\":\"object\",\"name\":\"LinearAxis\",\"id\":\"p1025\",\"attributes\":{\"ticker\":{\"type\":\"object\",\"name\":\"BasicTicker\",\"id\":\"p1027\",\"attributes\":{\"mantissas\":[1,2,5]}},\"formatter\":{\"type\":\"object\",\"name\":\"BasicTickFormatter\",\"id\":\"p1026\"},\"major_label_policy\":{\"type\":\"object\",\"name\":\"AllLabels\",\"id\":\"p1028\"}}}],\"below\":[{\"type\":\"object\",\"name\":\"LinearAxis\",\"id\":\"p1018\",\"attributes\":{\"ticker\":{\"type\":\"object\",\"name\":\"BasicTicker\",\"id\":\"p1020\",\"attributes\":{\"mantissas\":[1,2,5]}},\"formatter\":{\"type\":\"object\",\"name\":\"BasicTickFormatter\",\"id\":\"p1019\"},\"major_label_policy\":{\"type\":\"object\",\"name\":\"AllLabels\",\"id\":\"p1021\"}}}],\"center\":[{\"type\":\"object\",\"name\":\"Grid\",\"id\":\"p1024\",\"attributes\":{\"axis\":{\"id\":\"p1018\"}}},{\"type\":\"object\",\"name\":\"Grid\",\"id\":\"p1031\",\"attributes\":{\"dimension\":1,\"axis\":{\"id\":\"p1025\"}}}]}}],\"callbacks\":{\"type\":\"map\"}}};\n",
       "  const render_items = [{\"docid\":\"94d2e61d-8368-4d7e-b498-0607a0971987\",\"roots\":{\"p1001\":\"e4a9b01e-acfb-4409-9598-6c2aeac37a08\"},\"root_ids\":[\"p1001\"]}];\n",
       "  root.Bokeh.embed.embed_items_notebook(docs_json, render_items);\n",
       "  }\n",
       "  if (root.Bokeh !== undefined) {\n",
       "    embed_document(root);\n",
       "  } else {\n",
       "    let attempts = 0;\n",
       "    const timer = setInterval(function(root) {\n",
       "      if (root.Bokeh !== undefined) {\n",
       "        clearInterval(timer);\n",
       "        embed_document(root);\n",
       "      } else {\n",
       "        attempts++;\n",
       "        if (attempts > 100) {\n",
       "          clearInterval(timer);\n",
       "          console.log(\"Bokeh: ERROR: Unable to run BokehJS code because BokehJS library is missing\");\n",
       "        }\n",
       "      }\n",
       "    }, 10, root)\n",
       "  }\n",
       "})(window);"
      ],
      "application/vnd.bokehjs_exec.v0+json": ""
     },
     "metadata": {
      "application/vnd.bokehjs_exec.v0+json": {
       "id": "p1001"
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "from bokeh.plotting import figure, show \n",
    "from bokeh.io import output_notebook \n",
    "\n",
    "output_notebook()\n",
    "p = figure() \n",
    "# p.legend.title = 'Bayesian Learning'\n",
    "\n",
    "for (y,label,color) in zip([y1,y2,y3,y4,y5], \n",
    "                     ['P(h1|d)', 'P(h2|d)','P(h3|d)', 'P(h4|d)', 'P(h5/d)'], \n",
    "                     ['red','green','blue','yellow', 'orange']): \n",
    "    p.line(x, y, line_width=2,color=color)\n",
    "    p.circle(x,y, color = color)\n",
    "\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot above shows how the posterior probabilities of each hypothesis get updated as we receive more data. To make a prediction \n",
    "we consider all hypotheses and weight them by their posterior probabilities. \n",
    "\n",
    "Let's first consider how to do a prediction without observing any data. Let's say we want the probability that the next candy is lime. \n",
    "If we only had one hypothesis we can directly use the corresponding probability. For example if we have an h3 bag then the probability the next candy is lime is 0.5 - if we have an h5 bag the probability the next candy is like is 1.0. \n",
    "\n",
    "Now if we have multiple hypothesis and their prior we simply weigh them by the prior. So for $h1, h2, h3, h4, h5$ and the prior \n",
    "$P(h1) = 0.1, P(h2 = 0.2) P(h3 = 0.4) P(h4 = 0.2) P(h1 = 0.1)$\n",
    "we would have: \n",
    " \n",
    "$ 0.1 * 0 + 0.25 * 0.2 + 0.5 * 0.4 + 0.75 * 0.2 + 0.1 * 1 = 0.5$\n",
    "\n",
    "If we know that d1 is lime then we use the posterior probabilities \n",
    "after d1 that we calculated above to weigh the probabilities. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posteriors: [0.1 0.2 0.4 0.2 0.1]\n",
      "Lime probabilities for each hypotheses: [0.0, 0.25, 0.5, 0.75, 1.0]\n",
      "P(next=lime): 0.5\n",
      "Posteriors: [0.  0.1 0.4 0.3 0.2]\n",
      "Lime probabilities for each hypotheses: [0.0, 0.25, 0.5, 0.75, 1.0]\n",
      "P(next=lime): 0.65\n"
     ]
    }
   ],
   "source": [
    "lime_probs = [h1[1], h2[1], h3[1], h4[1], h5[1]]\n",
    "pd0 = posteriors[0,:]\n",
    "print('Posteriors:',pd0)\n",
    "print('Lime probabilities for each hypotheses:', lime_probs)\n",
    "print('P(next=lime):', np.dot(pd0, lime_probs))\n",
    "\n",
    "pd1 = posteriors[1,:]\n",
    "print('Posteriors:',pd1)\n",
    "print('Lime probabilities for each hypotheses:', lime_probs)\n",
    "print('P(next=lime):', np.dot(pd1, lime_probs))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.5        0.65       0.73076923 0.79605263 0.84710744 0.88597561\n",
      " 0.91500344 0.93648928 0.95238693 0.96420148 0.72637342]\n"
     ]
    }
   ],
   "source": [
    "prob_next_lime_bayes = np.zeros(11)\n",
    "\n",
    "\n",
    "for n in range(0,11): \n",
    "    for (i,h) in enumerate([h1,h2,h3,h4,h5]):\n",
    "        if (n == 0):\n",
    "            posterior[i] = prior[i]\n",
    "        else: # n > 0\n",
    "            posterior[i] = (prior[i] * likelihood(data_book[:n], h))\n",
    "        lime_probs = [h1[1], h2[1], h3[1], h4[1], h5[1]]\n",
    "    posterior /= np.sum(posterior)\n",
    "    prob_next_lime_bayes[n] = np.dot(posterior, lime_probs)\n",
    "\n",
    "print(prob_next_lime_bayes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the probability that the next candy is lime in a way similar to Figure 20.1 of the textbook. \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "  <div id=\"f8ba568e-9834-4556-a61b-fb1a50cada8d\" data-root-id=\"p1526\" style=\"display: contents;\"></div>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "(function(root) {\n",
       "  function embed_document(root) {\n",
       "  const docs_json = {\"63de7880-012e-4c8e-81e5-7af7c3194325\":{\"version\":\"3.1.1\",\"title\":\"Bokeh Application\",\"defs\":[],\"roots\":[{\"type\":\"object\",\"name\":\"Figure\",\"id\":\"p1526\",\"attributes\":{\"x_range\":{\"type\":\"object\",\"name\":\"DataRange1d\",\"id\":\"p1527\"},\"y_range\":{\"type\":\"object\",\"name\":\"DataRange1d\",\"id\":\"p1528\"},\"x_scale\":{\"type\":\"object\",\"name\":\"LinearScale\",\"id\":\"p1539\"},\"y_scale\":{\"type\":\"object\",\"name\":\"LinearScale\",\"id\":\"p1541\"},\"title\":{\"type\":\"object\",\"name\":\"Title\",\"id\":\"p1530\"},\"renderers\":[{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1578\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1572\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1574\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1573\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAEAAAACAAAAAwAAAAQAAAAFAAAABgAAAAcAAAAIAAAACQAAAAoAAAA=\"},\"shape\":[11],\"dtype\":\"int32\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"/v//////3z/NzMzMzMzkPyh2Yid2Yuc/lNdQXkN56T++Vs8OgRvrP1mYnoXpWew/cMrvS7VH7T+OgrxeuPftPxIuzyb0ee4/M+WsDb3a7j/ZCap6cz7nPw==\"},\"shape\":[11],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1579\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1580\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1575\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"orange\",\"line_width\":2}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1576\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"orange\",\"line_alpha\":0.1,\"line_width\":2}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1577\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"orange\",\"line_alpha\":0.2,\"line_width\":2}}}},{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1587\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1581\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1583\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1582\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAEAAAACAAAAAwAAAAQAAAAFAAAABgAAAAcAAAAIAAAACQAAAAoAAAA=\"},\"shape\":[11],\"dtype\":\"int32\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"/v//////3z/NzMzMzMzkPyh2Yid2Yuc/lNdQXkN56T++Vs8OgRvrP1mYnoXpWew/cMrvS7VH7T+OgrxeuPftPxIuzyb0ee4/M+WsDb3a7j/ZCap6cz7nPw==\"},\"shape\":[11],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1588\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1589\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1584\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"orange\"},\"fill_color\":{\"type\":\"value\",\"value\":\"orange\"},\"hatch_color\":{\"type\":\"value\",\"value\":\"orange\"}}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1585\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"orange\"},\"line_alpha\":{\"type\":\"value\",\"value\":0.1},\"fill_color\":{\"type\":\"value\",\"value\":\"orange\"},\"fill_alpha\":{\"type\":\"value\",\"value\":0.1},\"hatch_color\":{\"type\":\"value\",\"value\":\"orange\"},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.1}}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Circle\",\"id\":\"p1586\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":{\"type\":\"value\",\"value\":\"orange\"},\"line_alpha\":{\"type\":\"value\",\"value\":0.2},\"fill_color\":{\"type\":\"value\",\"value\":\"orange\"},\"fill_alpha\":{\"type\":\"value\",\"value\":0.2},\"hatch_color\":{\"type\":\"value\",\"value\":\"orange\"},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.2}}}}}],\"toolbar\":{\"type\":\"object\",\"name\":\"Toolbar\",\"id\":\"p1533\",\"attributes\":{\"tools\":[{\"type\":\"object\",\"name\":\"PanTool\",\"id\":\"p1557\"},{\"type\":\"object\",\"name\":\"WheelZoomTool\",\"id\":\"p1558\"},{\"type\":\"object\",\"name\":\"BoxZoomTool\",\"id\":\"p1559\",\"attributes\":{\"overlay\":{\"type\":\"object\",\"name\":\"BoxAnnotation\",\"id\":\"p1560\",\"attributes\":{\"syncable\":false,\"level\":\"overlay\",\"visible\":false,\"left_units\":\"canvas\",\"right_units\":\"canvas\",\"bottom_units\":\"canvas\",\"top_units\":\"canvas\",\"line_color\":\"black\",\"line_alpha\":1.0,\"line_width\":2,\"line_dash\":[4,4],\"fill_color\":\"lightgrey\",\"fill_alpha\":0.5}}}},{\"type\":\"object\",\"name\":\"SaveTool\",\"id\":\"p1561\"},{\"type\":\"object\",\"name\":\"ResetTool\",\"id\":\"p1562\"},{\"type\":\"object\",\"name\":\"HelpTool\",\"id\":\"p1563\"}]}},\"left\":[{\"type\":\"object\",\"name\":\"LinearAxis\",\"id\":\"p1550\",\"attributes\":{\"ticker\":{\"type\":\"object\",\"name\":\"BasicTicker\",\"id\":\"p1552\",\"attributes\":{\"mantissas\":[1,2,5]}},\"formatter\":{\"type\":\"object\",\"name\":\"BasicTickFormatter\",\"id\":\"p1551\"},\"major_label_policy\":{\"type\":\"object\",\"name\":\"AllLabels\",\"id\":\"p1553\"}}}],\"below\":[{\"type\":\"object\",\"name\":\"LinearAxis\",\"id\":\"p1543\",\"attributes\":{\"ticker\":{\"type\":\"object\",\"name\":\"BasicTicker\",\"id\":\"p1545\",\"attributes\":{\"mantissas\":[1,2,5]}},\"formatter\":{\"type\":\"object\",\"name\":\"BasicTickFormatter\",\"id\":\"p1544\"},\"major_label_policy\":{\"type\":\"object\",\"name\":\"AllLabels\",\"id\":\"p1546\"}}}],\"center\":[{\"type\":\"object\",\"name\":\"Grid\",\"id\":\"p1549\",\"attributes\":{\"axis\":{\"id\":\"p1543\"}}},{\"type\":\"object\",\"name\":\"Grid\",\"id\":\"p1556\",\"attributes\":{\"dimension\":1,\"axis\":{\"id\":\"p1550\"}}}]}}],\"callbacks\":{\"type\":\"map\"}}};\n",
       "  const render_items = [{\"docid\":\"63de7880-012e-4c8e-81e5-7af7c3194325\",\"roots\":{\"p1526\":\"f8ba568e-9834-4556-a61b-fb1a50cada8d\"},\"root_ids\":[\"p1526\"]}];\n",
       "  root.Bokeh.embed.embed_items_notebook(docs_json, render_items);\n",
       "  }\n",
       "  if (root.Bokeh !== undefined) {\n",
       "    embed_document(root);\n",
       "  } else {\n",
       "    let attempts = 0;\n",
       "    const timer = setInterval(function(root) {\n",
       "      if (root.Bokeh !== undefined) {\n",
       "        clearInterval(timer);\n",
       "        embed_document(root);\n",
       "      } else {\n",
       "        attempts++;\n",
       "        if (attempts > 100) {\n",
       "          clearInterval(timer);\n",
       "          console.log(\"Bokeh: ERROR: Unable to run BokehJS code because BokehJS library is missing\");\n",
       "        }\n",
       "      }\n",
       "    }, 10, root)\n",
       "  }\n",
       "})(window);"
      ],
      "application/vnd.bokehjs_exec.v0+json": ""
     },
     "metadata": {
      "application/vnd.bokehjs_exec.v0+json": {
       "id": "p1526"
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = figure() \n",
    "p.line(x, prob_next_lime_bayes, line_width=2,color=color)\n",
    "p.circle(x,prob_next_lime_bayes, color = color)\n",
    "\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bayesian learning is very powerful but can be computationally expensive as we have to consider all possible hypotheses and this can become prohibitive in large problems. A common approximation is to make predictions based on the most probable hypothesis. This is called the maximum a posteriori or MAP hypothesis. Typically as we get more data the probability of competing hypothesis vanishes and therefore with enough data the MAP estimate tends to be the same as the Bayesian one.\n",
    "\n",
    "So in practical terms rather than taking a sum we take a maximum. Let's modify the code for plotting the prediction to follow this approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.5  0.5  0.75 1.   1.   1.   1.   1.   1.   1.   0.75]\n"
     ]
    }
   ],
   "source": [
    "prob_next_lime_map = np.zeros(11)\n",
    "lime_probs = [h1[1], h2[1], h3[1], h4[1], h5[1]]\n",
    "for n in range(0,11): \n",
    "    for (i,h) in enumerate([h1,h2,h3,h4,h5]):\n",
    "        if (n == 0):\n",
    "            posterior[i] = prior[i]\n",
    "        else: # n > 0\n",
    "            posterior[i] = (prior[i] * likelihood(data_book[:n], h))\n",
    "        \n",
    "    posterior /= np.sum(posterior)\n",
    "    # this is the previous weighted sum \n",
    "    #prob_next_lime[n] = np.dot(posterior, lime_probs)\n",
    "    \n",
    "    # instead we find the maximum posterior hypothesis and predict based on it \n",
    "    max_i = np.argmax(posterior)\n",
    "    prob_next_lime_map[n] = lime_probs[max_i]\n",
    "print(prob_next_lime_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'legend'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[32], line 2\u001b[0m\n\u001b[0;32m      1\u001b[0m p \u001b[38;5;241m=\u001b[39m figure() \n\u001b[1;32m----> 2\u001b[0m p\u001b[38;5;241m.\u001b[39mline(x, prob_next_lime_bayes, line_width\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m2\u001b[39m,color\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mblue\u001b[39m\u001b[38;5;124m'\u001b[39m, legend\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mBayes\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m      3\u001b[0m p\u001b[38;5;241m.\u001b[39mcircle(x,prob_next_lime_bayes, color \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mblue\u001b[39m\u001b[38;5;124m'\u001b[39m, legend\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mBayes\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m      4\u001b[0m p\u001b[38;5;241m.\u001b[39mline(x, prob_next_lime_map, line_width\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m2\u001b[39m,color\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mred\u001b[39m\u001b[38;5;124m'\u001b[39m, legend\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mMAP\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\bokeh\\plotting\\_decorators.py:87\u001b[0m, in \u001b[0;36mglyph_method.<locals>.decorator.<locals>.wrapped\u001b[1;34m(self, *args, **kwargs)\u001b[0m\n\u001b[0;32m     85\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcoordinates \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m     86\u001b[0m     kwargs\u001b[38;5;241m.\u001b[39msetdefault(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcoordinates\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcoordinates)\n\u001b[1;32m---> 87\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m create_renderer(glyphclass, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mplot, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\bokeh\\plotting\\_renderer.py:133\u001b[0m, in \u001b[0;36mcreate_renderer\u001b[1;34m(glyphclass, plot, **kwargs)\u001b[0m\n\u001b[0;32m    127\u001b[0m plot\u001b[38;5;241m.\u001b[39mrenderers\u001b[38;5;241m.\u001b[39mappend(glyph_renderer)\n\u001b[0;32m    129\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m legend_kwarg:\n\u001b[0;32m    130\u001b[0m     \u001b[38;5;66;03m# It must be after the renderer is added because\u001b[39;00m\n\u001b[0;32m    131\u001b[0m     \u001b[38;5;66;03m# if it creates a new `LegendItem`, the referenced\u001b[39;00m\n\u001b[0;32m    132\u001b[0m     \u001b[38;5;66;03m# renderer must already be present.\u001b[39;00m\n\u001b[1;32m--> 133\u001b[0m     update_legend(plot, legend_kwarg, glyph_renderer)\n\u001b[0;32m    135\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m glyph_renderer\n",
      "File \u001b[1;32m~\\miniconda3\\Lib\\site-packages\\bokeh\\plotting\\_legends.py:57\u001b[0m, in \u001b[0;36mupdate_legend\u001b[1;34m(plot, legend_kwarg, glyph_renderer)\u001b[0m\n\u001b[0;32m     54\u001b[0m legend \u001b[38;5;241m=\u001b[39m _get_or_create_legend(plot)\n\u001b[0;32m     55\u001b[0m kwarg, value \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlist\u001b[39m(legend_kwarg\u001b[38;5;241m.\u001b[39mitems())[\u001b[38;5;241m0\u001b[39m]\n\u001b[1;32m---> 57\u001b[0m _LEGEND_KWARG_HANDLERS[kwarg](value, legend, glyph_renderer)\n",
      "\u001b[1;31mKeyError\u001b[0m: 'legend'"
     ]
    }
   ],
   "source": [
    "p = figure() \n",
    "p.line(x, prob_next_lime_bayes, line_width=2,color='blue', legend='Bayes')\n",
    "p.circle(x,prob_next_lime_bayes, color = 'blue', legend='Bayes')\n",
    "p.line(x, prob_next_lime_map, line_width=2,color='red', legend='MAP')\n",
    "p.circle(x,prob_next_lime_map, color = 'red', legend='MAP')\n",
    "\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that in the case of uniform hypothesis prior then the MAP learning reduces to choosing the hypothesis that maximizes the likelihood. So to summarize we have three large families/types of statistical learning: \n",
    "\n",
    "1. **Bayesian Learning** is the most powerful and flexible case in which all hypotheses are considered weighted by their probabilities \n",
    "\n",
    "2. **Maximum a Posteriori (MAP)** learning is a common approximation that only considers the most probable hypothesis. It tends to be easier to compute that BL as it does not require a big summation or integration over the possible hypotheses. \n",
    "\n",
    "3. **Maximum Likelihood (ML)** learning assumes that there is a uniform prior among the hypotheses. This is the simplest approach. \n",
    "\n",
    "With enough data they all converge. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note1:** In machine learning each particular model/classifier can be considered a hypothesis. Overfitting can occur when the hypothesis space is very expressive and can capture a lot of variation in the data due to noise. \n",
    "Bayesian and MAP learning methods use the prior to penalize complexity. That way one can control the tradeoff between the complexity of a hypothesis and its degree of fit to the data. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note2:** As I have mentioned before a lot of probability calculation are basically combinations of sums and products and the Bayesian learning computations above are no exception. In MAP learning instead of doing a weighted summation (dot-product) of the probabilities for different hypothesis we select the max. A similar differentiation can be observed in Hidden Markov Models when considering the difference between maximum likelihood estimation and filtering. In a filtering operation we are interested to know the current hidden state \n",
    "given a sequence of observations. In maximum likelihood state estimation we are interested to know the entire sequence of hidden state that is most likely given a sequence of data. One question to ask is whether the result of filtering for a particular state would be always the same as the corresponding state in the maximum likelihood sequence of states. The answer is no. \n",
    "\n",
    "To understand this consider what is called a trellis diagram that shows the transitions between states in an HMM. \n",
    "To make it concrete let's say that the states are sunny and cloudy and we are interested in whether the state at step 3 is S1 (sunny). There are multiple paths through the trellis that would result in that outcome. For example \n",
    "\n",
    "1. sunny, sunny, sunny \n",
    "2. sunny, cloudy, sunny \n",
    "3. could, sunny, sunny \n",
    "4. cloudy, cloudy, sunny \n",
    "\n",
    "For filtering we would calculate the probabilities of each of these paths by taking into account both the transition and observation model and then we would sum them all up. So in this case all paths are considered. \n",
    "If you view each path as a hypotheses then this can be viewed as a type of Bayesian Learning. \n",
    "\n",
    "In contrast in maximum likelihood state estimation we compute the probability of each path as before but now at \n",
    "the end we select the path that is most probable. That gives us single path and associated probability so the summation we used in filtering becomes a max operator. This corresponds to MAP learning. \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/trellis.png\" width=\"50%\"/>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning from incomplete data  \n",
    "\n",
    "\n",
    "There are many scenarios in learning probabilistic models in which there are hidden variables. This means that we have data for the evidence variables and we are interested in doing inference for some query variables but we do **NOT** have data for the hidden variables. \n",
    "\n",
    "There are many scenarios where this is the case. For example the book describes: \n",
    "\n",
    "1. Unsupervised clustering with mixtures of Gaussians\n",
    "2. Learning Bayesian networks with hidden variables \n",
    "3. Learning HMMs from observations without associated states \n",
    "\n",
    "In all these cases, a common approach is to use the Expectation-Maximization (EM) algorithm. \n",
    "\n",
    "This is an iterative algorithm in which at each step we improve our estimation of the parameters of a probabilistic model. The basic idea is simple to describe. Our probabilistic model is characterized by several parameters that we will denote as $\\theta$. For example for a discrete Naive Bayes classifier this would be the probabilities of each feature given the class, for a continuous Naive Bayes classifier, this would be the mean and variance of each feature given the class, for a Bayesian network these would be the conditional probability tables, for a HMM these would be the transition model and sensor model probabilities. \n",
    "\n",
    "There are two steps: \n",
    "\n",
    "1. In the **expectation step (E-step)** we use the current iteration model characterized by $\\theta_i$ to compute the expected values of the hidden variables for each example. In other words we use our current probability model to infer reasonable choices for the hidden variables for each example. \n",
    "\n",
    "2. In the **maximization step (M-step)** we have now \"complete\" data in the sense that have the data values for the non hidden variables that we had from the beginning as well as expected data we obtained by using our current guess of the parameters of the model. In this step we can perform maximum-likelihood parameter estimation using complete data and obtained an updated model characterized by a new $\\theta_{i+}$. \n",
    "\n",
    "These two steps are repeated until convergence (i.e the parameters between iterations don't change - or more accurately change by a very small amount). \n",
    "\n",
    "Another way to think of it is that if we have a probabilistic model we can sample it to generate data as well as use it to fill in missing values in the data that we have. After filling in these missing values we can learn a new model and repeat the process. In several notebooks we have seen this approach where we use a model to generate some data and then use that data to estimate a model (typicaly the estimated model is close to the original model). The EM-algorithm essentially alternates between these two steps. \n",
    "\n",
    "Once you understand this basic principle, all sorts of variations and improvements can be used. For example the E-step can be done with approximate rather than exact inference or the M-step can be done with non-parametric estimation. The notations for specific problems can become relatively intimidating at first glance but if you really understand the basic principle you should be able to understand them for specific cases. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Simple binary classification example\n",
    "\n",
    "In this notebook I will show a very simple example of this idea. Hopefully this will give you some general intuition about this approach. Then you can review the specific book examples that are more complicated (learning Gaussian mixtures, Bayesian networks with hidden variables, and learning hmm parameters). We will end by showing the general mathematical notation. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's consider a simple binary classification problem with one continous attribute. For example this could be classifying whether some one is a professional basketball player or not based on their height. We can generate some synthetic data for this problem by simply sampling two Gaussian distribution. Let's say that professional basketball players have an average height of 190cm and the average height of other people is 175cm. For simiplicity we will consider they both have a standard deviation of 10cm. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[179.63715753 180.51570694 193.67540732 201.11090391 213.51312748\n",
      " 188.04117362 166.94883342 188.09248273 180.01055863 165.59677033\n",
      " 181.56870182 187.16729697 192.42486165 180.43945761 199.0209023\n",
      " 198.31499142 183.14509836 186.49230076 176.83897353 195.55258663]\n",
      "[164.63507833 175.2532028  190.69444229 172.36674335 179.55339745\n",
      " 176.36402984 172.71454488 170.32021415 173.44061735 168.50358507\n",
      " 176.63069433 167.03338012 185.03892738 187.90147743 160.53426168\n",
      " 175.33064704 163.66025875 166.75525748 152.09508458 185.51640941]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# generate twenty samples of each class \n",
    "bball_samples = np.random.normal(190, 10, 20)\n",
    "other_samples = np.random.normal(175, 10, 20)\n",
    "print(bball_samples)\n",
    "print(other_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate 1000 samples of each class and plot histogram \n",
    "\n",
    "bball_mean_height = 190 \n",
    "other_mean_height = 175 \n",
    "bball_samples = np.random.normal(bball_mean_height, 10, 1000)\n",
    "other_samples = np.random.normal(other_mean_height, 10, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAlHklEQVR4nO3de3DU1f3/8dcayJJoEhog2cSETByJighVYLiUSmAgko4IQi0KKmmRolxmgDpoRIbgfEuQKqJlwOpoClWKzhSRKRSJQxK0GMu1InUAJUKEhChCEgImgZzfH/2xdZOF3HZPspvnY+Yz4372fD7nopIXZz95r8MYYwQAAGDJdW09AAAA0LEQPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABY1amtB1BfXV2dTp06pYiICDkcjrYeDgAAaAJjjCorKxUfH6/rrrv23ka7Cx+nTp1SYmJiWw8DAAC0QHFxsRISEq7Zpt2Fj4iICEn/HXxkZGQbjwYAADRFRUWFEhMT3T/Hr6XdhY8rH7VERkYSPgAACDBNeWSCB04BAIBVhA8AAGAV4QMAAFjV7p75AAAgUBljdOnSJV2+fLmth+IXISEh6tSpU6tLYRA+AADwgZqaGpWUlOjChQttPRS/Cg8PV1xcnEJDQ1t8D8IHAACtVFdXp6KiIoWEhCg+Pl6hoaFBVyjTGKOamhp9++23KioqUq9evRotJnY1hA8AAFqppqZGdXV1SkxMVHh4eFsPx2/CwsLUuXNnHT9+XDU1NerSpUuL7sMDpwAA+EhLdwICiS/mGPyrBAAA2hXCBwAAsIpnPgAA8KOXco9Y7W/e6BSr/bUEOx8AAMAqwgcAALCK8AEAQAf17bffyuVyaenSpe5zn376qUJDQ7V9+3a/9cszH4Av5GU3PDci0/44AKAZevTooTfffFPjx49XWlqabr31Vj388MOaOXOm0tLS/NYv4QMAgA7sF7/4haZPn64pU6Zo4MCB6tKli5YtW+bXPvnYBQCADu6FF17QpUuX9O677+rtt99uceXSpiJ8AADQwR07dkynTp1SXV2djh8/7vf++NgFAIAOrKamRlOmTNGkSZN06623atq0aTp48KBiY2P91ic7HwAAdGALFy5UeXm5XnnlFS1YsEC33Xabpk2b5tc+2fkAAMCP2nPF0fz8fK1cuVJ5eXmKjIyUJP3lL39R3759tWbNGj3xxBN+6ZfwAQBAB5Wamqra2lqPcz179tS5c+f82i8fuwAAAKsIHwAAwCrCBwAAsIpnPoC2RFn2oFD/K9Pb8wOGQHvAzgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrfdgEAwJ+8/VabPwXAb8yx8wEAAKwifAAA0EGtW7dO3bp1U3V1tcf5iRMn6tFHH/Vbv4QPAAA6qAceeECXL1/W5s2b3ee+++47/f3vf9evf/1rv/VL+AAAoIMKCwvT5MmTlZOT4z739ttvKyEhQampqX7rlwdOAbQb9cuUSw1LlTelDYCmmz59ugYOHKiTJ0/qxhtvVE5OjjIyMuRwOPzWJzsfAAB0YHfeeaf69eundevWad++fTp48KAyMjL82ic7HwAAdHCPPfaYXnrpJZ08eVKjRo1SYmKiX/tj5wMAgA5uypQpOnnypF5//XX95je/8Xt/hA8AADq4yMhITZw4UTfccIPGjx/v9/742AUAAH8KgIqjklRSUqIpU6bI6XT6va9m7XysWbNGffv2VWRkpCIjIzVkyBD94x//cL9vjFFWVpbi4+MVFham1NRUHTp0yOeDBgAAvvH9999rw4YN2rFjh2bNmmWlz2aFj4SEBC1btkx79uzRnj17NHLkSI0bN84dMJYvX64VK1Zo1apV2r17t1wul0aPHq3Kykq/DB4AALTOXXfdpRkzZuj555/XLbfcYqXPZn3sMnbsWI/Xv//977VmzRoVFhaqd+/eWrlypRYuXKgJEyZIktauXavY2FitX79eM2bM8N2oAQCAT3z99dfW+2zxA6eXL1/Whg0bVFVVpSFDhqioqEilpaVKS0tzt3E6nRo+fLh27dp11ftUV1eroqLC4wAAAMGr2Q+cHjx4UEOGDNEPP/ygG264Qe+995569+7tDhixsbEe7WNjY3X8+PGr3i87O1tLlixp7jDQ0dX/impfPtDVlK+/DpAHyJrN29zbeK7eKpq2Zf9UUwVar9k7H7fccosOHDigwsJCPfHEE5o6dar+85//uN+vX47VGHPNEq2ZmZkqLy93H8XFxc0dEgAA7YIxpq2H4He+mGOzdz5CQ0N18803S5IGDBig3bt36+WXX9ZTTz0lSSotLVVcXJy7fVlZWYPdkB9zOp1Wfq0HAAB/6dy5syTpwoULCgsLa+PR+NeFCxck/W/OLdHqOh/GGFVXVys5OVkul0u5ubm68847JUk1NTUqKCjQ888/39puAABot0JCQtS1a1eVlZVJksLDw/36xWxtwRijCxcuqKysTF27dlVISEiL79Ws8PHMM88oPT1diYmJqqys1IYNG5Sfn69t27bJ4XBo7ty5Wrp0qXr16qVevXpp6dKlCg8P1+TJk1s8QAAAAoHL5ZIkdwAJVl27dnXPtaWaFT5Onz6tRx55RCUlJYqKilLfvn21bds2jR49WpK0YMECXbx4UTNnztTZs2c1aNAgbd++XREREa0aJAAA7Z3D4VBcXJxiYmJUW1vb1sPxi86dO7dqx+OKZoWPN95445rvOxwOZWVlKSsrqzVjAgAgYIWEhPjkB3Qw44vlAACAVYQPAABgFeEDAABYRfgAAABWtbrOB4COh5Lj1+atJDxrBPwPOx8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivLqgC152W09gqDlrZx5U7Sk5HlL+wLwP+x8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKuocIqOrSNVHfUy15cuTfR4PY8/EQBYwM4HAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMCqZoWP7OxsDRw4UBEREYqJidH48eN1+PBhjzYZGRlyOBwex+DBg306aAAAELiaFT4KCgo0a9YsFRYWKjc3V5cuXVJaWpqqqqo82o0ZM0YlJSXuY+vWrT4dNAAACFydmtN427ZtHq9zcnIUExOjvXv36u6773afdzqdcrlcvhkhAAAIKq165qO8vFySFB0d7XE+Pz9fMTExSklJ0fTp01VWVtaabgAAQBBp1s7HjxljNH/+fA0bNkx9+vRxn09PT9cDDzygpKQkFRUVadGiRRo5cqT27t0rp9PZ4D7V1dWqrq52v66oqGjpkAAAQABocfiYPXu2PvvsM3388cce5ydNmuT+5z59+mjAgAFKSkrSli1bNGHChAb3yc7O1pIlS1o6DCD45GU33mZEpm/uAwBtoEUfu8yZM0ebN29WXl6eEhISrtk2Li5OSUlJOnr0qNf3MzMzVV5e7j6Ki4tbMiQAABAgmrXzYYzRnDlz9N577yk/P1/JycmNXnPmzBkVFxcrLi7O6/tOp9PrxzEAACA4NWvnY9asWXrrrbe0fv16RUREqLS0VKWlpbp48aIk6fz583ryySf1ySef6Ouvv1Z+fr7Gjh2r7t276/777/fLBAAAQGBp1s7HmjVrJEmpqake53NycpSRkaGQkBAdPHhQ69at07lz5xQXF6cRI0bonXfeUUREhM8GDQAAAlezP3a5lrCwMH3wwQetGhAAAAhufLcLAACwivABAACsInwAAACrCB8AAMCqFlc4BfymfmXOplTzhP94rZQ60fowOoKXco80ODdvdEobjATwL3Y+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXl1REcvJUA91aW3Wup8ADkp3l8cuxMg3NDburWonvVLxXeHsuEeytn3pH6B9oKOx8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivLqaP9aWko8WEqp+9HgE6/Z66ypJfDRat7KtrfH8vbouNj5AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFZR4RTwFx9VWP3k2JkG54bc1M0n925p/+rZ+HX1q6d+4qXNkBEtG1NHRvVSBAN2PgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgVbPCR3Z2tgYOHKiIiAjFxMRo/PjxOnz4sEcbY4yysrIUHx+vsLAwpaam6tChQz4dNAAACFzNCh8FBQWaNWuWCgsLlZubq0uXLiktLU1VVVXuNsuXL9eKFSu0atUq7d69Wy6XS6NHj1ZlZaXPBw8AAAJPs4qMbdu2zeN1Tk6OYmJitHfvXt19990yxmjlypVauHChJkyYIElau3atYmNjtX79es2YMcN3IwcAAAGpVc98lJeXS5Kio6MlSUVFRSotLVVaWpq7jdPp1PDhw7Vr1y6v96iurlZFRYXHAQAAgleLy6sbYzR//nwNGzZMffr0kSSVlpZKkmJjYz3axsbG6vjx417vk52drSVLlrR0GGhP6pcTH5HZ/GsQEOqXTpdeaNF9vJUKBxD8WrzzMXv2bH322Wf661//2uA9h8Ph8doY0+DcFZmZmSovL3cfxcXFLR0SAAAIAC3a+ZgzZ442b96snTt3KiEhwX3e5XJJ+u8OSFxcnPt8WVlZg92QK5xOp5xOZ0uGAQAAAlCzdj6MMZo9e7Y2btyoHTt2KDk52eP95ORkuVwu5ebmus/V1NSooKBAQ4cO9c2IAQBAQGvWzsesWbO0fv16vf/++4qIiHA/4xEVFaWwsDA5HA7NnTtXS5cuVa9evdSrVy8tXbpU4eHhmjx5sl8mAAAAAkuzwseaNWskSampqR7nc3JylJGRIUlasGCBLl68qJkzZ+rs2bMaNGiQtm/froiICJ8MGAAABLZmhQ9jTKNtHA6HsrKylJWV1dIxAQCAIMZ3uwAAAKsIHwAAwCrCBwAAsIrwAQAArGpxeXWgQ6MsPCyxWYLeW1/zRqdY6x8dBzsfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCoqnAJAgLNZBRXwBXY+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXl1QG0mrfy3oPbYBztWVuXQK/f/7zRKW00EoCdDwAAYBnhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhFeXX4T152w3MjMu2PAwDQrrDzAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsanb42Llzp8aOHav4+Hg5HA5t2rTJ4/2MjAw5HA6PY/Dgwb4aLwAACHDNDh9VVVXq16+fVq1addU2Y8aMUUlJifvYunVrqwYJAACCR7PrfKSnpys9Pf2abZxOp1wuV4sHBQAAgpdfnvnIz89XTEyMUlJSNH36dJWVlV21bXV1tSoqKjwOAAAQvHxe4TQ9PV0PPPCAkpKSVFRUpEWLFmnkyJHau3evnE5ng/bZ2dlasmSJr4cBf/NWvRTWfHLsjMfrITd1a/Y1vjT4xGs+ua6w529b1Ab+81LuEY/X80antFnftvuH//g8fEyaNMn9z3369NGAAQOUlJSkLVu2aMKECQ3aZ2Zmav78+e7XFRUVSkxM9PWwAABAO+H373aJi4tTUlKSjh496vV9p9PpdUcEAAAEJ7/X+Thz5oyKi4sVFxfn764AAEAAaPbOx/nz5/Xll1+6XxcVFenAgQOKjo5WdHS0srKyNHHiRMXFxenrr7/WM888o+7du+v+++/36cABAEBganb42LNnj0aMGOF+feV5jalTp2rNmjU6ePCg1q1bp3PnzikuLk4jRozQO++8o4iICN+NGgAABKxmh4/U1FQZY676/gcffNCqAQEAgODGd7sAAACrCB8AAMAqwgcAALCK8AEAAKzye5ExwANl2QGgw2PnAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFZRXr2jqV/efERm868BgGt4KfdIg3PzRqe0wUjQXrHzAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKyiwimAgDL4xGstuq6w5299PJLA5q0KaVv3TxXUjoOdDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhFeXWgnfnk2Jl2dR/AF9q6nDvaF3Y+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVzQ4fO3fu1NixYxUfHy+Hw6FNmzZ5vG+MUVZWluLj4xUWFqbU1FQdOnTIV+MFAAABrtnho6qqSv369dOqVau8vr98+XKtWLFCq1at0u7du+VyuTR69GhVVla2erAAACDwNbvOR3p6utLT072+Z4zRypUrtXDhQk2YMEGStHbtWsXGxmr9+vWaMWNG60YLAAACnk+f+SgqKlJpaanS0tLc55xOp4YPH65du3b5sisAABCgfFrhtLS0VJIUGxvrcT42NlbHjx/3ek11dbWqq6vdrysqKnw5JAAA0M745bddHA6Hx2tjTINzV2RnZysqKsp9JCYm+mNIAACgnfBp+HC5XJL+twNyRVlZWYPdkCsyMzNVXl7uPoqLi305JAAA0M74NHwkJyfL5XIpNzfXfa6mpkYFBQUaOnSo12ucTqciIyM9DgAAELya/czH+fPn9eWXX7pfFxUV6cCBA4qOjlbPnj01d+5cLV26VL169VKvXr20dOlShYeHa/LkyT4dOAAACEzNDh979uzRiBEj3K/nz58vSZo6dar+/Oc/a8GCBbp48aJmzpyps2fPatCgQdq+fbsiIiJ8N2oAABCwmh0+UlNTZYy56vsOh0NZWVnKyspqzbgAAECQ4rtdAACAVYQPAABgFeEDAABYRfgAAABW+bS8OizKy268zYhM/48DTfbJsTMNzg25qVsbjARXM/jEaw3OFfb8bRuMJLC8lHukXd0H7R87HwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8uptzVuZdMqidxjeSq635X3amrfy5mhfKIEOX2DnAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFZRXj2YeSvd7s/rOpCmlDMfclM3CyOBbU0pAV/Y87cWRtIx1S/vPm90il+vg3+w8wEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsosKpbVQPDThNqWaK9q8plUn92RdVT4H/YecDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFb5PHxkZWXJ4XB4HC6Xy9fdAACAAOWXOh+33367PvzwQ/frkJAQf3QDAAACkF/CR6dOndjtAAAAXvnlmY+jR48qPj5eycnJevDBB3Xs2LGrtq2urlZFRYXHAQAAgpfPdz4GDRqkdevWKSUlRadPn9b//d//aejQoTp06JC6devWoH12draWLFni62EENpsl2Cn33oCvyqlTlh2NqV+GnRLszfdS7hG/3Wve6BSf3RuefL7zkZ6erokTJ+qOO+7QqFGjtGXLFknS2rVrvbbPzMxUeXm5+yguLvb1kAAAQDvi9y+Wu/7663XHHXfo6NGjXt93Op1yOp3+HgYAAGgn/F7no7q6Wl988YXi4uL83RUAAAgAPg8fTz75pAoKClRUVKRPP/1Uv/zlL1VRUaGpU6f6uisAABCAfP6xyzfffKOHHnpI3333nXr06KHBgwersLBQSUlJvu4KAAAEIJ+Hjw0bNvj6lgAAIIjw3S4AAMAqwgcAALCK8AEAAKwifAAAAKv8XmSsQ6FUOYBWqF9uXaLkOoITOx8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKiqcthTVTAPSJ8fOtPUQAAS5l3KPNNpm3ugUCyNpv9j5AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWUV5colW5R/fLmQ27q1qI2gC2DT7zWru7jq3sX9vxti+7d1OuCQaCWSa8/7vY4RnY+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYJXDGGPaehA/VlFRoaioKJWXlysyMtJOpwFYXr1+CXKpaWXIvV3XEr7qqynl1Ztyna/mBeDabJZX91Y2vqX92ywT39K+6pdB91bevSml0tuqLHxzfn6z8wEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArPJb+Fi9erWSk5PVpUsX9e/fXx999JG/ugIAAAHEL+HjnXfe0dy5c7Vw4ULt379fP//5z5Wenq4TJ074ozsAABBA/BI+VqxYoWnTpumxxx7TbbfdppUrVyoxMVFr1qzxR3cAACCAdPL1DWtqarR37149/fTTHufT0tK0a9euBu2rq6tVXV3tfl1eXi7pv8VKrKn6wV5fPlJ1sbrBuYomzMPbdS3hq7683acl1/lqXgCu7Yeq89b68vb/dUv7r38vf86jpX3V/7nn7bqm/GxsSn/++Bl75Z5Nql1qfOzkyZNGkvnnP//pcf73v/+9SUlJadB+8eLFRhIHBwcHBwdHEBzFxcWNZgWf73xc4XA4PF4bYxqck6TMzEzNnz/f/bqurk7ff/+9unXr5rV9a1RUVCgxMVHFxcX2Sre3Ix19/hJr0NHnL7EGHX3+Emvgr/kbY1RZWan4+PhG2/o8fHTv3l0hISEqLS31OF9WVqbY2NgG7Z1Op5xOp8e5rl27+npYHiIjIzvkf3BXdPT5S6xBR5+/xBp09PlLrIE/5h8VFdWkdj5/4DQ0NFT9+/dXbm6ux/nc3FwNHTrU190BAIAA45ePXebPn69HHnlEAwYM0JAhQ/Taa6/pxIkTevzxx/3RHQAACCB+CR+TJk3SmTNn9Nxzz6mkpER9+vTR1q1blZSU5I/umszpdGrx4sUNPubpKDr6/CXWoKPPX2INOvr8JdagPczfYUxTficGAADAN/huFwAAYBXhAwAAWEX4AAAAVhE+AACAVQEfPnbu3KmxY8cqPj5eDodDmzZt8ng/IyNDDofD4xg8eLBHm+rqas2ZM0fdu3fX9ddfr/vuu0/ffPONxVm0TmNrIElffPGF7rvvPkVFRSkiIkKDBw/2+JbhQF6DxuZf/9//leMPf/iDu00gz19qfA3Onz+v2bNnKyEhQWFhYbrtttsafNFjIK9BY/M/ffq0MjIyFB8fr/DwcI0ZM0ZHjx71aBPI88/OztbAgQMVERGhmJgYjR8/XocPH/ZoY4xRVlaW4uPjFRYWptTUVB06dMijTaCuQVPmv3HjRt1zzz3q3r27HA6HDhw40OA+gTp/qfE1qK2t1VNPPaU77rhD119/veLj4/Xoo4/q1KlTHvextQYBHz6qqqrUr18/rVq16qptxowZo5KSEvexdetWj/fnzp2r9957Txs2bNDHH3+s8+fP695779Xly5f9PXyfaGwNvvrqKw0bNky33nqr8vPz9e9//1uLFi1Sly5d3G0CeQ0am/+P/92XlJTozTfflMPh0MSJE91tAnn+UuNrMG/ePG3btk1vvfWWvvjiC82bN09z5szR+++/724TyGtwrfkbYzR+/HgdO3ZM77//vvbv36+kpCSNGjVKVVVV7naBPP+CggLNmjVLhYWFys3N1aVLl5SWluYxv+XLl2vFihVatWqVdu/eLZfLpdGjR6uystLdJlDXoCnzr6qq0s9+9jMtW7bsqvcJ1PlLja/BhQsXtG/fPi1atEj79u3Txo0bdeTIEd13330e97G2Bj74Lrl2Q5J57733PM5NnTrVjBs37qrXnDt3znTu3Nls2LDBfe7kyZPmuuuuM9u2bfPTSP3H2xpMmjTJPPzww1e9JpjWwNv86xs3bpwZOXKk+3Uwzd8Y72tw++23m+eee87j3F133WWeffZZY0xwrUH9+R8+fNhIMp9//rn73KVLl0x0dLR5/fXXjTHBNX9jjCkrKzOSTEFBgTHGmLq6OuNyucyyZcvcbX744QcTFRVlXn31VWNMcK1B/fn/WFFRkZFk9u/f73E+mOZvzLXX4Ip//etfRpI5fvy4McbuGgT8zkdT5OfnKyYmRikpKZo+fbrKysrc7+3du1e1tbVKS0tzn4uPj1efPn20a9euthiuT9XV1WnLli1KSUnRPffco5iYGA0aNMhjWzrY1+DHTp8+rS1btmjatGnucx1h/sOGDdPmzZt18uRJGWOUl5enI0eO6J577pEU3GtQXf3frzf/8U5fSEiIQkND9fHHH0sKvvmXl5dLkqKjoyVJRUVFKi0t9Zif0+nU8OHD3fMLpjWoP/+mCKb5S01bg/LycjkcDvf3qdlcg6APH+np6Xr77be1Y8cOvfjii9q9e7dGjhzp/gOptLRUoaGh+slPfuJxXWxsbIMvxwtEZWVlOn/+vJYtW6YxY8Zo+/btuv/++zVhwgQVFBRICv41+LG1a9cqIiJCEyZMcJ/rCPN/5ZVX1Lt3byUkJCg0NFRjxozR6tWrNWzYMEnBvQa33nqrkpKSlJmZqbNnz6qmpkbLli1TaWmpSkpKJAXX/I0xmj9/voYNG6Y+ffpIknsO9b/c88fzC5Y18Db/pgiW+UtNW4MffvhBTz/9tCZPnuz+cjmba+CX8urtyaRJk9z/3KdPHw0YMEBJSUnasmWLxw+g+owxcjgcNoboV3V1dZKkcePGad68eZKkn/70p9q1a5deffVVDR8+/KrXBssa/Nibb76pKVOmePwt+GqCaf6vvPKKCgsLtXnzZiUlJWnnzp2aOXOm4uLiNGrUqKteFwxr0LlzZ/3tb3/TtGnTFB0drZCQEI0aNUrp6emNXhuI8589e7Y+++wz967Oj9WfS1PmF2hrcK35t0SgzV9qfA1qa2v14IMPqq6uTqtXr270fv5Yg6Df+agvLi5OSUlJ7ifdXS6XampqdPbsWY92ZWVlDf6WEIi6d++uTp06qXfv3h7nb7vtNvdvuwT7Glzx0Ucf6fDhw3rsscc8zgf7/C9evKhnnnlGK1as0NixY9W3b1/Nnj1bkyZN0gsvvCAp+Negf//+OnDggM6dO6eSkhJt27ZNZ86cUXJysqTgmf+cOXO0efNm5eXlKSEhwX3e5XJJUoO/vf54fsGwBlebf1MEw/ylxtegtrZWv/rVr1RUVKTc3Fz3rodkdw06XPg4c+aMiouLFRcXJ+m/fyh17txZubm57jYlJSX6/PPPNXTo0LYaps+EhoZq4MCBDX7t7MiRI+4v+gv2NbjijTfeUP/+/dWvXz+P88E+/9raWtXW1uq66zz/dw8JCXHvjAX7GlwRFRWlHj166OjRo9qzZ4/GjRsnKfDnb4zR7NmztXHjRu3YscMdqq5ITk6Wy+XymF9NTY0KCgrc8wvkNWhs/k0RyPOXmrYGV4LH0aNH9eGHH6pbt24e71tdA58+vtoGKisrzf79+83+/fuNJLNixQqzf/9+c/z4cVNZWWl+97vfmV27dpmioiKTl5dnhgwZYm688UZTUVHhvsfjjz9uEhISzIcffmj27dtnRo4cafr162cuXbrUhjNrumutgTHGbNy40XTu3Nm89tpr5ujRo+aPf/yjCQkJMR999JH7HoG8Bo3N3xhjysvLTXh4uFmzZo3XewTy/I1pfA2GDx9ubr/9dpOXl2eOHTtmcnJyTJcuXczq1avd9wjkNWhs/u+++67Jy8szX331ldm0aZNJSkoyEyZM8LhHIM//iSeeMFFRUSY/P9+UlJS4jwsXLrjbLFu2zERFRZmNGzeagwcPmoceesjExcUFxZ+FTZn/mTNnzP79+82WLVuMJLNhwwazf/9+U1JS4m4TqPM3pvE1qK2tNffdd59JSEgwBw4c8GhTXV3tvo+tNQj48JGXl2ckNTimTp1qLly4YNLS0kyPHj1M586dTc+ePc3UqVPNiRMnPO5x8eJFM3v2bBMdHW3CwsLMvffe26BNe3atNbjijTfeMDfffLPp0qWL6devn9m0aZPHPQJ5DZoy/z/96U8mLCzMnDt3zus9Ann+xjS+BiUlJSYjI8PEx8ebLl26mFtuucW8+OKLpq6uzn2PQF6Dxub/8ssvm4SEBPefA88++6zHH7jGBPb8vc1dksnJyXG3qaurM4sXLzYul8s4nU5z9913m4MHD3rcJ1DXoCnzz8nJ8dpm8eLF7jaBOn9jGl+DK79i7O3Iy8tz38fWGjj+/6ABAACs6HDPfAAAgLZF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGDV/wPM/1k6omFkrwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot\n",
    "bins = np.linspace(150, 220, 100)\n",
    "\n",
    "pyplot.hist(bball_samples, bins, alpha=0.5, label='x')\n",
    "pyplot.hist(other_samples, bins, alpha=0.5, label='y')\n",
    "pyplot.legend(loc='upper right')\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can clearly see in the histogram the height-distribution and overlap. You can also see that there is an equal number of instances for each class and that the standard deviation is the same.\n",
    "\n",
    "Now suppose that you are just given the nba_samples and other_samples and told that these are labeled samples for training a Naive Bayes classifier. You also know that they both have a standard deviation of 10cm so we will keep that. In this case the only parameter we are trying to estimate is the mean of each class. So 𝜃=(𝜇𝑛𝑏𝑎,𝜇𝑜𝑡ℎ𝑒𝑟).\n",
    "\n",
    "Given this data the maximum-likelihood estimate for the means is easily obtained by taking the statistical mean of the samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "190.53401424928552 174.95131543223874\n"
     ]
    }
   ],
   "source": [
    "mu_bball = np.mean(bball_samples)\n",
    "mu_other = np.mean(other_samples)\n",
    "print(mu_bball, mu_other)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have \"learned\" a model we can use it to predict. Suppose you are given a test height - lets say 183cm. You can calcuate the $P(183/nba)$ and $P(183/other)$ by using the corresponding probability density functions characterized by $\\mu_{nba}$ and $\\sigma = 10$ and $\\mu_{other}$ and $\\sigma = 10$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.025323796182471114 0.033224871192996344\n",
      "181 is more likely NOT a professional basketball player\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "test_height = 181\n",
    "\n",
    "p_bball = norm(mu_bball, 10).pdf(test_height)\n",
    "p_other = norm(mu_other, 10).pdf(test_height)\n",
    "print(p_bball, p_other)\n",
    "\n",
    "if (p_bball > p_other): \n",
    "    print(str(test_height) + \" is more likely a professional basketball player\")\n",
    "else: \n",
    "    print(str(test_height) + \" is more likely NOT a professional basketball player\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unsupervised learning \n",
    "\n",
    "Now let's make a very simple change to the problem above. Suppose that we know that we have two classes and we have heights. We also know that the standard deviation for each class is 10cm. However we are not given the labels of the heights but just a dataset of heights. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[183.25779348 166.95106925 199.62750884 ... 195.05905986 193.4413794\n",
      " 175.45981639]\n"
     ]
    }
   ],
   "source": [
    "heights = np.hstack([bball_samples,other_samples])\n",
    "np.random.shuffle(heights)\n",
    "print(heights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we know the means of the two classes, we can predict the class of each instance (height)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bball 183.25779347707206\n",
      "other 166.95106924922794\n",
      "bball 199.6275088362795\n",
      "bball 189.4955703895145\n",
      "bball 187.08277775934175\n",
      "other 163.49455706588833\n",
      "bball 205.42239628975142\n",
      "other 181.86652521723863\n",
      "other 176.55549881153638\n",
      "other 177.02698911845664\n",
      "other 169.66429240948463\n",
      "other 177.63082459926736\n",
      "other 169.3762095288119\n",
      "other 163.72589689590777\n",
      "bball 187.81539181016674\n",
      "other 170.25697717252567\n",
      "bball 203.61129219160765\n",
      "other 179.8861067241242\n",
      "bball 184.332336695728\n",
      "other 174.50195421670395\n"
     ]
    }
   ],
   "source": [
    "estimated_bball_samples = [] \n",
    "estimated_other_samples = []\n",
    "\n",
    "for h in heights[:20]: \n",
    "    p_bball = norm(mu_bball, 10).pdf(h)\n",
    "    p_other = norm(mu_other, 10).pdf(h)\n",
    "    if (p_bball > p_other): \n",
    "        print('bball', h) \n",
    "    else: \n",
    "        print('other', h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the idea of the EM-algorithm would be let's start with a reasonable guess of the two means, then predict the heights, the re-estimate the means and keep repeating until convergence. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 170, 165)\n",
      "(1, 185.578963340917, 161.94314250629293)\n",
      "(2, 188.21621662003975, 166.62429638620307)\n",
      "(3, 189.97027789969755, 169.05099263505087)\n",
      "(4, 191.00096489821175, 170.30347101989184)\n",
      "(5, 191.65572670540232, 171.07128821893536)\n",
      "(6, 192.11653003800765, 171.58232369135806)\n",
      "(7, 192.5188718023523, 172.00217303721027)\n",
      "(8, 192.79631610824157, 172.27866046032432)\n",
      "(9, 192.9728514261467, 172.45091243841435)\n",
      "(10, 193.076897915669, 172.55210531406644)\n",
      "(11, 193.13957653616455, 172.61260379104624)\n",
      "(12, 193.2025456312121, 172.6728482996125)\n",
      "(13, 193.22363369620817, 172.69284259837065)\n",
      "(14, 193.24466663248413, 172.71289382364347)\n",
      "(15, 193.2551875178223, 172.72291666418914)\n",
      "(16, 193.2551875178223, 172.72291666418914)\n",
      "(17, 193.2551875178223, 172.72291666418914)\n",
      "(18, 193.2551875178223, 172.72291666418914)\n",
      "(19, 193.2551875178223, 172.72291666418914)\n"
     ]
    }
   ],
   "source": [
    "mu_bball = 170 \n",
    "mu_other = 165 \n",
    "\n",
    "\n",
    "for i in range(0,20): \n",
    "    estimated_bball_samples = [] \n",
    "    estimated_other_samples = [] \n",
    "    print((i,mu_bball, mu_other))\n",
    "\n",
    "    for h in heights:  \n",
    "        # E-step - use current model to estimate values for the hidden variable (class membership)\n",
    "        p_bball = norm(mu_bball, 10).pdf(h)\n",
    "        p_other = norm(mu_other, 10).pdf(h)\n",
    "        if (p_bball > p_other): \n",
    "            estimated_bball_samples.append(h)\n",
    "        else: \n",
    "            estimated_other_samples.append(h)\n",
    "   \n",
    "    # M-step - using the estimated class values re-caculate the parameters of the model i.e the means \n",
    "    mu_bball = np.mean(estimated_bball_samples)\n",
    "    mu_other = np.mean(estimated_other_samples)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a simple example of the EM-algorithm with only one parameter to estimate per class to show the basic principle. In more real-world scenarios the probabilistic model can be much more complex with many parameters. \n",
    "In addition the prediction step typically is done using likelihoods and weighted samples rather than simple \n",
    "prediction as in this example. Therefore the E-step and M-step tend to be more complicated and have complex update equations but the basic principle is the same. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Advanced**: A great tutorial article on EM is: \n",
    "\n",
    "Bilmes JA. A gentle tutorial of the EM algorithm and its application to parameter estimation for Gaussian mixture and hidden Markov models. International Computer Science Institute. 1998 Apr 21;4(510):126.\n",
    "\n",
    "\n",
    "http://www.leap.ee.iisc.ac.in/sriram/teaching/MLSP_18/refs/GMM_Bilmes.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Author: Ron Weiss <ronweiss@gmail.com>, Gael Varoquaux\n",
    "# Modified by Thierry Guillemot <thierry.guillemot.work@gmail.com>\n",
    "# License: BSD 3 clause\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from sklearn import datasets\n",
    "from sklearn.mixture import GaussianMixture\n",
    "from sklearn.model_selection import StratifiedKFold\n",
    "\n",
    "colors = [\"navy\", \"turquoise\", \"darkorange\"]\n",
    "\n",
    "\n",
    "def make_ellipses(gmm, ax):\n",
    "    for n, color in enumerate(colors):\n",
    "        if gmm.covariance_type == \"full\":\n",
    "            covariances = gmm.covariances_[n][:2, :2]\n",
    "        elif gmm.covariance_type == \"tied\":\n",
    "            covariances = gmm.covariances_[:2, :2]\n",
    "        elif gmm.covariance_type == \"diag\":\n",
    "            covariances = np.diag(gmm.covariances_[n][:2])\n",
    "        elif gmm.covariance_type == \"spherical\":\n",
    "            covariances = np.eye(gmm.means_.shape[1]) * gmm.covariances_[n]\n",
    "        v, w = np.linalg.eigh(covariances)\n",
    "        u = w[0] / np.linalg.norm(w[0])\n",
    "        angle = np.arctan2(u[1], u[0])\n",
    "        angle = 180 * angle / np.pi  # convert to degrees\n",
    "        v = 2.0 * np.sqrt(2.0) * np.sqrt(v)\n",
    "        ell = mpl.patches.Ellipse(\n",
    "            gmm.means_[n, :2], v[0], v[1], angle=180 + angle, color=color\n",
    "        )\n",
    "        ell.set_clip_box(ax.bbox)\n",
    "        ell.set_alpha(0.5)\n",
    "        ax.add_artist(ell)\n",
    "        ax.set_aspect(\"equal\", \"datalim\")\n",
    "\n",
    "\n",
    "iris = datasets.load_iris()\n",
    "\n",
    "# Break up the dataset into non-overlapping training (75%) and testing\n",
    "# (25%) sets.\n",
    "skf = StratifiedKFold(n_splits=4)\n",
    "# Only take the first fold.\n",
    "train_index, test_index = next(iter(skf.split(iris.data, iris.target)))\n",
    "\n",
    "\n",
    "X_train = iris.data[train_index]\n",
    "y_train = iris.target[train_index]\n",
    "X_test = iris.data[test_index]\n",
    "y_test = iris.target[test_index]\n",
    "\n",
    "n_classes = len(np.unique(y_train))\n",
    "\n",
    "# Try GMMs using different types of covariances.\n",
    "estimators = {\n",
    "    cov_type: GaussianMixture(\n",
    "        n_components=n_classes, covariance_type=cov_type, max_iter=20, random_state=0\n",
    "    )\n",
    "    for cov_type in [\"spherical\", \"diag\", \"tied\", \"full\"]\n",
    "}\n",
    "\n",
    "n_estimators = len(estimators)\n",
    "\n",
    "plt.figure(figsize=(3 * n_estimators // 2, 6))\n",
    "plt.subplots_adjust(\n",
    "    bottom=0.01, top=0.95, hspace=0.15, wspace=0.05, left=0.01, right=0.99\n",
    ")\n",
    "\n",
    "\n",
    "for index, (name, estimator) in enumerate(estimators.items()):\n",
    "    # Since we have class labels for the training data, we can\n",
    "    # initialize the GMM parameters in a supervised manner.\n",
    "    estimator.means_init = np.array(\n",
    "        [X_train[y_train == i].mean(axis=0) for i in range(n_classes)]\n",
    "    )\n",
    "\n",
    "    # Train the other parameters using the EM algorithm.\n",
    "    estimator.fit(X_train)\n",
    "\n",
    "    h = plt.subplot(2, n_estimators // 2, index + 1)\n",
    "    make_ellipses(estimator, h)\n",
    "\n",
    "    for n, color in enumerate(colors):\n",
    "        data = iris.data[iris.target == n]\n",
    "        plt.scatter(\n",
    "            data[:, 0], data[:, 1], s=0.8, color=color, label=iris.target_names[n]\n",
    "        )\n",
    "    # Plot the test data with crosses\n",
    "    for n, color in enumerate(colors):\n",
    "        data = X_test[y_test == n]\n",
    "        plt.scatter(data[:, 0], data[:, 1], marker=\"x\", color=color)\n",
    "\n",
    "    y_train_pred = estimator.predict(X_train)\n",
    "    train_accuracy = np.mean(y_train_pred.ravel() == y_train.ravel()) * 100\n",
    "    plt.text(0.05, 0.9, \"Train accuracy: %.1f\" % train_accuracy, transform=h.transAxes)\n",
    "\n",
    "    y_test_pred = estimator.predict(X_test)\n",
    "    test_accuracy = np.mean(y_test_pred.ravel() == y_test.ravel()) * 100\n",
    "    plt.text(0.05, 0.8, \"Test accuracy: %.1f\" % test_accuracy, transform=h.transAxes)\n",
    "\n",
    "    plt.xticks(())\n",
    "    plt.yticks(())\n",
    "    plt.title(name)\n",
    "\n",
    "plt.legend(scatterpoints=1, loc=\"lower right\", prop=dict(size=12))\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets\n",
    "digits = datasets.load_digits()\n",
    "# We have 1797 instances of digits stored as 8x8=64 pixel images. These images are directly used as features. \n",
    "print(digits.data.shape)\n",
    "# \n",
    "print(digits.data) \n",
    "# the 10th instance \n",
    "print(digits.data[10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can use classifier to train a model for predicting the digit \n",
    "# of an 8 by 8 image \n",
    "\n",
    "from sklearn import svm\n",
    "clf = svm.SVC(gamma=0.001, C=100.)\n",
    "\n",
    "# as a simple example we train using all the data except the last image \n",
    "clf.fit(digits.data[:-1], digits.target[:-1]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can predict the last instance and the result is the number 8 \n",
    "print(clf.predict(digits.data[-1:]))\n",
    "# predicit instance number gives as the number 1 as the prediction \n",
    "print(clf.predict([digits.data[1]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%pylab inline\n",
    "\n",
    "plt.imshow(digits.images[-1:][0], cmap=plt.cm.gray_r, interpolation='nearest')\n",
    "plt.figure()\n",
    "plt.imshow(digits.images[1:][0], cmap=plt.cm.gray_r, interpolation='nearest')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import librosa\n",
    "import numpy as np\n",
    "\n",
    "fnames = glob.glob(\"/home/gtzan/data/sound/genres44k/*/*.wav\")\n",
    "\n",
    "genres = ['classical', 'country', 'disco', 'hiphop', 'jazz', 'rock', 'blues', 'reggae', 'pop', 'metal']\n",
    "\n",
    "# allocate matrix for audio features and target \n",
    "audio_features = np.zeros((len(fnames), 40))\n",
    "target = np.zeros(len(fnames))\n",
    "\n",
    "# compute the features \n",
    "for (i,fname) in enumerate(fnames): \n",
    "    print(\"Processing %d %s\" % (i, fname))\n",
    "    for (label,genre) in enumerate(genres): \n",
    "        if genre in fname: \n",
    "            audio, srate = librosa.load(fname)\n",
    "            mfcc_matrix = librosa.feature.mfcc(y=audio, sr=srate)\n",
    "            mean_mfcc = np.mean(mfcc_matrix,axis=1)\n",
    "            std_mfcc = np.std(mfcc_matrix, axis=1)\n",
    "            audio_fvec = np.hstack([mean_mfcc, std_mfcc])\n",
    "            audio_features[i] = audio_fvec\n",
    "            target[i] = label\n",
    "\n",
    "print(audio_features.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import svm, metrics\n",
    "from sklearn.model_selection import cross_val_predict\n",
    "\n",
    "clf = svm.SVC(gamma='scale', kernel='linear')\n",
    "clf.fit(audio_features, target)\n",
    "predicted = cross_val_predict(clf, audio_features, target, cv=10)\n",
    "\n",
    "print(\"Confusion matrix:\\n%s\" % metrics.confusion_matrix(target, predicted))\n",
    "print(\"Classification report for classifier %s:\\n%s\\n\"\n",
    "      % (clf, metrics.accuracy_score(target, predicted)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform 10-fold cross-validation to calculate accuracy and confusion matrix \n",
    "predicted = cross_val_predict(clf, audio_features, target, cv=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Confusion matrix:\\n%s\" % metrics.confusion_matrix(target, predicted))\n",
    "print(\"\\nClassification accuracy:%s\\n\"\n",
    "      % (metrics.accuracy_score(target, predicted)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random \n",
    "i = random.randint(0, 1000)\n",
    "predicted_label = int(clf.predict([audio_features[i]]))\n",
    "print(\"Prediction for %d %s is %s\"% (i, fnames[i], genres[predicted_label]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|Before ending this notebook I would like to mention how impressive it is that with a few lines of Python code and the right libraries we were able to do a full music genre classification system. Twenty years ago building such a system required many lines of C++ code and formed a large part of my PhD thesis which was completed in 2002. The journal article describing this work called \"Musical Genre Classification of Audio Signal\" is considered a classic paper in MIR."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}