"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/javascript": [
"\n",
"(function(root) {\n",
" function now() {\n",
" return new Date();\n",
" }\n",
"\n",
" var 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",
" var JS_MIME_TYPE = 'application/javascript';\n",
" var HTML_MIME_TYPE = 'text/html';\n",
" var EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n",
" var CLASS_NAME = 'output_bokeh rendered_html';\n",
"\n",
" /**\n",
" * Render data to the DOM node\n",
" */\n",
" function render(props, node) {\n",
" var 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",
" var cell = handle.cell;\n",
"\n",
" var id = cell.output_area._bokeh_element_id;\n",
" var 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",
" var cmd = \"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, {\n",
" iopub: {\n",
" output: function(msg) {\n",
" var 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",
" var cmd = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n",
" cell.notebook.kernel.execute(cmd);\n",
" }\n",
" }\n",
"\n",
" /**\n",
" * Handle when a new output is added\n",
" */\n",
" function handleAddOutput(event, handle) {\n",
" var output_area = handle.output_area;\n",
" var output = handle.output;\n",
"\n",
" // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n",
" if ((output.output_type != \"display_data\") || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n",
" return\n",
" }\n",
"\n",
" var 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",
" var bk_div = document.createElement(\"div\");\n",
" bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n",
" var script_attrs = bk_div.children[0].attributes;\n",
" for (var i = 0; i < script_attrs.length; i++) {\n",
" toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\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",
" var 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",
" var 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",
" var events = require('base/js/events');\n",
" var 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",
"\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",
" var NB_LOAD_WARNING = {'data': {'text/html':\n",
" \"
\\n\"+\n",
" \"
\\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",
" \"
\\n\"+\n",
" \"
\\n\"+\n",
" \"
re-rerun `output_notebook()` to attempt to load from CDN again, or
\"}};\n",
"\n",
" function display_loaded() {\n",
" var el = document.getElementById(\"1001\");\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",
"\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() {\n",
" console.error(\"failed to load \" + url);\n",
" }\n",
"\n",
" for (var i = 0; i < css_urls.length; i++) {\n",
" var url = css_urls[i];\n",
" const element = document.createElement(\"link\");\n",
" element.onload = on_load;\n",
" element.onerror = on_error;\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 (var i = 0; i < js_urls.length; i++) {\n",
" var url = js_urls[i];\n",
" var element = document.createElement('script');\n",
" element.onload = on_load;\n",
" element.onerror = on_error;\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",
" };var element = document.getElementById(\"1001\");\n",
" if (element == null) {\n",
" console.error(\"Bokeh: ERROR: autoload.js configured with elementid '1001' but no matching script tag was found. \")\n",
" return false;\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",
" \n",
" var js_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-1.4.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-1.4.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-tables-1.4.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-gl-1.4.0.min.js\"];\n",
" var css_urls = [];\n",
" \n",
"\n",
" var inline_js = [\n",
" function(Bokeh) {\n",
" Bokeh.set_log_level(\"info\");\n",
" },\n",
" function(Bokeh) {\n",
" \n",
" \n",
" }\n",
" ];\n",
"\n",
" function run_inline_js() {\n",
" \n",
" if (root.Bokeh !== undefined || force === true) {\n",
" \n",
" for (var 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",
" var cell = $(document.getElementById(\"1001\")).parents('.cell').data().cell;\n",
" cell.output_area.append_execute_result(NB_LOAD_WARNING)\n",
" }\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": "\n(function(root) {\n function now() {\n return new Date();\n }\n\n var 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\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 var NB_LOAD_WARNING = {'data': {'text/html':\n \"
\\n\"+\n \"
\\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 \"
\\n\"+\n \"
\\n\"+\n \"
re-rerun `output_notebook()` to attempt to load from CDN again, or
\"}};\n\n function display_loaded() {\n var el = document.getElementById(\"1001\");\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\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() {\n console.error(\"failed to load \" + url);\n }\n\n for (var i = 0; i < css_urls.length; i++) {\n var url = css_urls[i];\n const element = document.createElement(\"link\");\n element.onload = on_load;\n element.onerror = on_error;\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 (var i = 0; i < js_urls.length; i++) {\n var url = js_urls[i];\n var element = document.createElement('script');\n element.onload = on_load;\n element.onerror = on_error;\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 };var element = document.getElementById(\"1001\");\n if (element == null) {\n console.error(\"Bokeh: ERROR: autoload.js configured with elementid '1001' but no matching script tag was found. \")\n return false;\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 \n var js_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-1.4.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-1.4.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-tables-1.4.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-gl-1.4.0.min.js\"];\n var css_urls = [];\n \n\n var inline_js = [\n function(Bokeh) {\n Bokeh.set_log_level(\"info\");\n },\n function(Bokeh) {\n \n \n }\n ];\n\n function run_inline_js() {\n \n if (root.Bokeh !== undefined || force === true) {\n \n for (var 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 var cell = $(document.getElementById(\"1001\")).parents('.cell').data().cell;\n cell.output_area.append_execute_result(NB_LOAD_WARNING)\n }\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": {
"application/javascript": [
"\n",
"(function(root) {\n",
" function now() {\n",
" return new Date();\n",
" }\n",
"\n",
" var 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",
" var JS_MIME_TYPE = 'application/javascript';\n",
" var HTML_MIME_TYPE = 'text/html';\n",
" var EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n",
" var CLASS_NAME = 'output_bokeh rendered_html';\n",
"\n",
" /**\n",
" * Render data to the DOM node\n",
" */\n",
" function render(props, node) {\n",
" var 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",
" var cell = handle.cell;\n",
"\n",
" var id = cell.output_area._bokeh_element_id;\n",
" var 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",
" var cmd = \"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, {\n",
" iopub: {\n",
" output: function(msg) {\n",
" var 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",
" var cmd = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n",
" cell.notebook.kernel.execute(cmd);\n",
" }\n",
" }\n",
"\n",
" /**\n",
" * Handle when a new output is added\n",
" */\n",
" function handleAddOutput(event, handle) {\n",
" var output_area = handle.output_area;\n",
" var output = handle.output;\n",
"\n",
" // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n",
" if ((output.output_type != \"display_data\") || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n",
" return\n",
" }\n",
"\n",
" var 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",
" var bk_div = document.createElement(\"div\");\n",
" bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n",
" var script_attrs = bk_div.children[0].attributes;\n",
" for (var i = 0; i < script_attrs.length; i++) {\n",
" toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\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",
" var 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",
" var 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",
" var events = require('base/js/events');\n",
" var 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",
"\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",
" var NB_LOAD_WARNING = {'data': {'text/html':\n",
" \"
\\n\"+\n",
" \"
\\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",
" \"
\\n\"+\n",
" \"
\\n\"+\n",
" \"
re-rerun `output_notebook()` to attempt to load from CDN again, or
\"}};\n",
"\n",
" function display_loaded() {\n",
" var el = document.getElementById(null);\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",
"\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() {\n",
" console.error(\"failed to load \" + url);\n",
" }\n",
"\n",
" for (var i = 0; i < css_urls.length; i++) {\n",
" var url = css_urls[i];\n",
" const element = document.createElement(\"link\");\n",
" element.onload = on_load;\n",
" element.onerror = on_error;\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 (var i = 0; i < js_urls.length; i++) {\n",
" var url = js_urls[i];\n",
" var element = document.createElement('script');\n",
" element.onload = on_load;\n",
" element.onerror = on_error;\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",
" \n",
" var js_urls = [];\n",
" var css_urls = [];\n",
" \n",
"\n",
" var inline_js = [\n",
" function(Bokeh) {\n",
" /* BEGIN bokeh.min.js */\n",
" /*!\n",
" * Copyright (c) 2012 - 2019, Anaconda, Inc., and Bokeh Contributors\n",
" * All rights reserved.\n",
" * \n",
" * Redistribution and use in source and binary forms, with or without modification,\n",
" * are permitted provided that the following conditions are met:\n",
" * \n",
" * Redistributions of source code must retain the above copyright notice,\n",
" * this list of conditions and the following disclaimer.\n",
" * \n",
" * Redistributions in binary form must reproduce the above copyright notice,\n",
" * this list of conditions and the following disclaimer in the documentation\n",
" * and/or other materials provided with the distribution.\n",
" * \n",
" * Neither the name of Anaconda nor the names of any contributors\n",
" * may be used to endorse or promote products derived from this software\n",
" * without specific prior written permission.\n",
" * \n",
" * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \"AS IS\"\n",
" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE\n",
" * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE\n",
" * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE\n",
" * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR\n",
" * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF\n",
" * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS\n",
" * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN\n",
" * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)\n",
" * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF\n",
" * THE POSSIBILITY OF SUCH DAMAGE.\n",
" */\n",
" (function(root, factory) {\n",
" root[\"Bokeh\"] = factory();\n",
" })(this, function() {\n",
" var define;\n",
" var parent_require = typeof require === \"function\" && require\n",
" return (function(modules, entry, aliases, externals) {\n",
" if (aliases === undefined) aliases = {};\n",
" if (externals === undefined) externals = {};\n",
"\n",
" var cache = {};\n",
"\n",
" var normalize = function(name) {\n",
" if (typeof name === \"number\")\n",
" return name;\n",
"\n",
" if (name === \"bokehjs\")\n",
" return entry;\n",
"\n",
" var prefix = \"@bokehjs/\"\n",
" if (name.slice(0, prefix.length) === prefix)\n",
" name = name.slice(prefix.length)\n",
"\n",
" var alias = aliases[name]\n",
" if (alias != null)\n",
" return alias;\n",
"\n",
" var trailing = name.length > 0 && name[name.lenght-1] === \"/\";\n",
" var index = aliases[name + (trailing ? \"\" : \"/\") + \"index\"];\n",
" if (index != null)\n",
" return index;\n",
"\n",
" return name;\n",
" }\n",
"\n",
" var require = function(name) {\n",
" var mod = cache[name];\n",
" if (!mod) {\n",
" var id = normalize(name);\n",
"\n",
" mod = cache[id];\n",
" if (!mod) {\n",
" if (!modules[id]) {\n",
" if (parent_require && externals[id]) {\n",
" try {\n",
" mod = {exports: parent_require(id)};\n",
" cache[id] = cache[name] = mod;\n",
" return mod.exports;\n",
" } catch (e) {}\n",
" }\n",
"\n",
" var err = new Error(\"Cannot find module '\" + name + \"'\");\n",
" err.code = 'MODULE_NOT_FOUND';\n",
" throw err;\n",
" }\n",
"\n",
" mod = {exports: {}};\n",
" cache[id] = cache[name] = mod;\n",
" modules[id].call(mod.exports, require, mod, mod.exports);\n",
" } else\n",
" cache[name] = mod;\n",
" }\n",
"\n",
" return mod.exports;\n",
" }\n",
"\n",
" var main = require(entry);\n",
" main.require = require;\n",
"\n",
" main.register_plugin = function(plugin_modules, plugin_entry, plugin_aliases, plugin_externals) {\n",
" if (plugin_aliases === undefined) plugin_aliases = {};\n",
" if (plugin_externals === undefined) plugin_externals = {};\n",
"\n",
" for (var name in plugin_modules) {\n",
" modules[name] = plugin_modules[name];\n",
" }\n",
"\n",
" for (var name in plugin_aliases) {\n",
" aliases[name] = plugin_aliases[name];\n",
" }\n",
"\n",
" for (var name in plugin_externals) {\n",
" externals[name] = plugin_externals[name];\n",
" }\n",
"\n",
" var plugin = require(plugin_entry);\n",
"\n",
" for (var name in plugin) {\n",
" main[name] = plugin[name];\n",
" }\n",
"\n",
" return plugin;\n",
" }\n",
"\n",
" return main;\n",
" })\n",
" ([\n",
" function _(n,o,r){n(1),function(n){for(var o in n)r.hasOwnProperty(o)||(r[o]=n[o])}(n(102))},\n",
" function _(n,c,f){n(2),n(11),n(14),n(21),n(49),n(52),n(87),n(94),n(100)},\n",
" function _(e,n,a){e(3)()||Object.defineProperty(Object,\"assign\",{value:e(4),configurable:!0,enumerable:!1,writable:!0})},\n",
" function _(r,t,o){t.exports=function(){var r,t=Object.assign;return\"function\"==typeof t&&(t(r={foo:\"raz\"},{bar:\"dwa\"},{trzy:\"trzy\"}),r.foo+r.bar+r.trzy===\"razdwatrzy\")}},\n",
" function _(t,r,n){var o=t(5),c=t(10),a=Math.max;r.exports=function(t,r){var n,f,h,i=a(arguments.length,2);for(t=Object(c(t)),h=function(o){try{t[o]=r[o]}catch(t){n||(n=t)}},f=1;f= 0\");if(!isFinite(r))throw new RangeError(\"Count must be < ∞\");for(n=\"\";r;)r%2&&(n+=t),r>1&&(t+=t),r>>=1;return n}},\n",
" function _(t,i,n){var r=t(18),a=Math.abs,o=Math.floor;i.exports=function(t){return isNaN(t)?0:0!==(t=Number(t))&&isFinite(t)?r(t)*o(a(t)):t}},\n",
" function _(n,t,i){t.exports=n(19)()?Math.sign:n(20)},\n",
" function _(n,t,o){t.exports=function(){var n=Math.sign;return\"function\"==typeof n&&(1===n(10)&&-1===n(-20))}},\n",
" function _(n,r,t){r.exports=function(n){return n=Number(n),isNaN(n)||0===n?n:n>0?1:-1}},\n",
" function _(e,r,a){e(22)()||Object.defineProperty(Array,\"from\",{value:e(23),configurable:!0,enumerable:!1,writable:!0})},\n",
" function _(n,o,r){o.exports=function(){var n,o,r=Array.from;return\"function\"==typeof r&&(o=r(n=[\"raz\",\"dwa\"]),Boolean(o&&o!==n&&\"dwa\"===o[1]))}},\n",
" function _(e,l,r){var n=e(24).iterator,t=e(44),a=e(45),i=e(46),u=e(47),o=e(10),f=e(8),c=e(48),v=Array.isArray,h=Function.prototype.call,y={configurable:!0,enumerable:!0,writable:!0,value:null},s=Object.defineProperty;l.exports=function(e){var l,r,A,g,p,w,b,d,x,j,O=arguments[1],m=arguments[2];if(e=Object(o(e)),f(O)&&u(O),this&&this!==Array&&a(this))l=this;else{if(!O){if(t(e))return 1!==(p=e.length)?Array.apply(null,e):((g=new Array(1))[0]=e[0],g);if(v(e)){for(g=new Array(p=e.length),r=0;r
\"],_default:[0,\"\",\"\"]};function ve(e,t){var n;return n=void 0!==e.getElementsByTagName?e.getElementsByTagName(t||\"*\"):void 0!==e.querySelectorAll?e.querySelectorAll(t||\"*\"):[],void 0===t||t&&N(e,t)?b.merge([e],n):n}function ye(e,t){for(var n=0,r=e.length;n-1)i&&i.push(o);else if(l=ie(o),a=ve(f.appendChild(o),\"script\"),l&&ye(a),n)for(c=0;o=a[c++];)he.test(o.type||\"\")&&n.push(o);return f}me=r.createDocumentFragment().appendChild(r.createElement(\"div\")),(xe=r.createElement(\"input\")).setAttribute(\"type\",\"radio\"),xe.setAttribute(\"checked\",\"checked\"),xe.setAttribute(\"name\",\"t\"),me.appendChild(xe),h.checkClone=me.cloneNode(!0).cloneNode(!0).lastChild.checked,me.innerHTML=\"\",h.noCloneChecked=!!me.cloneNode(!0).lastChild.defaultValue;var Te=/^key/,Ce=/^(?:mouse|pointer|contextmenu|drag|drop)|click/,Ee=/^([^.]*)(?:\\.(.+)|)/;function ke(){return!0}function Se(){return!1}function Ne(e,t){return e===function(){try{return r.activeElement}catch(e){}}()==(\"focus\"===t)}function Ae(e,t,n,r,i,o){var a,s;if(\"object\"==typeof t){for(s in\"string\"!=typeof n&&(r=r||n,n=void 0),t)Ae(e,s,n,r,t[s],o);return e}if(null==r&&null==i?(i=n,r=n=void 0):null==i&&(\"string\"==typeof n?(i=r,r=void 0):(i=r,r=n,n=void 0)),!1===i)i=Se;else if(!i)return e;return 1===o&&(a=i,(i=function(e){return b().off(e),a.apply(this,arguments)}).guid=a.guid||(a.guid=b.guid++)),e.each(function(){b.event.add(this,t,i,r,n)})}function De(e,t,n){n?(Y.set(e,t,!1),b.event.add(e,t,{namespace:!1,handler:function(e){var r,i,a=Y.get(this,t);if(1&e.isTrigger&&this[t]){if(a.length)(b.event.special[t]||{}).delegateType&&e.stopPropagation();else if(a=o.call(arguments),Y.set(this,t,a),r=n(this,t),this[t](),a!==(i=Y.get(this,t))||r?Y.set(this,t,!1):i={},a!==i)return e.stopImmediatePropagation(),e.preventDefault(),i.value}else a.length&&(Y.set(this,t,{value:b.event.trigger(b.extend(a[0],b.Event.prototype),a.slice(1),this)}),e.stopImmediatePropagation())}})):void 0===Y.get(e,t)&&b.event.add(e,t,ke)}b.event={global:{},add:function(e,t,n,r,i){var o,a,s,u,l,c,f,p,d,h,g,v=Y.get(e);if(v)for(n.handler&&(n=(o=n).handler,i=o.selector),i&&b.find.matchesSelector(re,i),n.guid||(n.guid=b.guid++),(u=v.events)||(u=v.events={}),(a=v.handle)||(a=v.handle=function(t){return void 0!==b&&b.event.triggered!==t.type?b.event.dispatch.apply(e,arguments):void 0}),l=(t=(t||\"\").match(P)||[\"\"]).length;l--;)d=g=(s=Ee.exec(t[l])||[])[1],h=(s[2]||\"\").split(\".\").sort(),d&&(f=b.event.special[d]||{},d=(i?f.delegateType:f.bindType)||d,f=b.event.special[d]||{},c=b.extend({type:d,origType:g,data:r,handler:n,guid:n.guid,selector:i,needsContext:i&&b.expr.match.needsContext.test(i),namespace:h.join(\".\")},o),(p=u[d])||((p=u[d]=[]).delegateCount=0,f.setup&&!1!==f.setup.call(e,r,h,a)||e.addEventListener&&e.addEventListener(d,a)),f.add&&(f.add.call(e,c),c.handler.guid||(c.handler.guid=n.guid)),i?p.splice(p.delegateCount++,0,c):p.push(c),b.event.global[d]=!0)},remove:function(e,t,n,r,i){var o,a,s,u,l,c,f,p,d,h,g,v=Y.hasData(e)&&Y.get(e);if(v&&(u=v.events)){for(l=(t=(t||\"\").match(P)||[\"\"]).length;l--;)if(d=g=(s=Ee.exec(t[l])||[])[1],h=(s[2]||\"\").split(\".\").sort(),d){for(f=b.event.special[d]||{},p=u[d=(r?f.delegateType:f.bindType)||d]||[],s=s[2]&&new RegExp(\"(^|\\\\.)\"+h.join(\"\\\\.(?:.*\\\\.|)\")+\"(\\\\.|$)\"),a=o=p.length;o--;)c=p[o],!i&&g!==c.origType||n&&n.guid!==c.guid||s&&!s.test(c.namespace)||r&&r!==c.selector&&(\"**\"!==r||!c.selector)||(p.splice(o,1),c.selector&&p.delegateCount--,f.remove&&f.remove.call(e,c));a&&!p.length&&(f.teardown&&!1!==f.teardown.call(e,h,v.handle)||b.removeEvent(e,d,v.handle),delete u[d])}else for(d in u)b.event.remove(e,d+t[l],n,r,!0);b.isEmptyObject(u)&&Y.remove(e,\"handle events\")}},dispatch:function(e){var t,n,r,i,o,a,s=b.event.fix(e),u=new Array(arguments.length),l=(Y.get(this,\"events\")||{})[s.type]||[],c=b.event.special[s.type]||{};for(u[0]=s,t=1;t=1))for(;l!==this;l=l.parentNode||this)if(1===l.nodeType&&(\"click\"!==e.type||!0!==l.disabled)){for(o=[],a={},n=0;n-1:b.find(i,this,null,[l]).length),a[i]&&o.push(r);o.length&&s.push({elem:l,handlers:o})}return l=this,u\\x20\\t\\r\\n\\f]*)[^>]*)\\/>/gi,qe=/"
],
"text/plain": [
":Scatter [Droplet Diameter (um)] (Spindle Length (um))"
]
},
"execution_count": 2,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1004"
}
},
"output_type": "execute_result"
}
],
"source": [
"df = pd.read_csv('../data/good_invitro_droplet_data.csv', comment='#')\n",
"\n",
"scatter = hv.Scatter(\n",
" data=df,\n",
" kdims=['Droplet Diameter (um)'],\n",
" vdims=['Spindle Length (um)']\n",
")\n",
"\n",
"scatter"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is also useful to extract the data sets as Numpy arrays for speed in the MLE and bootstrap calculations."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"d = df['Droplet Diameter (um)'].values\n",
"ell = df['Spindle Length (um)'].values"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Performing MLE by direct numerical optimization\n",
"\n",
"We can take the same approach as we did with the mRNA counts from the smFISH experiments. We write a function for the log likelihood and then use numerical optimization to find the parameters.\n",
"\n",
"The joint probability density function for the data set is the product of the probability densities of Normals. Therefore, the log likelihood is the *sum* of *log* probability densities of Normals. We can code that up using the convenient functions in the `scipy.stats` module."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def theor_spindle_length(gamma, phi, d):\n",
" \"\"\"Compute spindle length using mathematical model\"\"\"\n",
" return gamma * d / np.cbrt(1 + (gamma * d / phi)**3)\n",
" \n",
" \n",
"def log_likelihood(params, d, ell):\n",
" \"\"\"Log likelihood of spindle length model.\"\"\"\n",
" gamma, phi, sigma = params\n",
" \n",
" if gamma <= 0 or gamma > 1 or phi <= 0:\n",
" return -np.inf\n",
"\n",
" mu = theor_spindle_length(gamma, phi, d)\n",
" return np.sum(st.norm.logpdf(ell, mu, sigma))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In looking at the code for the log likelihood, the convenience of the `scipy.stats` module is clear. We more or less code up the model as we would say it!\n",
"\n",
"We can now code up the MLE calculation. For our initial guesses, we can look at the data an make estimates. The initial slope of the $l$ vs. $d$ curve is about one-half, so we will guess an initial value of $\\gamma$ of 0.5. The plateau seems to be around 40 µm, so we will guess $\\phi = 40$ µm. Eyeballing the standard deviation in the data points, I would guess it is about 5 µm."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def spindle_mle(d, ell):\n",
" \"\"\"Compute MLE for parameters in spindle length model.\"\"\"\n",
" with warnings.catch_warnings():\n",
" warnings.simplefilter(\"ignore\")\n",
"\n",
" res = scipy.optimize.minimize(\n",
" fun=lambda params, d, ell: -log_likelihood(params, d, ell),\n",
" x0=np.array([0.5, 35, 5]),\n",
" args=(d, ell),\n",
" method='Powell'\n",
" )\n",
"\n",
" if res.success:\n",
" return res.x\n",
" else:\n",
" raise RuntimeError('Convergence failed with message', res.message)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's run the calculation to get MLEs for the parameters."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0.86047455, 38.23124995, 3.75342168])"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mle_params = spindle_mle(d, ell)\n",
"\n",
"mle_params"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can take a quick look as to how the theoretical curve might look with respect to the data."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":Overlay\n",
" .Scatter.I :Scatter [Droplet Diameter (um)] (Spindle Length (um))\n",
" .Curve.I :Curve [x] (y)"
]
},
"execution_count": 7,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1114"
}
},
"output_type": "execute_result"
}
],
"source": [
"# Compute theoretical curve using MLE parameters\n",
"d_theor = np.linspace(0, 250, 200)\n",
"ell_theor = theor_spindle_length(*mle_params[:-1], d_theor)\n",
"\n",
"# Make plots\n",
"scatter * hv.Curve(\n",
" data=(d_theor, ell_theor)\n",
").opts(\n",
" color='orange',\n",
" line_width=2\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of course, we are not using the full generative model here, since we are not using the parameter $\\sigma$. We will discuss in much more detail next week how to visualize the entire generative model as it relates to the measured data."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Least squares\n",
"\n",
"The likelihood for this model is\n",
"\n",
"\\begin{align}\n",
"L( \\gamma, \\phi, \\sigma; \\mathbf{l}, \\mathbf{d}) = \\prod_{i} f(l_i ; d_i, \\gamma, \\phi, \\sigma) = \\frac{1}{(2\\pi \\sigma^2)^{n/2}}\\,\\exp\\left[-\\frac{1}{2\\sigma^2}\\sum_{i}(l_i - \\mu_i)^2\\right],\n",
"\\end{align}\n",
"\n",
"with \n",
"\n",
"\\begin{align}\n",
"\\mu_i = \\frac{\\gamma d_i}{\\left(1+(\\gamma d_i/\\phi)^3\\right)^{\\frac{1}{3}}}\n",
"\\end{align}.\n",
"\n",
"Notice that the parameters $\\gamma$ and $\\phi$ only appear in the parameter $\\mu$. Therefore, to find the parameter values of $\\gamma$ and $\\phi$ that maximize the likelihood, we need only to find the values that minimize\n",
"\n",
"\\begin{align}\n",
"\\text{RSS} \\equiv \\sum_{i}(l_i - \\mu_i)^2,\n",
"\\end{align}\n",
"\n",
"a quantity known as the residual sum of squares. This is true for whatever value $\\sigma$ takes. So, assume for a moment we find the MLEs $\\gamma^*$ and $\\phi^*$, which give us a residual sum of squares of $\\text{RSS}^*$. Then, the log-likelihood is\n",
"\n",
"\\begin{align}\n",
"\\ell(\\gamma^*, \\phi^*, \\sigma) = -\\frac{n}{2}\\,\\ln 2\\pi - n\\ln\\sigma - \\frac{\\text{RSS}^*}{2\\sigma^2}.\n",
"\\end{align}\n",
"\n",
"The value of $\\sigma$ that maximizes the log-likelihood is then found by differentiating and setting to zero.\n",
"\n",
"\\begin{align}\n",
"\\frac{\\partial \\ell}{\\partial \\sigma} = -\\frac{n}{\\sigma} + \\frac{\\text{RSS}^*}{\\sigma^3} = 0.\n",
"\\end{align}\n",
"\n",
"This is solved to give $\\sigma^* = \\sqrt{\\text{RSS}^*/n}$.\n",
"\n",
"So, we can split the optimization problem in two parts. First, find the values of parameters $\\gamma$ and $\\phi$ that minimize the residual sum of squares. Then, using these values, compute the MLE for $\\sigma$ using the analytical result we have just derived.\n",
"\n",
"It turns out that this is quite powerful because the optimization problem for minimizing the residual sum of squares has nice properties. Importantly, we can use the powerful **[Levenberg-Marquardt algorithm](https://en.wikipedia.org/wiki/Levenberg–Marquardt_algorithm)** to perform the optimization. A variant of this is implemented in `scipy.optimization.least_squares()` ([docs](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html)). To use that algorithm you need to define a function for the residual; that is a function that returns the quantity $l_i - \\mu_i$. The function has call signature `f(params, *args)`, where `params` is a Numpy array containing the parameters you will use in the optimization (in our case, $\\gamma$ and $\\phi$), and `args` is a tuple of any other arguments passed to the function, which almost always include the data (in our case `(d, ell)`). We can code up that function."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def resid(params, d, ell):\n",
" \"\"\"Residual for spindle length model.\"\"\"\n",
" return ell - theor_spindle_length(*params, d)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can call `scipy.optimize.least_squares()` to perform the MLE. We can also specify bounds as a tuple containing a list of the lower bounds of the respective parameters and a list of the upper bounds of the respective parameters. We bound $\\gamma$ between zero and one, and $\\phi$ just has to be positive."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([ 0.8601068 , 38.24676911]), 3.754999871544005)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res = scipy.optimize.least_squares(\n",
" resid, np.array([0.5, 35]), args=(d, ell), bounds=([0, 0], [1, np.inf])\n",
")\n",
"\n",
"# Compute residual sum of squares from optimal params\n",
"rss_mle = np.sum(resid(res.x, d, ell)**2)\n",
"\n",
"# Compute MLE for sigma\n",
"sigma_mle = np.sqrt(rss_mle / len(d))\n",
"\n",
"# Take a look\n",
"res.x, sigma_mle"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Indeed, we got the same results as before. Let's code this procedure up into a function. Looking ahead to using `bebi103.draw_bs_reps_mle()`, we need the function to take the data set as a single argument."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def spindle_mle_lstq(data):\n",
" \"\"\"Compute MLE for parameters in spindle length model.\"\"\"\n",
" # Unpack data\n",
" d = data[:, 0]\n",
" ell = data[:, 1]\n",
" \n",
" res = scipy.optimize.least_squares(\n",
" resid, np.array([0.5, 35]), args=(d, ell), bounds=([0, 0], [1, np.inf])\n",
" )\n",
"\n",
" # Compute residual sum of squares from optimal params\n",
" rss_mle = np.sum(resid(res.x, d, ell)**2)\n",
"\n",
" # Compute MLE for sigma\n",
" sigma_mle = np.sqrt(rss_mle / len(d))\n",
" \n",
" return tuple([x for x in res.x] + [sigma_mle])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The big advantages of doing this are the robustness and speed of the Levenberg-Marquardt algorithm. Let's do a speed test."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Timing for MLE by Powell's method:\n",
"29.6 ms ± 555 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
"\n",
"Timing for MLE by least_squares:\n",
"9.17 ms ± 895 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
]
}
],
"source": [
"print(\"Timing for MLE by Powell's method:\")\n",
"%timeit spindle_mle(d, ell)\n",
"\n",
"print(\"\\nTiming for MLE by least_squares:\")\n",
"data = df[['Droplet Diameter (um)', 'Spindle Length (um)']].values\n",
"%timeit spindle_mle_lstq(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That's a substantial speed boost of nearly 5x."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Very important caveat\n",
"\n",
"We just did something you may have heard of. We minimized the sum of the squares of the residuals. Notice that this depended on the model for the residuals being Normal. **This does not in general work for all models!** We will also find that this does not in general work in a Bayesian setting with non-Uniform priors, which is almost always the case."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Confidence intervals\n",
"\n",
"We can now proceed to compute confidence intervals for our parameters. There are many ways to do this, including pairs, residual, and wild bootstrap, which will be topics for recitation. We will instead do a parametric bootstrap here. We use the generative model parametrized by the MLE to regenerate data sets. For each value of $d$ (which are fixed), we draw values of $l$ using the generative distribution). Then we perform MLE of those to get confidence intervals. To use `bebi103.draw_bs_reps_mle()`, then, we just are left to supply the bootstrap sample generating function."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def gen_spindle_data(params, d, size, rg):\n",
" \"\"\"Generate a new spindle data set.\"\"\"\n",
" mu = theor_spindle_length(*params[:-1], d)\n",
" sigma = params[-1]\n",
"\n",
" return np.vstack((d, rg.normal(mu, sigma))).transpose()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that we were careful to make sure the data set has the appropriate dimensions; the row indexes the trial. We can not compute our replicates and confidence intervals."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 1666/1666 [00:23<00:00, 70.82it/s]\n",
"100%|██████████| 1668/1668 [00:23<00:00, 70.87it/s]\n",
"100%|██████████| 1666/1666 [00:23<00:00, 70.77it/s]\n"
]
},
{
"data": {
"text/plain": [
"array([[ 0.82874058, 37.54114396, 3.55172247],\n",
" [ 0.89512945, 39.00949111, 3.95081461]])"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bs_reps = bebi103.draw_bs_reps_mle(\n",
" spindle_mle_lstq,\n",
" gen_spindle_data,\n",
" data = df[['Droplet Diameter (um)', 'Spindle Length (um)']].values,\n",
" mle_args=(),\n",
" gen_args=(d,),\n",
" size=5000,\n",
" n_jobs=3,\n",
" progress_bar=True,\n",
")\n",
"\n",
"# Compute confidence intervals\n",
"np.percentile(bs_reps, [2.5, 97.5], axis=0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The confidence intervals are given in each column for $\\gamma$, $\\phi$, and $\\sigma$, respectively, and are quite tight. This is because we have lots of data points. This does *not* mean that the data are tightly distributed about the theoretical curve, only that the MLE estimates for the parameters will not vary much if we repeated the experiment and analysis."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Computing environment"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPython 3.7.5\n",
"IPython 7.1.1\n",
"\n",
"numpy 1.17.4\n",
"scipy 1.3.1\n",
"pandas 0.24.2\n",
"tqdm 4.40.0\n",
"bokeh 1.4.0\n",
"holoviews 1.12.7\n",
"bokeh_catplot 0.1.7\n",
"bebi103 0.0.46\n",
"jupyterlab 1.2.3\n"
]
}
],
"source": [
"%load_ext watermark\n",
"\n",
"%watermark -v -p numpy,scipy,pandas,tqdm,bokeh,holoviews,bokeh_catplot,bebi103,jupyterlab"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}