{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Name: **Your name here** \n", "UID: **Your student ID num here**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Optional Homework: MCMC \n", "In this homework you will create a loss function for a logistic regression. Unlike your previous homeworks, where you \"solved\" for the optimal regression parameters using gradient optimization, in this assignment you create a confidence interval for the slope of the separation line between two classes." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from utility import *\n", "import numpy as np\n", "from numpy.random import randn, rand\n", "import matplotlib.pyplot as plt\n", "np.random.seed(0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create a classification problem in two dimensions\n", "The two classes will be separated by the line\n", " $$w^Tx = 0$$\n", "where $w$ is a 2-vector. The slope of this line is given by $m=-w[0]/w[1]$." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The NLL of the initial guess is 111.4339285159692\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAerElEQVR4nO3df6wd5Z3f8ffH1/FKznYFOIY4gH8k8rZLpC0LV4TVbrZKG7bGVeukVVrQFXFptC4OtGS1ldaUfyJVrgjZJCUtgTpZVLO+DUu7SWNVZClB7faPXbK+sMThxzoYgo2DBQ5kIapXAcy3f8xcOPdw5vycOfPMnM9LOjpnfp3z3Llznu+Z5/nOM4oIzMzMellVdwHMzCxdDhJmZlbIQcLMzAo5SJiZWSEHCTMzK7S67gKU6T3veU9s3ry57mKYmTXKww8//OOIWN9rWauCxObNm1laWqq7GGZmjSLpWNEyNzeZmVkhBwkzMyvkIGFmZoUcJMzMrJCDhJmZFXKQMLNqLC7C5s2walX2vLhYd4lsDK1KgTWzRCwuwq5dcPp0Nn3sWDYNsLBQX7lsZD6TMLPy3Xzz2wFi2enT2XxrFAcJMyvf8eOjzbdkOUiYWfk2bhxtviXLQcLMyrd3L6xdu3Le2rXZfGsUBwkzK9/CAuzbB5s2gZQ979vnTusGKiVISNom6Yiko5L29FguSV/Olx+WdEk+/29KerTj8aqkz+TLPivpRx3LtpdRVjObkoUFePZZePPN7NkBopEmToGVNAfcDlwBnAAOSToYEU90rHYlsDV/fAi4A/hQRBwBLu54nx8B3+zY7ksR8XuTltHMzMZTxpnEZcDRiHgmIl4D7gF2dK2zA7g7Mg8BZ0na0LXO3wOejojCIWvNzGy6yggS5wPPdUyfyOeNus5VwNe75t2QN0/dJensXh8uaZekJUlLp06dGr30beQrXc2sJGUECfWYF6OsI2kN8I+A/9ax/A7gA2TNUSeBL/T68IjYFxHzETG/fn3PGyvNluUrXY8dg4i3r3R1oDCzMZQRJE4AF3ZMXwA8P+I6VwKPRMQLyzMi4oWIOBMRbwJfJWvWskF8pauZlaiMIHEI2CppS35GcBVwsGudg8An8yyny4FXIuJkx/Kr6Wpq6uqz+DjwWAllbT9f6WpmJZo4uyki3pB0A3A/MAfcFRGPS7ouX34ncB+wHTgKnAauXd5e0lqyzKh/2fXWt0q6mKxZ6tkey62XjRuzJqZe883MRlTKKLARcR9ZIOicd2fH6wCuL9j2NLCux/xryijbzNm7d+Xom+ArXc1sbL7ium18pavVzdl1reL7SbTRwoKDgtXD95FoHZ9JmFl5nF3XOg4SZrOk6qYgZ9e1joOE2ayYxoWWvo9E6zhI1MEde1aHaTQF+T4SreMgMW0eNsPqMo2mIGfXtY6ySxjaYX5+PpaWluouRn+bN/e+2G3TpmzMfbOq+NizApIejoj5Xst8JjFt7tizurgpyMbgIDFt7tizurgpyMbgIDFt/jVndfItRW1EDhLT5l9zZtYgHpajDh42w8wawmcSZmZWyEFiVvmCPjMbgoNElVKtiH1Bn5kNyUGiKilXxB6p08yG5CBRlZQrYl/QZ2ZDcpCoSsoVsS/oM7MhOUhUJeWK2Bf0mdmQHCSqknJF7Av6zGxIpQQJSdskHZF0VNKeHssl6cv58sOSLulY9qyk70t6VNJSx/xzJD0g6an8+ewyyjo1qVfEHp5hbKkmrZlVYeKhwiXNAT8ArgBOAIeAqyPiiY51tgP/CtgOfAi4LSI+lC97FpiPiB93ve+twMsRcUseeM6OiN/tV5ZGDBVujbactNaZk7B2bVrx32xUVQ8VfhlwNCKeiYjXgHuAHV3r7ADujsxDwFmSNgx43x3A/vz1fuBjJZTVbCIpJ62ZVaGMIHE+8FzH9Il83rDrBPC/JD0saVfHOudFxEmA/PncXh8uaZekJUlLp06dmuDPMBss5aQ1syqUESTUY153G1a/dX4tIi4BrgSul/Qbo3x4ROyLiPmImF+/fv0om5qNLOWkNbMqlBEkTgAXdkxfADw/7DoRsfz8IvBNsuYrgBeWm6Ty5xdLKKvZRFJOWjOrQhlB4hCwVdIWSWuAq4CDXescBD6ZZzldDrwSESclvVvS3wCQ9G7gN4HHOrbZmb/eCXyrhLKaTST1pDWzsk18P4mIeEPSDcD9wBxwV0Q8Lum6fPmdwH1kmU1HgdPAtfnm5wHflLRclv8aEX+cL7sFuFfSp4DjwCcmLatZGXw7EJslE6fApsQpsGYJWVzM0r6OH886bfbudXRNVNUpsDYqX41lTTbM8ZvyKMg2Ep9JTJuvxrImG/b43bw5CwzdNm3KrvC3pPQ7k3CQmDZ/eazJhj1+V63KziC6SdlQMJYUNzelxFdjWZMNe/z6gpLWcJCA6fYR+MvTLrPWvzTs8esLSlrDQWLaHWz+8rTHLHbODnv8+oKS9oiI1jwuvfTSGNmmTRHZV3zlY9Om0d9rWAcOZO8vZc8HDlT3WW2Vwj6s49hJQQr73koFLEVBveqOa3ewNU8qGWI+dqwl3HHdj/sImieV8bp97NgMcJBwH0HzpJIh5mPHZoCDRFs62GYpyyaVX/BtOXbM+nCfRBuk0kY/LbP295pVzH0SbZdKG32B0k9y/At+dLN0pmmlcpBog1Ta6Huo7FKChYVsGIg338ye6wwQn/40rF6dBazVq7PpmvSMBbN4PYeVpyg3tomPsa6TaIOE8/UTLtpbJkr737279x+4e3dFpS124EDE2rUri7F2bcRP121K/59gtcLXSbRcwm30qV9KMPGuW70azpx55/y5OXjjjdLKOYyisffOsIpV77jtPOn8E6x27pNou4Tb6FNJRCoycXdOrwDRb36FClsdSfyfUAF3wZTHQaItUmqj7zC1SwnGrBUm7s6ZmxttfoWK6vwvrput6zncBVOyonaoJj5mtk8icZUP9VPUGD/EB03cZ9KAPokDB2KmxltqQj9YaujTJ1F7xV7mw0FiRk1QK0wQX962e3fE3Fy28dxcLQFi2QzFgkJS78NBqrtk6eoXJErpuJa0DbgNmAO+FhG3dC1Xvnw7cBr45xHxiKQLgbuB9wJvAvsi4rZ8m88CvwWcyt/m30bEff3KMbMd17Nuwt7xxcWsD+L48azJZu/eZFrrbAy++ePoKu24ljQH3A5cCVwEXC3poq7VrgS25o9dwB35/DeA34mIXwIuB67v2vZLEXFx/ugbIGyGTdg7nmh3jo3JQ2qVq4yO68uAoxHxTES8BtwD7OhaZwdwd35m8xBwlqQNEXEyIh4BiIifAk8C55dQJpslk9YKToVplYST/RqpjCBxPvBcx/QJ3lnRD1xH0mbgV4Dvdsy+QdJhSXdJOrvXh0vaJWlJ0tKpU6d6rWJtN0mtsLgI1167MhXm2msdKBpqOd5fc002/Qd/4LPDSZURJNRjXncDcd91JP088EfAZyLi1Xz2HcAHgIuBk8AXen14ROyLiPmImF+/fv2IRbfWGLfN6MYb4fXXV857/fVsvjWKU1+rUUaQOAFc2DF9AfD8sOtIehdZgFiMiG8srxARL0TEmYh4E/gqWbOWpaqpTTYvvTTafEtW4uNcNlYZQeIQsFXSFklrgKuAg13rHAQ+qczlwCsRcTLPevp94MmI+GLnBpI2dEx+HHishLJaFfwTzhKQ8DiXjTZxkIiIN4AbgPvJOp7vjYjHJV0n6bp8tfuAZ4CjZGcFy8Nk/hpwDfB3JT2aP7bny26V9H1Jh4GPAL89aVmtInX9hCvj7GXdutHmj6mpJ1pNkvoQMI1VdAFFEx++mK4mdVy9VMpVcPn7rFmz8n3WrCn1KrSyito0076wb1b3cxnwFddWqTrGQSjzM6uszQ4ciOfmNsUZFD9kU1zNgZkYJqJXhS1VfzG6rzgfT78g4aHCbXJ1DFWe+hjk0HO//D/W8lvs4+ssJFXUshVd9SxlaalOSU2Lhwq3atVx9VITGqB79NW8m9P8e7K+mjqKOq2+kaLO4ghnGzWNg4SNrldNM+2xLZow9kJBTbmR47UUdZpJaP0CoLONmsVBwkaTSrprE8ZeKKgpn5/bWEtRx0lCG/fMY+/e7N/SS0onezaEos6KJj7ccT0FHqx/eIml24yahDZp8XfvfudnOtsoTfTpuPaZhI2mSVcs1X1xQmJnO6N240x6+ctXvpJ1Uify59uYnN1ko2nKYP11ZFwlbtRd0oQEMiuHs5usPE3oMAYP5NPDqCc2TUggs+o5SNhoEmtCKdSkZrEpGiUJrSm/B6xaDhI2uibcys0/gyfWlN8DVi0HCWsn/wwuRRN+D1i1HCSsPlVmH7X0Z3DdCVs2e5zdZPVw9tHIvMusKv2ymxwkrB5NSaVNiHeZVcUpsJYeZx+NzLvM6uAgYfVw9tHI2rzL3NeSLgeJtkr9Wzet7KPU98MI2pqwlcqYkVagaFCnJj48wF8usYHlCg26jdiktxlLfD+M8+e18c5rHjOyfvj2pTNmWt+6im/7OXEFv25dsrVP4vFrquq4Rbqt1C9IuLmpjabRw1l1G8GkYy8tLsJLL/Vedvx47c1QHlrqbW3ua2mDUoKEpG2Sjkg6KmlPj+WS9OV8+WFJlwzaVtI5kh6Q9FT+fHYZZZ0J0/jWVV3LTRro+pXjnHNqbwRPKVOp7m6btva1tEbRKcawD2AOeBp4P7AG+B5wUdc624FvAwIuB747aFvgVmBP/noP8LlBZXFzU24abRlVtxFM2mRWVD5IohkqlXb4VJq92tjX0iRU2ScB/Cpwf8f0TcBNXev8Z+DqjukjwIZ+2y6vk7/eABwZVBYHiQ5Vf+uqruUmrb2KyrduXRKN4KlUzqkEK6tXvyBRRnPT+cBzHdMn8nnDrNNv2/Mi4iRA/nxurw+XtEvSkqSlU6dOjf1HNFZRW0HVI7NV3UYw6dhLReW77bYkGsFTGVoqpWYvS1MZQaLX7c67x/ooWmeYbfuKiH0RMR8R8+vXrx9l0+arM8F8GrXcJIFuYQF27oS5uWx6bi6bXlhIphE8hRFWE4iXlrgygsQJ4MKO6QuA54dcp9+2L0jaAJA/v1hCWdul7hSZFGq5IouLsH8/nDmTTZ85k00vLqbzMz4BicTL0tTdCd9KRe1Qwz6A1cAzwBbe7nz+YNc6/4CVHdd/Pmhb4POs7Li+dVBZZq5PIoG29WQVNbbPzblXtEtbOo1T6edpIvr0SZQyCqyk7cB/IMtWuisi9kq6Lg9Cd0oS8J+AbcBp4NqIWCraNp+/DrgX2AgcBz4RES/3K8fMjQLrYUGLrVqV1RO9eHztVvLXYXweKrytfIOBYkU1xjLXHK1T9LtAylpErZiHCm8rt60X69XY3snpO63jTvhqOEhMUxW9ail3HtdpOYAuZzd1c83ROm3rhE+Fg8S0eDzk6VtYyDKaXHOUIvXMIZ9YV6SoR7uJj6Szm3xpa32mmL7Tlkyhbs4cSlgJBx0eKjwB/dJV21qzzJgqKtJeh0Ydh0udv3H89eijpIPOQSIF/cYS8k+0RimqtMquSHt9/9esiXjXu6Z/uNR1SY7PYAYo6aBzkEhB0dGewIiktWvQT8V+lVbZFWnR97+Ow6WuMwm30g5Q0kHXL0i443painrVXi64PnBWUjQb1qHfbySUslMwRzkEqj5c6soc8gCEA0wj77coejTxkfSZRJFZ/6nUsL9/UNdSmU0jk55JlH2CNmt9IY3gPokZCBKz3ujasPGnBlVa3RXp7t3jV6yT9Em05bBqy99RKWc3tTxIRDSqTb50DfupOEqlVUYFN252U8N2a1+z/PWYln5BwmM3Wb0aOP7U4mLWB3H8eNb0u3dv76LWOeCcxzGyUXjsJktXAy+THXYklDo7XT2O0XSlfjX6JBwkrH4tHX+qjop6ubI6diyLuZ08Gkk1ykzQSzHYOEhYsRSP2AaZdtpoZ2UFWYW1HCgacILWWGXdIDLZbPCizoomPhrbcZ0ip5WUYpqdrql0Vs9aR3NZCXp1/v9wx7WNzLf5apzu5qVO0/qaNzAPYWJlfVXqTDZwx7WNzpe6lmqclrtRtym6dUbR/CqU1fTSJGU1KyabbFB0itHEh5ubSpRK20ULjNNyN842/a7ILvNv6deU1LBrI0tTRhNbnS28+GK6xDSh0TbxPokm7MJl48TbaW0zimEOCf+2mExdx3VlQQI4B3gAeCp/PrtgvW3AEeAosKdj/ueBvwQOA98Ezsrnbwb+Gng0f9w5THkaESQSr3xXSLQmbtIujBjv1/U421S9X4YJAE373wyU6HegbFUGiVuXK31gD/C5HuvMAU8D7wfWAN8DLsqX/SawOn/9ueXt8yDx2KjlaUSQ8E+tiTVtF07zrKCsOq3X+wwbuFpTr7Yu4hWrMkgcATbkrzcAR3qs86vA/R3TNwE39Vjv48BitC1IdH9jihqN295oW6KmtXtPq0+i6vLO3K1PmvZrZAL9gsSk2U3nRcRJgPz53B7rnA881zF9Ip/X7V8A3+6Y3iLpLyT9iaQPT1jOevS6OqYoT7H2FIbmSDYLpMA4I4/UOVpJUYYS1HNPiSr1zSBzhl+mKHosP4DvAI/1eOwA/qpr3Z/02P4TwNc6pq8B/mPXOjeT9UksX7fxc8C6/PWlZEHmFwrKtwtYApY2btxYbbgdVdEvke6fwi09ha3KqL+yW9P8MSWzcjv2gceRzyTSaG4CdgJ/Bqzt8zn/B5gfVJ7kmpuKvm3LB1obvmk1GbaymqFm5dLMSt048O+coYOnyiDxeVZ2XN/aY53VwDPAFt7uuP5gvmwb8ASwvmub9cBc/vr9wI+AcwaVJ7kgMSvftoSl9i9owi/xWakbh+rbasI/rARVBol1wINkKbAPLlfkwPuA+zrW2w78gCzL6eaO+UfzpqRH6Uh1Bf4J8HgeUB4B/uEw5UkuSMzKty1hKXVyN+lwmIW6MbUfEHWqLEik9kguSETMxrctYeNUBFX9y1wppaVJQbtq/YKEx26qWl33SvAw38Do4+pUOVyzk2XS0sD7XdXCo8C20SwOxdnHsLcbhWoHv/XAupaqfqPAOki0kWujsVU5XLNjt6XKQ4XPGrdrjK3KC/VSbN5wq6QN4iDRRjVektz0SqfqW46mdDvvZG+XaUlxkGijad9cOdeGSmdav/ZTCKazeIMgG0NR2lMTH0mmwPZTZXpsDam3dad4NiXbOJXUy5SuIbF64XtcJ6iFvZh13qO3SbszlbyCVMph9XPHdYpaeK5f5+isTdqdqeQV1NQqaQ3jIFGXVGqKEg2qdKpsh2/S7kxlqPMUs60sPQ4SdUmlpihRv0qn6k7tJu3OlH7Bd2dbQf0d6paYos6KJj4a1XGdSu/llFTdqd203ZliJ3vT9qGVBw/wl6gUa4qKjJpJM86umaHdWYm6s9OsPv2ChLObbCpGyaRpUqbSslHGh0pVndlpVi9nN1ntRmmHL8pUuvHGNNvL23ARITSrX8emx0HCpmKUTJqijKSXXnpnRfzpT9cfOJqUfttPSh3qlg43N1lyipqmepFWNpHU0SzVpmaaNjSb2ejc3GSNsHwdxbFjWQU7jO7KuY5f8G1qpklpAEJLg4OEJaGzXR+yyn85UGzaBOvWDf9evnLZrDwOEnVLYTjQBPRq1494O/vpttveWREXnW34ymWz8qyuuwAzrTvXc7k3Fmauhhk0rMby7uhsL9++Hfbvf2eqbF1XLs/Yv8xmxERnEpLOkfSApKfy57ML1tsm6Yiko5L2dMz/rKQfSXo0f2zvWHZTvv4RSX9/knImqy1pMSUYpl2/u738K1/xL3izqk3a3LQHeDAitgIP5tMrSJoDbgeuBC4CrpZ0UccqX4qIi/PHffk2FwFXAR8EtgFfyd+nXZo0Kl3Fxm3Xn1ZHq1sFbVZNGiR2APvz1/uBj/VY5zLgaEQ8ExGvAffk2w1633si4mcR8UPgaP4+7dKmtJgJpdyu39SL5RzYrAyTBonzIuIkQP58bo91zgee65g+kc9bdoOkw5Lu6miuGrTNWyTtkrQkaenUqVPj/h31cFrMCqmmXzaxVbCpgc3SMzBISPqOpMd6PAadDbz1Fj3mLWe33wF8ALgYOAl8YYhtVs6M2BcR8xExv379+iGLlIiUfz7bW5rYKtjEwGZpGpjdFBEfLVom6QVJGyLipKQNwIs9VjsBXNgxfQHwfP7eL3S811eB/zlom9ZxWkzyNm7sfQV4yq2CTQxslqZJm5sOAjvz1zuBb/VY5xCwVdIWSWvIOqQPAuSBZdnHgcc63vcqST8naQuwFfjzCctqNpa9e2HNmpXz1qxJu1XQ3V1WlkmDxC3AFZKeAq7Ip5H0Pkn3AUTEG8ANwP3Ak8C9EfF4vv2tkr4v6TDwEeC3820eB+4FngD+GLg+Is5MWFazsXUP/5H6kGfu7rKyeIA/swFGuRdGSjxYnw2r3wB/DhJmA7RplFezXjwKrNkE3L5vs8xBwmwAt+/bLHOQMBvAl7PYLPMosGZD8OUsNqt8JmFmZoUcJMzMrJCDhJmZFXKQMDOzQg4SCfP9AMysbs5uSpRvf21mKfCZRKJ8PwAzS4GDRKJ8PwAzS4GDRKI8XpCZpcBBIlEeL8jMUuAgkSiPF2RmKXB2U8I8XpCZ1c1nEmZmVshBwszMCjlImJlZIQcJMzMrNFGQkHSOpAckPZU/n12w3jZJRyQdlbSnY/4fSno0fzwr6dF8/mZJf92x7M5JymlmZuOZNLtpD/BgRNySV/57gN/tXEHSHHA7cAVwAjgk6WBEPBER/6xjvS8Ar3Rs+nREXDxh+czMbAKTNjftAPbnr/cDH+uxzmXA0Yh4JiJeA+7Jt3uLJAH/FPj6hOUxM7MSTRokzouIkwD587k91jkfeK5j+kQ+r9OHgRci4qmOeVsk/YWkP5H04aICSNolaUnS0qlTp8b7K8zMrKeBzU2SvgO8t8eiYccjVY950TV9NSvPIk4CGyPiJUmXAv9D0gcj4tV3vFHEPmAfwPz8fPf7mpnZBAYGiYj4aNEySS9I2hARJyVtAF7ssdoJ4MKO6QuA5zveYzXwj4FLOz7zZ8DP8tcPS3oa+EVgaVB5zcysPJM2Nx0EduavdwLf6rHOIWCrpC2S1gBX5dst+yjwlxFxYnmGpPV5hzeS3g9sBZ6ZsKxmZjaiSYPELcAVkp4iy166BUDS+yTdBxARbwA3APcDTwL3RsTjHe9xFe/ssP4N4LCk7wH/HbguIl6esKxmZjYiRbSnGX9+fj6WltwiZWY2CkkPR8R8r2W+4trMzAo5SJiZWSEHCTMzK+QgYWZmhRwkzAZZXITNm2HVqux5cbHuEplNjW9fatbP4iLs2gWnT2fTx45l0+B7y9pM8JmEWT833/x2gFh2+nQ232wGOEiY9XP8+GjzzVrGQcKsn40bR5tv1jIOEmb97N0La9eunLd2bTbfbAY4SJj1s7AA+/bBpk0gZc/79rnT2maGs5vMBllYcFCwmeUzCTMzK+QgYWZmhRwkzMyskIOEmZkVcpAwM7NCrboznaRTwLEpfNR7gB9P4XOawvtjJe+Plbw/Vkpxf2yKiPW9FrQqSEyLpKWiW/3NIu+Plbw/VvL+WKlp+8PNTWZmVshBwszMCjlIjGdf3QVIjPfHSt4fK3l/rNSo/eE+CTMzK+QzCTMzK+QgYWZmhRwkhiDpHEkPSHoqfz67xzoXSvrfkp6U9LikG+so6zQMsz/y9e6S9KKkx6ZdxmmQtE3SEUlHJe3psVySvpwvPyzpkjrKOS1D7I+/JenPJP1M0r+po4zTNMT+WMiPi8OS/lTS366jnIM4SAxnD/BgRGwFHsynu70B/E5E/BJwOXC9pIumWMZpGmZ/APwXYNu0CjVNkuaA24ErgYuAq3v8v68EtuaPXcAdUy3kFA25P14G/jXwe1Mu3tQNuT9+CPydiPhl4N+RaIe2g8RwdgD789f7gY91rxARJyPikfz1T4EngfOnVcApG7g/ACLi/5JVDG10GXA0Ip6JiNeAe8j2S6cdwN2ReQg4S9KGaRd0Sgbuj4h4MSIOAa/XUcApG2Z//GlE/CSffAi4YMplHIqDxHDOi4iTkAUD4Nx+K0vaDPwK8N3qi1aLkfZHS50PPNcxfYJ3/igYZp22mKW/dRij7o9PAd+utERj8p3pcpK+A7y3x6KbR3yfnwf+CPhMRLxaRtnqUNb+aDH1mNedTz7MOm0xS3/rMIbeH5I+QhYkfr3SEo3JQSIXER8tWibpBUkbIuJk3lzwYsF67yILEIsR8Y2KijoVZeyPljsBXNgxfQHw/BjrtMUs/a3DGGp/SPpl4GvAlRHx0pTKNhI3Nw3nILAzf70T+Fb3CpIE/D7wZER8cYplq8PA/TEDDgFbJW2RtAa4imy/dDoIfDLPcroceGW5ma6Fhtkfs2Tg/pC0EfgGcE1E/KCGMg4nIvwY8ADWkWXxPJU/n5PPfx9wX/7618lOJw8Dj+aP7XWXva79kU9/HThJ1lF5AvhU3WUveT9sB34APA3cnM+7Drgufy2yDJenge8D83WXueb98d78OHgV+Kv89S/UXe4a98fXgJ901BdLdZe518PDcpiZWSE3N5mZWSEHCTMzK+QgYWZmhRwkzMyskIOEmZkVcpAwM7NCDhJmZlbo/wOHZhPxa1+iQwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Create a matrix of data points and a vector of labels\n", "X, y = create_classification_problem(100, 2, cond_number=3)\n", "\n", "# Define the logistic loss function, and its gradient\n", "nll = lambda w: logreg_objective(w,X,y)\n", "\n", "# An initial guess of the minimizer (may not be close to center of distribution)\n", "# Note: I'm choosing a \"bad\" initial guess to produce burn-in samples for instructional purposes\n", "w_guess = np.array([[-10],[10]]) \n", "\n", "# Test the negative log likelihood function\n", "f = nll(w_guess)\n", "print('The NLL of the initial guess is ', f)\n", "ind = y.ravel()==1\n", "plt.scatter(X[ind,0], X[ind,1], color='blue')\n", "plt.scatter(X[~ind,0], X[~ind,1], color='red')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generate many samples from the posterios distribution\n", "Note: the NLL function above generates $-\\log(p(w)).$ " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Accepted proposals: 3536\n", "Rejected proposals: 1464\n" ] } ], "source": [ "iters = 5000 # number of MCMC samples to draw\n", "sigma = 3 # sigma for the Guassian proposal distribution\n", "\n", "# Counters to keep track of how many rejected and accepted proposals there have been \n", "reject_count=0;\n", "accept_count=0;\n", "\n", "# Arrays to store all the iterates be produced\n", "samps = np.zeros((iters,2)) # The samples of w from the distribution\n", "slopes = np.zeros((iters,1)) # The slopes of the samples\n", "nlls = np.zeros((iters,1)) # The NLL values of the samples\n", "\n", "# Run the Metropolis sampler \n", "w = w_guess\n", "for i in range(iters):\n", " # Make a proposal\n", " wp = w+sigma*randn(2,1) \n", " \n", " # The acceptance probability\n", " alpha = np.exp(nll(w)-nll(wp))\n", " \n", " # Should you accept this sample?\n", " if rand()" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Samples\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnmUlEQVR4nO3df3Ac5Zkn8O+j8diMHWLZYEg8tiLjdcziCCzQYvm8twu+EBNcCQqEGC+uStVumdpc9u6ArO+ksy/AlSlrow1wf+xeFdSmausQxhDMBNZwCizkts6HRGQkI7xBwQYje5wFgy2StQd7NHruj5kez4/unu751T2t76fKZannR78aaZ5+53nf93lFVUFERMHU5HUDiIiodhjkiYgCjEGeiCjAGOSJiAKMQZ6IKMBmed2AXJdeeqm2trZ63QwiooZy4MCBj1V1kdltvgryra2tGB4e9roZREQNRUQ+sLqN6RoiogBjkCciCjAGeSKiAGOQJyIKMAZ5IqIA89XsmnLFRuLoGxjHickEFjdHsG3DSnS1R71uFhGR5xo+yMdG4ujZO4ZEMgUAiE8m0LN3DAAY6Iloxmv4dE3fwHg2wBsSyRT6BsY9ahERkX80fE/+xGTC1XErTPkQURA1fE9+cXPE1XEzRsonPpmA4kLKJzYSr1IriYi80fBBftuGlYiEQ3nHIuEQtm1Y6fg5mPIhoqBq+HSNkVKpJNVSrZQPEZHfNHyQB9KBvpL8+eLmCOImAd1NyoeIyI8aPl1TDdVI+RAR+VEgevKVqkbKh4jIjxjkMypN+RAR+RHTNUREAcYgT0QUYAzyREQBxiBPRBRgDPJERAHGIE9EFGAM8kREAcYgT0QUYAzyREQBxiBPRBRgDPJERAHGIE9EFGAM8kREAcYgT0QUYDOi1PCO2Bh2Dx1DShUhEWxesxQ7u9q8bhYRUc0FPsjviI3hicGJ7Pcp1ez3DPREFHSBT9fsHjrm6jgRUZAEPsinVF0dJyIKkpoHeRG5WUTGReSwiHTX+nyFQiKujhMRBUlNg7yIhAD8DYCvA7gKwGYRuaqW5yy0ec1SV8eJiIKk1j356wEcVtX3VPU8gKcA3Frjc+bZ2dWGLZ0t2Z57SARbOls46EpEM0KtZ9dEAeSOcB4HsCb3DiJyN4C7AaClpaUmjdjZ1ZYX1GMjcazrfRUnJhNY3BzBtg0r0dUercm5iYi8VOuevFniO2/EU1UfU9UOVe1YtGhRjZuTDvA9e8cQn0xAAcQnE+jZO4bYSLzm5yYiqrda9+SPA8hNfi8BcKLG5zQVG4mjb2Ac8clE0W2JZAp9A+PszRNR4NQ6yP8SwAoRWQYgDuBOAH9S43MWMXrviWTK8j4nTII/EVGjq2mQV9UpEfkLAAMAQgB+oqqHanlOM30D47YBHgAWN0fq1BoiovqpeVkDVX0RwIu1Po+dUr30SDiEbRtW1qk1RET1E/jaNQBwUbgJieS06W1Rzq4hogCbEUH+3JR5gG8SYH/3+jq3hoiofgJfuwYApi3K1FgdJyIKihnRkw+JmBYkC4lkp1ZyYRQRBdGM6Mlb1anpvGIBF0YRUaDNiCBvVb/m6CeJoqmVxsIoIqIgmBHpGqC4fg0ALOveZ3pfLowioqCYET15K1YLoLgwioiCYkYH+W0bViISDuUd48IoIgqSGZOuMWPMouHsGiIKqhkd5IF0oGdQJ6KgmtHpGiKioGOQJyIKMAZ5IqIAY5AnIgowBnkiogCb8bNr/GJHbAy7h44hpYqQCDavWVq0QjeoZvLPTlRrDPI+sCM2hicGJ7Lfp1Sz3wc92M3kn52oHpiu8UhsJI51va9iWfe+vCCXa/fQsTq3qv6sfsaZ8LMT1QN78h6IjcTRs3es5ObiZjXwg8bqZ5wJPztRPbAn74G+gfGSAR5AtjRykFn9jDPhZyeqB/bkPeC0lLHVZideqNXg6OY1S03TVX762YkaGXvyHihVytjY1MQvA4/G4KiRQjEGR3fExip+bqsNXfzysxM1OlEf5T47Ojp0eHjY62bUnFlOPhIOYddtbb4slra850XLPXKP7LrFgxYRUS4ROaCqHWa3MV3jgUYrcczBUaLGxSDvkUYqcRwSsQzoRi+fi5iI/Ik5eSrJbhC0Fnl6IqoeBnkqyWxw1GqCIxcxEfkL0zXkyM6utrxUTGv3PtP7MU9P5C/syVNZuIiJqDGwJ09lsVrE1HnFAtvHxUbidZlVVOl56tVOolpjT57KsrOrDeuWLyw6/ubEp4iNxE0fY6wPiE8moADikwn07B2zvH+5Kj1PvdpJVA8VBXkR6RORd0TkLRF5TkSac27rEZHDIjIuIhsqbin5ztFPisszJJIp9A2Mm97frGaP3f3LVel56tVOonqotCf/MoCvqOrVAH4NoAcAROQqAHcCWAXgZgB/KyKhCs9FPmNVg8fseGwkjriL+9erXbV4PJGfVBTkVfXnqjqV+XYQwJLM17cCeEpVz6nq+wAOA7i+knOR/1jV4Ck8bqQ/3D5PrdtVq8cT+Uk1c/J/CuClzNdRALkTpo9njlGAbNuwEpFw/ge0SDiEbRtW5h2zK61sdv96tatWjyfyk5Kza0TkFQBfMLlpu6r+LHOf7QCmAPQbDzO5v+kEahG5G8DdANDS0uKgyeQXTmvw2KU5rIqyVTK7pbBdzXPDUAXu3TOKvoHxks/l9OfKbeP8SBgiwOTZJGfjkK9UXIVSRL4L4M8B/DtVPZs51gMAqror8/0AgAdU9XW755opVShnmnW9r5rm46PNEezvXl903KpK5+3XRfHaOyddBf5Sz1XYrnXLF6J/61rL5+obGEd8MoEmAaZt3jp+ripKwWNXhbKiIC8iNwN4GMAfq+rJnOOrADyJdB5+MYB/BLBCVW23Q2KQDyarQLtkwUV496Mz2WNGgLW6KAjyPw46CaRWz2XHLNA73bIxV7Q5gtZLIth/5FTec9/R0cJPAFRVtQzyhwHMAfBJ5tCgqv555rbtSOfppwDco6ovmT/LBQzywVWYfpk7uykvwBvWLV+I/3fklHluz4TVpwHDsu59jp8r19HejfntKuNiYcXuU0C5nzJoZqtZPXlV/T2b2x4C8FAlz0/BUVha2ar2zf4jpxBtjjgOqKWmNS528VyVnMcNuzRPIplC/+CE6YVp/5FTuOvx1xnoyRWueCXfMZvdYqXUtMZyZ8Qs73kxr2yy2+mTTttvxu6TR27qh8gJBnnyna72KK5tmV/yfpFwCDdeuQjrel/Fsu59WNf7alHpga72KBbMDbtuQ259/NhIHGfOTZV+UIZIetYQkR8wyJMnzOre5B4ffO+05WMF6Vz87ddF8eyBeMkaM/d/Y1VRz9rpH37/4AR69o5hMpHMb4NNsU3V9MXF6meslNuNWWIjcdsLIQUbgzx5on/r2qIgmDuwaFeX/v3ejdjfvR6vvXPStMbMD54+mBfQutqj2HVbGyLhnD93SZ8v2hyx3AAFSKdOzGbULJ5fOn3Tv3VtTd5gbjZmYbE1Yqlh8ozdAKLVvrK59eqtBkONxxkBDQCGPziFRHI6e59pTee38wK/CycmE1gwN4zTZ5NFtzUJsheXaZPHltIcCWPenFmWA8ZuNmaxK7bGqZozA4M81YzdqtUdsTHsHjpWtAn4TQ//wnRqpeGicBOWde/D4uYI5kfCRWmUQkZA+5dPP7O43T4Mz5sdwpnzxT35+ZEw7v/GKmz76UEkU/lBd1phW6vHEG4CCk8fCYfwwDdXoas9mt0k3YzTFcEstkZM11BN2KUJdsTG8MTgRNEm4Fff/78tA7wACDUJzpxPZZ/vzPkphJtK70R1YjLhelvCJgG2dLbgW9ea93bPnE8PxPZ9+xrT3bCMi4vdoO/UNPDoptXZlFG0OZK3uMtuA/X79ow6SsGw2BoxyFNN2KUJ+k12lAKA356zXk26uDmCVMEE82RK8bmLZpXccnBxcwQOrgV55sxKD9Q+e8A8d51MKX7w9EEMf3DK8gJyYjKBjVd/0fIcCmRr6RjjDLm98Z1dbZbjBYWfP6zq3bPYGjFdQzVhlyYoZwWq1fOdPpu0HTg1AlrP3rdKpmZyJZIp0+0NcxmfQKzMj4QtLxKG3HEDs3SLm9fKLIfvtNiawSqNRo2LQZ5qwmqlabkrUO0eV1jTJpeR/rh3z6jrc1YiEg5BxHxmTiG7gVCrAWgzggsDvrkKVxtbMdJohtyLGAN942KQp5pI956Li5Jt27AS258bMx3MtDKrSXDjlYsse81W/fMmudCTrVZ5A6d23dbm6sJyYjJhOphqtWG6GQVwz55R3LNnFCERXLFoLt47eTavV97xpYWWvXqrqZm7h44xyDewiksNVxMLlAWL1QyQ2EgcP3jmYFGO3U4kHHJVAdIgSAf4G69chD2/PFY0E6YWjKJpboqamRUtawLw8KbVGP7glGU9G7cKz5NbydOqnhBQXLCN/KVmVSirjUF+5oiNxLHtmdGiKYSA+dTCSoVDglRKy5q37oZdFcnynq8Jt1+3xHFvvhzRzAX4nhKfPBbMTU8b5fx6/2GQJ9+66/HXi+qtuyk13KiMXHupzUfqwfi04+SiFA4J+r59DQO9z9Ss1DBRpcxWvdrtJDV59ryrfL5fGYOpXgd4wN14RTKlXC3bYDhPnnzHbm73Q99qq6iMLxVrvcTdwiiulm0s7MmT7ziZ2/3gC4dM68bQBZdfPBsf/2vSdgpmcyRsW/HTDFfLNhYGefKlUnO7P6v2yGzAbOlswc6uNsRG4rYDqiLuCp6FQ+J6tazTOjtUGwzy1DCMYFGdWSvlTclsFE8OTeC5N+Mlxy/cfhoyG3Q1Gzw3xloKN0AvtcKXqo85eWoIuQXPqiGRTGHe7ODm9qcVNRmgLhXggQt70QL2NYyoPtiTJ18q/Ih/9vxU1XveQZil4zWrPWeN4yx17D0GefKd9EKpg0hOX9j8g7xXznaGdjWMqD6YriHfeeD5Q9kAT/5x9JOE620DWerYewzy5DuldnsyRMIh2005wk0oWWuenItPJnDvntHsRuKxkTjmzDIPIeuWL8ym3BLJVPb3ULgxCtUe0zXUkATA7ddFbWu6pGdZ8hNBNSmQfc2fPRDHuaniqazrli/EHR0tebNqUqrZHjwDfH2xJ0++Y9c7N+QGG6q/J4cmTAfCo80R9G9dy1k1PsIgT75z/zdWIRximsXPrIZMjFkzTmfVxEbiWNf7KpZ178O63ldd5/ypNFahJF/KnULZ5GJ3JPJWSATTqhCL6poL5oYx8sOvASheKFX4PNx60DmWGqaGVmppPjUWo369k9XLRnkGsmcX5JmuId/rao+iOVI6T0/1F22OQJDeccopo7SBk/UPVlsSknMM8tQQHvjmKtP51pdfPNujFhEA7O9ej9+7bJ7ruvi50yrtGGk65u7LxyBPDaGrPYrbr4tmA0NIBLdfF8XQ9ptwtHcjolxB6YkV/3Uf3v3oTFmPTamiVJgPieTVLVJc+CTAQO8Mgzw1hNhIHM8eiGd7dilVPHsgnn2jsxaKNyqt+FzqA8DmNUs5HbNCVVkMJSJ/CaAPwCJV/ThzrAfAnwFIAfiPqjpQjXNR8JmVrj36ScL0jf6Dpw8CcLeFHflf7uya1u59pvfh79uZioO8iCwFcBOAiZxjVwG4E8AqAIsBvCIiX1ZVlv0jW1ala62kVNGzdwy3XxfFswfiga4RPxMc7d1YdCxkMYWWJSucqUa65hEA/xn5n7xuBfCUqp5T1fcBHAZwfRXORQFnF9CtJJIpvPbOSSxZcFENWkRes1ojwbUTzlQU5EXkmwDiqnqw4KYogNy5T8czx8ye424RGRaR4ZMnT1bSHJrB4pOJsgcAyd+sBtU52O5MySAvIq+IyNsm/24FsB3AD80eZnLM9LKrqo+paoeqdixatMhd64koMKymw7JccWVKBnlV/aqqfqXwH4D3ACwDcFBEjgJYAuBNEfkC0j33pTlPswTAieo3n4KmnI0pmJltLCGRooB++cWzMbT9JtP7W02fZTVLZ8pO16jqmKpepqqtqtqKdGC/VlX/BcDzAO4UkTkisgzACgBvVKXFFGj9W9cWBfpSqymZmW0snVcswND2m/DoptXZFbOzQiHLee+lps+SvZrUk1fVQyLyNIB/BjAF4PucWUNO9W9dm/d9bCSO//LsW5a1y49+kuB0ugZibPT95sSn2dlQxgInoHizcLt58uzNl1a1xVCZHv3HOd8/pKrLVXWlqr5UrfPQzNPVHsX4zq8X9fDXLV+I/q1rTXO25G/7j5xyvMCJm4FXhjtDke/siI1h99AxpFTzFsUU9vANRm+OlSobn1ng5mbglWFZA/KVHbExPDE4kZd/fWJwIruv6I7YGJb3vIjW7n1Y3vNi9nhXe5QDsAFgFrg5u6Yy7MmTr1iVln1icAJD732SNxfeuAAAwM6uNpY2aHBWgburPYrhD07lfbq7tmU++gbGce+eUSzO1Kdnft4ce/LkK3arGK0WOxkXhm0bVnLbwAY1q0kwZ1YT7t0zWlRK2Gx2zf4jp1iV0iEGefKVcuqRGG/+rvYo+r59jasNLMgfpqYVk4mkadA2m11TiFUprTHIk69sXrO09J0K5F4YutqjrjewIP/JDdpOZ9Hk3o+bjFzAIE++srOrDVs6W1w9ppwLA/mfMb7idBaNcT9uMpKPQZ58x2mgD4lwo+eA2xEbc7QOInfQlpuM5OPsGvIlI3AbMyoKRcIh7LqtzXRGxbrlC8sqWUz+s3voWPZv4cEXDuH02SQAYG64CbNnhfBpIlk0u6bcxVNmm9VYrc1oJOzJk2/t7GrDkV23mJaUteuZLVv0uVo3jeokpYrW7n247+lRfJoJ8ABwNjmNc1PTeGTTauzvXp93sbdK79ilfaw2q7nr8dcr/Am8x548+Z7V3HeznpmxmKpQJNyERKUbkpJnzAbTrerXbNuwEj17x/JSNqUWT1l98jMCfSP38BnkyddiI3EIzCtNmvXMrBZTMcAHU3wygXW9r+LEZKIobdM3MG563C2rHn6jBHoGefK1voFx0wAvgGnPjFvCzTzGJ73CSpa1XAHbSGM+zMmTr1kNlimKS9IC7hdTXX7xbBzt3VjWZiXkP4lkCveYrJotJci/fwZ58jWrwTKr/T2t5sxbhf4Pf3ceO2Jj6N+61vX8fPIvJ3PjcxdMHf0kUbRb1YrL5tW6mXXBIE++5rYCoTHHPjeoz5sdst09ysjjd3wpuL25mchqBlZsJI72//5z3LNnNG/B1Ie/O593v+OnP7MM9I3U82dOnnytnEG0ji8txLMH4tnZFWfO29c9MfL4D75wqEqtJr8oLHXwwPOHMJlI2jzigkQyhbPnp4vWXXB2DVGVuR1Ec1LQqtC63lezC20oOApLHbj9uzgxmcD+7vW1aFrdMF1DgVPOtnCsQx9MN165CEB5F34gGLtPMchT4AThjUnV8dyb6YHXci78Qdl9ikGeAsdqsLaRBsuoOs6cT2FHbKzkhb85EsaWzhZEmyMQpGdvWdVGajTMyVPg2A3W5m4STjPD7qFj+PF3rrHd6H30/q/Vr0F1JuqjP/aOjg4dHh72uhk0A7R27/O6CVRHR3s32v7Oj/ZurGNrqk9EDqhqh9lt7MkTFZgdEpxP+afzQ5Vxsgq6tXsfQiLYvGZp4PYnYE6eZiSrRS4rLpuHXz90i+UKWWo8xiroUmMyKVU8MTiBHbGx7LEgbCPIIE8z0sv33VAU6FdcNg8v33cDYiNxNJWxoTj5j7FzmFm9eCvGCuigbCPIdA3NWC/fd0PRMeONbTYwa1XymPzn83NCeOvBmwGYbwhix/jd220j2EizbtiTJ8pht2iGAb4xrLhsXjbAA+WVBY6NxMveRtBvGOSJcjTaG5jyCYB3PzqD5T0v5uXW3eobGLecW69AQ+Xnma4hyjE/EnZcwIr8x/i0ZQyilis+mcDn54Rsb8/doMTP2JMnyoiNxHHm/JTXzaAq2j10rOyVzr89Z1/rJreUsZ9n4TDIE2X0DYwjyfnxgZJSRf/WtTUraXFiMuH7WTgM8kQZzMcHU2wkjqOfJCAAFswNI9xUvemxi5sjtrNw/KDiIC8i/0FExkXkkIj8KOd4j4gczty2odLzENUaq1cGU24v+/TZJCDpgmTV0HpJxPezcCoK8iJyI4BbAVytqqsA/HXm+FUA7gSwCsDNAP5WRKxHMYh8wKx6JflbJBzCo5tW25YuKOxlJ1MKEVTld73/yClcFDYPo37pNFTak/8egF5VPQcAqvpR5vitAJ5S1XOq+j6AwwCur/BcRDXV1R7FrtvabHt5RlA52rvRcjNxqp/br0vvGua2qujps0nHm4hEmyM42rvR8kJybmra1T7E9VZpkP8ygH8rIkMi8n9E5A8yx6MAjuXc73jmWBERuVtEhkVk+OTJkxU2h6gyXe1RjN7/NTy6aXU2iBtv7sIa49s2rGSNG4/1D04gNhLHgrnmF+ZK0++5wdrqQjKtwK7b2nxbi77kPHkReQXAF0xu2p55/AIAnQD+AMDTInIFYPq3b/oKqepjAB4D0qWGnTWbqLac7ivLP1hvKWBbJ37OrCZMpRTJ6Qu/KaflKaIFm8aHREwDfUjE9T7E9VQyyKvqV61uE5HvAdir6aL0b4jINIBLke65L8256xIAJypsK5FvGNPmyN8SyWmEQ/l9TicB3qy+/OY1S00XWBlVLv2q0nRNDMB6ABCRLwOYDeBjAM8DuFNE5ojIMgArALxR4bmIfKPcjaGpvkIirtc+WI217Oxqw5bOlmz6LiSSrXLpZ5WWNfgJgJ+IyNsAzgP4bqZXf0hEngbwzwCmAHxfVfmOoIYXG4mjb2AccZ9MjyNrkXDI9YW41IDpzq423wf1QhUFeVU9D2CLxW0PAXiokucn8hMjRcMevP8J0jNvXnvnpOMLciTcZDpgWliqeN3yhejfuraaza0pFigjMmH2xj76SYIBvkEogNfeOYltG1YWXZgj4RCubZmPwfdOI6Vqu+2fWS36/UdO4a7HX2+YQM+yBkQFrN7YTNE0lhOTCTwzPJEX4OfMSvfW+7euxeY1S7MzZnYPHTMtTWxVi76cGvVeYU+eqEAjvYHJ2uxZTUW/y3NT03hmeALDH5zKmymTW5q40XLupbAnT0SBdG5q2vT4/iOnsvu4FrI63sgY5IlcqFXJWqovq9Wrhcetft+N9HfAIE9UwO6N3b91LbZ0tjh6ni2dLaxv41NOyx2Y1aJvtNk1oi4L+9RSR0eHDg8Pe90MItz08C/w7kdnst+vuGweXr7vhrz7mA3QAsibrREbidsuuydvRMJNSCTN0zmPblrt2xIFVkTkgKp2mN3GgVeiArGROI6f/izv2PHTnyE2Es978zvpzXW1RxnkfSbabF0DHgAefOFQwwV5O0zXEBWo9k4/l188uxrNoiqIhENovSRiW7/m9NlgbeTOIE9UoNo7/Qxtv6mS5lAVXdsyv6wpsn7eqLsUBnmiAlY7+pSz08+O2BiWde+rtElUBVs6WzD43umS9yvcNMbvG3WXwiBPVMBsG8BydvrZERvDE4MTlqmBcJN1xUNDbtVDqszuoWMld5AKNwke+OaqvGN+36i7FAZ5ogLGNoCV7vRTamHN1DSwv3u9ZaCPhJuws6vN9/XKG0WpAC8ANl2/tOj37PeNukvh7BoiE5Xs9GM1tbKQkf6xKqK167b08npjmb3REw2JoPOKBRg99inOnC8umHb5xbPx4e/Ol9X2mUyRfo37ByewOGdXqMXNEdO6RX7ZqLsUBnmiEowa8icmE3lvfjNOAzyAbPrHeC6zc5Q6947YWF7wN+bnu2lHI3l002oAwLafHizaDCQSbsK5qWlMV7D0x+jtG3l3wPoiXOlG3Va/u2rjYigiG2Y15I1etlmgb3U4yLpgbhgjP/xaWec26qQ7uejsiI2h32ZcoNEY2/LFRuJ48IVDRdMdI+EQPn9RqGqfZKLNEezvXu/qQu+EMV5TKF3T/mrXz83FUERlsht0K/dNHgmHcP83VpW8n9W5c4NDfDKBbT89CABF7bEKJI1sR2wMO7vSF9i+gfGiIJ9IprBw3uzsxWBd76sVlYg28u5O03dOLwZW4zWJ5HT2E0S1FmQxyBPZqOagmwCueoFOz5FMKe57ehT37hnNe/6gBXgAeGJwIrsZiNXrkxvUK90DwC7vXhjQb7xyEZ49EM9emHNTPoW/b7tB4Eo7EYUY5IlsuB10W7d8oWkuvJyiVlbnNmPkoY3AMvxB8PLxBuNnnB8JYzJRvDpVgGwJiiaBoxx9uAkoLGUTDoll3r0wlRafTJimxawCtrFZiZVqztzhFEoiG27nzFtVLQTS+Xrj312Pv17WuZ1IJFPo97AXHw7Vfl5/IpkyDfBAepbMD54+iB2xMceDsGa1ylI2DzZLpVnd2yxgl5oWW82ZO+zJE9mwm/lipbDHXu4+obnndpt2qNVA66wmwVSJyFk468ULKdWKL3TTCsu0iZuetlnANmbRPDk0UXQhqsbMnVwM8kQlVDJnHqhsn1Dj3JUOIFoplTbI9eim1Q1VUbMalxqrYO40lSaAZcDe2dWWLUddzZk7hRjkiRqA2VxtIB1ECoOZ2TErTgN8JNzUMMv4q8kqbWL1+yikKD1LptJORCnMyRM1gMJSC82RMMIhKQrmzZEw7upscZ3LL7VT0rmp6YZZxl8tdgOvhb8Pq/pCftgZjEGeqMaqtU9oV3sU+7vX4/3ejZg3Z5Zp7nvenFnY2dWGXbe1uSpsVmqAclq9WcYfKrj61CtgLZgbRt+3r7HtYef+Pn78nWuqUtSuFpiuIaqx/q1riwZfK90ntNT8fSM4OUkpOBESQeslzqd0VkO4Cei745qifPXfvPZu3taM1RASwY+/Yx/U7ZQzQF8vDPJEdVDtjZ+dzN+vZHZOoc1rluLJIfezVbZ0tjgq8Vso3CTou+Ma09tevu+GslbzRk0WLAH2ZSrcqHVuvVxM1xA1IKfz942UwpbOFtPnmTc7VLRJRqEtnS3Y2dVWVuGv/sEJzJ7lbt58tDmSDfBWm3Xs7GrLFitzan/3+mwqq9Iy0o2EPXmiBuQ2PWBWrtioehgbiePePaOmM3KizZGKKiMq0vVYmgAgs/o0JIIrFs3F4Y/O5J3TKC4Wn0xYTtXMXUHqZpP0z8+5cEH0a4+7VliFkohMq1UWpjF+/7+9hITZ0lCHos0RtF4SMV0fIAAunhPCb8+VHj8QAI9sWu06DWU2DlLrOer1YleFkkGeiACUDnixkTju2zOK3DDfBOBPOlvqWgzNzTqAQkZ1SsC8lLPx/HdlUlSNgqWGiaikUmkMuxTRPxz8jWUtmWqrVrfUrP6M8fzGRauRAr0VDrwSkWO5c8O3bViJvoFxLOveB5H0jBgr5RRac0sAywFmM6UWd5Xao7dRVBTkRWS1iAyKyKiIDIvI9Tm39YjIYREZF5ENlTeViPzCSHUYM19On00Ckl5xK0gvJjK+Nmaw1NrizCCx08VnpRZ3uZ326VeVpmt+BOBBVX1JRG7JfH+DiFwF4E4AqwAsBvCKiHxZVStflUFEnjNLdSRTinlzZmH0fvNtDZ8ZnqjpvrPG9NE7Olrw+nun8qZ8Nkn6eOH9zfaKNbhZMexnlaZrFMDnM1/PB3Ai8/WtAJ5S1XOq+j6AwwCuN3k8ETWgcnbMMqu1Xy254bhvYLxoTr9RNjhXV3sU82Zb93NL1XxvFJX25O8BMCAif430BePfZI5HAQzm3O945lgREbkbwN0A0NLiPJ9GRN5xu2OWIXcKo9XsllLMZtcokN1qz80F6FObweIgDLoCDnryIvKKiLxt8u9WAN8DcK+qLgVwL4C/Mx5m8lSmn4lU9TFV7VDVjkWLFpX7cxBRHbndMctMbiVHwDxoFAqJ4JFNq01TKcZCKasLjdlxq/v6oXpktZQM8qr6VVX9ism/nwH4LoC9mbs+gwspmeMAcj/rLMGFVA4RNbjCUrvllgcwZusc7d2IRzatLhlcN69Ziq72KKYtBkVPTCZcXYCqcbHyu0rTNScA/DGAXwBYD+DdzPHnATwpIg8jPfC6AsAbFZ6LiHyk2uUBcp8vNhLH9ufGcOZ8OpVTuEDJLl3kpuSDn6tHVktFK15F5A8B/A+kLxafAfj3qnogc9t2AH8KYArAPar6Uqnn44pXInLCLJ9frWqSjahmK15V9f8CuM7itocAPFTJ8xMRmZkJPfBqYVkDImpIM62aZLlY1oCIKMAY5ImIAoxBnogowBjkiYgCjEGeiCjAfLUzlIicBPBBzqFLAXzsUXPs+LFdfmwTwHa5xXY558c2Ad6060uqaloXxldBvpCIDFtN8PeSH9vlxzYBbJdbbJdzfmwT4L92MV1DRBRgDPJERAHm9yD/mNcNsODHdvmxTQDb5Rbb5Zwf2wT4rF2+zskTEVFl/N6TJyKiCjDIExEFmO+CvIjcISKHRGRaRDoKbusRkcMiMi4iGzxs4wMiEheR0cy/W7xqS6Y9N2dek8Mi0u1lW3KJyFERGcu8Rp5tFCAiPxGRj0Tk7ZxjC0XkZRF5N/P/Ah+0yfO/KxFZKiKvicivMu/D/5Q57vXrZdUuz14zEblIRN4QkYOZNj2YOe7pa1VEVX31D8DvA1iJ9G5THTnHrwJwEMAcAMsAHAEQ8qiNDwD4S69fq0xbQpnX4goAszOv0VVetyvTtqMALvVBO/4IwLUA3s459iMA3ZmvuwH8lQ/a5PnfFYAvArg28/XFAH6dee95/XpZtcuz1wzpDas+l/k6DGAIQKfXr1XhP9/15FX1V6o6bnLTrQCeUtVzqvo+gMO4sKfsTHY9gMOq+p6qngfwFNKvFWWo6j8BOFVw+FYAf5/5+u8BdPmgTZ5T1d+o6puZr38H4FcAovD+9bJql2c07V8z34Yz/xQev1aFfBfkbUQBHMv5/ji8/SX/hYi8lfnY7eXHMb+9LrkUwM9F5ICI3O11Ywpcrqq/AdIBBMBlHrfH4Je/K4hIK4B2pHuovnm9CtoFePiaiUhIREYBfATgZVX11WsFeBTkReQVEXnb5J9dD1RMjtVs/meJNv5PAMsBrAbwGwA/rlU7nDTV5Jhf5sWuU9VrAXwdwPdF5I+8bpDP+ebvSkQ+B+BZpPdn/q1X7Shk0i5PXzNVTanqagBLAFwvIl+p5/md8GT7P1X9ahkPOw5gac73SwCcqE6Lijlto4g8DuAfatUOB+r6urihqicy/38kIs8hnVr6J29blfWhiHxRVX8jIl9EuifmKVX90Pjay78rEQkjHUj7VXVv5rDnr5dZu/zymqnqpIj8AsDN8MFrlauR0jXPA7hTROaIyDIAKwC84UVDMr84w7cAvG113zr4JYAVIrJMRGYDuBPp18pTIjJPRC42vgbwNXj7OhV6HsB3M19/F8DPPGwLAH/8XYmIAPg7AL9S1YdzbvL09bJql5evmYgsEpHmzNcRAF8F8A789rfl5aivxYj1t5DunZ4D8CGAgZzbtiM9k2QcwNc9bOP/AjAG4C2kf6Ff9Pg1uwXp2QZHAGz3+neYadMVSM/0OQjgkJftArAb6Y/yyczf1p8BuATAPwJ4N/P/Qh+0yfO/KwB/iHS67y0Ao5l/t/jg9bJql2evGYCrAYxkzv02gB9mjnv6WhX+Y1kDIqIAa6R0DRERucQgT0QUYAzyREQBxiBPRBRgDPJERAHGIE9EFGAM8kREAfb/AVTmv0U7TrLDAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "print('NLL values')\n", "plt.plot(nlls)\n", "plt.show()\n", "\n", "print('Samples')\n", "plt.scatter(samps[:,0], samps[:,1])\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Remove the burn-in samples" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "NLL values\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABFOklEQVR4nO2dd5wcxZXHf29XOSOxEiAhFgQIRBJIYEBHzogDHDiwDca+wzhgG0dOGAzmztjY5kywAaMDDD5yzllCICRArLKEsrSSVmmVtauwad79MdMzPT3V3dXd1WF66vv5gGY7VFdXV7169erVK2JmaDQajab8qYo7AxqNRqNRgxboGo1GkxK0QNdoNJqUoAW6RqPRpAQt0DUajSYldIryYXvvvTfX1tZG+UiNRqMpe6ZPn76JmWvcrotUoNfW1qKuri7KR2o0Gk3ZQ0QrZa7TJheNRqNJCVqgazQaTUrQAl2j0WhSghboGo1GkxK0QNdoNJqUoAW6RqPRpAQt0DUajSYlaIFeJrw9bz02NbfEnQ2NRpNgpAQ6EdUT0VwimkVEdblj/YnoPSJakvt3r3CzWrk0t7Tj+49Px7cenhZ3VjQaTYLxoqGfwcwjmXl07u9xACYw8yEAJuT+1oRAR0d2E5LVW3fFnBONRpNkgphcLgHwWO73YwAuDZwbjRiKOwMajaYckBXoDOBdIppORNfmjg1i5nUAkPt3oOhGIrqWiOqIqG7jxo3Bc1zJ6N0ClZLJMF6bvRYdGV2wmnQgK9DHMPNxAC4AcB0RnSr7AGYez8yjmXl0TY1rsDCNANIaeig8P6MBP35qJh6dWh93VjQaJUgJdGZem/u3EcBLAE4AsIGI9gWA3L+NYWVSk0XrkWrZsrMVALBhx56Yc6LRqMFVoBNRTyLqbfwGcC6AeQBeBXB17rKrAbwSViYrHa2gh0N1buijTS6atCATD30QgJcoW/k7AXiSmd8mos8BPEtE/wFgFYDLwsumRqMew5SVYS3QNenAVaAz83IAxwiObwZwVhiZ0ohhLXiUUl2VlegZraFrUoJeKVoG5EZH2oauGMOUtXrr7ljzodGoQgv0MiDNNvRNzS24673FaG3PRP7sfj26AAB6do10J0aNJjS0QC8j0mhxue+DpbhnwhJ8uDj6NQqGDV2bsjRpQQt0Taw0NmUDju1p64j82dqUpUkbWqBrEoEWqhpNcLRALwPypoEUir00zw9oNFGjBbqmYsl3JunrJzUVihboZYSeu1NLmkc+mspEC/QyQosdtVBOR9cdpSYtaIGuqVh0FEtN2tACvQzQGmS46PLVpAUt0DWJII7FPYaCrm3omrSgBXo5kUK5QzHaPQorRWPLgkajFC3Qy4g0apLxLrvXK0U16UILdE3FoidFNWlDC/QyIo2mgThNLgZpLFdNZaIFekJw2mRBy5twiL8r0WjUogV6ApixaiuOue1dzG3Y7nidFuxhoUtWkw60QE8A73+xAU0t7fhoiXNMcB23Wy1JMPdoksPm5pa4sxAYLdATQBXJ7W2ZRnGeBJGq+0nNe19swKjfvY+pSzfFnZVApEqg/8ejn+O6J2bEnQ3PVOV3n3e+TgsetRQWFmkqnbqVWwAAs13MnkknVZspTljYCAC4L+Z8eMbQ0CUk9p62DnTrXB12jioCvQWdJm1Ia+hEVE1EM4no9dzfxxDRJ0Q0l4heI6I+4WUz3VS5CBbz8cN+8zZWbd4VRbZSTyF8rkaTpdwX73kxuVwPYIHp74cAjGPmowC8BOBXKjNWSeRt6JJ1admm5hBzEw9xKMmUCAu+RqMOKYFOREMAjEVWiBsMB/BR7vd7AL6qNmuVQyUHiUqCo4m2uGjSgqyGfjeAGwBkTMfmAbg49/syAPury1ZlUVXlcaOFFAmgBIRywYeLnd1FNeknLaM1V4FORBcBaGTm6ZZT/w7gOiKaDqA3gFab+68lojoiqtu4UTccESTp5aJRSzqasEYl23e1xZ2FQMho6GMAXExE9QCeBnAmET3OzAuZ+VxmHgXgKQDLRDcz83hmHs3Mo2tqapRlPE0YNvRK9LaI0+SiFxZpDLbuzOqjD360XFmaM1dtxfCb38KmCBcsuQp0Zr6RmYcwcy2AKwBMZOYriWggABBRFYCbAfw91JymGEOsyLgtppVKnD/QJIeW9g7laf7v5OVoac/gs+VblKdtR5CFRV8nosUAFgJYC+AfarJUebh5uWhRp9FoZPC0sIiZJwGYlPt9D4B71GepPHlz7jpUVxHOO2Ifz/caI/8oh2YabUPXpI9UrRSNkx/mQg7U3zHW873du2RXfspaXNJknohTqGoTusYgTIeEKNtrqmK5lCtdO2UFeqdqLWGiJC2uappkEkf90gI9AVSid0sSSNNIR5M84qhfWqBrNJqKJ0wPsyg1dS3QywBrXavflL7gXHqQokkr2oaukBWbdqIjZUsw73x3UdxZUEYSFvcM7tc97ixoYsZJQjzx2Ur86rnZntPUNnTFLNnQhDPunIRHPl4Rd1Yc8SrUtDarBqMce3bV8eU19tz00jw8N70h7mxIkWqBvnb7HgBw3auz3NCTeWrR3i6aMNCToopJq/eIFkBq0R2kJswqoCdFKwyvHU8aBVBK+16NRk+KVhI7W9oxddlm54ss9UELP7Xo8tSEIXTjGElXxNL/JHhS2PHrl+bilVlrPd2j5Y8ajHJc0pi+Lf003khLp6419JhZWeEbPie3q9VUElqga5RQXaVFmkajUYMW6BqNpuJJi6OBFugJp60jg1+/NLf4YDrqnkaTGMJYTL59d3Z/0hkrt6lP3IZUC/Q0yL2PFm/EG3PXFR1LizZhJo43Sus6BU08PFe3GrXj3sCu1nYA2T1FAeCRKdGtVE+1QDdIspU6yXmLhIovAE0SUNG3/3XiUgDAxqbszmNxzI9VhEAvZxLscVnWtLR3YM223XFnQ5MY1I3WjM4hDoGebj90PaLW2DD85rcjfd7Zf/kQV35pKL495sBIn6uRQ4WGblW+tIYeAZ8t34yVm3fGnQ1pRIuitOm3/Fja2IzfvvZF3NnQ2KCySRlpVcUwvE63hi7g8vGfAvC3mXMYuH1z0ek0ynM9QamJExX1z2irRlqdtIYeDnHbodO4yYZGkybCaJ3VMWz6Li3QiaiaiGYS0eu5v0cS0adENIuI6ojohPCyWb4sWLcDZ9w5CY9NrReeLwrgI6hVSY5Do9GkhTD0reoY2q4XDf16AAtMf/8JwG3MPBLALbm/NRbWbM16UkxZusnX/WkX5zq2uyYtGMqX0Td0qo7eACL1RCIaAmAsgIdMhxlAn9zvvgC8hQyMgHJYgNPU0u54XtTJa3tz+bBm2248MGlZ3NnQuKDShh4nspOidwO4AUBv07GfAniHiO5EtmM4WXQjEV0L4FoAGDp0qN98JpqfPTPL970L1u1wPC/SYOMW53MatiHDwMj9+8Wck+Rz26vz8e4XG+LOhsaFQX26KUsrTn3LVUMnoosANDLzdMupHwD4GTPvD+BnAB4W3c/M45l5NDOPrqmpCZzhJPLSzDVqEkpCFy/BxX+bgkvvm6IX5kiwYlP5uMhWMueMGAQAOOWQvf0nYmm/cYykZUwuYwBcTET1AJ4GcCYRPQ7gagAv5q55DkBiJ0XLRE7aTIpGnw1Z/M4LiIh71KGpbIxmptIJIY467SrQmflGZh7CzLUArgAwkZmvRNZmflrusjMBLAktlxVMguV5LAsnNJrkE596EmRh0XcB3ENEnQDsQc5OrlGMcFI0+myI0OLcHd3nlQcqmlRhYZHCRD3iSaAz8yQAk3K/PwYwSn2W0smEhY2+7kuyW58KYaUFniatJNLkoomXJAu8JOctDo645W38x6Ofx50NTQWTaoHuZJpobc/gPe1O5glmxhtzCpttJHn0EAc7Wzt8j8Q08aIm2mLxwqI4SLVANxDNXN/57iJ89591mLpMnadGGCRp4vG56Q247skZ+b8TlDVNApm3Zjvmr90edzYiw9oc4nBbrLhoiwart+wCAGzd2RZzTpxJktA0dmIxUBpnJka1pnvn6tDSruRRzEV//RhAciKbyqBCCBtJaBu6Ypy+jSGLgoYHUNkLv6hqgVJEeBVVv3h2NsbcMTFQGmGwf//ucWdBkwKs+k0c3mgVq6EbmlPQQlcdpS2TYWxo2oN9+2aFTBIEnoG18/KqoL8wo0FhbtSRFDdQTZwErwSLNzTnUoqvQqVaQ3ckr6EHQ7Wd7IEPl+GkP0zE+u17ACTL5GIlLeaEJMlzZkZreybubGgUEIdgrwiBLhI71t1F/KL6k702Oxu0csvO1tyR5ArNGDZkKWtkOudHp9bj0JvfQmPTnvAzpAmFRAfnKhe8NgBVE3qZkL9eojV0pXOi8bWCML0RzGUk85iXZ2U7dCOOviYawghxHIdgT41AP+H2CZ6uL1mm65NKsr+WvmuCexsPRPUJ+3Sr2CmrxDO7Qb17pRboinEqT8NcEFTDtt7+zvz1gdKz0iWGXU9kSfLoIYkE2cGmI8O49ZV5aNi6S2GOoqG9o7LmBLTJJQZUmVyspoLv/Z81bLw/kigsrfVURRYT8Z4RNcAgrzpj1VY89slK/PyZ2cryEwZrtu3GSX+YkF/nAQAtepI3MipWoBvE6bbopLmUgyknSatYg2AUdeOOcCci3Yprd2sHZq/eJjxXWKyS7IrxwvQGrNu+B8/Wrc4f6yiHyqwQ4xsldYOLskfUkFSJItFHs2uUVlp9DkW/M6bW131Bsb5qSuQ5AGDiwg044fcT8EGosVicC+zBj9Kz96i5rlSYPNcmlzgJWvYiDX3aii3Ca61btok+/LZd2VAETsKyc0Ls6mkR6MyMWauzk2KzJDtjP7iVV9j+51c9/BlemRXuamThK1aYQI+TZEiGOFCmoosOiWvw+u0WgS66RmLYH8dQToTKhUVxvtL23W148rNVABDqhs5upRW2CWvykk24/ulZoT7D4LnpBZNL0s1EYaF6FbkMlSvQcwQVjl68ZKyTQzLPTojsFqNig4sYXB+t5b51Vxs2NWcDjy1Yt6MkCJkdHRn2VH/c5HVaRjwAsGFHoQwTXYctqMyrXikaIYYgeXPuOpcrnfH0ydjxz8RjraAqNcpVW5Ljjic7t3HJfR/j5pfnSafr1nmlSJ4XUW71vJxJtUB30p4MWfTBoo2BnuFFQ7de6VcbSIrGo1IA3T9pme3cQ1KZt2YHnsiZamRIkwZuh+gdmRmbm1tQO+4NfLp8c/SZ8oCKb6QnRcsYLx+vRPhL3Jsk+2PYXi7LNjarTTBhuBaXpUCZGe99sQGZOIyxCmEA01duBQA8NHlFvJlxQYUwLrgtBk/LKxUr0K2NqyPD+LcHP8GUpd52MBKNAuw+pLVd+hXWSWneqifx2joyeVt2mETZ0MwL2NwWs1nPvjp7Lb77zzr8Y2q9+oxFSFJGlDJMXbZJmdOB3uBCMYUCdRc8m5tbMG3FFvz0mVk+nyFxraWiaJNLMbe8Mh+jf/e+4lS9E5YXkddJ0cbcxOLabeUdqCtJo0w3Mgy8OVdt+I4oSbVAd0KVchnIhi5zT4LaQkpjc0WGW527+/0lnq4vGxJUh2UI2oHmV/UmecciIqoGUAdgDTNfRETPABieO90PwDZmHqk8hyFh9ThwK/u2joxwQY+Xj1aqobtP2iYZFS6HdpNoSvcrTQhByispaw/cEH03RtnJ9LLFi4Z+PYAFxh/MfDkzj8wJ8RcAvKg4b8oQLv03Hdvd2lE4Lrh//fY9OPw3bwtX2XnS0D0Jf4dzcTWPgFvQJQW30gtLdvotrzKR5baUe/69woJfUSEl0IloCICxAB4SnCMA/wbgKbVZi44Ms2Olq9+8E+0ZFrqoBRHSUiYX+eTB7G2hS1DCWtmYJgFANr89pxNB77mpuQXjP1qmvA6Vkw0dCK6oGOWXZC+XuwHcAEC04uIUABuYeYngHIjoWiKqI6K6jRuD+XyrxPzRgsRE92ND37GnDY1Nexw/uFOlsrvv8vGf4sAb35TOj1fKq1kmD79COSqB+LNnZuH3by7EvDU7lKZrrq/lOqorF1wFOhFdBKCRme0CfX8dDto5M49n5tHMPLqmpsZnNv3hLGsLNSuT8V/Rdrd1lBwzP3bLrtb8b0P4n/anD3DC7ROKGuoJB/YvTiM/sSLfmMttYY4dae04ki7LduxpBwC0Z9QGCUvr93QjjveWmRQdA+BiIroQQDcAfYjocWa+kog6AfgKgFFhZjJs3EwuTrS1O994w/Nz8r+NZ2zNRVQ0f/FeXctvezIVQ3O7SdEwxV9cE4xelIY0CUEO0L7KEeNVExkPnZlvZOYhzFwL4AoAE5n5ytzpswEsZOaGEPMYCubGZQ7A71VT9zIc3mhZNON0ZyU1gErBq8mlHD19xB109PlIAtVV0X+/oH7oV6CMJ0MNMhn2v2rTw22dLB+4yLZoTddYPix8ZjwtJKrHpqn9l1VskEqVvIoxivGHpx8c+bM9jfOZeRKASaa/v602O9Ghqu/0Milq7bDLbfbfSnnnPnq8Kmzlp5+LcVJc0kn2hQ0zapSaekWsFHUrziCCyZtbofPfsufiIqoOqMNDMKrdrR34n3cXoaW9dHI6acQR+90TYbmhhlBv2joyePbz1Z7qSlD8jowjdSWO7EkJw67uem10dsG5mva0YeJC591vzHfa5UdUF9yqx98mLsGjU8KPahdWPa2r3yp97QOTluKvE5fiiU/lw9jGtbAoCGWzUlTQfvxmffGGJmzZ2So89/Tnq3HDC3Pw+Kcr/SUeIVF+ufJzrVCEueIFaSuiezPM+Pmzs/GeZTuzklguITXSO99dDAD49pgDlaYbRnbbOgQdoocmYOwC5XfD7SjxqgCX4ZyoEL/V5ty7PkL/nl0w4zfnlJzbnHMw2BxBdE4DZrlv4mUkrpqUa+jysVL8FrpoxMfMWLFpp+C40992E6biEUC58ZnNxgbWbfn8ktQyMdcz3wuLooptE1IhBlFc7DT0fJIKyuX56WXnpGdLygW6PeZqYNYIPbsteoiH7qR5Wm3ATtcmZTLVSz4uH/9pKHkwBF2Q1b5Wwipfr6KnqI5G2GMF6TyEboum/yeRXz43W2l6eT/0pMZySTv/8+5izFy1zde9Qg0dcoLe/PfkJd421tBk8SN74tLm02JC8UpSR092rLbZ39Z53UgyXrIiBLo42mLh4PPTG3DdkzN8pS3qhe01dPd7rWmI0orLW6IkKwrqsEizTrw3iE+CCPTtu7Ori1WZqJxQL5zUC7u8xUV5ysBjnwSbaNV7ioZMmAUsSntPeweWbSy1oVsvdg6Ra08cK9CSilESKoVQaOFzfYofBvDndxYBAOY0bFeYIwsxR89kZuxsaXe9bu223bh3QjYWYFgRP71ifsc4tfWKEOhe8Fo9RN/ugUnL5O71+CwDLwK9uaUdLyia9AmlngZM02jPCRnxlmAW4l5kTxVR0SiyI4oXzD3jSUGYaFlEryib8wc/Wo4jbn0HjTv2OF43YWFj/neUuo2ToGab31FTEQKdKLsQwWwbS8QWdKZ7O1eX7vguugcoDSHgxG9enodfPDcbM1bJ+3bLoqLiCk1WHlI2BKanBV4xNTkvVa7aUkGj1EOfU+z1IRs+94056wAA610EujmJqhhHq58t35zfHMfcluNULirGD/23r84v2qBi1WbxxIcsR976Dr7xpaEYc/De0vc4bXDRpboKbR0dwnNWvFTiDbnGYd6VKUkErfxJ19D9UmVStbL+z8kwLfhBtgM1lCM709Ss1dtw6X1TMLB31/yxvXt1CZ5BH6zesguXj/8UXz52MO66fGTROeN9d8XQ5ipCQweADxcXb65hjXzoleaWdoz/aLnHLejsbehezChJMaEnQYjmbegJdosz6NFFXn/KmlwKfyfkk/tC3oae/deu75q4ILtQr7Gp0HajnEA3v4YxSb1ofVPJOYPbXvsi/ExZqBiBHprw8ZCuyLhgYNXAnLxc4hKkYQhN0bs8/PEKzF69Ter+5YIFXEH53RvhNMTOneSbW9EkW0Sd1ZLG5lDSZZars65rhYQbUMfTGIz3MUZSRe+nvVySg3DXckdvFC8auv3fpY+Vm4Apd0TlN2nRRlxy35SiY5ubW7C0sank2tdzdte73xfugCh+pksBvr+g0fkCF/75ST3WbtsNINhcTdF+pAoV0Yatu4Rxhnp0qQ6ctnhhkVyNZReTi+hohLG5iuqNk3lIT4pGgKql/la87NZlfeQ/ptbnf9u1Vy8rUUOnZA7AX0aYGX98eyEatsrPY5x+5ySc/ZePfD3PzPG3v4/xHy0PnI4Tt7wyH996ZJrt+ddmr8Xn9c7bBZrLVvX3vuCeyfj3R+s83bO0sRm3vTbfl0ue7C2GkKyykUoiF8UMM37/5oLIt1/M59WYwwnxe3mhYgR6kEIuTLxxyeSiJ+8KSybM7mF2JhdhOmWuoy9Y14QHJi3DdU/MkP4uTXvc/ZNl2NjUgr+8t1hJWk5stYlBAgA/fmomLvv7J+6JhDQR6qcsr3nsc/xjSj3qAzoTOGFo23a+5Xa7IY3/aDn+7UGJ8lSIkVej3cZhIhNRMQLdiluhMzN+++p8zFy1Nf+xPq/fisNveRu/e71gYzV66hNq+4uSCZA/p7z5SE9BHSt1u/SXjlFmbR3l3jXZI/Ib9yKerbdHs4DG/hmyMXOChM8tmFzkcxflIp5iLbxYQy+6TmvoyaM9w3h0aj2+JtCkHvp4Rf638fGuO9N9uymnD21XiUWTfl4qcVTebs/VrcZNL831dE9ahTkg3qQjzvddsG5HoPurTKNUrxQFv3PoNAqTouJrRO665mI+9r/exbw1Ia6kNeGWVyszQ1gHIqJiBbpbvZSdkDIq+J42d59TPxOoN75YKiS9tCl/2jxj3fbdnu751fNzivz85Z/l+ZayIO+CF2828nz9f4NFuzRGCH7CzqtyWxTfU0h86642PBjy/IhByR7ARSaXUr58/9Qws5OnYgW6LMwsFXMl+AIZ8UpREX5m9r00kienrcJJf5iIOQ3bHPMUVBYzM95fYL+rk0wnmVTEQcfkeWpaoXNkBF9YFLRjMdZJ+AlTXOzR5+C95WZyEZywtoUwO1AnoS0yx8RBqgW6o4nD5cvLfhKjgst8xCD5KUrHhyj1UscMj4FlG8PxSzZYuL7UDdHMXe+HP3kZFkHjsz/x2arEaPeAvA1dtLerrIAzrrKbLxAd/8NbC4r+Vm1ilN2r1k1Dj4pUC3Qn3OqYcT7DbiaX7L+1e/d0f6Zk3tyu9SIr/FRwexdK69/hVt0dudV45YjKokmCWcowX7u56bYKwvtmRxjuz8j7dttcKzos2sJQFXX1WzD85re935jLUp9u0UdWqViBbme2sHONsk8ne7JztXtRqrJ9exGkQYRBWIJE2i/Zg7121AF7+ctMgol66b+T0DW0Y7fRoci7x6sN3Q6ZTsF7tFT7h34m4dteCD1RylmHD/KYm+BIC3QiqiaimUT0uunYj4loERHNJ6I/hZNF/zjVj4yNRG9pz6B23Bt4bvpqT8+Sia/i1Bg8ubR5uDafvocH2NlrVcl36ZWDCfWDeWxqfcn8gi0BbABGVM35a7cXJXP87e97TkvGBi9jEnSbvxGfN9uX3Z8vumbNtt1C7d+KyiBmUhtCG/+aoy0am+7FMLTyMia4HsACAH0AgIjOAHAJgKOZuYWIBoaQP2VYC9cuvrSxi/hdpsUnTh+2MEwM1mBKrpXYzcgL1nuenrYKow7YC4cM6i241pvN08yOPW3o062z9wwK8DL5G2XjufXV+a7XqMhN/57ZqIJ9u3cu2ix5Y1P4O91bN6YuLKBxfjORopS9RaZ95ASh4NyYOya63m88pa0jgyoiqYB3zq7E8ov9RNdFGZbAQEpDJ6IhAMYCeMh0+AcA7mDmFgBg5mABMEJAxlQic6/Ttarc00qeECDvZuz6mXEvzsU5dwVfSm/lh4/728pPhN0oCgCWbGgqim+fTF1eDSyYx3ETrNNXbkHtuDcwfaX8knhnk0v2X3cNvfSCpY3NeHKau0urqm94yE1v4YrxcitHnZ7pxcRT5MnDpceiQlZDvxvADQDM6tyhAE4hotsB7AHwS2b+3HojEV0L4FoAGDp0aKDM+sXP6jXz6dvfWGB7ndty5eJnqtG6PU2u+nJxtDG5lPpqlSATJVGFMm3tjJIwcViEID9erQEFjZWxp82bA/iHi7Objn+0eBNGHeB9FbO1E6mS1NCrBYFYxpnWUjiaQW08xjyNvnJ5/rx+q9TlqkZ2Ii+XOEwurho6EV0EoJGZp1tOdQKwF4ATAfwKwLMkkATMPJ6ZRzPz6JqaGhV59ozIfNGvh9gsIBpazl9rv8ou7zsrY29TJLT9CWnv95R4tUh0JTWmzQcC48mVM1mosP8XBENWyy065zaB6PK3zD1mZDX0Y4b0lXiSGLNmW29aIe1NnquxoXdkGLss+5s6flOzQHcwHYWNjMllDICLiagewNMAziSixwE0AHiRs0wDkAEgv31PBDh9AOsWX/6fkcXLBIpUuoq8XGTSsyLrtmhgXlbes6v7oE82Ky/OWCN5JZSr6Nt3teGG52djV2txo77uSXUmJSuH79vH132fLd+M2nFvYFPATVsczQ/wv7Co6BlyMhGn3znJV/qeR0E2x3/76nzcO3Gprzzk006ihs7MNzLzEGauBXAFgInMfCWAlwGcCQBEdCiALgA2hZdV/4h67aAV06CgoQedFJXPj5es78xFh3zeo9cO4J4jo8N8aWZB8M5ds72sVnhedeIBwuN/+2AJnq1rwOOfriw6bux7GQZ2O1qJvrf5kBFbaPrKUjMDC663e55d+gAi8ZvMa7YlI0N5vLstio8/73Ff1aKVoi5ph0kQP/RHABxERPOQ1dyv5jjXvApwWr3lNnSUfRHjGTJui80t8gtlnEYXXjojI4zry7PWSt9j1yrsPq/18vs/cNZsrFqvCvxUvD999WhcdZJYoOfTNSXsJd9+Jsvtvqub+aaw6KdwXV6/8DBvE/fKVBXSw6xX1Y57w/2ZNmW7R7ji1SEdLv2deIHOzJOY+aLc71ZmvpKZj2Tm45hZzq8oIajqe4ydcmRsd/87eYVDfpz/LjonlTPjWnW1qmRO1CbpHS7xtn//pv0ks4jGJudd4J3y4gQRcKjAbTN7zlhIU+CWV9zdFfP58Z4dwTs4CePCOWPC0qykBLUl27UPL84EXq8odGj+J0VV2dBlV5LbnM39P4Eml6TQ3pFx3eXFSqF6lBasnR96/l7Jb7F+R1bYOGnoPRVs7WVG1YpTEQvW7cj7PDttal2E5d1FoWPNbNjhzdZ7y8syft/eG4+TZ5LojLGtXFjYaeifLneu906rOPNzPIL7ZEvsxhfnuO4IxMxYuTnY/q5xmCpUOSq8Mbdgiku8H3oSuPv9Jbjs759ghqK4wvPWiD1XDGG03WscEQLOGC724vngl6cDAPr37CKdHMM+hnWYS/8vuGcyJi3aaJMnO5NLsbhQNT9hMCskV0jHhSf5+N+mQyHbJPyO0kSrOGXyunbbbvz9w2Wu9empaYX5F7t0n5/egNP+PAmfLt/smJas2aLouGOKxUQV/9/6nOUbTV45xr8xCPToo8f4ZOH6rHDb5GOVXIbVCxkrVUT4zwsOwwcCYdine2ecPGyA49Jla+7aOzK44J7JUtd65UkfccsBeYETh2bi9nkPvektHDm42IvEqfF79ep42LTpSTY/pfe5yRrzHTJzMoVr3X3ERWe++886LFzfhAuO3AcHDOjpen32GeLjM3Od7pINzhE61263N5/l/dAlnynCu6+/t+s9ph5m4kLKRkM3dmL3shVXYdacMcKnS5gsBGBg727Cc9VVBCJvn9fJbPHSzDUY98IcqXSsjfyVWWvwa487C+XTsjnudRVjGLg9sbUjgxmrthUdczS5eBQM/23altBvOubOo4pIei1CQUMXmFwcEjH2FnUzkXnB7X0dd05Sko3iDExeIh5t+kE6JIY2uTiz0+Tgb7cjuBMZjqav7N+zC74zprbkeDVRtoE6ehwUn9tp8ajoZFHZnv68MAz2Yh56Z/566Wutud3S3Op4XpYwYpH46US87tPpZcLNyM0qiU2VH/90Jf7+4bIigewla/lJUdMA0NfCItNFfvvkIH15YaVo8fEgcfmvenga3vvCfhMVvxOXJR5C5rKzCc4VhaJTFgK93eyO5WMWO2xzC1CYZJ2woDSkDREwc9W2Eg3RiZ89M7vob7t3mNOwDcfc9i5emy12S1T55m/bdAZxu7vZsbGpBc/W2fvfO0U8Lnj9BWvwm3e22l7TntvP7eaX5+GOtxYWCYkhe/VwsaGXCn/R5ao9LT5xs5EHeJ5xZ8PW4k7w58/OLr3YBlFH+COHxWCqRIPQbdFyzRaHuqCKshDoRY3K5xL2sGW6kf6qLaUaGRGhuSWY77Xd8M2Y3J2yVLymKwrrh5vHUFx87//qcMPzc2y9U5wWgxmT72FOin77H8Whj8zP+tdj9pMWjgW3RcbUpZvwjykrCkLecRLSyeYuPnfvhCWhaZpGsnMtGz17c1v0+EwP15pDgDhuemP8a0l86y4t0AEEt0VFoaHLbHDhhFsOD5LYEcmKdZWjZ1wyZTS0Bz8Md2NeGcEm+sQbc0vh22x2NnYyuRiugmHWnI8tnTAz47JRQ6TuNb+veWHRNx76DLe99oXSuOBOz3Y65hWjnTpF2XQjTC+XK8Z/ml8p7Pa+C9fvKNm+Lgq9p0wEureSMBqw04qtrp2qsF9f8SSmH2RdEu32KHR7xbMOLw03/9jUevz9w2W299z88rzQFzfMs2hTUdAuENCi93Qzzxl98M1jD8c1/3Kg8BpzkCg/eOmIM1zcychWe9HConwaEvd7rSFO7THIBKtxpzUJt71n3XDeQtJbfmXs+eu378b5d08uWT8QxSRp2Ql0mQ54+M1vldxvLcsMc6hajB2NgkU1150xzPU+kTZ566vzhSYeM0EqkUxnMFPCR1wl732xASNueafkuB//ZuP7X3PKQbj5ohHCa16cuaYk0mFYbN3V6lnDnLSosWiC3NvzspPpxvuZO0CjzEQT2E51SjZsrZBcukFMeF7n2FTJWPN329lio7RF4JpRFgK9aMJB4nqjwhkFKNIoMuzPYyYManq5h5x163zsTgcd5n2wsNF2uzW7pJ+tW42vPTAVAHDzy3Px1wlLgmXCxAvTG9Aq1NBLcZosBOS9XIJo6eZni9YomGlpzxTc+gRKiAirHd4Lu3NB1Jxi7/zx7YUlx8QaenBhpcLk4pWg7UO8B7E4US975PqlLBYWBbWBZzKlw/SODCuL++DEAQN6uF4j83ZusueDher8bQ2Yge88mhUY9XeMFV4zeXHpczMM1K3cil88OxsvzPAWtc4vTgt57BqYbAjlIAO5eo/L4Zc0NiuJrW8XuVCW56Y34KoTD8ASwegkLFtwweSSTBu6CFFW7fqjKObyEqKjOhO0w55WvwVTl5W6W3lZjecX2Ue42fLc8mrElBGk7HhfEO8bZuBdBx/fMIT5LpvQvGIN3bnQZL+/jMeIKsxZUuFNIqXnC8rpNy/PMzJRci4swZTX0AMk77lJe3yWbXgC0/E4vb7KQ6D7+MI/e2aWq+3T68ISP8g8o1+Pzq71ym9erUVnHZXILHyJiy6dSqvnbg8hbPMaut15WQ09l5LVKyUMRFEeRbgJ6pIwDU7PdHxOKUEF+tJG8SSnkWyQidXHPvHm2RXUrm1UIbOTm135aA09h59yeGnmGtz3gb0HCBDR8EziGZeOHOyejM/MlqxWiyG+hF9EUSpt64LguDFqadgq9kOX2RUeCFZPvNZdgv+FWsJv609BL9wuaVLY3Sq/qcnZfyndnNzsLRVl6AhVjzrxoAGmRMN9lhNlIdCderZn61b7DtsZhZeLdVGL9VUOGNBDKh9+zUNeND3vDSn6zsGDPMeunJB51sYLRLZMoxjJGcguv3fzBbeelvm2wgk+QcmKJok9baAiYKNp+zzVJoso5sq6diooH1pDdyFjU1OZGTc8Pwdfvn+qr3SjaKZuu7Xn8+DyrX2bXDwMX70M06OAiLCxqQXNpg0zbHf18bHqUbZD39TcgrfnycfAMfBjOiAiZBiO8Udk8CM7tu0qjQkkSieMfVWbir6x8uRtCeNRNuvYImlPZeLl4nzcb4yEGNzQBXmQy4RfDd1L4/Csn4dcQwnA8be/L/VMp6zY3SNrcvESS8TMfS5b8YkwAq0tXN+E1Q5rDESvdLMxkelynRWjFHYLJpxFZdewdTeen96Ar0muapXhgUkF82i0botqnmVWGvwoHaooCw3dXBAdGcY3H/oUU5ZuimQIEzZuE3f563z2PlYt0Tr8NBehtTzjLl6xCUCMNa9mN1W79wjby2mFwDSxXaAFmxncr3v+97y1wVfhGm1nq8NzneqWXXmL/NO98FwuaNoHCxtLJuajbNd+n+QlcmrhuM+HeaAsBLpZJm3b3YYpSzfjx0/NzH/4JGjastgN/0XakRm/oxAvjSNuAV6KcNWG8Epzue5qbcfBN70lvM5M2LZxUdm/MnuN4z1mzx4nl1I3bU/V5Lfdc4LGUP/V89l4/t959HOc8T+Tis6pVtDd2pYfnLJoNwegl/7nsKucRrklSZ4fsZ/zRhol3zqXebcGsm67v/0sPQl0j0IgDvlvV0zm12yybFJtb0NXkye74F+ivO5xES7mLAXpcOY0ZLV7mc8fd/vpyHBRx1FOK0VF6djXUW1yAVAs7KzmF8B/xd/crD6cpch32oxVwIbdmNo6gkyKxquyy3pdAMV5t942e7XYdCFrQ3fDzlYu6kzdPC6kvVxc8jRxYWlcfplnljzH5kFh1tsoF+ao2uDCjF2HlCgNnYiqiWgmEb2e+/u3RLSGiGbl/rswrEzaFV5Qk4vT5gN+cctKiYJeJvaiV2aVmgqimBT1+swPF2/E85YVqnaraFWZXGx3YBLk1e2R5vrQHkXwDxOnHVq6yXmYnbrdTlsfCsJJhIbP17OWi/kv6wjR7p4w8KKhXw9ggeXYXcw8MvffmwrzVcR/vVa6XyMz53u8IH6mokqsgr16dBYetw67kiTOnYTl9U/PiiwfBuLAR+Jr1+T8/a9+ZBr+9PYiqfStCvp/X3KEl+zlsRPSfib3zEk5yXPZpL3koKtgdBlmp33Mbe/mf5tD5IrcJ53o213c1sJARii/MXed+N6kaOhENATAWAAPhZsdMeZJDbMGwwqM6AcP7OX/ZgFG/gbv1b3oeJ9uWQ/RkqX4knmPojLEZWJ56Fujpa91EpJrbHYmssOqoQ/qoy4+PmBjcvEQNXNa/Rb7C0NAVLKJmycXEMR05vv9TDf+4tnZuFciomiSFhbdDeAGAFad4UdENIeIHiGivUQ3EtG1RFRHRHUbN/obStkJGpW7lavG+u3u+OrRAPx/VFV1oW5lsZBwmsxxjeinqLnbRaT0OvKybmLthlWg+3dhEx8XVU+3N5I2A8lq6FKTog5uizYJJMlSGGQqxG+7Mm4jEF6Y0YDpK7eG9iwvuAp0IroIQCMzT7ecegDAMAAjAawD8D+i+5l5PDOPZubRNTX+zBvmoedWk93baDBBPmgUERfNz7F+VFmhpUp4bhBssJF/hkuN+5eD91aSByt2WqsXk4vd9U6UCHTFDU4Y0jfq8K5SwVxy1wouTa7KVMBt1DPAYTcxc/kcN7Sf52d7aZdJ0dDHALiYiOoBPA3gTCJ6nJk3MHMHM2cA/C+AE8LKpLkgbn8za8bf05bJN5ggNvSwJyUPzG1BZjynxMslUSYXZ6wePKryRAQctk/v0uOCa1WahUo3OFFbyKLyaXfxOpKtjyrLoXO1w/I2m8dsam61dcEcVuN9/9sguCllxvaNV4z/BMN+bT/V98PTD5Z+pmzdP3fEIDx41Sj87RvH4tihQiOGUlyX/jPzjQBuBAAiOh3AL5n5SiLal5kN6/+XAZSuO1aEqPA6zJOiAWRyWOLcyPO3TjoAQEEb9OsqGYWmFFcsF4K8HdSpIbW0e/MKCVtDF3lrGAqJHarro8w7TVm62X51o8N9f/9wGY4e0hdD+xcL8CiDmQHyCp11j08giMkle6PbBG4VEc47Yh9/D/FBkFgufyKikch+83oA31ORIRGioUprewYTFmYDGAWqPi439+hSnY/aFyQ547hV6Bga2Qm1/R0nwcLS0Jta2jBx4QacedggV5NLWENGIkIngUAXPc0pD3e85SwsrVgFj98pGbuNPNp9JCi7LeLKzbvw5LRVUtfukowhby3ap6etEoYvMHh48go0CVazRi7Qg9jQA6Yj2tHJTNQmNk8CnZknAZiU+31VCPkRYteIb3opOyioCmAId+vdg0+aZDEaamlsFesPtxTV8rNnZmHDjhZ8/J9noHtnS/xxy8svsuy+ripHBPE3XLe91H/c6ZleNyi2Ck+/HZZbRE0vyGqbl9w3RTpN0abaVphLy3bci3Odb7LJatRCzGsHUjvuDVw2agj+fNkxnldvGpfL3hZ5WUT7OH+4lZ2dI78MbgXu1VZZ2K7Mait3tqG7ffewHHqMSdJ5a3Z4Nln85KmZSvJABKGGLkRhOVj9l91W+UZBkrxH3LDLatQauh+em146qmr1UP9lq2HUCwfLInxumBOCURW38Ry7YbhbIwg7DsT3H5+OS0buV/xMyzVhlVUVkbQQUFkK1k2izzpsoMLU5bCGyFVvQ5crMfZwrYGdsJI1G0UFgVy3vQOAHXvcFzR5lc9Rd21lIdDDdPdx+0CXjhyMp212vBGRX+tkSdgQWB0Zqw09d96lEUThcv+ZYNLITJjahvykqEIvF8v7dKqOXhKd8qcPig8oLuMw4+HbLd1Pog1dxiwmY+5izppsZIlaQ09YXyomVIHu8hGPO8Cbq5GhgRsmBOODGpXc6rZmPP9bJ9U6phuFD6vVvBTZ6IXkBbrKji3I3EtYqM6SbL1R2VFGJcT27tU19zz/aYTdrKKuYuUh0EOMUaS67uUjQFZZNfTi89bnn3fEPlh6+wW26Ubhh25ddFRicgmpchKRvIZeFktd/KNau/WkoSsq2qiEWE3vrlLXEQF/fkcivo9EvpNucikLgR4WvbvKWZx+ctYhjucfvGpU/vf5R2Z9TvfNxQXJV4Dcv1YbunkixnnnmOgFWcmqVpfa+eVjB/t6DqHUni2bp7QRpw1dFUn0Q7eL4GhuVzK59lr/tMlFQHj+zxLBkgD8/JxDHa8xC6Mfnj4M8247D7ddcgROPbQGl4zMCrmCDb34XcxR5pxy0qdbdBHlDLyW+81jD/f1HJHJxW6CsmGrv40+yoXvnnqQ0vS8fEKVbqhR4h6S2P5cURz9EISv1tAFhGtDD455QpOI0KtrJwzq0w3//PcT8q5xeRu6wxjYrj59/7RhsWim1p143Ddn8FeaVQKTS6fqqJtCMlDtaSNvQ1dnR6+TCFSlErea0twiXhjYuGNPUV7DqHHabVFAWB4eRKTELiwz5DPkVbvNdmVGfkQMq+mJFZucV6SFgRe/XMB/gxAt/Ve1m1C5oVoAJDggaWCMkurpYjpdb7N94wm/n+D5mV6LUy8sEhCWDzZRsMBe5nRkr7Fq6D88fVjR39edUfw3kNVgVTVMLxqgNa9hVU5GsQDvXE3YutPbJgdpQXU/9rJgpykRP392Fmat2qb24RFxvcscl8ow25599ZU9WY4yEejhpEuQsb+5fxIv19h5uRj86rzDSu6tqpIrgyMHO29Q/cVab6tBvW7W61fgMxfmIUbs2wdLbr8Qnyzf7Hpflxj8xsNGtYYuK8xen7NOGJelHOjRxVlDVzlK+etE8f6xdkQ9QVwmJpdwJLqqGBwyn8zOhi4zQiAQBvcLvpvOhfdO9nS9dbPesL5Dr26d8hq6l/pfBivMNRHgGr4jBk8fA21yEbDV4x6Dsuxu61AyxJXphQt+6OKVok7Ub96JGy88HPd/8zjH6+at2eGemAfu+2BZ0d9u7cLcOX39hKHSz+nVtVSg9+7mrmtIx3/RVDSzG7ZLXReGaVcL9IhRMcSVsqHDRkOXuHnBuh3o1rkaFx61r6/8qcK1vpteZVAfuUUfBsZCLKOcfnB66VyC3T1p4z/PLzW7lSO/cHH3DcpPzspuSFG7t5oNNUKaqQslVTsqUqCfM2KQ0vRkPlm/Hln3xY2W1ZhezDVx42UTZq/acyeLhi6z0Cgp5aIamc4s6dTfMRY/dpmsDMr5R+6L+jvGopfkAkFXQpDoeul/BBQvJnC+Vup7SFzUJ+ePbp14kpFJ5SK4zNms9hhyz3jHtlysGxm3RbvgUGFyoCJtUBOM2becqzzNMFZja5NLAL59ci0OHtjL9bo+JvusCrfFrhJxtO380KUmRctDnhchq6Ff8y8HAgCacx3dgnXZeQAZgb5f3+ATxVZ+/+WjHM/36R79it1yw82NUAV9e5THd1AhX7yQKoHerXO1lCAZo3j3+l5d3SuXIaA+WFQcU2JPu/v2dlGvNvOLOZeyWR6e2xx6t2XDYRmB3jmEDSm+8SXnydw4zPbfO01tOICw6d+zS7QPVKRYh+HEtc5mUVNYpEqgE8nFPDFbA9waqJxJROYa8UUiX+/jhvYr+jvszS1UYe54ZP2fjXKx2syTamaKI1dfO25IDE/1T9SfTpU7bRit7P0FjSGkak+6BDrkYoCYhYWSpf9SC4ts8iLoDYb271H0d3mI82KsPux2GJ2rdWQlo6Gv3LzL9RrVWEMMR0FC+zZbWhTusSpDFHsFlAtlIdDHm8LTOkEkN2wyC/Tja/s7XrtwvXjrKutzvTzTM2VSX81veNDe7nMZQKFcrB2bbDjdqPHi6aOOZJaFHU9OWxXp80TN47aLj/CeTpm0MyfKQqCfe8Q+uPCofVyvs5uAMMc9P3DvnkUC+Nihe2GJw8YSOyWWQ8sIazsBtW1Xq+u95tcyJhFFWDX7ODm+di+p64zRjVVDT6g8DxW7dw6jLG4ee7jv+PVuWPdJDZv+PbrgX48p3g/Xz4bfadg8pSwEOgDc/02xln6syd5sV/Hn3nZe/vcz3zux5HxnS0yQfh5n0GVs6HZ5O2Rgb/ebTfXsAoeO7WfnhO9d4IT5HQf0kltYZJTd4g3FI6FKjLYoUgwG9+sein5+zSkH4a7LR4aQcvTfrqqK8NevH1t8rPKqDwAPAp2IqoloJhG9bjn+SyJiIlLrOiKJ7KKC+795HO687BgM7N1NOLQq8qH2qBLJ5MHOzi7SJKzXymoOUU0k/uq84cLjfly0jDzPsET6S+qkaJiIhNCb159SdmUh67Lao0t1yDnxhgqTy12XHxM8kQB40dCvB7DAfICI9gdwDoBojWYmiic4CfPWiuM2XHjUvvjaqKy3gPHdxh5dWEpv/pijTeYCmbYkq42KEGkzl40u9mqQqWiD+3XHqYfUCM8N2au7r7zZcd0ZB5ccO/VQ8bPdsGv7bsv6+0jEeik3RIK7b/fOZWd+ktXQ73OJTRQ1AxS4W550UCx6bR4pgU5EQwCMBfCQ5dRdAG5AjNN25r0CCUDTnmKb9wVHlpooDDdAUbWbfMMZOOuw0tAAom3obrzgsJJ45l6ZvKR0r8OTh+2NHwmEph2XjtwPU8adib1MFfKnZxfML8fs3y9QHt24/5vH4aFvjfZlg7QbuSRlUvThq0cX/b1vCIuZDOw0cb+LU56+ttS86MbeLsrJtRJb5MnG2OnbvTMG91OrbBh4LbN/Gz0kvx9wEKyvHrVPvqyGfjeygjvvj0REFwNYw8yznW4komuJqI6I6jZuFG/UGgRzGxC1ByeXJpEwGbJXd5w+vKBtGhVDVEe/d9ow3BAwkJJdxbMXjqXXi2yh5g7hFMULqawcNbgvunSqQvfO3ofQnW3cTJNiAz3r8OLOXeVmCVbsNFu/fdtxQ+Umps08+d0vOQp1GXOK1zj6SeCcEfuoWcBnSeLeK44VXxcSrgKdiC4C0MjM003HegC4CcAtbvcz83hmHs3Mo2tq/A3LnXjnp6cW8goqqXCydetHZxycnYAiwsA+BS3M+MbGxz77cLWBvazmFQNzP+QatlZQETuZJnovP35/X3kTIbL5G7ZQIx+yHi7/PuZAnDxM3NkEjaR4rGVxloEfG6fZBa6xqdgP/asKF/2YX3nMwQNw2D4SE+YO+JFP3TtXY/QBxd/v1xcWlBYZu7f52zlp9D26VIe2d+w5Iwbh85vOlr7erqO+4vj9HUN7WE2NVURYavKa69k12nkCGQ19DICLiagewNMAzgTwfwAOBDA7d3wIgBlEFHzM4pFDBxUqPRFw1JC+RedFH6o9FwDKLPx/ed5wfHTDGSXXWoNFDatRE5zJcKW0G5Idtm9h9yGn+DRnSmwpJxL4Jx00oEiA7N/ffuhrjk3+rRMPKDr3k7MOKZpDmH3ruXj8mi+55uniY/bDLf86At1stHq3iUA3O+0lFjc2g9MPHYirTjwAl44Unxdx9cm1tufM2bzEJs37v3kcLjraPfSxeTTwxDUn4u2cstJqif/zwg9Ock0L8Bcvnqi00z56SL/87++MsXebNTDvJHWQTTCze64YieGDeuORbx/vOY8y7NWzC2p6y89ttef2KbCa1M46fBAW/U7s1vznrx2Nf1jyn13cWHj/qMcqrgKdmW9k5iHMXAvgCgATmfmrzDyQmWtzxxsAHMfM68PNrpivHJf1pyUA91x+bJGQFPltGw3EuoWZWUgY27kZml5hgwrGQA8VxY5eOSFp5y978TH75TWjUw4paLF9LcGh/MTO/sNXjsJT156I0QeIF1WdedjAonydPGxA/rc1nvtJBw0o+rtv987o2ikrpMcK4rdfOnI/PPXdE3Gvxc3M+l5H7Fe8nZ6143Pbdmx0bX90ribccH6xR04HM/770iNxt8+hsLUjMDbS/vbJtbaLqS48al/87RvuE4B//OrRwuNtFoE+6oD+eOenp5aUmRUiwm8uGpH/W2Z00qVTVcminP36Fjp7pw2Zf5ILyrW/aT2Eub6YO7VLRg4GEWFYTS+M2Nd568QjB/fBT846RFqZmmdyU5bF6Pwm/OI0vP/z0/LHzxhub1Xo2rka1VVUJGPi9kgqGz90M8YCGmNC5bJRWZPCScMGYOiAHpjxm3Ow7PcXov6OsThZYD8+fXgNOlURrjrpgJJzBteflZ0ENeyQvXMxYnp27YQPf3UG5v42WPjOp757Im65aASGD7IfVl976jDMvvVcfMkkNA8e2As3nD8c3zvtINTfMTYf3Mrgj189qqgRGxjxYfbt2w2Xj86WV5/u2cb5m4tG4I9fyQqT7516EB759vF48Qcn5+/tbYqPY2g9B+UaV01v+0mfP3z1qBKt9dcXHo6Thg0oufbFH56MLp2qMPvWbLnuZ5ks+8mZB+PSkfvlw9cevm+f/Ibaj37neNx1+TGYMu5MDKvpiTu+chSOHNwXX/zX+fjh6Qdj0e/Ox0/OzM4pyLiYGqYOs+no85vOxjs/PRV/+MrR+PuVo/Dc90/C4H7dMSpnnqjp3RXXnFKqPJjNFQcP7IWh/Xvg7MMH5jsoYxeqI/brY9u5Dx/UGzePPRwPfPO4vNAcvk9vTLvpLNswskYkyhNMK6HPGZEdQJvNRF86sLhTr+nVtWhy/fwj9sHQAT3w/dOG5QXypF+eXvK8nl2q8fUTsvXqvCMKIw2zELdTPl790ZiSY2Yl5olrTsTPzzk0XyfsNs4wRpJO37hLpyqs+MOFqB3QA6ccsne+PI2y6dGlU16+3PqvI4q0bSs9cqPL75scI3pYTCyqPcxcYebI/hs1ahQHYc7qbfzolBW8p62dW9o6uLW9I3+uzfQ7DNo7MvzQ5OW8p63ddxqrNu/kJRuaOJPJKMyZPeu27eblG5uLnm8upzVbd/HfJi7hjg5xfna1tPPDk5fzjt2tvG7bbr7l5bm8bVcrMzO3tnfw3IZtUvl4d/56/mDhBl6wbrun/DfvaeNXZ63hFaZ3yGQy/MCkpdywdZentLzS0tZhWy5mMpkMvzSjIV8XW9o6uHlPG09csIF/++o81zRa2zv49je+4M3NLYHyu6yxiaev3MKfLd/ML0xfzas278yfu/OdhXzv+4tt73148nJ+ddaaomP3vr+Yf/bMTNt7Pl+xmSctamRm5qY9bbyzpY2ZmbfvbnWs33X1m/mtuetKjr89bx2/N389T1m6kd//Yj0zM29pbuGpSzflr9mwYzf/5d1Ftulv2L67pI7Vb2rmRz5ezhMXbODn61bzssYm27w58cmyTfzntxfyU5+t5N2t7fzPT+qL2tJd7y3i12YXyvCz5Zv52c9X+XqWCAB1LCFjiSMMYDB69Giuq6uL7HkajUaTBohoOjOPdruuLE0uGo1GoylFC3SNRqNJCVqgazQaTUrQAl2j0WhSghboGo1GkxK0QNdoNJqUoAW6RqPRpAQt0DUajSYlRLqwiIg2Aljp8/a9AWxSmJ1yQb935VGp767f254DmNk1XG2kAj0IRFQns1Iqbej3rjwq9d31ewdHm1w0Go0mJWiBrtFoNCmhnAT6+LgzEBP6vSuPSn13/d4BKRsbukaj0WicKScNXaPRaDQOaIGu0Wg0KaEsBDoRnU9Ei4hoKRGNizs/QSGiR4iokYjmmY71J6L3iGhJ7t+9TOduzL37IiI6z3R8FBHNzZ27l0S7QScEItqfiD4gogVENJ+Irs8dT/V7AwARdSOiaUQ0O/fut+WOV8K7VxPRTCJ6Pfd36t8ZAIioPpfnWURUlzsW/rvLbGsU538AqgEsA3AQgC4AZgMYEXe+Ar7TqQCOAzDPdOxPAMblfo8D8Mfc7xG5d+4K4MBcWVTnzk0DcBKy+2O/BeCCuN/N4Z33RXYjcQDoDWBx7t1S/d65/BKAXrnfnQF8BuDECnn3nwN4EsDrlVDPTe9dD2Bvy7HQ370cNPQTACxl5uXM3ArgaQCXxJynQDDzRwC2WA5fAuCx3O/HAFxqOv40M7cw8woASwGcQET7AujDzJ9w9sv/03RP4mDmdcw8I/e7CcACAIOR8vcGAM7SnPuzc+4/RsrfnYiGABgL4CHT4VS/swuhv3s5CPTBAFab/m7IHUsbg5h5HZAVfgAG5o7bvf/g3G/r8cRDRLUAjkVWU62I986ZHmYBaATwHjNXwrvfDeAGABnTsbS/swEDeJeIphPRtbljob97JwUZDxuRzaiSfC3t3r8sy4WIegF4AcBPmXmHg0kwVe/NzB0ARhJRPwAvEdGRDpeX/bsT0UUAGpl5OhGdLnOL4FhZvbOFMcy8logGAniPiBY6XKvs3ctBQ28AsL/p7yEA1saUlzDZkBtiIfdvY+643fs35H5bjycWIuqMrDB/gplfzB1O/XubYeZtACYBOB/pfvcxAC4monpkzaRnEtHjSPc752Hmtbl/GwG8hKzpOPR3LweB/jmAQ4joQCLqAuAKAK/GnKcweBXA1bnfVwN4xXT8CiLqSkQHAjgEwLTckK2JiE7MzXx/y3RP4sjl8WEAC5j5L6ZTqX5vACCimpxmDiLqDuBsAAuR4ndn5huZeQgz1yLbZicy85VI8TsbEFFPIupt/AZwLoB5iOLd454NlpwxvhBZr4hlAG6KOz8K3ucpAOsAtCHbC/8HgAEAJgBYkvu3v+n6m3LvvgimWW4Ao3MVZRmAvyG38jeJ/wH4F2SHi3MAzMr9d2Ha3zuX36MBzMy9+zwAt+SOp/7dc3k+HQUvl9S/M7IeebNz/803ZFYU766X/ms0Gk1KKAeTi0aj0Wgk0AJdo9FoUoIW6BqNRpMStEDXaDSalKAFukaj0aQELdA1Go0mJWiBrtFoNCnh/wFe1dqO+xHeOQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Samples\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+iUlEQVR4nO2df3Bc1ZXnv6dbLdNSCLISG0KDkNEYKByDlaiwvdrZwQyEBG+CYiYxxN7J1KTMzBSzGQPjGTmwwST2WDNOgK3a2dmCnamaKRtiOyaKiZhxTBlqdz3IrBzJcZzYCw62cJsBJ7YcsNpWq3X3j+7bev363vt+tvrX+VS5rH79473X7/X3nnvu+UFCCDAMwzC1SaTcB8AwDMOUDhZ5hmGYGoZFnmEYpoZhkWcYhqlhWOQZhmFqmIZyH4CVj3/846K9vb3ch8EwDFNVHDx48FdCiDmq5ypK5Nvb2zE0NFTuw2AYhqkqiOik7jl21zAMw9QwLPIMwzA1DIs8wzBMDcMizzAMU8OwyDMMw9QwFRVdU0v0DyexZc8xnB5L4eqWONbdfSN6OhPlPiyGYeoMFvkS0D+cxPoXDyOVzgAAkmMprH/xMACw0DMMM6Owu6YEbNlzLC/wklQ6gy17jpXpiBiGqVdY5EvA6bGUp+0MwzClIpDIE9EWIjpKRD8loh8QUYvlufVE9BYRHSOiuwMfaRVxdUvc03aGYZhSEdSS3wvgk0KIWwD8PwDrAYCIbgZwP4AFAD4L4L8TUTTgvqqGdXffiHis8HTjsSjW3X1jmY6IYZh6JZDICyF+LISYzD0cBHBN7u97AXxPCHFJCPE2gLcA3BZkX9VET2cCm1csRKIlDgKQaIlj84qFvOjKMMyME2Z0zR8C2J77O4Gs6EtO5bYVQUQPAngQANra2kI8nPLS05lgUWcYpuw4ijwRvQLgKsVTjwkhfph7zWMAJgFsk29TvF7ZMVwI8SyAZwGgq6uLu4ozDMOEiKPICyHuND1PRF8F8B8B/K4QQor0KQDXWl52DYDTfg+SKS2cuMUwtUsgdw0RfRbAXwL4HSHEuOWp3QCeJ6KnAFwNYD6AN4LsK0ycRK0UolcOIXWzT07cYpjaJqhP/r8BmAVgLxEBwKAQ4o+FEEeIaAeAnyPrxnlICJExfM6M4SRqpRA9N58Z9iDg9jxMiVss8gxT/QQSeSHEbxme2wRgU5DP98qq517H/uNn84+7O1qxbc3Sgtc4iVopRE/3mY/uOJR/HPbA4vY8Spm4xW4ghik/NZPxahd4ANh//CxWPfd6wTYnUSuF6OnemxEC6188jCdfOqIdBPqHk6Hu0769VIlbciaRHEtBYHrg8ns+DMP4o2ZE3i7wuu1OolYK0TO9N5XO4Nx4WvmcHAT8CKPb8yhV4pab+j39w0l09+3DvN4BdPft4wGAYUpAzYi8W5xErRSip/pMt/gtbOb2PEqVuOU0k2BLn2FmhrorNSzFS+crdno+yD4f3XEIGVGcCtASj+HS5FSR5Svx4yrych6lSNy6uiWOpOK45UyCF3wZZmaoGZHv7mhVumy6O1qLtjmJWilET36edYEVyFrXG76wAIB+EPDrKipn1m1To3qSKLdzpU6GmRlqRuS3rVnqKrqmnLixrlWDwLq7b8Tj/YfxwoF3kBECUSI8sPhabOxZWPD5XqNZShn98ub7F4zbdZZ+U2MUHetfNp4nwzDuqRmRB1BRgq7DZF3rBoGhk2exdXA0/7qMEPnHUgC9xveXOwlq3d03Fg1o0QjhwsT0Y9V5MgzjjZoS+VpANQhY4+mtvHDgnbz4efVx+/WJh2X9qwY0navGep5hwTH8TL3AIl8m3LhfJCo/vX27Vx+3H5+4yvp/ePsIhk6eLTp2N2sk9gGtvXdAuV/d+ful3LMYhplJ6i6EslR4ifl+vP8wtg6O5sVLuiUe7z+sfH2UVEU9C7d7je/3kw+gsv4FgG2Do0Xnu23N0qJFb6c1EjfnGQbcg5epJ1jkQ8BrzPcLB97xtP2Bxdc6bvca3+8nH0Bn5QtAKZDb1izFib7l+X9OayZuzjMMOLKHqSfYXaPAq7/Wq3/b5H7p7ttXtF/pCjG5d7zG91tfnxxLIUpUYM2q3qeLiAHCEUh5PtZFZpUNH9Sf7hTDzzC1BFvyNrxa5f3DSc/CZ3I/WPf7yPaR/H439izE8c334ETfchzffI/Sf9/TmcD+3jvw9MpFAICHt48YXUc9nYm8RS8HHtP5rrv7RqXoAqUTSAEUuLLCyJTlHrxMPVETIh9mDRQv/lopODp0wufW/TAFYN3OEU/n5lUEvZxvT2cCq5a0FQl9mALp5MoKw5/OPXiZeqLq3TWqSIm120ewdvuIr2QoL/5aleBITMK3sWch3j7zobaompX0FPIzBVMUiHRhqGYVJteRV//0xp6F6LqutWThh06RRGH507kHL1MvVL3Im4RWlhr2IvRe/LUmYXGyDE/82p8PWyXY9oFORXIshXm9AwWi3D+cRITIdSkFuy9cuoW27DmGh7ePoKUpBiGA86m0b/GPao5HurjYn84w3qh6kXey4NxYy1ZUmZg6q1wnOImWuKO4BVmolIItkBW/y2IRo8BLrO6boZNnsetgUimoqvNVzZjW7TwEEJDOZD/DWjLZb+z5A4uvLVh4tW4H1NcHAMYnJtE/nGTrnGFsVL1PPmwLzuqvlcgGHvY49nV331j0BUZy250IetxSmjNCFJQCcEMqncELB95RDgxRIuUsRDVjSk+JvMDr9uM19nxjz0KsXtKWt9yjRFi9pK1gofmyWPFte248XZZSxY/3H0bH+pfR3juAjvUva3MdGKZcVL3IB6nVrqOnM4FlN80p2KZKWBo6eRZTtvdO5bY7UYrj9oLO9z0lhCffvRN+3qeLJJKzCV2TlZlOaPKa1MYw5aDqRV5a3grjDoC61LAb3CQseU1qsmKP8GiJxxAJN7HTF2F3xgpzpmVaf5HowlmB7CDR+a0fo713AO29A1j05I8DWf5Brj/DzBRV75MHpiMl3JYadpNM46ZejOk1bvZhjfDo7tuHsVSxhdrcGEVLU6NRvMLCKUPW7guPRajAJ2/6vP7hJJ586UjeCm+Jx7DhCws8+dDdzAoot6+hk2cLkseWXD8bb5w4V3CsY6l0dl0ByC9Eu40a6h9Wr2cA4dfaYZgg1ITIS9xE0bgtTqWL8pCf0dOZ0L6GAM9lf3UiPj6RwRc/NUe5GBkmCQ8ZslYRtG7TRdf0Dyex7vuHigR27faRgs+2oirgZsq4lQgg/7mSjBDaBfj0lMi7eNxeM6f8iLBr7TBMEEhUkNXR1dUlhoaGSrqP7r592oiY/b135B9Lf6sK+Vrda5obo8rFUPs+ACgF0P4ewOyGCIMTfcuLtlkt26bGKMYnMvmIHi/NPHTfOQDEYxH84tufK9im+15jESAjgKmQb1kC0NIUU/r6VdfMdD4A8gvFXM6YmSmI6KAQokv1XCBLnoi+DeBeZNcb3wfwB0KI07nn1gP4GoAMgK8LIfYE2VdYuE2m2dizUCvyybFUvsZMPBbBpckpTIlp8dumeZ8uocoUofLu+VToomZHZXnaZzx+m3mYZikAkEpPFbnZdKTtq9whoRN4QH3NTG6j7o5WvHr0DNp7B0CYjoLicsZMuQjqrtkihPgvAEBEXwfwTQB/TEQ3A7gfwAIAVwN4hYhuEEJ4i/UrAV6SaRKa1xKmLetUegrxWLQg7PDVo2dCSagCwrda1fsQaO8dQEs8BiJgbDytTZKysnVwFD869K7Wt+7k1pB4zWUIk1iEYDpN1TXT3UMt8Rh+Mno+PzDaP5YblTPlIFB0jRDiN5aHzZi+r+8F8D0hxCUhxNsA3gJwW5B9hYWX4lSq11qtM4k9dM/LPvxGn/gNv2xujBbVnpHnM5ZK49x4GgLuFw/l4mX/cLKohtCTLx1xlaQ1U3R3tGJ2Uyz/uCUew5Yv3apc8Ja4vS/isSiI4Hi+XM6YmWkCL7wS0SYAvw/gPIBluc0JAIOWl53KbVO9/0EADwJAW1tb0MNxxEtJXtVr3VSctJfxBbI//rXbR7BzaLRggXjd3TcaffIqZFlgP4xPZIzuCT+kp0TRYudMRAO5JULAVxa3KV1L/cNJ5cANZAcB+31hXRC2vo4Irr5TLr/AzDSOC69E9AqAqxRPPSaE+KHldesBXCaEeIKI/hbA60KIrbnn/h7Ay0KIXaZ9zcTCa1DcLtwC0Pqa7aGd9vDCMIjHIkgpnNizQxb4cjGrIbsWYhfoWJTQ3Njgun6O7noSgKdXLip4r25BOBohZFz41exuPYYJi0ALr0KIO13u53kAAwCeQNZyt9bTvQbAaZefU9F4qW2j8zXbt1vj5U1RPV64NJldK7Afp99gqsYoYcLDbKPUXJqcwom+5UUD6W3tsz0VpDN1u7KLsS7JySTwchByClFlmFIRyCdPRPMtD78A4Gju790A7ieiWUQ0D8B8AG8E2Ve5kf7mh7ePYFZDBLObYiWpRe42LNGJKQFlDR6T/9lEpkSRLUFo7x0oGjD3Hz/rqayAzn2SUGz3muSUyFXqPNG3HPt772CBZ8pCUJ98HxHdiGwI5UkAfwwAQogjRLQDwM8BTAJ4qBIia/xiDyccS6URj0ULpvOq5B0T9rK/VkyJWF7o6Uxg59BoKP7xUmZxdne0hhphs3VwFK8ePeOYsSrXTOwuHzkzs8e563z3Oto/5lyNNAw4Hp8xUXfJUG6w/2guXJpUWsBOSVFXXt6I9z6YMO5L5acNy2WzeklbyTJlvQqeiWdWLipauA0DnQ9cVX/f6lZZdtMc/OjQu75nPVZUSWZhojoX9v3XHyaffNUXKAsbVfs83Y9d+nN1vtpffZh2LJCmqpy4sWehtpeqF3RJWWEQpmlgEnhCdrD0g/27lS63tdtHiqKTrH7zXQeToQi8V1Y993q+eFp77wBWPfe643vCaIfI1DYs8jbcVDqUSH+uqVDVtjVLcaJvOZ7JdVFSoVr8W7UkWDhpBMGEWMbTl7sOy7/raMWvPvQvuPK7tQ7eptd6uf5hoorEkp3NTITVDpGpXVjkbbj9cVgjanRCaN1usqzsi3/9w0m8evSMq+NQEVTgVy9pw5FvfRZv9y3HkutnB/ik4Pxk9Hyg9QD53boR76tb4qGKo5cy124jseyEXRqaqT1qqgqlW0wLVbqEp9lNMTQ1Nijf49SyDjAPHrpyvH4RCOYz3zo4iq2Do5g/txlvvX8h0LEEJYhVHY9FseymOY4FxYDpUhVhLXr7aSLvBy8hvUx9Unci3z+cxLqdh5DOxTbne5UiG42i+9E88fni+iw6UVZVaTTVO5HleJ2acbulMZco5AbTYPBmmQU+KKl0xvXCs7WdoleiRPjul28ty0Knlwxupj6pO5HfsPtIXuAl6SmBDbuPFCQlOf1odCWCYxHCli8V/+B1DajHUml09+3D+MRkaL7gMASecYcpgseL8OrCSN24fKz3LcPYqbsQyvbeAe1zXsLdTC4AVYkDALjrqdcqyjp2m47PTEOY9t2bDADVgO7UDcttZzOGsVOyevL1jMnHbn3OatFVmpyywHtHAEqBd1MTfyyVVtaUt94jXP6ACZu6E3ldga7ZTTHXU+y7nnrNKNgysiFMPztTOcj8CbmWs3No1HXGrr2mvNt2lHY4y5VxS92FUD7x+QWIRQtDHmNRwvJbPlGUBLX+xcPoH04WvNbJ5RKLUEHvUxb42kWu5XgtyWCd6flJZlIl7KnuVYYB6lDkezoT2PJ7tyKRq0WSaIljy+/dilePnlH+2J586UjBNpPAyyYU0qLyEnPNrZ+rEz+ZsS2WxiV+kpk4y5XxQt25awB1NMLDmtT6c+Np9A8nXU2Fz6fS+R9aT2fC2GTESiwCrLytDbsOJtnyryBKFX1kdRd6aUcpcdO4Rge7eeqPurPkdZh+VG4tJPvUWdUmTkV6Ctg5dAqfarvC7eEyM0Apl6Xveuo1AN5aRQLTnaxUOGW5spunPmGRz2HKELRaSPPnNjt+lnVxTdZ0l64hHZcmp/CvZWxozcws0u2nukdMFSS37DmmHHwI5ntYvlfl5lm7fQTtvQPoWP+yp1r8THVQl+4aFT2dCWzYfUTpY7VaSHsfud1VvLscGOyuIVOcPgc0hk9YZQqckLkRC775L7gw4c3lpnIf6twqXjpZ2XFy52SEyGcIh9W8hik/LPIWNnxhgas6IHsfuT3/ty4pqqUphu6+fez7LDMzIfDyHnm8/7BngVdhCqs0lcdwut/cNnB/4cA7LPI1BLtrLHidOgNqnyqQXVyz+j7X7Tzkqj44U33Masj+jHR9BVSY3H5PvnREGz2ju9/GUmlHX7vb8W4mBkZm5qi7sgal4PH+w9g2OMruFsY1uhIa/cNJxyYqLU0xXExnkEqbaxTZy2vM6x1wdY9GiXB88z0uXslUClzWoMS8evQMCzzjGlNOxF/u+qnxvQLZWaKbvAq7D95tSK+uPzGHX1Yn7K4JAe7Cw3hBrtfM6x1Ad9++vFtl1XOvu64g6saosIdUqlw90QjlB4woEVYvaVP64zn8snphS16BV4vFrYXEMADw4cXJ/AKodVHVa3kEE6qAgSC1501ZtmzNVzYs8jb8FIzS1YpnGCsE4LJYpMiXLmPV/XyeUDy2VrJcvGkv3vtgIv+aKy9vxIHH7vK8L50Rw8ZN5cPuGht+6oKoonJWL2nLPy53M2ym/CRa4vituc2Oi6U6YhEqerzKco8lWuJ4euUinOhbjv29dygFHgDe+2ACizft9bx/N32MmcokFEueiP4cwBYAc4QQv8ptWw/gawAyAL4uhNgTxr6C4tSYQedfT46l0N23Tzu9NXXn4ZLD9UlDhPDWX2WjVNzUm9cxf24zTvx6vHAjAV3XtRrj2e0C77TdhC6sksMtK5/AIk9E1wK4C8CoZdvNAO4HsADA1QBeIaIbhBBlVTnVD23/8bNY9dzr+FJXmzZlXOLkunm8/zBeOPAOMkIU9Hn921ffZIGvQyanROBuYN0drTjx61RRm8l0RuDJl45g6ORZ5T0XNgnNupOpVAdTGYThrnkawF+g0D14L4DvCSEuCSHeBvAWgNtC2FcgdJbU/uNn85EDTuhcN4/3H8bWwdG8ZSNTxG954l+UP3L2k9UHYbR71M0uz42nlfdcKerPeC2kxlQOgbSGiL4AICmEOGR7KgHAmv53KrdN9RkPEtEQEQ2dOXMmyOEEwoulrfrRPX9gVPFK4DeX1J87Ba4hzziz//hZNDZ4+5nKzNsrL29UPq/bbsJPNjhTGTi6a4joFQBXKZ56DMA3AHxG9TbFNqUnRAjxLIBngWzGq9PxVALW+OP+4SQ27D4CP+1SOfSSaYgQJh1uHrex8xJp2R947K7QomsA87oTU7k4irwQ4k7VdiJaCGAegEOUXWG/BsBPiOg2ZC13a9rcNQBOBz7agHR3tCpdNrMaIq5/SNYpatAF1XV334h13z9U5G9l6oepEjRTt0a8+BV0pnbw7a4RQhwWQswVQrQLIdqRFfZPCSH+DcBuAPcT0SwimgdgPoA3QjniAGxbsxTdHa0F27o7WvHX993iqrkHgIIpqpseriaXzJY9x9AQYadNPeMvoNKMriwBU5+UJBlKCHGEiHYA+DmASQAPlTuyRmINl7QjMwFNtpV1uuqmnIFTtA7DBCVCwJRASaNrJFy/pvoITeRz1rz18SYAm8L6/FJj9TfqYprtswD2qTOVwCeuKKw2qSIMcfaTDc6UH47kU6Bz69hnAW57uDJMKUmOpdDeO4DOb/1YWTAsrOJifrLBmfLDtWs0mNw6Emm9rH/xp77T1RkmLM6Np7Hu+9loZqtlHVZxMZ17kquwVjZsyQekpzOBX3z7c64afDOMjngsArkGHyXyfT+lMwKP7jhUYKWHJc720sVO25nKgEU+JMYn2JJn/DMxKfCVxW040bccxzffg1+eGXd+k4aMEAXumLDEmbNeq5O6ctfoasuEQSVPWWc3xfDxjzSGkmLPlAZZkmDr4CgSLfHAhb+s7hhVKWw/4hykHj1TPupG5GVtGYn8UQFQCr3XAaGSI23Oj6fzTSqYyies+0gaHmGKM2e9Vh9108i7Y/3LRuvIKuT2AUGia40GlKaccDwW5eqVjG+iRPjul29lUa4DTI2860bk23sHXL1u9ZK2vAVvx6mLvYxFDssSe2blIvz5zkOOtU0YRiZEmVCFAXuBE6EqFxZ5OFvykiiR8XUn+pY7fkZYVr2b4lUM4wWr0HsRbdU9HY9FuRJlhWAS+bqJrnFbz0P64FW4bXXW05nARy8LliRFAAs8Ezoyk9trghQnQlUvdSPyG3sWYvWSNkehlr55FV4KP/lpsWaF5Z0pJV5FmxOhqpe6ia4BskIvF051i6vWKJqZaKtmpyUew1iKI2GY0uJVtHXRY34Todi/P3PUlchbcRJy64AwkzTPamCRZ0qGrMnkVbTDirUHuNDZTFO3Ig+UT8hNi7uVGmvP1BZeRTvMWPuwaukw7qhrkVcR1jRS14Wqu6MV8+Z8ROkqYphSI+9JVXG9S5MZDJ08q73fw0qEYv/+zFI3C69uCKskK2AuV7yxZ2HRcwwzk2Tv9cLqqVMC2Do4isf7D5d031zobGZhS95C2NNIXeJJ/3ASPxk97+sYmfomW62ScGHCfw7GtE9cXVTvhQPvGN2YQWe7Yfr3JfZGP0ETv2oJFnkLMzWNdNMblmFUpNJTiMeixsxsE90drY73n/xMVf2mrutaAy+ahl3oTNXJbf/xs1j13Oss9GCRLyCsMDGn4mbse2SCkEpn8OrRM5jyIfDb1izFPIcSHxHSF/TbdfBU0QzAz2w3zEJnqrUv1fZ6tfbZJ28hjHrZ8schrSH547D6OXWDRjwWySdrRYnQ3dHK7QUZJdICdkM8FsUzKxflBc3pfVMC2KYJDNC5eCrdcDFZ+7UOi3wO6WdMpTN5oU20xD3X5njhwDuO23WDyeYVt+D45nvyjSO2rVmKzSsWui6nwNQPV8RjyvtIdafYM1nd9Cb2mnFd6Yumbq39WoTdNShOzsgIkbfgvU4pdT5S6/ahk2dxadK66BTRDiZy27rvH0I6w8UOmCwfXJrEht1H8kZJRggkDD0NrJa21SfuNS8jQsCshqh20bQcmaymcGWGLXkA4RZfcipuJt051tpjqfQUhk7qLYqezgS2/N6tmN0U83w8TG2SmRL5zGirURLT/KIbbNt7OhPY33sHEh4t8K8sbsPmFQuRaImDUDjbDTME2QumcGUmYKlhItoAYA2AM7lN3xBCvJx7bj2ArwHIAPi6EGKP0+eVstSwiXm9A8rpKQF420VpYStODUdMJY/d1sjp7tvHmbFMESZLHlCXye4fTuLh7SPK+z8ei2BiUjjWb7IGGqhoicfQPKuhrHVqVD55oHYGA1Op4TDcNU8LIb5j2+HNAO4HsADA1QBeIaIbhBAVGTcYZvElp5o4ppA3p5aEkmU3zcG2wdFQKlXGY1Hc9+kEfnToXa6ZU+X4Wfzs6Uxg6OTZIsMkGiFsXnGLoxjrjBorY6l0/t4qV52abWuW1m10Tal88vcC+J4Q4hKAt4noLQC3AajIpeywkzNMNXGcmpIA5mSU/uEkdh1MagXeTYcgKzIcjwujVT9h9hnOTAk89gNnMdYFGpgoVZ0ap/WAehB0FWGI/J8S0e8DGALwqBDiHIAEgEHLa07lthVBRA8CeBAA2traQjgc78xkF/oHFl/raPmYBgG52KbDT5+RSg9/Y5yJRghNjeYltv7hZP6evuup1/Dm+xeMr78wkXG0ur0mY0nc3nNuF3K5sqUeR5EnolcAXKV46jEAfwfg28hGXH0bwHcB/CHUkVzKu0EI8SyAZ4GsT97VUZeAsLvQ625OuzvH62eWwtrmmJ3qJzMlHEVbWs9uBF7iZHWbZqaJljjGJyZxbrz4nnXjCvUi3FzZUo9jdI0Q4k4hxCcV/34ohHhPCJERQkwBeA5ZlwyQtdytbZSuAXA6/MOvTJyiDDb2LMTxzfcoR0JAPUIC4FZrTCCk9exW4O3vU7Hk+tnK7auXtGF/7x144vMLfCcYeol648qWegK5a4joE0KId3MPvwjgZ7m/dwN4noieQnbhdT6AN4Lsq5pwa1XoLGi53T4b4IgaJgh+E5aubokrZ6YAtIX2Xj16psA95NUV2j+cdBXzbz3GMDtX1RJBffJ/Q0SLkNWlEwD+CACEEEeIaAeAnwOYBPBQpUbWlAK3VoVuqhslUk5VCexaYfwhrWc/MetNjRGl22RWQ0S7PqRyrUihl5a4Tujlva9DJdylqGxZKwQSeSHEfzI8twnApiCfX624tSp0i7APLL5WORtggWf8snnFQmWopBtU7p1UOuNYSdXqWvGyKGqqkqkT7pkMnqg2uKxBCXBrVahi6psbI9w1igkdvwIfFGm5e1kUNfnR7/u0PkAi7OCJWoHLGpSAns6ENvXbjlyEPdG3HNfPacJvLtWNV4uZQfzEs7vBqXTe1S1xrWjrfO4mP/qug8mSl0moNVjkS4SsDfJ233Ls773DlYXhFPXAtSgZP3R3tPqOZ3dCAPn6N/b7Mx6Lov1jcaObUdVq0FQl029NqXqGRb5KSLTEsWpJm6/68rKe+DMrF4V/YEzFs23N0pKVq07kfN+JlqyYy/1EkBVkp1K+qhlGT2cCn2q7QvseDov0Bvvkq4T9vXcAALqua/VcIla6im567OVSHR5TwbT3DiAeiyCVDt+aX3bTnKIy3QCgbi1SjG6GMfjLc9r3cFikN1jkK4j5c5uVLpv5c5vzf1sXl+xhlioSLXEMnTyLtdtHQj9epnqwd3QKIxy3JR7Dq0fPBOpXrJthmNxLHBbpDXbXVBB7H7m9QNCBrMDvfeR25eudGjLHIoT2j8U5WocpIB6L4umVi9AS99+fIB6LYsMXFgR2nTyw+Frldp34R4hr0XiFLfkKQyfoKkw/sJZ4DBu+sACP7jjk6rPisWggi4wpP26t81Q6g0d3HEJGCN8WfSqdwdrtI4gQ4GdN16l3gi6H5CuLs0UMrTXs3fZhqFdY5KsYXdJVoiWe9+G7ddNsXrGQXTpVjhetle4Q63ucmo6o8Fr1VDbPccLUl8Few95tH4Z6hUW+inGTdOWmfj3Axc/qHYK8n35a5L8PA5PbUYeuL8PzB9TuR6c+DPWaDcsiX8W4SeV2U7+eoE9MYeoDgex9tHnFLXhk+4jr6BgT1hmlW5y6N/UPJ7WzB50xU++15lnkqxAvVolT/XouesZITo+l8vdRGK47r4uyqj6s+4+fxarnXs8LvWnGqVusrfda8yzyVYYbq0Q1CEix5/LFjA5rSeEwiOSqqboVUl3ilHW7aeDQRerUe615Fvkqw8kqcRoE7EWcuvv2uRJ62TuWLf/apf1jcTy8fSS065sRInS3iMkw6bqu1dN76iWpikW+ynCySnSDwNrtI/kpuPRz9g8nceHSpKv9Sj8oC3zt4lSCwA+pdAaP/eBwPmTTb7jjvN4BXN0Sx7Kb5mDb4KjyPly7fUQ5mHyQmijaVk+15jkZqsq4QpPAIq0SN1PQ/cfP4q6nXsP6Fw+XpGcsw1i5MJHJrwfJcEdVYbLuDrUlDiDfRnPXwaTR0Fi8aW/RY1VlVzE1hS17jmFe7wC6+/bVdGVLFvkqon84iQsTxZZ3LEJ5q8TtFPTN9y9w8lMdUKmVS1WFybatWWoUeiA7MzAVW3vvgwnjY8nFjND2YK412F1TRWzZcwzpTLEd85HLspdR+tfZb85IKvU+0IU7WsMl5/UOKI+/FGWTaznahkW+itC5Ys6Np7Fu5yGkp6azGP0KvZ+sR6Z2kQvufpjdFMPYeFp5H0ZcTDFMGd2luEd1v69qT6Rid00VYXLFpG2/RNnMQTf9vfLyRmWTB1kbnGEA4KOXxfIdzpob3fUySLTE8czKRRj+5mewakmb+kUCju4RVfMQeY9eeXmj8j327brXqVD9vmS0WjW7dljkqwhTxxwVp8dSSj/n/LnN+M3FTIGFRZjun+l1P0ztcj6Vznc4O/Ktz2K1TrSRFfcTtk5oG3sWIh4rlpkpOJfSMLXRPPDYXUUCHiHgrgVXobtvX35Bdf09Nxe97qOzotrBw44pZLlaYHdNFeE1G1FaJlY/J5D13dtvXAHg1aNnCvYjp6iV6tetd8hnBUgv2K3bjT0L0XVdq6tG9RJdLRy7y0XnFtG5Rtbfc3PBcUwJFJTwkFa3qr+yWxdMLSRSBRZ5IvrPAP4UwCSAASHEX+S2rwfwNQAZAF8XQuwJui8mK8BuOkOZfnRublz54+ofTgZKcV+9pA2vHj3Dfv4SUGqBB7JCeddTrxUUF3NTM8mKrkieNUrGbX0ZqzhHXBTf0y2omgYPK7WQSBVI5IloGYB7AdwihLhERHNz228GcD+ABQCuBvAKEd0ghOCYvRBQVZ+MRQnNjQ04n0o7/ui83LhBp6W7DiZx36cT2HUwySGbVcqb71/Ab33jZXznS7fm7ym3Ignoo2Gs293Ul7EPBG6jbJysblNtejeVXiudoJb8nwDoE0JcAgAhxPu57fcC+F5u+9tE9BaA2wC8HnB/DLxbUnZUNy4AXLg0WVRrxOkHcqJvubKwlCSVzmDr4CgUblmmipicypYo2Dk0in89fjbvwmtujGLTF4vdIVZ00TDWBX43s0unTmg6TFa3U236oL+1SiCoyN8A4LeJaBOAiwD+XAjxfwEkAAxaXncqt40JCS+WlOq9APDkS0dwbnw643UslS6aIptqhchWhV/qanNMiS9BiXJmhkmlM0XX+cJEBo/uzHYf092PbqxhN7NLP35wJ6tblZQlt0trPshvrRJwtK+I6BUi+pni373IDhKzASwBsA7ADiIiqBPtlHMrInqQiIaIaOjMmTMBToXxQk9nAk2NxWO8PXJg2U1zlBfzyssbsfeR2/NTaKZ+yUwJrN0+YiwPMKthWmpmN8WKFkNN4ZJA1uJ2cs4kWuJYvaRNGY2jPXYXrqRqx9GSF0LcqXuOiP4EwItCCAHgDSKaAvBxZC13a93PawCc1nz+swCeBYCurq7a+WarAKcpcv9w0rFWiN8pNFN76Mpe2634i4ppncktYnepqPDToARwtyhc7QT1lPYDuAMAiOgGAI0AfgVgN4D7iWgWEc0DMB/AGwH3xYRE/3AS3X37tOIt64CbBPy9DyaweNPeqgolY0qPrHgqC5C5jTM3hTTqXCpW/N6Huhr0uu3VSFCf/D8A+Aci+hmACQBfzVn1R4hoB4CfIxta+RBH1lQGKsvKTkYIrPv+IWWdHCvvfTDBZRCqjFjE3frI7KYYnvj8Alfhuiq2Do5i18FT2hh5qyg7hU+6cZ34DWk0NQyvFQKJvBBiAsBqzXObAGwK8vlM+Lh1rzgJvGTd3Te6GhCYysDtArgQwNDJs/i38xd97yuVntLWULKKslP4pFMz+qAhjbqG4bUCB7bVGSVxr7C+1xxjqTS2Do5qxTUei2L1kjbH8heyWJ79vVZRdlobMrlO3Cyu1jtc1qDOCLOv65WXN2bLH/stU8hUJQmLz7zrulbHjGhZLE8XZ+4UPlkPLpVSwiJf49gXtJbdNKco+zQWIYDcu2iArMAfeOwuzOsdKMVhM2UkGiFkDAO3VaSHTrpvGfj0ykVKi9tNHH2tu1RKCYkKigft6uoSQ0ND5T6MmsG0yCp9pdIq6rquNT8YtDTF8OHFSaWFHo9FC6bHbhuBM5WL1W8uM1hlT1YdV17eiLsWXOUY2mjFfu9YmYma7fbMbNnruBYgooNCiC7lcyzy1Y3pxnUrwBECll7fisFfnstPh5dcPzv/2I41JtlNtA5T2SRyM7xXj57Ji2z7x+IlaeztN549KLrSG7Ui9CaR54XXKkZ14+4/fharnsuWCHK7yDolsu+zNlu2PrZjr1a5eQVPo6uZ5FgKWwdHCxpj/GT0vGO/VRO6VKJy5VXoBqxSDGSVBvvkqxinGzfMRVYrESLM6x0omFYHKUfMlJ4GAiY9TNpT6QxO/DrlKw8iSoSrrrisJCV6TRUjGTVsydcwYXR4Ur0/I0TordBmNUS01h8THC8CL0mOpbDspjme3/fA4msda9FYWfXc62jvHcj/kzNRO7K8gXXGuXVwNJ9dy6hhka9hpCvFbW9OO1GigvZrqnoeYbVCuzQ5hUgN1QupFXYdTOKjs7zdP7JE7+YVC9ESj+W3X6aoN+3kcrRiqhjphM71FMQlVS2wyFcxbm7cns5Evjen16JLDyy+Fj2diXyPzymNj96vS0iWKpbUUuW/WiGVzuDyeGPRtdIh77H+4SQ27D6CsdR0Ketz4+mimZ8XX3mQipGqXse1sujqBPvkq5hta5a6Dguzxhmr3jNvzkccfZ1h+vhXL2lzZYExpUdG1+jCIZNjqfxsTq7DDJ08q3z9A4uvNUZc6drxuSFoxch6EHQVLPJVjp8bV/cepwUsXUcpr1BuX15irMuNkxBWK4Rsz4CNPQu1vXgJ07M1a3NsYDoLVfLCgXeMhcmAwlLWXnhg8bXagYXRw+4axjXSz5oIGCHR0pT101ZLzW65YPijQ++W+1BCRwDYNjiK/uGkcrFUVWBMWuMbexYWCWxGCKPAA9MRNqa1HJUrsuu6VqVgDfz03VAW/2sVFnnGE9JHH4QPL2Z7yfqxwAiFXYZmgvs+ncCWPccK/Mu1hADyLhTrQnuiJa6tPZccS6F/OOl5ZiMHzMf7Dxtdf6rZ5pY9x6AaPlS+fmYadtdUOTORDh426Vy7OC+WvD1TciYzbWvNRaNCulDsHZpMrPv+IU/7kDXqdf58iW6maDqeIL7+WodFvopxarZQSro7WgNnC3qJpkmOpdCx/uV8nR2n2ipMIRHKZjbrkC4UL4OnqaBdhIBPXKGuPPnoDv3gYKoN77TwzzWU1LC7popx21qtFIQdqSCt+igRujtaldacTH55dGdtCrw851KsVJgEnoC8sIbVs3fp9dM+9fGJSWzYfQTzegfQ3bfPeO1MteGdkvuqZY1npmFLvopxarZQalYvaVNOu7s7WvHGiXOeShcf33xP0baO9S8rBcFUBreasbva5vUOlLwfCwGYe3kj1m4fCbU0hfX6nxufXsswWdtRIqXAW0sZmGS8Fgf+MGCRr2Kcmi2Umq7rWgvC5SIEfGVxGzb2LCxYKwBl28np0FlglfqjjQDKBcCgPPnSkYL1lZamWIFAloK5lzfivQ8mQv9cP+0gVQvxspSBxPSpXqK+qnEtyy/srqlivNQHCZtp3+203M1qiKLruuw03Zop62SO6qJsKnH6HY9F8NTKRYiU4NDOjacLKkGeL7HAAyiJwPth9ZI2ZZ6G24Q5L/e9vHet33UtR+ewyFcxqpC3Uva7fLz/MDrWv4z23gGs3T7iej3ANLPQ/bgBg/iXQmFdsHpJG37x7c+hpzMxI2GcpZgtVCq6e8A0m/N735dzLascsLumyunpTMzINNM+bdahWg/QtXdz+mHqenvKLlZO0RTy9aUIgbzokPDDhIOplIHffI1yr2XNNCzyjCvcTptVVrs99tqLD1TX27OnM6EdeOyzg7fPfBhKc4htg6Pouq4VPZ2JktXqrxbkAKorhQBk4+LHxtNoaYpBCLhOJnOz0HpZLIL+4aQvA6fca1kzTSCRJ6LtAKQjrAXAmBBiUe659QC+BiAD4OtCiD1B9sWUFzeLoCa/aClmHDpL3z4oqAq5+cGaGbru7hurrlFKLEq+FkStuElKIwCrFG64u556DW++f6HoM60VLt0utF6YyPjOCXHTOLyWCCTyQoiV8m8i+i6A87m/bwZwP4AFAK4G8AoR3SCE4EagVYpu2ixJlClCQWfp27HG9QcRfGtmqFPmZqVgrR6pc3NJ8W7vHdB+jkoIvczS9j5ye5HQz5/bjL2P3J5/7KUyqTXL1Uu0TJCZZTUSiruGiAjAlwHIIf5eAN8TQlwC8DYRvQXgNgDqli9MxWEXwisNoXYElKU5s1+cGp2rinJJrFP6jT0L8fyBUWOiUSXw9MpF2LLnGB7ePqJ1fzQ1Oi8k69ZQvMzSrIKuwmvY7OlcDR2vmd8ztZZVCYTlk/9tAO8JId7MPU4AGLQ8fyq3rQgiehDAgwDQ1tYW0uHUL2HE/6os3fc+mNCmxofpy5zJ+GXVtD0WJTQ3Nij9x7Eo4cKlyYL+tpUu8AAKzlF3uCo3ip1SiqK87l65uiVujJapVCGfyV61jsM3Eb1CRD9T/LvX8rIHALxgfZvio5T3lxDiWSFElxCia84c7/0kmWnCiv/VuTKmRHHP1zB9mTMdv2wPQZ3dFANsC4TyRrY+Zz22MkVzesJLmQJdt7ErL2909f7+4SS6+/blSxi4uXbW6+4Fee9VW7SMrldtu4fvzAuOIi+EuFMI8UnFvx8CABE1AFgBYLvlbacAWIOcrwFwOswDZ4qZifjfUsbluzl+PyJixxrv/+iOQ1h20xy83bccTY0NSNtMc4HseaqeS6UzM172ONESx+olbdlBpwRsW7NUKejvfTDh2DDb7yDtt16OvPd0M8lKjZYxrTuUwrAJw11zJ4CjQohTlm27ATxPRE8hu/A6H8AbIeyLMTATFk0pfZlOx/94/2FsGxzNTwn9VN20R29IK8rN/lU4NcgICwKy2cPInkPY5Q6uvLxRWytI8sKBd4wuBb9uEz/35+ymWP4zZzpaJqhL0WndIWxXUxgifz8KXTUQQhwhoh0Afg5gEsBDHFlTesKK/9WVES51Z3vT8fcPJwsEXmK19N388HTRMFsHR9ESjyl98fL7CyMuXhelNLsphovpKa1Fay0FvC3kiB7ToroVJ3Hya2T4yTmwHspMRsuEUd7bKVINCNcwCzzXFEL8gRDifyi2bxJCdAghbhRC/HPQ/TDOhFXLplyd7XXHv+ymOXh0xyFjl6IwfPljqXTRD0J+f+vuvhGxgA741Uva8N0v36o8xyc+vwCbVyxES7zYDWO9hlv2HAu1MiUB+NWH7mYFBBhdZTpjIkJkvBbLbvK+FnfeNhhbayXt770jkMCbXIJhuETddEQL09XEGa81RJgWTTk626uOf9lNc7DrYNJo+USJHN0EbqM3pgC0xGM4n0oXfH/9w0nfhd5V0RMbdh/Jzxoui2WHFukKM7kDwl5M9GJFRyJU1NBbHjegb/SeEcJo7b569IxyfwnDDMoqgmFGZDlZ6mG4RO1JfHbCdjWRqKByrl1dXWJoaKjch8FUELpYdokppl36sb22CrRndbo5Dh2qAmyq43FTyyfIcWg/z0WHLzm2qb7neCyCiUmRDwVccv1sDP7ynHZQViXNmRKwZjfF8OHFyYJFb+t35SXj1g2671feE07P+yGMQYqIDgohulTPsSXP+MIeS18qd47JQpI/Zl39FGnteY3eUO3TyVJT1Zjv7mhVCk2QuG6dtRyPRbB5xS2uCrdZj9kk8KuXtKHrulbjAGldeM4I4Thg2C3j/uGkcaA+N55GLEr52dUV8RiIgIe3j2DLnmO4cGmy6NgECusMecFkqfcPJzE+MVn0XFDLu9SJWVxqmPGMKllq//GzWPVc+AnNOt9klAhPr1yEjT0LHdcivLo4VPs0+UjjsYiyLPCJX6v3qxNhN8epKi/9zMpF+MW3PwcAShHSYYoLkjOQsNoBWrEvljv5EtIZgeZZDXh65SJcmpzCufHpXAVd0TNZZ8gruut8RTyG9S8WRzW1xGMlLe8dBmzJM57RWWthVHq046ZMsdNahBe/s84qMx2HrlCZSrRNlmuEqCCbdujk2SK/rfTv210D9vDSoMgZiN81gGiEjG0aT1t8+244PZYKZUbmhO46E6mTyppnNVS0wAMs8kyF43Yx2TTl1bk47MxuimH5LZ/I13mx7kt3HEMn9QObyio0Wa5SzJNjqWyzcoVIWuP6pRDrwkv9Yu3IpRsgo0QQENqyDpkpgebGKC5MqL9zOaC5RVruXvAToaK7zg97GMgrDRZ5puIxCbibRSv7DzeiiVM+N54uiKO3+49Vx/HojkPa41bNCNyKglOzcmtiUthhldYQP9MMxqkK58X0FJ5ZuUgbcRMWTbEIUumpgu8giJ9cdZ11ax2VmlVrhX3yjGd0SVGlTpay4yWN3hpHPeVBYJxioE1ipRqYwhIF636doo+8Yl0sNrWY3NizEKuX6IsKZoTAlj3HcN+nE/n3u+3b6/Z18VgUf7XiFjy9clFJ22CWs59yUNiSZzyjasIxE8lSdvxGqXjNsLRa3/aZg64yp06kVJaxKbpEh/x8Uz0Za/TR6bEUQIWZoioSikHINJPa2LNQG+8NZAegXQeTedF146KRM4WHt48Yw2PtM7dS+saruQY9izzji3IkS9nxm5ji1kcvsZYUsCfK6KbCuqzGoZNncWnS6vqI4L5PX4NdB5OeFhXl55uKXVljxfuHk1i38xDSBpX3a5k69dG1DryNDRFcmiyO64kSYUqIAvF0anAy01RrDXoWeaZq8Vurx26VmYxbe0kBuxBLuZIWvak2uKonrYwz37xioeOaAVD8+SZ3kfUYtuw5VlRFE5h25wSxTJ0yOIHswLvqudeVAh8h4LtfvrVo3/XWpq9UsMgzVUsQEbBaZbqG4M2NUWz64rRv1zRDmNXgnLGqs7rlIqr1eOzRMrqMWF2xK7u7yHTssrplEGQbRl1GqIC5T4Gu4xRQnS6SSoJFnqlawhIBtw3BTb58N2sBOivXur1/OIldB5MFAk8A7vu02lWgc5XY3UVBKpR6Sbv36gpzolpdJJUEizxT1YQlAm4agjsJWBgx0yqXkMB0ES+V4ALOA5TfWY/X0rpymy5BjJl5WOQZxiVSwB7dcUhplYcRHqkbKJJjKXR+68cFxbqk4G5e4TxA+Z31+Ilg6ulMeBL5mQ69rTdY5BnGA1LY/FjFCY3LxBq2aHIJqbpBOQmuF1fLXU+9VtDQe/7cZt8RTG4aYwDlCb2tNzgZimE8YkoQMuEmoUb1Gid0guslWcwu8ADw5vsXtMlUTrMWN40xgOxibHvvgOfidmH0+q0X2JJnGB/4WQtw4zKxvsZtwpZdcFVVQiU6y98u8JIpZAci06xFNVtwE1ZpRVYxdWPVh9GCr55gkWeYGcTN4CBf46ZBiF1wTQIv8bpAfFksglkNkaJuWYBZcK2L2Y/3H3YUfLdVTIPU469H2F3DMBWKynUjG2jo3ERuhNLrAvG58TQuTU7h6ZWL8v1Tpbtk7fYRVz1PN/YsxPHN9+BECDH5psXpjvUvG0s91CNsyTNMhVKKZCBCVgy7+/YVfNb8uc1alw1QaCm7aadYyhK8psVpVSnmeodFnmHKjCkCJsxkIGshNLsfe+8jtysXX61I4XbTvMM0W9D1lXUbSukm4cpairneYXcNw5QRLxEwbtAJ5ayGSFGNnlQ6gydfOpKPUhmfyNZ/V1WiBKaF28lKJ6hr6Uu2rVmK+XObC7bNn9vsOpTSGt2kI8x69dVOIJEnokVENEhEI0Q0RES3WZ5bT0RvEdExIro7+KEyTO1hWkT0w7Y1S4uEvrujFROKwmBA1t9uH2CW3TTHGOrp5NMXMEe59A8ncercxYJtp85d9DSwyf4AEU2Mp9t69PVAUHfN3wB4Ugjxz0R0T+7x7UR0M4D7ASwAcDWAV4joBiFEuB2BGabK8ZtsZEJlEbuJ1AGyA8yPDr1b1FjlU21X5IXbyV1isrCB8KJj+oeT2kL8buP064Gg7hoB4KO5v68AcDr3970AvieEuCSEeBvAWwBuU7yfYeoanVUcdls5L0lWY6l0UUng/cfP5qNWpLtkdlOs6L1uMn/DGti27DkG1fwkHouwP95CUJFfC2ALEb0D4DsA1ue2JwBY66qeym0rgogezLl6hs6cORPwcBimupiptnKqLN2WeLFIm7CWSu7pTGD4m5/J+/AJQEs8hstiETy8fcSYhRrWwKYbFC6m1a6pesXRXUNErwC4SvHUYwB+F8DDQohdRPRlAH8P4E6oW0sqJ1ZCiGcBPAsAXV1dvFrC1BUzWTPdHqnjJhTSimoxU36mlyzUsJqBBCmfXE84irwQ4k7dc0T0TwD+LPdwJ4D/mfv7FACrU+waTLtyGIaxUK6a6aoBZnxiUlkIDTAvZnrxs4c1sHHnKHcEXXg9DeB3ALwG4A4Ab+a27wbwPBE9hezC63wAbwTcF8MwIaOy7h/deQgZRatA02KmVz97GAMbd45yR1CRXwPgvxJRA4CLAB4EACHEESLaAeDnACYBPMSRNQxT+UiBfOwHh3FhIvuTJRQ2BVdRLtcJd45yhkQFJQ10dXWJoaGhch8GwzAeUfn3dX1pmfAhooNCiC7Vc1zWgGGYwLDrpHJhkWcYJhTYdVKZcO0ahmGYGoZFnmEYpoZhkWcYhqlhWOQZhmFqGBZ5hmGYGqai4uSJ6AyAk7mHHwfwqzIezkzC51qb8LnWJpV4rtcJIeaonqgokbdCREO64P5ag8+1NuFzrU2q7VzZXcMwDFPDsMgzDMPUMJUs8s+W+wBmED7X2oTPtTapqnOtWJ88wzAME5xKtuQZhmGYgLDIMwzD1DAVKfJEdIKIDhPRCBHVVIF5IvoHInqfiH5m2dZKRHuJ6M3c/7PLeYxhoTnXDUSUzF3bESK6p5zHGBZEdC0RvUpEvyCiI0T0Z7ntNXdtDedac9eWiC4jojeI6FDuXJ/Mba+a61qRPnkiOgGgSwhRaQkHgSGi/wDgQwD/JIT4ZG7b3wA4K4ToI6JeALOFEH9ZzuMMA825bgDwoRDiO+U8trAhok8A+IQQ4idEdDmAgwB6APwBauzaGs71y6ixa0tEBKBZCPEhEcUA/B9k+1qvQJVc14q05GsZIcT/AnDWtvleAP+Y+/sfkf3BVD2ac61JhBDvCiF+kvv7AwC/AJBADV5bw7nWHCLLh7mHsdw/gSq6rpUq8gLAj4noIBE9WO6DmQGuFEK8C2R/QADmlvl4Ss2fEtFPc+6cip3m+oWI2gF0AjiAGr+2tnMFavDaElGUiEYAvA9grxCiqq5rpYp8txDiUwA+B+Ch3LSfqQ3+DkAHgEUA3gXw3bIeTcgQ0UcA7AKwVgjxm3IfTylRnGtNXlshREYIsQjANQBuI6JPlvmQPFGRIi+EOJ37/30APwBwW3mPqOS8l/NzSn/n+2U+npIhhHgv96OZAvAcauja5ny2uwBsE0K8mNtck9dWda61fG0BQAgxBuA1AJ9FFV3XihN5ImrOLeaAiJoBfAbAz8zvqnp2A/hq7u+vAvhhGY+lpMgfRo4vokaubW6B7u8B/EII8ZTlqZq7trpzrcVrS0RziKgl93ccwJ0AjqKKrmvFRdcQ0fXIWu9AttH480KITWU8pFAhohcA3I5sudL3ADwBoB/ADgBtAEYBfEkIUfULlppzvR3Z6bwAcALAH0nfZjVDRP8ewP8GcBjAVG7zN5D1VdfUtTWc6wOosWtLRLcgu7AaRdYo3iGE+BYRfQxVcl0rTuQZhmGY8Kg4dw3DMAwTHizyDMMwNQyLPMMwTA3DIs8wDFPDsMgzDMPUMCzyDMMwNQyLPMMwTA3z/wEHcykYzJlBOgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "samps = samps[100:,:]\n", "nlls = nlls[100:]\n", "\n", "print('NLL values')\n", "plt.plot(nlls)\n", "plt.show()\n", "\n", "print('Samples')\n", "plt.scatter(samps[:,0], samps[:,1])\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Make a 95% confidence interval for the slope of the decision line\n", "Recall that for a weight $w$ the corresponding slope is $-w[0]/w[1].$ Youn should hae to write about 5 lines of code to compute the lower and upper bounds." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The confidence interval is [0.23916648934392704, 0.500554618647888]\n" ] } ], "source": [ "slopes = -samps[:,0]/samps[:,1]\n", "slopes_sorted = np.sort(slopes)\n", "\n", "n = len(slopes)\n", "lower = slopes_sorted[int(n*.05)]\n", "upper = slopes_sorted[int(n*.95)]\n", "\n", "print(f'The confidence interval is [{lower}, {upper}]')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Try tinkering with the $\\sigma$ value in the scripts above\n", "Then answer the following..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### What happens is sigma is too big?\n", "\n", "The acceptance probability is too small, and the Markov chain has trouble moving. This causes a very small number of unique iterates to be generated." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### What happens is sigma is too small?\n", "\n", "The chain almost always accepts, but the distance travelled on each update is very small, causing the Markov chain to have trouble exploring the distribution." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Why did you remove the burn-in samples? \n", "\n", "The burn-in samples are biased by the initial guess (i.e., they do not behave like i.i.d. samples) and have very very low likelihood." ] }, { "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.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }