algorithm_documentation.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "PwcPFB386x5z",
        "scrolled": false
      },
      "outputs": [],
      "source": [
        "from random import randint\n",
        "from math import sqrt\n",
        "import scipy.stats as st\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YAinuzPL8oly"
      },
      "source": [
        "#### Create normally distributed sample with parameters $\\mu = 2, \\sigma = 7$\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "d7yW4s1v79xa"
      },
      "outputs": [],
      "source": [
        "W = st.norm(2,7)\n",
        "X = W.rvs(size=20000)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "YT_OhZW9OWkx",
        "outputId": "22c28934-3304-4ee3-9b6b-de41215d0069"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([  2.57507309,  -8.31634145,   6.35990358, ...,   8.0748093 ,\n",
              "       -12.71109677,   7.54867633])"
            ]
          },
          "execution_count": 24,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "X"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CzFXu8BR-Ay4"
      },
      "source": [
        "#### Let $m$ be a randomly selected natural number from $[1, n]$. Then suppose that: $Am=\\dfrac{1}{n+\\sqrt{n}}\\sum ^{n}_{i=1}X_{i}-\\dfrac{1}{1+\\sqrt{n}}X_{m}$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "rAWx5nIP9IYm"
      },
      "outputs": [],
      "source": [
        "n = len(X) \n",
        "m = randint(0,n-1) \n",
        "Xm = X[m]\n",
        "Am = np.sum(X) / (n + sqrt(n)) - Xm/(1+sqrt(n))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "3BWATA34OWk0",
        "outputId": "2b5bddb7-9f68-4e5b-a77e-33f7143abaac"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "1.9208903225894802"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Am"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kKuOSBSq_TP4"
      },
      "source": [
        "#### Let $Y_j$ be a new random variable with the following distribution: $Y_{j}=\\begin{cases}X_{j}-Am,j=1,..,m-1;\\\\ X_{j+1}-Am,j=m,..,n-1\\end{cases}$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "cq04IGbZAS8h"
      },
      "outputs": [],
      "source": [
        "Yj = [X[j] - Am if j < m - 1 else X[j+1] - Am for j in range(n-1)]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "9JAvZnQYOWk2",
        "outputId": "895fd6df-d780-43ef-d782-651b51e268c9"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([  0.65418277, -10.23723178,   4.43901326, ...,   6.15391898,\n",
              "       -14.63198709,   5.62778601])"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "np.array(Yj)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tmN-YFohAuPj"
      },
      "source": [
        "#### Moving from random variable $X$ to $Yj$ allows us to eliminate parameter $\\mu$. The next step will be checking if $Yj$ in fact follows normal distribution with parameters $\\mu = 0, \\sigma = 7$. This will be done in two ways:\n",
        "* Overlaying $PDF$ of $N(0,7)$ over the histogram with $Y_j$ distribution\n",
        "* Using Shapiro-Wilk test for normality"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "id": "UrBZldj8BAjV",
        "outputId": "8bfc9c14-dc5d-4187-a64f-4e693883b282"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxV9Z3/8dfnnpsECBA0hDVAQDYBN0AE3FCqBkVwwbVVa/FnOyPdnJmOrQtq25laZ2o7U1trq9Zq3YobIosL7qxhT4BA2CQhCWFLQvbc+/n9kQMTYyA3kJtzl8/z8bgPzz3nm9z3V+MnJ99zzvcrqooxxpj44PM6gDHGmPZjRd8YY+KIFX1jjIkjVvSNMSaOWNE3xpg44vc6QFPdu3fXjIwMr2MYY0xUWbVq1T5VTWupXcQV/YyMDLKysryOYYwxUUVEdoXSzoZ3jDEmjljRN8aYOGJF3xhj4khIRV9EMkUkV0TyROS+Zo4nicir7vHlIpLR6NiZIrJURHJEZIOIdGi7+MYYY1qjxaIvIg7wJDAFGAHcIiIjmjSbCRxU1cHAE8Bj7tf6gReB76nqSGASUNdm6Y0xxrRKKGf644A8Vd2uqrXAK8D0Jm2mA8+723OAySIiwOXAelVdB6Cq+1U10DbRjTHGtFYoRb8vsLvR+3x3X7NtVLUeKAVSgaGAisgiEVktIj85+cjGGGNOVLjv0/cDFwDnApXAhyKySlU/bNxIRO4G7gbo379/mCMZY0z8CuVMvwDo1+h9uruv2TbuOH4KsJ+Gvwo+VdV9qloJzAdGN/0AVX1aVceq6ti0tBYfKDPGGHOCQin6K4EhIjJQRBKBm4G5TdrMBe5wt2cAi7VhdZZFwBki0sn9ZXAxsLFtohvT/jLSeyMiX3tlpPf2OpoxIWlxeEdV60VkFg0F3AGeVdUcEXkUyFLVucAzwAsikgccoOEXA6p6UER+Q8MvDgXmq+q7YeqLMWG3q6AInd31a/vlkSIP0hjTeiGN6avqfBqGZhrve6jRdjVwwzG+9kUabts0xhjjMXsi1xhj4ogVfWOMiSNW9I0xJo5Y0TfGmDgScYuoGBMpVJXVXx7iky0lbC0up7ouQOqVP+T5+gIynRX0lENeRzSm1azoG9OMj3P38qsFm9lcVI5PICM1mU5JDh0GjmF2/WU8Un87031f8C8J/yBd9nkd15iQWdE3ppHK2noeeCubN1YXMCC1E7+67gyuOrM3XTokACByEXkPDuPVwCU8H7ic92rG8vOE52hyR7MxEcuKvjGukvIabn92BZuLyvjB5CHcc8lpJPmdr7U7zVfIz3wvcbv/Pe6t/Sfurftnul3UHVWlYXJZYyKXXcg1BigqreamPy1l574Knvv2udx72dBmC35j6bKPlxN/wa3OB6RMuJEH3sqmYfYRYyKXFX0T98qr6/j2cyvYW17DCzPHMWlYj5C/1hHll/5nKV32D/6+/Eue+GBrGJMac/Ks6Ju4Fggq97y0hry9h/njt0YzNuPUVn8PETj0yfPcODad//lwK3NW5YchqTFtw4q+iWu/X5zHp1tK+Pk1o7hwyP9N632s2TSPJcmBx28cTdXOddz70nIS0zJs9k0TkexCrolby7fv54n3NnN448fc+thUbm1yvPnZNMua/V41AdDZndmrT3FVzX8w+K77eTfxZ3R61GbfNJHFir6JSxU19fz41bXUHSpi55gX6Tz2qwX+WMW9JT2klN8lPMmtdQ/wm/oZwFNtkNaYtmPDOyYu/fd7W9hTWs3++b+hs1S36fee6Gzkm84HPBO4ksQ+w9r0extzsqzom7izIb+Uvy7ZwTfP609NweawfMZ9/pfpxQFSM39AfSAYls8w5kRY0TdxJRBUfvrmerp3TuInmcPD9jldpIqHE54nMW0Af1/+Zdg+x5jWsqJv4spbawrILijj/qtOJ6VjQlg/6zLfKqp2ruWJD7ZwqLI2rJ9lTKis6JuYd+T2S19CEj985n1qCrdwzTn9wj5lgghUfPxnDh6uZtC079tC6iYi2N07JuYdWcz8D/VX8+v6Hvwj9U9MmN0FOPG7dEJVUbyLHyZ8xCtjr2LNxMUM8O11P9du5TTesDN9ExdKNZk/1k9nsm81E5xN7frZP/C/gZ8A/xu4tl0/15jmWNE3ceG5wBWU04l/8b/W7p/dUw7xLecD3ghcyI5gr3b/fGMas6JvYp4kduTZ+ilc5stihM+bO2m+559LInX8rv46Tz7fmCOs6JuY12X0VMpI5gf+Nz3LkCZl3OG8x9vBiXa2bzxlRd/EtIqaerqeew2TfGs5w7fD0ywz/QtIIMBfAld6msPENyv6Jqa9unI3TqcUvu/hWf4RPeQQ1zmfMSdwEb5OKV7HMXEqpKIvIpkikisieSJyXzPHk0TkVff4chHJcPdniEiViKx1Xzb7lGk3gaDy1yU7qc7fyBhfZCxucpcznxoS6TJ6qtdRTJxqseiLiAM8CUwBRgC3iMiIJs1mAgdVdTDwBPBYo2PbVPVs9/W9NsptzNc0nQO/y/CJfHmgkvKst72OdtRg3x4u82XRZfRVVNbWex3HxKFQzvTHAXmqul1Va4FXgOlN2kwHnne35wCTxVaINu3syENYR15XX38dfSmhcstSr6N9xXf983A6duX11QVeRzFxKJSi3xfY3eh9vruv2TaqWg+UAqnusYEiskZEPhGRC5v7ABG5W0SyRCSrpKSkVR0wpjnZwQEs1xHc4X8PNLJmuRwjW6gpyuOFpTttIXXT7sJ9IbcQ6K+q5wD3Ai+JyNeWI1LVp1V1rKqOTUtL+9o3Maa1nq2fQiequcn5yOsoXyMC5avfZUvxYZbvOOB1HBNnQin6BUC/Ru/T3X3NthERP5AC7FfVGlXdD6Cqq4BtwNCTDW3M8RzQLswLTmCG8ykpUul1nGZVbvqUlI4JvLB0l9dRTJwJpeivBIaIyEARSQRuBuY2aTMXuMPdngEsVlUVkTT3QjAiMggYAmxvm+jGNO+NwIXUksA3nQ+8jnJMWl/DjWPTWZRTRHFZ267cZczxtFj03TH6WcAiYBPwmqrmiMijIjLNbfYMkCoieTQM4xy5rfMiYL2IrKXhAu/3VNX+njVhowovBS5ljOQyzJfvdZzj+tb4AQRUeXmFLbJi2k9IUyur6nxgfpN9DzXargZuaObrXgdeP8mMxoRsuQ5nu/bhvxL+6HWUFg1ITebioWm8vOJLZl0yGL9jz0qa8LOfMhNTXq6fTFcqmOpb5nWUkNwyrj/FZTV8utXuWjPtw4q+iRm+jl1ZEBzHdc5ndJA6r+OE5NLhPejeOZHXVkb2UJSJHVb0TcxIHnUptSRwi7PY6ygtSnJAREj0O2xb/AoL1ufjJHezZRRN2FnRNzFBVel8xmWcI1sj/gIuQE2Ao08Or7hgKeL4+dOPp7CrwJZRNOFlRd/EhJw9ZSSmDWCG84nXUVptiK+Ac2QrrwYmeR3FxAEr+iYmzFmVj9bXMtVZ7nWUE3Kj8zFbNZ3E3vbsogkvK/om6tUFgsxdt4fKvOWkSIXXcU7IVGcZHamm85mXeR3FxDgr+ibqfZxbwoGKWio2fOh1lBPWRarI9K2k0/ALqakPeB3HxDAr+ibqvbE6n+6dE6naucbrKCflGucLnA6d+Wiz3bNvwseKvolqhypr+XDTXqaf3ReC0X2GfL4vm8Dhg7y1xubZN+FjRd9EtXfWF1IbCHLd6KZLPEQfvwSp2PQJizfvpbQyOh4uM9HHir6Jam+szmd4ry6M7BMbC41X5HxEbSDIguxCr6OYGGVF30StL/dXsubLQ1xzTvSf5R9RW7yNQWnJvGlDPCZMrOibqPXO+j0ATD0ztqYuuPbsvizfcYCCQ1VeRzExyIq+iToZ6b0REX75wkKq8zfS79RkRMTrWG0iyYF7r29YSnpE5m2ICCJic/KYNmNF30SdXQVFbHnwdBJ7DORXGSuOzmETC2oCUPejKsZKLuMvmXy0bzYnj2krVvRNVHonMB4fQa6M0mkXWjLNWUKu9mdLMHauV5jIYEXfRKV5wQmM922kh5R6HSUsMp0V+AjybmC811FMjLGib6JOQo9BbNc+XO1b6nWUsOkhpZzn28S84HhUvU5jYokVfRN1kk+/CD/1ZDorvY4SVlf5lrFN+5Kr/byOYmKIFX0TVVSV5NMv4kLfBk6Rw17HCaspNsRjwsCKvokqq788hD+lB1c7sTu0c0SqlDPRl8O8oBV903as6Juo8s66PWh9LZf5VnkdpV1c5VvGDu1NQo+BXkcxMcKKvokagaDy7oZCKretpIvEx9Oqmc5KHAIkD7/I6ygmRljRN1Eja+cBSsprqNz8uddR2s0pcpjzfdl0Gn4BarfxmDZgRd9EjQXZRST6fVRtz/I6Srua6ltGwim9yS4o8zqKiQEhFX0RyRSRXBHJE5H7mjmeJCKvuseXi0hGk+P9ReSwiPxr28Q28SYYVBblFHHRkDS0Nj6Gdo64wslCA3XMcyeYM+ZktFj0RcQBngSmACOAW0RkRJNmM4GDqjoYeAJ4rMnx3wALTj6uiVfr8g9RWFrNlFG9vI7S7lKkgqqda3l3Q6EN8ZiTFsqZ/jggT1W3q2ot8AowvUmb6cDz7vYcYLK40x6KyDXADiCnbSKbeLQwuwi/T/jG6T29juKJytwl5B+sImePDfGYkxNK0e8L7G70Pt/d12wbVa0HSoFUEekM/DvwyPE+QETuFpEsEckqKbFFoU2DI1Moiwi/f/tzyvOy6Jac6HUsT1RtXYbjExZm22yb5uSE+0Luw8ATqnrcRydV9WlVHauqY9PS0sIcyUSLXQVF6Oyu5Dw4ioRT+vDboatjZgrl1gpWlzN+0Km2jKI5aaEU/QKg8eQf6e6+ZtuIiB9IAfYD5wG/FpGdwI+An4nIrJPMbOLMwsA4fAS53Imvu3aayhzZi20lFWwtLvc6iolioRT9lcAQERkoIonAzcDcJm3mAne42zOAxdrgQlXNUNUM4LfAf6jq79sou4kTC4LjGOfbRKrEd7G7YmQvRLAhHnNSWiz67hj9LGARsAl4TVVzRORREZnmNnuGhjH8POBe4Gu3dRpzIvKCfdiq6UzxxfaMmqHo0bUDo/ufwgIr+uYk+ENppKrzgflN9j3UaLsauKGF7/HwCeQzcW5h8FwArojxaZRDNWVUL37x7ia+3F9J/9ROXscxUcieyDURbUFgHOfIVnrJQa+jRIQrRjY8p7Awxy7omhNjRd9ELH9KT3J0IFOcFV5H8VySAyJC/9RkaorymP3ntxARMtJ7ex3NRBkr+iZidRo2EYApPiv6NQHQ2V3R2V15IH0VSX2HU/jQAHYV2Pi+aR0r+iZidRp6PiNlB/189sBeY1e4F7UXBc71OImJRlb0TUQqLK0iqe9wG9ppxmDfHoZIPguCVvRN61nRNxFpkXtbYqbdqtmsKb4VrAiejq9jfD6hbE6cFX0TkRZkF1G7bxeDfTadcHOucFYSxEenIbZ+rmkdK/om4uw7XMPKnQeozF3idZSINUJ20V+Kj17sNiZUVvRNxHkvp5igYkX/OEQahr46DDiL0qo6r+OYKGJF30SchTlFDEjtRF3JDq+jRLRMZwXiJLB4c7HXUUwUsaJvIkppZR1L8vaRGYcrZLXW2bKN+vJ9LNhg9+qb0FnRNxHlg03F1AeVzJFW9FviE6Vyy1I+3VpCZW2913FMlLCibyLKwpwieqd04Kz0bl5HiQqVuUuorgvyca49wGZCY0XfRIyKmno+3VLCFSN74fOJ13GiQk1+DqnJiTbHvgmZFX3juSNr4fY4axI19UF+NesWRKzohyLJF2Tnknd4a0Ue4k84uqawTcRmjsWKvvHckbVwv3PdxaRSSsWdX8btWritVROA18euwZfUiQ9+ev7RSdlsIjZzLFb0TUSo1gQWB8/hMmcVjqjXcaLKRF8OXag4uuCMMcdjRd9EhCXBkVTQ8egMkiZ0iRLgG77VvB8YQ506XscxEc6KvokIC4Pn0oVKzvdlex0lKmU6KzlEF5YHT/c6iolwVvSN98TH+4GxTPatJlECXqeJShf71tGRahviMS2yom8816HfKA7ShUxb/PyEdZA6LvGtY1FgLEG1O5/MsVnRN57rOHQiHajhIt96r6NEtUxnBSWcwiod6nUUE8Gs6BtPBYNKp6ETmORbRyep8TpOVLvUt4ZE6lhoyyia47Cibzy1Zvch/F1SbWinDXSWai70bbCib47Lir7x1KKcIjRQxyW+NV5HiQmZvhUUkEZir8FeRzERyoq+8YyqsjC7iOpd60iRSq/jxITLnFU4BOg01FbUMs0LqeiLSKaI5IpInojc18zxJBF51T2+XEQy3P3jRGSt+1onIte2bXwTzTYVlvPlgUpbIasNdZMKJvg20mnoRFTtyWbzdS0WfRFxgCeBKcAI4BYRGdGk2UzgoKoOBp4AHnP3ZwNjVfVsIBP4k4j42yq8iW4LswvxCVTmLfc6SkzJ9K0gITWdLcWHvY5iIlAoZ/rjgDxV3a6qtcArwPQmbaYDz7vbc4DJIiKqWqmqR1Z36ADYqYc5amFOEedmnEqwstTrKDHlcicL1aBNt2yaFUrR7wvsbvQ+393XbBu3yJcCqQAicp6I5AAbgO81+iVwlIjcLSJZIpJVUmKLQcSDbSWH2VJ82JZFDIMeUkpN/iYWZBd6HcVEoLBfyFXV5ao6EjgX+KmIdGimzdOqOlZVx6alpYU7kokAi3IazkKvsGURw6JyyxdsLipn574Kr6OYCBNK0S8A+jV6n+7ua7aNO2afAuxv3EBVNwGHgVEnGtbEjoXZRZzVrxt9unX0OkpMqsxdCjQMoRnTWChFfyUwREQGikgicDMwt0mbucAd7vYMYLGqqvs1fgARGQAMB3a2SXITtfIPVrI+v5QrRvb0OkrMCpSXcGZ6CgtsXN800WLRd8fgZwGLgE3Aa6qaIyKPisg0t9kzQKqI5AH3Akdu67wAWCcia4E3gX9W1X1t3QkTXRZsaChEV51hS/qFU+aoXqzbfYg9h6q8jmIiSEi3T6rqfGB+k30PNdquBm5o5uteAF44yYwmxszbUMiovl0ZkJrsdZSYljmyF79emMuinCLuPH+g13FMhLAnck272n2gknW7D3HVGX28jhLzBqV1ZljPLjbEY77Cir5pV0duI7ShnfBKckBEWP76n1i+fR9O8imICBnp9u893lnRN+3q3Q1FnNE3hf6pnbyOEtNqAqCzu7Ls0nWI+Pjbjyehs7uyq8DO+uOdFX3Tbo4M7VxpZ/ntZpjsZqAUMj843usoJkJY0TftxoZ22p8ITPUtY2lwBCXa1es4JgJY0Tft5j/+toCawq0M6J6MiBx9mfCa6iwliI+FgXFeRzERwGa8NO1i94FK6D6Q2f6X+N7sr55xyiNlHqWKD8N8+QyRfN4JTADe8DqO8Zid6Zt2MX+DO7Tjs2mUvXC1s5SVOgync6rXUYzHrOibdjF/QyE1hVvo57NZVL0w1bcUxUen4ed7HcV4zIq+CbvdBypZl19K5ebPvY4Stwb5ihghO0kefpHXUYzHrOibsHvXHdqpzP3C4yTxbaqzlKS+wxuur5i4ZUXfhN289Xs4Mz2F+tJir6PEtat9y4D/+yVs4pMVfRNWeXsPk11QxrSzbK4dr/XzlVCzJ5d56/d4HcV4yIq+Cau56/YgghX9CFGx+TOyC8rYYStqxS0r+iZsVJW31xYw8bRUenT92iqZxgNHLqbPW2dn+/HKir4Jm3X5pezaX8n0s/t6HcW4AuX7ODfjFOatt3H9eGVF34TNW2sKSPT7yBxli59HkqvP6kNucTmbCu1J6HhkRd+ERX0gyLz1hVw6rAddOyR4Hcc0MvXMPvh9wptrCryOYjxgRd+ExdLt+9l3uIZrzrELuJHm1OREJg1L4+21BQSC6nUc086s6JuweGvNHrok+Zk0rIfXUUwzrj0nneKyGpZu2+91FNPOrOibNpcxoD//WLqFwpUL6JjotymUI9Dk03vQJcnPG2vyvY5i2pkVfdPmSjr2x5fUibfHrEBndz36MpGjQ4LDlWf0ZmF2EZW19V7HMe3Iir5pc8kjLyGNg4z3bfQ6ijmOa0f3pbI2wPsbbXqMeGJF37SpfYdr6DhoLNc5n+OIXSSMNEkOR4fbJpyWRn3pXu7+j78gImSk2zKW8cBWzjJt6q01BYjj53rnU6+jmGbUBPjKUNvjdUv4Y8o0ih/qR89Hd3uYzLQXO9M3bUZVmbMqn5o9Wxjqs3vAo8G1zucE8TE3MNHrKKadhFT0RSRTRHJFJE9E7mvmeJKIvOoeXy4iGe7+y0RklYhscP95advGN5EkZ08Zm4vKOZz9gddRTIgG+/ZwhmznzcAFXkcx7aTFoi8iDvAkMAUYAdwiIiOaNJsJHFTVwcATwGPu/n3A1ap6BnAH8EJbBTeRZ86qfBIdH5UbP/E6immF651PydGBJPQY6HUU0w5COdMfB+Sp6nZVrQVeAaY3aTMdeN7dngNMFhFR1TWqemQ6vxygo4gktUVwE1lq64O8vbaAy0b2JFhj0/ZGk2ucL0ikls5nXu51FNMOQin6fYHGV3jy3X3NtlHVeqAUSG3S5npgtarWNP0AEblbRLJEJKukxBbOjkaLNxdzsLKOGWPSvY5iWqmbVJDpW0nyyEuorgt4HceEWbtcyBWRkTQM+Xy3ueOq+rSqjlXVsWlpae0RybSxOavy6dEliQsHd/c6ijkBNzkf43TozKKcIq+jmDALpegXAP0avU939zXbRkT8QAqw332fDrwJ3K6q2042sIk8e8ur+Si3hGtH98Xv2A1h0WiCbyN1h4p4daXdthnrQvk/dCUwREQGikgicDMwt0mbuTRcqAWYASxWVRWRbsC7wH2q+kVbhTaR5R9Z+QSCyo1j+7Xc2EQknygV699nybb9fLm/0us4JoxaLPruGP0sYBGwCXhNVXNE5FERmeY2ewZIFZE84F7gyG2ds4DBwEMistZ92bSLMSQYVF5e8SUTBqVyWlpnr+OYk3A4+wN8Aq9l2dl+LAvpb3FVna+qQ1X1NFX9pbvvIVWd625Xq+oNqjpYVcep6nZ3/y9UNVlVz2702hu+7pj2lJHem+TBY8k/WMXcJ35is2lGuUD5fi4emsacVfk2z34MswFYc8J2FRRx281Xk0opB2/Mttk0Y8BN5/ajqKyaT7bYuVmssqJvTpjT+VQ+DI5mhvMJiWK3+sWCyaf3JK1LEi8s3eV1FBMmVvTNCet85uUEcLjVWex1FNMGkhxI9DtsXfgcizcXk9Ctl82+GYOs6JsTEggqnc+6nAt96xngs6GAWHBkBs7Nl36BX5RH7rkWnd2VXQV2734ssaJvTshHm/fi79qDW50PvY5i2lgvOcgVvixeDUyiShO9jmPamBV9c0L+umQn9eX7+IZvtddRTBjc5rxHKZ15JzDB6yimjVnRN622pbicz/P2Ub76XRLsAm5MGu/bxFDZzd8CNglbrLGib1rtuS92kuT3cXjdIq+jmDARgduc98nWgST2GeZ1HNOGrOibVjlUWcuba/K55uy+BKvKvI5jwuha53M6U0nX0Vd7HcW0ISv6plVeXrGb6rogd16Q4XUUE2adpZqbnI/pdPqF7DlU5XUc00as6JuQ1QeCvLB0JxMGpTK8lz15Gw/u9C8E4LkvdnicxLQVK/qmRRnpvRERUkZdzJ7Sat5+/Ec2x06cSJd9VG7+jJdX7Kasus7rOKYNWNE3LdpVUETwoa5MuuYmBskeym/dbHPsxJGyFW9yuKaeV1Z86XUU0was6JuQfBY8gxwdyHedeThiMzDGk9ribYwfdCrPfbGTukDQ6zjmJFnRNyH5Y2AavdjPNc7nXkcxHrj7okEUllYzb/0er6OYk2RF37QosfdQlgZHcpd/AUlS73Uc44FJQ3swpEdn/vjxNoI2135Us6JvWpQy/gZSOMwtNs9O3PL5hFmXDmZL8WFbPD3KWdE3x7W1uJxOQydwh7OIZKnxOo7x0NQz+zCoezK/+3Crne1HMSv65rh+9+FWgrVVfNv/ntdRjMcc92x/c1E5H2wq9jqOOUFW9M0xbS4qY976Qsqz5nKqlHsdx0SAaWf1YUBqJ/5n8VZU7Ww/GlnRN8f0xPtb6JLkp2zlm15HMRHC7/i455LBZBeUsXizLZ4Tjazom2ZlF5SyKKeYmRcOJFh92Os4xkNJDojI0ddN5w2k7lARM5943cb2o5Df6wAmMv3m/S2kdEzgOxcM5MdehzGeOrKMYmNvBV7nR3X38M76PUw/u69HycyJsDN98zVZOw+wePNe7r5oEF07JHgdx0Sgab4l1BZv57/ey6W23p7SjSZW9M1XBIPKz9/dRM+uSdx5fobXcUyE8oly8JO/svtAFS8t3+V1HNMKVvTNURnpvekyahLrdh8i5+VfkZyUYLNpmmOq3rGaCYNS+d/FeRyusSe1o0VIRV9EMkUkV0TyROS+Zo4nicir7vHlIpLh7k8VkY9E5LCI/L5to5u2tqtoP2dPu5MRspOy67LQ2V1tNk1zTEkOvDH7dvZX1NLv8plHL/RmpPf2Opo5jhaLvog4wJPAFGAEcIuIjGjSbCZwUFUHA08Aj7n7q4EHgX9ts8QmbLqOnUYBaTzgfxGfzaRpWlATgJrvFnOd7zPSJlzLjgeHoLO7sqvApmmIZKGc6Y8D8lR1u6rWAq8A05u0mQ48727PASaLiKhqhap+TkPxNxGsuKyalAk3Mdm3monORq/jmChyX8LLJFLHo/W3eR3FhCCUot8X2N3ofb67r9k2qloPlAKpoYYQkbtFJEtEskpKSkL9MtOGfj5vI+L4ecj/N6+jmCjTQw7xA/+bLA6OZnHgbK/jmBZExIVcVX1aVceq6ti0tDSv48SdT7eUMG99IaVLX2OAz56yNK33bWchg2QPj9bfDo7d5hvJQin6BUC/Ru/T3X3NthERP5AC7G+LgCa8qusCPPh2NoO6J1O6fI7XcUyUSpQAj/j/yk7tRbeJN3sdxxxHKEV/JTBERAaKSCJwMzC3SZu5wB3u9gxgsdpsTFHhDx/lsWt/JT+/ZhQE7LY7c+IudLK53vcJXcfPILug1Os45hhaLD8ft98AAAynSURBVPruGP0sYBGwCXhNVXNE5FERmeY2ewZIFZE84F7g6G2dIrIT+A3wbRHJb+bOH+OR7IJS/vDxNq49py/nD+7udRwTAx5MeJFgZRk/mbPe1tONUCHNvaOq84H5TfY91Gi7GrjhGF+bcRL5TJjU1Ae497W1pHZO5OGrR3odx8SIblLB/vf/wMbO9/P0p9u555LBXkcyTUTEhVzTvjLSe9PrGzPZUnyY9c/cR7fkRHvy1rSZ4LalVGz+nF/Pzyax52n20FaEsaIfh4oklVPGX8/NzmKqbttqT96aNlUTgNyzXqCXU84F37mPiofS7KGtCGJFP84cqKil+7R/o7/s5X7/372OY2LUqVLOEwl/YIf24pH6272OYxqxoh9HgkHl3tfW4nTsyu8TfkcXqfI6kolhE52N/LMzl1cDlzAvcJ7XcYzLin4ceerTbXycW8KBD//MKJ9Nh2vC70f+1zlHtvLvdXeT0L2/13EMVvTjxkeb9/L4olymntmbw2sXeB3HxIkECfDHxN/SiWrSrn2A0so6ryPFPSv6cWBLcTnff3kNI/t05fEZZ3kdx8SZXnKQpxJ/iz8lje+/soaAravrKSv6MW7/4Rruej6LjokOf759LB0THa8jmTg0xreVA+8/xadbSnjo7WzsgX3vWNGPYeXVdXz7uZUUl1Xz9G1j6J3S0etIJo4dXreI7118Gn9f/iX/uzjP6zhxK6Qnck30qa4LcPffVrGpsIynbx/DOf1P8TqSMfx75jD2llfzm/e30L1zEreeZxd325ud6cegmvoAs15azdLt+yl6+3Emn97r6FOR9uSt8ZKI8Nj1ZzJpWBr3v7WB11bubvmLTJuyM/0YU10X4LsvrOKTLSXsX/Qkh29YDXz1aVt5pMybcMYACY6Pp741hrtfWMVPXl+Potx0rp3xtxc7048hFTX13PncSj7dWsJj159ht2aaiNUhweHp28Zw8dA0/v31DTzz+Q6vI8UNO9OPchnpvdlVUITTJZUe1z9EQloG+999gpt/9bHX0Yz5iiSHrw8vOgkMuPEBfj4PCg5W8cBVp+Pz2RBkOFnRj3K7CorY8OAZzKz9NyrowO8THmfSDeuBrjaMYyJKTYBmJ/aTRx/h4bfX8+wXO9hzqIr/vvEskpOsNIWLDe9EMVWl85mXMaP2YfwEmJP4CJOc9V7HMqZVknxBHp42igMfPs2CDQUMuefPJHTvZ1Mxh4kV/ShVUVPPva+tI3XKDxnj28pbSQ8y3Gd3Qpjoc+QvgLKrPublpP+kT/cUhtz1W/alDLeHuMLAin4UWrJtH1N+9xlvry3g0Gcv8reE/yRNbCjHRL+JzkbmJf2MkbKT7lf/K//04mpKymu8jhVTrOhHkdLKOn76xnpu/fNyfAIv/7/xlC55BUfsbMjEjl5ykJcTf8HBj55jce5eLnviE+asyidoc/a0CbtaEgVq64O8uGwX/7N4K2VVdXz34kH8+BtD6ZBg8+iY2OSXIDWrXmfHtuV0n/Ij/vUfdXz/d//gwIdP01sOsTO/0OuIUcuKfgSrCwR5e+0envwojx37KqjauYaDi5/lZ7/awc+8DmdMmNUEQGeVEdSf82bwAh7rczNJt/0XlVuWsiG/lDPSU7yOGJWs6EeYjPTefFlSSudRk+k67jr8KT2o3buDg588T+W3tiDDofETtnZbpol1PlGudz4j07eCvwSu5L/6T+Hq33/OpGFp3HXBICaelmr39reCFf0Ioaqszy+lfNQMhp89iUo6MEZyucf/Vy7ptxbf9jJEbPFyE7+SpYYf+t/kx398jf9duJ5nPt/Bt55ZzsDuydw6rj+/vPMKdm3/+uydA/r2suGgRqzoe0hVyS4oY0F2IQuzi9i+r4LkkZO4ylnGLc5izpE8bH40Y74qMVDFrEuHgJNA8rDzqT5nCr/cV4Fe/9/c6t/Elb5lXOFk0d29o00eKfI4cWSxot/OCkurWJK3nyXb9rNk2z4KS6txfML4Qacy88KB3HbRCB7/qf1nMeZYvvpk72pgNbnBdMZ/PJrC8y/i/vq7eKD+O4yUXVzg20CHAcupqg3YAkIuqy5htP9wDWMuuYqDCakk9RxMYq/B+FN6ABCsKqNq13qqtq+iausytleX89LRr7RhHGNaY5gvn0OfbeTDS+ayWfvxfnAsnwdG8UzgSnrePI1RDy9iaM8unN0vhTPTuzGqTwpDenaOyzvgQir6IpIJ/A5wgL+o6q+aHE8C/gaMAfYDN6nqTvfYT4GZQAD4gaouarP0Ee62Z1bA5B9zCpAhRYySbZztW8REXw4jH8tBZ3eBswGEI4XeLswac+JE4HTZzem+3fzA/yYVmkT3Fwfw6B9eYH1+KfM3FPHyioYn1x+cOoKZFwz0OHH7a7Hoi4gDPAlcBuQDK0VkrqpubNRsJnBQVQeLyM3AY8BNIjICuBkYCfQBPhCRoaoaaOuORKL7pgzn8m9MZu+de+kqVU2O2oMmxoRbstSgu7L4SebpR/f5u/UmscdAfvTsTu4q2fO1r4n1C7+hnOmPA/JUdTuAiLwCTAcaF/3pwMPu9hzg99Iwh+p04BVVrQF2iEie+/2Wtk380B2ZgripY/0HPlb7Tok+KmuDIe8H6Gp33Rjjma/P7lkBZCOPlDU762eHXxQdc4W5tqoXXv5ikZYmNBKRGUCmqt7lvr8NOE9VZzVqk+22yXffbwPOo+EXwTJVfdHd/wywQFXnNPmMu4G73bfDgFx3uzuw72Q6GKFitV9gfYtW1rfo1LhvA1Q1raUviIgLuar6NPB00/0ikqWqYz2IFFax2i+wvkUr61t0OpG+hTLhWgHQr9H7dHdfs21ExA+k0HBBN5SvNcYY005CKforgSEiMlBEEmm4MDu3SZu5wB3u9gxgsTaMG80FbhaRJBEZCAwBVrRNdGOMMa3V4vCOqtaLyCxgEQ23bD6rqjki8iiQpapzgWeAF9wLtQdo+MWA2+41Gi761gP3tPLOna8N+cSIWO0XWN+ilfUtOrW6by1eyDXGGBM7bBEVY4yJI1b0jTEmjkRc0ReRn4vIehFZKyLviUgfd7+IyP+ISJ57fLTXWVtLRB4Xkc1u/jdFpFujYz91+5YrIld4mfNEiMgNIpIjIkERGdvkWFT3DRqmInHz54nIfV7nORki8qyI7HWfrzmy71QReV9Etrr/PMXLjCdKRPqJyEcistH9efyhuz/q+yciHURkhYisc/v2iLt/oIgsd382X3VvuDk2VY2oF9C10fYPgKfc7SuBBTRMVDMeWO511hPo2+WA391+DHjM3R4BrAOSgIHANsDxOm8r+3Y6DQ/WfQyMbbQ/FvrmuLkHAYluf0Z4nesk+nMRMBrIbrTv18B97vZ9R342o+0F9AZGu9tdgC3uz2DU98+tfZ3d7QRguVsLXwNudvc/BfzT8b5PxJ3pq2rjGceS+b9JaqYDf9MGy4BuItK73QOeBFV9T1Xr3bfLaHhuARpNV6GqO4Aj01VEDVXdpKq5zRyK+r7RaCoSVa0FjkxFEpVU9VMa7rJrbDrwvLv9PHBNu4ZqI6paqKqr3e1yYBPQlxjon1v7DrtvE9yXApfSMP0NhNC3iCv6ACLySxHZDXwTeMjd3RfY3ahZvrsvWn2Hhr9cIPb61lgs9C0W+tCSnqp6ZDKYIqCnl2HagohkAOfQcEYcE/0TEUdE1gJ7gfdp+Av0UKOTyRZ/Nj0p+iLygYhkN/OaDqCq96tqP+DvwKzjf7fI0lLf3Db30/Dcwt+9S9p6ofTNRD9tGCeI6nu5RaQz8DrwoyajB1HdP1UNqOrZNIwSjAOGt/Z7eDL3jqp+I8SmfwfmA7OJkikdWuqbiHwbmApMdn/4IEb6dgxR0bcWxEIfWlIsIr1VtdAdNt3rdaATJSIJNBT8v6vqG+7umOkfgKoeEpGPgAk0DHX73bP9Fn82I254R0SGNHo7Hdjsbs8Fbnfv4hkPlDb6cy0qSMNiND8BpqlqZaNDsTxdRSz0LZSpSKJd46lU7gDe9jDLCRMRoWGGgE2q+ptGh6K+fyKSduSOPxHpSMMaJ5uAj2iY/gZC6ZvXV6SbuUL9OpANrAfeAfo2unL9JA1jWBtodIdItLxouIi5G1jrvp5qdOx+t2+5wBSvs55A366lYTyxBigGFsVK39w+XEnDnSDbgPu9znOSfXkZKATq3P9mM4FU4ENgK/ABcKrXOU+wbxfQMHSzvtH/Z1fGQv+AM4E1bt+ygYfc/YNoOJHKA/4BJB3v+9g0DMYYE0cibnjHGGNM+FjRN8aYOGJF3xhj4ogVfWOMiSNW9I0xJo5Y0TfGmDhiRd8YY+LI/wfF8Nmb4YIeqAAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Overlaying PDF of N(0,7) over the histogram of Y_j distribution\n",
        "\n",
        "W1 = st.norm(0,7)\n",
        "z = np.linspace(-21,21,1000)\n",
        "y = W1.pdf(z)\n",
        "plt.plot(z,y)\n",
        "plt.hist(Yj, bins=50, density=True, edgecolor='black')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ut2sO970CeVX",
        "outputId": "a0424b8b-a322-46b6-83cc-6ed8ee78a528"
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/scipy/stats/morestats.py:1676: UserWarning: p-value may not be accurate for N > 5000.\n",
            "  warnings.warn(\"p-value may not be accurate for N > 5000.\")\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "(0.9999418258666992, 0.9624379277229309)"
            ]
          },
          "execution_count": 10,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Shapiro-Wilk test\n",
        "\n",
        "st.shapiro(Yj)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hfNt8cyKCQ62"
      },
      "source": [
        "#### Elimination of parameter $σ$ will then be done by creating a new variable: \n",
        "$Z_{k}=Y_{k} / \\sqrt{B_{k}}$, where $B_{k}=\\frac{1}{n-k-1} \\sum_{j=k+1}^{n-1} Y_{j}^{2}, k=1, \\ldots, n-2 .$\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "wARSctDVFO_X"
      },
      "outputs": [],
      "source": [
        "Bk = [1 / (n - k - 1) * sum(list(map(lambda x: x**2, Yj[k:n-2]))) for k in range(n-2)]\n",
        "Zk = [Yj[k] / sqrt(Bk[k]) for k in range(n-2)]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "KvZJY1QXOWk5",
        "outputId": "5e7a78a8-9489-4d57-a319-6d931daa96d5"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([ 0.09394802, -1.47014545,  0.63749506, ...,  0.46613058,\n",
              "        0.67149322, -1.41421356])"
            ]
          },
          "execution_count": 12,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "np.array(Zk)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y2FtDdGBFs40"
      },
      "source": [
        "#### Let $F_k$ be a $CDF$ of $F_{n-k-1}(Z_k)$, which follows Student distribution $t_{n-k-1}$. If null hypothesis is correct, then $F_k$ should be independent and uniformly distributed on segment $[0,1]$. We will check it in two ways:\n",
        "* Overlaying $PDF$ of $U([0,1])$ over the histogram with $F_k$ distribution\n",
        "* Using Kolmogorov-Smirnov and chi-square tests"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "EuAitjp7F12N"
      },
      "outputs": [],
      "source": [
        "Fk = np.array([st.t.cdf(x=Zk[k],df=n-(k+1)-1) for k in range(len(Zk))])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "l3X4CdV5OWk6",
        "outputId": "f8f69e69-da28-4674-a7f3-686b4902e5f0"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([0.53742431, 0.07076905, 0.73809514, ..., 0.66355961, 0.71446112,\n",
              "       0.19591328])"
            ]
          },
          "execution_count": 14,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Fk"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "id": "b6mYZJl9GSgY",
        "outputId": "67732c59-df10-44c9-ae64-367297c87d5c",
        "scrolled": true
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAP9UlEQVR4nO3db4xc1XnH8e8Tu8ZqS2hVb1TkPyxRTRWLVAKtgChSoYJWxqrsF4kiW6ItlYuVtI4qEVWloiKJkxdNo6Z/VDfBUhFNpOA4eZFuEkdWS/nTIky9FEJiI0eLA3gdKIYSIhVhY3j6Yi50WM/s3N29M7Nz9vuRVsy99zD3OTNzf3Pm3DvjyEwkSaPvXcMuQJLUDANdkgphoEtSIQx0SSqEgS5JhVg5rB2vWbMmx8fHh7V7SRpJjz766IuZOdZp29ACfXx8nKmpqWHtXpJGUkQ8022bUy6SVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA11SEcbXXUxEnPc3vu7iYZc2MEP76r8kNemZU8+Tn3z3eevj088PoZrhcIQuSYUw0CWpEAa6JBXCQNfbPKkkjTZPiuptnlSSRpsjdKkBfrrRUuAIXWqAn260FDhCV08XrMDRpzQCHKGrpzNv4OizMr7uYp45Nbr97lb/JWt/madnnhtCRWqSgd5nHkBl6T618tMhVDN/Tg2VzUDvMw8gqWxzfWob9MBtWQT6UnrAm+YnAGm4ug3aYPADt2UR6EvpAW+anwAkvcWrXCQtSV7bP39FjdBH6QqEty4FnK2EqRKngXrzMeqt26fP1Z99vuOxoxqBHhF3Ab8NvJCZl3fYHsDfAluAV4GbM/O/mi603VzBPSpXIMz3UsBRerNyGqi3ph6j5fjG0P3YWXrH+aDVGaHfDfw98OUu228ENlZ/VwNfrP7bN6N+6dhcuo3codwX8XIMpab0+82zyQsKfJ77LzKzd6OIceDbXUbodwL3Z+Y91fJx4LrMnPMZmpiYyKmpqXkX/OlvHeWL93yLay9Zcd62B555Y17r39527bXnrX/k8MO8dubseetXX7CKq6/5QO16H3jggUZqXdD6Dv1qvNYO+5jzfhpoP0zDeowWVE8fn5u5/p9uxw7QXK19PhbmY67H6N+ffYM3O0Ts6v99nqf/+W8WtL+IeDQzJzpta2IOfS1wsm15plp3XqBHxC5gF8CGDRsa2HX/vHbmbJcXR+cXqkZPU2/aeqfux84bQ6hmuN7Mzm9i/3Lk1b7sr4kR+reBv8jM/6iW7wX+NDPnHH4vdIRe7aPr9MN81r+9rcNj0G0fqz/7U850eF12+9jYVK0LWl/juV10rfN47Jpq35SF7Lffj1Gjr68+PjcL3kdTtfb5WJiPno9RwzX1e4R+CljftryuWlckf9fk/80139/P+29q7nYp6vb68soO1dFEoE8CuyNiP62Toa/0mj9XGfp9tUFTb57zPYne7zeqhfDKDtVR57LFe4DrgDURMQN8EvgZgMz8EnCQ1iWL07QuW/z9fhXbD0vx4C3VqDzW3cITyg3QYT43o/K6GAU9Az0zd/TYnsAfNVbRgDnyGRwf66VrmM+Nr4vm+NV/qUDd/lGS5Wg5/QMtRX31X1LLIEa9ozJVspwuZDDQJS2IUyVLj1MuktSm2688/twFK5b8NJYj9EKV/GuOJfdNg7OQ301a6p9IDPRClfwFlVGaEx2VeeblqMQpIwO9IaNy4C7Ha6yHqcTQ0NJloDfEA1fSsHlSVJIKYaBLUiEMdEkqhHPoKsaonJiW+sVAVzE8Ma3lzikXSSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SClEr0CNic0Qcj4jpiLitw/YNEXFfRDwWEU9ExJbmS5UkzaVnoEfECmAvcCOwCdgREZtmNftz4EBmXgFsB/6h6UIlSXOrM0K/CpjOzBOZeRbYD2yb1SaBt/6pmIuAHzdXoiSpjjqBvhY42bY8U61r9yngpoiYAQ4CH+90RxGxKyKmImLq9OnTCyhXktRNUydFdwB3Z+Y6YAvwlYg4774zc19mTmTmxNjYWEO7liRBvUA/BaxvW15XrWu3EzgAkJkPA6uBNU0UKEmqp06gHwE2RsSlEbGK1knPyVltngWuB4iI99EKdOdUJGmAegZ6Zp4DdgOHgCdpXc1yNCL2RMTWqtkngFsi4nvAPcDNmZn9KlqSdL6VdRpl5kFaJzvb193RdvsY8MFmS5MkzYffFJWkQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhagV6RGyOiOMRMR0Rt3Vp85GIOBYRRyPiq82WKUnqZWWvBhGxAtgL/CYwAxyJiMnMPNbWZiPwZ8AHM/PliHhPvwqWJHVWZ4R+FTCdmScy8yywH9g2q80twN7MfBkgM19otkxJUi91An0tcLJteaZa1+4y4LKIeCgiDkfE5qYKlCTV03PKZR73sxG4DlgHPBgR78/Mn7Q3iohdwC6ADRs2NLRrSRLUG6GfAta3La+r1rWbASYz8/XM/BHwQ1oB/w6ZuS8zJzJzYmxsbKE1S5I6qBPoR4CNEXFpRKwCtgOTs9p8k9bonIhYQ2sK5kSDdUqSeugZ6Jl5DtgNHAKeBA5k5tGI2BMRW6tmh4CXIuIYcB/wJ5n5Ur+KliSdr9YcemYeBA7OWndH2+0Ebq3+JElD4DdFJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKUSvQI2JzRByPiOmIuG2Odh+KiIyIieZKlCTV0TPQI2IFsBe4EdgE7IiITR3aXQj8MfBI00VKknqrM0K/CpjOzBOZeRbYD2zr0O4zwOeA1xqsT5JUU51AXwucbFueqda9LSKuBNZn5nfmuqOI2BURUxExdfr06XkXK0nqbtEnRSPiXcAXgE/0apuZ+zJzIjMnxsbGFrtrSVKbOoF+CljftryuWveWC4HLgfsj4mngGmDSE6OSNFh1Av0IsDEiLo2IVcB2YPKtjZn5SmauyczxzBwHDgNbM3OqLxVLkjrqGeiZeQ7YDRwCngQOZObRiNgTEVv7XaAkqZ6VdRpl5kHg4Kx1d3Rpe93iy5IkzZffFJWkQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhagV6RGyOiOMRMR0Rt3XYfmtEHIuIJyLi3oi4pPlSJUlz6RnoEbEC2AvcCGwCdkTEplnNHgMmMvPXgG8Af9l0oZKkudUZoV8FTGfmicw8C+wHtrU3yMz7MvPVavEwsK7ZMiVJvdQJ9LXAybblmWpdNzuB73baEBG7ImIqIqZOnz5dv0pJUk+NnhSNiJuACeDznbZn5r7MnMjMibGxsSZ3LUnL3soabU4B69uW11Xr3iEibgBuB67NzDPNlCdJqqvOCP0IsDEiLo2IVcB2YLK9QURcAdwJbM3MF5ovU5LUS89Az8xzwG7gEPAkcCAzj0bEnojYWjX7PPDzwNcj4vGImOxyd5KkPqkz5UJmHgQOzlp3R9vtGxquS5I0T35TVJIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpELUCPSI2R8TxiJiOiNs6bL8gIr5WbX8kIsabLlSSNLeegR4RK4C9wI3AJmBHRGya1Wwn8HJm/grw18Dnmi5UkjS3OiP0q4DpzDyRmWeB/cC2WW22Af9U3f4GcH1ERHNlSpJ6icycu0HEh4HNmfkH1fLvAFdn5u62Nj+o2sxUy09VbV6cdV+7gF3V4q8CxxdY9xrgxZ6tymKflwf7vDwsps+XZOZYpw0rF17P/GXmPmDfYu8nIqYyc6KBkkaGfV4e7PPy0K8+15lyOQWsb1teV63r2CYiVgIXAS81UaAkqZ46gX4E2BgRl0bEKmA7MDmrzSTwe9XtDwP/lr3mciRJjeo55ZKZ5yJiN3AIWAHclZlHI2IPMJWZk8A/Al+JiGngf2iFfj8tetpmBNnn5cE+Lw996XPPk6KSpNHgN0UlqRAGuiQVYkkH+nL8yYEafb41Io5FxBMRcW9EXDKMOpvUq89t7T4UERkRI3+JW50+R8RHquf6aER8ddA1Nq3Ga3tDRNwXEY9Vr+8tw6izKRFxV0S8UH1Pp9P2iIi/qx6PJyLiykXvNDOX5B+tE7BPAe8FVgHfAzbNavOHwJeq29uBrw277gH0+TeAn61uf2w59LlqdyHwIHAYmBh23QN4njcCjwG/WC2/Z9h1D6DP+4CPVbc3AU8Pu+5F9vnXgSuBH3TZvgX4LhDANcAji93nUh6hL8efHOjZ58y8LzNfrRYP0/pewCir8zwDfIbWbwS9Nsji+qROn28B9mbmywCZ+cKAa2xanT4n8O7q9kXAjwdYX+My80FaV/11sw34crYcBn4hIi5ezD6XcqCvBU62Lc9U6zq2ycxzwCvALw2kuv6o0+d2O2m9w4+ynn2uPoquz8zvDLKwPqrzPF8GXBYRD0XE4YjYPLDq+qNOnz8F3BQRM8BB4OODKW1o5nu89zTQr/6rORFxEzABXDvsWvopIt4FfAG4ecilDNpKWtMu19H6FPZgRLw/M38y1Kr6awdwd2b+VUR8gNZ3Wy7PzDeHXdioWMoj9OX4kwN1+kxE3ADcDmzNzDMDqq1fevX5QuBy4P6IeJrWXOPkiJ8YrfM8zwCTmfl6Zv4I+CGtgB9Vdfq8EzgAkJkPA6tp/YhVqWod7/OxlAN9Of7kQM8+R8QVwJ20wnzU51WhR58z85XMXJOZ45k5Tuu8wdbMnBpOuY2o89r+Jq3RORGxhtYUzIlBFtmwOn1+FrgeICLeRyvQTw+0ysGaBH63utrlGuCVzHxuUfc47DPBPc4Sb6E1MnkKuL1at4fWAQ2tJ/zrwDTwn8B7h13zAPr8r8B/A49Xf5PDrrnffZ7V9n5G/CqXms9z0JpqOgZ8H9g+7JoH0OdNwEO0roB5HPitYde8yP7eAzwHvE7rE9dO4KPAR9ue473V4/H9Jl7XfvVfkgqxlKdcJEnzYKBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQvwfMHsgDnz9ISYAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Overlaying PDF of U([0,1]) over the histogram of F_k distribution\n",
        "\n",
        "W2 = st.uniform()\n",
        "z = np.linspace(0,1,1000)\n",
        "y = W2.pdf(z)\n",
        "plt.plot(z,y)\n",
        "plt.hist(Fk, bins=50, density=True, edgecolor='black')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "PiwqXuEWG6aU",
        "outputId": "6f4cad53-2fd5-4d4e-ad76-7615cc9fc941"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "KstestResult(statistic=0.003549961678092295, pvalue=0.9626416156745609)"
            ]
          },
          "execution_count": 18,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Kolmogorov-Smirnov test for uniform distribution\n",
        "\n",
        "st.kstest(Fk, 'uniform', args=(0, 1))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "dVG4JhLROWk8",
        "outputId": "3dc29030-e2f7-4f7e-99f9-4fc49cee76ca",
        "scrolled": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "Power_divergenceResult(statistic=3342.4016839312735, pvalue=1.0)"
            ]
          },
          "execution_count": 19,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Chi-square test for uniform distribution\n",
        "\n",
        "st.chisquare(Fk)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E-pM-0uCOWk8"
      },
      "source": [
        "#### All above-mentioned steps make up ${Parameter\\space Elimination\\space Method\\space (PEM)}$ test for normality:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6hog4rIaOWk9"
      },
      "outputs": [],
      "source": [
        "def pem_test(X: list):\n",
        "  \"\"\"\n",
        "    pem_test(X) carries out a test for normality based on Parameter Elimination Method\n",
        "\n",
        "    :param: X - sample for testing the method on\n",
        "\n",
        "    :returns: dictionary with test statistic ('stat') and p-value ('pvalue')\n",
        "  \"\"\"\n",
        "  n = len(X) \n",
        "  m = randint(0,n-1) \n",
        "  Xm = X[m]\n",
        "  Am = np.sum(X) / (n + sqrt(n)) - Xm/(1+sqrt(n))\n",
        "  Yj = [X[j] - Am if j < m - 1 else X[j+1] - Am for j in range(n-1)]\n",
        "  Bk = [1 / (n - k - 1) * sum(list(map(lambda x: x**2, Yj[k:n-2]))) for k in range(n-2)]\n",
        "  Zk = [Yj[k] / sqrt(Bk[k]) for k in range(n-2)]\n",
        "  Fk = np.array([st.t.cdf(x=Zk[k],df=n-(k+1)-1) for k in range(len(Zk))])\n",
        "  ks_test = st.kstest(Fk, 'uniform', args=(0, 1))\n",
        "  return {'stat': ks_test[0], 'pvalue': ks_test[1]}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FgoZd5EyOWk9"
      },
      "source": [
        "#### Performing $PEM$ test on normally distributed $X$ sample"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "1oqRI32UOWk9",
        "outputId": "34e1f15d-286c-40fc-f49c-c9d6f0f286eb"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "{'pvalue': 0.33730529232864115, 'stat': 0.0066618550455152525}"
            ]
          },
          "execution_count": 21,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pem_test(X)"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "name": "algorithm_documentation.ipynb",
      "provenance": []
    },
    "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.4"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}