{ "cells": [ { "cell_type": "code", "execution_count": 57, "id": "1ce7bbe3", "metadata": {}, "outputs": [], "source": [ "%matplotlib inline \n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "def gellipse(mu, cov, n=100, *args, **kwargs):\n", " \"\"\"\n", " Contour plot of 2D Multivariate Gaussian distribution.\n", "\n", " gellipse(mu, cov, n) plots ellipse given by mean vector MU and\n", " covariance matrix COV. Ellipse is plotted using N (default is 100)\n", " points. Additional parameters can specify various line types and\n", " properties. See description of matplotlib.pyplot.plot for more details.\n", " \"\"\"\n", " mu = mu.copy()\n", " if mu.shape == (2,):\n", " mu.shape = (2, 1)\n", " if mu.shape != (2, 1) or cov.shape != (2, 2):\n", " raise RuntimeError('mu must be a two element vector and cov must be 2 x 2 matrix')\n", "\n", " d, v = np.linalg.eigh(4 * cov)\n", " d = np.diag(d)\n", " t = np.linspace(0, 2 * np.pi, n)\n", " x = v @ np.sign(d) @ np.sqrt(np.abs(d)) @ np.array([np.cos(t), np.sin(t)]) + mu\n", " return plt.plot(x[0], x[1], *args, **kwargs)" ] }, { "cell_type": "markdown", "id": "6229cfc7", "metadata": {}, "source": [ "For 2D gaussian distribution\n", "\n", "$p\\left(\\begin{bmatrix} x \\\\ y \\end{bmatrix}\\right)=\\mathcal{N}\\left(\\begin{bmatrix} x \\\\ y \\end{bmatrix} \\middle| \\begin{bmatrix} \\mu_x \\\\ \\mu_y \\end{bmatrix}, \\begin{bmatrix} \\Sigma_{xx} & \\Sigma_{xy} \\\\ \\Sigma_{yx} & \\Sigma_{xx} \\end{bmatrix}\\right)$\n", "\n", "the conditional probability\n", "\n", "$p(x | y) = \\mathcal{N}\\left(x \\middle| \\mu_{x|y}, \\Lambda_{xx}^{-1}\\right)$\n", "\n", "where\n", "\n", "$\\mu_{x|y} = \\mu_{x} - \\Lambda_{xx}^{-1} \\Lambda_{xy}(y-\\mu_y)$\n", "\n", "and\n", "\n", "${\\begin{bmatrix} \\Sigma_{xx} & \\Sigma_{xy} \\\\ \\Sigma_{yx} & \\Sigma_{xx} \\end{bmatrix}}^{-1}\n", "={\\begin{bmatrix} \\Lambda_{xx} & \\Lambda_{xy} \\\\ \\Lambda_{yx} & \\Lambda_{xx} \\end{bmatrix}}$" ] }, { "cell_type": "code", "execution_count": 74, "id": "6dcef7c4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mean: [1.89401505 1.53989025]\n", "Cov:\n", " [[12.99098534 9.96611533]\n", " [ 9.96611533 8.692699 ]]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABhX0lEQVR4nO3dd3hUVfoH8O/MJJk0kgAhPSQBQu+9KE0BQVGRqi6CIgKiK7LuWtCf4Kos4rqKIIiyFLHgLoIIKGWlCtJBSggJBJKQhISSDGkzmZn7++MwuTPJpJKp+X6e5z65c+bcmTMp3JdT3qOQJEkCERERkYtQOroBRERERDXB4IWIiIhcCoMXIiIicikMXoiIiMilMHghIiIil8LghYiIiFwKgxciIiJyKQxeiIiIyKV4OLoBdc1oNCIjIwMNGjSAQqFwdHOIiIioGiRJwu3btxEREQGlsvK+FbcLXjIyMhAdHe3oZhAREVEtpKWlISoqqtI6bhe8NGjQAID48AEBAQ5uDREREVWHRqNBdHR06X28Mm4XvJiGigICAhi8EBERuZjqTPnghF0iIiJyKQxeiIiIyKUweCEiIiKXwuCFiIiIXAqDFyIiInIpDF6IiIjIpTB4ISIiIpfC4IWIiIhcCoMXIiIicim1Dl727t2LkSNHIiIiAgqFAhs3brR4fvLkyVAoFBZH7969q3zd9evXo23btlCr1Wjbti02bNhQ2yYSERGRG6p18FJQUIBOnTph8eLFFdZ54IEHkJmZWXps3bq10tc8ePAgxo8fj4kTJ+LUqVOYOHEixo0bh0OHDtW2mUREVF3p6cCuXeIrUUWc4PdEIUmSdNcvolBgw4YNePTRR0vLJk+ejNzc3HI9MpUZP348NBoNfv7559KyBx54AA0bNsS3335brdfQaDQIDAxEXl4e9zYiIqquFSuA554DjEZAqQQ+/xx49lnxnE4HlJQAHh6AWi1fU1Agvvr4iGsAUU+nA1QqwNu7dnULCwFJEmUqlSjT6wGtVlzr41O7ukVF4vOp1eKzAIDBABQX16yuQgH4+sp1i4vFc15egKdnzesajeL9AMDPT66r1YrP4ukp6te0riSJ7w8g2mDaM8j086xJXdPPvuzvyfLlwJQpqAs1uX/bdM7L7t27ERISgpYtW2Lq1KnIzs6utP7BgwcxdOhQi7Jhw4bhwIEDFV6j1Wqh0WgsDiIiqoH0dPmGBIiv06fL/7NeuBDw9wdeeMHyupAQUZ6aKpctWSLKyt7QYmNFeUKCXLZqlSibMMGybtu2ovz4cbls3TpR9vDDlnV79BDl+/bJZZs3i7L777es27+/KN+2TS779VdR1qePZd3hw0W5+dSF338XZZ06WdYdPVqUf/21XHb6tCiLj7esO3GiKF++XC67eFGURUZa1p02TZR/8olclpkpyoKCLOvOni3K339fLsvLE2X+/iKwMZkzR5TNmSOX6fVy3bw8ufz990XZ7NnWf0+mTXNID4zNgpfhw4fj66+/xq+//op//vOfOHLkCAYPHgytVlvhNVlZWQgNDbUoCw0NRVZWVoXXzJ8/H4GBgaVHdHR0nX0GIqJ6ISlJviGZGAxAcrJj2kPOyYl+T2w2bFRWZmYmYmJi8N133+Gxxx6zWsfLywurV6/G448/Xlr29ddfY8qUKSguLrZ6jVartQiINBoNoqOjOWxERFRd6elATIzljUmlAi5fBqKiOGzEYSPxPcjJqfz35C7VZNjI467frZrCw8MRExODpKSkCuuEhYWV62XJzs4u1xtjTq1WQ23+x0RERDUTFSWGMaZNEzdWlUrMeTHdkLy85JucOfObp4mnp3xTrm1d85u9iYeHHETUtq55cGKiUllvW03qmgdetamrVFqvq1ZbBos1ratQWK9r7edZnbpRUcCyZcCMGdZ/T+zIbnlebty4gbS0NISHh1dYp0+fPtixY4dF2fbt29G3b19bN4+IqH6bMkX8D3rXLvG1jiZhkpuZOtUpfk9q3fOSn5+PZLNxrpSUFJw8eRKNGjVCo0aNMHfuXIwePRrh4eG4fPky3njjDQQHB2PUqFGl1zz11FOIjIzE/PnzAQAvvfQS+vfvjwULFuCRRx7Bjz/+iJ07d2L//v138RGJiKhaoqIc8r9ocjFO8HtS6+Dl6NGjGDRoUOnj2bNnAwAmTZqEpUuX4vTp01izZg1yc3MRHh6OQYMGYd26dWjQoEHpNampqVAq5c6fvn374rvvvsObb76Jt956C82bN8e6devQq1ev2jaTiIiI6opOJ69+eukl68OJdlAnE3adCfO8EBER2UhBgVg6DQD5+dbnydSSU07YJSIiIhfn4QFMmiSfO6oZDntnIiIici1qtUgu6GDcVZqIiIhcCoMXIiIicikMXoiIiKh6CgrEvkpBQXLWZAfgnBciIiKqPvONGx2EwQsRERFVj48PcOGCfO4gDF6IiIioepRKID7e0a3gnBciIiJyLex5ISIiouopKRE7kAPAc89Z3xXcDhi8EBERUfXodMALL4jzyZMZvBAREZGTU6mAMWPkcwdh8EJERETV4+0N/Oc/jm4FJ+wSERGRa2HwQkRERC6FwQsRERFVT2EhEBkpjsJChzWDc16IiIioeiQJyMiQzx2EwQsRERFVj7c3cOKEfO4gDF6IiIioelQqoHNnR7eCc16IiIjItbDnhYiIiKqnpAT4+mtx/uSTzLBLRERETk6nA55+WpyPHcvghYiIiJycSgWMGCGfOwiDFyIiIqoeb29gyxZHt4ITdomIiMi1MHghIiIil8LghYiIiKqnsBCIjxcHtwcgIiIipydJQHKyfO4gDF6IiIioery9gf375XMHYfBCRERE1aNSAf36OboVnPNCREREroU9L0RERFQ9ej2wYYM4HzUK8HBMGMHghYiIiKpHqwXGjRPn+fkMXoiIiMjJKZXAgAHyuaOaUdsL9+7di5EjRyIiIgIKhQIbN24sfa6kpASvvvoqOnToAD8/P0REROCpp55CRkZGpa+5atUqKBSKckdxcXFtm0lERER1xccH2L1bHD4+DmtGrYOXgoICdOrUCYsXLy73XGFhIY4fP4633noLx48fxw8//IALFy7g4YcfrvJ1AwICkJmZaXF4O3A5FhERETmXWg8bDR8+HMOHD7f6XGBgIHbs2GFR9umnn6Jnz55ITU1F06ZNK3xdhUKBsLCw2jaLiIiI3JzdBqzy8vKgUCgQFBRUab38/HzExMQgKioKDz30EE6cOFFpfa1WC41GY3EQERGRDRQVAZ07i6OoyGHNsEvwUlxcjNdeew1PPPEEAgICKqzXunVrrFq1Cps2bcK3334Lb29v9OvXD0lJSRVeM3/+fAQGBpYe0dHRtvgIREREZDQCp06Jw2h0WDMUknT3mxMoFAps2LABjz76aLnnSkpKMHbsWKSmpmL37t2VBi9lGY1GdO3aFf3798eiRYus1tFqtdBqtaWPNRoNoqOjkZeXV6P3IiIioioYDMCvv4rzwYNFxt06otFoEBgYWK37t02XSpeUlGDcuHFISUnBr7/+WuNgQqlUokePHpX2vKjVaqjV6rttKhEREVVFpQKGDHF0K2w3bGQKXJKSkrBz5040bty4xq8hSRJOnjyJ8PBwG7SQiIiIXFGte17y8/ORbNoWG0BKSgpOnjyJRo0aISIiAmPGjMHx48exefNmGAwGZGVlAQAaNWoELy8vAMBTTz2FyMhIzJ8/HwAwb9489O7dG/Hx8dBoNFi0aBFOnjyJJUuW3M1nJCIiorqg1wPbtonzYcNcL8Pu0aNHMWjQoNLHs2fPBgBMmjQJc+fOxaZNmwAAnTt3trhu165dGDhwIAAgNTUVSrMMfbm5uXjuueeQlZWFwMBAdOnSBXv37kXPnj1r20wiIiKqK1ot8NBD4tyB2wPUyYRdZ1KTCT9ERERUA0VFQP/+4nzv3jrNsus0E3aJiIjIjfj4AEeOOLoV9ktSR0RERFQXGLwQERGRS2HwQkRERNVTVAT06ycOB24PwDkvREREVD1GI3DggHzuIAxeiIiIqHrUamDDBvncQRi8EBERUfV4eABW9jG0N855ISIiIpfCnhciIiKqHoMB2LdPnN97b53uKl0TDF6IiIioeoqLAdPWQPn5gJ+fQ5rB4IWIiIiqR6EA2raVzx2EwQsRERFVj68vcPaso1vBCbtERETkWhi8EBERkUth8EJERETVU1QEDBkiDm4PQERERE7PaAR27pTPHYTBCxEREVWPWg2sXSufOwiDFyIiIqoeDw/gyScd3QrOeSEiIjeXng7s2iW+0t1zgu8ngxciInJfK1YAMTHA4MHi6/Ll8nNGI1BQIA5zxcWirKSk6rparSjT6eQySapdXUmSy3U6UabVWr5Gbeqaz00pKRFlxcW1q/vFF5bfzxUr4AgMXoiIyD2lpwPPPSffkI1GYMYMuccgIQHw9wdiYy2vmzJFlC9ZIpelpoqykBDLui+8IMoXLpTLrl8XZf7+lnVffVWUzZsnlxUWynULC+XyefNE2auvWr6Gqe7163LZwoWi7IUXLOuGhIjy1FS5bMkSUTZlimXd2FhRnpAgl61aJcomTBCP09OB6dMtv5/TpjmkB4bBCxERuaekpPIrYoxGIDnZMe1xdda+nwaDQ76fCkky73tyfRqNBoGBgcjLy0NAQICjm0NERI6Sni6GNsxvuCoVcPkyEBUlyk25Ssw3GCwuFjdlLy/A01OUVVRXqwX0elHPy0uUSZLci1KTur6+8n5BOp0YtvHwsFzVYxqKqkldHx9AeaevoqRE1FepAG/vmtWt6vt5l2py/2bPCxERuaeoKDHHRaUSj1Uq4PPP5RutUimCi7I7I3t7izJT4FJZXbValJmCEUAEFbWpa77RoZeXKCu7HLk2dZVmt3pPT1FmHrhUt25V3087Ys8LERG5r8JCoGVL0eOxbx8QH+/oFrm+9HQxVNSiRZ0GLjW5fzPPCxERuS9JAq5eFecREY5ti7uIinJIb4s5Bi9EROS+vL2Bw4flc3ILDF6IiMh9qVRAjx6ObgXVMU7YJSIiIpfCnhciInJfej2wbp04Hz9eLCcml8efIhERuS+tFvjTn8T5o48yeHET/CkSEZH7UiqB+++Xz8ktMHghIiL35eMD7Njh6FZQHat1GLp3716MHDkSERERUCgU2Lhxo8XzkiRh7ty5iIiIgI+PDwYOHIizZ89W+brr169H27ZtoVar0bZtW2zYsKG2TSQiIiI3VOvgpaCgAJ06dcLixYutPv/BBx/go48+wuLFi3HkyBGEhYVhyJAhuH37doWvefDgQYwfPx4TJ07EqVOnMHHiRIwbNw6HDh2qbTOJiIjIzdTJ9gAKhQIbNmzAo48+CkD0ukRERGDWrFl49c523lqtFqGhoViwYAGmTZtm9XXGjx8PjUaDn3/+ubTsgQceQMOGDfHtt99Wqy3cHoCIiEoVFsp5Xo4cERsaklNy+MaMKSkpyMrKwtChQ0vL1Go1BgwYgAMHDlR43cGDBy2uAYBhw4ZVeo1Wq4VGo7E4iIiIAIjtAc6dE4d7beVXr9kkeMnKygIAhIaGWpSHhoaWPlfRdTW9Zv78+QgMDCw9oqOj76LlRETkVry9gV27xMHtAdyGTdeNKcy37IYYTipbdrfXvP7668jLyys90tLSat9gIiJyLyoVMHCgOFQqR7eG6ohNlkqHhYUBED0p4eHhpeXZ2dnlelbKXle2l6Wqa9RqNdRq9V22mIiIiFyFTXpe4uLiEBYWhh1ma+t1Oh327NmDvn37Vnhdnz59LK4BgO3bt1d6DRERUYX0emDjRnHo9Y5uDdWRWve85OfnIzk5ufRxSkoKTp48iUaNGqFp06aYNWsW3n//fcTHxyM+Ph7vv/8+fH198cQTT5Re89RTTyEyMhLz588HALz00kvo378/FixYgEceeQQ//vgjdu7cif3799/FRyQionpLqwVGjRLn+fncHsBN1PqnePToUQwaNKj08ezZswEAkyZNwqpVq/C3v/0NRUVFeP7553Hr1i306tUL27dvR4MGDUqvSU1NhdIsXXPfvn3x3Xff4c0338Rbb72F5s2bY926dejVq1dtm0lERPWZUgmYeu+5PYDbqJM8L86EeV6IiIhcj8PzvBARERHZCoMXIiIicikMXoiIyH0VFYntAXr0EOfkFjjtmoiI3JfRCBw9Kp+TW2DwQkRE7kutBjZvls/JLTB4ISIi9+XhATz4oKNbQXWMc16IiIjIpbDnhYiI3JfBAPz6qzgfPJibM7oJBi9EROS+iouBoUPFeX4+4Ofn2PZQnWDwQkRE7kupBDp1ks/JLTB4ISIi9+XjA5w86ehWUB1jGEpEREQuhcELERERuRQGL0RE5L6KioCBA8XB7QHcBue8EBGR+zIagT175HNyCwxeiIjIfanVwPffy+fkFhi8EBGR+/LwAMaOdXQrqI5xzgsRERG5FPa8EBGR+zIYgN9/F+e9e3N7ADfB4IWIiNxXcTFwzz3inNsDuA0GL0RE5L4UCqBFC/mc3AKDFyIicl++vkBSkqNbQXWME3aJiIjIpTB4ISIiIpfC4IWIiNxXcTHw4IPiKC52dGuojnDOCxERuS+DAdi6VT4nt8DghYiI3JeXF7BypXxOboHBCxERuS9PT2DyZEe3guoY57wQERGRS2HPCxERuS+DATh9Wpx36MDtAdwEgxciInJfxcVAly7inNsDuA0OGxERuaL0dGDXLvGVKqZQAKGhQOPGQEaGo1tDdYTBCxGRq1mxAoiJAQYPFl9XrJCfKy4GCgqAkhK5zGgUZQUFlq+j1Yoyna52dSVJritJcrlOV/O6Wq3l+5nqGo21q1tSIsrWrAFycoAbN4DWrS2/V+SybBq8xMbGQqFQlDtmzpxptf7u3but1j9//rwtm0lE5DrS04HnnpNv1EYjMG2a3AMzcSLg7w8sXy5fc/GiKIuMtHytadNE+SefyGWZmaIsKMiy7uzZovz99+WyvDxR5u8P6PVy+Zw5omzOHLlMr5fr5uXJ5e+/L8pmz7Z8v6AgUZ6ZKZd98okomzbNsm5kpCi/eFEuW75clM2YUfH3ilyWTee8HDlyBAazpEBnzpzBkCFDMHbs2EqvS0xMREBAQOnjJk2a2KyNREQuJSnJsocBEJNSk5OBqCjHtMmV8HvlFhSSZN5/Z1uzZs3C5s2bkZSUBIWVrcl3796NQYMG4datWwgqG/VXk0ajQWBgIPLy8iwCICIit5CeLoaKzAMYlQq4fFnckIuLxQ3ay0vkOAFE3aIicW4+YVWrFT0inp5yArea1JUkoLBQnPv6ivklgBjeKSmpWV0PD0Ctlt/PNGzl4wMolTWvW1ICpKQAbdpU/L0ip1KT+7fd5rzodDqsXbsWzzzzjNXAxVyXLl0QHh6O++67D7t27aq0rlarhUajsTiIiNxWVJQYEjEt+VWpgM8/l2/G3t4i6DAFLoC4ofv5lV9po1aLMvPMszWpq1DIdc3/Xffyqnld82AEkOsqlbWr6+kJtGxZ+feKXJbdlkpv3LgRubm5mFxJpsPw8HAsX74c3bp1g1arxVdffYX77rsPu3fvRv/+/a1eM3/+fMybN89GrSYickJPPSXmjWRnizkccXGObpHzmjIFGDZMDBW1aMHAxU3Ybdho2LBh8PLywk8//VSj60aOHAmFQoFNmzZZfV6r1UJrNvNco9EgOjqaw0ZE5L4KCsRkVIC5S8ht1GTYyC49L1euXMHOnTvxww8/1Pja3r17Y+3atRU+r1aroS7bhUhE5M5UKmDECPmcqJ6xS/CycuVKhISE4MEHH6zxtSdOnEB4eLgNWkVE5KK8vYEtWxzdCiKHsXnwYjQasXLlSkyaNAkeHpZv9/rrr+Pq1atYs2YNAODjjz9GbGws2rVrVzrBd/369Vi/fr2tm0lEREQuwubBy86dO5Gamopnnnmm3HOZmZlITU0tfazT6fDKK6/g6tWr8PHxQbt27bBlyxaMMHWPEhERUb1n1zwv9sA8L0Tk9goLgU6dxPmpUyJvCpGLc7oJu0REVIckSSz9NZ0T1TMMXoiIXI23N7B/v3xOVM8weCEicjUqFdCvn6NbQeQwdtsegIiIiKgusOeFiMjV6PXAhg3ifNQosVEhUT3C33giIlej1QLjxonz/HwGL1Tv8DeeiMjVKJXAgAHyOVE9w+CFiMjV+PgAu3c7uhVEDsOQnYiIiFwKgxciIiJyKQxeiIhcTVER0LmzOIqKHN0aIrvjnBciIldjNIo9jUznRPUMgxciIlfj7Q1s3y6fE9UzDF6IiFyNSgUMGeLoVhA5DOe8EBERkUthzwsRkavR64Ft28T5sGHMsEv1Dn/jiYhcjVYLPPSQOOf2AFQP8TeeiMjVKJVA9+7yOVE9w+CFiMjV+PgAR444uhVEDsOQnYiIiFwKgxciIiJyKQxeiIhcTVER0K+fOLg9ANVDnPNCRFSV9HQgKQmIjweiohzdGrElwIED4jwtDWjZ0rHtIbIz9rwQEVVmxQogJgYYPFh8/eIL+bmSEqCgACgutrymoEAc5vsOVVS3sFCUGwxymV4vysr2qpjqfvUVoFCIsjZtRBuJ6hEGL0REFUlPB557Tg5CjEZgxgxRDgBLlgD+/sCUKZbXxcaK8oQEuWzVKlE2YYJl3bZtRfnx43LZunWi7OGHLev26CHKZ84EJElu07RpcpuI6gEGL0REFUlKKr9rs8EAJCc7pj0mztgmIjtSSJIpfHcPGo0GgYGByMvLQ0BAgKObQ0SuLD1dDBWZBwsqFXD5spj7UlIC6HSizHx354IC8dXHR04iV1HdwkLRi+LtLZ4DxLCRViuu9fGxrJueLoaKKmoTkYuqyf2bPS9ERBWJigKWL5eDCpUK+PxzOUjw9AT8/CyDEUCU+flZZr+tqK6vryg3vQcg0v37+VkGLqa6LVtW3iaieoA9L0REldFqgaefBvLygEWLgObNHd0iIT1dDBW1aMHAhdxCTe7fDF6IiCpTUCAmyQJiE0Q/P8e2h8hN1eT+zTwvRESV8fQE3n1XPicih2PwQkRUGS8vYM4cR7eCiMxwwi4RERG5FJsGL3PnzoVCobA4wsLCKr1mz5496NatG7y9vdGsWTMsW7bMlk0kIqqcJAE5OeJwrymCRC7L5sNG7dq1w86dO0sfq8yXA5aRkpKCESNGYOrUqVi7di1+++03PP/882jSpAlGjx5t66YSEZVXWAiEhIhzTtglcgo2D148PDyq7G0xWbZsGZo2bYqPP/4YANCmTRscPXoUH374IYMXIiIiAmCHOS9JSUmIiIhAXFwcJkyYgEuXLlVY9+DBgxg6dKhF2bBhw3D06FGUlJRYvUar1UKj0VgcRER1xs9PDBdJEntdiJyETYOXXr16Yc2aNdi2bRu++OILZGVloW/fvrhx44bV+llZWQgNDbUoCw0NhV6vx/Xr161eM3/+fAQGBpYe0dHRdf45iIiIyHnYNHgZPnw4Ro8ejQ4dOuD+++/Hli1bAACrV6+u8BqFaZv3O0w59MqWm7z++uvIy8srPdLS0uqo9UREROSM7Jrnxc/PDx06dEBSUpLV58PCwpCVlWVRlp2dDQ8PDzRu3NjqNWq1Gmq1us7bSkQEQGwP8Oqr4nzBAoD/3hA5nF3zvGi1WiQkJCA8PNzq83369MGOHTssyrZv347u3bvDk5kticgR9Hrgk0/Eodc7ujVEBBsHL6+88gr27NmDlJQUHDp0CGPGjIFGo8GkSZMAiCGfp556qrT+9OnTceXKFcyePRsJCQn497//jRUrVuCVV16xZTOJiCrm6Qm88YY4+J8oIqdg02Gj9PR0PP7447h+/TqaNGmC3r174/fff0dMTAwAIDMzE6mpqaX14+LisHXrVrz88stYsmQJIiIisGjRIi6TJiLH8fIC3nvP0a0gIjPcVZqIiIgcjrtKExHVFUkSWXYBwNcXqGDlIxHZDzdmJCKqTGEh4O8vDlMQQ0QOxeCFiIiIXAqHjYiIKuPrKzZkNJ0TkcOx54WIaiY9Hdi1S3ytDxQK4NYt4PBh4OpVR7eGiMDghYhqYsUKICYGGDxYfF2xQpQXFQEFBZZJ3AwGUVZ2nkhN6hYXi3LzjVmNRlFWUFD7ulqtKNPp5DJJsl532TLrn5mIHIbBCxFVT3o68NxzIiAAxNdp00T5/feLCa2bN8v19+0TZT16WL7Oww+L8nXr5LLjx0VZ27aWdSdMEOWrVsllCQmiLDbWsu6UKaJ8yRK5LDVVlIWEWNZ94QVRvnChXHb9ujwx1/wzP/+89c9MRA7D4IWIqicpSb6JmxgMQHKyY9pjD0lJokfGnLt/ZiIXwCR1RFQ96eli2MQ8gFGpgMuXgcaNRblaDXjcWQdgMIihHIXCcqJrUVH16xYXi+e8vOTU/EajeA0A8POrXV2tVgxbeXqK+oBlPhdT3co+c1RUTb57RFSFmty/2fNCRNUTFQUsXy5u3oD4+vnnotzHR9zwPcwWMKpUouxOMGIwStAbjJC8vausW8pU13xPIaVSlJkHIzWtq1aLMlPgAojAqWzdyj4zETkMe16IqGbS08WwSYsWQFQUinQGXMzJR8r1AmTmFSEzrxiZucXI1BRDU1SCfK0e+cV6FJUYSl/CU6WAl0qJ4AZqhAV4IzzQG82b+KNrTEN0ig6Cv9rJsjiU+cxEVPe4PQAR1TlJkpCak48/dp/BGY2EpJQMJOVcQPqtonLTQqpSYpBQYjCg4EYhrtywXGGkVADdYxthTLcojOgQ7hyBTFQUgxYiJ8KeFyKyKl+rx9HLN3Hsyi2cTMvF6at5yC0ssVq3oa8nmjXxR2SQD8KDvBEe4I2wQB808vOCv9oD/moP+KlVUCkVKDFI0BuNKC4xIue2trS35myGBsev3MLV3KLS1/VXe+DFwS3wdL84eHlwlJvIndXk/s3ghYgAAJriEhxJuYlDKTdx6NINnMnQwGC0/OfBS6VAmxtp6KC5ilbTJyI+qhHiQ/zR2F9dZ+24mluEjSeu4r/H0pFyXeRcad7ED59P7I4WIf5VXE1ErorBC4MXoioZjRLOZORhT2IO9ibl4HhqbrlgJbqRD3rGNkbnpkHoFBWIVmENoPZQ2a19/z2ejg9+OY/r+To09PXE6md6omNUkF3en4jsi8ELgxciqzTFJdh1Phu/ns/GvqTruFmgs3g+LtgPveIaoVezRugV1xgRQT4OaqnsRr4WT686gj/S89DYzws7Zg9AIz+vqi8kIpfC4IXBC1GprLxi7Ei4hu1ns/D7pRsoMch/8v5qD/Rt3hgDWjVB//gmiG7knBsP5mv1GP3ZASReu42HO0Vg0eNdHN0kIqpjXG1EVM+l3SzEltOZ+PlMFk6l5Vo81yLEH/e3CcWgVk3QNaYhPFU1mAhbVAQMHy7Of/5Z5HexA3+1BxaO7YhHl/yGTacy8PKQlogL9qv6QiJySwxeiNxEVl4xtpzOxOY/MnAiNbe0XKEAukQHYWi7MAxpG4rmTe5i0qvRCOzZI5/bUceoIPRrEYx9Sdex+VQGXrwv3q7vT0TOg8ELkQvLLdRh8x+Z2HQyA0eu3CzNt6JQAL3jGuPBjuEY2jYUIQHedfOGajXw/ffyeV05ckRs5HjvveU3cjTzUMdw7Eu6jn3J1xm8ENVjDF6IXIxOb8SeCzn44Xg6/peQDZ1B7gHpHtMQD3UMx4gO4XUXsJjz8ADGjq3b15w8GVi9Wn48aZLlLtJmmt3pNcrKK67bNhCRS2HwQuQCJEnC6at5+OH4VWw6lWGxSqhNeABGdYnAQx0jnGJ1UI0cOWIZuADi8cyZVntgQhqI3p7s2wxeiOozBi9ETiyvsAQbT17Ft4dTcT7rdml5kwZqPNo5AqO6RKFtRA1X1aWnA0lJQHx8zVPeGwzA77+L89695Q0La2vfPuvlv/1mNXgp1In9kfy8+E8XUX3GfwGInIwkSTicchPrjqRhy+lMaPViWMjLQ4mhbUMxulsU7m0RDI+arBIyWbECeO45MdlWqQSWLhWPARGYFBeLCTPmuzsXF4vnvLwAnQ645x5Rfu1a+V2Yzeuadnc2GsUqJcCyrlYLdO9uvZ39+lktvpEvepwa+zPPC1F9xuCFyEnkFurw32Pp+OZwKi7lFJSWtwptgAk9ozGqSySCfO/ipp2eLgcugPg6YwbQpQuQnw8UFAAjRwIxMcDly/J1EyYAP/4ILF8OPPmk2Fk5ORkIDQWCg4GcHLnulCnAN98A//oXMGuWKEtNBeLixATf5GS5t+eFF4AvvyzfToUCCA+3+hESr4nep6iGzpmPhojsg8ELkYOduZqHNQcv48eTGaW9LL5eKjzcKQLje0Sjc3QQFArF3b9RUlL55c1GI9CrFyBJoicGEEFMQYHI4aI0690pKBCPk5KAs2eB9u2BwkJR7u1tOYRUUCB6Yby9gXXrRJlWKwKj5ctFkFMRSbIMcswcunQDANAjtlFtvgNE5CaYYZfIAYpLDNh6OhNrDl7BSbMkcm3CAzCxdwwe7hwBf3Ud/98iPV0ED5XlZ1GpRM+HXg9cuCDmxRQXA4sWAa++CowZA/znP/JQUGQkkJcHnDgBdO4s6q5YIXpVRowAPv+8/HuqVKJnp0kT4MoVoE0b68+XCV70BiO6v7cTuYUlWD+jL7rFNKzL7w4RORgz7BI5qWuaYqw5eBnfHk4rXTHkqVJgRIdwPNUnBl2bNqx5L0t1J+BGRYlej2nTxLwUpbJ8IGMwiHkper1c5u1tOVcFENea5rvk5VVc11pvj8Eg96y0bAl8+qlYXQSIwOXzz61+jt2JOcgtLEFjPy90jAqs5BtCRO6OPS9EdnA6PQ8r9l/C5j8yob+zc3NEoDee7B2Dcd2j0aRBLRO+lZ2A+/nnwLPPiud0OqCkRORmMSWUKywUAYNeL5LNDRwImP8ToFQCCQmiR8V82KikRH69J58UZevXi/eVJMthI1NdlQq4fr3inhdTgKLTAX/5C3DrFvDOO0CzZlY/6rSvjmLb2Wt49p44vPlQ29p9v4jIabHnhcgJGIwSdpzLwr/3X8bhyzdLy3vGNsIz98Ti/jahtVsxZGJtAu706cADD4jA4P33gXnzgOefB5YsEXUkCbh6teLXlCSx0qhsT4unpzgKCoCtW+98QEP5euZ1AdGOZctEu4xG6z0rXl6i96USWXnF+F9CNgBgbPfoSusSkftj8EJUxwp1enx/JA0rfktB2k2xRNhDqcDIThF4pl8cOtTVkEdVQzLWeJtl3c3IsOx1ASqdLAtABBorV8rn1TF1qtjMMTlZrFSqaW4ZAMv2XITeKKFnXCO0CmtQ4+uJyL0weCGqIzcLdFh14DLWHLyM3MISAECQryee7NUUT/WJRWhdp+uPjy8/b0WlEgECALzxBvDXv4phI/Pn8/PF+Y0blV9vjaenSOdfU1FRFQctkiSGswDR61Nmzs81TTG+OZwKAJjF/YyICAxeiO5a2s1CfLHvEr4/mobiEhEIxDT2xbP3NsOYrlHw8brLLLQVKTsBt+yQjJeX9d4R01CPn1/l19cVvR7YsEGcjxplGUwBInDxv7PTdX5+uaGoxb8mQ6c3okdsQ/Rp3rhu20ZELsmmwcv8+fPxww8/4Pz58/Dx8UHfvn2xYMECtGrVqsJrdu/ejUGDBpUrT0hIQOvWrW3ZXKIaOZehwbI9F7HldCYMdybhdogMxPQBzfFA+zColHeZm6U6q4gmThT7A926BcyfL092rehavR7YvFmcP/SQyLcybFj1h3QMBuD0aXHeoUP1tgfQaoFx48R5fn754KUSZzPy8PWhKwCA2UNa1U2+GyJyeTYNXvbs2YOZM2eiR48e0Ov1mDNnDoYOHYpz587Bz9pEPzOJiYkWs42bNGliy6YSVdvJtFws/jUJO+9MIAWA/i2bYHr/ZujTvHHd3GDLriL69FMx6bXs6p/iYtFbAgD//rfoxVi1CnjxRflaU1K4I0eAnTvFcBIgBxKVDemUVVwsMvKarq/i7xiAaEPv3qJtGRkioDLn6wskJooA6ubN0tc0GiW8/eNZGCXgoY7h7HUholI2DV5++eUXi8crV65ESEgIjh07hv79+1d6bUhICIKCgmzYOqKaOXTpBhbvSsa+pOsAxNSMBzuEY8bA5mgXUYd5R6ytIpo5U/R03HuvKFu+XCSCe/hh4JVXRJmnJ9C0qdhzyMRoFMNCv/wC/Pe/cnlwsAh+AHnei1IpllSb9iFSKORl0CUl8t5HgYGiriSJ1Ufmr6FWi3p6vVy2Zg1w+LA4b93aMsOuVisCNSvB1vrj6Th65RZ8vVSY82Cbuvv+EpHLu4t1mjWXdyeZVaNGVaf27tKlC8LDw3Hfffdh165dtm4akVWSJGHvhRyMW3YQ45f/jn1J16FSKjCmWxR2zh6AxU90rdvABbC+iggQ2WjL8vICFi4Uh5eXZXI5E4PBMnABRP6VRx8Vc02GDxdfH35Y5Hjx9wdiY8Wuzv7+YufnKVPE+ejRIildmzbiNfz9gZAQoH9/cb5tmwiqzF/3+ectA7Fp00SABojnZs4s93zm+Uv4++ZzAIA/3xeP8ECfWn0ricg92W3CriRJmD17Nu655x60b9++wnrh4eFYvnw5unXrBq1Wi6+++gr33Xcfdu/ebbW3RqvVQqvVlj7WaDQ2aT/VL5IkYc+FHHy8M6k0fb+XSomx3aMwfUBzRDey4caAFa0iMv/9f+45sepHpbKc33LwoOjdML/W1INS1s2b5cuqKy8PyMys3bXmy7lzc8s9bTQY8defEqEpNqJTdBCevSeu9u0kIrdktwy7M2fOxJYtW7B//35E1XA1w8iRI6FQKLBp06Zyz82dOxfz5s0rV84Mu1QbkiRhf/J1/GvHBRxPzQUAeHsq8UTPGDzXvxnCAut4uXNFli0TvRKSJK8CsraZ4Zdfip4M8yEXg0HsFm1KCvfaa8B775W/du9eoGtXeThIrQYaNJCHjYqKxGs0bCjqGAxiPs2f/2w5F2fSpIqHjdLSKt+76OJFkfHX7PlV3R/G3Pueg7enElv+fC+aN/Gv028tETmnGmXIl+zghRdekKKioqRLly7V6vp3331Xat26tdXniouLpby8vNIjLS1NAiDl5eXdTZOpHvotOUcas/Q3KebVzVLMq5ullnO2Su/8dFa6pimyf2Py8yVJhC6SlJhovU5amiQplXI9QJJUKlHf9HjMGFF30iTLej4+klRYKJ774ANRNmmS5esHBoryCxfk91Moyr9fWlrln+XLL0U9U/0vv6zw+XOhzaSWr/0kxby6WVr1W0r1v19E5PLy8vKqff+26bCRJEl48cUXsWHDBuzevRtxcbXr/j1x4gTCw8OtPqdWq6FW13JfGCIAh1Nu4p/bE3EoRQyjeHko8WSvppgxoDlC6jqxXHV5egIffCDOY2Ot16kow+7ly/Jj0/9eVq0Sc0t69hSPTb0qNZGUVH74qaqMvkDVy7HvPJ+XkIQZR/XQ5ukwoGUTTOwdU7P2EVG9YdPgZebMmfjmm2/w448/okGDBsjKygIABAYGwsdHTMB7/fXXcfXqVaxZswYA8PHHHyM2Nhbt2rWDTqfD2rVrsX79eqxfv96WTaV66MzVPCzclog9F3IAiDktj/eMxoyBLew3PFQRLy+RHbcyFc2NadtWrOIxbcpo0qOHmGOybZvYdNEU9L/0khiiKpt/xbQH0p2/1Soz+lamiuXYUmQkXvk1C5fzriEyyAcfj+8M5d3mySEit2XT4GXp0qUAgIEDB1qUr1y5EpPvpBjPzMxEampq6XM6nQ6vvPIKrl69Ch8fH7Rr1w5btmzBiBEjbNlUqkcu5uTjo+0XsOW0mHDqoVRgXI9ovDCoBSKCXGhVS3Uy7JYVGCgnjDOpKhNvdd/vLizbcwk7zl2Dl0qJz57sioZ+1dw3iYjqJbtN2LWXGk34oXolI7cIn+xMwn+Pp8NglKBQAI90isDLQ1oipnE1kq3Zk9Eor+YJD5eT01mTnn5Xmx7WWB2/3/8SruHZNUchScB7o9rjyV4cLiKqj2py/+beRuT2bhXosHhXMr46eAU6gxjyuL9NCP4ytBXahN9FgFud9P21VVQkv2ZVmWzLDskUFgKdOonzU6dEBlsTgwFYt04MCY0bB8TUIlCoSUbeKpzL0ODFb09AkoAJPaLxRM+mdfK6ROTeGLyQ2youMWDlb5fx2e5k3C4Wydt6xTXC3x5ohW4xVSdKrJS19P2TJon5IaZeEp1OnndiPqnclJXWvK4p3b9KJTLaAuI6vV7U9/aW9xGyVhcQQYskyZNoTXULCsT7+PiIYZ6ZM8Vzr71mme3WzrI1xZiy+ggKdQb0bd4Yf3+0PfcuIqJqsWuGXSJ7MBglrD+WjsEf7saCX87jdrEercMaYNXTPfDdc73vPnCpKH2/v79l4rZPPhFl06ZZXh8ZKcovXpTLli8XZRMnisd+fiLwiIgAQkPlzRAB4Ouv5Wy35jp1EuUnTwL794vj55/lbLfp6SINv0nZbLd2VKjT49k1R5GZV4zmTfyw9Mlu8FTxnyMiqh72vJBb2XshB/N/Po+ETJFpOSLQG38Z2gqPdom8+12eTSpK3+8sVCqgXz9x/p//yOUVLa2uaqlzHdPpjZix9jj+SM9DQ19P/HtyDwT6etrt/YnI9XHCLrmFxKzbeHfLudJNExt4e2DmoBaY3DcW3p6qun2z9HQxV6TskuFz58Qk1roaNgLkoaCaDBuZ19XrxbJppRK4ccN6u03Zbu3AaJTw0rqT+OlUBnw8Vfh6ai90bdrQLu9NRM6NE3ap3riRr8VHOy7g28OpMEqAp0qBib1j8eLgFrZbbhsVJea4mOaOmJYMt2xpWa+6S5ABkZTO06z3QasFZs8W5x99JAcj1uqamCbm6vVyj8uoUfL72XCpc3VIkoS5P53FT6cy4KlSYNnEbgxciKhW2PNCLkmrN2D1gcv49H/JuK0Vk3GHtw/Da8Nb22fZs04nJrzevAn83/8BzZrV7esXFIi5KkDVq41qeq29l1bf8dGOC1j0vyQoFMAnE7rg4U4RdntvInJ+7HkhtyVJEradvYb5Pyfgyo1CAEC7iAD830Nt0atZY/s1xMtL9IjYiqcn8Pbb8nlNKJXAgAHyeVl1uNS5upbsSsai/yUBAOY93I6BCxHdFfa8kMtIyNRg3k9n8fslsQdRkwZq/G1YK4zuGuWYVPJl87zYMu+LC1u6+yIW/HIeAPC3B1rh+YHV2E6AiOod9ryQW8krLMFHOxLx1e9XYJQAtYcSU+9thhkDm8NP7aBf4S+/FHNHTHleJk4E1qwRE2YVCuCLL+4+f4obBENf7L1UGrj8ZUhLBi5EVCfY80JOy2CU8P3RNCzcloibBToAwIgOYXhjRBtENfSt4mobsrbaqCyFAjh/XuwIbZq0K0liRRAgJteaErKZViV5esp1ywZHDkwmV1sr9qfg75vPAQBm3R+PWfe3rOIKIqrPanL/ZlYockrHU2/h0SW/4fUfTuNmgQ7xIf74+tle+OzJbrYLXNLTgV27qk7aVp08L5IEtGoFvP++XHbunJhI6+8P7Nwpv8+cOaLsz38W73/kiBy4ADVPJldUJHaWbtFCtNUBPtudXBq4/HlwCwYuRFSnOGxETiXnthb/+Pk81h8XN+oGag/MGtIST/WJsW0G1rLp/k09HaacLOY9JTExok51EtXl5gK3b4usuKal1QAwdKg8vGSyfLlYumzttU3J5MLCgEuXxNG+vTycVFQkrlGrgZUrgYQEUd66tV17bSRJwj+3X8DiXWJ7gj/fF4+X74+3y3sTUf3BYSNyCgajhG8Op2LhL+ehubMP0dhuUfjbA63RpIG6iqvvUkVJ5y5fBqKjxePsbKBJE3H+3nvAm2+K4EOSqh/IWKNQAPv2Af37V/4apvb84x/AkiWizDzI6tEDOHoUWLUKeOYZhySikyQJf9+cgH//lgIAeG14a0wf0Nym70lE7oPDRuRSzlzNw2Of/Ya3Np6BpliP9pEB2PB8Xywc28n2gQtQedr8ykyYIIZ5rlwBxo+v3XtLErB9e9XBz2uvia+ffSaXWRtOysio3We5SwajhDc2nC4NXN55pB0DFyKyGQ4bkcNoikvw0fYLWHPwMoySGCL6y9CWmNgntu72IaqO+PjyvScqlZgzkp8vHvuazbP561+BWbPktP/p6RXnYpk9G3jySaB7dxGoWBMaWr5MoQD27AEeegjQaIBHHhFBVtnXMAUme/eK9ufkiF4ha5/FRrR6A2avO4UtpzOhVAALRnfE2O7RNns/IiIGL2R3kiThpz8y8ffN55BzWwsAeLhTBN58sA1CAryruNoGapo23zztv/lcGWt27646mZ1pPo15YKJQAHFxYpdq035FmZkVB1k+PuKxn59dtwDQFJdg2ppjOHjpBjxVCvxrfGc81JEJ6IjIthi8kF2l3SzEnI1nsPdCDgCgWbAf3nmkPe6JD3Zsw6ZMAYYNq1na/PT0ygMXoHpzYYqKyveoGI2iLQMHymXVDbJq81lqIVtTjEkrjyAhUwN/tQc+n9gN/Vo4+OdIRPUCgxeyC73BiJW/XcZHOy6gqMQALw8lZg5sgekDm0HtUce7PteGJAENG4qJr77VXIpdnSXTTz8NPPCAWPVTNkBRKsWqILW64h6VsqZMATp2BPbvB+65R7TXGhtvAZByvQATVxxC+q0iBPurserpHmgfGWiz9yMiMscJu2RzZ67m4dHPfsN7WxNQVGJAr7hG+OWle/HS/fHOEbgAInmcKQeLKZFcVUxzZSoza5ZYTfTFF/JSa0BeKdSyJRAZKXapNu0cbepRadxYLNXW6+XrvvgC6N1bzKXp3VsMW9nZsSu3MGbpAaTfKkJMY1/8MKMvAxcisisulSabKdTp8fHOJKzYnwKDUUKgjyfmjGiDsd2joFA4YC+iytR2F+cVK+RhnIoolWJFEiB6S65cEXlfnnhClG3ZIibmduwIfPKJPNTTrx9w4ACwYQPw6KNimKppU8seHDstgzbZ/EcGZn9/Cjq9ER0iA/HvyT3ssyKMiNwe9zYih/st+Tpe++EPpN0sAgA81DEcb49s57w3Ol9f6yuLqmIaxunZs+I6RiNw8CAwdiwQFCSCF1POGHNeXpZzXMqqbLWRHXK4LN1zER/8kggAuL9NKBY93hm+XvwnhIjsj//yUJ3SFJdg/tbz+PZwKgAgMsgHf3+0HQa3trIc2JkoFNXvbTFXUgKsXVv9+gcPimDG22xV1bBhInAqOwS1c6ecNReofEm3DZUYjHhzwxmsO5oGAHimXxzmPNjGvsvZiYjMMHihOrMrMRtv/HAamXnFAICn+sTgbw+0hr+jdn62B50OWLSo8joKBdCnD1BcDDz1lCj76it5jouHhzjKMi1/Nqnpku46kFdYgpnfHMf+5OtQKoC3R7bDpL6xNns/IqLq4JwXumt5hSV4Z/O50v2IYhr74oPRHdGrWWP7NiQ9XQytxMeLG3p6upgzAgB9+1Z+k9fpgHnzxPnbb8t5XKp6r+ho4PXXxQ7SZ86Ur1d2n6TazKux9t42XgYNAMnZtzF1zTGkXC+Ar5cKi5/o4vw9aETksmpy/2bwQndlx7lreGPDaeTc1kKhEEMKrwxtBR8vO68iKrux4pNPiuEc81/vJUuAGTPkVT86nRj28fQUX02BxbVrIrAw34zRVNfDQ7xu2U0cJ0yQr09MFNeaAgxABDqxscDWreLxc89VnJXXCew6n40/f3sCt7V6RAb54IunuqNtBP+eiMh2GLwweLE5TXEJ5m2Se1uaN/HDB2M6oVtMQ/s3xtrGihW5dElkri0sBF55BVi6VHx9913gL38RPSKrV4u6t26JCbaFhaI35sMPxbDP2rXl551cuiTS+E+ZIl7/3/8WAUxFu1U7KUmS8PneS1jwy3lIEtAzthE++1NXBPs76URrInIbXG1ENvVb8nX89T+nkJFXDIUCeO7eZnh5SEt4ezooZ0t1ksWZ/P67CC7atpWXLwNiUmyfPsCf/lT+mh49gHPnxHlenvWNDy9dEvUOHxbHY4+JoSrzDLymjRSHDbPb0uaaKC4x4LX1f2DjyQwAwOM9m2Lew+3g5cF0UETkXBi8ULUV6Qz4x88JWH1Q3PRjGvvin2M7oXtsI8c2zNoqnIqoygRYu3aJIMNcnz7Ajh3ll0xv3iwy5f70U/n3euwxkYzOxNOz4t2qP/sMmD5d5GxxEmk3C/H818dx+moeVEoF5o5siz/1jnG+fDxEROCwEVXT8dRb+Mv3p5ByvQAAMLF3DF4b3hp+zrKSyDxZnEoFjBwJbNxoWceULC4qSgwFmTY8NAU0ej2g1Yp65it9ytb97DNg5kzxnPmGigoFMHq0yILr7w+kporAylpQ5URDSHsv5ODP351AbmEJGvp6YskTXdGXexQRkZ1xzguDlzpTYjDi0/8lYfGuZBglICzAGx+M6Yj+La0kWXOkggIgIEAEErt3A926ASEhcqp/07LiioKFmqwEMhqBY8fEMXNm+eAkMVGk/V+4EPjb3yp+HTtnxy3LaBSJ5z7cnghJAjpGBWLpn7ohMsin6ouJiOoY57xQnbh8vQAvrTuJU2m5AIBRXSIxd2Q7BPo64SqZ1avlIGLQINE7kp0NZGUBaWnysuIC0XMEHx85KVxJiVxuzlrvTEmJWHnUoYMIcqz1qpw8KXpczPckssZO2XGt0RSX4C/fn8KOc9cAABN6RGPuw+0cN2+JiKgGOBOPypEkCd8fScOIRftwKi0XAd4eWPxEF/xrfGf7Bi7p6WJOSnp6xXWOHAH+7/+AF1+Uy4xGsSTa3x/4+99Fyn1TgBAZKcovXpTrL18OhIYCAwaIgMc01yU+XtQ9fVqu+/XXomz0aLmnpqzx40XgU1UeFztkx7UmIVODRxb/hh3nrsFLpcQ/HuuAf4zuyMCFiFwGe17Iwq0CHd7YcBo/n8kCAPRu1ggfjeuMCHsPJZgvMVYogJdfFnNaoqPl+SiTJ8vLmsuyNhpaVFS+h8VgkMv8/ESiufh4saPz7dui/Nq18nW3bq16R+kbNyp+zg7ZccuSJAnrjqTh7U1nodUbERnkg8+e7IpO0UF2awMRUV2wS8/LZ599hri4OHh7e6Nbt27YZ74qw4o9e/agW7du8Pb2RrNmzbBs2TJ7NLPeO5B8HQ98shc/n8mCp0qB14a3xtfP9rZ/4JKebrnEWJKAjz4CWrUSc1kAsbKnosDFXH6+WN5cUCBWEen1wLp1Yrl0QYFYQfTqq2J+zC+/AIMHi5wxMTFy8DJihBiG+t//gBdekPcl2rzZ+nsePiwCoaefLh/gKJXA99+LuS52nKxboNVj9ven8NoPp6HVGzGgZRP89OI9DFyIyCXZvOdl3bp1mDVrFj777DP069cPn3/+OYYPH45z586hqZWloikpKRgxYgSmTp2KtWvX4rfffsPzzz+PJk2aYPTo0bZubr2kNxjx8c4kLNmdDEkCmjXxw6IJXdA+MtAxDTpwoOJlzwkJIriZPLl6r7V+PXD9OrBnj1xm6mHp0kUuy86Wz41GICfH8rFpdREg9iiqzPnzIudLbKz1vYjGjq1e2+tIYtZtPP/1MVzMKYBKqcBfhrbE9P7NoeTGikTkomy+2qhXr17o2rUrli5dWlrWpk0bPProo5g/f365+q+++io2bdqEhISE0rLp06fj1KlTOHjwYJXvx9VGNZORW4SXvjuBI5dvARATN98e2c7+6f1NzIeLKrJ9O/DSSyKQsaZJExGcvPuuCDROnhRzYz79FBg3TgQze/eKeTHe3mJ5dHX+DJRKkePlwQcrr6dQAAsWAH/9q3hsp72IrPnP0TS89eMZFJcYERqgxqePd0XPOAfn5SEissJpVhvpdDocO3YMr732mkX50KFDccC0YV4ZBw8exNChQy3Khg0bhhUrVqCkpASeZfaD0Wq10Gq1pY81Gk0dtd797Th3DX/97ynkFpbAX+2B9x/rgIc7RTiuQWWHiyqyaZP13o977gH27xeTcp9+Wt5HyCQqSkzMNVdVL4q5gACgY8eqE+JJklgirVCIrQeiouwetNwuLsFbG8+UZsu9Nz4Y/xrfmWn+icgt2HTOy/Xr12EwGBBa5oYRGhqKrKwsq9dkZWVZra/X63H9+vVy9efPn4/AwMDSIzo6uu4+gJvS6g2Y99NZTF1zFLmFJegQGYjNL97j2MAFqDrNvynba0CAZRI5k+nTxVeVCvj55/LPp6TcXftyc0UQsny5Zabefv2s13/11cpXStnIidRbeHDRfmw8mSGGiYa0xOqnezJwISK3YZcJu2VTjEuSVGnacWv1rZUDwOuvv468vLzSIy0trQ5a7L7SbhZi7LKDWPnbZQDAlHvisH5GX8QGV7Gs1x5Maf6tUSrFUE9+vtgk8cgReRUQIM7HjxfPv/uu9WGg6u5/VBHTnKspU+S9jgAxj8Xa77PRKIaL7MRglLBkVzLGLjuI1JuFiAzywffTeuPF++I5v4WI3IpNh42Cg4OhUqnK9bJkZ2eX610xCQsLs1rfw8MDjRs3LldfrVZDreb/KKvj1/PX8PK6U8grKkGQryc+HNMJ97e1/nNwiKgoYPFi4PnnLcuVSrHKp3dvwOPOr6yXl+j9mDRJPA4MFM95eADt25cf2lGpgP79y5ebgiWjUdT517/Exolt2pS//oMP5Mfx8SJQAkRemAULymfTtWMel6y8Yry87iQOXhLLsx/qGI73RnVAoI8TJhQkIrpLNu158fLyQrdu3bBjxw6L8h07dqBv2c3w7ujTp0+5+tu3b0f37t3LzXeh6jEYJXy4LRHPrDqKvKISdIoSw0ROFbiYzJghMuJu2yaCja5dRRCxaJGYWGtOrxfLpVevtsxmW3Zox7TKp0eP8uXLl4v9jnbtEsuXX3xRpPa3dn2zZvJ7KBRi1ZKfnzj/61/FdgCmYMiOeVx+OZOF4Z/sxcFLN+DrpcIHYzri08e7MHAhIrdl89VG69atw8SJE7Fs2TL06dMHy5cvxxdffIGzZ88iJiYGr7/+Oq5evYo1a9YAEEul27dvj2nTpmHq1Kk4ePAgpk+fjm+//bZaS6W52sjS9Xwt/vztCRy4KP5H/lSfGMx5sA3UHk6eTdV8r6EWLYCgIDFsZD7XpahIBDiA/Fx6upg7Ex8vyq2t8qnu6p/arBKy48qi28UleOenc/jPMTGvpn1kABZN6IJmTSrI/EtE5MScZrURAIwfPx43btzAO++8g8zMTLRv3x5bt25FTEwMACAzMxOpqaml9ePi4rB161a8/PLLWLJkCSIiIrBo0SLmeKmFY1du4vmvj+OaRgtfLxXmP9YBj3SOdHSzLJkHG+Y3+zvBLADg0iWxzLnsnBVJAt57TwwheXmJZdZTp4pyhULs7lw2EVxJiTyZt18/ESSFhIjH2dmWKf1rs0rITiuLDqfcxOzvTyL9VhEUCmBa/+Z4eUi88welRER1gLtKu6mvD13B3E1nUWKQ0CLEH0uf7Ir40AaObpZgCliOHhUrckzBxtKlwJgx4ti9u/x1AQEiWy4gNkd88EFg507xODERaN3acqKuQgGkploGE+Y9OomJ8l5HQNW7STsBnd6If+28gGV7LkKSgKiGPvhoXGfmbiEil+dUPS9kX1q9AXM3ncW3h8WqqxEdwrBwTCf4qe34o66oNwWoOAmdJImlzqblztZoNOL6KVOAf/xDDlwA4NCh8iuMJAk4eBB47DGRz0WhAL76Sn6+TRtg2TKR0ffyZbEXkRMHLxeu3cZL351EQqbIZTSmWxTeHtkWDbw5t4WI6hcGL27kmqYY09cew4nUXDGHdFgrzBjQvNJl6XXOPDhRKkVvypNPivMbN6qXhA4QgYa1TsFp08RqIPOemVatgJs3K36to0fFSqWGDeWeG0C0Y/p08T6SJNq4fLld9xyqDr3BiC/2peBfOy5AZzCioa8n5j/WAQ+0D3d004iIHMIueV7I9o5duYmHPt2PE6m5CPD2wMrJPfD8wBb2DVzS08WcE1NwYgoO/P3FxNqEhOoFLs2bix4RazlfDAbg9Glg4EC5LDERmD27fF2lUmzGaNpG4Nat8u9vNMpBktEogiMHJJarSHJ2PsYsO4gFv5yHzmDEoFZNsO3l/gxciKheY8+LG/jP0TS8seE0SgwSWob6Y/nE7o5JOnfggPWhG5MTJyq/3tTbolaL4ZuKAp0RI4AyW06U9vSYrlEoRC9KVJRYbg2IfYx0usoDKINBrBayczr/cs0wSvj3/hQs3J4Ind6IBmoP/N/IthjTLcq+ASkRkRNi8OLCjEYJC7adx+d7LgEAHmgXhg/HdYK/Pee3VMeaNWIS7tq1luWmm7ApwHntNeD118WuzOZ1VSp5V+aQECAzEwgOLv8+RqO4rlEjMXk3Lk6Ut28vJ5T77jt5l2elUh4yMn8vOyWWq8ilnHz89b9/4NgVsVnmgJZN8I/RHRAeaGVLBCKieoirjVxUgVaPWetOYsc5kSL/xcEt8PL9LR2bBj49HWjatOIVPyUlYujGlF05MVFkpz1zRiSAa95cTKz1L5OnJDERyMgQQUWjRuL1r18X15TNgnv5ctW9Jua5WLZtk4MZU2I5B815MfW2fLg9EVq9Ef5qD7z1UBuM6x7N3hYicntcbeTmMnKLMGX1USRkauDlocTCMR2dI39LVJTIrWI+Ydc0dAMAnp4ijf9LL4nHMTFiiKhssBEcLK7fu1cMHzVtKrLemvPzE5OBZ8yQU/tXN6OteS6WKVPEBGA7JZarSEKmBq+t/wOn0sWE4nvjg/GP0R0RGcTeFiKistjz4mJOpeXi2TVHkXNbi2B/L3w+sTu6xTR0dLNkxcXA448Dt2+LSbe1GYIxLbX29xfDPdaWXJvXdXDgcTe0egMW/5qMpbsvQm+U0MDbA28+yN4WIqp/2PPipv6XcA0vfHMCRSUGtA5rgC8ndUdUQ19HN8uSwQBs3CjOw2uxIsZaHpjKljDbKaOtLRy7chOvrj+N5GwxH2dYu1C880h7hAZ4O7hlRETOjcGLi1j7+xX8349nYJTEBM4lT3Z1vom5gAherJ1XR3q69TwwpiXMw4ZZBiqmZdMA0KGDvJGik8vX6vHhtkSsPngZkgQE+6vx90faYXgHLn8mIqoOJ7z7kTmjUcLC7YlYuvsiAGB892i8O6o9PFVOmqLHPICwFkyYp+cvm44/KaniZczWljAXFwNdulh/LSe17WwW3v7xLLI0xQCAsd2iMOfBNgjy9XJwy4iIXAeDFyem1Rvwt//+gR9PZgAAXr6/Jf58n50Tz9WUWi0m0mZmAjk5lQcUV69aTsSNj7fM1WLO2hJmhQKIiJDPndjV3CK8/eNZ7EwQq8OaNvLFe6Pa4974Jg5uGRGR63HS/75TgVaPZ1cfxY8nM+ChVGDhmI546f545w5cAGDlSmDmTOCdd8TS5xUrRLlWK3pdvvpKzpzbpo38PCB6VZYvL99jU9FKIl9fEQBdvSrOnZDeYMSX+y5hyEd7sDPhGjxVCrwwqAW2v9yfgQsRUS1xtZETyi3U4elVR3AiNRc+niosm9gNA1q6wI0uPV0sf7aWe2XePODLL8vvWWQtN4tpBZGfnwh4XHQl0am0XLz+w2mcu7ORYo/YhnhvVAe0dJbdvYmInAhXG7mwa5piPLXiMBKv3UagjydWPt0DXZs60VLoylibs2Kaq2JSNla2NpfFhVcQASL4XLgtEd8cToUkAYE+nnhjRGuM7Rbt2CSCRERugsGLE7lyowB/WnEIaTeLENJAja+m9EKrMBf6X7q1OSumuSqLFwN//asYKrL2fG0UFwMTJ4rzr74Sexc5kNEo4T/H0vCPn8/jVmEJAGBUl0jMebANgv3VDm0bEZE74ZwXJ5GcnY+xyw4i7WYRYhr7Yv2Mvs4TuKSnA7t2Vb3bclSUCFJMzOeqqNVicq75nJaaZMW1xmAA/vtfcaSm1u416siZq3l4bOkBvLr+NG4VlqBlqD++e643/jW+MwMXIqI6xjkvTuDCtdt44ovfcT1fh9ZhDbBmSk+ENHCSRGXmSeOUSuDTT4FJk8QEWdPkYZ1O7Fvk6Ql4eckZciMjxVG27uXL4mjbVg5cCgrEVx8feUKv6XU9PETwY2Kq+9VXwPPPi6GoyhLZ2VDazULM/zkBP5/JgiQB/moPzLo/HpP6xjrvcnYiIidUk/s3/3V1sIRMDSYsF4FL2/AAfDO1t/MELmWTxhmNYiWRvz+g18v15swRZXPmiMdRUcA99wCtWonyvDy57vvvi/Iff7TscQkKEnUzM+WyTz4RZdOmWbYrMlKUz5wpz6ExJbKrqneojhiMEp5dfQT3frALW0+LwOWRzhH4318G4Nl7mzFwISKyIc55caCzGXn405eHcKuwBB0iA/HVlJ7OlayssqRxzqCiycE2nux76NINzPvpXOkqIgBYMak77msTatP3JSIigcNGDnI2Iw9PfHEIeUUl6BQdhDXP9ESgj6ejm2WpoqXP586JybnWho2MRmD0aNEjYppEa22IqaKhoOoOG129an3yb9ll13Xoam4R3t+agC1/iN6hAG8P9IxrhE8f7wofL9fYmoCIyFlxqbSTS86+jadWHEZeUQm6Ng3Cqmd6IsDbyQIXQAQBn34qhmcAeYKteVZcQMxz8brTY1RQAGzdKs69vctn2DWva85aJt7K6pom/06bJnpc7nbybyWKdAYs23MRy/ZchFZvhFIBPN6zKWYPaYnGnIxLRGR3DF7sLPVGIZ788hBuFOjQPjLAeQMXkylTgEOHxLyVf/5TZM2tjJeXyLJrOrd124YNE0NFNkhkZzRK2HQqAx/8ch4ZeWIvol5xjfD2yHZoG+G8vXpERO6Ow0Z2lJlXhLHLDiL9VtGdpbR90MjPiea4UKlDl27gva0J+CNdTDaODPLBnAfbYHj7MOffooGIyAVx2MgJ3SzQ4ckvDyH9VhFiG/ti7ZReDFycUMr1AszfmoDt58QGiv5qD8wY2BxT7omDtyfntRAROQMGL3ZQpDNgyuojuJRTgMggH3w9tTdCApxkOXRdMxiA06fFeYcO5TdZdFK3CnT45H9JWPv7FeiNUum8lpeHtGSSOSIiJ8PgxcYMRgkvfXcCJ1JzEejjidXP9EBkkI+jm1V9BQUiBwsA5OZan1hrrrgY6NJFnOfnV13fwYpLDFh14DI+25UMTbHIXTOoVRO8MaIN4rmBIhGRU2LwYkOSJGHuprPYfu4avDyU+HJSd7QIccEbonlCuqooFEBEhHzupPQGI9YfT8e/diQhSyMm47YOa4A3H2yLe+KDHdw6IiKqDIMXG/pi3yV89fsVKBTAJ+M7o0dsI0c3qeZ8fOSstT7V6DHy9RU5WJyUJEnYfu4aFm5LRHJ2PgAxGfflIS0xqkskVNz1mYjI6TF4sZHdidmY//N5AMCbD7bF8A7hDm5RLSmVIuFcUpL4Wp3lyKa9jeLjbZ7ttiYOXbqBBb+cx/HUXABAkK8nXhjUAn/qHcPJuERELoTBiw1czMnHi9+egCQBj/eMxjP9Yh3dpNoruzGj+eaH1rLiLl8OzJhhvb6DnE7Pwz93JGJ3Yg4AwNtTiWfvaYbnBjRz7hw7RERkFfO81PX7F5fg0SW/4VJOAbrHNMQ3U3vDy8NFN+mraHsAUwr+oCCRvO7CBdHLkp4ONG0qb5ZYtr6dJWbdxr92XMAvZ7NEU5QKTOgRjZfui3ff1V5ERC6KeV4cRJIkvPL9KVzKKUB4oDeW/qmb6wYugPWNGSvb/NA0tFTd+jaScr0AH++8gE2nMiBJYt7wo50j8dJ98YgNdu7VT0REVDWb3VkvX76MKVOmIC4uDj4+PmjevDnefvtt6HS6Sq+bPHkyFAqFxdG7d29bNbNOffX7FWw/dw2eKgU+n9gNTRq4eH6Q+Hh5OMhEpRKp+AExMTc/X94yoKr6NpZ2sxCv/vcP3P/RHvx4UgQuIzqEYfus/vjX+M4MXIiI3ITNel7Onz8Po9GIzz//HC1atMCZM2cwdepUFBQU4MMPP6z02gceeAArTfvjAPCy9R45deBchgbvbkkAALw+vA06RgU5tkF1ISqq8s0Py+Zwqaq+jaTeKMRnu5Px32Pp0BtFz8/g1iGYPaQl2kcG2vS9iYjI/uw652XhwoVYunQpLl26VGGdyZMnIzc3Fxs3bqzVezhizkuhTo+HPt2PSzkFuK91CL6c1N299r9JT6/Z5oc1rV9LV24UYPGvyfjhxFUY7gQt/Vo0xuwhrdAtpqHN3peIiOqe0855ycvLQ6NGVec62b17N0JCQhAUFIQBAwbgvffeQ0hIiNW6Wq0WWq229LFGo6mz9lbXwm2JuJRTgNAANRaO7eRegQsgApCaBCE1rV9Dl3LysXhXMn48mVEatNwbH4yX7otHd1fMpUNERDVit+Dl4sWL+PTTT/HPf/6z0nrDhw/H2LFjERMTg5SUFLz11lsYPHgwjh07BrW6/ByS+fPnY968ebZqdpWOp97CqgOXAQAfjOnEzRZtKDHrNpbuTsamUxm4E7NgYKsm+PN98ejalD0tRET1RY2HjebOnVtlsHDkyBF079699HFGRgYGDBiAAQMG4Msvv6xRAzMzMxETE4PvvvsOjz32WLnnrfW8REdH22XYSKs34KFF+5GUnY/Hukbio3Gdbfp+9dXx1Fv4bFcydiZkl5bd1zoEL94Xj87RQY5rGBER1RmbDhu98MILmDBhQqV1YmNjS88zMjIwaNAg9OnTB8uXL6/p2yE8PBwxMTFISkqy+rxarbbaI2MPy/dcQlJ2PoL9vfDWg20d0gZ3JUkS9iVdx2e7k/H7pZsAxJLn4e3DMGNAC3SI4kRcIqL6qsbBS3BwMIKDq7dx3dWrVzFo0CB069YNK1euhLLsMtpquHHjBtLS0hAe7lzp9XNua7Fsz0UAwFsPtUVDDhfVCYNRwrazWVi6+yJOX80DAHgoFRjVJRLTBzZH8yb+Dm4hERE5ms3mvGRkZGDgwIFo2rQpPvzwQ+Tk5JQ+FxYWVnreunVrzJ8/H6NGjUJ+fj7mzp2L0aNHIzw8HJcvX8Ybb7yB4OBgjBo1ylZNrZWPd15Agc6ATtFBeLhThKOb4/IKtHr852ga/v3bZaTeLAQg0vg/3rMppt7bDBFB1dgUkoiI6gWbBS/bt29HcnIykpOTEVVm5Yn5NJvExETk5Yn/YatUKpw+fRpr1qxBbm4uwsPDMWjQIKxbtw4NGjSwVVNr7GJOPr47kgYAmDOijfutLrKjbE0xVh+8jLW/pyKvqASA2DBxYu8YTO4bi8b+Lp7oj4iI6hz3NqqF13/4A98eTsN9rUOwYnIPm7yHu0vMuo0v913CjyczoDOILQhiG/tiyj1xGN0tCr5e3LmCiKg+cdo8L+7gZoEOPxy/CgCYNqC5g1vjWgxGCTsTrmHVb5dx8NKN0vLuMQ3x7L3NMKRtKFRK9mIREVHlGLzU0DeHrkCrN6JDZCB6xDK3SHXkFZZg3dFUrDl4Bem3igCIHZ6HtQvFs/c2Y44WIiKqEQYvNbThhOh1mdw3lnNdqnA+S4PVB65gw4l0FJeIoaGGvp6Y0LMp/tQ7BpGchEtERLXA4KUGkrPzcTGnAJ4qBYa0C3V0c5xScYkBW09n4utDqTh25VZpeZvwADzdNxYPd46At6fKgS0kIiJXx+ClBrafywIA9G0ejABvTwe3xrlcysnHN4dS8d/j6cgtFKuGPJQKDG0Xikl9YtEzrhF7qoiIqE4weKmBU2m5AMQmgCR6WbadzcK6I2k4cFGegBsZ5IPHe0ZjXPdohAR4O7CFRETkjhi81MCFa/kAxBBIfSVJEv5Iz8P3R9Ow6VQGbhfrAQBKBTC4dQie6NUUA1qGcNUQERHZDIOXajIYJVy5UQAAaBFS/1LUX8/XYuOJq/j+aFppEAeIXpYx3aIwrkc0J+ASEZFdMHipJoNRgvFOOr/6MuG0QKvH9nNZ2HgiA/uTr8Nw5xug9lBiePswjOsejd7NGkPJXhYiIrIjBi/VZH5/drOkxBZKDEbsvZCDjSczsONcVukSZwDoHB2Esd2j8FDHCAT6cMIyERE5BoOXalIpFQjy9URuYQmu3ChEkK/77CKt0xtx4OJ1/Hw6C9vPZeHWndVCABAX7IdHOkfgkc6RiAv2c2AriYiIBAYv1aRQKNAlOgi7EnNw9MotdIoOcnST7kpxiQH7kq7j59OZ2JFwrXTiLQA0aaDGyI4ReLRLBDpEBnKJMxERORUGLzXQr0UwdiXm4OtDVzC5b6zLrajJua3FrsRs/JqQjX1JOSjQGUqfa9JAjQfahWF4+zD0atbY5T4bERHVHwxeamB8j2h8+msyLuUU4JvDqZjYO8bRTaqU0SjhXKYGv57Pxv/OZ5fmqTEJD/TG8PbhGN4hDN2aNuTEWyIicgkMXmqggbcnZg5qjve3nse8TWcR19gP9zhRwjpJkpB6sxD7k6/jQPINHLh43WL+CgB0iAzE4NYhGNw6BB0iAxmwEBGRy2HwUkNT722G01c1+OlUBiavPIxXH2iNKffEOSQIMBolJGXn49iVWzieegu/X7pRumuziZ+XCv1aBOO+NiEY1CqEGW+JiMjlMXipIYVCgYVjOsJolLDldCbe25qAjSevYsbA5hjePtxmc0UkSUL6rSKcy9TgXIYGJ9JycSL1lsVEWwDwVCnQJboh+rUIRr8WjdEpOgieKqVN2kREROQICsnNkpZoNBoEBgYiLy8PAQG2S+MvSRK+OZyK97YkoPDOxNdgfy8MaBmCQa2boENkIKIa+tY4mCkuMSD1ZiEuXy9A6s1CXLlRiMSs20jI0pQLVADA10uFztFB6Nq0IbrFNkTP2EbwUzMmJSIi11KT+zeDl7t0q0CH1QcvY/WBy+Xml6g9lIgL9kOwvxoBPh5ooPaEt6cSJUYJBoOEEqMRBVo9buTrcLNAhxsFOuQVlVTwTqJXJT6kAdqEB6BjVCC6xTRE67AG8GDPChERuTgGL3YMXkx0eiOOXrmJXeez8VvyDVzMyYdWb6z6QisCvD0Q09gPTRv7IraxL5o38UfbiAA0b+LPISAiInJLNbl/c3yhjnh5KNG3eTD6NherjwxGCVdvFeHi9XzkFupwu1gPTVEJtHojPJRKeKgU8FAq4OulQmN/NRr5eaGxnxeaNFC7VfZeIiKiusbgxUZUSgWaNvZF08a+jm4KERGRW+EYBBEREbkUBi9ERETkUhi8EBERkUth8EJEREQuhcELERERuRQGL0RERORSGLwQERGRS2HwQkRERC6FwQsRERG5FAYvRERE5FIYvBAREZFLYfBCRERELoXBCxEREbkUt9tVWpIkAIBGo3FwS4iIiKi6TPdt0328Mm4XvNy+fRsAEB0d7eCWEBERUU3dvn0bgYGBldZRSNUJcVyI0WhERkYGGjRoAIVC4ejmQKPRIDo6GmlpaQgICHB0c+o9/jycC38ezoU/D+dS334ekiTh9u3biIiIgFJZ+awWt+t5USqViIqKcnQzygkICKgXv3yugj8P58Kfh3Phz8O51KefR1U9LiacsEtEREQuhcELERERuRQGLzamVqvx9ttvQ61WO7opBP48nA1/Hs6FPw/nwp9Hxdxuwi4RERG5N/a8EBERkUth8EJEREQuhcELERERuRQGL0RERORSGLzYUWxsLBQKhcXx2muvObpZ9cZnn32GuLg4eHt7o1u3bti3b5+jm1QvzZ07t9zfQVhYmKObVW/s3bsXI0eOREREBBQKBTZu3GjxvCRJmDt3LiIiIuDj44OBAwfi7NmzjmlsPVDVz2Py5Mnl/l569+7tmMY6EQYvdvbOO+8gMzOz9HjzzTcd3aR6Yd26dZg1axbmzJmDEydO4N5778Xw4cORmprq6KbVS+3atbP4Ozh9+rSjm1RvFBQUoFOnTli8eLHV5z/44AN89NFHWLx4MY4cOYKwsDAMGTKkdN84qltV/TwA4IEHHrD4e9m6dasdW+ic3G57AGfXoEED/i/TAT766CNMmTIFzz77LADg448/xrZt27B06VLMnz/fwa2rfzw8PPh34CDDhw/H8OHDrT4nSRI+/vhjzJkzB4899hgAYPXq1QgNDcU333yDadOm2bOp9UJlPw8TtVrNv5cy2PNiZwsWLEDjxo3RuXNnvPfee9DpdI5uktvT6XQ4duwYhg4dalE+dOhQHDhwwEGtqt+SkpIQERGBuLg4TJgwAZcuXXJ0kwhASkoKsrKyLP5W1Go1BgwYwL8VB9q9ezdCQkLQsmVLTJ06FdnZ2Y5uksOx58WOXnrpJXTt2hUNGzbE4cOH8frrryMlJQVffvmlo5vm1q5fvw6DwYDQ0FCL8tDQUGRlZTmoVfVXr169sGbNGrRs2RLXrl3Du+++i759++Ls2bNo3Lixo5tXr5n+Hqz9rVy5csURTar3hg8fjrFjxyImJgYpKSl46623MHjwYBw7dqxeZ95l8HKX5s6di3nz5lVa58iRI+jevTtefvnl0rKOHTuiYcOGGDNmTGlvDNmWQqGweCxJUrkysj3zLvIOHTqgT58+aN68OVavXo3Zs2c7sGVkwr8V5zF+/PjS8/bt26N79+6IiYnBli1bSof26iMGL3fphRdewIQJEyqtExsba7XcNGM8OTmZwYsNBQcHQ6VSletlyc7OLvc/TLI/Pz8/dOjQAUlJSY5uSr1nmleRlZWF8PDw0nL+rTiP8PBwxMTE1Pu/FwYvdyk4OBjBwcG1uvbEiRMAYPGPBNU9Ly8vdOvWDTt27MCoUaNKy3fs2IFHHnnEgS0jANBqtUhISMC9997r6KbUe3FxcQgLC8OOHTvQpUsXAGLO2J49e7BgwQIHt44A4MaNG0hLS6v39w0GL3Zy8OBB/P777xg0aBACAwNx5MgRvPzyy3j44YfRtGlTRzfP7c2ePRsTJ05E9+7d0adPHyxfvhypqamYPn26o5tW77zyyisYOXIkmjZtiuzsbLz77rvQaDSYNGmSo5tWL+Tn5yM5Obn0cUpKCk6ePIlGjRqhadOmmDVrFt5//33Ex8cjPj4e77//Pnx9ffHEE084sNXuq7KfR6NGjTB37lyMHj0a4eHhuHz5Mt544w0EBwdb/EesXpLILo4dOyb16tVLCgwMlLy9vaVWrVpJb7/9tlRQUODoptUbS5YskWJiYiQvLy+pa9eu0p49exzdpHpp/PjxUnh4uOTp6SlFRERIjz32mHT27FlHN6ve2LVrlwSg3DFp0iRJkiTJaDRKb7/9thQWFiap1Wqpf//+0unTpx3baDdW2c+jsLBQGjp0qNSkSRPJ09NTatq0qTRp0iQpNTXV0c12OIUkSZKjAiciIiKimmKeFyIiInIpDF6IiIjIpTB4ISIiIpfC4IWIiIhcCoMXIiIicikMXoiIiMilMHghIiIil8LghYiIiFwKgxciIiJyKQxeiIiIyKUweCEiIiKXwuCFiIiIXMr/A8SxQOoe+63+AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mu = np.array([0, 0])\n", "Sigma = np.array([[2.0, 1.2], \n", " [1.2, 1.0]])\n", "\n", "Lambda=np.linalg.inv(Sigma)\n", "\n", "x=[10]\n", "y=[15]\n", "\n", "for i in range(50):\n", " mm = mu[0]-Lambda[0,1]/Lambda[0,0]*(y[-1]-mu[1])\n", " ss = 1.0/np.sqrt(Lambda[0,0])\n", " x.append(np.random.normal(mm, ss)) # sample new x\n", " y.append(y[-1]) # keep the same y\n", "\n", " mm = mu[1]-Lambda[1,0]/Lambda[1,1]*(x[-1]-mu[0])\n", " ss = 1/np.sqrt(Lambda[1,1])\n", " y.append(np.random.normal(mm, ss))\n", " x.append(x[-1])\n", "\n", "plt.plot(x,y,'r:.');\n", "gellipse(mu, 4*Sigma)\n", "\n", "print(\"Mean:\", np.mean(np.c_[x, y], axis=0))\n", "print(\"Cov:\\n\", np.cov(np.c_[x, y].T))\n" ] }, { "cell_type": "code", "execution_count": null, "id": "fa8ea70b", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "b639b80a", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.13" } }, "nbformat": 4, "nbformat_minor": 5 }