01_mle_map.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#default_exp mle_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#hide\n",
    "from nbdev.showdoc import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to MLE and MAP\n",
    "> Summary: Maximum Likelihood Estimation, Maximum a Posteriori Estimation, MLE, MAP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum Likelihood Estimation (MLE)\n",
    "\n",
    "Maximum Likelihood Estimation (MLE) is a principle that estimates the parameters of a statistical model, which makes the observed data most probable. In other words, MLE maximizes the data likelihood.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameter Estimation: Estimating the Probability of Heads\n",
    "\n",
    "Let's assume we have a random variable $X$ representing a coin. We can estimate the probability that it will turn up heads ($X = 1$) or tails ($X = 0$).\n",
    "\n",
    "> Task: Estimate the probability of heads $\\theta = P(X = 1)$\n",
    "\n",
    "\n",
    "Evidently, if $P(X=1)=\\theta$, then $P(X=0)=1-\\theta$. Since we do not know the \"true\" probability of heads, i.e. $P(X=1) = \\theta$, we will use $\\hat\\theta$ to refer to its estimate.\n",
    "\n",
    "**Question:** _What is the probability of $\\theta = P(X=1)?$_\n",
    "\n",
    "In general, _Maximum Likelihood Estimation_ principle asks to choose parameter $\\theta$ that maximizes $P(Data|\\theta)$, or in other words maximizes the probability of the observed data. We assume that $\\theta$ belongs to the set $\\Theta \\subset \\mathbb{R}^n$. Therefore,\n",
    "\n",
    "$$\\hat\\theta_{MLE} = \\underset{\\theta}{\\arg\\max}\\ P(Data|\\theta)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In regards to our coin flip example, if we flip the coin repeatedly, we observe that:\n",
    "\n",
    "- It turns up heads $\\alpha_1$ times\n",
    "- It turns up tails $\\alpha_0$ times\n",
    "\n",
    "Intuitively, we can estimate the $P(X=1)$ from our training data (number of tosses) as the fraction of flips that ends up heads:\n",
    "\n",
    "$$ P(X=1) = \\frac{\\alpha_1}{\\alpha_1 + \\alpha_0}$$\n",
    "\n",
    "For instance, if we flip the coin 40 times, seeing 18 heads and 22 tails, then we can estimate that:\n",
    "\n",
    "$$\\hat\\theta = P(X=1) = \\frac{18}{40} = 0.45$$\n",
    "\n",
    "And if we flip it 5 times, observing 3 heads and 2 tails, then we have:\n",
    "\n",
    "$$\\hat\\theta = P(X=1) = \\frac{3}{5} = 0.6$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to Calculate MLE?\n",
    "\n",
    "First step in calculating the maximum likelihood estimator $\\hat\\theta$ is to define $P(Data|\\theta)$. If we flip the coin once, then $P(Data|\\theta) = \\theta$ if the flip results in heads and $P(Data|\\theta) = 1 - \\theta$, if the flips turns tails. If we observe $D = \\{1,0,1,1,0\\}$ by tossing the coin 5 times, assuming the flips are independent and identically distributed (i.i.d), then we have:\n",
    "\n",
    "$$P(Data|\\theta) = \\theta\\cdot(1-\\theta)\\cdot\\theta\\cdot\\theta\\cdot(1-\\theta) = \\theta^3\\cdot(1-\\theta)^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, if we flip the coin $n$ times, observing $\\alpha_H$ heads and $\\alpha_T$ tails, then\n",
    "\n",
    "$$P(Data|\\theta) = \\theta^{\\alpha_H}\\cdot(1-\\theta)^{\\alpha_T}$$\n",
    "\n",
    "The next step is to find the value of $\\theta$ that maximizes the $P(Data|\\theta)$. When finding the MLE, it is often easier to maximize the log-likelihood function since,\n",
    "\n",
    "$$\\underset{\\theta}{\\arg\\max} \\log P(Data|\\theta) = \\underset{\\theta}{\\arg\\max}\\ P(Data|\\theta)$$\n",
    "\n",
    "Let's call $J(\\theta) = \\log P(Data|\\theta)$. Thus, in order to find the value of the $\\theta$ that maximizes the $J(\\theta)$, we calculate the derivative of $J(\\theta)$ with respect to $\\theta$, set it to zero and solve for $\\theta$.\n",
    "\n",
    "$$\\frac{\\partial J(\\theta)}{\\partial \\theta} = \\frac{\\partial[\\alpha_H \\log \\theta + \\alpha_T \\log (1-\\theta)]}{\\partial \\theta}= \\alpha_H \\frac{1}{\\theta} - \\alpha_T \\frac{1}{1-\\theta} = 0$$\n",
    "\n",
    "Solving this for $\\theta$ gives, $$\\hat{\\theta} = \\dfrac{\\alpha_H}{\\alpha_H + \\alpha_T}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** How good is this MLE estimation?\n",
    "\n",
    "If I fliped the coin 5 times: 3 heads, 2 tails: $\\hat{\\theta}_{MLE}=\\frac{3}{5}=0.6$. What if I flipped 30 heads and 35 tails?  $\\hat{\\theta}_{MLE}=\\frac{30}{65}=0.46$\n",
    "\n",
    "Which estimator should we trust more? Let's assume that the coin is a goverment minted coin, meaning it's a fair coin and is \"close\" to 50-50 chance of heads/tails. It tells us that $\\theta$ should be more likely about $0.5$. Therefore, we cannot quite rely on the estimate $\\hat{\\theta} = 0.6$.\n",
    "\n",
    "In general, if we have plenty of data, MLE works well, but if we have a few observations such as 5 coin flips, our estimates will be unreliable. This leads us to the second principle for estimating parameters. This principle allows us to integrate our prior assumptions along with observed data to develop our ultimate estimate. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Map a Posteriori Estimation (MAP)\n",
    "\n",
    "Considering our coin flip example, we assume that the coin is a goverment minted coin, meaning the $\\theta$ is close to $0.5$. What can do we now that we have prior knowledge? How can we estimate the probability of heads?\n",
    "\n",
    "We have prior knowledge about the coin, but particularly, at the beginning we do not have enough coin flips to estimate the probability of heads ($\\theta$). However, we can add a number of _imaginary_ coin flips. For instance, if we add 10 imaginary heads and 10 imaginary tails before we even start the first flip, $\\hat{\\theta}=\\frac{10}{20}=0.5$, which describes our prior knowledge about the coin. After the first flip, if it turns up head $\\hat{\\theta}=\\frac{1+10}{1+10+10}=0.52$. \n",
    "\n",
    "We can see that as the number of coin flips increases, our final estimate becomes better, but more importantly, when we don't have plenty of flips, our estimate is still reliable. The more confident we are about our prior assumptions, the higher number of imaginary flips we can consider. Thus, we have:\n",
    "\n",
    "$$\\hat{\\theta} = \\dfrac{\\alpha_H + \\lambda_H}{(\\alpha_H+\\lambda_H) + (\\alpha_T+\\lambda_T)}$$\n",
    "\n",
    "where $\\lambda_H$ and $\\lambda_T$ are imaginary (or virtual) heads and tails respectively.\n",
    "\n",
    "### Bayesian Approach\n",
    "\n",
    "We choose a Bayesian approach, therefore, rather than estimating a single $\\theta$, we obtain a distribution over possible values of $\\theta$. Then, choose the value of $\\theta$ that is most probable, given the observed data and prior belief.\n",
    "\n",
    "We need Bayes rule to proceed.\n",
    "\n",
    "Chain rule: $$P(X,Y)=P(X|Y)P(Y)=P(Y|X)P(X)$$\n",
    "\n",
    "Bayes rule: $$P(X|Y)=\\frac{P(Y|X)P(X)}{P(Y)}$$\n",
    "\n",
    "Using the Bayes rule, we have: $$P(\\theta|Data)=\\frac{P(Data|\\theta)P(\\theta)}{P(Data)}$$\n",
    "\n",
    "Or equivalently, $$P(\\theta|Data)\\propto P(Data|\\theta)P(\\theta)$$\n",
    "\n",
    "We can get rid of $P(Data)$, because it's independent of the parameter $\\theta$.\n",
    "\n",
    "$P(\\theta|Data)$ is called \"posterior\", $P(Data|\\theta)$ is called \"likelihood\" and $P(\\theta)$ is called \"prior\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MLE vs. MAP\n",
    "\n",
    "1- Maximum Likelihood estimation (MLE):\n",
    "  \n",
    "* Choose value of $\\theta$ that maximizes the probability of observed\n",
    "data.\n",
    "$$\\hat{\\theta}_{MLE}=\\underset{\\theta}{\\arg\\max}\\ P(Data|\\theta)$$\n",
    "\n",
    "2- Maximum a posteriori (MAP) estimation:\n",
    "\n",
    "* Choose value of $\\theta$ that is most probable given observed data and\n",
    "prior belief.\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\hat{\\theta}_{MAP}&=\\underset{\\theta}{\\arg\\max}\\ P(\\theta|Data)\\\\\n",
    "&=\\underset{\\theta}{\\arg\\max}\\ P(Data|\\theta) P(\\theta)\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Map Estimation for Binomial Distribution\n",
    "\n",
    "Likelihood is Binomial:$P(Data|\\theta)={n\\choose \\alpha_H}\\theta^{\\alpha_H}(1-\\theta)^{\\alpha_T}$\n",
    "\n",
    "If we assume prior is Beta distribution: $P(\\theta)=\\frac{\\theta^{\\beta_H-1}(1-\\theta)^\\beta_T-1}{B(\\beta_H,\\beta_T)}\\sim Beta(\\beta_H, \\beta_T)$\n",
    "\n",
    " $B(x,y)=\\int_o^1 t^{x-1}(1-t)^{y-1}dt$\n",
    "\n",
    "Then, posterior is Beta distribution: $P(\\theta|Data)\\sim Beta(\\beta_H+\\alpha_H, \\beta_T+\\alpha_T)$\n",
    "\n",
    "\n",
    "And,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\hat{\\theta}_{MAP}&=\\underset{\\theta}{\\arg\\max}\\ P(Data|\\theta) P(\\theta)\\\\\n",
    "&=\\frac{\\alpha_H+\\beta_H -1}{\\alpha_H+\\beta_H+\\alpha_T+\\beta_T -2}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "- **Conjugate prior:** $P(\\theta)$ is the conjugate prior for likelihood function $P(\\theta|Data)$ if $P(\\theta)$ and $P(\\theta|Data)$ have the same form.\n",
    "\n",
    "- Beta prior is equivalent to extra coin flips\n",
    "\n",
    "- As the number of samples (e.g. coin flips) increases, the effect of prior is \"washed out\". It means as $N\\rightarrow \\infty$, prior is \"forgotten\".\n",
    "\n",
    "- For small sample size, prior is important."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implemented Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "import altair as alt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import norm\n",
    "\n",
    "def flip_coin(num_of_experiments = 1000, num_of_flips = 30):\n",
    "    \"\"\"\n",
    "    Flip the coin `num_of_flips` times and repeat this experiment `num_of_experiments` times. And\n",
    "    return the number of heads grouped together in all the experiments.\n",
    "    \"\"\"\n",
    "    all_heads = []\n",
    "    for i in range(num_of_experiments):\n",
    "        heads = tails = 0\n",
    "        for counter in range(num_of_flips):\n",
    "            num = np.random.randint(0,2)\n",
    "            if num == 0:\n",
    "                heads += 1\n",
    "            else:\n",
    "                tails += 1\n",
    "        all_heads.append(heads)\n",
    "\n",
    "    # group the number of heads in all the experiments\n",
    "    flip_heads = []\n",
    "    for flip in range(num_of_flips + 1):\n",
    "        num_of_heads = 0\n",
    "        for h in all_heads:\n",
    "            if h == flip:\n",
    "                num_of_heads += 1\n",
    "        flip_heads.append(num_of_heads)\n",
    "    return flip_heads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize the variables\n",
    "num_of_flips = 40\n",
    "num_of_experiments = 3000\n",
    "head_counts = flip_coin(num_of_experiments,num_of_flips)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot a chart and see the the number of heads in 30 flips in all the experiments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<div id=\"altair-viz-10142b9ecd3f4c34bde38391ed6390dd\"></div>\n",
       "<script type=\"text/javascript\">\n",
       "  (function(spec, embedOpt){\n",
       "    const outputDiv = document.getElementById(\"altair-viz-10142b9ecd3f4c34bde38391ed6390dd\");\n",
       "    const paths = {\n",
       "      \"vega\": \"https://cdn.jsdelivr.net/npm//vega@5?noext\",\n",
       "      \"vega-lib\": \"https://cdn.jsdelivr.net/npm//vega-lib?noext\",\n",
       "      \"vega-lite\": \"https://cdn.jsdelivr.net/npm//vega-lite@4.0.2?noext\",\n",
       "      \"vega-embed\": \"https://cdn.jsdelivr.net/npm//vega-embed@6?noext\",\n",
       "    };\n",
       "\n",
       "    function loadScript(lib) {\n",
       "      return new Promise(function(resolve, reject) {\n",
       "        var s = document.createElement('script');\n",
       "        s.src = paths[lib];\n",
       "        s.async = true;\n",
       "        s.onload = () => resolve(paths[lib]);\n",
       "        s.onerror = () => reject(`Error loading script: ${paths[lib]}`);\n",
       "        document.getElementsByTagName(\"head\")[0].appendChild(s);\n",
       "      });\n",
       "    }\n",
       "\n",
       "    function showError(err) {\n",
       "      outputDiv.innerHTML = `<div class=\"error\" style=\"color:red;\">${err}</div>`;\n",
       "      throw err;\n",
       "    }\n",
       "\n",
       "    function displayChart(vegaEmbed) {\n",
       "      vegaEmbed(outputDiv, spec, embedOpt)\n",
       "        .catch(err => showError(`Javascript Error: ${err.message}<br>This usually means there's a typo in your chart specification. See the javascript console for the full traceback.`));\n",
       "    }\n",
       "\n",
       "    if(typeof define === \"function\" && define.amd) {\n",
       "      requirejs.config({paths});\n",
       "      require([\"vega-embed\"], displayChart, err => showError(`Error loading script: ${err.message}`));\n",
       "    } else if (typeof vegaEmbed === \"function\") {\n",
       "      displayChart(vegaEmbed);\n",
       "    } else {\n",
       "      loadScript(\"vega\")\n",
       "        .then(() => loadScript(\"vega-lite\"))\n",
       "        .then(() => loadScript(\"vega-embed\"))\n",
       "        .catch(showError)\n",
       "        .then(() => displayChart(vegaEmbed));\n",
       "    }\n",
       "  })({\"config\": {\"view\": {\"continuousWidth\": 400, \"continuousHeight\": 300}}, \"data\": {\"name\": \"data-c7ad18303904753003e679bfb19caa8f\"}, \"mark\": \"bar\", \"encoding\": {\"x\": {\"type\": \"quantitative\", \"field\": \"Number of Heads\"}, \"y\": {\"type\": \"quantitative\", \"field\": \"Number of Ways\"}}, \"title\": \"Distribution of Heads\", \"width\": 360, \"$schema\": \"https://vega.github.io/schema/vega-lite/v4.0.2.json\", \"datasets\": {\"data-c7ad18303904753003e679bfb19caa8f\": [{\"Number of Heads\": 0, \"Number of Ways\": 0}, {\"Number of Heads\": 1, \"Number of Ways\": 0}, {\"Number of Heads\": 2, \"Number of Ways\": 0}, {\"Number of Heads\": 3, \"Number of Ways\": 0}, {\"Number of Heads\": 4, \"Number of Ways\": 0}, {\"Number of Heads\": 5, \"Number of Ways\": 0}, {\"Number of Heads\": 6, \"Number of Ways\": 0}, {\"Number of Heads\": 7, \"Number of Ways\": 0}, {\"Number of Heads\": 8, \"Number of Ways\": 0}, {\"Number of Heads\": 9, \"Number of Ways\": 0}, {\"Number of Heads\": 10, \"Number of Ways\": 2}, {\"Number of Heads\": 11, \"Number of Ways\": 9}, {\"Number of Heads\": 12, \"Number of Ways\": 16}, {\"Number of Heads\": 13, \"Number of Ways\": 28}, {\"Number of Heads\": 14, \"Number of Ways\": 54}, {\"Number of Heads\": 15, \"Number of Ways\": 124}, {\"Number of Heads\": 16, \"Number of Ways\": 174}, {\"Number of Heads\": 17, \"Number of Ways\": 257}, {\"Number of Heads\": 18, \"Number of Ways\": 326}, {\"Number of Heads\": 19, \"Number of Ways\": 379}, {\"Number of Heads\": 20, \"Number of Ways\": 351}, {\"Number of Heads\": 21, \"Number of Ways\": 351}, {\"Number of Heads\": 22, \"Number of Ways\": 289}, {\"Number of Heads\": 23, \"Number of Ways\": 249}, {\"Number of Heads\": 24, \"Number of Ways\": 156}, {\"Number of Heads\": 25, \"Number of Ways\": 109}, {\"Number of Heads\": 26, \"Number of Ways\": 68}, {\"Number of Heads\": 27, \"Number of Ways\": 36}, {\"Number of Heads\": 28, \"Number of Ways\": 15}, {\"Number of Heads\": 29, \"Number of Ways\": 5}, {\"Number of Heads\": 30, \"Number of Ways\": 0}, {\"Number of Heads\": 31, \"Number of Ways\": 2}, {\"Number of Heads\": 32, \"Number of Ways\": 0}, {\"Number of Heads\": 33, \"Number of Ways\": 0}, {\"Number of Heads\": 34, \"Number of Ways\": 0}, {\"Number of Heads\": 35, \"Number of Ways\": 0}, {\"Number of Heads\": 36, \"Number of Ways\": 0}, {\"Number of Heads\": 37, \"Number of Ways\": 0}, {\"Number of Heads\": 38, \"Number of Ways\": 0}, {\"Number of Heads\": 39, \"Number of Ways\": 0}, {\"Number of Heads\": 40, \"Number of Ways\": 0}]}}, {\"mode\": \"vega-lite\"});\n",
       "</script>"
      ],
      "text/plain": [
       "alt.Chart(...)"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = range(num_of_flips + 1)\n",
    "source = pd.DataFrame({\n",
    "  'Number of Heads': x,\n",
    "  'Number of Ways': head_counts\n",
    "})\n",
    "#     Bar chart\n",
    "bar_chart = alt.Chart(source).mark_bar().encode(\n",
    "    x='Number of Heads',\n",
    "    y='Number of Ways',\n",
    ").properties(title='Distribution of Heads',width=360)\n",
    "bar_chart\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does the plot look like? Does the bell shape ring a bell?\n",
    "\n",
    "Now, we plot a Normal distribution with mean = `number_of_flips` and standard deviation = $\\sqrt{mean/2}$  and check how it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<div id=\"altair-viz-10772a3e796f4a2bab5b9f8799ac775f\"></div>\n",
       "<script type=\"text/javascript\">\n",
       "  (function(spec, embedOpt){\n",
       "    const outputDiv = document.getElementById(\"altair-viz-10772a3e796f4a2bab5b9f8799ac775f\");\n",
       "    const paths = {\n",
       "      \"vega\": \"https://cdn.jsdelivr.net/npm//vega@5?noext\",\n",
       "      \"vega-lib\": \"https://cdn.jsdelivr.net/npm//vega-lib?noext\",\n",
       "      \"vega-lite\": \"https://cdn.jsdelivr.net/npm//vega-lite@4.0.2?noext\",\n",
       "      \"vega-embed\": \"https://cdn.jsdelivr.net/npm//vega-embed@6?noext\",\n",
       "    };\n",
       "\n",
       "    function loadScript(lib) {\n",
       "      return new Promise(function(resolve, reject) {\n",
       "        var s = document.createElement('script');\n",
       "        s.src = paths[lib];\n",
       "        s.async = true;\n",
       "        s.onload = () => resolve(paths[lib]);\n",
       "        s.onerror = () => reject(`Error loading script: ${paths[lib]}`);\n",
       "        document.getElementsByTagName(\"head\")[0].appendChild(s);\n",
       "      });\n",
       "    }\n",
       "\n",
       "    function showError(err) {\n",
       "      outputDiv.innerHTML = `<div class=\"error\" style=\"color:red;\">${err}</div>`;\n",
       "      throw err;\n",
       "    }\n",
       "\n",
       "    function displayChart(vegaEmbed) {\n",
       "      vegaEmbed(outputDiv, spec, embedOpt)\n",
       "        .catch(err => showError(`Javascript Error: ${err.message}<br>This usually means there's a typo in your chart specification. See the javascript console for the full traceback.`));\n",
       "    }\n",
       "\n",
       "    if(typeof define === \"function\" && define.amd) {\n",
       "      requirejs.config({paths});\n",
       "      require([\"vega-embed\"], displayChart, err => showError(`Error loading script: ${err.message}`));\n",
       "    } else if (typeof vegaEmbed === \"function\") {\n",
       "      displayChart(vegaEmbed);\n",
       "    } else {\n",
       "      loadScript(\"vega\")\n",
       "        .then(() => loadScript(\"vega-lite\"))\n",
       "        .then(() => loadScript(\"vega-embed\"))\n",
       "        .catch(showError)\n",
       "        .then(() => displayChart(vegaEmbed));\n",
       "    }\n",
       "  })({\"config\": {\"view\": {\"continuousWidth\": 400, \"continuousHeight\": 300}}, \"hconcat\": [{\"data\": {\"name\": \"data-8b48908db5d271e2956e2a8272fc8e07\"}, \"mark\": {\"type\": \"line\", \"color\": \"red\"}, \"encoding\": {\"x\": {\"type\": \"quantitative\", \"field\": \"x\"}, \"y\": {\"type\": \"quantitative\", \"field\": \"y\"}}, \"title\": \"Normal Distribution\", \"width\": 360}, {\"data\": {\"name\": \"data-c7ad18303904753003e679bfb19caa8f\"}, \"mark\": \"bar\", \"encoding\": {\"x\": {\"type\": \"quantitative\", \"field\": \"Number of Heads\"}, \"y\": {\"type\": \"quantitative\", \"field\": \"Number of Ways\"}}, \"title\": \"Distribution of Heads\", \"width\": 360}], \"$schema\": \"https://vega.github.io/schema/vega-lite/v4.0.2.json\", \"datasets\": {\"data-8b48908db5d271e2956e2a8272fc8e07\": [{\"x\": 0, \"y\": 2.600281868827203e-10}, {\"x\": 1, \"y\": 1.8276568877457265e-09}, {\"x\": 2, \"y\": 1.1623567955302993e-08}, {\"x\": 3, \"y\": 6.688901526032106e-08}, {\"x\": 4, \"y\": 3.4828975312041506e-07}, {\"x\": 5, \"y\": 1.6409567867287274e-06}, {\"x\": 6, \"y\": 6.995586696268014e-06}, {\"x\": 7, \"y\": 2.6984954724388374e-05}, {\"x\": 8, \"y\": 9.418674667969561e-05}, {\"x\": 9, \"y\": 0.000297459915550561}, {\"x\": 10, \"y\": 0.000850036660252035}, {\"x\": 11, \"y\": 0.0021979480031862693}, {\"x\": 12, \"y\": 0.005142422126351767}, {\"x\": 13, \"y\": 0.010886507726916078}, {\"x\": 14, \"y\": 0.02085355003628301}, {\"x\": 15, \"y\": 0.03614447853363626}, {\"x\": 16, \"y\": 0.05668582612248957}, {\"x\": 17, \"y\": 0.0804410163156249}, {\"x\": 18, \"y\": 0.10328830949345566}, {\"x\": 19, \"y\": 0.12000389484301362}, {\"x\": 20, \"y\": 0.126156626101008}, {\"x\": 21, \"y\": 0.12000389484301362}, {\"x\": 22, \"y\": 0.10328830949345566}, {\"x\": 23, \"y\": 0.0804410163156249}, {\"x\": 24, \"y\": 0.05668582612248957}, {\"x\": 25, \"y\": 0.03614447853363626}, {\"x\": 26, \"y\": 0.02085355003628301}, {\"x\": 27, \"y\": 0.010886507726916078}, {\"x\": 28, \"y\": 0.005142422126351767}, {\"x\": 29, \"y\": 0.0021979480031862693}, {\"x\": 30, \"y\": 0.000850036660252035}, {\"x\": 31, \"y\": 0.000297459915550561}, {\"x\": 32, \"y\": 9.418674667969561e-05}, {\"x\": 33, \"y\": 2.6984954724388374e-05}, {\"x\": 34, \"y\": 6.995586696268014e-06}, {\"x\": 35, \"y\": 1.6409567867287274e-06}, {\"x\": 36, \"y\": 3.4828975312041506e-07}, {\"x\": 37, \"y\": 6.688901526032106e-08}, {\"x\": 38, \"y\": 1.1623567955302993e-08}, {\"x\": 39, \"y\": 1.8276568877457265e-09}], \"data-c7ad18303904753003e679bfb19caa8f\": [{\"Number of Heads\": 0, \"Number of Ways\": 0}, {\"Number of Heads\": 1, \"Number of Ways\": 0}, {\"Number of Heads\": 2, \"Number of Ways\": 0}, {\"Number of Heads\": 3, \"Number of Ways\": 0}, {\"Number of Heads\": 4, \"Number of Ways\": 0}, {\"Number of Heads\": 5, \"Number of Ways\": 0}, {\"Number of Heads\": 6, \"Number of Ways\": 0}, {\"Number of Heads\": 7, \"Number of Ways\": 0}, {\"Number of Heads\": 8, \"Number of Ways\": 0}, {\"Number of Heads\": 9, \"Number of Ways\": 0}, {\"Number of Heads\": 10, \"Number of Ways\": 2}, {\"Number of Heads\": 11, \"Number of Ways\": 9}, {\"Number of Heads\": 12, \"Number of Ways\": 16}, {\"Number of Heads\": 13, \"Number of Ways\": 28}, {\"Number of Heads\": 14, \"Number of Ways\": 54}, {\"Number of Heads\": 15, \"Number of Ways\": 124}, {\"Number of Heads\": 16, \"Number of Ways\": 174}, {\"Number of Heads\": 17, \"Number of Ways\": 257}, {\"Number of Heads\": 18, \"Number of Ways\": 326}, {\"Number of Heads\": 19, \"Number of Ways\": 379}, {\"Number of Heads\": 20, \"Number of Ways\": 351}, {\"Number of Heads\": 21, \"Number of Ways\": 351}, {\"Number of Heads\": 22, \"Number of Ways\": 289}, {\"Number of Heads\": 23, \"Number of Ways\": 249}, {\"Number of Heads\": 24, \"Number of Ways\": 156}, {\"Number of Heads\": 25, \"Number of Ways\": 109}, {\"Number of Heads\": 26, \"Number of Ways\": 68}, {\"Number of Heads\": 27, \"Number of Ways\": 36}, {\"Number of Heads\": 28, \"Number of Ways\": 15}, {\"Number of Heads\": 29, \"Number of Ways\": 5}, {\"Number of Heads\": 30, \"Number of Ways\": 0}, {\"Number of Heads\": 31, \"Number of Ways\": 2}, {\"Number of Heads\": 32, \"Number of Ways\": 0}, {\"Number of Heads\": 33, \"Number of Ways\": 0}, {\"Number of Heads\": 34, \"Number of Ways\": 0}, {\"Number of Heads\": 35, \"Number of Ways\": 0}, {\"Number of Heads\": 36, \"Number of Ways\": 0}, {\"Number of Heads\": 37, \"Number of Ways\": 0}, {\"Number of Heads\": 38, \"Number of Ways\": 0}, {\"Number of Heads\": 39, \"Number of Ways\": 0}, {\"Number of Heads\": 40, \"Number of Ways\": 0}]}}, {\"mode\": \"vega-lite\"});\n",
       "</script>"
      ],
      "text/plain": [
       "alt.HConcatChart(...)"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "x = np.arange(0,num_of_flips)\n",
    "mean = num_of_flips / 2\n",
    "stddev = np.sqrt(mean / 2)\n",
    "y = norm.pdf(x,mean,stddev)\n",
    "\n",
    "data = pd.DataFrame({\n",
    "    'x': x,\n",
    "    'y': y\n",
    "})\n",
    "normal_chart = alt.Chart(data).mark_line(color='red').encode(\n",
    "    alt.X('x'),\n",
    "    y = 'y'\n",
    ").properties(title='Normal Distribution',width=360)\n",
    "normal_chart | bar_chart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By comparing the charts we realize that\n",
    "\n",
    "- As the sample size become larger, the distribution of samples approximate a normal distribution.\n",
    "\n",
    "- As we flip the coin repeatedly, number of heads and tails are getting equal because the number of heads has been peaked around 20 flips out of 40, which indicates that the probability of heads is close to $0.5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#hide\n",
    "# import altair as alt\n",
    "# import pandas as pd\n",
    "# import numpy as np\n",
    "\n",
    "# rand = np.random.RandomState(42)\n",
    "\n",
    "# df = pd.DataFrame({\n",
    "#     'xval': range(100),\n",
    "#     'yval': rand.randn(100).cumsum()\n",
    "# })\n",
    "\n",
    "# slider = alt.binding_range(min=0, max=100, step=1, name='cutoff:')\n",
    "# selector = alt.selection_single(name=\"SelectorName\", fields=['cutoff'],\n",
    "#                                 bind=slider, init={'cutoff': 50})\n",
    "\n",
    "# alt.Chart(df).mark_point().encode(\n",
    "#     x='xval',\n",
    "#     y='yval',\n",
    "#     color=alt.condition(\n",
    "#         alt.datum.xval < selector.cutoff,\n",
    "#         alt.value('red'), alt.value('blue')\n",
    "#     )\n",
    "# ).add_selection(\n",
    "#     selector\n",
    "# )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "[1] Tom Mitchell [Estimating Probabilities](http://www.cs.cmu.edu/~tom/mlbook/Joint_MLE_MAP.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}