DS_L4_22_Apr_Assignment_ipynb_txt.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "DS_L4_22_Apr_Assignment.ipynb.txt",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NVMUAJt8gzbA",
        "colab_type": "text"
      },
      "source": [
        "# MLE and Hypothesis Testing"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "d1-kEc8hXx3N",
        "colab_type": "text"
      },
      "source": [
        "### Read the given csv file [```q4.csv```](https://raw.githubusercontent.com/Ipsit1234/SHALA/master/q4.csv) of samples. Assume that it is drawn from a exponential distribution. What is the difference of log likelihoods for $\\lambda = 1$ and $\\lambda = 2$ ? Which of the two is more strongly supported?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DHfUkvD3KcNN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "from scipy import stats\n",
        "from matplotlib import pyplot as plt\n",
        "import seaborn as sns\n",
        "import random\n",
        "from numpy.random import normal\n",
        "from scipy.stats import norm\n",
        "from matplotlib import pyplot"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UAphR126LEwa",
        "colab_type": "code",
        "outputId": "7d6ceb94-7ceb-460e-fdee-85276b4bf40f",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 204
        }
      },
      "source": [
        "col_Names=[\"Sequence\", \"Coverage\"]\n",
        "url = 'https://raw.githubusercontent.com/Ipsit1234/SHALA/master/q4.csv'\n",
        "data = pd.read_csv(url,names=col_Names)\n",
        "data.head(5)"
      ],
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>Sequence</th>\n",
              "      <th>Coverage</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>0</td>\n",
              "      <td>0.019404</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>1</td>\n",
              "      <td>0.043634</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>2</td>\n",
              "      <td>0.023797</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>3</td>\n",
              "      <td>0.444340</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>4</td>\n",
              "      <td>0.578160</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "   Sequence  Coverage\n",
              "0         0  0.019404\n",
              "1         1  0.043634\n",
              "2         2  0.023797\n",
              "3         3  0.444340\n",
              "4         4  0.578160"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 9
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nDpo9jpidQHk",
        "colab_type": "text"
      },
      "source": [
        "# PDF of exponential distribution is $\\lambda e^{-\\lambda x}$\n",
        "Define a function to find pdf. You can use any module to get pre-defined functions also. And then find the likelihood and log-likelihood"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vdFipQwCaE9Y",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "47f7eaae-0fcc-465b-99c3-7a7367239f64"
      },
      "source": [
        "def pdf_exp(lam,x):\n",
        "  return lam*np.exp(-lam*x)\n",
        "\n",
        "likelihood1=1\n",
        "likelihood2=1\n",
        "\n",
        "for a in data['Coverage'].values[:]:\n",
        "  likelihood1 = likelihood1 * pdf_exp(lam=2,x=a)\n",
        "  likelihood2 = likelihood2 * pdf_exp(lam=1,x=a)\n",
        "log_likelihood1=np.log(likelihood1)\n",
        "log_likelihood2=np.log(likelihood2)\n",
        "print('log_likelihood1-log_likelihood2:', log_likelihood1-log_likelihood2)"
      ],
      "execution_count": 26,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "log_likelihood1-log_likelihood2: 20.366416970604874\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "p0IVzZ_IPyLG",
        "colab_type": "code",
        "outputId": "f9c07eec-1518-4bdd-fe73-6e4d3441a737",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 298
        }
      },
      "source": [
        "from scipy.stats import expon\n",
        "# data_expon_1 = expon.rvs(scale=1,loc=0,size=1000)\n",
        "ax = sns.distplot(data['Coverage'].values,\n",
        "                  kde=True,\n",
        "                  bins=30,\n",
        "                  color='grey',\n",
        "                  hist_kws={\"linewidth\": 15,'alpha':1})\n",
        "ax.set(xlabel='Exponential Distribution', ylabel='Frequency')"
      ],
      "execution_count": 31,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[Text(0, 0.5, 'Frequency'), Text(0.5, 0, 'Exponential Distribution')]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 31
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3hV9Z3v8feHACFAAIEAykVA5aqCGoLXDq3a0epox3rtmVpHW6bt2Bnn9rRnzpxSOmfOM52eMzP11NZax1FrvdTLWNqC1FpbVFAuQUBQkHsCCqhIICCB8D1/7A2NMZcNZGXtZH9ez7Mf91rrt9b+rAT3N2v91votRQRmZla4uqQdwMzM0uVCYGZW4FwIzMwKnAuBmVmBcyEwMytwLgRmZgUusUIg6T5J2yW91kKbaZJelbRS0u+SymJmZs1TUvcRSPoYsAd4MCJOb2J5P2A+cFlEbJY0KCK2t7bdgQMHxsiRI9s8r5lZZ7ZkyZJ3IqKsqWVdk/rQiJgnaWQLTT4LPBURm7PtWy0CACNHjmTx4sXHH9DMrIBI2tTcsjT7CMYAJ0j6raQlkm5urqGk6ZIWS1q8Y8eOdoxoZtb5pVkIugLnAFcAfwj8T0ljmmoYEfdERHlElJeVNXlkY2ZmxyixU0M5qAbejYhaoFbSPGASsCbFTGZmBSfNI4KfARdK6iqpJzAVeD3FPGZmBSmxIwJJjwDTgIGSqoEZQDeAiLg7Il6X9AywHDgE3BsRzV5qamZmyUjyqqGbcmjzHeA7SWUwM7PW+c5iM7MC50JgZlbgXAjMzApcmpePWkJmzpzZbp81Y8aMdvssM0uGjwjMzAqcC4GZWYFzITAzK3AuBGZmBc6FwMyswLkQmJkVOBcCM7MC50JgZlbgXAjMzAqcC4GZWYFzITAzK3AuBGZmBc6FwMyswLkQmJkVuMQKgaT7JG2X1OJziCVNkXRQ0rVJZTEzs+YleURwP3BZSw0kFQHfBn6VYA4zM2tBYoUgIuYB77XS7KvAk8D2pHKYmVnLUusjkDQU+GPgBzm0nS5psaTFO3bsSD6cmVkBSbOz+N+Br0XEodYaRsQ9EVEeEeVlZWXtEM3MrHCk+czicuBRSQADgU9JOhgRT6eYycys4KRWCCJi1OH3ku4HfuEiYGbW/hIrBJIeAaYBAyVVAzOAbgARcXdSn2tmZkdHEZF2hqNSXl4eixcvTjtG6mbOnJl2hFbNmDEj7QhmliVpSUSUN7XMdxabmRU4FwIzswLnQmBmVuBcCMzMCpwLgZlZgXMhMDMrcC4EZmYFzoXAzKzAuRCYmRU4FwIzswLnQmBmVuBcCMzMCpwLgZlZgXMhMDMrcC4EZmYFzoXAzKzAuRCYmRU4FwIzswKXWCGQdJ+k7ZJea2b5f5O0XNIKSfMlTUoqi5mZNS/JI4L7gctaWL4B+IOIOAP4R+CeBLOYmVkzuia14YiYJ2lkC8vnN5h8GRiWVBYzM2tevvQR3AbMaW6hpOmSFktavGPHjnaMZWbW+aVeCCR9nEwh+FpzbSLinogoj4jysrKy9gtnZlYAEjs1lAtJZwL3ApdHxLtpZjEzK1SpHRFIGgE8BXwuItaklcPMrNAldkQg6RFgGjBQUjUwA+gGEBF3A98ABgDflwRwMCLKk8pjZmZNS/KqoZtaWf4F4AtJfb6ZmeUm9c5iMzNLlwuBmVmBcyGwxERE2hHMLAcuBJaYTZs2pR3BzHLgQmCJWbhwYdoRzCwHLgSWmDfeeIP3338/7Rhm1goXAkvUokWL0o5gZq1wIbDEjB8/nsrKSg4cOJB2FDNrgQuBJWbq1Kl88MEHLF++PO0oZtYCFwJLzPDhwxkyZAgLFy70paRmecyFwBIjialTp7J9+3Y2btyYdhwza4YLgSXq9NNPp2fPnrzyyitpRzGzZrgQWKK6du3KOeecw+rVq9m5c2faccysCS4Elrjy8nIk+VJSszzlQmCJ69OnDxMmTKCyspK6urq045hZIy4E1i4qKirYv38/K1asSDuKmTXiQmDtwpeSmuUvFwJrF5KoqKhg+/btHpXULM8k+czi+4Arge0RcXoTywV8F/gUsBe4JSIqk8pj7W/mzJlNzn/ggQeYMWNGO6cxs+YkeURwP3BZC8svB07LvqYDP0gwi+WZmpqatCOYWVZihSAi5gHvtdDkauDByHgZ6CfpxKTyWH5ZvHhx2hHMLCvNPoKhQFWD6ersvI+QNF3SYkmLd+zY0S7hLFmVlZUcPHgw7RhmRgfpLI6IeyKiPCLKy8rK0o5jbaC2tpZVq1alHcPMSLcQbAGGN5gelp1nBWDAgAF+lKVZnsipEEg6I4HPngXcrIxzgV0R8VYCn2N5qKKigi1btrBli2u/WdpyPSL4vqSFkr4iqW8uK0h6BFgAjJVULek2SV+S9KVsk9nAemAt8CPgK0cb3jquSZMm0b17d48/ZJYHcrqPICIuknQacCuwRNJC4D8j4tkW1rmplW0G8OdHE9Y6j+LiYiZNmkRlZSWXXnopvXr1SjuSWcHKuY8gIt4E/gH4GvAHwJ2S3pB0TVLhrHObMmUK9fX1VFb6PkKzNOXaR3CmpH8DXgc+AfxRRIzPvv+3BPNZJ1ZWVsbo0aNZvHgxhw4dSjuOWcHK9Yjg/wGVwKSI+PPDQ0FExFYyRwlmx2TKlCnU1NSwevXqtKOYFaxcC8EVwMMRsQ9AUhdJPQEi4sdJhbPOb8yYMfTt29eXkpqlKNdC8GugpMF0z+w8s+PSpUsXysvL2bhxI9u2bUs7jllByrUQ9IiIPYcnsu97JhPJCs0555xDt27dWLBgQdpRzApSroWgVtLZhycknQPsSyaSFZqSkhLOOussVqxY4VFJzVKQayG4A3hc0guSXgQeA25PLpYVmvPOO4+I4OWXX047ilnByfWGskWSxgFjs7NWR8SB5GJZoenXrx8TJ05kyZIlfOxjH6NHjx5pRzIrGEcz6NwU4EzgbOAmSTcnE8kK1fnnn09dXR1LlixJO4pZQcn1hrIfA/8HuJBMQZgClCeYywrQiSeeyKhRo3jllVeor69PO45Zwcj1mcXlwITs+EBmibngggt46KGHWLFiBZMnT047jllByPXU0GvAkCSDmAGMHj2awYMHM3/+fPx3h1n7yLUQDARWSZoradbhV5LBrDBJ4vzzz2fHjh28+eabaccxKwi5nhr6ZpIhzA6bOXPmkfePPPJIzuvNmDEjiThmBSHXy0d/J+lk4LSI+HV2nKGiZKOZmVl7yPWqoS8CTwA/zM4aCjydVCgzM2s/ufYR/DlwAVADRx5SMyipUGZm1n5yLQT7I6Lu8ISkrkCrl3RIukzSaklrJX29ieUjJD0vaamk5ZI+lXt0MzNrC7kWgt9J+nugRNKlwOPAz1taQVIRcBdwOTCBzN3IExo1+wfgpxFxFnAj8P2jCW9mZscv10LwdWAHsAL4M2A2rT+ZrAJYGxHrs0cTjwJXN2oTQJ/s+77A1hzzmJlZG8n1qqFDwI+yr1wNBaoaTFcDUxu1+SbwK0lfBXoBlzS1IUnTgekAI0aMOIoIZmbWmlyvGtogaX3jVxt8/k3A/RExDPgU8GNJH8kUEfdERHlElJeVlbXBx5qZ2WFHM9bQYT2A64D+rayzBRjeYHpYdl5DtwGXAUTEAkk9yNzFvD3HXGZmdpxyOiKIiHcbvLZExL+TeaB9SxYBp0kaJak7mc7gxsNSbAYuBpA0nkyR2XFUe2BmZsclpyOCho+pJFM8yltbNyIOSrodmEvmLuT7ImKlpG8BiyNiFvA3wI8k/RWZjuNbPMKpmVn7yvXU0P9t8P4gsBG4vrWVImI2mSuMGs77RoP3q8jcqGZmZinJ9aqhjycdxMzM0pHrqaG/bml5RPxr28SxXOzatSvtCGbWiRzNVUNT+H1n7x8BCwEPGJ+C6urqtCOYWSeSayEYBpwdEbsBJH0T+GVE/ElSwax5VVVVrTcyM8tRrkNMDAbqGkzXZedZCnxEYGZtKdcjggeBhZL+Kzv9aeCBZCJZSw4cOMBbb72Vdgwz60RyvWronyTNAS7KzvrTiFiaXCxrztatWzl06FDaMfJORCAp7RhmHVKup4YAegI1EfFdoFrSqIQyWQt8WqhpK1euTDuCWYeV66BzM4CvAf89O6sb8FBSoax5VVVV9O/f2jBPhefZZ5/lwIEDaccw65ByPSL4Y+AqoBYgIrYCpUmFsqZFBNXV1QwfPrz1xgWmpqaGF198Me0YZh1Srp3FdRERkgJAUq8EM1kzdu7cSW1tLcOGDWPZsmVpx8k78+bNY968eR+ZP2PGjBTSmHUcuR4R/FTSD4F+kr4I/Jqje0iNtYHD/QM+IjCzttTqEYEyl2I8BowDaoCxwDci4tmEs1kjVVVVdO/eHT+cx8zaUquFIHtKaHZEnAH4yz9FVVVVDBs2jC5djuZiLzOzluX6jVIpaUqiSaxF+/fvZ/v27QwbNiztKGbWyeTaWTwV+BNJG8lcOSQyBwtnJhXMPmzLli1EhPsHzKzNtVgIJI2IiM3AH7ZTHmvG4YHmhg4dmnISM+tsWjsieJrMqKObJD0ZEZ9pj1D2UVVVVQwaNIiSkpK0o5hZJ9NaH0HDwVtGH+3GJV0mabWktZK+3kyb6yWtkrRS0sNH+xmF4NChQ1RVVfm00DGqr69PO4JZXmutEEQz71slqQi4C7gcmADcJGlCozankRm24oKImAjccTSfUSi2b99OXV0dI0aMSDtKhzR//vy0I5jltdYKwSRJNZJ2A2dm39dI2i2pppV1K4C1EbE+IuqAR4GrG7X5InBXROwEiIjtx7ITnd3mzZsBXAiO0bx589i5c2faMczyVouFICKKIqJPRJRGRNfs+8PTfVrZ9lCg4aO0qrPzGhoDjJH0kqSXJV3W1IYkTZe0WNLiHTt2tLZPnU5VVRWlpaX07ds37SgdkiSeeeaZtGOY5a2070zqCpwGTANuAn4kqV/jRhFxT0SUR0R5Id5Vu3nzZkaMGOHx9o/RtGnTWLNmDatXr047illeSrIQbAEa9m4Oy85rqBqYFREHImIDsIZMYbCsXbt2UVNT447i4zB16lQGDRrEnDlzqKura30FswKTZCFYBJwmaZSk7sCNwKxGbZ4mczSApIFkThWtTzBTh+P+geNXVFTEFVdcwa5du5ocndSs0CVWCCLiIHA7MBd4HfhpRKyU9C1JV2WbzQXelbQKeB74u4h4N6lMHVFVVRXdunVj8ODBaUfp0EaMGMHkyZNZsGABhdjPZNaSRPsIImJ2RIyJiFMi4p+y874REbOy7yMi/joiJkTEGRHxaJJ5OiIPNNd2LrnkErp3784vfvELIo7qamizTs3fLnls//79bNu2zaeF2kivXr345Cc/yebNm1myZEnacczyhgtBHquurvZAc21s8uTJjBo1il//+tfU1LR2K4xZYXAhyGObN29GkoeebkOSuPLKK6mvr2f27Nk+RWSGC0Feq6qqYvDgwRQXF6cdpVPp378/06ZNY/Xq1bz++utpxzFLnQtBnqqvr6e6utr9Awk577zzGDJkCHPmzGHfvn1pxzFLlQtBntq2bRsHDhxw/0BCunTpwlVXXUVtbS3PPusnsFphy/UJZdbOOvuNZDNnzsybz1q6dClLly790LwZM2Yc0/aaW+9Y1ilUzf2s/HNKjo8I8lRVVRV9+/alT5/WxvYzMzs+LgR5KCKODDRnZpY0F4I89P7777Nnzx73D5hZu3AhyEMbN24E4OSTT043iJkVBBeCPLRp0yZ69uxJIT57IV/4RjMrJC4EeSYi2LBhAyNHjvSDaFK0fPnytCOYtRsXgjzz/vvvU1NTw8iRI9OOUtDmzJnD+++/n3YMs3bhQpBnNmzYAOBCkAcef/xx6uvr045hljgXgjyzadMmevXqxcCBA9OOUtCuvvpqtm7dyq9+9au0o5glzoUgj7h/IH+MHz+eqVOnsnDhQlatWpV2HLNEuRDkkZ07d7J7925fNponLr30UoYOHcrPfvYz3nvvvbTjmCUm0UIg6TJJqyWtlfT1Ftp9RlJIKk8yT747fP/AqFGj0g1iQOah99dddx1FRUU8/vjjHDhwIO1IZolIrBBIKgLuAi4HJgA3SZrQRLtS4C+BV5LK0lFs2LCB3r17M2DAgLSjWFbfvn359Kc/zdtvv80zzzyTdhyzRCR5RFABrI2I9RFRBzwKXN1Eu38Evg18kGCWvBcRrF+/ntGjR7t/IM+MGTOGCy+8kMrKShYtWpR2HLM2l2QhGApUNZiuzs47QtLZwPCI+GVLG5I0XdJiSYt37NjR9knzwNtvv83evXsZPXp02lGsCR//+McZM2YMc+bMSTuKWZtLrbNYUhfgX4G/aa1tRNwTEeURUd5Zh11Yt24dgAtBnurSpQvXXHONh/2wTinJQrAFaDh85rDsvMNKgdOB30raCJwLzCrUDuP169czePBgSktL045izSguLuamm25KO4ZZm0uyECwCTpM0SlJ34EZg1uGFEbErIgZGxMiIGAm8DFwVEYsTzJSXDhw4wObNm3000AH069cv7QhmbS6xQhARB4HbgbnA68BPI2KlpG9Juiqpz+2INm3aRH19PaecckraUcysACX6zOKImA3MbjTvG820nZZklny2bt06ioqK/EQyM0uF7yzOA+vWrePkk0+mW7duaUcxswLkQpCyXbt2sWPHDp8W6iQOHjyYdgSzo5boqSFr3Zo1a4DMTUvW8T388MPceOONdO/evU22N3PmzGaXzZgxo00+I43PamsdOXs+8BFBytasWUP//v09rEQnsXHjRh588EFqa2vTjmKWMxeCFNXV1bFhwwbGjBnjYSU6iRtuuIFt27Zx3333ecRS6zBcCFK0bt066uvrfVqoExk7diw333wz+/bt49577z3yxDmzfOZCkKI1a9ZQXFzsy0Y7meHDh3PbbbfRu3dvfvzjH7NgwYK0I5m1yIUgJRHBm2++yWmnnUZRUVHacayNDRgwgNtuu41x48b5cZeW91wIUrJlyxZqa2t9WqgTKy4u5rrrruMTn/hE2lHMWuRCkJLXX3+dLl26cOqpp6YdxRIkiYsuuijtGGYtciFIQUSwcuVKRo8eTUlJSdpxzKzAuRCkYMuWLezatYuJEyemHcXMzIUgDStXrqSoqIhx48alHcXywKxZs9i3b1/aMayAuRC0s4hg1apVnHLKKfTo0SPtOJYHXn31Ve666y5WrFhBRKQdxwqQC0E7q66upqamxqeF7Ijp06fTr18/nnrqKR566CHfkWztzoWgnb322msUFRUxduzYtKNYnhgyZAi33norl19+OdXV1fzgBz/g+eefp66uLu1oViA8+mg7qq+vZ+XKlYwZM4bi4uK041ge6dKlCxUVFYwbN45nn32WefPmUVlZmXYsKxA+ImhHa9asoba2lsmTJ6cdxfJUnz59+MxnPsOtt97q5yNbu0m0EEi6TNJqSWslfb2J5X8taZWk5ZKek3RyknnStnTpUnr37u2byKxVw4cP59Zbb007hhWIxAqBpCLgLuByYAJwk6QJjZotBcoj4kzgCeBfksqTtpqaGtauXcvkyZPp0sUHYtY6D01u7SXJb6QKYG1ErI+IOuBR4OqGDSLi+YjYm518GRiWYJ5ULVu2jIjgrLPOSjuKdRI/+clP2Lhxoy85teOWZGfxUKCqwXQ1MLWF9rcBc5paIGk6MB3okEM2RwRLly5l5MiR9O/fP+041kls3bqVBx54gBNPPJHzzz+fCRMm+GjTjkle/KuR9CdAOfCdppZHxD0RUR4R5WVlZe0brg2sX7+enTt3+mjA2tQdd9zBlVdeSV1dHU8++SR33nkn8+fPZ+/eva2vbNZAkkcEW4DhDaaHZed9iKRLgP8B/EFE7E8wT2rmz59P7969mTChcReJ2bHr1q0b55xzDmeffTZr1qxhwYIFPPvss/zmN79h4sSJnHPOOQwfPtx9DdaqJAvBIuA0SaPIFIAbgc82bCDpLOCHwGURsT3BLKl5++23Wb9+PRdffDFdu/q2DWt7khg7dixjx45l27ZtLFmyhGXLlrF8+XIGDRrEWWedxRlnnEGvXr3Sjmp5KrFvpog4KOl2YC5QBNwXESslfQtYHBGzyJwK6g08nv2rZXNEXJVUpjQsWLCA7t27U15ennYUOwozZ87skJ81ePBgFi1adGR6+/btzJ07l7lz57bZZ1jnk+ifqBExG5jdaN43Gry/JMnPT1tNTQ2vvfYaU6ZM8QBzZpa38qKzuLNasGABEcG5556bdhQzs2a5ECTk/fffZ9GiRUyaNMlDBZhZXnMhSMhvf/tbAKZNm5ZqDjOz1vgylgRs27aNZcuWcf7559O3b9+045gdt+9+97sMHTqUoUOHMmzYME488cRjvgouIti3bx+7d+8+8qqpqWH37t3s2bOn2fXuu+8++vTpQ9++fSkrK2PIkCGUlZVRVFR0rLtlWS4ECXjuuefo0aMHF154YdpRzNrESSedRFVVFStXrgQyl6z27duXE044gT59+tCzZ09KSkro2rUrXbp04cCBA9TV1R157d+//0Nf/PX19R/5jJKSEkpLS5vNUFRUxFtvvcUbb7xxZP3u3bt3yNEG8o0LQRtbtWoVb775JpdeeiklJSVpxzFrE9dddx0Au3fvZsuWLbz11lvs3LmT9957j40bN7Jnz54mv9y7d+9O9+7dKS4uprS0lBEjRtC7d29KS0spLS2lT58+lJaW0rt37yNHGM1dTvv5z38egEOHDvHuu+/y9ttvs3nzZjZu3JjMThcQF4I2tHfvXmbPns2JJ57oK4WsUyotLWXcuHGMGzfuQ/Mjgvr6+iOvbt260bVr10Tuau7SpQtlZWWUlZVxxhlnAC3fixERvru6Fe4sbkNz585l3759XH311R78ywqKJLp27UpxcTE9e/akW7duefPl+8Mf/pAVK1Zw6NChtKPkLX9btZFVq1axfPlyLrroIgYPHpx2HDPLOnToEE899RTf+973qKysbPIUVqFzIWgD27Zt4+mnn2bo0KFcdNFFaccxswa+/OUvc/3119OjRw9+/vOfc+edd/LKK69w4MCBtKPlDfcRHKe9e/fy6KOPUlxczA033OBL2czyjCTGjx/PuHHjWLduHS+88ALPPPMML7zwAueeey5TpkyhuLg47ZipciE4DnV1dTz22GPs3r2bW265pcVL38wsXZI49dRTOfXUU9m0aRMvvPACzz33HC+99BIVFRVMnTqVnj17ph0zFS4Ex+iDDz7g4Ycfprq6mmuuuYZhwzrtUzbNOp2TTz6Zk08+ma1bt/LCCy8wb948FixYQHl5Oeedd17B/VHnQnAM9uzZwyOPPMLbb7/Ntdde6wfOmHVQJ510EjfccAPbt2/nxRdf5OWXX2bhwoWcfvrpVFRUcNJJJ6UdsV24EBylDRs28OSTT7J//36uv/56xo4dm3YkMztOgwYN4pprrmHatGksWLCAZcuWsWzZMoYPH05FRQXjx4/v1P1/LgQ52r9//5HDx/79+3PzzTczaNCgtGOZWRvq378/V1xxBRdffDGvvvoqCxcu5Mknn6RXr15MnDiRM844g6FDh+bNPRJtxYWgFQcPHmTZsmU8//zz1NbWMnnyZC6//HK6d++edjQzS0iPHj0499xzmTp1Km+++SavvvoqS5YsYeHChZxwwglMnDiRCRMmMGTIkE5RFFwImvHOO++wfPlylixZwt69exk2bBg33XQTQ4cOTTuambUTSYwZM4YxY8bwwQcf8MYbb/Daa6/x0ksv8eKLL1JSUsKoUaOOvPr3798hC0OihUDSZcB3yTyz+N6I+OdGy4uBB4FzgHeBGyJiY5KZmhIR1NTUUF1dTVVVFWvXruXdd98FYOzYsVRUVDBq1KgO+Qs2s7bRo0cPJk+ezOTJk9mzZw/r1q1jw4YNbNiwgVWrVgGZEVQHDx585DVo0CBOOOEESkpK8vr7I7FCIKkIuAu4FKgGFkmaFRGrGjS7DdgZEadKuhH4NnBDEnkOf9HX1tayZ88eamtr2bt3L7t27eKdd96hrq4OgK5dux7pIBo3bhx9+vRJIo6ZdWC9e/dm0qRJTJo0iYjgvffeY8OGDbz11lts27aNysrKD9253K1btyPPUigtLaWkpIQePXrQo0cPSkpKKC4upqioiKKiIrp27UpRURGSiIgPvUpLSxN54mGSRwQVwNqIWA8g6VHgaqBhIbga+Gb2/RPA9yQpIqKtw1RVVfHEE08cme7Zsye9evWitLSUSZMmMXDgQIYOHcqQIUM69dUBZta2JDFgwAAGDBhwZF5EsHPnTrZt28auXbvYtWsXNTU1R/7w3Ldv3zENcXHBBRdwySWXtGV8AJTAd25mw9K1wGUR8YXs9OeAqRFxe4M2r2XbVGen12XbvNNoW9OB6dnJscDqZj52IPBOM8s6i0LYR/B+diaFsI+Q//t5ckSUNbWgQ3QWR8Q9wD2ttZO0OCLK2yFSagphH8H72ZkUwj5Cx97PJEcf3QIMbzA9LDuvyTaSugJ9yXQam5lZO0myECwCTpM0SlJ34EZgVqM2s4DPZ99fC/wmif4BMzNrXmKnhiLioKTbgblkLh+9LyJWSvoWsDgiZgH/AfxY0lrgPTLF4ni0evqoEyiEfQTvZ2dSCPsIHXg/E+ssNjOzjsFPKDMzK3AuBGZmBa7DFgJJ/SU9K+nN7H9PaKZdvaRXs6/GndV5S9JlklZLWivp600sL5b0WHb5K5JGtn/K45fDft4iaUeD3+EX0sh5PCTdJ2l79r6ZppZL0p3Zn8FySWe3d8bjlcM+TpO0q8Hv8RvtnbEtSBou6XlJqyStlPSXTbTpeL/Pxrcwd5QX8C/A17Pvvw58u5l2e9LOegz7VgSsA0YD3YFlwIRGbb4C3J19fyPwWNq5E9rPW4DvpZ31OPfzY8DZwGvNLP8UMAcQcC7wStqZE9jHacAv0s7ZBvt5InB29n0psKaJf7Md7vfZYY8IyAxP8UD2/QPAp1PM0taODM8REXXA4eE5Gmq4/08AFyufR7VqWi772eFFxDwyV8U152rgwch4Gegn6cT2Sdc2ctjHTiEi3oqIyuz73cDrQOMhiRkE7fcAAAaaSURBVDvc77MjF4LBEfFW9v3bwOBm2vWQtFjSy5I6SrEYClQ1mK7mo//YjrSJiIPALmAAHUsu+wnwmewh9hOShjexvKPL9efQ0Z0naZmkOZImph3meGVPx54FvNJoUYf7feb1EBOSfg0MaWLR/2g4EREhqbnrYE+OiC2SRgO/kbQiIta1dVZLzM+BRyJiv6Q/I3MU9ImUM9nRqyTz/+IeSZ8CngZOSznTMZPUG3gSuCMiatLOc7zyuhBERLPD7EnaJunEiHgre9i1vZltbMn+d72k35Kp4PleCI5meI7qDjw8R6v7GREN9+leMn1DnU0uv+8OreGXZUTMlvR9SQOj0QCTHYGkbmSKwE8i4qkmmnS432dHPjXUcHiKzwM/a9xA0gnZh98gaSBwAR8eBjtfFcrwHK3uZ6Nzq1eROSfb2cwCbs5ebXIusKvBac9OQdKQw31YkirIfPd0tD9cyO7DfwCvR8S/NtOsw/0+8/qIoBX/DPxU0m3AJuB6AEnlwJciM/z1eOCHkg6R+Yf3z/HhB+PkpUhneI52l+N+/oWkq4CDZPbzltQCHyNJj5C5amagpGpgBtANICLuBmaTudJkLbAX+NN0kh67HPbxWuDLkg4C+4AbO+AfLpD5Y/JzwApJr2bn/T0wAjru79NDTJiZFbiOfGrIzMzagAuBmVmBcyEwMytwLgRmZgXOhcDMrMC5EFiiGo3++mpTI4zmC0l3SOrZYHq2pH6trLMxe49KU/NXZF+rJP0vST2yy06S9EQL2+wn6SutfO787H+nSfpFa/vWaN1PS5rQYPpbkpq9edM6P18+aomStCcieqedIxeSNgLlR3O3a3PrNJyfHY7gHuBARHz+o1v5yDZHkhmp8/QmlnXNji11eHoa8LcRceVRZL4/u/1mi5EVFh8RWLuT1FeZZxCMzU4/IumL2fd7JP1bdqz35ySVZedPzg4cuFzSfyn7/AlJv5X0bUkLJa2RdFF2fpGk70halF3nz7Lzp2XXeULSG5J+kr0D9C+Ak4DnJT2fbXvkr31JT0taks01/Wj2NyL2AF8CPq3MczRGKjtuv6SJ2eyvZnOeRuZmyVOy876TzfyCMs/TWHX459TgI/pI+mX2Z3q3pC6N20i6VtL9ks4nc4f2d7LbPyU7/9psu4slLc0eydyn39+Zv1HSTEmV2WXjjuZnYPnNhcCSVtLo1NANEbELuB24X9KNwAkR8aNs+15k7iqeCPyOzB2qAA8CX4uIM4EVDeYDdI2ICuCOBvNvI3Nr/xRgCvBFSaOyy87Ktp1A5lkIF0TEncBW4OMR8fEm9uPWiDgHKCdzt/NRjfSaHWtnAx8daO1LwHcjYnJ229Vknq+xLiImR8TfZdudDfxlRIxpYvMVwFez+3MKcE0LOeaTGQLh77LbPzLuVvbU1f3ADRFxBpmRB77cYPV3IuJs4AfA3+a049YhuBBY0vZlv3AOvx4DiIhnyXyh3wU0fOrYIeCx7PuHgAsl9QX6RcTvsvMfIPMglMMOD/y1BBiZff9JMuO9vEpmmOAB/P5LeGFEVEfEIeDVBuu05C8kLQNeJjOg2LGMnNnU8yIWAH8v6WtkRufc18y6CyNiQwvL1kdEPfAIcOExZAMYC2yIiDXZ6Vx+ztYJuBBYKrKnL8aTGYulyceMZuXSibU/+996fj9+loCvNihAoyLiV43aN16nuazTgEuA8yJiErAU6JFDrobbKCXz5bmm4fyIeJjMqZp9wGxJzQ2xXdvC5hv/jKKJ+UeVtxlN/ZytE3AhsLT8FZmRRD8L/KcyQ/tC5t/ktdn3nwVezJ5K2nn4/D+ZQb9+R8vmkhnkrBuApDGSerWyzm4yjx9srC+wMyL2Zs+Nn9vKdj4k21n8feDpiNjZaNloYH321NTPgDNbyNGcCmVGcO0C3AC8mJ2/TdL47Pw/btC+ue2vBkZKOjU7ncvP2ToBV3VLWol+P0ojwDPAf5I5HVQREbslzQP+gcz5/VoyX2z/QOYZEzdk1/s8cLcyl3eup/URHe8l8xd4pSQBO2j9cab3AM9I2tqon+AZ4EuSXifzZflyK9s57PnsZ3cB/gv4xybaXA98TtIBMk/a+98R8Z6kl7IdynOAX7byOYuA7wGnAs9nPwsyfQ2/ILPvi4HDV289Cvwo20F+uOgSER9I+lPgcWWecbEIuDvHfbUOzJePWl5RB7rc1Kyz8KkhM7MC5yMCM7MC5yMCM7MC50JgZlbgXAjMzAqcC4GZWYFzITAzK3D/HzPwu/vAte8rAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MtOaxPgg1Hym",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "from numpy import exp\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "from scipy.special import factorial\n",
        "import pandas as pd\n",
        "from mpl_toolkits.mplot3d import Axes3D\n",
        "import statsmodels.api as sm\n",
        "from statsmodels.api import Poisson\n",
        "from scipy import stats\n",
        "from scipy.stats import norm\n",
        "from statsmodels.iolib.summary2 import summary_col"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Hio8uXrtLjl7",
        "colab_type": "code",
        "outputId": "b73fec62-bc96-4a74-9bc4-69fc96c75e72",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 299
        }
      },
      "source": [
        "# generate a sample\n",
        "import math\n",
        "import scipy.optimize as optimize\n",
        "\n",
        "random.seed(101)\n",
        "size=4000;\n",
        "sample = normal(loc=40, scale=6, size=4000)\n",
        "# calculate parameters\n",
        "sample_mean = sample.mean()\n",
        "sample_std = sample.std()\n",
        "print('Mean=%.3f, Standard Deviation=%.3f' % (sample_mean, sample_std))\n",
        "# define the distribution\n",
        "dist = norm(sample_mean, sample_std)\n",
        "# sample probabilities for a range of outcomes\n",
        "values = [value for value in range(10, 70)]\n",
        "probabilities = [dist.pdf(value) for value in values]\n",
        "# plot the histogram and pdf\n",
        "pyplot.hist(sample, bins=15, density=True)\n",
        "pyplot.plot(values, probabilities)\n",
        "pyplot.show()\n",
        "\n",
        "l1= size*math.log(2*math.pi*(sample_std**2))/2 + np.sum(((sample-sample_mean)**2)/(2 * (sample_std**2)))\n",
        "ll = size*math.log(2*math.pi*(sample_std**2))/2 + np.sum(((sample-sample_mean)**2)/(2 * (sample_std**2)))\n",
        "print('log likelihood:', ll)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Mean=40.097, Standard Deviation=6.029\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3xU9Z3/8ddnZpJwR7koyMWAXBQQRSJIFUEoFmorutUVbVerdm23dWu323Zxd2ut66+r3d3abrW/1q1abxWtt6WCUhVRQW7hpnLTEIKACOEuAklmzmf/OCcYQ0ImZCbfOWc+z8djHjNzzpnM+6vDJyefc853RFUxxhgTXTHXAYwxxmSXFXpjjIk4K/TGGBNxVuiNMSbirNAbY0zEJVwHqK9bt25aXFzsOoYxxoTKsmXLdqpq94bW5VyhLy4uprS01HUMY4wJFRHZ1Ng6a90YY0zEWaE3xpiIs0JvjDERZ4XeGGMizgq9McZEnBV6Y4yJOCv0xhgTcVbojTEm4qzQG2NMxOXclbHG5KLi6bMAEDy+Gn+VC2Lv8ufUGF7yziVF/JivrbjrktaIaEyjrNAbk6Z+so27C+5nVGw9+7Qdk+NL2apdeTh5MTNSF7GfDq4jGtMga90Y05RUkm/FZ/JS4XQGy2Z+UPNNRlTdzzeq/5EPvJP554InWFj09/wk8TBFVLtOa8xRbI/emGOpfA+e/QbTC1bxUupcflzzdSo5EYBXvJG84o1kSLKCGxIvcX1iDjUk+Fnyq45DG/NZVuiNaUyqBv50HXz8EX9XfQsveqMb3GyNFvODmm9RpQV8Iz6bOakSlungVg5rTOOsdWNMY5b+HnasgUv/u9EiX9fPktfwIV35j4Lf0YaqVghoTHqs0BvTkAM74LWfwWkT4PQvpfWST2jLD2u+Sf/YR/ww8VSWAxqTPiv0xjTklduh5hBM+TmIpP2yhd5QHklO4vr4S5wr67KXz5hmsEJvTH2bl8LKx2HMt6HbwGa//K7k1WzRbvy84He05XAWAhrTPFbojanLS8HsH0DHnnDhD4/rRxykDT9KfpN+se38KPFkhgMa03xW6I2pa/kjsG0lXHwnFHU87h+zyBvCH5IXc31iDmxeksGAxjSfFXpjah3cDa/eAaeeD8O+0uIfd3dyGnu0Ayy8NwPhjDl+VuiNqTXvLji8t9kHYBtziDb8KTUO1s2Cjz/KQEBjjo8VemMADu+HFY/CWddAj2EZ+7F/TE0AL+m3hIxxxAq9MQCrn4OagzDy6xn9sRXa0z8Xv/QhSCUz+rONSZcVemPA35vvfjr0Lsn8zy65ET7+EN57KfM/25g0WKE3Zsc62LIURvxNRnrzRxk0GTr1gtIHMv+zjUlDWoVeRCaLyHoRKROR6Q2sLxKRJ4P1i0WkuM664SKyUERWi8g7ItImc/GNyYAVj0KsAM6alp2fH0/4LaENc2HXhuy8hzHH0GShF5E4cB8wBRgCXC0iQ+ptdiOwR1UHAPcAdwevTQCPAd9S1aHAeKAmY+mNaalkNax6AgZPgfbdsvc+51wLsQSUPpi99zCmEens0Y8CylS1XFWrgRnA1HrbTAUeDh4/DUwUEQEuBt5W1VUAqrpLVVOZiW5MBrz3Ihzc5bdtsqljD39ytJWP+3PoGNOK0in0vYDNdZ5vCZY1uI2qJoF9QFdgEKAiMkdElovIjxp6AxG5SURKRaS0srKyuWMw5vgtfxQ6ngIDJmb/vc69EQ7t8c/wMaYVZftgbAK4APhqcH+5iBz1L0pV71fVElUt6d69e5YjGRPYtxU2vApnXwOxY3/Bd0YUj4Vug2CpHZQ1rSudQr8V6FPnee9gWYPbBH35zsAu/L3/N1R1p6oeBGYD57Q0tDEZsfKPoB6M+FrrvJ+If6rl1lL4cGXrvKcxpFfolwIDRaSfiBQC04CZ9baZCVwXPL4CmKuqCswBzhSRdsEvgHHAmsxEN6YFPM8/26Z4LHTp13rve9Y0SLT1e/XGtJImC33Qc78Zv2ivBZ5S1dUicoeIXBps9gDQVUTKgO8D04PX7gF+gf/LYiWwXFVnZX4YxjTTpvmwd1P2D8LW1/YE/3jAulmg2rrvbfJWWl8Orqqz8dsudZfdVufxYeDKRl77GP4plsY4VTz9032MewruY2KsHec+UUjVE62873H6l2DdC/Dhcug1snXf2+QluzLW5J02VDE5tpQ/p8ZQRWHrBxj0BZC4v1dvTCuwQm/yzvmxd2kr1bzojXIToF0XKD4f1r7g5v1N3kmrdWNMlEyMLedjbcti74xWeb+6LaNa18b7cUfBG0y49X8o11MafF3FXZdkO5rJE7ZHb/KK4PH5+Ape94ZT43A/5+WUP0vmxbFSZxlM/rBCb/LKcCnnJNnLKym3B0G30ZVVXn++ELdCb7LPCr3JKxPjy0mpMM87y3UU/pIqYUSsjJPZ7TqKiTgr9CavTIotp1QHs5eOrqMwx/PbN5PiyxwnMVFnhd7kjV5UckbsA1523LapVaa92OD1tD69yTor9CZvTIwvB+AVL1emWxJe9koYE1tDJw64DmMizAq9yRufjy1ng9fT/8LuHDEnVUKBpLgoZpOcmeyxQm/yw+H9nBdbw8tebrRtaq3U09iuJ9jZNyarrNCb/LDhVQolxSupXGnb+JQYL6dGMj62iiKqXccxEWWF3uSH9S+yWzuwXAe5TnKUOd65tJMqLoi94zqKiSgr9Cb6Ukl4bw6veSPwcvAjv8gbwn5tx8UxO83SZEfufeqNybTNi+Hw3pw5rbK+GhK84Z3JuPgqwOaoN5lnhd5E3/rZEC/kTe9M10kaNd87kx6yhwFS/1s6jWk5K/Qm+ta/CMVj+YS2rpM0an7wS2is9elNFlihN9G2awPs3gCDJrtOckxbtDvlXg8r9CYrrNCbaCuf598PmOg0Rjrme2cyOraWApKuo5iIsUJvoq18HnTuA136u07SpPneMNpLFefI+66jmIixQm+iy0vBxjeg/zgQcZ2mSQu9oSQ1xtj4266jmIixQm+ia9sqOLwX+l/kOklaPqYdK3WAXThlMi6tQi8ik0VkvYiUicj0BtYXiciTwfrFIlIcLC8WkUMisjK4/Taz8Y05htr+fL8LncZojvneMIbLRk7gY9dRTIQ0WehFJA7cB0wBhgBXi8iQepvdCOxR1QHAPcDdddZtUNWzg9u3MpTbmKaVz4OTh0GHk1wnSdsbqeHERPlcbLXrKCZC0tmjHwWUqWq5qlYDM4Cp9baZCjwcPH4amCgSgqaoia6aQ/DBIug/3nWSZlmlp7Ff21r7xmRUOoW+F7C5zvMtwbIGt1HVJLAP6Bqs6yciK0TkdREZ29AbiMhNIlIqIqWVlZXNGoAxDfpgEaSqQlfoU8RZ5A3hwvg7oDYdgsmMbB+M3Qb0VdURwPeBP4pIp/obqer9qlqiqiXdu3fPciSTF8rnQawA+o5xnaTZ3vCG01t2wu5y11FMRKRT6LcCfeo87x0sa3AbEUkAnYFdqlqlqrsAVHUZsAHIvXliTfSUz4M+o6Cog+skzTbfG+Y/2DDXbRATGekU+qXAQBHpJyKFwDRgZr1tZgLXBY+vAOaqqopI9+BgLiLSHxgI2G6Kya6Du/1TK/uPd53kuFRoD7ZoN9jwmusoJiKaLPRBz/1mYA6wFnhKVVeLyB0icmmw2QNAVxEpw2/R1J6CeSHwtoisxD9I+y1V3Z3pQRjzGRvfADS0hR6EN1JnQsWb/lz6xrRQIp2NVHU2MLvestvqPD4MXNnA654BnmlhRmOap3weFHaEU3LrawObY753JtdUvQZbl0Hf0a7jmJCzK2NN9JTPg35jIZ7WfkxOessbCoj16U1GWKE30bKnAvZsDHHbxreXjnDKCCi3Pr1pOSv0JlrKX/fv+493mSIz+o/zWzdVB1wnMSFnhd5ES/k86NgTukXgLN7iseAlYfMi10lMyFmhN9HhebDxdX9vPgozcPQZDbEEbHzTdRITclboTXRsfxcO7oJ+41wnyYyiDtBrJFTMd53EhJwVehMdFcGeb4imJW5S8Vj4cAVU2bTF5vhZoTfRUTHf/8rAzvXn3Aux4gtAU7BpoeskJsSs0Jto8FJQscAvjFHSZ7Q/OVuF9enN8bNCb6Lho3egah8UR6htA1DYDnqfa4XetIgVehMNtQcsi893myMbii8Ivv92n+skJqSs0JtoqJgPXU6DTqe4TpJ5/caCetanN8fNCr0JPy8Fm97yC2IU9R4F8SJr35jjFt5Zn0xeKp4+66hlw6ScF4r28d2FHZi54Oj1oVfQxu/Tb3zDdRITUrZHb0LvvNhaABZ5ZzhOkkX9xvoHnA/tcZ3EhJAVehN6Y2Jr2OD1ZAcnuo6SPcVjAfVbVMY0kxV6E2pxUpwbW8cib4jrKNnVuwQSbWw6BHNcrNCbUBsim+gkh6LdtgFIFPlfdm4TnJnjYIXehNp5sTVAxPvztYovhO3v+F9+bkwzWKE3oTYmtoYy7xQqo9yfr1U7vcOmBW5zmNCxQm9Cy+/Pr8+PvXnwpywuaGftG9NsVuhNaA2VCjrKoegfiK2VKPQnObMLp0wzpVXoRWSyiKwXkTIRmd7A+iIReTJYv1hEiuut7ysiB0TkB5mJbUzd/nyeFHrw2zc71sAnO10nMSHSZKEXkThwHzAFGAJcLSL1/2XdCOxR1QHAPcDd9db/Anix5XGN+dSY2Bre93qxk86uo7Se2i9VsT69aYZ09uhHAWWqWq6q1cAMYGq9baYCDwePnwYmivhf2ikilwEbgdWZiWxMHvbna50ywu/T2/n0phnSKfS9gM11nm8JljW4jaomgX1AVxHpAPwT8NNjvYGI3CQipSJSWllZmW52k8eGyUY6yGEW5lPbBiBeAH3Ps0JvmiXbB2NvB+5R1QPH2khV71fVElUt6d69e5YjmSgYE/Tnl+TbHj1Yn940WzqzV24F+tR53jtY1tA2W0QkAXQGdgGjgStE5OfACYAnIodV9d4WJzd57bzYWt7Lt/58reJgOuaK+TD0MrdZTCikU+iXAgNFpB9+QZ8GXFNvm5nAdcBC4ApgrqoqcGSCcBG5HThgRd60VIIkJbH1PJOK6PzzgYamZAZ//G8XFfGnJx7jJ8mCBrepuOuSbEYzIdNkoVfVpIjcDMwB4sCDqrpaRO4ASlV1JvAA8KiIlAG78X8ZGJMVZwb9+bw6rbKOJAlKvUFHpmc2pilpffGIqs4GZtdbdludx4eBK5v4GbcfRz5jjlJb4BbnY38+sNAbyj8VzKAr+9iVj+0r0yx2ZawJnTGx1azz+rCbTq6jOFN7Wumo2DrHSUwYWKE3oeL359/Lv/Pn63lH+/GJFh25OtiYY7FCb0JluJTTTqry7/z5evw+/WDr05u0WKE3oXLekfPnT3ecxL1F3hAGx7bQlX2uo5gcZ4XehMqY2BrWen3Zk8f9+Vq17avRtldvmmCF3oRHstr683V82qe3Qm+OzQq9CY8Pl9NWqvO+P18rSYKl3ul2QNY0yQq9CY+Nb+Kp5PX58/Ut8s5gUGyr9enNMVmhN+FR8SZrtS/76OA6Sc6ovTrY+vTmWKzQm3BIVsHmxXk77UFj3tViDmgb69ObY7JCb8Jh6zJI5uH880349Hx669ObxlmhN+Gw8U1AWOINdp0k59T26btZn940wgq9CYeKN6HHmey3/vxRav/KGROzb+s0DbNCb3JfzWHYvOTTL8Y2n/GO9me/tuNzVuhNI6zQm9y3tRRSVf5X6JmjeMRY5J3B+bF3XUcxOcoKvcl9G98EiUHfMa6T5Kz53jD6xirpI9tdRzE5yAq9yX0bX4eeZ0HbE1wnyVkLvGEAnG/tG9MAK/Qmt1V9DFuWQv+LXCfJaRv0FLbrCda+MQ2yQm9y26a3wEtC//Guk+Q4YYE3jDGxNQie6zAmx1ihN7mtfB4k2kCf0a6T5Ly3vKF0k/0Mli2uo5gcY4Xe5Lbyef5B2II2rpPkvAWp2j69tW/MZ1mhN7nr4+2wY421bdK0ja6Uez3sfHpzlLQKvYhMFpH1IlImItMbWF8kIk8G6xeLSHGwfJSIrAxuq0Tk8szGN5G28XX/vv94lylC5S1vqD+TZarGdRSTQ5os9CISB+4DpgBDgKtFpP7MUjcCe1R1AHAPcHew/F2gRFXPBiYDvxORRKbCm4grnwdtT4Qew10nCY0F3jA6yGHYutx1FJND0tmjHwWUqWq5qlYDM4Cp9baZCjwcPH4amCgioqoHVTUZLG8DaCZCmzyg6hf6fuMgZh3GdC30huCpfPrXkDGkV+h7AZvrPN8SLGtwm6Cw7wO6AojIaBFZDbwDfKtO4T9CRG4SkVIRKa2srGz+KEz07CqD/VutbdNMe+nIGj0Vyq3Qm09lfVdJVRer6lDgXOBWETnq9AlVvV9VS1S1pHv37tmOZMKgfJ5/33+8wxDhtMAbCluWQPVB11FMjkin0G8F+tR53jtY1uA2QQ++M7Cr7gaquhY4AAw73rAmj5TPgxNOhS79XCcJnQXeMEhVwwcLXUcxOSKdQr8UGCgi/USkEJgGzKy3zUzguuDxFcBcVdXgNQkAETkVOB2oyEhyE12ppD+RWf/xrpOE0lJvMMQKrE9vjmjyDBhVTYrIzcAcIA48qKqrReQOoFRVZwIPAI+KSBmwG/+XAcAFwHQRqQE84NuqujMbAzERsm0lVO2zQn+cDtEG+oyyPr05Iq1THVV1NjC73rLb6jw+DFzZwOseBR5tYUaTb8pf8+/7jXObI8z6jYN5/w4Hd0O7Lq7TGMfsvDWTe8pf98+db9/VdZLw6j8OUNj4huskJgfYxUum1RVPn9XourYcZmXRQh5KTeauY2xnmtCrBIo6Q9krMPQy12mMY7ZHb3LKubH1FEnyyBdpmOMUT8BpF/mFXu06xXxnhd7klPNj71KlCf/MEdMyAyfBx9tgu81mme+s0JucMi72Nsu8QRymyHWU8Bvwef/+/Zfd5jDOWaE3OeMUdnJ6bDNzvRGuo0RDxx7+QW0r9HnPCr3JGRPiKwB4zTvbcZIIGTgJNi+GQ3tdJzEOWaE3OeOi2Eo2eSexQU9xHSU6Bl4Mmvr02gSTl6zQm5xQRDWfi60O2jbiOk509CqBNp3h/VdcJzEOWaE3OWFMbDVtpdraNpkWT8BpE6DsZTvNMo9ZoTc5YUJsJQe1iMXeGa6jRM/Ai+HAdvjobddJjCNW6E0OUCbEV7DAG0YVha7DRI+dZpn3rNAb5wbKVnrLTuZa2yY7OpwEPc/yr5I1eckKvXFuQiw4rTJlhT5rBl4cnGa5x3US44AVeuPchPgK1nin8hE2W2XWDJgE6sEGO80yH1mhN0514gAj5T1r22Rb7xJoe6L16fOUFXrj1LjY2yTEY27Kpj3Iqlg8OM3yFfA812lMK7NCb5y6KL6SXdqRlTrAdZToGzAJPtlhp1nmISv0xpkYHuNjK3ndOwvPPorZN+DzgMB7c1wnMa3M/nUZZ86WMrrIATvbprV06A59RsPama6TmFZmhd44c1F8JUmN8bo33HWU/DFkqv9FJLs2uE5iWpEVeuPMhNgKlukg9tPBdZT8ccaX/fs1/+s2h2lVaRV6EZksIutFpExEpjewvkhEngzWLxaR4mD5JBFZJiLvBPcTMhvfhNWp8hFDY5t4OTXSdZT8ckIf6DXSCn2eabLQi0gcuA+YAgwBrhaRIfU2uxHYo6oDgHuAu4PlO4Evq+qZwHXAo5kKbsLtktgiAGalznOcJA8NmQrbVsKeTa6TmFaSSGObUUCZqpYDiMgMYCqwps42U4Hbg8dPA/eKiKjqijrbrAbaikiRqla1OLkJtS/FF1PqDWKbXQ2bFcXTZzW6ro905M0iuPM/7+b3qUuOWl9x19HLTLil07rpBWyu83xLsKzBbVQ1CeyDo/4FfwVYbkXe9JcPGRLbxKzUaNdR8tJmPZl3vGK+GF/sOoppJa1yMFZEhuK3c77ZyPqbRKRUREorKytbI5JxqLZtM9sKvTMvpkZxTqyMHuxyHcW0gnQK/VagT53nvYNlDW4jIgmgM/ifIBHpDTwHXKuqDZ7Tpar3q2qJqpZ07969eSMwoXNJfDFLvMFsp4vrKHnrRc//JTs5vtRxEtMa0in0S4GBItJPRAqBaUD9Ky5m4h9sBbgCmKuqKiInALOA6aq6IFOhTYjtWMfpsc28YAdhndqoPVnr9WFKfInrKKYVNFnog577zcAcYC3wlKquFpE7ROTSYLMHgK4iUgZ8H6g9BfNmYABwm4isDG4nZXwUJjxWP4enwoupUa6T5L2XUqM4V9bTHZujPurSOesGVZ0NzK637LY6jw8DVzbwujuBO1uY0USFKqx+jiV6OpWc6DpN3pvtjeYf5Bm+EC/lsdQk13FMFtmVsab17FgLO9db2yZHvK+9KPNOYUrM2jdRZ4XetJ7Vz4LEeMnaNjlCeNEbxXmxNXRhv+swJous0JvWEbRtKL6AnXR2ncYEXkyNIi7KxfFS11FMFlmhN61j+7uwqwyGXu46ialjjZ7KRu9kpsbech3FZJEVetM63n0WJA5nXNr0tqYVCU+nxjEmvoa+st11GJMlVuhN9tW2bfpdCO27uU5j6nk6dSEpFf46Ps91FJMlVuhN9n2wEPZshGFfcZ3ENGA7XZjnnc2V8deJk3Idx2SBFXqTfaUPQVEnGPZXrpOYRjyVGs/JspdxsVWuo5gssEJvsuvgbv9LLoZfBYXtXacxjXjVG0GlduIqa99EkhV6k10r/wipKii53nUScwxJEjyTupCJseXwsR2UjRor9CZ7VGHZH6D3KDh5qOs0pgl/So0jIR6sesJ1FJNhVuhN9lTMh13v2958SGzQXiz1BsGKR/1f0iYyrNCb7Fn2ELTpbBdJhciTqYv8C9s+WOQ6iskgK/QmOz7ZCWtmwllXQ0Fb12lMmmalRkNhR1j+iOsoJoOs0JvsWPk4eDUw0to2YXKINv5psGueh8M20VlUWKE3med5/kHYvp+Dk053ncY01znXQc1BePcZ10lMhlihN5lX8QbsLreDsGHV6xw4aQiUPmgHZSPCCr3JvNKHoG0Xm8AsrETgvL+Dj96GDa+6TmMywAq9yayPt8O6F+Dsa6Cgjes05ngNnwadesGbv3CdxGSAFXqTWW/9t//nfskNrpOYlkgUwuf+HjYtgE0LXacxLWSF3mTOx9th6QP+vDZdT3OdxrTUOddCu67w5n+5TmJaKOE6gAmv4umzPvP8XxKPcX28iolLSti0eFYjrzKhUdgezvs2zP032LYKep7lOpE5TlboTUZ0Yx9fi7/C894FbNIeruOYFqj7C7wTfZlf1JY3fvNDbq65pcnXVtx1STajmeOUVutGRCaLyHoRKROR6Q2sLxKRJ4P1i0WkOFjeVUReE5EDInJvZqObXHJT4gUKqeHe5FTXUUwG7ac9j6Ym8cXYEvrLh67jmOPUZKEXkThwHzAFGAJcLSJD6m12I7BHVQcA9wB3B8sPAz8GfpCxxCbndGMffxN/mee9C6jQnq7jmAx7MDmFKgr4u/hM11HMcUpnj34UUKaq5apaDcwA6u+2TQUeDh4/DUwUEVHVT1R1Pn7BNxFVuzf/6+RlrqOYLNhFZ2akLuKy+AJ6Uek6jjkO6RT6XsDmOs+3BMsa3EZVk8A+oGu6IUTkJhEpFZHSykr7IIWJ7c3nh/uTXwLgbxN2kD2McuL0SlW9X1VLVLWke/furuOYZvhb25vPC9voyrOpsVwdf42+Yt9AFTbpFPqtQJ86z3sHyxrcRkQSQGdgVyYCmtzVlX1cG3+Z573zbW8+D/xX8kqqSXBH4g+AzYETJukU+qXAQBHpJyKFwDSg/lGZmcB1weMrgLmqNhtS1N1a8AQJUtybtC8WyQc7OJFfJr/C+PgqvhBb6jqOaYYmC33Qc78ZmAOsBZ5S1dUicoeI1M5a9QDQVUTKgO8DR07BFJEK4BfA10VkSwNn7Jgw2jCXK+Jv8NvUl9loe/N54w+pL7DW68tPCh6hnZ1jERppXTClqrOB2fWW3Vbn8WHgykZeW9yCfCYXVX8Cf/4eG7ye3Gu9+bySIs6/1lzPM0U/5buJ57grebXrSCYNOXEw1oTMaz+DvZu4teYbVFHoOo1pZct0ME8mx3NjfDYDZYvrOCYNVuhN82xdDot+AyOvZ4me4TqNceSu5DQO0JY7Cx7EDszmPiv0Jn2pGpj5XWh/Ekz6qes0xqE9dOLu5DRGx9ZxeWy+6zimCVboTfre+jVsfwcu+S9o09l1GuPYk6nxrPAG8M8Fj9OVfa7jmGOwQm/Ss2sDzLvL/3rAM77kOo3JAUqMW2u+QUcO8ZvCX5Eg6TqSaYQVetO0w/vhya/5Xw34xf9wncbkkHXal+k1f8vo2Dr+JfG46zimEVbozbGlkvD0DVC5Hq58GDraXPPms573LuD3ySlcn5gDK6zY5yL74pE8V/9bouq7PfEHvp54mVtrbuSJ/zkI2KRW5mj/nryGM+QDzn/hH+Ck06HXSNeRTB22R28adW18Dl9P/IX7k5fwRGqi6zgmh6WIc3PN30OHk2HG1+DADteRTB1W6E2DxsdW8JPEI/wlNdKufjRp2UMnmPY4HNoDT10LyWrXkUzACr05ymD5gF8X3MtaPZVbar6DZx8Tk66ew2HqvfDBQnjiKqg64DqRwQq9qWeUrGVG4Z18QhturP4Bh2jjOpIJmzOvgKn3Qfk8eORS+MRmLHfNCr054rLYfB4r/Bm7tBNXVt/Gdrq4jmTCasTX4KrHYftqeGgy7N3c9GtM1lihN4DyvcTT/LLwN5R6g/mr6tvZrCe7DmXC7vQvwteehY+3w4Nf8E/RNU5Yoc9zhdRwT8Fv+F7iWf6UvJDraqaznw6uY5moKD4frp8FXtIv9mWvuk6Ul6zQ57Ntb/NU4U+5PL6An9f8NT9MfpMau7TCZFqPM+GGOf5keI/9FTz/Hf/MHNNqJNe+8a+kpERLS0tdx4i2w/tg7v+Dpf/DTq8DP665nhe90a5TmYgroprvJp7lm/EX2E0nflxzPXO8c4+sr7jrEofpwk9ElqlqSUPrbI8+n6jCqifh1yWw5H4ouYEJVf9pRd60iioK+Y/kNKZW30mlduZ3hfdwXxhFAScAAAflSURBVMEv6YGdlZNt9nd6PkglYf1sWHgfbF7kX57+1afglBHsf9OmNDCta7UWM7X637gpPotbEs8yqWgZf/bGwLY+/nn4JuOsdRMBjc1XcxJ7mBZ/jWsSr9JD9rBVu3Jv8jJmpC5C7Y85kwN6SyU3xF/kqvhrtJcq6HchjLkZBkyCmH1Gm+NYrRvbo4+Yk9jD2Ng7fD6+jM/HllMgKV5PDedfUzcw1xthV7manLJFu3NH8lp+mfwKb395Gyz6Lfzxr6FTbxg8GQZNhuKx/hTZ5rhZoQ+7g7sZG3ubsbF3GBt7mzNi/oUpldqZh1KTeTw1kU1qUwub3Laf9nD+LXDet2HN/8Lq52DlH2Hp76GgHfS/CAZMgJ5nw8lDoaCt68ihklbrRkQmA78C4sDvVfWueuuLgEeAkcAu4CpVrQjW3QrcCKSA76rqnGO9l7VuGpCsho+3wf4PYU8F7FgN29f4Vx0e+AiAKk2w1BvMm95w3vCGs077WHvGhFoR1ZwXW8uE2HImxlfQW3YCkNQYG/QUVmsx67w+/PM1k+GEPtC5L7TvBiKOk7txrNZNk4VeROLAe8AkYAuwFLhaVdfU2ebbwHBV/ZaITAMuV9WrRGQI8AQwCjgFeAUYpKqpxt4vVIVeNbh5wS0FXurTey8FXo3/pdqpGkhV+7fkYag5CDWHgttB/u2ZRXSWT+jMJ3SSg3TmE7rKfnrKbrqxj5h8+v+pSgt4X3uxXvuwzuvDWj2VUm8Qhyly+B/DmGxSestOhkoFQ2IVDJMKhsYq6CH1zsdPtIXOvaBdN2jXBdp2gXYn+veFHaCwnf8XQmF7/76gLcQLIF5Y51YAsQRIzL+PxUHi/vMjN8m5Xygt7dGPAspUtTz4YTOAqcCaOttMBW4PHj8N3CsiEiyfoapVwEYRKQt+3sLjGcgxfbgCHkr3PNxGfrmpfrruyC9AbWB5UNwz6McF4Kmwn3bs0/bsoz17tCNrvFP5SLuwjS58pF3Yqt2o0B6kiGf0/Y3JbcIW7c4W7f6Zc+87cYC3bzkD9m3259PZtxn2bYGDu2DvB/DhSv/irOShLMWKAbVFv879kfW1j5tYVmvIpXD5bzMeM51C3wuoOyPRFqD+iddHtlHVpIjsA7oGyxfVe22v+m8gIjcBNwVPD4hISybF6AbsbMHrHdtX90nIx3KUKI0nSmOBEI9H7j5qUWjHAr8Lbp+R7nhObWxFThyMVdX7gfsz8bNEpLSxP1/CJkpjgWiNJ0pjgWiNJ0pjgcyMJ52jdVuBPnWe9w6WNbiNiCSAzvgHZdN5rTHGmCxKp9AvBQaKSD8RKQSmATPrbTMTuC54fAUwV/2jvDOBaSJSJCL9gIHAksxEN8YYk44mWzdBz/1mYA7+6ZUPqupqEbkDKFXVmcADwKPBwdbd+L8MCLZ7Cv/AbRL4zrHOuMmQjLSAckSUxgLRGk+UxgLRGk+UxgIZGE/OTYFgjDEms+yKGmOMiTgr9MYYE3GhLfQi8qCI7BCRd+ss6yIiL4vI+8H9iS4zNoeI9BGR10RkjYisFpFbguWhG5OItBGRJSKyKhjLT4Pl/URksYiUiciTwcH9UBCRuIisEJEXgudhHkuFiLwjIitFpDRYFrrPWS0ROUFEnhaRdSKyVkTGhHU8IjI4+P9Se9svIt9r6XhCW+iBPwCT6y2bDryqqgOBV4PnYZEE/lFVhwDnAd8JppAI45iqgAmqehZwNjBZRM4D7gbuUdUBwB78OZDC4hZgbZ3nYR4LwEWqenad87PD+Dmr9SvgJVU9HTgL//9TKMejquuD/y9n488ddhB4jpaOR1VDewOKgXfrPF8P9Awe9wTWu87YgrH9L/78QqEeE9AOWI5/NfVOIBEsHwPMcZ0vzTH0Dv5xTQBewL92PZRjCfJWAN3qLQvl5wz/mp2NBCeWhH089cZwMbAgE+MJ8x59Q05W1W3B44+Ak12GOV4iUgyMABYT0jEFrY6VwA7gZWADsFdVk8EmDU6HkaN+CfwIqJ3gqCvhHQv4Ezf9RUSWBdOPQEg/Z0A/oBJ4KGit/V5E2hPe8dQ1DX9SSGjheKJW6I9Q/1df6M4dFZEOwDPA91R1f911YRqTqqbU//OzN/5Edqc7jnRcRORLwA5VXeY6SwZdoKrnAFPwW4QX1l0Zps8Z/rVA5wD/X1VHAJ9Qr60RsvEAEBzzuRT4U/11xzOeqBX67SLSEyC43+E4T7OISAF+kX9cVZ8NFod6TKq6F3gNv71xQjBFBoRnOozzgUtFpAKYgd+++RXhHAsAqro1uN+B3/8dRXg/Z1uALaq6OHj+NH7hD+t4ak0Blqvq9uB5i8YTtUJfdyqG6/D73KEQTOv8ALBWVX9RZ1XoxiQi3UXkhOBxW/xjDWvxC/4VwWahGIuq3qqqvVW1GP9P6bmq+lVCOBYAEWkvIh1rH+P3gd8lhJ8zAFX9CNgsIoODRRPxr8QP5XjquJpP2zbQ0vG4PuDQggMVTwDbgBr83+o34vdOXwXex/+Sky6uczZjPBfg/zn2NrAyuH0xjGMChgMrgrG8C9wWLO+PP9dRGf6fpEWuszZzXOOBF8I8liD3quC2GviXYHnoPmd1xnQ2UBp83p4HTgz5eNrjTwrZuc6yFo3HpkAwxpiIi1rrxhhjTD1W6I0xJuKs0BtjTMRZoTfGmIizQm+MMRFnhd4YYyLOCr0xxkTc/wEATHOimyhlrwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        },
        {
          "output_type": "stream",
          "text": [
            "log likelihood: 12862.206950784697\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q214ixOGeR35",
        "colab_type": "text"
      },
      "source": [
        "### Read the csv files [```q4_1.csv```](https://raw.githubusercontent.com/Ipsit1234/SHALA/master/q4_1.csv) and [```q4_2.csv```](https://raw.githubusercontent.com/Ipsit1234/SHALA/master/q4_2.csv) to get samples from two unknown distributions. Choose a significance level and implement the the most general two-sample t-test"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bAsMqWiGdb52",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "col_Names=[\"Sequence\", \"Coverage\"]\n",
        "sample1 = pd.read_csv('https://raw.githubusercontent.com/Ipsit1234/SHALA/master/q4_1.csv',names=col_Names)\n",
        "sample2 = pd.read_csv('https://raw.githubusercontent.com/Ipsit1234/SHALA/master/q4_2.csv',names=col_Names)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rb29JEqw-Tth",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Student's t-test for independent samples\n",
        "from math import sqrt\n",
        "from numpy.random import seed\n",
        "from numpy.random import randn\n",
        "from numpy import mean\n",
        "from scipy.stats import sem\n",
        "from scipy.stats import t\n",
        "from numpy.random import seed\n",
        "from numpy.random import randn\n",
        "from scipy.stats import ttest_ind"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yIvOH3yWiTU9",
        "colab_type": "text"
      },
      "source": [
        "Find the parameters required to compute the T statistic"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_3TBrHZifG2u",
        "colab_type": "code",
        "outputId": "992f7985-e6bf-4f72-a7c7-0aa8425b9e81",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 68
        }
      },
      "source": [
        "# function for calculating the t-test for two independent samples\n",
        "def independent_ttest(data1, data2, alpha):\n",
        "\t# calculate means\n",
        "\tmean1, mean2 = mean(data1), mean(data2)\n",
        "\t# calculate standard errors\n",
        "\tse1, se2 = sem(data1), sem(data2)\n",
        "\t# standard error on the difference between the samples\n",
        "\tsed = sqrt(se1**2.0 + se2**2.0)\n",
        "\t# calculate the t statistic\n",
        "\tt_stat = (mean1 - mean2) / sed\n",
        "\t# degrees of freedom\n",
        "\tdf = len(data1) + len(data2) - 2\n",
        "\t# calculate the critical value\n",
        "\tcv = t.ppf(1.0 - alpha, df)\n",
        "\t# calculate the p-value\n",
        "\tp = (1.0 - t.cdf(abs(t_stat), df)) * 2.0\n",
        "\t# return everything\n",
        "\treturn t_stat, df, cv, p\n",
        "\n",
        "# seed the random number generator\n",
        "seed(101)\n",
        "# generate two independent samples from the csv files\n",
        "data1 = sample1['Coverage'].values[:]\n",
        "data2 = sample2['Coverage'].values[:]\n",
        "# calculate the t test\n",
        "alpha = 0.05\n",
        "t_stat, df, cv, p = independent_ttest(data1, data2, alpha)\n",
        "print('t=%.3f, df=%d, cv=%.3f, p=%.3f' % (t_stat, df, cv, p))\n",
        "# interpret via critical value\n",
        "if abs(t_stat) <= cv:\n",
        "\tprint('Accept null hypothesis that the means are equal.')\n",
        "else:\n",
        "\tprint('Reject the null hypothesis that the means are equal = YES')\n",
        "# interpret via p-value\n",
        "if p > alpha:\n",
        "\tprint('Accept null hypothesis that the means are equal.')\n",
        "else:\n",
        "\tprint('Reject the null hypothesis that the means are equal = YES')"
      ],
      "execution_count": 35,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "t=-16.821, df=148, cv=1.655, p=0.000\n",
            "Reject the null hypothesis that the means are equal = YES\n",
            "Reject the null hypothesis that the means are equal = YES\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SrOP9BwSiwVc",
        "colab_type": "text"
      },
      "source": [
        "Plot the distribution of both the samples to visualize"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "omfVDZqxfyAM",
        "colab_type": "code",
        "outputId": "2f9df978-20af-4415-b53f-a29f33beb932",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        }
      },
      "source": [
        "sns.distplot(sample1[1].values)\n",
        "sns.distplot(sample2[1].values)\n",
        "plt.legend(['first sample','second sample'])\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU1d348c/JZCP7vgeSALJDgLCJLO4LKi5Y1wpWRVweu9pq7c+2Pu3ztLV9bK1WStXiWlGkFutSUUFkEQhhJyxhTwhZyb5Ocn5/3EkMISETMpM7y/f9es1rZu49c+83N8k3J+ece47SWiOEEML9+ZgdgBBCCMeQhC6EEB5CEroQQngISehCCOEhJKELIYSH8DXrxDExMTotLc2s0wshhFvaunVrqdY6tqt9piX0tLQ0srOzzTq9EEK4JaXUse72SZOLEEJ4CEnoQgjhISShCyGEhzCtDb0rzc3N5Ofn09DQYHYoohcCAwNJSUnBz8/P7FCE8GouldDz8/MJDQ0lLS0NpZTZ4Qg7aK0pKysjPz+f9PR0s8MRwqu5VJNLQ0MD0dHRkszdiFKK6Oho+a9KCBfgUgkdkGTuhuR7JoRrcLmELoQQ4vy4VBt6Z29tOu7Q490xZWCPZZ577jlefPFFJkyYwK233srevXt5/PHH7Tr+0aNH2bBhA3fccUdfQ+2zthu3YmJizA5FCNFPXDqhm+Evf/kLn332GSkpKQBcf/31Z5WxWq34+p596Y4ePcpbb73lEgldCJeU/ffelc+6xzlxeChpculg0aJFHD58mKuvvppnn32WpUuX8sgjjwCwYMECFi1axJQpU/jxj3/Ml19+SWZmJpmZmYwfP57q6moef/xxvvrqKzIzM3n22WfPOHZhYSEzZ84kMzOT0aNH89VXXwHw4IMPkpWVxahRo/j5z3/eXj4tLY0nnniCzMxMsrKyyMnJ4corr2Tw4MEsXrwYgDVr1jBz5kzmzJnDsGHDWLRoEa2trWd9XW+88QaTJ08mMzOTBx54gJaWFmddQiGEiaSG3sHixYv55JNPWL16NTExMSxduvSM/fn5+WzYsAGLxcJ1113HCy+8wPTp06mpqSEwMJDf/OY3/P73v+ff//73Wcd+6623uPLKK3nyySdpaWmhrq4OgF//+tdERUXR0tLCpZdeys6dOxk7diwAAwcOZPv27Xz/+99nwYIFrF+/noaGBkaPHs2iRYsA2Lx5M3v37mXQoEFcddVVrFixgnnz5rWfNzc3l2XLlrF+/Xr8/Px46KGHePPNN7n77ruddBWFEGaxO6ErpSxANlCgtb62074A4DVgIlAG3Kq1PurAOF3CLbfcgsViAWD69On84Ac/4M477+Smm25qb6LpzqRJk/jOd75Dc3MzN9xwA5mZmQC88847LFmyBKvVSmFhIXv37m1P6G3NPWPGjKGmpobQ0FBCQ0MJCAigoqICgMmTJ5ORkQHA7bffzrp1685I6J9//jlbt25l0qRJANTX1xMXF+fAqyKEcBW9aXL5LpDbzb57gdNa6yHAs8Bv+xqYKwoODm5//fjjj/PSSy9RX1/P9OnT2bdv3zk/O3PmTNauXUtycjILFizgtdde48iRI/z+97/n888/Z+fOncyZM+eM8dwBAQEA+Pj4tL9ue2+1WoGzhwx2fq+1Zv78+Wzfvp3t27ezf/9+fvGLX5zX1y+EcG12JXSlVAowB3ipmyJzgVdtr5cDlyoPH5x86NAhxowZw09+8hMmTZrEvn37CA0Npbq6usvyx44dIz4+nvvvv5/77ruPnJwcqqqqCA4OJjw8nKKiIj7++ONex7F582aOHDlCa2sry5Yt46KLLjpj/6WXXsry5cspLi4GoLy8nGPHup19Uwjhxuxtcvkj8GMgtJv9ycAJAK21VSlVCUQDpR0LKaUWAgvBaB/uiT3DDM3yxz/+kdWrV+Pj48OoUaO4+uqr8fHxwWKxMG7cOBYsWMD3v//99vJr1qzhmWeewc/Pj5CQEF577TXS09MZP348w4cPJzU1lenTp/c6jkmTJvHII4+Ql5fHxRdfzI033njG/pEjR/KrX/2KK664gtbWVvz8/HjhhRcYNGhQn6+BEMK1KK31uQsodS1wjdb6IaXUbOBHXbSh7wau0lrn294fAqZorUvPOqBNVlaW7rzARW5uLiNGjDivL8QbrVmzpttO2P4m3zthFxm22GdKqa1a66yu9tnT5DIduF4pdRR4G7hEKfVGpzIFQKrtZL5AOEbnqBBCiH7SY0LXWj+htU7RWqcBtwFfaK3v6lRsJTDf9nqercy5q/6iz2bPnu0StXMhhGs473HoSqmngWyt9UrgZeB1pVQeUI6R+IUQQvSjXiV0rfUaYI3t9VMdtjcAtzgyMCGEEL0jt/4LIYSHkIQuhBAewrXncuntEKeeuMEQKFcZihgSEkJNTY2pMQghekdq6EII4SEkoXdQW1vLnDlzGDduHKNHj2bZsmUAbN26lVmzZjFx4kSuvPJKCgsLAcjLy+Oyyy5j3LhxTJgwgUOHDqG15rHHHmP06NGMGTOm/Rhr1qxh9uzZzJs3j+HDh3PnnXfSNrLzk08+Yfjw4UyYMIEVK1Z0GduePXvap8AdO3YsBw8eBOCGG25g4sSJjBo1iiVLlrSXDwkJ4bHHHmPUqFFcdtllbN68mdmzZ5ORkcHKlSsBWLp0KXPnzmX27NkMHTqUX/7yl12e+5lnnmHSpEmMHTv2jCl+hRCuxbWbXPrZJ598QlJSEh9++CEAlZWVNDc381//9V/861//IjY2lmXLlvHkk0/yyiuvcOedd/L4449z44030tDQQGtrKytWrGD79u3s2LGD0tJSJk2axMyZMwHYtm0be/bsISkpienTp7N+/XqysrK4//77+eKLLxgyZAi33nprl7EtXryY7373u9x55500NTW1z2n+yiuvEBUVRX19PZMmTeLmm28mOjqa2tpaLrnkEp555hluvPFGfvazn7Fq1Sr27t3L/Pnz22dy3Lx5M7t37yYoKIhJkyYxZ84csrK+uQnt008/5eDBg2zevBmtNddffz1r165t/5qEEK5DEnoHY8aM4Yc//CE/+clPuPbaa5kxYwa7d+9m9+7dXH755QC0tLSQmJhIdXU1BQUF7XOnBAYGArBu3Tpuv/12LBYL8fHxzJo1iy1bthAWFsbkyZPbp9nNzMzk6NGjhISEkJ6eztChQwG46667zqhpt5k2bRq//vWvyc/P56abbmov/9xzz/HPf/4TgBMnTnDw4EGio6Px9/fnqquuav+6AgIC8PPzY8yYMRw9erT9uJdffjnR0dEA3HTTTaxbt+6shP7pp58yfvx4AGpqajh48KAkdCFckCT0Di644AJycnL46KOP+NnPfsall17KjTfeyKhRo9i4ceMZZbubVfFcOk6Ba7FY2qfAtccdd9zBlClT+PDDD7nmmmv461//io+PD5999hkbN24kKCiI2bNnt0+/6+fn1z6VbsfpdztOvQv2Tb/7xBNP8MADD/TuixVC9DtpQ+/g5MmTBAUFcdddd/HYY4+Rk5PDsGHDKCkpaU/ozc3N7Nmzh9DQUFJSUnj//fcBaGxspK6ujhkzZrBs2TJaWlooKSlh7dq1TJ48udtzDh8+nKNHj3Lo0CEA/vGPf3RZ7vDhw2RkZPDoo48yd+5cdu7cSWVlJZGRkQQFBbFv3z6+/vrrXn/Nq1atory8nPr6et5///2zZny88soreeWVV9pHvBQUFLRPxSuEcC2uXUPv52GGu3bt4rHHHsPHxwc/Pz9efPFF/P39Wb58OY8++iiVlZVYrVa+973vMWrUKF5//XUeeOABnnrqKfz8/Hj33Xe58cYb2bhxI+PGjUMpxe9+9zsSEhK6XQAjMDCQJUuWMGfOHIKCgpgxY0aXtf933nmH119/HT8/PxISEvjpT39KcHAwixcvZsSIEQwbNoypU6f2+muePHkyN998M/n5+dx1111nNLcAXHHFFeTm5jJt2jTA6Gx94403ZNUjIVxQj9PnOotMn2u+pUuXkp2dzfPPP9/nY8n3TthFps/ts75OnyuEEMINuHaTi3CqBQsWsGDBArPDEEI4iMvV0GUadfcj3zMhXINLJfTAwEDKysokQbgRrTVlZWXt4/CFEOZxqSaXlJQU8vPzKSkpMTsU0QuBgYHtN0wJIczTY0JXSgUCa4EAW/nlWuufdyqzAHgGY21RgOe11i/1Nhg/Pz/S09N7+zEhhBDYV0NvBC7RWtcopfyAdUqpj7XWne9iWaa1fsTxIQohhLBHjwndtthz28TYfraHNHILIYSLsatTVCllUUptB4qBVVrrTV0Uu1kptVMptVwpldrNcRYqpbKVUtnSTi6EEI5lV0LXWrdorTOBFGCyUmp0pyIfAGla67HAKuDVbo6zRGudpbXOio2N7UvcQgghOunVsEWtdQWwGriq0/YyrXWj7e1LwETHhCeEEMJePSZ0pVSsUirC9noAcDmwr1OZxA5vrwdyHRmkEEKIntkzyiUReFUpZcH4A/CO1vrfSqmngWyt9UrgUaXU9YAVKAcWOCtgIYQQXbNnlMtOYHwX25/q8PoJ4AnHhiaEEKI3XOpOUSGEB9Iati6FrX+HpjoYegXEDDU7Ko/kUnO5CCE80Kqn4N/fAxTUlcGmF6Foj9lReSRJ6EII59nzPmx4DrLuhYVrYOaPISQBdi6Dplqzo/M4ktCFEM7RUAkf/QgSx8HVvwWlwC8QMm+Hxmo49IXZEXocSehCCOfY8GeoLYHrngOL3zfbw1MhYSwc3wDWBvPi80CS0IUQjldXDhv/AqNvhqTMs/cPvhia6yE/++x94rxJQhdCOF72K9BcCzMf63p/xCAITYCCrf0bl4eTYYtCiG9k/7135bPuOXubtRE2L4HBl0LciK4/pxQkT4R9HxojX4Kiex+rOIvU0IUQjrV7BdQUwbSHz10uaYLxfGqn82PyEpLQhRCOozVsfAFiR8DgS85dNigaQhOhaG//xOYFJKELIRzn2Hoo2gVTHzSaVXoSNwLKD0OzjHZxBEnoQgjH2bQYBkTC2G/ZVz5uFOgWKD3g3Li8hCR0IYRjVJwwOjknzAe/AfZ9JnIQ+PhBeZ5zY/MSktCFEI6x5SXjedK99n/Gxxci04xmF9FnktCFEH3XXA85r8KwayBiYO8+G5UBlQXSju4AktCFEH23aznUn4Ypi3r/2agMQMPpIw4Py9vYswRdoFJqs1Jqh1Jqj1Lql12UCVBKLVNK5SmlNiml0pwRrBDCBWkNm/4KcSMh7aLefz4yDZSPNLs4gD019EbgEq31OCATuEopNbVTmXuB01rrIcCzwG8dG6YQwmXlfda7oYqd+QZAWAqUH3J8bF6mx4SuDTW2t362h+5UbC7wqu31cuBSpc7nOyuEcCtaw5r/NdrNx952/seJzoCK49BidVxsXsiuNnSllEUptR0oBlZprTd1KpIMnADQWluBSuCsyRmUUguVUtlKqeySkpK+RS6EMF/eZ8YEWzN+BL7+53+cqAxotULlccfF5oXsSuha6xatdSaQAkxWSo0+n5NprZdorbO01lmxsbHncwghhKvQrbD6f4za+bjb+3asyHTj+fSxvsflxXo1ykVrXQGsBq7qtKsASAVQSvkC4UCZIwIUQrioE5vhZA7MerxvtXOAgFDjDlOpofeJPaNcYpVSEbbXA4DLgX2diq0E5ttezwO+0Fp3bmcXQniKxmrIXQkDL+x77bxNeCpU5jvmWF7Knhp6IrBaKbUT2ILRhv5vpdTTSqnrbWVeBqKVUnnAD4DHnROuEMIl7H3fmPf8uj+Cj4NuZwlPNZasa65zzPG8UI8LXGitdwLju9j+VIfXDcAtjg1NCOGSTmw2OkKHXgmxwxx33IhU47kyH2IucNxxvYjcKSqEsF/lCdj1DkQPhaFXOPbY4baEXnHCscf1IpLQhRD2aaox1gr1DzFmVPSxOPb4/sHGoheVktDPlyR0IUTPWlsg5zWjMzTrOxAQ4pzzhKdIQu8DSehCiJ7t/8hYhGL0vN7Pptgb4anGotFNtc47hweThC6EOLfSg3Docxg4DQZ2nsbJwdr+WMjwxfMiCV0I0T1rI+x8G4JiYOQNzj9fWLLxXFXg/HN5IEnoQoju5X4AdeXGzUO+Ac4/n38wBEYYC16IXpOELoToWlUhHFtvzHEePbj/zhuWBNWS0M+HJHQhRNf2f2jUyi/oPHWTk4UlQ00xtDT373k9gCR0IcTZqgqhaDdkzDaaQfpTWLIxk2PNqf49rweQhC6EONuRNeDjB2kz+v/cYUnGc9XJ/j+3m5OELoQ4U1OdMVdL6qT+r50DBMeAxV9GupyHHifnEkJ4mZM5xupBAy/suWz23x1/fuUDoYlSQz8PUkMXQpzpxCaj2SM8xbwYwpKMhC7LKvSKJHQhxDdqS425VJInmRtHWLIxL7o0u/SKJHQhxDcKtxvPiePMjaOtY/TULnPjcDP2LEGXqpRarZTaq5Tao5T6bhdlZiulKpVS222Pp7o6lhDCxRXugPCBEBRlbhztCX23uXG4GXs6Ra3AD7XWOUqpUGCrUmqV1npvp3Jfaa2vdXyIQoh+UX3KaG4ZPsfsSMA30JgbvUhq6L3RYw1da12otc6xva4GcoFkZwcmhOhneZ8Zz3EjzY2jTViy1NB7qVdt6EqpNIz1RTd1sXuaUmqHUupjpdQoB8QmhOhPBz+FgHAITTI7EkNYEpQflrnRe8HuhK6UCgHeA76nta7qtDsHGKS1Hgf8GXi/m2MsVEplK6WyS0pKzjdmIYSjtbbA4TUQNxyUMjsaQ1gyoKGoc+uu6I5dCV0p5YeRzN/UWq/ovF9rXaW1rrG9/gjwU0rFdFFuidY6S2udFRsb28fQhRAOU7QHGiqNxZ9dRVvHaJE0u9jLnlEuCngZyNVa/183ZRJs5VBKTbYdt8yRgQohnOjYeuO5P6fJ7cmAKPAPlYTeC/aMcpkOfBvYpZSyDVLlp8BAAK31YmAe8KBSygrUA7dpLbd4CeE2jq03ln8bEGl2JN9QCuJHGf89CLv0mNC11uuAczaqaa2fB553VFBCiH6kNRzbAEOvNDuSs8WPgl3vGjG6Stu+C5M7RYXwdiX7oK4M0qabHcnZEkZDY5UxPl70SBK6EN6urf18kB2zK/a3+NHGs4xHt4skdCG83dH1xtjzyHSzIzlb3AjjWdrR7SIJXQhvd3wjDJrmmm3UAaEQmSYjXewkCV0Ib1Z1EqoLIcXk6XLPJX601NDtJAldCG9WkGM8J080N45ziR8N5YeMpfHEOUlCF8KbncwBH19IGGN2JN2LHwW6FUpyzY7E5UlCF8KbFWw1Zlf0G2B2JN2Lt831J80uPZKELoS3am2Fk9sgeYLZkZxbZDr4BUtCt4MkdCG8VflhY0KuJBdP6D4+ED9SErodJKEL4a1OukGHaJv4Ucb6ojJF1DlJQhfCWxXkgO8AiB1udiQ9ix8NDRXGMEvRLUnoQnirgq2QlAkWeyZdNVnbFADS7HJOktCF8EatLUYTRmKm2ZHYJ962zqncMXpOktCF8EZlh8BaD4ljzY7EPoHhED5Qaug9kIQuhDc6tdN4duUbijqLHyU19B702HimlEoFXgPiAQ0s0Vr/qVMZBfwJuAaoAxZorXMcH64QwiFO7QSLP8QMMzuSc8v++zevW61QegA2LQGLX/efybrH+XG5KHtq6Fbgh1rrkcBU4GGl1MhOZa4GhtoeC4EXHRqlEMKxTu0yRrf4+psdif3Cko0pAGqKzI7EZfWY0LXWhW21ba11NZALJHcqNhd4TRu+BiKUUokOj1YI0XdaQ+FOSHCT9vM2YbaUIkMXu9WrNnSlVBowHtjUaVcy0HGNqHzOTvpCCFdQfQrqSt2nQ7RNcCz4+EG1JPTu2J3QlVIhwHvA97TWVedzMqXUQqVUtlIqu6Sk5HwOIYToq1O7jGd36hAFUD4QmiA19HOwK6ErpfwwkvmbWusVXRQpAFI7vE+xbTuD1nqJ1jpLa50VGxt7PvEKIfqqbYRL28067iQsGaoKZAqAbvSY0G0jWF4GcrXW/9dNsZXA3cowFajUWhc6ME4hhKOc2mnMYBgYZnYkvReWCE210FhtdiQuyZ57fqcD3wZ2KaW227b9FBgIoLVeDHyEMWQxD2PYoveOGxLC1Z3a5X7NLW1Ck4zn6pPu+QfJyXpM6FrrdcA5V4/VWmvgYUcFJYRwksZqY9rccXeYHcn5CbMl9KqT7jGpWD+TO0WF8CanbHdautsIlzb+wRAYIR2j3ZCELoQ3cdcRLh2FJUlC74YkdCG8yamdEBQNoW58319oonG3aKvV7EhcjiR0IbxJW4eoOme3mGsLSwLdAjXFZkficiShC+EtWpqheK/73fLfWZjtJnRpdjmLJHQhvEXpAWhpcv+EHhwLPhZJ6F1wg7WnhBDnreP0s/lbjOeyvDO3uxsfC4QkypwuXZAauhDeojLfmNwqJM7sSPouLFFq6F2QhC6Et6g6aSRC5QG/9mFJ0FgFjTVmR+JSPOA7K4TokdbGpFZhHjKrddvXIc0uZ5CELoQ3qD8NzXUQnmJ2JI4R2mEKANFOOkWFV3hr03GnHfuOKQOddmyHqbLNZu0pNfSAEAgIk4TeidTQhfAGVQWAcu87RDsLk5EunUlCF8IbVOZDSCz4BpgdieOEJhvL6bW2mB2Jy5CELoQ38KQO0TZhScZ8LrWynGUbSehCeLqmOqNT1BMTOkg7egc9dooqpV4BrgWKtdZnLUKolJoN/As4Ytu0Qmv9tCODFKK3tNYcKa1l2/EKCivr+epgKQABvj4EB/iSEB5IYvgAQgK8YFxAW4eop4xwaRMSB8pia0efYHY0LsGen+alwPPAa+co85XW+lqHRCREH+QVV/P6xmOs2lvEycqG9u0D/CwoBU3WVqyt3ywwnBgeyIjEMEYnh5MQFmhGyM7naSNc2vj4Qki81NA7sGcJurVKqTTnhyLE+dtxooJn/rOfdXml+Ft8uHh4LA9dPISpGVGkRAaxIqegvWxdk5XCygbyT9ez71QVq/cV88W+YgZGBTElPYqxKRFYfNx4etnOqgqMIX4BoWZH4nhhScbcNAJw3Dj0aUqpHcBJ4Eda6z0OOq7wMr0dL15V38wne06x/UQFwQG+XDEynqy0qPamlM1HTrP5yOkzPhPk78vg2BAGx4Yw64JYahqtbDt+mi1Hy3l3az6f5RYx+4I4xg+KwNfHA7qZKvM9r3beJiwJCrKhqdZYns7LOSKh5wCDtNY1SqlrgPeBoV0VVEotBBYCDBzoBjdjCJe2q6CS97cV0NzSyqwLYpl1QSyBfpZeHyckwJcZQ2O5aEgM+4uq+WJfMf/cXsDq/cXMGhbLxIGR+FrcNLFbG42hfe685Ny5tI2rrzoJMV2mHa/S54Suta7q8PojpdRflFIxWuvSLsouAZYAZGVl6c77hbBHk7WVf20vYNuJClIiB3DLxFRiQ/s+vlopxfCEMIbFh3KwuIbPc4v41/aTrNlfwqXD4xg/MNL9mmKqCgANER5agWrr6K0qkISOAxK6UioBKNJaa6XUZIyhkGV9jkyILpyubeKNTcc4VdnAJcPjuHhYnMOTrFKKC+JDGRoXQl5JDav2FrFiWwFrD5Zw2Yh4RieH4+MuS7hV2Jqwwj00oQeE2qYAKOi5rBewZ9jiP4DZQIxSKh/4OeAHoLVeDMwDHlRKWYF64DattdS+hcMdKa3lzU3HaNWau6elMSzBuZ18SimGxoUyJDaE3MJqPt17ire3nCDpQAmXj0zggvgQlKsn9orjEBgBgWFmR+I8YcmS0G3sGeVyew/7n8cY1iiE0+wqqOSd7BNEBvlz97RBxIT03y3sSilGJoUxPDGUHScq+Cy3iFc3HmVQdBBXjEzotzjOS8Vxz21uaROeDIf2Q4sVLF5wX8E5ePdXL9zC+rxSPtpVyMCoIL49bRBB/ub82PooxfiBkYxJCSf76GlW7y/mb18dZvuJ0zwwczCXDI/Dx5Xa2OtPQ10pDJxidiTOFZYMuhVqTnnezVO9JAlduCytNav3F/NZbjEjE8O4dVIqfi4w2sTXx4epGdFMGBjJlqPlbD9RwX2vZTM4NpgHZg5m7vgkAnx7P9rG4U5uN549tf28TduQzKoCr0/o5v92CNEFrTWf7DnFZ7nFjE+N4PbJA10imXfk7+vD9CExrHlsNn+6LZMAXws/fm8n03/zBf/7cS6HS0xeHu1kjvEcnmpuHM4WHAMWf6iUdnSpoQuX06o1H+w4yaYj5UxJj+K6cUkuParEz+LD3Mxkrh+XxLq8Ul7feIyXvjrCX788zOS0KG6dlMo1YxIZ4N/PtfaCHAiKAf+g/j1vf1M+to7RfLMjMZ0kdOFSWlo1K3Ly2XaighlDY7hqVILrjySxUUoxY2gsM4bGUlzdwHtbC1i25Tg/fHcHP1+5h6tGJ3Dj+GSmZkQ7fzy71pCfDRGDnHseVxGWDAVbjLZ0LyYJXbiMJmsry7YcZ/fJKi4bYYwxd5dk3llcaCAPzh7MolkZbD5Szns5+Xy86xTLt+aTEBbI3PFJ3DQ+xXlDL08fNToJ02c45/iuJiwZjq2DunKzIzGVJHThEhqaW3j4zRx2n6zimtEJXDQ01uyQHEIpxZSMaKZkRPP03NGs2lvE+9sKeNnWJDMyMYybJhjNNXGdZnvsyzqo6QUfMg3YWRdN/ZGzk9yU9KjzPrZLCu/QMerFJKEL09U1WVn42lbW5ZVy/bgkpmZEmx2SUwT6WbhuXBLXjUuirKaRD3ac5J/bCvjVh7n8z0e5TB8Sw7enDuLSEfF9bpKJPZ1Dk28o9QGe8YexR6EJRlu6JHQhzFPd0My9S7PJPlbOM/PG0tziHTcZR4cEsGB6Ogump5NXXMP72wpYkZPPwte3MjAqiPtmpKM15z2yJ7Y8h5LI8eCmTVa9ZvE3Frzw8pEurjUOTHiVyrpm7np5M1uPn+ZPt43nliwPH17XjSFxIfzoymGs/fHFvHDHBGJC/HnqX3v4w6f7+fpwGS2tvfsjF9BYTnjtESOhexOZAkBq6MIcZTWNfPvlzeQV1/DinRO4YpSL30J/Dn1p6+7KzS/V/fkAABsTSURBVBNSyEyN5PN9RazccZKNh8u4dkwiQ+Pt60CNqTBuKCqJnEBYzSGHxubSwlKgYCvUlkGwZzbb9URq6KLfFVc1cNuSrzlUUsPf5me5dTJ3BqUUQ+JCWDgjg29PHURLq+bvG47yj83HqW209vj52NPbaFF+lIWP6odoXUhbx+ipnebGYSKpoYt+daS0lvmvbKa0ppGl90xm2mDvrEnZQynFiMQwhsaFsPZgKav3FXO4tJabxiczIrH72RPjyrdSHj6KVkv/TWDmEsKSjOdTu2DwxebGYhKpoYt+s/1EBTe/uIGaRitv3T9VkrmdfC0+XDI8jocuHkxYoC+vf32Mj3YVdtm27tdcSVTlHk5Fe/iEXF3xDzGmCj61y+xITCMJXfSLz3OLuH3J14QE+PLegxeSmRphdkhuJzF8AA/OGszUjCjW5ZXy8rrD1HRqgkko3YQPrRTGTjcpSpOFJ3t1k4skdOF0b28+zv2vZTMkLoT3HryQ9BhZzPd8+Vp8uH5cMt/KSqWgop4X1+RRUt3Yvj+xdANNviGUhXvoGqI9CUuG0gPQXG92JKaQhC6cRmvNHz87wOMrdnHR0FjeXjjVIWt/CshMjeC+izJosray+MtDHCurBa1JLN1AUfQUtI+Xdo+1zY1etMfsSEzRY0JXSr2ilCpWSu3uZr9SSj2nlMpTSu1USk1wfJjC3VhbWnlixS7++NlBbp6QwsvzswgO8NIk4ySpUUE8OHsIQf4WXll/hIoTewhuKKQw5kKzQzNP21TBJ7eZG4dJ7KmhLwWuOsf+q4GhtsdC4MW+hyXcWV2TlQde38rbW07w8MWD+f0tY11uLnNPERXsz8KZGUQG+VO64xMA707oAyKNKYPbFvfwMj3+lmmt1wLnmsJsLvCaNnwNRCilEh0VoHAvZTWN3PG3TXyxv5j/njuKx64c7rYzJrqL0EA/7puRwaX+uziiE9hZ68UdzkpB0vhvFvfwMo6oNiUDJzq8z7dtO4tSaqFSKlsplV1SUuKAUwtXcrysjnmLN5JbWMWLd07k29PSzA7Ja0T61DNV72KjJYvXNx7jRHmd2SGZJ3kClOyDplqzI+l3/fp/sNZ6idY6S2udFRvrJbPAeYndBZXc9OIGymubePO+KVw1Wu7+7E/JxV9i0c0MyLyZkEBflm44SlFVg9lhmSNpvNEx6oXj0R3RS1UAdJxVKcW2TXiwjvOXHCyq5s3Nxwnys3DPhWkcKKrhQJHJ62l6mdSiVdQFxFEfN4HvTLfy17WHeHXjUUYNsRDu12J2eP0ryTYpWUEODJxqbiz9zBE19JXA3bbRLlOBSq11oQOOK9zA9hMVvLrxKFFB/iyaNfisRRqE8/laa0kqWcfxhCtA+RAV7M+3pw6ittHK7w8l09TqZX0YoQkQmuSVI13sGbb4D2AjMEwpla+UulcptUgptchW5CPgMJAH/A14yGnRCpeyPq+Ud7JPMCg6mIUzMwgb4Gd2SF4puXgtltYmTiRc3r4tJTKIWyamcqA2iMVHE9DeMc38N7y0Y7THJhet9e097NfAww6LSLg8rTX/2XOKLw+UMCopjG9lpcqwxH4y+Pi7Z20bevwdmnxDCKvOO2O63MGATo7iHwVxJAY2cUtSWT9GarLk8bD/Q2iohMBws6PpN/JbKHrF2tLKT97byZcHSpiUFsXtkwdKMjeRr7WWiJoDxlS5XQwPnRtfzqzoSpYXxrK+3EkLUruitnb0wh3mxtHP5DdR2K3J2sojb23jnex8Lhkexw2ZSfjIGHNTxVTswke3UhLR9epESsHCgYWMCKnjxaOJHKz1kj6OJNsN6wXe1ewiCV3YpaG5hUVvbOWTPad46tqRXDYiXm4YMpvWxFZso3pAMvWBcd0W8/WBHw7OJ9LPyh8OJVPRbOnHIE0SFAURg7yuY1QSuuhRfVML97+WzRf7ivn1jaP5zkXpZockgJD6fIIaSyiJ7Hn6pFDfVn44uIBaq4X/O5RMU2s/BGg2L+wYlYQuzqm20cqCv29mfV4pz8wby51TBpkdkrCJPb2NFh8/ysJG2lU+LaiRRWmF7K8N4untIU6OzgUkT4CK48Yao15CErroVnVDM3e/spnsY6d59tZMbslK7flDol/4WmuJqdxFafiYXi01Nz2qmuviy3jjcBDLjnh4e3ryROO5INvcOPqRJHTRpdpGK/f8fQs7TlTw/O3jmZvZ5fQ8wiRx5Vvx0S3ntdTcHcklzIhr4v9tCyWnzIOnNE6aAD6+cPxrsyPpN5LQxVnqm1q499Ut5Bw/zXO3j+fqMTJ5pitRrVbiT2+hImQIDQG9nxPJR8Gfp1YSP6CVBzeGU9zgoWnAPwgSxsKJTWZH0m88+M+z6Djfir2sLa28/vUx8opruCUrhYq65vM6jnCe6Mrd+FtrOdSHhaAj/DVLplVw0+ooHtoYxluzKvD3xLw+cCpkvwLWJvD1Nzsap/PEb6E4T9bWVt7afJyDxTXcOD6ZzNRIs0MSnWlNYtlG6gLiqArO6NOhRkS08ExWFdll/vzSUztJU6eAtcFrFo6WhC4AaGnVvLPlBPtOVXP9uCSy0qLMDkl0IbJ6P0GNJZyMmd7lnaG9dW1qIw9cUMubh4N42xM7SdtmW/SSdnRJ6IJWrXkvJ5/dJ6u4ZkwiUzOizQ5JdEVrkkvW0uAfZdzq7yA/HlPLjPhGnvLETtLQBOMGoxOS0IUXaNWa97cVsP1EBVeMjOeiITFmhyS6kVTyFcENpyiIuQiU4351LQr+PKWKhAGtLNoYTnG9h6WFgVPh+Ca8YcpJD/vOid7QWvPBjpNkHzvNxcPimD2s+9vHhcm0ZvShv9LoF05ZxBiHHz7CX7Pkwgqqm3148Oswz7qTNHUK1BZD+WGzI3E6SeheSmvNR7sK2XSknBlDY7hshCRzVxZftomYip2cjJmOVs6Zi2V4uNFJurXMn19s96CZGdNmGM9H1pobRz+QhO6FtNas2lvE+kNlTMuI5qpRCTLRlosbfWgJdQFxlERkOvU816Y2smhYLW8dHsA/DntIJ2nMUGMFo8NrzI7E6exK6Eqpq5RS+5VSeUqpx7vYv0ApVaKU2m573Of4UIWjfL6vmDUHSpiUFsmcsYmSzF1cbHkO8eVbyM24B+3j/E7Lx0bXMtPWSbrVEzpJlYKM2UYNvdWT2pLOZs8SdBbgBeBqYCRwu1Kqq9mAlmmtM22Plxwcp3CQz/cV8cW+YiYOjGRuZrLMZ+4GRh1aQr1/FHmpN/fL+SwKnptSRWKQcSdpkSd0kmbMgvpyKNpldiROZc93ajKQp7U+rLVuAt4G5jo3LOEMq/cX83luMeNTI7hxgiRzdxBdsYuk0vXsS59Pi2VAv523rZO0plnx4MZwGlv67dTOkT7LePbwZhd7EnoycKLD+3zbts5uVkrtVEotV0p1OS2fUmqhUipbKZVdUlJyHuGK8/XlgRJW7S0iMzWCmyemSDJ3E6MOLaHRL5yDA2/t93MPD2/hmUnV5JT78WROqHuP+gtLhNjhktDt9AGQprUeC6wCXu2qkNZ6idY6S2udFRvb+0mFxPn58kAJ/9lzirEp4cyTZO42oip2k1K8hn1pd2P1DTYlhjkpjXx3RC3Ljw3gd7vNicFhMmbDsY1gbTQ7EqexJ6EXAB1r3Cm2be201mVa67ar9BIw0THhib7QWvPJ7sL2ZH7LxFRJ5m5kTN5faPCLYH/anabG8b2RtdyZUceL+4N56UD/Nfs43OBLwFoPR78yOxKnsSehbwGGKqXSlVL+wG3Ayo4FlFId51e9Hsh1XIjifLS0av65rYC1B0uZkh7Ft7JSsfhIMncX0ad3kFzyFbnpC0yrnbdRCp4eX8M1yQ38amcobxxy0+GM6bPALxhy/212JE7TY0LXWluBR4D/YCTqd7TWe5RSTyulrrcVe1QptUcptQN4FFjgrIBFzxqtLTzyVo7tDtBYrh+XJDVzNzM27y80+EdxcNDtZocCGCNfnp1cxaWJjfxsWxhvuuMYdb9AGHo57PsQWt29l7drdg0y1Vp/BHzUadtTHV4/ATzh2NDE+ahptLLo9a2syytlzphEpsvcLG4n5vQ2Eks3kDPsh1h9g8wOp12ABf4ytZIHN4bzZE4YLa2Ku4fUmx1W74y4Dva+D/lbvpmJ0YN4wABT0aagop55L25g4+Ey/nDLOEnmbmrswReo94/m4KD+H9nSkwALvDitkssSG3lqeyjP7gl2r9EvQ68Aiz/kfmB2JE4hCd1DbD12mrnPr6Ogop5XFkzi5okpZockzkNc2RYSyjaxN+Pefh133hsBFlg8rZJ5g+r5U24wP9sWQrO73IAZGGaMdsn9wCNnX5SE7ua01ry64Si3LdlIcIAv/3zoQmZdIENC3ZLWjD34Z+oCYskbeIvZ0ZyTrw88k1XNomHG4hjzv4rgdKOb9NOMuA4qjkHBVrMjcTgPmKjBe1U1NPPkP3fzwY6TXDo8jj98axwRQZ6/bqKnSin6grjT29g0+ue0WJzb6bjpSLlDjnNxSDk+aWEsOZbAVf8J5+UZNYyOtDrk2E4z8gb4+HHYuhRSssyOxqGkhu6mvj5cxtV//IoPd57ksSuH8be7sySZuzHV2kzm/mepCBnM4eQbzA6nV2ZFV/HzC47TpBU3rY7k5YMDXLs1IzAMxtwMu9+Dhiqzo3EoSehuprbRytMf7OX2v32Nn0Wx/MELefjiIfjIGHO3NvT4O4TVHWP7sB/0y4yKjnZBSAO/G3mUmfFN/PeOUO7+KoITtS6cXiYugOY62PWu2ZE4lAtfcdFR2xzml//fl7yy/gh3TB7Ih4/OYMLASLNDE30U2FjK2IMvUBg9lZOxM8wO57yF+bbwtwsr+e/x1eSU+3L5p9H8dX+Qa65+lDQBEsbA1r97VOeoJHQ3sPdkFXe9vIn7X8smNNCP9x6cxq9vHENwgPvV5MTZJuT+DktrI9mjfmbclunGlIJvD65n1RXlXBTXxP/uCuHy/0TxYX6Aa+VNpWDiPXBqFxzbYHY0DiMJ3YUdKqnhe29vY86fv2LPySp+ft1I/v3oRUwcFGV2aMJBEkvWk1b4MXsy7qM6eJDZ4ThMUlArL02v5NWLKgi0wMNfhzPn80g+zA+gxVUSe+YdEBwHa39ndiQOI1U8F7Qzv4K/fXWED3eeJMDXwsIZGTw0ewjhQX5mhyYcKKCxnCm7nqIyOJ29GfeaHY5TzEpo4qL4clYcC+TFfUE8/HU4aSFW7sqo55a0BsL9TczufgNg+qPw6c/g8JfGIhhuThK6i6hrsvLxrlO8veU4W46eJjTAl/tnZnD/jAxiQgLMDk84mm5l2q4nCWiuYE3WC7RaPHeEkkXBLWkN3DSogY/zA1iaN4Bf7Qzlmd0hXJzYyLUpjVyS2EiQGdlo0v2waQn856ew8EuwuHdKdO/o3Vxrq2bL0XJ+98l+dp2spMnaSlSwP3PGJDJxUCSBfhY+3VNkdpjCCUYefpmkknVsGflTKsKGmx2Ow/Q0vj0WeCwNjsYFsLo0nK+LwvikIJAAn1YmhNcwJbKaMaG1hPie2ZM6Jd1JzYx+gXDFf8O782Hjn+Gi7zvnPP1EEno/a7K2svlIOZ/lFrFqbxEFFfX4+/owNjmcCQMjGRQdJIs2e7jBJ94j88BzHE28moMDbzM7HFOkBTVyz8Bi5qcWs69mABtOh7HpdCgbT4eh0KQHNTAmtI7RYbUMD3HyBGAj58KI62H1/0D6TEh23+UclDap6zkrK0tnZ2ebcu7+dqK8jk1Hylm9v5i1+0uobrQS4OvD9CExXDcukco6K/6+0j/tDdIK/s3UnU9yKmYaayf+mVaf3vWLDD7uWeOmO2rRcLB2ALurgthVHczB2gG0aIUFzahIK+Ojmhkf3cz4KCsDg1u6HxCUdU/vT15bBktmg26Bez+FcNedC0kptVVr3eUtrlJDd7DWVs3h0lpyjp3m6yNlbDpcTkGFUcOICfHnmjGJXDYynulDognyNy7/W5uOmxmy6AdKtzD2wPOMOvwSRZET+Wr8//U6mXs6i4LhIfUMD6lnHmU0tChya4LYVzOAopYw3j0ayKuHjOmEw/xaGRFuZWSElRERVkaGWxkaZiXAcp4nD46G296EpXPg1evgrhUQle64L66fSELvg/qmFg6X1nCgqJpd+VXsLqhkz8lKapuMyfOjgv2Zkh7F/TPSmZIRzbD4ULmj0wtFVu5hYu5viTu9jYOp89g68qeSzO0QaNGMD69lfHgtU9JbsbbCgSpfcsp82VPhR26lL28fGUB9i/E7ZVGaIaEtXHBwG+kxwQyODSY9Jpi0mGDCAu243olj4a734M1b4G+XwHV/NJpi3KgJ1K6ErpS6CvgTYAFe0lr/ptP+AOA1jLVEy4BbtdZHHRtq/9NaU1HXTEFFvfE4XU/+6XoOldRwqKSGgor69pslAnx9GJkUxryJKYxODmdcagRD40KkPdxLqdZmEks3Mjh/BSlFX9DoH8nGMb/iSMpcs0NzW74+MDLCqJVDA2A00xyrsbC3wpfcSl9yK3zZfuI0/9558owbmWJCAsiICWZgdBBJ4YEkhA8gITyAhLABJIYHEhHkZ/yupk6G+z6H5ffAO3cb7ekX/QCGXGZ0oLq4HhO6UsoCvABcDuQDW5RSK7XWezsUuxc4rbUeopS6Dfgt4LTZ+bXWtLRqWrSmtRVabO9b27fpDtvA2tpKfXMLDc0t1De1UtdkbX9f19RCRV0zp+uaKK9toqKumfLapvb3jdYze9sH+FnIiA1mwsBIbpmYypC4EIbEhTA4Nhhfi7SDe5vguhOE1xzBt6WOgKbThNadIKL6ADEVO/FtqafBP4q9GfeyN+M7NPuFmh2ux7EoyAhtISO0hWtTbevUZ91EQ3MLJ8rrOFxay5HSWo6UGM9fHSyhpLqR1k5dhwG+PsSFBRAZ5E/4AD8iQv/AJb6fcknxG4QvuxOrZQBlMVnUhl9Ac/gg1IAofIIisA68iEB/f3wtCl8fH3x8wKK+ed1xm8VHOb2CZ08NfTKQp7U+DKCUehuYC3RM6HOBX9heLweeV0op7YQe1w93FvLwWzmOPizhA/yICvYnMsiPpIhARiWFERnsT3xYIMkRA4xH5AAi2/6SCwFkFKxkTN7i9vdWn0CqgtM4lHIjp6KnUBg7Q5pXTBDoZ2FofChD48/+I2ptaaWkppHCygaKKhuM5yrjUVHfTGV9MwUVzWyom8IT9WOZyk4u8dnO5MJ9pJ/aSID6ZnrgjIY3aD2PG+4XzRrM41c7frhqj6NclFLzgKu01vfZ3n8bmKK1fqRDmd22Mvm294dsZUo7HWshsND2dhiw31FfSCcxQGmPpVyLu8XsbvGCxNwf3C1ecL+YB2mtu1zFpl87RbXWS4Alzj6PUiq7u2E9rsrdYna3eEFi7g/uFi+4Z8zdsed/hQIgtcP7FNu2LssopXyBcIzOUSGEEP3EnoS+BRiqlEpXSvkDtwErO5VZCcy3vZ4HfOGM9nMhhBDd67HJRWttVUo9AvwHY9jiK1rrPUqpp4FsrfVK4GXgdaVUHlCOkfTN5PRmHSdwt5jdLV6QmPuDu8UL7hlzl0y79V8IIYRjycBpIYTwEJLQhRDCQ7htQldKRSmlVimlDtqez1otWSl1sVJqe4dHg1LqBtu+pUqpIx32ZbpCzLZyLR3iWtlhe7pSapNSKk8ptczWSW1qvEqpTKXURqXUHqXUTqXUrR329ds1VkpdpZTab7s2j3exP8B2zfJs1zCtw74nbNv3K6WudFaMvYz3B0qpvbZr+rlSalCHfV3+fLhAzAuUUiUdYruvw775tp+jg0qp+Z0/a1K8z3aI9YBSqqLDPlOucZ9prd3yAfwOeNz2+nHgtz2Uj8LosA2yvV8KzHPFmIGabra/A9xme70YeNDseIELgKG210lAIRDRn9cYo7P+EJAB+AM7gJGdyjwELLa9vg1YZns90lY+AEi3HcfiAvFe3OFn9cG2eM/18+ECMS8Anu/is1HAYdtzpO11pNnxdir/XxgDPky7xo54uG0NHWO6gVdtr18Fbuih/DzgY611nVOjOrfextxOGfMNXIIxtUKvP3+eeoxXa31Aa33Q9vokUIyxME1/ap+eQmvdBLRNT9FRx69lOXCp7ZrOBd7WWjdqrY8AebbjmRqv1np1h5/VrzHu/zCTPde4O1cCq7TW5Vrr08Aq4Conxdmmt/HeDvzDyTE5nTsn9HitdaHt9Skgvofyt3H2N+zXtn9pn1XGjJHOZm/MgUqpbKXU121NREA0UKG1bptIIh9IdmKs0MtrrJSajFEbOtRhc39c42TgRIf3XV2b9jK2a1iJcU3t+ayj9fac9wIfd3jf1c+Hs9kb88227/dypVTbDYkufY1tzVnpwBcdNptxjfvMpedDV0p9BiR0sevJjm+01lop1e34S6VUIjAGYyx9mycwkpQ/xjjUnwBPu0jMg7TWBUqpDOALpdQujATkcA6+xq8D87XWbVNUOuUaexOl1F1AFtBxSfqzfj601oe6PkK/+gD4h9a6USn1AMZ/RJeYHJM9bgOWa61bOmxz1Wt8Ti6d0LXWl3W3TylVpJRK1FoX2pJJ8TkO9S3gn1rr5g7Hbqt5Niql/g78yFVi1loX2J4PK6XWAOOB94AIpZSvrYbZ1RQMpsSrlAoDPgSe1Fp/3eHYTrnGXejN9BT56szpKez5rKPZdU6l1GUYf1hnaa0b27Z38/Ph7GTTY8xa647TfbyE0QfT9tnZnT67xuERnqk339fbgIc7bjDpGveZOze5dJxuYD7wr3OUPat9zJag2tqmbwB2OyHGznqMWSkV2dY0oZSKAaYDe7XRU7Maoy+g28+bEK8/8E/gNa318k77+usa92V6ipXAbbZRMOnAUGCzk+K0O16l1Hjgr8D1WuviDtu7/Plwcrz2xpzY4e31QK7t9X+AK2yxRwJXcOZ/y6bEa4t5OEZH7cYO28y6xn1ndq/s+T4w2j8/Bw4CnwFRtu1ZGKsqtZVLw/jL7NPp818AuzCSzBtAiCvEDFxoi2uH7fneDp/PwEg2ecC7QIALxHsX0Axs7/DI7O9rDFwDHMCoRT1p2/Y0RkIECLRdszzbNczo8NknbZ/bD1zdTz+/PcX7GVDU4Zqu7OnnwwVi/l9gjy221cDwDp/9ju3a5wH3uEK8tve/AH7T6XOmXeO+PuTWfyGE8BDu3OQihBCiA0noQgjhISShCyGEh5CELoQQHkISuhBCeAhJ6EII4SEkoQshhIf4/wKIwkRWBFi+AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}