{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Conjugacy\n",
    "\n",
    "<hr>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "nbsphinx": "hidden",
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Colab setup ------------------\n",
    "import os, sys, subprocess\n",
    "if \"google.colab\" in sys.modules:\n",
    "    cmd = \"pip install --upgrade bebi103 watermark\"\n",
    "    process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "    stdout, stderr = process.communicate()\n",
    "    data_path = \"https://s3.amazonaws.com/bebi103.caltech.edu/data/\"\n",
    "else:\n",
    "    data_path = \"../data/\"\n",
    "# ------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "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=\"a6726bf3-d8e4-4105-afad-89dff38ac75a\">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",
       "    function drop(id) {\n",
       "      const view = Bokeh.index.get_by_id(id)\n",
       "      if (view != null) {\n",
       "        view.model.document.clear()\n",
       "        Bokeh.index.delete(view)\n",
       "      }\n",
       "    }\n",
       "\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",
       "\n",
       "    // Clean up Bokeh references\n",
       "    if (id != null) {\n",
       "      drop(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",
       "            drop(id)\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(\"a6726bf3-d8e4-4105-afad-89dff38ac75a\");\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.3.0.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.3.0.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.3.0.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.3.0.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-mathjax-3.3.0.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(\"a6726bf3-d8e4-4105-afad-89dff38ac75a\")).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(\"a6726bf3-d8e4-4105-afad-89dff38ac75a\");\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.3.0.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.3.0.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.3.0.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.3.0.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-mathjax-3.3.0.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(\"a6726bf3-d8e4-4105-afad-89dff38ac75a\")).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"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "\n",
    "import bokeh.plotting\n",
    "import bokeh.io\n",
    "bokeh.io.output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "We have talked about Bayes theorem as a model for learning. The idea there was that we know something before (*a priori*) acquiring data, and then we update our knowledge after (*a posteriori*). So, we come in with the prior and out with the posterior after acquiring data. It might make sense, then, that the prior and the posterior distributions have the same functional form. That is, the prior and the posterior have the same distribution, and the parameters of the distribution are updated from the prior to the posterior by the likelihood. When the prior and posterior have the same functional form, the prior is said to be **conjugate** to the likelihood. This *seems* pleasing: the likelihood serves to update the prior into the posterior, so it should determine the functional form of the prior/posterior such that they are the same.\n",
    "\n",
    "Conjugate priors are especially useful because they enable analytical calculation of the posterior, typically without have to do any mathematics yourself. This is perhaps best seen through example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Beta-Binomial conjugate pair\n",
    "\n",
    "As a motivating example of use of a conjugate prior, we will use data from one of the experiments from a paper by [Mossman, et al., 2019](https://doi.org/10.1093/jhered/esz009), in which the authors investigated the age of *Drosophila* parents on the viability of the eggs of their offspring. In one vial, they mated young (≤ 3 days old) males with old females (≥ 45 days old). Of the 176 eggs laid by their offspring over a 10-day period, 94 of them hatched, while the remainder failed to hatch. In another vial, they mated young males with young females. Of the 190 eggs laid by their offspring over a 10-day period, 154 hatched. (They took many more data than these, but we're using this experiment as a demonstration of conjugacy.)\n",
    "\n",
    "The story behind the data matches that of the Binomial distribution. Each egg can by thought of as a Bernoulli trial, with success being hatching and failure being failure to hatch. The number $n$ of hatches out of $N$ eggs is Binomially distributed.\n",
    "\n",
    "\\begin{align}\n",
    "n \\mid N, \\theta \\sim \\text{Binom}(N, \\theta).\n",
    "\\end{align}\n",
    "\n",
    "Writing out the PMF, this is\n",
    "\n",
    "\\begin{align}\n",
    "f(n\\mid N, \\theta) = \\begin{pmatrix}N\\\\n\\end{pmatrix} \\theta^n(1-\\theta)^{N-n}.\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "This identifies a single parameter to estimate, $\\theta$, which is the probability of hatching. Our goal, then, is to compute $g(\\theta\\mid n,N)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding the conjugate\n",
    "\n",
    "We seek a prior that is conjugate to the likelihood. That is, we want to choose $g(\\theta)$ such that $g(\\theta \\mid n, N)$ and $g(\\theta)$ have the same functional form. Intuitively, we might try a prior that has a similar functional form as the likelihood. The Beta distribution is one such distribution. Its PDF is\n",
    "\n",
    "\\begin{align}\n",
    "g(\\theta \\mid \\alpha, \\beta) = \\frac{\\theta^{\\alpha-1}(1-\\theta)^{\\beta-1}}{B(\\alpha,\\beta)},\n",
    "\\end{align}\n",
    "\n",
    "where $B(\\alpha,\\beta)$ is a [beta function](https://en.wikipedia.org/wiki/Beta_function),\n",
    "\n",
    "\\begin{align}\n",
    "B(\\alpha,\\beta) = \\frac{\\Gamma(\\alpha)\\Gamma(\\beta)}{\\Gamma(\\alpha+\\beta)}.\n",
    "\\end{align}\n",
    "\n",
    "Note that there are two parameters, $\\alpha$ and $\\beta$ that parametrize the prior. If in fact the Beta distribution is conjugate to the Binomial, it is these two parameters that get updated by the data.\n",
    "\n",
    "You should explore the Beta distribution in the [Distribution Explorer](https://distribution-explorer.github.io/continuous/beta.html) to get a feel for its shape and how the parameters $\\alpha$ and $\\beta$ affect its shape. Importantly, if $\\alpha = \\beta = 1$, we get a Uniform distribution on the interval [0, 1].\n",
    "\n",
    "We can now check to see if the Beta distribution is conjugate to the Binomial. If we insert a Beta distribution for the prior, we have\n",
    "\n",
    "\\begin{align}\n",
    "g(\\theta \\mid n, N, \\alpha, \\beta) &= \\frac{f(n\\mid N, \\theta)\\, g(\\theta\\mid \\alpha, \\beta)}{f(n \\mid N)} \\\\[1em]\n",
    "&= \\frac{1}{f(n\\mid N)}\\,\\begin{pmatrix}N\\\\n\\end{pmatrix}\\,\\theta^n(1-\\theta)^{N-n}\\,\\frac{\\theta^{\\alpha-1}(1-\\theta)^{\\beta-1}}{B(\\alpha,\\beta)} \\\\[1em]\n",
    "&= \\frac{1}{f(n\\mid N)\\,B(\\alpha,\\beta)}\\,\\begin{pmatrix}N\\\\n\\end{pmatrix}\\,\\theta^{n+\\alpha-1}(1-\\theta)^{N-n+\\beta-1}.\n",
    "\\end{align}\n",
    "\n",
    "In looking at this expression, the only bit that depends on $\\theta$ is $\\theta^{n+\\alpha-1}(1-\\theta)^{N-n+\\beta-1}$, which is exactly the $\\theta$-dependence of a Beta distribution with parameters $n+\\alpha$ and $N-n+\\beta$. Thus, the posterior must also be a Beta distribution. Therefore,\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{1}{f(n\\mid N)\\,B(\\alpha,\\beta)}\\,\\begin{pmatrix}N\\\\n\\end{pmatrix} = \\frac{1}{B(n+\\alpha, N-n+\\beta)}.\n",
    "\\end{align}\n",
    "\n",
    "We have just normalized the posterior without doing any nasty integrals! So, the posterior is\n",
    "\n",
    "\\begin{align}\n",
    "g(\\theta \\mid n, N, \\alpha, \\beta) = \\frac{\\theta^{n+\\alpha-1}(1-\\theta)^{N-n+\\beta-1}}{B(n+\\alpha, N-n+\\beta)}.\n",
    "\\end{align}\n",
    "\n",
    "Thus, we have prior and posterior distributions\n",
    "\n",
    "\\begin{align}\n",
    "&\\theta \\sim \\text{Beta}(\\alpha, \\beta),\\\\[1em]\n",
    "&\\theta \\mid n, N \\sim \\text{Beta}(n+\\alpha, N-n+\\beta).\n",
    "\\end{align}\n",
    "\n",
    "So, we can see that conjugacy is useful. For a given likelihood, if we know its conjugate prior, we can just immediately write down the posterior in a clear form. The [Wikipedia page on conjugate priors](https://en.wikipedia.org/wiki/Conjugate_prior) has a useful table of likelihood-conjugate pairs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plots of the posteriors\n",
    "\n",
    "We can now rapidly exactly plot the posteriors for the egg hatch probability for old mothers and young mothers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "  <div id=\"d6d966f0-a855-4894-8431-bd453806859d\" 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 = {\"9e75b192-2a7b-4314-8639-4b0a8ccf9b0c\":{\"version\":\"3.3.0\",\"title\":\"Bokeh Application\",\"roots\":[{\"type\":\"object\",\"name\":\"Figure\",\"id\":\"p1001\",\"attributes\":{\"x_range\":{\"type\":\"object\",\"name\":\"Range1d\",\"id\":\"p1010\"},\"y_range\":{\"type\":\"object\",\"name\":\"DataRange1d\",\"id\":\"p1003\"},\"x_scale\":{\"type\":\"object\",\"name\":\"LinearScale\",\"id\":\"p1011\"},\"y_scale\":{\"type\":\"object\",\"name\":\"LinearScale\",\"id\":\"p1012\"},\"title\":{\"type\":\"object\",\"name\":\"Title\",\"id\":\"p1008\"},\"renderers\":[{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1036\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1030\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1031\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1032\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAAAAACASAEiBYhkP4BIASIFiHQ/wOwBswfMfj+ASAEiBYiEP6CagWoGqok/wOwBswfMjj9wH8F9BPeRP4BIASIFiJQ/kHFBxgUZlz+gmoFqBqqZP7DDwQ4HO5w/wOwBswfMnj/oCqErhK6gP3AfwX0E96E/+DPhz4Q/oz+ASAEiBYikPwhdIXSF0KU/kHFBxgUZpz8YhmEYhmGoP6CagWoGqqk/KK+hvIbyqj+ww8EOBzusPzjY4WCHg60/wOwBswfMrj+kAJECRAqwP+gKoSuErrA/LBWxVMRSsT9wH8F9BPexP7Qp0aZEm7I/+DPhz4Q/sz88PvH4xOOzP4BIASIFiLQ/xFIRS0UstT8IXSF0hdC1P0xnMZ3FdLY/kHFBxgUZtz/Ue1HvRb23PxiGYRiGYbg/XJBxQcYFuT+gmoFqBqq5P+SkkZNGTro/KK+hvIbyuj9subHlxpa7P7DDwQ4HO7w/9M3RN0ffvD842OFgh4O9P3zi8YnHJ74/wOwBswfMvj8E9xHcR3C/P6QAkQJECsA/xgUZF2RcwD/oCqErhK7APwoQKUCkAME/LBWxVMRSwT9OGjlp5KTBP3AfwX0E98E/kiRJkiRJwj+0KdGmRJvCP9YuWbtk7cI/+DPhz4Q/wz8aOWnkpJHDPzw+8fjE48M/XkN5DeU1xD+ASAEiBYjEP6JNiTYl2sQ/xFIRS0UsxT/mV5lfZX7FPwhdIXSF0MU/KmKpiKUixj9MZzGdxXTGP25subHlxsY/kHFBxgUZxz+ydsnaJWvHP9R7Ue9Fvcc/9oDZA2YPyD8YhmEYhmHIPzqL6Syms8g/XJBxQcYFyT9+lflV5lfJP6CagWoGqsk/wp8Jfyb8yT/kpJGTRk7KPwaqGahmoMo/KK+hvIbyyj9KtCnRpkTLP2y5seXGlss/jr45+uboyz+ww8EOBzvMP9LISSMnjcw/9M3RN0ffzD8W01lMZzHNPzjY4WCHg80/Wt1pdafVzT984vGJxyfOP57neZ7nec4/wOwBswfMzj/i8YnHJx7PPwT3EdxHcM8/JvyZ8GfCzz+kAJECRArQPzUD1QxUM9A/xgUZF2Rc0D9XCF0hdIXQP+gKoSuErtA/eQ3lNZTX0D8KEClApADRP5sSbUq0KdE/LBWxVMRS0T+9F/Ve1HvRP04aOWnkpNE/3xx9c/TN0T9wH8F9BPfRPwEiBYgUINI/kiRJkiRJ0j8jJ42cNHLSP7Qp0aZEm9I/RSwVsVTE0j/WLlm7ZO3SP2cxncV0FtM/+DPhz4Q/0z+JNiXalGjTPxo5aeSkkdM/qzut7rS60z88PvH4xOPTP81ANQPVDNQ/XkN5DeU11D/vRb0X9V7UP4BIASIFiNQ/EUtFLBWx1D+iTYk2JdrUPzNQzUA1A9U/xFIRS0Us1T9VVVVVVVXVP+ZXmV9lftU/d1rdaXWn1T8IXSF0hdDVP5lfZX6V+dU/KmKpiKUi1j+7ZO2StUvWP0xnMZ3FdNY/3Wl1p9Wd1j9ubLmx5cbWP/9u/bv179Y/kHFBxgUZ1z8hdIXQFULXP7J2ydola9c/Q3kN5TWU1z/Ue1HvRb3XP2V+lflV5tc/9oDZA2YP2D+Hgx0OdjjYPxiGYRiGYdg/qYilIpaK2D86i+ksprPYP8uNLTe23Ng/XJBxQcYF2T/tkrVL1i7ZP36V+VXmV9k/D5g9YPaA2T+gmoFqBqrZPzGdxXQW09k/wp8Jfyb82T9Tok2JNiXaP+SkkZNGTto/dafVnVZ32j8GqhmoZqDaP5esXbJ2ydo/KK+hvIby2j+5seXGlhvbP0q0KdGmRNs/27Zt27Zt2z9subHlxpbbP/279e/Wv9s/jr45+ubo2z8fwX0E9xHcP7DDwQ4HO9w/QcYFGRdk3D/SyEkjJ43cP2PLjS03ttw/9M3RN0ff3D+F0BVCVwjdPxbTWUxnMd0/p9WdVnda3T842OFgh4PdP8naJWuXrN0/Wt1pdafV3T/r361/t/7dP3zi8YnHJ94/DeU1lNdQ3j+e53me53nePy/qvaj3ot4/wOwBswfM3j9R70W9F/XeP+LxiccnHt8/c/TN0TdH3z8E9xHcR3DfP5X5VeZXmd8/JvyZ8GfC3z+3/t36d+vfP6QAkQJECuA/7AGzB8we4D81A9UMVDPgP34E9xHcR+A/xgUZF2Rc4D8OBzsc7HDgP1cIXSF0heA/oAl/JvyZ4D/oCqErhK7gPzAMwzAMw+A/eQ3lNZTX4D/CDgc7HOzgPwoQKUCkAOE/UhFLRSwV4T+bEm1KtCnhP+QTj088PuE/LBWxVMRS4T90FtNZTGfhP70X9V7Ue+E/BhkXZFyQ4T9OGjlp5KThP5YbW25sueE/3xx9c/TN4T8oHp94fOLhP3AfwX0E9+E/uCDjgowL4j8BIgWIFCDiP0ojJ42cNOI/kiRJkiRJ4j/aJWuXrF3iPyMnjZw0cuI/bCivobyG4j+0KdGmRJviP/wq86vMr+I/RSwVsVTE4j+OLTe23NjiP9YuWbtk7eI/HjB7wOwB4z9nMZ3FdBbjP7Ayv8r8KuM/+DPhz4Q/4z9ANQPVDFTjP4k2JdqUaOM/0jdH3xx94z8aOWnkpJHjP2I6i+kspuM/qzut7rS64z/0PM/zPM/jPzw+8fjE4+M/hD8T/kz44z/NQDUD1QzkPxZCVwhdIeQ/XkN5DeU15D+mRJsSbUrkP+9FvRf1XuQ/OEffHH1z5D+ASAEiBYjkP8hJIyeNnOQ/EUtFLBWx5D9aTGcxncXkP6JNiTYl2uQ/6k6rO63u5D8zUM1ANQPlP3xR70W9F+U/xFIRS0Us5T8MVDNQzUDlP1VVVVVVVeU/nlZ3Wt1p5T/mV5lfZX7lPy5Zu2TtkuU/d1rdaXWn5T/AW/9u/bvlPwhdIXSF0OU/UF5DeQ3l5T+ZX2V+lfnlP+Jgh4MdDuY/KmKpiKUi5j9yY8uNLTfmP7tk7ZK1S+Y/BGYPmD1g5j9MZzGdxXTmP5RoU6JNieY/3Wl1p9Wd5j8ma5esXbLmP25subHlxuY/tm3btm3b5j//bv279e/mP0hwH8F9BOc/kHFBxgUZ5z/YcmPLjS3nPyF0hdAVQuc/anWn1Z1W5z+ydsnaJWvnP/p369+tf+c/Q3kN5TWU5z+Mei/qvajnP9R7Ue9Fvec/HH1z9M3R5z9lfpX5VebnP65/t/7d+uc/9oDZA2YP6D8+gvsI7iPoP4eDHQ52OOg/0IQ/E/5M6D8YhmEYhmHoP2CHgx0Odug/qYilIpaK6D/yiccnHp/oPzqL6Syms+g/gowLMi7I6D/LjS03ttzoPxSPTzw+8eg/XJBxQcYF6T+kkZNGThrpP+2StUvWLuk/NpTXUF5D6T9+lflV5lfpP8aWG1tubOk/D5g9YPaA6T9YmV9lfpXpP6CagWoGquk/6Jujb46+6T8xncV0FtPpP3qe53me5+k/wp8Jfyb86T8KoSuErhDqP1OiTYk2Jeo/nKNvjr456j/kpJGTRk7qPyyms5jOYuo/dafVnVZ36j++qPei3ovqPwaqGahmoOo/Tqs7re606j+XrF2ydsnqP+Ctf7f+3eo/KK+hvIby6j9wsMPBDgfrP7mx5caWG+s/ArMHzB4w6z9KtCnRpkTrP5K1S9YuWes/27Zt27Zt6z8kuI/gPoLrP2y5seXGlus/tLrT6k6r6z/9u/Xv1r/rP0a9F/Ve1Os/jr45+ubo6z/Wv1v/bv3rPx/BfQT3Eew/aMKfCX8m7D+ww8EOBzvsP/jE4xOPT+w/QcYFGRdk7D+Kxycen3jsP9LISSMnjew/GsprKK+h7D9jy40tN7bsP6zMrzK/yuw/9M3RN0ff7D88z/M8z/PsP4XQFUJXCO0/ztE3R98c7T8W01lMZzHtP17Ue1HvRe0/p9WdVnda7T/w1r9b/27tPzjY4WCHg+0/gNkDZg+Y7T/J2iVrl6ztPxLcR3Afwe0/Wt1pdafV7T+i3ot6L+rtP+vfrX+3/u0/NOHPhD8T7j984vGJxyfuP8TjE49PPO4/DeU1lNdQ7j9W5leZX2XuP57neZ7nee4/5uibo2+O7j8v6r2o96LuP3jr361/t+4/wOwBswfM7j8I7iO4j+DuP1HvRb0X9e4/mvBnwp8J7z/i8YnHJx7vPyrzq8yvMu8/c/TN0TdH7z+89e/Wv1vvPwT3EdxHcO8/TPgz4c+E7z+V+VXmV5nvP976d+vfre8/JvyZ8GfC7z9u/bv179bvP7f+3fp36+8/AAAAAAAA8D8=\"},\"shape\":[400],\"dtype\":\"float64\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAAAAABwcG8H+EpUGJndEDOdgjAehuj+286ZmiGmCT77y9IFJAjkYFcrP+UlukKOaW+EbCecPM8yEqq1KAa0rInQ8dIpq3nkFGwkzirE1dHkTtKtK4JX5UzX7XYskLZpQSAnMC1la1DcCqHXLecmGBeIxnMugixaDUN0BC+tFh+5TdKLL6SJct/iLQowzvuctSfGgTBNCfnZ6gnyMBfyQ3TYMVwxSLSlFZRnwTGH3K7YQFkhMqxpbXpjdXwyB7unweOG0zLh4kjb174mM0DAgPtix3YzUvtRRkvWwzPMoYPWmFkONPu1hm59lFQ0SVW7aS3xmDQ4IQPFejbbNOp7pjXs5xo1PHCbHd5AWDWjrCuMkwqUNXTJNCPUg841jw8knkqABTZ+bU9zYSk8NkU88yopNXE2Thr5JNKwozYWerOF+SrVNmOO/IqDcQU3mSHqKheHNDcyG7QYe55iN0uyxo1KCpA3VpmbAU1Pujf+s4kzsZXkN4kAmdluyg44g7+HLOYNNjgsXBZgVU9eODG9uErCA4Q45qL9xqZwqTj0mkhuvirPOFP1c5Q1bfI4CAq6DosOFTkqP+z/cEg3Ocx1Qo578Fg5w5fwQRfoeTnn/lHHHSCaObDY8Nkbmrk53p6T7jVn2DlbOxGeDaX2Od5G9mNceRQ69xRSvBoNMjrUEiwo+g9POsQfOP67GWo6/NdtOtFuhTo6sdVlIjahOsrVBabKDLs6yfG8ZxXQ1DqAAK1FuWDvOocbrleVLwc7iebeRKPNIDsHVSLbOeY3Ox8xivaZr1A7kGbq55biZjtlL6lf1dd+O/8upW29bpQ78Ma32vyfqjvYfwifQxHBO1kXv/WCidU77RlcRnPC6jvXu4si2F8APD+lFvyxvhM8O+PEfJx3Jzxvw2rE2n87PFRJBeQMyE88L1OOojQeYjzSIN2moWJ0PDUXM93to4Y8x2C2T7jTmDzBTy2QPuOqPBGMxyoKxLw800CzBqBozjynpwGzJcXfPM4vi634Z/A86k9YbP3AAD3BPk1wk+sQPecuHyo05yA9NGiKgtu0MD0aecIu71ZAPe0LZCQpok89nGqps/FPXj3l++NbH8JsPTYN9D6aBHs9iqmetb4jiT2DKDpv2CuXPRVnuxCsKKU9rZULWBQlsz32XZ3itSrBPf7g/6+Vg849WzD0EAzi2j3tnp5QF3vnPYGA7YX5VfQ90qWUqCN3AT4br+mUAsENPoAMGS2mIxk+vNdWuSESJT4SKk22mYUxPueKVnlM6jw++n2rOgStRz4x82gTSj1TPmM6HgsXCV8+RnaqFA3ZaD7H3xII1r9zPtHUkDo7K38+L1g4ROZriD5vzNxrOQCTPmCl/nPLXJ0+y3USk0OIpj79pHyMliyxPj34LOh1Abo+xSYRUDiPwz53dscv8DrNPgcZaDRJs9U+cn2QBOMB4D45FWRje3fnPtKoh6r0F/E+cfMb65+/+D4vbiB4Xs4BP0sfiumUdwk/wKJ4blIaEj+UX600SpUZPyk0LiaJ+CE/J49LLpYZKT9jqac0ZW0xP5+W2aTGDzg//V2iQTiEQD/GIM4EDYxGP2Wm9C2gm04/awskjw+pVD+WE41eFb1bP3rGTm4DhWI/xOIo6XeYaD/U1s0nuj5wP5s46aJbWHU/8UlweQLmez+kgLk66yKCP4xMtAIzdYc/RABN0hMvjj/f1mzi1lGTPwy15HwSm5g/QDrxi/Utnz+o1m0ty6ejPznlUVFWqKg/0iomNm7Hrj/U55zvfx2zP7si10JZoLc/mM3nlj0PvT+EJHxnx8jBPz4B8bKkqcU/RMYrAJNCyj+dK1hSUq7PP1/fGir/BNM/eO5mE3u61j+OhDTIWAjbP4METE4gAOA/oaBR+U3a4j9qiZe/LhzmPwymlVBez+k/pWVvxHH97T+jgtgT5FfxP4CWE2Sq9/M/r9Zx5THi9j9OF5rrLBv6P9iPRPS3pf0/oDTjJxvCAEBlwoYrF9wCQLr10g0TIQVAXDZGzMCQB0Bo3UG/MSoKQGnfZuHI6wxAGucWjS/TD0AsVrOgpm4RQA7LTHshAxNAJxCjPbWkFEAicUYfrVAWQKyh0Oj2AxhAKczmpSm7GUBhhnXMjnIbQKlIAsotJh1AdWrn09nRHkAnh43hoDggQOmIQtgAACFAyt9Y+tq8IUB7Z9SOB20iQH1UTDBzDiNAH4nCASmfI0A26uOaXB0kQF8OqmZzhyRATPRDMQ3cJEBJ7XamCxolQNgxO4iYQCVANrX/bSpPJUCcB3bnh0UlQL/Vq+nIIyVApPscelbqJECrGtSY6JkkQISHHniCMyRAj2VIHW24I0AKesmRMCojQPsO0NaLiiJA7obK1WvbIUDiMSSP4R4hQAoIqMwXVyBAyUSQPZEMH0Bub8LWZV0dQHWtHaMgpRtA4JfwQRfoGUDD3to3cCoYQEV9G2YTcBZA8Crfvpy8FEBkwXprURMTQDDD44MXdxFAEfW7wuDUD0CqwAIE694MQAoOQsCvDwpASZqdIMNpB0BDrlPW8u4EQLb4RvtQoAJAxe0r5kF+AEDRM6yhGBH9P+MaG1rWfPk/uFvTD0E99j/cw4HIEU/zP+Ug9QJdrvA/dVVLUHCt7D8rrsgduoboP7IzuNGV3uQ/tsHUL8iq4T+fLBdKIsLdP8a6LI6y7tg/45XR0K/H1D/o90+1YzrRP/oW+LREasw/uOkRfNlOxz8q8VBlDgTDP0dRswIV3L4/sIcghzHnuD9dgFNEPPyzP7zYrhgN5q8/0QaAQlJQqT8IcMYYbPmjPyleRYuDV58/1WAbFIZymD+AaZ81U/WSP3d2xy/wOo0/4C26P1tmhj+QMYng+g+BP+uoZzb71Xk/iAqwJLpwcz9nNhfiHBNtP1N1NiIHm2U/D6NEjDfoXz8c6OR+xWhXPwoHDfs5EFE/aySlMWa2SD8hmT7ezcZBP4r4rcsqZzk/vExbIv4GMj8JLV16CGkpP3hvKmehyCE/BkgCGnq3GD9hmpXtsg0RP0ex9S+uXAc/qn9fevTE/z4AZG8vv3D1PtifP/KUuOw+rFc9hP4W4z7GQacUAy7ZPlUN0c65edA+auJh+FtjxT6H7AUNXYq7PloPynTZlbE+CuA4Pg5Fpj6x5YrZavabPmU2bxNnZ5E+YKz6A3x5hT6y0Ix64EJ6PhcINsr00m8+bU5y7uIaYz40UoZAULlWPv6+UNpHxUo+4FsrOG47Pz57CfnZ6gkyPpkaH4C2oCQ+NGHyos9YFz4es4jyGyYKPqRxxgXN+fw9wn1GPazC7z2FVe36zDbhPTGtk64fc9I9L5fH0faLwz0P+rMvUHe0PUjPZx8KLKU95l5umW+ilT37mdBSutSFPX0/dMd/v3U9Ndei2f5hZT0DmXyRRb5UPTru1n4p2UM9QsxADRG6Mj1dHXBckGohPX4vMNzA6w8978yfZHbQ/Dy14XBmU5zpPNc93fGlZ9Y8Et+5CbRIwzw2+VHyrlKwPCd1Z3zhKZs8UwbvuB01hjy/I1T4htRxPEioXwecGlw8JakueK66RTwZB6hxIHkwPEuHqJ3leRg8eYBt9OfOATxHlhd8z1zpO99qwkkzqtE72+03kPoMuDszBeUwOfufO/toq81lwIQ74MYrAJNCajvogozyLjBQO2b8Fg6QbTM7S7/xoqisFjsjIqk3jLX5OiQwTc7+Sdw6BlY0OJ0svjrlRiFTcimfOvdszbWCHn86cbcTn4UCXjoP0q8MHOg7Ogw1ivVu/Bg6gVdQlRyB9Tk+6hw3gsLROUl+vJPGGKw5tzADyLo+hTl9kOO9YaVeOYOmtBUQCTU51kzMiwtqCzmnwwKY3+ngONJkG34PtLM4Ry9tnVSahTiFT20xhTdWOIpcLMttWiU4luQxm44a8zelZ8X6Aq6/Nx991HVcOog32eB6skEAUTfFKLTZ/8MVN72l13INQtk28lu6LVVhmjYv0CpiOJpYNufy59DQThQ2s+EiQ35hzTW/oUJSO2uCNdxOhWrxwjM1mq51Lvzh4TS05ZZLqtmKNKwxi634ZzA0SDk5fmTrzzOvCsyaPxVoM8ulkukRVfsyfxPjs1R9hjJ1ZtUGSasJMkCH0VuzP4MxLBrVp+S78TC47SoKA3FSMEr8hT66aKMv6+C6kZno4S4pUR2K/9kHLvMKnBpMdBEtkfzprBXU8itZeQojXR6gKochoiqlC/8of0oOI8A94ybzgdEtW/ToI6oC0AaknM8eAAAAAAAAAAA=\"},\"shape\":[400],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1037\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1038\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1033\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"#1f77b4\",\"line_width\":2}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1034\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"#1f77b4\",\"line_alpha\":0.1,\"line_width\":2}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1035\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"#1f77b4\",\"line_alpha\":0.2,\"line_width\":2}}}},{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1047\",\"attributes\":{\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1041\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1042\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1043\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAAAAACASAEiBYhkP4BIASIFiHQ/wOwBswfMfj+ASAEiBYiEP6CagWoGqok/wOwBswfMjj9wH8F9BPeRP4BIASIFiJQ/kHFBxgUZlz+gmoFqBqqZP7DDwQ4HO5w/wOwBswfMnj/oCqErhK6gP3AfwX0E96E/+DPhz4Q/oz+ASAEiBYikPwhdIXSF0KU/kHFBxgUZpz8YhmEYhmGoP6CagWoGqqk/KK+hvIbyqj+ww8EOBzusPzjY4WCHg60/wOwBswfMrj+kAJECRAqwP+gKoSuErrA/LBWxVMRSsT9wH8F9BPexP7Qp0aZEm7I/+DPhz4Q/sz88PvH4xOOzP4BIASIFiLQ/xFIRS0UstT8IXSF0hdC1P0xnMZ3FdLY/kHFBxgUZtz/Ue1HvRb23PxiGYRiGYbg/XJBxQcYFuT+gmoFqBqq5P+SkkZNGTro/KK+hvIbyuj9subHlxpa7P7DDwQ4HO7w/9M3RN0ffvD842OFgh4O9P3zi8YnHJ74/wOwBswfMvj8E9xHcR3C/P6QAkQJECsA/xgUZF2RcwD/oCqErhK7APwoQKUCkAME/LBWxVMRSwT9OGjlp5KTBP3AfwX0E98E/kiRJkiRJwj+0KdGmRJvCP9YuWbtk7cI/+DPhz4Q/wz8aOWnkpJHDPzw+8fjE48M/XkN5DeU1xD+ASAEiBYjEP6JNiTYl2sQ/xFIRS0UsxT/mV5lfZX7FPwhdIXSF0MU/KmKpiKUixj9MZzGdxXTGP25subHlxsY/kHFBxgUZxz+ydsnaJWvHP9R7Ue9Fvcc/9oDZA2YPyD8YhmEYhmHIPzqL6Syms8g/XJBxQcYFyT9+lflV5lfJP6CagWoGqsk/wp8Jfyb8yT/kpJGTRk7KPwaqGahmoMo/KK+hvIbyyj9KtCnRpkTLP2y5seXGlss/jr45+uboyz+ww8EOBzvMP9LISSMnjcw/9M3RN0ffzD8W01lMZzHNPzjY4WCHg80/Wt1pdafVzT984vGJxyfOP57neZ7nec4/wOwBswfMzj/i8YnHJx7PPwT3EdxHcM8/JvyZ8GfCzz+kAJECRArQPzUD1QxUM9A/xgUZF2Rc0D9XCF0hdIXQP+gKoSuErtA/eQ3lNZTX0D8KEClApADRP5sSbUq0KdE/LBWxVMRS0T+9F/Ve1HvRP04aOWnkpNE/3xx9c/TN0T9wH8F9BPfRPwEiBYgUINI/kiRJkiRJ0j8jJ42cNHLSP7Qp0aZEm9I/RSwVsVTE0j/WLlm7ZO3SP2cxncV0FtM/+DPhz4Q/0z+JNiXalGjTPxo5aeSkkdM/qzut7rS60z88PvH4xOPTP81ANQPVDNQ/XkN5DeU11D/vRb0X9V7UP4BIASIFiNQ/EUtFLBWx1D+iTYk2JdrUPzNQzUA1A9U/xFIRS0Us1T9VVVVVVVXVP+ZXmV9lftU/d1rdaXWn1T8IXSF0hdDVP5lfZX6V+dU/KmKpiKUi1j+7ZO2StUvWP0xnMZ3FdNY/3Wl1p9Wd1j9ubLmx5cbWP/9u/bv179Y/kHFBxgUZ1z8hdIXQFULXP7J2ydola9c/Q3kN5TWU1z/Ue1HvRb3XP2V+lflV5tc/9oDZA2YP2D+Hgx0OdjjYPxiGYRiGYdg/qYilIpaK2D86i+ksprPYP8uNLTe23Ng/XJBxQcYF2T/tkrVL1i7ZP36V+VXmV9k/D5g9YPaA2T+gmoFqBqrZPzGdxXQW09k/wp8Jfyb82T9Tok2JNiXaP+SkkZNGTto/dafVnVZ32j8GqhmoZqDaP5esXbJ2ydo/KK+hvIby2j+5seXGlhvbP0q0KdGmRNs/27Zt27Zt2z9subHlxpbbP/279e/Wv9s/jr45+ubo2z8fwX0E9xHcP7DDwQ4HO9w/QcYFGRdk3D/SyEkjJ43cP2PLjS03ttw/9M3RN0ff3D+F0BVCVwjdPxbTWUxnMd0/p9WdVnda3T842OFgh4PdP8naJWuXrN0/Wt1pdafV3T/r361/t/7dP3zi8YnHJ94/DeU1lNdQ3j+e53me53nePy/qvaj3ot4/wOwBswfM3j9R70W9F/XeP+LxiccnHt8/c/TN0TdH3z8E9xHcR3DfP5X5VeZXmd8/JvyZ8GfC3z+3/t36d+vfP6QAkQJECuA/7AGzB8we4D81A9UMVDPgP34E9xHcR+A/xgUZF2Rc4D8OBzsc7HDgP1cIXSF0heA/oAl/JvyZ4D/oCqErhK7gPzAMwzAMw+A/eQ3lNZTX4D/CDgc7HOzgPwoQKUCkAOE/UhFLRSwV4T+bEm1KtCnhP+QTj088PuE/LBWxVMRS4T90FtNZTGfhP70X9V7Ue+E/BhkXZFyQ4T9OGjlp5KThP5YbW25sueE/3xx9c/TN4T8oHp94fOLhP3AfwX0E9+E/uCDjgowL4j8BIgWIFCDiP0ojJ42cNOI/kiRJkiRJ4j/aJWuXrF3iPyMnjZw0cuI/bCivobyG4j+0KdGmRJviP/wq86vMr+I/RSwVsVTE4j+OLTe23NjiP9YuWbtk7eI/HjB7wOwB4z9nMZ3FdBbjP7Ayv8r8KuM/+DPhz4Q/4z9ANQPVDFTjP4k2JdqUaOM/0jdH3xx94z8aOWnkpJHjP2I6i+kspuM/qzut7rS64z/0PM/zPM/jPzw+8fjE4+M/hD8T/kz44z/NQDUD1QzkPxZCVwhdIeQ/XkN5DeU15D+mRJsSbUrkP+9FvRf1XuQ/OEffHH1z5D+ASAEiBYjkP8hJIyeNnOQ/EUtFLBWx5D9aTGcxncXkP6JNiTYl2uQ/6k6rO63u5D8zUM1ANQPlP3xR70W9F+U/xFIRS0Us5T8MVDNQzUDlP1VVVVVVVeU/nlZ3Wt1p5T/mV5lfZX7lPy5Zu2TtkuU/d1rdaXWn5T/AW/9u/bvlPwhdIXSF0OU/UF5DeQ3l5T+ZX2V+lfnlP+Jgh4MdDuY/KmKpiKUi5j9yY8uNLTfmP7tk7ZK1S+Y/BGYPmD1g5j9MZzGdxXTmP5RoU6JNieY/3Wl1p9Wd5j8ma5esXbLmP25subHlxuY/tm3btm3b5j//bv279e/mP0hwH8F9BOc/kHFBxgUZ5z/YcmPLjS3nPyF0hdAVQuc/anWn1Z1W5z+ydsnaJWvnP/p369+tf+c/Q3kN5TWU5z+Mei/qvajnP9R7Ue9Fvec/HH1z9M3R5z9lfpX5VebnP65/t/7d+uc/9oDZA2YP6D8+gvsI7iPoP4eDHQ52OOg/0IQ/E/5M6D8YhmEYhmHoP2CHgx0Odug/qYilIpaK6D/yiccnHp/oPzqL6Syms+g/gowLMi7I6D/LjS03ttzoPxSPTzw+8eg/XJBxQcYF6T+kkZNGThrpP+2StUvWLuk/NpTXUF5D6T9+lflV5lfpP8aWG1tubOk/D5g9YPaA6T9YmV9lfpXpP6CagWoGquk/6Jujb46+6T8xncV0FtPpP3qe53me5+k/wp8Jfyb86T8KoSuErhDqP1OiTYk2Jeo/nKNvjr456j/kpJGTRk7qPyyms5jOYuo/dafVnVZ36j++qPei3ovqPwaqGahmoOo/Tqs7re606j+XrF2ydsnqP+Ctf7f+3eo/KK+hvIby6j9wsMPBDgfrP7mx5caWG+s/ArMHzB4w6z9KtCnRpkTrP5K1S9YuWes/27Zt27Zt6z8kuI/gPoLrP2y5seXGlus/tLrT6k6r6z/9u/Xv1r/rP0a9F/Ve1Os/jr45+ubo6z/Wv1v/bv3rPx/BfQT3Eew/aMKfCX8m7D+ww8EOBzvsP/jE4xOPT+w/QcYFGRdk7D+Kxycen3jsP9LISSMnjew/GsprKK+h7D9jy40tN7bsP6zMrzK/yuw/9M3RN0ff7D88z/M8z/PsP4XQFUJXCO0/ztE3R98c7T8W01lMZzHtP17Ue1HvRe0/p9WdVnda7T/w1r9b/27tPzjY4WCHg+0/gNkDZg+Y7T/J2iVrl6ztPxLcR3Afwe0/Wt1pdafV7T+i3ot6L+rtP+vfrX+3/u0/NOHPhD8T7j984vGJxyfuP8TjE49PPO4/DeU1lNdQ7j9W5leZX2XuP57neZ7nee4/5uibo2+O7j8v6r2o96LuP3jr361/t+4/wOwBswfM7j8I7iO4j+DuP1HvRb0X9e4/mvBnwp8J7z/i8YnHJx7vPyrzq8yvMu8/c/TN0TdH7z+89e/Wv1vvPwT3EdxHcO8/TPgz4c+E7z+V+VXmV5nvP976d+vfre8/JvyZ8GfC7z9u/bv179bvP7f+3fp36+8/AAAAAAAA8D8=\"},\"shape\":[400],\"dtype\":\"float64\",\"order\":\"little\"}],[\"y\",{\"type\":\"ndarray\",\"array\":{\"type\":\"bytes\",\"data\":\"AAAAAAAAAAAAAAAAAAAAAIrTwhIEAAAAQp09wAWMjwQYPQzgPCuLCLFN9LUHf6IL2Y5HB+b9Jw6+1JeEowBKEJdEikd51SIS9ZcwLGlNwxNl3kSaoVs3Fd0BZ3KoC4gWH4bx9tmTuxex6XoepZ/VGB/PB6VqL9sZmpqeAOAVzxobMju1ceWxGz3ELQjnjIYctMVMFl2oUB2JZ/gbtJEOHlym8PsQS8Iem+zYGx7GbR9euHDHtBMRIB7kpNMIfKwgOjlQvCS/QSEuqpIfzurQIXTqXValMlkihvVTwmzf3SJJF1qEz6hcI74wTgJGldYjNJnE9befTST7gvx9RV7AJGgnrOwNzy4l3oFWpQrwmCWi5xjrhYQBJrsxpvLAiGUm6fevtUBWxyZMYgSHK3MmJ3OwVfbuSoMn6MQXAqHL3SdYERkFMck0KJ/lj6/6VIooi8vPT6Bt3ihuyhfsLRswKa/Vq9rEXX8pL7GJo8syzCkU3k2GpXwXKh+DifcRL2IqsU/NapVAqirOeIrfPbjxKtxQNDiDbTYrJlmWSceueit75k6sIum9K4PSKePCqf8rL+eJW2O4Pyyffu48KCJ+LPxZ0Y5VMrssgO19Rchc9yytA/iWOiIzLSeK/7zO7W0tx6y4pidjpi1mLAPvKQrgLSP1FDt7CxYuJbN6v1wZTS7qMzj16XeCLl923D8JkrYu/cDYnWyW6i6CP/nuyzgeL30+Q8Z5l1Avt6CMDHCdgS+vLW6LoBmyL2tgqTNgBOIvdsEkcJRjETCBlKYHOElAMGF23dp9oG0wvRn34hoxmjCHLCcxmoXGMKV2/nqg2fIwSOUS7u67HjFQP+aKnGtIMQxoeWaM7HIxD2PBbTufnDFqPccHoCLFMZlXQbn1fu4xvr201M2BFTI0LiBkb6s9MqsEGsW0BmQye7sv2qZ3ijLqdQR5rCGxMtjqyHqPu9UyslBNUU4H+zLGHF1gFXwgMwFI9mvjukMzeWq0lvstZzMYfxYazr2KM9WRiAtGTq4zbBNYqR3g0DNNJcxbo3nyM7linv3n4xM0DcWLuZIQNTQA3ZYQtfNVNAA3dJOVhHY0p6Fskju+ljRCs0oUsJ+2NLb2xJPwK9Y0d5ADNZhp9TQ3iSdJS2IUNdt9/zb2ITM1ZNWTwPO1UTXNS0cALyxwNUoSXjanJI01INrQFjbqqTVI891ugr/GNZdTDf9Kt+M1fieNRzHgADaS3nGLT4kcNpioJApL1jc2ShHE1RytUzaX1vvlfA1wNgSP/E6b5Ik2/OejurilpDaidLkTIEjANlGGN6fYZdk2hvhhAI2Y8zab5KbxGOsNN1y2HajMmSY3RW3GTr3lQDfue4QFuQJZN61Yx9kqU3I3WktN1aCWijda6KWdHhqjN81qymRDL7s3Gwa8Y0Ip0zcB3Zd6qsLqN9avQhztgwI4do9klrFjGTh4OaCtnkAxOEuuSeJsPUc4NVV3RYoIXzh+oaQr2op0OOez9S8y94o45B9KifmMoTgV7NjFRqi2ONKsrePvAs04nRnqcqBs4jj1gYPOKDf3OBkfon7YBQ05CNFBWDAAIjm4QH7sbig2OS5/EHOUEUs56Z8rdv1oYDksy9Dm779zObCIx9VgmIc5ji3jDVT8mzkuHPhsVnqwOerXmV5hRMM5ZXR5nvle1jl7jolVmsvpORgScCxXiv05w1gMudbMEDo2gdNZLvsiOgCZmYukTTU6MvgmzLTARzqwtpn54U9aOi5K58279Ww6FQQQSuyrfzpC9gCWpjWROtXzfZoClqI6huB2V9byszrpyzAlvEfFOjBXRJA0kNY6gRQ1McDH5zrL52CF+un4Oun7cOO08gk70Zf7nhDeGjtGVCdzl6grO5gUvlpRTzw7a3+EFNfPTDtEU5uzYChdOw1teMPPV2071Auas7RdfTuzjTRuTzqNO9dNaiyL7pw7zb/Wx/V7rDuR0ULysuS7O81wcN9rK8s7mCupCDxT2jvAxy/Fm1/pO+cIaH1JVPg7GlB+QDI1BzzHjgd+WgYWPN5KUZfHyyQ8Jp4R6GqJMzxC8zjNDkNCPOfbVRRG/FA8sjykP7xwXzwEM53bp/RsPFoi0SOYiXo8zpXcLRs0iDxgAd9K//eVPFNyuG9T2KM8Sof1KWzXsTyKBE/42O2/PA7s89yjb8w8yIymNQY12Tw9wrn7vz3mPG2Yxb2miPM8BFisXMkTAT3uuWAyJ7kNPcwn/Pzgvxk97DstxNU0Jj1dF5K8lhAzPc2g769mS0A9Z/LYJNO6Sz0W6Je8mn1XPW3UT2bYz2M9tNqz1s+icD04fo0U/tB7PQJahnV7J4c9yw5WP6Uwkz3bBWjx/qufPRScrCEcBqo92Eub8aJKtT2UzJv8cVjBPc3X4NU0JMw9g9cwIlK71j2/xoIaz0jiPcoMQktiSu09gO19Rchc9z0gH2GZXY4CPiha75KqWg0+2/LXAwgfFz67dLdx2CIiPmDnOR9zVSw+gcOYhY0KNj7nZSeaTRNBPgcM42sHWUo+BxXVN3o+VD6kkZfGHvtePvrLZw7om2c+XO9XhebqcT5HMnXGzBV7PnyFgaxWY4Q+9y5lioWRjj7hniogmdKWPnR0XQBZ+KA+YhpWdUYiqT4qPq/6YomyPkW74JsWO7s+Wj57dXzrwz4NwSjTVgbNPoNGsJM2D9U+lJzjWnZv3j61ZyTOLuflPgsAnTiFZe8+DLW8jg5p9j6nXbogf9z/PgDSNWhZjgY/sVlF/T/ODz/9vwXx81QWP1UQzH0HOx8/Yc6QSVS/JT+ltvCtfCkuP6E+pm1W1DQ/m5rORzCmPD+7kv5quJ5DP5rcl3qqwko/GtEkFkwsUj+kyvpBHZRYPyH5TJHzjGA/bOqJhN0xZj88QBXRBaNtP0O9y9fJs3M/lXDTMEUVej9c/bmKxTCBP3DDMfCZj4Y/c8+TedN6jT+lkYJU9yyTP2rJMYYL1pg/0nQOhiYDoD+G7au0+o2kP6FKemJNRKo/dJORhg+1sD9ZZzl0JSi1PwMOAeQhq7o/MbSFuve6wD/j8uPOx+TEP1rFWuF7+Mk/PgjVrHYQ0D8LI2tFVMfTP7OUqjUUPNg/CArZVviM3T/Re+8Cn+3hP1n7M/F8peU/BZZqNWkB6j9IyprbyRXvP6aGSLjle/I/Z7TcE37e9T8SmKFMZr35PzzqJIP6Iv4/AFscxm+MAUC2scBWzVMEQLGxYswOawdAau2wcqTUCkBMaPgbtJEOQJMXfhDrUBFA2qpN7mmBE0BHk/NRNdgVQHknfUwLUhhAwj3Ge3TqGkCjJ/LJtpsdQOsugMtoLyBAUjbnesKVIUD65JmvMvwiQKKrgWt+XSRAnUM+/PyzJUAmnKhetvkmQH0D2QyIKChA2mGgYE86KUB4NJhPGCkqQErKBOZO7ypAjWHKkfGHK0CBldwNwu4rQNrYqJdyICxANbTYEM0aLEBlTJfM0dwrQN4gRgrMZitA2YyXh1q6KkChVcYda9opQISoiwIpyyhApeo7/N2RJ0ACTIuVxzQmQFXd4BPiuiRAwGtSlasrI0Cy8o084o4hQKEpMEyB2B9AkKKsBniWHEAv91JBjmUZQBn0Y5k6UhZAqsxLKS9nE0CXMYk7HK0QQClH5m0aVQxA8xQ1McDHB0Dvg4NYrLYDQPL1CHN2IgBAyiY3mHwQ+j8csXnKdMT0PwYiVBmuUPA/Q7kHp21D6T+zmm9SKkTjPzTPE6XB7dw/VhJ/L1Be1T+O6XUVHQrPP1d6hP/ZJsY/2dY3qcELvz/Y7dngjFi1P19Wts3qxqw/NOYXkBYAoz/zU5RIvo2YP8chQ7o0BY8/5g8GVNUigz9KxQlqoQd3P8KX3QFRAGs/Zrg/RaTMXj8+zpExhRBRP/rm3L0WV0I/oJRXhKEVMz9VzxJL3TAjP0yGLO8NnBI/oFYorb9cAT+E3kLCTRfvPkDXsuH5o9o+KsmLfhDHxT7CxEAK1uywPoCL+z416pg+UzCNBLdKgT4zzil6soVmPuLCZSSMX0s+Br6QVjTZLj62GuQBSAAQPuiSiAtGTu498nslUa7xyT2C5S8HmtmjPfdXEYZ+yHo9zHxjLQpbTz0SwulDvjwfPYe5oo3V2+k8vmtSuT9FsTwRxAeZpe5xPDr9uL3Tlys8YFs03OeJ3TsAQY2E2DKEO4+u02bYTh87lEa9m7UTpzrtJ3MBdacYOoR14KK+NWg5QmdWD+2bgjgYl9UDHFk6NwMN7BpfaAM1AAAAAAAAAAA=\"},\"shape\":[400],\"dtype\":\"float64\",\"order\":\"little\"}]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1048\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1049\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1044\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"orange\",\"line_width\":2}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Line\",\"id\":\"p1045\",\"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\":\"p1046\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"y\"},\"line_color\":\"orange\",\"line_alpha\":0.2,\"line_width\":2}}}}],\"toolbar\":{\"type\":\"object\",\"name\":\"Toolbar\",\"id\":\"p1009\",\"attributes\":{\"tools\":[{\"type\":\"object\",\"name\":\"PanTool\",\"id\":\"p1023\"},{\"type\":\"object\",\"name\":\"WheelZoomTool\",\"id\":\"p1024\",\"attributes\":{\"renderers\":\"auto\"}},{\"type\":\"object\",\"name\":\"BoxZoomTool\",\"id\":\"p1025\",\"attributes\":{\"overlay\":{\"type\":\"object\",\"name\":\"BoxAnnotation\",\"id\":\"p1026\",\"attributes\":{\"syncable\":false,\"level\":\"overlay\",\"visible\":false,\"left_units\":\"canvas\",\"right_units\":\"canvas\",\"top_units\":\"canvas\",\"bottom_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\":\"p1027\"},{\"type\":\"object\",\"name\":\"ResetTool\",\"id\":\"p1028\"},{\"type\":\"object\",\"name\":\"HelpTool\",\"id\":\"p1029\"}]}},\"left\":[{\"type\":\"object\",\"name\":\"LinearAxis\",\"id\":\"p1018\",\"attributes\":{\"ticker\":{\"type\":\"object\",\"name\":\"BasicTicker\",\"id\":\"p1019\",\"attributes\":{\"mantissas\":[1,2,5]}},\"formatter\":{\"type\":\"object\",\"name\":\"BasicTickFormatter\",\"id\":\"p1020\"},\"axis_label\":\"g(\\u03b8 | n, N)\",\"major_label_policy\":{\"type\":\"object\",\"name\":\"AllLabels\",\"id\":\"p1021\"}}}],\"below\":[{\"type\":\"object\",\"name\":\"LinearAxis\",\"id\":\"p1013\",\"attributes\":{\"ticker\":{\"type\":\"object\",\"name\":\"BasicTicker\",\"id\":\"p1014\",\"attributes\":{\"mantissas\":[1,2,5]}},\"formatter\":{\"type\":\"object\",\"name\":\"BasicTickFormatter\",\"id\":\"p1015\"},\"axis_label\":\"\\u03b8\",\"major_label_policy\":{\"type\":\"object\",\"name\":\"AllLabels\",\"id\":\"p1016\"}}}],\"center\":[{\"type\":\"object\",\"name\":\"Grid\",\"id\":\"p1017\",\"attributes\":{\"axis\":{\"id\":\"p1013\"}}},{\"type\":\"object\",\"name\":\"Grid\",\"id\":\"p1022\",\"attributes\":{\"dimension\":1,\"axis\":{\"id\":\"p1018\"}}},{\"type\":\"object\",\"name\":\"Legend\",\"id\":\"p1039\",\"attributes\":{\"location\":\"top_left\",\"items\":[{\"type\":\"object\",\"name\":\"LegendItem\",\"id\":\"p1040\",\"attributes\":{\"label\":{\"type\":\"value\",\"value\":\"old mothers\"},\"renderers\":[{\"id\":\"p1036\"}]}},{\"type\":\"object\",\"name\":\"LegendItem\",\"id\":\"p1050\",\"attributes\":{\"label\":{\"type\":\"value\",\"value\":\"young mothers\"},\"renderers\":[{\"id\":\"p1047\"}]}}]}}],\"frame_width\":300,\"frame_height\":200}}]}};\n",
       "  const render_items = [{\"docid\":\"9e75b192-2a7b-4314-8639-4b0a8ccf9b0c\",\"roots\":{\"p1001\":\"d6d966f0-a855-4894-8431-bd453806859d\"},\"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": [
    "# Observed data\n",
    "n_old = 94\n",
    "N_old = 176\n",
    "n_young = 154\n",
    "N_young = 190\n",
    "\n",
    "# Assume uniform prior\n",
    "alpha = 1\n",
    "beta = 1\n",
    "\n",
    "# Instantiate the plot\n",
    "p = bokeh.plotting.figure(\n",
    "    frame_width=300,\n",
    "    frame_height=200,\n",
    "    x_axis_label=\"θ\",\n",
    "    y_axis_label=\"g(θ | n, N)\",\n",
    "    x_range=[0, 1],\n",
    ")\n",
    "\n",
    "# Theoretical theta values\n",
    "theta = np.linspace(0, 1, 400)\n",
    "\n",
    "# Old posterior\n",
    "p.line(\n",
    "    theta,\n",
    "    st.beta.pdf(theta, n_old + alpha, N_old - n_old + beta),\n",
    "    line_width=2,\n",
    "    legend_label=\"old mothers\",\n",
    ")\n",
    "\n",
    "# Young posterior\n",
    "p.line(\n",
    "    theta,\n",
    "    st.beta.pdf(theta, n_young + alpha, N_young - n_young + beta),\n",
    "    line_width=2,\n",
    "    color=\"orange\",\n",
    "    legend_label=\"young mothers\",\n",
    ")\n",
    "\n",
    "p.legend.location = \"top_left\"\n",
    "\n",
    "bokeh.io.show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In plotting the posteriors, we see that the young mothers clearly are more likely to have their offsprings' eggs hatch than the old mothers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comments on conjugates\n",
    "\n",
    "Use of conjugate priors have advantages and disavantages. First, the advantages.\n",
    "\n",
    "1. They are conceptually pleasing in that the prior and posterior have the same form, so the update of the prior by the data via the likelihood simply updates the parameters of the distribution describing the parameter values. (See, however, the criticisms of this approach below.)\n",
    "2. They make evaluation of the posterior analytical tractable, and in fact easy. For example, you now know how to immediately write down the posterior after an experiment consisting of a series of Bernoulli trials (a Binomial likelihood) using a Beta distribution for the prior and posterior.\n",
    "\n",
    "There are disadvantages and criticisms, however.\n",
    "\n",
    "1. Finding conjugates for a given likelihood is typically very challenging for all but simple distributions. The [table on Wikipedia](https://en.wikipedia.org/wiki/Conjugate_prior) is pretty complete for known useful conjugates and it contains a paltry few distributions. For any sort of complex model, it is hopeless to find conjugates, thereby restricting their utility.\n",
    "2. Our prior information about a parameter might not actually match the form of a conjugate prior. Of example, if someone tells you a coin is slightly biased, but not in which direction, you might choose a prior that is bimodal. This cannot be described with a Beta distribution. Sivia comments on this in his book: \"While we might expect our initial understanding of the object of interest to have a bearing on the experiment we conduct, it seems strange that the choice of the prior pdf should have to wait for, and depend in detail upon, the likelihood function.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python implementation: CPython\n",
      "Python version       : 3.11.5\n",
      "IPython version      : 8.15.0\n",
      "\n",
      "numpy     : 1.26.2\n",
      "scipy     : 1.11.4\n",
      "bokeh     : 3.3.0\n",
      "jupyterlab: 4.0.10\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -v -p numpy,scipy,bokeh,jupyterlab"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.5"
  },
  "widgets": {
   "state": {
    "85aaad5108ce4f5196619e4a454fb60e": {
     "views": [
      {
       "cell_index": 27
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}