chapter-02/2.5-local-methods-in-high-dimensions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.5 Local Methods in High Dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Mean squared error for estimating f(0):**\n",
    "\n",
    "*Assume that the relationship between X and Y is: $Y = f(X)$*\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    " \\text{MSE}(x_0) & = E_\\tau[f(x_0) - \\hat{y_0}]^2\\\\\n",
    " & = E_\\tau[(f(x_0) - E_\\tau(\\hat{y_0})) + (E_\\tau(\\hat{y_0}) - \\hat{y_0})]^2\\\\\n",
    " & = E_\\tau[(E_\\tau(\\hat{y_0}) - \\hat{y_0})^2  + 2(f(x_0) - E_\\tau(\\hat{y_0}))(E_\\tau(\\hat{y_0}) - \\hat{y_0})+ (f(x_0) - E_\\tau(\\hat{y_0}))^2]\\\\\n",
    " & = E_\\tau[(E_\\tau(\\hat{y_0}) - \\hat{y_0})^2] + E_\\tau[(E_\\tau(\\hat{y_0}) - f(x_0))^2]\\\\\n",
    " & = E_\\tau[\\hat{y_0} - E_\\tau(\\hat{y_0})]^2 + [E_\\tau(\\hat{y_0}) - f(x_0)]^2\\\\\n",
    " & = Var_\\tau(\\hat{y_0}) + Bias^2(\\hat{y_0})\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "We have broken down the MSE into two components: variance and squared bias. Such decomposition is always possible and is known as *the bias-variance decomposition*.\n",
    "\n",
    "*(2.26) Suppose that the relationship between Y and X is linear with some noise:*\n",
    "\n",
    "$$Y = X^T\\beta + \\varepsilon $$\n",
    "\n",
    "where $\\varepsilon \\sim N(0, \\sigma^2)$ and we fit the model by least squares to the training data. For a test point $x_0$ we have $\\hat{y_0}=x_0^T\\hat{\\beta}$ which can be written as $\\hat{y_0} = x_0^T\\beta + \\sum_{i=1}^N {l_i(x_0)\\varepsilon_i}$ where $l_i(x_0)$ is the $i$th element of $\\mathbf{X}(\\mathbf{X}^T\\mathbf{X})^{-1}x_0$\n",
    "\n",
    "Proof:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    " \\hat{\\beta}=(\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T\\mathbf{y} \\\\\n",
    " \\hat{\\beta}=(\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T(\\mathbf{X}\\beta + \\varepsilon) \\\\\n",
    " \\hat{\\beta}=(\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T(\\mathbf{X}\\beta + \\varepsilon) \\\\\n",
    " \\hat{\\beta}=\\beta + (\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T\\varepsilon \\\\\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "and by plugging $\\hat{B}$ into the linear model:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    " \\hat{y_0} = x_0^T(\\beta+(\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T\\varepsilon) \\\\\n",
    " \\hat{y_0} = x_0^T\\beta+x_0^T(\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T\\varepsilon\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "we can get $\\mathbf{X}(\\mathbf{X}^T\\mathbf{X})^{-1}x_0$ from $(x_0^T(\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T)^T$ by using two matrix properties:\n",
    "\n",
    "1. $(\\mathbf{AB})^T=\\mathbf{B}^T\\mathbf{A}^T$\n",
    "\n",
    "2. $(\\mathbf{A}^{-1})^T = (\\mathbf{A}^T)^{-1}$\n",
    "\n",
    "\n",
    "\n",
    "Under this model the least square estimates are unbiased, so the expected prediction error will be:\n",
    "$$\n",
    "\\begin{align}\n",
    "\\text{EPE}(x_0) & = E_{y_0|x_0}E_\\tau(y_0-\\hat{y_0})^2\\\\\n",
    "& = \\text{Var}(y_0|x_0) + Var_\\tau(\\hat{y_0}) + \\text{Bias}^2(\\hat{y_0})\\\\\n",
    "& = \\sigma^2 + E_{\\tau}x_0^T(\\mathbf{X}^T\\mathbf{X})^{-1}x_0\\sigma^2 + 0^2\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "*Proof*:\n",
    "$$\n",
    "\\begin{align}\n",
    "\\text{EPE}(x_0) & = E_{y_0|x_0}E_\\tau(y_0-\\hat{y_0})^2\\\\\n",
    "& = E_{y_0|x_0}E_\\tau((y_0 - f(x_0)) + (f(x_0) - \\hat{y_0}))^2\\\\\n",
    "& = E_{y_0|x_0}E_\\tau(y_0 - f(x_0))^2 + 2E_{y_0|x_0}E_\\tau(y_0 - f(x_0))(f(x_0) - \\hat{y_0}) + E_{y_0|x_0}E_\\tau(f(x_0) - \\hat{y_0})^2\\\\\n",
    "& = U_1 + U_2 + U_3\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "There are three components $U_1$, $U_2$, $U_3$ and we're going to expand them as well. \n",
    "\n",
    "$U_1 = E_{y_0|x_0}E_\\tau(y_0 - f(x_0))^2 = E_{y_0|x_0}(y_0-f(x_0))^2 = \\sigma^2$ \n",
    "\n",
    "*Note*: $f(x_0) = E_{y_0|x_0}(y_0) $\n",
    "\n",
    "\n",
    "$U_2 = 2E_{y_0|x_0}E_\\tau(y_0 - f(x_0))(f(x_0) - \\hat{y_0}) = 0$ \n",
    "\n",
    "*Note*: $E_{y_0|x_0}(y_0-f(x_0)) = 0$\n",
    "\n",
    "\n",
    "*$U_3$*:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "U_3 & = E_{y_0|x_0}E_\\tau(f(x_0) - \\hat{y_0})^2\\\\\n",
    "& = E_{y_0|x_0}E_\\tau((\\hat{y_0} - E_\\tau(\\hat{y_0})) + (E_\\tau(\\hat{y_0}) - f(x_0)))^2\\\\\n",
    "& = E_{y_0|x_0}E_\\tau(\\hat{y_0} - E_\\tau(\\hat{y_0}))^2 + 2E_{y_0|x_0}E_\\tau[(\\hat{y_0} - E_\\tau(\\hat{y_0}))(E_\\tau(\\hat{y_0}) - f(x_0))] + E_{y_0|x_0}E_\\tau(E_\\tau(\\hat{y_0}) - f(x_0))^2\\\\\n",
    "& = E_\\tau(\\hat{y_0} - E_\\tau(\\hat{y_0}))^2 + (E_\\tau(\\hat{y_0}) - f(x_0))^2\\\\\n",
    "& = \\text{Var}_\\tau(\\hat{y_0}) + \\text{Bias}_\\tau^2(\\hat{y_0}) \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Finally if we sum all $U_i$ we get: $$\\text{EPE}(x_0) = U_1+U_2+U_3 = \\sigma^2 + 0 + (\\text{Var}_\\tau(\\hat{y_0}) + \\text{Bias}_\\tau^2(\\hat{y_0}))$$\n",
    "\n",
    "$E_\\tau(\\hat{y_0}) = E_\\tau(x_0^T\\beta + \\sum_{i=1}^N {l_i(x_0)\\varepsilon_i})=x_0^T\\beta + E(\\sum_{i=1}^N {l_i(x_0)\\varepsilon_i}) = x_0^T\\beta + 0$ thus $\\text{Bias}_\\tau{\\hat{y_0}} = 0$ \n",
    "\n",
    "(2.27) and we can find variance: \n",
    "$$\n",
    "\\begin{align}\n",
    "\\text{Var}_\\tau(\\hat{y_0}) & = E_\\tau(\\hat{y_0} - E_\\tau(\\hat{y_0})) ^ 2\\\\\n",
    "& = E_\\tau(x_0^T(\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T\\varepsilon)\\\\\n",
    "& = E_\\tau(x_0^T(\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T\\varepsilon\\varepsilon^T\\mathbf{X}(\\mathbf{X}^T\\mathbf{X})^{-1}x_0)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $\\varepsilon\\varepsilon^T=\\sigma^2\\mathbf{I}_n$, so we can simplify further:\n",
    "$$\\text{Var}_\\tau(\\hat{y_0}) = \\sigma^2x_0^{T}E_\\tau[(\\mathbf{X}^T\\mathbf{X})^{-1})]x_0$$\n",
    "\n",
    "(2.28) if N is large and $\\tau$ were selected at random, and assuming E(X) = 0, then $\\mathbf{X}^T\\mathbf{X}$->$NCov(\\mathbf{X})$.\n",
    "\n",
    "*Proof*: By definition of covariance $\\text{Cov}(X) = E[(X-E(X))(X-E(X))^T] = E(XX^T) = \\frac{\\mathbf{X}^T\\mathbf{X}}{N}$\n",
    "\n",
    "and we can derive that:\n",
    "$$\n",
    "\\begin{align}\n",
    "E_{x_0}\\text{EPE}(x_0) & = E_{x_0}x_0^{T}\\text{Cov}^{-1}(X)x_0\\sigma^2/N+\\sigma^2\\\\\n",
    "& = \\text{trace}[\\text{Cov}^{-1}(X)\\text{Cov}(x_0)]\\sigma^2/N+\\sigma^2\\\\\n",
    "& = \\sigma^2(p/N)+\\sigma^2\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "**FIGURE 2.9**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAHwCAYAAAC/hfaiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XeYk1X6xvHvYRjKANKL1EGKKAgKoyKgYu+9O4ruqii2tayVXctvxbL2shZcFcuoa8W1txUQ6SAiKgJKR3pnGJxyfn88GRmGaUCSN3lzf64r10ySN8mTUO6c857ivPeIiIhI8qsWdAEiIiISHQp1ERGRkFCoi4iIhIRCXUREJCQU6iIiIiGhUBcREQkJhbpIKc65C51zoxPhtZ1zG5xzu+3A82Q75z6LbnXJwzl3l3NuhXNuSdC1iMSTQl12mnNurnNuUySAii9PBFjPCOfcxTF67kznnC/xPuc6526OxWsBeO/reu9/rWJN1Us8Lsd7f2S063HO9XfOFZX6s97gnDsg2q+1o5xzbYHrgT299y3KuL+Gc+6tyJ+dd871L3X/sMjt+5W4raNzzpe4PsI5l+eca1PitsOdc3Nj8Z4qEvki+L1zLtc5t8Q595RzrkGJ++9wzr1SzmP7OefGOOfWOudWOee+cc7tG7/qJdoU6hItJ0QCqPhyZdAFxVgD731d4BzgNufc0aUPKBmyIbO41J91Xe/92NIHOVOtstsqswOfY1tgpfd+WQXHjAbOA8prya8C7qrkdTYCf9/O2qLKOXc9cB9wA1Af6A20Az53ztWo5LG7AB8AjwONgFbAncDmWNYssaVQl5iKtBreLnH9Pufcl5H/3Ps75xY6526NdJXOdc5llzi2pnPuAefcfOfcUufc08652iXuP8k5N9U5t84594tz7mjn3BDgQOCJkj0GzrkuzrnPI62Rn51zZ5Z4nsbOuf9GnmcC0KGq7y8SZj8A3SLP5Z1zVzjnZgGzdva1I8/XMfJ7befcg865eZGW1ejI5zEqcvia4laz27Ybv49zbmLkcROdc31K3DfCOfePSCttvXPuM+dck6p+BqXqHeGcG+Kc+wbIBXYr57aWkfe9yjk32zl3SYnnuMNZS/oV59w64MIyXqe+c+4l59zyyOfxN+dcNefc4cDnQMvIZzGs9GO997977x/x3o8GCst5Ky8C3Z1zB1fwdh8DznHOVfr3JfLv4IFSt73nnLsu8vtNzrlFkc//Z+fcYVV4zl2wEL7Ke/+J9z7fez8XOBPIxL60VKQzgPf+Ne99ofd+k/f+M+/9tMpeWxKY914XXXbqAswFDi/nvgxgJvYf84HACqB15L7+QAHwEFATOBhr/eweuf9h4L9YK6Ie8D5wT+S+/YC1wBHYl9NWQJfIfSOAi0vUUAdYAPwJqA7sE6ljz8j9rwNvRI7rBiwCRpfzfjIBH3keB/TFguqwyP0eC5VGQO2dfe3I83WM/P6vyHtrBaQBfSKf2x81lXjchcXPE6llNXB+pIZzItcbl/i8fsH+k68duX5vOe+/P7Cwgr8LI4D5QNfIa6WXc9so4EmgFrA3sBw4NPIcdwD5wMmRP9vaZbzOS8B7kb8XmdjfsYuqUmOp51kI9C912zCslX51ic+wI+BLvc+Lsb+7r0RuOxyYW87rHBT5e+Ai1xsCm4CWwO6R+1qW+DvWoQq1H439+6lexn0vAq+V+DxfKeOYXYCVkWOPARoG/X+JLjt/UUtdomW4c25NicslAN77XCxMHgJewVoVC0s99u/e+83e+5HAh8CZzjkHDASu9d6v8t6vB+4Gzo485iLgee/95977Iu/9Iu/9jHJqOx77z/YF732B9/5b4G3gDOdcGnAacJv3fqP3fjr2n1xlVmBdtP8Gbvbef1nivnsiNW+K1mtHuqz/DPwl8l4LvfdjvPdV6So9DpjlvX85UsNrwAzghBLHvOC9nxmp+Q0saMvTstSf9RrnXJ0S9w/z3v8Qea380rcBLbAvQzd57/O891Mjn+OAEs8x1ns/PPJnu6nUZ5GG/T24xXu/3lvr9EHs71k0PQO0dc4dU8Ex9wAnOOe6VvJcX2NfvA6MXD8de4+Lsd6CmsCezrl07/1c7/0vVaivCbAi8pmW9lvk/nJ579cB/SJ1PQssj/SeNK/Ca0uCUqhLtJzsvW9Q4vJs8R3e+/HAr1jL9o1Sj1vtvd9Y4vo8rPXSFGvlTy4ODuCTyO0AbbDWZVW0A/YvGUJANhYuTbHW44JSNVSmife+ofd+D+/9Y6XuK/lc0XrtJlirtqrvuaSWZTzvPKzFX6zkueVcoG4Fz7e41J91g1J/hgvKeEzJ21oCxV/UyqunrOco1gRr7Zd8T6Ufv9MiX5j+EbmUd8xy4Ang/yp5Lo/1ypwTuelcICdy32zgGqxFvcw597pzrmUVSlwBNHFljznYNXJ/hbz3P3nvL/Tet8Z6iloCj1ThtSVBKdQl5pxzV2AtkcXAjaXubliqldc2ctwKrHuya4ngqO9tcBrYf/rlncssvfXgAmBkqRCq670fhHX7FmBfEkrWsDNKvn60XnsFkEfZ77myrRYXY18uSmqLdfXHQln1lLxtMdDIOVevgnoqek8rsO75ku8pVu/nBaABcGoFx9wPHAL0quS5XgNOd861A/bHemwA8N6/6r3vh70njw1+q8xYbFDbVrU55+pi3elflvWg8kR6uoYRGR8iyUmhLjHlnOuMnZ88D+sevdE5V7pr905n04wOxLqr3/TeF2Fdgg8755pFnquVc+6oyGOeA/7knDssMkCqlXOuS+S+pUDJud0fAJ2dc+c759Ijl32dc3t47wuBd4A7nHMZzrk9gQui+BFE5bUjn8fzwEORQWZpkQFxNbEvB0Wl3nNJH0VqONc5V905dxawZ6S2uPPeLwDGAPc452o557pjp1PKnHZVxuMLsR6fIc65epGQvK6qj4c/BmHWilytEanDlfFaBcDtwE0V1LMG6/4v/YW19HHfYl9I/g18GnkczrndnXOHRv4s87Avs0WVvQfv/VpsoNzjzgaJpjvnMrHPZiHwconDq0XeY/GlprMBnNc751pH6miD9SSMq+y1JXEp1CVa3ndbz1t+N9It+Apwn/f+O+/9LOBW4OXIf2Bg3b6rsdZbDnBZiXPjNwGzgXGRUdBfYIOK8N5PwAafPYwNmBvJlpbbo1iLaLVz7rFIN++R2HnYxZHXvA/rPQC4EutuXoK1VF6I1ocS5df+K/A9MBE7n38fUC0ybmEI8E2ki793qRpWYl+WrscGRt0IHO+9r7R7thwt3bbz1E/bzuc4BxsQthh4F7jde//Fdjz+KmxQ5a/Y9LRXsS89VfUzFp6tgE8jv5fuzSj2GnaOuiKPUv5I+pJexQbUvVritprAvVjgLwGaAbfAH4sI/VDek3nv/4n9m3oAWAeMx3qHDis13uIc7D0WX34B1mM9BuOdcxuxMJ+O/T2RJFU8ElMk7pwt+vFK5HyeiIjsJLXURUREQkKhLiIiEhLqfhcREQkJtdRFRERCQqEuIiISEkm3i1STJk18ZmZm0GWIiIjEzeTJk1d475tWdlzShXpmZiaTJk0KugwREZG4cc5VZflqdb+LiIiEhUJdREQkJBTqIiIiIZF059TLkp+fz8KFC8nLywu6lEDVqlWL1q1bk56eHnQpIiISgFCE+sKFC6lXrx6ZmZmUsdFSSvDes3LlShYuXEj79u2DLkdERAIQiu73vLw8GjdunLKBDuCco3HjxinfWyEikspCEepASgd6MX0GIiKpLTShHrS6detuc9vTTz/NSy+9FEA1IiKSikJxTn27zcmB7wZD7nzIaAs9hkD77Ki/zGWXXRb15yzJe4/3nmrV9N1MRERSsaU+JwcmDITceYC3nxMG2u1Rdscdd/DAAw8A0L9/f2666Sb2228/OnfuzNdffw1AYWEhN9xwA/vuuy/du3fnmWeeAWDDhg0cdthh9OzZk7322ov33nsPgLlz57L77rszYMAAunXrxoIFC6Jet4iIJKfwtdQnXwOrp5Z//4pxULR569sKc2H8RfDLs2U/puHe0OuRnS6toKCACRMm8NFHH3HnnXfyxRdf8Nxzz1G/fn0mTpzI5s2b6du3L0ceeSRt2rTh3XffZZdddmHFihX07t2bE088EYBZs2bx4osv0rt3752uSUREwiN8oV6Z0oFe2e1RdOqppwLQq1cv5s6dC8Bnn33GtGnTeOuttwBYu3Yts2bNonXr1tx6662MGjWKatWqsWjRIpYuXQpAu3btFOgiIrKN8IV6ZS3q4ZmRrvdSMtrB4SNiUdEfatasCUBaWhoFBQWAnRd//PHHOeqoo7Y6dtiwYSxfvpzJkyeTnp5OZmbmH9PV6tSpE9M6RUQkOaXeOfUeQyAtY+vb0jLs9gAcddRRPPXUU+Tn5wMwc+ZMNm7cyNq1a2nWrBnp6el89dVXzJtXpQ16REQkhYWvpV6Z4lHuUR79npubS+vWrf+4ft1111XpcRdffDFz586lZ8+eeO9p2rQpw4cPJzs7mxNOOIG99tqLrKwsunTpslP1iYhI+DnvfdA1bJesrCxfej/1n376iT322COgihKLPgtJKnGaXiqS7Jxzk733WZUdl3otdRFJDMXTSwtz7Xrx9FJQsIvsoNQ7py4iwfAeNv0Giz+FH++HCZduCfRihbnWcheRHaKWuohEX8EmWPcjrJ4Ga0pcNq+o/LG582DVZGjUK/Z1ioSMQl1Edpz3dj68dHivnwm+yI5Jqw0N9oLWJ0OD7pHLXvBxz7Knl+Lgkyxo3Bs6Xwltz4C0GnF9W6GjsQspQ6EuIlWTvx7WTN86vNdMg/x1W46pu5uFdtsztwR43d2gWtq2z9djyNbn1MGml/Z61G6b+QSMPQ++vR46DoSOl0JGq9i/z7DR2IWUolAXka35Ilj/y7bhveHXLcek72KBnXleidZ3N0ivV/XXqWx6aecr4bfPLdyn3wU/3ANtTrXbm/YDbTVcNd8NLn/sgkI9dBTqUbRkyRKuueYaJk6cSIMGDWjevDmPPPIInTt3LvP4zMxMJk2aRJMmTba6/emnnyYjI4MBAwbEo2xJZb+v3rbrfM30LSHgqkG9TnZ+e7c/bQnwOu2iE6rts8sPFlcNWh5llw2/wswn4ZfnYP4b0KCHhXvmuVA9o+zHp7pNS2Dh8HJOcWBfpLzXl6OQSclQz8mBwYNh/nxo2xaGDIHsnfzC6r3nlFNO4YILLuD1118H4LvvvmPp0qXlhnp5Yr1lq4RURedNi/Jh3cxtW9+5C7c8vkYjaNgDOl6yJbzr75kYoVl3N+j5AHT/P5ibAzMfhwmXwNQbocNF0GmQHZPqNs6HBe/Agrdh+TeAB1cdfEEZB3v4pBd0uhwyz4HqWn46DFIu1HNyYOBAyI00RObNs+uwc8H+1VdfkZ6evlUg9+jRgxEjRnD88cfzwQcfAHDllVeSlZXFhRdeCMA///lPPv74Y2rXrs2rr75Kx44dueOOO6hbty5//etfmT17NpdddhnLly8nLS2NN998kw4dOux4oRJOZZ03HfcnmPU0FG6AtT9C0e92n6sO9feAZgeX6DrvDrV3TfxWW/UM+9LR4WJYPtrCfcbD8NOD0Op4a723ONxa+ali3SwL8QVvw6rIwlwN9oK9boc2p8Hq78oYu1Ab2p4Fqyfbl6Nv/wrtB9iXo/pavCqZhS7Ur7kGplaw8+q4cbC51IZsublw0UXwbDk7r+69NzxSyT4x06dPp1ev7Z+CU79+fb7//nteeuklrrnmmj/Cv1h2djY333wzp5xyCnl5eRQVFW33a0gKKOu8qc+HlWOgxRF2KQ7vXbok/2hy56DZgXbJXWRfXn4ZCl+9D7vsDp2ugN0usHP/YeM9rJ0O89+Ghe/Amu/t9kb7wt73QutTYZdOW45v0M1+ltWL4z2sGGOnNmY/Y1+Smh1srffWJyf/35MUFLpQr0zpQK/s9lg755xz/vh57bXXbnXf+vXrWbRoEaeccgoAtWrVint9kiRy55d9u/dwyCfxrSXeMlpBj39At7/B/LcsmCZfDd/dCu0vgM5XJH/r03trhS9427rX188CnA0Y7PkItDkF6rQt//HljV1wDpr2tUveI/Dr8zDrGfjmLKjV3HpEOg6s+LlliwSYOhi6UK+sRZ2ZaV3upbVrByNG7Pjrdu3a9Y890UuqXr36Vq3r4u1Ti7kS3Z0u0bs+JTHNfxsoZw+HjBT6zzit5pbwWjkRZv4LfnkWZv3LuuQ7Xwktjy97el0iKiq0VnRxkOcusFMnzQ+BLtdbS7p28+i9Xq2msOdNsMcN8Nun1nr/4W748R5oeZy13nc9MrVObWyPBJk6mHJ/OkOGQEapcT8ZGXb7zjj00EPZvHkzQ4cO/eO2adOm4b3nxx9/ZPPmzaxZs4Yvv/xyq8f95z//+ePnAQccsNV99erVo3Xr1gwfPhyAzZs3k5tbqotVUpf3ttzq6DOgbkc7T1pSgFsKB67xvnDAMDh5IfS4G9bNgFEnw/sd4Md/wuaVQVdYtqJ8m8Y3YRAMbwVfHGSnFhruDb2HwalL4dDPoNOl0Q30klw1aHkM9H8fTpoDe94MK8fDiGPg/U72+eVVYWXAVFKUb+spJMCyxykX6tnZMHSotcyds59Dh+786HfnHO+++y5ffPEFHTp0oGvXrtxyyy20aNGCM888k27dunHmmWeyzz77bPW41atX0717dx599FEefvjhbZ735Zdf5rHHHqN79+706dOHJUuW7FyhEg5F+TDxMhv93fYMOHYa7PcsZLQDnP3cb6jmIddqCl1vgRPnwIFv2wj5qTfB8NYw7iJY9W3QFUJhHix83wY2vtMCvjoS5r4MzQ6Cvq/Dacvh4P/aGIGajeJbW5129sXwpAVWS0brLZ/fmPNh+Rj7cplqvIc1P8CMR2HECfBWY8hbWvax5Z0aixFtvRoy+ixSQP46+PoMWPIZ7HkL9LhLXaLbY81065qf85K1pJr0gc5X2cI28RoYlr8BfvvYutUXfQAFGyC9PrQ6EdqeBi2OhOq1K3+eIKz5AWY/bZ9f/jobfNlpEGRmb9/iQ8kmdyEs+QKWfGk/8yINrLod7fTOgrfK3tsgox2cPHenX76qW68q1ENGn0XIbZwPI4+HtT/Bfs9Ahz8HXVHy+n0N/DrMAn7DbKjVAjpdZgPDau8am9db9L4F+W+fWAu9ZlM7N97mNDtXnkyjzfM3wLzXYNaTsHoqVK8H7c+3gC8ecZ/Mfl8DS7+yEF/6Baz72W6v2dRCvMXh0OIw682Abc+pg50Ci1KPmUI9RemzCLFVk2HE8VC4ybqSWxwWdEXh4IsiA8OegMUf2WC0tqdb673JATs3dz9vOSx8zwa7Lf3STpvUbmW9Am1Os9HryTJwrzze2zn3WU/BvP9A0WZ7X50G2XtMqxl0hVVTmAcrxkZa41/YbANfZIvyNDvYQrz5YfaFpbyesRiOfleopyh9FiG18D345lyo1Qz6f2grvUn0rZ9to75/fR7y10LDfSzc251d9e7w3EWw4N3Iqm6jLBjq7mYB1+Y0G8QX1tMlm1da78esp2DDL9aq7XCRbcZTNzPo6rbmi6yHoTjEl39twe7SoElvC/AWh0Pj/ROiByXlQr1Lly4pPyXMe8+MGTMU6mHiPfz8KEy5zsLgoP/GbtSzbJG/IbIc7RO20EvNxjZnu9MgWDZ629ZY0wOsW33+27BynD1H/T23BHmD7om/Wl80+SILyllPwaL/2t/jlsfa57fr0cH0TnhvewgUh/jS/8Hvq+y++l23dKk3OyghFy1KqVCfM2cO9erVo3Hjxikb7N57Vq5cyfr162nfvn3Q5Ug0FBXA5GtsnnWb0+CAlxJjHfZU4j0sG2nhvnA4+EJryfnCEgc5/lgnoGFPG+jW+lSo3yWIihPPxgW2XsDsZ21wWZ121nLvcJH1PMVS3jJY8j87J77kS9g4127PaB3pTj8cWhwamzEUUZZSoZ6fn8/ChQu3Wdgl1dSqVYvWrVuTnp4edCmys/LXwzdn2znePW6w5T/D2mWbLDYugI+6bb1/fLH0hnDMZKirL9TlKsq300iznrJWcrV0+7La6fLobaVbsBGWjdoyQn3Nd3Z7en1ofuiWwW31Oiddz0lKhfoOSYDl/ETKlLvQBsStnQ5Z/7KFRiQxvFqNslfvc3Cu9mWosrUzbFrcr8Ns7EL9rtY13/787ev6Lsq31QOLR6ivGGu3VathXxSKB7c16pX0AxIV6hWJ8dQDkR226lubspa/Hvq9aXuJS+IYnln2/uRRmouccgpyYd7r1npfNclGmmeeZwG/Zvq2Da/Mc2HdTyXOi4+AgvWAs0GNxefFm/YN3akqhXpF9A9TEtGiD6zLvUYjG+HeYK+gK5LS1CCInZUTI9PiXrNR6FQDSvR+uDSoXtda9gB1O2wJ8eaH2GDGEKtqqIduQ5cqKW/Zvjgv5yfyh58fhynXWGvj4PeTYuBOSioObp26i77G+9ql54PwXgfIX731/b7Qutb3/7d1qSfaFLkEkZqhntG2nJZ6Cu1oJYmhqNA2gvj5UWh9EvTJsS5ISVzlbWMq0VGjIeSvKfu+wk02al7KlZrDaXsMsS6zrTjodlsg5UiKyt8AX59igd7lOuj3tgJdBMpvYKnhVanUDPX22XYOrHhHq1rNAQ+rxgddmaSK3MXwxcGw+EMb4d7zwaQfnSsSNWU1vFJ5K+HtkJrd77BtF9q3N8JP90PL46H1CcHVJeG3ehqMPM42jDjofWh1bNAViSQWjV3YYak5+r0shZvh0/1h02I49nstxSmxsfgTGH2GLYbR/wNouHfQFYlIEqjq6PfU7H4vS1pN6POKrRY1/mJbHlIkmmY9ZXPQ63WEo8Yr0EUk6hTqJTXoBnvfB4s/gNlDg65GwqKoEKZcDxMvt80sDv8aMloFXZWIhJBCvbTdr7LFDKZcB+tmBl2NJLuCjTD6dJjxkG3hedB7kF436KpEJKQU6qW5atB7mHXHjznPFjsQ2RGblsAX/W3ryV6PQtZjGuEuIjGlUC9LRivY7xlYNRGm/yPoapLfnBxbmvfVavZzTk7QFcXemuk28HLtj3DgcNj96qArEpEUoFAvT9szoP0A+GEILB8bdDXJq3it7Nx5gLefEwaGO9h/+xw+7ws+H474WlMkRSRuFOoVyXrc5keOPc92zZLt993grTe/ALv+3eBg6om12c/CiGOgTjs4cjw06hl0RSKSQhTqFUnfBQ54CTbMgSnXBl1Ncip385x58NtntlRqGPgimHqz9UK0OAKOGA112gRdlYikmNRdUa6qmh0Ie94MP94DLY+DNqcEXVFyqdEQfl9V9n1fHQWuuu3M1Kw/NO8PTfok3+jwgk0wdgAseAs6XmY9PNX0T0tE4k//81TFXnfAb5/ChEugSW9ti1lViz+B31ezzb7IaRnQ63Go0xqWjoBlI2yJ3h/vsZBvlGUB36w/NO2b2CGftwxGnggrJ8A+D0KXa8G5oKsSkRSlZWKrau1P8ElPC5r+H+k/7sqs+ha+OMhWT+t0pc0iqGgN54KNsHyMBfzSERaSvgBcGjTaF5ofXCLk6wXwhsqw9kcYcRzkLbUtU9WLIyIxUtVlYhXq22Pmv2DSlZD1BHS+IpgaksHG+fBZb3DpcORYyGi5/c9RYciXbskHEPJLvoSvT4O0WnDw+3YKQUQkRhTqseA9jDjWguboKVB/j2DqSGS/r4bP+0HuIjjiG2jQNTrPW7ARVozd0l2/coItDOTSoFGvLefkm/aLfcj/8jxMuBR22R36f2gj3UVEYkihHiubfoOP9rK92I8cC2k1gqsl0RRutsFvK8bAIZ9ZyMZKQW6pkB9fTsj3tVkM0eCLYNrf4Ye7bYR7vzehRv3oPLeISAUU6rG04F34+lToequdHxYLvDHnwbzX7Pxy5rnxff1yQ74aNOy1pbu+Wb8dC/nCPBh7Icz/D3S4GPZ9EqqlR/UtiIiUR6Eea+Mugl9fgMNH2rS3VDf1FvjxXuhxD3S9OehqIiE/rsQ5+fFQ9HuJkC8eeNev8tZ23nIYdbL1QOx9H+xxgwZKikhcBR7qzrk2wEtAc8ADQ733j5Y6xgGPAscCucCF3vspFT1vwoR6/nr4eG/whXDstOh18SajWU/DxEHQaRBk/SsxA69gE6wct6Ulv2JciZDvWWLgXSTk5+TYqne5861LH6Dva9D29ODeg4ikrKqGeiznqRcA13vvpzjn6gGTnXOfe+9/LHHMMUCnyGV/4KnIz8SXXg8OeAW+6AeTroYDhgVdUTAWvg+TroCWx0OvxxIz0AGq14bmh9gFtg35nx+Dnx6wkM9oC7kLbbQ92M9qtWzMgIhIAotZqHvvfwN+i/y+3jn3E9AKKBnqJwEveesuGOeca+Cc2zXy2MTX9ADoOtjmYLc6PvVacSsnwjdnW0u33+vJtYpamSE/3kL+x3u3BHqxojxruZeeXy8ikkDisva7cy4T2AcYX+quVsCCEtcXRm4r/fiBzrlJzrlJy5cvj1WZO6bb321xlAmX2jSuVLHhVxh5PNRqDgd/ANXrBF3Rzqle27rgu99h3fJlKW8dexGRBBHzUHfO1QXeBq7x3q/bkefw3g/13md577OaNm0a3QJ3VrV06POKjY4e9ycbBR52m1fCV8dAUQEc8jHUbh50RdGV0Xb7bhcRSRAxDXXnXDoW6Dne+3fKOGQRUHIrq9aR25LLLp2h50Ow5HP4+fGgq4mtgk0w6iTYOA8Oes8WYAmbHkNsffqS0jI0fVFEEl7MQj0ysv054Cfv/UPlHPZfYIAzvYG1SXM+vbSOA22w2NSbYM0PQVcTG77IdiNbPgb6vGxzvsOofTbsN9QWGMLZz/2G6ny6iCS8WE5p6wd8DXzPli26bgXaAnjvn44E/xPA0diUtj957yucr5YwU9rKsmmprTZXuyUcNR7SagZdUXRNuR5mPGS7ke1xXdDViIikjMCntHnvRwMVzm+KjHoPz84otZvD/s/BqBNh2m2wz31BVxQ9Pz9mgd75atteVEREEk5cRr+nlNYnWFf8T/fD0pFBVxMdC96ByddA61Ns7ECizkUXEUlxCvVY6PmQ7SM+dgD8viZ9VbQKAAAgAElEQVToanbO8jEwJhsa729ruldLC7oiEREph0I9FqrXsdXmNi2y/deT1bpZdiqhdms4+L82l1tERBKWQj1WmuwH3W6DuTkw9/Wgq9l+ectgxDGAs7notRJsfQAREdmGQj2Wut4KjXvbZicbF1R+fKIoyIWRJ8CmxbZaXL2OQVckIiJVoFCPpWrVbbU5nw/jLkiO1eaKCmHMubaue59XoUly7K8jIiIK9dir1wF6PQpLv4IZjwRdTcW8h8l/gYXv2Y5rbU4OuiIREdkOCvV42O3P0Ppk+O4WWD0t6GrKN+NBmPUv2OOvsHsSD/ATEUlRCvV4cM6WGa3R0KaHFeYFXdG25v0Hvr0B2p4Je4do0RwRkRSiUI+XWk2h9wuwdrrty51Ilo2yOfVND4QDXgSnvxYiIslI/3vHU8tjoNPlttzqki+Drsas/QlGngR1d4ODhkNaraArEhGRHaRQj7d97rftSsdeAL+vDraWTUtsLnpaTej/MdRsFGw9IiKyUxTq8VY9w5ZbzVsKEwbZiPMg5G+AEcdB3nLo/yHUzQymDhERiRqFehAa9YLud8L8/8DcV+P/+kUF8M1ZsGYq9HvD6hERkaSnUA/KHjdB074w6XLYOC9+r+s9TLwcFn8E+z4FrY6L32uLiEhMKdSDUi0NDnjZQnbsAFvJLR5+vAd+edaWsO04MD6vKSIicaFQD1Ld9pD1uE0pm/Fg7F9vzis2nS7zPOh+V+xfT0RE4kqhHrT2A6DNaTDtb7Dq29i9zpL/wfg/Q/NDYP/nbEEcEREJFYV60JyD/Z6Bmk1g7HlQsCn6r7Hme/j6FKjXGQ58B9JqRP81REQkcAr1RFCzMfQeBmt/hKk3R/e5cxfBiGOhel2bi16jQXSfX0REEoZCPVHseiR0vhpmPga/fRad58xfZ4H++1ro/xHUaROd5xURkYSkUE8ke98L9feEcRfC5pU791xF+fD16db6P/AtaNgjKiWKiEjiUqgnkuq1bbW5zStgwqU7vtqc9zD+EljyOez/rPUCiIhI6CnUE03DvW262YK3Yc5LO/Yc398Jc16Eve6A3S6MZnUiIpLAFOqJqMv10OxgmHQVbJizfY/95XmYfifs9ifodlts6hMRkYSkUE9E1dIi+5o7GHt+1VebW/wpTBgILY60aXKaiy4iklIU6omqTjvIehKWfwM/3Vf58au+hdGnQ/1ucOCbUC099jWKiEhCUagnssxzoe1ZMO12WDmp/OM2zoeRx0GNhjZ1LX2X+NUoIiIJQ6GeyJyD/Z6C2i0iq83lbnvM72tsLnpBri0uk9Ey/nWKiEhCUKgnuhoNofeLsO5n+PaGre8r3AyjToH1M+Ggd6FB12BqFBGRhKBQTwYtDoUu18GsJ2HRR3abL4Jxf4ZlI2D/F2yjFhERSWnVgy5AqqjHEFtM5ptz7Jz5poV2e5szoH12sLWJiEhCUEs9WaTVgnbnQsG6LYEOsPhDmJMTXF0iIpIwFOrJZNbT295WmAvfDY5/LSIiknAU6skkd/723S4iIilFoZ5MMtpu3+0iIpJSFOrJpMcQSMvY+ra0DLtdRERSnkI9mbTPhv2GQkY7wNnP/YZq9LuIiACa0pZ82mcrxEVEpExqqYuIiISEQl1ERCQkFOoiIiIhoVAXEREJCYW6iIhISCjURUREQkKhLiIiEhIKdRERkZBQqIuIiISEQl1ERCQkFOoiIiIhoVAXEREJCYW6iIhISCjURUREQkKhLiIiEhIKdRERkZBQqIuIiISEQl1ERCQkFOoiIiIhoVAXEREJCYW6iIhISCjURUREQkKhLiIiEhIKdRERkZBQqIuIiISEQl1ERCQkFOoiIiIhoVAXEREJCYW6iIhISCjURUREQkKhLiIiEhIKdRERkZBQqIuIiIREzELdOfe8c26Zc256Ofc3dM6965yb5pyb4JzrFqtaREREUkEsW+rDgKMruP9WYKr3vjswAHg0hrWIiIiEXsxC3Xs/ClhVwSF7Av+LHDsDyHTONY9VPSIiImEX5Dn174BTAZxz+wHtgNZlHeicG+icm+Scm7R8+fI4ligiIpI8ggz1e4EGzrmpwFXAt0BhWQd674d677O891lNmzaNZ40iIiJJo3pQL+y9Xwf8CcA554A5wK9B1SMiIpLsAmupO+caOOdqRK5eDIyKBL2IiIjsgJi11J1zrwH9gSbOuYXA7UA6gPf+aWAP4EXnnAd+AC6KVS0iIiKpIGah7r0/p5L7xwKdY/X6IiIiqUYryomIiISEQl1ERCQkFOoiIiIhoVAXEREJCYW6iIhISCjURUREQkKhLiIiEhIKdRERkZBQqIuIiISEQl1ERCQkFOoiIiIhoVAXEREJCYW6iIhISCjURUREQkKhLiIiEhIKdRERkZBQqIuIiISEQl1ERCQkFOoiIiIhoVAXEREJCYW6iIhISCjURUREQkKhLiIiEhIKdRERkZBQqIuIiISEQl1ERCQkFOoiIiIhoVAXEREJCYW6iIhISCjURUREQkKhLiIiEhIKdRERkZBQqIuIiISEQl1ERCQkFOoiIiIhoVAXEREJCYW6iIhISCjURUREQkKhLiIiEhLVKzvAOdcXuANoFzneAd57v1tsSxMREZHtUWmoA88B1wKTgcLYliMiIiI7qiqhvtZ7/3HMKxEREZGdUpVQ/8o5dz/wDrC5+Ebv/ZSYVSUiIiLbrSqhvn/kZ1aJ2zxwaPTLERERkR1Vaah77w+JRyEiIiKycyqd0uacq++ce8g5NylyedA5Vz8exYmIiEjVVWWe+vPAeuDMyGUd8EIsixIREZHtV5Vz6h2896eVuH6nc25qrAoSERGRHVOVlvom51y/4iuRxWg2xa4kERER2RFVaakPAl6MnEd3wCrgwlgWJSIiItuvKqPfpwI9nHO7RK6vi3lVIiIist3KDXXn3Hne+1ecc9eVuh0A7/1DMa5NREREtkNFLfU6kZ/1yrjPx6AWERER2Qnlhrr3/pnIr194778peV9ksJyIiIgkkKqMfn+8ireJiIhIgCo6p34A0AdoWuq8+i5AWqwLExERke1T0Tn1GkDdyDElz6uvA06PZVEiIiKy/So6pz4SGOmcG+a9nxfHmkRERGQHVGXxmdzIfupdgVrFN3rvtfWqiIhIAqnKQLkcYAbQHrgTmAtMjGFNIiIisgOqEuqNvffPAfne+5He+z8DaqWLiIgkmKp0v+dHfv7mnDsOWAw0il1JIiIisiOqEup3RTZzuR6bn74LcG1MqxIREZHtVpUNXT6I/LoWOATAOVen/EeIiIhIECo8p+6ca+Wcy3LO1Yhcb+acuxuYFZfqREREpMrKDXXn3DXAVKzLfZxz7mLgJ6A20Cs+5YmIiEhVVdT9PhDY3Xu/yjnXFpgJ9PXeT45PaSIiIrI9Kup+z/PerwLw3s8Hflagi4iIJK6KWuqtnXOPlbi+a8nr3vurY1eWiIiIbK+KQv2GUtfVShcREUlgFW3o8mI8CxEREZGdU5VlYkVERCQJxCzUnXPPO+eWOeeml3N/fefc+86575xzPzjn/hSrWkRERFJBLFvqw4CjK7j/CuBH730PoD/wYPEiNyIiIrL9Klp85o0Sv99X6r7PKnti7/0oYFVFhwD1nHMOqBs5tqCy5xUREZGyVdRS71Ti9yNK3dc0Cq/9BLAHtuvb98BfvPdFZR3onBvonJvknJu0fPnyKLy0iIhI+FQU6n4H76uqo7BlaFsCewNPOOd2KfPFvB/qvc/y3mc1bRqN7xMiIiLhU9E89Qzn3D5Y8NeO/O4il9pReO0/Afd67z0w2zk3B+gCTIjCc4uIiKScikJ9CfBQGb8XX99Z84HDgK+dc82B3YFfo/C8IiIiKamixWf678wTO+dew0a1N3HOLQRuB9Ijz/008A9gmHPue6z1f5P3fsXOvKaIiEgqKzfUnXOdgPuBjthAtr967xdV9Ym99+dUcv9i4MiqPp+IiIhUrKKBcs8DHwKnAVOwfdVFREQkQVV0Tr2e9/7ZyO/3O+emxKMgERER2TEVhXqtEiPeYesR8HjvFfIiIiIJpKqj30tf98ChsSpKREREtl/MRr+LiIhIfFW09vuNJX4/o9R9d8eyKBEREdl+FY1+P7vE77eUuq+i3ddEREQkABWFuivn97Kui4iISMCquqFL6Q1corGhi4iIiERRRaHewzm3zjm3Huge+b34+l5xqk9ERCQp5ORAZiZUq2Y/c3LiX0NFo9/T4lmIiIhIssrJgYEDITfXrs+bZ9cBsrPjV0dFLXURERGpgsGDtwR6sdxcuz2eFOoiIiI7Yc0aa5mXZf78+NaiUBcREdkBU6bAJZdAq1blH9O2bfzqgYqXiRUREZESNm2CN96Ap56C8eOhdm07Z56ZCXffvXUXfEYGDBkS3/oU6iIiIpWYPRuefhpeeAFWrYLdd4dHH4UBA6BBAzsmM9POoc+fby30IUPiO0gOFOoiIiJlKiiADz+0Vvmnn0L16nDyyXD55dC/P7hSy7BlZ8c/xEtTqIuIiJSwZAn8+98wdCgsWGDnzO+8Ey6+GFq2DLq6iinURUQk5XkPo0bBk0/CO+9YK/2II6yL/YQTrJWeDJKkTBERkehbuxZeftm62H/80c6PX3UVXHYZdO4cdHXbT6EuIiIpZ+pUC/KcHNi4EbKy4Pnn4ayzbNR6slKoi4hISsjLg7fesi72sWOhVi045xwYNAj23Tfo6qJDoS4iIqH266/wzDPWEl+xAjp1goceggsugEaNgq4uuhTqIiISOoWF8PHH1ir/5BPbOe2kk6xVfuihdj2MFOoiIhIay5bBc89Zy3zePNh1V7jttsqXcw0LhbqIiCQ17+Gbb6xV/tZbkJ9vrfEHH4QTT4T09KArjB+FuoiIJKX16+GVVyzMp0+H+vVttbfLLoMuXYKuLhgKdRERSSrff2/T0V5+GTZsgH32sRXgzj4b6tQJurpgKdRFRCTh5ORsvTnKnXdCjRrWKh89GmrWtBAfNAj222/bddhTlUJdREQSSk4ODBy4ZRvTefPgwgvt9w4d4IEH7HrjxkFVmLgU6iIiElfew5o1sHRp2ZecHNu3vLRmzWDmzPBOR4sGhbqIiOy0oiJYubL8oC55WbbMRqiXVq0aNG1adqADLF+uQK+MQl1EJORKn58eMqRq+34XFFiQVhbQS5facYWF2z5Herq1sJs3t0v37lt+L3lp1sy609PSIDPTutxLa9t2pz+K0FOoi4iEWFnnpy+5xJZO7dmz/JBeutRa3t5v+5y1am0J47Ztbd30kuFcMqwbNtz+QWxDhmxdM9gmK0OG7PjnkCqcL+tPLIFlZWX5SZMmBV2GiEhSKK/VW1rdutu2nMtqUTdvDvXqxX60+Y72LoSVc26y9z6rsuPUUhcRCbH588u+3TkYM2ZLUCfadqPZ2akd4jtKoS4iEmKNG9vOZKW1bQu9e8e/HoktjSMUEQmpceNs6ljpEeM6Px1eCnURkRCaPRtOOAHatYN//ct+Omc/hw5V13ZYqftdRCRkVqyAY46xkesffwydOtkmJxJ+CnURkRDZtMm2G124EP73Pwt0SR0KdRGRkCgqgvPPt3Ppb74JBxwQdEUSbwp1EZGQuOEGePtteOghOO20oKuRIGignIhICDz2mIX51VfDNdcEXY0ERaEuIpLkhg+3ID/5ZAt27S2euhTqIiJJbPx4OPdc2G8/W1o1LS3oiiRICnURkST1yy82F33XXeG//028pV4l/hTqIiJJqHguemGhzUVv1izoiiQRaPS7iEiSycuz8+fz58OXX0LnzkFXJIlCoS4ikkSKimDAAPjmG3jjDejbN+iKJJGo+11EJIncdJMtLPPAA3DGGUFXI4lGoS4ikiSeeMLC/Mor4brrgq5GEpFCXUQkCfz3v/CXv9i67o88ornoUjaFuogEJicHMjNtv+/MTLsu25o4Ec4+G3r1gldf1Vx0KZ8GyolIIHJyYOBAyM216/Pm2XXQXt8l/forHH88tGgB778PdeoEXZEkMrXURSQQgwdvCfRiubl2u5iVK+HYY6GgwOaiN28edEWS6NRSF5G4W7nSWuZlmT8/vrUkquK56HPmwBdfwO67B12RJAO11EUkbnJz4Z57YLfdyj+mXj3YvDl+NSWioiK48EIYPRpeegkOPDDoiiRZKNRFJOYKCuDZZ6FTJ7j1VujfH+69d9u1ytPSYN062Hdf+PbbQEpNCLfcAv/5D/zzn3DWWUFXI8lEoS4iMeM9vPsudOtmg+AyM+Hrr+G992wRlaFDoV07m57Vrh28+CJ88AEsX267jt15J+TnB/0u4uvJJy3ML78c/vrXoKuRZOO890HXsF2ysrL8pEmTgi5DRCoxejTceCOMHQtduli3+0knVW1+9apVcNVVNn2rZ08L+27dYl9z0N5/386jH3usfRmqrlFPEuGcm+y9z6rsOLXURSSqfvjBFkg58EAbDPfss/D99xZWVV0wpVEjm/L29tuwYIHNz773XuvGD6tJk2wu+j77wOuvK9BlxyjURSQqFiyAP/8ZuneHUaOsZT5rFlx88Y4H1Kmn2peEE06w88z9+sHPP0e37kQwZw4cd5xtn/rBB5qLLjtOoS4iO2X1ajs/3rmzta6vvRZ++QVuvnnbgXA7omlT28DktdfsS8Lee8PDD9sI8TBYtcq62/Pz4aOPbJEZkR2lUBeRHbJpE9x/v01Pu/9+OPNMmDnTNhxp3Di6r+WcdU1Pnw5HHGGbmfTvb18ektnmzXDKKbZq3PDhsMceQVckyU6hLiLbpbAQXnjBWuY33gh9+sDUqTaYrV272L72rrvayPlhw2DaNOvqf/LJ5Gy1F89FHzXK3s9BBwVdkYSBQl1EqsR7G53do4edO2/ZEr76Cj780MI1XpyDCy6wVnu/fnDFFXDkkcm3Et3gwTYg7t574Zxzgq5GwkKhLiKVGjPGWpInnmjnft96C8aNsy7woLRuDZ98As88A+PH25S3556zLx+J7umnLcwvu8x6O0SiRaEuIuWaMcPO+fbtC7NnWxhNnw6nnZYY+3k7Z4vafP+9TXu7+GLb0Wzx4qArK9+HH1rvwnHHweOPJ8bnKOGhUBeRbSxaZGHZtSt8+SXcdZeF+qWXQnp60NVtKzPT6nzsMTsl0LUrvPJK4rXaJ0+2ZV/33ltz0SU2FOoi8oc1a2xt9k6dbPDWVVfZCPPBgxN/7nS1albvd9/BnnvC+efbPPelS4OuzMyda70ITZpYa71u3aArkjBSqIsIeXnw0EPQoYMtGnPqqbbIyyOP2DzxZNKpk40ov/9+24O8a1eb5x6k1attLnpenuaiS2zFLNSdc88755Y556aXc/8Nzrmpkct051yhc65RrOoRkW0VFtrWnrvvDtdfb7ujTZliXdft2wdd3Y5LS7PNUKZMsXn0Z55p89xXrIh/LcVz0WfPtrnoe+4Z/xokdcSypT4MOLq8O73393vv9/be7w3cAoz03q+KYT0iEuG9tWJ79rTpYU2bwhdf2GjyffYJurro2XNPG7l/113wzjs2Qv699+L3+t7b9L+RI+10xsEHx++1JTXFLNS996OAqob0OcBrsapFRLaYMAEOPdS6gzdutAFbEybAYYcFXVlsVK9uYwImTrRu75NPhgEDrEs81v72N9tp7u674dxzY/96IoGfU3fOZWAt+rcrOGagc26Sc27S8uXL41ecSIjMnAlnnAH77w8//ghPPGE/zzrLBpmFXY8e9uXl73+3oO3WzXomYmXoUAvzSy6xdfBF4iER/imfAHxTUde7936o9z7Le5/VNNlG7YgE7LffYNAg64r++GO44w47v3vFFVCjRtDVxVeNGvB//2cL5zRoAMccY1P31q+P7ut8/DFcfjkcfbQtY6u56BIviRDqZ6Oud5GdlpNj87WrVbOf//433HYbdOxovw8aZNPTbr8d6tULutpgZWXZnPGbbrJV6PbaC/73v+g895Qp1iPSvTu88Ybmokt8OR/D1Rmcc5nAB977buXcXx+YA7Tx3m+synNmZWX5SZMmRa1GkTDIybEWZ27utveddZYNFOvYMf51JYOxY21jlZkzrffivvt2fE7+vHnQu7f1CIwda+vji0SDc26y9z6rsuNiOaXtNWAssLtzbqFz7iLn3GXOuctKHHYK8FlVA11EyjZ4cNmB3qKFDYRToJfvgAPg22/hmmusq7xHDxg9evufZ80aG3y4aZPNRVegSxBi2lKPBbXURbZVrVrZS6I6l5zbkgZl1Chrtc+dC9deaz0ctWtX/rjff7fz56NH2+C7Qw+NdaWSagJvqYtIfCxcWP567G3bxreWZHfQQbZP+6WX2gp7PXvaDnAV8R4uusjWnH/+eQW6BEuhLpLERoyw3cmc23Yke0YGDBkSSFlJrW5deOop+Owzm8ffp4+th795c9nH33abrcB3111w3nnxrVWkNIW6SBLyHh54AA4/HBo1snPCzz8P7dpZwLdrZ/Oks7ODrjR5HXGEbel64YW2Hv6++9rnXHKWQePGFuYXX2zBLxI0TbYQSTLr19vSo2+9Zfuav/CCTVHbYw+FeLTVr29T3k491RaR6dXLpqjl59v9q1ZZuB94oOaiS2JQS10kicyYAfvtZ+uY33+/7T6W6nPO4+G442D6dBs0VxzoxYqKrAteJBGopS6SJN5+27qCa9e2zVcOOSToilJLo0Y2Xa0s8+fHtxaR8qilLpLgCgrgxhvh9NNtb/ApUxToQSlvNoFmGUiiUKiLJLBly+Coo6yrfdAg28Kzdeugq0pdQ4bYrIKSNMtAEolCXSRBjR9vA7PGjLG9uJ98EmrWDLqq1JadbbMKNMtAEpXOqYskGO8tKK6+2pYaHTMG9tkn6KqkWHa2QlwSl1rqIglk0yabrnbZZbYy2eTJCnQRqTqFukiCmDMH+va1rvbbboMPPrAR1yIiVaXud5EE8MkncO651vX+/vtw/PFBVyQiyUgtdZEAFRXBP/5hW3a2aQOTJinQRWTHqaUuEpA1a+D8862b/bzz4Jlntp0uJSKyPRTqIgGYNs3WE583D554Ai6/XGuHi8jOU/e7SJy98gr07m0j3UeOhCuuUKCLSHQo1EXi5Pff4aqrrMt9331tulqfPkFXJSJholAXiYPFi2299ieegOuusw1ZWrQIuioRCRudUxeJsVGj4MwzYcMGeP11OOusoCsSkbBSS10kRryHhx+2leHq17e13BXoIhJLCnWRGNiwAc45x7raTzgBJkywbVNFRGJJoS4SZTNn2uj2N9+Ee++Fd96xlrqISKzpnLpIFA0fDgMG2Bapn34Khx8edEUikkrUUheJgsJCuOUWOOUU6NLFpqsp0EUk3tRSF9lJy5fbZixffAGXXAKPPQa1agVdlYikIoW6yE6YOBFOOw2WLYN//xsuuijoikQklan7XWQHPfss9OtnS7yOHq1AF5HgKdRFtlNeHlx8MQwcCAcfbOfPs7KCrkpERKEusl3mzrXW+XPPweDB8PHH0KRJ0FWJiBidUxepos8+swVlCgrgvffgxBODrkhEZGtqqYtUoqgI7r4bjj4aWraESZMU6CKSmBTqImXIyYHMTKhWDerWta72s8+GceOgU6egqxMRKZu630VKycmxQXC5uXZ90yZIT4fjjoM6dYKtTUSkImqpi5QyePCWQC+Wn2+3i4gkMoW6SCnz5pV9+/z58a1DRGR7qftdJMJ7GxBXnrZt41eLiMiOUEtdBJumNmgQ/O1v0LcvZGRsfX9GBgwZEkxtIiJVpVCXlJeba+u3P/MM3HQTjBoFQ4dCu3a2BGy7dnY9OzvoSkVEKqbud0lpK1bACSfA+PHw+ONw5ZV2e3a2QlxEko9CXVLWnDm2oMy8efDWW3DqqUFXJCKycxTqkpKmTIFjj4Xff7d90Pv1C7oiEZGdp3PqknI++8x2V6tZE775RoEuIuGhUJeU8tJLtjLcbrvB2LGwxx5BVyQiEj0KdUkJxXPQL7gADjrIRri3bBl0VSIi0aVQl9ArLLRR7YMHw7nn2h7o9esHXZWISPQp1CXUNm2C00+HJ5+EG26Al1+GGjWCrkpEJDY0+l1Ca+VK2/d87Fh49FG4+uqgKxIRiS2FuoTS3Lk2B33uXHjjDWuti4iEnUJdQmfqVDjmGMjLs+lrBx0UdEUiIvGhc+oSKp9/biGeng6jRyvQRSS1KNQlNF55xVaJy8y08+hduwZdkYhIfCnUJel5D/fdB+efDwceCF9/Da1aBV2ViEj8KdQlqRUW2qj2m2+Gs8/WHHQRSW0KdUlamzbBmWfCE0/AX/8KOTm2nruISKrS6HdJSqtWwUkn2YYsDz8M11wTdEUiIsFTqEvSmTfPpqz98gu8/rq11kVERKEuSea77yzQc3O3bKEqIiJG59QlaXz5pY1uT0uzOegKdBGRrSnUJSm8+qq10Nu1szno3boFXZGISOJRqEtC8x7uvx+ys6FPH5uD3rp10FWJiCQmhbokrMJCG9V+4402GO7TT6FBg6CrEhFJXAp1SUh5ebaYzGOPwbXXwmuvaQ66iEhlNPpdEs7q1TYH/euv4cEH4brrgq5IRCQ5KNQlocyfbwPiZs+21vnZZwddkYhI8lCoS8KYNs0CfcMG+OQTOOSQoCsSEUkuOqcuCeGrr2wOunM2B12BLiKy/RTqErjXX4ejjoI2bWwO+l57BV2RiEhyStlQz8mBzEyoVs1+5uQEXVFqeughOOccOOAAGxjXpk3QFYmIJK+UPKeekwMDB9r64WAbhAwcaL9nZwdXVyopKoLrr4dHHoEzzoCXXoJatYKuSkQkuaVkS33w4C2BXiw3126X2MvLs9b5I4/AX/5i3e8KdBGRnZeSLfX587fvdomeNWvg5JNh5Eh44AGbg+5c0FWJiIRDzFrqzrnnnXPLnHPTKzimv3NuqnPuB+fcyFjVUlrbtmXfXqMG/PxzvKpIHSXHLzyZRToAAAsoSURBVDRrZufOX33Vut8V6CIi0RPL7vdhwNHl3emcawA8CZzove8KnBHDWrYyZAhkZGx9W40aFjrdu8Pf/w6bNsWrmnArHr8wb55tzpKfD+npdk5dRESiK2ah7r0fBayq4JBzgXe89/Mjxy+LVS2lZWfD0KG2jadz9vP552HOHNs45K67bGvPTz6JV0Xhdcst245f2LxZ4xdERGIhyIFynYGGzrkRzrnJzrkB5R3onBvonJvknJu0fPnyqLx4djbMnWstxrlz7Xrz5vDyy/Dll9aaPOYYC/lFi6LykimloMC+OC1YUPb9Gr8gIhJ9QYZ6daAXcBxwFPB351znsg703g/13md577OaNm0a88IOPRS++85a7O+/D126wKOPWlBJxbyHd96Brl3h0kvL31mtvHENIiKy44IM9YXAp977jd77FcAooEeA9WylZk3rIp4+Hfr1s329990Xxo8PurLENXKkLSJz2mlQvTq89x4899y24xcyMmxcg4iIRFeQof4e0M85V905lwHsD/wUYD1l6tABPvoI3nwTli2z0Bo0yLYHFTNtGhx3HPTvb6cqnnvOejpOPLHs8QtDh2qRHxGRWHDe+9g8sXOvAf2BJsBS4HYgHcB7/3TkmBuAPwFFwL+9949U9rxZWVl+0qRJMam5MuvWwe23w2OPQZMmttd3dnbqTsuaNw9uu83GIdSvD7feCldeCbVrB12ZiEi4OOcme++zKj0uVqEeK0GGerFvv4XLLoMJE2w3sSeftPPuqWLlSrj7bnjiCftC85e/wM03Q8OGQVcmIhJOVQ31lFwmdmftsw+MGQNPPWUBnypz2zdutDDfbTdb4vW882DWLLjvPgW6iEgiUKjvoLQ0a63PmAFnnRXuue3F09M6dbLBg4ccAt9/b+fOtauaiEjiUKjvpLLmtp9xRjjmthdPT+vWzaantW8Po0fD8OGw555BVyciIqUp1KOk5Nz2Dz5I/rntJaenpaXZ9LTRo6Fv36ArExGR8ijUoygMc9tLTk9buHDr6WmpOspfRCRZKNRjIBnnts+bBxdcAHvvbYMA77vPBsH9+c+2kIyIiCQ+hXqMOAennw4//WRTvoYOtS75V16xc9WJYuVK2wK1c2f4z3/ghhvg11/hxhs131xEJNko1GNsl13g4Ydh0iQbaHb++XDYYTZqPki5uXDPPZqeJiISJgr1OCme2/7008HObS8ogGefhY4dbQW4/v3tPLqmp4mIJD+FehxVq2ZTw2bMgLPPju/cdu/h3Xft9QYOtF6Dr7+2Ue1du8b+9UVEJPYU6gFo3hxeegn+97/4zG0fNQr69IFTT7UvFsOH2/S0fv1i83oiIhIMhXqADjkktnPbv/+e/2/v7mLtqMowjv8fWxJbICgRGgVprWkgfFiQowGRj4AaI4QaQsSkEkgkxsQPMBCC8UJvMAhGPjTRGKwSKXBRMRpRlCAiopG0CC0FtAGlrVKLGK0YUgp9vZg5WmsFeyx7zl7n/0tO9px19p55z7rYz57Zs9bijDPg5JNhw4buEvvq1bBkicPTJKlFhvrAJse2r10LJ564Z8a2r18P558PixfDvfc6PE2SZgpDfZpYuBBuuw1WrJj62Pann4ZLLumGp91yS7f92GMOT5OkmcJQn0aSblrWRx/dvbHtOw5Pu/rqbo33devgyith//1HV78kaViG+jS0777/fWz78uWwYEF3w9v8+XDBBd3qaQ5PkySlptP0Zv+DiYmJWrly5dBljMz27d248ssugy1bujDf+Ua6RYtg2TLvZpekViVZVVUTL/U8z9SnuR3Hts+Zs+s747duNdAlSYb62Jg3r/vufFc2bBhtLZKk6clQHyOHHLJ77ZKkmcVQHyOXXw5z5/5729y5XbskSYb6GFm6tBvmNn9+N/xt/vzu96VLh65MkjQdOL/YmFm61BCXJO2aZ+qSJDXCUJckqRGGuiRJjTDUJUlqhKEuSVIjDHVJkhphqEuS1AhDXZKkRhjqkiQ1wlCXJKkRhrokSY0w1CVJaoShLklSIwx1SZIaYahLktSIVNXQNeyWJE8BTwxdxzTwGuBPQxcxA9jPo2E/j4b9PBovRz/Pr6oDXupJYxfq6iRZWVUTQ9fROvt5NOzn0bCfR2PIfvbyuyRJjTDUJUlqhKE+vr46dAEzhP08GvbzaNjPozFYP/uduiRJjfBMXZKkRhjqYyTJ65PcleThJGuTXDh0TS1LMivJr5J8b+haWpbkVUlWJHk0ySNJjh+6phYl+UT/vvFQkpuTvHLomlqQZFmSzUke2qFt/yR3JFnXP756VPUY6uPleeDiqjocOA74SJLDB66pZRcCjwxdxAxwLXB7VR0GLMY+3+OSHAR8HJioqiOBWcD7h62qGd8A3r1T22XAnVW1CLiz/30kDPUxUlVPVtX9/fbf6N78Dhq2qjYlORg4Hbh+6FpalmQ/4CTgawBV9VxV/WXYqpo1G5iTZDYwF/jDwPU0oap+Cvx5p+YlwA399g3Ae0dVj6E+ppIsAI4BfjlsJc26BrgU2D50IY17A/AU8PX+q47rk+w9dFGtqarfA58H1gNPAn+tqh8NW1XT5lXVk/32JmDeqA5sqI+hJPsA3wIuqqotQ9fTmiRnAJuratXQtcwAs4E3A1+uqmOAvzPCS5UzRf+d7hK6D1GvA/ZO8oFhq5oZqhtiNrJhZob6mEmyF12gL6+qW4eup1EnAGcm+R1wC3BqkhuHLalZG4GNVTV5xWkFXchrz3oH8NuqeqqqtgG3Am8buKaW/THJawH6x82jOrChPkaShO67x0eq6gtD19OqqvpkVR1cVQvobib6cVV5VvMyqKpNwIYkh/ZNpwEPD1hSq9YDxyWZ27+PnIY3JL6cvguc12+fB3xnVAc21MfLCcC5dGeOD/Q/7xm6KOn/9DFgeZLVwNHAZweupzn9lZAVwP3AGrr3fmeX2wOS3Az8Ajg0ycYkHwSuAN6ZZB3dVZIrRlaPM8pJktQGz9QlSWqEoS5JUiMMdUmSGmGoS5LUCENdkqRGGOpSA5K80A9xXJvkwSQXJ3lF/7eJJNcNVNfPhziuNFM5pE1qQJJnqmqffvtA4Cbg3qr69LCVSRolz9SlxlTVZuBDwEfTOWVyTfgkn0lyQ5J7kjyR5KwkVyZZk+T2fhpikhyb5O4kq5L8cIcpL3+S5HNJ7kvymyQn9u1H9G0PJFmdZFHf/kz/mCRX9Wt5r0lyTt9+Sr/PyfXUl/cznkmaAkNdalBVPU63ZvaBu/jzG4FTgTOBG4G7quoo4Fng9D7YvwicXVXHAsuAy3d4/eyqeitwETB5JeDDwLVVdTQwQTen+47OopstbjHdDFtXTX5QoFtt8CLgcGAh3cyJkqZg9tAFSBq5H1TVtiRr6IL/9r59DbAAOBQ4ErijP2meRbdc56TJhYRW9c+HbprMT/Xr0N9aVet2OubbgZur6gW6xS7uBt4CbAHuq6qNAEke6Pf5sz3yn0ozjGfqUoOSLAReYNerQ20FqKrtwLb614012+k+6AdYW1VH9z9HVdW7dn59v//Z/b5uojvzfxb4fpJTd6PcrTts/3OfknafoS41JskBwFeAL9XU7oT9NXBAkuP7/e2V5IiXOOZC4PGquo5uRao37fSUe4Bzkszq6zsJuG8KtUl6EX4iltowp790vRfwPPBNYErL81bVc0nOBq5Lsh/d+8Q1wNoXedn7gHOTbAM28Z8rrX0bOB54ECjg0qralOSwqdQoadcc0iZJUiO8/C5JUiMMdUmSGmGoS5LUCENdkqRGGOqSJDXCUJckqRGGuiRJjTDUJUlqxD8ApRo3/DvvH/YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "simulations = 10000\n",
    "sample_size = 500\n",
    "\n",
    "def least_square_error(x_0, y_0, train_x, train_y):\n",
    "    X = np.c_[np.ones((len(train_x), 1)), train_x]\n",
    "    beta = np.linalg.lstsq(X, train_y, rcond = None)[0]\n",
    "    return (np.dot(np.array([1, *x_0]), beta) - y_0) ** 2\n",
    "\n",
    "# 1-nearest neighbor error\n",
    "def nn_error(x_0, y_0, train_x, train_y):\n",
    "    X = (train_x * train_x).sum(axis=1)\n",
    "    return (y_0 - train_y[X.argmin()]) ** 2\n",
    "\n",
    "cubic_epe_ratio = []\n",
    "linear_epe_ratio = []\n",
    "\n",
    "for p in range(1, 11):\n",
    "    least_square_epe = [0, 0]\n",
    "    nn_epe = [0, 0]\n",
    "    for _ in range(simulations):\n",
    "        error = np.random.standard_normal(sample_size)\n",
    "        train_x = np.random.uniform(-1, 1, size=(sample_size, p))\n",
    "        train_y = [train_x[:, 0] + error, \n",
    "                   0.5 * (train_x[:, 0] + 1) ** 3 + error]\n",
    "        x_0 = np.zeros(p)\n",
    "        y_0 = [np.random.standard_normal(),\n",
    "               0.5 + np.random.standard_normal()]\n",
    "        for i in range(2):\n",
    "            least_square_epe[i] += least_square_error(x_0, y_0[i], \n",
    "                                                      train_x, train_y[i])\n",
    "            nn_epe[i] += nn_error(x_0, y_0[i], \n",
    "                                  train_x, train_y[i])\n",
    "    for i in range(2):\n",
    "        least_square_epe[i] /= simulations\n",
    "        nn_epe[i] /= simulations\n",
    "    linear_epe_ratio.append(nn_epe[0] / least_square_epe[0])\n",
    "    cubic_epe_ratio.append(nn_epe[1] / least_square_epe[1])\n",
    "\n",
    "# plot\n",
    "fig = plt.figure(1, figsize = (8, 8))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.set_title(\"Expected Prediction Error of 1NN vs. OLS\")\n",
    "\n",
    "axes.plot(np.arange(1, 11), linear_epe_ratio, '-o',\n",
    "          color = 'orange', label = \"Linear\")\n",
    "\n",
    "axes.plot(np.arange(1, 11), cubic_epe_ratio, '-o',\n",
    "          color = 'blue', label = \"Cubic\")\n",
    "\n",
    "axes.legend()\n",
    "axes.set_xlabel(\"Dimension\")\n",
    "axes.set_ylabel(\"EPE Ration\")\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}