TP1 Numpy and Maximum Likelihood in linear regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "import scipy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9147718188549299"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.normal()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 8.65696056, 11.66274356, 10.88337469])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu = 10\n",
    "variance = 2\n",
    "np.random.normal(mu, variance, size=(3, ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8, 1)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 0, 1, 1])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.binomial(1, p, size=(5, ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_841398/1429000599.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     22\u001b[0m     \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\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[0mgenerate_wishart\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meye\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m3\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/tmp/ipykernel_841398/1429000599.py\u001b[0m in \u001b[0;36mgenerate_wishart\u001b[0;34m(n, covariance, k)\u001b[0m\n\u001b[1;32m     20\u001b[0m         \u001b[0mk\u001b[0m \u001b[0mWishart\u001b[0m \u001b[0mmatrices\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     21\u001b[0m     \"\"\"\n\u001b[0;32m---> 22\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     23\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     24\u001b[0m \u001b[0mgenerate_wishart\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meye\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m3\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": [
    "from scipy import linalg\n",
    "\n",
    "def generate_wishart(n, covariance, k=1):\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",
    "    raise NotImplementedError\n",
    "    \n",
    "generate_wishart(5, np.eye(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.choice(6) + 1 # roll a die"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5, 2, 3, 1, 2, 3, 2, 5, 2, 3])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 5, 6, 5, 4, 6, 6, 6, 2, 6])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 19,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_841398/360823744.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     13\u001b[0m     \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 15\u001b[0;31m \u001b[0mbiased_die_histogram\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m]\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[0m\u001b[1;32m     16\u001b[0m \u001b[0mbiased_die_histogram\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\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/tmp/ipykernel_841398/360823744.py\u001b[0m in \u001b[0;36mbiased_die_histogram\u001b[0;34m(p, k)\u001b[0m\n\u001b[1;32m     11\u001b[0m         \u001b[0mNumber\u001b[0m \u001b[0mof\u001b[0m \u001b[0mrolls\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m     \"\"\"\n\u001b[0;32m---> 13\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     14\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0mbiased_die_histogram\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m]\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;31mNotImplementedError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "def biased_die_histogram(p, k):\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",
    "    raise NotImplementedError\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": {},
   "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": 20,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_841398/3220792203.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0mtest_singular_values_eigenvalues\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[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mtest_pinv\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;32m/tmp/ipykernel_841398/3220792203.py\u001b[0m in \u001b[0;36mtest_singular_values_eigenvalues\u001b[0;34m()\u001b[0m\n\u001b[1;32m      2\u001b[0m     \"\"\"Test that the square of the singular values of X gives the eigenvalue of X^T.dot(X)\n\u001b[1;32m      3\u001b[0m     \"\"\"\n\u001b[0;32m----> 4\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      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mtest_singular_values_eigenvalues\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 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",
    "    raise NotImplementedError\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",
    "    raise NotImplementedError\n",
    "\n",
    "test_pinv()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_841398/1688407605.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      9\u001b[0m     \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m \u001b[0mplot_boxplot_normal_increasing_mean_and_variance\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/tmp/ipykernel_841398/1688407605.py\u001b[0m in \u001b[0;36mplot_boxplot_normal_increasing_mean_and_variance\u001b[0;34m()\u001b[0m\n\u001b[1;32m      7\u001b[0m     \u001b[0mthe\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0maxis\u001b[0m \u001b[0mrepresenting\u001b[0m \u001b[0mtheta\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m     \"\"\"\n\u001b[0;32m----> 9\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     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0mplot_boxplot_normal_increasing_mean_and_variance\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 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",
    "    raise NotImplementedError\n",
    "\n",
    "plot_boxplot_normal_increasing_mean_and_variance()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": 22,
   "metadata": {},
   "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": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkkElEQVR4nO3df5xcdX3v8dd7ZrIb8juQ8CsEEiBYoraAMdZH1fahqIG2BCtWqLRYf3Bt5dqWaosXL/qgt63Va39KVay5/mgRf7ZNa7wUq2i9imZBBAIiIQSS8GtJwo9kk92dmc/945zZnN3M7Jzd7O7sWd7Px2Mfe+Z7fsx3z0ze883nnDlHEYGZmc1cpU53wMzMJpeD3sxshnPQm5nNcA56M7MZzkFvZjbDOejNzGY4B70dRtIHJP1jh577Fklvyzz+X5KelPTYOLb1C5Lul7RP0oUT2tEjIGm7pHMnadsh6fTJ2Ha6/ZdLum+ytm+Tw0H/HJQGX+OnLulA5vGbOt2/BkknA38IrI6I48exiWuBj0bEvIj4lxbP8RuSetK//VFJX5f0siPo9owWEf8VEc9rPJ7MDy2bOA7656A0+OZFxDzgYeBXM23/1On+ZZwM7I6IJ8a5/inAllYzJV0J/DXwZ8Bx6fP9PbB+nM9nNi056K2VLkmflfSspC2S1jRmSDpR0lck9Up6UNK7Wm1E0qclfVzSzem2vi3plMz8V0v6iaSnJX0UUNp+LnAzcGI62v50i+2/XdJWSXskbZR0Ytr+AHAq8G/p+t0j1ltIMuJ/Z0R8NSL2R8RgRPxbRLwnXWatpO9Leiod7X9UUlc6b0VaJqlktjmy7PR2Sfemf/c9ks7JdOEsSXemf/cXJM1O17lb0q9mtjErLV2d3eLvf0/at0ckvWXEvG5J/1vSw5IeT1+Ho9J5vyRpp6Q/lPREuo3fzqx7ftrnZyXtkvTu7Hrp9OdIPhwb+/iPJH1N0n8f0Y87Jb2uWf9tikSEf57DP8B24NwRbR8ADgLnA2Xgz4Fb03kl4DbgGqCLJEy3Aa9tsf1PA88CrwC6gb8BvpvOW5LOuwiYBfwBUAXels7/JWDnKH1/JfAkcE667b8DvjPa35aZty59rsoo238R8PNABVgB3Av8fjpvBRDZ9YFbMn1/A7ALeDHJh9fpwCmZfv0QOBE4Ot3uO9J5fwR8IbPN9cBdo/wNjwMvAOYCN6R9Oj2d/1fAxvQ55gP/Bvx5Zt9WST7sZqWvdR+wOJ3/KPDydHoxcE6z12TkPgZ+HfhB5vHPAbuBrk6/15/LPx7RWyvfjYhNEVEDPkfyDxaS4FoaEddGxEBEbAM+CVw8yra+FhHfiYh+4GrgpZKWk4TLloj4ckQMkpRRxnLQ9U3Ahoi4Pd32e9Ntr8ix7jHAkxFRbbVARNwWEbdGRDUitgOfAH4xZ9/eBnwoIjZHYmtEPJSZ/7cR8UhE7CEJ4LPS9n8Ezpe0IH38myT7v5lfB/5PRNwdEftJPqABkCTgcuAPImJPRDxLUqLKvk6DwLWR/E9mE7APeF5m3mpJCyJib0TcnvPv3gicIWlVpv9fiIiBnOvbJHDQWyvZwO0DZqdlilNIyilPNX6A/0FS425lR2MiIvYBe0hGsyeOmBfZxzmcCAyFZ7rt3cCyHOvuBpZkSy8jSTpD0r9LekzSMyRBuSRn35YDD4wyf+T+nQcQEY8A/w94vaRFwHlAq+Mmw/YfmX0BLAXmALdlXqf/m7Y37B7xQTfUD+D1JB/ED6XltpeO8rcMiYiDwBeASyWVgEto/UFlU8RBb2O1A3gwIhZlfuZHxPmjrLO8MSFpHkkp4RGS8kB2nrKPc3iE5IOnsf5ckpH6rhzrfh/oBy4cZZmPAT8BVkXEApIPNKXz9qe/52SWz54ZtAM4LUc/mvkMcClJ+ef7EdHq7xm2/0jq5Q1PAgeA52dep4WRHIBvK/2fyHrgWOBfgC+2WrRF/98EvAroi4jv53lOmzwOehurHwLPSvpjSUdJKkt6gaQXj7LO+ZJelh7I/BOSev8O4GvA8yX9WjqyfhfDw7KdzwO/Lems9GDrn5HUh7e3WzEiniY5znCdpAslzUkPfJ4n6UPpYvOBZ4B9kn4G+J3M+r0kHyiXpvvgLQwP9n8A3i3pRUqcnj0I3ca/kBx3+D3gs6Ms90XgzZJWS5oDvD/TvzpJSe2vJB0LIGmZpNe2e3JJXZLeJGlhWlJ7Bqi3WPxxkuM0Q9JgrwMfwaP5acFBb2OS1ux/haSm/CDJyPEfgIWjrHYDSQjtITnAeWm6rSdJRq0fJCmlrCIpW+TtyzeA/wl8hWR0exqjHysYuf5HgCuB9wG9JKPwK0iCFuDdwG+QHDD+JElJIuvtwHvSvj8f+F5m218C/pTkb3823ebROft1IP2bVgJfHWW5r5Mc1/gmsDX9nfXHafutaenpGxyqwbfzm8D2dL13kIzQm/lz4H1peejdmfbPAi8kOeZgHaakLGo2OZScFrkzIt7X6b4UiaRrgDMi4tJO92U8JP0WcHlE+Mtn00DLA1Fm1hmSjgbeSjKqLpy0jPS7JF8+s2nApRuzaUTS20lKSF+PiO90uj9jlR4D6CWp3d/Q4e5YyqUbM7MZziN6M7MZbtrV6JcsWRIrVqzodDfMzArltttuezIiljabN+2CfsWKFfT09HS6G2ZmhSLpoVbzXLoxM5vhHPRmZjOcg97MbIZz0JuZzXAOejOzGc5Bb2Y2wznozcxmOAe9WU6bt+/hJ4890+lumI3ZtPvClNl09YaPJzdK2v7BX+5wT8zGxiN6M7MZzkFvZjbDOejNzGY4B71ZDvW679tgxeWgN8uhv1rvdBfMxs1Bb5bDgcFap7tgNm4OerMcHPRWZLmCXtI6SfdJ2irpqibz3yHpLkl3SPqupNWZee9N17svvXGwWeEMunRjBdY26CWVgeuA84DVwCXZIE/dEBEvjIizgA8Bf5muuxq4GHg+sA74+3R7ZoVS9cFYK7A8I/q1wNaI2BYRA8CNwPrsAhGR/V74XKDxr2I9cGNE9EfEg8DWdHtmhVKPQ0Ef4dC3YslzCYRlwI7M453AS0YuJOmdwJVAF/DKzLq3jlh32bh6atZB1dqhcB+sBV0VdbA3ZmMzYQdjI+K6iDgN+GPgfWNZV9Llknok9fT29k5Ul8wmTK2eDXrX661Y8gT9LmB55vFJaVsrNwIXjmXdiLg+ItZExJqlS5fm6JLZ1KplyjXZ0b1ZEeQJ+s3AKkkrJXWRHFzdmF1A0qrMw18G7k+nNwIXS+qWtBJYBfzwyLttNrVq9UOj+AGP6K1g2tboI6Iq6QrgJqAMbIiILZKuBXoiYiNwhaRzgUFgL3BZuu4WSV8E7gGqwDsjwickW+FkR/HVuoPeiiXX9egjYhOwaUTbNZnp3xtl3T8F/nS8HTSbDrKlm8GqSzdWLP5mrFkO2YOxNZ9eaQXjoDfLIfuFqZpLN1YwDnqzHLKXKfa3ZK1oHPRmOQwf0TvorVgc9GY51Bz0VmAOerMcai7dWIE56M1y8IjeisxBb5aDa/RWZA56sxzqDnorMAe9WQ5V1+itwBz0ZjlkvyTlL0xZ0TjozXIYfjC2gx0xGwcHvVkOvgSCFZmD3iwHn0dvReagN8she8VKn3VjReOgN8uhVnPQW3E56M1y8OmVVmQOerMc6i7dWIE56M1y8IjeisxBb5ZDdhRfd9BbwTjozXLw6ZVWZA56sxxq9aBSUjrtL0xZsTjozXKo1ut0VUrptEf0Viy5gl7SOkn3Sdoq6aom86+UdI+kOyX9p6RTMvNqku5IfzZOZOfNpkqtzlDQu0ZvRVNpt4CkMnAd8GpgJ7BZ0saIuCez2I+ANRHRJ+l3gA8Bb0znHYiIsya222ZTq1av01X2iN6KKc+Ifi2wNSK2RcQAcCOwPrtARHwrIvrSh7cCJ01sN806q1aHWeUSks+jt+LJE/TLgB2ZxzvTtlbeCnw983i2pB5Jt0q6sNkKki5Pl+np7e3N0SWzqVWPoFSCSkke0VvhtC3djIWkS4E1wC9mmk+JiF2STgW+KemuiHggu15EXA9cD7BmzRr/K7JpJznrpkS5JNforXDyjOh3Acszj09K24aRdC5wNXBBRPQ32iNiV/p7G3ALcPYR9NesI2r1oCSolEoe0Vvh5An6zcAqSSsldQEXA8POnpF0NvAJkpB/ItO+WFJ3Or0E+AUgexDXrBBq9aBcEiXX6K2A2pZuIqIq6QrgJqAMbIiILZKuBXoiYiPwYWAe8CVJAA9HxAXAmcAnJNVJPlQ+OOJsHbNCqEVQkqiUSw56K5xcNfqI2ARsGtF2TWb63BbrfQ944ZF00Gw6qNWDSlmUfTDWCsjfjDXLoVYPyhKVknwJBCscB71ZDsnplaIkj+iteBz0Zjk0LmpWKfv0SiseB71ZDtV6cjDWNXorIge9WQ719PTKpEbvoLdicdCb5VCLxnn0HtFb8TjozXJofGHKNXorIge9WQ6N0yvLvgSCFZCD3iyHWj05vdI1eisiB71ZDvVITq8sS1T9hSkrGAe9WQ7VdESfXKa4070xGxsHvVkO9cYlEMoe0VvxOOjNcqg1Sjeu0VsBOejNcqjV0tKNz6O3AnLQm+VQi8bplR7RW/E46M1yqNVJTq8sO+iteBz0ZjkMnV5Z8h2mrHgc9GY5VGt1yiVRFq7RW+E46M1yqAfpZYo9orficdCb5dC4Z6wvgWBF5KA3y6GW3nik5BuPWAE56M1ySK5Hj28OboWUK+glrZN0n6Stkq5qMv9KSfdIulPSf0o6JTPvMkn3pz+XTWTnzabKocsUu3RjxdM26CWVgeuA84DVwCWSVo9Y7EfAmoj4WeDLwIfSdY8G3g+8BFgLvF/S4onrvtnka9xopFwquUZvhZRnRL8W2BoR2yJiALgRWJ9dICK+FRF96cNbgZPS6dcCN0fEnojYC9wMrJuYrptNjepQ0OObg1sh5Qn6ZcCOzOOdaVsrbwW+PpZ1JV0uqUdST29vb44umU2deiTBXvJFzaygJvRgrKRLgTXAh8eyXkRcHxFrImLN0qVLJ7JLZkesEeyVxh2mwkFvxZIn6HcByzOPT0rbhpF0LnA1cEFE9I9lXbPprFGqaXxhKgLfINwKJU/QbwZWSVopqQu4GNiYXUDS2cAnSEL+icysm4DXSFqcHoR9TdpmVhiHDsaKcvovxnV6K5JKuwUioirpCpKALgMbImKLpGuBnojYSFKqmQd8SRLAwxFxQUTskfQnJB8WANdGxJ5J+UvMJkmjVFNOL2oGuE5vhdI26AEiYhOwaUTbNZnpc0dZdwOwYbwdNOu07Ii+UhKA6/RWKP5mrFkbQ6dXpl+YguSOU2ZF4aA3a6NRpmmcXgn4BuFWKA56szYa59FXMkHvGr0ViYPerI2qa/RWcA56szbqw86jT0s3rtFbgTjozdoYfnqlSzdWPA56szZq9cOD3l+YsiJx0Ju1UcucXllJvzBVd43eCsRBb9ZG0xG9a/RWIA56szbqrtFbwTnozdpojN59eqUVlYPerI1GqGdPr/QNwq1Icl3UzOy5rJHp5ZKIcI3eisdBb9ZG47o25ZKoh2v0VjwOerM2sgdjGyN61+itSBz0Zm3UGqUbifAXpqyAHPRmbWTPo2+cv+Dr0VuROOjN2sgGfZBMe0RvReKgN2vj0EXNoDGi9yUQrEgc9GZt1IbOuimBR/RWQA56szYa58xXSkJpm78wZUXioDdrozF6r5R1qM0HY61Acl0CQdI6SfdJ2irpqibzXyHpdklVSReNmFeTdEf6s3GiOm42Varp+ZWVUmnoEgiu0VuRtB3RSyoD1wGvBnYCmyVtjIh7Mos9DLwZeHeTTRyIiLOOvKtmnTE0os9eAsE1eiuQPKWbtcDWiNgGIOlGYD0wFPQRsT2d58KlzThDNfqyL4FgxZSndLMM2JF5vDNty2u2pB5Jt0q6sNkCki5Pl+np7e0dw6bNJl9j9D6rXPKNR6yQpuIyxadExBrgN4C/lnTayAUi4vqIWBMRa5YuXToFXTLLr1Gjz954xDV6K5I8Qb8LWJ55fFLalktE7Ep/bwNuAc4eQ//MOm4wU6Nv3DPWNXorkjxBvxlYJWmlpC7gYiDX2TOSFkvqTqeXAL9AprZvVgS1ej05h16+laAVU9ugj4gqcAVwE3Av8MWI2CLpWkkXAEh6saSdwBuAT0jakq5+JtAj6cfAt4APjjhbx2zaq9ZiKOBdo7ciyvWFqYjYBGwa0XZNZnozSUln5HrfA154hH0066jBWjArudANac77evRWKL5nrFkbtXp9aCQvJTcI9yUQrEgc9GZtDNaDWZnLH5RK8sFYKxQHvVkbtVoMnW0Dydk3vvGIFYmD3qyNwUzpBpIDsq7RW5E46M3aqNaGl26SGr2D3orDQW/WRq0eh43oXaO3InHQm7UxWKsPnV4JaenGNXorEAe9WRu1egy76UilVHKN3grFQW/WxmA90vvFJsqu0VvBOOjN2qjW6sxyjd4KzEFv1kZ1ROmm7G/GWsE46M3aqNbqh39hyiN6KxAHvVkbzUf0DnorDge9WRvVWlBxjd4KzEFv1sZArU5XxWfdWHE56M3aGKzV6Sq7Rm/F5aA3a2Ogevg3Y126sSJx0Ju1MejSjRWcg96sjcNH9CWP6K1QHPRmbQzU6nRXhtfo6w56KxAHvVkb2ZuDQxL0gzV/M9aKw0FvNopaPajVY1iNvntWmYGqg96KI1fQS1on6T5JWyVd1WT+KyTdLqkq6aIR8y6TdH/6c9lEddxsKjRG7tkRfVe5RL+D3gqkbdBLKgPXAecBq4FLJK0esdjDwJuBG0asezTwfuAlwFrg/ZIWH3m3zaZGI9CHj+hL9FdrneqS2ZjlGdGvBbZGxLaIGABuBNZnF4iI7RFxJzBymPNa4OaI2BMRe4GbgXUT0G+zKdEY0XdlrnXTXfGI3oolT9AvA3ZkHu9M2/LIta6kyyX1SOrp7e3NuWmzyTfQbERfKTvorVCmxcHYiLg+ItZExJqlS5d2ujtmQ5rW6CslBqp1wrcTtILIE/S7gOWZxyelbXkcybpmHdd8RJ9Me1RvRZEn6DcDqyStlNQFXAxszLn9m4DXSFqcHoR9TdpmVggDTUb0jaAf8Ln0VhBtgz4iqsAVJAF9L/DFiNgi6VpJFwBIerGkncAbgE9I2pKuuwf4E5IPi83AtWmbWSE0HdHPKgPQP+igt2Ko5FkoIjYBm0a0XZOZ3kxSlmm27gZgwxH00axjBmtJHT57meLucqN041MsrRimxcFYs+mq+YjeNXorFge92SgODiaj9tmV8lDbUI3eQW8F4aA3G8XBtDwze9bw0yvBI3orDge92SgODDSCPjuibxyMdY3eisFBbzaKg+mofXjQe0RvxeKgNxtFY9SeLd00RvSu0VtROOjNRtGsdOMavRWNg95sFAerNcolNf1mrM+jt6Jw0JuN4uBgnaMyo3nwefRWPA56s1EcHKwNq8/DoXPqD/qsGysIB73ZKA4M1oYOvjbM6U4e9w046K0YHPRmo+gfrB82ou8ql6iUxP7+aod6ZTY2DnqzURwcrHFU1/ARvSTmdJU9orfCcNCbjaJvoHbYwViAud0Vj+itMBz0ZqPY119l/uxZh7V7RG9F4qA3G8W+/ipzuw+/bcO87gr7Bzyit2Jw0JuNYl9/lXlNgn5OV4W+fo/orRgc9Gaj2HewyrzuZjX6skf0VhgOerMWavXgwGCNed3NavQ+GGvF4aA3a2FfGuTzZh9euklG9C7dWDE46M1aGAr6ZqWbrgp9HtFbQTjozVrYPxT0TUo33RX6BmvU6zHV3TIbs1xBL2mdpPskbZV0VZP53ZK+kM7/gaQVafsKSQck3ZH+fHyC+282aZ49mAT93CYj+vndFSJgnw/IWgEcXnwcQVIZuA54NbAT2CxpY0Tck1nsrcDeiDhd0sXAXwBvTOc9EBFnTWy3zSZfo3Qzv0mNftGcZJT/dN8gC5p8ocpsOskzol8LbI2IbRExANwIrB+xzHrgM+n0l4FXSdLEddNs6u072Lp0s2hOFwB7+wamtE9m45En6JcBOzKPd6ZtTZeJiCrwNHBMOm+lpB9J+raklzd7AkmXS+qR1NPb2zumP8BssuxJQ3zx3MODfnE6on+qb3BK+2Q2HpN9MPZR4OSIOBu4ErhB0oKRC0XE9RGxJiLWLF26dJK7ZJbPnn1p0Kej96xG6cYjeiuCPEG/C1ieeXxS2tZ0GUkVYCGwOyL6I2I3QETcBjwAnHGknTabCnv7BlgwuzLsfrENjdKNR/RWBHmCfjOwStJKSV3AxcDGEctsBC5Lpy8CvhkRIWlpejAXSacCq4BtE9N1s8m1e/8AR889fDQPsOgol26sONqedRMRVUlXADcBZWBDRGyRdC3QExEbgU8Bn5O0FdhD8mEA8ArgWkmDQB14R0TsmYw/xGyi7R0l6CvlEvO7Ky7dWCG0DXqAiNgEbBrRdk1m+iDwhibrfQX4yhH20awjdu8fYNmi2S3nL57bxZ79Dnqb/vzNWLMWRhvRAxy3oJvHnzk4hT0yGx8HvVkTtXrw5L5+lszrbrnMcQtmO+itEBz0Zk30PttPtR6cuOiolsucsHA2jz59kAhf78amNwe9WRO7njoAwLJRgv64BbPpr9Z5+oDPvLHpzUFv1sQjadCPNqI/fmFyoPYxl29smnPQmzVxKOhbn3XTGO0/vLtvSvpkNl4OerMmHt7Tx4LZFeaPcmXKU5fOA2Dbk/unqltm4+KgN2ti6xP7WHXc/FGXWXjULJbO7+aBJ/ZNUa/MxsdBb9bE/U/s44zj5rVd7vSl83ig10Fv05uD3myE3fv62bN/gNOPHX1ED3DasXPZ+sQ+n2Jp05qD3myEO3c9DcCZx7cP+jNPWMAzB6s8vMcHZG36ctCbjfCjh/ZSEvzc8kVtl33RKYsB6Nm+d5J7ZTZ+DnqzEXoe2svzjl/A3O721/w749j5zJ9doechB71NXw56s4xnDg6yefseXrFqSa7lSyWxdsXRfHdrr+v0Nm056M0yvn1fL4O14NzVx+Ve57XPP54dew6w5ZFnJrFnZuPnoDfL+OrtO1k6v5tzTl6ce51Xrz6Ockn86x0j77BpNj046M1S25/czy0/7eWStSdTLin3eovndrHuBcdz4+Yd7OuvTmIPzcbHQW+W+vBN9zG7UubSl5w85nUvf/mpPHuwysdu2ToJPTM7Mg56M+Brdz7K1+56lP/2i6dy7ILWFzJr5eeWL+LXzlnG9d/Zxl07n56EHpqNn4PenvNu3bab93z5x5y1fBG/+0unj3s7V59/JsfOn83bPrvZl0WwacVBb89ZtXqw4bsP8lsbfsiJi47i+t98EV2V8f+TOGZeN5968xqqteD1H/se/37nIz7l0qYFTbc34po1a6Knp6fT3bAZbO/+ATbd/Sif+q8H2fbkfl75M8fy4Yt+lmNGuT/sWDy0ez9X3PAj7tr1NGctX8RbXraSV/3Msbm+gGU2XpJui4g1TeflCXpJ64C/AcrAP0TEB0fM7wY+C7wI2A28MSK2p/PeC7wVqAHvioibRnsuB71NpKcPDLJjTx8PPrmfO3c+xe0PP8UdO56iVg/OPGEB73rl6ax7wfFI+c+yyaNaq/P5zTv45He28fCeProrJdauPJqzly/iBcsWsmLJXE4+eg6zZ5Un9HntueuIgl5SGfgp8GpgJ7AZuCQi7sks87vAz0bEOyRdDLwuIt4oaTXweWAtcCLwDeCMiKi1ej4H/cwSEdTqQS2Ceh1q6eP6UFtk2g7NjzjUPlgLDg7W6K/WOThYG5ruH6xxcLBOf7XGs/1Vnto/yN6+AZ7qS3737uvnqb5D93PtKpd4wbIFvOTUY/jlF57A809cMOEBP1KtHvzwwT3ctOUxfvDgHn76+LPU6of+zS2Z182SeV0cPbeLxXO7OGZuF3O6KszpKjOnq8xRjd+zKsyeVWJWuUSlJCrlErPKolJKf6fts8olKmVRlihJqAQCSo3HakyDMr+t+EYL+jz/l1wLbI2IbenGbgTWA/dkllkPfCCd/jLwUSXvnvXAjRHRDzwoaWu6ve+P5w8ZzVN9A1z08WSz2Q+vYR9j0XRyaPnhbdllo3l7k8/IVs/dbnutlqXVczd5zpbP12Z/NNvWkfSfgHoa1FNVGewql1g0ZxaL53SxcM4sTl06l7Urj+bko+ckP8fM4fRj59FdmdoRdLkkXnraMbz0tGMA6Buo8tPH9/HQ7v08vLuPXU8dYPf+AfbuH+DeR55hT98A+/urDNamtqRa0qEPAzT8cdo0pPHBkP180Mh5w5YfudShNh22DIj2288a2tYRbGPY9kad236BduuP9vxnnrCAv7vk7HY9GLM8Qb8M2JF5vBN4SatlIqIq6WngmLT91hHrLhv5BJIuBy4HOPnksZ/DDMk/qOdl7wjU5AVOn6tFe/5lyfmGOrz98Df68OkW88ewvVb9H9b98fa/zf7I9qek5DUpSZRLykzTpC0dgZYOzW+0lyS6KmJ2pUz3rBLdlTKzh34nbbMrZWaVVYiR6ZyuCmctX8RZba6MOVir0zdQ48BAjb6BKn0Dyf9iqrU61XowWKtTrQXVep3B7O90ul4P6ukHbkTy4Zx93JjfaI8I6nFomcaHdT0Y9j+QhmYDgkMf+nkHMjHi8di3EU1GK4f607qPrbT7eG1XAWn78dxmgeWLW9+M/khMi6NDEXE9cD0kpZvxbGP+7Flc96ZzJrRfZp0yq1xi4VElFh7V+p61ZnnlOZdsF7A88/iktK3pMpIqwEKSg7J51jUzs0mUJ+g3A6skrZTUBVwMbByxzEbgsnT6IuCbkfwfZyNwsaRuSSuBVcAPJ6brZmaWR9vSTVpzvwK4ieT0yg0RsUXStUBPRGwEPgV8Lj3Yuofkw4B0uS+SHLitAu8c7YwbMzObeP7ClJnZDDDa6ZW+BIKZ2QznoDczm+Ec9GZmM5yD3sxshpt2B2Ml9QIPHcEmlgBPTlB3JpL7NTbu19i4X2MzE/t1SkQsbTZj2gX9kZLU0+rIcye5X2Pjfo2N+zU2z7V+uXRjZjbDOejNzGa4mRj013e6Ay24X2Pjfo2N+zU2z6l+zbgavZmZDTcTR/RmZpbhoDczm+EKG/SS3iBpi6S6pDUj5r1X0lZJ90l6baZ9Xdq2VdJVU9DHL0i6I/3ZLumOtH2FpAOZeR+f7L6M6NcHJO3KPP/5mXlN990U9evDkn4i6U5J/yxpUdre0f2V9mFK3zuj9GO5pG9Juid9//9e2t7yNZ3i/m2XdFfah5607WhJN0u6P/29eIr79LzMfrlD0jOSfr8T+0zSBklPSLo709Z0/yjxt+l77k5J47+zUkQU8gc4E3gecAuwJtO+Gvgx0A2sBB4gubxyOZ0+FehKl1k9hf39CHBNOr0CuLuD++4DwLubtDfdd1PYr9cAlXT6L4C/mCb7q6PvnRF9OQE4J52eD/w0fd2avqYd6N92YMmItg8BV6XTVzVe1w6+lo8Bp3RinwGvAM7Jvp9b7R/gfODrJDfr/HngB+N93sKO6CPi3oi4r8msoRuSR8SDQOOG5EM3OY+IAaBxk/NJl94o/deBz0/F8x2BVvtuSkTEf0RENX14K8kdyaaDjr13RoqIRyPi9nT6WeBemtyHeZpZD3wmnf4McGHnusKrgAci4ki+fT9uEfEdknt2ZLXaP+uBz0biVmCRpBPG87yFDfpRNLuZ+bJR2qfCy4HHI+L+TNtKST+S9G1JL5+ifmRdkf53cEPmv9Kd3EcjvYVkNNPQyf01nfbLEEkrgLOBH6RNzV7TqRbAf0i6TdLladtxEfFoOv0YcFxnugYkN0XKDrimwz5rtX8m7H03rYNe0jck3d3kpyOjqWZy9vEShr+5HgVOjoizgSuBGyQtmMJ+fQw4DTgr7ctHJvK5j6BfjWWuJrkj2T+lTZO+v4pG0jzgK8DvR8QzdPA1HeFlEXEOcB7wTkmvyM6MpCbRkXO6ldwK9QLgS2nTdNlnQyZr/7S9lWAnRcS541httBuST/iNytv1UcnN0n8NeFFmnX6gP52+TdIDwBnAhN1aK+++k/RJ4N/Th5N+M/cc++vNwK8Ar0rf9FOyv9qYVje5lzSLJOT/KSK+ChARj2fmZ1/TKRURu9LfT0j6Z5Ky1+OSToiIR9PSwxOd6BvJh8/tjX01XfYZrffPhL3vpvWIfpxa3ZA8z03OJ8O5wE8iYmejQdJSSeV0+tS0j9umoC+N58/W+V4HNM4A6OjN3CWtA/4IuCAi+jLtHd1fdO69c5j0eM+ngHsj4i8z7a1e06ns21xJ8xvTJAfX7ybZV5eli10G/OtU9y017H/W02GfpVrtn43Ab6Vn3/w88HSmxDM2U33UewKPXr+OpGbVDzwO3JSZdzXJWRL3Aedl2s8nOUvhAeDqKernp4F3jGh7PbAFuAO4HfjVKd53nwPuAu5M30wntNt3U9SvrSQ1yTvSn49Ph/3VqfdOi368jOS/9ndm9tP5o72mU9i3U0nOSPpx+npdnbYfA/wncD/wDeDoDvRtLrAbWJhpm/J9RvJB8ygwmObXW1vtH5Kzba5L33N3kTm7cKw/vgSCmdkMNxNLN2ZmluGgNzOb4Rz0ZmYznIPezGyGc9Cbmc1wDnozsxnOQW9mNsP9f0TPp8jChvvRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgsUlEQVR4nO3deZgc9X3n8fene07d6LAEkkDCFmCwzWEF7PVFHmMbcIziOLZh4yROWBNvQkzWORYHL8vifXaDrxxPSBxs49jENsFHHCWB4NjG69gxRgJzSICMkASSkIQOdM7R093f/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+ld7a0uxeyYHOhmx1CpBmu27uPcRTNbXYrZkBzoZsewYedBekoVXrVoVqtLMRuSA93sGB7bug+AV/oI3XLAgW52DE/vPEixIJbOndrqUsyG5EA3O4ZNu3pYfFK3/3eL5YJ7qdkxbNx1iCU+OreccKCbDSIi2LT7EEvmONAtHxzoZoPYeaCfnlLF4+eWGw50s0Fs3dsLwKKTultciVlzHOhmg9ixvx+A+TO6WlyJWXMc6GaDeP5AH+BAt/xwoJsNYsf+PooFMWdqR6tLMWuKA91sEDv29zNvWieFglpdillTHOhmg9ixv4/5MzpbXYZZ0xzoZoPYc6jEnGkOdMuPpgJd0qWS1klaL+n6BstPlXSfpJ9KelTS5SNfqtnY2tszwKwp/h/olh9DBrqkInArcBlwNnCVpLPrVvsocFdEnA9cCfzVSBdqNtb29Q4wq9snRC0/mjlCvxBYHxEbIqIE3AmsqFsngBnp9EzguZEr0WzsDVSqHOwv+wjdcqWZQF8IbM7Mb0nbsm4C3idpC3A38LuNnkjSNZJWS1q9c+fO4yjXbGzs6x0AcKBbrozUSdGrgL+NiEXA5cAdkl703BFxW0Qsj4jl8+b5Duo2fu3tSQJ9ZrcD3fKjmUDfCizOzC9K27KuBu4CiIgfA13A3JEo0KwV9vWWAJg1xWPolh/NBPoqYJmkpZI6SE56rqxb51ngzQCSXk4S6B5TsdyqHaHP8hG65ciQgR4RZeBa4F7gCZKrWdZKulnSFelqvw98QNIjwFeB90dEjFbRZqPNQy6WR23NrBQRd5Oc7My23ZiZfhx43ciWZtY6h0plAKZ1NfUjYjYu+JOiZg0c6q8AMLXDgW754UA3a+BQf5mCoKvdPyKWH+6tZg0cKpWZ2tGG5P+0aPnhQDdroKe/wpTOYqvLMBsWB7pZA7UjdLM8caCbNdBT8hG65Y8D3ayBg/0+Qrf8caCbNdBTKjO104Fu+eJAN2ugp7/ClA4PuVi+ONDNGvBJUcsjB7pZA4f6Kx5ysdxxoJs10FMq093hHw/LF/dYszrlSpVqQGebx9AtXxzoZnVKlSoAHW3+8bB8cY81q1Mqp4Fe9I+H5Yt7rFmdw4HuI3TLGfdYszr9DnTLKfdYszq1MfROB7rljHusWR2PoVteucea1fEYuuWVe6xZHV+2aHnlHmtWx0MullfusWZ1PORieeUea1bHly1aXrnHmtXxZYuWV+6xZnWOjKH7n3NZvjjQzep4DN3yyj3WrE6pXAEc6JY/7rFmdXwduuWVe6xZHV+HbnnlHmtWpxbo7UW1uBKz4Wkq0CVdKmmdpPWSrh9knfdIelzSWklfGdkyzcZOf6VKR1sByYFu+TLkbc0lFYFbgbcAW4BVklZGxOOZdZYBHwFeFxEvSHrJaBVsNtpK5SqdHm6xHGqm114IrI+IDRFRAu4EVtSt8wHg1oh4ASAinh/ZMs3GTqlc9QlRy6Vmeu1CYHNmfkvalnUGcIakH0m6X9KljZ5I0jWSVktavXPnzuOr2GyUOdAtr0aq17YBy4CLgauAz0qaVb9SRNwWEcsjYvm8efNG6KXNRlap4kC3fGqm124FFmfmF6VtWVuAlRExEBEbgZ+RBLxZ7pTKVV+yaLnUTK9dBSyTtFRSB3AlsLJunW+RHJ0jaS7JEMyGkSvTbOz0e8jFcmrIXhsRZeBa4F7gCeCuiFgr6WZJV6Sr3QvslvQ4cB/whxGxe7SKNhtNHkO3vBryskWAiLgbuLuu7cbMdAAfTh9mueYhF8sr91qzOv0+KWo55V5rVqdUrvrmFpZL7rVmdUrlio/QLZfca83qlCoeQ7d8cq81q+OrXCyv3GvN6jjQLa/ca83qJJct+gbRlj8OdLM6/l8ullfutWYZ1WowUAkHuuWSe61ZRu0G0b4O3fLIvdYsoxbovmzR8si91iyjdoNoD7lYHrnXmmU40C3P3GvNMg4HuodcLIfca80yDo+h+wjdcsi91izDQy6WZ+61Zhn9DnTLMfdas4zaEXqnx9Ath9xrzTI8hm555l5rluExdMsz91qzDAe65Zl7rVlGqVIBfB265ZN7rVmGj9Atz9xrzTIc6JZn7rVmGf2HL1v0HYssfxzoZhmH/x96u380LH/ca80y/M+5LM/ca80ySuUqbQVRKKjVpZgNmwPdLKNU9g2iLb/cc80yShUHuuWXe65ZRqlc9fi55VZTPVfSpZLWSVov6fpjrPcuSSFp+ciVaDZ2PORieTZkz5VUBG4FLgPOBq6SdHaD9aYD1wE/GekizcZKv4dcLMea6bkXAusjYkNElIA7gRUN1vsYcAvQN4L1mY2pUrlKZ5s/VGT51EygLwQ2Z+a3pG2HSboAWBwR/3KsJ5J0jaTVklbv3Llz2MWajbZ+D7lYjp1wz5VUAD4N/P5Q60bEbRGxPCKWz5s370Rf2mzE9Q9U6HSgW04103O3Aosz84vStprpwCuA70vaBLwGWOkTo5ZHpUrVgW651UzPXQUsk7RUUgdwJbCytjAi9kXE3IhYEhFLgPuBKyJi9ahUbDaK+gcc6JZfQ/bciCgD1wL3Ak8Ad0XEWkk3S7pitAs0G0v95YpPilputTWzUkTcDdxd13bjIOtefOJlmbVGf9lH6JZf7rlmGaVy1f8613LLPdcso98f/bccc881y+gvV+hs9xi65ZMD3SwVEeknRf1jYfnknmuWKleDavhuRZZf7rlmqcM3iPZJUcsp91yzVP9ABcDXoVtuOdDNUqVKeoTuMXTLKfdcs1T/QBLo/m+LllfuuWapw2PoHnKxnHKgm6VKZQ+5WL6555ql+svJSVEPuVheueeapWpDLg50yyv3XLNUTyk5Qp/S4TF0yycHulmqp1QGHOiWXw50s1RveoTe3dHUbQLMxh0Hulnq8JCL/9ui5ZQD3SzVO1A7QnegWz450M1SPaUyBfk6dMsv91yzVE+pwpSONiS1uhSz4+JAN0v1liq+wsVyzYFulupxoFvOOdDNUj2lii9ZtFxzoJulegfKPkK3XHOgm6U85GJ550A3S/WWKnT7Q0WWYw50s5SP0C3vHOhmqYP9ZaZ3tbe6DLPj5kA3AyKCfb0DzOj2VS6WXw50M5Lhlko1mOEjdMsxB7oZsK93AICZ3Q50y6+mAl3SpZLWSVov6foGyz8s6XFJj0r6rqTTRr5Us9Gzvy8J9BkOdMuxIQNdUhG4FbgMOBu4StLZdav9FFgeEa8Cvg58fKQLNRtN+3p8hG7518wR+oXA+ojYEBEl4E5gRXaFiLgvInrS2fuBRSNbptno2t+X3H7OY+iWZ80E+kJgc2Z+S9o2mKuBexotkHSNpNWSVu/cubP5Ks1GmcfQbSIY0ZOikt4HLAc+0Wh5RNwWEcsjYvm8efNG8qXNTsj+3toYui9btPxqpvduBRZn5helbUeRdAlwA/CmiOgfmfLMxsbenhIS/mCR5VozR+irgGWSlkrqAK4EVmZXkHQ+8DfAFRHx/MiXaTa6dh4sMWdqB8WC71Zk+TVkoEdEGbgWuBd4ArgrItZKulnSFelqnwCmAV+T9LCklYM8ndm4tPNAP3Ondba6DLMT0tSAYUTcDdxd13ZjZvqSEa7LbEztPNjPvOkOdMs3f1LUDNh1wIFu+edAt0kvInyEbhOCA90mvf19ZUrlKvM8hm4550C3Se+5vb0AnDyzu8WVmJ0YB7pNes/uSf5rxeLZDnTLNwe6TXqba4F+0pQWV2J2YhzoNulteaGX6Z1tzJriT4lavjnQbdJ7ZvchFs2eguRPiVq+OdBt0vvZjoOcMX9aq8swO2EOdJvU9vUOsHVvL2ctmNHqUsxOmAPdJrV12w8AcNbJ01tcidmJc6DbpPbEtv0AnLXAgW7550C3SW31My8wf0YnC2Z0tboUsxPmQLdJKyL4yYbdXLR0jq9wsQnBgW6T1sZdh3j+QD8XnT671aWYjQgHuk1a//b4DgDeuMz3t7WJwYFuk9Y9a7bzioUzWDzbH/m3icGBbpPSUzsO8PDmvbz9lae0uhSzEeNAt0npiz/eREdbgfcsX9TqUsxGjAPdJp3t+/r4+oNbuOLcU5jjm1rYBOJAt0nn4/c+SbUK1715WatLMRtRDnSbVL69djvffGgr/+UNS30y1CYcB7pNGuu2H+APvvYI55wyg+su8dG5TTwOdJsUnti2n1/9/E/oai/ymfe9ms62YqtLMhtxDnSb0CKClY88xy//9X8gwR1XX+ShFpuw2lpdgNloeWrHAW7513V854kdnLd4Fp9536tZMNP/hMsmLge6TSiVavCj9bu4c9Wz3LNmO1Pai3zksrO4+vVLaSv6D1Kb2BzolnsH+8v8+OndfH/d83zvyefZtq+Pk6a088E3vZRr3nA6J03taHWJZmPCgW65ERHsPlRi465D/GzHAR7dvI9HtuzlZzsOUA2Y2lHkdS+byw1vfzlvOXu+T3zapONAt3EhIugbqLLzQD/b9/exI/PYvr+fZ/f0sGHnQQ70lQ9vM2tKO+cumsXbzlnAhUtn83NLZtPR5mEVm7yaCnRJlwJ/DhSBz0XEn9Qt7wS+BLwa2A28NyI2jWypNt5EBKVKld5ShZ70kUyX6RlIpntLlXS6zKH+Cvv7BtjXM8C+3uSxt/fIdKlcfdFrdLYVmD+ji1NnT+EXz1vI0rlTOX3eVF46bxqLTur2jSnMMoYMdElF4FbgLcAWYJWklRHxeGa1q4EXIuJlkq4EbgHeOxoFj6WIIAKqEQTp1+DFbVUIgmok29S+Zrep37a2fW39ciVpK1eDSjWdzrRV0/Zy9cVtlWpQGaStUq17pG0DlSoDlSqlcpWBSlAqVykd1VZN26JB25Hpagxvn07vbGPmlHZmdiePM+ZPY2Z3OzPS+XnTOlkws4v5M7qYP72LGd1tDm2zJjVzhH4hsD4iNgBIuhNYAWQDfQVwUzr9deAvJSkihvnjPrS7Vm3mtn/fQDUCDodjGqjVRkF6dFu1WgvSI2GaDdfIPOdEUxC0FQoUCtBeKNDRVqC9WPsq2osFOjNtUzraMm06vH52ve6OIt3tRaZ0FOnuKDKloy0znSyrtXe3FykWHM5mo6WZQF8IbM7MbwEuGmydiChL2gfMAXaNRJFZs6a0c+b86SAoSBQEIpmWhETaJgoFkjYy66briGS+UEiWZ7ct6EhbIfuch6eVec3aemkdBR31mkdtW6ujkLx+dtu2QvJabcX0axq8bYUCxQJNtRUKUBykrZjWZWYT15ieFJV0DXANwKmnnnpcz/HWcxbw1nMWjGRZZmYTQjOXBGwFFmfmF6VtDdeR1AbMJDk5epSIuC0ilkfE8nnzfB9HM7OR1EygrwKWSVoqqQO4ElhZt85K4NfT6V8Gvjca4+dmZja4IYdc0jHxa4F7SS5bvD0i1kq6GVgdESuBzwN3SFoP7CEJfTMzG0NNjaFHxN3A3XVtN2am+4B3j2xpZmY2HP5YnZnZBOFANzObIBzoZmYThAPdzGyCUKuuLpS0E3jmODefyyh8CnUEuK7hGa91wfitzXUNz0Ss67SIaPhBnpYF+omQtDoilre6jnqua3jGa10wfmtzXcMz2erykIuZ2QThQDczmyDyGui3tbqAQbiu4RmvdcH4rc11Dc+kqiuXY+hmZvZieT1CNzOzOg50M7MJYtwHuqR3S1orqSpped2yj0haL2mdpLdl2i9N29ZLun4Mavx7SQ+nj02SHk7bl0jqzSz7zGjXUlfXTZK2Zl7/8syyhvtujOr6hKQnJT0q6R8kzUrbW7q/0hrGtO8co47Fku6T9Hja/69L2wd9T8ewtk2SHktff3XaNlvSv0l6Kv160hjXdGZmnzwsab+k32vF/pJ0u6TnJa3JtDXcP0r8RdrfHpV0wQm9eHIfzfH7AF4OnAl8H1ieaT8beAToBJYCT5P8e99iOn060JGuc/YY1vsp4MZ0egmwpoX77ibgDxq0N9x3Y1jXW4G2dPoW4JZxsr9a2nfqajkZuCCdng78LH3fGr6nY1zbJmBuXdvHgevT6etr72kL38ftwGmt2F/AG4ELsn15sP0DXA7cQ3IHy9cAPzmR1x73R+gR8URErGuwaAVwZ0T0R8RGYD3JDa0P39Q6IkpA7abWo07JTTvfA3x1LF7vBAy278ZERHw7Isrp7P0kd8EaD1rWd+pFxLaIeCidPgA8QXLv3vFqBfDFdPqLwC+2rhTeDDwdEcf7SfQTEhE/ILkvRNZg+2cF8KVI3A/MknTy8b72uA/0Y2h08+qFx2gfC28AdkTEU5m2pZJ+Kun/SXrDGNWRdW36p9ztmT+DW7mP6v0myRFKTSv313jaL4dJWgKcD/wkbWr0no6lAL4t6UEl9wkGmB8R29Lp7cD8FtRVcyVHH1S1en/B4PtnRPvcuAh0Sd+RtKbBoyVHR400WeNVHN2RtgGnRsT5wIeBr0iaMYZ1/TXwUuC8tJZPjeRrn0BdtXVuAMrAl9OmUd9feSNpGvAN4PciYj8tfE8zXh8RFwCXAb8j6Y3ZhZGMJbTkemglt8m8Avha2jQe9tdRRnP/NHXHotEWEZccx2bHunn1UDe1HrahalRyc+xfAl6d2aYf6E+nH5T0NHAGsPpE62m2rkx9nwX+OZ1t5sbfo1qXpPcDvwC8Oe3gY7K/hjDq+2U4JLWThPmXI+KbABGxI7M8+56OmYjYmn59XtI/kAxV7ZB0ckRsS4cMnh/rulKXAQ/V9tN42F+pwfbPiPa5cXGEfpxWAldK6pS0FFgGPEBzN7UeDZcAT0bEllqDpHmSiun06WmNG8agltrrZ8fi3gnUzroPtu/Gqq5LgT8CroiInkx7S/cXres7L5Kej/k88EREfDrTPth7OlZ1TZU0vTZNcoJ7DUffKP7XgX8cy7oyjvorudX7K2Ow/bMS+LX0apfXAPsyQzPDN5Znf4/zjPE7ScaV+oEdwL2ZZTeQXJWwDrgs0345yVUBTwM3jFGdfwt8sK7tXcBa4GHgIeAdY7zv7gAeAx5NO87JQ+27MaprPcm44cPp4zPjYX+1qu8MUsfrSf4sfzSzny4/1ns6RnWdTnL1zyPpe3VD2j4H+C7wFPAdYHYL9tlUYDcwM9M25vuL5BfKNmAgza6rB9s/JFe33Jr2t8fIXMl3PA9/9N/MbILI85CLmZllONDNzCYIB7qZ2QThQDczmyAc6GZmE4QD3cxsgnCgm5lNEP8f/3t2Mk/uyNkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEICAYAAAC6fYRZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiBUlEQVR4nO3deXxcdb3/8dcn+942Tehe2kpZCsgWWdyvgAK/K+WK1wcKiorww4v7jt6fC8pVr9frT396VVSuiBtcUKkIVxb1B/qjQFlEGpaWUmia0KYpyWSZLJP5/P44J+k0TdK0kzOTnHk/H495zJnzPcvnJO17vvmeM2fM3RERkXgqyncBIiISHYW8iEiMKeRFRGJMIS8iEmMKeRGRGFPIi4jEmEJe9mJmnzezn+a7DgAz+5KZ7TKzFw5i3VeY2SYz6zGz8yIo76CY2VYzOyOibbuZHRbFtsPtv8rMnopq+xINhXyBCUNv5JE2s2TG6wvzXd8IM1sOfBRY4+4LD2ITVwHfdvcad//NBPt4m5ltCI+9zcxuN7NXZlF2rLn7ve5+xMjrKN+wZPoo5AtMGHo17l4DPA+8MWPez/JdX4blQIe77zzI9Q8FNk7UaGYfAf438C/AgnB//wGsPcj9icxICnkZT5mZ/cTMus1so5k1jTSY2WIzu9nM2s3sWTP7wEQbMbMfm9n3zOzOcFv/18wOzWh3M/uAmW0Jh2W+ZmZFYe/wTmBx2Mv+8QTbv9TMNpvZbjNbZ2aLw/nPAKuA34brl49Zbw5BT/8Kd/+Vu/e6+5C7/9bdPx4uc7KZ3WdmnWEv/9tmVha2rQhrL8nY5p/M7D1jansiPO5mMzsxo4TjzewxM+sysxvMrCJc53Eze2PGNkrDn8sJExz/x8PaWs3s3WPays3s38zseTPbEf4eKsO215pZi5l91Mx2htt4V8a654Q1d5vZdjP7WOZ64fT1BG+MIz/jT5jZ78zs/WPqeMzM/mG8+iVH3F2PAn0AW4Ezxsz7PNAPnAMUA18G1odtRcBDwGeBMoIg3QK8YYLt/xjoBl4NlAPfBP6c0e7AH4F6gsB4GnhP2PZaoGWS2l8H7AJODLf9f4B7Jju2jLazgBRQMsn2TwJOBUqAFcATwIfCthVh7SUZy/8po/Z/BLYDLwMMOAw4NKOuB4DF4XE/AVwetn0CuCFjm2uBv01yDDuAY4Bq4OdhTYeF7d8A1oX7qAV+C3w542ebInijKw1/133AvLC9DXhVOD0POHG838nYnzHwFuD+jNfHAR1AWb7/rRfyQz15Gc+f3f02dx8Grif4zwpBaDW6+1XuPujuW4AfABdMsq3fufs97j4AfAY4zcyWZbR/1d13u/vzBMMnb51ijRcC17r7w+G2rwy3vWIK684Hdrl7aqIF3P0hd1/v7il33wp8H3jNFGt7D/Cv7v6gBza7+3MZ7d9y91Z3300QvseH838KnGNmdeHrtxP8/MfzFuA/3f1xd+8leHMGwMwMuAz4cPiz7SYYlsr8PQ0BV3nwF8xtQA9wREbbGjOrc/cX3f3hKR73OuBwM1udUf8N7j44xfUlAgp5GU/m1Sx9QEU4NHEowRBK58gD+DTBmPZEto1MuHsPsJugF7tPO/DcmLbJLA6Xz9x2B7BkCut2AA2Zwy1jmdnhZnarmb1gZgmCkGyYYm3LgGcmaR/7860BcPdW4C/A+WY2FzgbmOg8yWL2/dmNaASqgIcyfk//Hc4f0THmTW60DuB8gt79c+EQ22mTHMsod+8HbgAuMrMigjfsid6kJEcU8nIgtgHPuvvcjEetu58zyTqjvXYzqyEYPmgdr51gyCazbTKtBG86I9uuJuihb5/CuvcBA8B5kyzzXeBJYLW71xG8mVnY1hs+V2Usn3kF0DbgJVOoYzzXARcRDPnc5+4THU8b+/7sRuwCksDRGb+nOR6cbN+v8C+QtcAhwG+AGydadIL6LwROB/rc/b6p7FOio5CXA/EA0G1mnzSzSjMrNrNjzOxlk6xzjpm9Mjxp+UWC8f3MHujHzWxeOITzQYKe4FT8AniXmR0fnlj9F4Lx4K37W9HduwjOK3zHzM4zs6rwJOfZZvav4WK1QALoMbMjgfdmrN9O8GZyUfgzeDd7h/oPgY+Z2UkWOCzzhPN+/IbgPMMHgZ9MstyNwDvNbI2ZVQGfy6gvTTCM9g0zOwTAzJaY2Rv2t3MzKzOzC81sjrsPEfwM0hMsvoPgvMyoMNTTwNdRL35GUMjLlIVj9H9PMIb8LEGP8YfAnElW+zlBAO0mOJl50Zj2WwhO5j4K/A740RRruQv4X8DNBL3alzD5uYGx638d+Ajwz0A7Qe/7fQQhC/Ax4G0EJ45/wL5vPpcCHycY+jka+H8Z2/4v4GqCY+8Ot1k/xbqS4TGtBH41yXK3E5zD+AOwOXzO9Mlw/vpwuOku9oy578/bga3hepcT9MzH82Xgn8MhoY9lzP8JcCzBOQbJM3PXl4ZINCy49LHF3f95gnYnGA7ZnNPCZjgz+yxwuLuPfUOcFczsHcBl7q4Pls0AE554EpHcM7N64BKC3vSsEw4d/RPBB8tkBtBwjcgMYWaXEgwb3e7u9+S7ngMVjvm3E4zV/zzP5UhIwzUiIjGmnryISIzNqDH5hoYGX7FiRb7LEBGZVR566KFd7t44XtuMCvkVK1awYcOGfJchIjKrmNlzE7VpuEZEJMYU8iIiMaaQFxGJMYW8iEiMKeRFRGJMIS8iEmMKeRGRGFPIi4jk2Tfv2sS9m9oj2bZCXkQkj9ydb/1hE/dv2R3J9hXyIiJ51Ds4zHDaqauM5gYECnkRkTxKJIcAqKsojWT7WYe8mVWY2QNm9lcz22hmXwjnrzSz+81ss5ndEH7Hp4iIZEj0hyFfOUNDnuBb71/n7scRfPfnWWZ2KvBV4BvufhjwIsG33YiISIZEMgXM4J68B3rCl6Xhw4HXATeF868Dzst2XyIicTM6XDOTx+TNrNjMHgV2AncCzwCd7p4KF2kBlkyw7mVmtsHMNrS3R3MJkYjITDU6XDNTe/IA7j7s7scDS4GTgSMPYN1r3L3J3ZsaG8e9572ISGzt6cnP4JAf4e6dwB+B04C5Zjby98dSYPt07ktEJA4S/cGAR23FDB2uMbNGM5sbTlcCZwJPEIT9m8PFLgZuyXZfIiJxk0gOUVVWTGlxNFe0T8dbxyLgOjMrJnjTuNHdbzWzZuCXZvYl4BHgR9OwLxGRWEn0DzEnoqEamIaQd/fHgBPGmb+FYHxeREQmkEimIjvpCvrEq4hIXiX6hyK7fBIU8iIiedWVHFJPXkQkroKevEJeRCSWgjF5DdeIiMROOu10qycvIhJPvYMp0h7dLQ1AIS8ikjcjn3bV1TUiIjEU9ReGgEJeRCRvor45GSjkRUTyZnS4Rj15EZH4ifoLQ0AhLyKSN1F/YQgo5EVE8qYr7MlHdS95UMiLiORNIpmipryEkojuJQ8KeRGRvEn0D0V6SwNQyIuI5E0iGe0tDUAhLyKSN0FPXiEvIhJLiWQq0ssnQSEvIpI36smLiMSYxuRFRGIqnXa6B6L9whBQyIuI5EXPYAr3aG9OBtMQ8ma2zMz+aGbNZrbRzD4Yzq83szvNbFP4PC/7ckVE4iEXtxmG6enJp4CPuvsa4FTgCjNbA3wKuNvdVwN3h69FRITgyhqI9uZkMA0h7+5t7v5wON0NPAEsAdYC14WLXQecl+2+RETiomsW9eRHmdkK4ATgfmCBu7eFTS8ACyZY5zIz22BmG9rb26ezHBGRGWv0DpQzfUx+hJnVADcDH3L3RGabuzvg463n7te4e5O7NzU2Nk5XOSIiM9rImPyc2RDyZlZKEPA/c/dfhbN3mNmisH0RsHM69iUiEge5+FYomJ6rawz4EfCEu/97RtM64OJw+mLglmz3JSISFyM9+ZqIr5Ofjq2/Ang78DczezSc92ngK8CNZnYJ8BzwlmnYl4hILCT6h6gtL6G4yCLdT9Yh7+5/Biaq8vRsty8iEkfBzcmiHaoBfeJVRCQvEv1DkX7t3wiFvIhIHuTi5mSgkBcRyYtEfyryK2tAIS8ikhdBT17DNSIisZSLLwwBhbyISM4Np53ufl1dIyISS3tuM6zhGhGR2Nm0sweAlzTWRL4vhbyISI41t3YBsGZxXeT7UsiLiORYc1uC+dVlHFJbHvm+FPIiIjm2sTXBmsV1BPd3jJZCXkQkhwZTaTbt6GHNouiHakAhLyKSU8+09zA4nM7JeDwo5EVEcqq5NfjivKMV8iIi8dPclqCitIiVDdFfPgkKeRGRnGpuTXDEwrrIvyxkhEJeRCRH3J2NrV05O+kKCnkRkZzZ3pkk0Z/K2UlXUMiLiOTMyElX9eRFRGKouS2BGRy1qDZn+1TIi4jkSHNrgpUN1VSVRX/3yREKeRGRHGluS+R0qAamKeTN7Foz22lmj2fMqzezO81sU/g8bzr2JSIyG3X1DdHyYjKnJ11h+nryPwbOGjPvU8Dd7r4auDt8LSJSkO7Z1A7ACcty29+dlpB393uA3WNmrwWuC6evA86bjn2JiMxGNz3UwuI5FZy8sj6n+41yTH6Bu7eF0y8AC8ZbyMwuM7MNZrahvb09wnJERPLjha5+7t3UzptOXJqzT7qOyMmJV3d3wCdou8bdm9y9qbGxMRfliIjk1K8f2U7a4fyTluZ831GG/A4zWwQQPu+McF8iIjOSu3PTQ9toOnQeKxuqc77/KEN+HXBxOH0xcEuE+xIRmZEe3dbJM+29vDkPvXiYvksofwHcBxxhZi1mdgnwFeBMM9sEnBG+FhEpKDc91EJFaRH/46WL8rL/afnYlbu/dYKm06dj+yIis1H/0DDr/trK2ccsoraiNC816BOvIiIR+cUDz9Pdn8rbUA0o5EVEIvFCVz9fv+NpXn14Iy9/yfy81aGQFxGJwBd+u5Gh4TRfWnsMZrm9Nj6TQl5EZJrd1byD2x9/gQ+cvprl86vyWotCXkRkGvUOpPjcuo0cvqCGS1+1Kt/lTM/VNSIiAqnhNB++4VG2dya56fLTKCvJfz86/xWIiMRAOu184ubHuKN5B18492iaVuT2RmQTUciLiGTJ3bnq1mZ+9fB2PnLm4Vz88hX5LmmUhmtERLIwmErzxVubuX79c1zyypW8/3WH5bukvSjkRUQO0vbOJFf87GEe3dbJZa9exZVnH5nXyyXHo5AXETlA7s7vN+7gyl89xtCw890LT+TsY/Nzb5r9UciLiByAp3d0c9Vvm/nz5l0ctaiO77ztBFY11uS7rAkp5EVEpuCZ9h5+eO8WbtzQQk15CV8492guPGU5JcUz+/oVhbyIyATSaWf9lg6u/ctW7npiB+UlRVx0ynI+eMbh1FeX5bu8KVHIi4hkcHee3tHDLY9u5zePbKe1q5/66jI+ePpq3n7aoTTUlOe7xAOikBeRgpccHObBrbv5w5M7ufvJHWzbnaS4yHjV6gY+dc5RvH7NAipKi/Nd5kFRyItIwenoGeDRbZ088nwn67d08NeWToaGnfKSIl5xWAOXv+YlnLlmAYfUVuS71Kwp5EUktobTTsuLfTy9o4cn2hI80Zbg8dYutu1OAlBcZByzZA6XvHIVp66q55SV86ksm5099oko5EVkVhtOOzsS/bS8mOS5jl62dvSydVcfW3b1sqW9h4FUGgAzWDG/mmOXzOGiUw7l+GVzOXbpHKrK4h2D8T46EZnVBlNp2nsGaO8eYEeif/TR1tVPa2eS1s5+2rqSDA376DrFRcayeZWsaKjmlYfN57BDajjskFqOXFhLdXnhRV7hHbGI5EU67XQPpOjqG6IzOciLfUN09g3yYu8gu/uG2N07wO7eQXb1DNLRM0BH7yCdfUP7bKekyFhQV8GiORUct2wu5xy7iGX1lSybV8Wy+iqWzqukdIZfu55LCnkR2a+h4TR9A8P0DKbo6U/RM5CidyB47ulP0T2Qort/iO7+4DmRTNE9EDx3JYdI9A+RSA6R9vG3bwZzK0uZV11GQ005RyysZX51OY214aOmnAV1FSyYU05DdTlFRTPr/jAzWeQhb2ZnAd8EioEfuvtXot6nSKFxdwZSaZKDwySHwsfgMP0Z0yPPfftMp+gbHKZ3IJjuHRimb3BkXorewWEGw3Ht/akuK6a2opS6yhJqK0qZX1PGqsZq5lSWjj7mVpUxp7KUeVXB9MhzsYI7EpGGvJkVA98BzgRagAfNbJ27N0e5X5F8Sg2nGUil6R8aZiCVDh/DDAztPW9kun9oeM/8veal6U8N75keGqY/XGa89gNVZFBVVkJlWTFVZcVUlhZTXV5CbUUJC+sqqCoP5leXl1BTVkJVeQm15SVUl5dQXV5MbUUJNeWlwXR5KTUVJQrqGSjqnvzJwGZ33wJgZr8E1gIKeYnMSK92YCgM1zBk+4f2hO1o8IbL9Y/OzwjmcUJ67Hp7bTdcdniiMYkpqigtorykmIrSIipKiykvCZ4rSoqpqyihorY8fF1EZRjO5aXB8lWlxUFbafFo29jnqrLi0e3OtNviyvSLOuSXANsyXrcAp2QuYGaXAZcBLF++POJyJB/cncHh9F5DBsmM3ujer0em9+3NZobqXr3gkaANe7pTHVqYiBlUlBRTVlI0GrgjQVteUkRZSRG1FaWj0xUlxZSXFlFWvGeZzNcVpRNvayTMy0dfK3hleuX9xKu7XwNcA9DU1JRdF0iyNpx2ejNOro2cYAtOsgVjtb1jnvsGM5/Dsd6hFMnBNMnBFMmh4QlPuE2muMioKNkTlOWlQaCO9HDrKkupGBOmQcAW79UbHgnWsUE9drqitJiy4iJKi01BK7ERdchvB5ZlvF4azpOIuAeXqXX2BpepdfYN0ZUcojMZXN2QGLnSoT9FIrnnaojuMNT7BoentB8zqCotpqq8hOqyYirLSqgqC8ZpD6ktD8Z4w3kjwwMjz5WjQwpFwZhw6Z7grsiY1mVwItmLOuQfBFab2UqCcL8AeFvE+4ydgdQwOxMDox8Kae8eoKNnkF09A3T0BtO7ewd5sS8I9dQk3ebycKihrrKEuopSaitKWDK3MjyJVkLNyHM4XR1OV5cFJ9uqyoLXFaUaVhCZDSINeXdPmdn7gN8TXEJ5rbtvjHKfs427094zwLbdSVpe7KPlxSRtXUnaOvtp7Qo+3be7d3DcdedWlVJfXUZDdTmrGqupr65nXlUp86rKmJvxPHLpWjC8Ea/7cojI5CIfk3f324Dbot7PTJccHOaZ9h427exm044ent3Vy9aOPp7r6N1niGRuVSmL51SyeE4FJyyfy8K6Cg6pLeeQunIOqa2gsbac+uoyDWeIyH7l/cRrHPUMpHhsWyePtnTS3JqguS3Bs7t68XAUpaTIWF5fxYqGak5bNZ9D51exrL6SpfOqWDK3siDvryEi0VCaTIPdvYOs39LB+i0d3L9lN0/v7B4N9GX1laxZVMe5xy3miAW1rF5Qw/L6aspK1AsXkegp5A+Cu9PcluDuJ3Zy95M7+eu2TgCqyoo56dB5nH3sQk5YPo/jls5hbtXs+B5IEYknhfwB2La7j988sp1fP7KdLbt6MYPjls7lo2cezssPa+ClS+donFxEZhSF/H64O/c908E1927hT0+1A3DyynouffUqzjhqAY21s+tLfUWksCjkJ+Du3NG8g2/dvYmNrQkaasr40BmrOf/EpSyrr8p3eSIiU6KQH8dTL3Rz1a0b+cvmDlY1VPOVNx3LeScs0TXmIjLrKOQzDA2n+bc7nuKH9z5LTXkJXzj3aC48ZTklGmcXkVlKIR/akejnfT9/mAe3vsgFL1vGJ846kvpqXRkjIrObQh6475kO3v+Lh+kbHOabFxzP2uOX5LskEZFpUfAhv35LB+/8zwdYOq+SX1x6KqsX1Oa7JBGRaVPQIf9YSyfvuW4Dy+uruOF/nqbhGRGJnYI9o7h5ZzcXX/sAc6tKuf6SUxTwIhJLBRnyvQMpLr72QYqLivjpJaewcE5FvksSEYlEQQ7XfOPOp9nemeSmy09jRUN1vssREYlMwfXkH9/exbV/eZa3nbKcphX1+S5HRCRSBRXyw2nn07/+G/XV5XzyDUfmuxwRkcgVVMhff99WHmvp4rNvXMOcqtJ8lyMiErmCCfmegRRfv+NpXrW6gTe+dFG+yxERyYmCCfnb/tZG90CKD51xOGaW73JERHKiYEL+podaWNVYzYnL5+a7FBGRnCmIkH+uo5cHnt3Nm09aql68iBSUrELezP7RzDaaWdrMmsa0XWlmm83sKTN7Q3ZlZufmh1ooMnjTCUvzWYaISM5l+2Gox4E3Ad/PnGlma4ALgKOBxcBdZna4uw9nub8Dlk47Nz+8nVeubtQnW0Wk4GTVk3f3J9z9qXGa1gK/dPcBd38W2AycnM2+Dtb6LR1s70zy5pPUixeRwhPVmPwSYFvG65Zw3j7M7DIz22BmG9rb26e9kJseaqG2ooTXr1kw7dsWEZnp9hvyZnaXmT0+zmPtdBTg7te4e5O7NzU2Nk7HJkf1DqS47fE23njcYn0/q4gUpP2Oybv7GQex3e3AsozXS8N5ObWxNUH/UJozj1IvXkQKU1TDNeuAC8ys3MxWAquBByLa14Q2tnYBcPTiulzvWkRkRsj2Esp/MLMW4DTgd2b2ewB33wjcCDQD/w1ckY8ra5pbEzTUlNFYW57rXYuIzAhZXULp7r8Gfj1B29XA1dlsP1vNbQmOWlSnD0CJSMGK7SdeB1NpNu3oYY2GakSkgMU25J9p72FwOM3Ri+fkuxQRkbyJbchvbE0AsGaRevIiUrhiG/LNrQkqSotYqe9wFZECFt+Qb+viyIV1FBfppKuIFK5Yhry709ya0ElXESl4sQz57Z1JEv0pjceLSMGLZciPnHTVJ11FpNDFMuSbWxMUGRy5UCEvIoUtniHflmBlQzWVZbrzpIgUtniGfGuCNfoQlIhI/EK+q2+I7Z1JnXQVESGGIb+5vQeAIxbW5LkSEZH8i13IdyUHAZhfrdsLi4jELuQTyRQAdZWlea5ERCT/4hfy/UMA1FVkdat8EZFYiF/IJ4OQr61QT15EJHYh35UcorK0mLKS2B2aiMgBi10SJpIp6io1VCMiAnEM+f4h6jRUIyICxDXkdWWNiAgQx5BPpnRljYhIKKuQN7OvmdmTZvaYmf3azOZmtF1pZpvN7Ckze0PWlU6RevIiIntk25O/EzjG3V8KPA1cCWBma4ALgKOBs4D/MLOc3BIykRxijkJeRATIMuTd/Q53T4Uv1wNLw+m1wC/dfcDdnwU2Aydns68p1kOiP6UTryIioekck383cHs4vQTYltHWEs7bh5ldZmYbzGxDe3t7VgX0DQ4znHZdQikiEtpvGprZXcDCcZo+4+63hMt8BkgBPzvQAtz9GuAagKamJj/Q9TPtuaWBevIiIjCFkHf3MyZrN7N3An8PnO7uIyG9HViWsdjScF6kdHMyEZG9ZXt1zVnAJ4Bz3b0vo2kdcIGZlZvZSmA18EA2+5qKrqR68iIimbIdvP42UA7caWYA6939cnffaGY3As0EwzhXuPtwlvvar5Gbk2lMXkQkkFUauvthk7RdDVydzfYPlMbkRUT2FqtPvO7pySvkRUQgbiHfH5x4rdVtDUREgLiFfHKIqrJiSotjdVgiIgctVmmo2wyLiOwtXiGvLwwREdlLvEJePXkRkb3EL+R1ZY2IyKh4hby+MEREZC+xCvku3UteRGQvsQn5dNrp1nCNiMheYhPyvYMp0q5bGoiIZIpNyI982lWXUIqI7BGfkNdthkVE9hG/kNeYvIjIqPiE/MhwjXryIiKj4hPy+sIQEZF9xCfk9YUhIiL7iE/IJ3UveRGRseIT8v1DVJcVU6J7yYuIjIpNInYl9WlXEZGxYhPyiaRuMywiMlZ8Qr5/SFfWiIiMkVXIm9kXzewxM3vUzO4ws8XhfDOzb5nZ5rD9xOkpd2LBbYbVkxcRyZRtT/5r7v5Sdz8euBX4bDj/bGB1+LgM+G6W+9mvRL9uMywiMlZWIe/uiYyX1YCH02uBn3hgPTDXzBZls6/9SejEq4jIPrIexDazq4F3AF3A34WzlwDbMhZrCee1jbP+ZQS9fZYvX35QNaTTTveAvhVKRGSs/fbkzewuM3t8nMdaAHf/jLsvA34GvO9AC3D3a9y9yd2bGhsbD/wIgJ7BFO66OZmIyFj77fq6+xlT3NbPgNuAzwHbgWUZbUvDeZHQbYZFRMaX7dU1qzNergWeDKfXAe8Ir7I5Fehy932GaqbLyC0NdAmliMjesk3Fr5jZEUAaeA64PJx/G3AOsBnoA96V5X4mpZuTiYiML6uQd/fzJ5jvwBXZbPtAdOkLQ0RExhWLT7zOry7j7GMW0lhbnu9SRERmlFgMYjetqKdpRX2+yxARmXFi0ZMXEZHxKeRFRGJMIS8iEmMKeRGRGFPIi4jEmEJeRCTGFPIiIjGmkBcRiTEL7kAwM5hZO8E9cKaqAdgVUTkzWaEeNxTuseu4C8uBHveh7j7uvdpnVMgfKDPb4O5N+a4j1wr1uKFwj13HXVim87g1XCMiEmMKeRGRGJvtIX9NvgvIk0I9bijcY9dxF5ZpO+5ZPSYvIiKTm+09eRERmYRCXkQkxmZFyJvZWWb2lJltNrNPjdNebmY3hO33m9mKPJQ57aZw3B8xs2Yze8zM7jazQ/NR53Tb33FnLHe+mbmZxeISu6kct5m9JfydbzSzn+e6xqhM4d/6cjP7o5k9Ev57PycfdU4nM7vWzHaa2eMTtJuZfSv8mTxmZice1I7cfUY/gGLgGWAVUAb8FVgzZpl/Ar4XTl8A3JDvunN03H8HVIXT7y2U4w6XqwXuAdYDTfmuO0e/79XAI8C88PUh+a47h8d+DfDecHoNsDXfdU/Dcb8aOBF4fIL2c4DbAQNOBe4/mP3Mhp78ycBmd9/i7oPAL4G1Y5ZZC1wXTt8EnG5mlsMao7Df43b3P7p7X/hyPbA0xzVGYSq/b4AvAl8F+nNZXISmctyXAt9x9xcB3H1njmuMylSO3YG6cHoO0JrD+iLh7vcAuydZZC3wEw+sB+aa2aID3c9sCPklwLaM1y3hvHGXcfcU0AXMz0l10ZnKcWe6hOBdf7bb73GHf7Yuc/ff5bKwiE3l9304cLiZ/cXM1pvZWTmrLlpTOfbPAxeZWQtwG/D+3JSWVweaAeOKxRd5FzozuwhoAl6T71qiZmZFwL8D78xzKflQQjBk81qCv9ruMbNj3b0zn0XlyFuBH7v7183sNOB6MzvG3dP5Lmymmw09+e3AsozXS8N54y5jZiUEf8515KS66EzluDGzM4DPAOe6+0COaovS/o67FjgG+JOZbSUYq1wXg5OvU/l9twDr3H3I3Z8FniYI/dluKsd+CXAjgLvfB1QQ3MQrzqaUAfszG0L+QWC1ma00szKCE6vrxiyzDrg4nH4z8AcPz1zMYvs9bjM7Afg+QcDHZXx20uN29y53b3D3Fe6+guBcxLnuviE/5U6bqfw7/w1BLx4zayAYvtmSwxqjMpVjfx44HcDMjiII+facVpl764B3hFfZnAp0uXvbgW5kxg/XuHvKzN4H/J7gLPy17r7RzK4CNrj7OuBHBH++bSY4kXFB/iqeHlM87q8BNcB/heeZn3f3c/NW9DSY4nHHzhSP+/fA682sGRgGPu7us/0v1qke+0eBH5jZhwlOwr5ztnfkzOwXBG/aDeG5hs8BpQDu/j2Ccw/nAJuBPuBdB7WfWf5zEhGRScyG4RoRETlICnkRkRhTyIuIxJhCXkQkxhTyIiIxppAXEYkxhbyISIz9f5quit3r8VeTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "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": {},
   "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": 24,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_841398/3108092858.py\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/tmp/ipykernel_841398/3108092858.py\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/tmp/ipykernel_841398/3108092858.py\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": 27,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_841398/3077871989.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\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      2\u001b[0m \u001b[0;31m# Answer 1.\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": [
    "raise NotImplementedError\n",
    "# Answer 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_841398/2453162538.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# Answer 2.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\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[0;31mNotImplementedError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# Answer 2.\n",
    "raise NotImplementedError"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": 29,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_841398/1719254416.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     54\u001b[0m     \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtesting\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0massert_allclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb_true\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     55\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 56\u001b[0;31m \u001b[0mtest_nointercept\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[1;32m     57\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     58\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mgenerate_dataset_with_intercept\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msigma\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.7\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignal_strength\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\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;32m/tmp/ipykernel_841398/1719254416.py\u001b[0m in \u001b[0;36mtest_nointercept\u001b[0;34m()\u001b[0m\n\u001b[1;32m     51\u001b[0m     \u001b[0mCheck\u001b[0m \u001b[0mthat\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb_true\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0my\u001b[0m \u001b[0mwhen\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mnoise\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     52\u001b[0m     \"\"\"\n\u001b[0;32m---> 53\u001b[0;31m     \u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb_true\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgenerate_dataset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msigma\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignal_strength\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\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     54\u001b[0m     \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtesting\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0massert_allclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb_true\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     55\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/tmp/ipykernel_841398/1719254416.py\u001b[0m in \u001b[0;36mgenerate_dataset\u001b[0;34m(n, p, sigma, signal_strength)\u001b[0m\n\u001b[1;32m     45\u001b[0m     \u001b[0mb_true\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mnp\u001b[0m \u001b[0marray\u001b[0m \u001b[0mof\u001b[0m \u001b[0mshape\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mp\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[1;32m     46\u001b[0m     \"\"\"\n\u001b[0;32m---> 47\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     48\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     49\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mtest_nointercept\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 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 NotImplementedError\n",
    "\n",
    "def generate_dataset(n, p, sigma=0.7, signal_strength=2):\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",
    "    raise NotImplementedError\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 accout 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",
    "    raise NotImplementedError\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": {},
   "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": 33,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_841398/166140221.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     22\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_lsqr\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[1;32m     25\u001b[0m \u001b[0mtest_same_estimate\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;32m/tmp/ipykernel_841398/166140221.py\u001b[0m in \u001b[0;36mtest_lsqr\u001b[0;34m()\u001b[0m\n\u001b[1;32m     14\u001b[0m     \"\"\"Test that all methods yield the correct beta when there is no noise\n\u001b[1;32m     15\u001b[0m     \"\"\"\n\u001b[0;32m---> 16\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     17\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mtest_same_estimate\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 compute_least_squares_using_SVD(X, y):\n",
    "    raise NotImplementedError\n",
    "    \n",
    "def compute_least_squares_using_pseudo_inverse(X, y):\n",
    "    raise NotImplementedError\n",
    "    \n",
    "def compute_least_squares_using_method3(X, y):\n",
    "    raise NotImplementedError\n",
    "\n",
    "def compute_least_squares_using_scipy(X, y):\n",
    "    raise NotImplementedError\n",
    "    \n",
    "def test_lsqr():\n",
    "    \"\"\"Test that all methods yield the correct beta when there is no noise\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "        \n",
    "def test_same_estimate():\n",
    "    \"\"\" Test that all methods yield the same estimate of beta\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "        \n",
    "test_lsqr()\n",
    "test_same_estimate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": 34,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_841398/3955241513.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m     \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0merror_residual_plot\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/tmp/ipykernel_841398/3955241513.py\u001b[0m in \u001b[0;36merror_residual_plot\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0merror_residual_plot\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;32m----> 2\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      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0merror_residual_plot\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 error_residual_plot():\n",
    "    raise NotImplementedError\n",
    "    \n",
    "error_residual_plot()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": 35,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_841398/4148399174.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# Question 1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\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      3\u001b[0m \u001b[0;31m# Question 2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNotImplementedError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# Question 1\n",
    "raise NotImplementedError\n",
    "# Question 2\n",
    "raise NotImplementedError"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}