GOF - 4. Goodness of fit.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as st\n",
    "import pylab\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from fitter import Fitter, get_common_distributions, get_distributions\n",
    "from scipy.stats import entropy as kl_div, kstest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Fitting a distribution with Fitter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fitter is a python library that utilizes the Scipy library's distribution to graph and calculate metrics for goodness of fit.\n",
    "\n",
    "For the source code, please refer to the [github repository](https://github.com/cokelaer/fitter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('data/sample_claims.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# f = Fitter(df)\n",
    "# f.fit()\n",
    "# f.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the Fitter API to see how well an exponential distribution fits the observed data.\n",
    "\n",
    "All calculations below can be replicated and validated using the Fitter source code provided in the github link above. The Fitter package abstracts all the logic away, but it  is good to understand the underlying concepts and implementation of the various tests using the exponential distribution as an example.\n",
    "\n",
    "Note that Fitter comes with some evaluation metrics for goodness of fit which will be covered later. These include:\n",
    "\n",
    "- sum of squared errors\n",
    "- information criterion\n",
    "- kullback-leibler information\n",
    "- kolmogorov-smirnov test\n",
    "- chi square test (not included)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scipy.stats.expon\n"
     ]
    },
    {
     "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>sumsquare_error</th>\n",
       "      <th>aic</th>\n",
       "      <th>bic</th>\n",
       "      <th>kl_div</th>\n",
       "      <th>ks_statistic</th>\n",
       "      <th>ks_pvalue</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>expon</th>\n",
       "      <td>8.627424e-08</td>\n",
       "      <td>3575.446841</td>\n",
       "      <td>-1990.559271</td>\n",
       "      <td>inf</td>\n",
       "      <td>0.18499</td>\n",
       "      <td>0.002388</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       sumsquare_error          aic          bic  kl_div  ks_statistic  \\\n",
       "expon     8.627424e-08  3575.446841 -1990.559271     inf       0.18499   \n",
       "\n",
       "       ks_pvalue  \n",
       "expon   0.002388  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD4CAYAAAAUymoqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfUklEQVR4nO3de5RU5Znv8e/TN5GL3JQWQW0kJLEF46UFNMZ0QnRAHZtMNIoaiXGGsNSZyTnJaHs0mVxHknPMxaNLBhMTTGYUYxLlRKIhmDIxCaB4R0QaRGnlIihiAS1093P+2LuxKLq79i66ui79+6y1V+3L++56Hxb4c19qb3N3RERE4ijL9wBERKT4KDxERCQ2hYeIiMSm8BARkdgUHiIiEltFvgfQEw4//HCvqanJqu/OnTsZMGBAzw4oj1RP4Su1mlRPYeuunhUrVmx19yOy2W9JhEdNTQ1PPvlkVn0TiQT19fU9O6A8Uj2Fr9RqUj2Frbt6zOzVbPer01YiIhKbwkNERGJTeIiISGwlcc1DRKQre/fupbm5mZaWlkjtBw8ezKpVq3I8qt4zePBgXnnlFUaPHk1lZWWP7VfhISIlrbm5mUGDBlFTU4OZZWz/7rvvMmjQoF4YWe/YsWMHe/bsobm5mTFjxvTYfnXaSkRKWktLC8OHD48UHKXIzBg+fHjkI6+oFB4iUvL6anB0yEX9Cg8REYlN4SEiIrHpgjlQ0/jQvvn1c87L40hERIqDjjxERHLsF7/4BRMnTuSkk07ii1/8IsuWLePEE0+kpaWFnTt3csIJJ/DCCy+QSCQ466yz+PSnP01tbS2zZ8+mvb0dgHvuuYcJEyYwfvx4rr/++n37HjhwIDfeeCMf+chHmDx5Mps3b+6VmnTkISJ9x9cHZ2yS1U26X3+ny02rVq1iwYIF/OUvf6GyspKrr76a1atXc8EFF3DTTTexe/duLr/8csaPH08ikWD58uW8+OKLHHvssUydOpVf//rXnHHGGVx//fWsWLGCoUOHcs455/DAAw8wffp0du7cyeTJk/nOd77Dddddx5133slNN92UTRWxKDxERHJoyZIlrFixgtNOOw2A3bt3M2LECL72ta9x2mmn0a9fP2699dZ97SdOnMhxxx0HwIwZM3j88ceprKykvr6eI44IHoB72WWX8ac//Ynp06dTVVXF+eefD8Cpp57K4sWLe6UuhYeI9B3dHCF06OkfCbo7M2fO5Oabb95v/aZNm0gmk+zdu5eWlpZ9j01Pv63WzHD3LvdfWVm5r095eTmtra09Nvbu6JqHiEgOTZkyhfvvv58tW7YA8NZbb/Hqq68ya9YsvvWtb3HZZZftdw1j+fLlvPLKK7S3t7NgwQLOPPNMJk2axGOPPcbWrVtpa2vjnnvu4eMf/3i+SgJ05CEiklO1tbV8+9vf5pxzzqG9vZ3KykoaGhqoqKjg0ksvpa2tjTPOOINHH32UsrIyTj/9dBobG3n++ef3XTwvKyvj5ptv5hOf+ATuzrnnnktDQ0Ne61J4iIjk2MUXX8zFF1/c6bby8nKWLVsGBC9u6t+/PwsWLDig3aWXXsqll156wPpkMrlv/sILL+TCCy/soVF3T6etREQkNh15iIgUiPr6+qJ5Ba6OPESk5HV3t1JfkIv6FR4iUtL69evHtm3b+myAuDvbtm2jX79+PbpfnbYSkZI2evRompubefPNNyO1b2lp6fH/0OZTS0sLQ4YMYfTo0T26X4WHiJS0ysrKWG/QSyQSnHzyyTkcUe/KVT06bSUiIrEpPEREJDaFh4iIxKbwEBGR2BQeIiISm8JDRERiU3iIiEhsCg8REYktUniY2VQzW21mTWbW2Ml2M7Nbw+3Pmdkpmfqa2TAzW2xma8LPoeH6GjPbbWbPhNPcnihURER6TsbwMLNy4HZgGlALzDCz2rRm04Bx4TQLuCNC30ZgibuPA5aEyx3WuvtJ4TQ72+JERCQ3ohx5TASa3H2du+8B7gXSX2HVANztgaXAEDMbmaFvAzA/nJ8PTD+4UkREpLdEebbVKGBDynIzMClCm1EZ+la7+0YAd99oZiNS2o0xs6eBHcBN7v7n9EGZ2SyCoxyqq6tJJBIRSjlQMpnkyxPa9i1nu59CkUwmi76GVKVWD5ReTaqnsOWqnijhYZ2sS3+2cVdtovRNtxE4xt23mdmpwANmdoK779hvJ+7zgHkAdXV1nu0LVBKJBLc8vnPf8vrLsttPoUgkEkXzMpkoSq0eKL2aVE9hy1U9UU5bNQNHpyyPBt6I2Ka7vpvDU1uEn1sA3P09d98Wzq8A1gIfjFKMiIj0jijh8QQwzszGmFkVcAmwMK3NQuCK8K6rycA74Smp7vouBGaG8zOBBwHM7IjwQjtmdhzBRfh1WVcoIiI9LuNpK3dvNbNrgUeAcuAud19pZrPD7XOBRcC5QBOwC7iyu77hrucA95nZVcBrwEXh+rOAb5pZK9AGzHb3t3qkWhER6RGRXgbl7osIAiJ13dyUeQeuido3XL8NmNLJ+l8Bv4oyLhERyQ/9wlxERGJTeIiISGwKDxERiU3hISIisSk8REQkNoWHiIjEpvAQEZHYFB4iIhKbwkNERGJTeIiISGwKDxERiU3hISIisSk8REQkNoWHiIjEpvAQEZHYFB4iIhKbwkNERGJTeIiISGwKDxERiU3hISIisSk8REQkNoWHiIjEpvAQEZHYFB4iIhKbwkNERGJTeIiISGwKDxERiS1SeJjZVDNbbWZNZtbYyXYzs1vD7c+Z2SmZ+prZMDNbbGZrws+hafs8xsySZvaVgylQRER6XsbwMLNy4HZgGlALzDCz2rRm04Bx4TQLuCNC30ZgibuPA5aEy6l+APwui5pERCTHohx5TASa3H2du+8B7gUa0to0AHd7YCkwxMxGZujbAMwP5+cD0zt2ZmbTgXXAyqyqEhGRnKqI0GYUsCFluRmYFKHNqAx9q919I4C7bzSzEQBmNgC4Hjgb6PKUlZnNIjjKobq6mkQiEaGUAyWTSb48oW3fcrb7KRTJZLLoa0hVavVA6dWkegpbruqJEh7WyTqP2CZK33TfAH7g7kmzzrqHO3GfB8wDqKur8/r6+gy77VwikeCWx3fuW15/WXb7KRSJRIJs/ywKUanVA6VXk+opbLmqJ0p4NANHpyyPBt6I2Kaqm76bzWxkeNQxEtgSrp8EXGhm3wOGAO1m1uLut0UYq4iI9IIo1zyeAMaZ2RgzqwIuARamtVkIXBHedTUZeCc8JdVd34XAzHB+JvAggLt/zN1r3L0G+CHwHwoOEZHCkvHIw91bzexa4BGgHLjL3Vea2exw+1xgEXAu0ATsAq7srm+46znAfWZ2FfAacFGPViYiIjkT5bQV7r6IICBS181NmXfgmqh9w/XbgCkZvvfrUcYnIiK9S78wFxGR2BQeIiISm8JDRERiU3iIiEhsCg8REYlN4SEiIrEpPEREJDaFh4iIxKbwEBGR2BQeIiISm8JDRERiU3iIiEhsCg8REYlN4SEiIrEpPEREJDaFh4iIxKbwEBGR2BQeIiISm8JDRERiU3iIiEhsCg8REYlN4SEiIrEpPEREJDaFh4iIxKbwEBGR2BQeIiISm8JDRERiU3h4O2PtderLnsn3SEREikak8DCzqWa22syazKyxk+1mZreG258zs1My9TWzYWa22MzWhJ9Dw/UTzeyZcHrWzD7dE4V2WZs7v6tq5GdV3+NQWnL5VSIiJSNjeJhZOXA7MA2oBWaYWW1as2nAuHCaBdwRoW8jsMTdxwFLwmWAF4A6dz8JmAr8p5lVZFtgJl5Wzqt+JABjbFOuvkZEpKREOfKYCDS5+zp33wPcCzSktWkA7vbAUmCImY3M0LcBmB/OzwemA7j7LndvDdf3Azy70qJb5yMBOM425vqrRERKQpT/ox8FbEhZbgYmRWgzKkPfanffCODuG81sREcjM5sE3AUcC3wuJUxIaTOL4CiH6upqEolEhFIOlEwmGXl4NbwNnzuyOev9FIpkMln0NaQqtXqg9GpSPYUtV/VECQ/rZF360UBXbaL0PbCB+zLgBDM7HphvZr9z95a0NvOAeQB1dXVeX1+fabedSiQS3LtlFKdWwqY3N9OQ5X4KRSKRINs/i0JUavVA6dWkegpbruqJctqqGTg6ZXk08EbENt313Rye2iL83JL+xe6+CtgJjI8wzqytbT8KgOMsvSwREelMlPB4AhhnZmPMrAq4BFiY1mYhcEV419Vk4J3wlFR3fRcCM8P5mcCDAGHbinD+WOBDwPpsC4yi45rHGNsEnvNLLCIiRS/jaSt3bzWza4FHgHLgLndfaWazw+1zgUXAuUATsAu4sru+4a7nAPeZ2VXAa8BF4fozgUYz2wu0A1e7+9YeqbYL2xnENh/EcHsX3t0Ihx2Vy68TESl6kW6BdfdFBAGRum5uyrwD10TtG67fBkzpZP3PgZ9HGVdPWucjg/DY+rLCQ0QkA/3CPLQuvO7B1jX5HYiISBFQeIQ6rnuwrSm/AxERKQIKj9C+8NCRh4hIRgqP0PtHHgoPEZFMFB6hV72avV4O2zfA3t35Ho6ISEFTeIRaqeA1HwE4bFub7+GIiBQ0hUeKjlNXV9+6gJrGh/I8GhGRwqXwSKGn64qIRKPwSLHOw2dclSk8RES6o/BIsba948hDD0gUEemOwiNFx5HHWNtIL7yDSkSkaCk8UrzFILb7AAbZbqp5O9/DEREpWAqP/RirPXj9yPFlr+V5LCIihUvhkebF9mMBON4UHiIiXVF4pHnRg/CoLVuf34GIiBQwhUeaF9trAKi1V/M7EBGRAqbwSLPGR7HXy4NX0u7Zme/hiIgUJIVHmj1U0uRHUWYOW1blezgiIgVJ4dGJjusebHouvwMRESlQCo9OdNxxxabn8zsQEZECpfDoxCpXeIiIdEfh0YlV7ccEM5tXQntbfgcjIlKAFB6deJvDeMOHwd5d8NYr+R6OiEjBUXh04f3rHrpoLiKSTuHRBV33EBHpmsKjC7rjSkSkawqPLryoIw8RkS5FCg8zm2pmq82sycwaO9luZnZruP05MzslU18zG2Zmi81sTfg5NFx/tpmtMLPnw89P9kShcb3mI6BqICQ3QfLNfAxBRKRgZQwPMysHbgemAbXADDOrTWs2DRgXTrOAOyL0bQSWuPs4YEm4DLAV+Ht3nwDMBH6edXUHwSmDkR8JFl5/Mh9DEBEpWFGOPCYCTe6+zt33APcCDWltGoC7PbAUGGJmIzP0bQDmh/PzgekA7v60u3e8RHwl0M/MDsmuvIN09MTgc8OyvHy9iEihqojQZhSwIWW5GZgUoc2oDH2r3X0jgLtvNLMRnXz3Z4Cn3f299A1mNovgKIfq6moSiUSEUg6UTCb58oTOfwj4/Pb+TAC2P/97nqmoz2r/vS2ZTGb9Z1GISq0eKL2aVE9hy1U9UcLDOlnnEdtE6dv5l5qdAHwXOKez7e4+D5gHUFdX5/X19VF2e4BEIsEtj3f+6PV//upV8MJ3GLJzLfUf+yiUV2b1Hb0pkUiQ7Z9FISq1eqD0alI9hS1X9UQ5bdUMHJ2yPBp4I2Kb7vpuDk9tEX5u6WhkZqOB3wBXuPvaCGPMiZpvLWVt+0hobdGPBUVEUkQJjyeAcWY2xsyqgEuAhWltFgJXhHddTQbeCU9Jddd3IcEFccLPBwHMbAjwEHCDu/8l+9J6xor2DwYzG5bndyAiIgUkY3i4eytwLfAIsAq4z91XmtlsM5sdNlsErAOagDuBq7vrG/aZA5xtZmuAs8NlwvYfAL5qZs+EU2fXQ3rFCu8ID100FxHpEOWaB+6+iCAgUtfNTZl34JqofcP124Apnaz/NvDtKOPqDSvaxwUzOvIQEdlHvzDPYK0fxXYfADteh+0bMncQEekDFB4ZOGU8te/oQ6euRERA4RGJLpqLiOxP4RHBU64jDxGRVAqPCJ5tHwtWHjxhd0/nPygUEelLFB4R7KIfHDkBvA1eW5rv4YiI5J3CI6qxnwg+1z6a33GIiBQAhUdUY8PXiig8REQUHpEdPRkqB8CWF2FH+qO9RET6FoVHVBVVMOZjwbyOPkSkj1N4xDE2fJpK05L8jkNEJM8UHnF8IAyPdX+E9s5fICUi0hcoPOIYdhwMOQZ2vw0bn8n3aERE8kbhEYdZyqkrXfcQkb5L4RFXx6mrtbruISJ9V6T3eQjUND4EwCBaeP7Q8uAhiS3vUPP1x/e1WT/nvHwNT0SkV+nII6Z36Q9HTwweVbJmcb6HIyKSFwqPbNQ2BJ8rf5PfcYiI5InCIxsd4bFmMQPYnd+xiIjkgcIjG4cdBcecDm3vMaXsqXyPRkSk1yk8snXCPwBwfrke0S4ifY/CI1u1FwDGx8ueZSC78j0aEZFepfDIQk3jQ9R8ZwV/azueQ6yVs8tW5HtIIiK9SuFxEB5qnwTAeTp1JSJ9jMLjIDzcNpE2N84qe47D0LvNRaTvUHgchK0M5m/ttVRZGxeU/zXfwxER6TUKj4O0oC14t/ll5X8A9zyPRkSkd0QKDzObamarzazJzBo72W5mdmu4/TkzOyVTXzMbZmaLzWxN+Dk0XD/czP5oZkkzu60nisylh9sn8qYfxvFlG2DDsnwPR0SkV2QMDzMrB24HpgG1wAwzq01rNg0YF06zgDsi9G0Elrj7OGBJuAzQAnwV+Er2ZfWevVRwX1t9sPDET/I6FhGR3hLlyGMi0OTu69x9D3Av0JDWpgG42wNLgSFmNjJD3wZgfjg/H5gO4O473f1xghApCve0TaHdDV58AHZuzfdwRERyLsoj2UcBG1KWm4FJEdqMytC32t03Arj7RjMbEWPcmNksgqMcqqurSSQScbrvk0wm+fKEg32l7FBeff1Exux6lrW/+iYbjvmHg9xf9pLJZNZ/FoWo1OqB0qtJ9RS2XNUTJTysk3XpV4a7ahOlb1bcfR4wD6Curs7r6+uz2k8ikeCWxw/+Ntunys7hp1XPMvatBGMv/yGU5edehEQiQbZ/FoWo1OqB0qtJ9RS2XNUT5b9wzcDRKcujgTcitumu7+bw1Bbh55bowy48j7V/BAYfA9tfhTWP5Hs4IiI5FSU8ngDGmdkYM6sCLgEWprVZCFwR3nU1GXgnPCXVXd+FwMxwfibw4EHWklftlMGkLwYLj31Xt+2KSEnLGB7u3gpcCzwCrALuc/eVZjbbzGaHzRYB64Am4E7g6u76hn3mAGeb2Rrg7HAZADNbD3wf+LyZNXdyd1dhqvsCDDgC3nga1vw+36MREcmZSO8wd/dFBAGRum5uyrwD10TtG67fBkzpok9NlHEVnKr+8NEvwe9vhMTNMO4csM4u+4iIFDf9wryn1X0BBozQ0YeIlDSFR0+r6g9nfimYT9ysax8iUpIUHrlw6pXvH3289Nt8j0ZEpMcpPHKhqj+c9W/B/MP/C/boce0iUloUHrlS9wU4cgK88xr86f/kezQiIj1K4ZEr5RVw3g+C+b/+X3hzdX7HIyLSgxQeuXT0aXDKTGjfCw99WRfPRaRkKDxy7VNfh0OHwfo/w9O/yPdoRER6hMIj1/oPg6k3B/O/u06nr0SkJCg8esOJFwfT3l3wy8/D3t35HpGIyEFRePQGMzjvFhj+AdjyIjx8wJt8RUSKisKjtxwyCC76GZQfAit+pusfIlLUFB696cgJMO27wfzCf4GX9ewrESlOCo/eVnclnPk/wNvglzOh+cl8j0hEJDaFRz5M+Xc46fLgAvp/XQRbVuV7RCIisSg88sEM/v5HMO7vYPdb8NNpsGF5vkclIhKZwiNfyivgs/PhQ+fC7rdh/gXwst59LiLFQeGRT5WHwmd/DidfDq274Z4ZsPxOPcZERApepNfQysGpaXxo3/z6Oeftv7G8Ai64LXj/x+Pfh0Vf4aH/dx+Ne2fxLv0PbC8iUgB05FEIzOBT/w4X3gVVgzivfDkPVd1Anb2U75GJiHRK4VFIxn8GvvgYz7fXcEzZm9x/yDfhwWtg57Z8j0xEZD8Kj0IzfCyf2fMNftT6ad7ziuCX6LedCkvv0DOxRKRgKDwK0B4q+UHrRUzbMwfGfDy4G+vhRvjRSbB0rkJERPJO4VHA1vlRcMWDcMk9cOSJkNwED18Pt3wYHrkRtq3N9xBFpI/S3VY9KPWuqoNpsx8z+PC58KFp/NON3+Saigc4qWUd/O02+NttLGv/MJPO/0eobYCBI7Icec/o9q6yGG3ypZDHJpKqEP6uKjyKhRmL2+tYvKeOE20tl5f/gQvK/8qkspdg0VeCF02Nnsix5cdB8yA46iQoK8/3qEWkRCk8itBzPpbrWsfyzdbP8amyp/jh+HXQtAQ2LGUMS+HH/w1Vg2B0HRwzGUadGjzRd2B1cCQjInKQFB5FLEl/Hmg/kx9eejO07IBX/sTrf/45o3atgu2vwro/BlOHAUfAER+Gwz8Ih4+DYWNhyDHBVNU/f4WISNGJFB5mNhX4EVAO/Njd56Rtt3D7ucAu4PPu/lR3fc1sGLAAqAHWA59197fDbTcAVwFtwL+4ux76lEm/w+D481mzeSCj6uthx0bYsCx44OKm54Jp55vBtP7PB/bvPxwGjQyn6iBo+h8OAw6HQ4dCvyFw6BA45LDguyr76yhGpA/LGB5mVg7cDpwNNANPmNlCd38xpdk0YFw4TQLuACZl6NsILHH3OWbWGC5fb2a1wCXACcBRwB/M7IPu3tYzJfcRh42EE6YHEwTPy9r+Gmx9OZzWwNuvBOu2b4Bd24Jp8wvR9m9lUDUQqgYEU+Wh3F/1Hi1eyXtUwYJ7oeKQ4M2JFVXBZ3kF/1bxKnupoNXL4fGXoawyuDZTVhHss6wcrHzf54jNL8EL24JtqROWsmxhkHXzCRnnT7GX369vwxEdhabVnbrcRXhmyNRBO9bA64d13+hg9WKwD3y3Cd4Y3Gvfl2vFUM94WwfA63543sYQ5chjItDk7usAzOxeoAFIDY8G4G53d2CpmQ0xs5EERxVd9W0A6sP+84EEcH24/l53fw94xcyawjH8LfsyBTMYemwwjTt7/23tbbBzK7z7RnDEsnNLcISSfDN4ZPzu7cFvTVq2w3vvBqfIWnfDezuCKVSXeuP3qqc7HcY1qX/j/vDLjMOuBeil1538+pCUhZ/k7ntOBXgqd/vvbXUAK/I9ip5TDPX8Nvy7et3efwJm5GUMUcJjFLAhZbmZ4OgiU5tRGfpWu/tGAHffaGYd95mOApZ2sq/9mNksYFa4mDSz1RFq6czhwNYs++acfTd2m4KuJwulVg+UXk2qJ29u4X//xy2ZGnVXz7HZfnOU8Ojs+Df9meFdtYnSN5vvw93nAfMy7CsjM3vS3esOdj+FQvUUvlKrSfUUtlzVE+UX5s3A0SnLo4E3Irbpru/m8NQW4eeWGN8nIiJ5FCU8ngDGmdkYM6siuJi9MK3NQuAKC0wG3glPSXXXdyEwM5yfCTyYsv4SMzvEzMYQXITXO1pFRApIxtNW7t5qZtcCjxDcbnuXu680s9nh9rnAIoLbdJsIbtW9sru+4a7nAPeZ2VXAa8BFYZ+VZnYfwUX1VuCaHN9pddCnvgqM6il8pVaT6ilsOanHXK88FRGRmPRUXRERiU3hISIisfXp8DCzqWa22syawl+5Fwwzu8vMtpjZCynrhpnZYjNbE34OTdl2Q1jHajP7u5T1p5rZ8+G2W8NHyRDekLAgXL/MzGpyWMvRZvZHM1tlZivN7F+LvJ5+ZrbczJ4N6/lGMdeTMpZyM3vazH5b7PWY2fpwHM+Y2ZMlUM8QM7vfzF4K/x2dnvd63L1PTgQX8NcCxwFVwLNAbb7HlTK+s4BTgBdS1n0PaAznG4HvhvO14fgPAcaEdZWH25YDpxP8fuZ3wLRw/dXA3HD+EmBBDmsZCZwSzg8CXg7HXKz1GDAwnK8ElgGTi7WelLr+J/DfwG+L+e9b+B3rgcPT1hVzPfOBfwznq4Ah+a4np38ZC3kK/wAfSVm+Abgh3+NKG2MN+4fHamBkOD8SWN3Z2Anubjs9bPNSyvoZwH+mtgnnKwh+gWq9VNeDBM87K/p6gP4EDxuZVMz1EPyeagnwSd4Pj2KuZz0HhkdR1gMcBrySvv9819OXT1t19UiVQrbfI12A1Ee6dPV4mOZO1u/Xx91bgXeA4TkbeSg8HD6Z4P/Wi7ae8BTPMwQ/bl3s7kVdD/BD4DqgPWVdMdfjwO/NbIUFjzKC4q3nOOBN4KfhacUfm9kA8lxPXw6PbB6dUqiyeTxMr9dvZgOBXwFfcvcd3TXtZF1B1ePube5+EsH/sU80s/HdNC/oeszsfGCLu0d9HGBB1xP6qLufQvDE72vM7Kxu2hZ6PRUEp7DvcPeTgZ0Ep6m60iv19OXwKMbHoMR9pEtzOJ++fr8+ZlYBDAbeytXAzaySIDj+y91/Ha4u2no6uPt2gidCT6V46/kocIGZrQfuBT5pZr+geOvB3d8IP7cAvyF4Mnex1tMMNIdHtwD3E4RJXuvpy+ER5bErhSbWI13CQ9l3zWxyeFfFFWl9OvZ1IfCohyc8e1r43T8BVrn790ugniPMbEg4fyjwKeClYq3H3W9w99HuXkPw7+BRd7+8WOsxswFmNqhjHjgHeKFY63H3TcAGM/tQuGoKwRM48ltPri5YFcNE8EiVlwnuRrgx3+NJG9s9wEZgL8H/FVxFcA5yCbAm/ByW0v7GsI7VhHdQhOvrCP7hrAVu4/2nCvQDfknwSJnlwHE5rOVMgkPg54BnwuncIq7nRODpsJ4XgK+F64uynrTa6nn/gnlR1kNwjeDZcFrZ8W+7WOsJv+8k4Mnw79wDwNB816PHk4iISGx9+bSViIhkSeEhIiKxKTxERCQ2hYeIiMSm8BARkdgUHiIiEpvCQ0REYvv/cRZ+ERAfp8gAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "exp_fit = Fitter(df, distributions=[\"expon\"])\n",
    "exp_fit.fit()\n",
    "exp_fit.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy histogram function returns a tuple of the observed frequencies and the edges of the bin\n",
    "\n",
    "We then adjust the x values to the mid-point of the edges using a list comprehension\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "y, x = np.histogram(df, bins=100, density=True)\n",
    "x = [(bin_edge_value + x[i + 1]) / 2.0 for i, bin_edge_value in enumerate(x[0:-1])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we've defined our bin values and frequencies, we proceed to fit the data with an exponential distribution using the [Scipy library](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.expon.html). The underlying object [rv_continuous.fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.fit.html#scipy.stats.rv_continuous.fit) utilizes the maximum likelihood estimate (MLE) algorithm detailed in the previous notebook as a default, but method of moments (MOM) is available as well.\n",
    "\n",
    "The fit method returns a tuple containing general optimized location and scale parameters, which depending on the chosen distribution, can be solved in a system of equations to obtain the distributional parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "location parameter: 24.0\n",
      "scale parameter: 2965.8333333333335\n",
      "lambda: 0.00033717336330429894\n"
     ]
    }
   ],
   "source": [
    "params = st.expon.fit(df)\n",
    "\n",
    "fitted_pdf = st.expon.pdf(x, *params)\n",
    "\n",
    "print(f'location parameter: {params[0]}')\n",
    "print(f'scale parameter: {params[-1]}')\n",
    "print(f'lambda: {1 / params[-1]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the exponential distribution is parameterized as follows:\n",
    "\n",
    "$ f_{x}(X) = \\lambda e^{-\\lambda x} $\n",
    "\n",
    "Then following from notebook 2, the `scale` parameter would be $ \\dfrac{1}{\\beta} $ and given our choice of fit, we can say that we think the underlying distribution generating the data follows $~Exp(\\lambda = 0.00034)$.\n",
    "\n",
    "If we subscribe to this assumption, we can start to make very powerful inferences about the data as follows below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "probability of observing a claim < $5000 is: 0.81\n",
      "probability of observing a claim < $3000 is: 0.63\n",
      "probability of observing a claim between $3000 and $5000 is: 0.18\n"
     ]
    }
   ],
   "source": [
    "print(f'probability of observing a claim < $5000 is: {np.round(st.expon.cdf(5000, *params), 2)}')\n",
    "print(f'probability of observing a claim < $3000 is: {np.round(st.expon.cdf(3000, *params), 2)}')\n",
    "print(f'probability of observing a claim between $3000 and $5000 is: {np.round(st.expon.cdf(5000, *params) - np.round(st.expon.cdf(3000, *params), 2), 2)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1ff963c99d0>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD4CAYAAAAUymoqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZ50lEQVR4nO3df5RX9X3n8eeLGQYDEREZEBAdVDRC0qohiEn12BgTME0hacjRdgPr2sOy0bNpc9JkbE+2m3a3JdndJrWxopvaRRtjWI11tuK6Btv0l6hD/AmKDIgwgjDiT1DBYd77x/2gX8dh5t75db9feD3O+Z577+d+Pvf7eSfgi3u/93u/igjMzMyKGFH2BMzMrPY4PMzMrDCHh5mZFebwMDOzwhweZmZWWH3ZExgMEyZMiKamprKnYWZWU9atW/diRDT2Z+wRER5NTU20traWPQ0zs5oi6bn+jvVlKzMzK8zhYWZmhTk8zMysMIeHmZkV5vAwM7PCHB5mZlaYw8PMzApzeJiZWWEODzMzK+yI+Ib5QDU13/3O+tblny1xJmZmtcFnHmZmVpjDw8zMCnN4mJlZYQ4PMzMrzOFhZmaFOTzMzKwwh4eZmRXm8DAzs8IcHmZmVpjDw8zMCnN4mJlZYQ4PMzMrzOFhZmaFOTzMzKwwh4eZmRXm8DAzs8IcHmZmVpjDw8zMCnN4mJlZYQ4PMzMrzOFhZmaF5QoPSfMkbZTUJqm5h/2SdG3a/7ikc/saK2m8pPskbUrL41N7k6Q3JT2aXisGo1AzMxs8fYaHpDrgOmA+MBO4XNLMbt3mAzPSaylwfY6xzcCaiJgBrEnbh2yOiLPTa1l/izMzs6GR58xjDtAWEVsi4gBwG7CgW58FwM2RWQuMkzS5j7ELgJVpfSWwcGClmJnZcMkTHlOB7RXb7aktT5/exk6KiJ0AaTmxot90SY9I+rmkC3qalKSlkloltXZ0dOQow8zMBkue8FAPbZGzT56x3e0ETo6Ic4CvAbdKGvu+g0TcGBGzI2J2Y2NjH4c0M7PBlCc82oFpFdsnATty9ult7K50aYu03A0QEfsjYk9aXwdsBs7IU4yZmQ2PPOHxMDBD0nRJDcBlQEu3Pi3A4nTX1Vzg1XQpqrexLcCStL4EuAtAUmP6oB1Jp5J9CL+l3xWamdmgq++rQ0R0SroauBeoA26KiPWSlqX9K4DVwKVAG/AGcEVvY9OhlwOrJF0JbAMWpfYLgT+S1AkcBJZFxEuDUq2ZmQ2KPsMDICJWkwVEZduKivUArso7NrXvAS7uof0O4I488zIzs3L4G+ZmZlaYw8PMzApzeJiZWWEODzMzK8zhYWZmhTk8zMysMIeHmZkV5vAwM7PCHB5mZlaYw8PMzApzeJiZWWEODzMzK8zhYWZmhTk8zMysMIeHmZkV5vAwM7PCHB5mZlaYw8PMzApzeJiZWWEODzMzK8zhYWZmhTk8zMysMIeHmZkV5vAwM7PCHB5mZlaYw8PMzApzeJiZWWG5wkPSPEkbJbVJau5hvyRdm/Y/LuncvsZKGi/pPkmb0vL4bsc8WdJeSV8fSIFmZjb4+gwPSXXAdcB8YCZwuaSZ3brNB2ak11Lg+hxjm4E1ETEDWJO2K30PuKcfNZmZ2RDLc+YxB2iLiC0RcQC4DVjQrc8C4ObIrAXGSZrcx9gFwMq0vhJYeOhgkhYCW4D1/arKzMyGVJ7wmApsr9huT215+vQ2dlJE7ARIy4kAksYA3wS+3dukJC2V1CqptaOjI0cZZmY2WPKEh3poi5x98ozt7tvA9yJib2+dIuLGiJgdEbMbGxv7OKSZmQ2m+hx92oFpFdsnATty9mnoZewuSZMjYme6xLU7tZ8HfFHSd4FxQJektyLiBznmamZmwyDPmcfDwAxJ0yU1AJcBLd36tACL011Xc4FX06Wo3sa2AEvS+hLgLoCIuCAimiKiCfg+8CcODjOz6tLnmUdEdEq6GrgXqANuioj1kpal/SuA1cClQBvwBnBFb2PToZcDqyRdCWwDFg1qZWZmNmTyXLYiIlaTBURl24qK9QCuyjs2te8BLu7jff9znvmZmdnw8jfMzcysMIeHmZkV5vAwM7PCHB5mZlaYw8PMzApzeJiZWWEODzMzK8zhYWZmhTk8zMysMIeHmZkV5vAwM7PCHB5mZlaYw8PMzApzeJiZWWEODzMzK8zhYWZmhTk8zMysMIeHmZkV5vAwM7PCHB5mZlaYw8PMzApzeJiZWWEODzMzK8zhYWZmhTk8zMysMIeHmZkV5vAwM7PCHB5JA2+XPQUzs5qRKzwkzZO0UVKbpOYe9kvStWn/45LO7WuspPGS7pO0KS2PT+1zJD2aXo9J+vxgFNqbvxz5fX7S8MdD/TZmZkeMPsNDUh1wHTAfmAlcLmlmt27zgRnptRS4PsfYZmBNRMwA1qRtgCeB2RFxNjAPuEFSfX8LzGNPjOVU7QBiKN/GzOyIkefMYw7QFhFbIuIAcBuwoFufBcDNkVkLjJM0uY+xC4CVaX0lsBAgIt6IiM7UfgzD8F/0LTGZ4/QGJ/DaUL+VmdkRIU94TAW2V2y3p7Y8fXobOykidgKk5cRDnSSdJ2k98ASwrCJMqOizVFKrpNaOjo4cZRze5pgCwGnaMaDjmJkdLfKEh3po6342cLg+eca+v0PEgxExC/gYcI2kY3roc2NEzI6I2Y2NjX0dsldbUnicOmLngI5jZna0yBMe7cC0iu2TgO7/RD9cn97G7kqXtkjL3d3fOCKeAvYBH84xz357Pk7grRjJqXJ4mJnlkSc8HgZmSJouqQG4DGjp1qcFWJzuupoLvJouRfU2tgVYktaXAHcBpL71af0U4Exga38LzCMYwbNxYvrQ3MzM+tLnXUwR0SnpauBeoA64KSLWS1qW9q8AVgOXAm3AG8AVvY1Nh14OrJJ0JbANWJTafwVolvQ20AV8JSJeHJRqe7ElJnOWtg3125iZHRFy3QIbEavJAqKybUXFegBX5R2b2vcAF/fQfgtwS555DabNMYXPjGiFzgNQ3zDcb29mVlP8DfNkS9cU6tUFLz9b9lTMzKqewyPZEpOzlRc3lTsRM7Ma4PBI3gmPPQ4PM7O+ODySvYxmd4yDF9vKnoqZWdVzeFTYEpPhxWfKnoaZWdVzeFTY3DXFl63MzHJweFTYEpPhzZc5p/nHNDXfXfZ0zMyqlsOjwub0obkfU2Jm1juHRwU/INHMLB+HR4X2aGR/1PvR7GZmfXB4VOhiBM/FJE7zZSszs145PLrZHFP8dF0zsz44PLrZGNOYrhf4AG+VPRUzs6rl8OhmQ9cpjFBwptrLnoqZWdVyeHSzIZoAmDniuXInYmZWxRwe3bTHBF6L0czU1rKnYmZWtRwe7yOeipM5a4R/VdDM7HAcHj3Y0HUKH9I26DpY9lTMzKqSw6MHG+IUxmg/vORfFTQz64nDowcbupqylRceL3UeZmbVyuHRg00xlbejDl54ouypmJlVJYdHDw4wkraYArueLHsqZmZVyeFxGBviFJ95mJkdhsPjMDZ0nQKv74S9HWVPxcys6jg8DuPQN83Z5bMPM7PuHB6H8VTXydnKC/7cw8ysO4fHYbzCsTB2qj/3MDPrQa7wkDRP0kZJbZKae9gvSdem/Y9LOrevsZLGS7pP0qa0PD61XyJpnaQn0vKTg1Fov5z4EX/Xw8ysB32Gh6Q64DpgPjATuFzSzG7d5gMz0mspcH2Osc3AmoiYAaxJ2wAvAp+LiI8AS4Bb+l3dQE05Bzo2wluvlTYFM7NqlOfMYw7QFhFbIuIAcBuwoFufBcDNkVkLjJM0uY+xC4CVaX0lsBAgIh6JiEM/5bceOEbSqP6VN0DT5gABz7eW8vZmZtUqT3hMBbZXbLentjx9ehs7KSJ2AqTlxB7e+zeARyJif/cdkpZKapXU2tExRLfTTp0NCLY/NDTHNzOrUXnCQz20Rc4+ecb2/KbSLOA7wL/vaX9E3BgRsyNidmNjY55DFnfMWJg0C7Y/ODTHNzOrUXnCox2YVrF9ErAjZ5/exu5Kl7ZIy92HOkk6CbgTWBwRm3PMcUg0Nd/NLc+fyGttD/jx7GZmFfKEx8PADEnTJTUAlwEt3fq0AIvTXVdzgVfTpajexraQfSBOWt4FIGkccDdwTUT8S/9LGxzrus5grN6EjqfLnoqZWdXoMzwiohO4GrgXeApYFRHrJS2TtCx1Ww1sAdqA/wl8pbexacxy4BJJm4BL0jap/+nAtyQ9ml49fR4yLNbFjGxl29qypmBmVnUUkesjiKo2e/bsaG3t/x1RTc1397I3eHjUV2j85fnwhRv6/R5mZtVG0rqImN2fsf6GeZ/Euq4z/KG5mVkFh0cO67pmwMvPwt7dfXc2MzsKODxyWNd1Rrbi73uYmQEOj1yejOlQ1+BLV2ZmicMjhwOMzJ5z5TMPMzPA4ZHftPNgxy/gwL6yZ2JmVjqHR16n/SocPABbS//eoplZ6RweeZ38caj/AGxeU/ZMzMxK5/DIa+Qx0PQJaHN4mJk5PIo47ZOwZxO8sq3smZiZlcrhUcRpF2fLzfeXOw8zs5I5PIpoPBPGTvWlKzM76jk8ipCyS1dbfg4HO8uejZlZaRweRZ1+Mex/FZ5fV/ZMzMxK4/Ao6tSLQCN8y66ZHdXqy55Araj8zY+tp380+9zjV3//ve3LP1vG1MzMhp3PPPrj9E9ll61e31X2TMzMSuHw6I+ZC4CAp7r/lLuZ2dHB4dEfE8+Cxg/B+r8teyZmZqVwePTXrC/Ac//CRF4ueyZmZsPO4dFfsxYCwfw6/8aHmR19HB791XgmTJzFZ+vWlj0TM7Nh51t1++HQ7blX1c3i90auYhIvsYvxJc/KzGz4+MxjAFZ3nQfApXX+bXMzO7o4PAbg2ZjM+q5T+DVfujKzo4zDY4D+z8Hz+eiITTRpZ9lTMTMbNg6PAbrj4IW8HXX8Vp2fdWVmR49c4SFpnqSNktokNfewX5KuTfsfl3RuX2MljZd0n6RNaXl8aj9B0t9L2ivpB4NR5FDqYBz3dn2MRXU/h7ffLHs6ZmbDos/wkFQHXAfMB2YCl0ua2a3bfGBGei0Frs8xthlYExEzgDVpG+At4FvA1/tf1vD6m4OfYpz2wfo7y56KmdmwyHPmMQdoi4gtEXEAuA1Y0K3PAuDmyKwFxkma3MfYBcDKtL4SWAgQEfsi4p/JQqQmrO06i01dU+Hhvyp7KmZmwyJPeEwFtldst6e2PH16GzspInYCpOXE/NMGSUsltUpq7ejoKDJ0CIgfHbwYnm+FHY+WPBczs6GXJzzUQ1vk7JNnbL9ExI0RMTsiZjc2Ng7GIQfkpwcvgJGjodVnH2Z25MsTHu3AtIrtk4AdOfv0NnZXurRFWu7OP+3q8xpj4MO/AU/cDm+8VPZ0zMyGVJ7weBiYIWm6pAbgMqD7D1m0AIvTXVdzgVfTpajexrYAS9L6EuCuAdZSvvOvyu64euC6smdiZjak+ny2VUR0SroauBeoA26KiPWSlqX9K4DVwKVAG/AGcEVvY9OhlwOrJF0JbAMWHXpPSVuBsUCDpIXApyNiw8DLHWITz8qetvvgDVmQjPbzrszsyJTrwYgRsZosICrbVlSsB3BV3rGpfQ9w8WHGNOWZV1W68BvZj0Q9cB1c/K2yZ2NmNiT8DfPBNmnmu2cf/uzDzI5QDo+hcOE34MBef/ZhZkcsh8dQqDz72Fv2d1DMzAafw2OoXPT70PkW/OwPy56Jmdmgc3gMlcYz4ONXw6M/guceKHs2ZmaDyuExlC78PThuGtz9NTj4dtmzMTMbNA6PodQwBuZ/B3ZvyD7/MDM7Qjg8htqZl8IZ8+Af/hRe3lr2bMzMBoXDY6hJcOl/A9XB7f8OOg+UPSMzswFzeAyHcSfDr18Lz6+D+/+o7NmYmQ2Yw2O4zFoIH/tt+Ne/gGfuLXs2ZmYD4vAYTp/+rzDpI3DnMn/+YWY1zeExnEYeA4v+F0QX3PIF2Pdi2TMyM+sXh8dwm3A6/OYqeG0H/GgR7N9b9ozMzApzeJTh5PNg0V/Dzsdg1WLfgWVmNcfhUZYz58Pn/hw2r4FbF8H+18uekZlZbg6PMp37ZVjwl/DsP8HKz/kJvGZWMxweZTvnt+CyW2H303DTp6HjmbJnZGbWp1w/Q2sD09R89zvrW5d/9v0dzpwHi++C2y6HGy/id/ct5s6uCw7f38ysZD7zqBYnnwfL/hmmnM33Gq7nu/U3MJq3yp6VmVmPHB7VZOwUWNzCtZ0L+WLdP/KzUV+HDS0QUfbMzMzew+FRberq+bPOL/HFA3/IK3EsrPoy3Pol6NhY9szMzN7h8KhSv4gz+NyB/wKf+RN47l/huvPg9iv9gbqZVQWHRxU7SB2cfxV89XH4xFdh4z1w3Rz48eXQ9jPo6ip7imZ2lPLdVoOo8q6qgfR5nzEnwCXf5tw1H+KK+v/L5U/fz4SNq9naNYmfHryAr/1uM0yY0Y8ZD64+7yrL2acs1Tw3s0rV8GfV4VFDXmIs/6PzS/xF5+eZN+JhfrN+Db9Tfwf84HaY9OHsVwtP/xRM/SjU+f9aMxs6/i9MDTrASFq6Pk7LgY8ziZd4cOFeWH8n/NN/h3/8LhxzHJz8cZg2B06eC5N/Ofs9dTOzQeLwqHG7GA9zvwxz/wO88RI8+3NoWwPb1sIz96ReghNOhxM/AhPPyi5xTTgDjp8ODaNLnb+Z1aZc4SFpHvDnQB3ww4hY3m2/0v5LgTeAfxsRv+htrKTxwE+AJmAr8KWIeDntuwa4EjgI/MeI8E/v5TF6PMz6fPYC2LcH2h/Knt77whPwfCus/+l7x4xpzH4m97iT4NjJ6XUijJmQ7Rt9AnzgeBg5Ovs9djMzcoSHpDrgOuASoB14WFJLRGyo6DYfmJFe5wHXA+f1MbYZWBMRyyU1p+1vSpoJXAbMAqYAP5N0RkQcHJySjyJjTsie3nvm/HfbDuyDPZvhxWeyXzN85Tl4+TnYtQHa7ocDh3m6b11Ddjls1Fg4ZiyMOhYaPpheo2HkaL5e/zxvRQP7GQkPPg/1DVA36t1l3UguGPE4b1NPZ4yA9kkwoh5G1GVL1WXrGpGWlesj0kvZEh2mTe9fOvTMBl2eM485QFtEbAGQdBuwAKgMjwXAzRERwFpJ4yRNJjurONzYBcBFafxK4B+Ab6b22yJiP/CspLY0hwf6X6a9o2EMTP6l7NWTt16DfR3Zrxzu64A3XoQ3X4E3X85eB/Zmffa/Dm9uz8Jo/17ofItldXupV7p9+J5bezz8LQ0VGz8c1MpySkFyKFwq1p8Zlc09EPzxiPf2f2d45fZhQilXWA1xoDkwj2hPjuoEYP6BPy1tDnnCYyqwvWK7nezsoq8+U/sYOykidgJExE5JEyuOtbaHY72HpKXA0rS5V1J/v4I9Aaja34PVdwr3qep6+uFIqweOvJpcT2muRMuv7KtTb/Wc0t93zhMePf0TpvvDlg7XJ8/Y/rwfEXEjcGMfx+qTpNaImD3Q41QL11P9jrSaXE91G6p68nzDvB2YVrF9ErAjZ5/exu5Kl7ZIy90F3s/MzEqUJzweBmZImi6pgezD7JZufVqAxcrMBV5Nl6R6G9sCLEnrS4C7KtovkzRK0nSyD+Ef6md9ZmY2BPq8bBURnZKuBu4lu932pohYL2lZ2r8CWE12m24b2a26V/Q2Nh16ObBK0pXANmBRGrNe0iqyD9U7gauG+E6rAV/6qjKup/odaTW5nuo2JPUo/FsRZmZWkJ+qa2ZmhTk8zMyssKM6PCTNk7RRUlv6lnvVkHSTpN2SnqxoGy/pPkmb0vL4in3XpDo2SvpMRftHJT2R9l2bHiVDuiHhJ6n9QUlNQ1jLNEl/L+kpSeslfbXG6zlG0kOSHkv1fLuW66mYS52kRyT9Xa3XI2lrmsejklqPgHrGSbpd0tPp79H5pdcTEUfli+wD/M3AqUAD8Bgws+x5VczvQuBc4MmKtu8CzWm9GfhOWp+Z5j8KmJ7qqkv7HgLOJ/v+zD3A/NT+FWBFWr8M+MkQ1jIZODetHws8k+Zcq/UI+GBaHwk8CMyt1Xoq6voacCvwd7X85y29x1ZgQre2Wq5nJfDbab0BGFd2PUP6h7GaX+l/wHsrtq8Bril7Xt3m2MR7w2MjMDmtTwY29jR3srvbzk99nq5ovxy4obJPWq8n+waqhqmuu8ied1bz9QCjgV+QPTmhZush+z7VGuCTvBsetVzPVt4fHjVZDzAWeLb78cuu52i+bHW4R6pUs/c80gWofKTL4R4P095D+3vGREQn8CpwwpDNPEmnw+eQ/Wu9ZutJl3geJfty630RUdP1AN8HvgFU/rZxLdcTwP+TtE7Zo4ygdus5FegA/jpdVvyhpDGUXM/RHB79eXRKterP42GGvX5JHwTuAH4nIl7rrWsPbVVVT0QcjIizyf7FPkfSh3vpXtX1SPo1YHdErMs7pIe2qqkn+UREnEv2xO+rJF3YS99qr6ee7BL29RFxDrCP7DLV4QxLPUdzeNTiY1CKPtKlPa13b3/PGEn1wHHAS0M1cUkjyYLjRxFx6EdFaraeQyLiFbInQs+jduv5BPDrkrYCtwGflPQ31G49RMSOtNwN3En2ZO5aracdaE9ntwC3k4VJqfUczeGR57Er1abQI13SqezrkuamuyoWdxtz6FhfBO6PdMFzsKX3/ivgqYj4syOgnkZJ49L6B4BPAU/Xaj0RcU1EnBQRTWR/D+6PiH9Tq/VIGiPp2EPrwKeBJ2u1noh4Adgu6czUdDHZEzjKrWeoPrCqhRfZI1WeIbsb4Q/Knk+3uf0Y2Am8TfavgivJrkGuATal5fiK/n+Q6thIuoMitc8m+4uzGfgB7z5V4Bjgf5M9UuYh4NQhrOVXyE6BHwceTa9La7ieXwIeSfU8Cfyn1F6T9XSr7SLe/cC8Jush+4zgsfRaf+jvdq3Wk97vbKA1/Zn7W+D4suvx40nMzKywo/mylZmZ9ZPDw8zMCnN4mJlZYQ4PMzMrzOFhZmaFOTzMzKwwh4eZmRX2/wFqr7lV3S063AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pylab.hist(df, bins=100, density=True)\n",
    "pylab.plot(x, fitted_pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But, how do we know if our assumptions are 'good' or not? This is where goodness of fit testing would come in."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluating the goodness of fit\n",
    "\n",
    "While all these tests (excluding chi square) comes out of the box with the Fitter library, it is good to get a conceptual understanding of the logic and how to implement it manually.\n",
    "\n",
    "Note that all these metrics are relative measures and provide information when comparing between different fits/models. Absolute standalone values from these metrics rarely provide enough information to make informed decisions regarding model selection.\n",
    "\n",
    "One common pattern to note is that in many of the tests and formulae, a $ log(\\tilde{L}) $ is usually seen in the calculations. This is because when multiplying small numbers repeatedly, the results tend to get rounded to 0. Thus, a log transformation is applied to allow for additions during algebraic manipulation without compromising the extrema of the distribution functions. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sum of squared errors\n",
    "\n",
    "Sum of square errors is a measure of quadratic loss, and is a commonly used metric which is also referred to as the L2 loss.\n",
    "\n",
    "$ SSE = \\Sigma (x_{pred} - x_{actual})^2$\n",
    "\n",
    "\n",
    "Note that SSE metric would be impacted materially with outlier values, which would suggest that more attention be put into cleaning up the data before fitting when using this method. Obviously when using SSE, a lower results suggests a better fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sum of squared errors: 8.627423773657961e-08\n"
     ]
    }
   ],
   "source": [
    "sq_error = pylab.sum((fitted_pdf - y) ** 2)\n",
    "\n",
    "print(f'Sum of squared errors: {sq_error}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Information criterion\n",
    "\n",
    "The `Akaike` and `Bayesian` information criterion metrics (AIC and BIC) aim to help users select from a set of models, and uses the likelihood function as a starting point and adding penalizing terms based on the number of parameters of the model and observations.\n",
    "\n",
    "$ AIC = -2ln(\\tilde{L}) + 2k $\n",
    "\n",
    "$ BIC = -2ln(\\tilde{L}) + kln(N) $\n",
    "\n",
    "Where k is the model degrees of freedom, and N is the number of observations. In both these measures, a lower result represents a better fit. While these measures look similar, the penalization terms results in different motivations behind the selections.\n",
    "\n",
    "Generally, AIC presents the danger that a model overfits the data as more degrees of freedom in the observed model would result in a lousier score. BIC on the other hand presents the danger of underfitting the data given the errors from the likelihood are added to the size of the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AIC: 3575.446841073396\n",
      "BIC: -1990.5592711125862\n"
     ]
    }
   ],
   "source": [
    "logLik = np.sum(st.expon.logpdf(x, *params))\n",
    "k = len(params)\n",
    "n = len(df)\n",
    "aic = 2 * k - 2 * logLik\n",
    "bic = n * np.log(sq_error / n) + k * np.log(n)\n",
    "\n",
    "print(f'AIC: {aic}')\n",
    "print(f'BIC: {bic}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kullback Leiblier divergence\n",
    "\n",
    "The KL divergence is a measure that aims to quantify `statistical distance` between 2 probability distributions. Given a random variable X that represents our observations, a starting point to measure difference between 2 distributions p and q would be:\n",
    "\n",
    "$ \\Delta_{x_i}(p, q) = log(p_{\\theta}(x_i)) - log(q_{\\phi}(x_i)) $\n",
    "\n",
    "$ \\Delta_{x_i}(p, q) =  log\\bigg(\\dfrac{p_{\\theta}(x_i)}{q_{\\phi}(x_i)}\\bigg) = log(likelihood ratio) $\n",
    "\n",
    "Now that we have a measure of difference between 2 distributions for a single data point, we need to extend this concept to get the expectation (or average difference) across all data points.\n",
    "\n",
    "$ \\Delta(p, q) =  \\int p_{\\theta}(x_i) log\\bigg(\\dfrac{p_{\\theta}(x_i)}{q_{\\phi}(x_i)}\\bigg) dx $\n",
    "\n",
    "Note that depending on whether p or q is in the denominator, the KL divergence measure would result in a different answer, making it asymmetric. This leads to controversy around actually calling this a distance measure as opposed to a divergence instead.\n",
    "\n",
    "The choice of reference and approximate distributions lead to different results and application of reverse KL (q in the numerator) can be observed in ML models when evaluating cross entropy.\n",
    "\n",
    "For more information please refer to the [Scipy kl_div documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.kl_div.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kl divergence is: 1.743687567240447\n"
     ]
    }
   ],
   "source": [
    "# Adjust 0 values to a small number to prevent div 0 errors when calculating KL divergence\n",
    "\n",
    "y[y == 0] = 0.0001\n",
    "\n",
    "kullback_leibler = kl_div(fitted_pdf, y)\n",
    "\n",
    "print(f'kl divergence is: {kullback_leibler}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, note that a standalone number as above doesn't provide too much context. This must be used in conjunction with other KL divergence figures on different distributions for comparison."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kolmogorov-Smirnov test\n",
    "\n",
    "The KS test is used when wanting to know how well the distribution of the observed data conforms to some theoretical distribution, and uses the cumulative  distribution functions (CDF) as a basis of comparison. The hypotheses of the KS test can be summarized as follows:\n",
    "\n",
    "$ H_0 $ : Two samples drawn from populations with the same distribution\n",
    "\n",
    "$ H_1$ : Two samples drawn from populations with different distributions\n",
    "\n",
    "The test statistic is commonly known as D and can be defined as follows:\n",
    "\n",
    "$ D = sup|F_{S}(X) - F_{T}(X)| $\n",
    "\n",
    "The T represents the theoretical distribution and S is the observed sample. This means that the test statistic represents the greatest distance between the 2 CDFs, and is assumed to follow the Kolmogorov distribution. Depending on the value of D, we can get p-values to assess the probability of observing the max distance between our data's CDF and the theoretical one given the null hypothesis is true.\n",
    "\n",
    "In the example below, we see that D=0.18 and the probability of observing this given the 2 distributions are the same is ~0.2%. At a 95% confidence level, we would reject the null hypothesis that the 2 distributions are the same since 0.2% < 5%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ks statistic D: 0.18499042546927535\n",
      "probability of observing this value given H0 is true: 0.0023880598948170517\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1ef34297340>"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD4CAYAAADrRI2NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkyUlEQVR4nO3de3xU9Z3/8ddncuUSbiFcw1XRioC3CCpKsVYFtFqtPn6Kumq3dd3CT9vd3yru+qO6ttbrr64/tSxa1631wnpD63oparG1XrgZ5A7hFkIgJBDInSQz3/1jhjAkk2SAZCZn5v18PMYzc86Zmc83wptvvuec7zHnHCIikjx88S5ARERiS8EvIpJkFPwiIklGwS8ikmQU/CIiSSY1Xl/cv39/N3LkyHh9vYiIJy1fvrzMOZdzPJ8Rt+AfOXIky5Yti9fXi4h4kpltP97P0FCPiEiSUfCLiCQZBb+ISJJR8IuIJBkFv4hIkmk3+M3seTPbY2arW9luZvakmRWY2TdmdmbHlykiIh0lmh7/C8C0NrZPB8aEHrcBvzn+skREpLO0ex6/c+7PZjayjV2uBH7ngvM7f2lmfcxssHNuV0cVKV1AIAD+g9BYB4314A9/NECgAfyNEGgMPg80Bt8TaATnh4A/tAyAC3/4m712wQfu8GtcK0vaeR7egLAXbU5F3sq2qKYv7yJTnGuq9S4r4BzLtpcz4NSpjJz0vbjV0REXcA0FdoS9LgqtaxH8ZnYbwd8KGD58eAd8tUTF3wDVpVBdBjVlULMPasuhdj/U7Ye6A3CwMviorw49qqChFhprg0t/fbxbIUfF4l2ANHPon+M8B1811Hs++CP9CYvY5XDOzQfmA+Tl5alb0lEa66F8G+wtgH1b4MAO2F8IB4qgcncw9Fvrjab1gMxekJEVfKT3hO79IK07pHcPLlMzIa0bpGZASgakpgeXKemQkhZ8+NIgJTW49KWGPXxgKeBLCVv6Qg87vA5rtj40CmkW2tbKEthUWsWBWn9o/6b/gBnuiD+eYc+trWBsZVub7xFp3a4DdfzrH9ZQdbCRf71iHNfm5ca1no4I/iJgWNjrXKC4Az5XIqkth6LlsCsfSlbD7tXBsHf+w/uk94Q+w6F3Lgw5A7IGQ88B0CMHevSH7tnQrR9k9g6GuEfVNfj513fX8vJXhfEuRaRdJw3sycs/PoeTBmbFu5QOCf53gNlm9iowCTig8f0OVLMPtn4Km/8EhV9C2YbD2/oMh4HjYeyV0H8MZJ8I/UZDt74J3zst2FPJ7Je/Zv3uSv5uymguGHNcc1aJdCqfD84c3pfMtJR4lwJEEfxm9gowFehvZkXAz4E0AOfcPOA9YAZQANQAt3ZWsUmjfDusewfWvg1FywAHGb1g+Dkw/loYdnawJ5/Zu9NL+WhtCXPfXs2+mq41xl/fGKBv93ReuPVspp48IN7liHhKNGf1XN/OdgfM6rCKklV9Nax+A1b8DoqWBtcNPg2+fTec8B0YelZwDD1W5TQGeOj99Tz/162cMrgX3zttSMy+OxoZqT5uOGcEA3tlxrsUEc+J27TMErK/ED7//7DyVThYATnfgu/eFxy+6Te6U7/668JynvtsK/WNgRbbtpVVs2lPFbecN5J7ZnyLjNSu8SuqiBw/BX+87NsCf3k8GPgYjLsazro1OJzTyePzgYDjuc+28MgHG+jVLS1ir7l7Rir/ftNZXHrqoE6tRURiT8EfI3UNft79ZhcNNfuZUDCPUwpfwZmPjbnXsmbUrdRkDgyeC1V83PdYaNcn6/eweEMp004dxMM/mEDv7mmd/p0i0nUo+GNgS2kVs19awZg97/MvaS/TnwO86p/KrxuvoXRjX9hYBpTFrJ70VB8PXHkqN54zAkvws39EpCUFfydb+PVOHn7rC37hm89F6V/ROOgMKi56iEsHn8GlcaqpW3oK3dP1v14kWelvfyepqW9k7ttr2L7iI/7Q7Tdksx8ueoDUc2fTx6fZsEUkfhT8nWDdrgpmv7yCyfveYkHG77A+I7Ef/BcM1YzVIhJ/Cv4O8PG6En753+uo9wdPiyytqOWfM/6Lm9MWwskz4Or5wXlwRES6AAX/cdpTWcc/vraSft3TmTiqHymugZtKfsOEfX+EvL+FGY+GJiETEekaFPzHwTnHP7+5itp6P8/+fR4nZHeD12+FfX+Ei+bC+f+Q8HPmiIj36CjjcXhjxU4+WreHf7r0ZE7o3wPe/Wlwfp1LfgkX/KNCX0S6JAX/MSreX8v9f1jDxJH9+OHkUfDRfcF5di74RzhvdrzLExFplYL/GHy1ZS9XP/M5/oDj0Wsn4Fv2HPz1ieCUC9/5v/EuT0SkTQr+o+APOP7to01c/+yXZKb5+K+/O5cRtevhg3tgzKVw2eMa3hGRLk8Hd6O0p6KOO1/N54ste/n+6UP4xVXj6emvgH+/JXiHq6vm6ewdEfEEBX8UPt1Yyj8syKem3s8j10zg2rNyMefgjb+Hyl3www+D96kVEfEABX8bGvwBHv/jRuZ9upmTB2bx1MwzGHPofplf/gY2fgDTH4Xcs+JbqIjIUVDwt2LHvhruePVrvi7cz/UTh/Pz7409fL/M/YXwyQNw0jSY+OP4FioicpQU/BH8af0e7nz1a5yDp2aeweUTwm476Bz89/8BDGY8poO5IuI5Cv5mduyrYfbLKxie3YN5N57JiOweR+6w7h3Y9GHwIq0+w+JTpIjIcVDwhwkEHP/0+krMjPk3ncWwft2P3KHuALx3FwwaD5Nuj0+RIiLHScEf5ndfbOPLLft46OrxLUMfYPHDUFUC170MKfrRiYg36QKukK1l1Tz0wXqmnpzD/zo7whDOgZ2w9Dk4fabO4hERT1PwEwz9219cTnqKj4eunhD5PrR/eRycH759V+wLFBHpQEk/XvF2/k7++c1VpKb4eGrmmQzqndlyp/LtwQnYzrgJ+o6MeY0iIh0paYO/tt7Pfe+sYcGyHeSN6Mu/XX8GQ/t0i7zznx8B88GUf4ptkSIinSApg3/D7kpmv7yCgtIqZl14Aj/77kmkprQy6rV3M+S/ErxQq/fQ2BYqItIJki7431xRxD1vriIrM5Xf/XAiF4zJafsNf3kcUtKDd9MSEUkASRX89Y0B5ryxivG5vfnNjWcyICvCeH646jJY9Rqc+TeQNTA2RYqIdLKkOqtna1k19f4AN50zov3QB1jxn+Cvh7M1H4+IJI6kCv6NJZUAnHRohs22+Bth2X/AqCkw4FudXJmISOwkVfBvKqnEZzA6p0f7O2/8AA7sgIm3dX5hIiIxlFTBv6GkkpHZPQ5Pr9yWJfOhVy6cNL3zCxMRiaGkCv5NJVXRDfOUboCtn0LerZqTR0QSTtIEf12Dn217qzlpYM/2d176XPAUzjNv7vzCRERiLKrgN7NpZrbBzArMbE6E7b3N7A9mttLM1pjZrR1f6vHZXFpFwHH41omt8TcET+E85XvQs51z/EVEPKjd4DezFOBpYDowFrjezMY2220WsNY5dxowFXjczNI7uNbjsqmkCoCTB7UT/Jv/BLXlMP7aGFQlIhJ70fT4JwIFzrktzrl64FXgymb7OCDLgtNa9gT2AY0dWulx2lhSSarPGNn8jlrNrX4DMnvDCd+JTWEiIjEWTfAPBXaEvS4KrQv3FHAKUAysAu50zgWaf5CZ3WZmy8xsWWlp6TGWfGw2llQyqn8P0lPbaHJDLax/NzjMk5oRu+JERGIomuCPdDdx1+z1pUA+MAQ4HXjKzHq1eJNz851zec65vJyc2I6fbyyp4qT2hnk2/RHqq2DcNbEpSkQkDqIJ/iIg/JZUuQR79uFuBd50QQXAVqDLXO5aU9/IjvIaThrQTvCveh165MDIC2JTmIhIHEQT/EuBMWY2KnTA9jrgnWb7FAIXAZjZQOBkYEtHFno8CvZU4Rxtn8pZVxHs8Z96lc7dF5GE1m7COecazWw28CGQAjzvnFtjZreHts8DHgBeMLNVBIeG7nbOlXVi3UdlY+iMnjaHeja8B411MO4HMapKRCQ+ouraOufeA95rtm5e2PNi4JKOLa3jbCqpJD3Fx4h+3VvfafWb0HsY5E6MXWEiInGQFFfubiypZHROj9bvslVfA1sWw7cuB19S/EhEJIklRcptLKlq+8KtbZ+B/yCM+W7sihIRiZOED/7KugZ27q9te3K2go8gtRuMOD92hYmIxEnCB//a4goAxg5pcVnBYQWLYNQFkBbFXblERDwu4YN/dSj4T20t+Pduhn1b4EQN84hIckj44F9TfICcrIzW77Fb8HFwqeAXkSSR+MG/s4Jx7Q3z9BsN2SfErigRkThK6OCva/BTUFrFqUN6R96hoQ62/gVOvDi2hYmIxFFCB//63ZX4A45xQ1vp8W//KzTWaphHRJJKQgf/muIDAK33+As+gpQMGKnTOEUkeSR08K/eWUGvzFRy+3aLvEPBxzByMqS3MZWDiEiCSejgX1t8gFOH9CZ4Y7BmqsugbAOMmhL7wkRE4ihhg7/BH2Dd7srWx/cLvwwuh58bu6JERLqAhA3+zaVV1DcGWh/f3/FlcHx/yBmxLUxEJM4SNvjX7Axesdtmj3/IGbq3rogknYQN/tXFB+iWlsKo/hHuutVQC8X5MPycmNclIhJvCRv8a4orOGVwFim+CAd2d66AQIPG90UkKSVk8AcCjrXFFa2P7xd+EVwO0922RCT5JGTw7yivoepgY+tTMRd+CTnfgu79YluYiEgXkJDBf+jm6hHvuhUIwI4lMGxSjKsSEekaEjT4KwEYMyDCgd3SdXDwgMb3RSRpJWzwD+3TjazMtJYbD43v64weEUlSCRr8VYwZGKG3D1D4FfQcCH1HxrQmEZGuIuGCv9EfYPOeqtZvrl74ZbC3H2n+HhGRJJBwwb99Xw31/kDk4K/eCwcKYWhe7AsTEekiEi74N4UO7J4Uaahn98rgcvCEGFYkItK1JFzwHzqV88RIZ/Ts+ia4HKTgF5HklXDBv6GkkmH9utE9PbXlxt3fQO9hunBLRJJawgX/ppJKTm7twO6ub9TbF5Gkl1DBX98YYEtpNWMiBf/BKthboPF9EUl6CRX82/ZW0xhwkQ/slqwBnHr8IpL0Eir4Nzad0ROhx787dGBXPX4RSXIJFvxV+AxOyIl0Rs9K6NYPeg2NfWEiIl1IQgX/ppJKRmT3IDMtpeXG3d8Ee/u6YldEklxUwW9m08xsg5kVmNmcVvaZamb5ZrbGzD7t2DKjs6GkMvL4vr8B9qzT+L6ICFEEv5mlAE8D04GxwPVmNrbZPn2AZ4ArnHOnAtd2fKltO9joZ/vemsjj+6XrwV8Pg0+LdVkiIl1OND3+iUCBc26Lc64eeBW4stk+M4E3nXOFAM65PR1bZvu2ldXgDzhdsSsi0o5ogn8osCPsdVFoXbiTgL5mttjMlpvZ33RUgdHasa8GgBHZPVpu3P0NpHWH7BNiXJWISNcTYV6DFiIdDXURPucs4CKgG/CFmX3pnNt4xAeZ3QbcBjB8+PCjr7YNReXB4M/t263lxl3fwMBx4Itw0FdEJMlE0+MvAoaFvc4FiiPs84Fzrto5Vwb8GWgxoO6cm++cy3PO5eXk5BxrzRHt3F9LZpqP7B7pR24IBGD3Kp2/LyISEk3wLwXGmNkoM0sHrgPeabbP28AFZpZqZt2BScC6ji21bUXltQzt0w1rfrpmxU6or4SBp8ayHBGRLqvdoR7nXKOZzQY+BFKA551za8zs9tD2ec65dWb2AfANEACec86t7szCm9u5v5ahfbu33FAWGm3qf3IsyxER6bKiGePHOfce8F6zdfOavX4UeLTjSjs6ReW1jBvau+WGsk3BZf+TYluQiEgXlRBX7tbUN7Kvup6hfSIc2C3bAJl9oEf/mNclItIVJUTw7yyvBVo5o6dsU7C3r6kaRESABAn+ov1tBf9GDfOIiIRJjOBv6vE3O7hbux+qSiBHwS8ickhCBP/O8lrSU3zk9Mw4csPeguBSPX4RkSYJEfxF5TUM7pOJz9dsHL90Q3Cp4BcRaZIQwb9zf23r4/u+NOgzIvZFiYh0UQkR/Ieu2m2hbFNwYraUqC5XEBFJCp4P/roGP6WVB1se2IXQGT1jYl+UiEgX5vngLw6dytmix+9vgPKtmqpBRKQZzwf/ztbO4d+3FQKNOrArItKM54P/0Dn8Q5sHf9mhM3o01CMiEs7zwb+zvJYUnzGoV+aRG5pm5VTwi4iE83zwF5XXMKhXJqkpzZpStgmyhkBGhJuvi4gkMc8Hf5vn8GuqBhGRFjwf/EXltS3H9507PCuniIgcwdPBX98YoKSiruU5/FUlcLACsjW+LyLSnKeDf/eBOgIOcpufw1++LbjsNyrmNYmIdHWeDv5D5/C3GOop3x5cao4eEZEWPB38eyrrABiQ1Ww65v2Hgn94jCsSEen6PB38ZVX1APRvPg9/+XboOQjSMiO8S0QkuXk8+A+S6jN6d0s7csP+7dBXwzwiIpF4O/grD5LdM73lDVjKt2t8X0SkFd4O/qqDLYd5/A1QUaQev4hIKzwe/PXkND+we6AIXEA9fhGRVng8+CP0+A+d0aMev4hIRJ4Nfucce6vqI5/RA+rxi4i0wrPBX1HbSL0/QP+e6Udu2L8dLAV6DY1PYSIiXZxng7+06iBAyzH+8u3QO1c3WBcRaYVng78sFPwRx/g1vi8i0qrEC36dwy8i0ibPBn9p5aHgDxvjr6+B6j3q8YuItMGzwV9WdZAUn9G3e1jw7y8MLvuMjEtNIiJe4N3gr6ynX49m0zXoHH4RkXZ5N/gjXbx16AYsGuMXEWlVVMFvZtPMbIOZFZjZnDb2O9vM/GZ2TceVGFkw+Judw1++HVK7Qc8Bnf31IiKe1W7wm1kK8DQwHRgLXG9mY1vZ72Hgw44uMpKI8/Ts3x68+YpZ5DeJiEhUPf6JQIFzbotzrh54Fbgywn7/G3gD2NOB9UXknKO06iA5kU7l1Pi+iEibogn+ocCOsNdFoXVNzGwocBUwr60PMrPbzGyZmS0rLS092lqbVB5spL4xcOQYv3OhHr+CX0SkLdEEf6RxE9fs9RPA3c45f1sf5Jyb75zLc87l5eTkRFliS2WHzuHPChvjry2HgxXq8YuItCOaCW2KgGFhr3OB4mb75AGvWnBsvT8ww8wanXMLO6LI5iLea/dA6JcS3WBdRKRN0QT/UmCMmY0CdgLXATPDd3DOjTr03MxeAN7trNCHVqZrqAj9W6RZOUVE2tRu8DvnGs1sNsGzdVKA551za8zs9tD2Nsf1O0ObwZ81ONbliIh4SlRzFzvn3gPea7YuYuA75245/rLaVlZ5EJ9Bvx5hY/yVu8B80HNgZ3+9iIinefLK3dKq4HQNKeHTNVTsCoa+5uEXEWmTN4O/MsJ0DZXFGuYREYmCJ4M/4jw9Fbug15D4FCQi4iGeDf4W0zWoxy8iEhXPBb9zruUEbfXVUHcAein4RUTa47ngr673U9fQbLqGil3BZZaGekRE2uO54G+ariE8+CsPXbylHr+ISHu8F/yHLt7KUo9fRORYeDf4w8f41eMXEYma54L/pIFZ3HvZKeT27X54ZcUuyOgFGVnxK0xExCM8d5nr6JyejM7peeRKncopIhI1z/X4I6rYpWEeEZEoJUbwV+7SgV0RkSh5P/gDfqjcrR6/iEiUvB/81aXg/BrjFxGJkveDv+nOWxrqERGJhveDv/LQxVvq8YuIRMP7wa8ev4jIUUmM4PelQo+ceFciIuIJ3g/+yl3QcxD4UuJdiYiIJ3g/+CuKdSqniMhR8H7wV+7SgV0RkaPg/eDXvXZFRI6Kt4P/YCXUV6rHLyJyFLwd/IduwKIev4hI1Lwd/E0Xbw2Kbx0iIh7i7eCvLg0uew6Mbx0iIh6SGMGvi7dERKLm/eD3pUJmn3hXIiLiGd4O/qo90L0/+LzdDBGRWPJ2YlaXaZhHROQoeTz4S6FH/3hXISLiKd4P/p4D4l2FiIinpMa7gONSXaqhHpFO1tDQQFFREXV1dfEuJalkZmaSm5tLWlpah3+2d4O/vhoaajTUI9LJioqKyMrKYuTIkZhZvMtJCs459u7dS1FREaNGjerwz49qqMfMppnZBjMrMLM5EbbfYGbfhB6fm9lpHV5pczqHXyQm6urqyM7OVujHkJmRnZ3dab9ltRv8ZpYCPA1MB8YC15vZ2Ga7bQW+7ZybADwAzO/oQluoLgsue2iMX6SzKfRjrzN/5tH0+CcCBc65Lc65euBV4MrwHZxznzvnykMvvwRyO7bMCKr2BJca6hEROSrRBP9QYEfY66LQutb8LfB+pA1mdpuZLTOzZaWlpdFXGYmGekTkOJ133nnHtK0t9913H4899tixlhQT0QR/pN83XMQdzS4kGPx3R9runJvvnMtzzuXl5BxnYCv4ReQ4ff755y3W+f3+VrclimjO6ikChoW9zgWKm+9kZhOA54Dpzrm9HVNeG6rLIKMXpGV2+leJSND9f1jD2uKKDv3MsUN68fPvndrufr///e958sknqa+vZ9KkSTzzzDP07t2bWbNm8dFHH9G3b18efPBB7rrrLgoLC3niiSe44ooreOGFF3jrrbc4ePAgW7duZebMmfz85z8HoGfPnlRVVbF48WLuv/9+Bg8eTH5+PmvXrm3aBvDII4/w4osv4vP5mD59Og899BDPPvss8+fPp76+nhNPPJEXX3yR7t27d+jPprNE0+NfCowxs1Fmlg5cB7wTvoOZDQfeBG5yzm3s+DIjqN6j8X2RJLFu3ToWLFjAX//6V/Lz80lJSeGll16iurqaqVOnsnz5crKysrj33ntZtGgRb731FnPnzm16/5IlS3jppZfIz8/ntddeY9myZS2+Y8mSJfzyl79k7dq1R6x///33WbhwIV999RUrV67krrvuAuDqq69m6dKlrFy5klNOOYXf/va3nftD6EDt9vidc41mNhv4EEgBnnfOrTGz20Pb5wFzgWzgmdCR6EbnXF7nlY0u3hKJg2h65p3h448/Zvny5Zx99tkA1NbWMmDAANLT05k2bRoA48ePJyMjg7S0NMaPH8+2bdua3n/xxReTnZ0NBAP7s88+Iy/vyIiaOHFixHPmP/roI2699dam3ny/fv0AWL16Nffeey/79++nqqqKSy+9tMPb3VmiuoDLOfce8F6zdfPCnv8I+FHHltaO6jLoNzqmXyki8eGc4+abb+ZXv/rVEesfe+yxptMefT4fGRkZTc8bGxub9mt+amSkUyV79OjR6ndH2v+WW25h4cKFnHbaabzwwgssXrz4qNoUT96dq0cTtIkkjYsuuojXX3+dPXuCp3Hv27eP7du3R/3+RYsWsW/fPmpra1m4cCGTJ0+O+r2XXHIJzz//PDU1NU3fDVBZWcngwYNpaGjgpZdeOorWxJ83p2wI+KFmry7eEkkSY8eO5Re/+AWXXHIJgUCAtLQ0nn766ajff/7553PTTTdRUFDAzJkzWwzztGXatGnk5+eTl5dHeno6M2bM4MEHH+SBBx5g0qRJjBgxgvHjx1NZWXksTYsLcy7imZmdLi8vz0U6wBKVqlJ47ESY/ihMuq1jCxORI6xbt45TTjkl3mUcsxdeeIFly5bx1FNPxbuUoxbpZ29my4/3GKo3h3qazuHXUI+IyNHy5lCPLt4SkSjdcsst3HLLLfEuo0vxdo9fN2ERETlq3g5+9fhFRI6ad4PfUiCzT7wrERHxHO8Gf4/+4PNm+SIi8eTN5Kwu0zn8ItLl7N+/n2eeeabpdXFxMddcc02b79m2bRvjxo3r7NKO4M3gr9IEbSLS9TQP/iFDhvD666/HsaLIvHs6p+bpEYm99+fA7lUd+5mDxsP0h9rdrfm0zD/84Q/58Y9/zJIlS/D7/UycOJEFCxZQVlbG3Llzyc7OZsOGDUyZMoVnnnkGn8/HK6+8woMPPohzjssuu4yHH34YCE7PfOedd/Luu+/SrVs33n77bQYOHEhpaSm33347hYWFADzxxBNMnjyZ++67j8LCQrZs2UJhYSE//elPueOOO5gzZw6bN2/m9NNP5+KLL2bWrFlcfvnlrF69mm3btnHTTTdRXV0NwFNPPXXMN3s5Xt7s8VeX6YwekSQSaVrmDRs2cMUVV3Dvvfdy1113ceONNzYNmSxZsoTHH3+cVatWsXnzZt58802Ki4u5++67+eSTT8jPz2fp0qUsXLgQgOrqas455xxWrlzJlClTePbZZwG48847+dnPfsbSpUt54403+NGPDs9FuX79ej788EOWLFnC/fffT0NDAw899BAnnHAC+fn5PProo0e0YcCAASxatIgVK1awYMEC7rjjjtj88CLwXo+/vhoaqjXUIxIPUfTMO0Nr0zLPnTuXs88+m8zMTJ588smm/SdOnMjo0cFRgeuvv57PPvuMtLQ0pk6dyqG7/91www38+c9/5vvf/z7p6elcfvnlAJx11lksWrQICE7JHD4/f0VFRdOcPJdddhkZGRlkZGQwYMAASkpK2mxDQ0MDs2fPbvqHa+PG2Ny6JBLvBb8u3hJJOq1Ny7x7926qqqpoaGigrq6uaWrlSNMwtzUvWVpaWtN7UlJSmqZ0DgQCfPHFF3Tr1q3Few5NAd38Pa359a9/zcCBA1m5ciWBQIDMzPjdPdB7Qz3VZcGlhnpEkkZr0zLfdtttPPDAA9xwww3cfffhW30vWbKErVu3EggEWLBgAeeffz6TJk3i008/paysDL/fzyuvvMK3v/3tNr/3kksuOWJyt/z8/Db3z8rKanWWzgMHDjB48GB8Ph8vvvhi071948G7PX4N9YgkjUjTMl955ZWkpqYyc+ZM/H4/5513Hp988gk+n49zzz2XOXPmsGrVKqZMmcJVV12Fz+fjV7/6FRdeeCHOOWbMmMGVV17Z5vc++eSTzJo1iwkTJtDY2MiUKVOYN29eq/tnZ2czefJkxo0bx/Tp05k1a1bTtp/85Cf84Ac/4LXXXuPCCy9s9cYvseC9aZkLv4QvnoYZj0LWoI4vTESO4LVpmRcvXsxjjz3Gu+++G+9SjltnTcvsvR7/8HOCDxEROSbeC34RkTZMnTqVqVOnxruMLs17B3dFJObiNSSczDrzZ67gF5E2ZWZmsnfvXoV/DDnn2Lt3b6ed8qmhHhFpU25uLkVFRZSWlsa7lKSSmZlJbm5up3y2gl9E2pSWlsaoUaPiXYZ0IA31iIgkGQW/iEiSUfCLiCSZuF25a2alwPZjeGt/oKyDy4mnRGsPJF6bEq09oDZ5QWvtGeGcO67JyuIW/MfKzJYd7+XKXUmitQcSr02J1h5Qm7ygM9ujoR4RkSSj4BcRSTJeDP758S6ggyVaeyDx2pRo7QG1yQs6rT2eG+MXEZHj48Uev4iIHAcFv4hIkvFM8JvZNDPbYGYFZjYn3vWEM7PnzWyPma0OW9fPzBaZ2abQsm/YtntC7dhgZpeGrT/LzFaFtj1pobs/m1mGmS0Irf/KzEbGoE3DzOxPZrbOzNaY2Z1ebpeZZZrZEjNbGWrP/V5uT1gtKWb2tZm9myDt2RaqJd/MliVIm/qY2etmtj709+ncuLfJOdflH0AKsBkYDaQDK4Gx8a4rrL4pwJnA6rB1jwBzQs/nAA+Hno8N1Z8BjAq1KyW0bQlwLmDA+8D00PqfAPNCz68DFsSgTYOBM0PPs4CNodo92a7Qd/cMPU8DvgLO8Wp7wtr1D8DLwLsJ8uduG9C/2Tqvt+k/gR+FnqcDfeLdpk5tcAf+4M4FPgx7fQ9wT7zralbjSI4M/g3A4NDzwcCGSLUDH4baNxhYH7b+euDfw/cJPU8leDWfxbh9bwMXJ0K7gO7ACmCSl9sD5AIfA9/hcPB7tj2h79lGy+D3bJuAXsDW5t8R7zZ5ZahnKLAj7HVRaF1XNtA5twsgtBwQWt9aW4aGnjdff8R7nHONwAEgu9Mqbyb0q+MZBHvJnm1XaFgkH9gDLHLOebo9wBPAXUAgbJ2X2wPggD+a2XIzuy20zsttGg2UAv8RGpJ7zsx6EOc2eSX4LcI6r56H2lpb2mpj3NpvZj2BN4CfOucq2to1wrou1S7nnN85dzrBnvJEMxvXxu5duj1mdjmwxzm3PNq3RFjXZdoTZrJz7kxgOjDLzKa0sa8X2pRKcBj4N865M4BqgkM7rYlJm7wS/EXAsLDXuUBxnGqJVomZDQYILfeE1rfWlqLQ8+brj3iPmaUCvYF9nVZ5iJmlEQz9l5xzb4ZWe75dzrn9wGJgGt5tz2TgCjPbBrwKfMfMfo932wOAc644tNwDvAVMxNttKgKKQr9dArxO8B+CuLbJK8G/FBhjZqPMLJ3gAYx34lxTe94Bbg49v5ngGPmh9deFjsSPAsYAS0K/7lWa2Tmho/V/0+w9hz7rGuATFxrQ6yyhGn4LrHPO/b+wTZ5sl5nlmFmf0PNuwHeB9V5tj3PuHudcrnNuJMG/D5845270ansAzKyHmWUdeg5cAqz2cpucc7uBHWZ2cmjVRcDauLepMw/UdPBBkhkEzyzZDPxLvOtpVtsrwC6ggeC/vn9LcIztY2BTaNkvbP9/CbVjA6Ej86H1eQT/oG8GnuLwldWZwGtAAcEj+6Nj0KbzCf66+A2QH3rM8Gq7gAnA16H2rAbmhtZ7sj3N2jaVwwd3PdseguPhK0OPNYf+nnu5TaHvPB1YFvqztxDoG+82acoGEZEk45WhHhER6SAKfhGRJKPgFxFJMgp+EZEko+AXEUkyCn4RkSSj4BcRSTL/AxguuxgfloGgAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# weird issue occurring without converting explicitly to int from the pandas dataframe\n",
    "# to check out scipy ktest source code https://github.com/scipy/scipy/blob/v1.9.1/scipy/stats/_stats_py.py#L7889-L8039\n",
    "\n",
    "_data = np.array([x for x in df['claims']]).astype(int)\n",
    "\n",
    "dist_fitted = st.expon(*params)\n",
    "ks_stat, ks_pval = kstest(_data, dist_fitted.cdf)\n",
    "\n",
    "print(f'ks statistic D: {ks_stat}')\n",
    "print(f'probability of observing this value given H0 is true: {ks_pval}')\n",
    "\n",
    "# Plotting the ks test to visualize\n",
    "\n",
    "data_pdf = y / sum(y)\n",
    "data_cdf = np.cumsum(data_pdf)\n",
    "\n",
    "plt.plot(x, data_cdf, label=\"empirical\")\n",
    "plt.plot(x, dist_fitted.cdf(x), label='exponential')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Chi square test\n",
    "\n",
    "The chisquare test is used to determine whether or not there is a statistically significant difference between the expected frequencies and the observed. The test statistic is defined as:\n",
    "\n",
    "$ \\chi^2 = \\Sigma \\dfrac{(O_i - E_i)^2}{E_i} $\n",
    "\n",
    "Which is approximately distributed $ \\chi^2(df) $ where df is the degrees of freedom in the model, or number of bins - 1 - number of parameters in model. \n",
    "\n",
    "- $k$ represents the total number of bins\n",
    "- $O_i$ represents the number of observed data points in the bin i\n",
    "- $E_i$ represents the number of expected data points in bin i\n",
    "\n",
    "We note that once again, the choice of bins here would directly affect the results, and would suggest referring back to the general guidelines in notebook 1 for bin selection. Once a set of bins are selected, calculating $O_i$ is quite straight forward given our data. Calculating $E_i$ on the other hand would depend on the assumed distribution. Following on from our exponential distribution assumption, given that `a` and `b` are some numerical end points of a selected bin:\n",
    "\n",
    "$ Pr_{\\theta}(a < X \\leq b) = \\int_{a}^{b} \\lambda e^{-\\lambda} dx = e^{-a \\lambda} - e^{-b \\lambda} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chisquare test statistic: 96.00216105934356\n",
      "chisquare p-value: 0.5095574753231276\n"
     ]
    }
   ],
   "source": [
    "f_obs, bin_edges = np.histogram(df, bins=100, density=True)\n",
    "cdf = getattr(st, 'expon').cdf(bin_edges, loc=params[0], scale=params[1])\n",
    "expected_values = len(df) * np.diff(cdf)\n",
    "results = st.chisquare(f_obs, expected_values, ddof=len(params))\n",
    "\n",
    "\n",
    "print(f'chisquare test statistic: {results[0]}')\n",
    "print(f'chisquare p-value: {results[1]}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given our selection of bins, the probability of observing a chisquare statistic of 96 is ~50%. Thus we don't have enough evidence to reject the exponential distribution."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.5 ('base')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "7f3e2235ce0c3d3b0e4a04e7f9c82d922b04c8831f13f1d259fbff6721a8fb0e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}