{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Linear Regression" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Preliminaries\n", "\n", "- Goal \n", " - Maximum likelihood estimates for various linear regression variants\n", "- Materials \n", " - Mandatory\n", " - These lecture notes\n", " - Optional\n", " - Bishop pp. 140-144 \n", " - [G. Deng et al., _A model-based approach for the development of LMS algorithms_', ISCAS-05 symposium, 2005](./files/Deng-2005-A-model-based-approach-for-the-development-of-LMS-algorithms.pdf)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Regression - Illustration\n", "\n", "\n", "\n", "\n", "Given a set of (noisy) data measurements, find the 'best' relation between an input variable $x \\in \\mathbb{R}^D$ and input-dependent outcomes $y \\in \\mathbb{R}$\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Regression vs Density Estimation\n", "\n", "\n", "- Observe $N$ IID data **pairs** $D=\\{(x_1,y_1),\\dotsc,(x_N,y_N)\\}$ with $x_n \\in \\mathbb{R}^D$ and $y_n \\in \\mathbb{R}$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- [Q.] We could try to build a model for the data by density estimation, $p(x,y)$, but what if we are interested only in (a model for) the responses $y_n$ for **given inputs** $x_n$?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- [A.] We will build a model only for the conditional distribution $p(y|x)$. \n", " - Note that, since $p(x,y)=p(y|x)\\, p(x)$, this is a building block for the joint data density.\n", " - In a sense, this is density modeling with the assumption that $x$ is drawn from a uniform distribution." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Next, we discuss model (1) specification, (2) ML estimation and (3) prediction for the linear regression model. \n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### 1. Model Specification for Linear Regression\n", "\n", "\n", "- In a _regression_ model, we try to 'explain the data' by a purely deterministic term $f(x,w)$, plus a purely random term $\\epsilon_n$ for 'unexplained noise',\n", "\n", " $$\n", " y_n = f(x_n,w) + \\epsilon_n\n", " $$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In _linear regression_, we assume that \n", "\n", "$$f(x,w)=w^T x \\,.$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In _ordinary linear regression_, the noise process $\\epsilon_n$ is zero-mean Gaussian with constant variance $\\sigma^2$, i.e.\n", "\n", "$$\n", "y_n = w^T x_n + \\mathcal{N}(0,\\sigma^2) \\,,\n", "$$\n", "or equivalently, the likelihood model is \n", "$$\n", "p(y_n|\\,x_n,w) = \\mathcal{N}(y_n|\\,w^T x_n,\\sigma^2) \\,.\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- For full Bayesian learning we should also choose a prior $p(w)$; In ML estimation, the prior $p(w)$ is uniformly distributed (so it can be ignored).\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### 2. ML Estimation for Linear Regression Model\n", "\n", "- Let's work out the log-likelihood for multiple observations\n", "$$\\begin{align*}\n", "\\log p(D|w) &\\stackrel{\\text{IID}}{=} \\sum_n \\log \\mathcal{N}(y_n|\\,w^T x_n,\\sigma^2) \\propto -\\frac{1}{2\\sigma^2} \\sum_{n} {(y_n - w^T x_n)^2}\\\\\n", " &= -\\frac{1}{2\\sigma^2}\\left( {y - \\mathbf{X}w } \\right)^T \\left( {y - \\mathbf{X} w } \\right)\n", "\\end{align*}$$\n", "where we defined $N\\times 1$ vector $y = \\left(y_1 ,y_2 , \\ldots ,y_N \\right)^T$ and $(N\\times D)$-dim matrix $\\mathbf{X} = \\left( x_1 ,x_2 , \\ldots ,x_n \\right)^T$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Set the derivative $\\nabla_{w} \\log p(D|w) = \\frac{1}{\\sigma^2} \\mathbf{X}^T(y-\\mathbf{X} w)$ to zero for\n", "the maximum likelihood estimate\n", "$$\\begin{equation*}\n", "\\boxed{\\hat w_{\\text{ML}} = (\\mathbf{X}^T \\mathbf{X})^{-1} \\mathbf{X}^T y}\n", "\\end{equation*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The matrix $\\mathbf{X}^\\dagger \\equiv (\\mathbf{X}^T \\mathbf{X})^{-1}\\mathbf{X}^T$ is also known as the **Moore-Penrose pseudo-inverse** (which is sort-of-an-inverse for non-square matrices)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that size ($N\\times D$) of the data matrix $\\mathbf{X}$ grows with number of observations, but the size ($D\\times D$) of $\\mathbf{X}^T\\mathbf{X}$ is independent of training data set." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### 3. Prediction of New Data Points \n", "\n", "- Now, we want to apply the trained model. New data points can be predicted by\n", "$$\\begin{equation*}\n", "p(y_\\bullet \\,|\\, x_\\bullet,\\hat w_{\\text{ML}}) = \\mathcal{N}(y_\\bullet \\,|\\, \\hat w_{\\text{ML}}^T x_\\bullet, \\sigma^2 ) \n", "\\end{equation*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that the expected value of a predicted new data point\n", "\n", "$$\n", "\\mathrm{E}[y_\\bullet] = \\hat w_{\\text{ML}}^T x_\\bullet = x_\\bullet^T \\hat{w}_{\\text{ML}} = \\left( x_\\bullet^T \\mathbf{X}^\\dagger \\right) y\n", "$$\n", "\n", "can also be expressed as a linear combination of the observed data points \n", "\n", "$$y = \\left( {y_1 ,y_1 , \\ldots ,y_N } \\right)^T \\,.$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Deterministic Least-Squares Regression\n", "\n", "- (You may say that) we don't need to work with probabilistic models. E.g., there's also the deterministic **least-squares** solution: minimize sum of squared errors,\n", "$$\\begin{align*} \\hat w_{\\text{LS}} &= \\arg\\min_{w} \\sum_n {\\left( {y_n - w ^T x_n } \\right)} ^2 \n", " = \\arg\\min_{w} \\left( {y - \\mathbf{X}w } \\right)^T \\left( {y - \\mathbf{X} w } \\right)\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Setting the gradient \n", "$ \\frac{\\partial \\left( {y - \\mathbf{X}w } \\right)^T \\left( {y - \\mathbf{X}w } \\right)}{\\partial w} = -2 \\mathbf{X}^T \\left(y - \\mathbf{X} w \\right)\n", "$ to zero yields the **normal equations** \n", "$\\mathbf{X}^T\\mathbf{X} \\hat w_{\\text{LS}} = \\mathbf{X}^T y$ and consequently\n", "$$\n", "\\boxed{\\hat w_{\\text{LS}} = (\\mathbf{X}^T \\mathbf{X})^{-1} \\mathbf{X}^T y} \n", "$$\n", "which is the same answer as we got for the maximum likelihood weights $\\hat w_{\\text{ML}}$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- $\\Rightarrow$ Least-squares regression ($\\hat w_{\\text{LS}}$) corresponds to (probabilistic) maximum likelihood ($\\hat w_{\\text{ML}}$) if \n", " 1. **IID samples** (determines how errors are combined), and\n", " 1. Noise $\\epsilon_n \\sim \\mathcal{N}(0,\\,\\sigma^2)$ is **Gaussian** (determines error metric) \n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Probabilistic vs. Deterministic Approach\n", "\n", "- The (deterministic) least-squares approach assumed IID Gaussian distributed data, but these assumptions are not obvious from looking at the least-squares (LS) criterion." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- If the data were better modeled by non-Gaussian assumptions (or not IID), then LS might not be appropriate." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The probabilistic approach makes all these issues completely transparent by focusing on the **model specification** rather than the error criterion." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Next, we will show this by two examples: (1) samples not identically distributed, and (2) few data points." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Not Identically Distributed Data\n", "\n", "- What if we assume that the variance of the measurement error varies with the sampling index, $\\epsilon_n \\sim \\mathcal{N}(0,\\sigma_n^2)$?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Let's make the log-likelihood again (use $\\Lambda \\triangleq \\mathrm{diag}[1/\\sigma_n^2]$): \n", "$$\\begin{align*}\n", "\\mathrm{L(w)} &\\triangleq \\log p(D|w) \n", " \\propto -\\frac{1}{2} \\sum_n \\frac{(y_n-w^T x_n)^2}{\\sigma_n^2} = -\\frac{1}{2} (y- \\mathbf{X}w)^T \\Lambda (y- \\mathbf{X} w)\\,.\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Set derivative $\\partial \\mathrm{L(w)} / \\partial w = -\\mathbf{X}^T\\Lambda (y-\\mathbf{X} w)$\n", "to zero to get the **normal equations** \n", "$\\mathbf{X}^T \\Lambda \\mathbf{X} \\hat{w}_{\\text{WLS}} = \\mathbf{X}^T \\Lambda y$ \n", "and consequently \n", "$$ \\boxed{\\hat{w}_{\\text{WLS}} = \\left(\\mathbf{X}^T \\Lambda \\mathbf{X}\\right)^{-1} \\mathbf{X}^T \\Lambda y}$$ " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- This is also called the **Weighted Least Squares** (WLS) solution. (Note that we just stumbled upon it, the crucial aspect is appropriate model specification!)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note also that the dimension of $\\Lambda$ grows with the number of data points. In general, models for which the number of parameters grow as the number of observations increase are called **non-parametric models**." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### CODE EXAMPLE\n", "\n", "We'll compare the Least Squares and Weighted Least Squares solutions for a simple linear regression model with input-dependent noise:\n", "\n", "$$\\begin{align*}\n", "x &\\sim \\text{Unif}[0,1]\\\\\n", "y|x &\\sim \\mathcal{N}(f(x), v(x))\\\\\n", "f(x) &= 5x - 2\\\\\n", "v(x) &= 10e^{2x^2}-9.5\\\\\n", "\\mathcal{D} &= \\{(x_1,y_1),\\ldots,(x_N,y_N)\\}\n", "\\end{align*}$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "scrolled": false, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAG2CAYAAACd5Zf9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XlclOXaB/DfsMywiLgiqCCIC5q7pGkuue8KHTVNTT1p9aaVlalpooHGqWNleSrfXCs1W8EttU6Kkembmpq5gAsKua8gIgzMXO8fjwyMDDgss/L7fj7z0bnneZ65BsHn4r6v+75VIiIgIiIicjIutg6AiIiIyBKY5BAREZFTYpJDRERETolJDhERETklJjlERETklJjkEBERkVNikkNEREROiUkOEREROSUmOUREROSUmOQQERGRU3KYJCcvLw9vvPEGQkJC4OnpiYYNGyI6Ohp6vd7WoREREZEdcrN1AOZ6++23sXTpUnz22Wd46KGHsH//fkycOBG+vr546aWXbB0eERER2RmHSXL27NmDYcOGYdCgQQCA4OBgfPnll9i/f7+NIyMiIiJ75DBJTpcuXbB06VIkJyejSZMmOHz4MH799VcsXry42HNycnKQk5NjeK7X63Hjxg3UrFkTKpXKGmETERFROYkIbt++jbp168LFpRSVNuIg9Hq9zJo1S1Qqlbi5uYlKpZK33nqrxHPmzZsnAPjggw8++OCDDyd4pKWllSp3UImIwAGsX78er732Gv7973/joYcewqFDhzBt2jS89957GD9+vMlz7u/JSU9PR1BQENLS0lC1alVrhU5ERHbs7bffxltvvYXZs2dj5syZRZ6T7WVkZCAwMBC3bt2Cr6+v2ec5TJITGBiIWbNmYcqUKYa2BQsWYM2aNThx4oRZ18jIyICvry/S09OZ5BARkUFMTAyioqKgVquh1WoRHR2NuXPn2josuqes92+HmUKelZVVZBzO1dWVU8iJiKjc5s6da0hw1Go1Exwn4TBJzpAhQ7Bw4UJs2bIFZ8+eRVxcHN577z1ERkbaOjQiInJwMTExhgRHq9UiJibG1iFRBXCYJGfJkiUYPnw4nn/+eTRr1gzTp0/Hs88+y29EIiIql/yhqujoaOTk5CA6OhpRUVG8vzgBh5lC7uPjg8WLF5c4Zbwi6PV6aLVai74H0f3c3d3h6upq6zCIKp3CCU7+EFX+n1FRUUbPyfE4TJJjDVqtFikpKazzIZuoVq0a/P39uYYTkRXpdDqTRcb5z3U6nS3CogriMLOrKkJJ1dkigtTUVOTm5pZ+sSGichARZGVl4cqVK6hWrRoCAgJsHRIRkV0p6+wq9uTck5eXh6ysLNStWxdeXl62DocqGU9PTwDAlStX4Ofnx6ErIqIKwO6Ke/K7JNVqtY0jocoqP7nOzc21cSRERM6BSc59WA9BtsLvPSKiisUkh4iIiJwSkxwHJyJ45plnUKNGDahUKhw6dAjXr1+Hn58fzp49a9Y1cnJyEBQUhAMHDlg2WCIiIitikuPgtm3bhtWrV2Pz5s24ePEiWrRogdjYWAwZMgTBwcFmXUOj0WD69OnciI6IiJwKkxwHd/r0aQQEBKBz587w9/dHbm4uVqxYgUmTJpXqOmPGjEFiYiKOHz9uoUiJiIisi0mOA5swYQJeeOEFpKamQqVSITg4GFu3boWbmxs6depkOC46Ohp169bF9evXDW1Dhw5Ft27dDAsf1qxZE507d8aXX35p9c9BRERkCVwn5wHu3LlT7Guurq7w8PAw61gXFxfDWiglHevt7W12bB988AFCQ0Px6aefYt++fXB1dcXChQsRHh5udNycOXOwbds2TJo0CXFxcVi6dCl++eUXHD582GjRww4dOiAxMdHs9yciIrJnTHIeoEqVKsW+NnDgQGzZssXw3M/PD1lZWSaP7d69OxISEgzPg4ODce3atSLHlWYBal9fX/j4+MDV1RX+/v4AgLNnz6Ju3bpGx7m6umLNmjVo06YNZs2ahSVLluDTTz9FgwYNjI6rV6+e2cXKRM5s/vz5cHV1NblnUUxMDHQ6HebPn2/9wIioVDhc5WTu3r1r1LuUr2HDhli0aBHefvttDBkyBGPGjClyjKenZ7FJGlFl4urqanIX6vzNHLkiNZFjYE/OA2RmZhb72v3/0V25cqXYY+/fC8tSPSa1atXCzZs3Tb72yy+/wNXVFWfPnkVeXh7c3Iz/+W/cuIHatWtbJC4iR2JqF2pTu1UTkX1jkvMApamRsdSxpdG2bVusWbOmSPtXX32F77//HgkJCXjiiScQExODN9980+iYv/76C23btrVIXESOpnCis2DBAmi1WiY4RA6Gw1VOpl+/fjh69KhRb87ff/+N//mf/8Hbb7+NLl26YPXq1YiNjcXevXuNzk1MTETfvn2tHTKR3Zo7dy7UajW0Wi3UajUTHCIHwyTHybRs2RLh4eH4+uuvASiFzBMmTECHDh0wdepUAECfPn0wdepUjB071jAct2fPHqSnp2P48OE2i53I3sTExBgSHK1WW6RGh4jsG5McBzdt2rQi9T1z587FBx98AL1eD5VKhf/+97/Ytm2b0QaQ7733Hk6dOmWYPfbee+/htddeM5rmTlSZFa7BycnJQXR0tMliZCKyX6zJcUIDBw7EyZMncf78eQQGBj7w+JycHLRu3Rovv/yyFaIjsn+mioxNFSMTkX1jkuOkXnrpJbOP1Wg0eOONNywYDZFj0el0JouM85/rdDpbhEVEpcQkh4joPiUt9MceHCLHwZocIiIickpMcoiIiMgpMckhIiIip8Qkh4iIiJwSkxwiIiJySkxyiIiIyCkxySEiIiKnxCTHwU2YMAEqlQoqlQru7u6oU6cO+vTpg5UrV0Kv19s6PCIiIpthklNB5s+fX+yeNjExMSUuLlZe/fv3x8WLF3H27Fls3boVPXr0wEsvvYTBgwcjLy/PYu9LRERkz5jkVBBXV1eTm/fl74Hj6upqsffWaDTw9/dHvXr10K5dO8yePRsbNmzA1q1bsXr1aou9LxERkT3jtg4VxNTmfaY2+bOWnj17onXr1vj+++8xadIkq743ERGRPWCSU4EKJzoLFiyAVqu1SYKTLywsDH/++adN3puIiMjWOFxVwebOnQu1Wg2tVgu1Wm3TzfxEBCqVymbvT0REZEtMcipYTEyMIcHRarXFFiNbw/HjxxESEmKz9yciIrIlJjkVqHANTk5ODqKjo00WI1vDjh07cOTIEfzjH/+w+nsTERHZA9bkVBBTRcamipEtIScnB5cuXYJOp8Ply5exbds2xMbGYvDgwXjqqacs8p5ERET2jklOBdHpdCaLjPOf63Q6i733tm3bEBAQADc3N1SvXh2tW7fGhx9+iPHjx8PFhZ11RERUOalERGwdhLVkZGTA19cX6enpqFq1qtFr2dnZSElJQUhICDw8PGwUIVVm/B4kIjKtpPt3SfhrPhERETklJjlERETklJjkEBERkVNikkNEREROiUkOEREROSUmOUREROSUmOQQERGRU3KoJOf8+fMYO3YsatasCS8vL7Rp0wYHDhywdVhERERkhxxmxeObN2/i0UcfRY8ePbB161b4+fnh9OnTqFatmq1DIyIiIjvkMD05b7/9NgIDA7Fq1Sp06NABwcHB6NWrF0JDQ20dmtMJDg7G4sWLDc9VKhXi4+NtGBEREVHpOUySs3HjRoSHh2PEiBHw8/ND27ZtsWzZshLPycnJQUZGhtHD2UyYMAERERHFvn7w4EEMHjwYfn5+8PDwQHBwMJ544glcu3bN7Pe4ePEiBgwYUBHhEhERWY3DJDlnzpzBJ598gsaNG2P79u147rnn8OKLL+Lzzz8v9pzY2Fj4+voaHoGBgVaM2PauXLmC3r17o1atWti+fTuOHz+OlStXIiAgAFlZWWZfx9/fHxqNxoKRPpiIIC8vz6xjc3NzLRaHJa9NREQVy2GSHL1ej3bt2uGtt95C27Zt8eyzz2Ly5Mn45JNPij3n9ddfR3p6uuGRlpZmxYht77fffkNGRgaWL1+Otm3bIiQkBD179sTixYsRFBRk9nUKD1edPXsWKpUK33//PXr06AEvLy+0bt0ae/bsKfLe3bp1g6enJwIDA/Hiiy/izp07htfXrFmD8PBw+Pj4wN/fH08++SSuXLlieD0hIQEqlQrbt29HeHg4NBoNEhMTi8SWH8/XX3+Nxx57DB4eHlizZo1ZMVy8eBGDBg2Cp6cnQkJCsG7dOpNDdUuXLsWwYcPg7e2NBQsWmP11IyIi23KYJCcgIADNmzc3amvWrBlSU1OLPUej0aBq1apGD3OJAHfu2OZRUfvC+/v7Iy8vD3FxcajozebnzJmD6dOn49ChQ2jSpAlGjx5t6Gk5cuQI+vXrh8cffxx//vknvvrqK/z666+YOnWq4XytVouYmBgcPnwY8fHxSElJwYQJE4q8z4wZMxAbG4vjx4+jVatWxcYzc+ZMvPjiizh+/Dj69etnVgxPPfUULly4gISEBHz33Xf49NNPjRKtfPPmzcOwYcNw5MgR/POf/yzHV42IiKxKHMTo0aOlS5cuRm3Tpk2TTp06mX2N9PR0ASDp6elFXrt7964cO3ZM7t69KyIimZkiSrph/Udmpvlfl/Hjx8uwYcOKfX327Nni5uYmNWrUkP79+8s777wjly5dKvGaDRo0kPfff9/wHIDExcWJiEhKSooAkOXLlxteP3r0qACQ48ePi4jIuHHj5JlnnjG6ZmJiori4uBi+vvf7/fffBYDcvn1bRER27twpACQ+Pr7EWPPjWbx4sVH7g2I4fvy4AJB9+/YZXj958qQAKPLZp02bVmIMFeX+70EiIlKUdP8uicP05Lz88svYu3cv3nrrLZw6dQrr1q3Dp59+iilTptg6NLu2cOFCXLp0CUuXLkXz5s2xdOlShIWF4ciRI+W6buFelYCAAAAw9IIcOHAAq1evRpUqVQyPfv36Qa/XIyUlBYBSED1s2DA0aNAAPj4+eOyxxwCgSM9ceHi4WfHcf9yDYkhKSoKbmxvatWtnOKdRo0aoXr36A69NRESOwWHWyXn44YcRFxeH119/HdHR0QgJCcHixYsxZswYi7yflxeQmWmRS5v13hWpZs2aGDFiBEaMGIHY2Fi0bdsWixYtwmeffVbma7q7uxv+rlKpACh1U/l/Pvvss3jxxReLnBcUFIQ7d+6gb9++6Nu3L9asWYPatWsjNTUV/fr1g1arNTre29vbrHjuP+5BMSQlJZm8jpgY1jM3BiIisi8Ok+QAwODBgzF48GCrvJdKBTjjvU2tViM0NNSoALeitWvXDkePHkWjRo1Mvn7kyBFcu3YN//rXvwwz3vbv32/VGMLCwpCXl4eDBw+iffv2AIBTp07h1q1bFRoHERHZjsMMV1Hx0tPTcejQIaNHamoqNm/ejLFjx2Lz5s1ITk5GUlISFi1ahB9++AHDhg2zWDwzZ87Enj17MGXKFBw6dAgnT57Exo0b8cILLwBQelLUajWWLFmCM2fOYOPGjYiJibFqDGFhYejduzeeeeYZ/P777zh48CCeeeYZeHp6GnqmiqPT6RAWFoZNmzZVaMxERFSxHKonh0xLSEhA27ZtjdrGjx+PqKgoeHl54dVXX0VaWho0Gg0aN26M5cuXY9y4cRaLp1WrVti1axfmzJmDrl27QkQQGhqKJ554AgBQu3ZtrF69GrNnz8aHH36Idu3aYdGiRRg6dKjVYgCAzz//HE8//TS6desGf39/xMbG4ujRo/Dw8Cjx2iKCpKQkpKenV1i8RERU8VRiqgjBSWVkZMDX1xfp6elFppNnZ2cjJSUFISEhD7zJkXP6+++/ERgYiP/+97/o1auX1d+f34NERKaVdP8uCXtyqNLasWMHMjMz0bJlS1y8eBEzZsxAcHAwunXrZuvQiIioAjDJoUorNzcXs2fPxpkzZ+Dj44POnTtj7dq1RjPHiIjIcTHJoUqrX79+6Nevn63DICIiC+HsKiIiInJKTHKIiIjIKTHJISIiIqfEJIeIiIicEpMcIiIickpMcoiIiMgpMcmhCnH27FmoVCocOnQIgLLVhEql4oaXRERkM0xyHNjSpUvh4+ODvLw8Q1tmZibc3d3RtWtXo2MTExOhUqmQnJwMAAgODsbixYuLvfZ3332Hjh07wtfXFz4+PnjooYfw6quvmh1b586dcfHiRfj6+pbyUxEREVUMJjkOrEePHsjMzMT+/fsNbYmJifD398e+ffuQlZVlaE9ISEDdunXRpEmTB173v//9L0aNGoXhw4fj999/x4EDB7Bw4UJotVqzY1Or1fD393/gjt6Wlpuba9ZxOp0Oer3eIjGIiFEiSkRE1sEkx4E1bdoUdevWRUJCgqEtISEBw4YNQ2hoKH777Tej9h49eph13c2bN6NLly547bXX0LRpUzRp0gQRERFYsmSJ2bHdP1y1evVqVKtWDdu3b0ezZs1QpUoV9O/fHxcvXjQ6b9WqVWjWrBk8PDwQFhaGjz/+2Oj1mTNnokmTJvDy8kLDhg0xd+5co0Rm/vz5aNOmDVauXImGDRtCo9HA1B60+fFs3rwZzZs3h0ajwblz58yK4bfffkObNm3g4eGB8PBwxMfHmxyq2759O8LDw6HRaJCYmGj2146IiCoGt3UojghQqCfEqry8ADN7QB577DHs3LkTs2bNAgDs3LkTM2bMgF6vx86dO9G7d29otVrs2bPH7CTF398f69atw19//YUWLVqU+WPcLysrC4sWLcIXX3wBFxcXjB07FtOnT8fatWsBAMuWLcO8efPwn//8B23btsXBgwcxefJkeHt7Y/z48QAAHx8frF69GnXr1sWRI0cwefJk+Pj4YMaMGYb3OXXqFL7++mt89913cHV1LTGe2NhYLF++HDVr1oSfn98DY7h9+zaGDBmCgQMHYt26dTh37hymTZtm8vozZszAokWL0LBhQ1SrVq3Cvo5ERGQmqUTS09MFgKSnpxd57e7du3Ls2DG5e/eu0pCZKaKkOtZ/ZGaa/Zk+/fRT8fb2ltzcXMnIyBA3Nze5fPmyrF+/Xjp37iwiIrt27RIAcvr0acN5DRo0kPfff9/kNTMzM2XgwIECQBo0aCBPPPGErFixQrKzs4uNIyUlRQDIwYMHRURk586dAkBu3rwpIiKrVq0SAHLq1CnDOR999JHUqVPH8DwwMFDWrVtndN2YmBjp1KlTse/7zjvvSPv27Q3P582bJ+7u7nLlypVizykcz6FDh4zaHxTDJ598IjVr1iz4PhGRZcuWmfzs8fHxJcZwvyLfg0REJCIl379Lwp4cB9ejRw/cuXMH+/btw82bN9GkSRP4+fmhe/fuGDduHO7cuYOEhAQEBQWhYcOGZl3T29sbW7ZswenTp7Fz507s3bsXr776Kj744APs2bMHXl5eZYrVy8sLoaGhhucBAQG4cuUKAODq1atIS0vD008/jcmTJxuOycvLMype/vbbb7F48WKcOnUKmZmZyMvLQ9WqVY3ep0GDBqhdu/YD41Gr1WjVqpXhuTkxJCUloVWrVvDw8DC83qFDB5PXDw8Pf2AMRERkOUxyiuPlBWRm2u69zdSoUSPUr18fO3fuxM2bN9G9e3cAypBTSEgIdu/ejZ07d6Jnz56lDiM0NBShoaGYNGkS5syZgyZNmuCrr77CxIkTS30tAHB3dzd6rlKpDPUy+UW/y5YtQ8eOHY2Oyx9y2rt3L0aNGoU333wT/fr1g6+vL9avX493333X6Hhvb2+z4vH09DQqjDYnBhEpUkyd/xnuZ24cRERkGUxyiqNSAQ5yk+rRowcSEhJw8+ZNvPbaa4b27t27Y/v27di7d2+ZE5N8wcHB8PLywp07d8obrkl16tRBvXr1cObMGYwZM8bkMbt370aDBg0wZ84cQ1t+sbC1YggLC8PatWuRk5MDjUYDAEaz24iIyH4wyXECPXr0wJQpU5Cbm2voyQGUJOd//ud/kJ2dbXJm1fnz5w0zgvIFBQXhww8/RFZWFgYOHIgGDRrg1q1b+PDDD5Gbm4s+ffpY7HPMnz8fL774IqpWrYoBAwYgJycH+/fvx82bN/HKK6+gUaNGSE1Nxfr16/Hwww9jy5YtiIuLs2oMTz75JObMmYNnnnkGs2bNQmpqKhYtWgQAD5wuv2fPHkycOBG7du1CnTp1KjRuIiIqilPInUCPHj1w9+5dNGrUyOjm2b17d9y+fRuhoaEIDAwsct6iRYvQtm1bo8fGjRvRvXt3nDlzBk899RTCwsIwYMAAXLp0CT/++COaNm1qsc8xadIkLF++HKtXr0bLli3RvXt3rF69GiEhIQCAYcOG4eWXX8bUqVPRpk0b/Pbbb5g7d65VY6hatSo2bdqEQ4cOoU2bNpgzZw6ioqIAwKhOx5Q7d+4gKSnJ7LV7iIiofFRSXEGBE8rIyICvry/S09OLFKtmZ2cjJSUFISEhD7xZERW2du1aTJw4Eenp6fD09Czzdfg9SERkWkn375JwuIqolD7//HM0bNgQ9erVw+HDhzFz5kyMHDmyXAkOERFVPCY5RKV06dIlREVF4dKlSwgICMCIESOwcOFCW4dFRET3YZJDVEozZswwWmGZiIjsEwuPiYiIyCkxyblPJarDJjvD7z0ioorFJOee/BVttVqtjSOhyirr3oaw968MTUREZcOanHvc3Nzg5eWFq1evwt3dHS4uzP/IOkQEWVlZuHLlCqpVq1bizulERGQ+Jjn3qFQqBAQEICUlpUK3CiAyV7Vq1eDv72/rMIiInAaTnELUajUaN27MISuyOnd3d/bgEBFVMCY593FxceFqs0RERE6AhSdERETklJjkEBERkVNikkNEREROiUkOEREROSUmOUREpTB//nzExMSYfC0mJgbz58+3bkBEVCwmOUREpeDq6oqoqKgiiU5MTAyioqK4FACRHeEUciKiUpg7dy4AICoqyvA8P8GJjo42vE5EtqeSSrQrYEZGBnx9fZGeno6qVavaOhwicmD5iY1arYZWq2WCQ2RBZb1/M8khIiojjUYDrVYLtVqNnJwcW4dD5LTKev9mTQ4RURnExMQYEhytVltsMTIR2Q6THCKiUipcg5OTk4Po6GiTxchEZFssPCYiKgVTRcamipGJyPYcticnNjYWKpUK06ZNs3UoRFSJ6HQ6k0XGc+fORXR0NHQ6nY0iI6L7OWTh8b59+zBy5EhUrVoVPXr0wOLFi806j4XHREREjqfSFB5nZmZizJgxWLZsGapXr27rcIiIiMhOOVySM2XKFAwaNAi9e/d+4LE5OTnIyMgwehAREVHl4FCFx+vXr8cff/yBffv2mXV8bGws3nzzTQtHRURERPbIYXpy0tLS8NJLL2HNmjXw8PAw65zXX38d6enphkdaWpqFoyQiIqrE7Kzw3mEKj+Pj4xEZGWm0+Z1Op4NKpYKLiwtycnIeuDEeC4+JiIgqmAjw++/AsmXArl3AsWOAu3uFvkVZ798OM1zVq1cvHDlyxKht4sSJCAsLw8yZM7nzLxERkTXdvAmsWaMkN4Xvzz/9BAwcaLu4CnGYJMfHxwctWrQwavP29kbNmjWLtBMREZEFiAC//KIkNt9+C+Tv2ebhAYwcCUyeDDz6qG1jLMRhkhwiIiKykStXgM8+A5YvB5KTC9pbt1YSmzFjgGrVbBdfMRw6yUlISLB1CERERM5Jrwf++1+l12bDBiA3V2mvUgUYPVpJbsLDAZXKtnGWwKGTHCIiIqpg588Dq1YBK1YAZ88WtHfooCQ2TzwB+PjYLLzSYJJDRERU2eXlAVu3Kr02W7YovTiAMgQ1dqyS3LRqZdsYy4BJDhERUWV19qzSY7NyJXDhQkF7165KYjN8OODpabPwyotJDhERUWWi1QIbNyq9Nj/9pMyYAoBatYDx44FJk4CwMNvGWEGY5BAREVUGycnK7KjVq4GrVwvae/dWem2GDQM0GpuFZwlMcoiIiJxVdjbw3XcFqxHnCwgAJk4Enn4aaNjQdvFZGJMcIiIiZ/PXX0pi88UXysrEAODioqxEPHmy8qeb86cAzv8JiYiIKoM7d4CvvlKSm717C9obNFB6bCZOBOrXt118NsAkh4iIyFGJAAcOKInNl18Ct28r7W5uSo3N5MlKzU0l3d+RSQ4REZGjSU8H1q5VkptDhwraGzVSEpvx44E6dWwXn51gkkNEROQIRIDfflMSm6+/Bu7eVdo1GuAf/1CSm+7d7XqbBWtjkkNERGTPrl8HPv9cmf597FhB+0MPKYnN2LFAzZq2i8+OMckhIiKyN3o9kJCg9Np8/72ygB8AeHkpe0dNngw88gh7bR6ASQ4REZG9uHRJWaxv+XLg9OmC9nbtlMRm9GjA19fsy82fPx+urq6YO3dukddiYmKg0+kwf/788sdtp1xsHQAREVGlptMpm2M+/jgQGAi8/rqS4Pj4AM89p8yeOnBA+XspEhwAcHV1RVRUFGJiYozaY2JiEBUVBVcnn3XFnhwiIiJbSEtTNsZcuRJITS1o79RJ6bUZORLw9i7XW+T34ERFRRme5yc40dHRJnt4nIlKJH9nLueXkZEBX19fpKeno2rVqrYOh8hpVPYucSKz5eYCW7YotTbbtim1NwBQowYwbpyyOWaLFhX+tvmJjVqthlardbgEp6z3bw5XEVG5VfYucaIHOn0amD0bCAoCIiOBH35QEpzHHlPWuzl/Hli82CIJDqD04OQnOGq12qESnPLgcBURlVtl7xInMiknB4iPV3ptfv65oN3Pr2BzzMaNrRJKTEyMIcHRarWIiYmpHD+XUomkp6cLAElPT7d1KEROKTo6WgCIWq0WABIdHW3rkIis79gxkZdfFqlZU0RZwk9EpRLp31/ku+9EtFqrhpP/c5n/83j/c0dQ1vs3a3KIqEJpNBrDb4w5OTm2DofIOrKygG++UXptdu8uaK9fH/jnP5VHgwZWD6u4HtWK7GnNysrChQsX0KhRo/KGW6yy3r85XEVEFabSdolT5XXokJLYrF2r7CcFKJthDh6szJDq37/CNscsS4G/TqczmcjkP9fpdGWK5caNG9i8eTPi4uKwfft2tGzZEv/3f/9XpmtZlEX6lewUh6uILMcZusSJzJJiv2SLAAAgAElEQVSRIfK//ysSHl4wHAWIhISILFwocv68Rd62uJ8pa/2spaWlyZIlS6RXr17i6uoqAAyPkJAQyczMtNh7l/X+zSSHiMrN1v/5ElmcXi+yd6/I00+LeHsXJDbu7iIjR4r89JOITmfxMGz5y8Tw4cONEpuWLVtKVFSUHDx4UPR6vUXfu6z3bw5XEVG5WapLnMjmbt4E1qxRhqSOHClob9pUGY566imgdm2rhVN4JuOCBQsqfM0bvV6Pffv2IS4uDnFxcdiwYQPCwsIAAI8//jguXryIiIgIREZGIjQ0tELe05JYeExERFSYCJCYqCQ2334LZGcr7R4ewIgRSnLTpYtNN8esyAL/3NxcJCQkGJKaCxcuGF6LjY3FrFmzyhtuubHwmIiIqDyuXgU++0zZHDMpqaC9VSslsRkzBqhe3Xbx3VORBf4HDx5Ez549cevWLUObj48PBg4ciIiICAwcOLCiwrYJJjlERFR56fXKQn3LlikL9+XmKu3e3sqO35MnAw8/bNNem8Lun/qd/xzAAxOda9euYdOmTVCr1RgzZgwAoFmzZsjNzYWfnx+GDRuGyMhI9OzZExqNxuKfxRqY5BARUeVz4QKwahWwYgWQklLQ/vDDSmIzapSyC7gdMbW2janVxgs7d+4c4uPjERcXh8TEROj1eoSFhRmSHA8PD+zfvx+NGzd2yu1XmOQQEVHlkJenbIq5bJmySWZ+QbyvLzB2rJLctG5t2xhLUJoC/yVLlmD16tX4448/jI5t06YNIiIikJeXBzc3JQXILyx2RkxyiIjIuZ09C6xcqTzOny9o79pVSWyGDwc8PW0WnrnuX+gvn16vR58+fdCxY0dD28GDB/HHH3/AxcUFXbp0QUREBCIiIhASEmKlaO0DkxwiInI+Wi2wcaPSa/PTT8qMKQCoVQsYPx6YNAlw4B4MrVaLHTt2ID4+Hhs2bMClS5fw+++/4+GHHwYAPPvss3j00UcxdOhQ1LbiFHd7wySHiIicR3KyMjvqs8+AK1cK2nv3Vnpthg0DHLSoNjMzE1u3bkVcXBy2bNmCjIwMw2tVq1bF6dOnDUlOx44djXp2KismOUTklMqyzw85qOxs4LvvlF6bXbsK2gMCgIkTgaefBho2tF185SAiUN2b2bV//36MHDnS8Jq/v79hRlSPHj2gVqttFabdYpJDRE7J1dXV5IyTwjNUyMH99ZeS2HzxhbIyMQC4uAADBii9NoMGAW7m3+bsJTFOSUkxzIhq164dFi9eDADo0qULHnnkEXTt2hWRkZHo2LEjXFxcLB6PQ7PEHhP2intXETm2efPmFbtHT3R0tMybN69IG7hpqHPJzBRZsULkkUeMN8cMChJ5802R1NQyX9pWe7Dp9Xo5dOiQzJ8/X1q3bm20P1T9+vUtvi+UI+AGnWZgkkPk2MpyE8p/Ta1WM8FxZPv3izz7rIiPT0Fi4+Ym8vjjIlu3iuTlVcjb2CIx7tGjh1Fi4+LiIo899ph88MEHcu7cOYu9ryNhkmMGJjlEjq8sN6H8BEetVlsrTKoIt26JfPyxSNu2xr02jRqJ/OtfIhcvWuRtLZUYZ2dny5YtW2Tq1KmSm5traH/ppZfEw8NDhg4dKqtWrZKrV69WyPs5EyY5ZmCSQ+QcSnMTYk+Og9HrRXbvFpkwQcTLqyCxUatFRo8W2bFDRKezeBgVlRinp6fLl19+KSNHjpQqVaoYemt+/vlnwzFXrlyRzMzM8obs1JjkmIFJDpHzMOcmxJocB3Ltmsj774s0b27ca9O8udJ+7ZrVQqmIxHjfvn3Sv39/cXd3NxqKqlu3rjz//PNy9OhRC0TuvJjkmIFJDpFzMOcmZKsiUioFvV7pmRk9WumpyU9sPD2Vnpzdu5VjrKisifGpU6ckOTnZ8PzAgQOGxKZp06Yya9Ys2bt3r+is0AvljJjkmIFJDpHjM/cmVNqZWGRFFy8qNTWNGhn32rRtq9Tg3Lplk7BKkxjr9Xr5448/ZO7cudKyZUsBIOPHjzd6/f3335djx45ZK3ynVtb7N9fJISKHUZpdmEtaz8TUOihkYTod8OOPyro2mzYpm2UCyk7fY8Yo69q0a2fjEEveADMvLw+7du1CXFwc4uPjce7cOcMxrq6uyMrKMjxXqVSYNm2adQKnYqlE8jf0cH4ZGRnw9fVFeno6qlatautwiKzKXhY6Kw9n+AyVTlpaweaYqakF7Z06KYnNyJGAt7ft4nsAnU4HV1dXw/MWLVrg6NGjAABPT0/069cPkZGRGDx4MGrUqGGrMJ1eWe/f7MkhqiScYQVg9s44iNxcYMsWpddm2zZAr1faq1cHnnpK2RyzRQvbxliCW7du4YcffkBcXBx++eUXpKSkwMvLCwAwYcIEHDlyBJGRkejbt6+h3Vk43S8SFhk8s1OsyaHKjrONyKJOnxZ5/XURf3/jWpvHHhNZu1bk7l1bR1is8+fPyyeffCJ9+/YtMiNq06ZNtg7Pauy1YN/pC4/feustCQ8PlypVqkjt2rVl2LBhcuLEiVJdg0kOEdeNoQqWnS2yfr1Ir17GiY2fn8iMGSKFZhzZq88++8woqQEgzZo1k9mzZ8u+ffsq3bYK9vjLUFnv3w5Tk9O/f3+MGjUKDz/8MPLy8jBnzhwcOXIEx44dg7eZ47msySFSaDQaaLVaqNVq5OTk2DocckQnTijDUZ9/Dly7prSpVEDfvkqtzZAhgJ3tii0i+OOPPxAXF4dOnTph0KBBAICTJ0+iSZMm6NixIyIjIxEREYGmTZvaOFrbyh/GVqvV0Gq1Jguyrams92+HSXLud/XqVfj5+WHXrl3o1q2bWecwySGyv/+8yIFkZQHffqskN7/+WtBerx7wz38qj+Bgm4VnSl5eHhITEw0zotLS0gAAkZGR+P777w3HXb58GXXq1LFVmHbJnn4ZqnSFx+np6QBQYjV7Tk6O0T9MRkaGxeMismf3T8HOfw6wcJdKcPiwktisWQPc+78Xrq7AoEFKr03//oCbfd1OdDodJk+ejA0bNuDGjRuGdi8vLwwYMACjRo0yOp4JjrGYmBhDgqPVahETE+OY/0dU/MiZ5en1ehkyZIh06dKlxOPmzZtXZJwVrMmhSspeCwrJTmVkiHz6qcjDDxvX2gQHiyxYIHL+vK0jNHLjxg3Zvn27UdsjjzwiAKRmzZoyceJE2bhxo2RlZdkoQsfhTDU59pV6m2nq1Kn4888/8Wvh7lITXn/9dbzyyiuG5xkZGQgMDLR0eER26UELnel0OluERfZEBNi3T+m1+fJL4M4dpd3dHYiIUHptevUCXFxsG+c958+fR3x8POLj45GQkAAAuHLlCqpXrw4AWLBgAdzc3PDoo4/Czc56muxVaRbcdAQO96/+wgsvYOPGjfjll19Qv379Eo/VaDTQaDRWiozIvnGNGSrWzZvKUNSyZcCRIwXtTZsqa9qMHw/Urm27+ApJSUnB+vXrER8fj99//93otRYtWiA1NdWQ5PTq1csWITo0Z/tlyGEKj0UEL7zwAuLi4pCQkIDGjRuX+hosPCai8nCqhdJEgMREJbH59lsgO1tp9/AAhg9Xem26dlVmTNmQXq9HXl4e1Pdman3yySd4/vnnAShbJ3Tq1MkwI6pRo0a2DJUsyOkLj6dMmYJ169Zhw4YN8PHxwaVLlwAAvr6+8PT0tHF0RFQZOMOq0bh6FfjsM2D5ciApqaC9VSslsRkzRlmZ2IZyc3Oxa9cuw1DU7NmzDYlNREQENmzYgMjISAwbNgz+/v42jZXsnCUKhCwBJgqIAciqVavMvgYXAySi8rLHoswH0ulEfvxRZMQIEXf3giJib2+RSZNE/u//RGy84F1mZqZ89913Mm7cOKlevbrR//ODBg2yaWxke05feCyOMapGRE6ucBHmggUL7HutoQsXgFWrgBUrgJSUgvaHH1Z6bUaNUnYBt7Hs7GzUrVvXaJmP2rVrY+jQoYiMjGRtDZWZw9TkVATW5FBZOFUdBlUYe1oozUhenrIp5rJlyiaZ+YWivr7A2LFKctO6tc3CS01NRXx8PJKSkvDRRx8Z2vv27YuTJ08iMjISkZGR6Ny5s9Hu31S5OX1NDpGtOEUdBlUou1wo7dw5pcdm5Urg/PmC9i5dlMRm+HDABjtmiwiOHTuG+Ph4xMXF4cCBA4bXZs+ejXr16gEAvv76a/j6+kJl40JncjIWGDqzW6zJobJyyDoMsgi7+l7QakW+/VakXz8Rlaqg1qZmTZFXXhE5dsz6MRWydu1aady4sVF9jUqlkq5du8q7774r169ft2l8+ebNm1fsv190dLTMmzfPugFREU6/C3lFYJJD5cHdu8luVo1OTlZ2+PbzM16NuFcvZUfw7GzrxFFITk6ObNu2TdLS0gxta9asMfzMDBo0SJYvXy6XL1+2emwPYjf/rlQsqyU548ePl127dpX2NLvAJIfKKz/BUavVtg6FbMCmv/HfvSuydq3IY48ZJzb+/iKvvy5y6pTl3rsYt2/flm+++UaefPJJ8fX1FQASGxtreP3WrVvy1VdfSUZGhtVjKy276qGjIqyW5Dz++OOi0WikUaNGsnDhQvn7779LewmbYZJD5cGeHLKJv/4SeeklkRo1ChIbFxeRgQNF4uKUISsrysrKkhUrVsjgwYNFo9EYDUXVqVNH3nnnHavGU5H4M14+er3I1auWubZVh6uuXbsmixcvljZt2oibm5v0799fvvnmG9Fa+YettJjkUFnxtzyyqsxMkRUrRB55xLjXJjBQZP58kdRUq4Zz+/Ztw9+zsrLE29vbkNiEhobK9OnTZffu3aLT6awalyWwt/bBtFqRpCSRDRtE3n5bZOJEkc6dC/JwSyQ6NqvJ+eOPP2Tq1Kni4eEhtWrVkmnTpklycnJ5L2sRTHKoLDheT1azf7/Is8+K+PgUJDZubiKRkSI//CCSl2eVMPR6vfz555/y5ptvSps2baR58+ZGr8+aNUuio6PlyJEjorfxIoIViT05xm7cENmzR2TVKpFZs0QiIkTCwpRvycK5d+GHSqWsLVnRbLIY4MWLF/Hjjz/ixx9/hKurKwYOHIijR4+iefPmeOedd/Dyyy+X5/JEdsHZNqwjO5OeDqxbp6xrc/BgQXtoqLI55oQJgBW2LtDr9dizZw/i4uIQHx+P06dPG15zcXFBWloaAgMDAQCxsbElXssR15a6f/ft/OeAc29gq9Mpqw8kJQEnThQ8kpKAy5eLP8/LS9m/NSxMeeT/vXFjm6xUULzSZlNarVa+/fZbGTRokLi7u0v79u3lk08+MSos+/LLL6VatWqlvbTFsSeHiOyCXi+ye7fIhAkiXl4Fvwar1SKjR4vs2KFsxWBFU6dONaqv0Wg0MmTIEFm5cqVcLeX4g6P1fjpavGVx+7bSUbhmjcgbbyg7fLRsKaLRFN8rA4jUq6dM2nv+eZEPP1R2B0lNtfq3p/V6cgICAqDX6zF69Gj8/vvvaNOmTZFj+vXrh2rVqpUn9yIicj7XrwNffKFsjnn0aEF78+bKgn3jxgE1a1o0hIyMDGzduhVxcXGYOXMm2rZtCwDo06cPvvjiCwwePBiRkZHo168fqlSpUqb3KLz1Rf7z+3tK7Imz9NaKAH//XdATU7hnpvD6kPfTaJQemPxemfyemaZN7WLXj3Ip9bYOX3zxBUaMGAEPDw9LxWQx3NaBiKxp/vz5cHVxwdyuXZXhqO+/B+5tAZHr7o6/HnoIbT/6COjUCbDgSr+XL1/Gxo0bERcXh59//hlarRaAsuLwwoULlXhycyEiUKvVFfa++YlN/srQ9pjgOKK7d4GTJ4sOLyUlAXfuFH+en1/R4aWwMKBBA8Ded9Ao6/2be1cREVnC5cv4edw4BP30ExoXbm/bFj/Uq4fRmzdjuoVv+hcuXMDIkSPx22+/GW1y3KRJE0RGRmLUqFEme+Mrkt3u8WXnRJSamMJJTP7fz51TXjfFzQ1o1Mg4iclPaqpXt+5nqEjcu4qIyNZ0OuCnn5Rem40b0SsvDwCQAeBkeDjaL12KmB9+sMiwjYjg8OHDuHDhAgYOHAgA8PPzw4kTJyAiCA8PR2RkJCIiItCsWTOr7BFll3t82RmtFjh1qmgik5Sk1KQXp3p14yQmP5Fp2BBwd7de/PaOSQ4RUXn9/beyMeaKFUBqakH7I48Akyfj4zNn8PrChVB37lyhwzY6nQ67d+82zIg6e/YsAgMDce7cOahUKri5uWH9+vVo2rSpYWaUtVTW2UrFuXbN9AymM2cKNoq/n4sLEBJieoipVi2LjnA6DSY5RERlkZcHbNmi9Nps3Qro9Up79epKAfGkSUDLlgCAWQDm/fvfhl6N8t7kd+7ciTVr1mDjxo24du2aod3DwwPt2rVDenq6YfJH7969y/VeZWGqyNhUMbKzycsDUlJMDzFdv178eT4+RZOYsDBl2EmjsV78zohJDhFRaZw5o/TYrFoFXLxY0N69uzJD6vHHAU9Po1PKO2yTnp4Ob29vuLkp/2XHx8dj5cqVAIDq1atjyJAhiIiIQN++feHt7V3+z1hOzjJbqTi3bhUU+hbumTl1CsjNLf68oKCiw0thYUBAAHtlLIVJDlEl5IiLtdlUTg4QH6/02vz8c0F77drKYn2TJgFNmpg8tazDNhcvXjTMiNqxYwe2bduGnj17AgBGjRoFnU6HyMhIdOvWDe52VoRR0veOo/Tg6PXKyOP9w0snTgCXLhV/nqdnQfJSuGemcWPADvLPSodJDlEl5OrqavJGW/iGTFDuaMuWAZ9/rhRVAMqv3H36KL02Q4cCJUy5Lu2wzcmTJw31NXv37jWaEZWYmGhIcjp16oROnTpV6EetrDIzgeTkosNLyclAdnbx59Wta3oGU2CgUktD9oFJDlEl5GiLtVnV3bvAN98oyc2vvxa0160L/POfwNNPA8HBZl2qNMM2R44cQatWrYyO69Chg2FGVFhYWNk+D0FEWQzP1CJ5f/9d/HlqtfEieYV7aLgKiWNgkkNUSRVOdBYsWMDF2g4fVhKbNWsK5u66ugKDBinDUQMGKIuQlIKpYZu8vDz8+uuvuHbtGrwKbfLTokULhIaGIiQkBJGRkRg6dCjq169fnk9U6WRnF79IXmZm8efVrm268Dc42P4XyaOScTFAokquUi/Wdvs2sH69ktzs21fQHhxcsDlmvXrlfpu7d+/ip59+QlxcHDZt2oTr96ba1KhRA5cvXzYUFOfm5tpdfY29EQGuXDE9g+ns2eIXyXN1VfY8NbV1gYV30qAKwMUAiajUKuVibSJKQrNsmZLg5P+K7+4OREQotTa9elVYYcVrr72Gjz/+GFlZWYa2GjVqYOjQoYiIiDA6lglOAa0WOH3a9Noyt24Vf161aqZnMDVsWGL5VBEszncOTHKIKqlKt1jbrVvKUNSyZcCffxa0N2miJDZPPaVs7lMO58+fx8aNGzFhwgR43ptGrtFokJWVhcDAQEN9TdeuXQ29N5Xd9etFV/o9cUJJcIqbaa5SFSySd/8QU+3aFTMdm8X5TqIit0K3d2Xdqp3I2URHRwsAiY6ONqvdYen1Ir/8IjJunIiHh4jSjyOi0YiMHSuya5dyTDmcOHFCYmNjpWPHjgJAAEh8fLzh9XPnzsn+/ftFX873cWS5uSInT4ps2iSyaJHIpEkiXbqI1KpV8E9i6lGlikj79iJjxojExIh8843In3+K3L1rnbjv/3lwup8PB1LW+zd/lSCqhJx9sTZcvapM+16+XOkWyNeypdJrM3ZsuXYr/Pvvv/Hxxx8jLi4OJwpfH8r07sI7eQcFBSEoKKjM7+VIMjJMDy+dPKkMPxUnMND01gV169p2kTwW5zs+Fh4TkXPQ64EdO5ThqLi4gqVnvb2BUaOU5KZDhzLdNXNzc3Hz5k343RvOSk5ORtOmTQEodTQ9e/Y0zIgKCAiosI9kj/R6IC3NdOFv4QWg7+fhUVDoW3h4qUkT+18kr1IX59sJFh4TUeV08aKyxcKKFcqWC/nCw5XEZtSoMi1qkpWVhR9//NEwI6pHjx747rvvAABNmjTBSy+9hA4dOmDgwIGGfaKcyZ07yoJ49/fMJCcrSwkVJyDAOJHJ/3tQkGMuklcpi/OdCJMcIjtm7gyPSjcTRKdTNsVctkzZJDN/eM3XFxgzRklu2rQp9WVv3LiBzZs3Iy4uDtu3b8fdQnfzP/74AzqdDq73Fk5ZvHhxhXwUWxIBLlwwPcRUeDP1+7m7F79Inq+v9eK3tEpXnO+EmOQQ2TFzZ3hUlpkg70+bhvaHDqHbqVPKErb5Hn0UG/z88FezZpizcGGZrz906FDs3r3b8Dw4ONgwI+rRRx81JDiOJjtb2TzS1BBTSYvk1apV/CJ5zj45rLLupO5snPzblMixmbv9glNv05CbC2zcCCxbhmk//ghVfhlhzZrKtO9JkxDz3XfKZ23b1qxLHj9+3DAM9cMPP6D6vSLkoUOHIiMjA5GRkYiMjETr1q2hcpDtoUWUemtTWxecPavU0pji6qqsIWNqkbxataz6EeyK0xfnVxYWmetlpziFnBxV/tRVtVpd4hRWc49zCMnJIjNmiPj5Gc0rPhMSIk8AsjAqSkTMm9ar0+lk7969MnPmTGnatKlhqjcA+eKLL4yOs3darciJEyLx8SL/+pfIhAkinTqJVK9e8nRsX1+Rjh1Fxo8Xeestke+/Fzl2TCQnx9afiOjBynr/5uwqIgdh7gyP4o5ziLqd7Gzg+++VWpuEhIJ2f39g4kRlc8zQUEMvVX4xaEm9VYmJiRg1ahQuXLhgaFOr1ejduzciIiIQERGB2rVrW/iDld7Nm6aHl06fBvLyTJ+jUilDSaZ2x65Tx7bTsYnKg7OriJyYuTM8SjrOrut2jh5VEpsvvgBu3FDaVCplU8zJk5VNMgtteTB37lzDuiVqtdrweTIzM7Ft2zZUq1YNvXv3BgCEhobiwoUL8PHxwcCBAxEZGYkBAwbYxS86Op0ylHR/IpOUpOzPVBxvb9MzmBo3Bu4ttExEAIeriOyduauumnOcXa3gmpkpsnKlMtZSeFwlMFBk/nyR1NRiT71/WC4iIkKGDBkiHh4eAkD69OljdHxiYqJkZ2db+hMVKyNDZN8+kS++EJkzR+Qf/xBp0UJZeLmkIab69UV69xaZMkVkyRKRn34SSUsr9yLNRA6nrPdvJjlEdqw8CU1x7Tav2zlwQOS550SqVi24m7u6ikRGivzwg0heXomn58c/ePBg6d69u6hUKqMam4YNG8qsWbOsvo2CTidy7pzI9u0iH3wg8vzzIj17itStW3Iio9GItGwpMmKEyNy5ImvXiuzfryRGjmrevHkl1o3NmzfPugGRw+O2DkROyNwZHqWZCVLcUI9FZWQA69YpQ1J//FHQHhoKTJoETJig1N0UQ0Rw+vRpfPnll4ahtcTEROzatQsAEBAQgIsXL2LKlClYsmSJRWdEZWUVv0heoY3Gi6hTx/Tu2EFBygwnZ2LXQ6NUuVgm57JP7MkhsmJPjl4v8ttvIhMninh5FXRdqNUio0aJ/Pyz0v1RDJ1OJ7t375bp06dLaGioAJCXX37ZEO+GDRvk/fffl5SUFMPnqqgeAr1e5MIFkR07RD7+WOTFF0X69RNp0KDkXhk3N5FmzZROqVmzRFavFtm7V+TmzQoJy6HY1dAoOTwOV5mBSY7jYHe3ZVjlxnP9usjixSIPPWScATRrJvLeeyJXrxZ7ak5OjmzdulWeeeYZqVOnjtEwlEajkU2bNlVcnCKSnS3y118i334rsnChsln5ww+L+PiUnMzUqCHSubPIP/8p8s47Ihs3iiQlKdO7qYDNh0bJaTDJMQOTHMdRmhoTMo9Fv6Z6vcjOnSJPPmlcTevpqSzM8uuvZlXLrlmzxiix8fX1lSeffFK++eYbuX37dplDu3JFJDFRZNkykVdfFRk0SKRRIxEXl+ITGRcX5ZjBg0WmT1fOTUwsMUcjE/ITHLVabetQyIGxJoecilOv4GsjFlnB9fJlYPVqYPlyZd+AfG3aKFO/n3wSMLF55ZUrV7Bx40bExcWhT58+mDZtGgBg4MCBCAoKMkz1fuyxx6BWq80KJTcXSEkxvbZM/qx0U6pWNb11QWgooNGU5otB9+PmlmRzFkq67BJ7chwPu7vtUF6eyNatIo8/rhSh5Hd9VKki8swzylxpE702Z86ckXfffVe6du0qLi4uht6aRx55xOi4B82KunlTZM8epd5l1iyRiAiRsDARd/fie2VUKpHgYKWu5qWXRD75RKm3uXCB07EthTU5tuOMw/0crjIDkxzHxO5uO5GWJvLmm0Wrbzt2FFm+XKSY4SS9Xi+PPvqo0TAUAGnfvr0sWLBA/vrrryLn5OWJnD6tzCh/7z0ld+rWTaROnZJrZby8RNq2FRk9Wgl1/XqRQ4dE7tyx8NeGjHC42bac8evP4SpySuzutrG8PGDLFmXq99atBbs8VqsGjBunDEm1bGk4XKfT4bfffsOvv/6K119/HQCgUqlQp04duLq6olu3boiMjMSwYcMQFBSEzExlWGntWuMhpuRkoISdK1Cvnukhpnr1ABcXS35ByBzc3NK2ONxfiIWSLrvEnhzHwu5uGzpzRmT2bJGAAOOuku7dlWV7s7IMh2ZnZ8uWLVtk0qRJ4ufnZ+ipSUpKEhFlOCgx8ax8+226LFmirN7bq5dIvXoPXiSvRQuR4cNF3nhDZM0aZSTMkRfJI7ImZxrurzTDVR999JEEBweLRqORdu3ayS+//GL2uUxyHIczdrfez+7GzXNyRL76StlHoHC2Ubu2yGuvKXOkC9mzZ4+MHDlSqlSpci+x8RCgpXh5TZBWrb6VQYNuSbt2It7eJSczfn7KUNQzz4i8+67Ili3KUNUDFj4mIjM4y3B/pRiu+uqrrzBt2jR8/PHHePTRR/G///u/GDBgAI4dO4agoCBbh0cVqDJ0d9vNqrBJScpw1IygUXQAACAASURBVGefAdeuKW0qFdCnjzIcNXQooFbj0qVLwKXLAOrgxAngq6888fXXjwCYAFfXh6DT1Qfggqws4M8/lUc+NzegUSPTu2NXr26dj0lU2XC4H441XNWhQwd57rnnjNrCwsJk1qxZZp3PnhyyNzYbksvKEvn8c5GuXY27VerWFXnjDck5cUaOHhX5z38uSP/+u6R27R8E2CsaTVaJvTLVqyv7bU6cKPL22yLx8SInTnCRPCJrc7bhfqcfrsrJyRFXV1f5/vvvjdpffPFF6datm1nXYJJD9qgs4+ZlHuo6fFhk6lSRatUMmYnexUXOth4iKyI3ytCBWgkKuisqVV6Ji+SFhioL6r3yisinn4r88ouy4B6nYxPZnjMO9zv9cNW1a9eg0+lQp04do/Y6deoo3egm5OTkIKfQFI2MjAyLxkhUFg/aMHP+/PlwdXU1as8f6kpISEDXrl0xf/58AKaHuvJuZeL6R+uh/nwZqif/bmhPdWmAT/WTsEo/ERcO1wMO57/ifu/PDABJ8Pe/hfBwHwwZ0hidO9dEo0aAh0fFfx2IqGJUhuF+czlMkpPv/t2FRaTYHYdjY2Px5ptvWiMsojJ70Li5qdqduXPnIiEhATt27DAcN3v2O4iN3YzHH4/H3ayhmNFjH8IPLcOAW1+iDjIBALlwQzwisAyT8V99bwhc4O5+EX17CMLCVGjaFPjxxw+Rk3MYo0c/hsGDB6FGjRpW/GoQUXnl/9JjCmty7FRZhquys7MlPT3d8EhLS+NwFdkVc8fN89vnz4+WlBSRcePWCvCS1K27UYAdAlwQQMQXN+V5/EcOorXRGFOSqom86/8vGdB+nzRr9pVoNOMFaC2AlwCQffv22eDTExGZp6zDVSoREVsmWaXRsWNHtG/fHh9//LGhrXnz5hg2bBhiY2MfeH5GRgZ8fX2Rnp6OqlWrWjJUchKmhoryxcTEQKfTlfhbU0mKW5xLaX8bzz33Prp2nWxYIG/Xrsu4fLkqAM/7riTogl8xGcswAt/AE9kAgDw3Da50Gw6XZybjm8tH8MqrLyMvL89wVv369REZGYmIiAh069YNbm4O17FLRJVEWe/fDvW/2iuvvIJx48YhPDwcnTp1wqefforU1FQ899xztg6NnFRZpnmbmxjl5enw6qvvo1Onafjoo4LNJE+cmAtgLpYuBZYuLXx2fj1aDlSqUxg/qA7Cj81DzzNr0AwF9WY5TZogMSwMNV98EW179QIAtNipR15eHpo3b25IbNq3b1/sUC8RkVOwSL+SBX300UfSoEEDUavV0q5dO9m1a5fZ53J2FZVFaadi3v/63bsif/4pMnLkNwK8Ia1a/Snt2yv7WZY0Hbt2bWWG96RJIosWiYwZs05UCJF+ru6yHhCtSmU4OM/DQzbXqSMdC+0N9cILLxhiys3NNaxATETkaCrFcFV5cbiKyiq/5ya/OPj+ISYR4MqVgv2XVq/eiz17bqB69U64das6ivspc3UFQkOLLpDXtClQs2bBce/PmIGr//43plevjho3bxra/3R3xzpvb3x86xZuF7puaGgoYmJiMHr06Ar+ShARWV9Z799McojMpNFooNUK3N2b4euvDxuGl/JrZm7dKunsm6hf/w56965vtLFkw4aAWl3MKTodsG0bkqZPR+iJEwVjy1Wr4vcmTbAcwLL9+wEA3t7eGDBgACIiInDs2DG4u7uXuVaIiMjeVIqaHCJruXGjcI0MsGlTErTaPwGEIjfXDZGRRc9RqYCQEOPdsadM6YW8vCNwd09HWloJ22oXlpoKrFgBrFwJ/P03mt5r/s3FBWs9PfFBWho6VK2KDgBujhgBNzc3rFy5Ep6e9xckExFVbkxyqNLKywPOnjXujcl/5G/hVKCp4W9qdQ602iNo1UqDESNaGpKaxo2NF8mLiYlBXt4O8/aNyc0FNm0Cli2DbN8O1b0O1msAPgewHMBxvR4NatXCK1evIvTebzLffPNNBX01yFIsOUOPiErGJIecXkZG0SQmKQk4eRLQaos/LzAQUKtP4/TpLRg0qDFefnkAwsKAunU1WLBgK6KiojB8eDSGDzd98yo8PTz/OXDfYlynTgHLl0NWr4bq8mUAgArAzwCWAYgD0LRlS4y4NyOqTZs2nBHlYOxmI1aiSohJDjkFvR5IS0OROpkTJ4CLF4s/z8OjoNC3cPFvkyaAtzcwf/4X934LH2B0XknLo5ta/yb/z6ioKLjm5mJW06bIfP99VD1wAICS2KBOHWDiRPzVsSOi/v1vREZGYkFEBBo1alT+LxDZTOF/+/znxa2RREQVi4XH5FDu3AGSk4v2zCQnA3fvFn9eQIBxIpP/96AgwMWlYmMsbngi9/Bh7Bg9Gh2Tk1HtXnKkB7ANwJmePTF12zbA3b3oBckpPGiGHhEVj7OrzMAkxzGIABcumB5iSk0t/jx3d6UupnASk/93X1/rxW/kzh3gm2+Q+/HHcN+3z9CcCmCNuzvO9+uHrmPGYODAgfyerASUGXrKPmWFNw8mopJxdhU5nOxspSTF1BBTZmbx59WqVXR4qWlTZWaTPexMcP36dexesgSB27ej7bFjQEYG3AHkAfhRrcbJHj3QeMoUvNq3LzQaja3DJSt50EasRFTx7OCWQM5MBLh61fQMprNnlVoaU1xdlTVk7h9eatpUSXLsTWpqKn5Yvx53V6xAt+RkDC38YsOGwKRJuNynD/q1bYuBrq62CpNsxOxCdCKqUExyqELk5gJnzhQdXjpxAii0QG8Rvr5FE5mwMGUV4GIXybMjn3/2GXa89Ra6JSdjHADve+1alQrJzZsj+K23UGXwYMDFBfVsGSjZzP+3d+9hUdX5H8Dfw8CAFzRvICorC5lYZLqYZmZFq2kXCbJ1STI1c9fWflltW1brQFCU3R7LLXvW1PyZm/tb3KHykqihWVqmiY+thqt4AdGUVCAvMzDz+f1xYmRkwJlxbufM+/U889ScOXP4zleZ8/Z7vt/PudRE9KbPici7GHLILadOOb+8tH+/UnfGGZ0OSEhwfokpNlZ5XQ1sNhu+/vprpKSkoENDA7BkCe54+WU8+MvSbwD4KSYGePhhdHniCaQEwZATa7QEntVqdTrJuLUVekTkHQw51IzV6rxIXlmZcn+mlrRr53wFU58+gFqL8VosFpSUlMBkMuHjoiL0/fFHLBo6FB2++w4wm9ENQIPBgPqMDLT5n/9Bl2HDgiq1sUZL4LUWIjmCQ+RbDDkhrK6u5SJ5rS386NWr+eWl5GSgZ8+gOr977Ny5c1ixYgVMJhNWrlyJyNpaTAKwEcBVALBli7LjgAHA1KkIHz8e4VdcEbD2toY1WogolHEJucbZbEBlpfNLTFVVLb8vMlIpiHfx5aWrrgKio/3Xfn+xWq3Q/zIh+MiRI4jv1QsjAUwFcA+Axuo10r49dOPHA1OnAqmpqkl1rNFCRGrGOjku0HLIOXu2eZG8sjLlcfZsy++LjW0eZBqL5Kl1EZCr81AOHjwIk8kEk8mEdu3aYfXq1UoiXLQIJ155Bd2adtyQIUqw+f3vgfbt/fhpvIc1WohIrVgnJwSIAMeOOV/BdOhQy+8LD2+5SF6QXmW5LC3NQ8nLy0NOTg7S0tIwcOBAlJaWKvsDSNfr0XDHHQgvLgZsNnQDlM6ZMAF4+GGgf3//fxAvYo0WIgpFDDlByGy+UCTv4jkzdXUtv69zZ8dRmaZF8kLpbgEtzUPJyckBAJSUlAAAEnU65MbHI7OmBu1raoDPPlMOcPPNyqjN2LHqnTHdBGu0EFGoYsgJEBGgutr5xN/y8paL5IWFOS+Sl5wcnEXyAsFsNiM1NRWpqakwGo148cUXYbFYcPfdd2Pj2rV4/tpr8cC5c+j5n/9cuE9E167ApEnKqE3fvgFtvzexRgsRhTKGHB+rrwcOHHA+8ffkyZbf16GD8xVMSUnKpGBfUmNtldraWqxevRomkwmrVq1C3S9DXuHh4bBYLEiJiMDyxEREREdDt23bhTeOHKmM2txzjzqqD7qJNVqIKJQx5HjJ6dPOLy/t368EHWd0OqB3b+e1Zbp3D9zCHTXVVtm3bx8ee+wxrF+/HhaLxb49Li4OST16IGH7dkzV6XBzfT3w9tvKiz16AJMnA1OmKNfyNIw1WogolDHkeEHfvsrKppa0bdtykby2bf3XTlcFc22V8vJyVFdXY/DgwQCALl26YO3atWhoaMBVV12FzMxMjE9Jwfm5c9Fn61Z0AgAR2HQ6rBTB2fHj8fvFi4PjTp5EAaDGkVoiT/Gb3gs6dlT+27Nn88tLjUXywsIC20Z3NQ06jXNaAhFwRAQ7d+60L/XetWsXbrjhBmz5pSBfp06d8MEHHyC1b18kl5YC8+cDs2dfOEDv3sCUKQibPBmlixbBaDRib3IyRzEoZKlppJboskkIqampEQBSU1Pj1eMeOCBSW+vVQwYNg8EgAMRgMPj153711Vfy+OOPS0JCggCwP/R6vYwYMULMZrOIzSaydavI1Kki7duLKPO5pSEsTHZffbXIZ5+JWK0Ox83Ly5OcnBy/fhaiYJOXlycAJC8vz+lzomDj6fmbIYda1PjF1xh0fPkFeP78eYfn9957rz3YtGnTRjIyMmTx4sXy008/iZw6JfK3v4lcd5092Agg0qePyOzZIseO+aydRFrhz99vosvFkOMChhzX+eNfeqdPn5alS5fKfffdJ+3atZOysjL7a8uXL5cHH3xQTCaTnDlzRhm12bRJ5MEHRdq0uRBsIiNFsrNFNmxQ9iEilwVqpJbIXQw5LmDIcU1LgcYbQaeqqkrmzZsnt99+u0RERDhcipozZ07zN5w4IfLGGyL9+jmO2qSkiLz1lshPP3ncFqJQxpEcUhNPz9+ceEzN+Kq2ysaNG3Hrrbc6bOvXrx8yMzORkZGBQYMGKRttNqCkRJlEbDIBjUvD27YFsrKUujZDhqjm5phEwYZVsClU8Aad5HUigu+++w5FRUXo0aMHHnnkEQDAuXPn0K1bN6SkpNiDTd+m1YWPHgU++AB4/32l7HOj1FQl2Nx/v1IlkYg81lI5iGApE0HkDG/QSQHV0NCATZs2wWQyoaioCBUVFQCAa665xh5y2rRpg8rKSlzR9K6gViuwZo0yavPpp8pzQAkz2dlKuBk40N8fh0izWAWbQglHcuiyPfnkk1i8eDFONrlPRdu2bXHHHXcgIyMD2dnZ0F18aenwYWDhQuXxSyACANx4oxJsfvc7oF07P30CIiIKZhzJIb84deoUiouLMW7cOHtwOXnyJE6ePIkuXbogPT0dGRkZGDlyJNpcfAfv+npgxQpl1Oazz5QpxADQuTO29OmDnddfj2lz5zb7mazCSkREnlBZHV4KhCNHjuDdd9/FyJEjERMTg6ysLHz33Xf215988kmUlJTg2LFjWLhwIdLT0x0Dzr59wLPPAvHxwL33AqtXKwEnLQ34xz+AI0ew7q678Mjf/ob8/HyHn904T0Cv1/vr4xIRkUZwJIecqqysxIcffgiTyYStW7c6vHbNNdfg1KlT9uf9+/dvfgCzGfj3v5VRm5KSC9tjY4FJk4CHHwauvNK+OZjvl0VEROrEOTkEQFkRdfbsWbT7ZR5MSUkJbrvtNgCATqfDDTfcYF8R1adPn5YPtHu3Emz+93+Bxjk6Oh0wapQy12bMGCAiosW3NwYbg8EQsPtlERFRcPH0/M2QE8Lq6+vxxRdf2FdEjR07Fm+99RYAZbVUVlYWRo4cifT0dMTFxbV8oLNngf/7PyXcbN58YXuvXsBDDymP3r1dbldkZCQsFgsMBgPMZrOnH4+IiDSCE4/JJWfOnEFxcTFMJhNWrFjhcNlp/fr19v8PDw9HYWFh6wfbsUMJNkuXArW1yja9Hrj7bmXUZvRo5bkb8vPz7QHHYrEgPz+fIzlEROQRhpwQk5qairKyMvvzbt26IT09HZmZmfjtb3976QPU1gIffaSEm+3bL2xPTFTm2UyaBLQ26tMKVmElIiJvYsjRqIqKChQVFWH9+vUoLCxEeLjyRz1y5EiYzWZkZmYiMzMTN95446VXLokA33yjBJtly5TLU4Ayt+bee5VRm7Q0IMzzxXrOJhk7m4xMRETkKoYcjRAR7NmzByaTCSaTCdubjLJs2rQJaWlpAIBXX30Vb7/9dvPifM6cPAl8+KESbr7//sL25GQl2Dz4INC1q1fazyqsRETkbZx4rAFr1qzBY489hr1799q36XQ6DBs2DJmZmbj//vtbnzjclAjwxRdKsCksVJaCA0BUFDBunBJuhg3jzTGJiMhvOPE4RFgsFmzYsAGxsbG47rrrAABdunTB3r17YTAYMGLECGRkZCA9PR2xsbGuH/j4cWDxYuXmmE3CEq67Tgk22dlA03tOeVlubi70er3TS1KseExEFBzU9l3Niscq8PPPP6OwsBAPPPAAYmJiMGrUKMyZM8f+empqKpYvX44TJ05g5cqVmDp1qmsBx2YDiouV+0T16gU8/bQScNq3V4LN1q3KCqrp030acABAr9fDaDSy4jERURBT3Xe1hJCamhoBIDU1NYFuyiVZrVZZsGCBjBkzRiIjIwWA/REbGyvPPPOM5wevrBTJzxdJSBBRLlApj8GDRebPF6mr894HcUNeXp4AkLy8PKfPiYgo8ALxXe3p+ZshJ4j89NNPDs/79etnDzZJSUny1FNPyZdffikNDQ3uH7y+XuTjj0XGjBEJC7sQbK64QuTRR0V27vTSp7g8jb8sBoOBAYeIKEj5+7va0/M3Jx4HkIjg+++/R1FREUwmE/bu3YsTJ07Yb2753nvv4fjx48jMzERKSoprK6IuduAAsGABsGgRUFV1Yfvw4colqfvuAy6+W3iAseIxEVHw8+d3NSceq4TNZsOWLVvswWb//v3218LCwrBt2zYMHz4cADBt2jTPfojFAnz8sbJCat06ZcwGUJZ7T5yoFO1LTr7cj+ITrHhMRBT8VPNd7YthJW87cOCAPPTQQ5KQkCBRUVGSmJgoRqNRzGazW8cJhstVb7zxhsP8msjISBkzZowsXLhQjh8/fnkH/+EHkaeeEunWzXGuzYgRIv/8p8j58975ED7COTlERMGPc3K8bPXq1TJp0iRZs2aN7N+/Xz7++GOJiYmRP//5z24dx58hp7a2VpYtWyZZWVlSWFho3/7DDz9Ix44dJTs7WwoLC6Xucif5nj0rsmSJyM03OwabuDiR554T2b//Mj+Jf7T0S8KgQ0QUPAL1Xe3p+VsVl6tGjx6N0aNH258nJiairKwM8+bNw+uvvx7Aljk6fvw4PvnkE5hMJqxbtw4WiwWAcrfvsWPHAgD69u2LEydOICIi4vJ+2K5dyuWoJUuA06eVbWFhwJ13KnNt7rwTCL/wxxvstQ1Y8ZiIKPip7btaFSHHmZqaGnTu3LnVfcxms8NkqNrGO2V7mcViwciRI7Fp0yZIk3ncffr0QWZmJu677z6H/T0OOD//DPzzn0q4+eabC9t79wamTAEmT1bq3TjRWNsAcLwHVNN7RgVSawErKK/zEhGFILV9V6sy5Ozfvx9z587FG2+80ep+L7/8Ml544QWft6dxZrmIYNCgQcjIyEBmZib69evX6oool0ZXcnKUu33Pnw/84x9K0AGUUZp77lFGbUaMAC5RgMnZzS6d3RSTiIhIM3xy8cxFOTk5DpNwnT2+/fZbh/ccOXJErrzySpkyZcolj3/+/HmpqamxPyoqKnw2J2fbtm1y+PBht97T2rXNDoCsvPtukQEDHOfaXHmlyOzZIseOedRO1qEhIiK1UWWdnOrqalRXV7e6T0JCAqKiogAAVVVVSEtLw5AhQ/DBBx8gLMy9u1IEW50cwPFy0ay//hWLpk6FbsECZEdEIKK+XtkpMhIYO1YZtbnllsu+OSbr0Lgv2Oc0ERFpmSrr5HTt2hVdu3Z1ad8jR44gLS0NqampWLRokdsBJ1jNmjULbc6cwRGjEXuMRkxufKG+HrjmGiXYTJgAXGL+kau8XdsgVE7+wT6niYiInPDJuJKXNV6iuu2226SyslKOHj1qf7gjGOrk2FmtIuvXi2RliRgM9stRPwMikyeLbN4sYrN59Uf6orZBKC39Zh0fIqLA0HSdnEWLFrU4Z8cdQRFyjh4VKSgQSUpymGuzDZDper1E++ik6cswEkonf85pIvKunJycFn+P8vLyJCcnx78NoqCk6ZDjLQELOQ0NIitXimRkiOj1F8JNdLR8e/31MtAPAcHXXyShdPJv/IwGgyHQTSFSvVAaDSbPMeS4wO8h59AhkZwckfh4xxVSQ4eKLFwoL//1r5r65Q6Fk38ohTkifwml0WDyDEOOC/wSciwWkX//W+SOO0R0ugvBpnNnkRkzRL7/3r6rloZpQ+Hkzy9iIt8Jhe8Q8hxDjgt8GnL27ROZOVOke3fHUZtbbxVZulTk3Dnv/8wgEQonfw6pE/leKIwGk2c0fe+qoHbuHHD33cDnn1/YFhOj3GJhyhSgT5/Atc0PnFVNdlZdWe3Udr8WIrXxdnkLIkClt3UIKm3aAOfPKwX6Ro1S6tqMGQNc7g04VSJUTv5qu18LkZpc/I+lxucAf7/o8gS04rG/+azi8Y4dSrG+3r29d0wiohDQ0j30eG89akqVFY81Y+DAQLeAiEKIliqNh8poMAUGQw4Rkcpo6TYjvBRMvsSQQ0SkMs4m9/PyDlFznJNDRKRSjcGmcUUSAw5plafnb4YcIiIV0+v1sNlsMBgMMJvNDq+pbX4OUUs8PX+H+bBNRETkQ/n5+bDZbABgry3T9DWj0Qi9Xh+o5hEFHOfkEBGp0MWTjI1Go32OTuNzXr6iUMeQQ0SkMi1NMm4adBhwiDgnh4hIdVqqkxMZGQmLxQK9Xo+GhoYAtY7I+zgnh4goROTm5jYLOE3v/WS1Wh3m5xCFKoYcIiKVa3r5ymw2Iy8vD0ajkUGHQh7n5KiQlkq6E9HlcTY/x1mxQKJQxJCjQloq6U5El4f3fiJqGSceq9TF/3pjSXciItIqVjx2gZZCDsCS7kREFBoYclygtZADXFgy6qykOxERkRZwCXkIarpk9OKS7kRERKGOIUeluGSUiIiodVxdpUJcMkpERHRpDDkqxCWjREREl8aJx0RERBTUOPGYiIiIqAmGHCIiItIkhhwiIiLSJIYcIiIi0iSGHCIiItIkhhwiIiLSJIYcIiIi0iSGHCIiItIkhhwiIiLSJIYcIiIi0iSGHCIiItIkhhwiIiLSJIYcIiIi0iSGHCIiItIkhhwiIiLSJIYcIiIi0iSGHCIiItIkhhwiIiLSJIYcIiIi0iTVhRyz2YwBAwZAp9OhtLQ00M0hIiKiIKW6kPP000+jR48egW4GERERBTlVhZzVq1ejuLgYr7/+eqCbQkREREEuPNANcNWPP/6IqVOnoqioCG3btnXpPWazGWaz2f68trbWV80jIiKiIKOKkRwRwaRJkzBt2jQMGjTI5fe9/PLL6Nixo/0RHx/vw1YSERFRMAloyMnNzYVOp2v1sW3bNsydOxe1tbV49tln3Tr+s88+i5qaGvujoqLCR5+EiIiIgo1ORCRQP7y6uhrV1dWt7pOQkICsrCx8+umn0Ol09u1WqxV6vR7Z2dlYvHixSz+vtrYWHTt2RE1NDTp06HBZbSciIiL/8PT8HdCQ46rDhw87zKepqqrCqFGjUFhYiCFDhqBXr14uHYchh4iISH08PX+rYuLxr371K4fn7du3BwAkJSW5HHCIiHJzc6HX6zFr1qxmr+Xn58NqtSI3N9f/DSMin1DFxGMiIm/Q6/UwGo3Iz8932J6fnw+j0Qi9Xh+glhGRL6hiJOdiCQkJUMFVNiIKMo0jOEaj0f68MeDk5eU5HeEhIvVSxZwcb+GcHCICLozcGAwGWCwWBhyiIKfpicfewpBDRI0iIyNhsVhgMBgcioYSUfDx9PzNOTlEFHLy8/PtAcdisTSbo0NE2sCQQ0QhpekcHLPZjLy8PKeTkYlI/VQ58ZiIyBPOJhk7m4xMRNrAkENEIcNqtTqdZNz43Gq1BqJZROQjnHhMREREQY0Tj4mIiIiaYMghIiIiTWLIISIiIk1iyCEiIiJNYsghIiIiTWLIISIiIk1iyCEiIiJNYsghIiIiTWLIISIiIk1iyCEiIiJNYsghIiIiTWLIISIiIk1iyCEiIiJNYsghIgoyubm5yM/Pd/pafn4+cnNz/dsgIpUKD3QD/ElEACi3bCciClb19fV44YUXcP78eTzzzDP27bNnz0ZBQQGee+45fo9RSGn8+954HneVTtx9h4pVVlYiPj4+0M0gIiIiD1RUVKBXr14u7x9SIcdms6GqqgrR0dHQ6XReO25tbS3i4+NRUVGBDh06eO245Ij97D/sa/9gP/sH+9k/fNnPIoK6ujr06NEDYWGuz7QJqctVYWFhbiVAd3Xo0IG/QH7AfvYf9rV/sJ/9g/3sH77q544dO7r9Hk48JiIiIk1iyCEiIiJN0udyLaJX6PV63HrrrQgPD6krgH7HfvYf9rV/sJ/9g/3sH8HWzyE18ZiIiIhCBy9XERERkSYx5BAREZEmMeQQERGRJjHkEBERkSYx5Ljo3Xffxa9//WtERUUhNTUVmzZtanX/5cuX4+qrr0ZkZCSuvvpqmEwmP7VU3dzp5/nz52P48OHo1KkTOnXqhBEjRmDr1q1+bK16ufv3udGyZcug0+mQkZHh4xZqh7t9ffr0aUyfPh1xcXGIiopCv379sGrVKj+1Vr3c7ec5c+agb9++aNOmDeLj4/HEE0/g/PnzfmqtOn3xxRcYM2YMevToAZ1Oh6Kioku+Z+PGjUhNTUVUVBQSExPx3nvv+aGlTQhd0rJlyyQiIkLmz58vu3fvlhkzZki7du3k0KFDTvffvHmz6PV6KSgokD179khBQYGEh4fL119/7eeWq4u7/Tx+/Hh55513ZMeOHbJnzx6ZPHmydOzYUSorVmxgMAAAB1pJREFUK/3ccnVxt58bHTx4UHr27CnDhw+Xe+65x0+tVTd3+9psNsugQYPkzjvvlC+//FIOHjwomzZtktLSUj+3XF3c7ecPP/xQIiMjZenSpXLgwAFZs2aNxMXFyeOPP+7nlqvLqlWr5Pnnn5fly5cLADGZTK3uX15eLm3btpUZM2bI7t27Zf78+RIRESGFhYV+arEIQ44LBg8eLNOmTXPYlpycLDNnznS6/7hx42T06NEO20aNGiVZWVk+a6MWuNvPF2toaJDo6GhZvHixL5qnGZ70c0NDgwwbNkzef/99mThxIkOOi9zt63nz5kliYqJYLBZ/NE8z3O3n6dOny2233eaw7cknn5SbbrrJZ23UGldCztNPPy3JyckO2/74xz/KDTfc4MumOeDlqkuwWCzYvn07br/9doftt99+OzZv3uz0PVu2bGm2/6hRo1rcnzzr54udPXsW9fX16Ny5sy+aqAme9nNeXh66deuGKVOm+LqJmuFJX3/yyScYOnQopk+fjtjYWKSkpKCgoABWq9UfTVYlT/r5pptuwvbt2+2Xt8vLy7Fq1SrcddddPm9vKGnpXLht2zbU19f7pQ3BUZIwiFVXV8NqtSI2NtZhe2xsLI4dO+b0PceOHXNrf/Ksny82c+ZM9OzZEyNGjPBFEzXBk37+6quvsGDBApSWlvqjiZrhSV+Xl5fj888/R3Z2NlatWoX//ve/mD59OhoaGmA0Gv3RbNXxpJ+zsrJw4sQJ3HTTTRARNDQ04JFHHsHMmTP90eSQ0dK5sKGhAdXV1YiLi/N5GxhyXKTT6Ryei0izbZezPyk87bdXX30VH330ETZs2ICoqChfNU8zXO3nuro6PPDAA5g/fz66du3qr+Zpijt/p202G2JiYvD3v/8der0eqampqKqqwmuvvcaQcwnu9POGDRvw0ksv4d1338WQIUOwb98+zJgxA3FxcZg1a5Y/mhsynP25ONvuKww5l9C1a1fo9fpm/yI4fvx4s4TaqHv37m7tT571c6PXX38dBQUFWLduHfr37+/LZqqeu/28f/9+HDx4EGPGjLFvs9lsAIDw8HCUlZUhKSnJt41WKU/+TsfFxSEiIgJ6vd6+rV+/fjh27BgsFgsMBoNP26xGnvTzrFmzMGHCBDz88MMAgGuvvRZnzpzBH/7wBzz//PMIC+NMDm9o6VwYHh6OLl26+KUN/JO8BIPBgNTUVKxdu9Zh+9q1a3HjjTc6fc/QoUOb7V9cXNzi/uRZPwPAa6+9hvz8fHz22WcYNGiQr5upeu72c3JyMnbt2oXS0lL7Iz09HWlpaSgtLUV8fLy/mq46nvydHjZsGPbt22cPkgCwd+9exMXFMeC0wJN+Pnv2bLMgo9frIcpiHJ+1NdS0dC4cNGgQIiIi/NMIv01xVrHG5YkLFiyQ3bt3y+OPPy7t2rWTgwcPiojIhAkTHGbxf/XVV6LX6+WVV16RPXv2yCuvvMIl5C5wt59nz54tBoNBCgsL5ejRo/ZHXV1doD6CKrjbzxfj6irXudvXhw8flvbt28ujjz4qZWVlsmLFComJiZEXX3wxUB9BFdzt55ycHImOjpaPPvpIysvLpbi4WJKSkmTcuHGB+giqUFdXJzt27JAdO3YIAHnzzTdlx44d9qX6M2fOlAkTJtj3b1xC/sQTT8ju3btlwYIFXEIerN555x3p3bu3GAwG+c1vfiMbN260v3bLLbfIxIkTHfb/17/+JX379pWIiAhJTk6W5cuX+7nF6uROP/fu3VsANHvk5OT4v+Eq4+7f56YYctzjbl9v3rxZhgwZIpGRkZKYmCgvvfSSNDQ0+LnV6uNOP9fX10tubq4kJSVJVFSUxMfHy5/+9Cc5depUAFquHiUlJU6/cxv7duLEiXLLLbc4vGfDhg0ycOBAMRgMkpCQIPPmzfNrm3UiHJsjIiIi7eGcHCIiItIkhhwiIiLSJIYcIiIi0iSGHCIiItIkhhwiIiLSJIYcIiIi0iSGHCIiItIkhhwiIiLSJIYcIiIi0iSGHCIiItIkhhwiUq0TJ06ge/fuKCgosG/75ptvYDAYUFxcHMCWEVEw4L2riEjVVq1ahYyMDGzevBnJyckYOHAg7rrrLsyZMyfQTSOiAGPIISLVmz59OtatW4frr78eO3fuxLfffouoqKhAN4uIAowhh4hU79y5c0hJSUFFRQW2bduG/v37B7pJRBQEOCeHiFSvvLwcVVVVsNlsOHToUKCbQ0RBgiM5RKRqFosFgwcPxoABA5CcnIw333wTu3btQmxsbKCbRkQBxpBDRKr2l7/8BYWFhdi5cyfat2+PtLQ0REdHY8WKFYFuGhEFGC9XEZFqbdiwAXPmzMGSJUvQoUMHhIWFYcmSJfjyyy8xb968QDePiAKMIzlERESkSRzJISIiIk1iyCEiIiJNYsghIiIiTWLIISIiIk1iyCEiIiJNYsghIiIiTWLIISIiIk1iyCEiIiJNYsghIiIiTWLIISIiIk1iyCEiIiJNYsghIiIiTfp/CkmJqTTUJpAAAAAASUVORK5CYII=", "text/plain": [ "Figure(PyObject
)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "using PyPlot, LinearAlgebra\n", "\n", "# Model specification: y|x ~ 𝒩(f(x), v(x))\n", "f(x) = 5*x .- 2 \n", "v(x) = 10*exp.(2*x.^2) .- 9.5 # input dependent noise variance\n", "x_test = [0.0, 1.0]\n", "plot(x_test, f(x_test), \"k--\") # plot f(x)\n", "\n", "# Generate N samples (x,y), where x ~ Unif[0,1]\n", "N = 50\n", "x = rand(N)\n", "y = f(x) + sqrt.(v(x)) .* randn(N)\n", "plot(x, y, \"kx\"); xlabel(\"x\"); ylabel(\"y\") # Plot samples\n", "\n", "# Add constant to input so we can estimate both the offset and the slope\n", "_x = [x ones(N)]\n", "_x_test = hcat(x_test, ones(2))\n", "\n", "# LS regression\n", "w_ls = pinv(_x) * y\n", "plot(x_test, _x_test*w_ls, \"b-\") # plot LS solution\n", "\n", "# Weighted LS regression\n", "W = Diagonal(1 ./ v(x)) # weight matrix\n", "w_wls = inv(_x'*W*_x) * _x' * W * y\n", "plot(x_test, _x_test*w_wls, \"r-\") # plot WLS solution\n", "ylim([-5,8]); legend([\"f(x)\", \"D\", \"LS linear regr.\", \"WLS linear regr.\"],loc=2);" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Too Few Training Samples\n", "\n", "- If we have fewer training samples than input dimensions, $\\mathbf{X}^T\\mathbf{X}$ will not be invertible. (Why?)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- As a general recipe, in case of (expected) problems, **go back to full Bayesian!** Do proper model specification, Bayesian inference etc. Let's do this next. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **Model specification**. Let's try a Gaussian prior for $w$ (why is this reasonable?)\n", "\n", "$$\n", "p(w) = \\mathcal{N}(w|0,\\Sigma) = \\mathcal{N}(w|0,\\varepsilon I)\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **Learning**. Let's do Bayesian inference,\n", "\n", "$$\\begin{align*}\n", "\\log p(w|D) &\\propto \\log p(D|w)p(w) \\\\\n", " &\\stackrel{IID}{=} \\log \\sum_n p(y_n|x_n,w) + \\log p(w)\\\\\n", " &= \\log \\sum_n \\mathcal{N}(y_n|\\,w^Tx_n,\\sigma^2) + \\log \\mathcal{N}(w|0,\\varepsilon I)\\\\\n", " &\\propto \\frac{1}{2\\sigma^2}\\left( {y - \\mathbf{X}w } \\right)^T \\left( {y - \\mathbf{X}w } \\right) + \\frac{1}{2 \\epsilon}w^T w\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **Done!** The posterior $p(w|D)$ specifies all we know about $w$ after seeing the data." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Too Few Training Samples, cont'd: the MAP estimate\n", "\n", "- As discussed, for practical purposes, you often want a point estimate for $w$, rather than a posterior distribution." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- For instance, let's take a **Maximum A Posteriori (MAP) estimate**. Set derivative \n", "$$\\nabla_{w} \\log p(w|D) = -\\frac{1}{\\sigma^2}\\mathbf{X}^T(y-\\mathbf{X}w) + \\frac{1}{\\varepsilon} w\n", "$$ \n", "to zero, yielding\n", "$$\n", "\\boxed{ \\hat{w}_{\\text{MAP}} = \\left( \\mathbf{X}^T\\mathbf{X} + \\frac{\\sigma^2}{\\varepsilon} I \\right)^{-1}\\mathbf{X}^T y }\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that, in contrast to $\\mathbf{X}^T\\mathbf{X}$, the matrix $\\left( \\mathbf{X}^T\\mathbf{X} + (\\sigma^2 / \\varepsilon) I \\right)$ is always invertible! (Why?)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note also that $\\hat{w}_{\\text{LS}}$ is retrieved by letting $\\varepsilon \\rightarrow \\infty$. Does that make sense?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Adaptive Linear Regression\n", "\n", "- What if the data arrives one point at a time?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Two standard _adaptive_ linear regression approaches: RLS and LMS. Here we shortly recap the LMS approach." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **Least Mean Squares** (LMS) is gradient-descent on a 'local-in-time' approximation of the square-error cost function." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Define the cost-of-current-sample as \n", "$$\\begin{equation*}\n", "E_n(w) = \\frac{1}{2}(y_n - w^Tx_n)^2\n", "\\end{equation*}$$ \n", "and track the optimum by gradient descent (at each sample index $n$):\n", "$$\\begin{equation*}\n", "w_{n+1} = w_n - \\eta \\, \\left. \\frac{\\partial E_n}{\\partial w} \\right|_{w_n}\n", "\\end{equation*}$$\n", "which leads to the LMS update:\n", "$$\n", "\\boxed{ w_{n+1} = w_n + \\eta \\, (y_n - w_n^T x_n) x_n }\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- (OPTIONAL) This is not a probabilistic modelling derivation. Is there also a Bayesian treatment of LMS? Sure, e.g., have a look at [G. Deng et al., _A model-based approach for the development of LMS algorithms_', ISCAS-05 symposium, 2005](./files/Deng-2005-A-model-based-approach-for-the-development-of-LMS-algorithms.pdf)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "open(\"../../styles/aipstyle.html\") do f\n", " display(\"text/html\", read(f, String))\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Julia 1.1.0", "language": "julia", "name": "julia-1.1" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.1.0" } }, "nbformat": 4, "nbformat_minor": 1 }