kolmogorov_smirnov.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "2d208ae2-5439-4173-8c25-2565646348e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as st\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43ef3304-5f38-433e-a5b6-a90e29bb5249",
   "metadata": {},
   "source": [
    "* https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test\n",
    "* https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kstwo.html\n",
    "* https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kstest.html\n",
    "* https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ks_2samp.html#id1\n",
    "* https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kstwobign.html\n",
    "* https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.kolmogorov.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38a98dcd-3b83-4bcf-ad75-4a05e2d9740d",
   "metadata": {},
   "source": [
    "## Kolmogorov-Smirnov Test"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33ab90b8-f596-451d-bac0-049d46d8dfd4",
   "metadata": {},
   "source": [
    "Empirical CDF $F_n(x)$ computed using $n$ sample observation $X_1, X_2,\\dots, X_n$ is defined as \n",
    "\n",
    "$$\n",
    "F_n(x) = \\frac{1}{n}\\sum_{i=1}^{n}\\mathbb{1}_{X_i \\leq x}\n",
    "$$\n",
    "\n",
    "where $\\mathbb{1}_{X_i \\leq x}$ is an indicator function that is equal to 1 for $X_i \\leq x$ and 0 otherwise."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afbd6b86-2f2f-4638-a979-2e2f1805fef5",
   "metadata": {},
   "source": [
    "## One-Sample KS Test\n",
    "\n",
    "One sample KS test can be used to compare a sample with a reference distribution. Is the sample drawn from the reference disribution? The test statistic $D_n$ is defined as\n",
    "\n",
    "$$\n",
    "D_n = \\sup_x|F_n(x) - F(x)|\n",
    "$$\n",
    "\n",
    "where $F_n(x)$ is the emprical CDF of the sample and $F(x)$ is the reference CDF. $\\sup$ is the supremum of the set of distances (intuitively the largest absolute difference)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "db82bd03-1f9f-465f-9792-0ff72de6dd8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def empirical_cdf_factory(sample):\n",
    "    def empirical_cdf(x):\n",
    "        return (1 / len(sample)) * sum([int(X_i < x) for X_i in sample])\n",
    "    \n",
    "    return np.vectorize(empirical_cdf)\n",
    "  \n",
    "np.random.seed(0)\n",
    "sample = np.random.randn(100) + 6\n",
    "empirical_cdf = empirical_cdf_factory(sample)\n",
    "cdf = st.norm(loc=6, scale=1).cdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "58195f26-c258-4e4f-a1a2-fc5d72b98008",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAEWCAYAAAB2c65HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA2CElEQVR4nO3deXgUVdbH8e8BAsgumwsg4IaiaNTI5qtmBBR9RRRcUcfojOC8bjiiI4wyrrjhiLtGUXQUFcEFFQVcAo6AChKRRRBRMSg7iBC2JOf9ozqhCSEkkE6lO7/P85TdXXW7+lSIOX2XutfcHREREYk/VcIOQERERHaPkriIiEicUhIXERGJU0riIiIicUpJXEREJE4piYuIiMQpJXGRXTCz9WZ2YDHHnzaz2/bwM1LNLGtPzhEPzGyEmd1dAeLYy8zeNbPfzeyNsOMR2V1K4hKXzOwnM9sYSbD52+Ox+Cx3r+Pui4o5fpW73xWLz85ngevMbLaZbTCzLDN7w8zaRY6PMLMtZvZHZJttZveaWf2oc6SZWW55/MyiPu+/MTx/hpltilzHSjN708z2K+HbzwX2ARq5+3mxilEk1pTEJZ71iCTY/O2a8g7AzKqW00c9AlwPXAc0BA4F3gb+N6rMA+5eF2gCXA50BD43s9pRZaaG/TMrY9e4ex3gYKAOMLSE72sJLHD3nNJ+oJlVK+17RGJFSVwSTqQG+LmZPWxma81skZl1juz/xcyWm9llUeVHRJrEJ0ZqsZPMrGXUcTezg6PKPmVm48xsA/Cnwk3EZtbTzDLNbJ2Z/WBm3SP7LzezeZHPWGRm/Up4PYcAVwMXufsn7r7Z3bPd/RV3v69weXff5O5fAWcBjQgSeml+fheY2fRC+24ws7GR52eY2dzIdSwxswFFnONw4GmgU6SmvDbq8N5m9n7k/V+Y2UFR7zss8u+w2szmm9n5JYnZ3dcSfKlJ3tW5zOwOYDBwQSS2v0T2XxH591ljZuOL+B242sy+B76P7Dsz8u+81symmNlRUeV/MrMBZjbLgib7182sZtTxnf2O1Dez4Wb2W+Rne3c5flGUOKQkLomqAzCLIImNBF4DjieosV0CPG5mdaLKXwzcBTQGMoFXijl3H+AeoC6wXXOxmbUHXgJuAhoAJwE/RQ4vB84E6hEk1ofN7NgSXEsXIMvdvyxB2QLu/gcwETixNO8D3gXaRL485OtD8HMEGA70i9T6jwQ+KeKz5wFXsa3m3yDq8IXAHcDewEKCnyWRFoOJkc9pGin3pJm13VXAZtYI6BU5X7Hncvd/AUOA1yOxDTeznsCgyDmaAJ8Brxb6mLMJfq/amtkxwPNAP4LfsWeAsWZWI6r8+UB3oDVwFJAWia2435ERQA7B7+kxwKnAX3d1/VJ5KYlLPHs7UgvK366MOvaju7/g7rnA60AL4M5ILXYCsIXgD2W+9919srtvBv5JUINssZPPfcfdP3f3PHffVOjYX4Dn3X1i5PgSd/8OwN3fd/cfPDAJmEDJEmwj4LcSlCvKrwTN7/k6FvqZdSz8BnfPBt4BLoKCloDDgLGRIlsJElk9d1/j7l+XMqa33P3LSFP2K2yrPZ8J/BT5d8tx95nAGKC4PutHzex3YCXBF7Brd/NcVwH3uvu8SFxDgOTo2njk+Gp33wj0BZ5x9y/cPdfdXwQ2E3RhFMTm7r+6+2qCL0b511nk74iZ7QOcAfR39w3uvhx4mOALiEiRlMQlnp3t7g2itmejji2Ler4RwN0L74uuif+S/8Td1wOrgf138rm/7GQ/BF8WfijqgJmdbmbTIs27awn+YDcu5lz5VgElHbBVWDOCa8k3rdDPbNpO3jeSSBInqIW/HUnuAL0JYv850vXQqZQxLY16ns22f4eWQIfoLxkELST7FnOu69y9PkFNd2+g+W6eqyXwSFTZ1YAR/Pzy/VKo/I2Fzt+C7X9ndnadO/sdaQkkAb9FnfMZgpYEkSJpgIZIoKDWHWlmb0hQiy1KcUv//QIcVHhnpJl1DPBngpr8VjN7myBR7MrHwBNmluLu03dZettn1gG6EmmuLqWJQBMzSyZI5jfkH4j0t/c0syTgGmAUUT+/KKVdIvEXYJK7dyttsO7+rQXjEp6IdFGU9ly/APe4e3HdKNHXk19+d362Rf6ORPZvBhrvzoA7qZxUExcJnGFm/2Nm1Qn6xqe5e3E17p0ZDlxuZl3MrIqZNTOzw4DqQA1gBZBjZqcT9Hfukrt/DzwJvGrB/eTVzaymmV1oZrcULm9mNczsOIKBXmuAF0p7Ee6+FXgDeJDgC83EyLmrm9nFZlY/UmYdkLeT0ywDmkd+piXxHnComV1qZkmR7XgLBsmVxIsEt42dtRvnehoYaGZHRK6zvpkV14z/LHCVmXWwQG0z+18zq1uCOIv8HXH33wi6WB4ys3qRYweZ2cklvH6phJTEJZ69a9vf8/zWHpxrJPAvgmbU4wgGv5VaZPDZ5QR9mb8Dk4CWkUFm1xHUWtcQNFGP3dl5inAd8DjwBLCWoDn2HIK+1nw3m9kfBM3vLwEzgM7uvmF3roXgZ9IVeKNQzfBS4CczW0fQl3zxTt7/CTAHWGpmK3f1YZGf0akEfcC/EjRH30/w5WeX3H0Lwa14t5X2XO7+VuT4a5Hrmg2cXsxnTQeuJPg3WUMwoC6thHEW+TsSOfxngi98cyPnHc3ud6VIJWDupW3xEkksZjaCYPT3rWHHIiJSGqqJi4iIxCklcRERkTil5nQREZE4pZq4iIhInIq7+8QbN27srVq1CjsMERGRcjFjxoyV7t6kqGNxl8RbtWrF9Oklnu9CREQkrpnZzzs7puZ0ERGROKUkLiIiEqeUxEVEROJU3PWJF2Xr1q1kZWWxaVPhVSGlIqlZsybNmzcnKSkp7FBERBJCQiTxrKws6tatS6tWrTAryaJQUt7cnVWrVpGVlUXr1q3DDkdEJCHErDndzJ43s+VmNnsnx83MHjWzhWY2K7J84G7ZtGkTjRo1UgKvwMyMRo0aqbVERKQMxbJPfATQvZjjpwOHRLa+wFN78mFK4BWf/o1ERMpWzJrT3X2ymbUqpkhP4CUP5n2dZmYNzGy/yJq6IiJSmaSnQ69e0LgxjBgRbLsyevT25TMygv1Dh8J77+36/dHlp06FMWOC1wMHBq8jHHAHz4O8PMhzyK3fiN8eH8PWrdBo6EBs9SoW3ZLO1q1gBp07l/TC90yYfeLNgF+iXmdF9u2QxM2sL0FtnQMOOKBcgouVRx99lKeeeopjjz2WV155JexwirRgwQL69+/P999/T926dTn44IN57LHHmDdvHj179uTAAw8kOzubffbZh5tvvpkzzzwTgNtvv51nn32WJk2CiYW6d+/OfffdF+aliEi8GDkSFiwIEupucoeNG2HrOqiWDbk5kJO/5QaPuTmQmwu5eXDXJZCdDafPgUNWwd+PDV5f9xskb9yWsPPydvysVcC5BwfPhwCNgH5vBa8bNIA1a3b7MkolpgugRGri77n7kUUcew+4z93/G3n9MfAPdy92OraUlBQvPGPbvHnzOPzww8ss7j3h7rg7VaoU3VNx2GGH8dFHH9G8efMSnS8nJ4dq1crvu9amTZto164d//73v+nRowcAGRkZNG7cmJUrVzJ06FDei3zDzczM5Oyzz2b48OF06dKF22+/nTp16jBgwICdnr8i/VuJSAXy7rvBY+TvTrT16+GXX2DxYli6FJYvD7YVK7Y9z982by7+Y6pUgdq1g61WrR2f16oFNWsGW40a27bo19WrQ1ISVKsWPOZv+a9r1oSTTiq7H42ZzXD3lKKOhVkTXwK0iHrdPLIv7vz000+cdtppdOjQgRkzZjBu3DhGjRrFqFGj2Lx5M+eccw533HEHV111FYsWLeL000/niiuuoG/fvlx77bXMnj2brVu3cvvtt9OzZ09GjBjBm2++yfr168nNzWXcuHE7LTd27Fiys7P54YcfOOecc3jggQcA+PDDDxk0aBC5ubk0btyYjz/+mA0bNhR5nmgjR46kU6dOBQkcIDU1FQiSebTk5GQGDx7M448/TpcuXWL6MxaRxLVxI8xv0YMFC+DnB4NkHb2tXr3je2rWhKZNg22ffeDII4PnjRoFNeH69YOt8PNatYLm7kQRZhIfC1xjZq8BHYDfy6I/vH9/yMzc07NsLzkZhg0rvsz333/Piy++SMeOHZkwYQLff/89X375Je7OWWedxeTJk3n66af58MMP+fTTT2ncuDGDBg3ilFNO4fnnn2ft2rW0b9+erl27AvD1118za9YsGjZsWGy5zMxMZs6cSY0aNWjTpg3XXnstNWvW5Morr2Ty5Mm0bt2a1ZH/A+65554iz1O7du2C65g9ezbHHXdciX82xx57LA8++GDB64cffpiXX34ZgPvvv5/TTjutxOcSkTiWnh5k0s6dYcoUGDQInnkG2rQJatkPPcTWnKC5OntD8LghO3jctAl+pDX3cQsLaEODBnDAAcHWufO25y1awP77Q5MmUKdOYiXj3RWzJG5mrwKpQGMzywL+BSQBuPvTwDjgDGAhkA1cHqtYykPLli3p2LEjABMmTGDChAkcc8wxAKxfv57vv/+ekwq1r0yYMIGxY8cyNNIHtGnTJhYvXgxAt27daNiw4S7LdenShfr16wPQtm1bfv75Z9asWcNJJ51UcD/2rs6zJ83bhbtjbrjhhmKb00UkQY0cCfXqwdixQNDv/PVXMOldyH0bTp4Om6KauqtYUCuuVxf23Rfa1vqRbr0nUW9AG+rVC+cS4lEsR6dftIvjDlxd1p+7qxpzrETXZt2dgQMH0q9fv2Lf4+6MGTOGNm3abLf/iy++2OF8OytXo0aNgtdVq1YlJyen1J8X7YgjjmDSpEnFxh1t5syZ6uMWqeTcgxr1ut/WMfhK+PLLzsyenUHef4PjrVv3YObZPTjuOGjbFg4/HFq2hKpVw407EWju9Bg47bTTeP7551m/fj0AS5YsYfny5UWWe+yxxwpqszNnztzp+UpSLl/Hjh2ZPHkyP/74I0BBc3pJztOnTx+mTJnC+++/X7Bv8uTJzJ6945w9s2bN4q677uLqq8v8u5iIVHCrV8OoUfCXvwRN3V99BfMXwJtvBk3et94a3OW1fDksWgSvvQY33QT/+79w4IFK4GUlIaZdrWhOPfVU5s2bR6dOnQCoU6cOL7/8Mk2bNt2u3G233Ub//v056qijyMvLo3Xr1gUjv3enXL4mTZqQnp5Or169yMvLo2nTpkycOLFE59lrr71477336N+/P/379ycpKYmjjjqKRx55hJUrV/LZZ59xzDHHkJ2dTdOmTXn00Uc1qE2kEsjJgS+/hPHjg+2rr4Jbr+rXh65d4dBawcCxldPUV12eYnqLWSxU9FvMpHj6txKJH3l5MHly0N09enRw73OVKtC+PZx6Kpx2WvC8WjUgchcLhe5ikT1XUW8xExGRCsY9uMNn5Eh49VVYsiS4f/rss+Hi4+fT9Y1+JD0wZNso9K6DgjdmZga38ki5UhIXEREWLQoS98iRMG9eULs+/fRgArUePYJEznzgrZ2cIDkZ+vQpv4AFUBIXEam08vKC/u2HH4aJE4N9J54ITz8N554bTJxSIH9Gtejm8s6d1XweMiVxEZFKJjsb/vOf4Jbc776D/faDu+6CP/85GGlepIceCh6LmBZVwqMkLiKS6NLTYeRINm+BX5fAQ6vSeGJDGqcctZKMw86l4T03ktSrB8yfD6k7md9Cfd4VkpK4iEiC2/DsSPgmk+k5ybjDocfApGFw4uFg5xGZS3MX1OddISmJV3CdO3dmypQppT5WnOJWG3vppZd44IEHMDOqVavGxRdfzIABA0hLS2PSpEnUq1ePjRs30rFjR4YMGVKwGlurVq2oW7cuVSMzODz55JN0Lq8FdUWkSIsXw+DBkDYdqlZJZsw1GVx3HaQenF+i8fZ92m3aqI87ziiJV3BFJen85Ul3J4EX54MPPmDYsGFMmDCB/fffn82bN/PSSy8VHH/wwQc599xzcXeGDRvGKaecwuzZs6levTpAwcIuIhKuNWvg3nvh0UeD17c1hwNawomPhhuXlD1Nu1pGXn75Zdq3b09ycjL9+vUjNzcXCGZru+mmmzjiiCPo2rUrX375JampqRx44IGMjSwUMGLECHr27ElqaiqHHHIId9xxR8F569SpAwTLgJ544omcddZZtG3bdrtjEKwY1q5dO44++mhuueUWAJ599lmOP/54jj76aHr37k12dnax13DvvfcydOhQ9t9/fwBq1KjBlVdeuUM5M+OGG25g33335YMPPtjdH5mIlLFNm4Jbwg46KHi84IKgm/uggyBJVbaElJhJPDV111tkJa+C8iNGBM9Xrtyx7C7MmzeP119/nc8//5zMzEyqVq3KK6+8AsCGDRs45ZRTmDNnDnXr1uXWW29l4sSJvPXWWwwePLjgHF9++SVjxoxh1qxZvPHGGxSelQ6C5UkfeeQRFixYsN3+Dz74gHfeeYcvvviCb775hptvvhmAXr168dVXX/HNN99w+OGHM3z48GKvY3eWIf3uu+8KXv/pT38iOTmZDh06lPgcIrLncnPhpZeC1vCbboIOHWDmTHjxxWChEUlc+m5WBj7++GNmzJjB8ccfD8DGjRsL5kmvXr063bt3B6Bdu3bUqFGDpKQk2rVrx08//VRwjm7dutEoclNmr169+O9//0tKyvaz7LVv375gedFoH330EZdffjm1atUCti09Onv2bG699VbWrl3L+vXry3xt78JT9qo5XaT8zZgBf/1rMHg8JQVeeAFOOSXsqKS8JGYSL+3AjOjyjRuX+v3uzmWXXca99967w7GkpCQsshpAlSpVCpYOrVKlynbLhlqhFQMKv4btlzstibS0NN5++22OPvpoRowYQcYuruuII45gxowZnFLCvwAzZ87U4iciIdm4Ee64I2hUbNo0mCL1/PODuc13MHp0uccn5SMxm9PLWZcuXRg9enTBcqOrV6/m559/LtU5Jk6cyOrVq9m4cSNvv/02J5xwQonf261bN1544YWCPu/8pUf/+OMP9ttvP7Zu3VrQvF+cgQMHctNNN7F06VIAtmzZwnPPPbdDOXfn0Ucf5bfffitoZRCR8vPZZ3D00XD//XD55TB3Llx44U4SOASVE7WSJaTErImXs7Zt23L33Xdz6qmnkpeXR1JSEk888QQtS9EZ1b59e3r37k1WVhaXXHLJDk3pxenevTuZmZmkpKRQvXp1zjjjDIYMGcJdd91Fhw4daNKkCR06dOCPP/4o9jxnnHEGy5Yto2vXrrg7ZsYVV1xRcPymm27irrvuIjs7m44dO/Lpp58WjEwXkdj74w+45RZ48klo3Ro++ghK1BiWP+YnLS2G0UkYtBRpBTBixAimT5/O448/HnYoMRfv/1YiYVn6P715bnYnBq8bwPXXw9DpqUSmZdi5M8+EAQO0TGic01KkIiJx6vff4brr4KzPoUZT+Pxz6NQJSC3FSZS8E5aSeAWQlpZGmpq5RKSQb74JVhP78Udo8c8x3HYbRMbGKjELkEAD2+KtW6Ay0r+RSMmNGAEdO8KGDfDjhQO5O3fgtgQuEpEQSbxmzZqsWrVKSaICc3dWrVpFzZo1ww5FpELbuDG47/vyy4PlumfOhBZZU2Hq1LBDkwooIZrTmzdvTlZWFitWrAg7FClGzZo1CxZMEZEdLVwYNJ9/8w3ceivcfju7HrwmlVpCJPGkpKQiZzITEYkXb70V3AFWtSq8/z6ccUbYEUk8SIjmdBGReLV1azDfea9ewdznM2cqgUvJJURNXEQkHq1bFzSfT5wI/23blw4nN6Jay8j0zb17w6pVwfPMTEhODitMqcCUxEVEQvDrr0GNe84ceP55OGEqO/+LnJwMffqUY3QSL5TERUTK2dy50L07rFkDi7r1pcVUID19+0JjxoQSm8QX9YmLiJSjyZPhhBOCvvDJk6FF9gJYsCDssCROqSYuIlIe+vZl+ecL8HkwviYc1QZq3oD6u2WPKImLiJSDWd/C6rlQvx4c2Q6S8v/6qr9b9kBMk7iZdQceAaoCz7n7fYWOHwC8CDSIlLnF3cfFMiYRkfKUd2Vfpk2DE2an07s3vPwyJGniQikjMesTN7OqwBPA6UBb4CIza1uo2K3AKHc/BrgQeDJW8YiIlLetW2Hs542YNLsR/fvDqFGgmYelLMVyYFt7YKG7L3L3LcBrQM9CZRyoF3leH/g1hvGIiJSbrVuDVvJz5t1LtQfu5eGHoYqGEksZi2VzejPgl6jXWUCHQmVuByaY2bVAbaBrUScys75AX4ADDjigzAMVESlL+Ql89Gj497/hhhvCjkgSVdjfCy8CRrh7c+AM4D9mtkNM7p7u7inuntKkSZNyD1JEpKS2boWLLopK4P/tHcy+JhIDsUziS4AWUa+bR/ZF+wswCsDdpwI1gcYxjElEJGbyE/iYMVE18FWrtk2fKlLGYpnEvwIOMbPWZladYODa2EJlFgNdAMzscIIkrvVERSTuRCfwhx9WE7qUj5glcXfPAa4BxgPzCEahzzGzO83srEixG4Erzewb4FUgzd09VjGJiMTC1q1w4YXbEnj//mFHJJVFTO8Tj9zzPa7QvsFRz+cCJ8QyBhGRWMp5Kp0Zd4/nzV/HMGwYXL90IKRO3VZAM7JJDIU9sE1EJG7l5MArT66jwa9zggR+fRGFNCObxJDFW+t1SkqKT58+PewwRKSSc4e//jVYRlRN6BJLZjbD3VOKOqaauIjIbrj11iCBDx6sBC7hURIXESmlxx6DIUOgb1+4fVIqpKaGHJFUVkriIiKl8PrrQd/32WfDE0+AhR2QVGpK4iIiO5Oevu350KGsOTqVfS9KZUa9VEavSqVa19Rg9LlISJTERUR2ZuRIuC9YQTkrC2bPgVp7wZFHQtX8v54afS4hiul94iIicS0tDYAffoDOrw6gZrMBTJkCSc3CDUskn5K4iMjOpKWxbBmcdkJwT/j48dBMCVwqECVxEam80tODJvOd2HJUCn0m3cKvvzbmk0/gsMPKMTaRElCfuIhUXiNH7nRgmjvMf2U6h3z7JqNHQ8eO5RuaSEmoJi4ilVekz7vgMcrNN8HQycFtZGecUa5RiZSYkriIVF5FJG+AF16AoUPh6qvh//6vfEMSKQ01p4tI5bVyZbBF+ewz6NcPunaFYcPCCUukpFQTF5HK69xzg8eMDAB+/BF69YLWrWHUKKimv5BSwelXVEQSV3Gjz595Bm68seDlunXQo0dwK9m778Lee5dTjCJ7QElcRBJX/ujz5OSij/foAUBuLlx0EXz3XXAv+KGHlluEIntESVxEEltyckFz+c784x8wbhw8+SR06VIuUYmUCQ1sE5FKbfhweOihYCT63/4WdjQipaOauIgkrqg+76J89lmQuDUSXeKVkriIJK5In3dRsrKCwekaiS7xTL+2IpK45s8PHtu02W735s1w3nmwYQN8+qlGokv8UhIXkcTVr1/wWGhg29//DtOmBTXwtm3LPyyRsqIkLiKJa8iQHXa99FIwCn3AgKA2LhLPlMRFJHF17rzdy5kzg8p5airce284IYmUJd1iJiKJa8qUYANWr4bevaFRI3j9dQ1kk8SgX2MRSVyDBgGQ90kGF18cjEifPBmaNg05LpEyoiQuIgnvjjvgww/hqaegY8ewoxEpO2pOF5GEtmoV3HknXHbZtsHqIolCSVxEEtbGjTBvHhxzTFALNws7IpGyFdMkbmbdzWy+mS00s1t2UuZ8M5trZnPMbCdrBoqIlM7GjTB7DmAwZgzstVfYEYmUvZj1iZtZVeAJoBuQBXxlZmPdfW5UmUOAgcAJ7r7GzDTcRETKxHXXwcUboF07aNQ67GhEYiOWNfH2wEJ3X+TuW4DXgJ6FylwJPOHuawDcfXkM4xGRSuLll+G55+CAA6BRw7CjEYmdWCbxZsAvUa+zIvuiHQocamafm9k0M+te1InMrK+ZTTez6StWrIhRuCKSCObNCwawnXgitG4VdjQisRX2wLZqwCFAKnAR8KyZNShcyN3T3T3F3VOaNGlSvhGKSNzIzobzz4fateHVVzWQTRJfLO8TXwK0iHrdPLIvWhbwhbtvBX40swUESf2rGMYlIgnqmmtgzpzgnvBmzShy7nSRRBLLmvhXwCFm1trMqgMXAmMLlXmboBaOmTUmaF5fFMOYRCSRpKcHE6GnprL0sFRWvPAu//wnnNpyfrAfdpg/XSSRxCyJu3sOcA0wHpgHjHL3OWZ2p5mdFSk2HlhlZnOBT4Gb3H1VrGISkQQzciRkZrIhGxZ8D+2OhNtvDzsokfJj7h52DKWSkpLi06dPDzsMEakIUlPJzYWjVmewciVkZsJ++4UdlEjZMrMZ7p5S1DHNnS4iccsJauDzlsPEiUrgUvmEPTpdRGS3LV0Ky5bB4MHQpUvY0YiUP9XERSQuzZ4Nq7+HvRvAbbeFHY1IOFQTF5G4k50NF1wA1arC4YdD1aphRyQSjlLXxM2sNrDJ3XNjEI+IyC5df30wMxvDn6G67iCTSmyXSdzMqhDc430xcDywGahhZiuB94Fn3H1hTKMUEYl47bVgXvSBA6Hz5W3CDkckVCVpTv8UOIhgtbF93b2FuzcF/geYBtxvZpfEMEYREQAWLoS+fYP5W+68E3j33WATqaRK0pzeNTIt6nbcfTUwBhhjZkllHpmIVG7p6cFkLhF5Dr9/DeM2weG90qhWLQ0mTYLp06FHj/DiFAnRLmvi+QnczLoWPmZml0WXEREpM5HZ2PItWgR/rIfD2kCjRpGdt9wCffqEEp5IRVCagW2Dzaw3MACoAzxH0D/+YiwCExEhORkyMnjvvaCyfe21kPpo1PHGjYP2dZFKqjS3mJ0M/ABkAv8FRrr7ubEISkQkX1YWXHZZkM8feCDsaEQqltLUxPcG2hMk8uZASzMzj7fJ10UkPoweTU4O9DkfNm+G11+HmjXDDkqkYilNTXwa8KG7dye41Wx/4POYRCUi0rgxdz3VmM8+g6efhkMPDTsgkYqnNDXxru6+GMDdNwLXmdlJsQlLRCq7ef8YweIHIC0tjUt0E6tIkUoy2Usrd/8pP4FHc/fJZmZAM3fPikmEIlLpLF8Oa4aN4G97wRGPp4UdjkiFVZKa+IORWdveAWYAK4CawMHAn4AuwL8AJXER2WN5ecFAtk8tgy+nQe3aYUckUnHtMom7+3lm1pZg2tUrgP2AbGAeMA64x903xTRKEak0HnoIPvwQnnwSjjoq7GhEKrYS9Ym7+1zgnzGORUQquWnTYNAg6N0brlo/FIYCAwaEHZZIhbXL0elmNiTqebfYhiMildXatXDRRdC8ebDAib3/Hrz3XthhiVRoJbnFrHvU8/tjFYiIVF7u8Ne/BhO7vPYaNGgQdkQi8aHU64mLiJS1Z56BMWOCGdk6dAg7GpH4UZIk3tTM/g5Y1PMC7v7vmEQmIpXCrFnQvz+cfjrceGPY0YjEl5Ik8WeBukU8FxHZIxs2wAUXQMOGMGIEVCnNHJIiUqJbzO4ws6ZAtruvN7O9gBsJVjJ7JNYBikjiuuYauPu73hzVrxNNm0ZGoaemBo+ZmcGqJyKyUyX93vsa0DDy/A7gIGANMDIWQYlI4nv55aD2fdhhcMjBRRRITtZa4SK7UJJpVy8jSNqpkSlWLwAeANYTrGT2ZyDT3WfFNFIRSRjffQdLrxjIf5pDm2/HbP+XKCMjrLBE4k5J+sQzgA3ALKARsAx4l2Cg29WR47/HJjwRSTTZ2XDeefA0U0lpDtV0j4zIbitJn/jPZvYYMB7IA65098VmdgCwqqiFUURECvTuDatWFbxcPB8eWwrta2eSVCM5vLhEEkBJp119ysz+A+S5e3Zk9yrgophFJiIJZ+kyWLoUWh4ASa2T1ectsodK3JDl7usLvd5Q9uGISMIYODB4HDMGgLlz4fjj4fiT4aOP0FRTImVA/xuJSGxENaFv2BD0g9epA6++qn5wkbIS06kVzKy7mc03s4Vmdksx5XqbmZtZSizjEZFylJ4ebMDVV8O8efDKK7DffiHHJZJAYpbEzawq8ARwOtAWuCiyLnnhcnWB64EvYhWLiITnhRfgxRfhttuga9ewoxFJLLGsibcHFrr7InffQjBhTM8iyt1FsDraphjGIiLlrW9fVp/bl6uvhlNOgcGDww5IJPHEMok3A36Jep0V2VfAzI4FWrj7+8WdyMz6mtl0M5u+YsWKso9URMpc7rwFLBy3gHr1gmb0qlXDjkgk8YS23ICZVQH+TTAPe7HcPd3dU9w9pUmTJrEPTkT2iDssWADZG2HkSNh337AjEklMsUziS4AWUa+bR/blqwscCWSY2U9AR2CsBreJxL+nn4Zly6FVq6ApXURiI5Y3enwFHGJmrQmS94VAwcwO7v470Dj/tZllAAPcfXoMYxKRsjJwYHAbWWQEOn37woIF/L4O2mbCcdUyqd0yOcwIRRJezGri7p4DXEMwXes8YJS7zzGzO83srFh9roiEZ8sWmDMHatSAGh2SMc3IJhJT5u5hx1AqKSkpPn26KusiFU1ODnTrBtOmwdSpWgpcpKyY2Qx3L7KrWfMmiUjp9e4dPEamVAW45ZZgFdEXX1QCFykvSuIiUnpRU6oCjBoFDz0E//d/8Oc/hxSTSCUU2i1mIpIY5s6FK66ATp3g4YfDjkakclESF5Hdtm4dnHMO1K4Nb7wB1auHHZFI5aLmdBHZLQ6kpcEPP8DHH0OzZrt6h4iUNSVxEdktvyyGtyYFfeEnnxx2NCKVk5rTRaTUVq+BRT/CBRfADTeEHY1I5aWauIiUynffwco5QT/4c8+BWdgRiVReqomLSImtXg09eoBVgXZHQp06YUckUrmpJi4iJbJ1K5x3HixeDC0v6ETN5mFHJCJK4iKyS+5w7bXwySfBjGzN/3xv2CGJCGpOF5ESePxxeOYZ+Mc/NCObSEWiJC4ixZowAfr3h7POgiFDIjt79942f7qIhEbN6SKyo/R0GDmS7Gyo9TVM2wuOWQ1VTiFY5aRTp7AjFBGUxEWkKCNHkjczk2+3JBeMRK9WNer4gAGhhSYi2yiJi8gO8hxmVUnmpLwMPs2Amp3DjkhEiqIkLiLbcYfvF8DatfDcS9BZCVykwlISF5HtDB4MG5eeSZcucOmlYUcjIsVREheRAk88AXffDX/5ywC6Pxt2NCKyK0riIpVRZPR5tHcuHc211zZmWPIIrjlwKWa3hBSciJSU7hMXqYxGjoTMzIKXa9bC3/4W9H9fdRVUnfBhkOhFpEJTTVykskpOhowMMjPhpJOgxSEwdizUaJgG/dLCjU1ESkQ1cZFK7Mcf4fTToUEDGD8eGjYMOyIRKQ3VxEUqqS1b4dRTYfPmYGGT5lqVTCTuKImLVEK5ufDtLFiSCx9/DIcfHnZEIrI7lMRFKplNm+CxNWnM2wCjxmoadJF4pj5xkUSUng6pqdu2+fMB2DLmXRbsn8oDc87k5BfSOPPMMIMUkT2lJC6SiArdQgZBDfyee2D1Ghg2DC67LJTIRKQMqTldJFFFbiED2LgRzu4JEzN70HJ4D664ItTIRKSMKImLJLiNG+Hss2HiRBg+HC6/POyIRKSsxLQ53cy6m9l8M1toRczhaGZ/N7O5ZjbLzD42s5axjEekstm4EXr2VAIXSVQxS+JmVhV4AjgdaAtcZGZtCxWbCaS4+1HAaOCBWMUjUtnk5gUJ/KOPlMBFElUsa+LtgYXuvsjdtwCvAT2jC7j7p+6eHXk5DdB0EyK7Ej3yfMqUYN+UKduNQs/dvJXZs4ME/vzzSuAiiSqWfeLNgF+iXmcBHYop/xfgg6IOmFlfoC/AAQccUFbxicSn/JHnyclFHt6wAe7fdAvL1vzG8y9AWlp5Bici5alCDGwzs0uAFODkoo67ezqQDpCSkuLlGJpIxRQ18hwIlh/LyGD5cuhxCkyf1YPhSuAiCS+WSXwJ0CLqdfPIvu2YWVfgn8DJ7r45hvGIJLSFC6F7d/j1V3jzzaA/XEQSWyyT+FfAIWbWmiB5Xwj0iS5gZscAzwDd3X15DGMRSWhffEHB7GuffAIdO4Ybj4iUj5gNbHP3HOAaYDwwDxjl7nPM7E4zOytS7EGgDvCGmWWa2dhYxSOSqMaOhT/9CerXD8a3KYGLVB4x7RN393HAuEL7Bkc97xrLzxdJdE8/DVdfDccdB++9B02bhh2RiJQnzZ0uEocc+Oc/4W9/gzPOgE8/VQIXqYyUxEXizNY7hnBXzSEMGQJ9+8Jbb0Ht2mFHJSJhqBC3mIlIySxZAucP6syUKXD33TBoEJiFHZWIhEVJXKSiSU8PJnSJduONfLRXD247fz63rbuJGnfeQpd/dg4nPhGpMNScLlLRFFoL3IHXXoNTT4VGjSD12HV02Wd2aOGJSMWhmrhIRRSZkW3lSrjkEhg/Mnh8+uk21KqdEXZ0IlJBKImLVFBffAHnnQfLlsEzz8CVV6r/W0S2p+Z0kfKWng4jRgTPV64MVh97993g9fz5eGYmWUvgxBOhatVgApe+fZXARWRHqomLlLeRI2Hp0iJXJ1m5EpYmJfPowj6cdia89BLsvXf5hygi8UFJXCQM++4bPDZuDBkZuMPL/4H+/duwYUMGdz0AT98IVdRWJiLFUBIXCdnixXDVVfDBB8GKosOHw2GHhR2ViMQDfc8XCUleHjz5JBxxBEyeDI8+GjwqgYtISakmLhKC7I1w+p+CpN2tWzDWrVWrsKMSkXijJC5SjjZvhmWL4aefYFZ9eP75YHybRp6LyO5Qc7pIOcjLg5dfhjZtYNGP0LAhzJ0Ll1+uBC4iu09JXKQspacH930PHVqwa/XRqdzeagSXXgoH1lvJCbUzOfJI2G+/8MIUkcSgJC5SlqLmPZ8xI+jvnjULsrPhlVfgo48gKSUZ+vQJNUwRSQzqExfZXampO+7LzGTjYcn85esBvHpTsGDJN8MyuOcqqFEDILgvXESkLCiJi5RGfjP5gAE7HFq/HrJqJvPI9D68PStY6/vmm6F+/XKOUUQqDSVxkdKYOnXb84wM8vLg/ffh4Yfh0xlQqxZcfhV8PxCaNQsvTBGpHJTERUpjzBgANmyAF1+ERx6BBQugeXO4//5gpTHNdS4i5UVJXCRaenowOG0n1rbtxH317iU9Hdasgfbt4dVXoXdvSEoqxzhFRFASF9le/ujy5OSCXVu2worlsGw5jJ8ED1aBXr3ghhugUyfd5y0i4VESF4nWqBF06cL6F8fwzjvBbWETJkBuLhx5JFx8MfzYBw44IOxARUSUxEUKrFsHn1w6htdfh3eawsaNQbIeMCBI3u3ahR2hiMj2lMSl0srLg6+/hvHjg23qVMjJCaZEveyyYD6WE07Qmt4iUnEpiUul8ttvQfP4+PEwcSKsXBnsP+aYoMZ91c8Dad4cqj5wb7iBioiUgJK4JKzs7GDq0y+/hC++CB5P/TmdPozk6iS4bW9oeFhwS1j1Iw+Fe9Oh7ypYG3bkIiIloyQuCWHdOvjuO/j2220Je/bsYEAaBGt1d+gAt1Ubyf7LMqlyXDJFDipPTy/HqEVE9oySuMQNd1i2DObN23H79Ve4knROYzyj9x5D+/bwaO2BtKy7ir1eSqdpU6BvX1iZCccla/5yEUkISuJSYWzZAkuWwOLF228//7zt+YYN28rXrQuHHw5duwaPf315JI3nTKLXqsi92wOBVUDTqA9JTtYKYiKSMMzdY3dys+7AI0BV4Dl3v6/Q8RrAS8BxBH9uL3D3n4o7Z0pKik+fPj02Acvui5rpzAlGfq94cjRrqzUm6ZUR7P3OCCb+M4MVK+CQd4bSes57bNkKW7dQ8JiTu/0p/0QGTZvCoOpD6ZA3ldfPH8OBB8KZUway/09TqV6D7ZvE8ydpUS1bRBKImc1w95SijsWsJm5mVYEngG5AFvCVmY1197lRxf4CrHH3g83sQuB+4IJYxZTI8vKC/t+cnOCx8JaTA1u3bnssvOXkwObNwbZp07bn0dvGjUFNeMOGYNBY9PNhmSM5dGMm31ZJJic3aPo+t13wzewyIA248MIg1gEGZ1eD6kmQVB3q1gkek5KC5Tpr1oAaNSE7A/baCxgKTIWOD0cu9lfgtyJ+CKpli0glE7OauJl1Am5399MirwcCuPu9UWXGR8pMNbNqwFKgiRcTVFnWxMeMgRaXphZbxoGM2mfyfMMBuMN/slJ5s14ab9ZNY+/clTz227nblS3K6NppjK6TRoOclTy18lzS69zIR3v1oPWW+dy3tt92b/bCz4HHq9/IB9V6cFDOfIZt6se/qg1hCp3pkDuFO3MH7fI6BzGEqXSmE1MYwiD68QwLaMOZvMuNPLTL9+eX71nlXW60h7iq0Wg2123MRZtHcM7vI6haFQ5an8nihsk83juD+vWD5TcbNKDgef7rpk2D+7CrVt3lx4qICCHVxIFmwC9Rr7OADjsr4+45ZvY70AhYGV3IzPoCfQEOKMP5Lhs2hDp1dl2uWTM49rDged3foXUrOP5AqLsZ6m6g6FHOUHDg0IPg5IOD8ntnQPKRUKUl7PM7NPp8W1kr9D6L/OeUY2C/g2GftdDsIzjvZOjYAlr9Ci0ztr3XLNIXHHnM3zegJ6xqA/v/BIe/CQ/3hU0toflMOPidbe+rUiXYLPJYJbLv68egejuo9gHYv2HOaKAxMCKyAZBM2z59eLLvrn+eIiJSNmJZEz8X6O7uf428vhTo4O7XRJWZHSmTFXn9Q6TMyqLOCeoTFxGRyqW4mngsJ5RcArSIet08sq/IMpHm9PoE3agiIiKyC7FM4l8Bh5hZazOrDlwIjC1UZizBuCeAc4FPiusPFxERkW1i1ice6eO+BhhPcIvZ8+4+x8zuBKa7+1hgOPAfM1sIrCZI9CIiIlICMZ3sxd3HAeMK7Rsc9XwTcF4sYxAREUlUWmRRREQkTimJi4iIxCklcRERkTilJC4iIhKnYroASiyY2Qrg5zI8ZWMKzRCXoHSdiUXXmVh0nYmlrK+zpbs3KepA3CXxsmZm03c2E04i0XUmFl1nYtF1JpbyvE41p4uIiMQpJXEREZE4pSQO6WEHUE50nYlF15lYdJ2Jpdyus9L3iYuIiMQr1cRFRETilJK4iIhInKq0SdzMWpjZp2Y218zmmNn1YccUC2ZW08y+NLNvItd5R9gxxZKZVTWzmWb2XtixxIqZ/WRm35pZpplNDzueWDGzBmY22sy+M7N5ZtYp7JjKmpm1ifw75m/rzKx/2HHFgpndEPkbNNvMXjWzmmHHFAtmdn3kGueUx79lpe0TN7P9gP3c/WszqwvMAM5297khh1amzMyA2u6+3sySgP8C17v7tJBDiwkz+zuQAtRz9zPDjicWzOwnIMXdE3rSDDN7EfjM3Z8zs+pALXdfG3JYMWNmVYElQAd3L8sJrUJnZs0I/va0dfeNZjYKGOfuI8KNrGyZ2ZHAa0B7YAvwIXCVuy+M1WdW2pq4u//m7l9Hnv8BzAOahRtV2fPA+sjLpMiWkN/czKw58L/Ac2HHInvGzOoDJwHDAdx9SyIn8IguwA+JlsCjVAP2MrNqQC3g15DjiYXDgS/cPdvdc4BJQK9YfmClTeLRzKwVcAzwRcihxESkiTkTWA5MdPeEvE5gGHAzkBdyHLHmwAQzm2FmfcMOJkZaAyuAFyLdI8+ZWe2wg4qxC4FXww4iFtx9CTAUWAz8Bvzu7hPCjSomZgMnmlkjM6sFnAG0iOUHVvokbmZ1gDFAf3dfF3Y8seDuue6eDDQH2keafBKKmZ0JLHf3GWHHUg7+x92PBU4Hrjazk8IOKAaqAccCT7n7McAG4JZwQ4qdSHfBWcAbYccSC2a2N9CT4MvZ/kBtM7sk3KjKnrvPA+4HJhA0pWcCubH8zEqdxCN9xGOAV9z9zbDjibVIc+SnQPeQQ4mFE4CzIv3FrwGnmNnL4YYUG5FaDe6+HHiLoP8t0WQBWVGtRqMJknqiOh342t2XhR1IjHQFfnT3Fe6+FXgT6BxyTDHh7sPd/Th3PwlYAyyI5edV2iQeGfA1HJjn7v8OO55YMbMmZtYg8nwvoBvwXahBxYC7D3T35u7eiqBZ8hN3T7hv+mZWOzIQk0jz8qkETXgJxd2XAr+YWZvIri5AQg06LeQiErQpPWIx0NHMakX+9nYhGIeUcMysaeTxAIL+8JGx/LxqsTx5BXcCcCnwbaS/GGCQu48LL6SY2A94MTLytQowyt0T9varSmAf4K3g7yDVgJHu/mG4IcXMtcArkabmRcDlIccTE5EvY92AfmHHEivu/oWZjQa+BnKAmSTuFKxjzKwRsBW4OtYDMivtLWYiIiLxrtI2p4uIiMQ7JXEREZE4pSQuIiISp5TERURE4pSSuIiISJxSEhcREYlTSuIiIiJxSklcRIplZp+aWbfI87vN7LGwYxKRQGWesU1ESuZfwJ2R6SSPIVioQ0QqAM3YJiK7ZGaTgDpAqrv/EXY8IhJQc7qIFMvM2hHMwb9FCVykYlESF5GdMrP9gFcI1oJeb2aJuIytSNxSEheRIplZLYJ1n29093nAXQT94yJSQahPXEREJE6pJi4iIhKnlMRFRETilJK4iIhInFISFxERiVNK4iIiInFKSVxERCROKYmLiIjEqf8H4NJa1y/ZMm0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(2, 9, 50)\n",
    "    \n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 4))\n",
    "ax.plot(x, cdf(x), \"b\", label=\"reference CDF\")\n",
    "ax.step(x, empirical_cdf(x), \"r-.\", label=\"empirical CDF\")\n",
    "ax.set_xlabel(\"$x$\")\n",
    "ax.set_ylabel(\"$F(x)\")\n",
    "ax.set_title(\"Empirical CDF vs the Reference\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "id": "e9cf677d-8aae-4c65-b577-3dcf9c5ea433",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ks_stat(cdf_1, cdf_2, x):\n",
    "    return np.abs(cdf_1(x) - cdf_2(x)).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "b3f43993-b5f6-4af0-9e36-c6013f8bda70",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.05824863872383251\n"
     ]
    }
   ],
   "source": [
    "# under the null hypothesis that the sample comes from the reference distribution,\n",
    "# sqrt(n) * D_n statisic follows kolmogorov distribution\n",
    "\n",
    "D_n = ks_stat(cdf, empirical_cdf, sample)\n",
    "print(D_n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "id": "ea24c49e-f044-40a2-b11e-2440dc94aec4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8667717341286234"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# KS two-sided test statistic distribution for finite sample size n\n",
    "1 - st.kstwo.cdf(D_n, len(sample))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "id": "5ead5940-0b66-4ee8-822b-192251cf26c0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8865884365301926"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 - st.kstwobign.cdf(np.sqrt(len(sample)) * D_n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "id": "2b1c1e99-db37-4fb7-a5ec-0b467d13bd56",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KstestResult(statistic=0.05824863872383251, pvalue=0.8667717341286234)"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "st.kstest(sample, cdf, alternative=\"two-sided\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af46f95f-22b1-402f-ad3c-5e98cebf7f39",
   "metadata": {},
   "source": [
    "## Two-Sample KS Test\n",
    "\n",
    "Two-sample KS Test is used to compare two samples using their empirical CDFs. Are both samples drawn from the same\n",
    "unknown distribution?\n",
    "The statistic is defined as\n",
    "\n",
    "$$\n",
    "D_{n,m} = \\sup_x|F_n(x) - F_m(x)|\n",
    "$$\n",
    "\n",
    "where $F_n(x)$ and $F_(m)$ are the empirical CDFs of the first and second sample respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "id": "22a905cc-4185-49dc-9010-1cfcbbf03cbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_1 = np.random.randn(100) + 6\n",
    "sample_2 = np.random.randn(100) + 6\n",
    "\n",
    "empirical_cdf_n = empirical_cdf_factory(sample_1)\n",
    "empirical_cdf_m = empirical_cdf_factory(sample_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "id": "1a061d0f-5eda-49a9-9386-d433aada0c5f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.10999999999999999\n"
     ]
    }
   ],
   "source": [
    "D_nm = ks_stat(empirical_cdf_n, empirical_cdf_m, sample_1)\n",
    "print(D_nm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "id": "50e3216e-7263-4899-9a94-763519af9382",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5439007965907301"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# https://github.com/scipy/scipy/blob/7f7db53324bc6c9b7e8eac0a42441529da7b8b0c/scipy/stats/_stats_py.py#L7786\n",
    "n = (len(sample_1) * len(sample_2)) / (len(sample_1) + len(sample_2))\n",
    "1 - st.kstwo.cdf(D_nm, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "id": "71ae6028-548a-46d8-a2aa-bd93808a79cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KstestResult(statistic=0.10999999999999999, pvalue=0.5439007965907301)"
      ]
     },
     "execution_count": 174,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# st.ks_2samp can also be used\n",
    "# asymp, because auto/exact mode is using some strange algorithm\n",
    "st.kstest(sample_1, sample_2, alternative=\"two-sided\", mode=\"asymp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "746ed05e-8948-4b02-af1a-7818cbf27d9d",
   "metadata": {},
   "source": [
    "For large samples, the null hypothesis is rejected at level $\\alpha$ if\n",
    "\n",
    "$$\n",
    "D_{n,m} > \\sqrt{-\\log\\left(\\frac{\\alpha}{2}\\right)\\frac{1 + \\frac{m}{n}}{2m}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "id": "75f49bf6-a454-408d-888a-0dc4d6494822",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ks_two_sample_ppf(alpha, n, m):\n",
    "    return np.sqrt(-np.log(alpha / 2) * ((1 + m / n) / 2 * m))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "id": "08939ebd-7a8e-4ace-bc36-40ac5630b43d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 178,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D_nm > ks_two_sample_ppf(0.5, 100, 100)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}