{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Regression" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Preliminaries\n", "\n", "- Goal \n", " - Introduction to Bayesian (Linear) Regression\n", "- Materials \n", " - Mandatory\n", " - These lecture notes\n", " - Optional\n", " - Bishop pp. 152-158 \n", " - In this and forthcoming lectures, we will make use of some elementary matrix calculus. The most important formulas are summarized at the bottom of this notebook in an [OPTIONAL SLIDE on matrix calculus](#matrix-calculus). For derivations, see Bishop Appendix C. \n", " " ] }, { "attachments": {}, "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}^M$ 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}^M$ and $y_n \\in \\mathbb{R}$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Assume that we are interested in (a model for) the responses $y_n$ for **given inputs** $x_n$?, I.o.w. we are interested in building a model for the conditional distribution $p(y|x)$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that, since $p(x,y)=p(y|x)\\, p(x)$, building a model $p(y|x)$ is similar to density estimation with the assumption that $x$ is drawn from a uniform distribution." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Bayesian Linear Regression\n", "\n", "- Next, we discuss (1) model specification, (2) Inference and (3) a prediction application for a Bayesian linear regression problem. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "##### 1. Model Specification\n", "\n", "\n", "- In a traditional _regression_ model, we try to 'explain the data' by a purely deterministic function $f(x_n,w)$, plus a purely random term $\\epsilon_n$ for 'unexplained noise':\n", "$$\n", " y_n = f(x_n,w) + \\epsilon_n\n", "$$ " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In a _linear regression_ model, i.e., linear w.r.t. the parameters $w$, we assume that \n", "$$f(x_n,w)= \\sum_{j=0}^{M-1} w_j \\phi_j(x_n) = w^T \\phi(x_n)$$\n", "where $\\phi_j(x)$ are called basis functions.\n", " - For notational simplicity, from now on we will assume $f(x_n,w) = w^T x_n$, with $x_n \\in \\mathbb{R}^M$." ] }, { "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, i.e.\n", "$$\n", "\\epsilon_n \\sim \\mathcal{N}(0,\\beta^{-1}) \\,.\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Hence, given a data set $D=\\{(x_1,y_1),\\dotsc,(x_N,y_N)\\}$, the likelihood for an ordinary linear regression model is \n", "\n", "$$\\begin{align*}\n", "p(y\\,|\\,\\mathbf{X},w,\\beta) &= \\mathcal{N}(y\\,|\\,\\mathbf{X} w,\\beta^{-1} \\mathbf{I}) \\\\\n", " &= \\prod_n \\mathcal{N}(y_n\\,|\\,w^T x_n,\\beta^{-1}) \\tag{B-3.10}\n", "\\end{align*}$$\n", "\n", "where $w = \\left(\\begin{matrix} w_1 \\\\ w_2 \\\\ \\vdots \\\\ w_{M} \\end{matrix} \\right)$, the $(N\\times M)$-dim matrix $\\mathbf{X} = \\left(\\begin{matrix}x_1^T \\\\ x_2^T \\\\ \\vdots \\\\ x_N^T \\end{matrix} \\right) = \\left(\\begin{matrix}x_{11},x_{12},\\dots,x_{1M}\\\\ x_{21},x_{22},\\dots,x_{2M} \\\\ \\vdots \\\\ x_{N1},x_{N2},\\dots,x_{NM} \\end{matrix} \\right) $ and $y = \\left(\\begin{matrix} y_1 \\\\ y_2 \\\\ \\vdots \\\\ y_N \\end{matrix} \\right)$.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- For full Bayesian learning we should also choose a prior $p(w)$:\n", "\n", "$$\\begin{equation}\n", "p(w\\,|\\,\\alpha) = \\mathcal{N}(w\\,|\\,0,\\alpha^{-1}\\mathbf{I}) \\tag{B-3.52}\n", "\\end{equation}$$\n", "\n", "- For simplicity, we will assume that $\\alpha$ and $\\beta$ are fixed and known. \n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "##### 2. Inference\n", "\n", "- We'll do Bayesian inference for the parameters $w$. \n", "\n", "$$\\begin{align*}\n", "p(w|D) &\\propto p(D|w)\\cdot p(w) \\\\\n", " &= \\mathcal{N}(y\\,|\\,\\mathbf{X} w,\\beta^{-1} \\mathbf{I}) \\cdot \\mathcal{N}(w\\,|\\,0,\\alpha^{-1}\\mathbf{I}) \\\\\n", " &\\propto \\exp \\big( -\\frac{\\beta}{2} \\big( {y - \\mathbf{X}w } \\big)^T \\big( {y - \\mathbf{X}w } \\big) - \\frac{\\alpha}{2}w^T w \\big) \\qquad \\text{(B-3.55)} \\\\\n", " &= \\exp\\big( -\\frac{1}{2} w^T\\big(\\underbrace{\\beta \\mathbf{X}^T\\mathbf{X} + \\alpha \\mathbf{I}}_{\\Lambda_N}\\big)w + \\big(\\underbrace{\\beta \\mathbf{X}^Ty}_{\\eta_N}\\big)^T w - \\frac{\\beta}{2}y^T y \\big) \\\\\n", " &\\propto \\mathcal{N}_c\\left(w\\,|\\,\\eta_N,\\Lambda_N \\right)\n", "\\end{align*}$$\n", "\n", "with natural parameters (see the [natural parameterization of Gaussian](https://nbviewer.jupyter.org/github/bertdv/BMLIP/blob/master/lessons/notebooks/The-Gaussian-Distribution.ipynb#natural-parameterization)):\n", "\n", "$$\\begin{align*}\n", "\\eta_N &= \\beta\\mathbf{X}^Ty \\\\\n", "\\Lambda_N &= \\beta \\mathbf{X}^T\\mathbf{X} + \\alpha \\mathbf{I}\n", "\\end{align*}$$\n", "\n", "- Or equivalently (in the moment parameterization of the Gaussian):\n", "$$\\begin{align*}\n", "p(w|D) &= \\mathcal{N}\\left(w\\,|\\,m_N,S_N \\right) \\qquad &&\\text{(B-3.49)} \\\\\n", "m_N &= \\beta S_N \\mathbf{X}^T y \\qquad &&\\text{(B-3.53)}\\\\\n", "S_N &= \\left(\\alpha \\mathbf{I} + \\beta \\mathbf{X}^T \\mathbf{X}\\right)^{-1} \\qquad &&\\text{(B-3.54)}\n", "\\end{align*}$$\n", "\n", "- Note that B-3.53 and B-3.54 combine to\n", "$$\n", "m_N = \\left(\\frac{\\alpha}{\\beta}\\mathbf{I} + \\mathbf{X}^T \\mathbf{X} \\right)^{-1} \\mathbf{X}^T y\\,.\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ " - (Bishop Fig.3.7) Illustration of sequential Bayesian learning for a simple linear model of the form $y(x, w) =\n", "w_0 + w_1 x$. (Bishop Fig.3.7, detailed description at Bishop, pg.154.)\n", "\n", "

" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "##### 3. Application: Predictive Distribution\n", "\n", "- Assume we are interested in the distribution $p(y_\\bullet \\,|\\, x_\\bullet, D)$ for a new input $x_\\bullet$. This can be worked out as (exercise B-3.10)\n", "\n", "$$\\begin{align*}\n", "p(y_\\bullet \\,|\\, x_\\bullet, D) &= \\int p(y_\\bullet \\,|\\, x_\\bullet, w) p(w\\,|\\,D)\\,\\mathrm{d}w \\\\\n", "&= \\int \\mathcal{N}(y_\\bullet \\,|\\, w^T x_\\bullet, \\beta^{-1}) \\mathcal{N}(w\\,|\\,m_N,S_N)\\,\\mathrm{d}w \\\\\n", "&= \\int \\mathcal{N}(y_\\bullet \\,|\\, z, \\beta^{-1}) \\underbrace{\\mathcal{N}(z\\,|\\,x_\\bullet^T m_N,x_\\bullet^T S_N x_\\bullet)\\,\\mathrm{d}z}_{=\\mathcal{N}(w\\,|\\,m_N,S_N)\\,\\mathrm{d}w} \\quad \\text{(sub. }z=x_\\bullet^T w \\text{)} \\\\\n", "&= \\int \\underbrace{\\underbrace{\\mathcal{N}(z \\,|\\, y_\\bullet, \\beta^{-1})}_{\\text{switch }z \\text{ and }y_\\bullet} \\mathcal{N}(z\\,|\\,x_\\bullet^T m_N,x_\\bullet^T S_N x_\\bullet)}_{\\text{Use Gaussian product formula SRG-6}}\\,\\mathrm{d}z \\\\\n", "&= \\int \\mathcal{N}\\left(y_\\bullet\\,|\\, m_N^T x_\\bullet, \\sigma_N^2(x_\\bullet) \\right) \\underbrace{\\mathcal{N}\\left(z\\,|\\, \\cdot, \\cdot\\right)}_{\\text{integrate this out}} \\mathrm{d}z\\\\\n", "&= \\mathcal{N}\\left(y_\\bullet\\,|\\, m_N^T x_\\bullet, \\sigma_N^2(x_\\bullet) \\right)\n", "\\end{align*}$$\n", "with\n", "$$\\begin{align*}\n", "\\sigma_N^2(x_\\bullet) = \\beta^{-1} + x^T_\\bullet S_N x_\\bullet \\tag{B-3.59}\n", "\\end{align*}$$\n", "\n", "- So, the uncertainty $\\sigma_N^2(x_\\bullet)$ about the output $y_\\bullet$ contains both uncertainty about the process ($\\beta^{-1}$) and about the model parameters ($x^T_\\bullet S_N x_\\bullet$). \n", "\n", "- (See the [OPTIONAL SLIDE](#change-of-variable-derivation) below for the step in this derivation where $\\mathcal{N}(w\\,|\\,m_N,S_N)\\,\\mathrm{d}w$ gets replaced $\\mathcal{N}(z\\,|\\,x_\\bullet^T m_N,x_\\bullet^T S_N x_\\bullet)\\,\\mathrm{d}z$.)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Example Predictive Distribution\n", "\n", "- As an example, let's do Bayesian Linear Regression for a synthetic sinusoidal data set and a model with 9 Gaussian basis functions \n", "$$\\begin{align*}\n", "y_n &=\\sum_{m=1}^9 w_m \\phi_m(x_n) + \\epsilon_n \\\\\n", "\\phi_m(x_n) &= \\exp\\left( \\frac{x_n-\\mu_m}{\\sigma^2}\\right) \\\\\n", "\\epsilon_n &\\sim \\mathcal{N}(0,\\beta^{-1})\n", "\\end{align*}$$\n", "\n", "

\n", "\n", "- The predictive distributions for $y$ are shown in the following plots (Bishop, Fig.3.8)\n", "

\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "- And some plots of draws of posteriors for the functions $w^T \\phi(x)$ (Bishop, Fig.3.9) \n", "\n", "

" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Maximum Likelihood Estimation for Linear Regression Model\n", "\n", "- Recall the posterior mean for the weight vector\n", "$$\n", "m_N = \\left(\\frac{\\alpha}{\\beta}\\mathbf{I} + \\mathbf{X}^T \\mathbf{X} \\right)^{-1} \\mathbf{X}^T y\n", "$$\n", "where $\\alpha$ is the prior precision for the weights." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The Maximum Likelihood solution for $w$ is obtained by letting $\\alpha \\rightarrow 0$, which leads to \n", "\n", "$$\\begin{equation*}\n", "\\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 if we have fewer training samples than input dimensions, i.e., if $Ngradient \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 so-called **normal equations** \n", "$\\mathbf{X}^T\\mathbf{X} \\hat w_{\\text{LS}} = \\mathbf{X}^T y$ and consequently\n", "\n", "$$\n", "\\hat w_{\\text{LS}} = (\\mathbf{X}^T \\mathbf{X})^{-1} \\mathbf{X}^T y\n", "$$\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 the (probabilistic) maximum likelihood solution ($\\hat w_{\\text{ML}}$) if the probabilistic model includes the following assumptions:\n", " 1. The observations are independently and identically distributed (**IID**) (this determines how errors are combined), and\n", " 1. The noise signal $\\epsilon_n \\sim \\mathcal{N}(0,\\,\\beta^{-1})$ is **Gaussian** distributed (determines the error metric) \n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- If you use the Least-Squares method, you cannot see (nor modify) these assumptions. The probabilistic method forces you to state all assumptions explicitly! " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Not Identically Distributed Data\n", "\n", "- Let's do an example regarding changing our assumptions. What if we assume that the variance of the measurement error varies with the sampling index, $\\epsilon_n \\sim \\mathcal{N}(0,\\beta_n^{-1})$?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The likelihood is now (using $\\Lambda \\triangleq \\mathrm{diag}(\\beta_n)$ )\n", "$$\n", "p(y\\,|\\,\\mathbf{X},w,\\Lambda) = \\mathcal{N}(y\\,|\\,\\mathbf{X} w,\\Lambda^{-1} ) \\,.\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Combining this likelihood with the prior $p(w) = \\mathcal{N}(w\\,|\\,0,\\alpha^{-1}\\mathbf{I})$ leads to a posterior\n", "\n", "$$\\begin{align*}\n", "p(w|D) &\\propto p(D|w)\\cdot p(w) \\\\\n", " &= \\mathcal{N}(y\\,|\\,\\mathbf{X} w,\\Lambda^{-1} \\mathbf{I}) \\cdot \\mathcal{N}(w\\,|\\,0,\\alpha^{-1}\\mathbf{I}) \\\\\n", " &\\propto \\exp \\left\\{ \\frac{1}{2} \\left( {y - \\mathbf{X}w } \\right)^T \\Lambda \\left( {y - \\mathbf{X}w } \\right) + \\frac{\\alpha}{2}w^T w \\right\\} \\\\\n", " &\\propto \\mathcal{N}\\left(w\\,|\\,m_N,S_N \\right) \n", "\\end{align*}$$\n", "with\n", "$$\\begin{align*}\n", "m_N &= S_N \\mathbf{X}^T \\Lambda y \\\\\n", "S_N &= \\left(\\alpha \\mathbf{I} + \\mathbf{X}^T \\Lambda \\mathbf{X}\\right)^{-1} \n", "\\end{align*}$$\n", "\n", "- And maximum likelihood solution \n", "$$\n", "\\hat{w}_{\\text{ML}} = \\left. m_N\\right|_{\\alpha \\rightarrow 0} = \\left(\\mathbf{X}^T \\Lambda \\mathbf{X}\\right)^{-1} \\mathbf{X}^T \\Lambda y\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- This maximum likelihood solution 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: Least Squares vs Weighted Least Squares\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": 1, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAG2CAYAAACH2XdzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABkq0lEQVR4nO3deXhTZdoG8DtNm5YW2sraspcdQaAsBZFd9q2Nihsq4IqDA26AwlBKwyLKIMqofKKCM+qgo6bsiwhFVtl3FNnLWsDSlhaaNnm+P45JG5rSpE2z9f5dV66ZnHNy8uYQe+685znvqxIRAREREZGX83N3A4iIiIicgaGGiIiIfAJDDREREfkEhhoiIiLyCQw1RERE5BMYaoiIiMgnMNQQERGRT2CoISIiIp/AUENEREQ+gaGGiIiIfILXhJq8vDz84x//QFRUFCpUqIAGDRogMTERJpPJ3U0jIiIiD+Dv7gbYa/bs2ViwYAG+/PJLtGjRArt378aoUaMQFhaGcePGubt5RERE5GZeE2q2b9+O2NhYDBo0CABQv359/Pe//8Xu3bvd3DIiIiLyBF4Tarp06YIFCxbg+PHjaNKkCQ4cOIAtW7Zg3rx5Rb4mJycHOTk5lucmkwl//vknqlSpApVK5YJWExERUWmJCDIzM1GzZk34+d2lcka8hMlkkrfeektUKpX4+/uLSqWSmTNn3vU1U6dOFQB88MEHH3zwwYcPPFJSUu563leJiMALLFmyBOPHj8d7772HFi1aYP/+/Xj11Vcxd+5cjBgxwuZr7uypSU9PR926dZGSkoLQ0FBXNZ2IyKlmz56NmTNnYtKkSZg4cWKxy4m8XUZGBurUqYMbN24gLCysyO285vLT+PHj8dZbb+Hxxx8HANx33304e/YsZs2aVWSoCQwMRGBgYKHloaGhDDVE5LUCAgKQmJiIKVOmWC2fMWMGgoKCYDQa+TeOfFJxpSNeE2qys7MLXUdTq9W8pZuIyp2EhIQi190ZdIjKE68JNUOGDMGMGTNQt25dtGjRAvv27cPcuXPx7LPPurtpRERE5AG8pqYmMzMTU6ZMgV6vR2pqKmrWrIknnngC8fHx0Gg0du0jIyMDYWFhSE9PZ9csERGRl7D3/O01ocYZ7D0oRqMRubm5LmwZlWcBAQFQq9XubgYRkcey9/ztNZefXEFEcPnyZdy4ccPdTaFyJjw8HBERERw/iYioFBhqCjAHmurVqyM4OJgnGCpzIoLs7GykpqYCACIjI93cIiIi78VQ8xej0WgJNFWqVHF3c6gcqVChAgAgNTUV1atX56UoIqIS8ppZusuauYYmODjYzS2h8sj8vWMtFxFRyTHU3IGXnMgd+L0jIio9hhoiIiLyCQw1PkBE8OKLL6Jy5cpQqVTYv38/rl+/jurVq+PMmTN27SMnJwd169bFnj17yraxREREZYShxgesWbMGixcvxooVK3Dp0iW0bNkSs2bNwpAhQ1C/fn279hEYGIg333yTk+AREZHXYqjxASdPnkRkZCQ6d+6MiIgI5Obm4vPPP8fzzz/v0H6GDx+OzZs349ixY2XUUiIiorLDUOPlRo4cib///e84d+4cVCoV6tevj9WrV8Pf3x/333+/ZbvExETUrFkT169ftywbOnQounXrZpkUtEqVKujcuTP++9//uvxzEBERlRbHqSlGVlZWkevUajWCgoLs2tbPz88yHsndtg0JCXGofR988AEaNmyITz/9FLt27YJarcaMGTPQvn17q+0mT56MNWvW4Pnnn4der8eCBQvwyy+/4MCBA1azn8fExGDz5s0OtYGIiMgTMNQUo2LFikWuGzhwIFauXGl5Xr16dWRnZ9vctnv37khOTrY8r1+/Pq5du1ZoO0en4goLC0OlSpWgVqsREREBADhz5gxq1qxptZ1arcZXX32FNm3a4K233sL8+fPx6aefol69elbb1apVy+7iYiIiIk/CUOODbt26ZdWDZNagQQPMmTMHL730Eh577DEMHz680DYVKlQoMpgRERF5MoaaYty8ebPIdXcOZ2+ev8eWgpd4AJRpb0jVqlWRlpZmc90vv/wCtVqNM2fOIC8vD/7+1l+BP//8E9WqVSuzthEREZUVFgoXIyQkpMjHnb0hd9u2YD3N3bZ1hujoaBw9erTQ8m+//RY//vgjkpOTkZKSAp1OV2ibw4cPIzo62intICIiciWGGh/Ur18/HDlyxKq35vz583j55Zcxe/ZsdOnSBYsXL8asWbOwY8cOq9du3rwZffv2dXWTiYiISo2hxgfdd999aN++Pb777jsASvHxyJEjERMTg1deeQUA0KdPH7zyyit46qmnLJfYtm/fjvT0dDzyyCNuazsREVFJqcTR2228WEZGBsLCwpCeno7Q0FCrdbdv38bp06cRFRVls8jW26xatQpvvvkmDh8+XKiepyjDhg1DdHQ0Jk2aVMatozv52vePiMiZ7nb+LoiFwj5q4MCB+OOPP3DhwgXUqVOn2O1zcnLQunVrvPbaay5oHRERkfMx1PiwcePG2b1tYGAg/vGPf5Rha4iIiMoWa2qIiIjIJzDUEBERkU9gqCEiIiKfwFBDREREPoGhhoiIiHwCQw0RERH5BIYaIiIi8gkMNUREROQTGGp8wMiRI6FSqaBSqRAQEIAaNWqgT58++OKLL2AymdzdPCIiIpdgqHGShIQE6HQ6m+t0Oh0SEhLK9P379++PS5cu4cyZM1i9ejV69uyJcePGYfDgwcjLyyvT9yYiIvIEDDVOolarER8fXyjY6HQ6xMfHQ61Wl+n7BwYGIiIiArVq1ULbtm0xadIkLF26FKtXr8bixYvL9L2JiMjzufvHtysw1DjJlClTkJiYaBVszIEmMTERU6ZMcXmbevXqhdatW+PHH390+XsTEblbeTiJO8LdP75dgRNaOpE5uMTHx2P69OkwGAxuCzRmzZo1w8GDB932/kRE7mI+iQOw+jtc8AdneVLwHGV+7u4f387GUONkU6ZMsQQajUbj9i+JiEClUrm1DURE7lAeTuKO8sQf387Ey09OptPpLIHGYDAU2fXpKseOHUNUVJRb20BE5C4FSwMCAwPLdaAxmzJliuUc5Qk/vp2JocaJCv4CyMnJKVRj42obNmzAoUOH8PDDD7vl/YmIPIEvn8RLwtN+fDsTQ42T2OrStFU8XFZycnJw+fJlXLhwAXv37sXMmTMRGxuLwYMH45lnninT9yYi8mS+fBJ3lKf9+HY21tQ4idFotNmlaX5uNBrL9P3XrFmDyMhI+Pv745577kHr1q3x4YcfYsSIEfDzY3YlovLpzh+c5ucAyl2PTVE/vgH4zDFhqHGSu90aWNZfksWLF3MsGiKiO5SHk7gj3P3j2xUYaoiIyCeVh5O4I9z549tVGGqIiMgnlYeTOFljsQURERH5BIYaIiIi8gkMNUREROQTGGqIiIjIJzDUEBERkU/wqlBz4cIFPPXUU6hSpQqCg4PRpk0b7Nmzx93NIiIiIg/gNbd0p6Wl4YEHHkDPnj2xevVqVK9eHSdPnkR4eLi7m0ZEREQewGt6ambPno06depg0aJFiImJQf369fHggw+iYcOG7m6az6lfvz7mzZtnea5SqZCUlOS29hAREdnDa0LNsmXL0L59ewwbNgzVq1dHdHQ0Fi5ceNfX5OTkICMjw+rha0aOHIm4uLgi1+/btw+DBw9G9erVERQUhPr16+Oxxx7DtWvX7H6PS5cuYcCAAU5oLRERUdnxmlBz6tQpfPLJJ2jcuDHWrl2L0aNHY+zYsfj3v/9d5GtmzZqFsLAwy6NOnToubLH7paamonfv3qhatSrWrl2LY8eO4YsvvkBkZCSys7Pt3k9ERAQCAwPLsKXFExHk5eXZtW1ubq7b20BERK7nNaHGZDKhbdu2mDlzJqKjo/HSSy/hhRdewCeffFLka95++22kp6dbHikpKS5ssftt27YNGRkZ+OyzzxAdHY2oqCj06tUL8+bNQ926de3eT8HLT2fOnIFKpcKPP/6Inj17Ijg4GK1bt8b27dsLvXe3bt1QoUIF1KlTB2PHjkVWVpZl/VdffYX27dujUqVKiIiIwJNPPonU1FTL+uTkZKhUKqxduxbt27dHYGAgNm/eXKht5vZ899136NGjB4KCgvDVV18BABYtWoTmzZsjKCgIzZo1w8cff1yojW3atEFQUBDat2+PpKQkqFQq7N+/36E2EBGRZ/CaUBMZGYl7773Xalnz5s1x7ty5Il8TGBiI0NBQq4e9RICsLPc8REp8mKxEREQgLy8Per0e4qyd/mXy5Ml48803sX//fjRp0gRPPPGEpRfj0KFD6NevHx566CEcPHgQ3377LbZs2YJXXnnF8nqDwQCdTocDBw4gKSkJp0+fxsiRIwu9z4QJEzBr1iwcO3YMrVq1KrI9EydOxNixY3Hs2DH069cPCxcuxOTJkzFjxgwcO3YMM2fOxJQpU/Dll18CADIzMzFkyBDcd9992Lt3L3Q6HSZOnGhz3/a2gYiI3Ey8xBNPPCFdunSxWvbqq6/K/fffb/c+0tPTBYCkp6cXWnfr1i05evSo3Lp1S0REbt4UUeKF6x83b9p/XEaMGCGxsbFFrp80aZL4+/tL5cqVpX///vLuu+/K5cuX77rPevXqyfvvv295DkD0er2IiJw+fVoAyGeffWZZf+TIEQEgx44dExGRp59+Wl588UWrfW7evFn8/Pwsx/dOO3fuFACSmZkpIiIbN24UAJKUlHTXtprbM2/ePKvlderUkW+++cZqmU6ns3xfPvnkE6lSpYpVexYuXCgAZN++fQ61wRnu/P4REVG+u52/C/KanprXXnsNO3bswMyZM3HixAl88803+PTTTzFmzBh3N82jzZgxA5cvX8aCBQtw7733YsGCBWjWrBkOHTpUqv0W7LGIjIwEAMvloz179mDx4sWoWLGi5dGvXz+YTCacPn0agFLAHBsbi3r16qFSpUro0aMHABTqeWvfvr1d7Sm43dWrV5GSkoLnnnvOqg3Tp0/HyZMnAQC///47WrVqhaCgIMvrYmJiit03ERF5Lq8Zp6ZDhw7Q6/V4++23kZiYiKioKMybNw/Dhw8vk/cLDgZu3iyTXdv13s5UpUoVDBs2DMOGDcOsWbMQHR2NOXPmWC7FlERAQIDl/6tUKgBK3ZP5f1966SWMHTu20Ovq1q2LrKws9O3bF3379sVXX32FatWq4dy5c+jXrx8MBoPV9iEhIXa1p+B25nYsXLgQHTt2tNpOrVYDUIp+ze02kyIu0dnbBiIici+vCTUAMHjwYAwePNgl76VSAb54LtNoNGjYsKFV0a6ztW3bFkeOHEGjRo1srj906BCuXbuGd955x3JH2u7du532/jVq1ECtWrVw6tSpIkNvs2bN8PXXXyMnJ8dyZ5cz20BERK7nVaGGbEtPT7fcsWNWuXJlHDx4EEuWLMHjjz+OJk2aQESwfPlyrFq1CosWLSqz9kycOBGdOnXCmDFj8MILLyAkJATHjh3DTz/9hPnz56Nu3brQaDSYP38+Ro8ejcOHD0On0zm1DQkJCRg7dixCQ0MxYMAA5OTkYPfu3UhLS8Prr7+OJ598EpMnT8aLL76It956C+fOncOcOXMAoFAPTkH/+te/oNfr8fPPPzu1vUREVHoMNT4gOTkZ0dHRVstGjBiB+Ph4BAcH44033kBKSgoCAwPRuHFjfPbZZ3j66afLrD2tWrXCpk2bMHnyZHTt2hUigoYNG+Kxxx4DAFSrVg2LFy/GpEmT8OGHH6Jt27aYM2cOhg4d6rQ2PP/88wgODsZ7772HCRMmICQkBPfddx9effVVAEBoaCiWL1+Ol19+GW3atMF9992H+Ph4PPnkk1Z1Nne6du2apS6HiIg8i0qKKiTwQRkZGQgLC0N6enqh27tv376N06dPIyoq6q4nNfJdX3/9NUaNGoX09HRUqFDBpe/N7x8RUdHudv4uiD01VG79+9//RoMGDVCrVi0cOHAAEydOxKOPPuryQENERM7BUEPl1uXLlxEfH4/Lly8jMjISw4YNw4wZM9zdLCIiKiGGGiq3JkyYgAkTJri7GURE5CReM/geERER0d0w1BAREZFPYKghIiIin8BQQ0RERD6BoYaIiIh8AkMNERER+QSGGnKKM2fOQKVSWeagSk5Ohkqlwo0bN9zaLiIiKj8YarzYggULUKlSJeTl5VmW3bx5EwEBAejatavVtps3b4ZKpcLx48cBAPXr18e8efOK3PcPP/yAjh07IiwsDJUqVUKLFi3wxhtv2N22zp0749KlSwgLC3PsQxEREZUQQ40X69mzJ27evIndu3dblm3evBkRERHYtWsXsrOzLcuTk5NRs2ZNNGnSpNj9rl+/Ho8//jgeeeQR7Ny5E3v27MGMGTNgMBjsbptGo0FERMRdZ7x2hdzcXLu2MxqNMJlMbm0DERGVDkONF2vatClq1qyJ5ORky7Lk5GTExsaiYcOG2LZtm9Xynj172rXfFStWoEuXLhg/fjyaNm2KJk2aIC4uDvPnz7e7bXdeflq8eDHCw8Oxdu1aNG/eHBUrVkT//v1x6dIlq9ctWrQIzZs3R1BQEJo1a4aPP/7Yav3EiRPRpEkTBAcHo0GDBpgyZYpVaEhISECbNm3wxRdfoEGDBggMDIStOVvN7VmxYgXuvfdeBAYG4uzZszAYDJgwYQJq1aqFkJAQdOzY0er4AsDChQtRp04dBAcHQ6vVYu7cuQgPD3e4DUTkWRISEqDT6Wyu0+l0SEhIcG2DyGEMNUURAbKy3PNw4ATYo0cPbNy40fJ848aN6NGjB7p3725ZbjAYsH37drtDTUREBI4cOYLDhw87dsyKkZ2djTlz5uA///kPfvnlF5w7dw5vvvmmZf3ChQsxefJkzJgxA8eOHcPMmTMxZcoUfPnll5ZtKlWqhMWLF+Po0aP44IMPsHDhQrz//vtW73PixAl89913+OGHHyw1PkW1Z9asWfjss89w5MgRVK9eHaNGjcLWrVuxZMkSHDx4EMOGDUP//v3xxx9/AAC2bt2K0aNHY9y4cdi/fz/69Oljc74oe9tARJ5DrVYjPj6+ULDR6XSIj4+HWq12U8vIblKOpKenCwBJT08vtO7WrVty9OhRuXXrlrLg5k0RJV64/nHzpt2f6dNPP5WQkBDJzc2VjIwM8ff3lytXrsiSJUukc+fOIiKyadMmASAnT560vK5evXry/vvv29znzZs3ZeDAgQJA6tWrJ4899ph8/vnncvv27SLbcfr0aQEg+/btExGRjRs3CgBJS0sTEZFFixYJADlx4oTlNR999JHUqFHD8rxOnTryzTffWO1Xp9PJ/fffX+T7vvvuu9KuXTvL86lTp0pAQICkpqYW+ZqC7dm/f79l2YkTJ0SlUsmFCxestn3wwQfl7bffFhGRxx57TAYNGmS1fvjw4RIWFuZwGwoq9P0jIrdITEwUAJKYmGjzObnH3c7fBXFCSy/Xs2dPZGVlYdeuXUhLS0OTJk1QvXp1dO/eHU8//TSysrKQnJyMunXrokGDBnbtMyQkBCtXrsTJkyexceNG7NixA2+88QY++OADbN++HcHBwSVqa3BwMBo2bGh5HhkZidTUVADA1atXkZKSgueeew4vvPCCZZu8vDyrYuPvv/8e8+bNw4kTJ3Dz5k3k5eUhNDTU6n3q1auHatWqFdsejUaDVq1aWZ7v3bsXIlKo7ignJwdVqlQBAPz+++/QarVW62NiYrBixYoStYGIPMuUKVMAAPHx8Zg+fToMBgMSExMty8mzMdQUJTgYuHnTfe9tp0aNGqF27drYuHEj0tLS0L17dwDKJaSoqChs3boVGzduRK9evRxuRsOGDdGwYUM8//zzmDx5Mpo0aYJvv/0Wo0aNcnhfABAQEGD1XKVSWWpNzEW6CxcuRMeOHa22M3f57tixA48//jimTZuGfv36ISwsDEuWLME///lPq+1DQkLsak+FChWsCplNJhPUajX27NlTqJu5YsWKAAARKVT8bP4MJWkDEXmeKVOmWAKNRqNhoPEiDDVFUakALzkx9ezZE8nJyUhLS8P48eMty7t37461a9dix44dJQ4iZvXr10dwcDCysrJK21ybatSogVq1auHUqVMYPny4zW22bt2KevXqYfLkyZZlZ8+edVoboqOjYTQakZqaWuiWeLNmzZph586dVssK3n1GRN5Pp9NZAo3BYIBOp2Ow8RIMNT6gZ8+eGDNmDHJzcy09NYASal5++WXcvn3bZpHwhQsXChWx1q1bFx9++CGys7MxcOBA1KtXDzdu3MCHH36I3Nxc9OnTp8w+R0JCAsaOHYvQ0FAMGDAAOTk52L17N9LS0vD666+jUaNGOHfuHJYsWYIOHTpg5cqV0Ov1Tnv/Jk2aYPjw4XjmmWfwz3/+E9HR0bh27Ro2bNiA++67DwMHDsTf//53dOvWDXPnzsWQIUOwYcMGrF69uthb15955hnUqlULs2bNclp7icj5zEXB5ktO5ucAGGy8AO9+8gE9e/bErVu30KhRI9SoUcOyvHv37sjMzETDhg1Rp06dQq+bM2cOoqOjrR7Lli1D9+7dcerUKTzzzDNo1qwZBgwYgMuXL2PdunVo2rRpmX2O559/Hp999hkWL16M++67D927d8fixYsRFRUFAIiNjcVrr72GV155BW3atMG2bduc/kdm0aJFeOaZZ/DGG2+gadOmGDp0KH799VfL8XvggQewYMECzJ07F61bt8aaNWvw2muvISgo6K77PXfuXKHb14nIs9wZaAAlyCQmJtq8K4o8j0psFQT4qIyMDISFhSE9Pb1Qcent27dx+vRpREVFFXuCIirohRdewG+//YbNmzeXeB/8/hG5X0JCAtRqtc0fSzqdDkajkWPVuMndzt8F8fITkYPmzJmDPn36ICQkBKtXr8aXX35ZaJBAIvI+dwssvPTkHRhqiBy0c+dOvPvuu8jMzESDBg3w4Ycf4vnnn3d3s4iIyj2GGiIHfffdd+5uAhER2cBCYSIiIvIJDDV3KEd10+RB+L0jIio9hpq/mEe7zc7OdnNLqDwyf+/uHHWZiIjsx5qav6jVaoSHh1vmIgoODi52QDWi0hIRZGdnIzU1FeHh4ZwFmIioFBhqCoiIiAAAS7AhcpXw8HDL94+IiEqGoaYAlUqFyMhIVK9eHbm5ue5uDpUTAQEB7KEhInIChhob1Go1TzJERERehoXCRERE5BMYaoiIiMgnMNQQERGRT2CoISIiIp/AUENEREQ+gaGGiIiIfAJDDREREfkEhhoiIiLyCQw1RERE5BMYaoiIiMgnMNQQERGRT2CoISIiIp/gtaFm1qxZUKlUePXVV93dFCIiIvIAXhlqdu3ahU8//RStWrVyd1OIiIjIQ3hdqLl58yaGDx+OhQsX4p577nF3c4iIiMhDeF2oGTNmDAYNGoTevXsXu21OTg4yMjKsHkREROSb/N3dAEcsWbIEe/fuxa5du+zaftasWZg2bVoZt4qIiIg8gdf01KSkpGDcuHH46quvEBQUZNdr3n77baSnp1seKSkpZdxKIiIicheViIi7G2GPpKQkaLVaqNVqyzKj0QiVSgU/Pz/k5ORYrbMlIyMDYWFhSE9PR2hoaFk3mYiIiJzA3vO311x+evDBB3Ho0CGrZaNGjUKzZs0wceLEYgMNERER+TavCTWVKlVCy5YtrZaFhISgSpUqhZYTERFR+eM1NTVERKRISEiATqezuU6n0yEhIcG1DSKn4b9t6Xh1qElOTsa8efPc3QwiIpdSq9WIj48vdPLT6XSIj4/n5Xgvxn/b0vGay09ERKSYMmUKACA+Pt7y3HzSS0xMtKwn78N/29LxmrufnIF3PxGRLzGf7DQaDQwGA096PoT/ttbsPX8z1BARebHAwEAYDAZoNBrk5OS4uzmllpCQALVabfMErtPpYDQay01dia/925aGvedvr66pISIqz3Q6neWkZzAYiiww9SasKVH44r+tS0g5kp6eLgAkPT3d3U0hIiqVxMREASCJiYk2n3szX/5s9ijvn98We8/fDDVERF6mqJOcL538zJ9Fo9H4zGeyR3n4ty0Je8/fvPuJiMjLGI1Gm4Wj5udGo9EdzXKqKVOmYPr06ZZLMOWlSLY8/NuWJRYKExGRx+HdP1QQC4WJiMgrFRyXJScnB4mJiTaLh4nuxMtPRETkMWwNNGdrQDoiWxhqiIjIY7CmhEqDNTVERETk0VhTQ0REROUKQw0RERH5BIYaIiIi8gkMNUREROQTGGqIiIjIJzDUEBERkU9gqCEiIiKfwFBDREREPoGhhoiIiErMYDDgt99+c3czADDUEBERkYNu3ryJ77//Hk899RSqV6+Onj17wmQyubtZnPuJiIiIinft2jUsW7YMSUlJWLduHXJycizrgoKCcPbsWURFRbmxhQw1REREZIfZs2djzpw5lucNGzaEVquFVqtFp06d4Ofn/os/DDVEREQEABARHD58GHq9Hnq9Hu+88w769esHAIiLi8PPP/+MuLg4aLVatGzZEiqVys0ttsZQQ0REVI6ZTCZs374der0eSUlJOHnypGWdXq+3hJoHHngAe/fudVcz7cJQQ0REVE5dvnwZbdq0wZUrVyzLAgMD0bdvX8TFxWHIkCFubJ3jGGqIqFQSEhKgVqsxZcqUQut0Oh2MRiMSEhJc3zAispKRkYHVq1fj6tWreOWVVwAANWrUQMWKFXH79m0MHjwYcXFx6N+/PypWrOjm1pYMQw0RlYparUZ8fDwAWAUbnU6H+Ph4JCYmuqtpROXelStXLHcsrV+/HgaDAaGhoXjxxReh0WigUqmwZs0a1K1bFxqNpsTvc/48sGwZsGIFsGQJEBrqxA/hAIYaIioVc5ApGGwKBhpbPThEVLb++9//4uOPP8bWrVshIpbljRs3hlarRXZ2tiXENGrUyOH9iwCHDwNLlyqP3bvz161ZAzz6aKk/Qokw1BBRqRUMNtOnT4fBYGCgIXIREcGBAwfQuHFjhISEAAD++OMPbNmyBQDQvn17yx1LzZs3L/EdS3l5wNatQFKSEmROn85fp1IBnToBsbFAx46l/UQlp5KCEc7HZWRkICwsDOnp6Qh1V98YkQ8LDAyEwWCARqOxGpiLiJzLaDRi69atSEpKgl6vx5kzZ/Ddd99h2LBhAJRQs2bNGsTFxaFOnTolfp+sLGDtWiXErFgB/Pln/rrAQKBPHyXIDBkC1KhR2k9VNHvP3+ypISKn0Ol0lkBjMBig0+nYU0PkRDk5OVi/fj30ej2WLVuGq1evWtYFBQUhJSXF8rxx48Zo3Lhxid7nyhVg+XKlR2b9eqDg75PKlYHBg5Ug07cv4Gn1xAw1RFRqd9bQmJ8DYLAhKgURsVwuunTpEgYPHmxZFx4ejiFDhkCr1aJv376WS08l8dtv+fUxO3YoNTNmUVFAXJwSZB54APD34OTgwU0jIm9gqyjYVvEwEdnn0qVLWLZsGfR6PUJCQvDDDz8AAOrXr49Bgwahfv360Gq16NatGwICAkr0HkYj8OuvSohJSgKOH7de3769EmJiY4GWLZWaGW/AUENEpWI0Gm0WBZufG41GdzSLyKucOHHCMjXBjh07LHcsBQYGIisry9ILs2LFihK/x61byuWkpUuVy0upqfnrAgKAnj2VEDN0KFC7dqk+jtuwUJiIiMiNRo0ahcWLF1sti4mJgVarRVxcHJo1a1bifV+7BqxcqQSZtWuB7Oz8dWFhwMCBSpDp31957qlYKExERORB8vLysGXLFiQlJWHy5MmoVq0aAKB169bw9/dHjx49oNVqMXToUNQuRVfJyZP59TFbtgAmU/662rXzLyt17w6UYrw9j8SeGiIiojJy69Yt/PTTT9Dr9Vi+fDmuX78OAPj888/x7LPPAgDS09NhMplwzz33WF7nyPQjJhOwZ09+kDl82Hr7Vq3yC32jo72nPqYg9tQQERG5ycmTJzFhwgSsWbMG2QWu+VSuXBlDhw7Fvffea1kWZuO6T3HTj8THz8DatUqR77JlwMWLBV8LdOuWXx8TFeX8z+epGGqIiIhK6eLFi7h27RpatWoFAAgNDUVSUhJMJhPq1KljGdG3a9eu8LfjnmhbdxBOmvQuZs06hpYtD+P991ug4LRqISFKXUxcnFInU7my0z+iV+DlJyIiohL4/fffodfrkZSUhF9//RXdu3dHcnKyZf2nn36Kdu3aoW3btiWemuCNNz7A3Ll/QKXSQqQbgPxbuCMilJ6Y2FigVy8gKKiUH8iD8fITEfkcR+oMiMrCnj178MMPPyApKQnHjh2zWmc0GpGXl2fpiXnxxRcd3r8IcOBAfn3Mvn3jLMsBoHnz/ELfmBjAz690n8fXMNQQkdcors4gsWB/PJETGI1GqNVqy/OpU6di5cqVAAB/f3/06tULWq0WsbGxiIyMLNF75OYCmzfnB5mzZ/PXqVQCkS1Qq1fCaPwBTzzxDAezvAuGGiLyGrbqDGyNaExUGtnZ2Vi3bh30ej1WrFiBvXv3ol69egCAJ554AhUqVIBWq8XAgQMRHh5eovfIzISl0HflSuDGjfx1FSooE0X6+S1DUtLzSEz8O6ZMeQc6XQhH6S6OlCPp6ekCQNLT093dFCIqhcTERAEgGo1GAEhiYqK7m0Re7vr16/Lll19KXFycVKhQQQBYHvPnz3fKe1y8KLJggciAASIajYhyUUl5VK0qMmqUSFKSSFZW/nf8zu92Uct9nb3nb68pFJ41axZ+/PFH/Pbbb6hQoQI6d+6M2bNno2nTpnbvg4XCRL4jMDDQMit4TsFphIkctHHjRvTp08dqSo969epBq9VCq9Wic+fOdt2xdCcR4NgxpTdm6VJg507r9Y0a5dfHdO6s3IptVpL6MV+uOfO5QuFNmzZhzJgx6NChA/Ly8jB58mT07dsXR48eLdXMpETkfXQ6nSXQGAwG6HQ6dseTXY4dOwa9Xo+aNWti5MiRAID27dtDrVbj3nvvtQSZ1q1bl+iOJaMR2LYtvz7mxAnr9TEx+QPhNW9e9EB4dwsfRX3XWXMG7738lJqaKgBk06ZNdr+Gl5+IvN+d3e/ltTue7GM0GmXHjh0yceJEadq0qeWSUnR0tNV2ly9fLvF7ZGUpl41GjRKpVs36spJGo1xuWrBA5MKF0n6a4vnqfx/2nr+9NtT88ccfAkAOHTpU5Da3b9+W9PR0yyMlJYWhhsiLsc6AHDFx4kSpWbOmVX2MRqORAQMGyMKFC8VkMpV436mpIp9/LjJ0qEiFCtZBJjxcZPhwkf/9TyQjw4kfyE6+WHPmczU1BYkIYmNjkZaWhs2bNxe5XUJCAqZNm1ZoOWtqiLyTL9cMUOlkZWUhOTkZAwcOtFw2euKJJ7BkyRJUqlQJAwcOhFarxYABA0r89/+PP/IvK23dmj92DADUq5dfH9O1KxAQUPR+XMHXas7sranxylAzZswYrFy5Elu2bLnrTKY5OTlW/5gZGRmoU6cOQw2RF2GQoaJcu3YNy5cvR1JSEtatW4fbt2/j0KFDaNmyJQBg165duHr1Kh588EEEBgY6vH+TCdi1SwkxSUlK0W9B0dH5QaZ1a8+ZKNJcQ2OuOfOF4Q7svtGn7DuNnOuVV16R2rVry6lTpxx+LWtqiLwPLzlRQZcuXZJ58+ZJjx49xM/Pz+rSUoMGDWTdunWl2v+tWyIrV4q8+KJIZKT1ZSV/f5HevUXmzxc5e9ZJH8jJyrKmZurUqUXuJzExUaZOnVrq9yiKz9XUmEwmGTNmjNSsWVOOHz9eon0w1BB5J18tfqTimUwmuX37tuX5qlWrrIJMmzZtZNq0aXLgwIES18j8+afIf/4j8sgjIhUrWgeZSpVEHn1U5OuvRdLSnPShykhZ/wBw5w8Mnws1L7/8soSFhUlycrJcunTJ8sjOzrZ7Hww1RN7LF4sfyTaj0Shbt26V8ePHS6NGjWT8+PGWdbdv35Z+/frJ3LlzS9Rjb3b6tMi8eSI9e4qo1dZBpmZNkZdfFlmzRqRAnvJ4ruhJcdcPDJ8rFC5qvIBFixZZxhooDgffI/Juvlb8SPkMBgM2bNiApKQkLF26FJcvX7asa9GiBQ4fPlyq/YsA+/blF/oeOGC9vmXL/PqYdu04UeTduKNmx6cLhUuKoYbIe/li8aO3KauibRFBs2bNcPz4ccuy0NBQDBo0CFqtFv3790elSpUc3q/BAGzapISYZcuAlJT8dX5+QJcuykB4Q4cCDRs6vPtyzdU/MOw9fzOLEpHHKzgiak5ODhITExEfHw+dTufuppUr5hFr7zzu5n+fgrNZFyU1NRWfffYZhg8fDpPJBEDpie/evTsiIiLw0ksvYc2aNbh69Sq++eYbDBs2zKFAk5EBfPst8OSTQPXqQN++wEcfKYEmOBjQaoHFi4ErV5TA89prDDSOsjWit8co04tgHoY1NUTeh3c/eZaS1FScOnVK5s6dK127drW6Y2nbtm2WbdLT08VoNJaoTefPi3z8sUjfviIBAdb1MdWrizz3nMiyZSIOlGBSETy9poahhog8mjtvIyXb7C3aXrdunbRu3drqbiUA0q5dO9HpdJKSklKi9zeZRA4eFNHpRNq3tw4xgEiTJiITJohs3SqSl1eaT0oF2QwwOTny6ejR8hIg+6KjRVq0EClFAXdRfK5Q2BlYU0NE5Bx31lQYjUZs374dVapUQfPmzQEAW7duRZcuXaBWq9GtWzfExcUhLi4OdevWdfj98vKUUXzNA+GdPp2/TqUCOnXKL/Rt1sxJH5KsTJs6FVVv3MCYDh2UKcd37VKqr++sqVmyBHjsMae+d5nN0j1y5Eg8++yz6NatW6kaSERE3unOmor27dsjJSUFqampeOmll7BgwQIAwP33349///vfGDhwIKpUqeLw+2RlAWvXKkFm5Urg+vX8dYGBQO/eSqHv4MFARISTPhzlu3gxP7zs3Impu3YB6emFt7vnHqBDB2UK8pgY4IEHXN/WvzgcajIzM9G3b1/UqVMHo0aNwogRI1CrVq2yaBsREXmYxMRETJ06FS1atMDZs2dhMBiwZ88eAEB4eDiCg4Mt2/r5+eHpp592aP9XrgDLlytBZv164Pbt/HWVKysBJjZWKQCuWNEpH4kAJazs3q2EGHOQuXCh8HZBQUDbtkp4MQeZhg09Zo6IEl1+un79Or766issXrwYhw8fRu/evfHcc88hNjYWAe6exesuePmJiMhxN2/eRMWKFS13Od1zzz1IS0sDANSsWRO1a9fGzp07MXXq1BLd0v3778olpaVLgR07rCeKjIpSQkxcnNIB4O/wT3Eq5PZtZaCeAr0w+P33wtv5+QEtWuT3wHTooAzo44bzfJldfgKAKlWqYNy4cRg3bhz27duHL774Ak8//TQqVqyIp556Cn/729/QuHHjEjeeiIjc6+TJk9Dr9dDr9Th27BguX74Mo9GIxMREVKxYEampqdBqtWjfvj38/Pws49TYw2RSwot5ILw7z6ft2+fXx7Rs6TGdAN7JaFQOsLkHZudO4OBBIDe38LZRUdY9MG3bAiEhrm9zKZQq8166dAnr1q3DunXroFarMXDgQBw5cgT33nsv3n33Xbz22mvOaicREZUhEcH+/fuRlJQEvV6PQ4cOWa3ftWvXXXthihsI8dYt4OeflR6Z5cuB1NT8dQEBQM+eSogZOhSoXbsUH6Q8E1EG5CnYA7N7N3DzZuFtq1WzroNp315Z5uUcDjW5ublYtmwZFi1ahHXr1qFVq1Z47bXXMHz4cMsASUuWLMHLL7/MUENE5CXmzp2LN9980/JcrVajR48e0Gq1iI2NRe0SJI3r14EVK5TemLVrgezs/HWhocCgQUqQ6d8fCAtzxqcoZ/78Mz+8mIPMlSuFtwsJUeZ+KNgLU6+eT3aBORxqIiMjYTKZ8MQTT2Dnzp1o06ZNoW369euH8PBwJzSPiIic6datW1i/fj30ej0efvhhDBo0CADQp08fVKhQAf3794dWq8WgQYNQuXJlh/d/6lT+ZaXNm5VLTWa1a+dfVureHdBonPWpyoHsbOX26YK9MCdPFt7O3x9o1cq6F6Z5c8CO0Z59gcOh5v3338ewYcMQFBRU5Db33HMPThccRICIiNzmxo0bWLlyJfR6PdasWYOsrCwAQE5OjiXU3Hfffbh27ZrV3Uv2EFGucJiDzJ3zTrZqpRT5xsYC0dE+2TngfHl5wJEj1r0whw8r9TF3atzYupC3TRugQgWXN9lTOBxqHL09j4iI3OP27duIjY3Fhg0bkJeXZ1leu3ZtxMXF4dFHH7UsU6lUdgcagwHYuDF/osiCd/6q1UC3bvn1MVFRTvs4vklEGUmw4CWkPXuUIqQ7RUTkBxhzHcw997i+zR6MN8eRxyirGYCJyovjx4/j4MGDeOSRRwAAQUFBSE1NRV5eHu69915otVrExcWhXbt2UDnYZXLjBrB6tVLou3o1kJmZvy4kRKmLiYsDBg5UxpOhIqSmWl9C2rXLelRBs0qV8i8hmf+3Vi12dRWDoYY8hnkGYMD6ToqCMzSTZ2MwdS0RwZ49e6DX65GUlISjR48iMDAQ/fv3R8W/Rqb717/+hWrVqqFJkyYO7z8lJf+yUnKyclXELCJC6YmJjQV69VLGZKM73Lyp9LoUDDFnzxbeTqNRLhsVrINp0kQZJ4YcwlBDHsN8IiwYbAoGmuJuGSX3YzB1jT179mDx4sVISkrC+fPnLcv9/f3RrVs3pKamWkLNAw4MWS+iDGFinl9p3z7r9c2b5xf6xsTwnGvFYAAOHbKugzl2zLpSGlB6Wpo1s66DadVKmfeBSs/pU2l6MM7S7R3snQGYPNOdM/nanNmXHJKVlSXZ2dmW53PnzrXMeB0SEiKPPPKIfPXVV5KWlubwvg0GkZ9/Fhk7VqRePevZrlUqkQceEHnvPZHjx533ebye0Sjy++8i//mPyN//LtKpk0hgYOHpwgGROnVEHn5Y5J13RDZsEOH5p0Q4S7cNnCbBe9w5AzB5F3PPjHnCQ/a0Oe7PP/+03LG0du1azJ8/H88++ywA4PTp09DpdNBqtejduzcqOHi3S2am9USRf814AEC5jNS3r9IbM3gwUL26Mz+Vl7p4sfB4MPZM7NihA2fadBJ7z98MNeRxeEL0DQymjjt//jyWLl0KvV6P5ORkq2kHRo4ciUWLFpV435cuKXcqLV2qjOxrMOSvq1pVCTBxcUCfPoCDd3X7loITO5qDjBdO7OhrynTuJ6KycmcNjfk5UPww7OQ5dDqdJdAYDAbodDr++xUjLS0N9evXtwoyLVu2tNyxFB0d7dD+RJSSDnOh76+/Wq9v1Ci/PqZz53IzNps188SOBXthPHxiR7o7hhryGLaKgm0VD5NnYzC9O5PJhN27d0Ov1yM1NRWff/45AGXQ0i5duiA3N9cSZBo1auTQvo1GYPv2/ELfEyes18fE5M943bx5OetUuHNix127lEDjoxM7llcMNeQxzDMA33niMz+3dwZgch8GU9tyc3OxadMm6PV6LF26FBf+upzh5+eHd955B9X+mkjwp59+QoCDv/6zs4GfflKCzIoVwNWr+es0GuDBB5UgM2QIULOm0z6SZ3NkYseqVQsPaOcDEzuWVww15DFKMwMweQYG08Lef/99JCYm4saNG5ZlFStWxIABA6DVahFSoAfA3kBz9aoSYJKSlEBTcPDZ8HDriSL/mmfYt9k7sWNwsBJaChbz+ujEjuUVQw0ROU15D6bXr1/H8uXL0adPH9SqVQsAUKlSJdy4cQPVqlVDbGws4uLi8OCDD951/jxb/vgjvz5m2zbr4U/q1cuvj+na1cdLPRyZ2PG++6zrYJo3V5aTz+K/LhFRKZw7dw5JSUlISkrCL7/8AqPRiHnz5mHcuHEAgIceeghNmzZF586doXagGtdkUs7Z5iBz9Kj1+ujo/CDTurWPdjbk5Skf3NwDY8/EjuZemHI+sWN5xVBDROSgtLQ0fPTRR9Dr9di7d6/VutatW6NygcmPKleujK5du9q139u3gQ0blBCzfLlyG7aZvz/QvbtS5Dt0KFC3rjM+iQcpOLGjuQdm716lZ+ZOnNiRisBQQ0RUDJPJhCtXriAyMhKAMh2E+bZ1lUqFLl26IC4uDnFxcWjQoIFD+05LUwbAW7oUWLPGupa1UiVgwAClN2bgQKVexmekphaug+HEjlRKDDVERDYYDAYkJydb7liKjIzEnj17AAChoaGYNGkSatWqhaFDh6K6g8Punj2bf9v1L79YX02pWVPpiYmLA3r08JEpgTixI7kIQw15Bc7+TK5w8+ZNrFmzBnq9HitXrkR6gaHws7KykJaWhnv+uswxdepUu/crotS2mutjDhywXt+yZX59TLt2Xn4O58SO5EYMNeQVOPszucLo0aPx9ddfW55HRERY7ljq2bMnAh044ebmAps25QeZlJT8dX5+QJcu+UGmYUNnfgoXMpmUEf4KXkLatw+wNS1GnTrWPTDt2gGcroacjKGGvIKtAdxsDfRGZI8zZ85Y7lj617/+hZYtWwIAhg4dih07dkCr1UKr1aJTp07wc6DbJCMDWL1aCTGrVlnPeRgcDPTrp4SYQYOUMd+8jr0TO4aHW/fAdOgA/FWPRFSWOKEleRVOdkklISI4fPgw9Ho9kpKSsG/fPss6nU6Hf/zjHwCUgmCVSgWVA0WoFy7kTxS5YYP1qPvVqysj+cbGAr17e9kdxo5O7FiwF4YTO5KTcZZuGxhqfANnfyZHHD9+HAMHDsTJAgO0+fn5oWvXrtBqtXjooYdQp04du/cnogyVYr6stHu39fomTfLnV+rY0UsmiuTEjuThOEs3+STO/kx3k5OTgw0bNiA7OxsPP/wwAKBevXpITU1FYGAg+vbtC61WiyFDhqCqA9d/8vKArVvzg8ypU/nrVCqgU6f8+phmzZz9qZyMEzuSD2OoIa/B2Z/JloyMDKxevRp6vR6rVq1CZmYmGjVqhIceeggqlQqBgYH46aef0KJFC1SsWNHu/WZlAevW5U8UWXAIlcBA5XJSXBwweLAyFpxHMk/sWLAHhhM7kg9jqCGvwNmf6U7ffvstvvzyS/z8888wGAyW5ZGRkejTpw9ycnIs8yt17NjRrn1euaKM5Lt0KbB+vXJVxqxyZSXAxMYCffsCDuQj1yk4saP5fzmxI5UjDDXkFTj7M506dQr169e33I30888/Y/Xq1QCAJk2aWO5Y6tChg0N3LP3+e/5lpe3blc4Ns6io/MtKXbp42FyI5okdC/bCcGJHKudYKExEHklEcODAAej1euj1ehw6dAjbtm3D/fffDwDYtm0bkpOTodVq0bx5c7v3azIBO3bkB5k762Hbt88PMi1bekjnxZ0TO+7apQxwx4kdqZxgoTAReR2j0YitW7dabr0+c+aMZZ1arcaBAwcsoaZz587o3LmzXfu9dQv4+ef8iSILXpEJCAB69lRCzNChQO3azvxEJcCJHYlKjKGGiDzGjh070L17d8vzChUqoF+/ftBqtRg8eLDV7NfFuX5dmSgyKQlYu9Y6E4SGKhNExsUB/fsDYWHO+wwOc3Rix4J1MJzYkcgKQw0RuVx6ejpWrlwJvV6P+vXr47333gMAdOrUCc2bN0eHDh2g1WrRt29fBAcH273fU6fyLytt2WJ9daZ27fzLSt27K3MnupwjEzu2bm1dB9O0qZdPCkVU9lhTQ0QucenSJSxbtgx6vR4bNmxA7l/jokRGRuL8+fOW4l4RsXtEXxElI5hnvD582Hp9q1b5A+FFR7u4UyM3V6l7MffA2DOxo7kXhhM7EllhTQ0ReYwnn3wSS5YsQcHfUM2bN4dWq0VcXJxViCku0BgMwMaNSpBZtsx65H61GujWLb8+JirK6R/FtoITO5p7YDixI5HLMdQQkdOICPbu3Yvly5dj0qRJ0Px1jScyMhIigo4dO1qCTNOmTe3e740b+RNFrl6tTBxpFhKi1MWYJ4p0oOym5O6c2HH3bqWRd+LEjkQuxctPXiIhIQFqtdrmAHM6nQ5GoxEJCQmubxiVe3l5ediyZYvljqVz584BANauXYu+ffsCAC781Z1Sq1Ytu/ebkpJfH5OcrNzVbBYRofTExMYCvXopcyqWGU7sSOR2Pnv56eOPP8Z7772HS5cuoUWLFpg3bx66du3q7maVObVabXPk3IIj7RK50m+//YbZs2dj+fLluF7gbp3g4GD0798fYQVuKbInzIgABw/mB5m9e63XN2+eX+gbE1NGNbOc2JHIq3lVqPn222/x6quv4uOPP8YDDzyA//u//8OAAQNw9OhR1K1b193NK1O2pgSwNXWAJ2Cvkm9KS0tDZmam5b81g8GAxYsXAwCqVKmCIUOGQKvVok+fPqhg52BveXnA5s1Kke+yZUCBYWmgUgGdOytFvrGxyphyTuXoxI4Fe2Cio23Ok8DvPpGbiReJiYmR0aNHWy1r1qyZvPXWW3a9Pj09XQBIenp6WTTPJRITEwWAaDQaASCJiYnublIh5jbe2bailpPnOn/+vHz00UfSu3dv8ff3l+HDh1vWmUwmiY+Pl40bN0pubq7d+8zIEPnf/0SeekrknntElD4a5REUJDJ0qMjnn4tcueLED2IyiZw9K/L99yITJoj06CFSsaL1m5sfVauKDBwokpAgsnKlSGqq3W/D7z5R2bD3/O01oSYnJ0fUarX8+OOPVsvHjh0r3bp1s2sfvhBqRMQSaDQajbubUqQ7/4jzj7r3OHbsmMyaNUtiYmIEgNXjgQceEJPJ5PA+L14U+b//U7KCRlM4Q4wcKZKUJJKV5aQPcf26yJo1IomJIkOGiNSoYTvABAeLdOsm8sYbIt9+K3L6tBKASoHffSLn87lQc+HCBQEgW7dutVo+Y8YMadKkic3X3L59W9LT0y2PlJQUrw813tBTY+ZNbaV8rVu3tgoy999/v7z77rty/Phxu/dhMokcPSoyc6ZIx46Fs0TDhkqO+OUXkby8UjY4K0tkyxaR998XeeIJZee2Aoy/v0h0tMhLLyldQQcPijjQw+QIfveJnMtnQ822bduslk+fPl2aNm1q8zVTp04t9EvTm0ONN/4C9IZepfLIYDDI+vXrZcyYMdKoUSPJzMy0rJs5c6b069dPFixYIBcvXrR7n3l5Ips3i7z5pkjjxoUzRUyMyIwZIocPl6IzJDdX5MABkYULRV54QaRNGxG12naIadxYZPhwkXnzRLZtE8nOLuGblgy/+0TO43OhpiSXn3ypp8Ybr9Xz16pnycrKkh9//FGefvppueeee6yC/v/+978S7lO5bDRqlEi1ataZQqMRGTBAZMECkQsXSrBzk0nk5EmRJUtEXn9dpEsX5XKRrQATEaEU40yfLrJunciff5bo8zgLv/tEzuVzoUZEKRR++eWXrZY1b968XBQKT506tcg/jImJiTJ16lTXNqgY3tir5MuSkpKkQoUKVkGmatWq8uyzz8ry5cvl1q1bdu8rNVXkiy9EYmNFKlSwzhbh4UrnyHffKQXBDrlyRWTFCpH4eCUNValiO8BUqiTSq5fIxIkiP/wgkpLicNdPWf73xO8+kfP5ZKhZsmSJBAQEyOeffy5Hjx6VV199VUJCQuTMmTN2vd6bQ4038cZeJV9y7tw5mT9/vqxZs8ay7NSpUwJA6tevL6+99pps2rRJ8hwoZvnjD5E5c0S6dhXx87POGHXrivz97yI//yxiMNi5w8xMkeRkkffeExk2TKRePdsBRqMR6dBBZMwYkS+/VAp1jEbHDogNZfUd5XefqGzYe/72qnFqHnvsMVy/fh2JiYm4dOkSWrZsiVWrVqFevXrubhoVYDQabY6dY35uLDh1MpWaiODYsWOWEX13794NABgyZAj69esHAIiKisLRo0fRrFkzuyaLNJmUYVvMA+EdPWq9Pjo6fyC81q2LGTT3zokdd+1SdujGiR3LatwnfveJ3IvTJBB5KRHBP/7xD3z//fc4fvy4ZblKpULnzp3x+OOP45VXXrF7fzk5wIYN+RNFXrqUv87fH+jePX+iyCJ/RzgysWPt2vmD2blpYkdzkNFoNDAYDB43kCURKew9fzPUEHkJg8GA/fv3IyYmxrKsW7du2Lx5MzQaDXr37o24uDgMHToUNWrUsGufaWnAypVKkFmzBrh5M39dxYrAwIFKkBkwALjnHhs7cHRiR3MPjAdN7BgYGAiDwQCNRoMcW+GLiNzOZ+d+IipPbt68ibVr10Kv12PFihXIzMzElStXULVqVQDAW2+9hTFjxmDAgAF2B/WzZ/MvK23apMwWYFazZv5EkT173nHlxzyxY8EQU9TEjtHR1vMiNWrkkRM76nQ6S6AxGAzQ6XTsqSHyYgw1RB7m+vXrWLZsGfR6PX766Sfcvn3bsq5GjRo4fvy4JdQMHDiw2P2JAPv3K/MrLV2qTG9UUMuW+fUx7dr9NVFkTo6yYcE6mN9+K7zzghM7mnthvGRixztraMzPATDYEHkphppyjhPweQYRsRTw6vV6vPDCC5Z1DRs2hFarRVxcHDp16gS1Wl3s/nJzlV4Yc33MuXP56/z8gC5d8oNMw/oFJnZc/FcvTCkndvR0toqCbRUPE5F3Yagp59Rqtc0/4gX/6JPziQiOHDliuWNpxIgRGDt2LADlrqW2bdsiNjYWWq0WLVu2tOuOpYwMpS5m6VKlTiY9PX9dcDDQty8QFysY0iYFlU/+FV6e3wns2QNkZhbeYdWq1peQOnQAqlVz1iFwK96lROSbWChMRXbD804Q5zKZTNixYwf0ej30ej1OnjxpWdezZ09s2LDB4X1euKD0xCxdqty5VLBzpXp14PG+f+KJRrvQzrgTAfv/CjJXrhTeUXAw0L69dS9MvXoeWQdDROUP736ygaGmaLy1tWzl5eWhUaNGOHv2rGVZYGAg+vTpA61WiyFDhqCaHb0gIsCRI/mFvrt25a+rgGzE1tmHJxrtwv3+O1H11E6oCgQnC7VaGf+lYC9M8+bKfdtERB6IocYGhpq7462tzpGZmYnVq1dj//79mDlzpmX5oEGDsGXLFgwePBharRb9+/dHRTvqUfLygG3b8gt9T50C1MjDvTiKjtiJwdV3opN6F6qnHoLK1mWTxo2te2DatAEqVHDeByYiKmO8pZscwltbSyc1NRXLli1DUlIS1q9fbwmFo0ePRt26dQEACxcuRNWqVaHRaIrdX1YWsG6dEmJWLBeE/nkaHbALf8NOdFLtRDu/vQgyZv/15gVeGBFh3QPTvj1QubKzP65NLDonIndjqCGvubXVE0+aq1atwjvvvIOtW7fCVGDY/8aNG0Or1VrdqVSzZs277is1FVi+HNj4bSqyknehTe5OPIpdmIOdqIrr+RsKACOASpXyC3jNQaZWLbfVwbDonHyRJ/7dobsos9mnPBAntCzMmybgc3dbTSaT7N+/Xy5cuGBZtmTJEsus1+3atZPp06fL4cOHxWTnrNG/78mUb15Klvn13pNvMUxOo57NiR1NZTSx492UZCZrzlBNvsbdf3dI4ZOzdJcWQ01hJTlxuZOrT5p5eXnyyy+/yGuvvSZRUVECQHQ6nWV9enq6fPDBB3L27Nnid2YwiHHXHjk5/hPZed+z8ntgS8mDX+EAo1LJrQbNxTRihMi//iWyc6fI7dtl8vnupqR/zM3rNRoN/+iTT2BYdz+GGhsYanxDWZ80c3NzZcWKFfLcc89JtWrVLD0xACQoKEjGjx9f/E6MRpHffxf56ivJ/dtY+bNpJzGoA232wqQG1pYTbR6SG2+/I7Jhg4gHfT9L+sfc/G+j0Whc0UyiMsew7l72nr959xN5JWffqWU0Gi31Lzk5OahWrRoy/xqQLjw8HEOGDIFWq0Xfvn0REhJSeAeXLllNKWDauQt+6TcKbZaGcOxVxyCjaQdUHRiD1s93QGhTz5jYsSiO3u7P4QHIV/EOUfex+/ztkojlIdhT4xuc9Yvp0qVLsmDBAunXr5+0aNHCqg5m3LhxMmbMGFm/fr0YDAbrF964IbJ+vcisWSJarUitWjZ7YLIRJFtxv7yPcfJK5a8lYfhxWbfWJDk5pfn07mFvzwu76clXsafGvXj5yQaGGu9X2pPmH3/8Ie+995507txZVCqV1aWl3377rfALbt8W+fVXkfnzRZ5+WqRZM5sBJg9+cgD3yUI8Jy9igbTBXml7n0GmTBHZvVvEzrphj2TvH3MWVJKvYlh3P4YaGxhqvFtpT5pTpkyxCjEAJCYmRmbOnCnHjh0TycsTOXJEZNEikb/9TaR9e5GAAJsh5mql+rK0wqPyOuZIV2ySEGSKWi3Ss6fIvHkip06V4YFwIUf+mHtb0TmRPRjWPYO952+OU0New95JCPPy8rBlyxYkJSXhmWeeQdu2bQEAnTp1gr+/P3r06IG42Fg81KEDIs+fV+pgXn65yIkdTVWq4lKdGGzP7YAlp2Kw6VYHXMtUpjQICQH69wc+iQUGDXLZOHcu4ehM1ncbq2PKlClISEgoclBHjvdBnoqTn3oXhhovwMGfFHf7jG+++SbWr1+PZ599FsuWLcP168pgdYGBgUqo+fNP9DGZkDFhAiocOgRMn170xI7t2iGjeQy258Xg6+Md8N/t9ZF3PX9Auxo1gBeGAnFxQK9eQFCQkz+oh3D2H3MOzkfeqLiwTh7GRT1HHsFbLz+x+7No165dk4cfflhCQkLyb7sGpF+lSvJ1hw5yqVcvkUaNbF5CErVaJDpa5KWXxLTwM/n9+4OSGJ8rbdsW3rR5c5G33hLZvr3Mx7zzaaxNIKKS4OUnH2Kry9/WpYHy4OLFizhx4gS6deuGhIQE+JlMuLZhAx7PykLPkBB0r1ABNf/8E36ZmdZTWAOFJnbMa9kGm3dXUGa8ngGcOZO/qUoFdO6s9MbExiovpdIr+F2ePn06b/kmIqfiODVepLyO/3H8+HHo9Xrof/wRqTt3om94OD4ZNQopP/yAKufOwcaoMbgEIKNZMzR96imriR1v3gTWrlVmvF65EkhLy39NUBDQt68SYgYPBqpXd9EHLIc43gcROcLe8zdDjZcpLyeDgwcPYuWiRTj3ww+ISElBBwAxAKra2DYDwPWoKEQ9+ij+d/YsXluyBC9Nm4Ypf/VsXb4MLFumzHj9889AwcNWtaoSYGJjgT59lMJfKlvlNZwTUclx8D0bvLWmxsyXB3/KTUuTvJ9/FnnvPZFhw+R6aKjtiR0DAgpN7JiYkGB1XKZNS5SjR5Wx8Tp1Krybhg1F3nhD5JdflLu4yXVYU0NEJcFxamzw5lDjUycDg0Fkzx7J+eADOfPgg3IuPFzybAQYIyA3ataUnMcfL3Zix4CAIAEeED+/f0rjxoWDTEyMyIwZIocPe/dAeN6MBe9EVFIsFPYhjo4X4lFMJuDECaVod+dO5G3bBtWBA1Dn5kIDoF6BTdMqVsQ9ffsqhbwdOsCvXTuEhYUVuevsbGD9emDatH3IzT0LoDpMJuCPPwCNRrndOjYWGDoUqFmzrD8oFYfjfRBRWWNNjRfwqnFq7pjYEbt2ATduFNosDcBOAMdDQ1Ghe3e0GDkSHYYOhb//3XP21avAihVKfcy6dcCtW/nrwsOBunUP4eDBREye3B7Tp0905icjIiI3YaGwDd4aajxWeroyCq85xOzcCVy4UHi7oCAgOhrSoQNeW7IEJ8LD0e6xxxCn1aJNmzZQqVSFX1PAiRNKiFm6FNi6Ven8yXcWHTtewcyZMejaFQgIsN2zRURE3sve8zcvP5F9cnKAAwfye2B27gR++63QZiaVCn8EBGCzwYBdAI5VrIifr1xBQHAwVAASpk1DeHj4Xd/KZAJ271Zuu166FDh61Hp9dLRyWenSpQWoWfMq4uN5OYOIiBhqyBaTSQks5vCyc6cSaHJzC29bvz5+Cw3Fd2fOYH1GBvaKIMtgQEBAAHr37o2n4uJg9PNDwF+bFxVocnKADRuUELNsmXIVy8zfH+jePb8+pp6lEGd0kR+BPTRE5Am8qnzABzDUlHcigHlSR3MvzO7dNid2RNWqyGvXDn+EhaH+o4+iQrduQLVqWJKQgGnTpqFixYoYNHAgtFotBg4cWOwlvrQ0YNUqJcisXg3cvJm/rmJFYOBAJcgMGADcc4+TP7eD+IeJiEqCc565FkONk3jNSe/PP/MLeM1B5i4TOyImBhnNmmFdWhr+s3kz1v30E27fvo3vH30UD1dTZqoeMWIEYmJi8OCDDyIwMPCub3/2bH59zC+/AHl5+etq1lR6YmJjgZ49gWJ25VL8w0REJcFpblzMBbeXe4yyHKfGI8fgyM4W2bpV5P33RZ54oviJHV98UeSzz0QOHpTrV67IBx98ID169BA/Pz/LZJEAJCoqSr7++mu7mmAyiezdKzJ1qkibNoXfukULkUmTlCFoPH2iSJ8aK4iIXMqXB091BQ6+Z0NZD77n1pNebq7IgQMiCxcq4aRNGyWs2AoxjRqJPPmkyLx5SujJzhaTySSZmZmW3Z08edIqyLRu3VoSEhLkwIEDYipm9DqDQWT9epFXXhGpW9f6rf38RLp1E/nnP0VOnCjrg5Jv6tSpRf47JCYmytSpU+3aD/8wEVFJmf9uaDQadzfF6zDU2OCKEYVdctIzmUROnhRZskTk9ddFunQRCQ62HWBq1BAZOlREpxNZu1bk+nXLboxGo2zbtk3Gjx8vjRo1kqFDh1q9zahRo2Tu3Lly6tSpYpuUni7y7bdKVgoPt25CcLBIXJzI4sUiV686/WjYxZk9afzDRESO4g+i0mGoscFV0yQ4/aR35YrIihXKNZwBA0SqVLEdYCpVEunZU2TiRJEffhA5d67QnAA5OTmyZs0aeemllyQiIsKqNyY8PFxycnLsbtaFCyKffCLSr5+IRmPdlOrVRZ57TmTZMuUqmCdwRk8a/zARkaN46br0GGps8OSeGsvlkcxMkU2bLBM7Sv36tgNMwYkdFy8WOXrUrqKUfv36WQWZ0NBQeeKJJ+S7776TjIyMu77WZBI5dEhk+nTlre9sUpMmIuPHi2zZ4rkTRZYmlPAPExE5yiPrLb0QQ40NHldT89fEjvLJJ7K3bVs5CIhRpbI5sWNqtWoizzxT7MSOZqmpqfLZZ5/JkCFD5GqBaz5z5syRiIgIeemll2TNmjXF9szk5ioZ6/XXldmtCzZNpRK5/36Rd94ROXbMsWPlTiXpSeMfJiIqCWfV85V3DDU2uPXup4QEkd9/F/nqK5GxY0U6dRIJCrLZC3MjNFTkoYdkfZ8+0hOQ2ZMm2dWG06dPy9y5c6Vbt25Wdyx98cUXlm1u3bolxmJ6dG7eFPnxR5ERIwpf6QoMFBk0SOTTT0UuXXL4MLldqXvSitgn/zAREZUdhhobyjLUFDrpXbwokpQkMmmSnGjYULKLCDASHi7Sp4/I5MkiS5fK3PHjHTrpTp06VV5++WVp06aN1WUlABIZGSk9e/aUE3bcZnTlinI395AhhbNW5cpKJ9EPPyhXx7wVLx8REXkne8/fHHzPGdLTkdCtmzKQ3UMPFZrYsaH5/wQGAm3bAjExQIcOyv82agQUmNDxtaFD8dYHH8BgMECj0RQamMloNGL79u0QEXTt2hVqtRqffPIJAMDPzw9du3aFVqvFuXPnMHfuXLz88sto2LAhbDl+XBkELykJ2L5diTBmUVHKIHixsUCXLspUBd7M1mBXtgbFIiIi7+XlpyoP0b27MjdSQX5+QIsW+eElJgZo2VKZRvoudDqdJdAYDAbodDpMmDABP//8M5KSkrB06VKkpqaiR48e2Lhxo9WJecKECZg1axZ0Oh3mzp1baLRKkwn49df8EX3vnI+yXTslxMTFKU0tZvJsr2I0Gm2O3snJL4mIfIhrOo48Q5ldfnruOeUupUcfFZkzR6msLcF1mjsvhzz22GNWl6LMj/DwcBkxYoRVbUxRtSK3bil3g7/wgjJkzZ03UPXtK/LRRyIpKc45FERERM5m7/lbJVLwooNvy8jIQFhYGNLT04udbNEheXmlvj5j6/LIwIEDsXr1agBAxYoV8fTTT0Or1aJHjx4IsNHjExgYCIPBgICACHz22SUsXQqsXQtkZeVvExpqPVFkWFipmk1ERFTm7D1/8/KTM5Qi0Jw6dQp6vR4LFy6ESqXCiBEjLOuef/55tGrVCmlpaYiIiMC0adOK3M/rr8+HwfAyVKo45OZ2QYHdoHbt/PqY7t0Bjca+tnnNJJ1ERERgqHE5EcGBAweg1+uh1+tx6NAhq/VbtmzBk08+CQB46KGH8NBDDxWxH2DvXqXId+HCK7hy5e+W5YoD6N49Hf/8Zze0bVuy+hjOTE1ERN7EK0LNmTNnoNPpsGHDBly+fBk1a9bEU089hcmTJ0Njb7eDh9Dr9Xj44Yctz9VqNXr06AGtVouhQ4eiTp06Rb7WYACSk5Ui32XLgPPnzWtqwM/PhG7d/BAXBwwdCnz11TLEx8dj1apEtGvn2F09d/bQFAw2Dz74IDZs2GCz6JaIiMidvCLU/PbbbzCZTPi///s/NGrUCIcPH8YLL7yArKwszJkzx93Ns+n27dtYv3499Ho92rRpg7//XelJ6d27N8LDw9GzZ09otVoMGjQIlStXLnI/6enA6tVKj8zq1UBGRv66kBCgbt2jaN78OBYujEPB3ZTmrh5bPTTx8fGYNm0ajEYjevXq5ZZAw8thRER0V66oWi4L7777rkRFRTn0mrKeJuHGjRvy9ddfyyOPPCIhISGWu5XatWtntZ3BYLjrfs6dU2ZD6NNHuUPpzkm3X3hBZOVK5c6msnLnnVhqtVoAiFqtLrs3dbBNxS0nBUdDJiJv5/OD76Wnp9+1hwMAcnJykJOTY3meUbCbw8keffRRJCUlITc317Ksdu3aiIuLg1artdr2zjuXRIBDh/IHwtu713rfzZvnF/rGxChD4JQ1Wz00arUaRqMROp3OLT01Bdu0YcMG9OrVy/K84OUw9tpYY20UEZUbLgpZTnXixAkJDQ2VhQsX3nW7qVOnFpo6AGXUU2MeU+bee++VSZMmya5du8RkMhW5fW6uyIYNIuPGFZ6IW6USeeABkXffVaaLcqc7e2g8oVfE3Abzo2BbPKF9nohTRHgH9qoR2eYVcz8VFToKPnbt2mX1mgsXLkijRo3kueeeK3b/t2/flvT0dMsjJSWlzELNkSNH5PdiEkhmpsj334s8/bQyn1LBIBMUJDJ0qMjnnyvzMHmCXr16WQKNp50QCw5I6Ent8mQlncyTXIeXWIls84pQc/XqVTl27NhdH7cKFI5cuHBBmjRpIk8//XSxM03bUtY1NbZcuqTMaD1woDLDdcEgU7WqyMiRInq9MjO2JzH/Ee3Vq5fV84IBwl2/Gu88OfNEbT/zcdJoNO5uChWBvWpEhXlFqHHE+fPnpXHjxvL4449LXl5eifbhilBjMokcPSoya5ZIp07KpaSCQaZhQ5E33hD55ReREn6MMufJvxaL+oPPE3Xx2FPjPfhvRWTNp0KN+ZJTr1695Pz583Lp0iXLwxFlFWry8kS2bBF5802Rxo2tQwwgEhMjMmOGyOHDSujxdJ56Xd9WqLpbfQ3l469/78NeNaJ8PhVqFi1aVGTNjSPKKtRER1uHGI1GpH9/kU8+EblwwalvVa7dGbYKnpgTExOlR48ePFHb4Mk9b2Qbe2qIrPlUqHGWsgo1o0eLhIeLDB8u8t13IhkZTt092cATtf08teeNbGOvGlFhnKXbhrKapfvPP4FKlQAbE2dTGeHowuSLCo4dVNSYQpyehMoje8/fDDVERB6CYZ3INoYaGxhqiIiIvI+9528XDLhPREREVPYYaoiIiMgnMNQQAOVavk6ns7lOp9PxOj4REXk8hhoCkD+T853BxnzXhVqtdlPLiIiI7OPv7gaQZzDfbREfH295zttIiYjIm/DuJ7JiDjIajQYGg4GBhoiI3I63dNvAUGOfwMBAGAwGaDQa5OTkuLs5RERUzvGWbioRnU5nCTQGg6HI4mEiIiJPw1BDFgVraHJycpCYmGizeJi8E+9wIyJfx0JhAmB7bhlbxcPkvcx3uAEocl4hIiJvxlBDAACj0WizKNj83Gg0uqNZ5ES8w42IfB0LhYnKGd7hRkTehnc/2cBQQ6TgHW5E5E1495MXYiEnuQLvcCMiX8VQ40E4VQGVNd7hRkS+jIXCHoSFnFSWeIcbEfk6hhoPU/AkM336dBZyktPwDjci8nUsFPZQLOQkIiJSsFDYi7GQk4iIyHEMNR6GhZxEREQlw5oaD8JCTiIiopJjqPEgLOQkIiIqORYKExERkUdjoTARERGVKww1RB6AU2QQEZUeQw2RB+AUGUREpcdCYSIPwCkyiOyTkJAAtVpt878JnU4Ho9HIns1yjKGGyENwigyi4pl7NQHrIS4K/gig8ot3PxF5GE6RQXR3d/ZislfT99l7/mZPDZEHsTVFBv9IE1ljryYVhYXCRB6CU2QQ2W/KlCmW8K/RaBhoCAB7aog8AqfIIHIMezXJFoYaIg/AKTKI7FdUTQ3A8F/esVCYiIi8RlFFwSwW9m0sFCYiIp/DXk26G/bUEBERkUfjhJZERERUrjDUEBERkU9gqCEiIiKfwFBDREREPoGhhoiIiHwCQw0RERH5BIYaIiIi8gleF2pycnLQpk0bqFQq7N+/393NISIiIg/hdaFmwoQJqFmzprubQURERB7Gq0LN6tWrsW7dOsyZM8fdTSEiIiIP4zVzP125cgUvvPACkpKSEBwcbNdrcnJykJOTY3mekZFRVs0jIiIiN/OKnhoRwciRIzF69Gi0b9/e7tfNmjULYWFhlkedOnXKsJVERETkTm4NNQkJCVCpVHd97N69G/Pnz0dGRgbefvtth/b/9ttvIz093fJISUkpo09CRERE7ubWWbqvXbuGa9eu3XWb+vXr4/HHH8fy5cuhUqksy41GI9RqNYYPH44vv/zSrvfjLN1ERETex97zt1tDjb3OnTtnVQ9z8eJF9OvXD99//z06duyI2rVr27UfhhoiIiLvY+/52ysKhevWrWv1vGLFigCAhg0b2h1oiIiIyLd5RaEwERERUXG8oqfmTvXr14cXXDUjIiIiF2JPDREREfkEhhoiIiLyCQw1RERE5BMYaoiIiMgnMNQQERGRT2CoISIiIp/AUENEREQ+gaGGiIiIfAJDDREREfkEhhoiIiLyCQw1RERE5BMYaoiIiMgnMNQQERGRT2CoISIiIp/AUENEREQ+gaGGiMiDJSQkQKfT2Vyn0+mQkJDg2gYReTCGGiIiD6ZWqxEfH18o2Oh0OsTHx0OtVrupZUSex9/dDSAioqJNmTIFABAfH295bg40iYmJlvVEBKhERNzdCFfJyMhAWFgY0tPTERoa6u7mEBHZzRxkNBoNDAYDAw2VK/aevxlqiIi8RGBgIAwGAzQaDXJyctzdHCKXsff8zZoaIiIvoNPpLIHGYDAUWTxMVJ6Vq5oac6dURkaGm1tCRL5q1qxZ8PPzw8SJEwutmz17NkwmE95++22H9jl79mzMnDkTkyZNwsSJEzF79mzEx8fj9u3bNt+HyNeYz9vFXVwqV5efzp8/jzp16ri7GURERFQCKSkpqF27dpHry1WoMZlMuHjxIipVqgSVSuW0/WZkZKBOnTpISUlhrU4Z47F2DR5n1+Bxdg0eZ9coy+MsIsjMzETNmjXh51d05Uy5uvzk5+d314RXWqGhofwPxkV4rF2Dx9k1eJxdg8fZNcrqOIeFhRW7DQuFiYiIyCcw1BAREZFPYKhxgsDAQEydOhWBgYHuborP47F2DR5n1+Bxdg0eZ9fwhONcrgqFiYiIyHexp4aIiIh8AkMNERER+QSGGiIiIvIJDDVERETkExhq7PTxxx8jKioKQUFBaNeuHTZv3nzX7Tdt2oR27dohKCgIDRo0wIIFC1zUUu/myHH+8ccf0adPH1SrVg2hoaG4//77sXbtWhe21ns5+n0227p1K/z9/dGmTZuybaAPcfRY5+TkYPLkyahXrx4CAwPRsGFDfPHFFy5qrfdy9Dh//fXXaN26NYKDgxEZGYlRo0bh+vXrLmqtd/rll18wZMgQ1KxZEyqVCklJScW+xuXnQqFiLVmyRAICAmThwoVy9OhRGTdunISEhMjZs2dtbn/q1CkJDg6WcePGydGjR2XhwoUSEBAg33//vYtb7l0cPc7jxo2T2bNny86dO+X48ePy9ttvS0BAgOzdu9fFLfcujh5nsxs3bkiDBg2kb9++0rp1a9c01suV5FgPHTpUOnbsKD/99JOcPn1afv31V9m6dasLW+19HD3OmzdvFj8/P/nggw/k1KlTsnnzZmnRooXExcW5uOXeZdWqVTJ58mT54YcfBIDo9fq7bu+OcyFDjR1iYmJk9OjRVsuaNWsmb731ls3tJ0yYIM2aNbNa9tJLL0mnTp3KrI2+wNHjbMu9994r06ZNc3bTfEpJj/Njjz0m//jHP2Tq1KkMNXZy9FivXr1awsLC5Pr1665ons9w9Di/99570qBBA6tlH374odSuXbvM2uhr7Ak17jgX8vJTMQwGA/bs2YO+fftaLe/bty+2bdtm8zXbt28vtH2/fv2we/du5ObmlllbvVlJjvOdTCYTMjMzUbly5bJook8o6XFetGgRTp48ialTp5Z1E31GSY71smXL0L59e7z77ruoVasWmjRpgjfffBO3bt1yRZO9UkmOc+fOnXH+/HmsWrUKIoIrV67g+++/x6BBg1zR5HLDHefCcjWhZUlcu3YNRqMRNWrUsFpeo0YNXL582eZrLl++bHP7vLw8XLt2DZGRkWXWXm9VkuN8p3/+85/IysrCo48+WhZN9AklOc5//PEH3nrrLWzevBn+/vyTYa+SHOtTp05hy5YtCAoKgl6vx7Vr1/C3v/0Nf/75J+tqilCS49y5c2d8/fXXeOyxx3D79m3k5eVh6NChmD9/viuaXG6441zInho7qVQqq+ciUmhZcdvbWk7WHD3OZv/973+RkJCAb7/9FtWrVy+r5vkMe4+z0WjEk08+iWnTpqFJkyauap5PceQ7bTKZoFKp8PXXXyMmJgYDBw7E3LlzsXjxYvbWFMOR43z06FGMHTsW8fHx2LNnD9asWYPTp09j9OjRrmhqueLqcyF/dhWjatWqUKvVhRJ/ampqoQRqFhERYXN7f39/VKlSpcza6s1KcpzNvv32Wzz33HP43//+h969e5dlM72eo8c5MzMTu3fvxr59+/DKK68AUE68IgJ/f3+sW7cOvXr1cknbvU1JvtORkZGoVasWwsLCLMuaN28OEcH58+fRuHHjMm2zNyrJcZ41axYeeOABjB8/HgDQqlUrhISEoGvXrpg+fTp7053EHedC9tQUQ6PRoF27dvjpp5+slv/000/o3Lmzzdfcf//9hbZft24d2rdvj4CAgDJrqzcryXEGlB6akSNH4ptvvuH1cDs4epxDQ0Nx6NAh7N+/3/IYPXo0mjZtiv3796Njx46uarrXKcl3+oEHHsDFixdx8+ZNy7Ljx4/Dz88PtWvXLtP2equSHOfs7Gz4+Vmf/tRqNYD8ngQqPbecC8usBNmHmG8X/Pzzz+Xo0aPy6quvSkhIiJw5c0ZERN566y15+umnLdubb2N77bXX5OjRo/L555/zlm47OHqcv/nmG/H395ePPvpILl26ZHncuHHDXR/BKzh6nO/Eu5/s5+ixzszMlNq1a8sjjzwiR44ckU2bNknjxo3l+eefd9dH8AqOHudFixaJv7+/fPzxx3Ly5EnZsmWLtG/fXmJiYtz1EbxCZmam7Nu3T/bt2ycAZO7cubJv3z7LrfOecC5kqLHTRx99JPXq1RONRiNt27aVTZs2WdaNGDFCunfvbrV9cnKyREdHi0ajkfr168snn3zi4hZ7J0eOc/fu3QVAoceIESNc33Av4+j3uSCGGsc4eqyPHTsmvXv3lgoVKkjt2rXl9ddfl+zsbBe32vs4epw//PBDuffee6VChQoSGRkpw4cPl/Pnz7u41d5l48aNd/2b6wnnQpUI+9qIiIjI+7GmhoiIiHwCQw0RERH5BIYaIiIi8gkMNUREROQTGGqIiIjIJzDUEBERkU9gqCEiIiKfwFBDREREPoGhhoiIiHwCQw0RERH5BIYaIvJaV69eRUREBGbOnGlZ9uuvv0Kj0WDdunVubBkRuQPnfiIir7Zq1SrExcVh27ZtaNasGaKjozFo0CDMmzfP3U0jIhdjqCEirzdmzBisX78eHTp0wIEDB7Br1y4EBQW5u1lE5GIMNUTk9W7duoWWLVsiJSUFu3fvRqtWrdzdJCJyA9bUEJHXO3XqFC5evAiTyYSzZ8+6uzlE5CbsqSEir2YwGBATE4M2bdqgWbNmmDt3Lg4dOoQaNWq4u2lE5GIMNUTk1caPH4/vv/8eBw4cQMWKFdGzZ09UqlQJK1ascHfTiMjFePmJiLxWcnIy5s2bh//85z8IDQ2Fn58f/vOf/2DLli345JNP3N08InIx9tQQERGRT2BPDREREfkEhhoiIiLyCQw1RERE5BMYaoiIiMgnMNQQERGRT2CoISIiIp/AUENEREQ+gaGGiIiIfAJDDREREfkEhhoiIiLyCQw1RERE5BMYaoiIiMgn/D9gKXxSr7aWZgAAAABJRU5ErkJggg==", "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": [ "##
OPTIONAL SLIDES
" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Some Useful Matrix Calculus\n", "\n", "When doing derivatives with matrices, e.g. for maximum likelihood estimation, it will be helpful to be familiar with some matrix calculus. We shortly recapitulate used formulas here. \n", "\n", "- We define the **gradient** of a scalar function $f(A)$ w.r.t. an $n \\times k$ matrix $A$ as\n", "$$\n", "\\nabla_A f \\triangleq\n", " \\begin{bmatrix}\n", "\\frac{\\partial{f}}{\\partial a_{11}} & \\frac{\\partial{f}}{\\partial a_{12}} & \\cdots & \\frac{\\partial{f}}{\\partial a_{1k}}\\\\\n", "\\frac{\\partial{f}}{\\partial a_{21}} & \\frac{\\partial{f}}{\\partial a_{22}} & \\cdots & \\frac{\\partial{f}}{\\partial a_{2k}}\\\\\n", "\\vdots & \\vdots & \\cdots & \\vdots\\\\\n", "\\frac{\\partial{f}}{\\partial a_{n1}} & \\frac{\\partial{f}}{\\partial a_{n2}} & \\cdots & \\frac{\\partial{f}}{\\partial a_{nk}}\n", " \\end{bmatrix}\n", "$$\n", " \n", "\n", " \n", "- The following formulas are useful (see Bishop App.-C)\n", "$$\\begin{align*}\n", "|A^{-1}|&=|A|^{-1} \\tag{B-C.4} \\\\\n", "\\nabla_A \\log |A| &= (A^{T})^{-1} = (A^{-1})^T \\tag{B-C.28} \\\\\n", "\\mathrm{Tr}[ABC]&= \\mathrm{Tr}[CAB] = \\mathrm{Tr}[BCA] \\tag{B-C.9} \\\\\n", "\\nabla_A \\mathrm{Tr}[AB] &=\\nabla_A \\mathrm{Tr}[BA]= B^T \\tag{B-C.25} \\\\\n", "\\nabla_A \\mathrm{Tr}[ABA^T] &= A(B+B^T) \\tag{B-C.27}\\\\\n", " \\nabla_x x^TAx &= (A+A^T)x \\tag{from B-C.27}\\\\\n", "\\nabla_X a^TXb &= \\nabla_X \\mathrm{Tr}[ba^TX] = ab^T \\notag\n", "\\end{align*}$$" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Derivation Predictive Distribution\n", "\n", "- In the [derivation of the predictive distribution](#predictive-distribution), we replaced $\\mathcal{N}(w\\,|\\,m_N,S_N)\\,\\mathrm{d}w$ with $\\mathcal{N}(z\\,|\\,x_\\bullet^T m_N,x_\\bullet^T S_N x_\\bullet)\\,\\mathrm{d}z$. Here we discuss why that is allowed. \n", "\n", "- Since $z = x^T w$ (drop the bullet for notational simplicity), we have\n", "$$p(z) = \\mathcal{N}(z|m_z,S_z)$$\n", "with\n", "$$\\begin{aligned} m_z &:= E[z] = E[x^T w] = x^T E[w] = x^T m_N \\\\ \\Sigma_z &:= E[(z-m_z)(z-m_z)^T] \\\\   &= E[(x^T w - x^T m_N)(x^T w - x^T m_N)^T] \\\\   &= x^T E[(w - m_N)(w - m_N)^T]x \\\\   &= x^T S_N x \\end{aligned}$$\n", "\n", "- Then we equate probability masses in both domains:\n", "$$ \\mathcal{N}(z|m_z,\\Sigma_z)\\mathrm{d}z = \\mathcal{N}(w|m_N,S_N)\\mathrm{d}w$$\n", "\n", "$$ \\Rightarrow \\mathcal{N}(z|x^T m_N,x^T S_N x)\\mathrm{d}z = \\mathcal{N}(w|m_N,S_N)\\mathrm{d}w $$" ] }, { "cell_type": "code", "execution_count": 1, "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\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "anaconda-cloud": {}, "celltoolbar": "Raw Cell Format", "kernelspec": { "display_name": "Julia 1.8.2", "language": "julia", "name": "julia-1.8" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.8.2" } }, "nbformat": 4, "nbformat_minor": 4 }