2018-01-09-bayesian-linreg-posterior.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Linear Regression part 3: Posterior\n",
    "\n",
    "\n",
    "Now I have [priors on the weights]({% post_url 2018-01-03-bayesian-linreg %}) and [observations]({% post_url 2018-01-08-bayesian-linreg-sample %}). In this post, I'll show a formula for finding the posterior on the weights, and show one plot using it. The [next post]({% post_url 2018-01-10-bayesian-linreg-plots %}) will have more plots.\n",
    "\n",
    "My goal is to find the distribution of the weights given the data. This is given by\n",
    "\n",
    "$$p(\\mathbf w | \\mathcal D) = \\frac{p(\\mathcal D | \\mathbf w)p(\\mathbf w)}{p(\\mathcal D)}.$$\n",
    "\n",
    "Since I'm assuming the prior is Gaussian and likelihood is a combination of Gaussians, the posterior will also be Gaussian. That means there is a [closed form expression](https://en.wikipedia.org/wiki/Conjugate_prior)\n",
    "for the mean and covariance of the posterior. Skipping ahead, I can use the equations from [\"Computing the Posterior\" in the class notes](http://www.inf.ed.ac.uk/teaching/courses/mlpr/2017/notes/w6c_bayesian_regression.html#computing-the-posterior):\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    V_N &= \\sigma_y^2(\\sigma_y^2V_0^{-1} + \\Phi^\\top \\Phi)^{-1}\\\\\n",
    "    \\textbf w_N &= V_N V_0^{-1}\\textbf w_0 + \\frac{1}{\\sigma_y^2}V_N \\Phi^\\top\\textbf y.\n",
    "\\end{aligned}$$\n",
    "\n",
    "## Code\n",
    "\n",
    "I'll convert this to code. Heads up, I know this isn't the most efficient way to do this. I'll try to update this when I find more tricks.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports!\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# helper functions you can skip over :D\n",
    "SAVE = True\n",
    "def maybe_save_plot(filename):\n",
    "    if SAVE:\n",
    "        plt.tight_layout()\n",
    "        plt.savefig('images/' + filename, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variables\n",
    "\n",
    "\\\\( w_0 \\\\) and \\\\( V_0 \\\\) are the prior's mean and variance, which I defined back in [priors on the weights](2018-01-03-bayesian-linreg.ipynb). The code for that was"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_w = 0\n",
    "mu_b = 0\n",
    "\n",
    "sigma_w = 0.2\n",
    "sigma_b = 0.2\n",
    "\n",
    "w_0 = np.hstack([mu_b, mu_w])\n",
    "V_0 = np.diag([sigma_b, sigma_w])**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\\\(V_0^{-1}\\\\) is the inverse of the prior's variance. It shows up a few times, so I'll\n",
    "compute it once. It doesn't look like I can use `np.linalg.solve` on it, so I'll use\n",
    "`np.linalg.inv`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "V0_inv = np.linalg.inv(V_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\\\(\\Phi\\\\) is the augmented input matrix. In this case, it's the `x` values of the [observations](2018-01-08-bayesian-linreg-sample.ipynb), with the column of 1s I add to deal with the bias term. So, from [the last post](2018-01-08-bayesian-linreg-sample.ipynb), I had `x` as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_in = 2 * np.random.rand(21, 1) - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Then \\\\(\\Phi\\\\) is\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "Phi_X_in = np.hstack((\n",
    "    np.ones((X_in.shape[0], 1)),  # pad with 1s for the bias term\n",
    "    X_in\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\\\(\\textbf y\\\\) is also from [the last post]({% post_url 2018-01-08-bayesian-linreg-sample %}). It's the vector containing all the observations.\n",
    " The code used there was"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "true_w = np.array([[2, 0.3]]).T\n",
    "\n",
    "def f(x):\n",
    "    x_bias = np.hstack((\n",
    "        np.ones((x.shape[0], 1)),  # pad with 1s for the bias term\n",
    "        x\n",
    "    ))\n",
    "\n",
    "    return x_bias @ true_w\n",
    "\n",
    "true_sigma_y = 0.1\n",
    "noise = true_sigma_y * np.random.randn(X_in.shape[0], 1)\n",
    "\n",
    "y = f(X_in) + noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "But since I already have \\\\(\\Phi\\\\), I'll skip the function and just use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = Phi_X_in @ true_w + noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\\\\(\\sigma_y\\\\) is my guess of `true_sigma_y`. On a real dataset, I might not know the true \\\\(\\sigma_y\\\\), so I keep separate `true_sigma_y` and `sigma_y` constants that I can use to explore what happens if my guess is off. I'll start with imagining I know it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_y = true_sigma_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The rest is a matter of copying the equation over correctly and hoping I got it right!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "V_n = sigma_y**2 * np.linalg.inv(sigma_y**2 * V0_inv + (Phi_X_in.T @ Phi_X_in))\n",
    "w_n = V_n @ V0_inv @ w_0 + 1 / (sigma_y**2) * V_n @ Phi_X_in.T @ y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Complete code\n",
    "\n",
    "Putting it all together, I get"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the prior\n",
    "mu_w = 0\n",
    "mu_b = 0\n",
    "\n",
    "sigma_w = 0.2\n",
    "sigma_b = 0.2\n",
    "\n",
    "w_0 = np.hstack([mu_b, mu_w])[:, None]\n",
    "V_0 = np.diag([sigma_b, sigma_w])**2\n",
    "\n",
    "# Get observations\n",
    "true_sigma_y = 0.1\n",
    "true_w = np.array([[2, 0.3]]).T\n",
    "\n",
    "X_in = 2 * np.random.rand(11, 1) - 1\n",
    "\n",
    "Phi_X_in = np.hstack((\n",
    "    np.ones((X_in.shape[0], 1)),  # pad with 1s for the bias term\n",
    "    X_in\n",
    "))\n",
    "\n",
    "true_sigma_y = 0.05\n",
    "noise = true_sigma_y * np.random.randn(X_in.shape[0], 1)\n",
    "\n",
    "y = Phi_X_in @ true_w + noise\n",
    "\n",
    "# Compute the posterior\n",
    "sigma_y = true_sigma_y  # I'm going to guess the noise correctly\n",
    "\n",
    "V0_inv = np.linalg.inv(V_0)\n",
    "V_n = sigma_y**2 * np.linalg.inv(sigma_y**2 * V0_inv + (Phi_X_in.T @ Phi_X_in))\n",
    "w_n = V_n @ V0_inv @ w_0 + 1 / (sigma_y**2) * V_n @ Phi_X_in.T @ y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sweet! Now it seems like after doing all that code and math, I should be rewarded with pretty graphs!\n",
    "\n",
    "## Plotting\n",
    "\n",
    "In this post, I'll just show one graph. `w_n` is the mean guess of the weights, so I can plot that function. I can also compare it to the weights from\n",
    "least squares and the true weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f7007b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGoCAYAAABL+58oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlYl2Xe///nxQ6yKe6I4K6IKAqKKSZZaWWajmGLtm/O\nzD3TfXfbbOnM1DS/xvw290wzYU5NjWgl7fs6UWoCgrjjvnwQUcSFXbbP5/z9cRmj5QKGsvh6HIdH\nwudazgvLXsd5vc/zbRljEBEREWlL3Jp7ACIiIiJNTQFHRERE2hwFHBEREWlzFHBERESkzVHAERER\nkTZHAUdERETaHAUcERERaXMUcERERKTNUcARERGRNsejuW7csWNHExER0Vy3FxERkVZo7dq1R4wx\nnc53XLMFnIiICLKzs5vr9iIiItIKWZblaMhxekUlIiIibY4CjoiIiLQ5CjgiIiLS5jRbDY6IiEhL\nU1tbS35+PlVVVc09lMuej48PPXr0wNPT84LOV8ARERE5KT8/n4CAACIiIrAsq7mHc9kyxnD06FHy\n8/Pp1avXBV1Dr6hEREROqqqqIiQkROGmmVmWRUhIyA+aSVPAEREROYXCTcvwQ/8cFHBERESkzVHA\nERERuQALFiwgLS3ttO+lpaWxYMGCJr3Pvn37iIqKatJrXoj169fz0Ucf1X/93nvv8dRTTzXjiM5N\nAUdEROQCxMXFkZSUVB9y0tLSSEpKIi4urplHdn51dXWNPue7AWfKlCn88pe/bMphNSkFHBERkQuQ\nmJhIamoqSUlJzJ8/n6SkJFJTU0lMTPxB133mmWeIiooiKiqK//u//wPsQHL77bczaNAgZsyYQWVl\nJQC//OUviYyMJDo6mv/93/8FoKioiB/96EfExcURFxfHN998A8Dvfvc7Zs+ezZgxY5g9ezbx8fFs\n2bKl/r7jx48nOzubNWvWMHr0aGJiYrjiiivYvn07NTU1zJ8/n+XLlzNs2DCWL1/Oyy+/zE9/+lPA\nnmW66qqriI6OZsKECeTl5QFw11138bOf/YwrrriC3r1788YbbwBw8OBBxo0bx7Bhw4iKimLlypU/\n6Gd2RsaYZvk1YsQIIyIi0pLk5uY2+px58+YZwMybN+8H3z87O9tERUWZ8vJyU1ZWZiIjI01OTo4B\nzKpVq4wxxtx9993m6aefNkeOHDH9+/c3LpfLGGPM8ePHjTHG3HrrrWblypXGGGMcDocZOHCgMcaY\n3/72t2b48OGmsrLSGGPMM888Y+bPn2+MMaagoMD079/fGGNMSUmJqa2tNcYY8/nnn5vp06cbY4x5\n6aWXzE9+8pP6sZ769eTJk83LL79sjDHmxRdfNFOnTjXGGHPnnXeaGTNmGKfTabZs2WL69OljjDFm\n4cKF5g9/+IMxxpi6ujpTWlp6xp/Hmf48gGzTgJyhGRwREZELlJaWRnJyMvPmzSM5Ofl7NTmNtWrV\nKqZNm0a7du3w9/dn+vTprFy5krCwMMaMGQPArFmzWLVqFUFBQfj4+HDvvffy1ltv4efnB8AXX3zB\nT3/6U4YNG8aUKVMoLS2lvLwcsF8r+fr6ApCUlFQ/o5KamsqMGTMAKCkp4eabbyYqKor//u//Pm2W\n52zS09O57bbbAJg9ezarVq2q/+ymm27Czc2NyMhICgsLAfv13ksvvcTvfvc7Nm3aREBAwA/6uZ2J\nAo6IiMgF+LbmJjU1lccff7z+ddUPDTln8t0l05Zl4eHhwZo1a5gxYwYffPABkyZNAsDlcpGRkcH6\n9etZv349Bw4cwN/fH4B27drVXyM0NJSQkBA2btzI8uXLmTlzJgDz5s0jMTGRzZs38/777//gXZ29\nvb3rf29PwMC4ceNYsWIFoaGh3HXXXSxZsuQH3eNMFHBEREQuQFZW1mk1N9/W5GRlZV3wNRMSEnjn\nnXeorKykoqKCt99+m4SEBPLy8khPTwfglVdeYezYsZSXl1NSUsL111/Pn//8ZzZs2ADAtddey7PP\nPlt/zfXr15/1fjNnzmTBggWUlJQQHR0N2DM4oaGhALz88sv1xwYEBFBWVnbG61xxxRW89tprACxb\ntoyEhIRzPqfD4aBLly7cf//93HfffeTk5JznJ9N4CjgiIiIX4NFHH/1eQXFiYiKPPvroBV9z+PDh\n3HXXXYwcOZJRo0Zx33330b59ewYMGMDf//53Bg0axPHjx5kzZw5lZWVMnjyZ6Ohoxo4dyzPPPAPA\nX//6V7Kzs4mOjiYyMpJFixad9X4zZszgtddeIykp6bTn+tWvfkVMTMxpq60SExPJzc2tLzI+1bPP\nPstLL71EdHQ0KSkp/OUvfznnc3711VcMHTqUmJgYli9fzs9//vML+XGdk/XtdNGlFhsba7Kzs5vl\n3iIiImeydetWBg0a1NzDkJPO9OdhWdZaY0zs+c7VDI6IiIi0OQo4IiIi0uYo4IiIiEibo4AjIiIi\nbY4CjoiIiLQ5CjgiIiLS5ijgiIiItBDFxcU899xzzT2M7ykoKKhv5XAu3+6Y/F3vvPMOubm5TT2s\nc1LAERERaSHOFXBO3XTvUuvevXt936oLoYAjIiLSgrz88vd/fduJobb2zJ9/2xmhsvL7n53PL3/5\nS3bv3s2wYcOYO3cuX331FQkJCUyZMoXIyEj27dtHVFRU/fELFy7kd7/7HQC7d+9m0qRJjBgxgoSE\nBLZt2/a96w8ZMoTi4mKMMYSEhNT3gLrjjjv4/PPPcTqdzJ07l7i4OKKjo3n++ecBTrtvZWUlSUlJ\nREZGMm3aNEaNGsWpG/f+5je/YejQocTHx1NYWMjq1at57733mDt3LsOGDWP37t389a9/JTIykujo\naG655Zbz/2AugMdFuaqIiIg02lNPPcXmzZvr+0d99dVX5OTksHnzZnr16sW+ffvOeu4DDzzAokWL\n6NevH5mZmfz4xz/myy+/PO2YMWPG8M033xAeHk7v3r1ZuXIld9xxB+np6SQnJ/Piiy8SFBREVlYW\n1dXVjBkzhmuvvfa0Zp/PPfcc7du3Jzc3l82bNzNs2LD6zyoqKoiPj+fJJ5/k0Ucf5R//+AePPfYY\nU6ZMYfLkyfWvuZ566in27t2Lt7c3xcXFTfgT/A8FHBERkbO4666zf+bpee7P/fzO/XlDjRw5kl69\nep3zmPLyclavXs3NN99c/73q6urvHZeQkMCKFSsIDw9nzpw5LF68mAMHDtC+fXvatWvHZ599xsaN\nG+tfR5WUlLBz50769+9ff41Vq1bV946Kioqqb9IJ4OXlxeTJkwEYMWIEn3/++RnHGx0dze23385N\nN93ETTfd1MCfROMo4IiIiLRg7dq1q/+9h4cHLper/uuqqioAXC4XwcHB5+wcDjBu3Dj+/ve/k5eX\nx5NPPsnbb7/NG2+8Ud/92xjDs88+y8SJE08771wzR6fy9PSsn+1xd3c/a93Qhx9+yIoVK3j//fd5\n8skn2bRpEx4eTRtJVIMjIiLSQgQEBFBWVnbWz7t06cLhw4c5evQo1dXVfPDBBwAEBgbSq1cvXn/9\ndcAOKhs2bPje+WFhYRw5coSdO3fSu3dvxo4dy8KFCxk3bhwAEydOJDk5mdraWgB27NhBRUXFadcY\nM2YMqampAOTm5rJp06ZGPZfL5WL//v0kJibypz/9iZKSEsrLy897jcZSwBEREWkhQkJCGDNmDFFR\nUcydO/d7n3t6ejJ//nxGjhzJNddcw8CBA+s/W7ZsGS+++CJDhw5l8ODBvPvuu2e8x6hRo+pfOSUk\nJHDgwAHGjh0LwH333UdkZCTDhw8nKiqKBx988HuzMD/+8Y8pKioiMjKSxx57jMGDBxMUFHTO57rl\nllt4+umniYmJYefOncyaNYshQ4YQExPDz372M4KDgxv1c2oIyxjT5BdtiNjYWHNq1bWIiEhz27p1\nK4MGDWruYbRoTqeT2tpafHx82L17N1dffTXbt2/Hy8urye91pj8Py7LWGmNiz3euanBERESkwSor\nK0lMTKS2thZjDM8999xFCTc/lAKOiIiINFhAQACt4Q3MeWtwLMsKsywrzbKsXMuytliW9fNzHBtn\nWVadZVnn389ZRERE5CJpyAxOHfCIMSbHsqwAYK1lWZ8bY07bc9myLHfgT8BnF2GcIiIiIg123hkc\nY8xBY0zOyd+XAVuB0DMc+l/Am8DhJh2hiIiISCM1apm4ZVkRQAyQ+Z3vhwLTgOTznP+AZVnZlmVl\nFxUVNW6kIiIiIg3U4IBjWZY/9gzNw8aY0u98/H/AL4wxru+f+R/GmMXGmFhjTGynTp0aP1oREZE2\nzt3dnWHDhjF06FCGDx/O6tWrL8l977vvvkve8ftiatAqKsuyPLHDzTJjzFtnOCQWeO3k9swdgest\ny6ozxrzTZCMVERG5DPj6+ta3XPj000/51a9+xddff33R7/vCCy9c9HtcSucNOJadWl4EthpjnjnT\nMcaYXqcc/zLwgcKNiIi0Zr9/fwu5Bd99YfHDRHYP5Lc3Dm7w8aWlpbRv3x6wG2pOnTqV48ePU1tb\nyx/+8AemTp3K/Pnz6dChAw8//DAAv/nNb+jcuTM///nPefrpp0lNTaW6uppp06bx+9//noqKCpKS\nksjPz8fpdDJv3jxmzpzJ+PHjWbhwIbGxscyZM4esrCxOnDjBjBkz+P3vfw9AREQEd955J++//z61\ntbW8/vrrp+2m3JI0ZAZnDDAb2GRZ1rddvH4N9AQwxiy6SGMTERG57Jw4cYJhw4ZRVVXFwYMH+fLL\nLwHw8fHh7bffJjAwkCNHjhAfH8+UKVO45557mD59Og8//DAul4vXXnuNNWvW8Nlnn7Fz507WrFmD\nMYYpU6awYsUKioqK6N69Ox9++CFgdwz/rieffJIOHTrgdDqZMGECGzdurO8a3rFjR3JycnjuuedY\nuHBhi535OW/AMcasAqyGXtAYc9cPGZCIiEhL0JiZlqZ06iuq9PR07rjjDjZv3owxhl//+tesWLEC\nNzc3Dhw4QGFhIREREYSEhLBu3ToKCwuJiYkhJCSEzz77jM8++4yYmBjAngHauXMnCQkJPPLII/zi\nF79g8uTJ9Z3ET5WamsrixYupq6vj4MGD5Obm1gec6dOnAzBixAjeeutMVSstg3YyFhERaaFGjx7N\nkSNHKCoq4qOPPqKoqIi1a9fi6elJREQEVVVVgF0g/PLLL3Po0CHuuecewO4o/qtf/YoHH3zwe9fN\nycnho48+4rHHHmPChAnMnz+//rO9e/eycOFCsrKyaN++PXfddVf9fQC8vb0Buxj6u404WxJ1ExcR\nEWmhtm3bhtPpJCQkhJKSEjp37oynpydpaWk4HI7646ZNm8Ynn3xCVlYWEydOBGDixIn885//pLy8\nHIADBw5w+PBhCgoK8PPzY9asWcydO5ecnJzT7llaWkq7du0ICgqisLCQjz/++NI9cBPSDI6IiEgL\n8m0NDtizMP/6179wd3fn9ttv58Ybb2TIkCHExsaeVtzr5eVFYmIiwcHBuLu7A3DttdeydetWRo8e\nDYC/vz9Lly5l165dzJ07Fzc3Nzw9PUlOPn0Lu6FDhxITE8PAgQMJCwtjzJgxl+jJm5ZljGmWG8fG\nxprW0KxLREQuH1u3bmXQoEHNPYxGc7lcDB8+nNdff51+/fo193CazJn+PCzLWmuMiT3fuXpFJSIi\n0orl5ubSt29fJkyY0KbCzQ+lV1QiIiKtWGRkJHv27GnuYbQ4msERERGRNkcBR0RERNocBRwRERFp\ncxRwREREpM1RwBEREWlB/P39z/n5H//4x0s0ktZNAUdERKQVUcBpGC0TFxEROZOPfwmHNjXtNbsO\ngeueatChBw8eZObMmZSWllJXV0dycjIffvhh/U7HgwcPZvHixSQlJZGfn4/T6WTevHnMnDmTTz75\nhIcffhg/Pz/Gjh3Lnj17+OCDD5r2WVo4BRwREZEW6JVXXmHixIn85je/wel0UllZSUJCAn/729/q\nu42/+eabdO/enQ8//BCAkpISqqqquP/++/nyyy/p27cvM2fObM7HaDYKOCIiImfSwJmWiyUuLo57\n7rmH2tpabrrppvr+VKcaMmQIjzzyCL/4xS+YPHkyCQkJrF+/nl69etXvajxr1iwWL158qYff7FSD\nIyIi0gKNGzeOFStWEBoayl133cWSJUu+d0z//v3JyclhyJAhPPbYYzz++OPNMNKWSQFHRESkBXI4\nHHTp0oX777+f++67j5ycHAA8PT2pra0FoKCgAD8/P2bNmsXcuXPJyclh4MCB7Nu3j927dwPw6quv\nNtszNCe9ohIREWmBvvrqK55++mk8PT3x9/evn8F54IEHiI6OZvjw4dxxxx3MnTsXNzc3PD09SU5O\nxsfHh8WLF3PDDTfg5+dHQkICZWVlzfw0l55ljGmWG8fGxprs7OxmubeIiMiZbN26lUGDBjX3MJrU\nV199xcKFC1vlKqoz/XlYlrXWGBN7vnP1ikpERETaHL2iEhERacPGjx/P+PHjm3sYl5xmcERERE7R\nXKUbcrof+ueggCMiInKSj48PR48eVchpZsYYjh49io+PzwVfQ6+oRERETurRowf5+fkUFRU191Au\nez4+PvTo0eOCz1fAEREROcnT05NevXo19zCkCegVlYiIiLQ5CjgiIiLS5ijgiIiISJujgCMiIiJt\njgKOiIiItDkKOCIiItLmKOCIiIhIm6OAIyIiIm2OAo6IiIg0ifXroaV0uVDAERERkQtWWfmfUFNS\nAvv3N+94vqVWDSIiItJoZWWwejVkZ8O0aRAZCePGgWU198hsCjgiIiLSYMXF8M03sG4duFwwZAh0\n6WJ/1lLCDSjgiIiISAMZA8uWwbFjMGwYjB0L7ds396jOTAFHREREzqqoCDIzYeJE8PSEKVMgMBCC\ngpp7ZOemgCMiIiLfc/AgrFwJubng5QXR0dCzJ4SFNffIGkYBR0REROrV1MAbb8COHeDjA1deCaNG\ngZ9fc4+scRRwRERELnPG2MXD7dvbr6Hc3GDCBIiLs0NOa6SAIyIicpkyBnbvhhUroKAAHn4Y/P3h\nllsu4GLV5bDhVSg7BBPmNflYG0sBR0RE5DJjDGzf/p9gExgI114L3t4XcLHjDlizGHJSoLoEwkaB\nywlu7k0+7sZQwBEREbnMHDsGy5dDcLC9KmroUHBvTB4xBvatgsxFsP0jwILIKTBqDoSNbBEb4ijg\niIiItHFOJ2zeDIWF9kxNSAjceae9KsqtMU2baqtg0+uQ+TwUbgLfDjDmYYi7D4JCL9r4L4QCjoiI\nSBtVVwcbNsCqVXD8OHTtan/PwwMiIhpxodICyHoB1r4MlUehcyTc+FeITgJP34s0+h9GAUdERKQN\nysuzl3uXlkJoKEyaBP37N/Lt0f4syEyG3HftupoB10P8QxCR0CJeQ52LAo6IiEgbUV0NFRXQoYP9\nGqpTJ5g6FXr3bkQeqauB3Hfs+poDa8E7EEY+CCPvhw69Lur4m5ICjoiISCtXVWW3U8jIsIPNvfdC\nu3Ywe3YjLlJeBGtfgqwXofwQhPSF6xfC0FvB2/+ijf1iUcARERFppSoq7FCzZo09ezNgACQ09u3R\nwQ120fCm18FZA30mwNS/2f9sVAVyy6KAIyIi0kpt2WIXEEdG2sGma9cGnuisg+0fQsYiyFsNnn4w\n/A77VVSn/hc0ll2Hy0hJd3CsspZnb425oGs0JQUcERGRVqKkBL75Brp3h2HDICYGevWya20apPIY\nrEuBNf+Akv0Q3BOu/QPEzAbf4EaPp9bp4vPcQlLSHaTvOYqXuxuTo7vhdBnc3Zq3CFkBR0REpIU7\ndsyeqdmwwf7a9+TKbE/PBoabw9vsouGNy6G20l4FNen/s1dFXcCOw4WlVby6Jo9X1+RRWFpNaLAv\nj04aQFJsGB39L2Q75KangCMiItKCff21/cvNDUaMgDFjICioASe6XLDrc8hIhj1p4O4N0TfDqIeg\n65BGj8MYQ8aeY6Rk7OPTLYU4XYYr+3fiyZvCSRzYudlnbL5LAUdERKSFOXTIDjG+vnZdTXw8jB4N\nAQENOLmqFNa/Amueh2N7IKAbXPUYjLgb2nVs9FjKqmp5K+cAKRkOdh0uJ8jXk3vGRHD7qHAiOrZr\n/MNdIgo4IiIiLcSBA3YDzO3bITERrrzSXhk1YEADTj62BzIXw7qlUFMGPeIg8TcQORXcPRs9lm2H\nSklJd/D2ugNU1jgZEhrEghnRTBnaHR/P5m2k2RAKOCIiIs3M4bCDze7d9qzNVVfByJENONEY2POV\nXV+z41O7nmbwNLvpZY8RjR5HTZ2LT7YcYmm6gzX7juHl4caUod2ZHR/O0LDGFyE3JwUcERGRZvbN\nN/9phBkbC15e5zmhptIuGM58Hoq2gl9HGDcXYu+BwG6Nvn9B8Qleyczjtaw8jpTX0LODH7++fiA3\njwijfbvzDaZlUsARERG5hIyBHTvsVVHTptltFSZPtmduPM/3Jql4P2T9A9b+C6qK7WLhqc9B1I/A\n06dR43C5DN/sPkJKuoMvthZigKsGdGbW6HCu7NcJtxZWNNxYCjgiIiKXgMsFW7fCypV2EXH79lBW\nZgecwMBznGgM5GXYTS+3fgAYGDgZ4udAz9GNbnpZUlnLGzn5LMtwsOdIBR3aefHglX24bWRPwjr4\n/aBnbEkUcERERC4ylwsWL7aDTceO9szNkCHn6YRQVw2b37Traw5uAJ8gGP0Tu+llcM9Gj2HzgRKW\nZjh4Z/0BqmpdDO8ZzJ9nDuW6qG6tomi4sRRwRERELgKnE3btsldAublBVJTdTmHQoPMEm7JCyH4R\nsv8JFUXQaSBM/jNEzwSvxi3Lrqp18tGmg6RkOFiXV4yPpxs3DQtlVnw4UaEN2Uyn9VLAERERaUK1\ntbBunV1jU1oKDz4I3brB2LHnOfHAWrtoePNb4KqFfhMh/iHondjo11D7j1WyLDOP1Oz9HKuooXfH\ndsyfHMmPRvQgyLfxS8ZbIwUcERGRJlBbC1lZsHo1lJdDz54wZcp5GmA6a2Hre3bTy/w14BUAcffC\nyAcgpE+j7u9yGb7eWcTSdAdfbj+MBVwT2YXZ8RFc0Sek1RcNN5YCjoiIyA9gjD3B4nLZBcTdusGM\nGRAefo6Jl4qjkPMyrHkBygqgfS+Y9BQMux18zlVx/H3HK2p4fe1+lmbkkXesko7+3vw0sS+3juxJ\n92DfH/x8rZUCjoiIyAWorISMDNi7F+65B7y94Sc/AX//c5xUuMXuDbXpdairgt7j7fqaftc0uunl\nhv3FLEl38P7GAmrqXIyM6MD/ThzApMFd8fI4V5HP5UEBR0REpBHKy+3XUNnZUFMDkZFQVWXvY3PG\ncONywvaP7dVQ+1aChy8MvcVuetl5UKPuXVXr5L0NBSzNcLAxvwQ/L3duHtGD2aPDGdi1cTM/bZ0C\njoiISAMVFMA//2mvkBoyxC4c7tz5LAefKLb7Qq1ZDMUOCOwBV/8eht8Bfh0add99RypYlukgNTuf\nkhO19Ovsz+NTBzMtJpQAn8ujaLixzhtwLMsKA5YAXQADLDbG/OU7x9wO/AKwgDJgjjFmQ9MPV0RE\n5NI6fhyOHoW+fe2C4VGjYPhwCAk5ywlHdtqroda/ArUV9mZ81zxub87n3vB5BafLkLbtMEsyHKzY\nUYSHm8XEwV2ZFR9OfO8OWI1cWXW5achPug54xBiTY1lWALDWsqzPjTG5pxyzF7jSGHPcsqzrgMXA\nqIswXhERaSYLFiwgLi6OxMTE+u+lpaWRlZXFo48+2owjuziKiuyl3ps2QUAA/Pzn9v4111xzhoNd\nLtj9pb3b8K4vwN0LombAqAeh+7BG3fdIeTXLs/bzSmYeB4pP0CXQm4ev7setI3vSJbBx7RguZ+cN\nOMaYg8DBk78vsyxrKxAK5J5yzOpTTskAejTxOEVEpJnFxcWRlJREamoqiYmJpKWl1X/dlhw5Amlp\nkJsLHh72jM0VV5xlc77qctjwqj1jc3Qn+HeB8b+G2LvB/2zvrr7PGENO3nFS0h18tOkQNU4Xo3uH\n8NgNg7g6sgue7ioabqxG1eBYlhUBxACZ5zjsXuDjs5z/APAAQM+ejd9mWkREmk9iYiKpqakkJSUx\nZ84ckpOT68NOW+By2SGmrMzegXjsWIiPh3Zn2jz4+D5Y8w/ISYHqEugeA9MWw+Bp4NHw7tuVNXW8\nu76AlHQHuQdLCfD24LZRPZkV35O+nQOa7NkuR5YxpmEHWpY/8DXwpDHmrbMckwg8B4w1xhw91/Vi\nY2NNdnZ2I4crIiLNbf78+TzxxBPMmzePxx9/vLmH84Pl5cGKFXZNzXXX2fva1NTYy75PYwzsW2Wv\nhtr+EWBB5FS76WWPuEbtNry7qJyUdAdv5uRTVlXHwK4BzB4dzk3DQmnnrfU/52JZ1lpjTOz5jmvQ\nT9GyLE/gTWDZOcJNNPACcN35wo2IiLROaWlpJCcnM2/ePJKTk0lMTGyVMzjG2PvXrFgB+/bZszR9\n+9qfWdZ3wk3tCXvfmsznoXAz+HaAsf8NsfdCUGiD71nndPHF1kJSMhx8s+sonu4W10V1447R4YwI\nb6+i4SbWkFVUFvAisNUY88xZjukJvAXMNsbsaNohiohIS3Bqzc23webUr1uTr7+Gr76CwEB71mb4\ncPD87mrr0gLIegGyX4ITx6DzYJjyLAy5GTwbvkPw4dIqXjtZNHyotIruQT7MnTiApNgwOgV8d5pI\nmkpDZnDGALOBTZZlrT/5vV8DPQGMMYuA+UAI8NzJBFrXkOkjERFpPbKysk4LM9/W5GRlZbX4gGMM\nbN0KHTva+9ZERdmb8g0bZhcSn3Zgfpa92/DW9+xN+gZcbze9jEho8GsoYwxr9h5jSYaDTzcfos5l\nSOjXkcenDuaqgZ3xUNHwRdfgGpymphocERG52Fwu2LzZ7hFVVAQjR8L115/hwLoayH3HDjYFOeAd\nBMNnQ9x90KFXg+9XXl3H2zn5pGQ42FFYTqCPBzfHhjErPpxeHc9UrSyN1aQ1OCIiIq3Nxo32a6hj\nx+xZmxkz7LYKpykvgux/QvaLUF4IIX3h+oUw9FbwPldTqdPtKCwjJd3BWzn5VNQ4iQoNZMGPorlx\naHd8vRoa9MB2AAAgAElEQVTXY0qahgKOiIi0GXV14O5uv0k6fBh8fOCWW2DAgO+8XTq4ATIWweY3\nwFkDfa+GUX+HPhPOsuHN99XUufgs9xAp6Q4y9x7Dy8ONydHdmB0fzrCwYBUNNzMFHBERafVqamDt\nWvjmG5g6Ffr1g/HjYcKEU4KNsw62fWAv885LB892dl+okQ9Cp/4NvtfBkhO8mpnHq1n7KSqrpkd7\nX3553UCSYsPo0K7he+DIxaWAIyIirVZVFWRlQXo6VFZCr17g52d/Vl88XHkMcpbYK6JK9kNwT7j2\nSYiZBb7BDbqPMYbVu4+Sku7g862FuIxhfP9OzB4dzpX9O+PuptmalkYBR0REWiVj4KWXoLDQnrEZ\nNw7Cwk454PBWe7Zmw3KoO2Gvgpr0FAy4DtwaVhdTcqKWt04WDe8pqqC9nyf3JfTi9pHh9AzxuzgP\nJk1CAUdERFqN8nL7VdTYsXatzdVX25v0de9+8gCXC3Z+Zje93PMVuHtDdJLd9LLrkAbfJ7eglJQM\nB++sO8CJWifDwoL5fzcP5Ybobvh4qmi4NVDAERG5BC63TtxNrbTUrq9ZuxacTnumpndve+YGgKpS\nWL/M3m34+F4I6AZXzYMRd0O7kAbdo7rOySebD7Ek3cFax3G8PdyYOqw7s+MjGNIj6OI9nFwUCjgi\nIpfA5dKJu6nV1sInn8D69fYrqehoSEiw+0YBcHQ3rFkM65ZBTRn0GAlXPWb3iHL/7tbEZ5Z/vJJX\nMvNYnrWfoxU1RIT48dgNg5gxogfBfioabq0UcERELoG23om7qVVV2Uu8PTzs5d4xMfZrqeBg7KSz\nO82erdnxKbh52F284x+C0BENur7LZVi56wgp6Q6+3FYIwIRBXZgdH87Yvh1xU9Fwq6eAIyJyiSQm\nJjJnzpz6TtwKN99XWGjvOrxrF/z85+DrC3fffXJrmppKyH7NDjZF28CvI4ybC3H3QkDXBl2/uLKG\nN9bmszTDwb6jlYS082LO+D7cOrInPdqraLgtUcAREblE2kon7ouhoMDu7L1tG3h52S0VvuVWuh+y\n/gFr/wVVxdA1Gm5KhsHTwdOnQdfflF/CkvR9vLehgOo6F7Hh7fnva/ozKaor3h4qGm6LFHBERC6B\nttSJu6kdOwaLF9uvpMaPh1GjwNfH2JvxZSTbm/MBDJwM8XOg5+gGNb2sqnXywcaDpGQ42LC/GD8v\nd340ogezRoUT2T3w4j6UNDsFHBGRS6A1d+JuasbA3r1w8CCMGQMdOth9ovr1A2/3atj8ph1sDm0E\nn2C44r8g7n4IDjv/xYG8o5Usy3SwPHs/xZW19OnUjt/dGMn0ET0I9GlY4bG0fuomLiIil4Qxdm3N\n119Dfj4EBcFPfwqenkDZIch6Eda+BBVF0GmgvXdN9EzwOn8XbqfL8PWOwyxJd/D1jiLcLItrI+2i\n4dF9QtQXqg1RN3EREWkxDh2Cd9+1Z22Cg2HyZBg2DDwK19pNL7e8Da466D8RRj0Evcc36DXUsYoa\nlmftZ1mmg/zjJ+gU4M1/XdWP20b2pGtQw+pzpG1SwBERkYvC5bL7Q/n72/2hnE67EWb04Frct78L\nLz8P+WvAK8BeCTXyAQjpc97rGmNYt7+YpekOPth0kJo6F/G9O/Cr6wZx7eAueLo3rBu4tG0KOCIi\n0qScTti4EVatgsBAuPNO+59z7jyKtfYlePZFKCuA9r1g0p9g2G3gc/6i3xM1Tt7bcIAl6Q62FJTi\n7+3BLXFhzIoPp3+XgEvwZNKaKOCIiEiTqKuDdevslgrFxdCtm73c2xzajJWZjLXxdXBW26+fJv8Z\n+l17coObc9tTVM7SjDzeWLuf0qo6BnQJ4ImbopgWE4q/t/43JmemfzNERC6iy6kHVXa23VYhLAxu\nuM5JX9fHWJmLYN9K8PC1Z2pGPQidB533WnVOF//edpilGQ5W7jyCh5vFpKiu3DE6griI9ioalvNS\nwBERuYjacg+q6mrIyrL7Qg0aZLdT6Na+mJ5HUrA+XgzFeRDYA67+PQy/A/w6nPeaRWXVLM/K45XM\nPApKqugW5MMj1/Rn5sgwOgeoaFgaTgFHROQiaos9qE6cgMxMyMiwe0aNHAmDOu3EO3MR4etfhdoK\n6HkFXPOEvTmf+7n/V2OMIdtxnCXpDj7ZfJBap2Fs347Mv3EwVw/qjIeKhuUCKOCIiFxkbakH1Zo1\n8MUXUFMDAwe4uKbXvwnZtQj+9gW4e0HUDPs1VPdh571WRXUdb687wNIMB9sOlRHg48Gs+HBmxYfT\np5P/JXgaacsUcERELrLW3oOqtNRuo+DlZTe/HNSnnAkdXyUw93nYvhP8u8D4X0Ps3eDf+bzX21lY\nxtIMB2/mHKC8uo7IboE8NX0IU4Z1x89L/1uSpqF/k0RELqLW3IOquNhe6r1uHVx1FYwZtJchBS8w\nZE8KbC2B7sNh+j8g8ibw8DrntWqdLj7PLSQl3UH6nqN4ubtxQ3Q3ZsWHM7xnsIqGpckp4IiIXESt\nsQfV0aOwcqW9l42F4eq+KxmxexF8/hFYbhA51W562SPuvLsNF5ZW8eqaPF5dk0dhaTWhwb48OmkA\nM2PDCPH3vkRPJJcj9aISEZHTLFkCBY4T3BD+OpGlz+N+ZDP4drBfQcXeC0Gh5zzfGEP6nqMszXDw\n6ZZCnC7Dlf07ccfocMYP6Iy7m2Zr5MKpF5WIiDRIQYH9KmrSJAg0B5jR/gV8D76MtecYdB4MU56F\nITeDp+85r1NWVctbOQdIyXCw63A5wX6e3Du2F7eP6kl4yPkbZoo0JQUcEZHL1P79sGIF7Nxp6O2V\nhfVmMuS9i59xwcAb7KaXEWPP+xpq26FSUtIdvL3uAJU1Tob2COLpGdHcOLQ7Pp7ul+hpRE6ngCMi\ncplxuWDZMti3u4YYr7d5JHARAaU5cCjIrq0ZeT+0jzjnNWrqXHyy5RAp6fvI2nccbw83bhzandnx\n4QwNC74kzyFyLgo4IiKXAWPsV1GhoeBWeZjR1S8x0+tFvGoKIaAfXL8Qht4K3ufef+ZA8Qlezczj\ntaw8jpTXEB7ix2+uH8SMET1o3+7cK6lELiUFHBGRNswY2LbNfhXFwfXc0f95fHe/QV9nDfS9GkY9\nB32uOmfTS5fL8M3uI6SkO/hiayEGmDCwM7PiwxnXrxNuKhqWFkgBR0SkDXK5YMsWWLWijpCiD5js\nsYhQ0jF728HwO+3dhjv2O+c1SipreSMnn2UZDvYcqaBDOy8evLIPt43sSVgHv0v0JCIXRgFHRKQN\nqi4+RtE7S5hlXiCA/Rj/njDySayYWeB77hqZzQdKSEl38O6GA1TVuhjeM5g/zxzK9UO64e2homFp\nHRRwRETagLo6WL8einK3Mil4Eb4bl3OV8wQmIgHi/4TVfxK4nT2cVNU6+WjTQVIyHKzLK8bX051p\nMaHcPiqcqNCgS/gkIk1DAUdEpBWrrYW12S4OrfiUIScWEctXGA8fe9+aUQ9hdY065/n7j1WyLDOP\n1Oz9HKuooXfHdsyfHMmPRvQgyNfzEj2FSNNTwBERaaUK80rZvHQZMTXPE89e6vy6Y+LnY424C9qF\nnPU8l8vw9c4ilqY7+HL7YSzgmsguzI6PYEzfEPWFkjZBAUdEpBU5cQJK9u6mq2MxndctY0JNGVWd\nRsKV8/AYNAXczz7rcryihtTs/SzLzCPvWCUd/b35aWJfbh3Zk+7B596lWKS1UcAREWkFKsoNOz5J\nIyB3EX1cn2HcPLCipsOoB/EJHXHOc9fvLyYl3cH7GwuoqXMxslcH5k4cwMTBXfHyOPvycJHWTAFH\nRKQFKztWQd4Hy+m853li2EaVeycqRjyK//h7IKDrWc+rqnXy3oYClmY42JhfQjsvd5JiezArPpyB\nXQMv4ROINA8FHBGRlqg4D9b8A7/sJQyuKabYN5qyK5IJGP0j8PA+62n7jlSwLNNBanY+JSdq6dfZ\nn8enDmZaTCgBPioalsuHAo6ISEthDKWb0yn/dzLdij/AssBt0I2UDZ5DcGT8WZteOl2GL7cdJiXD\nwYodRXi4WUyM6srs+HBG9eqgomG5LCngiIg0t9oqStPfxLl6Ee2rNuJJMI4e/0XEjPuxgsMIOMtp\nR8qrWZ61n1cy8zhQfIIugd7899X9uWVkGF0CfS7pI4i0NAo4IiLNpewQZL1Ibfo/Caw9QpE1kNw+\n/0fPyTOJaH/mVgjGGHLyjpOS7uCjTYeocbq4ok8Ij90wiKsju+Dpfv6i4QULFhAXF0diYmL999LS\n0sjKyuLRRx9tsscTaU4KOCIil1r+WirTkvHd+zaWy4kzfBJbgh6k/8TxdGp35tdJlTV1vLu+gJR0\nB7kHSwnw9uC2UT2ZFd+Tvp3PNsdzZnFxcSQlJZGamkpiYiJpaWn1X4u0FQo4IiKXgrMWk/su1SsW\n4VOUhTsB7O18P71n3o9PSB+GneW0XYfLWZrh4M21+ZRV1zGwawB/nDaEqcO60877wv4KT0xMJDU1\nlaSkJObMmUNycnJ92BFpKxRwREQupoojsPYl6tJfxOPEQSrozTdef8J/7G3ExAeC1/dPqXO6+GJr\nISkZDr7ZdRRPd4vrh3Rjdnw4I8LbN0nRcGJiInPmzOGJJ55g3rx5CjfS5ijgiIhcDIc2YzKSYdPr\nWM5qigMTWclfCB1/DeOGu+F5hhXbh0ureO1k0fCh0iq6B/kwd+IAkmLD6BRw9qXhFyItLY3k5GTm\nzZtHcnIyiYmJCjnSpijgiIg0FZcTtn+EyViE5VhFHb5U9b+NgKsfwj9wIFM8wf07Db2NMazZe4wl\nGQ4+3XyIOpchoV9HnrgpiqsGdsbdremXeJ9ac/NtsDn1a5G2QAFHROSHOlEM61IwaxZjFedR5hZG\nBo/j6DCbCfEdCOgM3120XV5dx9s5+aRkONhRWE6gjwd3XRHB7fHh9OrY7qIONysr67Qw821NTlZW\nlgKOtBmWMaZZbhwbG2uys7Ob5d4iIk2iaAdkLoINr0JtJQe9rmBFzUMc73wDCVd6MGgQuH1n1fb2\nQ2UszXDwVk4+FTVOokIDuSM+ghuHdsfXy/3M9xGRepZlrTXGxJ7vOM3giIg0hssFu/8NGcmw+9+4\n3LywhtyMFf8gR44MZZgX9O9/+qbDNXUuPss9xJJ0B2v2HsPLw43J0d24Y3QEQ3sEaadhkYtAAUdE\npCGqy2H9K7DmeTi6ixqvLqzx/A3ptXdzU1Qn+nWDId1OP+VgyQlezczj1az9FJVVE9bBl19dN5Cb\nY8Po0O4My6dEpMko4IiInMuxvbDmH7AuBapLKQ0YwVeeL7ChZioRfbxIGgfh4f853BjD6t1HSUl3\n8PnWQlzGkDigM7PjwxnXv9NFKRoWke9TwBER+S5jYN9KyFhkr4pyc8eKnIoZNYeU9+Lo0B3uToAe\nPf5zSsmJWt46WTS8p6iC9n6e3JfQi1mjwgnrcOa2CyJy8SjgiIh8q/YEbEyFzOfh8BZcviHs6/EI\nX5Tcy503dsfbG+6/H7xOebuUW1BKSsY+3llXwIlaJ8PCgnkmaSjXD+mGj6eKhkWaiwKOiEjJAch6\nAda+DCeO4ewYxZZef+OjvBlUH/BlyBCorQVvbzvcVNc5+XjTIVIyHKx1HMfH042pQ0OZFR/OkB5B\nzf00IoICjohcroyB/WsgMxly3wMMDLiessEP8X9vj4XjFkOHwtix0KGDfUr+8UpeycxjedZ+jlbU\nEBHix2M3DOLmEWEE+Z1ha2IRaTYKOCJyeamrgS1v28GmYB14B1E5dA77Ot1P5JgIAoBrKmDQIAgK\nApfL8PWOI6SkO/hyWyEAEwZ1YXZ8OGP7dsRNRcMiLZICjohcHsoPQ/Y/7V/lhRDSj5KxC/n3kVvZ\ntMEfHx/oG2e/goqPh+LKGv6xIp+lmQ4cRyvp6O/Fj8f35dZRPQkN9m3upxGR81DAEZG2rWC9vdvw\n5jfBWQN9r6E08iE+3HYV21e54e0NCQl2qPHygo35xaSkO3hvQwHVdS7iItrzP9f0Z1JUV7w9VDQs\n0loo4IhI2+Osg23v26uh8tLBsx0Mv5Oa4Q/i1a0fdccg/99w1VUwciTg7uSDjQdJyXCwYX8xfl7u\n/GhED2aNCieye2BzP42IXAD1ohKRtqPyGOT8C9a8AKX5EByOGfkg+9rP4quMIHx84NZb7UOdTjhQ\nXMmyTAfLs/dTXFlLn07tmB0fzvQRPQj0OXvR8IIFC4iLizutMWVaWhpZWVk8+uijF/spRS5r6kUl\nIpePwly7hcKG5VB3AiISMNctYAeTWLHKnQMHIDDQXhFV5zSs2HmYJekOvt5RhJtlMXFwF2bFhzO6\nd0iD+kLFxcWRlJRU35E7LS2t/msRaRkUcESkdXK5YOendtPLvV+Dhw9EJ8Goh6DLYNZkwscfQ/v2\ncOON0KNvNW+ty+d/FjrIP36CzgHe/Oyqftw6siddg3wadevExERSU1NJSkpizpw5JCcn14cdEWkZ\nFHBEpHWpKoV1S2HNYji+FwK6w4T5uGLuYtPuEAIqoDcQHQ3e3gZncDHLMh188NZBapwu4nt34FfX\nDeLawV3wdHe74GEkJiYyZ84cnnjiCebNm6dwI9LCKOCISOtwdLddNLx+GdSUQ9gomDAfZ/8b2bDZ\nk5UvwvHjdrDp1sPJe5sPsCTdwZaCUvy9PbhlZBiz48Pp1yWgSYaTlpZGcnIy8+bNIzk5mcTERIUc\nkRZEAUdEWi5jYPeXdrDZ+Rm4eUDUdPs1VOhwNm2Cz/8GpaXQvTsMji9n9eE8fvfH/ZRW1TGgSwB/\nuCmKm2JC8fduur/uTq25+TbYnPq1iDQ/BRwRaXlqKmDDa3awObId2nWCK38BsXdT490VNzf7L6+a\nGggMchESfZjP9zn449tH8HS3mBTVjdnx4cRFtG9Q0XBjZWVlnRZmvq3JycrKUsARaSG0TFxELooL\nWkpdnAdr/mEv9a4qgW5DYdQciJpOVZ03a9ZAejokJkKvQdW8lpXHq5l5FJRU0S3Ih9tG9mTmyDA6\nBzSuaFhEWo8mWyZuWVYYsAToAhhgsTHmL985xgL+AlwPVAJ3GWNyLmTgItI2NHgptTHgWG33htr2\nIWDBoBvt11A946k8YZGxAjIzoarK4N3jOIs3O/j63YPUOg0J/Try2ymDmTCwMx4/oGhYRNqWhryi\nqgMeMcbkWJYVAKy1LOtzY0zuKcdcB/Q7+WsUkHzynyJymTrvUuraKrt9QmYyHNoEPsFwxc8g7j4I\nDqu/zhtvwPbddVR2PsBGDwe7d5cR4OPB7PgIbo/vSZ9O/s30hCLSkp034BhjDgIHT/6+zLKsrUAo\ncGrAmQosMfb7rgzLsoIty+p28lwRuUydcSl16UHIfhGyX4LKI9BpENz4FxiSBF5+lJRA+id2f6iC\n8jI2ejj42BygYn8dkd0CeWr6EKYM646fl0oIReTsGvU3hGVZEUAMkPmdj0KB/ad8nX/ye6cFHMuy\nHgAeAOjZs2fjRioirc6pS6kz3vw7h7quo+uR1eByQv9JEP8Q9LoSLIvjx2HlJ7BuvYt9dYU8t3Uf\nGwuP4eXuxg3R3Zg9OpyYsOCLUjQsIm1PgwOOZVn+wJvAw8aY0gu5mTFmMbAY7CLjC7mGiLQOaWlp\n3HbLzXz1t/9icOkKuLmOsvwV7O91I2HTfw8hfQB7Q+J334HMDVVsd+axx8qjpKaa0GpffjFpIEmx\nPQjx927mpxGR1qZBAceyLE/scLPMGPPWGQ45AISd8nWPk98TkctRxRHcV/+ZfT8PwHvrM9ChN1y3\ngHWlYWSs28yjIX0oKYHAQEPG3qO8vNPBpupCjDFcOaATs+PDGT+gM+5umq0RkQvTkFVUFvAisNUY\n88xZDnsP+KllWa9hFxeXqP5G5DJ0aBNkLIJNrzPOWQ19rrJXQ/W9BtzcGAf0ibyel5bW8sn2AxwO\ncLD3WDnBvp7cl9CL20f1JDykXXM/hYi0AQ2ZwRkDzAY2WZa1/uT3fg30BDDGLAI+wl4ivgt7mfjd\nTT9UEWmRXE7Y/pEdbByrwNMPYm6HkQ9C54H1h+XlQepnpXy218Ee5wHqcDLEJ4iFNw9lcnQ3fDzd\nm/EhRKStacgqqlXAOeeJT66e+klTDUpEWoETxyEnxd6YryQPgsLgmsdh+B3g277+sJo6F+/mHGTh\nOw4KXcfxdHNjSkx37rwinKFhwc34ACLSlmmdpYg0TtEOyFwEG16F2koIHwMTn4QB14O7/VeKMfDN\nuhMsTc8j63geR8tr6B7oxy+vGMQtI3sQ7OfVzA8hIm2dAo6InJ/LBbu+sDfl2/0luHvBkJth1IN2\nO4WTnE7Dq2lH+NdqB7sqCzHA+H6duTshgoS+HXFT0bCIXCIKOCJydtVlsP5VWPM8HN0F/l0h8TEY\ncRf4d6o/rKSyln+tzOellQ6O11Xg6+bF9EF9+NkNPQnv6Nd84xeRy5YCjoh837G9sGYxrFsK1aUQ\nOgKmvwCRU8HjP6+XNuSV8M+VDj7ddoCqWhdhvu25Z1Q/7pvYFV8vFQ2LSPNRwBERmzGwd4VdX7P9\nY3Bzh8ibIH4O9PhP496qWifvrz/I4jQHO48V42G5MyM2lNmjwxncPagZH0BE5D8UcEQud7UnYONy\nyHweDueCXwgkPAJx90Jg9/rD9h+rZMnqPF5bs5+ymhoCrXZM7BTJnEk9GBrpiTooiEhLooAjcrkq\nyYesF2Dty/aS7y5DYOrfIWoGePoA4HIZvt5ZREq6g7Tth7GwCLO6ML1nOHddF0JEhKVgIyItkgKO\nyOXEGNifab+Gyn0PMPby7vg59nLvk2nleEUNqdn7WZqRx/7jlQR5e/NfiX2ZGdcTKn0JDW3exxAR\nOR8FHJHLQV01bHkbMpLh4HrwCYLRP4a4+6F9eP1h6/cXk5Lu4P0NBdQ4XXRz78CVngOYGtOV6de6\n2Qe1P8s9RERaEAUckbasrBDWvgRZL0LFYejYH274fxB9C3j7A3Cixsn7GwtYmuFgY34JPh7u9HEL\nY4BHOAnRASQkQOfOzfwcIiKNpIAj0hYVrLN7Q215C5w10O9ae1O+3leBmz0Ts+9IBUszHLy+Np+S\nE7X06ejPE1MHE9c5lG2bPElIgJCQZn4OEZELpIAj0lY462Db+3aw2Z8BXv72hnwjH4SOfe1DXIYv\ncwtJyXCwYkcR7m4WQ4K70s0Vzg0DOzB5tF2DM7BPMz6HiEgTUMARae0qj9krobJegNIDEBwOE/8I\nMbPsWhvgSHk1y7P280pmHgeKT9CpnTfXdu9P+2NhBNb6MGI0XHFF8z6GiEhTUsARaa0Kc+3VUBtT\noe4E9BoH1y+E/hPBzR1jDDmOY6SkO/ho0yFqnC7G9A1h3uRBnNjVha25bowcC6NHQ7t2zf0wIiJN\nSwFHpDVxOWHHp3bTy70rwMMHopNg1EPQZTAAlTV1vLMuj5QMB1sPlhLg7cGNkT3pURXOzRP96dED\nyiNg0kTw9W3ex7kQCxYsIC4ujsTExPrvpaWlkZWVxaOPPtqMIxORlkQBR6Q1qCqBdcvsppfH90Fg\nKEz4rV1j49cBgF2Hy1ma4eDNtfmUVdcxsGsA/5MwBO/C7hzc6UGtH5SW2pfz92+2J/nB4uLiSEpK\nIjU1lcTERNLS0uq/FhH5lgKOSEt2ZJcdata/AjXlEDbKDjaDbgR3T+qcLr7YfJAl6Q5W7z6Kp7vF\n9UO6ccfocLaubM+uLIuAAJg0CYYPBy+v89+ypUtMTCQ1NZWkpCTmzJlDcnJyfdgREfmWAo5IS2MM\n7P7Srq/Z+Rm4eULUj+xl3qHDAThcWsVrWft4JTOPQ6VVhAb7MnfiAOJCwogb4o1lgbMQBg6EYcPA\no439l56YmMicOXN44oknmDdvnsKNiHxPG/trT6QVq6mADa9C5v/f3p2HV1Udeh//rswEQgghYc5h\nnmcSSBjEIIqigIgCWtBqta+0T997O7x2ury9re1rr5372gapVRScKOU61FmJikggYZ4EwxBGCTMh\nIeNZ94910IiEBDiZdn6f5zkP5+y9z95rZZ9Dftl7DQvg2A5ongDjfgTJ90FMW6y1rN59nEVZeby1\n5TPK/ZZreiXwi6kDSCxP5OOVhtfzoVUU9OoFI0bUd4VqT2ZmJhkZGcybN4+MjAzS09MVckTkSxRw\nROrbqX2wZgGse8a1tWk/BKY9Dv2nQVgkBcVlvLRqL4uy8th55CyxzcL5+qgu3DXSx+kDzfloOXx8\n3I02PH069OhR3xWqXZXb3JwPNpVfi4iAAo5I/bAW8la621CfvAYY164mda5rZ2MMOz4rYHHWpyxb\nd4DC0goGdozl0dsHccvADkRHhlJRAY89B1FRMHOmux3VFGb2zs7O/lKYOd8mJzs7WwFHRD5nrLX1\ncuDk5GSbk5NTL8cWqTdlxbBlqRtt+MhmaBbnekKl3A+xnSgt9/P2ts94ZlUea/acICIshMmDOjAn\nzUf/dq1YuxY2bYJ774XwcNcrKiamaQQbEREAY8xaa21yddvpCo5IXThzGHL+DjlPQdExSOgLk/8E\nA2dARDSHT5/j+bd38Hz2fo4WlNC5dTN+fFMfZiR3Jjosguxs+OMLUFgIXbtCURHExkLLlvVdMRGR\nhkkBR6Q2HciBrAzY9pIbpK/3TW5Qvq7XYIGPdx1n0artvLP9CH5rSe+dyJw0H+N6JhASYjh1Cv74\nGJw7Bz17wtixkJRU35USEWn4FHBEgq28FLa/4oLNwRyIbAkjvgkjHoDW3Th9roxlH7tGw7uPFhIX\nHc4DY7vxtZFJdG4dTWEh5Oa6nlCxsTB8OPTrBx061HfFvE+jJIt4hwKOSLAUHnO3oLKfgLOfQevu\ncNNvYMidEBnD1kOnWbxsEy+tP8S5sgqGJrXi9zMGM2lge6LCQzlzBt58E9audW1qfvADNzDfhAn1\nXbGmQ6Mki3iHAo7I1Tq8CVY/Dpv/ARUl0H08TPn/0GMCJX7LG5s/Y1HWZtbmnSQqPISpgzsyJ83H\ngEv9ohUAAB/QSURBVI5upu8zZ+Ddt2D9ete5atAgGDPGG6MONzYaJVnEOxRwxLNq9XaDv8J17149\n33X3Do+GobPdaMMJvTlwsojn3t7Ji9n7OV5YStc2zZl3Sz9uH9aJ2OhwwIUZY6C4GDZsgKFDYfRo\niIu7uqLJ1dEoySLeoIAjnlUrtxvOnYR1i2DN3+D0PohNgusfhmFz8Ee2YkXuMRa9ns3yT/IBmNC3\nLXPSfIzu3oaQENeXOz8fVqxwu5s+3Q3Q9/3vN86Zvb1IoySLeIMCjnhWUG83HN3hrtZsfAHKisA3\nBib+CnpP4lSJn39kH2Dx6g3kHS+iTYsIvnVtD+4amUSHVl+klkOHXLDZvt3dfhox4ourOAo3DYNG\nSRbxDgUc8bSrut3g90Puu7A6w01+GRoJA+9wt6HaD2LTgVMsWraVVzYeoqTcT0qXOL5/Q29u7N+O\niLCQL+1q3Tp45RU36vC4cTByJERHB7myctU0SrKId2gkY/G083+RX9YVnJIC2PCcazh8Yhe0aOdG\nGk6+l+KIOP616TCLsvLYuP8U0RGhTBvakdmpPvq2/2LUPWth71432nCnTnD2rGtEnJLiQo6IiFwZ\njWQsTd5l3244sdu1rVm/GErOQMdkmP536DuFfafLWfxBHkty1nGqqIweiS34+ZT+TBvWkZZR4Z/v\nwlo3hs2HH8L+/W5+qFmzoEULN0ifiIjUDQUc8awa3W6wFvZ84OaG2vkmhIS6WbxHPkhFh+G8vyOf\nRYs28MHOo4QYw8T+bZmd6iOtWzzmggmgcnPhvffg8GE3QN/NN7ueUSIiUvd0i0qaptIi2LzE3YbK\n3wbR8ZB8HyR/g+MhrVmSc4BnV+dx4OQ5EmMiuWtkErNSkmgX++X7S36/+zckBLKyYM0ad6Vm0CAI\nDa2HeomIeJxuUYlczOkD7jbUuqddl++2A2HqX7ADprP+cDGL38jjX5s2UFrhJ61bPD+Z1Jfr+7Ul\nPPTLjYYrKmDzZtcraswYd6UmJcX1jAoJqeLYIiJSZxRwxPushf2r3dxQ218FLPS5GUY+yLn2qby8\n8RCLMnLYeugMLSLDuHNEZ2an+ujZNuYruyovd4PyffQRnDoF7dpBTGAzXbEREWk4FHDEu8pLYMsy\nN37N4Q0QFQtp34KUB9hdHs/irH0sXfgeZ4rL6dMuhl/eOoBpQzvSPLLqr8WLL8Knn7qeUZMmuRm+\nL2iKIyIiDYACjnhPwRHIedI9CvOhTS+4+XeUD5jJe7sLWbwsjxWfbiE81HDTgPbMSfOR7Iv7SqNh\ngJISN/nlsGGue/eoUZCWBl27KtiIiDRkCjjiHQfXuas1W5aBvwx63gAjH+Ro4mhezNnPc39cw6HT\nxbSPjeIHN/RiRkpnEmMuPijNuXOuwXBWlnvevDkMHuyCjYiINHwKONK4VZS5djWr57t2NhEtIPk+\n7IgHyC6IZ1FWHm9uWU5ZhWVszzb8bEp/ruuTSFjoxVsCWwvLl7twU1LixrEZOxY6dqzjeomIyFVR\nwJHGqegErF0I2U/AmYMQ1wUmPsLZfrN4aXsBixfn8clnO4mJCmNOahdmpybRLaFFlbsrKYHISHfb\n6cgR17Zm7Fho27bOaiQiIkGkgCONy5Gt7mrNpiVQXgxdx8Gk3/Jp7CgWrznAP9/M5mxJOf07tOS/\npg9k8uAOREdU/TE/dQpWroSNG+Fb34JWrdzIw+rqLSLSuCngSMPnr4Cdb7lJL/d8CGFRMGgmZSnf\n5J1j8Tzz4V6ydq8kIjSEWwa1Z3aaj6GdW1200fB5x4+7rt4bN7qrNkOGfNHNW+FGRKTxU8CRhqv4\ntJsXas0COLkXWnaE635Gfq9ZPLvpLM8/uY/8gjw6xTXjhzf2YUZyJ+JbRFa723PnICPDPU9JgdGj\noWXLS79HREQaFwUcaXiO5cKax92M3qVnoXMq9rr/ZFVkGovXHOStNzbgt5ZxvRL4dZqPcb0SCQ25\ndJ/tw4fd+DXXXAPNmsG0aeDzuUkwRUTEexRwpGGwFna95ya9zH0HQsJhwHTODruffx5KYNHbeeTm\nr6VVdDj3j+nKXSOT8MU3r3a3Bw64mb137nTj2Awb5kJN//51UCcREak3CjhSv0oLYePzbtLLYzuh\neSJc+2N2dr6dhZvO8dKTBykqzWdw51b89o7B3DKoPVHh1c+JcPIkvPoq7N4N0dEwfrybJyrq4sPe\niIiIxyjgSP04mefa1qxf5NratB9C+ZQM3rBpPJN9mOw3dxAZFsKUwR2Yk+ZjUKdW1e7SWigqcoPy\nNWsGZ87ADTdAcjJERNRBnUREpMFQwJG6Yy3krXSTXu54HTDQbwpH+9/H0/sSeeH1/Rw7u40u8dH8\nx819uX14J1pFV59MrIUdO9ytqPJymDvXXan59rc1nYKISFOlgCO1r6wYNv/D3YY6shmaxWFH/Rtr\nEqbxxKYy3lt0BChgfJ+2zEnzMbZHG0KqaTQM4PfDtm2wYoUbnC8uzg3OZ60LNgo3IiJNlwKO1J4z\nh91Iw2ufgqLjkNiPool/4MWSkTyTc5Q9xw4Q3zyCB8d1566RSXSKi76s3W/fDkuXQps2cNttMGCA\nxrARERFHAUeCb3+2G21420tukL7eN7Gn+xzm53Xk5dcPUVy2l+G+OP59Qk9uHNCOyLDqGw2Du/20\ncaMbkG/IEDdP1MyZ0Lu3go2IiHyZAk4T9+ijj5KSkkJ6evrnyzIzM8nOzuahhx6q+Y7KS2Hby260\n4YNrIbIl5ckP8F7MFOZvtqxfdopm4YeZNrQTs1OT6N8htsa7LiuDdevclApnzrhAc37k4b59L6e2\nIiLSVCjgNHEpKSnMmDGDJUuWkJ6eTmZm5ueva+TsUXcLKvvvcPYzaN2dk+N+xZOFqSzOOcHJopN0\nS2jOzyb3Y/rwTrSMCr+s8m3dCq+/DoWFbmC+qVOhW7crqKiIiDQpCjhNXHp6OkuWLGHGjBnMnTuX\njIyMz8POJR3e5G5DbV4KFSXY7texadjD/DnPx/K3jxFi8rm+b1vuTvOR1j3+kvNCXai42DUUbtbM\nPdq1cyMQ+3xXWVkREWkyFHCE9PR05s6dy8MPP8y8efOqDjcV5bDjNdcbKm8lhEdTPPAuXo64mce2\nhLJ/6zkSYgr4zvie3DmiM+1jm11WOYqKICsLVq92Iw5PnOiu1uiKjYiIXC4FHCEzM5OMjAzmzZtH\nRkYG6enpXw45507Cumdgzd/g9H5sbGcOjfgJGadHsySngNLyUkZ2bc0Pb+zDDf3aERF2eS1+Cwpg\n1SrIznYNifv1g8GDg1xJERFpUhRwmrjKbW7OB5vPX/dv7ya93PgClBVRkTSarB7f59E93dj44Vma\nR5xlZnJnZqf66N0u5orL8O67sHkzDBwIY8ZAQkIQKygiIk2SsdbWy4GTk5NtTk5OvRxbvvCVXlR+\nP5uW/Y6Ybc/S1b8HQiMp6HUrL5hJPLY9mtPnyujVtgVz0rowbWhHWkRefkY+cQI++ghGjoS2bd28\nUdZC69ZBrpyIiHiOMWattTa5uu10BaeJ+7wreEkBbHgOVj/OoBO7sDHtyfV9lz+cHMVr68sICzHc\nOKANc1J9jOja+rIaDZ939KgbdXjzZtfFOynJBZy4uCBXSkREmjwFnKbuxG5YvQDWL4bSAsraDyez\n76/45e5e7Mspo13LUL53fVdmpXQmseWVT8X98suwYQOEh0NamnvEXPldLRERkUtSwGmKrIU9H0DW\nfNj5JjYklBNdbubpiolk5MZRtscyukdLfjLZx4S+bQkLvbJhgj/7zF2hMQZatXLzRKWmQvTlzcgg\nIiJy2RRwmpLSItj0ouvmfXQ7NroNW7s/wCPHRrNyWzgxUWHMTu3E10b66JHY4ooPk5fnZvbetQtm\nz4YePWDcuCDWQ0REpBoKOE3B6QOui/e6p+HcSUra9OcN3095eG9fjm8NoW/7ljxym4+pQzoQHXFl\nHwlrYfduF2zy8qBFC7jhBtfORkREpK4p4HiVtbAvy402vP1VLJYjHa5jQclEnjzQnojQUCYNbMec\nNB/DkuKuqNFwZRUV8NJL7nbUpEkwdKhrbyMiIlIfFHC8prwEtixzk14e3og/MpZ1He7iV/mjWb+r\nJR1bNeP/TExiZkpn2rSIvOLD+P2wfbub3XvmTAgLgzlzXFfvMH2qRESknulXkVcUHIGcJ92jMJ+i\n2B68nPhd/t+BQRScjmRcrwSemOYjvU8ioSFXfrXG73fdvFesgGPHoE0bOH3aBZvExCDWR0RE5Coo\n4DR2B9e521BbloG/jAMJY3nMfpsXjvQgtlkEs0a7RsNd2jS/6kOdOQNPPeUG5mvbFu64A/r2hZAr\n62QlIiJSa6oNOMaYJ4FbgHxr7YCLrI8FFgNJgf391lr7VLALKpVUlMH2V12w2b8af3hzVreeyi/y\nx7B9fyIDO8by6O0+pgzuQFR46FUdqqwMjhyBTp3cuDVJSW4SzN69XXsbERGRhqgmV3AWAo8Bz1Sx\n/tvANmvtZGNMArDDGPOstbY0SGWU84pOwNqnIPvvcOYghc2TWNLyQX6fn0zJuRZMHtSBX6f5GNy5\n1VUfqrQUcnLg449dyPne9yAyEqZNC0I9REREalm1Acda+6ExpsulNgFijOuG0wI4AZQHpXTiHNnq\nrtZsWgLlxeTFpvCn0Dm8dHwAnVq34DuTkrhjeGfimkdc9aGKi2HNGje797lz0K0bXHMNRFz9rkVE\nROpMMNrgPAa8AhwCYoCZ1lr/xTY0xnwT+CZAkgZIuTR/Bex8E7IyYO8KKkKjWNn8Oh45Po5P8jsx\nvncif0/zMa5nAiFX0Wj4QidOwPLl0KuXG3m4c+eg7VpERKTOBCPgTAQ2AOOB7sA7xpgV1tozF25o\nrV0ALAA3m3gQju09xafdvFBrFsDJvZyNbMdzEXfz1zOjCSGeGWM6s2BkEp1bB2e+g7Nn3W2o8nI3\nfk2HDvCd70B8fFB2LyIiUi+CEXDuBX5trbVArjFmD9AHWBOEfTcdx3LdbagNz0FZIXuiB/Gnin/n\n1dPDGZwUz89u8nHTgPZX3Wj4vNOnYeVKWLfODdI3eLAbG9AYhRsREWn8ghFw9gHXASuMMW2B3sDu\nIOzX+/x+2L3cTXqZ+w4VIeF8GHENvzubTm5FD24d0pGXU30M6Bgb1MNu3uxGHbbWBZsxYxRqRETE\nW2rSTfx54FqgjTHmAPAzIBzAWjsfeBhYaIzZDBjgh9baY7VWYi8oOQsbn3e3oY7t5Gx4axabGTxR\nlE5MdAdm3+zj9mGdiI0O3lwHx465PJWY6Lp6DxsGo0e7Wb5FRES8xrg7S3UvOTnZ5uTk1Mux683J\nPFizALvuGUzJGfZE9OLPhRN43Z/KuL4duTutC6O6xwe10fCRI24CzG3bXMPhO+8M2q5FRETqnDFm\nrbU2ubrtNJJxbbMW9n4Eq+djd7yOxfB+SCqPlVzPvtAB3DkuicwRSXRo1Syohz10CD74AHbscOPX\njBkDqalBPYSIiEiDpYBTW8qKYfM/YPXjcGQzhaGxLKqYzMLSCSR16cnX03zc2L8dEWHBnefgfEPh\nXbtg3z5IT4cRI6BZcPOTiIhIg6aAE2xnDkH2E9i1CzFFx9kb6iOj7AHe9l/DpKFdWZjmo0+7lkE9\npLWwZ4+7FZWcDAMGwMiRLthEXvmE4SIiIo2WAk6w7M+G1RnYbS+Dv4L3SWZB6YMcjR/B3RO68B9D\nOxITFbxGw+CCzaefumBz4AC0rJSbNPKwiIg0ZQo4V6O8FLa9hM2ajzm0liITzXNl1/OsfyL9+g/m\nf6f6SO3WGlNLs1IuXQpbt7qeUJMnuy7fYTqjIiIiCjhX5OxRWPsU/jVPEFJ4hP2mA38r+zoroidw\n6+jevDAiibYto4J+WL/fBZrevd0VmkGDoGdPGDgQQoMz/p+IiIgnKOBcjsMbsavnYzctJcRfykf+\nwTxZfg+lvnTmjOrK/+3XlvDQ4DYaBjfS8MaN8NFHbq6oKVPcODa9ewf9UCIiIp6ggFOdinLY8RoV\nqzII3b+KEiJZUj6OpaGTGJY8kv9ITaJHYkytHNpayM52UyqcPu3miZo1S8FGRESkOgo4VSk6AesX\nUbbqccLPHuQwCSws+xrr42/httH9eX5IR5pH1s6Pz++HkBDX3XvrVtd4ePJk6N7dLRMREZFLU8C5\nUP4n+LMy8G98kbCKc+RU9GORnUl4v5uZPaobP/XF1Vqj4eJid8UmOxseeABiYuCuu1x7GwUbERGR\nmvNkwHn00UdJSUkhPT3982WZmZlkZ2fz0EMPffUNfj/kvkPpyr8QkfcBZYTzUvloXms2mZHjruXn\nyZ1JiKm9AWWKimD1avcoLnYNh8vK3DqNYyMiInL5PBlwUlJSmDFjBkuWLCE9PZ3MzMzPX39J8Rns\nhmcpWTmfqIK9nLBxLCqfwR7f7dw6ejBP9kkkrBYaDX+pCMXw5z+7f/v2hbFjXVsbERERuXKeDDjp\n6eksWbKEGTNmMHfuXDIyMj4POwCc2E3px/Nh/WIiKgrZ5u/BCyH/Ruzw27kzrRvdElrUavnOnIHc\nXNcTKioKxo+HLl3cTN8iIiJy9TwZcMCFnLlz5/Lwww8zb9480q+9FnZlcvbDx4jOew9jQ3jNP5IP\n4qaTOvYGfj64I80iancwmZMnXVfvDRvc6549XTubESNq9bAiIiJNjmcDTmZmJhkZGfxi3o84+eFf\nOfrIMhJK91NsW7LQTuNY768xdexwft+5Va01Gj6voADeew82bXKNhYcOdbN7x9RO73IREZEmz5MB\nJzMzk+/edwcf/fJmfEcXE3VtGVuLQ/hLyIN0HHcPd43sQevmtT9ZU3m5mzohNNTNGTViBIwa9eU5\no0RERCT4PBlwSlb8ibX3VMChl3nbn8L6DrOIbp1E1N71PDC+X60f/+BBWLHCXbm5/36Ijobvflfz\nRImIiNQVT/7KLepxE4t2xlA45F4mX5PKTfHRgTXX1+px9+1zM3vn5rrGw6mprgd6aKjCjYiISF3y\n5K/dCbfdR0To/USF190MlNu3w4svuqs1EyZASorGsBEREakvngw4LaPCa/0Y1rorNRUV0KeP6xE1\naRIMGeJGHhYREZH648mAU5ushU8+cbeiDh8Gn88FnLAwdfcWERFpKBRwLsOuXfDWW5CfD61bw9Sp\nMGhQfZdKRERELqSAU42KCveIiHDdvgGmT4f+/d2M3yIiItLwKOBUobwc1q+HlSvdVZrx46FXL/fQ\nzN4iIiINmwLOBcrKYO1aF2wKCqBTJ9fOBhRsREREGgsFnAv861+wcSN07Qq33eYmwVSwERERaVya\nfMA5dw5Wr4aBAyE+3s0RNXw4JCXVd8lERETkSjXZgFNYCKtWwZo1UFrqBuiLj4eEhPoumYiIiFyt\nJhlw3n3XXbUpL3e9oa65BhIT67tUIiIiEixNJuAUFEBMjHt+PtiMGQNt2tRvuURERCT4PB9wjh93\nM3tv2gT33ON6RE2cqIbDIiIiXubZgJOf74LNli1uNu+UFDf6MCjciIiIeJ0nA05FBTz9tBvTZtQo\nSEuDFi3qu1QiIiJSVzwZcEJDYcYM1yMqOrq+SyMiIiJ1zZMBB74YfVhERESaHk0XKSIiIp6jgCMi\nIiKeo4AjIiIinqOAIyIiIp6jgCMiIiKeo4AjIiIinqOAIyIiIp6jgCMiIiKeo4AjIiIinqOAIyIi\nIp6jgCMiIiKeo4AjIiIinqOAIyIiIp6jgCMiIiKeo4AjIiIinqOAIyIiIp6jgCMiIiKeo4AjIiIi\nnqOAIyIiIp6jgCMiIiKeo4AjIiIinqOAIyIiIp6jgCMiIiKeo4AjIiIinqOAIyIiIp6jgCMiIiKe\no4AjIiIinqOAIyIiIp6jgCMiIiKeo4AjIiIinqOAIyIiIp6jgCMiIiKeo4AjIiIinqOAIyIiIp5T\nbcAxxjxpjMk3xmy5xDbXGmM2GGO2GmM+CG4RRURERC5PTa7gLARurGqlMaYV8FdgirW2P3BHcIom\nIiIicmWqDTjW2g+BE5fY5C5gmbV2X2D7/CCVTUREROSKBKMNTi8gzhjzvjFmrTHm7qo2NMZ80xiT\nY4zJOXr0aBAOLSIiIvJVwQg4YcBw4GZgIjDPGNPrYhtaaxdYa5OttckJCQlBOLSIiIjIV4UFYR8H\ngOPW2kKg0BjzITAY2BmEfYuIiIhctmBcwXkZGGOMCTPGRAMjge1B2K+IiIjIFan2Co4x5nngWqCN\nMeYA8DMgHMBaO99au90Y8yawCfADT1hrq+xSLiIiIlLbqg041to7a7DNb4DfBKVEIiIiIldJIxmL\niIiI5yjgiIiIiOco4IiIiIjnKOCIiIiI5yjgiIiIiOco4IiIiIjnKOCIiIiI5yjgiIiIiOco4IiI\niIjnKOCIiIiI5yjgiIiIiOco4IiIiIjnKOCIiIiI5yjgiIiIiOco4IiIiIjnKOCIiIiI5yjgiIiI\niOco4IiIiIjnKOCIiIiI5yjgiIiIiOco4IiIiIjnKOCIiIiI5yjgiIiIiOco4IiIiIjnKOCIiIiI\n5yjgiIiIiOco4IiIiIjnKOCIiIiI5yjgiIiIiOco4IiIiIjnKOCIiIiI5yjgiIiIiOco4IiIiIjn\nGGtt/RzYmKNAXi0eog1wrBb339A0tfpC06uz6ut9Ta3Oqq/31UadfdbahOo2qreAU9uMMTnW2uT6\nLkddaWr1haZXZ9XX+5panVVf76vPOusWlYiIiHiOAo6IiIh4jpcDzoL6LkAda2r1haZXZ9XX+5pa\nnVVf76u3Onu2DY6IiIg0XV6+giMiIiJNlAKOiIiIeE6jDjjGmDuMMVuNMX5jTJXd0IwxNxpjdhhj\nco0xP6q0vKsxZnVg+YvGmIi6KfmVMca0Nsa8Y4z5NPBv3EW2STfGbKj0KDbG3BpYt9AYs6fSuiF1\nX4vLU5M6B7arqFSvVyot9+I5HmKMWRX47G8yxsystK5RnOOqvpOV1kcGzldu4Px1qbTux4HlO4wx\nE+uy3FeqBvX9njFmW+B8vmeM8VVad9HPdkNWg/p+3RhztFK97q+07p7A5/9TY8w9dVvyK1eDOv+h\nUn13GmNOVVrXGM/xk8aYfGPMlirWG2PMnwM/j03GmGGV1tXNObbWNtoH0BfoDbwPJFexTSiwC+gG\nRAAbgX6BdUuAWYHn84G59V2naur7KPCjwPMfAf9VzfatgRNAdOD1QuD2+q5HbdQZOFvFcs+dY6AX\n0DPwvANwGGjVWM7xpb6Tlbb5FjA/8HwW8GLgeb/A9pFA18B+Quu7TkGob3ql7+nc8/UNvL7oZ7uh\nPmpY368Dj13kva2B3YF/4wLP4+q7TsGo8wXbfwd4srGe40CZrwGGAVuqWD8JeAMwQCqwuq7PcaO+\ngmOt3W6t3VHNZiOAXGvtbmttKfACMNUYY4DxwNLAdk8Dt9ZeaYNiKq6cULPy3g68Ya0tqtVS1a7L\nrfPnvHqOrbU7rbWfBp4fAvKBakf1bEAu+p28YJvKP4elwHWB8zkVeMFaW2Kt3QPkBvbXkFVbX2tt\nZqXvaRbQqY7LGEw1Ob9VmQi8Y609Ya09CbwD3FhL5Qymy63zncDzdVKyWmKt/RD3B3RVpgLPWCcL\naGWMaU8dnuNGHXBqqCOwv9LrA4Fl8cApa235BcsbsrbW2sOB558BbavZfhZf/RL9KnC58A/GmMig\nlzD4alrnKGNMjjEm6/wtOZrAOTbGjMD9xbir0uKGfo6r+k5edJvA+TuNO581eW9Dc7ll/gbuL9/z\nLvbZbshqWt/pgc/pUmNM58t8b0NT43IHbj92BZZXWtzYznFNVPUzqbNzHFYbOw0mY8y7QLuLrPqp\ntfblui5PbbtUfSu/sNZaY0yVffwDSXkg8FalxT/G/dKMwI1N8EPgF1db5qsVpDr7rLUHjTHdgOXG\nmM24X4oNTpDP8SLgHmutP7C4QZ5jqRljzGwgGRhXafFXPtvW2l0X30Oj8SrwvLW2xBjzv3BX68bX\nc5nqyixgqbW2otIyL57jetfgA461dsJV7uIg0LnS606BZcdxl8zCAn8hnl9ery5VX2PMEWNMe2vt\n4cAvt/xL7GoG8N/W2rJK+z5/ZaDEGPMU8IOgFPoqBaPO1tqDgX93G2PeB4YC/8Sj59gY0xJ4DRf0\nsyrtu0Ge4wtU9Z282DYHjDFhQCzuO1uT9zY0NSqzMWYCLuSOs9aWnF9exWe7If/yq7a+1trjlV4+\ngWt7dv69117w3veDXsLgu5zP5Szg25UXNMJzXBNV/Uzq7Bw3hVtU2UBP43rTROA+XK9Y19opE9dO\nBeAeoKFfEXoFV06ovrxfuccb+IV5vm3KrcBFW783MNXW2RgTd/5WjDGmDTAa2ObVcxz4HP837v72\n0gvWNYZzfNHv5AXbVP453A4sD5zPV4BZxvWy6gr0BNbUUbmvVLX1NcYMBR4Hplhr8ystv+hnu85K\nfmVqUt/2lV5OAbYHnr8F3BCodxxwA1++Ct1Q1eQzjTGmD65h7apKyxrjOa6JV4C7A72pUoHTgT/A\n6u4c10bL5bp6ANNw9+9KgCPAW4HlHYDXK203CdiJS8Q/rbS8G+4/x1zgH0BkfdepmvrGA+8BnwLv\nAq0Dy5OBJypt1wWXkkMueP9yYDPul95ioEV91ykYdQZGBeq1MfDvN7x8joHZQBmwodJjSGM6xxf7\nTuJupU0JPI8KnK/cwPnrVum9Pw28bwdwU33XJUj1fTfwf9j58/lKYHmVn+2G/KhBfR8BtgbqlQn0\nqfTe+wLnPRe4t77rEqw6B17/J/DrC97XWM/x87genGW438PfAB4EHgysN8BfAj+PzVTq6VxX51hT\nNYiIiIjnNIVbVCIiItLEKOCIiIiI5yjgiIiIiOco4IiIiIjnKOCIiIiI5yjgiIiIiOco4IiIiIjn\n/A+86AIctXJf3AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f700be0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "grid_size = 0.01\n",
    "x_grid = np.arange(-1, 1, grid_size)[:, None]\n",
    "\n",
    "Phi_X = np.hstack((\n",
    "    np.ones((x_grid.shape[0], 1)),  # pad with 1s for the bias term\n",
    "    x_grid\n",
    "))\n",
    "\n",
    "lsqsq_w = np.linalg.lstsq(Phi_X_in, y)[0]\n",
    "\n",
    "plt.clf()\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(X_in, y, 'xk', label='observations')\n",
    "plt.plot(x_grid, Phi_X @ true_w, '--b', alpha=0.5, label='true weights')\n",
    "plt.plot(x_grid, Phi_X @ w_n, label='Bayesian')\n",
    "plt.plot(x_grid, Phi_X @ lsqsq_w, label='lstsq')\n",
    "plt.legend()\n",
    "plt.ylim([1.5, 2.5])\n",
    "maybe_save_plot('2018-01-09-mean')  # Graph showing three lines, all near each other, going through the points.\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "These are pretty similar! They are different at least in part due to the prior, which are centered at 0, meaning that it expects most lines to go through the origin and have a slope of 0.\n",
    "\n",
    "This is more obvious if I make the true bias is very far away from 0. Then the Bayesian fit might not even go through the points! This might remind you of the effects of regularization, which makes extreme values less likely, at the cost of sometimes having poorer fits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "true_w = np.array([[100, 0.3]]).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Graph showing three lines, two near each other going through the points, and one closer to 0.](images/2018-01-09-mean-far-away.png)\n",
    "\n",
    "\n",
    "### [Next](2018-01-10-bayesian-linreg-plots.ipynb)\n",
    "\n",
    "## See Also\n",
    "\n",
    " - Still thanks to [MLPR](http://www.inf.ed.ac.uk/teaching/courses/mlpr/2017/notes/)!\n",
    " - Wikipedia on [Conjugate Priors](https://en.wikipedia.org/wiki/Conjugate_prior)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}