TP1 Numpy and Maximum Likelihood in linear regression.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 95,
      "metadata": {
        "id": "vzqL8DwIn5tK"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import statsmodels.api as sm\n",
        "import scipy\n",
        "import math as m\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "p_Npo1BOn5tR"
      },
      "source": [
        "# Random numbers\n",
        "\n",
        "Numpy let you generate random numbers. Using the random number generators given as examples, implement the functions below."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "umjn0-5Kn5tW",
        "outputId": "77efcfdc-4ad5-4581-eefa-cce263fc717e"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "-0.022480621054556123"
            ]
          },
          "metadata": {},
          "execution_count": 13
        }
      ],
      "source": [
        "np.random.normal()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "NXE9KSFFn5tZ",
        "outputId": "3f1f8486-eca8-4b20-ed8d-11cbc330a482"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([ 9.28050528,  9.02849548, 11.44746155])"
            ]
          },
          "metadata": {},
          "execution_count": 14
        }
      ],
      "source": [
        "mu = 10\n",
        "variance = 2\n",
        "np.random.normal(mu, variance, size=(3, ))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "_3dQWUyEn5tc",
        "outputId": "274707b8-4fd0-423e-8051-3bfb6a7f46a8"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(9, 1)"
            ]
          },
          "metadata": {},
          "execution_count": 15
        }
      ],
      "source": [
        "p = 0.8\n",
        "n = 10\n",
        "np.random.binomial(n, p), np.random.binomial(1, p) # Binomial(n,p), Bernoulli(p)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ojGQ7wsAn5tg",
        "outputId": "4d12e06f-d63e-4c44-e475-8c0e539677e5"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([1, 1, 1, 1, 1])"
            ]
          },
          "metadata": {},
          "execution_count": 16
        }
      ],
      "source": [
        "np.random.binomial(1, p, size=(5, ))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WrslsPVVn5tk",
        "outputId": "8e0661ba-54ae-48bf-904e-4836f0f47b0f"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[array([[ 3.04133182, -0.7618349 , -2.09535178],\n",
              "        [-0.7618349 ,  3.31089133, -0.61026798],\n",
              "        [-2.09535178, -0.61026798,  3.42569592]]),\n",
              " array([[10.59244253, -3.97024679,  5.86676138],\n",
              "        [-3.97024679,  5.21469601, -2.65869232],\n",
              "        [ 5.86676138, -2.65869232,  5.70905089]]),\n",
              " array([[10.26221634,  0.30680997,  0.84552063],\n",
              "        [ 0.30680997,  9.68695397, -2.95650254],\n",
              "        [ 0.84552063, -2.95650254,  6.40178509]])]"
            ]
          },
          "metadata": {},
          "execution_count": 17
        }
      ],
      "source": [
        "from scipy import linalg\n",
        "\n",
        "def generate_wishart(n, covariance, k=1):\n",
        "\n",
        "    p = np.eye(3).shape[0]\n",
        "    Y = [np.random.multivariate_normal(np.zeros(p), covariance, size = n) for i in range(k)]\n",
        "    X = [np.matmul(y.T, y) for y in Y]\n",
        "    return X\n",
        "    \"\"\"\n",
        "    Return k Wishart matrices with degrees of freedom n\n",
        "    and given covariance.\n",
        "    \n",
        "    Parameters\n",
        "    ----------\n",
        "    n: int\n",
        "        degrees of freedom\n",
        "    covariance: np array of shape (p, p)\n",
        "        covariance matrix\n",
        "    k: int or None (optional)\n",
        "        number of wishart matrices to return\n",
        "    \n",
        "    Return\n",
        "    ------\n",
        "    W: np array of shape (k, p, p)\n",
        "        k Wishart matrices\n",
        "    \"\"\"    \n",
        "generate_wishart(5, np.eye(3), 3)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "JF20-d62n5tn",
        "outputId": "52f025dd-25e0-489b-a3d3-e7ce39772466"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "5"
            ]
          },
          "metadata": {},
          "execution_count": 18
        }
      ],
      "source": [
        "np.random.choice(6) + 1 # roll a die"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Y7H8eJqFn5tp",
        "outputId": "038199d3-4c81-4353-972f-7094e1c6e009"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([3, 5, 2, 3, 4, 6, 2, 3, 4, 1])"
            ]
          },
          "metadata": {},
          "execution_count": 19
        }
      ],
      "source": [
        "1 + np.random.choice(6, p =[1/6, 1/6, 1/6, 1/6, 1/6, 1/6], size=(10, ))  # roll a die 10 times"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ADpNm_rSn5tr",
        "outputId": "ee933c75-5665-4377-a8f1-ba2511967037"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([3, 6, 3, 5, 6, 6, 1, 3, 6, 6])"
            ]
          },
          "metadata": {},
          "execution_count": 20
        }
      ],
      "source": [
        "# roll a biased die 10 times\n",
        "1 + np.random.choice(6, p =[1/10, 1/10, 1/10, 1/10, 1/10, 1/2], size=(10, )) "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 513
        },
        "id": "SctiOWa9n5tv",
        "outputId": "9a13458d-edfd-4dd7-a0c9-acd080746a6e"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAARL0lEQVR4nO3dfYxldX3H8fengH1AKiBTRB5c0xISNOUhN4sGSlAEgRKwjWkhrUVLsmqwgdTEoE2k1X9sGrWpGMkWKNjiahVRUhHYIAmSKDK7XeS5IMGwK7KDqzyoiVn99o85mwzDHXbmnjtzd368X8nNPed3fuec79md/czZ3z33nFQVkqR2/dakC5AkLS+DXpIaZ9BLUuMMeklqnEEvSY3be9IFDHPQQQfVmjVrJl2GJK0amzZterqqpoYt2yODfs2aNUxPT0+6DElaNZL8cKFlDt1IUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxu026JMcnuT2JA8kuT/JxV37gUk2Jnmkez9ggfUv6Po8kuSCcR+AJOmlLeaMfifwwao6GngTcFGSo4FLgduq6kjgtm7+BZIcCFwGnACsBS5b6BeCJGl57Dboq+rJqtrcTT8HPAgcCpwLXNt1uxZ4x5DV3w5srKodVfVTYCNwxjgKlyQtzpK+GZtkDXAccBdwcFU92S36MXDwkFUOBZ6YM7+1axu27XXAOoAjjjhiKWVpT5eszH58iI401KI/jE3ySuB64JKqenbuspp9TFWvf2VVtb6qBlU1mJoaersGSdIIFhX0SfZhNuSvq6qvds1PJTmkW34IsH3IqtuAw+fMH9a1SZJWyGKuuglwFfBgVX1qzqIbgV1X0VwAfH3I6rcApyc5oPsQ9vSuTZK0QhZzRn8i8C7grUm2dK+zgE8ApyV5BHhbN0+SQZIrAapqB/Bx4O7u9bGuTZK0QlJ74AdYg8GgvE1xQ/wwVlp2STZV1WDYMr8ZK0mNM+glqXEGvSQ1zqCXpMbtkc+MVWclPsT0A0ypeZ7RS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS49r7ZqzfJpWkF/CMXpIat9sz+iRXA2cD26vqjV3bl4Cjui77Az+rqmOHrPs48Bzwa2DnQjfFlyQtn8UM3VwDXA58fldDVf3lrukknwSeeYn131JVT49aoCSpn90GfVXdkWTNsGXdg8P/AnjreMuSJI1L3zH6PwGeqqpHFlhewK1JNiVZ13NfkqQR9L3q5nxgw0ssP6mqtiX5A2Bjkoeq6o5hHbtfBOsAjjjiiJ5lSZJ2GfmMPsnewJ8DX1qoT1Vt6963AzcAa1+i7/qqGlTVYGpqatSyJEnz9Bm6eRvwUFVtHbYwyb5J9ts1DZwO3Ndjf5KkEew26JNsAL4DHJVka5ILu0XnMW/YJslrk9zUzR4M3JnkHuB7wDeq6ubxlS5JWozFXHVz/gLt7x7S9iPgrG76MeCYnvVJknpq7xYIkl5eVuK2J7Cqb33iLRAkqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY1bzKMEr06yPcl9c9r+Mcm2JFu611kLrHtGkoeTPJrk0nEWLklanMWc0V8DnDGk/dNVdWz3umn+wiR7AZ8FzgSOBs5PcnSfYiVJS7fboK+qO4AdI2x7LfBoVT1WVb8CvgicO8J2JEk99Bmj/0CS73dDOwcMWX4o8MSc+a1d21BJ1iWZTjI9MzPToyxJ0lyjBv3ngD8EjgWeBD7Zt5CqWl9Vg6oaTE1N9d2cJKkzUtBX1VNV9euq+g3w78wO08y3DTh8zvxhXZskaQWNFPRJDpkz+2fAfUO63Q0cmeT1SV4BnAfcOMr+JEmj23t3HZJsAE4BDkqyFbgMOCXJsUABjwPv7fq+Friyqs6qqp1JPgDcAuwFXF1V9y/LUUiSFpSqmnQNLzIYDGp6enq0lZPxFjPMSv2ZtXIsK3EcsHJ/L9qz+PMFQJJNVTUYtsxvxkpS4wx6SWqcQS9JjTPoJalxBr0kNW63l1dKmqOVK6H0suIZvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1LjdBn2Sq5NsT3LfnLZ/SfJQku8nuSHJ/gus+3iSe5NsSTLiI6MkSX0s5oz+GuCMeW0bgTdW1R8D/wd8+CXWf0tVHbvQI64kSctrt0FfVXcAO+a13VpVO7vZ7wKHLUNtkqQxGMcY/d8C31xgWQG3JtmUZN1LbSTJuiTTSaZnZmbGUJYkCXoGfZJ/AHYC1y3Q5aSqOh44E7goyckLbauq1lfVoKoGU1NTfcqSJM0xctAneTdwNvBXVcOflFBV27r37cANwNpR9ydJGs1IQZ/kDOBDwDlV9YsF+uybZL9d08DpwH3D+kqSls9iLq/cAHwHOCrJ1iQXApcD+wEbu0snr+j6vjbJTd2qBwN3JrkH+B7wjaq6eVmOQpK0oN0+M7aqzh/SfNUCfX8EnNVNPwYc06s6SVJvPhxcernyQecvG94CQZIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuEUFfZKrk2xPct+ctgOTbEzySPd+wALrXtD1eSTJBeMqXJK0OIs9o78GOGNe26XAbVV1JHBbN/8CSQ4ELgNOANYCly30C0GStDwWFfRVdQewY17zucC13fS1wDuGrPp2YGNV7aiqnwIbefEvDEnSMuozRn9wVT3ZTf8YOHhIn0OBJ+bMb+3aXiTJuiTTSaZnZmZ6lCVJmmssH8ZWVQG9ngJcVeuralBVg6mpqXGUJUmiX9A/leQQgO59+5A+24DD58wf1rVJklZIn6C/Edh1Fc0FwNeH9LkFOD3JAd2HsKd3bZKkFbLYyys3AN8BjkqyNcmFwCeA05I8ArytmyfJIMmVAFW1A/g4cHf3+ljXJklaIZkdXt+zDAaDmp6eHm3lZLzFDLNSf2atHMtKHAe0cyz+fC1NSz9fPSTZVFWDYcv8ZqwkNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1buSgT3JUki1zXs8muWRen1OSPDOnz0f7lyxJWoq9R12xqh4GjgVIshewDbhhSNdvV9XZo+5HktTPuIZuTgV+UFU/HNP2JEljMq6gPw/YsMCyNye5J8k3k7xhoQ0kWZdkOsn0zMzMmMqSJPUO+iSvAM4Bvjxk8WbgdVV1DPAZ4GsLbaeq1lfVoKoGU1NTfcuSJHXGcUZ/JrC5qp6av6Cqnq2q57vpm4B9khw0hn1KkhZpHEF/PgsM2yR5TZJ002u7/f1kDPuUJC3SyFfdACTZFzgNeO+ctvcBVNUVwDuB9yfZCfwSOK+qqs8+JUlL0yvoq+rnwKvntV0xZ/py4PI++5Ak9eM3YyWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxvYM+yeNJ7k2yJcn0kOVJ8m9JHk3y/STH992nJGnxej1KcI63VNXTCyw7Eziye50AfK57lyStgJUYujkX+HzN+i6wf5JDVmC/kiTGE/QF3JpkU5J1Q5YfCjwxZ35r1/YCSdYlmU4yPTMzM4ayJEkwnqA/qaqOZ3aI5qIkJ4+ykapaX1WDqhpMTU2NoSxJEowh6KtqW/e+HbgBWDuvyzbg8Dnzh3VtkqQV0Cvok+ybZL9d08DpwH3zut0I/E139c2bgGeq6sk++5UkLV7fq24OBm5IsmtbX6iqm5O8D6CqrgBuAs4CHgV+Abyn5z4lSUvQK+ir6jHgmCHtV8yZLuCiPvuRJI3Ob8ZKUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS40YO+iSHJ7k9yQNJ7k9y8ZA+pyR5JsmW7vXRfuVKkpaqz6MEdwIfrKrN3QPCNyXZWFUPzOv37ao6u8d+JEk9jHxGX1VPVtXmbvo54EHg0HEVJkkaj7GM0SdZAxwH3DVk8ZuT3JPkm0ne8BLbWJdkOsn0zMzMOMqSJDGGoE/ySuB64JKqenbe4s3A66rqGOAzwNcW2k5Vra+qQVUNpqam+pYlSer0Cvok+zAb8tdV1VfnL6+qZ6vq+W76JmCfJAf12ackaWn6XHUT4Crgwar61AJ9XtP1I8nabn8/GXWfkqSl63PVzYnAu4B7k2zp2j4CHAFQVVcA7wTen2Qn8EvgvKqqHvuUJC3RyEFfVXcC2U2fy4HLR92HJKk/vxkrSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9Jjev7cPAzkjyc5NEklw5Z/ttJvtQtvyvJmj77kyQtXZ+Hg+8FfBY4EzgaOD/J0fO6XQj8tKr+CPg08M+j7k+SNJo+Z/RrgUer6rGq+hXwReDceX3OBa7tpr8CnJrkJZ8zK0kar5EfDg4cCjwxZ34rcMJCfapqZ5JngFcDT8/fWJJ1wLpu9vkkD49Y10HDtj9WK/e7qpVjWf7jgHaOxZ+vpfLna9brFlrQJ+jHqqrWA+v7bifJdFUNxlDSxLVyLK0cB3gse6JWjgOW71j6DN1sAw6fM39Y1za0T5K9gVcBP+mxT0nSEvUJ+ruBI5O8PskrgPOAG+f1uRG4oJt+J/Ctqqoe+5QkLdHIQzfdmPsHgFuAvYCrq+r+JB8DpqvqRuAq4D+TPArsYPaXwXLrPfyzB2nlWFo5DvBY9kStHAcs07HEE2xJapvfjJWkxhn0ktS4ZoI+ydVJtie5b9K19JHk8CS3J3kgyf1JLp50TaNK8jtJvpfknu5Y/mnSNfWRZK8k/5vkfyZdSx9JHk9yb5ItSaYnXU8fSfZP8pUkDyV5MMmbJ13TKJIc1f197Ho9m+SSsW2/lTH6JCcDzwOfr6o3TrqeUSU5BDikqjYn2Q/YBLyjqh6YcGlL1n0Let+qej7JPsCdwMVV9d0JlzaSJH8PDIDfr6qzJ13PqJI8Dgyqavm/ZLTMklwLfLuqruyu/vu9qvrZpOvqo7u9zDbghKr64Ti22cwZfVXdweyVPataVT1ZVZu76eeAB5n9hvGqU7Oe72b36V6r8swiyWHAnwJXTroWzUryKuBkZq/uo6p+tdpDvnMq8INxhTw0FPQt6u72eRxw12QrGV033LEF2A5srKrVeiz/CnwI+M2kCxmDAm5Nsqm79chq9XpgBviPbkjtyiT7TrqoMTgP2DDODRr0e6gkrwSuBy6pqmcnXc+oqurXVXUss9+cXptk1Q2rJTkb2F5VmyZdy5icVFXHM3vn2Yu6Yc/VaG/geOBzVXUc8HPgRbdLX0264adzgC+Pc7sG/R6oG8++Hriuqr466XrGofsv9e3AGZOuZQQnAud0Y9tfBN6a5L8mW9Loqmpb974duIHZO9GuRluBrXP+l/gVZoN/NTsT2FxVT41zowb9Hqb7APMq4MGq+tSk6+kjyVSS/bvp3wVOAx6abFVLV1UfrqrDqmoNs/+t/lZV/fWEyxpJkn27D/nphjlOB1bllWpV9WPgiSRHdU2nAqvuooV5zmfMwzawB929sq8kG4BTgIOSbAUuq6qrJlvVSE4E3gXc241tA3ykqm6aYE2jOgS4truK4LeA/66qVX1pYgMOBm7oHguxN/CFqrp5siX18nfAdd2Qx2PAeyZcz8i6X7ynAe8d+7ZbubxSkjScQzeS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXu/wEIdsZ9oMGwTgAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAMoUlEQVR4nO3dX4il9X3H8fcnrsHUJDXW6bK40hUiFilEw2ASDKHVGkwjuhdBlDYsRdhepEVJITW5KYFeJDdJelEKi2u7pcY/1YgiJY1sDGmgNZlVU/+sQSNKdlF30ihqLxo0317MY11mZ52z58+c/U7eLxjmPM/5931Q3vvMb845k6pCktTPu+Y9gCRpPAZckpoy4JLUlAGXpKYMuCQ1tWUjn+yss86qHTt2bORTSlJ7Bw4c+HlVLazev6EB37FjB0tLSxv5lJLUXpLn19rvEookNWXAJakpAy5JTRlwSWrKgEtSUwZckpoa6WWESZ4DXgPeBN6oqsUkZwJ3ADuA54Brqurl2YwpSVrtRM7A/6CqLqyqxWH7JmB/VZ0H7B+2JUkbZJIllKuBfcPlfcDOyceRJI1q1IAX8J0kB5LsHvZtraoXhssvAlvXumOS3UmWkiwtLy9POK6kXyvJ7L8aG/Wt9B+vqsNJfht4IMlTR19ZVZVkzT/tU1V7gD0Ai4uL/vkfSZqSkc7Aq+rw8P0IcA9wMfBSkm0Aw/cjsxpSknSsdQOe5PQk73vrMvBJ4HHgPmDXcLNdwL2zGlKSdKxRllC2AvdkZa1oC/DNqvp2kh8Bdya5HngeuGZ2Y0qSVls34FX1LPChNfb/N3DZLIaSJK3Pd2JKUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWpq5IAnOSXJI0nuH7bPTfJQkmeS3JHk3bMbU5K02omcgd8AHDxq+6vA16vqg8DLwPXTHEyS9M5GCniS7cCngZuH7QCXAncNN9kH7JzFgJKktY16Bv4N4AvAr4bt3wJeqao3hu1DwNlr3THJ7iRLSZaWl5cnGlaS9LZ1A57kSuBIVR0Y5wmqak9VLVbV4sLCwjgPIUlaw5YRbnMJcFWSPwJOA94P/C1wRpItw1n4duDw7MaUJK227hl4VX2xqrZX1Q7gWuC7VfXHwIPAZ4ab7QLundmUkqRjTPI68L8CPp/kGVbWxPdOZyRJ0ihGWUL5f1X1PeB7w+VngYunP5IkaRS+E1OSmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlPrBjzJaUl+mOTHSZ5I8uVh/7lJHkryTJI7krx79uNKkt4yyhn4/wKXVtWHgAuBK5J8FPgq8PWq+iDwMnD97MaUJK22bsBrxevD5qnDVwGXAncN+/cBO2cyoSRpTSOtgSc5JcmjwBHgAeCnwCtV9cZwk0PA2bMZUZK0lpECXlVvVtWFwHbgYuB3R32CJLuTLCVZWl5eHnNMSdJqJ/QqlKp6BXgQ+BhwRpItw1XbgcPHuc+eqlqsqsWFhYWJhpUkvW2UV6EsJDljuPwe4HLgICsh/8xws13AvbMaUpJ0rC3r34RtwL4kp7AS/Dur6v4kTwK3J/kb4BFg7wznlCStsm7Aq+q/gIvW2P8sK+vhkqQ58J2YktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1NS6AU9yTpIHkzyZ5IkkNwz7z0zyQJKnh+8fmP24kqS3jHIG/gbwl1V1AfBR4HNJLgBuAvZX1XnA/mFbkrRB1g14Vb1QVQ8Pl18DDgJnA1cD+4ab7QN2zmpISdKxTmgNPMkO4CLgIWBrVb0wXPUisPU499mdZCnJ0vLy8gSjSpKONnLAk7wXuBu4sapePfq6qiqg1rpfVe2pqsWqWlxYWJhoWEnS20YKeJJTWYn3rVX1rWH3S0m2DddvA47MZkRJ0lpGeRVKgL3Awar62lFX3QfsGi7vAu6d/niSpOPZMsJtLgE+CzyW5NFh35eArwB3JrkeeB64ZjYjSpLWsm7Aq+oHQI5z9WXTHUeSNCrfiSlJTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNbVl3gOMLJn9c1TN/jlgcx2LpLnxDFySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlPrBjzJLUmOJHn8qH1nJnkgydPD9w/MdkxJ0mqjnIH/I3DFqn03Afur6jxg/7AtSdpA6wa8qr4P/GLV7quBfcPlfcDOKc8lSVrHuGvgW6vqheHyi8DW490wye4kS0mWlpeXx3w6SdJqE/8Ss6oKOO4nJ1XVnqparKrFhYWFSZ9OkjQYN+AvJdkGMHw/Mr2RJEmjGDfg9wG7hsu7gHunM44kaVSjvIzwNuA/gPOTHEpyPfAV4PIkTwN/OGxLkjbQun/QoaquO85Vl015FknT4B8M+bXhOzElqSkDLklNGXBJaqrPHzXWyWcj1lrB9VbpODwDl6SmDLgkNWXAJakp18AlcD1fLXkGLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNTRTwJFck+UmSZ5LcNK2hJEnrGzvgSU4B/g74FHABcF2SC6Y1mCTpnU1yBn4x8ExVPVtVvwRuB66ezliSpPVsmeC+ZwM/O2r7EPCR1TdKshvYPWy+nuQnYz7fWcDPx7zvaJKZPvxRNsuxzP44wGM5cZvl/y/YPMcy6XH8zlo7Jwn4SKpqD7Bn0sdJslRVi1MYae42y7FsluMAj+VktVmOZVbHMckSymHgnKO2tw/7JEkbYJKA/wg4L8m5Sd4NXAvcN52xJEnrGXsJpareSPLnwL8BpwC3VNUTU5vsWBMvw5xENsuxbJbjAI/lZLVZjmUmx5GqmsXjSpJmzHdiSlJTBlySmjrpA57kliRHkjw+71kmkeScJA8meTLJE0lumPdM40pyWpIfJvnxcCxfnvdMk0pySpJHktw/71kmkeS5JI8leTTJ0rznGVeSM5LcleSpJAeTfGzeM40jyfnDf4u3vl5NcuPUHv9kXwNP8gngdeCfqur35j3PuJJsA7ZV1cNJ3gccAHZW1ZNzHu2EJQlwelW9nuRU4AfADVX1n3MebWxJPg8sAu+vqivnPc+4kjwHLFbV7N+UNENJ9gH/XlU3D69y+42qemXec01i+PiRw8BHqur5aTzmSX8GXlXfB34x7zkmVVUvVNXDw+XXgIOsvJu1nVrx+rB56vB1cp8JvIMk24FPAzfPexZBkt8EPgHsBaiqX3aP9+Ay4KfTijc0CPhmlGQHcBHw0HwnGd+w5PAocAR4oKraHgvwDeALwK/mPcgUFPCdJAeGj7Ho6FxgGfiHYVnr5iSnz3uoKbgWuG2aD2jAN1iS9wJ3AzdW1avznmdcVfVmVV3IyjtwL07ScnkryZXAkao6MO9ZpuTjVfVhVj4l9HPDEmQ3W4APA39fVRcB/wO0/rjqYRnoKuBfpvm4BnwDDevFdwO3VtW35j3PNAw/2j4IXDHvWcZ0CXDVsHZ8O3Bpkn+e70jjq6rDw/cjwD2sfGpoN4eAQ0f9VHcXK0Hv7FPAw1X10jQf1IBvkOEXf3uBg1X1tXnPM4kkC0nOGC6/B7gceGq+U42nqr5YVduragcrP+J+t6r+ZM5jjSXJ6cMvyBmWHD4JtHv1VlW9CPwsyfnDrsuAdr/sX+U6prx8AhvwaYSTSnIb8PvAWUkOAX9dVXvnO9VYLgE+Czw2rB0DfKmq/nWOM41rG7Bv+K36u4A7q6r1y+82ia3APSvnCmwBvllV357vSGP7C+DWYenhWeBP5zzP2IZ/TC8H/mzqj32yv4xQkrQ2l1AkqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpv4P1sTlIXjORu4AAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "source": [
        "def biased_die_histogram(p1, k):\n",
        "    assert len(p1)==6, \"Il faut que la taille de p soit égale à 6\"\n",
        "    results = 1 + np.random.choice(6, p =p1, size=k) \n",
        "    plt.hist(results, bins  = [1,2,3,4,5,6,7], color = \"red\", width=0.5, align='mid')\n",
        "    plt.show()\n",
        "    \"\"\"\n",
        "    Draw a biased die with probabilities p, k times, and plot the resulting histograms.\n",
        "    There should be six bars.\n",
        "    \n",
        "    Parameters\n",
        "    ----------\n",
        "    p : np array of shape (6,)\n",
        "        The probabilities associated with each outcome of the dice.\n",
        "    k: int\n",
        "        Number of rolls\n",
        "    \"\"\"\n",
        "biased_die_histogram([1/6, 1/6, 1/6, 1/6, 1/6, 1/6], 100)\n",
        "biased_die_histogram([1/10, 1/10, 1/10, 1/10, 1/10, 1/2], 100)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bXcGMli9n5t0"
      },
      "source": [
        "# SVD, Pseudo-Inverse in numpy\n",
        "\n",
        "Generate a random matrix A of size (3, 6) with iid N(0,1) entries, and verify that the singular values squared are the eigenvalues of the matrix $A^TA$.\n",
        "\n",
        "Verify that the pseudo-inverse computed by ``np.linalg.pinv`` is the same as the pseudo-inverse computed from its definition using the SVD."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "U2uF2uiSn5t3",
        "outputId": "666f1ee4-1a88-4b07-fa2b-5855b49b4be1"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mêmes valeurs propores ? :  True\n",
            "Même inverse ? :  True\n"
          ]
        }
      ],
      "source": [
        "def test_singular_values_eigenvalues():\n",
        "    \"\"\"Test that the square of the singular values of X gives the eigenvalue of X^T.dot(X)\n",
        "    \"\"\"\n",
        "    A = np.random.normal(0, 1, size = (3,6))\n",
        "    U,S,V = np.linalg.svd(A)\n",
        "    eig_vals_A_AT = np.sqrt(np.linalg.eigvals(np.matmul(A,A.T)))\n",
        "    eig_vals_A_AT = np.sort(eig_vals_A_AT)[::-1] #Sorting them in descending order\n",
        "    print(\"Mêmes valeurs propores ? : \", np.min(np.isclose(S, eig_vals_A_AT)))\n",
        "\n",
        "test_singular_values_eigenvalues()\n",
        "\n",
        "def test_pinv():\n",
        "    \"\"\"Test that computing the pseudo inverse with pinv or via svd gives the same result\n",
        "    \"\"\"\n",
        "    A = np.random.normal(0, 1, size = (3,6))\n",
        "    U,S,V = np.linalg.svd(A, full_matrices=False)\n",
        "    diagonal = np.zeros((S.shape[0], S.shape[0]))\n",
        "    inverse_svd = np.matmul(V.T, np.matmul(np.diag(1/S), U.T))\n",
        "    print(\"Même inverse ? : \", np.min(np.isclose(inverse_svd, np.linalg.pinv(A))))\n",
        "\n",
        "test_pinv()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IwiVMNLyn5t5"
      },
      "source": [
        "# Boxplots\n",
        "\n",
        "Using matplotlib, for each $\\theta = 1,2,...20$, generate 100 independent copies of a random variable with $N(\\theta, \\theta^{-1})$ distribution. Plot the generated random variables with boxplots, one boxplot for each value of $\\theta$ and the x-axis representing $\\theta$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {
        "id": "YWKSahS9n5t6",
        "outputId": "b2d19f02-2001-47cf-9e68-31dbcd3b2c2e",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 268
        }
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD7CAYAAABzGc+QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAYcklEQVR4nO3dfXBcZ3XH8e/RRpaISIJEXmqiuKYUMkJbCKBmaDEvciCTpEwCNDBoKDEj1R7RRhjS1g5WhzjDmEkEDvWoHdwEqWCIFygkJDAUHCpBqhkI2CEvMgqEl0CUOrGp06SYkS07p3/slVjLu9rdu2/37v4+Mzu7e1f37JF0dfTsc5/nPubuiIhI/DTVOgEREQlHBVxEJKZUwEVEYkoFXEQkplTARURiSgVcRCSm8hZwM7vAzCbN7Mdmtt/MNgbbO8zsHjN7NLhvr3y6IiKywPKNAzezlcBKd7/fzM4A9gFvA94HHHb3m8zseqDd3TdXOmEREUnLW8BP2cHsLuCfg9ub3P1AUOS/4+4XLrfv2Wef7atXrw6bq4hIQ9q3b99v3P2cpdtPKyaIma0GXgXcB5zn7geCl54Ezsu3/+rVq9m7d28xbyki0vDM7FfZthd8EtPMng98Bfiguz+b+Zqnm/FZm/JmtsHM9prZ3kOHDhWRsoiILKegAm5mzaSL9+3ufkew+amg62Shn/xgtn3d/VZ373H3nnPOOeUTgIiIhFTIKBQDxoAZd78l46W7gXXB43XAXeVPT0REcimkD/x1wHuBh83sgWDbFuAm4EtmNgD8CnhXZVIUEZFs8hZwd58CLMfLl5Q3HRERKZRmYoqIxJQKuIhITKmAi4jEVFETeUREpDzSA/xOVuzMeBVwEZFq2noWAH7DmTlfY+szBYVSARcRqSK78dllX29vb+fw1sJiqYCLiFTRQjeJulBERGKq2GKdjUahiIjElAq4iEgRUqkUyWSSRCJBMpkklUrVLBd1oYiIFCiVSjE8PMzY2Bhr1qxhamqKgYEBAPr6+qqeT9Er8pSip6fHtaCDiMRVMplkdHSU3t7exW2Tk5MMDQ0xPT1dsfc1s33u3nPKdhVwEZHCJBIJ5ubmaG5uXtw2Pz9Pa2srJ06cqNj75irg6kIRESlQV1cXK1asqHUai1TARUQKsfUspq8Grs4ygzJ4PX1f2CzKclABFxEpRFCYs03AgfKM6y6WCriISBFqUahz0ThwEZGYUgEXEYkpFXARkZhSARcRiSkVcBGRmFIBFxGJKRVwEalrHR0dmNmyt46OjlqnGYrGgYtIXTv8gRNAjtmTiyp3HZNKUgEXkfqWZwZld3d3Ra8kWEkq4CLSEJqamnJeSTCu1AcuIg2hq6uLqampk7ZNTU3R1dVVo4xKpwIuIg1heHiYgYEBJicnmZ+fZ3JykoGBAYaHh2udWmjqQhGRhrCw5NnQ0BAzMzN0dXWxbdu2miyFVi5akUdEJOJyrcijLhQRkZhSF4qIRF62IYBRui53raiAi0h0BcuU+Q1ZJuLUYAmzqFEXiohE19ZnaL3pKLeceQNsfQa78VnY+gy3nHkDrTcdbejiDTqJKSIRl2sGJTRON4pOYopILLW0tLB9+3bcffG2fft2Wlpaap1azakPXEQibf369WzevBmAwcFBdu7cyebNmxkcHKxxZrWnFriIVEwqlSKZTJJIJEgmk6RSqaJjjI6OMjg4yJYtW2hra2PLli0MDg4yOjpagYzjJW8BN7NxMztoZtMZ27aa2RNm9kBwu6KyaYpI3KRSKYaHhxkdHWVubo7R0VGGh4dDF/G5uTncfTGWFNYC/wxwWZbtn3T3i4LbN8qblojE3bZt2xgbG6O3t5fm5mZ6e3sZGxtj27ZttU6tbuTtA3f3e81sdeVTEZF6MjMzw9q1a2udRl0rpQ/8WjN7KOhiaS9bRiISex0dHTz33HM5X4/zMmZREnYUyqeAjwIe3G8H+rN9oZltADYArFq1KuTbiUic1PMyZlESqoC7+1MLj83sNuDry3ztrcCtkJ7IE+b9RCRmtj5DMplk//79WV9ulAk4lRaqC8XMVmY8fTsQzwXlRKRiZmZmOHbs2EkTcI4dO0ZTk0Yvl0shwwhTwPeAC81s1swGgBEze9jMHgJ6gQ9VOE8RiZl6XMIsavIWcHfvc/eV7t7s7p3uPubu73X3P3H3V7j7le5+oBrJikj1lDoJpx6XMIsaTaUXkVMsTMIZGxtjzZo1TE1NMTAwAFDwEmT1uIRZ1OhqhCJyimQyyejoKL29vYvbJicnGRoaYnpap7yqLdfVCFXAReQUiUSCubk5mpubF7fNz8/T2trKiRMa/ldtuQq4ulBE5BRdXV2sWLHilO3d3d01yEZyUQEXkd8Llimbvhq4OttEnMfTX9PgK+FEhQq4iCyyG5/N+zXt7e0c3lr5XCQ/FXARWZR5TkwrwUefCriIZKViHX2a0yoiElMq4CIiMaUCLiISUyrgIiIxpZOYInVGo0cah1rgInWio6Mja/GGdFHXMmb1Ry1wkTqhZcwajwq4SJ3QLMrGowIuUifcnUQikXU1+KamJl1FsA6pD1ykjnR1dTExMXHSOpQTExNaxqxOqYCL1BEtY9ZY1IUiUke0jFlj0Yo8IiIRl2tFHnWhiIjElLpQRCIi1yQczaKUXNQCF6mx5WZQAppBKTmpBS5SY5pBKWGpgIvUWrBAcK5WeHd3N9PT09XMSGJCBVwkIpqampibm6O5uXlx2/z8PK2trTXMSqJMfeAiEdHV1cXU1NRJ26ampjSLUnJSAReJCM2ilGKpC0UkIjSLUoqlmZgiIhGXayamWuAiZaBlzKQW1AcuUga7d++mu7sbSA/72717d40zkkagLhSRUmw9q4CveabyeUhdUxeKSAXkW8ZMS5hJJakLRaQE+T7BHj58uEqZSCNSARcpgz179py0jNmePXtqnZI0ABVwkRJ1dnZyzTXXnDQB55prrqGzs7PWqUmdUwEXKdHIyAgnTpygv7+flpYW+vv7OXHiBCMjI7VOTepc3gJuZuNmdtDMpjO2dZjZPWb2aHDfXtk0RaKrr6+PHTt20NbWhpnR1tbGjh07NINSKi7vMEIzewPwW2CXuyeDbSPAYXe/ycyuB9rdfXO+N9MwQokarYIjcRB6TUx3vxdYeir9KuCzwePPAm8rOUORGlg46Zj5WMVb4iJsH/h57n4gePwkcF6Z8hGpioVlzBZuwEnPtYyZxEHJE3nc3c0sZ5PFzDYAGwBWrVpV6tuJlIWWMZN6ELaAP2VmK939gJmtBA7m+kJ3vxW4FdJ94CHfT6Ss8s2gBM2ilOgL24VyN7AueLwOuKs86YhUR2Z/d7Y+cHfXLEqJvEKGEaaA7wEXmtmsmQ0ANwFvMbNHgTcHz0ViJ1cfuEgc5O1Ccfdcg1kvKXMuIlWnEScSZ5qJKSISUyrgIiIxpeuBS2xpGTNpdCrgElsLxdrMVLilIamAS7xkWcLMbzjz1O1axkwagAq4xIom4Ij8ngq4xEq2rhJ1oUijUgGX2Mo8ibnwWIVcGokKuMSWirU0Oo0DFxGJKbXApeo0flukPFTApeo0flukPNSFIiISUyrgUjVaxkykvNSFIlWjZcxEyksFXKpGsyhFyksFXKpm6QlLncQUKY0KuFSdZlCKlIcKuFSdirVIeWgUiohITKmAi4jElLpQpGDZpsCDukREakUFXAqmKfAi0aIuFMlLMyhFokktcMlLMyhFokkFXPLSDEqRaFIXSoNIpVIkk0kSiQTJZJJUKlXwvu5+0i3btsOHD1cqdRHJQS3wBpBKpdi4cSNtbW24O0eOHGHjxo0A9PX1FRwn2wxK0CgUkVpRC7wBbNq0iUQiwfj4OEePHmV8fJxEIsGmTZuKirO01Z3ZIheR6lMBbwCzs7Ps2rWL3t5empub6e3tZdeuXczOztY6NREpgbpQGsSll15a6xREpMzUAm8AnZ2drFy5komJCQAmJiZYuXIlnZ2dNc5MREqhAt4ARkZGOH78OP39/QD09/dz/PhxRkZGapyZiJTCqnkSqqenx/fu3Vu192t0HR0dPP3008t+TXt7u4YAikScme1z956l29UHXsc0g1KkvqmAR1wqlWLbtm3MzMzQ1dXF8PBwwWO3NYNSpL6pgEdYKpVieHiYsbEx1qxZw9TUFAMDA0BhE3CydY/pSoIi9UN94BGWTCYZHR2lt7d3cdvk5CRDQ0NMT08XFSvbtbxVyEXiQX3gMTQzM8PatWtP2d7UVPzgIRVrkfqjYYQR1tXVxcTExEkXkJqYmKCrq6vGmYlIFJRUwM3sMTN72MweMDP1jZTZ8PAwAwMDTE5OAunuk4GBAYaHh2ucmYhEQTm6UHrd/TdliCNLLJyoHBoaWrzftm1bUVcQFJH6VdJJTDN7DOgptIDrJGbhNAlHRBZU6iSmA3vMzIF/dfdbs7zxBmADwKpVq0p8u8ahSTgikk+pLfDz3f0JMzsXuAcYcvd7c329WuCFyzbsbym1wEUaQ0Va4O7+RHB/0MzuBC4GchZwKdzSf6yagCMiS4UehWJmbWZ2xsJj4FKguNklIiISWikt8POAO4OP+qcBu939m2XJqg7k6gIpthWdbR1KtcRFBEoo4O7+C+CVZcylrmQW2VK6P1SsRSQXzcQUEYkpFXARkZhSAS+zjo4OzOykG3DS846OjhpnKSL1QFcjzGJoaIjbbruNo0eP0tLSwvr16xkdHS1oX03AEZFqUQFfYmhoiJ07d3LzzTczODjIzp072bx5M0BBRdxufDbviUczw7eWI1sRaWRa0GGJ1tZWPvaxj3HdddctbrvlllvYsmULc3NzeffXDEoRKbdcMzFVwJdYrgCH/VlpFqWIlCJXAddJzCVaWlrYvn37SYsobN++nZaWlhpnJiJyMvWBL7F+/frFPm9Id59s3ryZwcHBGmYlInIqFfAlFk5UbtmyZfF+cHCw4FEoC5Z2xWgavIiUW131gZfr+iOZ8VRwRaTWGqIP3N1P6rvOfF6IpZNwQBNwRCS6ItOFkq31XO3W79NPP73sexYyRFBEpFoiU8AXCmctuy38hjNh61nLvy4iEhGRKeBRkG8WpWZQikiUqIAvsVw3SXt7exUzERFZngp4Bq1DKSJxUlejUEREGokKuIhITKmAi4jEVF0U8HwTcIqdhJMtjohI1NTFScx8E3CguEk4OnEpInFQFy1wEZFGVBct8HwzKBe/RkSkjtRFAdc6lCLSiNSFIiISUyrgIiIxVdMCvnT4XynX4F4aZ+lN1zERkXpT0z7wcg3/0zVMRKQRqQtFRCSmVMBFRGKqpl0oGr8tIhJeTQu4xm+LiISnLhQRkZhSARcRiSkVcBGRmKqLa6EsyBwznvlYY8JFpB7VVQFXoRaRRlLzAp5vpqWmwIuIZFdSH7iZXWZmPzGzn5nZ9cXu7+6n3JZuP3z4cCkpiojUrdAF3MwSwL8AlwMvB/rM7OXlSkxERJZXSgv8YuBn7v4Ldz8GfAG4qjxpiYhIPqUU8POBxzOezwbbTmJmG8xsr5ntPXToUAlvJyIimSo+Dtzdb3X3HnfvOeeccyr9diIiDaOUAv4EcEHG885gm4iIVEEpBfyHwEvN7MVmtgJ4N3B3edISEZF8Qo8Dd/fjZnYt8C0gAYy7+/6yZSYiIssqaSKPu38D+EaZchERkSLoYlYiIjGlAi4iElMq4CIiMVXzi1ktyHYpWF1dUEQkt8i0wHfv3k13dzdNTU10d3eze/fuWqckIhJpkSjgqVSKjRs3cuTIEdydI0eOsHHjRlKpVK1TExGJrEgU8E2bNpFIJBgfH+fo0aOMj4+TSCTYtGlTrVMTEYmsSBTw2dlZdu3aRW9vL83NzfT29rJr1y5mZ2drnZqISGRFooCLiEjxIlHAOzs7WbduHZOTk8zPzzM5Ocm6devo7OysdWoiIpEViQI+MjLC8ePH6e/vp7W1lf7+fo4fP87IyEitUxMRiaxIFPC+vj527NhBW1sbAG1tbezYsYO+vr4aZyYiEl1WzckyPT09vnfv3qq9n4hIPTCzfe7es3R7JFrgIiJSPBVwEZGYUgEXEYkpFXARkZhSARcRiamqjkIxs0PAr/J82dnAb0p8q3LEiFoc5VLZOFHKpVxxlEtl41Qzlz9093NO2erukboBe6MQI2pxlIu+J+USrThRyEVdKCIiMaUCLiISU1Es4LdGJEbU4iiXysaJUi7liqNcKhun5rlU9SSmiIiUTxRb4CIiUoDIFHAzGzezg2Y2XUKMC8xs0sx+bGb7zWxjyDitZvYDM3swiHNjCTklzOxHZvb1EmI8ZmYPm9kDZhbqamBm9gIz+7KZPWJmM2b2ZyFiXBjksHB71sw+GDKfDwU/22kzS5lZa4gYG4P99xeTR7Zjzcw6zOweM3s0uG8PGeedQT7PmdkpFx8qMMbHg9/TQ2Z2p5m9IGScjwYxHjCzPWb2ojBxMl77OzNzMzs7RC5bzeyJjGPnirC5mNlQ8PPZb2bLXnM6Ry5fzMjjMTN7IEwuZnaRmX1/4e/SzC4OGeeVZva94G/8a2Z2Zr44i8oxDKZMQ2neALwamC4hxkrg1cHjM4CfAi8PEceA5wePm4H7gNeGzOk6YDfw9RK+r8eAs0v8+X4W+Ovg8QrgBSXGSwBPkh6fWuy+5wO/BJ4XPP8S8L4iYySBaeB04DTg28Afhz3WgBHg+uDx9cDNIeN0ARcC3wF6Qsa4FDgteHxzCbmcmfH4A8DOMHGC7RcA3yI9j2PZYzFHLluBvy/yd5wtTm/wu24Jnp8b5vvJeH078JGQuewBLg8eXwF8J2ScHwJvDB73Ax8t9GcUmRa4u98LHC4xxgF3vz94/H/ADOliUWwcd/ffBk+bg1vRJwvMrBP4C+DTxe5bTmZ2FukDZwzA3Y+5+/+WGPYS4Ofunm9iVi6nAc8zs9NIF+H/LnL/LuA+d/+dux8Hvgu8o5AdcxxrV5H+J0dw/7Ywcdx9xt1/Ukgey8TYE3xPAN8H8i5NlSPOsxlP2yjgGF7m7/CTwKYSYxQlR5z3Aze5+9Hgaw6GzcXMDHgXkAqZiwMLreWzKOAYzhHnZcC9weN7gL/MF2dBZAp4uZnZauBVpFvPYfZPBB+tDgL3uHuYOP9E+qB/LkwOGRzYY2b7zGxDiP1fDBwC/i3ozvm0mbWVmNO7KeDAz8bdnwA+AfwaOAA84+57igwzDbzezF5oZqeTbgFdECafwHnufiB4/CRwXgmxyqkf+I+wO5vZNjN7HHgP8JGQMa4CnnD3B8PmEbg26NIZL6SLKoeXkf6932dm3zWzPy0hn9cDT7n7oyH3/yDw8eDn+wngwyHj7CfdgAB4J0Ucx3VZwM3s+cBXgA8uaYUUzN1PuPtFpFs/F5tZssgc3gocdPd9Yd5/iTXu/mrgcuBvzewNRe5/GumPbZ9y91cBR0h3E4RiZiuAK4F/D7l/O+kD9sXAi4A2M/urYmK4+wzp7oU9wDeBB4ATYfLJEtsJ8Ymr3MxsGDgO3B42hrsPu/sFQYxrQ+RwOrCFkMU/w6eAlwAXkf6nvT1knNOADuC1wD8AXwpa0mH0EbIREng/8KHg5/shgk+4IfQDf2Nm+0h3/R4rdMe6K+Bm1ky6eN/u7neUGi/oapgELity19cBV5rZY8AXgLVm9vmQOTwR3B8E7gTynixZYhaYzfgU8WXSBT2sy4H73f2pkPu/Gfilux9y93ngDuDPiw3i7mPu/hp3fwPwNOlzHmE9ZWYrAYL7ZT+aV5qZvQ94K/Ce4B9KqW6niI/mGV5C+h/tg8Gx3Ancb2Z/UEwQd38qaBQ9B9xG8cfwglngjqCb8wekP90ue1I1m6Dr7h3AF0PmAbCO9LEL6cZMqO/J3R9x90vd/TWk/6H8vNB966qAB/+Jx4AZd7+lhDjnLJz5N7PnAW8BHikmhrt/2N073X016e6GCXcvqpUZvH+bmZ2x8Jj0Ca6iRuq4+5PA42Z2YbDpEuDHxeaSodSWy6+B15rZ6cHv7BLS5yuKYmbnBverSP8x7i4hp7tJ/0ES3N9VQqySmNllpLvernT335UQ56UZT6+iyGMYwN0fdvdz3X11cCzPkh4o8GSRuazMePp2ijyGM3yV9IlMzOxlpE/Ih7mg1JuBR9x9NmQekO7zfmPweC0Qqism4zhuAv4R2FnwzoWe7az0jXRBOADMkz5IBkLEWEP6o+9DpD9SPwBcESLOK4AfBXGmKeAsdZ54byLkKBTgj4AHg9t+YDhknIuAvcH39FWgPWScNuB/gLNK/JncSLqgTAOfIxhVUGSM/yL9j+hB4JJSjjXghcB/kv4j/DbQETLO24PHR4GngG+FiPEz4PGMY7iQ0SPZ4nwl+Pk+BHwNOD9MnCWvP0b+USjZcvkc8HCQy93AypDf0wrg88H3dT+wNsz3A3wGGCzxmFkD7AuOv/uA14SMs5H0p8efAjcRTLAs5KaZmCIiMVVXXSgiIo1EBVxEJKZUwEVEYkoFXEQkplTARURiSgVcRCSmVMBFRGJKBVxEJKb+H6izAO0vwMwtAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "source": [
        "def plot_boxplot_normal_increasing_mean_and_variance():\n",
        "    \"\"\"Using matplotlib, for each theta = 1,2,...,20, \n",
        "    generate 100 independent copies of a random variable \n",
        "    with N(theta, theta^{-1}) distribution. \n",
        "    Plot the generated random variables with boxplots, \n",
        "    one boxplot for each value of theta and \n",
        "    the x-axis representing theta.\n",
        "    \"\"\"\n",
        "    data = [list() for i in range(1, 20)]\n",
        "    for theta in range(1,20):\n",
        "      for i in range(100):\n",
        "        N = np.random.normal(theta, 1/theta)\n",
        "        data[theta-1].append(N)\n",
        "    plt.boxplot(data)\n",
        "\n",
        "plot_boxplot_normal_increasing_mean_and_variance()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AxGpJ1_xn5t9"
      },
      "source": [
        "# Quantile quantile plots\n",
        "\n",
        "In order to check visually that a collection of random variables follow a certain distribution, quantile quantile plots are convenient."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "id": "vWum8p6Rn5t-"
      },
      "outputs": [],
      "source": [
        "n = 10\n",
        "X = np.random.randint(1000, size=n) # generate n numbers\n",
        "X_sort = np.sort(X) # sort them\n",
        "scores  = np.arange(1, n + 1) / (n+1) # their corresponding scores (many other methods are possible)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "id": "Sh3Seghun5uA",
        "outputId": "5f84ea75-c681-42a8-f3a4-e6a1f10a05c7",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 827
        }
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df5xcdX3v8dd7ZrIb8juQ8CsEEiBYoraAMdZHlfahqIG2RCtWqLRYf3Bt5dqWaosXL/qgt63Va39KVay5/mgRf7ZNa7wUq2i9imZBBAIiIQSS8GtJwo9kk92dmc/945zZnN3M7Jzd7O7sWd7Px2Mfe+Z7fsx3z0ze883nnDlHEYGZmc1cpU53wMzMJpeD3sxshnPQm5nNcA56M7MZzkFvZjbDOejNzGY4B70dRtIHJP1jh577Fklvyzz+X5KelPTYOLb1C5Lul7RP0msntqfjJ2m7pPMmadsh6fTJ2Ha6/ZdLum+ytm+Tw0H/HJQGX+OnLulA5vGbOt2/BkknA38IrI6I48exiWuBj0bEvIj4lxbP8RuSetK//VFJX5f0siPp90wWEf8VEc9rPJ7MDy2bOA7656A0+OZFxDzgYeBXM23/1On+ZZwM7I6IJ8a5/inAllYzJV0J/DXwZ8Bx6fP9PbB+nM9nNi056K2VLkmflfSspC2S1jRmSDpR0lck9Up6UNK7Wm1E0qclfVzSzem2vi3plMz8V0n6iaSnJX0UUNp+HnAzcGI62v50i+2/XdJWSXskbZR0Ytr+AHAq8G/p+t0j1ltIMuJ/Z0R8NSL2R8RgRPxbRLwnXWatpO9Leiod7X9UUlc6b0VaJqlktjmy7PR2Sfemf/c9ks7JdOEsSXemf/cXJM1O17lb0q9mtjErLV2d3eLvf0/at0ckvWXEvG5J/1vSw5IeT1+Ho9J5vyRpp6Q/lPREuo3fzqx7QdrnZyXtkvTu7Hrp9OdIPhwb+/iPJH1N0n8f0Y87Jb2uWf9tikSEf57DP8B24LwRbR8ADgIXAGXgz4Fb03kl4DbgGqCLJEy3Aa9psf1PA88C5wLdwN8A303nLUnnXQTMAv4AqAJvS+f/ErBzlL6/AngSOCfd9t8B3xntb8vMW5c+V2WU7b8I+HmgAqwA7gV+P523Aojs+sAtmb6/AdgFvJjkw+t04JRMv34InAgcnW73Hem8PwK+kNnmeuCuUf6Gx4EXAHOBG9I+nZ7O/ytgY/oc84F/A/48s2+rJB92s9LXug9YnM5/FHh5Or0YOKfZazJyHwO/Dvwg8/jngN1AV6ff68/lH4/orZXvRsSmiKgBnyP5BwtJcC2NiGsjYiAitgGfBC4eZVtfi4jvREQ/cDXwUknLScJlS0R8OSIGScooYzno+iZgQ0Tcnm77vem2V+RY9xjgyYiotlogIm6LiFsjohoR24FPAL+Ys29vAz4UEZsjsTUiHsrM/9uIeCQi9pAE8Flp+z8CF0hakD7+TZL938yvA/8nIu6OiP0kH9AASBJwOfAHEbEnIp4lKVFlX6dB4NpI/iezCdgHPC8zb7WkBRGxNyJuz/l3bwTOkLQq0/8vRMRAzvVtEjjorZVs4PYBs9MyxSkk5ZSnGj/A/yCpcbeyozEREfuAPSSj2RNHzIvs4xxOBIbCM932bmBZjnV3A0uypZeRJJ0h6d8lPSbpGZKgXJKzb8uBB0aZP3L/zgOIiEeA/we8XtIi4Hyg1XGTYfuPzL4AlgJzgNsyr9P/Tdsbdo/4oBvqB/B6kg/ih9Jy20tH+VuGRMRB4AvApZJKwCW0/qCyKeKgt7HaATwYEYsyP/Mj4oJR1lnemJA0j6SU8AhJeSA7T9nHOTxC8sHTWH8uyUh9V451vw/0A6Oddvkx4CfAqohYQPKBpnTe/vT3nMzy2TODdgCn5ehHM58BLiUp/3w/Ilr9PcP2H0m9vOFJ4ADw/MzrtDCSA/Btpf8TWQ8cC/wL8MVWi7bo/5uAVwJ9EfH9PM9pk8dBb2P1Q+BZSX8s6ShJZUkvkPTiUda5QNLL0gOZf0JS798BfA14vqRfS0fW72J4WLbzeeC3JZ2VHmz9M5L68PZ2K0bE0yTHGa6T9FpJc9IDn+dL+lC62HzgGWCfpJ8Bfiezfi/JB8ql6T54C8OD/R+Ad0t6kRKnZw9Ct/EvJMcdfg/47CjLfRF4s6TVkuYA78/0r05SUvsrSccCSFom6TXtnlxSl6Q3SVqYltSeAeotFn+c5DjNkDTY68BH8Gh+WnDQ25ikNftfIakpP0gycvwHYOEoq91AEkJ7SA5wXppu60mSUesHSUopq0jKFnn78g3gfwJfIRndnsboxwpGrv8R4ErgfUAvySj8CpKgBXg38BskB4w/SVKSyHo78J60788HvpfZ9peAPyX5259Nt3l0zn4dSP+mlcBXR1nu6yTHNb4JbE1/Z/1x2n5rWnr6Bodq8O38JrA9Xe8dJCP0Zv4ceF9aHnp3pv2zwAtJjjlYhykpi5pNDiWnRe6MiPd1ui9FIuka4IyIuLTTfRkPSb8FXB4R/vLZNNDyQJSZdYako4G3koyqCyctI/0uyZfPbBpw6cZsGpH0dpIS0tcj4jud7s9YpccAeklq9zd0uDuWcunGzGyG84jezGyGm3Y1+iVLlsSKFSs63Q0zs0K57bbbnoyIpc3mTbugX7FiBT09PZ3uhplZoUh6qNU8l27MzGY4B72Z2QznoDczm+Ec9GZmM5yD3sxshnPQm5nNcA56M7MZzkFvltPm7Xv4yWPPdLobZmM27b4wZTZdveHjyY2Stn/wlzvcE7Ox8YjezGyGc9Cbmc1wDnozsxnOQW+WQ73u+zZYcTnozXLor9Y73QWzcXPQm+VwYLDW6S6YjZuD3iwHB70VWa6gl7RO0n2Stkq6qsn8d0i6S9Idkr4raXVm3nvT9e5LbxxsVjiDLt1YgbUNekll4DrgfGA1cEk2yFM3RMQLI+Is4EPAX6brrgYuBp4PrAP+Pt2eWaFUfTDWCizPiH4tsDUitkXEAHAjsD67QERkvxc+F2j8q1gP3BgR/RHxILA13Z5ZodTjUNBHOPStWPJcAmEZsCPzeCfwkpELSXoncCXQBbwis+6tI9ZdNq6emnVQtXYo3AdrQVdFHeyN2dhM2MHYiLguIk4D/hh431jWlXS5pB5JPb29vRPVJbMJU6tng971eiuWPEG/C1ieeXxS2tbKjcBrx7JuRFwfEWsiYs3SpUtzdMlsatUy5Zrs6N6sCPIE/WZglaSVkrpIDq5uzC4gaVXm4S8D96fTG4GLJXVLWgmsAn545N02m1q1+qFR/IBH9FYwbWv0EVGVdAVwE1AGNkTEFknXAj0RsRG4QtJ5wCCwF7gsXXeLpC8C9wBV4J0R4ROSrXCyo/hq3UFvxZLrevQRsQnYNKLtmsz0742y7p8CfzreDppNB9nSzWDVpRsrFn8z1iyH7MHYmk+vtIJx0JvlkP3CVM2lGysYB71ZDtnLFPtbslY0DnqzHIaP6B30ViwOerMcag56KzAHvVkONZdurMAc9GY5eERvReagN8vBNXorMge9WQ51B70VmIPeLIeqa/RWYA56sxyyX5LyF6asaBz0ZjkMPxjbwY6YjYOD3iwHXwLBisxBb5aDz6O3InPQm+WQvWKlz7qxonHQm+VQqznorbgc9GY5+PRKKzIHvVkOdZdurMAc9GY5eERvReagN8shO4qvO+itYBz0Zjn49EorMge9WQ61elApKZ32F6asWBz0ZjlU63W6KqV02iN6K5ZcQS9pnaT7JG2VdFWT+VdKukfSnZL+U9IpmXk1SXekPxsnsvNmU6VWZyjoXaO3oqm0W0BSGbgOeBWwE9gsaWNE3JNZ7EfAmojok/Q7wIeAN6bzDkTEWRPcb7MpVavX6Sp7RG/FlGdEvxbYGhHbImIAuBFYn10gIr4VEX3pw1uBkya2m2adVavDrHIJyefRW/HkCfplwI7M451pWytvBb6eeTxbUo+kWyW9ttkKki5Pl+np7e3N0SWzqVWPoFSCSkke0VvhtC3djIWkS4E1wC9mmk+JiF2STgW+KemuiHggu15EXA9cD7BmzRr/K7JpJznrpkS5JNforXDyjOh3Acszj09K24aRdB5wNXBhRPQ32iNiV/p7G3ALcPYR9NesI2r1oCSolEoe0Vvh5An6zcAqSSsldQEXA8POnpF0NvAJkpB/ItO+WFJ3Or0E+AUgexDXrBBq9aBcEiXX6K2A2pZuIqIq6QrgJqAMbIiILZKuBXoiYiPwYWAe8CVJAA9HxIXAmcAnJNVJPlQ+OOJsHbNCqEVQkqiUSw56K5xcNfqI2ARsGtF2TWb6vBbrfQ944ZF00Gw6qNWDSlmUfTDWCsjfjDXLoVYPyhKVknwJBCscB71ZDsnplaIkj+iteBz0Zjk0LmpWKfv0SiseB71ZDtV6cjDWNXorIge9WQ719PTKpEbvoLdicdCb5VCLxnn0HtFb8TjozXJofGHKNXorIge9WQ6N0yvLvgSCFZCD3iyHWj05vdI1eisiB71ZDvVITq8sS1T9hSkrGAe9WQ7VdESfXKa4070xGxsHvVkO9cYlEMoe0VvxOOjNcqg1Sjeu0VsBOejNcqjV0tKNz6O3AnLQm+VQi8bplR7RW/E46M1yqNVJTq8sO+iteBz0ZjkMnV5Z8h2mrHgc9GY5VGt1yiVRFq7RW+E46M1yqAfpZYo9orficdCb5dC4Z6wvgWBF5KA3y6GW3nik5BuPWAE56M1ySK5Hj28OboWUK+glrZN0n6Stkq5qMv9KSfdIulPSf0o6JTPvMkn3pz+XTWTnzabKocsUu3RjxdM26CWVgeuA84HVwCWSVo9Y7EfAmoj4WeDLwIfSdY8G3g+8BFgLvF/S4onrvtnka9xopFwquUZvhZRnRL8W2BoR2yJiALgRWJ9dICK+FRF96cNbgZPS6dcAN0fEnojYC9wMrJuYrptNjepQ0OObg1sh5Qn6ZcCOzOOdaVsrbwW+PpZ1JV0uqUdST29vb44umU2deiTBXvJFzaygJvRgrKRLgTXAh8eyXkRcHxFrImLN0qVLJ7JLZkesEeyVxh2mwkFvxZIn6HcByzOPT0rbhpF0HnA1cGFE9I9lXbPprFGqaXxhKgLfINwKJU/QbwZWSVopqQu4GNiYXUDS2cAnSEL+icysm4BXS1qcHoR9ddpmVhiHDsaKcvovxnV6K5JKuwUioirpCpKALgMbImKLpGuBnojYSFKqmQd8SRLAwxFxYUTskfQnJB8WANdGxJ5J+UvMJkmjVFNOL2oGuE5vhdI26AEiYhOwaUTbNZnp80ZZdwOwYbwdNOu07Ii+UhKA6/RWKP5mrFkbQ6dXpl+YguSOU2ZF4aA3a6NRpmmcXgn4BuFWKA56szYa59FXMkHvGr0ViYPerI2qa/RWcA56szbqw86jT0s3rtFbgTjozdoYfnqlSzdWPA56szZq9cOD3l+YsiJx0Ju1UcucXllJvzBVd43eCsRBb9ZG0xG9a/RWIA56szbqrtFbwTnozdpojN59eqUVlYPerI1GqGdPr/QNwq1Icl3UzOy5rJHp5ZKIcI3eisdBb9ZG47o25ZKoh2v0VjwOerM2sgdjGyN61+itSBz0Zm3UGqUbifAXpqyAHPRmbWTPo2+cv+Dr0VuROOjN2sgGfZBMe0RvReKgN2vj0EXNoDGi9yUQrEgc9GZt1IbOuimBR/RWQA56szYa58xXSkJpm78wZUXioDdrozF6r5R1qM0HY61Acl0CQdI6SfdJ2irpqibzz5V0u6SqpItGzKtJuiP92ThRHTebKtX0/MpKqTR0CQTX6K1I2o7oJZWB64BXATuBzZI2RsQ9mcUeBt4MvLvJJg5ExFkT0Fezjhga0WcvgeAavRVIntLNWmBrRGwDkHQjsB4YCvqI2J7Oc+HSZpyhGn3Zl0CwYspTulkG7Mg83pm25TVbUo+kWyW9ttkCki5Pl+np7e0dw6bNJl9j9D6rXPKNR6yQpuIyxadExBrgN4C/lnTayAUi4vqIWBMRa5YuXToFXTLLr1Gjz954xDV6K5I8Qb8LWJ55fFLalktE7Ep/bwNuAc4eQ//MOm4wU6Nv3DPWNXorkjxBvxlYJWmlpC7gYiDX2TOSFkvqTqeXAL9AprZvVgS1ej05h16+laAVU9ugj4gqcAVwE3Av8MWI2CLpWkkXAkh6saSdwBuAT0jakq5+JtAj6cfAt4APjjhbx2zaq9ZiKOBdo7ciyvWFqYjYBGwa0XZNZnozSUln5HrfA154hH0066jBWjArudANac77evRWKL5nrFkbtXp9aCQvJTcI9yUQrEgc9GZtDNaDWZnLH5RK8sFYKxQHvVkbtVoMnW0Dydk3vvGIFYmD3qyNwUzpBpIDsq7RW5E46M3aqNaGl26SGr2D3orDQW/WRq0eh43oXaO3InHQm7UxWKsPnV4JaenGNXorEAe9WRu1egy76UilVHKN3grFQW/WxmA90vvFJsqu0VvBOOjN2qjW6sxyjd4KzEFv1kZ1ROmm7G/GWsE46M3aqNbqh39hyiN6KxAHvVkbzUf0DnorDge9WRvVWlBxjd4KzEFv1sZArU5XxWfdWHE56M3aGKzV6Sq7Rm/F5aA3a2Ogevg3Y126sSJx0Ju1MejSjRWcg96sjcNH9CWP6K1QHPRmbQzU6nRXhtfo6w56KxAHvVkb2ZuDQxL0gzV/M9aKw0FvNopaPajVY1iNvntWmYGqg96KI1fQS1on6T5JWyVd1WT+uZJul1SVdNGIeZdJuj/9uWyiOm42FRoj9+yIvqtcot9BbwXSNugllYHrgPOB1cAlklaPWOxh4M3ADSPWPRp4P/ASYC3wfkmLj7zbZlOjEejDR/Ql+qu1TnXJbMzyjOjXAlsjYltEDAA3AuuzC0TE9oi4Exg5zHkNcHNE7ImIvcDNwLoJ6LfZlGiM6Lsy17rprnhEb8WSJ+iXATsyj3embXnkWlfS5ZJ6JPX09vbm3LTZ5BtoNqKvlB30VijT4mBsRFwfEWsiYs3SpUs73R2zIU1r9JUSA9U64dsJWkHkCfpdwPLM45PStjyOZF2zjms+ok+mPaq3osgT9JuBVZJWSuoCLgY25tz+TcCrJS1OD8K+Om0zK4SBJiP6RtAP+Fx6K4i2QR8RVeAKkoC+F/hiRGyRdK2kCwEkvVjSTuANwCckbUnX3QP8CcmHxWbg2rTNrBCajuhnlQHoH3TQWzFU8iwUEZuATSParslMbyYpyzRbdwOw4Qj6aNYxg7WkDp+9THF3uVG68SmWVgzT4mCs2XTVfETvGr0Vi4PebBQHB5NR++xKeahtqEbvoLeCcNCbjeJgWp6ZPWv46ZXgEb0Vh4PebBQHBhpBnx3RNw7GukZvxeCgNxvFwXTUPjzoPaK3YnHQm42iMWrPlm4aI3rX6K0oHPRmo2hWunGN3orGQW82ioPVGuWSmn4z1ufRW1E46M1GcXCwzlGZ0Tz4PHorHge92SgODtaG1efh0Dn1B33WjRWEg95sFAcGa0MHXxvmdCeP+wYc9FYMDnqzUfQP1g8b0XeVS1RKYn9/tUO9MhsbB73ZKA4O1jiqa/iIXhJzusoe0VthOOjNRtE3UDvsYCzA3O6KR/RWGA56s1Hs668yf/asw9o9orcicdCbjWJff5W53YfftmFed4X9Ax7RWzE46M1Gsa+/yrwmQT+nq0Jfv0f0VgwOerNR7DtYZV53sxp92SN6KwwHvVkLtXpwYLDGvO5mNXofjLXicNCbtbAvDfJ5sw8v3SQjepdurBgc9GYtDAV9s9JNV4U+j+itIBz0Zi3sHwr6JqWb7gp9gzXq9ZjqbpmNWa6gl7RO0n2Stkq6qsn8bklfSOf/QNKKtH2FpAOS7kh/Pj6x3TebPM8eTIJ+bpMR/fzuChGwzwdkrQAOLz6OIKkMXAe8CtgJbJa0MSLuySz2VmBvRJwu6WLgL4A3pvMeiIizJrjfZpOuUbqZ36RGv2hOMsp/um+QBU2+UGU2neQZ0a8FtkbEtogYAG4E1o9YZj3wmXT6y8ArJWniumk29fYdbF26WTSnC4C9fQNT2iez8cgT9MuAHZnHO9O2pstERBV4GjgmnbdS0o8kfVvSy5s9gaTLJfVI6unt7R3TH2A2WfakIb547uFBvzgd0T/VNzilfTIbj8k+GPsocHJEnA1cCdwgacHIhSLi+ohYExFrli5dOsldMstnz7406NPRe1ajdOMRvRVBnqDfBSzPPD4pbWu6jKQKsBDYHRH9EbEbICJuAx4AzjjSTptNhb19AyyYXRl2v9iGRunGI3orgjxBvxlYJWmlpC7gYmDjiGU2Apel0xcB34yIkLQ0PZiLpFOBVcC2iem62eTavX+Ao+cePpoHWHSUSzdWHG3PuomIqqQrgJuAMrAhIrZIuhboiYiNwKeAz0naCuwh+TAAOBe4VtIgUAfeERF7JuMPMZtoe0cJ+kq5xPzuiks3Vghtgx4gIjYBm0a0XZOZPgi8ocl6XwG+coR9NOuI3fsHWLZodsv5i+d2sWe/g96mP38z1qyF0Ub0AMct6ObxZw5OYY/MxsdBb9ZErR48ua+fJfO6Wy5z3ILZDnorBAe9WRO9z/ZTrQcnLjqq5TInLJzNo08fJMLXu7HpzUFv1sSupw4AsGyUoD9uwWz6q3WePuAzb2x6c9CbNfFIGvSjjeiPX5gcqH3M5Rub5hz0Zk0cCvrWZ900RvsP7+6bkj6ZjZeD3qyJh/f0sWB2hfmjXJny1KXzANj25P6p6pbZuDjozZrY+sQ+Vh03f9RlFh41i6Xzu3ngiX1T1Cuz8XHQmzVx/xP7OOO4eW2XO33pPB7oddDb9OagNxth975+9uwf4PRjRx/RA5x27Fy2PrHPp1jatOagNxvhzl1PA3Dm8e2D/swTFvDMwSoP7/EBWZu+HPRmI/zoob2UBD+3fFHbZV90ymIAerbvnexumY2bg95shJ6H9vK84xcwt7v9Nf/OOHY+82dX6HnIQW/Tl4PeLOOZg4Ns3r6Hc1ctybV8qSTWrjia727tdZ3epi0HvVnGt+/rZbAWnLf6uNzrvOb5x7NjzwG2PPLMJPbMbPwc9GYZX719J0vnd3POyYtzr/Oq1cdRLol/vWPkHTbNpgcHvVlq+5P7ueWnvVyy9mTKJeVeb/HcLta94Hhu3LyDff3VSeyh2fg46M1SH77pPmZXylz6kpPHvO7lLz+VZw9W+dgtWyehZ2ZHxkFvBnztzkf52l2P8t9+8VSOXdD6Qmat/NzyRfzaOcu4/jvbuGvn05PQQ7Pxc9Dbc96t23bzni//mLOWL+J3f+n0cW/n6gvO5Nj5s3nbZzf7sgg2rTjo7TmrVg82fPdBfmvDDzlx0VFc/5svoqsy/n8Sx8zr5lNvXkO1Frz+Y9/j3+98xKdc2rSg6fZGXLNmTfT09HS6GzaD7d0/wKa7H+VT//Ug257czyt+5lg+fNHPcswo94cdi4d27+eKG37EXbue5qzli3jLy1byyp85NtcXsMzGS9JtEbGm6bw8QS9pHfA3QBn4h4j44Ij53cBngRcBu4E3RsT2dN57gbcCNeBdEXHTaM/loLeJ9PSBQXbs6ePBJ/dz586nuP3hp7hjx1PU6sGZJyzgXa84nXUvOB4p/1k2eVRrdT6/eQef/M42Ht7TR3elxNqVR3P28kW8YNlCViyZy8lHz2H2rPKEPq89dx1R0EsqAz8FXgXsBDYDl0TEPZllfhf42Yh4h6SLgddFxBslrQY+D6wFTgS+AZwREbVWz+egn1kiglo9qEVQr0MtfVwfaotM26H5EYfaB2vBwcEa/dU6BwdrQ9P9gzUODtbpr9Z4tr/KU/sH2ds3wFN9ye/eff081Xfofq5d5RIvWLaAl5x6DL/8whN4/okLJjzgR6rVgx8+uIebtjzGDx7cw08ff5Za/dC/uSXzulkyr4uj53axeG4Xx8ztYk5XhTldZeZ0lTmq8XtWhdmzSswql6iURKVcYlZZVErp77R9VrlEpSzKEiUJlUBAqfFYjWlQ5rcV32hBn+f/kmuBrRGxLd3YjcB64J7MMuuBD6TTXwY+quTdsx64MSL6gQclbU239/3x/CGjeapvgIs+nmw2++E17GMsmk4OLT+8LbtsNG9v8hnZ6rnbba/VsrR67ibP2fL52uyPZts6kv4TUE+Deqoqg13lEovmzGLxnC4WzpnFqUvnsnbl0Zx89Jzk55g5nH7sPLorUzuCLpfES087hpeedgwAfQNVfvr4Ph7avZ+Hd/ex66kD7N4/wN79A9z7yDPs6Rtgf3+VwdrUllRLOvRhgIY/TpuGND4Ysp8PGjlv2PIjlzrUpsOWAdF++1lD2zqCbQzb3qhz2y/Qbv3Rnv/MExbwd5ec3a4HY5Yn6JcBOzKPdwIvabVMRFQlPQ0ck7bfOmLdZSOfQNLlwOUAJ5889nOYIfkH9bzsHYGavMDpc7Voz78sOd9Qh7cf/kYfPt1i/hi216r/w7o/3v632R/Z/pSUvCYliXJJmWmatKUj0NKh+Y32kkRXRcyulOmeVaK7Umb20O+kbXalzKyyCjEyndNV4azlizirzZUxB2t1+gZqHBio0TdQpW8g+V9MtVanWg8Ga3WqtaBarzOY/Z1O1+tBPf3AjUg+nLOPG/Mb7RFBPQ4t0/iwrgfD/gfS0GxAcOhDP+9AJkY8Hvs2oslo5VB/WvexlXYfr+0qIG0/ntsssHxx65vRH4lpcXQoIq4HroekdDOebcyfPYvr3nTOhPbLrFNmlUssPKrEwqNa37PWLK8855LtApZnHp+UtjVdRlIFWEhyUDbPumZmNonyBP1mYJWklZK6gIuBjSOW2Qhclk5fBHwzkv/jbAQultQtaSWwCvjhxHTdzMzyaFu6SWvuVwA3kZxeuSEitki6FuiJiI3Ap4DPpQdb95B8GJAu90WSA7dV4J2jnXFjZmYTz1+YMjObAUY7vdKXQDAzm+Ec9GZmM5yD3sxshnPQm5nNcNPuYKykXuChI9jEEuDJCerORHK/xsb9Ghv3a2xmYr9OiYilzWZMu6A/UpJ6Wh157iT3a2zcr7Fxv8bmudYvl27MzGY4B72Z2Qw3E4P++k53oAX3a2zcr7Fxv8bmOdWvGR3N52EAAASQSURBVFejNzOz4WbiiN7MzDIc9GZmM1xhg17SGyRtkVSXtGbEvPdK2irpPkmvybSvS9u2SrpqCvr4BUl3pD/bJd2Rtq+QdCAz7+OT3ZcR/fqApF2Z578gM6/pvpuifn1Y0k8k3SnpnyUtSts7ur/SPkzpe2eUfiyX9C1J96Tv/99L21u+plPcv+2S7kr70JO2HS3pZkn3p78XT3GfnpfZL3dIekbS73din0naIOkJSXdn2pruHyX+Nn3P3Slp/HdWiohC/gBnAs8DbgHWZNpXAz8GuoGVwAMkl1cup9OnAl3pMqunsL8fAa5Jp1cAd3dw330AeHeT9qb7bgr79Wqgkk7/BfAX02R/dfS9M6IvJwDnpNPzgZ+mr1vT17QD/dsOLBnR9iHgqnT6qsbr2sHX8jHglE7sM+Bc4Jzs+7nV/gEuAL5OcrPOnwd+MN7nLeyIPiLujYj7mswauiF5RDwING5IPnST84gYABo3OZ906Y3Sfx34/FQ83xFote+mRET8R0RU04e3ktyRbDro2HtnpIh4NCJuT6efBe6lyX2Yp5n1wGfS6c8Ar+1gX14JPBARR/Lt+3GLiO+Q3LMjq9X+WQ98NhK3AosknTCe5y1s0I+i2c3Ml43SPhVeDjweEfdn2lZK+pGkb0t6+RT1I+uK9L+DGzL/le7kPhrpLSSjmYZO7q/ptF+GSFoBnA38IG1q9ppOtQD+Q9Jtki5P246LiEfT6ceA4zrTNSC5KVJ2wDUd9lmr/TNh77tpHfSSviHp7iY/HRlNNZOzj5cw/M31KHByRJwNXAncIGnBFPbrY8BpwFlpXz4ykc99BP1qLHM1yR3J/iltmvT9VTSS5gFfAX4/Ip6hg6/pCC+LiHOA84F3Sjo3OzOSmkRHzulWcivUC4EvpU3TZZ8Nmaz90/ZWgp0UEeeNY7XRbkg+4Tcqb9dHJTdL/zXgRZl1+oH+dPo2SQ8AZwATdmutvPtO0ieBf08fTvrN3HPsrzcDvwK8Mn3TT8n+amNa3eRe0iySkP+niPgqQEQ8npmffU2nVETsSn8/IemfScpej0s6ISIeTUsPT3SibyQfPrc39tV02We03j8T9r6b1iP6cWp1Q/I8NzmfDOcBP4mInY0GSUslldPpU9M+bpuCvjSeP1vnex3QOAOgozdzl7QO+CPgwojoy7R3dH/RuffOYdLjPZ8C7o2Iv8y0t3pNp7JvcyXNb0yTHFy/m2RfXZYudhnwr1Pdt9Sw/1lPh32WarV/NgK/lZ598/PA05kSz9hM9VHvCTx6/TqSmlU/8DhwU2be1SRnSdwHnJ9pv4DkLIUHgKunqJ+fBt4xou31wBbgDuB24FeneN99DrgLuDN9M53Qbt9NUb+2ktQk70h/Pj4d9len3jst+vEykv/a35nZTxeM9ppOYd9OJTkj6cfp63V12n4M8J/A/cA3gKM70Le5wG5gYaZtyvcZyQfNo8Bgml9vbbV/SM62uS59z91F5uzCsf74EghmZjPcTCzdmJlZhoPezGyGc9Cbmc1wDnozsxnOQW9mNsM56M3MZjgHvZnZDPf/AUTPp8hr8SiwAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZgc9X3n8fene07d6LAEkkDCFmCwzWEF7PVFHmMbcIziOLZh4yROWBNvQkzWORYHL8vifXaDrxxPSBxs49jENsFHHCWB4NjG69gxRgJzSICMkASSkIQOdM7R093f/aOqpVLTo+mRZqanZj6v5+lnqn5V1f2d6t98puZX1VOKCMzMLP8KrS7AzMxGhgPdzGyCcKCbmU0QDnQzswnCgW5mNkE40M3MJggH+iQm6SZJf9eC190k6ZJ0WpK+IOkFSQ8cx3O9U9JmSQclnT/y1R4fSSHpZaPwvEvS524b6efOvMavSPr2aD2/jR4H+gSWhlztUZXUm5n/lVbXl3o98BZgUURceBzbfxK4NiKmRcRP6xemvzA+JGmNpEOStkj6mqRXnmjhE1VEfDki3lqbH61fTjbyHOgTWBpy0yJiGvAs8I5M25dbXV/qNGBTRBw6ge3XHmP5nwPXAR8CZgNnAN8C3n6cr2c2bjnQrUPSlyQdkLRW0vLaAkmnSPqGpJ2SNkr60GBPIqlb0qckPSNpn6QfSupOl/1q2r5b0g2Zba4GPge8Nv2r4X81eN6CpI+m2z+f1jpTUqekg0AReETS0w22XQb8DnBVRHwvIvojoic9Av2TdJ23S/qppP3p0M1Nme0vlrSl7jmzw0VFSX8s6el0/z0oaXFm9UskPSVpr6Rb078WOiTtyf6FIOklknokzWvwPRQlfVLSLkkbqPtFlO6Lz0vaJmmrpP8tqZgue3/6PnwyHdLaKOmyzLbvl7QhrX1j7a+22nbp9A/S1R9J36P3pn/tvCPzPO1pfeNmyGvSigg/JsED2ARcUtd2E9AHXE4SjP8XuD9dVgAeBG4EOoDTgQ3A2wZ5/luB7wML0+f6T0AncDZwEHhjOv9poFyrBXg/8MNj1P2bwPr09acB3wTuyCwP4GWDbPtB4Jkh9svFwCvT7/dVwA7gFzPLtgy2H4E/BB4DzgQEnAvMydT1z8As4FRgJ3BpuuyvgFsyz3kd8E/H+B6eBBaT/IVxX/rcbenyfwD+BpgKvAR4APitzL4dAD6Qvif/FXgurXUqsB84M133ZOCcRu9J/T4G/gj4+8z8CuCxVvdxP8KBPlkeDB7o38nMnw30ptMXAc/Wrf8R4AsNnrsA9ALnNlh2I3BnZn4qUBpGoH8X+O3M/JlpSNUC7ViBfgPpL6hh7Kc/A/40nR4q0NcBKwZ5ngBen5m/C7g+u28BpfOrgfcM8jzfAz6YmX9rLdCB+UA/0J1ZfhVwX2bfrs8sm5JuuyB9H/YC78pu3+g9aRDopwAHgBnp/NeBP2p1H/cjPORibM9M9wBd6RUUpwGnpMMFeyXtBf6YJETqzQW6gBcNe5D88G+uzUQyVr57GPWdAjyTmX+GI2E2lN0kR56DknSRpPvSYaV9JEfEc5usbTGNv+ea+n07DSAifpLOXyzpLOBlwMpBnuOo/cfR++I0oB3YlnmP/obkSP1FNURETzo5LX0f3kvy/W6T9C9pLUOKiOeAHwHvkjQLuAwYL+dkJjUHug1mM7AxImZlHtMj4vIG6+4iGbp5aYNl20iCDwBJU4A5w6jjOZLgqjmVZMhmRxPbfhdYlD0v0MBXSMJ0cUTMBD5DMiQBcIjkqBZIxrOB7Dj3Zhp/z834IvA+4FeBr0dE3yDrHbX/SL7/7Ov3A3Mz79GMiDinmQIi4t6IeAvJL70ngc8eR/3vBn4cEVuHsa2NEge6DeYB4ICk/56e8CxKeoWkn6tfMSKqwO3Ap9MTqUVJr5XUSfLn+C9Ier2kDuBmhtfvvgr8N0lLJU0D/g/J+G15qA0j4imS8eqvpic4OyR1SbpS0vXpatOBPRHRJ+lC4D9nnuJnJH+xvF1SO/BRkvMANZ8DPiZpWXrC81WSmv1l9XfAO0lC8UvHWO8u4EOSFkk6CajVTURsA74NfErSjPQE8kslvWmoF5c0X9IKSVNJfikcBKqDrL6D5BxG1reAC0jG/49Vv40hB7o1FBEV4BeA84CNJEfhnwNmDrLJH5CcIFwF7AFuAQoRsZbkSpOvkBxtvgBsGeQ5GrkduAP4QVpHH/C7w9j+Q8Bfkpy03UsyRPJO4J/S5b8N3CzpAMl4/121DSNiX7r8c8BWkiP2bO2fTtf/NskJxs8D3c0UFRGbgYdIxqf//Rirfha4F3gkXf+bdct/jeSk9eMk+/brDDHMlCoAHyb5C2gP8CaSk6aN3AR8MR3WeU9afy/wDWBpg5qsRWonZcxsjEm6HXguIj7a6lqOh6QbgTMi4n2trsUSo/bxYTMbnKQlwC8Bubx2W9Js4GqScwA2TnjIxWyMSfoYsAb4RERsbHU9wyXpAyQnZO+JiB8Mtb6NHQ+5mJlNED5CNzObIFo2hj537txYsmRJq17ezCyXHnzwwV0R8aL/+wMtDPQlS5awevXqVr28mVkuSXpmsGUecjEzmyAc6GZmE4QD3cxsgnCgm5lNEEMGuqTbldwpZs0gyyXpLyStl/SopAtGvkwzMxtKM0fofwtceozllwHL0sc1wF+feFlmZjZcQwZ6+tHePcdYZQXwpUjcD8yS1Mx/ezMzsxE0EtehL+ToO6psSdu2jcBzm9kEERFUA8rVKtXq0V8r1aASQbkSVNP1qult1WrT1SoEQdTm69epJl+P2iYard/MOtnl6frVICDZntrtO4/UFBx5/uz3e9Q6kTzHJS+fz7mLZ434Ph7TDxZJuoZkWIZTTz11iLXNrFkRQX+5Sm+pQu9Ahb6BzNdSld7afKlCX7lCqVxloBIMVKrpdJVSJfk6UI7D87VlA5U4srxSPRy85WoSdOVqJKFce0RQqaTtcaTdQIIFM7vGbaBv5ehbZC1K214kIm4DbgNYvny5312zOqVyld2H+tl1oMSuQ/3sOtDP7kMl9vYMcKBvgAN9ZQ70DbA//Xqgr8z+3gF6BiqHjwyHqyDoaCvQXizQUUy+trfp8HRtWXtRTO9qp70gigXRVhQFibaCKBSSr8XaskIhWXbMddRwnaKS5RIUpPSRBKEy8wU1WiddVsjOH1k/+X5FoUDdNsd+zlqbOFKHRDqftNe2JV3nyPpHbzeaRiLQVwLXSrqT5G7m+9JbY5lZnYhg+/4+1m0/wDO7e9i8p4fNL/Sw5YVetu7tZW/PQMPtamE6o6uN6V3tTO9qY+7cqczoamd6VztTO4t0dxTpbk8fHUW6MtPd7el8R5HOtiSka4FdLIxuyNjYGTLQJX0VuBiYK2kL8D9J7jRORHwGuBu4HFhPcifz3xitYs3yZs+hEqs27eGBjXt4bMs+nty+n/19R26H2tVeYNFJU1h8UjfnnzqLl0zvYu60TuZM62DutE7mpdNTOoqjfnRn+TdkoEfEVUMsD5J7RpoZ8PTOg9zz2DbuWbOdtc/tB6CzrcA5p8zgHeeewlkLpnPmghmcPm8qc6Z2OKhtxPgWdGYjoFoN7l27nS/8xyYe2Jhc5XvBqbP4w7edyUVLZ/PKRTPpbCu2uEqb6BzoZifo/g27+dg/P87a5/az6KRuPnLZWaw4byELZna1ujSbZBzoZsepv1zhE/+6js/9cCMLZ3XzZ+89j3ece4pPMlrLONDNjsPB/jK/dcdqfrR+N7/22tP4yGUvp7vDQyrWWg50s2HqL1f4jS88wEPP7uWT7z6XX371olaXZAY40M2G7X98aw2rNr3AX1x1Plece0qryzE7zP8P3WwY/u3xHdy1egvX/vzLHOY27jjQzZrUW6pw4z+u4awF07nukmWtLsfsRTzkYtakO1c9y7Z9ffzpe8+jvehjIRt/3CvNmlCuVLntBxu4cOlsXnP6nFaXY9aQA92sCT9cv4tt+/r4zdctaXUpZoNyoJs14RsPbWXWlHZ+/qyXtLoUs0E50M2GUCpX+e4TO7j8lSf7/7HYuOZANxvCg8+8QE+pwsVnzGt1KWbH5EA3G8IPntpJW0G89qU+GWrjmwPdbAgPPvMCr1g4k+ld7a0uxeyYHOhmx1CpBmu27uPcRTNbXYrZkBzoZsewYedBekoVXrVo5O/QbjbSHOhmx/DY1n0AvNJH6JYDDnSzY3h650GKBbF07tRWl2I2JAe62TFs2tXD4pO6/b9bLBfcS82OYeOuQyzx0bnlhAPdbBARwabdh1gyx4Fu+eBANxvEzgP99JQqHj+33HCgmw1i695eABad1N3iSsya40A3G8SO/f0AzJ/R1eJKzJrjQDcbxPMH+gAHuuWHA91sEDv291EsiDlTO1pdillTHOhmg9ixv5950zopFNTqUsya4kA3G8SO/X3Mn9HZ6jLMmuZANxvEnkMl5kxzoFt+NBXoki6VtE7SeknXN1h+qqT7JP1U0qOSLh/5Us3G1t6eAWZN8f9At/wYMtAlFYFbgcuAs4GrJJ1dt9pHgbsi4nzgSuCvRrpQs7G2r3eAWd0+IWr50cwR+oXA+ojYEBEl4E5gRd06AcxIp2cCz41ciWZjb6BS5WB/2UfolivNBPpCYHNmfkvalnUT8D5JW4C7gd9t9ESSrpG0WtLqnTt3Hke5ZmNjX+8AgAPdcmWkTopeBfxtRCwCLgfukPSi546I2yJieUQsnzfPd1C38WtvTxLoM7sd6JYfzQT6VmBxZn5R2pZ1NXAXQET8GOgC5o5EgWatsK+3BMCsKR5Dt/xoJtBXAcskLZXUQXLSc2XdOs8CbwaQ9HKSQPeYiuVW7Qh9lo/QLUeGDPSIKAPXAvcCT5BczbJW0s2SrkhX+33gA5IeAb4KvD8iYrSKNhttHnKxPGprZqWIuJvkZGe27cbM9OPA60a2NLPWOVQqAzCtq6kfEbNxwZ8UNWvgUH8FgKkdDnTLDwe6WQOH+ssUBF3t/hGx/HBvNWvgUKnM1I42JP+nRcsPB7pZAz39FaZ0FltdhtmwONDNGqgdoZvliQPdrIGeko/QLX8c6GYNHOz3EbrljwPdrIGeUpmpnQ50yxcHulkDPf0VpnR4yMXyxYFu1oBPiloeOdDNGjjUX/GQi+WOA92sgZ5Sme4O/3hYvrjHmtUpV6pUAzrbPIZu+eJAN6tTqlQB6Gjzj4fli3usWZ1SOQ30on88LF/cY83qHA50H6FbzrjHmtXpd6BbTrnHmtWpjaF3OtAtZ9xjzep4DN3yyj3WrI7H0C2v3GPN6viyRcsr91izOh5ysbxyjzWr4yEXyyv3WLM6vmzR8so91qyOL1u0vHKPNatzZAzd/5zL8sWBblbHY+iWV+6xZnVK5QrgQLf8cY81q+Pr0C2v3GPN6vg6dMsr91izOrVAby+qxZWYDU9TgS7pUknrJK2XdP0g67xH0uOS1kr6ysiWaTZ2+itVOtoKSA50y5chb2suqQjcCrwF2AKskrQyIh7PrLMM+Ajwuoh4QdJLRqtgs9FWKlfp9HCL5VAzvfZCYH1EbIiIEnAnsKJunQ8At0bECwAR8fzIlmk2dkrlqk+IWi4102sXApsz81vStqwzgDMk/UjS/ZIubfREkq6RtFrS6p07dx5fxWajzIFueTVSvbYNWAZcDFwFfFbSrPqVIuK2iFgeEcvnzZs3Qi9tNrJKFQe65VMzvXYrsDgzvyhty9oCrIyIgYjYCPyMJODNcqdUrvqSRculZnrtKmCZpKWSOoArgZV163yL5OgcSXNJhmA2jGCdZmOm30MullND9tqIKAPXAvcCTwB3RcRaSTdLuiJd7V5gt6THgfuAP4yI3aNVtNlo8hi65dWQly0CRMTdwN11bTdmpgP4cPowyzUPuVheudea1en3SVHLKfdaszqlctU3t7Bccq81q1MqV3yEbrnkXmtWp1TxGLrlk3utWR1f5WJ55V5rVseBbnnlXmtWJ7ls0TeItvxxoJvV8f9ysbxyrzXLqFaDgUo40C2X3GvNMmo3iPZ16JZH7rVmGbVA92WLlkfutWYZtRtEe8jF8si91izDgW555l5rlnE40D3kYjnkXmuWcXgM3UfolkPutWYZHnKxPHOvNcvod6BbjrnXmmXUjtA7PYZuOeRea5bhMXTLM/daswyPoVueudeaZTjQLc/ca80ySpUK4OvQLZ/ca80yfIRueeZea5bhQLc8c681y+g/fNmi71hk+eNAN8s4/P/Q2/2jYfnjXmuW4X/OZXnmXmuWUSpXaSuIQkGtLsVs2BzoZhmlsm8QbfnlnmuWUao40C2/3HPNMkrlqsfPLbea6rmSLpW0TtJ6SdcfY713SQpJy0euRLOx4yEXy7Mhe66kInArcBlwNnCVpLMbrDcduA74yUgXaTZW+j3kYjnWTM+9EFgfERsiogTcCaxosN7HgFuAvhGsz2xMlcpVOtv8oSLLp2YCfSGwOTO/JW07TNIFwOKI+JdjPZGkayStlrR6586dwy7WbLT1e8jFcuyEe66kAvBp4PeHWjcibouI5RGxfN68eSf60mYjrn+gQqcD3XKqmZ67FVicmV+UttVMB14BfF/SJuA1wEqfGLU8KlWqDnTLrWZ67ipgmaSlkjqAK4GVtYURsS8i5kbEkohYAtwPXBERq0elYrNR1D/gQLf8GrLnRkQZuBa4F3gCuCsi1kq6WdIVo12g2VjqL1d8UtRyq62ZlSLibuDuurYbB1n34hMvy6w1+ss+Qrf8cs81yyiVq/7XuZZb7rlmGf3+6L/lmHuuWUZ/uUJnu8fQLZ8c6GapiEg/KeofC8sn91yzVLkaVMN3K7L8cs81Sx2+QbRPilpOueeapfoHKgC+Dt1yy4FulipV0iN0j6FbTrnnmqX6B5JA939btLxyzzVLHR5D95CL5ZQD3SxVKnvIxfLNPdcs1V9OTop6yMXyyj3XLFUbcnGgW16555qlekrJEfqUDo+hWz450M1SPaUy4EC3/HKgm6V60yP07o6mbhNgNu440M1Sh4dc/N8WLacc6Gap3oHaEboD3fLJgW6W6imVKcjXoVt+ueeapXpKFaZ0tCGp1aWYHRcHulmqt1TxFS6Waw50s1SPA91yzoFuluopVXzJouWaA90s1TtQ9hG65ZoD3SzlIRfLOwe6Waq3VKHbHyqyHHOgm6V8hG5550A3Sx3sLzO9q73VZZgdNwe6GRAR7OsdYEa3r3Kx/HKgm5EMt1SqwQwfoVuOOdDNgH29AwDM7HagW341FeiSLpW0TtJ6Sdc3WP5hSY9LelTSdyWdNvKlmo2e/X1JoM9woFuODRnokorArcBlwNnAVZLOrlvtp8DyiHgV8HXg4yNdqNlo2tfjI3TLv2aO0C8E1kfEhogoAXcCK7IrRMR9EdGTzt4PLBrZMs1G1/6+5PZzHkO3PGsm0BcCmzPzW9K2wVwN3NNogaRrJK2WtHrnzp3NV2k2yjyGbhPBiJ4UlfQ+YDnwiUbLI+K2iFgeEcvnzZs3ki9tdkL299bG0H3ZouVXM713K7A4M78obTuKpEuAG4A3RUT/yJRnNjb29pSQ8AeLLNeaOUJfBSyTtFRSB3AlsDK7gqTzgb8BroiI50e+TLPRtfNgiTlTOygWfLciy68hAz0iysC1wL3AE8BdEbFW0s2SrkhX+wQwDfiapIclrRzk6czGpZ0H+pk7rbPVZZidkKYGDCPibuDuurYbM9OXjHBdZmNq58F+5k13oFu++ZOiZsCuAw50yz8Huk16EeEjdJsQHOg26e3vK1MqV5nnMXTLOQe6TXrP7e0F4OSZ3S2uxOzEONBt0nt2T/JfKxbPdqBbvjnQbdLbXAv0k6a0uBKzE+NAt0lvywu9TO9sY9YUf0rU8s2BbpPeM7sPsWj2FCR/StTyzYFuk97PdhzkjPnTWl2G2QlzoNuktq93gK17ezlrwYxWl2J2whzoNqmt234AgLNOnt7iSsxOnAPdJrUntu0H4KwFDnTLPwe6TWqrn3mB+TM6WTCjq9WlmJ0wB7pNWhHBTzbs5qKlc3yFi00IDnSbtDbuOsTzB/q56PTZrS7FbEQ40G3S+rfHdwDwxmW+v61NDA50m7TuWbOdVyycweLZ/si/TQwOdJuUntpxgIc37+Xtrzyl1aWYjRgHuk1KX/zxJjraCrxn+aJWl2I2YhzoNuls39fH1x/cwhXnnsIc39TCJhAHuk06H7/3SapVuO7Ny1pditmIcqDbpPLttdv55kNb+S9vWOqToTbhONBt0li3/QB/8LVHOOeUGVx3iY/ObeJxoNuk8MS2/fzq539CV3uRz7zv1XS2FVtdktmIc6DbhBYRrHzkOX75r/8DCe64+iIPtdiE1dbqAsxGy1M7DnDLv67jO0/s4LzFs/jM+17Ngpn+J1w2cTnQbUKpVIMfrd/Fnaue5Z4125nSXuQjl53F1a9fSlvRf5DaxOZAt9w72F/mx0/v5vvrnud7Tz7Ptn19nDSlnQ++6aVc84bTOWlqR6tLNBsTDnTLjYhg96ESG3cd4mc7DvDo5n08smUvP9txgGrA1I4ir3vZXG54+8t5y9nzfeLTJh0Huo0LEUHfQJWdB/rZvr+PHZnH9v39PLunhw07D3Kgr3x4m1lT2jl30Szeds4CLlw6m59bMpuONg+r2OTVVKBLuhT4c6AIfC4i/qRueSfwJeDVwG7gvRGxaWRLtfEmIihVqvSWKvSkj2S6TM9AMt1bqqTTZQ71V9jfN8C+ngH29SaPvb1Hpkvl6oteo7OtwPwZXZw6ewq/eN5Cls6dyunzpvLSedNYdFK3b0xhljFkoEsqArcCbwG2AKskrYyIxzOrXQ28EBEvk3QlcAvw3tEoeCxFBBFQjSBIvwYvbqtCEFQj2ab2NbtN/ba17WvrlytJW7kaVKrpdKatmraXqy9uq1SDyiBtlWrdI20bqFQZqFQplasMVIJSuUrpqLZq2hYN2o5MV2N4+3R6Zxszp7Qzszt5nDF/GjO725mRzs+b1smCmV3Mn9HF/OldzOhuc2ibNamZI/QLgfURsQFA0p3ACiAb6CuAm9LprwN/KUkRMcwf96HdtWozt/37BqoRcDgc00CtNgrSo9uq1VqQHgnTbLhG5jknmoKgrVCgUID2QoGOtgLtxdpX0V4s0Jlpm9LRlmnT4fWz63V3FOluLzKlo0h3R5EpHW2Z6WRZrb27vUix4HA2Gy3NBPpCYHNmfgtw0WDrRERZ0j5gDrBrJIrMmjWlnTPnTwdBQaIgEMm0JCTSNlEokLSRWTddRyTzhUKyPLttQUfaCtnnPDytzGvW1kvrKOio1zxq21odheT1s9u2FZLXaiumX9PgbSsUKBZoqq1QgOIgbcW0LjObuMb0pKika4BrAE499dTjeo63nrOAt56zYCTLMjObEJq5JGArsDgzvyhta7iOpDZgJsnJ0aNExG0RsTwils+b5/s4mpmNpGYCfRWwTNJSSR3AlcDKunVWAr+eTv8y8L3RGD83M7PBDTnkko6JXwvcS3LZ4u0RsVbSzcDqiFgJfB64Q9J6YA9J6JuZ2Rhqagw9Iu4G7q5ruzEz3Qe8e2RLMzOz4fDH6szMJggHupnZBOFANzObIBzoZmYThFp1daGkncAzx7n5XEbhU6gjwHUNz3itC8Zvba5reCZiXadFRMMP8rQs0E+EpNURsbzVddRzXcMzXuuC8Vub6xqeyVaXh1zMzCYIB7qZ2QSR10C/rdUFDMJ1Dc94rQvGb22ua3gmVV25HEM3M7MXy+sRupmZ1XGgm5lNEOM+0CW9W9JaSVVJy+uWfUTSeknrJL0t035p2rZe0vVjUOPfS3o4fWyS9HDavkRSb2bZZ0a7lrq6bpK0NfP6l2eWNdx3Y1TXJyQ9KelRSf8gaVba3tL9ldYwpn3nGHUslnSfpMfT/n9d2j7oezqGtW2S9Fj6+qvTttmS/k3SU+nXk8a4pjMz++RhSfsl/V4r9pek2yU9L2lNpq3h/lHiL9L+9qikC07oxZP7aI7fB/By4Ezg+8DyTPvZwCNAJ7AUeJrk3/sW0+nTgY50nbPHsN5PATem00uANS3cdzcBf9CgveG+G8O63gq0pdO3ALeMk/3V0r5TV8vJwAXp9HTgZ+n71vA9HePaNgFz69o+DlyfTl9fe09b+D5uB05rxf4C3ghckO3Lg+0f4HLgHpI7WL4G+MmJvPa4P0KPiCciYl2DRSuAOyOiPyI2AutJbmh9+KbWEVECaje1HnVKbtr5HuCrY/F6J2CwfTcmIuLbEVFOZ+8nuQvWeNCyvlMvIrZFxEPp9AHgCZJ7945XK4AvptNfBH6xhbW8GXg6Io73k+gnJCJ+QHJfiKzB9s8K4EuRuB+YJenk433tcR/ox9Do5tULj9E+Ft4A7IiIpzJtSyX9VNL/k/SGMaoj69r0T7nbM38Gt3If1ftNkiOUmlbur/G0Xw6TtAQ4H/hJ2tToPR1LAXxb0oNK7hMMMD8itqXT24H5Lair5kqOPqhq9f6CwffPiPa5cRHokr4jaU2DR0uOjhppssarOLojbQNOjYjzgQ8DX5E0Ywzr+mvgpcB5aS2fGsnXPoG6auvcAJSBL6dNo76/8kbSNOAbwO9FxH5a+J5mvD4iLgAuA35H0huzCyMZS2jJ9dBKbpN5BfC1tGk87K+jjOb+aeqORaMtIi45js2OdfPqoW5qPWxD1ajk5ti/BLw6s00/0J9OPyjpaeAMYPWJ1tNsXZn6Pgv8czrbzI2/R7UuSe8HfgF4c9rBx2R/DWHU98twSGonCfMvR8Q3ASJiR2Z59j0dMxGxNf36vKR/IBmq2iHp5IjYlg4ZPD/WdaUuAx6q7afxsL9Sg+2fEe1z4+II/TitBK6U1ClpKbAMeIDmbmo9Gi4BnoyILbUGSfMkFdPp09MaN4xBLbXXz47FvROonXUfbN+NVV2XAn8EXBERPZn2lu4vWtd3XiQ9H/N54ImI+HSmfbD3dKzqmippem2a5AT3Go6+UfyvA/84lnVlHPVXcqv3V8Zg+2cl8Gvp1S6vAfZlhmaGbyzP/h7nGeN3kowr9QM7gHszy24guSphHXBZpv1ykqsCngZuGKM6/3xvtzMAAACxSURBVBb4YF3bu4C1wMPAQ8A7xnjf3QE8BjyadpyTh9p3Y1TXepJxw4fTx2fGw/5qVd8ZpI7Xk/xZ/mhmP11+rPd0jOo6neTqn0fS9+qGtH0O8F3gKeA7wOwW7LOpwG5gZqZtzPcXyS+UbcBAml1XD7Z/SK5uuTXtb4+RuZLveB7+6L+Z2QSR5yEXMzPLcKCbmU0QDnQzswnCgW5mNkE40M3MJggHupnZBOFANzObIP4//3t2MpF4X+4AAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEICAYAAAC6fYRZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxcdb3/8dcn+942Tehe2kpZCsgWWdyvgAK/K+Wq1wcKiorww4v7jt6fC8pVr9frT396VVSuiIhwcaEiXFnUH+iPAmURaVhaSqFpQpumJJNlskzm8/vjnKTTNEnTTs5Mcub9fDzmMWfO9yyfk7Tv+eZ7zpwxd0dEROKpKN8FiIhIdBTyIiIxppAXEYkxhbyISIwp5EVEYkwhLyISYwp52YuZfcHMfpbvOgDM7MtmtsvMXjiIdV9hZpvMrMfMzouivoNhZlvN7IyItu1mdlgU2w63/yozeyqq7Us0FPIFJgy9kUfazJIZry/Id30jzGw58DFgjbsvPIhNXAl8x91r3P03E+zj7Wa2ITz2NjO73cxemU3dcebu97r7ESOvo3zDkumjkC8wYejVuHsN8Dzwxox51+e7vgzLgQ5333mQ6x8KbJyo0cw+Cvxv4F+ABeH+/gNYe5D7E5mRFPIynjIz+6mZdZvZRjNrGmkws8Vm9kszazezZ83sgxNtxMx+YmbfN7M7w239XzM7NKPdzeyDZrYlHJb5upkVhb3DO4HFYS/7JxNs/xIz22xmu81snZktDuc/A6wCfhuuXz5mvTkEPf3L3f1X7t7r7kPu/lt3/0S4zMlmdp+ZdYa9/O+YWVnYtiKsvSRjm38ys/eOqe2J8LibzezEjBKON7PHzKzLzG40s4pwncfN7I0Z2ygNfy4nTHD8nwhrazWz94xpKzezfzOz581sR/h7qAzbXmtmLWb2MTPbGW7j3RnrnhPW3G1m283s45nrhdPXEbwxjvyMP2lmvzOzD4yp4zEz+4fx6pcccXc9CvQBbAXOGDPvC0A/cA5QDHwFWB+2FQEPAZ8DygiCdAvwhgm2/xOgG3g1UA58C/hzRrsDfwTqCQLjaeC9YdtrgZZJan8dsAs4Mdz2/wHumezYMtrOAlJAySTbPwk4FSgBVgBPAB8O21aEtZdkLP+njNr/EdgOvAww4DDg0Iy6HgAWh8f9BHBZ2PZJ4MaMba4F/jbJMewAjgGqgZ+HNR0Wtn8TWBfuoxb4LfCVjJ9tiuCNrjT8XfcB88L2NuBV4fQ84MTxfidjf8bAW4H7M14fB3QAZfn+t17ID/XkZTx/dvfb3H0YuI7gPysEodXo7le6+6C7bwF+CJw/ybZ+5+73uPsA8FngNDNbltH+NXff7e7PEwyfvG2KNV4AXOPuD4fbviLc9ooprDsf2OXuqYkWcPeH3H29u6fcfSvwA+A1U6ztvcC/uvuDHtjs7s9ltH/b3VvdfTdB+B4fzv8ZcI6Z1YWv30Hw8x/PW4H/dPfH3b2X4M0ZADMz4FLgI+HPtptgWCrz9zQEXOnBXzC3AT3AERlta8yszt1fdPeHp3jc64DDzWx1Rv03uvvgFNeXCCjkZTyZV7P0ARXh0MShBEMonSMP4DMEY9oT2TYy4e49wG6CXuw+7cBzY9omszhcPnPbHcCSKazbATRkDreMZWaHm9mtZvaCmSUIQrJhirUtA56ZpH3sz7cGwN1bgb8AbzazucDZwETnSRaz789uRCNQBTyU8Xv673D+iI4xb3KjdQBvJujdPxcOsZ02ybGMcvd+4EbgQjMrInjDnuhNSnJEIS8HYhvwrLvPzXjUuvs5k6wz2ms3sxqC4YPW8doJhmwy2ybTSvCmM7LtaoIe+vYprHsfMABMdmnl94AngdXuXkfwZmZhW2/4XJWxfOYVQNuAl0yhjvFcC1xIMORzn7tPdDxt7PuzG7ELSAJHZ/ye5nhwsn2/wr9A1gKHAL8Bbppo0QnqvwA4Hehz9/umsk+JjkJeDsQDQLeZfcrMKs2s2MyOMbOXTbLOOWb2yvCk5ZcIxvcze6CfMLN54RDOhwh6glNxA/BuMzs+PLH6LwTjwVv3t6K7dxGcV/iumZ1nZlXhSc6zzexfw8VqgQTQY2ZHAu/LWL+d4M3kwvBn8B72DvUfAR83s5MscFjmCef9+A3BeYYPAT+dZLmbgHeZ2RozqwI+n1FfmmAY7ZtmdgiAmS0xszfsb+dmVmZmF5jZHHcfIvgZpCdYfAfBeZlRYaingW+gXvyMoJCXKQvH6P+eYAz5WYIe44+AOZOs9nOCANpNcDLzwjHttxCczH0U+B3w4ynWchfwv4BfEvRqX8Lk5wbGrv8N4KPAPwPtBL3v9xOELMDHgbcTnDj+Ifu++VwCfIJg6Odo4P9lbPu/gKsIjr073Gb9FOtKhse0EvjVJMvdTnAO4w/A5vA506fC+evD4aa72DPmvj/vALaG611G0DMfz1eAfw6HhD6eMf+nwLEE5xgkz8xdXxoi0bDg0scWd//nCdqdYDhkc04Lm+HM7HPA4e4+9g1xVjCzdwKXurs+WDYDTHjiSURyz8zqgYsJetOzTjh09E8EHyyTGUDDNSIzhJldQjBsdLu735Pveg5UOObfTjBW//M8lyMhDdeIiMSYevIiIjE2o8bkGxoafMWKFfkuQ0RkVnnooYd2uXvjeG0zKuRXrFjBhg0b8l2GiMisYmbPTdSm4RoRkRhTyIuIxJhCXkQkxhTyIiIxppAXEYkxhbyISIwp5EVEYkwhLyKSZ9+6axP3bmqPZNsKeRGRPHJ3vv2HTdy/ZXck21fIi4jkUe/gMMNpp64ymhsQKORFRPIokRwCoK6iNJLtZx3yZlZhZg+Y2V/NbKOZfTGcv9LM7jezzWZ2Y/gdnyIikiHRH4Z85QwNeYJvvX+dux9H8N2fZ5nZqcDXgG+6+2HAiwTfdiMiIhkSyRQwg3vyHugJX5aGDwdeB9wczr8WOC/bfYmIxM3ocM1MHpM3s2IzexTYCdwJPAN0unsqXKQFWDLBupea2QYz29DeHs0lRCIiM9XocM1M7ckDuPuwux8PLAVOBo48gHWvdvcmd29qbBz3nvciIrG1pyc/g0N+hLt3An8ETgPmmtnI3x9Lge3TuS8RkThI9AcDHrUVM3S4xswazWxuOF0JnAk8QRD2bwkXuwi4Jdt9iYjETSI5RFVZMaXF0VzRPh1vHYuAa82smOBN4yZ3v9XMmoFfmNmXgUeAH0/DvkREYiXRP8SciIZqYBpC3t0fA04YZ/4WgvF5ERGZQCKZiuykK+gTryIieZXoH4rs8klQyIuI5FVXckg9eRGRuAp68gp5EZFYCsbkNVwjIhI76bTTrZ68iEg89Q6mSHt0tzQAhbyISN6MfNpVV9eIiMRQ1F8YAgp5EZG8ifrmZKCQFxHJm9HhGvXkRUTiJ+ovDAGFvIhI3kT9hSGgkBcRyZuusCcf1b3kQSEvIpI3iWSKmvISSiK6lzwo5EVE8ibRPxTpLQ1AIS8ikjeJZLS3NACFvIhI3gQ9eYW8iEgsJZKpSC+fBIW8iEjeqCcvIhJjGpMXEYmpdNrpHoj2C0NAIS8ikhc9gynco705GUxDyJvZMjP7o5k1m9lGM/tQOL/ezO40s03h87zsyxURiYdc3GYYpqcnnwI+5u5rgFOBy81sDfBp4G53Xw3cHb4WERGCK2sg2puTwTSEvLu3ufvD4XQ38ASwBFgLXBsudi1wXrb7EhGJi65Z1JMfZWYrgBOA+4EF7t4WNr0ALJhgnUvNbIOZbWhvb5/OckREZqzRO1DO9DH5EWZWA/wS+LC7JzLb3N0BH289d7/a3ZvcvamxsXG6yhERmdFGxuTnzIaQN7NSgoC/3t1/Fc7eYWaLwvZFwM7p2JeISBzk4luhYHqurjHgx8AT7v7vGU3rgIvC6YuAW7Ldl4hIXIz05Gsivk5+Orb+CuAdwN/M7NFw3meArwI3mdnFwHPAW6dhXyIisZDoH6K2vITiIot0P1mHvLv/GZioytOz3b6ISBwFNyeLdqgG9IlXEZG8SPQPRfq1fyMU8iIieZCLm5OBQl5EJC8S/anIr6wBhbyISF4EPXkN14iIxFIuvjAEFPIiIjk3nHa6+3V1jYhILO25zbCGa0REYmfTzh4AXtJYE/m+FPIiIjnW3NoFwJrFdZHvSyEvIpJjzW0J5leXcUhteeT7UsiLiOTYxtYEaxbXEdzfMVoKeRGRHBpMpdm0o4c1i6IfqgGFvIhITj3T3sPgcDon4/GgkBcRyanm1uCL845WyIuIxE9zW4KK0iJWNkR/+SQo5EVEcqq5NcERC+si/7KQEQp5EZEccXc2tnbl7KQrKORFRHJme2eSRH8qZyddQSEvIpIzIydd1ZMXEYmh5rYEZnDUotqc7VMhLyKSI82tCVY2VFNVFv3dJ0co5EVEcqS5LZHToRqYppA3s2vMbKeZPZ4xr97M7jSzTeHzvOnYl4jIbNTVN0TLi8mcnnSF6evJ/wQ4a8y8TwN3u/tq4O7wtYhIQbpnUzsAJyzLbX93WkLe3e8Bdo+ZvRa4Npy+FjhvOvYlIjIb3fxQC4vnVHDyyvqc7jfKMfkF7t4WTr8ALBhvITO71Mw2mNmG9vb2CMsREcmPF7r6uXdTO286cWnOPuk6IicnXt3dAZ+g7Wp3b3L3psbGxlyUIyKSU79+ZDtphzeftDTn+44y5HeY2SKA8HlnhPsSEZmR3J2bH9pG06HzWNlQnfP9Rxny64CLwumLgFsi3JeIyIz06LZOnmnv5S156MXD9F1CeQNwH3CEmbWY2cXAV4EzzWwTcEb4WkSkoNz8UAsVpUX8j5cuysv+p+VjV+7+tgmaTp+O7YuIzEb9Q8Os+2srZx+ziNqK0rzUoE+8iohE5IYHnqe7P5W3oRpQyIuIROKFrn6+ccfTvPrwRl7+kvl5q0MhLyISgS/+diNDw2m+vPYYzHJ7bXwmhbyIyDS7q3kHtz/+Ah88fTXL51fltRaFvIjINOodSPH5dRs5fEENl7xqVb7LmZ6ra0REBFLDaT5y46Ns70xy82WnUVaS/350/isQEYmBdNr55C8f447mHXzx3KNpWpHbG5FNRCEvIpIld+fKW5v51cPb+eiZh3PRy1fku6RRGq4REcnCYCrNl25t5rr1z3HxK1fygdcdlu+S9qKQFxE5SNs7k1x+/cM8uq2TS1+9iivOPjKvl0uORyEvInKA3J3fb9zBFb96jKFh53sXnMjZx+bn3jT7o5AXETkAT+/o5srfNvPnzbs4alEd3337CaxqrMl3WRNSyIuITMEz7T386N4t3LShhZryEr547tFccMpySopn9vUrCnkRkQmk0876LR1c85et3PXEDspLirjwlOV86IzDqa8uy3d5U6KQFxHJ4O48vaOHWx7dzm8e2U5rVz/11WV86PTVvOO0Q2moKc93iQdEIS8iBS85OMyDW3fzhyd3cveTO9i2O0lxkfGq1Q18+pyjeP2aBVSUFue7zIOikBeRgtPRM8Cj2zp55PlO1m/p4K8tnQwNO+UlRbzisAYue81LOHPNAg6prch3qVlTyItIbA2nnZYX+3h6Rw9PtCV4oi3B461dbNudBKC4yDhmyRwufuUqTl1Vzykr51NZNjt77BNRyIvIrDacdnYk+ml5MclzHb1s7ehl664+tuzqZUt7DwOpNABmsGJ+NccumcOFpxzK8cvmcuzSOVSVxTsG4310IjKrDabStPcM0N49wI5E/+ijrauf1s4krZ39tHUlGRr20XWKi4xl8ypZ0VDNKw+bz2GH1HDYIbUcubCW6vLCi7zCO2IRyYt02ukeSNHVN0RncpAX+4bo7Bvkxd5BdvcNsbt3gN29g+zqGaSjZ4CO3kE6+4b22U5JkbGgroJFcyo4btlczjl2EcvqK1k2r4pl9VUsnVdJ6Qy/dj2XFPIisl9Dw2n6BobpGUzR05+iZyBF70Dw3NOfonsgRXf/EN39wXMimaJ7IHjuSg6R6B8ikRwi7eNv3wzmVpYyr7qMhppyjlhYy/zqchprw0dNOQvqKlgwp5yG6nKKimbW/WFmsshD3szOAr4FFAM/cvevRr1PkULj7gyk0iQHh0kOhY/BYfozpkee+/aZTtE3OEzvQDDdOzBM3+DIvBS9g8MMhuPa+1NdVkxtRSl1lSXUVpQyv6aMVY3VzKksHX3MrSpjTmUp86qC6ZHnYgV3JCINeTMrBr4LnAm0AA+a2Tp3b45yvyL5lBpOM5BK0z80zEAqHT6GGRjae97IdP/Q8J75e81L058a3jM9NEx/uMx47QeqyKCqrITKsmKqyoqpLC2muryE2ooSFtZVUFUezK8uL6GmrISq8hJqy0uoLi+huryY2ooSaspLg+nyUmoqShTUM1DUPfmTgc3uvgXAzH4BrAUU8hKZkV7twFAYrmHI9g/tCdvR4A2X6x+dnxHM44T02PX22m647PBEYxJTVFFaRHlJMRWlRVSUFlNeEjxXlBRTV1FCRW15+LqIyjCcy0uD5atKi4O20uLRtrHPVWXFo9udabfFlekXdcgvAbZlvG4BTslcwMwuBS4FWL58ecTlSD64O4PD6b2GDJIZvdG9X49M79ubzQzVvXrBI0Eb9nSnOrQwETOoKCmmrKRoNHBHgra8pIiykiJqK0pHpytKiikvLaKseM8yma8rSife1kiYl4++VvDK9Mr7iVd3vxq4GqCpqSm7LpBkbTjt9GacXBs5wRacZAvGanvHPPcNZj6HY71DKZKDaZKDKZJDwxOecJtMcZFRUbInKMtLg0Ad6eHWVZZSMSZMg4At3qs3PBKsY4N67HRFaTFlxUWUFpuCVmIj6pDfDizLeL00nCcRcQ8uU+vsDS5T6+wbois5RGcyuLohMXKlQ3+KRHLP1RDdYaj3DQ5PaT9mUFVaTFV5CdVlxVSWlVBVFozTHlJbHozxhvNGhgdGnitHhxSKgjHh0j3BXZExrcvgRLIXdcg/CKw2s5UE4X4+8PaI9xk7A6lhdiYGRj8U0t49QEfPILt6BujoDaZ39w7yYl8Q6qlJus3l4VBDXWUJdRWl1FaUsGRuZXgSrYSakedwujqcri4LTrZVlQWvK0o1rCAyG0Qa8u6eMrP3A78nuITyGnffGOU+Zxt3p71ngG27k7S82EfLi0naupK0dfbT2hV8um937+C4686tKqW+uoyG6nJWNVZTX13PvKpS5lWVMTfjeeTStWB4I1735RCRyUU+Ju/utwG3Rb2fmS45OMwz7T1s2tnNph09PLurl60dfTzX0bvPEMncqlIWz6lk8ZwKTlg+l4V1FRxSW84hdeUcUltBY2059dVlGs4Qkf3K+4nXOOoZSPHYtk4ebemkuTVBc1uCZ3f14uEoSkmRsby+ihUN1Zy2aj6Hzq9iWX0lS+dVsWRuZUHeX0NEoqE0mQa7ewdZv6WD9Vs6uH/Lbp7e2T0a6MvqK1mzqI5zj1vMEQtqWb2ghuX11ZSVqBcuItFTyB8Ed6e5LcHdT+zk7id38tdtnQBUlRVz0qHzOPvYhZywfB7HLZ3D3KrZ8T2QIhJPCvkDsG13H795ZDu/fmQ7W3b1YgbHLZ3Lx848nJcf1sBLl87ROLmIzCgK+f1wd+57poOr793Cn55qB+DklfVc8upVnHHUAhprZ9eX+opIYVHIT8DduaN5B9++exMbWxM01JTx4TNW8+YTl7Ksvirf5YmITIlCfhxPvdDNlbdu5C+bO1jVUM1X33Qs552wRNeYi8iso5DPMDSc5t/ueIof3fssNeUlfPHco7nglOWUaJxdRGYphXxoR6Kf9//8YR7c+iLnv2wZnzzrSOqrdWWMiMxuCnngvmc6+MAND9M3OMy3zj+etccvyXdJIiLTouBDfv2WDt71nw+wdF4lN1xyKqsX1Oa7JBGRaVPQIf9YSyfvvXYDy+uruPF/nqbhGRGJnYI9o7h5ZzcXXfMAc6tKue7iUxTwIhJLBRnyvQMpLrrmQYqLivjZxaewcE5FvksSEYlEQQ7XfPPOp9nemeTmy05jRUN1vssREYlMwfXkH9/exTV/eZa3n7KcphX1+S5HRCRSBRXyw2nnM7/+G/XV5XzqDUfmuxwRkcgVVMhfd99WHmvp4nNvXMOcqtJ8lyMiErmCCfmegRTfuONpXrW6gTe+dFG+yxERyYmCCfnb/tZG90CKD59xOGaW73JERHKiYEL+5odaWNVYzYnL5+a7FBGRnCmIkH+uo5cHnt3NW05aql68iBSUrELezP7RzDaaWdrMmsa0XWFmm83sKTN7Q3ZlZueXD7VQZPCmE5bmswwRkZzL9sNQjwNvAn6QOdPM1gDnA0cDi4G7zOxwdx/Ocn8HLJ12fvnwdl65ulGfbBWRgpNVT97dn3D3p8ZpWgv8wt0H3P1ZYDNwcjb7Oljrt3SwvTPJW05SL15ECk9UY/JLgG0Zr1vCefsws0vNbIOZbWhvb5/2Qm5+qIXaihJev2bBtG9bRGSm22/Im9ldZvb4OI+101GAu1/t7k3u3tTY2DgdmxzVO5DitsfbeONxi/X9rCJSkPY7Ju/uZxzEdrcDyzJeLw3n5dTG1gT9Q2nOPEq9eBEpTFEN16wDzjezcjNbCawGHohoXxPa2NoFwNGL63K9axGRGSHbSyj/wcxagNOA35nZ7wHcfSNwE9AM/DdweT6urGluTdBQU0ZjbXmudy0iMiNkdQmlu/8a+PUEbVcBV2Wz/Ww1tyU4alGdPgAlIgUrtp94HUyl2bSjhzUaqhGRAhbbkH+mvYfB4TRHL56T71JERPImtiG/sTUBwJpF6smLSOGKbcg3tyaoKC1ipb7DVUQKWHxDvq2LIxfWUVykk64iUrhiGfLuTnNrQiddRaTgxTLkt3cmSfSnNB4vIgUvliE/ctJVn3QVkUIXy5Bvbk1QZHDkQoW8iBS2eIZ8W4KVDdVUlunOkyJS2OIZ8q0J1uhDUCIi8Qv5rr4htncmddJVRIQYhvzm9h4AjlhYk+dKRETyL3Yh35UcBGB+tW4vLCISu5BPJFMA1FWW5rkSEZH8i1/I9w8BUFeR1a3yRURiIX4hnwxCvrZCPXkRkdiFfFdyiMrSYspKYndoIiIHLHZJmEimqKvUUI2ICMQx5PuHqNNQjYgIENeQ15U1IiJAHEM+mdKVNSIioaxC3sy+bmZPmtljZvZrM5ub0XaFmW02s6fM7A3Zlzo16smLiOyRbU/+TuAYd38p8DRwBYCZrQHOB44GzgL+w8xyckvIRHKIOQp5EREgy5B39zvcPRW+XA8sDafXAr9w9wF3fxbYDJyczb6mWA+J/pROvIqIhKZzTP49wO3h9BJgW0ZbSzhvH2Z2qZltMLMN7e3tWRXQNzjMcNp1CaWISGi/aWhmdwELx2n6rLvfEi7zWSAFXH+gBbj71cDVAE1NTX6g62fac0sD9eRFRGAKIe/uZ0zWbmbvAv4eON3dR0J6O7AsY7Gl4bxI6eZkIiJ7y/bqmrOATwLnuntfRtM64HwzKzezlcBq4IFs9jUVXUn15EVEMmU7eP0doBy408wA1rv7Ze6+0cxuApoJhnEud/fhLPe1XyM3J9OYvIhIIKs0dPfDJmm7Crgqm+0fKI3Ji4jsLVafeN3Tk1fIi4hA3EK+PzjxWqvbGoiIAHEL+eQQVWXFlBbH6rBERA5arNJQtxkWEdlbvEJeXxgiIrKXeIW8evIiInuJX8jryhoRkVHxCnl9YYiIyF5iFfJdupe8iMheYhPy6bTTreEaEZG9xCbkewdTpF23NBARyRSbkB/5tKsuoRQR2SM+Ia/bDIuI7CN+Ia8xeRGRUfEJ+ZHhGvXkRURGxSfk9YUhIiL7iE/I6wtDRET2EZ+QT+pe8iIiY8Un5PuHqC4rpkT3khcRGRWbROxK6tOuIiJjxSbkE0ndZlhEZKz4hHz/kK6sEREZI6uQN7MvmdljZvaomd1hZovD+WZm3zazzWH7idNT7sSC2wyrJy8ikinbnvzX3f2l7n48cCvwuXD+2cDq8HEp8L0s97NfiX7dZlhEZKysQt7dExkvqwEPp9cCP/XAemCumS3KZl/7k9CJVxGRfWQ9iG1mVwHvBLqAvwtnLwG2ZSzWEs5rG2f9Swl6+yxfvvygakinne4BfSuUiMhY++3Jm9ldZvb4OI+1AO7+WXdfBlwPvP9AC3D3q929yd2bGhsbD/wIgJ7BFO66OZmIyFj77fq6+xlT3Nb1wG3A54HtwLKMtqXhvEjoNsMiIuPL9uqa1Rkv1wJPhtPrgHeGV9mcCnS5+z5DNdNl5JYGuoRSRGRv2abiV83sCCANPAdcFs6/DTgH2Az0Ae/Ocj+T0s3JRETGl1XIu/ubJ5jvwOXZbPtAdOkLQ0RExhWLT7zOry7j7GMW0lhbnu9SRERmlFgMYjetqKdpRX2+yxARmXFi0ZMXEZHxKeRFRGJMIS8iEmMKeRGRGFPIi4jEmEJeRCTGFPIiIjGmkBcRiTEL7kAwM5hZO8E9cKaqAdgVUTkzWaEeNxTuseu4C8uBHveh7j7uvdpnVMgfKDPb4O5N+a4j1wr1uKFwj13HXVim87g1XCMiEmMKeRGRGJvtIX91vgvIk0I9bijcY9dxF5ZpO+5ZPSYvIiKTm+09eRERmYRCXkQkxmZFyJvZWWb2lJltNrNPj9NebmY3hu33m9mK3Fc5/aZw3B81s2Yze8zM7jazQ/NR53Tb33FnLPdmM3Mzi8UldlM5bjN7a/g732hmP891jVGZwr/15Wb2RzN7JPz3fk4+6pxOZnaNme00s8cnaDcz+3b4M3nMzE48qB25+4x+AMXAM8AqoAz4K7BmzDL/BHw/nD4fuDHfdefouP8OqAqn31coxx0uVwvcA6wHmvJdd45+36uBR4B54etD8l13Do/9auB94fQaYGu+656G4341cCLw+ATt5wC3AwacCtx/MPuZDT35k4HN7r7F3QeBXwBrxyyzFrg2nL4ZON3MLIc1RmG/x+3uf3T3vvDlemBpjmuMwlR+3wBfAr4G9OeyuAhN5bgvAb7r7i8CuPvOHNcYlakcuwN14fQcoDWH9UXC3e8Bdk+yyFrgpx5YD8w1s0UHup/ZEPJLgG0Zr1vCeeMu4+4poAuYn5PqojOV4850McG7/my33+MO/wIymF0AAAIFSURBVGxd5u6/y2VhEZvK7/tw4HAz+4uZrTezs3JWXbSmcuxfAC40sxbgNuADuSktrw40A8YViy/yLnRmdiHQBLwm37VEzcyKgH8H3pXnUvKhhGDI5rUEf7XdY2bHuntnXqvKjbcBP3H3b5jZacB1ZnaMu6fzXdhMNxt68tuBZRmvl4bzxl3GzEoI/pzryEl10ZnKcWNmZwCfBc5194Ec1Ral/R13LXAM8Ccz20owVrkuBidfp/L7bgHWufuQuz8LPE0Q+rPdVI79YuAmAHe/D6gguIlXnE0pA/ZnNoT8g8BqM1tpZmUEJ1bXjVlmHXBROP0W4A8enrmYxfZ73GZ2AvADgoCPy/jspMft7l3u3uDuK9x9BcG5iHPdfUN+yp02U/l3/huCXjxm1kAwfLMll0VGZCrH/jxwOoCZHUUQ8u05rTL31gHvDK+yORXocve2A93IjB+ucfeUmb0f+D3BWfhr3H2jmV0JbHD3dcCPCf5820xwIuP8/FU8PaZ43F8HaoD/Cs8zP+/u5+at6GkwxeOOnSke9++B15tZMzAMfMLdZ/tfrFM99o8BPzSzjxCchH3XbO/ImdkNBG/aDeG5hs8DpQDu/n2Ccw/nAJuBPuDdB7WfWf5zEhGRScyG4RoRETlICnkRkRhTyIuIxJhCXkQkxhTyIiIxppAXEYkxhbyISIz9f5quit1M6IRJAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "20%, 50% and 60% quantiles: [-1.37638192  0.          0.3249197 ]\n"
          ]
        }
      ],
      "source": [
        "# Getting theoretical quantities on distributions\n",
        "cauchy = scipy.stats.cauchy # many others are available\n",
        "# the probability density function\n",
        "X = np.arange(-100, 100, 0.1)\n",
        "plt.figure()\n",
        "plt.title(\"The pdf of Cauchy density\")\n",
        "plt.plot(X, cauchy.pdf(X))\n",
        "plt.show()\n",
        "# the cumulative distribution function\n",
        "X = np.arange(-100, 100, 0.1)\n",
        "plt.figure()\n",
        "plt.title(\"The cdf of Cauchy density\")\n",
        "plt.plot(X, cauchy.cdf(X))\n",
        "plt.show()\n",
        "# the percent point function (inverse of cdf)\n",
        "X = np.arange(0, 1, 0.01)\n",
        "plt.figure()\n",
        "plt.title(\"The ppf of Cauchy density\")\n",
        "plt.plot(X, cauchy.ppf(X))\n",
        "plt.show()\n",
        "# Theoretical 20%, 50% and 60% quantiles of cauchy distribution\n",
        "quantiles = cauchy.ppf((0.2, 0.5, 0.6))\n",
        "print(\"20%, 50% and 60% quantiles:\", quantiles)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C2ay7KWSn5uC"
      },
      "source": [
        "1. Generate 100 independent copies of N(1,3) random variables and plot the QQ-plot\n",
        "against the theoretical quantiles of the N(0,1), N(5, 3) and N(1,3) distribution. Superpose on the plot the line y=x.\n",
        "\n",
        "2. Generate 100 independent copies of chi-square random variables with 5 degrees of freedom,\n",
        "and compare those aginst the theoretical quantiles of the N(5,1) distribution, the $\\chi^2_4$ distribution, the $\\chi^2_5$ distribution, the $\\chi^2_6$ distribution.Superpose on the plot the line y=x."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "id": "aoYHmjIYn5uD",
        "outputId": "f53ce641-7bf2-476f-d07d-24fbd7552758",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 354
        }
      },
      "outputs": [
        {
          "output_type": "error",
          "ename": "NotImplementedError",
          "evalue": "ignored",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
            "\u001b[0;32m<ipython-input-26-fa407db4868c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     22\u001b[0m     \u001b[0msm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mqqplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mline\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'45'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 24\u001b[0;31m \u001b[0mtest_qqplot_sm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
            "\u001b[0;32m<ipython-input-26-fa407db4868c>\u001b[0m in \u001b[0;36mtest_qqplot_sm\u001b[0;34m()\u001b[0m\n\u001b[1;32m     19\u001b[0m     \"\"\"\n\u001b[1;32m     20\u001b[0m     \u001b[0mX\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnormal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 21\u001b[0;31m     \u001b[0mqqplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstats\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnorm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     22\u001b[0m     \u001b[0msm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mqqplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mline\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'45'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-26-fa407db4868c>\u001b[0m in \u001b[0;36mqqplot\u001b[0;34m(data, distribution)\u001b[0m\n\u001b[1;32m     12\u001b[0m         \u001b[0mTheoretical\u001b[0m \u001b[0mdistribution\u001b[0m \u001b[0mto\u001b[0m \u001b[0mcompare\u001b[0m \u001b[0mto\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m     \"\"\"\n\u001b[0;32m---> 14\u001b[0;31m     \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     15\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mtest_qqplot_sm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mNotImplementedError\u001b[0m: "
          ]
        }
      ],
      "source": [
        "def qqplot(data, distribution):\n",
        "    \"\"\"\n",
        "    Generate a qqplot using the quantiles of the data and the theoretical data of the distribution.\n",
        "    Superpose on the plot the line y=x.\n",
        "    \n",
        "    Parameters\n",
        "    ----------\n",
        "    data: np array of shape (100,)\n",
        "        Empirical data\n",
        "    \n",
        "    distribution: distribution from scipy.stats\n",
        "        Theoretical distribution to compare to\n",
        "    \"\"\"\n",
        "    raise NotImplementedError\n",
        "    \n",
        "def test_qqplot_sm():\n",
        "    \"\"\"\n",
        "    Test the plot is the same as in sm.qqplot\n",
        "    \"\"\"\n",
        "    X = np.random.normal(0, 1, size=100)\n",
        "    qqplot(X, scipy.stats.norm)\n",
        "    sm.qqplot(X, line='45')\n",
        "    \n",
        "test_qqplot_sm()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6EVqMY6Pn5uD"
      },
      "outputs": [],
      "source": [
        "raise NotImplementedError\n",
        "# Answer 1."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZTYG9Dhnn5uE"
      },
      "outputs": [],
      "source": [
        "# Answer 2.\n",
        "raise NotImplementedError"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "42H_l4tEn5uF"
      },
      "source": [
        "# Generate a simple regression model\n",
        "\n",
        "Given the sample size n and dimension p, generate a dataset (y, X) where y has size n and X has size n times p,\n",
        "such that\n",
        "$y = X \\beta^* + \\beta_0^* + \\epsilon$ where $\\epsilon$ has iid $N(0,\\sigma^2)$ entries, $X$ has iid $N(0,1)$ entries, and $\\beta^*\\in R^p$ has norm equal to the signal strength parameter, and $\\beta_0^*$ is the intercept."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "metadata": {
        "id": "gZc_ejfZn5uG"
      },
      "outputs": [],
      "source": [
        "def linear_regression(X, beta_true, sigma):\n",
        "    \"\"\"\n",
        "    Applies the linear regression model\n",
        "    \n",
        "    Parameters\n",
        "    -----------\n",
        "    \n",
        "    X: np array of shape (n, p)\n",
        "    input features\n",
        "\n",
        "    beta_true: np array of shape (p)\n",
        "    true coefficients\n",
        "    \n",
        "    sigma: float\n",
        "    noise standard deviation \n",
        "\n",
        "    Return\n",
        "    ------\n",
        "    y: np array of shape (n,)\n",
        "    responses generated accroding to the linear regression model\n",
        "    y= X beta_star + eps\n",
        "    \"\"\"\n",
        "    raise np.matmul(X,b_true)+np.random.normal(0, sigma**2, size = (n))\n",
        "\n",
        "def generate_dataset(n, p, sigma=0.7, signal_strength=2):\n",
        "\n",
        "    X = np.random.normal(0, 1, size = (n,p))\n",
        "    b_true = np.random.uniform(0, 1, size = (p))\n",
        "    b_true = (b_true / np.linalg.norm(b_true))*signal_strength\n",
        "    eps = np.random.normal(0, sigma**2, size = (n))\n",
        "    y = np.matmul(X,b_true)+eps\n",
        "\n",
        "    \"\"\"\n",
        "    Generate a data set (with intercept = 0)\n",
        "    Take beta_star uniformly in the sphere of length: signal_strength\n",
        "    Parameters\n",
        "    ----------\n",
        "    n: int\n",
        "        number of samples\n",
        "    p: int\n",
        "        dimension\n",
        "    sigma: float\n",
        "        noise standard deviation\n",
        "    signal_strength: float\n",
        "        norm of the true parameter b_true\n",
        "    \n",
        "    Return\n",
        "    X: np array of shape (n, p)\n",
        "        Features\n",
        "    y: np array of shape (n,)\n",
        "        Target\n",
        "    b_true: np array of shape (p,)\n",
        "    \"\"\"\n",
        "    return X, y, b_true\n",
        "\n",
        "def test_nointercept():\n",
        "    \"\"\"\n",
        "    Check that X.dot(b_true) = y when the noise is 0\n",
        "    \"\"\"\n",
        "    X, y, b_true = generate_dataset(100, 10, sigma=0, signal_strength=1)\n",
        "    np.testing.assert_allclose(X.dot(b_true), y)\n",
        "    \n",
        "test_nointercept()\n",
        "\n",
        "def generate_dataset_with_intercept(n, p, sigma=0.7, signal_strength=2):\n",
        "    \"\"\"\n",
        "    Generate a data set, same as before, but here the returned array X should have\n",
        "    p+1 columns, and the first column contains only 1,1,...,1 to account for\n",
        "    possible intercept.\n",
        "    Parameters\n",
        "    ----------\n",
        "    n: int\n",
        "        number of samples\n",
        "    p: int\n",
        "        dimension\n",
        "    sigma: float\n",
        "        noise standard deviation\n",
        "    signal_strength: float\n",
        "        norm of the true parameter b_true\n",
        "    \n",
        "    Return\n",
        "    X: np array of shape (n, p + 1)\n",
        "        Features, last column contains only ones\n",
        "    y: np array of shape (n,)\n",
        "        Target\n",
        "    b_true: np array of shape (p,)\n",
        "        true coefficients (last one is the intercept)\n",
        "    \"\"\"\n",
        "    X = np.ones((n,p+1))\n",
        "    X[:,:-1] = np.random.normal(0,1,size = (n,p))\n",
        "    b_true = np.ones((p+1))\n",
        "    b_true[:-1] = np.random.uniform(0, 1, size = (p))\n",
        "    b_true[:-1] = (b_true[:-1] / np.linalg.norm(b_true[:-1]))*signal_strength\n",
        "    eps = np.random.normal(0, sigma**2, size = (n))\n",
        "    y = np.matmul(X,b_true)+eps\n",
        "    return X, y, b_true\n",
        "\n",
        "def test_intercept():\n",
        "    \"\"\"\n",
        "    Check that X.dot(b_true) = y when the noise is 0 and intercept = 0\n",
        "    \"\"\"\n",
        "    X, y, b_true = generate_dataset_with_intercept(100, 10, sigma=0, signal_strength=1)\n",
        "    np.testing.assert_allclose(X[:, :-1].dot(b_true[:-1]) + b_true[-1], y)\n",
        "    \n",
        "test_nointercept()\n",
        "test_intercept()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wQTrTyMpn5uI"
      },
      "source": [
        "# Least-squares\n",
        "\n",
        "Given a design matrix X and response y,\n",
        "implement the least-squares estimate $\\hat\\beta = (X^TX)^\\dagger X^Ty$\n",
        "using several methods:\n",
        "\n",
        "1. using the SVD of X `np.linalg.svd`\n",
        "2. using the pseudo-inverse of X `np.linalg.pinv`\n",
        "3. using the inverse of $X^TX$ if full-rank, and its pseudo-inverse otherwise\n",
        "4. using `scipy.linalg.lstsq`\n",
        "\n",
        "For $p=10, n=20$, generate a dataset using the previous functions and verify that all methods give the correct estimate of $\\beta$ when there is no noise.\n",
        "\n",
        "For $p=10, n=20$, generate a dataset using the previous functions and verify that all methods give the same estimate.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 69,
      "metadata": {
        "id": "h3O4tMEqn5uJ",
        "outputId": "9831ea38-6253-49d7-fa2d-d8d7038f48b4",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Méthode 1: Le resultat est-il correct?  True\n",
            "Méthode 2: Le resultat est-il correct?  True\n",
            "Méthode 3: Le resultat est-il correct?  True\n",
            "Méthode 4: Le resultat est-il correct?  True\n",
            "Le même β est-il retourné? : True\n"
          ]
        }
      ],
      "source": [
        "def compute_least_squares_using_SVD(X, y):\n",
        "    U, S, Vt = np.linalg.svd(X, full_matrices = False)\n",
        "    VinvSUt = np.matmul(np.matmul(Vt.T, np.diag(1/S)), U.T)\n",
        "    return np.matmul(VinvSUt, y)\n",
        "    \n",
        "def compute_least_squares_using_pseudo_inverse(X, y):\n",
        "    return np.dot(np.linalg.pinv(X), y)\n",
        "    \n",
        "def compute_least_squares_using_method3(X, y):\n",
        "    return np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), y)\n",
        "\n",
        "def compute_least_squares_using_scipy(X, y):\n",
        "    return scipy.linalg.lstsq(X, y)[0]\n",
        "    \n",
        "def test_lsqr():\n",
        "    n = 20\n",
        "    p = 10\n",
        "    X, y, b_true = generate_dataset_with_intercept(n, p, sigma=0, signal_strength=2)\n",
        "\n",
        "    b_hat = compute_least_squares_using_SVD(X, y)\n",
        "    print(\"Méthode 1: Le resultat est-il correct? \", np.min(np.isclose(b_true, b_hat)))\n",
        "\n",
        "    b_hat = compute_least_squares_using_pseudo_inverse(X, y)\n",
        "    print(\"Méthode 2: Le resultat est-il correct? \", np.min(np.isclose(b_true, b_hat)))\n",
        "\n",
        "    b_hat = compute_least_squares_using_method3(X, y)\n",
        "    print(\"Méthode 3: Le resultat est-il correct? \", np.min(np.isclose(b_true, b_hat)))\n",
        "\n",
        "    b_hat = compute_least_squares_using_scipy(X, y)\n",
        "    print(\"Méthode 4: Le resultat est-il correct? \", np.min(np.isclose(b_true, b_hat)))\n",
        "    \n",
        "    \"\"\"Test that all methods yield the correct beta when there is no noise\n",
        "    \"\"\"\n",
        "    \n",
        "\n",
        "def test_same_estimate():\n",
        "    \"\"\" Test that all methods yield the same estimate of beta\n",
        "    \"\"\"\n",
        "\n",
        "    flag = True \n",
        "    n = 20\n",
        "    p = 10\n",
        "    X, y, b_true = generate_dataset_with_intercept(n, p, sigma=0, signal_strength=2)\n",
        "\n",
        "    b_hat1 = compute_least_squares_using_SVD(X, y)\n",
        "    b_hat2 = compute_least_squares_using_pseudo_inverse(X, y)\n",
        "    b_hat3 = compute_least_squares_using_method3(X, y)\n",
        "    b_hat4 = compute_least_squares_using_scipy(X, y)\n",
        "\n",
        "    flag = np.min(np.isclose(b_hat1, b_hat2))\n",
        "    flag = np.min(np.isclose(b_hat1, b_hat3))\n",
        "    flag = np.min(np.isclose(b_hat1, b_hat4))\n",
        "\n",
        "    print(\"Le même β est-il retourné? :\", flag)\n",
        "\n",
        "        \n",
        "test_lsqr()\n",
        "test_same_estimate()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6wT-6tVZn5uK"
      },
      "source": [
        "# The distribution of the MLE, the residuals, and the loss\n",
        "\n",
        "By computing 100 independent datasets with $p=10, n=30$, verify with histograms and\n",
        "QQ-plots that the squared norm $\\|y-X\\hat\\beta\\|^2/\\sigma^2$ has chi square distribution, and recall its degrees of freedom."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 90,
      "metadata": {
        "id": "tDAdgh5Hn5uL",
        "outputId": "f185796d-f316-4c04-a488-7d857ebec754",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        }
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAQHUlEQVR4nO3df4xlZX3H8fenrGhFK4tMtsgSFyuxQVMrmVCs1hAxitSwtDEEY+qqJBsTbLW2UayJ+KfUVqtNa7MV6toQhKIWYrS6pRjTP6AdEPmtLCiym4Udi4LWpLr67R/3bHsd7uzOvWfmzp3H9yu5uec85zlzvjl75rPPPPfec1NVSJLa8kvrXYAkafUZ7pLUIMNdkhpkuEtSgwx3SWrQpvUuAODEE0+sbdu2rXcZkrSh3Hrrrd+tqrlR22Yi3Ldt28bCwsJ6lyFJG0qSh5bb5rSMJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1aCY+oaoNJBl/H78QRpo6R+6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBh013JNcmeRgkruG2j6U5L4kdyT5XJLjh7a9N8neJN9I8pq1KvwXWjLZQ9IvjJWM3D8JnLukbQ/woqr6DeCbwHsBkpwOXAS8sNvnb5Mcs2rVSpJW5KjhXlVfBR5b0vblqjrUrd4MbO2WtwOfrqr/qapvAXuBM1exXknSCqzGnPtbgS92yycDDw9t29e1SZKmqFe4J3kfcAi4aoJ9dyZZSLKwuLjYpwxJ0hITh3uSNwOvA95Y9X+3/dsPnDLUbWvX9iRVtauq5qtqfm5ubtIyNAlfjJWaN1G4JzkXeDdwflX9aGjTDcBFSZ6a5FTgNOA/+pcpSRrHUe/nnuRq4GzgxCT7gMsYvDvmqcCeDEZ1N1fV26rq7iTXAvcwmK65pKp+ulbFS5JGS83AFynMz8/XwsLCepexcUw6TXL437rPF274ZR3SzEhya1XNj9rmJ1QlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KCj3hVSmgnesEwaiyN3SWqQ4S5JDTLcJalBzrlrepw3l6bGkbskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ06argnuTLJwSR3DbWdkGRPkvu7581de5J8LMneJHckOWMti5ckjbaSkfsngXOXtF0K3FhVpwE3dusArwVO6x47gY+vTpmSpHEcNdyr6qvAY0uatwO7u+XdwAVD7Z+qgZuB45OctFrFSpJWZtI59y1VdaBbfgTY0i2fDDw81G9f1/YkSXYmWUiysLi4OGEZkqRRer+gWlUFjH0DkKraVVXzVTU/NzfXtwxJ0pBJw/3Rw9Mt3fPBrn0/cMpQv61dmyRpiiYN9xuAHd3yDuD6ofY3de+aOQt4fGj6RpI0JUe95W+Sq4GzgROT7AMuAz4IXJvkYuAh4MKu+xeA84C9wI+At6xBzZKkozhquFfVG5bZdM6IvgVc0rcoSVI/fkJVkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGHfV97lojyfj71Ni38JH0C8qRuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3qFe5J/jjJ3UnuSnJ1kqclOTXJLUn2JrkmybGrVawkaWUmDvckJwN/BMxX1YuAY4CLgMuBj1TV84HvARevRqGSpJXrOy2zCfjlJJuApwMHgFcC13XbdwMX9DyGJGlME4d7Ve0H/gL4DoNQfxy4Ffh+VR3quu0DTh61f5KdSRaSLCwuLk5ahiRphD7TMpuB7cCpwHOA44BzV7p/Ve2qqvmqmp+bm5u0DEnSCH2mZV4FfKuqFqvqJ8BngZcBx3fTNABbgf09a5QkjalPuH8HOCvJ05MEOAe4B7gJeH3XZwdwfb8SJUnj6jPnfguDF05vA+7sftYu4D3Au5LsBZ4NXLEKdUqSxrDp6F2WV1WXAZctaX4QOLPPz5Uk9eMnVCWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJalCv97lLG0Iy2X5Vq1uHNEWO3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDeoV7kmOT3JdkvuS3JvkpUlOSLInyf3d8+bVKlaStDJ9R+4fBf6lqn4deDFwL3ApcGNVnQbc2K1LkqZo4nBP8izgFcAVAFX146r6PrAd2N112w1c0LdISdJ4+ozcTwUWgX9I8rUkn0hyHLClqg50fR4BtozaOcnOJAtJFhYXF3uUIUlaqk+4bwLOAD5eVS8B/pslUzBVVcDIbxmuql1VNV9V83Nzcz3KkCQt1Sfc9wH7quqWbv06BmH/aJKTALrng/1KlCSNa+Jwr6pHgIeTvKBrOge4B7gB2NG17QCu71WhJGlsm3ru/4fAVUmOBR4E3sLgP4xrk1wMPARc2PMYkqQx9Qr3qrodmB+x6Zw+P1eS1I+fUJWkBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWpQ31v+Su1Lxt+nRn4BmTQ1jtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBvcM9yTFJvpbk8936qUluSbI3yTVJju1fpiRpHKsxcn8HcO/Q+uXAR6rq+cD3gItX4RiSpDH0CvckW4HfBT7RrQd4JXBd12U3cEGfY0iSxtd35P5XwLuBn3Xrzwa+X1WHuvV9wMmjdkyyM8lCkoXFxcWeZayTZPyHJE3BxOGe5HXAwaq6dZL9q2pXVc1X1fzc3NykZUiSRuhzy9+XAecnOQ94GvArwEeB45Ns6kbvW4H9/cuUJI1j4pF7Vb23qrZW1TbgIuDfquqNwE3A67tuO4Dre1cpSRrLWrzP/T3Au5LsZTAHf8UaHEOSdASr8k1MVfUV4Cvd8oPAmavxcyVJk/ETqpLUIL9DVVpLfv+q1okjd0lqkOEuSQ0y3CWpQc65OycqqUGO3CWpQYa7JDXIcJekBhnuktQgX1CVZpUv9qsHR+6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJatDE4Z7klCQ3Jbknyd1J3tG1n5BkT5L7u+fNq1eupBVJJnuoGX1G7oeAP6mq04GzgEuSnA5cCtxYVacBN3brkqQpmjjcq+pAVd3WLf8AuBc4GdgO7O667QYu6FukJGk8qzLnnmQb8BLgFmBLVR3oNj0CbFlmn51JFpIsLC4urkYZkqRO73BP8gzgM8A7q+qJ4W1VVcDIbw+oql1VNV9V83Nzc33LkCQN6RXuSZ7CINivqqrPds2PJjmp234ScLBfiZKkcfV5t0yAK4B7q+rDQ5tuAHZ0yzuA6ycvT5I0iT7fofoy4A+AO5Pc3rX9GfBB4NokFwMPARf2K1GSNK6Jw72q/h1Y7o2x50z6cyVJ/fkJVUlqUJ9pGUktmvSTqjXyjXFaJ47cJalBG3/k7ihDkp7EkbskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDNv5dISXNlknu1OpdWledI3dJapDhLkkNMtwlqUGGuyQ1yBdUJc0OX4xdNY7cJalBaxbuSc5N8o0ke5NculbHkaR1l4z/WGNrEu5JjgH+BngtcDrwhiSnr8WxJElPtlYj9zOBvVX1YFX9GPg0sH2NjiVJWmKtXlA9GXh4aH0f8FvDHZLsBHZ2qz9M8o0j/LwTge+uaoV9/ixaft+j17k2xx1n/8nO5fTrHtS5/ufrSJY/l+tV9+h9V/5vvn51T/47PoUpjiGrl0WrU/dzl9uwbu+WqapdwK6V9E2yUFXza1xSbxuhzo1QI2yMOjdCjbAx6twINcLGqRPWblpmP3DK0PrWrk2SNAVrFe7/CZyW5NQkxwIXATes0bEkSUusybRMVR1K8nbgS8AxwJVVdXePH7mi6ZsZsBHq3Ag1wsaocyPUCBujzo1QI2ycOkn56S5Jao6fUJWkBhnuktSgmQn3JKckuSnJPUnuTvKOEX3OTvJ4ktu7x/vXqdZvJ7mzq2FhxPYk+Vh364U7kpwx5fpeMHSObk/yRJJ3LumzLucyyZVJDia5a6jthCR7ktzfPW9eZt8dXZ/7k+yYco0fSnJf9+/5uSTHL7PvEa+NKdT5gST7h/5dz1tm36ncHmSZGq8Zqu/bSW5fZt+pnMvlsmfWrsuxVdVMPICTgDO65WcC3wROX9LnbODzM1Drt4ETj7D9POCLQICzgFvWsdZjgEeA587CuQReAZwB3DXU9ufApd3ypcDlI/Y7AXiwe97cLW+eYo2vBjZ1y5ePqnEl18YU6vwA8KcruCYeAJ4HHAt8fenv2lrWuGT7XwLvX89zuVz2zNp1Oe5jZkbuVXWgqm7rln8A3Mvgk64b0XbgUzVwM3B8kpPWqZZzgAeq6qF1Ov7PqaqvAo8tad4O7O6WdwMXjNj1NcCeqnqsqr4H7AHOnVaNVfXlqjrUrd7M4LMb62qZc7kSU7s9yJFqTBLgQuDqtTj2Sh0he2bquhzXzIT7sCTbgJcAt4zY/NIkX0/yxSQvnGph/6+ALye5tbuNwlKjbr+wXv9RXcTyvzyzcC4BtlTVgW75EWDLiD6zdE7fyuAvs1GOdm1Mw9u76aMrl5lKmJVz+TvAo1V1/zLbp34ul2TPRrsuf87MhXuSZwCfAd5ZVU8s2Xwbg+mFFwN/DfzztOvrvLyqzmBw18tLkrxineo4ou4DZOcD/zRi86ycy59Tg791Z/b9uUneBxwCrlqmy3pfGx8Hfg34TeAAg2mPWfUGjjxqn+q5PFL2zPp1OcpMhXuSpzA4uVdV1WeXbq+qJ6rqh93yF4CnZHDDoamqqv3d80Hgcwz+zB02K7dfeC1wW1U9unTDrJzLzqOHp62654Mj+qz7OU3yZuB1wBu7X/YnWcG1saaq6tGq+mlV/Qz4+2WOPwvnchPw+8A1y/WZ5rlcJns2xHW5nJkJ927+7Qrg3qr68DJ9frXrR5IzGdT/X9OrEpIcl+SZh5cZvNB215JuNwBv6t41cxbw+NCfd9O07MhoFs7lkBuAw+8y2AFcP6LPl4BXJ9ncTTW8umubiiTnAu8Gzq+qHy3TZyXXxppa8trO7y1z/Fm4PcirgPuqat+ojdM8l0fInpm/Lo9ovV/RPfwAXs7gz547gNu7x3nA24C3dX3eDtzN4NX9m4HfXoc6n9cd/+tdLe/r2ofrDIMvK3kAuBOYX4c6j2MQ1s8aalv3c8ngP5sDwE8YzE9eDDwbuBG4H/hX4ISu7zzwiaF93wrs7R5vmXKNexnMrR6+Nv+u6/sc4AtHujamXOc/dtfcHQzC6aSldXbr5zF4V8gDa1nnqBq79k8evhaH+q7LuTxC9szUdTnuw9sPSFKDZmZaRpK0egx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KD/BXjLBgnCvUPlAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "source": [
        "def error_residual_plot():\n",
        "    p = 10\n",
        "    n = 30\n",
        "    sigma=0.7\n",
        "    data = list()\n",
        "    for i in range(1000):\n",
        "      X, y, b_true = generate_dataset_with_intercept(n, p, sigma, signal_strength=2)\n",
        "      b_hat = compute_least_squares_using_SVD(X, y)\n",
        "      data.append(np.linalg.norm(y - np.dot(X, b_hat))**2 / sigma**2)\n",
        "    plt.hist(data, 20, color = \"red\", width=0.8, align='mid')\n",
        "    \n",
        "error_residual_plot()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kKCs8TcNn5uM"
      },
      "source": [
        "# The curve of the risk of the MLE in high dimension\n",
        "\n",
        "1. Fix $n=200$, noise level $\\sigma^2=0.7$, signal strength $2$, and a grid of values for $p/n$, say $0.1, 0.2, 0.3, 0.4, ..., 2.9, 3.0$. For each for of $n, p$, generate $10$ independent datasets and plot the boxplots corresponding the squared error $\\|\\hat\\beta - \\beta^*\\|^2$ for each given ratio $p/n$. If $p>n$, make sure that $\\hat\\beta$ is the solution $(X^TX)^\\dagger X^Ty$. Plot the corresponding theoretical curves for the expected squared error found using the formula for expected inverse Wishart matrices.\n",
        "\n",
        "2. Bonus: Plot also the boxplots for the the bias and variance of the squared error. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 100,
      "metadata": {
        "id": "f2_Viox5n5uN",
        "outputId": "83121dff-9725-47a7-9cc8-c3caf096515d",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        }
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1376: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n",
            "  X = np.atleast_1d(X.T if isinstance(X, np.ndarray) else np.asarray(X))\n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "{'boxes': [<matplotlib.lines.Line2D at 0x7f01f32e5310>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3380110>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3583610>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3593b90>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f367c050>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36845d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34b3b50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f347c0d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37a2590>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f379fa10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3859f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f384c4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f375ba50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3737f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f371d4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f372ba10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3725f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35bc4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35c3a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f332af50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f32f74d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f31f8a50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f321ef50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36134d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35f8a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35fcf50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34db4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34c8a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f365ef90>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36444d0>],\n",
              " 'caps': [<matplotlib.lines.Line2D at 0x7f01f339bb50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f33780d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35980d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3598610>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35947d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3594cd0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f357fb90>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f357ffd0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36b4050>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36b4550>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34aa5d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34aab10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34a0b10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34ae090>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37b2090>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37b25d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37b5550>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37b5a90>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f383e290>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f383e850>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3842f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f385d4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f375f4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f375fa10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3766a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3766f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3769f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f376ff90>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37264d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3726a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f370ea10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f370ef50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35d8f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35be4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35b74d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35b7a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3315a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3315f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3320f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f32fa4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f331d4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f331da10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f320aa10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f320af50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3226f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f322f4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36094d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3609a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35fba10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35fbf50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35faf50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34be4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34cc4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34cca10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34c1a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34c1f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36439d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36524d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36624d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3662a10>],\n",
              " 'fliers': [<matplotlib.lines.Line2D at 0x7f01f3378b90>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35830d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35935d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36a4a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3684050>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34b35d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34aeb10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37a2090>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f379f510>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3859a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f385df50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f375b4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3737a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f376fa10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f372b4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3725a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35bef50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35c34d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f332aa10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f32faf50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f31f84d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f321ea10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f322ff50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35f84d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35fca10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34bef50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34c84d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f365ea10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3652f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36594d0>],\n",
              " 'means': [],\n",
              " 'medians': [<matplotlib.lines.Line2D at 0x7f01f3378650>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3598b50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3593110>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36a4f10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f36b4a90>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34b3090>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34ae5d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37b2b10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37b5fd0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f38594d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f385da10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f375ff50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37374d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f376f5d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3726f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37254d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35bea10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35b7f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f332a4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f32faa10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f331df50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f321e4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f322fa10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3609f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35fc4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34bea10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34ccf50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f365e4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3652a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3662f50>],\n",
              " 'whiskers': [<matplotlib.lines.Line2D at 0x7f01f32f4bd0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f339b610>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3380650>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3380b50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3583b50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3594110>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f357f110>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f357f650>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f367c5d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f367cad0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3684b10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34aa090>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34a0090>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34a05d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f347c5d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f347cb10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37a2a90>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37a2fd0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f379ff50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f383e550>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f38424d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3842a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f384ca10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f384cf50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f375bf50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37664d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f37694d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3769a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f371da10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f371df50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f372bf50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f370e4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35d84d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35d8a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35bca10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35bcf50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35c3f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f33154d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f33204d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3320a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f32f7a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f32f7f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f31f8f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f320a4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f32264d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3226a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3613a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3613f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35f8f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35fb4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35fa4d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f35faa10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34dba10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34dbf50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34c8f50>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f34c14d0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3643050>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3643cd0>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3644a10>,\n",
              "  <matplotlib.lines.Line2D at 0x7f01f3644f50>]}"
            ]
          },
          "metadata": {},
          "execution_count": 100
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAe6UlEQVR4nO3df5Ac5X3n8ffXy4LANrCL1linH5FT4Nz69hz/2AO7vHZYOK4wodA5hgPunDPJuFTnMzJOnMvZmSuQSO3FJDnHtsiZUrHEgC9LfMRxZB8+wMfaeKoOzIoIWWJtR06USEQ2MhKSKSJuJX/vj+5ZRjPdPT2/e3o+r6qpnen5zjPPzDPz3Wee5+luc3dERKT/varXFRARkfZQQhcRyQkldBGRnFBCFxHJCSV0EZGcOK1XT7xy5Upfv359r55eRKQv7dix4yfuPhZ1X+qEbmZDwALwrLtfVXXfGcC9wNuB54Hr3H1fUnnr169nYWEh7dOLiAhgZn8Xd18jQy43A4sx9xWAI+5+AfBHwO0NlCsiIm2QKqGb2Rrgl4G7YkI2APeE1x8ALjMza716IiKSVtoe+meA3wZ+FnP/amA/gLufAI4C51UHmdlGM1sws4VDhw41UV0REYlTN6Gb2VXAc+6+o9Unc/dt7j7p7pNjY5Fj+iIi0qQ0PfR3AVeb2T7gfuBSM/tiVcyzwFoAMzsNOIdgclRERLqkbkJ390+6+xp3Xw9cDzzq7h+oCtsOfDC8fk0Yo6N+iYh0UdPr0M3sNmDB3bcDs8B9ZrYXOEyQ+EVEpIsaSuju/k3gm+H1Wyq2HweubWfFRESkMT3bU1R6o3I1qUbFRPJFx3IZMOUkrmQukj9K6CIiOaEhF2mL6h2D9QtApPvUQ5e2qBzKUTIX6Q0ldBGRnFBCFxHJCSV0EZGcUEIXEckJJXQRkSbMzc0xMTHB0NAQExMTzM3N9bpKWrYoItKoubk5isUis7OzTE1NUSqVKBQKANxwww09q5d66CIiDZqZmWF2dpbp6WmGh4eZnp5mdnaWmZmZntbLerVmeHJy0nWS6N4ws46sFe9UuSJZMzQ0xPHjxxkeHl7etrS0xIoVKzh58mRHn9vMdrj7ZNR96qGLiDRofHycUql0yrZSqcT4+HjsY8zslEsnKKGLiDSoWCxSKBSYn59naWmJ+fl5CoUCxWIx9jGN7E3dbPLXpKiISIU0h5guT3xu2rSJxcVFxsfHmZmZaduEqLs3NYSpMfQBpDF0kWRZ+I7ExbY0hm5mK8zsO2b2tJntMbMtETE3mtkhM9sZXj6UqsYiItI2aYZcXgYudfcXzWwYKJnZ19398aq4P3P3m9pfRRERSaNuD90DL4Y3h8OLflfLstHR0eVxx/IkzujoaI9rJdJ5nV610qhUq1zMbMjMdgLPAY+4+xMRYe83s11m9oCZrY0pZ6OZLZjZwqFDh1qotmTJkSNHlmfuy5cjR470uloiQGeXC6Y5pWO5w9ONTk+qhO7uJ939LcAa4CIzm6gK+Sqw3t3fDDwC3BNTzjZ3n3T3ybGxsVbqLSKSSuUywV5M2kd1eDrV6WloHbq7vwDMA1dUbX/e3V8Ob94FvL091RMR6Y5u9qQ7pe6kqJmNAUvu/oKZnQlcDtxeFbPK3Q+GN68GFtteUxGRDjr80ZPA2RH3dHZX/mqjo6PLvffyP5eRkREOHz5c97FpVrmsAu4xsyGCHv2X3P1rZnYbsODu24GPmtnVwAngMHBjE69DRKRnbMuxuHXf+Obu1aM8RFNdhzS0Y9EAavdOE1HlaScjyZp6n8mEHXmaetyyzeck3Hc0VXmV25J2LNKu/yIiDaocFoHkoZFGe/5+69k1/wT81qihoFpK6CKSCdXDCln+hRc1LALph0aSRP0DSDvso4QuIpnQ7AGp+kFUoh8ZGWn78+jwuSKSW1lYili59rzydppVK41SD11EGpKFoZG0dWh0aCRtTzpqnHt5extU1yNtb14JXUQakoWhkYoVH22rR2U5dVe2dHCJYyuvTQldRDomC735QaKELiIdk4XefKd0a6KzEZoUFRlwc3NzTExMMDQ0xMTEBHNzc72uUtssj3VXXVod6+70RGf1JG5a6qGLDLC5uTmKxSKzs7NMTU1RKpUoFAoAbTs/Zi81M9ZdmUzbvUd15d+kspt9XvXQRQbYzMwMs7OzTE9PMzw8zPT0NLOzs8zMzPS6arGiliK2cxlide+7E+V2aghKx3IZQDqWi5QNDQ1x/PhxhoeHl7ctLS2xYsUKTp6MP8pgO052XKl6V3qIP8JgI5+3Zo/P0ohuf9ZbOkm0iOTX+Pg4pVLplG2lUonx8fGu1qOTZ72qPmORmbVt8rJ6GKXXlNBFBlixWKRQKDA/P8/S0hLz8/MUCgWKxWLLZWfhXLOdnrzs9BBKozQpKjLAyhOfmzZtYnFxkfHxcWZmZtoyIdrIcb1bOcKgvEJj6ANIY+jSqjTt2+pYd+xzxB1fPOLY4o3WuR+0dDx0M1sBPAacEcY/4O63VsWcAdxLcC7R54Hr3H1fi/UWkQxp5BjgndTK4WXzLs2Qy8vApe7+opkNAyUz+7q7P14RUwCOuPsFZnY9wTlHr+tAfUWkRzp5DHBpj7qToh54Mbw5HF6qW3UDcE94/QHgMlMri4h0VapJ0fAE0TuAC4A/dvcnqkJWA/sB3P2EmR0FzgN+UlXORmAjwLp161qruYhkmiY6uy9VQnf3k8BbzOxc4C/MbMLddzf6ZO6+DdgGwaRoo48Xkc6o/EEdN3HY6DHAOznW3ezxwvOuoWWL7v6Cmc0DVwCVCf1ZYC1wwMxOA84hmBwVkT6Q5qiInTwGeLmcSnFJupHjlleXm5eVLnHqjqGb2VjYM8fMzgQuB75XFbYd+GB4/RrgUc/zuyYibdWtnX/ynpbS9NBXAfeE4+ivAr7k7l8zs9uABXffDswC95nZXuAwcH3HaiwiIpHqJnR33wW8NWL7LRXXjwPXtrdqItJpja4t7+RJHRo5vKxE067/IgOskbXljY5dN0oJvHVK6CLSMVqN0l1K6AMiK7ttS/MaOeFy2thGlyI2opWz10tzlNAHhHbb7n+NJMi0sZ1eiijdpeOhi4jkhBK6iHRMs2evl+ZoyEVEOkZj592lHrpIH+j0me4lH9RDl5bpqHqd18jp3LK0s1CnDNLxWRqhHrq0zLYcC07/VXGxLcdi4+fm5piYmGBoaIiJiQnm5ua6WNv8Kyf/6ktlki+LOr5JO4+j0imDdHyWRqiHLl01NzdHsVhkdnaWqakpSqUShUIBoC0nJs6rRn4FNbO2XEkxH3SS6AHRyMl521F2XLkTExNs3bqV6enp5W3z8/Ns2rSJ3bsbPsT+wGj1hMtJ25upSyX9M+iupJNEK6EPiKwk9KGhIY4fP87w8PDytqWlJVasWMHJkydbqkcndSqJpS03SwldeispoWsMXbpqfHycUql0yrZSqcT4+HiPapRO9bG621luZdkirVBCl64qFosUCgXm5+dZWlpifn6eQqFAsVjsddUyr7xssXzJ+koU6T5NikpXlSc+N23axOLiIuPj48zMzAzkhGjU8sK4pYWdPnSt5IPG0AdEVsbQ+1m7X1Oz71vdA25pDD3XksbQ6/bQzWwtcC9wPuDANnf/bFXMJcBfAn8bbvqyu9/WSqVF+lUrh65t1w5Z/bizkLQuzZDLCeDj7v6Umb0W2GFmj7j7M1Vx33b3q9pfRZH+4u7petoRh65tx2FrNTwzuOpOirr7QXd/Krz+U2ARWN3piolId2jP3fxoaFLUzNYTnDD6iYi732lmTwP/APyWu++JePxGYCPAunXrGq2riLSZ9tzNl9STomb2GuBbwIy7f7nqvrOBn7n7i2Z2JfBZd78wqTxNinaXJkWbU70Spazeqfuandxsx6RoI7Hac7f/tLynqJkNA18DHnL3T6eI3wdMuvtP4mKU0LtLCb05jb5vUf8Ako5y2Mj71szeqvXaoV/33B1kra5yMWAWWIxL5mb2euDH7u5mdhHB2PzzLdRZ2qyTJwMGnd29rJHD3DaqkX+QaQ8vW95zt7KH3g977kq0NGPo7wJ+Ffiume0Mt/0OsA7A3e8ErgE+bGYngH8Ervc8dM9ypJMnA9bZ3bMnbTuU99ytHkOfmZnpcA2lE+omdHcvAYldDHe/A7ijXZUSke7Qnrv5oj1FB0Q39h7MYw+9HUc0bEesSFlLY+gi7Zbn42nrdHzSS0ro0nV5HnNvdO9PTSZLOymhi/SIdtGXdlNCF0nQ6eWeIu2khC4Dq3K4I3Znng4u9xRpNyV0yY1GJ1vTHhWx1XpoXFy6RQldciMLk61ZqIMMLp1TVEQkJ5TQJdOqT4zcDqOjo6eUV74+OjralvIbUV0PkVZoyEUyrRNDGFEH0Co/R7dpWEbaST10EZGcUA9dBo7WlkteKaHLwNHacskrDblILmRpolOkV5TQpauiEm87km55orP6EnU+UJG80pCLdFUnT9HWKVH1S9r7s3oZolaySLfU7aGb2VozmzezZ8xsj5ndHBFjZvY5M9trZrvM7G2dqa5Id1X29itvR530OeoxSubSTWmGXE4AH3f3NwHvAD5iZm+qinkvcGF42Qh8vq21lIHUqeEZkbxKc07Rg8DB8PpPzWwRWA08UxG2Abg3PDH042Z2rpmtCh8r0pR+HJ4R6aWGJkXNbD3wVuCJqrtWA/srbh8It1U/fqOZLZjZwqFDhxqrqWRa9eoSEem+1AndzF4D/DnwMXc/1syTufs2d59098mxsbFmipCM0rixSO+lWuViZsMEyfx/uPuXI0KeBdZW3F4TbhPpikb3/mx05YpIP6ib0C345M8Ci+7+6Ziw7cBNZnY/cDFwVOPn0k2N7P2pc3lKXqXpob8L+FXgu2a2M9z2O8A6AHe/E3gQuBLYC7wE/Fr7qyqSTL1uGXRpVrmUgMRZrnB1y0faVSnJr6ihkXYcFEu9bhHtKSpdFjU0ooNiibSHjuUiIpIT6qFLZnVqeEYkr5TQJbNsS+3uDiMjIxze3IO66IBb0geU0CWzsjTRqQQu/UBj6DKwqnvdIv1OPXTpuuoE2qu14up1S94ooUtXdXIYpfIfRa+HaER6QQldckMJXAadxtBFRHJCCV1EJCc05DJAdPAqkXxTQh8QWVrTLSKdoSEXEZGcUEIXEckJJXQRkZyom9DN7G4ze87Mdsfcf4mZHTWzneHllvZXU0RE6kkzKfoF4A7g3oSYb7v7VW2pkYiINKVuD93dHwMOd6EuMiDMLPWBsRqJFRl07RpDf6eZPW1mXzezfxYXZGYbzWzBzBYOHTrUpqeWfuPup1zaFSsy6NqR0J8Cfs7dfxHYCnwlLtDdt7n7pLtPjo2NteGpRUSkrOWE7u7H3P3F8PqDwLCZrWy5ZiIi0pCWE7qZvd7CwU0zuygs8/lWy5XO0Hi0SH7VXeViZnPAJcBKMzsA3AoMA7j7ncA1wIfN7ATwj8D1rsHOzFLTiORX3YTu7jfUuf8OgmWNIiLSQ9pTVEQkJ5TQRURyQgldRCQnlNBFRHJCJ7jIgcoliFrFIjK41EPPgXISVzIXGWxK6CIiOaGELiKSE0roIiI5oYQuIpITSugiIjmhhC4ikhNK6CIiOaGE3sdGR0drzrlpZoyOjva4ZiLSC9pTtI8dOXIkcmcinbxCZDCphy4ikhNK6CIiOVE3oZvZ3Wb2nJntjrnfzOxzZrbXzHaZ2dvaX00REaknTQ/9C8AVCfe/F7gwvGwEPt96tUREpFF1E7q7PwYcTgjZANzrgceBc81sVbsqKCIi6bRjDH01sL/i9oFwWw0z22hmC2a2cOjQoTY8tYiIlHV1UtTdt7n7pLtPjo2NdfOpRURyrx3r0J8F1lbcXhNukw7zW8+GzedEbxeRgdOOhL4duMnM7gcuBo66+8E2lCt12JZjsTsW+ebu10dEeqtuQjezOeASYKWZHQBuBYYB3P1O4EHgSmAv8BLwa52qrIiIxKub0N39hjr3O/CRttVIRESaoj1FRURyQgfn6nNRB+IaGRnpQU1EpNeU0PtY5YSomUVOkIrI4NCQi4hITiihZ9Dc3BwTExMMDQ0xMTHB3Nxcr6skIn1AQy4ZMzc3R7FYZHZ2lqmpKUqlEoVCAYAbbkhccCQiA0499IyZmZlhdnaW6elphoeHmZ6eZnZ2lpmZmV5XTUQyzno1kTY5OekLCws9ee4sGxoa4vjx4wwPDy9vW1paYsWKFZw8eTL2cZoUFRkMZrbD3Sej7lMPPWPGx8cplUqnbCuVSoyPj/eoRiLSL5TQM6ZYLFIoFJifn2dpaYn5+XkKhQLFYrHXVRORjNOkaMaUJz43bdrE4uIi4+PjzMzMJE6Ilncu0rCLyGDTGLqISB/RGLqIyADQkEsGVR+fRcMoIpKGeugZ5O7LSVzJXETSUkLvEu3OLyKdliqhm9kVZvZ9M9trZp+IuP9GMztkZjvDy4faX9X+Vd6df+vWrRw/fpytW7dSLBYjk/ro6ChmdsrKFTNjdHS029UWkT5Td5WLmQ0BPwAuBw4ATwI3uPszFTE3ApPuflPaJx6kVS4TExNs3bqV6enp5W3z8/Ns2rSJ3bt3nxIbt/RQSxJFBJJXuaSZFL0I2OvufxMWdj+wAXgm8VGybHFxkampqVO2TU1Nsbi4WBPrt54Nm8+J3i4ikiBNQl8N7K+4fQC4OCLu/Wb2HoLe/G+4+/7qADPbCGwEWLduXeO17VPl3fkre+hxu/PblmPxPfTNnayliPS7dk2KfhVY7+5vBh4B7okKcvdt7j7p7pNjY2Nteursa3R3/vK4eeVFp5UTkXrS9NCfBdZW3F4Tblvm7s9X3LwL+P3Wq5YfjezOr9PKiUiz0kyKnkYwjHIZQSJ/Evi37r6nImaVux8Mr78P+M/u/o6kcgdpUhSa21lICV1EqrU0KeruJ8zsJuAhYAi42933mNltwIK7bwc+amZXAyeAw8CNbat9DkQtORwdHeXw4cOR8ZXJX0ldRNLSwbm6ICopK1GLSDNaXbYoLYpaiqhliCLSbkroXWBbjtVsGxkZ4fDm7tdFRPJLCb0LtHJFRLpBB+fqkurjs4iItJt66F2iXrmIdJoSegt0IgoRyRINuTQpbm25iEivqIfepCNHjkSuLRcR6RX10EVEckI99CZpZyERyRol9GZtPgpoXbmIZIeGXEKV5/KsvGiiU0T6hXroocMfPQlEDZmcjH1M9Y5C6qmLSC/lOqGPjo5y5MiRmu0jIyM1h66NOt7Kcuzm6PKVwEUkS3Kd0BvpdZeTc/XSw7hjlouIZE3fJvTKxBvXU1avW0QGSd9NipYnLyvFTV66e02Cdnf1ukUkl1L10M3sCuCzBKegu8vdP1V1/xnAvcDbgeeB69x9X3urGmhm8lK9bhEZBHUTupkNAX8MXA4cAJ40s+3u/kxFWAE44u4XmNn1wO3AdWkr0enJSxGRQZCmh34RsNfd/wbAzO4HNgCVCX0DsDm8/gBwh5mZp+waNzN5KSIip0qT0FcD+ytuHwAujotx9xNmdhQ4D/hJZZCZbQQ2Aqxbt+6VO8K9LkVEpHldnRR1923uPunuk2NjY918ahGR3EuT0J8F1lbcXhNui4wxs9OAcwgmR0VEpEvSJPQngQvN7A1mdjpwPbC9KmY78MHw+jXAo2nHz0VEpD3qjqGHY+I3AQ8RLFu82933mNltwIK7bwdmgfvMbC9wmCDpi4hIF6Vah+7uDwIPVm27peL6ceDa9lZNREQa0Xd7ioqISDQldBGRnFBCFxHJCevVYhQzOwT8XcRdK6naISlBv8VmpR5ZiM1KPbIQm5V69FtsVurR7difc/foHXnKRyTMyoVg5UwuY7NSjyzEZqUeWYjNSj36LTYr9chCbPmiIRcRkZxQQhcRyYksJvRtOY7NSj2yEJuVemQhNiv16LfYrNQjC7FADydFRUSkvbLYQxcRkSYooYuI5EWjy2I6dQHuBp4DdqeIXQvME5w1aQ9wc0LsCuA7wNNh7JYU5Q8BfwV8rU7cPuC7wE7qLDECziU4m9P3gEXgnTFxvxCWV74cAz6WUO5vhK9rNzAHrEiIvTmM2xNVZlQbAKPAI8Bfh39HEmKvDcv+GTBZp9w/CN+LXcBfAOfWif/dMHYn8DDwT+p9boCPAw6sTCh3M8Hhn8vv95VJ5QKbwnrvAX4/odw/qyhzH7AzIfYtwOPlzxFwUULsLwL/N/zcfRU4O+k7EdV+CbFx7RcXX9OGCbE17RcXG9V+CeXWtF9SudXtl1BuXPvFxde0YUJsTRsSk6eANwBPAHvDOp2emGcaTbydugDvAd5GuoS+CnhbeP21wA+AN8XEGvCa8Ppw+Oa8o075vwn8KekS+sqUr+8e4EPh9dOpSGAJjxkCfkSwI0HU/auBvwXODG9/CbgxJnaCIJmfRXBQtm8AF9Rrg/BD/4nw+ieA2xNixwn+IX2TUxNCVOy/Ak4Lr99eLjch/uyK6x8F7kz63IRfpocIdl5bmVDuZuC30nwegenwfTsjvP26NJ9d4L8BtySU+zDw3vD6lcA3E2KfBH4pvP7rwO8mfSei2i8hNq794uJr2jAhtqb94mKj2i+h3Jr2S4itab+kOsS0X1zZNW2YEFvThsTkKYLv9PXh9juBDyfljMwMubj7YwSH3k0Te9Ddnwqv/5Sgx7s6Jtbd/cXw5nB4iZ0JNrM1wC8Dd6WvfTIzO4fgyzkb1un/ufsLKR56GfBDd4/ao7bsNODM8MQiZwH/EBM3Djzh7i+5+wngW8CvVAbEtMEGgn9GhH//dVysuy+6+/ernzgm9uGwHhD0bNbUia88O/irCdsw4XPzR8BvU9HWDX7GomI/DHzK3V8OY56rV66ZGfBvCH49xcU6r5xU9xzCNoyJfSPwWHj9EeD9YWzcd6Km/eJiE9ovLr6mDRNia9qvzvf4lPZr8DsfF1vTfvXKjWi/uPiaNkyIrWnDhDx1KcEve6j4/sXJTEJvlpmtB95K8B8tLmbIzHYS/Hx9xN1jY4HPEHyQfpbi6R142Mx2hOdLjfMG4BDwJ2b2V2Z2l5m9OkX51xN+kCKf3P1Z4A+BvwcOAkfd/eGY8N3Au83sPDM7i1d+mtZzvrsfDK//CDg/xWMa9evA1+sFmdmMme0H/h1wS0LcBuBZd3865fPfZGa7zOxuMxtJiHsjwXv4hJl9y8z+RYqy3w382N3/OiHmY8AfhK/tD4FPJsTuIUjSEAyR1LRh1Xcisf3SfH9Sxte0YXVsUvtVxtZrv4g6xLZfVWxi+8W8ttj2q4pPbMOq2Mg2rM5TwA+BFyr+aR4g5p/YsqTue7cvwHpSDLlUxL8G2AH8Ssr4cwnGtCZi7r8K+O/h9UuoP+Sy2l/56fY08J6YuEngBHBxePuzhD+VE8o+neA4DucnxIwAjwJjBP/RvwJ8ICG+EL5fjwGfBz5Trw3CD1Tl/UfqtRdVP9nrxBYJxl8t7WeB4MuyJSqW4FfKE8A54e19VAyLRby+8wmGtl4FzBCcwCUudjewleDn8UUEw11W5/V9Hvh4nff4cwS9NAh6g99IiP2nBD/vdwC3As8nfSfqtF/k9yeq/erE17RhXGxM+y3Hpmi/6teX1H7VsUntF/faatovpuykNqyOrdeG5Tw1Beyt2L426jN2ymOT7uz2Je5LERM7TDDG9psNPsctRIyZhvf9HsF/wX0EvZmXgC+mLHdzQrmvB/ZV3H438L/qlLcBeLhOzLXAbMXtf0/4DylFff8r8B/rtQHwfWBVeH0V8P167UXKhA7cSDA5dFYjnwVgXVUdl2OBf07Qw9kXXk4Q/IJ5fYpyq1979e3/DUxX3P4hMJbw+k4DfkwwFJH0PEd5JbEYcCzl+/BG4DtJ34m49ouKrdN+kfFRbZhUdnX7VccmtV+KctfHlZvUfgmvLa79osqObMMUdT6lDSu23wL8J4JOXXme4p3AQ0nf674ccgnHtWaBRXf/dJ3YMTM7N7x+JnA5wSx3DXf/pLuvcff1BMMdj7r7B2LKfbWZvbZ8nWCCaHdMuT8C9pvZL4SbLiOY+U5yAwnDLaG/B95hZmeF78llBON0kczsdeHfdQS9oT+tUz6cer7YDwJ/meIxdZnZFQRDW1e7+0sp4i+suLmB+Db8rru/zt3Xh+14gGBi6kcx5a6quPk+Ytow9BWCiTXM7I288isqzr8EvufuBxJiIBgz/6Xw+qUEK1IiVbThq4D/QjBRlvSdqGm/Rr4/SWVHtWFCbE37RcXGtR9BYo0qt6b9El5fXPvFvRc17ZdQdk0bJrwXNW0Yk6cWCXrq14QPrf/9S8r23bwQJK+DwBJBIxYSYqcIxq/Ly6CWl5tFxL6ZYAniLoIv6y0p63MJCUMuwM8TDLOUlxkV65T3FoLlTLsIPlgjCbGvBp4n/NlZp9wtBMltN3Af4Qx+TOy3Cf6RPA1clqYNgPOA/0OQZL4BjCbEvi+8/jLBF/ChhNi9wP6K9ruzTj3+PHyNuwiWeq2Oi616Tft4ZZVLVLn3ESwf20WQ/FYlxJ4OfDGsx1PApUl1AL4A/IcU7/EUwc/vpwmGG96eEHszwWqJHwCf4pVeYeR3Iqr9EmLj2i8uvqYNE2Jr2i8uNqr9Esqtab+E2Jr2S6pDTPvFlV3ThgmxNW1ITJ4iyDPfCd/r/0nC99vdteu/iEhe9OWQi4iI1FJCFxHJCSV0EZGcUEIXEckJJXQRkZxQQhcRyQkldBGRnPj/YzfLLgRgV+AAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "source": [
        "# Question 1\n",
        "n = 200\n",
        "sigma = m.sqrt(0.7)\n",
        "signal_strength = 2\n",
        "p_values = [int(0.1*i*n) for i in range(1,31)]\n",
        "data = [list() for p in p_values]\n",
        "\n",
        "for k in range(len(p_values)):\n",
        "  for i in range(10):\n",
        "    if (p_values[k]/n == 1):\n",
        "      continue\n",
        "    X, y, b_true = generate_dataset_with_intercept(n, p_values[k], sigma**2, signal_strength)\n",
        "    b_hat = compute_least_squares_using_SVD(X, y)\n",
        "    data[k].append(np.linalg.norm(b_hat - b_true)**2)\n",
        "plt.boxplot(data)\n",
        "# Question 2\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "putpCxQEn5uO"
      },
      "outputs": [],
      "source": [
        ""
      ]
    }
  ],
  "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.8.10"
    },
    "colab": {
      "name": "TP1 Numpy and Maximum Likelihood in linear regression.ipynb",
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}