{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "
" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import HTML\n", "\n", "HTML('''\n", "
''')" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import scipy.interpolate\n", "import scipy.stats " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import bokeh.plotting\n", "import bokeh.models # LinearAxis \n", "from bokeh.palettes import Category10_10 as palette" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "
\n", " \n", " Loading BokehJS ...\n", "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "\n", "(function(global) {\n", " function now() {\n", " return new Date();\n", " }\n", "\n", " var force = true;\n", "\n", " if (typeof (window._bokeh_onload_callbacks) === \"undefined\" || force === true) {\n", " window._bokeh_onload_callbacks = [];\n", " window._bokeh_is_loading = undefined;\n", " }\n", "\n", "\n", " \n", " if (typeof (window._bokeh_timeout) === \"undefined\" || force === true) {\n", " window._bokeh_timeout = Date.now() + 5000;\n", " window._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", " \"\\n\"+\n", " \"from bokeh.resources import INLINE\\n\"+\n", " \"output_notebook(resources=INLINE)\\n\"+\n", " \"\\n\"+\n", " \"
\"}};\n", "\n", " function display_loaded() {\n", " if (window.Bokeh !== undefined) {\n", " var el = document.getElementById(\"780c6343-d35d-488e-a691-c732a2e0244d\");\n", " el.textContent = \"BokehJS \" + Bokeh.version + \" successfully loaded.\";\n", " } else if (Date.now() < window._bokeh_timeout) {\n", " setTimeout(display_loaded, 100)\n", " }\n", " }\n", "\n", " function run_callbacks() {\n", " window._bokeh_onload_callbacks.forEach(function(callback) { callback() });\n", " delete window._bokeh_onload_callbacks\n", " console.info(\"Bokeh: all callbacks have finished\");\n", " }\n", "\n", " function load_libs(js_urls, callback) {\n", " window._bokeh_onload_callbacks.push(callback);\n", " if (window._bokeh_is_loading > 0) {\n", " console.log(\"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.log(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n", " window._bokeh_is_loading = js_urls.length;\n", " for (var i = 0; i < js_urls.length; i++) {\n", " var url = js_urls[i];\n", " var s = document.createElement('script');\n", " s.src = url;\n", " s.async = false;\n", " s.onreadystatechange = s.onload = function() {\n", " window._bokeh_is_loading--;\n", " if (window._bokeh_is_loading === 0) {\n", " console.log(\"Bokeh: all BokehJS libraries loaded\");\n", " run_callbacks()\n", " }\n", " };\n", " s.onerror = function() {\n", " console.warn(\"failed to load library \" + url);\n", " };\n", " console.log(\"Bokeh: injecting script tag for BokehJS library: \", url);\n", " document.getElementsByTagName(\"head\")[0].appendChild(s);\n", " }\n", " };var element = document.getElementById(\"780c6343-d35d-488e-a691-c732a2e0244d\");\n", " if (element == null) {\n", " console.log(\"Bokeh: ERROR: autoload.js configured with elementid '780c6343-d35d-488e-a691-c732a2e0244d' but no matching script tag was found. \")\n", " return false;\n", " }\n", "\n", " var js_urls = [\"https://cdn.pydata.org/bokeh/dev/bokeh-0.12.6dev4.min.js\", \"https://cdn.pydata.org/bokeh/dev/bokeh-widgets-0.12.6dev4.min.js\"];\n", "\n", " var inline_js = [\n", " function(Bokeh) {\n", " Bokeh.set_log_level(\"info\");\n", " },\n", " \n", " function(Bokeh) {\n", " \n", " },\n", " \n", " function(Bokeh) {\n", " \n", " document.getElementById(\"780c6343-d35d-488e-a691-c732a2e0244d\").textContent = \"BokehJS is loading...\";\n", " },\n", " function(Bokeh) {\n", " console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/dev/bokeh-0.12.6dev4.min.css\");\n", " Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/dev/bokeh-0.12.6dev4.min.css\");\n", " console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/dev/bokeh-widgets-0.12.6dev4.min.css\");\n", " Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/dev/bokeh-widgets-0.12.6dev4.min.css\");\n", " }\n", " ];\n", "\n", " function run_inline_js() {\n", " \n", " if ((window.Bokeh !== undefined) || (force === true)) {\n", " for (var i = 0; i < inline_js.length; i++) {\n", " inline_js[i](window.Bokeh);\n", " }if (force === true) {\n", " display_loaded();\n", " }} else if (Date.now() < window._bokeh_timeout) {\n", " setTimeout(run_inline_js, 100);\n", " } else if (!window._bokeh_failed_load) {\n", " console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n", " window._bokeh_failed_load = true;\n", " } else if (force !== true) {\n", " var cell = $(document.getElementById(\"780c6343-d35d-488e-a691-c732a2e0244d\")).parents('.cell').data().cell;\n", " cell.output_area.append_execute_result(NB_LOAD_WARNING)\n", " }\n", "\n", " }\n", "\n", " if (window._bokeh_is_loading === 0) {\n", " console.log(\"Bokeh: BokehJS loaded, going straight to plotting\");\n", " run_inline_js();\n", " } else {\n", " load_libs(js_urls, function() {\n", " console.log(\"Bokeh: BokehJS plotting callback run at\", now());\n", " run_inline_js();\n", " });\n", " }\n", "}(this));" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "bokeh.plotting.output_notebook()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question:\n", "Given two normal distributions of different mean and standard deviation, consider a measurement with a value between the two means, what is the probability for it belonging to each of the two distributions?\n", "\n", "As an example, consider two groups of runners. The first group of runners trains regularly and is in good shape. The second group of runners does not train but still enjoy running occaisionally. You measure how fast everyone can run the mile then calculate the mean and standard deviation for each group. Then you find you have the time for one person but forgot to label what group they came from. Based on their time for running the mile, what is the probability they are in the group that trains regularly?\n", "\n", "To put some numbers to this example, lets say the group that trains regularly have a mean time of 6.5 minutes with a standard deviation of 0.75 minutes and that the occasional runners have a mean time of 10 minutes with a standard devaition of 2 minutes. Additionally, lets say there are 20 people in the faster group and 80 people in the slower group.\n", "\n", "Mean times:
\n", "$ m_1 = 6.5\\,\\text{min}$
\n", "$ m_2 = 12.0\\,\\text{min}$
\n", "\n", "Standard deviations:
\n", "$ s_1 = 0.75\\,\\text{min}$
\n", "$ s_2 = 3.0\\,\\text{min}$
\n", "\n", "Group size:
\n", "$ n_1 = 20\\,\\text{people}$
\n", "$ n_2 = 80\\,\\text{people}$
\n", "\n", "## Caveat:\n", "Note that often a Gaussian distribution is a simplified view of a more complicated genuine distribution. For this situation, it is much more likely that both groups of runners would have an asymmetric distribution. We could safely expect both groups of runners, but especially the group that does not train regularly, would exhibit a longer tail for slower times. In the following analysis, we proceed under the optimized assumption that each group can be described by a single Gaussian distribution." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "m1 = 6.5\n", "m2 = 12.0\n", "\n", "s1 = .75\n", "s2 = 3.0\n", "\n", "n1 = 20\n", "n2 = 80\n", "n = n1 + n2" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x_min = 2.0\n", "x_max = 25.0\n", "x = np.linspace(x_min, x_max, 200)\n", "y1 = scipy.stats.norm.pdf(x, loc=m1, scale=s1) * n1 / n\n", "\n", "y2 = scipy.stats.norm.pdf(x, loc=m2, scale=s2) * n2 / n\n", "\n", "y = y1 + y2\n", "\n", "p = bokeh.plotting.figure(width=400, height=300, \n", " x_axis_label='Time (m)',\n", " y_axis_label='% of People')\n", "p.line(x, y, legend='combined', color=palette[0], line_width=2)\n", "p.line(x, y1, legend='fast group', color=palette[1])\n", "p.line(x, y2, legend='slow group', color=palette[2])\n", "\n", "p.legend.location = 'top_right'\n", "\n", "bokeh.plotting.show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we calculate the area under the curves, we should get the percent of people in each group." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "dx = x[1] - x[0]" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.999674635352\n" ] } ], "source": [ "a = dx * y.sum()\n", "print(a)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.199999999882\n" ] } ], "source": [ "a1 = dx * y1.sum()\n", "print(a1)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.79967463547\n" ] } ], "source": [ "a2 = dx * y2.sum()\n", "print(a2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now lets consider the runner who was not labelled. \n", "\n", "My initial approach was to compare his time, $t$ to the means and standard distributions, calculating the difference between the means as\n", "$$\n", "\\begin{align}\n", " d_1 &= | t - m_1|\\\\\n", " d_2 &= | t - m_2|\n", "\\end{align}\n", "$$ \n", "then calculate the probability of being in the fast group as\n", "$$\n", "p_\\text{fast} =\\left\\{ \\begin{array}{cl}\n", " 1 &\\mbox{ if $t < m_1$} \\\\\n", " 1 - \\dfrac{s_1 d_1}{s_1 d_1 + s_2 d_2} &\\mbox{ if $m_1 <= t < m_2$}\\\\\n", " 0 &\\mbox{ if $m_2 <= t$}\n", " \\end{array} \\right.\n", "$$\n", "\n", "This probability will have the following distribution." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "p1_v1 = np.zeros_like(x)\n", "p1_v1[x < m1] = 1\n", "\n", "d1 = np.abs(x - m1)\n", "d2 = np.abs(x - m2)\n", "\n", "mask = np.logical_and(m1 < x, x < m2)\n", "p1_v1[mask] = 1 - s1 * d1[mask] / (s1 * d1[mask] + s2 * d2[mask])" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = bokeh.plotting.figure(width=400, height=300, toolbar_location='above',\n", " x_axis_label='Time (m)',\n", " y_axis_label='% of People')\n", "p.extra_y_ranges = {'prob': bokeh.models.Range1d(start=0, end=1.02)}\n", "p.add_layout(bokeh.models.LinearAxis(y_range_name='prob', \n", " axis_label='Probability'), 'right')\n", "\n", "p.y_range = bokeh.models.Range1d(start=0, end=0.15)\n", "p.line(x, y, legend='combined', color=palette[0], line_width=2)\n", "p.line(x, y1, legend='fast group', color=palette[1])\n", "p.line(x, y2, legend='slow group', color=palette[2])\n", "p.line(x, p1_v1, legend='p fast', color=palette[3], y_range_name='prob')\n", "p.line(x, 1-p1_v1, legend='p slow', color=palette[4], y_range_name='prob')\n", "\n", "p.legend.location = 'top_right'\n", "\n", "bokeh.plotting.show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is clearly wrong as the crossover from being more likely to be in fast group is much too close to the slow group. Lets see what happens when we switch the standard deviations to use\n", "$$\n", "p_\\text{fast} =\\left\\{ \\begin{array}{cl}\n", " 1 &\\mbox{ if $t < m_1$} \\\\\n", " 1 - \\dfrac{s_2 d_1}{s_2 d_1 + s_1 d_2} &\\mbox{ if $m_1 <= t < m_2$}\\\\\n", " 0 &\\mbox{ if $m_2 <= t$}\n", " \\end{array} \\right.\n", "$$\n", "\n", "With this change, the probability will have the following distribution." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [], "source": [ "p1_v1[mask] = 1 - s2 * d1[mask] / (s2 * d1[mask] + s1 * d2[mask])\n", "p2_v1 = 1 - p1_v1" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = bokeh.plotting.figure(width=400, height=300, toolbar_location='above',\n", " x_axis_label='Time (m)',\n", " y_axis_label='% of People')\n", "p.extra_y_ranges = {'prob': bokeh.models.Range1d(start=0, end=1.02)}\n", "p.add_layout(bokeh.models.LinearAxis(y_range_name='prob', \n", " axis_label='Probability'), 'right')\n", "\n", "p.y_range = bokeh.models.Range1d(start=0, end=0.15)\n", "p.line(x, y, legend='combined', color=palette[0], line_width=2)\n", "p.line(x, y1, legend='fast group', color=palette[1])\n", "p.line(x, y2, legend='slow group', color=palette[2])\n", "p.line(x, p1_v1, legend='p fast', color=palette[3], y_range_name='prob')\n", "p.line(x, p2_v1, legend='p slow', color=palette[4], y_range_name='prob')\n", "\n", "p.legend.location = 'top_right'\n", "\n", "bokeh.plotting.show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This seems more reasonable but I did not take into account the relative difference in the size of the groups. Lets add a scale factor related to this difference.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_\\text{fast} =\\left\\{ \\begin{array}{cl}\n", " 1 &\\mbox{ if $t < m_1$} \\\\\n", " 1 - \\dfrac{s_2 d_1 \\frac{n_1}{n}}{s_2 d_1 \\frac{n_1}{n} + s_1 d_2 \\frac{n_2}{n}} &\\mbox{ if $m_1 <= t < m_2$}\\\\\n", " 0 &\\mbox{ if $m_2 <= t$}\n", " \\end{array} \\right.\n", "$$" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [], "source": [ "p1_v1[mask] = 1 - s2 * d1[mask] * n1 / n / (s2 * d1[mask] * n1 / n + s1 * d2[mask] * n2 / n)\n", "p2_v1 = 1 - p1_v1" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = bokeh.plotting.figure(width=400, height=300, toolbar_location='above',\n", " x_axis_label='Time (m)',\n", " y_axis_label='% of People')\n", "p.extra_y_ranges = {'prob': bokeh.models.Range1d(start=0, end=1.02)}\n", "p.add_layout(bokeh.models.LinearAxis(y_range_name='prob', \n", " axis_label='Probability'), 'right')\n", "\n", "p.y_range = bokeh.models.Range1d(start=0, end=0.15)\n", "p.line(x, y, legend='combined', color=palette[0], line_width=2)\n", "p.line(x, y1, legend='fast group', color=palette[1])\n", "p.line(x, y2, legend='slow group', color=palette[2])\n", "p.line(x, p1_v1, legend='p fast', color=palette[3], y_range_name='prob')\n", "p.line(x, p2_v1, legend='p slow', color=palette[4], y_range_name='prob')\n", "\n", "p.legend.location = 'top_right'\n", "\n", "bokeh.plotting.show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This looks like a linear switch between the two means. How does this change if we change one of the mean values? Let make $m_2=15.0$ minutes." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m2 = 15.0\n", "\n", "y1 = scipy.stats.norm.pdf(x, loc=m1, scale=s1) * n1 / n\n", "y2 = scipy.stats.norm.pdf(x, loc=m2, scale=s2) * n2 / n\n", "y = y1 + y2\n", "\n", "p1_v1 = np.zeros_like(x)\n", "p1_v1[x < m1] = 1\n", "\n", "d1 = np.abs(x - m1)\n", "d2 = np.abs(x - m2)\n", "\n", "mask = np.logical_and(m1 < x, x < m2)\n", "p1_v1[mask] = 1 - s1 * d1[mask] / (s1 * d1[mask] + s2 * d2[mask])\n", "p1_v1[mask] = 1 - s2 * d1[mask] * n1 / n / (s2 * d1[mask] * n1 / n + s1 * d2[mask] * n2 / n)\n", "p2_v1 = 1 - p1_v1" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = bokeh.plotting.figure(width=400, height=300, toolbar_location='above',\n", " x_axis_label='Time (m)',\n", " y_axis_label='% of People')\n", "p.extra_y_ranges = {'prob': bokeh.models.Range1d(start=0, end=1.02)}\n", "p.add_layout(bokeh.models.LinearAxis(y_range_name='prob', \n", " axis_label='Probability'), 'right')\n", "\n", "p.y_range = bokeh.models.Range1d(start=0, end=0.15)\n", "p.line(x, y, legend='combined', color=palette[0], line_width=2)\n", "p.line(x, y1, legend='fast group', color=palette[1])\n", "p.line(x, y2, legend='slow group', color=palette[2])\n", "p.line(x, p1_v1, legend='p fast', color=palette[3], y_range_name='prob')\n", "p.line(x, p2_v1, legend='p slow', color=palette[4], y_range_name='prob')\n", "\n", "p.legend.location = 'top_right'\n", "\n", "bokeh.plotting.show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This still looks like a linear switch between the two means. This seem a littles strange. Lets see how this changes if we change the relative size of the distributions? Let change the mean of the slow group back to the original value, $m_2 = 12.0$, and adjust the size of the groups so $n_1=40$ and $n_2=60$ people." ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m2 = 12.0\n", "n1 = 40.0\n", "n2 = 60.0\n", "n = n1 + n2\n", "\n", "y1 = scipy.stats.norm.pdf(x, loc=m1, scale=s1) * n1 / n\n", "y2 = scipy.stats.norm.pdf(x, loc=m2, scale=s2) * n2 / n\n", "y = y1 + y2\n", "\n", "p1_v1 = np.zeros_like(x)\n", "p1_v1[x < m1] = 1\n", "\n", "d1 = np.abs(x - m1)\n", "d2 = np.abs(x - m2)\n", "\n", "mask = np.logical_and(m1 < x, x < m2)\n", "p1_v1[mask] = 1 - s1 * d1[mask] / (s1 * d1[mask] + s2 * d2[mask])\n", "p1_v1[mask] = 1 - s2 * d1[mask] * n1 / n / (s2 * d1[mask] * n1 / n + s1 * d2[mask] * n2 / n)\n", "p2_v1 = 1 - p1_v1" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = bokeh.plotting.figure(width=400, height=300, toolbar_location='above',\n", " x_axis_label='Time (m)',\n", " y_axis_label='% of People')\n", "p.extra_y_ranges = {'prob': bokeh.models.Range1d(start=0, end=1.02)}\n", "p.add_layout(bokeh.models.LinearAxis(y_range_name='prob', \n", " axis_label='Probability'), 'right')\n", "\n", "p.y_range = bokeh.models.Range1d(start=0, end=0.15)\n", "p.line(x, y, legend='combined', color=palette[0], line_width=2)\n", "p.line(x, y1, legend='fast group', color=palette[1])\n", "p.line(x, y2, legend='slow group', color=palette[2])\n", "p.line(x, p1_v1, legend='p fast', color=palette[3], y_range_name='prob')\n", "p.line(x, p2_v1, legend='p slow', color=palette[4], y_range_name='prob')\n", "\n", "p.legend.location = 'top_right'\n", "\n", "bokeh.plotting.show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This shows that adjusting the verifies that the linear switch was a coincidence of our selection. I am still considering if this implementation of the relative size is appropriate. \n", "\n", "In a more general sense, while the approach considered so far provides a switch between 0 and 1 for the times that occur between the two distributions, it does not apply what is know about a normal distribution and the probability of being a certain distance from the mean. As a more concrete example, there is a finite probability, small as it may be in this case, that a runner from the slow group would run faster than the mean for the fast group. We can apply what know about a Gaussian distribution to account for this.\n", "\n", "Lets reset the group sizes to the original $n_1 = 20$ and $n_2 = 80$." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n1 = 20.0\n", "n2 = 80.0\n", "n = n1 + n2\n", "\n", "y1 = scipy.stats.norm.pdf(x, loc=m1, scale=s1) * n1 / n\n", "y2 = scipy.stats.norm.pdf(x, loc=m2, scale=s2) * n2 / n\n", "y = y1 + y2\n", "\n", "p1_v1 = np.zeros_like(x)\n", "p1_v1[x < m1] = 1\n", "\n", "d1 = np.abs(x - m1)\n", "d2 = np.abs(x - m2)\n", "\n", "mask = np.logical_and(m1 < x, x < m2)\n", "p1_v1[mask] = 1 - s1 * d1[mask] / (s1 * d1[mask] + s2 * d2[mask])\n", "p1_v1[mask] = 1 - s2 * d1[mask] * n1 / n / (s2 * d1[mask] * n1 / n + s1 * d2[mask] * n2 / n)\n", "p2_v1 = 1 - p1_v1" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = bokeh.plotting.figure(width=400, height=300, toolbar_location='above',\n", " x_axis_label='Time (m)',\n", " y_axis_label='% of People')\n", "p.extra_y_ranges = {'prob': bokeh.models.Range1d(start=0, end=1.02)}\n", "p.add_layout(bokeh.models.LinearAxis(y_range_name='prob', \n", " axis_label='Probability'), 'right')\n", "\n", "p.y_range = bokeh.models.Range1d(start=0, end=0.15)\n", "p.line(x, y, legend='combined', color=palette[0], line_width=2)\n", "p.line(x, y1, legend='fast group', color=palette[1])\n", "p.line(x, y2, legend='slow group', color=palette[2])\n", "p.line(x, p1_v1, legend='p fast', color=palette[3], y_range_name='prob')\n", "p.line(x, p2_v1, legend='p slow', color=palette[4], y_range_name='prob')\n", "\n", "p.legend.location = 'top_right'\n", "\n", "bokeh.plotting.show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now lets calculate the probability a runner would have a given time for each distribution. To do this, we will use [`scipy.stats.cdf`](https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.stats.norm.html), the [cumulative distribution function](https://en.wikipedia.org/wiki/Cumulative_distribution_function), taking the difference of probabilities between a small $\\delta t$." ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x_lower = x - dx / 2\n", "x_upper = x + dx / 2\n", "pd1 = (scipy.stats.norm.cdf(x_upper, loc=m1, scale=s1) - \n", " scipy.stats.norm.cdf(x_lower, loc=m1, scale=s1)) * n1 / n\n", "pd2 = (scipy.stats.norm.cdf(x_upper, loc=m2, scale=s2) - \n", " scipy.stats.norm.cdf(x_lower, loc=m2, scale=s2)) * n2 / n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Summing each of these probability distributions should give us the percent of the total population in each group." ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [], "source": [ "atol = 0.001 # absolute tolerance\n", "assert np.allclose(pd1.sum(), n1 / n, atol=atol)\n", "assert np.allclose(pd2.sum(), n2 / n, atol=atol)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can consider the unlabelled runner, comparing the likelihoods he is in each group. This will take the form\n", "$$\n", "p_\\text{fast} = \\dfrac{pd_1(t)}{pd_1(t) + pd_2(t)}\n", "$$\n", "where $pd_1(t)$ and $pd_2(t)$ are the probabilities a member of the fast and slow group would respectively have a time, $t$." ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [], "source": [ "p1_v2 = pd1 / (pd1 + pd2)\n", "p2_v2 = 1 - p1_v2" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = bokeh.plotting.figure(width=400, height=300, toolbar_location='above',\n", " x_axis_label='Time (m)',\n", " y_axis_label='% of People')\n", "p.extra_y_ranges = {'prob': bokeh.models.Range1d(start=0, end=1.02)}\n", "p.add_layout(bokeh.models.LinearAxis(y_range_name='prob', \n", " axis_label='Probability'), 'right')\n", "\n", "p.y_range = bokeh.models.Range1d(start=0, end=0.15)\n", "p.line(x, y, legend='combined', color=palette[0], line_width=2)\n", "p.line(x, y1, legend='fast group', color=palette[1])\n", "p.line(x, y2, legend='slow group', color=palette[2])\n", "p.line(x, p1_v2, legend='p fast', color=palette[3], y_range_name='prob')\n", "p.line(x, p2_v2, legend='p slow', color=palette[4], y_range_name='prob')\n", "\n", "p.legend.location = 'top_right'\n", "\n", "bokeh.plotting.show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While at first this seems a bit strange (the oscillatory nature of the probability of being in the slow group), on closer inspection we see this does make sense. A time close to mean values of the distributions would get grouped with that distribution. The reason the probability of the slow group goes back up for the quickest times, $t < $ 4.18 min, is that group has such a wide standard deviation that the probability decays much slower than for the narrow distribution of the fast group.\n", "\n", "While the larger standard deviation makes it more likely the fastest runner did not train, we could impose a constraint to have \"p fast\" smoothly transition to 1 at the quickest times. This would coorespond to the slow group having a more assymmetric distribution, more likely exhibiting a longer tail toward the slower times." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets compare my original probability and this corrected one." ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p_per = bokeh.plotting.figure(width=400, height=300,\n", " y_axis_label='% of People')\n", "p_per.line(x, y, legend='combined', color=palette[0], line_width=2)\n", "p_per.line(x, y1, legend='fast group', color=palette[1])\n", "p_per.line(x, y2, legend='slow group', color=palette[2])\n", "\n", "p_prob = bokeh.plotting.figure(width=400, height=300,\n", " y_axis_label='Probability',\n", " x_range=p_per.x_range)\n", "p_prob.line(x, p1_v1, legend='p_v1 fast', color=palette[3])\n", "p_prob.line(x, p2_v1, legend='p_v1 slow', color=palette[4])\n", "p_prob.line(x, p1_v2, legend='p_v2 fast', color=palette[3], line_dash='dashed')\n", "p_prob.line(x, p2_v2, legend='p_v2 slow', color=palette[4], line_dash='dashed')\n", "\n", "p_diff = bokeh.plotting.figure(width=400, height=100,\n", " x_axis_label='Time (m)',\n", " y_axis_label='Residuals',\n", " x_range=p_per.x_range)\n", "p_diff.line(x, p1_v2-p1_v1, color=palette[3])\n", "p_diff.line(x, p2_v2-p2_v1, color=palette[4])\n", "\n", "p = bokeh.layouts.gridplot([[p_per], [p_prob], [p_diff]])\n", "bokeh.plotting.show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This shows a general agreement though there are some important differences. In particular, the time at which the probability switches from one group to the other is around 7.5 minutes rather than just over 9 minutes. This agrees much better with what we would expect looking at the two distributions." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python [conda env:p36]", "language": "python", "name": "conda-env-p36-py" }, "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.6.0" } }, "nbformat": 4, "nbformat_minor": 2 }