robustness_of_the_two_sample_t_test_interactive.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Two-sample t-test and robustness\n",
    "**Posted on 11 May 2018 by John**\n",
    "source:\n",
    "https://www.johndcook.com/blog/2018/05/11/two-sample-t-test/ site visited on 28.05.2020\n",
    "\n",
    "This notebook conists the explanations and code samples from John D. Cook shown on his blog expanded with some interactivity. One can play with different settings and explore the effect when one of the distributions is not symmetric. \n",
    "Thus, the questions in the post \"...How far can you be from a normal distribution and still do OK? Can you have any distribution as long as it’s symmetric? Does a little asymmetry ruin everything? If something does go wrong, how does it go wrong?\" can be explored interactivly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# interaktive controls\n",
    "#\n",
    "import numpy as np\n",
    "from scipy.stats import norm, t, gamma, uniform, ttest_ind\n",
    "import matplotlib.pyplot as plt\n",
    "import ipywidgets as widgets\n",
    "\n",
    "%matplotlib widget\n",
    "\n",
    "#\n",
    "#\n",
    "style = {'description_width': 'initial'}\n",
    "#\n",
    "SldSampleSize = widgets.IntSlider(min=50, max=5000, value=1000, description=\"sample size\", style=style, continious_update=False)\n",
    "ITXTNTrials = widgets.IntText(min=1, max= 5000, value=1000, description=\"number of trials\", style=style, continuous_update=False, layout = {'width': '20%'})\n",
    "ITXTNPerGroup = widgets.IntText(min=1, max=5000, value=36, description=\"number per group\",  style=style, continuous_update=False, layout = {'width': '20%'})\n",
    "#\n",
    "TGBTNDistr1 = widgets.ToggleButtons(\n",
    "    options=['gauss', ''],\n",
    "    description='',\n",
    "    disabled=False,\n",
    "    button_style='', \n",
    "    tooltips=['Normal distribution', ''],\n",
    ")\n",
    "distrib_1 = widgets.VBox(\n",
    "    [\n",
    "        widgets.Label(value=\"1st distribution:\"),\n",
    "        TGBTNDistr1\n",
    "    ]\n",
    ")\n",
    "#\n",
    "TGBTNDistr2 = widgets.ToggleButtons(\n",
    "    options=['gauss', 'gamma'],\n",
    "    description='',\n",
    "    disabled=False,\n",
    "    button_style='', \n",
    "    tooltips=['Normal distribution', 'Gamma distribution'],\n",
    ")\n",
    "distrib_2 = widgets.VBox(\n",
    "    [\n",
    "        widgets.Label(value=\"2nd distribution:\"),\n",
    "        TGBTNDistr2\n",
    "    ]\n",
    ")\n",
    "#\n",
    "# Erwartungswert und Standrdabweichung der Standardverteilung\n",
    "#\n",
    "SldMean1 = widgets.FloatSlider(min=-50, max=150, value=100.0, description=\"\", continuous_update=False)\n",
    "SldSig1 = widgets.FloatSlider(min=0.1, max=50, value=15, description='', continuous_update=False)\n",
    "#\n",
    "set1 = widgets.VBox(\n",
    "    [\n",
    "        widgets.Label(value=\"$µ_1$ (mean value):\"),\n",
    "        SldMean1,\n",
    "        widgets.Label(value=\"$\\sigma_2$ (standard deviation):\"),\n",
    "        SldSig1\n",
    "        ]\n",
    ")\n",
    "#\n",
    "SldMean2 = widgets.FloatSlider(min=-50, max=150, value=100.0, description=\"\", continuous_update=False)\n",
    "SldSig2 = widgets.FloatSlider(min=0.1, max=50, value=15, description='', continuous_update=False)\n",
    "SldShift = widgets.IntSlider(min= 0, max=200, value=90, continuous_update=False)\n",
    "#\n",
    "set2 = widgets.VBox(\n",
    "    [\n",
    "        widgets.Label(value=\"$µ_2$ (mean value):\"),\n",
    "        SldMean2,\n",
    "        widgets.Label(value=\"$\\sigma_2$ (standard deviation):\"),\n",
    "        SldSig2,\n",
    "        widgets.Label(value=\"shift parameter for gamma distribution\"),\n",
    "        SldShift\n",
    "])\n",
    "#\n",
    "# diagram range selectors\n",
    "xRange = widgets.IntRangeSlider(min=-20, max=250, value=[40,140], style = style)\n",
    "xRng = widgets.VBox([\n",
    "    widgets.Label(value=\"diagram x-range:\"),\n",
    "    xRange\n",
    "])\n",
    "# yRange = widgets.IntRangeSlider(min=0, max=20, value=[0, 0.03], step=0.01, description = \"diagram y-range\", style = style)\n",
    "# buttons\n",
    "# \n",
    "BtnReload = widgets.Button(description=\"test\",button_style='info', tooltip=\"do an independent 2-sample t-test of the 2 distributions\")\n",
    "BtnReset = widgets.Button(description=\"Reset\", button_style='warning', tooltip=\"reset to default values\")\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulate_trials:\n",
    "# according to : https://www.johndcook.com/blog/2018/05/11/two-sample-t-test/ site visited on 28.05.2020\n",
    "# \n",
    "#\n",
    "# two-sample t-test with a certain amount of values of each group\n",
    "# \n",
    "#\n",
    "# parameter_in:  num_trials:  number of repetitions\n",
    "#                numPerGroup: number of random values selected for the test\n",
    "#                group1:      distribution 1\n",
    "#                group2:      distribution 2\n",
    "#\n",
    "# parameter_out: num_rejects: number of tests rejecting the null hypothesis\n",
    "#\n",
    "def simulate_trials(num_trials, numPerGroup, group1, group2):\n",
    "    num_reject = 0\n",
    "    for _ in range(num_trials):\n",
    "        a = group1.rvs(numPerGroup)\n",
    "        b = group2.rvs(numPerGroup)\n",
    "        stat, pvalue = ttest_ind(a, b)\n",
    "        if pvalue < 0.05:\n",
    "            num_reject += 1\n",
    "    return(num_reject)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize global variables\n",
    "#\n",
    "nSampleOld = 0\n",
    "mu1Old = 0\n",
    "sig1Old = 0\n",
    "mu2Old = 0\n",
    "sig2Old = 0\n",
    "n1 = 0\n",
    "n2 = 0\n",
    "x1 = 0\n",
    "x2 = 0\n",
    "rejects = 0\n",
    "dist2_old = ''\n",
    "shiftOld = 90\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot function\n",
    "#\n",
    "def plot_diagrams(dictDistributions, ax):\n",
    "    #\n",
    "    ax.clear()\n",
    "    x1 = dictDistributions['x1']\n",
    "    ax.plot(x1, dictDistributions['n1'].pdf(x1), color='C1', label='distribution 1')\n",
    "    #\n",
    "    x2 = dictDistributions['x2']\n",
    "    #\n",
    "    ax.plot(x2, dictDistributions['n2'].pdf(x2), color='blue', ls='-.', label='distribution 2')\n",
    "    #\n",
    "    ax.set_xlim(xRange.value[0], xRange.value[-1])\n",
    "#     ax.set_ylim(yRange.value[0], yRange.value[-1])\n",
    "    ax.legend(loc='best')\n",
    "    ax.set_facecolor('beige')\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_anotation(dictDistributions, ax):\n",
    "    # get informations about mean and rejects due to two-sample t.test\n",
    "    #\n",
    "    # clear current texts\n",
    "    for txt in ax1.texts:\n",
    "#         txt.set_visible(False)\n",
    "        txt.remove()\n",
    "    #\n",
    "    textstr = '\\n'.join((\\\n",
    "                         r'$\\mathrm{µ_1}=%.2f$' % (dictDistributions['mu1'], ),\n",
    "                         r'$\\mathrm{\\sigma_1}=%.2f$' % (dictDistributions['sig1'], ),\n",
    "                         r'$\\mathrm{µ_2}=%.2f$' % (dictDistributions['mu2'], ),\n",
    "                         r'$\\mathrm{\\sigma_2}=%.2f$' % (dictDistributions['sig2'], ),\n",
    "                         r'number of rejects: %3d' % (dictDistributions['rejects'], ),\n",
    "                         'ratio: {:.2f}%'.format(dictDistributions['ratio'], ),\n",
    "                        ))\n",
    "    props = dict(boxstyle='round', facecolor='khaki', alpha=0.5)\n",
    "    ax.text(0.05, 0.75, textstr, transform=ax.transAxes, fontsize=10, verticalalignment='top', bbox=props)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# callback function when button 'BtnReset' was pressed\n",
    "#\n",
    "def reset(*args):\n",
    "    set1.children[1].value = 100         # µ_1\n",
    "    set1.children[3].value = 15          # \\sigma_1\n",
    "    set2.children[1].value = 100         # µ_2\n",
    "    set2.children[3].value = 15          # \\sigma_2\n",
    "    #\n",
    "    distrib_1.children[1].value = 'gauss' # toggle-button_1\n",
    "    distrib_2.children[1].value = 'gauss' # toggle_button_2\n",
    "    #\n",
    "    xRange.value = [40,140]\n",
    "    SldSampleSize.value = 1000\n",
    "    ITXTNTrials.value = 1000\n",
    "    ITXTNPerGroup.value = 36\n",
    "    #\n",
    "    update_controls(True)\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# callback function when button 'BtnReload' was pressed\n",
    "#\n",
    "def do_test(*args):\n",
    "    update_controls(\"do test\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_controls(*args):\n",
    "    #\n",
    "    global nSampleOld, mu1Old, sig1Old, mu2Old, sig2Old, n1, n2, x1, x2, rejects, dist2_old, ax1, anTxt, shiftOld\n",
    "    #\n",
    "    #\n",
    "    dist1 = distrib_1.children[1].value\n",
    "    dist2 = distrib_2.children[1].value\n",
    "    #\n",
    "    # presets for distribution 2\n",
    "    #\n",
    "    signalChange = False\n",
    "    #\n",
    "    # normal distribution:\n",
    "    if (dist2 == 'gauss') and (dist2_old != dist2):\n",
    "        set2.children[1].value = 100\n",
    "        set2.children[3].value = 15\n",
    "        set2.children[5].disabled = True\n",
    "        set2.children[0].value = \"$µ_1$ (mean value):\".format(set1.children[1].value)\n",
    "        set2.children[2].value = \"$\\sigma_1$ (standard deviation):\".format(set2.children[3].value)\n",
    "        signalChange = True\n",
    "    #\n",
    "    # gamma distribution:\n",
    "    if (dist2 == 'gamma') and (dist2_old != dist2):\n",
    "        set2.children[1].value = 10\n",
    "        set2.children[3].value = 3.87\n",
    "        set2.children[5].disabled = False\n",
    "        set2.children[5].value = 90\n",
    "        set2.children[0].value = \"$µ_2$: (a = {:.2f})\".format((10 / 3.87)**2)\n",
    "        set2.children[2].value = \"$\\sigma_2$: (b = {:.2f})\".format(3.87**2/10)\n",
    "        signalChange = True\n",
    "    #\n",
    "    #\n",
    "    # just recalc distributions if parameters changed or BtnReload was pressed\n",
    "    #\n",
    "    mu1 = SldMean1.value\n",
    "    sig1 = SldSig1.value\n",
    "    mu2 = set2.children[1].value\n",
    "    sig2 = set2.children[3].value\n",
    "    shift = set2.children[5].value\n",
    "    nSize = SldSampleSize.value\n",
    "    nTrials = ITXTNTrials.value\n",
    "    #\n",
    "    reCalc = False\n",
    "    if (nSampleOld != nSize):\n",
    "        reCalc = True\n",
    "    #\n",
    "    if (mu1Old != mu1) or (sig1Old != sig1):\n",
    "        reCalc = True\n",
    "    #\n",
    "    if (mu2Old != mu2) or (sig2Old != sig2):\n",
    "        reCalc = True\n",
    "    #\n",
    "    if (shiftOld != shift):\n",
    "        reCalc = True\n",
    "    #\n",
    "    if reCalc or signalChange:\n",
    "        #\n",
    "        # 1st distribution\n",
    "        n1 = norm(mu1, sig1)\n",
    "        x1 = np.arange(n1.ppf(0.01), n1.ppf(0.99), 1/200)\n",
    "        #\n",
    "        # 2nd distribution\n",
    "        if dist2 == 'gauss':\n",
    "            n2 = norm(mu2, sig2)\n",
    "            x2 = np.arange(n2.ppf(0.01), n2.ppf(0.99), 1/200)\n",
    "        if dist2 == 'gamma':\n",
    "            #\n",
    "            # clc mu to a(shape) and sigma to b(scale)\n",
    "            # a = (mu/sigma)**2\n",
    "            # b = sigma**2/mu\n",
    "            # \n",
    "            a = (mu2 / sig2) ** 2\n",
    "            b = (sig2**2) / mu2\n",
    "            n2 = gamma(a = a, scale=b, loc=shift)\n",
    "            x2 = np.arange(n2.ppf(0.01), n2.ppf(0.99), 1/200)\n",
    "            #\n",
    "            set2.children[0].value = \"$µ_2$: (a = {:.2f})\".format(a)\n",
    "            set2.children[2].value = \"$\\sigma_2$: (b = {:.2f})\".format(b)\n",
    "            #\n",
    "        dictPlot = {\n",
    "            'x1': x1,\n",
    "            'n1': n1,            \n",
    "            'x2': x2,\n",
    "            'n2': n2,\n",
    "        }\n",
    "        #\n",
    "        plot_diagrams(dictPlot, ax1)\n",
    "        #\n",
    "        rejects = simulate_trials(nTrials, ITXTNPerGroup.value, n1, n2)\n",
    "        dictAnotate = {\n",
    "            'mu1': SldMean1.value,\n",
    "            'sig1': sig1,\n",
    "            'mu2': mu2,\n",
    "            'sig2': sig2,\n",
    "            'rejects': rejects,\n",
    "            'ratio': rejects / nTrials * 100\n",
    "        }\n",
    "        #\n",
    "        add_anotation(dictAnotate, ax1)\n",
    "    #\n",
    "    # do t-test when button 'test' was pressed\n",
    "    #\n",
    "    if args[0] == \"do test\":\n",
    "        rejects = simulate_trials(nTrials, ITXTNPerGroup.value, n1, n2)\n",
    "    #\n",
    "    dictAnotate = {\n",
    "            'mu1': mu1,\n",
    "            'sig1': sig1,\n",
    "            'mu2': mu2,\n",
    "            'sig2': sig2,\n",
    "            'rejects': rejects,\n",
    "            'ratio': rejects / nTrials * 100\n",
    "        }\n",
    "        #\n",
    "    add_anotation(dictAnotate, ax1)\n",
    "    #\n",
    "    nSampleOld = nSize\n",
    "    mu1Old = mu1\n",
    "    sig1Old = sig1\n",
    "    mu2Old = mu2\n",
    "    sig2Old = sig2\n",
    "    dist2_old = dist2\n",
    "    shiftOld = shift\n",
    "#\n",
    "SldSampleSize.observe(update_controls, 'value')\n",
    "ITXTNTrials.observe(update_controls, 'value')\n",
    "ITXTNPerGroup.observe(update_controls, 'value')\n",
    "SldMean1.observe(update_controls, 'value')\n",
    "SldMean2.observe(update_controls, 'value')\n",
    "SldSig1.observe(update_controls, 'value')\n",
    "SldSig2.observe(update_controls, 'value')\n",
    "SldShift.observe(update_controls, 'value')\n",
    "TGBTNDistr1.observe(update_controls, 'value')\n",
    "TGBTNDistr2.observe(update_controls, 'value')\n",
    "xRange.observe(update_controls, 'value')\n",
    "# yRange.observe(update_controls, 'value')\n",
    "BtnReload.on_click(do_test)\n",
    "BtnReset.on_click(reset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4083a186d86240d6a63acb8222fd3a65",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bec7cb86d0b34334af0ba9b4b498888f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HBox(children=(IntSlider(value=1000, description='sample size', max=5000, min=50, style=SliderS…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# interactive controls and diagram\n",
    "#\n",
    "# controls for distribution parameters\n",
    "#\n",
    "# diagram:\n",
    "#  - prob. density function of both distributions\n",
    "#  - mean-values, number of rejected hypotheses\n",
    "#\n",
    "fig1, ax1 = plt.subplots()\n",
    "#\n",
    "# arange controls\n",
    "intWidgets = widgets.VBox([\n",
    "    widgets.HBox([\\\n",
    "                  SldSampleSize,\n",
    "                  ITXTNTrials,\n",
    "                  ITXTNPerGroup,\n",
    "                  BtnReload,\n",
    "                  BtnReset\n",
    "                 ]),\n",
    "    widgets.HTML(value=\"<hr>\"),\n",
    "    widgets.HBox([\\\n",
    "                  widgets.VBox([\\\n",
    "                                distrib_1,\n",
    "                                set1\n",
    "                               ]),\n",
    "                  widgets.HTML(value=\"<div>\"),\n",
    "                  widgets.VBox([\\\n",
    "                                distrib_2,\n",
    "                                set2,\n",
    "                               ]),\n",
    "                  xRng\n",
    "                      ]),\n",
    "    widgets.HTML(value=\"<hr>\"),\n",
    "#     yRange\n",
    "])\n",
    "#\n",
    "# draw diagram for the 1st time with default settings\n",
    "update_controls(True)\n",
    "#\n",
    "intWidgets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Overview interactive controls:\n",
    "- distribution paramerters (µ, $\\sigma$)\n",
    "- sample size\n",
    "- number of t-tests\n",
    "- number of randomly picked numbers\n",
    "\n",
    "Overview of scipy.stats functions:\n",
    "- normal distribution: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html#scipy.stats.norm\n",
    "- gamma distribution: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gamma.html#scipy.stats.gamma\n",
    "- t-test : https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.html#scipy.stats.ttest_ind"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright © 2020 IUBH Internationale Hochschule"
   ]
  }
 ],
 "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.8.2"
  },
  "nbTranslate": {
   "displayLangs": [
    "de"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "en",
   "targetLang": "de",
   "useGoogleTranslate": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}