81_normal_t_comparison.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Valuation Under Independent t-Distributed Assumptions\n",
    "\n",
    "This notebook contains the empricial extensions described in Section 8.1 of the paper \"What is the value of experimentation and measurement?\"\n",
    "\n",
    "> We extend the model to incorporate the heavy tail in the distribution of uplifts shown in related work by modelling the true value of the propositions, as well as the estimation noise, with Generalized Student’s t-distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rulu.normal_normal_model as nnm\n",
    "import rulu.t_t_model as ttm\n",
    "from rulu.utils import get_test_params\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Case study - comparison of t-distributed / normal assumptions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_normal_t_improvement_against_M(\n",
    "    N, mu_X, sigma_X, mu_epsilon, sigma_1, sigma_2, M_list, verbose=False):\n",
    "    \n",
    "    # Constants\n",
    "    NUM_RUNS = 10000\n",
    "    T_DF = 3\n",
    "    CI_HIGH_PERCENTILE = 95\n",
    "    CI_LOW_PERCENTILE = 5    \n",
    "    \n",
    "    # Some functions take sigma, some take sigma squared\n",
    "    sigma_sq_1 = sigma_1 ** 2\n",
    "    sigma_sq_2 = sigma_2 ** 2\n",
    "    sigma_sq_X = sigma_X ** 2\n",
    "    \n",
    "    # Initialise empty list to hold the results\n",
    "    normal_E_Ds = []\n",
    "    normal_CI_lows = []\n",
    "    normal_CI_highs = []\n",
    "    t_E_Ds = []\n",
    "    t_CI_lows = []\n",
    "    t_CI_highs = []\n",
    "    \n",
    "    print(\"Parameters: N = {}, mu_X = {}, sigma_X = {}, \"\n",
    "      \"mu_epsilon = {}, sigma_1 = {}, sigma_2 = {}\"\n",
    "      .format(N, mu_X, sigma_X, mu_epsilon, sigma_1, sigma_2))\n",
    "\n",
    "    for M in M_list:\n",
    "        # Normal assumptions\n",
    "        normal_improvement = (\n",
    "            nnm.get_samples(\n",
    "                n_samples=NUM_RUNS, N=N, M=M, \n",
    "                mu_X=mu_X, mu_epsilon=mu_epsilon, \n",
    "                sigma_sq_X=sigma_sq_X, sigma_sq_1=sigma_sq_1,\n",
    "                sigma_sq_2=sigma_sq_2, verbose=False)['improvement']\n",
    "        )\n",
    "\n",
    "        normal_E_D = np.mean(normal_improvement)\n",
    "        normal_CI_low = np.percentile(normal_improvement, CI_LOW_PERCENTILE)\n",
    "        normal_CI_high = np.percentile(normal_improvement, CI_HIGH_PERCENTILE)\n",
    "\n",
    "        normal_E_Ds.append(normal_E_D)\n",
    "        normal_CI_lows.append(normal_CI_low)\n",
    "        normal_CI_highs.append(normal_CI_high)\n",
    "        \n",
    "        # t-assumptions\n",
    "        t_improvement = (\n",
    "            ttm.get_samples(\n",
    "                n_samples=NUM_RUNS, N=N, M=M, \n",
    "                mu_X=mu_X, mu_epsilon=mu_epsilon, \n",
    "                sigma_sq_X=sigma_sq_X, sigma_sq_1=sigma_sq_1,\n",
    "                sigma_sq_2=sigma_sq_2, nu=T_DF, verbose=False)['improvement']\n",
    "        )\n",
    "        \n",
    "        t_E_D = np.mean(t_improvement)\n",
    "        t_CI_low = np.percentile(t_improvement, 5)\n",
    "        t_CI_high = np.percentile(t_improvement, 95)\n",
    "\n",
    "        t_E_Ds.append(t_E_D)\n",
    "        t_CI_lows.append(t_CI_low)\n",
    "        t_CI_highs.append(t_CI_high)\n",
    "        \n",
    "        if verbose:\n",
    "            print(\"M = {}: normal improvement = {} ([{}, {}]), t improvement {} ([{}, {}])\"\n",
    "                  .format(M, np.round(normal_E_D, 3), \n",
    "                          np.round(normal_CI_low, 3), np.round(normal_CI_high, 3),\n",
    "                          np.round(t_E_D, 3), np.round(t_CI_low, 3), np.round(t_CI_high, 3)))\n",
    "            \n",
    "    # Prepare result in a dict\n",
    "    result = {\n",
    "        'normal_E_Ds': normal_E_Ds,\n",
    "        'normal_CI_lows': normal_CI_lows,\n",
    "        'normal_CI_highs': normal_CI_highs,\n",
    "        't_E_Ds': t_E_Ds,\n",
    "        't_CI_lows': t_CI_lows,\n",
    "        't_CI_highs': t_CI_highs,\n",
    "    }\n",
    "            \n",
    "    return result\n",
    "\n",
    "\n",
    "def generate_normal_t_improvement_against_M_plot(M_list, result):\n",
    "    fig, ax = plt.subplots()\n",
    "    \n",
    "    ax.errorbar(result['normal_E_Ds'], M_list, \n",
    "                xerr=[np.array(result['normal_E_Ds']) - np.array(result['normal_CI_lows']), \n",
    "                      np.array(result['normal_CI_highs']) - np.array(result['normal_E_Ds'])],\n",
    "                fmt='x', capsize=3)\n",
    "    \n",
    "    # Slightly displace the mean and error bar under t-distributed assumptions\n",
    "    # so that it will be visible\n",
    "    ax.errorbar(result['t_E_Ds'], \n",
    "                (np.array(M_list) * 1.1 if params['N'] > 500 else np.array(M_list) + 2), \n",
    "                xerr=[np.array(result['t_E_Ds']) - np.array(result['t_CI_lows']), \n",
    "                      np.array(result['t_CI_highs']) - np.array(result['t_E_Ds'])],\n",
    "                fmt='^', capsize=3, alpha=0.85)\n",
    "    \n",
    "    ax.axvline(0, c='black')\n",
    "    \n",
    "    # The zero line has to be visible\n",
    "    ax.set_xlim(min(min(result['normal_CI_lows'])-0.01, min(result['t_CI_lows'])-0.01, -0.01) * 1.05, \n",
    "                max(max(result['normal_CI_highs']), max(result['t_CI_highs']), 0.01) * 1.05)\n",
    "    \n",
    "    if params['N'] > 1000:\n",
    "        ax.set_yscale('log')\n",
    "    \n",
    "#     ax.set_xlabel('Value gained')\n",
    "#     ax.set_ylabel('M')\n",
    "    \n",
    "    return fig, ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_list = [\n",
    "    {'N': 6700, 'mu_X': 0, 'mu_eps': 0, 'sigma_X': 0.6, 'sigma_1': 1, 'sigma_2': 0.4,\n",
    "     'M': [10, 20, 50, 100, 200, 500, 1000], 'name':'ecomm'},\n",
    "    {'N': 6700, 'mu_X': 0, 'mu_eps': 0, 'sigma_X': 0.6, 'sigma_1': 0.6, 'sigma_2': 0.4,\n",
    "     'M': [10, 20, 50, 100, 200, 500, 1000], 'name':'ecomm'},\n",
    "    {'N': 184, 'mu_X': 19.9, 'mu_eps': 0, 'sigma_X': 10, 'sigma_1': 2, 'sigma_2': 0.4, \n",
    "     'M': [10, 20, 40, 60, 80, 100], 'name':'marketing'},\n",
    "    {'N': 184, 'mu_X': 19.9, 'mu_eps': 0, 'sigma_X': 10, 'sigma_1': 0.8, 'sigma_2': 0.4, \n",
    "     'M': [10, 20, 40, 60, 80, 100], 'name':'marketing'},\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters: N = 6700, mu_X = 0, sigma_X = 0.6, mu_epsilon = 0, sigma_1 = 1, sigma_2 = 0.4\n",
      "M = 10: normal improvement = 0.617 ([0.3, 0.935]), t improvement 2.19 ([1.095, 3.306])\n",
      "M = 20: normal improvement = 0.579 ([0.357, 0.807]), t improvement 1.661 ([1.034, 2.291])\n",
      "M = 50: normal improvement = 0.527 ([0.388, 0.667]), t improvement 1.119 ([0.83, 1.417])\n",
      "M = 100: normal improvement = 0.481 ([0.384, 0.578]), t improvement 0.806 ([0.643, 0.971])\n",
      "M = 200: normal improvement = 0.433 ([0.366, 0.498]), t improvement 0.567 ([0.474, 0.66])\n",
      "M = 500: normal improvement = 0.36 ([0.32, 0.4]), t improvement 0.341 ([0.298, 0.384])\n",
      "M = 1000: normal improvement = 0.297 ([0.271, 0.322]), t improvement 0.223 ([0.199, 0.246])\n",
      "Parameters: N = 6700, mu_X = 0, sigma_X = 0.6, mu_epsilon = 0, sigma_1 = 0.6, sigma_2 = 0.4\n",
      "M = 10: normal improvement = 0.242 ([-0.025, 0.513]), t improvement 0.876 ([-0.115, 1.904])\n",
      "M = 20: normal improvement = 0.228 ([0.042, 0.418]), t improvement 0.659 ([0.099, 1.225])\n",
      "M = 50: normal improvement = 0.207 ([0.089, 0.324]), t improvement 0.441 ([0.191, 0.693])\n",
      "M = 100: normal improvement = 0.189 ([0.11, 0.268]), t improvement 0.316 ([0.178, 0.456])\n",
      "M = 200: normal improvement = 0.17 ([0.117, 0.223]), t improvement 0.223 ([0.146, 0.299])\n",
      "M = 500: normal improvement = 0.142 ([0.11, 0.173]), t improvement 0.133 ([0.099, 0.168])\n",
      "M = 1000: normal improvement = 0.117 ([0.097, 0.137]), t improvement 0.087 ([0.068, 0.106])\n",
      "Parameters: N = 184, mu_X = 19.9, sigma_X = 10, mu_epsilon = 0, sigma_1 = 2, sigma_2 = 0.4\n",
      "M = 10: normal improvement = 0.375 ([0.0, 0.913]), t improvement 0.382 ([0.0, 1.383])\n",
      "M = 20: normal improvement = 0.32 ([0.051, 0.673]), t improvement 0.349 ([0.02, 0.89])\n",
      "M = 40: normal improvement = 0.249 ([0.079, 0.456]), t improvement 0.278 ([0.072, 0.564])\n",
      "M = 60: normal improvement = 0.206 ([0.078, 0.36]), t improvement 0.23 ([0.082, 0.421])\n",
      "M = 80: normal improvement = 0.168 ([0.067, 0.289]), t improvement 0.189 ([0.074, 0.336])\n",
      "M = 100: normal improvement = 0.136 ([0.055, 0.231]), t improvement 0.154 ([0.06, 0.273])\n",
      "Parameters: N = 184, mu_X = 19.9, sigma_X = 10, mu_epsilon = 0, sigma_1 = 0.8, sigma_2 = 0.4\n",
      "M = 10: normal improvement = 0.049 ([-0.048, 0.211]), t improvement 0.047 ([-0.046, 0.253])\n",
      "M = 20: normal improvement = 0.041 ([-0.026, 0.143]), t improvement 0.044 ([-0.036, 0.192])\n",
      "M = 40: normal improvement = 0.032 ([-0.011, 0.092]), t improvement 0.039 ([-0.02, 0.131])\n",
      "M = 60: normal improvement = 0.026 ([-0.006, 0.071]), t improvement 0.034 ([-0.011, 0.104])\n",
      "M = 80: normal improvement = 0.022 ([-0.004, 0.056]), t improvement 0.029 ([-0.008, 0.082])\n",
      "M = 100: normal improvement = 0.017 ([-0.003, 0.044]), t improvement 0.023 ([-0.007, 0.065])\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM8AAACZCAYAAACIVm61AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAMuUlEQVR4nO3dfYwU9R3H8feXu4Mzel4rKKVASwU0NRZ7qYFW0icf0dZqjGkrtjaNEWlKpfWPxqZtmv5hyx99EjGRU6w1FbWtxqC2EUp9SMUeD4JEpQoYo0eN9rQgkgMO7ts/du/Y29udm52Z3dmd/bwSIjO7M/M7L1/nNz9/n9+YuyMilRuXdgNEGpWKRyQiFY9IRCoekYhUPCIRtabdAIC2tjY/66yz0m6GyChbtmzpc/eTS31m9TBU3dHR4fv370+7GSKjmNkWdz+71GfqtolEVBfdttCe/jW89QJ86BPw2RvTbo00uca682y4FXauhWeWp90SkeSLx8w+bma3m9lfzOw7iZ34QB8c3As+yP6Dh/npvU8Of/Tbda+wYXcftz+1O7HLiYwlVPGY2V1m9raZvVC0f4GZvWxmu8zsJgB33+Hui4GvAvMTa2nPyvxFWxhvA0x86Q9s2N0HwC3rd7Jk9VbmTOtM7HIiYwl757kbWFC4w8xagNuAi4EzgKvM7Iz8Z18BHgP+mkgrD/TB86vBWqBjMhMY5MqWp1i66u/MX7YegBULuzhn5qRELicSRqjicfengXeLds8Fdrn7q+5+GLgfuCz//TXufjFwdblzmtkiM9tsZpsHBgaCG9CzEgaPgtnQwbQwyOLj1rNn70EAFY7UXJxnnqnAGwXbvcBUM/uCmS03s5UE3Hncvdvdz3b3s9va2oKv9OY2GByA1vG5557W8bTZUU4/upMbzp0FMNyFE6mVxIeq3f1J4MlET3r1n0dsbtjdx8I7elh93TxunDmJ5f/YxZLVW9V1k5qKUzx7gOkF29Py+0Izs0uBS9vb2yu68PbefVzRNXW4UJaeN5t5p57E9t59Kh6pmdDTc8xsBvCou5+Z324FXgHOI1c0m4CF7v5ipY3Q9BypV7Gn55jZfcCzwOlm1mtm17r7EWAJ8DiwA/hTlMIRaVSpTgwt6LZd19/fn1o7RMoJuvNoVrVIAM2qFqmCVGdVRx1tE6kH6raJBMh2t61nZW7u25+/DQfeOTaBVKTKGr/btrE7Vzw718LGlfDCgzDv+sTaKFJO43fbbv4wUPAznDAZlm5LpF0ime62bT48nb1HWnNxBR+E/r3c+MA2BeOk6hq7eA70cea41/jf0Xb62ydBxxQOHurnia07FIyTqmvs4ulZSTsDDNLCm/sO8vreQ+w72s63WtZqgqhUXWMPGLy5Ddra6Rw8wDgzBgedtnFHmeOvJtpOkVIaf8AAmHHTY5x0/Hi+Me8j/LHndd49cJjXln0pwRZKs8r0gMFQgnTFwi5uvPB0VizsGrFfpFoavniKg3HnzJzEFV1T2d67L+WWSdZlotsmUi1B3bbGHjAQSZHuPCIBMj1gIJIWFY9IRCoekYhUPCIRqXhEItJQtUhEGqoWCaChapEqUPGIRJT94ileTadwpZ1y3xEJIfvFs7F75HbPymMr7ZT7jkgIqY62Je3qO//Fd784a0QEe8P/Ornt5t9z78efhSMHYee6/Ae3wn9fhtYJKbVWGl2m7jzP7HqHJau3DgfhNuzuY0n/Ip7Zf0ruC307cyvsDK2007czxdZKo8tU8QCc0jGBa1ZtZP6y9VyzaiM/n3Bv7oOLfgGH3oOOKdAxOffPQ+/BRb9Mt8HSsFItHjO71My6jxw5ktg5O49rY/KJE9iz9yCTT5zAp1ry67cNvVF7XEtue1wLDB4Z+ewjUoFUi8fdH3H3Ra2tyT16LT1/Nv0Dg9xw7iz6BwZ5reuHuQ+G3qh9cO+xP4MD8J+tMHdRYteX5pGpAQNgxFuxPz1zIktWAxwe9UZtkbgy9cwzf9bEEa+TP2fmJFYs7GL+rIkpt0yySHPbRAJobptIFah4RCJS8YhEpOIRiUjFIxKRikckIhWPSEQqnqAgXKngXCXHS6YlXjxmdrmZ3WFmD5jZhUmfP3FBQbhSwblKjpdMCzW3zczuAr4MvO3uZxbsXwDcArQAd7r7Mnd/GHjYzD4I/ApYm3yzy/vtulf4wQWncftTu5kzrXNkMG53H9t799F/+Cg/uOC03M79b8FDJSaGKjgnYwh757kbWFC4w8xagNuAi4EzgKvM7IyCr/wk/3lN3bI+F3CbM61zdDBu9VbmTOsc/k4gBedkDKHuPO7+tJnNKNo9F9jlnnt7rpndD1xmZjuAZcDf3P25BNsa2tdWPgscC8ZNPnECb713iFmnnMAtfy8qgo7JcEVR1+tAH3R/PheYG9eSywENBeeOL5pkeuunqviTSD2L88wzFXijYLs3v+97wPnAlWa2uNzBZrbIzDab2eaBgYEYzSivOBjXeVxbuAMVnJMQEs/zuPtyYHmI73UD3ZCbVZ1kGx64/jPAsa7aDefO4o89r7P0/NmcM3MSM2567NiXSwXhCoNzhf6zdfR3FaRrWnGKZw8wvWB7Wn5faNVcq3qocEYG47YOvy172LzrRx9cSXCu1PHSFOJ02zYBs83sY2Y2Hvg6sKaSE1Qlhn3ebCD3luxSwbjtvfuGvyMSR6gwnJndB3wBmAS8BfzM3VeZ2SXA78gNVd/l7jdHaYTCcFKvgsJwqSZJC7pt1/X396fWDpFy6rZ4hujOI/VKMWyRKtCb4UQiUrdNJIC6bSJVoOIJI0pmJ0wWKOlrSk1lbqH3qoiS2QmTBUr6mlJTqQ4YuPsjwCMdHR3XJXXOoTzPkKBcz+LPzyx73AjlMj/lKAvUFDLXbSvO6gTleoKOi0VZoKaQubckwLE8z5ByuZ5R2Z5ySmV+yqkkCxREOaG61xTPPJFzPVEoC9Q0MvfMA8fyPEPK5XoKjcj4FKsks1NJFiiIckJ1L5PdtkJBuZ7iAiqrksxOUi/RUk6o7mVuwKA4qxOU6wk6TmQsmp4jEiBoeo4mhopEpDuPSABNDBWpAhWPSEQqHpGIVDwiEal4RCLSULVIRBqqFgmgoWqRKlDxiESk4hGJSMWTpDRWvIm7Sk9cTbzKj4onSWmseBN3lZ64mniVn0yF4YJWwAm7ik4l5xyl0lV24tIqPanK1J0naAWcsKvoVHLO1GmVnlRl6s4Do1fOKRR7FZ2xVLLKTlxJrdITVxOv8tMUq+cMqekqOtWmVXpSl7nVc4pXzikUZhWdYoGr6hSr5Yo3Sa3SE1cTr/KTuW5bOYmsojOWWq54k9QqPXE18So/mRowCFoBJ+wqOpWcU5qbJoaKBNDEUJEqUPGIRKTiEYlIxSMSkYpHJCIVj0hEKp5qyVrOJe3cULXE+D2peKolazmXtHND1RLj95T49BwzOxX4MdDp7lcmff5iFeVtiJ7rqfi6tc72VJNyQyWFuvOY2V1m9raZvVC0f4GZvWxmu8zsJgB3f9Xdr61GY0upNG8TNdcT97oNTbmhksLeee4GVgD3DO0wsxbgNuACoBfYZGZr3P2lpBs5lqAMTylVz/VAbbM91VQvuaFqiZFHCnXncfengXeLds8FduXvNIeB+4HLwl7YzBaZ2WYz2zwwMBC6wUnIVK6n2pQbKivOM89U4I2C7V5gnplNBG4GuszsR+7+y1IHu3s30A25iaEx2hGY4SklSq6n2Jg5n6zkXOolN1QtMX5PiQ8YuPs7wOKkz5uUmuR6IDs5l3rJDVVLjN9TnKHqPcD0gu1p+X2hJRHDrjRvEzXXE/e6kj2h8zxmNgN41N3PzG+3Aq8A55Ermk3AQnd/sdJGKM8j9Sp2nsfM7gOeBU43s14zu9bdjwBLgMeBHcCfohSOSKNKNUk69H4e4JtAoxfeJKAv7UbEoPaX9lF3P7nUB3URwzazzeVujY2i0X8Gtb9ymtsmEpGKRySieimeDMxjafifQe2vUF0884g0onq584g0HBWPSESpF0+pTFCjKJdzahRmNt3MnjCzl8zsRTNbmnabKmFm7Wa20cyez7f/5zW9fsr/k7SF3BSf4UwQcFUamaAozOxzwPvAPUPTlhqJmU0Bprj7c2bWAWwBLm+gf/8GHO/u75tZG/BPYKm7/6sW10/7zhMrE5S2MjmnhuHub7r7c/m/7yc3zWpquq0Kz3Pez2+25f/U7G6QdvGUygQ1zC8vS/ITf7uAnnRbUhkzazGzbcDbwDp3r1n70y4eqQNmdgLwIPB9d38v7fZUwt2PuvsnyUVi5ppZzbrPaRdP7EyQxJN/VngQuNfdH0q7PVG5+17gCWBBra6ZdvFsAmab2cfMbDzwdWBNym1qGvkH7lXADnf/TdrtqZSZnWxmH8j//ThyA0//rtX1Uy2eRs8Elco5pd2mCs0nFwc518y25f9cknajKjAFeMLMtpP7D/E6d3+0VhfX9ByRiNLutok0LBWPSEQqHpGIVDwiEal4RCJS8YhEpOIRiej/rIkCjKQlDmMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 216x162 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM8AAACZCAYAAACIVm61AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAOQUlEQVR4nO3df4wc9XnH8ffjtX04jn2OuNghNuolPnDjIlcmKwOmak78SExS45bSAqYCUlrbpFcQ/QNRtWqUSm34oxUCjHQ22CK0NRisCBkwIpDKodI5Z59jY/FDtu8oDWccuMWcbcz5x909/WP3Lnvr3b252V+zc5+XZOl2Zm720Wkfz8x35/Mdc3dEZOKm1LoAkXql5hEJSc0jEpKaRyQkNY9ISFNrXQBAU1OTNzc3F93m4MGDACxatKgKFYmk7d27N+XuX863LhLN09zcTFdXV9FtWltbAdi5c2flCxLJMLP/K7ROp20iIcWjed74d3j++3DqE+jcUOtqZJKIR/N0PAaHfwa7N8DujbWuRiaJsjePmX3DzNrNbJuZ3VPu/Z/nVIqTA6c5PWSc+J92GB4C4O+27qf9Fz0Vf3uZvAI1j5ltNrOPzeytnOUrzOygmXWb2YMA7v6uu68D/hy4uvwl5+jcwHQb5DeDMzk1OAUG+unoSfHTfUdYsqCx4m8vk1fQI89TwIrsBWaWAB4HbgAWA7eZ2eLMuhuBl4EdZas0n1MpeHMLDebMa5zBaabRP3CW+za9DsDyhU0VfXuZ3AI1j7u/ARzLWbwM6Hb399z9LPAssCqz/XZ3vwG4vdA+zWyNmXWZWVdfX1+46js3pE/TzJgxLcEwCQZ8Outm/Dzc/kQmoJTveeYDH2S97gWuMLNW4CaggSJHHnffCGwESCaT4XIRR/fD8DmYOp2zpz6l0c4xI+EsGjocanciE1H2L0ndfSews9z7zev25wHo6EnRtmUfx86c5f2HvseUnhQ80UlHT0qnblIxpYy2HQEuznq9ILMsMDNbaWYbjx8/XkIZcKD3OOtXL+W+ay8B0tc6Ny2dz4He0vYrUowFTZKaWTPwkrtflnk9FTgEXEu6afYAq9397YkWkUwmXbfnSBSZ2V53T+ZbF3So+hlgF7DIzHrN7G53HwTagFeBd4HnwjSOSL0KdM3j7rcVWL6DEoajzWwlsLKlpSXsLkRqpqa357j7i+6+prFRX2ZK/YnHvW0iNVDT5inXaJtILei0TSSkeJ22jWR5TqWU75GKi9dp20iWp3OD8j1ScfE6bTv5ETx3B+xaDz6cDsll8j0i5Rar07b2M9/m1cMn041jCfBhPvjMuP3JX9a6NImhWDXP5V9McdmZNxm4YC7MmsfABXNJnO3nne7/rXVpEkOxap5lQweYwjBHT57j18c+5+jJc3zJPuPOxM9qXZrEULwGDGbNZZoNMcc+5wvDnzHHPgdLsGTKe+XZv0iWeA0Y/M1ukmfauS6xiaevfp3rEpvY9xdv8ZfnHijP/kWyRGLG0HLp6EkBsH71UpYvbOLKhRfStmVfjauSuIrVNc+B3uPctHT+aHp0+cIm1q9eytUtF9a4MomjWB151n1r4XnLli9sUhRbKiJeAwYiVRSvAQORKorVNY9INal5REJS84iEpOYRCUnNIxKShqpFQtJQtUhIOm0TCUnNIxJSvJsn38w52TPrBNlepIB4N0++mXOyZ9YJsr1IAbG6qxrg4dcOcf/1l9L+ix6WfNrI8p+uGV23+9OZfOPINoY8wZyOx6DvIExtqGG1Us9id+R55OfpRyouWdBI28AaOk6k4wgdJ5rY/36K6TbIJ0NfSM+wk9LjFyW82B15AG7ZsAuAHzX8F3ccuod5sxs4d6KPlxqe41Ofnd5o1kVw5gR858cwMxOWe+ybNapY6lGsvyT9ZqKHebMbONJ/mnUzXidhwwxbIr1ySgKGB/Nf+4gEEMsvSbeuvYqta6/i/aUPMHBumHuvaeHSoW5mTYevXnCGRjsFp/vTT9L+MGuOg2VrCu9UJEcsT9sg84TsffNHJwPpWLiNK7fsY/2fLWX1E528/6Pvnf9LV6ytfqFSt2I3YDDyROyRJ2TnTgZyoPf46DYipYjdkef+6y8FNBmIVF7sjjwi1aLmEQlJzSMSkppHJCQ1j0hIah6RkNQ8IiFNjuYJGnIrFpQrZb8SS2VvHjP7YzN7wsy2mtm3y73/UIKG3IoF5UrZr8RSoDsMzGwz8EfAx+5+WdbyFcAjQAJ40t0fcvcXgBfM7EvAvwFVfSDoSBgOSAfiFjSy/ORHkAnFdZxo4sCpOay7qDu9fe8i7l9wEAZPw+HX0jtRUE4CCHrkeQpYkb3AzBLA48ANwGLgNjNbnLXJP2bWV9VIGA4ygbgt++gYTDdTx4km2rqTLJnZ/9vtP/zd9A+pw2MeQa+gnIwn0JHH3d8ws+acxcuAbnd/D8DMngVWmdm7wEPAK+7+qzLWGthIGA5g7qwGfvCbtcw81MRHJ87QMveLPHJmJY/0jWxxDL7zr7DxW+mA3JQEDA+dH5TLR+G5Sa2Ua575wAdZr3szy/4WuA642czWFfplM1tjZl1m1tXX11dos5I1zpjGn0zbxZH+08yb3UDjjGnnb9S5Id0wUxSUk+DKfle1uz8KPBpgu43ARoBkMunlrGHr2qtGf+7oSdH2k2u495oW/rPz19x33SVj7qxufvBlOLo/HYw73T92Rx+O8zBghecmtVKa5whwcdbrBZllgZnZSmBlS0tLCWUU1tGTom3LPtbfuXzM07Gzcz4A3P58uDdQeG5SK+W0bQ9wiZl9zcymA7cC2yeyg0rEsLODbsUCcfm2F5mIoEPVzwCtQJOZ9QI/dPdNZtYGvEp6qHqzu79dsUoDGhmmhmCBuOztRSYi6GjbbQWW7wB2hH3zSp+2iVRSLGfPEamGyXFvm0gFxHrSQ5FK0mmbSEg6bRMJafI0TzmzNxPN/QShbFDdmTzXPOXM3kw09xOEskF1p6Yzhrr7i8CLyWTyr0vdV3aOJ1e+B13l5nrG7Gsk45OPcj+SEZvTtuwcT658D7rKzfWM2ddIxicf5X4kI1ZzVWfneHJlP+gqf64n2zG4Kc9p1KlUuNxPEMoG1Z1Jc82T/aCrgrme8Sj3I1lic80DY3M8uTpeeoCBfcMFcz3Zmh98Of9OwuZ+glA2qO7E6rStkNwHXRXM9YwnbO4nCGWD6k5sBgyK5XKC5HqC7ktkhLmXNQEdSjKZ9K6urqLbtLa2ArBz587KFySSYWZ73T2Zb92kGTAQKTfdGCoSUmyueUSqTc0jEpKaRyQkNY9ISGoekZA0VC0SkoaqRULSaZtISGoekZDUPCIhxat56mkGmkrMwFMN9fQ3rrB4NU89zUBTiRl4qqGe/sYVVtdhuPNmzMl66nWu9qMtLJnZz/LZqdFlxWbQKfiexWbWCUoz8MRCXR95is2Yk2vJzH7aupOBZ9Ap+J7FZtYJSjPwxEJdH3lg7Iw5Dw/N5v6+7xfcdm7TOe449AcBZ9AppMDMOkFVcgaeatAsP6Mm1R0GjTOmlT6DTqk0A09s1P3sOWNmzOm8j61XFJlBJ/OA3yAz6BRScGadoCo5A081aJafUXV/2jZGkRloRp+MXeoMOqWq5Aw81aBZfkbV9YDBRGa5megMOuV4T4k3zZ4jUkRkZ88RqWdqHpGQ1DwiIal5REJS84iEpOYRCSk6zaOcyMTVayYoakJ+9qLTPMqJTFy9ZoKiJuRnr+y355jZ14F/ABrd/ebAv1gki/Nw7yJIHUq/KLBNEOXK9BRTlrxPEMoE1VygI4+ZbTazj83srZzlK8zsoJl1m9mDAO7+nrvfXc4iy5KhoXyZnmLKVeu4lAmquaBHnqeA9cDTIwvMLAE8DlwP9AJ7zGy7u78TqpJZ8wrnZHa/zDtnvwLALUXyOkGUJ9NTTIl5nyDqPRMUNSEzSoGOPO7+BnAsZ/EyoDtzpDkLPAusCvrGZrbGzLrMrKuvr2yf3HFFItNTKmWCIqGUa575wAdZr3uBK8zsQuBfgKVm9vfu/uN8v+zuG4GNkL4xlGXFIz2LvzobKP7E6yDKkekppuS8TxD1ngmKmpAZpbIPGLj7J8C6Cf9iFXIikcn0lKreM0FRE/KzV8pQ9RHg4qzXCzLLAgsawy5XhqZcmZ5ilPeZPALnecysGXjJ3S/LvJ4KHAKuJd00e4DV7v72RItQnkeiquQ8j5k9A+wCFplZr5nd7e6DQBvwKvAu8FyYxhGpVzVNkprZSmAlcAsQ5IuKJiA17la1E+X6VFs4v+PuX863IhIx7KDMrKvQITQKolyfaiu/6NzbJlJn1DwiIdVb80T91uso16fayqyurnlEoqTejjwikaHmEQkpks2TLyeUs77BzLZm1ndm7n6ISm13mVmfme3P/PurKtaWN3eVtd7M7NFM7QfM7PII1dZqZsez/m7/VK3aQnP3SP0DEkAP8HVgOvAmsDhnmx8A7ZmfbwW2Rqi2u4D1Nfrb/SFwOfBWgfXfBV4BDLgS6IxQba2kb/+q+Wcw6L8oHnmC5IRWAT/J/LwNuNbMLCK11Yznz11lWwU87Wm/BOaY2UURqa3uRLF58uWE5hfaxtP32B0HqhGhDFIbwJ9mTou2mdnFedbXStD6a+UqM3vTzF4xs9+rdTHjiWLz1LsXgWZ3XwK8xm+PkFLcr0jfR/b7wGPACzWuZ1xRbJ4gOaHRbTLRiEagGpOXjVubu3/i7mcyL58EovQQz5IzWJXi7ifc/bPMzzuAaWYW6YRiFJtnD3CJmX3NzKaTHhDYnrPNduDOzM83A//tmavOWteWcw1xI+m4RlRsB+7IjLpdCRx396O1LgrAzL4yct1qZstIfzYjPZtj5B6r6O6DZjaSE0oAm939bTP7Z6DL3bcDm4D/MLNu0heht0aotnvN7EZgMFPbXdWoDUZzV61Ak5n1Aj8EpmVqbwd2kB5x6wY+B0qbiqi8td0M3GNmg8AAcGuV/kMMTbfniIQUxdM2kbqg5hEJSc0jEpKaRyQkNY9ISGoekZDUPCIh/T/5tTy9eDDbbgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 216x162 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM8AAACZCAYAAACIVm61AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAO4klEQVR4nO3df5BV5X3H8ff37k9YFxBRQiEWBIaM06FVHIuYZJzQkdQamaFiLAajEwc7DdVgnIp2Okk606mZ6dTSkBFptIWkUQx1WjomTSyGSRyatRCNsbUIGtFFWHbD8puFvbvf/nHOwt3l7u69z7l777m7n9fMDveec+45z2Hny3POw3M/x9wdESleptINEKlWKh6RQCoekUAqHpFAKh6RQLWVbgBAXV2dNzY2Mm/evEo3RaSf3bt3d7j75fnWpaJ4GhsbWbBgATt27Kh0U0T6MbP9g63TZZtIoOoonpanCtvuVAd871449euR2b9Ijuoonlc3FrZdy1Ow90fwapHFUOj+RXIMe89jZs8AtwKH3f234mWTgS3ATOA94A537zQzA9YBtwCngXvc/edJG3mks5PJL6y6aPmGg3OY33SURRM6INsFe1+iq7eG7p8+SXP7HqhtyLu/J1rnsWbGnqTNkjGukJ7nn4BPD1i2Ftju7nOB7fF7gN8H5sY/q4AnS9HIzmx93uXzm46yet917Dw+BTr2crYXDmWbqLcsdOwddH/rPvxYKZolY9ywPY+7/8TMZg5YvBS4KX69CdgBPBIv3+zRbNOfmdkkM5vm7geTNvSz7ffmXX7FlG6+9PZ8ttVvocsvobamhk6fQE1HO4/aMk5kJub51BFYlnOp9o0FSZsnY1DoPc/UnII4BEyNX08HPsjZrjVedhEzW2Vmu8xsV3d3d2AzYOK4Ov543H/ivT1kMrVkDHqthgw9LDm1LXi/IsNJ/P887u5mVvT3Gtx9I7ARoLm5ecjPb+65mS3335B33c53Ouj99j4mNEBX90ma6+qorzHAWD6tneV3Xfy5mWtf7L/g+ovvp0SGE1o8bX2XY2Y2DTgcLz8AfDRnuxnxskQ29Szha3mW73yng9XffY31K7fSNHsKv+h7v/waFs2eUvgBfvf+pE2UMSj0sm0b8Pn49eeBf8tZfrdFFgLHSnG/8+DiuXmXv9F6jPUrLhTKotlTWL/iGt5oPRa0P5Fi2HDfJDWzZ4kGB6YAbcBXgH8FngeuBPYTDVUfiYeq1xONzp0G7nX3XcM1orm52TU9R9LIzHa7+3X51hUy2vZHg6xanGdbB75YXPNEqlN1zDAQSSEVj0ggFY9IIBWPSCAVj0ggFY9IIBWPSCAVj0ggFY9IIBWPSCAVj0ggFY9IIBWPSKDqKp5i89VCc9ySHlfGhOoqnmLz1UJz3JIeV8aERBkGZrYGuA9w4JfAvcA04DngMmA3sNLdzyVsZ+REG+TJb8srznEDYOc3YIgcN5EQwT2PmU0HHgCui8MQa4A7ga8DT7j7HKAT+EIpGvrES28Pu82Gg3OiDDeIctu8lzNez8mezJA5bueP0aqnNEjhkqbn1ALjzKwbGA8cBD4FrIjXbwK+SgnCD9dt38ua6VP7560NMD8OANmw7Equ37+MM41X0HYyy9RLauHscVjy19B02eDHWPsiax748sUrlOsmeQQXj7sfMLO/Ad4HzgA/IrpMO+ru2XizQXPbQnx49AxrnvqvIbe5ormBlmcfZ1btaU54A/U1GTrP9HDau3j5yb9ga/PKUjVHxrgkl22XEiWEzgJ+A2ji4ljeoT5fdOjhD5tuG3abiePqWFD/HuZZJtlpmjlJk5+kliyzzw1/6ZeXct0kjySXbb8H/Mrd2wHM7AXgRmCSmdXGvc+guW3FhB72ufehx8kfuntBlOX2GJ+78Uq+0/J+v2iqa4nS6QdzURhiH+W6SR5JhqrfBxaa2fg4cmox8L/Aj4Hb421yM91G3PkQxBXX8NDN81i/4hpWf/c1dr7TUa4myBgSXDzu3gJsBX5ONEydIepJHgEeMrN9RMPVT5egnQUFFYaGIBZzDJE+w4YeloNCDyWthgo9rK4ZBiIpouIRCaTiEQmk4hEJpOIRCaTiEQmk4hEJpOIRCaTiEQmk4hEJpOIRCaTiEQmk4hEJVH3FkzRDrVRZbgMp223Mqb7iSZqhVqost4GU7TbmJE3PKb9istsGUpablFCinsfMJpnZVjP7PzN7y8xuMLPJZvaSme2N/7y0VI0tJLstV78cN4COvZzthaO9jeC9BWW59Tu+ct0kR9KeZx3wH+5+u5nVE2W3PQZsd/fHzWwtsJboq9mJFZLdlmt+TqbBoo9A1zc/weGeZi6fMB5qJhaU5dbv+IPluoGy3cag4OIxs4nAJ4F7AOJI3XNmthS4Kd5sE7CDEhUPFJbdluuK5gbufvpVHhv3Ard0nyFTO57OU+foBBqU5SYJJOl5ZgHtwD+a2W8TBR4+CEx194PxNoeAqfk+bGargFUADQ2F33cUkt2Wa+K4OqZOaGDmqb3U1/QyjpNRsnYsOMttIGW7jTlJiqeWKArtT929xczWEV2inefubmZ5E0ZCctugsOy2XH1xVK9/YiMPD8hxg+Gz3HINmusGynYbg5IMGLQCrXEEFUQxVNcCbWY2DSD+83CyJoZTjpuMpCS5bYeAD8ysbwiqL/RwG1HYIZQ49LDYXLWkOW5Jjy+jW6LcNjP7HeBbQD3wLtHzeTLA88CVwH7gDnc/MtR+lNsmaTVUbluioWp3fx3It+PFSfYrUg2qb3qOSEqoeEQCqXhEAql4RAKpeEQCqXhEAql4RAKpeEQCqXhEAql4RAKpeEQCqXhEAql4RAJVV/GMVLDgSAUhFkvBiVWluopnpIIFRyoIsVgKTqwqiUMPzawG2AUccPdbzWwW8BxwGVEoyMo4WSexI52dTA4NPIxtODiH+U1HWTQh/ip2touuPds55zVMKGEQ4hOt81gzY0/i/Uh6laLneRB4K+f914En3H0O0Al8oQTHAKAzW594H/ObjrJ633XnwxA/OHCA49la6jIWFIQ4mHUffqwk+5H0StTzmNkM4A+AvwIeMjMDPgWsiDfZBHwVeDLJcXJ9tr2Y7Jz8rpjSzd1vf5x5zV1s7HqYbKaRTq/hmDdS09HOo7aME5mJCY9ypOBwxvMUnFhVkvY8fwf8GdAbv78MOOru2fh9KzA93wfNbJWZ7TKzXd3d3QmbUZy+LLebT26j3hwyNQD0Wg0ZelhyaltZ2yPVKUli6K3AYXffbWY3Ffv5kNy2zT03s+X+G4o91EX6Iqk+c3kbmaNZptR1UV9j8Vpj+bR2lt+V7DhDZrwNRsGJVSXJZduNwG1mdgvQCEwgyq6eZGa1ce8zAziQvJmRTT1L+FrCfeRmuV01+wcX3i/vH4ZYEQpOrCpJctsedfcZ7j4TuBN42d3vAn4M3B5vVtHctnxKneU2GGW8jX6JctvO7yS6bHs4Hqq+imioejLwGvA5dz871OeV2yZpNWK5bX3cfQfR0xBw93eB60uxX5E0q64ZBiIpouIRCaTiEQmk4hEJpOIRCaTiEQmk4hEJpOIRCaTiEQmk4hEJpOIRCaTiEQmk4hEJlM7iSVt+WVpy3Sohbb+LFEln8aQtvywtuW6VkLbfRYokyTD4KLAZmAo4sNHd15nZZGALMBN4D7jD3TuL2vmJNojz2cqdf1auXLdSUDZcZSXpebLAl939amAh8EUzuxpYC2x397nA9vh9sHLnn5Ur160UlA1XWcE9j7sfBA7Gr0+Y2VtEMVNLgZvizTYRfcP0kaJ23jz1QubZqy+WJKutGOXJdSuFgGy4YilLblAluecxs5nANUALMDUuLIBDRJd1+T5Tsdy24SjXTQpRiqzqS4B/Ab7k7sej0NCIu7uZ5U0YGTK3bUB+WSmy2opRjly3UgjKhiuWsuQGlTRut46ocP7Z3V+IF7eZ2TR3P2hm04DDRe+4gvllqc51qwRlyQ0q+LItzqV+GnjL3f82Z9U2orw2KEFuW7nzz8qV61YKyoarrODcNjP7OPBT4JdcyKp+jOi+53ngSmA/0VD1kaH2pdw2SasRyW1z91cAG2T14tD9ilSLdM4wEKkC6Sme4x9WugWFG8tz3UajwPl76SmeEweH3yYtxvJct9EocP5eSbKqSyGbzZ6fz5ZG5+e8jW+FvS8B0PXKN3nuzXPcM72Kes0ENJeuv9T0PN0+2NhDOvTNefvgwAHwXs54PceztSzOvlLpppWN5tL1l5qeB0rzvNGRNHdyB5mjv6LVJnLOMzRmnMYT+7mv7eGUzHUbaWWYS1cJgfP3UtPzVIPlvd+n3pyzvRlqMwYZzXUby1LT87T5ZFrKPIetWJ3/8Di9J7N8pMHp6u6hua6O+pr0zHUbaWWZS1cJgfP3UtPzHPJLK92EIe18p4PFh1azZ+VrNP35e+xZ+RoLu59i57IWuOt7lW6eJBE4fy81xTPj0vGVbsKQqmnO20jRXLr+SvJM0qQ0t03Saqi5banpeUSqTSp6HjNrB04BHZVuSxlMYfSf52g6x99098vzrUhF8QCY2a7BusfRZCyc51g4R9Blm0gwFY9IoDQVzyic95HXWDjPsXCO6bnnEak2aep5RKqKikckUNmLx8w+bWZ7zGyfmV2UY21mDWa2JV7fEqeRVpUCzvEeM2s3s9fjn/sq0c4kzOwZMztsZm8Ost7M7O/jv4M3zOzacrdxxLl72X6AGuAd4CqgHvgFcPWAbf4E2BC/vhPYUs42lukc7wHWV7qtCc/zk8C1wJuDrL8F+AFRwtJCoKXSbS71T7l7nuuBfe7+rrufA54jCobPtZQoIB5gK7DYcjN806+Qc6x67v4TYKg8vqXAZo/8DJgUJ8iOGuUununABznvW+Nlebdx9yxwDLisLK0rjULOEeAP48uZrfGzjkabQv8eqpYGDCrj34GZ7j4feIkLPa1UkXIXzwEg91/ZGfGyvNuYWS0wEaimgLRhz9Hdf+3uZ+O33wJG40NwCvldV7VyF89/A3PNbJaZ1RMNCAwMAMgNir8deNnjO9AqMew5Drj2vw14q4ztK5dtwN3xqNtC4JhfeG7TqFDWDAN3z5rZauCHRKNSz7j7/5jZXwK73H0b0ZMXvm1m+4huSO8sZxuTKvAcHzCz24geTXmEaPStqpjZs0RPAJxiZq3AV4A6AHffAHyfaMRtH3AaSHc0UgBNzxEJpAEDkUAqHpFAKh6RQCoekUAqHpFAKh6RQCoekUD/D7luuUGXDwGXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 216x162 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM8AAACZCAYAAACIVm61AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAPgklEQVR4nO3de3Bc5XnH8e+zq4tBSHawQBUYEL7EGZpxa6yCcZhYxSlQQkLHYwM1NYmHDNDGhZhkahOmk0x6G9oGx8Ezvkxo47RJbOKSVlMPU4jB0zQmdmVMnMYgW1ATS5ZtqZZ8l6WVnv5xjuyVvJJW7zm7e470fGZ22HPZc94X7TPv2eN3fyuqijFm9BKFboAxcWXFY4wjKx5jHFnxGOPIiscYR0WFbgBAZWWl1tTU5P28jY2NAMycOTPv5zbxsGfPnnZVvSbTtkgUT01NDQ0NDXk/b11dHQA7duzI+7lNPIjIh0Nts8s2YxyNz+LZtcHtdWfb4UfL4Oz/hdsecG+TKZjxWTy7N7q9btcGOPga7M7BG921TaZgRvzMIyL/ANwPHFfVj/vrrga2ADXAIeBBVe0QEQHWAPcB54DPq+rbuWl6AKePsfrb34T2A97yK49ftsv61unMKutkXkW7tyLVRVfjdro1ScXOF6GtEYpKhz3N6uaZrJjSGHbrTURkM/J8F7h30LpVwHZVnQFs95cBfh+Y4T8eB9aF08zwrTnysWG3zyrrZHlTLTtPVQJwuKWFU6kiihMC2gftBwOfw8TbiCOPqv6niNQMWv0AUOc/3wTsAFb667+n3mzTn4vIJBGpVtXWsBocivIqOAX7u38DgIfalmXc7drKHh49cCczy7vY2PUVUokJdGiSkzqBZHsbz8pCTicmDnOiE7Awy8uxF+eMshOm0Fw/81SlFcRRoMp/fj1wOG2/Zn/dZUTkcRFpEJGGtrY2x2bk1sQriqmqKOXuM/WUiEIiCUCfJEnQyz1n6wvcQlNIgf+dR1VVREb9vQZV3QhsBKitrc3v9yJuexx+DLdcVwHAlifuyLjbzvfbWf6DvXzmmmMkOlNUFndRkhR/q7C4uo3Fj2R+LUDNqm2ja5OJFdfiOdZ/OSYi1cBxf30LcEPaflP8ddFy+xPw4+Hf2P2Fs3bJbKZOe/XS8uLZzJtWmZs2mVhxvWyrBz7nP/8c8G9p6x8Vz1zgZOQ+7/ieXjBj2O37mk+ydsmlQpk3rZK1S2azr/lkaOcw8SYjfZNURH6Id3OgEjgGfA34V+Bl4EbgQ7xb1Sf8W9Vr8e7OnQOWqeqI825qa2vVpueYKBKRPapam2lbNnfb/nCITQsy7KvAF0fXPGPiaXzOMDAmBFY8xjiy4jHGkRWPMY6seIxxZMVjjCMrHmMcWfEY48iKxxhHVjzGOLLiMcaRFY8xjqx4jHE0/oonrHy0XGa4ZWK5bpEz/oonrHy0XGa4ZWK5bpETKMNARFYAXwAU+CWwDKgGNgOTgT3AUlXtDtjO8Jw+dimnbZjctmGluuDg697zLDPczNjjPPKIyPXAU0CtH4aYBB4GngdWq+p0oAN4LIyGhm11c3a/jLC+dfrF7LZ+h1taONObAElmneHmcm4TbUHTc4qAK0SkB7gSaAXuApb42zcBXydK4YflVbBwI2tWbaOm8qPeumGy1WalBYHMm1bJ7l81csO++0hedTWUlkBfL1w4Bff8DZRNzqoJa1ZtY8VTXx5duy3XLXKci0dVW0Tk74FfA+eB1/Au0zpVNeXvNmRuWxTsP3IKgIc2vDXsfteWl/LoS7upqijl4TP/zNTiPjrO99Jx/jwApdrFG+v+nK3lS3PeZhMdQS7bPoKXEHozcB1QxuWxvMO9vjChhw75aP3hhy2dXcwpOUSx9FKmZy4+ikgxrftADhqbxnLdIifIZdungP9V1TYAEXkF+AQwSUSK/NFnyNy2goUepuWjjRR62K8/s+2pu6azfNdXB0RS9bsVL/k+G6MKQ+xnuW6RE+RW9a+BuSJypR85tQDYD7wJLPL3Sc90i6X08MNn7p7J2iWzWf6Dvex8v73QTTMF5lw8qroL2Aq8jXebOoE3kqwEnhGRJrzb1S+F0M7QZRtIGEb4oeu5TbSNGHqYDxZ6aKJquNDD8TfDwJiQWPEY48iKxxhHVjzGOLLiMcaRFY8xjqx4jHFkxWOMIyseYxxZ8RjjyIrHGEdWPMY4suIxxtH4KZ5c5p7lO8NtJJbxlhfjp3hymXuW7wy3kVjGW14ETc+Jj/S8tn6uuW3pLMNt3Ao08ojIJBHZKiLvici7InKHiFwtIq+LyEH/vx8Jq7Fhcc1Ny0WGW1htMwWgqs4PvFy2L/jPS4BJwN8Cq/x1q4DnRzrOnDlzNOe+fevFpzet/HdVVZ0/f77Onz8/60P8rKlNZ3/jNf1ZU5uqqu76n/f0yNem6rm/+7h3/G/9luoLv6l6pt25mf1tCyStryYYoEGHeN86X7aJyETgk8Dn/SLsBrpF5AGgzt9tE7ADL9cgUh7a8FbWuW3pLMPN9Aty2XYz0Ab8o4jsFZHviEgZUKWqrf4+R4GqTC/Oe25bSLlnkchwG4llvOXHUEPSSA+gFkgBt/vLa4C/wEsMTd+vY6Rj5eWyLY3rZZvqpUu3b/7HewMu4cJum4kGhrlsCzLyNAPN6kVQgRdDdStwTESqAfz/Hg9wjkixDDeTLkhu21HgsIj03x7qDz2sxws7hIiGHrrmpuUiwy2stpn8C5TbJiK/DXwH707bB3i/z5MAXgZuBD4EHlTVE8Mdx3LbTFQNl9sW6B9JVfUdvM8+gy0Iclxj4mD8TM8xJmRWPMY4suIxxpEVjzGOrHiMcWTFY4wjKx5jHFnxGOPIiscYR1Y8xjiy4jHGkRWPMY6seIxxNHaLpxDBf1ELPxwtC0sclbFbPIUI/ota+OFoWVjiqAQOPRSRJNAAtKjq/SJyM7AZmAzsAZaql6yTX4NCDlc3z2TFlMaB+wQIPVzfOp1ZZZ3Mq/C/gp3qoqtxO92apCJP4YcZ+2TyJoyR52ng3bTl54HVqjod6AAeC+Ecga058rFQjzerrJPlTbUXQxAPt7RwKlVEcUJCCT/MRth9MqMTaOQRkSnAp4G/Ap4REQHuApb4u2wCvg6sC3IeJ+VVsDDtMmT3Nh5qWzZgl/3d3le/B6/P1rWVPTx64E5mlnexsesrpBIT6NAkJ3UCyfY2npWFnE5MdO7CyE4M7GNQL84J71jjQNCR51vAnwF9/vJkvOiplL/cDFyf6YV5z23Lgf4Mt7vP1FMiCokkAH2SJEEv95ytL3ALTS4FSQy9HziuqntEpG60r1fVjcBG8AJAXNsxpAzBf1ueuGPAct0PKzKuz1Z/FNVnrjlGojNFZXEXJUnxtwqLq9tY/IjbsbNRs2pbuAe0sMRRCXLZ9gngsyJyHzABqMALPpwkIkX+6DMFaAneTAe3P5HTw6dnuE2d9uql5cWXoqliJ8f/z8aaILltz6rqFFWtAR4G3lDVR4A3gUX+bpHJbQs7Dy0fGW4jsYy3wsrF7/OsBDaLyF8Ce4GXcnCOUVvxex8N9XhPzp922bp50yrzOuqE3SczOqEUj6ruwPs1BFT1A+C2MI5rTJSN3RkGxuSYFY8xjqx4jHFkxWOMIyseYxxZ8RjjyIrHGEdWPMY4suIxxpEVjzGOrHiMcWTFY4wjKx5jHMW3eKKeMRb3DLdCiPrfdJD4Fk/UM8binuFWCFH/mw4SJMPgBuB7QBWgwEZVXSMiVwNbgBrgEPCgqnYEb+ogg3LZwCHHLEBu22ADctxSXXDwdbr6kvT8dB3lechwyyXLh8ssyMiTAr6sqrcAc4EvisgtwCpgu6rOALb7y3lRyByzATlu7Qe50AdHU2WUSCovGW65ZPlwmTmPPKraCrT6z0+LyLt4MVMPAHX+bpvwvmG6MlArMxmcywYZs9mGEzS3bbBrK3v40oFZ1JdsoUuvoiiZpEMr8pThlksh58MNJWa5caF85hGRGmA2sAuo8gsL4CjeZV2m18Q+t22wiVcU8+QVP0H7ekkkikiIZbiNZWFkVV8F/AvwJVU95YWGelRVRSRjJlvg3LYhMsZGk8EWNLdtsJ3vt9P3T01UlEJXzxnKi4v9HLfcZ7jlUuj5cEOJWW5c0LjdYrzC+b6qvuKvPiYi1araKiLVwPGgjcwoYhljF3Pblm6lbFolvxgLOW75FrG/6UicL9v8XOqXgHdV9YW0TfV4eW2Q59y2QuaYRSHHLVcsHy4zUXVLuhWRO4GfAr/kUlb1V/E+97wM3Ah8iHer+sRwx6qtrdWGhgandgRRV1cHwI4dO/J+bhMPIrJHVWszbQtyt+2/ABli8wLX4xoTF/GdYWBMgcWjeGI252kAm+MWDw7vsXgUT8zmPA1gc9ziweE9loug91Ctfv0AKzLMYwtFiHPbBlvfOp1bJ7Ry29HN3oqdL7L7gzbe7qrmyeqm0M83XkRpnl3kR5412+M5L2xWWSfvHGrnQh8gSS70wd5D7cwq6yx002ItSvPsIj/yABzprWBFSPPP0oU9ty1deW8nf130p7SmriLRV0RfX4o/KHqL584/wpoLcZ3jFgU5mmfnMK8u8iNPXN17rp6k9JFIFNHTqyQSRRRJn81xG0NiMfJcd/fTbLk9/HlhYc9tG+D7L9B9GOg6w4TSJF09vZSXFsd6jlsU5GyencO8uniMPDGb8wSwc+465vZsoHHpXsqeO0Tj0r3M7dnAzrnrCt00k4nDeyzyxRPXeVVjea5bIUXp/eA8ty1MNrfNRNVwc9siP/IYE1WRGHlEpA1vBnYhVALtBTp32MZSXyAa/blJVa/JtCESxVNIItIw1LAcN2OpLxD9/thlmzGOrHiMcWTF44eQjBFjqS8Q8f6M+888xriykccYR1Y8xjgaN8UjIveKSKOINInIZfnZIlIqIlv87bv8FNRIyqIvnxSRt0UkJSKLCtHGbGXRl2dEZL+I7BOR7SJyUyHamZGqjvkHkATeB6YCJcAvgFsG7fMnwHr/+cPAlkK3O0BfaoBZeL9isajQbQ7Yl98FrvSf/3GU/i7jZeS5DWhS1Q9UtRvYjBdIn+4BvGB6gK3AAknPDo6OEfuiqodUdR+X8vSiKpu+vKmq5/zFnwNT8tzGIY2X4rkeOJy23Oyvy7iPqqaAk8DkvLRudLLpS1yMti+PAa/mtEWjEIsvwxkjIn8E1ALzC92WfuNl5GkBbkhbnuKvy7iPiBQBE4Eohq1l05e4yKovIvIp4Dngs6p6IU9tG9F4KZ7/BmaIyM0iUoJ3Q2BwmEB6QP0i4A31P6VGTDZ9iYsR+yIis4ENeIWTm1/ccFXoOxZ5vLNzH3AA7+7Oc/66b+D9UQAmAD8CmoDdwNRCtzlAX34H7/PDWbzR81eFbnOAvvwEOAa84z/qC93m/odNzzHG0Xi5bDMmdFY8xjiy4jHGkRWPMY6seIxxZMVjjCMrHmMc/T9biXWcru57OgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 216x162 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for params in params_list:\n",
    "\n",
    "    result = (\n",
    "        get_normal_t_improvement_against_M(\n",
    "            N=params['N'], mu_X=params['mu_X'], sigma_X=params['sigma_X'],\n",
    "            mu_epsilon=params['mu_eps'], sigma_1=params['sigma_1'], \n",
    "            sigma_2=params['sigma_2'], M_list=params['M'], verbose=True)\n",
    "    )\n",
    "\n",
    "    fig, ax = generate_normal_t_improvement_against_M_plot(params['M'], result)\n",
    "    \n",
    "    fig.set_size_inches(3, 2.25)\n",
    "    fig.savefig(\"./output/fig/t_case_study_{}_{}_{}.pdf\"\n",
    "                .format(params['name'],\n",
    "                        str(params['sigma_1']).replace('.', '-'),\n",
    "                        str(params['sigma_2']).replace('.', '-')),\n",
    "                transparent=True, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# General runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_run_setting(count, N, M, mu_X, sigma_X, mu_eps, sigma_1, sigma_2):\n",
    "     print(\"Cycle {}: N = {}, M = {}, mu_X = {}, sigma_X = {}, \"\n",
    "           \"mu_epsilon = {}, sigma_1 = {}, sigma_2 = {}\"\n",
    "           .format(count, N, M, np.round(mu_X, 4), \n",
    "                   np.round(sigma_X, 4), np.round(mu_eps, 4),\n",
    "                   np.round(sigma_1, 4), np.round(sigma_2, 4)))\n",
    " \n",
    "\n",
    "def calculate_norm_and_t_improvements(cycles=100):\n",
    "    # Constants\n",
    "    NUM_RUNS = 1000\n",
    "    T_DF = 3\n",
    "    CI_HIGH_PERCENTILE = 95\n",
    "    CI_LOW_PERCENTILE = 5  \n",
    "    \n",
    "    # Initialise empty list to hold the results\n",
    "    normal_E_Ds = []\n",
    "    normal_CI_lows = []\n",
    "    normal_CI_highs = []\n",
    "    t_E_Ds = []\n",
    "    t_CI_lows = []\n",
    "    t_CI_highs = []\n",
    "    \n",
    "    for cycle in range(1, cycles + 1, 1):\n",
    "        params = get_test_params()\n",
    "        N = int(params['N'])\n",
    "        M = int(params['M'])\n",
    "        mu_X = params['mu_X']\n",
    "        mu_epsilon = params['mu_epsilon']\n",
    "        sigma_sq_X = params['sigma_sq_X']\n",
    "        sigma_sq_1 = params['sigma_sq_1']\n",
    "        sigma_sq_2 = params['sigma_sq_2']\n",
    "\n",
    "        # Reconciling the use of sigmas by numpy\n",
    "        # and sigma_sqs in the theoretical calculations\n",
    "        sigma_X = np.sqrt(sigma_sq_X)\n",
    "        sigma_1 = np.sqrt(sigma_sq_1)\n",
    "        sigma_2 = np.sqrt(sigma_sq_2)\n",
    "    \n",
    "        print_run_setting(cycle, N, M, mu_X, sigma_X,\n",
    "                          mu_epsilon, sigma_1, sigma_2)\n",
    "    \n",
    "        # Normal assumptions\n",
    "        normal_improvement = (\n",
    "            nnm.get_samples(\n",
    "                n_samples=NUM_RUNS, N=N, M=M, \n",
    "                mu_X=mu_X, mu_epsilon=mu_epsilon, \n",
    "                sigma_sq_X=sigma_sq_X, sigma_sq_1=sigma_sq_1,\n",
    "                sigma_sq_2=sigma_sq_2, verbose=False)['improvement']\n",
    "        )\n",
    "\n",
    "        normal_E_D = np.mean(normal_improvement)\n",
    "        normal_CI_low = np.percentile(normal_improvement, CI_LOW_PERCENTILE)\n",
    "        normal_CI_high = np.percentile(normal_improvement, CI_HIGH_PERCENTILE)\n",
    "\n",
    "        normal_E_Ds.append(normal_E_D)\n",
    "        normal_CI_lows.append(normal_CI_low)\n",
    "        normal_CI_highs.append(normal_CI_high)\n",
    "        \n",
    "        # t-assumptions\n",
    "        t_improvement = (\n",
    "            ttm.get_samples(\n",
    "                n_samples=NUM_RUNS, N=N, M=M, \n",
    "                mu_X=mu_X, mu_epsilon=mu_epsilon, \n",
    "                sigma_sq_X=sigma_sq_X, sigma_sq_1=sigma_sq_1,\n",
    "                sigma_sq_2=sigma_sq_2, nu=T_DF, verbose=False)['improvement']\n",
    "        )\n",
    "        \n",
    "        t_E_D = np.mean(t_improvement)\n",
    "        t_CI_low = np.percentile(t_improvement, 5)\n",
    "        t_CI_high = np.percentile(t_improvement, 95)\n",
    "\n",
    "        t_E_Ds.append(t_E_D)\n",
    "        t_CI_lows.append(t_CI_low)\n",
    "        t_CI_highs.append(t_CI_high)\n",
    "    \n",
    "        print(\"Cycle {}: Normal E(D) = {}; t E(D) = {}\"\n",
    "             .format(cycle, normal_E_D, t_E_D))\n",
    "        \n",
    "    # Prepare result in a dict\n",
    "    result = {\n",
    "        'normal_E_Ds': normal_E_Ds,\n",
    "        'normal_CI_lows': normal_CI_lows,\n",
    "        'normal_CI_highs': normal_CI_highs,\n",
    "        't_E_Ds': t_E_Ds,\n",
    "        't_CI_lows': t_CI_lows,\n",
    "        't_CI_highs': t_CI_highs,\n",
    "    }\n",
    "            \n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cycle 1: N = 743, M = 329, mu_X = 3.7665, sigma_X = 9.221, mu_epsilon = -0.0438, sigma_1 = 8.863, sigma_2 = 2.9277\n",
      "Cycle 1: Normal E(D) = 1.8975043738303834; t E(D) = 0.696125145049461\n",
      "Cycle 2: N = 2576, M = 830, mu_X = 1.2699, sigma_X = 2.1169, mu_epsilon = -4.2157, sigma_1 = 2.0734, sigma_2 = 0.5917\n",
      "Cycle 2: Normal E(D) = 0.5857895746468357; t E(D) = 0.2244220866944466\n",
      "Cycle 3: N = 558, M = 382, mu_X = -4.3298, sigma_X = 5.2404, mu_epsilon = 5.2075, sigma_1 = 5.5045, sigma_2 = 5.3395\n",
      "Cycle 3: Normal E(D) = 0.030380646754874355; t E(D) = 0.008654740933419353\n",
      "Cycle 4: N = 2515, M = 1366, mu_X = -3.8029, sigma_X = 4.7719, mu_epsilon = 8.8466, sigma_1 = 8.294, sigma_2 = 7.5481\n",
      "Cycle 4: Normal E(D) = 0.11917893642517761; t E(D) = 0.044937624633430574\n",
      "Cycle 5: N = 73, M = 37, mu_X = -1.5309, sigma_X = 3.0347, mu_epsilon = -0.5342, sigma_1 = 5.5688, sigma_2 = 4.7152\n",
      "Cycle 5: Normal E(D) = 0.1466374567223039; t E(D) = 0.06419475746459284\n",
      "Cycle 6: N = 18, M = 5, mu_X = 5.8109, sigma_X = 6.5851, mu_epsilon = -7.1895, sigma_1 = 5.403, sigma_2 = 4.0285\n",
      "Cycle 6: Normal E(D) = 0.6283847114051593; t E(D) = 0.19595238046646699\n",
      "Cycle 7: N = 3002, M = 1097, mu_X = 4.7965, sigma_X = 9.5829, mu_epsilon = -2.4087, sigma_1 = 6.501, sigma_2 = 3.6759\n",
      "Cycle 7: Normal E(D) = 1.0481782855579986; t E(D) = 0.40322032109700273\n",
      "Cycle 8: N = 74, M = 41, mu_X = 0.477, sigma_X = 0.4327, mu_epsilon = 5.9466, sigma_1 = 4.8978, sigma_2 = 1.9317\n",
      "Cycle 8: Normal E(D) = 0.039081474760881744; t E(D) = 0.02535315344571938\n",
      "Cycle 9: N = 811, M = 515, mu_X = -1.4213, sigma_X = 5.4496, mu_epsilon = -7.8215, sigma_1 = 6.0916, sigma_2 = 2.6024\n",
      "Cycle 9: Normal E(D) = 0.7601320850848359; t E(D) = 0.26489588393905356\n",
      "Cycle 10: N = 3025, M = 1412, mu_X = 7.3243, sigma_X = 7.5958, mu_epsilon = -7.3031, sigma_1 = 1.4666, sigma_2 = 0.9571\n",
      "Cycle 10: Normal E(D) = 0.06622997105372495; t E(D) = 0.040066118286449445\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'normal_E_Ds': [1.8975043738303834,\n",
       "  0.5857895746468357,\n",
       "  0.030380646754874355,\n",
       "  0.11917893642517761,\n",
       "  0.1466374567223039,\n",
       "  0.6283847114051593,\n",
       "  1.0481782855579986,\n",
       "  0.039081474760881744,\n",
       "  0.7601320850848359,\n",
       "  0.06622997105372495],\n",
       " 'normal_CI_lows': [1.554117754211027,\n",
       "  0.5298329567107,\n",
       "  -0.16024247684374232,\n",
       "  -0.032312773304483705,\n",
       "  -0.4696363601393994,\n",
       "  -1.760415791234639,\n",
       "  0.8917353550097207,\n",
       "  -0.06414662529003909,\n",
       "  0.6087667813758106,\n",
       "  0.04708183768946075],\n",
       " 'normal_CI_highs': [2.233280704687485,\n",
       "  0.639956496556603,\n",
       "  0.22335257300775102,\n",
       "  0.2692829857359111,\n",
       "  0.7766349608135799,\n",
       "  3.3206495632028967,\n",
       "  1.2078238967546737,\n",
       "  0.14388196509369147,\n",
       "  0.9118124909525646,\n",
       "  0.08670451321290215],\n",
       " 't_E_Ds': [0.696125145049461,\n",
       "  0.2244220866944466,\n",
       "  0.008654740933419353,\n",
       "  0.044937624633430574,\n",
       "  0.06419475746459284,\n",
       "  0.19595238046646699,\n",
       "  0.40322032109700273,\n",
       "  0.02535315344571938,\n",
       "  0.26489588393905356,\n",
       "  0.040066118286449445],\n",
       " 't_CI_lows': [0.5444811133134807,\n",
       "  0.19928961307274634,\n",
       "  -0.08712052410933795,\n",
       "  -0.01859980882732748,\n",
       "  -0.1909419133231386,\n",
       "  -0.9106605204285834,\n",
       "  0.32492230856235826,\n",
       "  -0.025687620288312573,\n",
       "  0.19891169552666424,\n",
       "  0.02530376824108407],\n",
       " 't_CI_highs': [0.8546100585570586,\n",
       "  0.24928362963168985,\n",
       "  0.09618357344935632,\n",
       "  0.1094187199509189,\n",
       "  0.31772775760422534,\n",
       "  1.402083201760571,\n",
       "  0.48665784310366184,\n",
       "  0.08420068219401601,\n",
       "  0.3308180081321436,\n",
       "  0.05532707711032829]}"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calculate_norm_and_t_improvements(cycles=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}