From a7afa4e4503a2f253919b9476617f4d5ab4d377d Mon Sep 17 00:00:00 2001 From: sakshikakde Date: Sun, 8 May 2022 23:53:47 -0400 Subject: [PATCH] homework 9 --- homework9/mcmc_SakshiKakde.ipynb | 405 +++++++++++++++++++++++++++++++ homework9/utility.py | 140 +++++++++++ 2 files changed, 545 insertions(+) create mode 100644 homework9/mcmc_SakshiKakde.ipynb create mode 100644 homework9/utility.py diff --git a/homework9/mcmc_SakshiKakde.ipynb b/homework9/mcmc_SakshiKakde.ipynb new file mode 100644 index 0000000..33cf25e --- /dev/null +++ b/homework9/mcmc_SakshiKakde.ipynb @@ -0,0 +1,405 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "Name: **Sakshi Kakde** \n", + "UID: **117472448**" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "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." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 1, + "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)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "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]$." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 2, + "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()" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "The NLL of the initial guess is 111.4339285159692\n" + ] + }, + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAe4klEQVR4nO3df4xdZ53f8ffHY+zFqxLiiUmtOLaDklUVxDaIaej+AbsVvwJiMRURJBqCgYB3Q2lZoa0U5FKkbN0SViqwKtvFBFgTmyU07Qq3lEb5sVSoKmkmizdNlmZtAk6cDWRiZ0OEEcTJt3+cM/Gd4Zw79849P55z7uclHc29555755kz5z7fc57n+zxHEYGZmVmRdW0XwMzM0uUgYWZmpRwkzMyslIOEmZmVcpAwM7NS69suQJXOO++82LlzZ9vFMDPrlHvvvfeJiNhS9FqvgsTOnTtZWFhouxhmZp0i6XjZa25uMjOzUg4SZmZWykHCzMxKOUiYmVkpBwkzMyvlIGFm9Th0CHbuhHXrsp+HDrVdIluDXqXAmlkiDh2CPXvg9Ons+fHj2XOA+fn2ymVj85WEmVVv796zAWLJ6dPZeusUBwkzq97DD4+33pLlIGFm1du+fbz1liwHCTOr3r59sGnT8nWbNmXrrVMcJMysevPzsH8/7NgBUvZz/353WndQJUFC0hWSHpR0TNL1Ba9vlHRL/vrdknbm6+clHRlYnpN0Wf7at/LPXHrtJVWU1cwaMj8PP/whPPdc9tMBopMmDhKSZoDPAm8CLgWulnTpis2uBZ6MiIuBTwE3AkTEoYi4LCIuA64BfhARRwbeN7/0ekQ8PmlZzcxsPFVcSVwOHIuIhyLiF8BXgV0rttkFHMgf3wq8VpJWbHN1/l4zM0tEFUHiAuCRgecn8nWF20TEGeApYHbFNu8E/mzFui/lTU0fKwgqAEjaI2lB0sLi4uJa/4b+8ChXM6tQEh3Xkl4FnI6I+wdWz0fEy4FX58s1Re+NiP0RMRcRc1u2FN5YaXosjXI9fhwizo5ydaAwszWqIkg8Clw48Hxbvq5wG0nrgXOAkwOvX8WKq4iIeDT/+TTwFbJmLRvGo1zNrGJVBIl7gEskXSRpA1mFf3jFNoeB3fnjK4G7IiIAJK0D3sFAf4Sk9ZLOyx+/AHgLcD82nEe5mlnFJp7gLyLOSPoQcBswA3wxIh6QdAOwEBGHgS8AN0s6BpwiCyRLXgM8EhEPDazbCNyWB4gZ4A7g85OWtfe2b8+amIrWm5mtgfIT+l6Ym5uLhYWFtovRnpUzb0I2ytWDmMxsCEn3RsRc0WtJdFxbRTzK1VLgDLte8f0k+mZ+3kHB2uP7SPSOryTMrDrOsOsdBwmzaVJ3U5Az7HrHQcJsWjQx2NL3kegdB4k2uGPP2tBEU5DvI9E7DhJN89QZ1pYmmoKcYdc7HifRtJ07iwe87diRzblvVhcfe1bC4yRS4o49a4ubgmwNHCSa5o49a4ubgmwNHCSa5rM5a5NvKWpjcpBoms/mzKxDPC1HGzx1hpl1hK8kzMyslIPEtPKAPjMbgYNEXVKuhD2gz8xG5CBRh9QrYc/UaWYjcpCoQ+qVsAf0mdmIHCTqkHol7AF9ZjYiB4k6pF4Je0CfmY3IQaIOqVfCHtBnZiOqJEhIukLSg5KOSbq+4PWNkm7JX79b0s58/U5JP5N0JF/+ZOA9r5T0f/P3/JEkVVHWRnShEvb0DGuWcuKaWdUmHnEtaQb4LPB64ARwj6TDEfHXA5tdCzwZERdLugq4EXhn/tr3I+Kygo/+j8AHgLuB/w5cAXxz0vI2xqOqe2kpcW0pL2EpcQ3877Z+quJK4nLgWEQ8FBG/AL4K7FqxzS7gQP74VuC1w64MJG0FXhQR34nshhdfBt5WQVnNJpJ64ppZ1aoIEhcAjww8P5GvK9wmIs4ATwGz+WsXSfqupP8p6dUD259Y5TMBkLRH0oKkhcXFxcn+ErNVpJ64Zla1tjuuHwO2R8QrgI8AX5H0onE+ICL2R8RcRMxt2bKllkKaLUk9cc2salUEiUeBCweeb8vXFW4jaT1wDnAyIn4eEScBIuJe4PvAr+Xbb1vlM80al3rimlnVqggS9wCXSLpI0gbgKuDwim0OA7vzx1cCd0VESNqSd3wj6aXAJcBDEfEY8BNJ/zjvu3g38PUKymo2kS4krplVaeLspog4I+lDwG3ADPDFiHhA0g3AQkQcBr4A3CzpGHCKLJAAvAa4QdIzwHPA70bEqfy1DwJ/CryQLKupO5lN1mtOXLNpoix5qB/m5uZiYWGh7WKYGWT5wnv3Zr3627dnbXKOrkmSdG9EzBW91nbH9fTxSCzrulGO4dRnQraR+UqiSStHYkHW6+lGbeuKUY/hnTuzwLDSjh3ZCH9LyrArCQeJJvmLY1036jG8bl12BbGSlE0FY0lxc1MqPBLLum7UY9gDSnrDQaLJPgJ/cfpn2vqYRj2GPaCkN6Y7SDTdueYvTr9MY+fsqMewB5T0R0T0ZnnlK18ZY9mxIyL7ei9fduwY73PGcfBg9vlS9vPgwfp+V5+lsB/bOH5SkMK+t0qRjWkrrFenu+PanWvdlEqWmI8f6wl3XJdxH0E3pTJft48fmwLTHSTcR9BNqWSJ+fixKTDdQaJPnWvTlGWTyhl8n44fsxLT3SfRF6m00Tdl2v5es5q5T6LvUmmjL1H5RY7P4MczTVeZVjkHiT5IpY2+QG1DCebns2kgnnsu+9lmgPjgB2H9+ixgrV+fPW9BYSyYxrEcVq2y3NguLmOPk+iLhPP1Ey7a8yZK+7/uuuI/8LrraiptsYMHIzZtWl6ETZsinp7dkf4/wFqHx0n0XMJt9KkPJZh4161fD88++8vrZ2bgzJnKyrmasnn3nmUd60j4H2BJcJ9E3yXcRp9KIlKZibtzigLEsPU1KW1xJPF/QE3cDVMdB4m+SKmNfkBjQwnWWCtM3J0zMzPe+pqU1fn/fnb6xnK4G6ZiZe1QXVymtk8icbVP9VPWID/CL5q4zyTxPomDB2Pq5lrqQj9YahjSJ9F6xV7l4iAxpSaoFSaIL2ddd13EzEz25pmZxgPEkimLBaWk4sNBartk6RoWJCrpuJZ0BfAZYAa4KSI+seL1jcCXgVcCJ4F3RsQPJb0e+ASwAfgF8C8j4q78Pd8CtgI/yz/mDRHx+LByTG3H9bSbsHf80KGsD+Lhh7Nmm337kmmtszXwDSDHV2vHtaQZ4LPAm4BLgaslXbpis2uBJyPiYuBTwI35+ieA346IlwO7gZtXvG8+Ii7Ll6EBwqbYhL3jiXbn2Bp5Sq1qVdFxfTlwLCIeiohfAF8Fdq3YZhdwIH98K/BaSYqI70bE3+brHwBemF91mI1u0lrBqTC9knCyXydVESQuAB4ZeH4iX1e4TUScAZ4CZlds83bgLyPi5wPrviTpiKSPSVLRL5e0R9KCpIXFxcVJ/g7rqklqhUOH4L3vXZ4K8973OlB01FK8v+aa7PnNN/vqcFJJpMBKehlZE9TvDKyez5uhXp0v1xS9NyL2R8RcRMxt2bKl/sJamtbaZvThD8Mzzyxf98wz2XrrFKe+1qOKIPEocOHA8235usJtJK0HziHrwEbSNuDPgXdHxPeX3hARj+Y/nwa+QtasZanqapPNyZPjrbdkJT7PZWdVESTuAS6RdJGkDcBVwOEV2xwm65gGuBK4KyJC0ouBbwDXR8T/WtpY0npJ5+WPXwC8Bbi/grJaHXwKZwlIeJ7LTps4SOR9DB8CbgO+B3wtIh6QdIOkt+abfQGYlXQM+Ahwfb7+Q8DFwL/O+x6OSHoJsBG4TdJ9wBGyK5HPT1pWq0lbp3BVXL3MruwaW2X9GnT1IqtrUp8CprPKBlB0cfFgupa0MXqpklFw+eds2LD8czZsqGwkWlXF7KKmB/dN876eFB5xbbVqYx6EKn9njbXZP589GD9gRzyL4gfsiKs5OBVTRBRV2FL9g9E96nxthgUJTxVuk2tjqvLU5yAHOHSIn75rD7/K2f3yUzbxAfbzVc0nU8w6lI16lrK0VKekpsVThVu92hi91IUG6L17lwUIgF/lNP+Wva0Vs6n+kbLO4ghnG3WNg4SNr6imaXpuiy7MvVBSU27n4VaK2WQS2rAg6GyjbnGQsPGkku7ahbkXSmrK07PbWynmWpLQ1nrlsW9f9m8pktLFno2grLOii4s7rhvgyfpHl1i6zbhJaJMW/7rrfvl3OtsoTQzpuPaVhI2nSyOW2h6gkNjVzrjdOJMOf/njP846qRP5822NnN1k4+nKZP1tZFwlbtxd0oUEMquGs5usOl3oMAZP5FNg3AubLiSQWf0cJGw8iTWhlOpSs1iDxklC68r5gNXLQcLG14Vbufk0eGJdOR+wejlIWD/5NLgSXTgfsHo5SFh76sw+6ulpcNsJWzZ9nN1k7XD20di8y6wuw7KbHCSsHV1JpU2Id5nVxSmwlh5nH43Nu8za4CBh7XD20dj6vMvc15IuB4m+Sv1b11T2Uer7YQx9TdhKZc5IK1E2qVMXF0/wl0tsYrlSw24jVsUtxhLfD2v5E/t45zXPGdk+fPvSKdPUt66uGquKyv3gwYiZmWRrn8TjV6PauEW6LecgMW2a+NbVWctNGuSKyrZyP7R8Su6z57O8L9o3LEhU0ich6QpJD0o6Jun6gtc3Srolf/1uSTsHXvtovv5BSW8c9TNtiCZ6OOucQG/SNJ6isg3avLn1RvCUMpXa7rbpa19Lb5RFj1EXYAb4PvBSYAPwV8ClK7b5IPAn+eOrgFvyx5fm228ELso/Z2aUzyxafCWRa6Ito86rlUlPLcvKtrQfZmdbP3VN5ew5lWavPva1dAl1NjcBvwHcNvD8o8BHV2xzG/Ab+eP1wBOAVm67tN0on1m0OEgMqPtbV2ctN2nNVVa2mZnsMxJoBE+lck4lWFm7hgWJKpqbLgAeGXh+Il9XuE1EnAGeAmaHvHeUzwRA0h5JC5IWFhcXJ/gzOqqsraDumdnqbCOYdN6lsrIdOJB9RgIDDlKZWiqlZi9LU+fHSUTE/oiYi4i5LVu2tF2cZrWZYF53LTdJkFsq2+zs2XUvfOHZx4k0gqcww2oC8dISV0WQeBS4cOD5tnxd4TaS1gPnACeHvHeUz7S2776WQi03zM9+dvbxyZNnA2gqp/EJSCReVqbtTvheKmuHGnUh62N4iKzjeamT+WUrtvlnLO+4/lr++GUs77h+iKzTetXPLFqmrk8igbb1ZK3WL2HP60uncSr9PF3EkD6JSmaBlfRm4NN5Bf/FiNgn6Yb8Fx+W9CvAzcArgFPAVRHxUP7evcD7gDPA70XEN8s+c7VyTN0ssJ4WtNy6dVk9UcTza/eSvw5r56nC+8o3GChXVmMscc3RO2XnBVLWImrlPFV4X7ltvVxRY/sgp+/0jjvh6+Eg0aQ6etVS7zxuy1IAnZkpft01R+/0rRM+FQ4STfF8yM2bn8/GRrjmqETqmUO+sK5JWY92F5eks5s8tLU9DaXv9CVLqIgzhxJWwYGHZ4FNwLB01T7XLlOirkp05aFx3XXtHCptnuP46zFERQeeg0QKyr5ls7M+ReuYokqrjkp0tRnPmzxU2hqS4yuYVVR04DlIpKDsaE9gRtLWdehUsezfWFaJT1KJln3/2zhU2rqScCvtKiqK3sOChDuum1LWq3bqVPH205Ki2bEO/bKZUOpIohr1EGjiUGkrc8gTEK6igbxfB4kmFaWrTntyd9vzT42prHJ69tnqK9FRD4Gy7arMRmorc2javx6raiJ6l11idHFJurmpzLQ3unZs/qlhzR9FrWaTtKRN0ifRl8OqL39HrZzd1PMgEdGpNvnKdazReZxKq4oKbq3ZTR3brUNN89ejKcOChOdusnZ1cP6pQ4ey1rCHH86aPfbtKy5qmxPOeR4jG4fnbrJ0dXCY7KgzobTZ6eq2/GalPhp9Eg4S1r6ezj/VRkW9VFkdP57F3EGejaQeVSbopRhsHCSsWIpHa8c0nTY6WFlBVmEtBYoOXKB1VlUJeslmg5d1VnRx6WzHdWqcUlKZJjtdU+isnsZO5qoS9Nr8/+GOaxuLb/HVSSublwY18TXvYA5CJar6urSZbOCOaxuPh7lWbtzWu7W09pWN+i5bX7WOjYusTFXNiskmG5RdYnRxcXNTRVJot+iRcVvv1traN2zQXZV/S1lzUsfGRVaqima2Nlt58WC6xKTecNuBPonUd+GgcWPuWmN03bF9tcPC5xaTa+u4ri1IAJuB24Gj+c9zS7bbnW9zFNidr9sEfAP4f8ADwCcGtn8PsAgcyZf3j1KeTgSJDlTAEZF0LdyVXbhk3DPstZ6R171fVgsCXfu/jCTh70GV6gwSnwSuzx9fD9xYsM1m4KH857n543PzIPFP8m02AN8G3hRng8R/GLc8nQgSPt2aWNd2YVNXEhHVNXsUfcYowatXdWovo16xOoPEg8DW/PFW4MGCba4GPjfw/HPA1QXbfQb4QPQtSKz81pQ1Gk9Dw21Futb23VSfRN1l7VpwntgU/cHDgsSk2U3nR8Rj+eMfAecXbHMB8MjA8xP5uudJejHw28CdA6vfLuk+SbdKunDCcrajaHRMWZ5i6ykM3ZFsFkiJcWceaXOmkmEZSm3dU6JOQ7PInOWXKYseSwtwB3B/wbIL+LsV2z5Z8P7fB/7VwPOPAb8/8Hw98E3g9wbWzQIb88e/A9w1pHx7gAVgYfv27fWG23GVnYmsPBXu6SVsXdZyZt6bJpCarXaV1qd9uepx5CuJNJqbgC8CfzTkd8wAT41SnuSam8q+cUsHWh++aS0ZtbKaomblSkxRvbj63zpFB0+dQeIPWd5x/cmCbTYDPyDrrD43f7w5f+3fAP8ZWLfiPVsHHv9T4DujlCe5IDFN37hEpfYvSP1MfIrqxdH6tlL/h1WkziAxS9aPcDRvllqq/OeAmwa2ex9wLF/em6/bBgTwPVakugL/jiwt9q+AvwD+wSjlSS5ITNM3LlEpdXJ35XCYknoxuROINtUWJFJbkgsSEdPzjUvUWiqCuv5lrpTS0pWg3YRhQcJzN9WtjXsleJrv542bkVPndM1OlklLB+931QrPAts30zoV5xCj3m4U6p0A15PrWqqGzQLrINE3rokmUud0zY7flipPFT5N3KYxkToH6qXYvOGWSVuNg0TftDwcueuVTt2jilO6nXeyt8u0pDhI9E2Lcyf0odJp6mw/hWA6rTcJsjGVpT11cUkyBbZMnamxLaXdtp3i2ZVs41RSL1MaQ2Ltwve4TkxPezDbvEdvl3ZpKrkFqZTD2ueO69T09Dq/ze6QLu3SVHIL+jirq1XPQaINqdQSFRtW6dTdBt+lXZrKVOcpZltZehwk2pBKLVGxskoH6u/Q7tIuTekMfjDbat++7Mqrq5lpVpOyzoouLp3puE6l57IhTXRod22XptbJ3rX9Z9XCE/wlKLVaokZryaJZy+6Zol1aubYz06xdw4KEs5usduNm0XQpU2nJOPNDpajNzDRrn7ObrFXjtsGXZSp9+MPtD0Ar0odBhF3q07FmOUhY7cbNoinLSDp5cnlF/L73wXnntR80upR+WyalznRLi5ubLDllzVOraatJqi9NNV1vMrO1c3OTdcLSWIrjx7MKdlxtnb33pakmpckHLR0OEpaEwXZ9yM7MlwLFjh0wOzva57QxeM5NNdZnDhJtS2E60AQUtetHnM2A+sxnfrkiLtLG2btHLlufrW+7AFNtZa7nUloMTF0Ns9q0Gku7Y6nNfPNm+MlP4Jlnzm7b5tn7/PzU/ctsSkx0JSFps6TbJR3Nf55bst3ufJujknYPrP+WpAclHcmXl+TrN0q6RdIxSXdL2jlJOZPVh7SYiozSrj/YZv7EE/ClL/ns3axukzY3XQ/cGRGXAHfmz5eRtBn4OPAq4HLg4yuCyXxEXJYvj+frrgWejIiLgU8BN05YzjR1aVa6mq2lXb+pjla3CNo0mzRI7AIO5I8PAG8r2OaNwO0RcSoingRuB64Y43NvBV4rrSXfJXF9SYupQKrt+l0eKOfgZlWYNEicHxGP5Y9/BJxfsM0FwCMDz0/k65Z8KW9q+thAIHj+PRFxBngKKMxvkbRH0oKkhcXFxQn+lBY4LWaZFFMwu9oi2OXgZmlZNUhIukPS/QXLrsHt8kmixh2ZNx8RLwdenS/XjPl+ImJ/RMxFxNyWLVvGfXu7Uj19tud1tUWwq8HN0rNqdlNEvK7sNUk/lrQ1Ih6TtBV4vGCzR4HfGni+DfhW/tmP5j+flvQVsj6LL+fvuRA4IWk9cA5wcpQ/qHOcFpO07duLR3+n3iLY1eBm6Zm0uekwsJSttBv4esE2twFvkHRu3mH9BuA2SeslnQcg6QXAW4D7Cz73SuCu6NP8IdYZRS2CErz5ze2UZ1Tu7rKqTBokPgG8XtJR4HX5cyTNSboJICJOAX8A3JMvN+TrNpIFi/uAI2RXD5/PP/cLwKykY8BHKMiaMmvC/Dzs3r18mpAIOHAg7fZ9d3dZVTzBn9kqxr0fRio8YZ+NatgEfw4SZqvoyyyvZmU8C6zZBNy+b9PMQcJsFW7ft2nmIGG2Cg9nsWnmWWDNRuDhLDatfCVhZmalHCTMzKyUg4SZmZVykDAzs1IOEgnz/QDMrG3ObkqUb39tZinwlUSifD8AM0uBg0SifD8AM0uBg0SiPF+QmaXAQSJRni/IzFLgIJEozxdkZilwdlPCPF+QmbXNVxJmZlbKQcLMzEo5SJiZWSkHCTMzKzVRkJC0WdLtko7mP88t2W53vs1RSbvzdX9P0pGB5QlJn85fe4+kxYHX3j9JOc3MbG0mvZK4HrgzIi4B7syfLyNpM/Bx4FXA5cDHJZ0bEU9HxGVLC3Ac+C8Db71l4PWbJiynmZmtwaRBYhdwIH98AHhbwTZvBG6PiFMR8SRwO3DF4AaSfg14CfDtCctjZmYVmjRInB8Rj+WPfwScX7DNBcAjA89P5OsGXUV25RAD694u6T5Jt0q6sKwAkvZIWpC0sLi4uIY/wczMyqwaJCTdIen+gmXX4HZ5BR8lH7Oaq4A/G3j+X4GdEfHrZFceBwrflf3e/RExFxFzW7ZsWeOvNzOzIquOuI6I15W9JunHkrZGxGOStgKPF2z2KPBbA8+3Ad8a+Ix/CKyPiHsHfufJge1vAj65WjnNzKx6kzY3HQZ25493A18v2OY24A2Szs2zn96Qr1tyNcuvIsgDzpK3At+bsJxmZrYGk87d9Anga5KuJctOegeApDngdyPi/RFxStIfAPfk77khIk4NfMY7gDev+Nx/IemtwBngFPCeCctpZmZroOV9xd02NzcXCwsLbRfDzKxTJN0bEXNFr3nEtZmZlXKQMDOzUg4SZmZWykHCzMxKOUiYrebQIdi5E9aty34eOtR2icwa49uXmg1z6BDs2QOnT2fPjx/PnoPvLWtTwVcSZsPs3Xs2QCw5fTpbbzYFHCTMhnn44fHWm/WMg4TZMNu3j7ferGccJMyG2bcPNm1avm7Tpmy92RRwkDAbZn4e9u+HHTtAyn7u3+9Oa5sazm4yW838vIOCTS1fSZiZWSkHCTMzK+UgYWZmpRwkzMyslIOEmZmV6tWd6SQtkt1GtW7nAU808Hu6wvtjOe+P5bw/lktxf+yIiC1FL/QqSDRF0kLZrf6mkffHct4fy3l/LNe1/eHmJjMzK+UgYWZmpRwk1mZ/2wVIjPfHct4fy3l/LNep/eE+CTMzK+UrCTMzK+UgYWZmpRwkRiBps6TbJR3Nf55bsM1lkv63pAck3SfpnW2UtQmj7I98u/8h6e8k/bemy9gESVdIelDSMUnXF7y+UdIt+et3S9rZQjEbM8L+eI2kv5R0RtKVbZSxSSPsj49I+uu8vrhT0o42yrkaB4nRXA/cGRGXAHfmz1c6Dbw7Il4GXAF8WtKLmytio0bZHwB/CFzTWKkaJGkG+CzwJuBS4GpJl67Y7FrgyYi4GPgUcGOzpWzOiPvjYeA9wFeaLV3zRtwf3wXmIuLXgVuBTzZbytE4SIxmF3Agf3wAeNvKDSLibyLiaP74b4HHgcIRjD2w6v4AiIg7gacbKlPTLgeORcRDEfEL4Ktk+2XQ4H66FXitJDVYxiatuj8i4ocRcR/wXBsFbNgo++MvIuJ0/vQ7wLaGyzgSB4nRnB8Rj+WPfwScP2xjSZcDG4Dv112wloy1P3rqAuCRgecn8nWF20TEGeApYLaR0jVvlP0xTcbdH9cC36y1RGvkO9PlJN0B/P2Cl/YOPomIkFSaNyxpK3AzsDsiOnvGVNX+MLPhJL0LmAN+s+2yFHGQyEXE68pek/RjSVsj4rE8CDxest2LgG8AeyPiOzUVtRFV7I+eexS4cOD5tnxd0TYnJK0HzgFONlO8xo2yP6bJSPtD0uvITrx+MyJ+3lDZxuLmptEcBnbnj3cDX1+5gaQNwJ8DX46IWxssWxtW3R9T4B7gEkkX5f/7q8j2y6DB/XQlcFf0d/TqKPtjmqy6PyS9Avgc8NaISPdEKyK8rLKQtSPfCRwF7gA25+vngJvyx+8CngGODCyXtV32tvZH/vzbwCLwM7I22Te2XfaK98Obgb8h63vam6+7gexLD/ArwH8CjgH/B3hp22VueX/8o/w4+CnZFdUDbZe55f1xB/DjgfricNtlLlo8LYeZmZVyc5OZmZVykDAzs1IOEmZmVspBwszMSjlImJlZKQcJMzMr5SBhZmal/j92fxuplrqOgQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "### Generate many samples from the posterios distribution\n", + "Note: the NLL function above generates $-\\log(p(w)).$ \n", + "\n", + "**You will have to fill in the formula for the acceptance probability, alpha.**" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 26, + "source": [ + "iters = 5000 # number of MCMC samples to draw\n", + "sigma = 30 # 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(wp)) / np.exp(-nll(w))######## FiLL IN THIS LINE OF CODE #######\n", + " \n", + " # Should you accept this sample?\n", + " if rand()" + ] + }, + "metadata": { + "needs_background": "light" + } + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Samples\n" + ] + }, + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoI0lEQVR4nO2df5BV5Znnv48NQperSyh/0tBpgoilwcGkR6SY2TWODITeWlhnN0LhTn5Mwcwu2cTJlpkmuDFu7KIrZBxnaiqmYCZbY0n8MRundQdWg2Pc2VCA1QQC6sDYFERtXTVFGY1BDO2zf9xz2tOnz3t+vufn/X6quvre99w+5z333v6e9zzv834fUVUQQghpJueU3QFCCCH5QZEnhJAGQ5EnhJAGQ5EnhJAGQ5EnhJAGM6XsDni58MILtaenp+xuEEJIrThw4MDPVfWioG2VEvmenh4MDw+X3Q1CCKkVIvIz0zaGawghpMFQ5AkhpMFQ5AkhpMFQ5AkhpMFQ5AkhpMFUKruGtDc9/TsntZ0c7CuhJ4Q0B47kSSUIEviwdkJIPCjyhBDSYCjyhBDSYDKJvIhsFZGjInJYRP5ORGZ4tm0SkREROSYiyzP3lBBCSGKyjuR3A/i4ql4D4J8BbAIAEbkKwBoAVwNYAeA7ItKR8ViEEEISkknkVfWHqnrWeboPwGzn8SoAD6nqGVU9AWAEwHVZjkWajSmLhtk1hGTDZgrlFwA87DzuQkv0XV5x2iYhIhsAbACA7u5ui90hdYOCToh9IkVeRJ4CcGnAps2q+pjzms0AzgLYkbQDqroNwDYA6O3tZVVxQgixSKTIq+pNYdtF5HMA/g2A31FVV6RHAczxvGy200YIIaRAsmbXrADwVQD/VlV/5dn0OIA1IjJNROYCmA/g2SzHIoQQkpysMfm/BDANwG4RAYB9qvpHqvq8iDwC4AW0wjgbVXUs47FIDVi3fS/2HD81/nzpvJnYsX5JiT0ipL2RDyMs5dPb26usDFVf/ALvQqEnJF9E5ICq9gZt44pXYo0ggQ9rJ4TkD10oyQToBElIs+BInoxDJ0hCmgdFnlhj6byZidoJIflDkSfW2LF+ySRB56QrIeXCmDyxSlUE/fJNO3HWkzg2RYCRLZxbIO0HR/IEADB0sDkLkv0CDwBntdVOSLvBkbwFhg6OYuuTx/DqW6cxa0Ynbl++AKuvDfRjqyRDB0ex6dEjxu11y67xC3xUOyFNhiKfEVcgT/+6taB39K3T44JZF6Hf+uSx8f576ZrRiT39N5bQI0KILSjyGQkSyNO/HsPWJ4/VRuRffet0ovYw6n5XQ0jTaGuRt7Hwx6ZAlsWsGZ0YDejvrBmdifZTlbuaKRIcmpkihXWBkMrQthOvthb+mITQ375u+1709O8c/1m3fW+i4+RFT//OQIHvnNqB25cvSLQv013NbQ8fytLFxIxs6QsU9LMKXLl5V6F9IaRs2lbkbXH78gXonDqxfK1fIIOMu/YcP1Wa0A8dHMXSwadDL2hbbl6YePQddvdS9KrZkS19mN4xWenfG1MKPWkr2jpcYwNXCMPi0FUy7vKHVEwkFfi5/TtRteSV98aCe2RqJ6SJUOQtsPrartpMLpoyabJQRYEnhLSgyLcZeUwI5yHwy+55Bi++8e748/kXn4fdX7khhyMR0mzaNiZvyqLJY+FPlYy7kmbMlIFf4AHgxTfexbJ7nkm0n6CYfFg7IU2krUfyRa3k3LF+SayyeElSOq/cvGtCbHl6h+DowMrIvty+fEFkTD6P9yXJPv0CH9Vu4ujAytTvEyFNoa1FvkiijLvCUjr9AukXLuDDrJEoAXPnDmymNQqCQzYC4ETJlggUdNLuZArXiMg3ReSwiBwSkR+KyCynXUTkL0RkxNn+CTvdJUD2rJHV13ZZDWWcGOyD/6/SCjwLlBBil6wj+a2q+t8AQES+BODrAP4IwKcBzHd+FgO4z/lNKoLtUEbWEXsccZ9/8XmZjkFIO5JJ5FX1bc/T8/DhXfsqAPerqgLYJyIzROQyVX0ty/GIXaoSyog7emd2DSHJyRyTF5EBAL8P4BcAPuU0dwF42fOyV5y2SSIvIhsAbACA7u7urN1pC6Z3SGBoJmmoxSSudbMWJoSYkdZgO+QFIk8BuDRg02ZVfczzuk0ApqvqnSLy9wAGVfXHzrZ/APAnqjocdqze3l4dHg59SaMpIrsm7Fhxjmv62ywXhrgj+TwvPrbPiZAiEZEDqtobtC1yJK+qN8U8zg4AuwDcCWAUwBzPttlOGwkhiaiUFWpJkgVkk6IF3m2n0JO6kylcIyLzVfVF5+kqAEedx48D+KKIPITWhOsvGI+vBq7fe10o8o6CkCaSNSY/KCILAHwA4GdoZdYArRH9SgAjAH4F4PMZj0MsENecrGhODvYlFuwko29eDEg7kzW75vcM7QpgY5Z9E/vkYU5mwl9Me4q07H9N5CW6RYdiWBmLVA2ueK0pacQkrjlZWvELm0A9qy3hDxP6sP1UdeQdds51rPdLmgdFvoakLbNnKvNnKtgd5LdjCq3EIagknx9bFbuSkCZcBMTrU93q/ZLmQZGvIWmLhweZk5nK/IVVs4oT864CUa6Vy+55ZnyBVZ53Cu4dVByTOkJsQ5GvIWmLh8epYuUSt5pV0lJ67gUhKkafFq9YR7lWuvbFea+knTWj03jR9F8g6ZJJbEORryGmsEscr3ibVayC3DDj4o/Rp70byDoCT2pfnBT3Timu62dcN9E03DF0BA/ufxljqugQwdrFc3D36oXWj0OqRdsWDakzcYqHF0HWWqlujN5muKenf+f4T5kIWnMdaQqi51GD9o6hI3hg30sYc1a4j6nigX0v4Y6hI9aPRaoFR/I1JEnYJS1L580MDNkkrWZ1yfnn4vV33rfVrVDSCnvabJ60E7Zl8OD+l43tHM03G4p8Tcm7eLhJ4JNOFEYJfNkjbhNx8+jjCrrpolkUYwaPKlM7aQ4M15BJmITXL1Jl10r1CmxU2OHkYF+pfvQ71i+JfReUx/vaIcH7NLWT5sCRfBtiK8QQVHjExBSJlycfh6CqU6ZwhEtP/87Si45474JMF9K8smvWLp6DB/a9FNhOmg1F3hJz+3dOqHNahfqmQdhe5v/Jno/ECkN40yWzhmi8WURJ9pU0kyYv64OwPueVPunG3Zld035Q5C3gF3igVSJrbv/OSgq9TeII/CXnn5ton66wBhmqebOIiojnF203fMfQkdyE9+7VCynqbQhF3gKmKESVprTKmuC85PxzsX/zsgltF0zrwNtnoo3SisgiqhoP7HsJvR+d2ehzJMVCkW8Dkgp8mD+Ntz1qhGvafviuFbH7ZMoianJ+N03NiE0o8iQQV6DDxLinf6e1fPokXHPnE7HuBEwk9d4pOmRDUzNiE6ZQWsCUhFan5LQgEYsz2jalBrq+LLZXn2YV+CoQ54IRxxZ66OAolg4+jbn9O7F08GkMHWSFTTIZirwFTgz2TRL0qmbXBJF1lLpj/RKcHOxLVMnJ9NqovsQR+KTnc8G0jugXWebkYB9uvb7buD3Kh8idlB596zQUH9pNU+iJH4ZrLFEXQY9LERO1UWK87J5nJqQ9xslzd4VzyjmCsx9Mnvqecs7k+6vDd60o5Q7h7tUL0fvRmbHtn72ktZsm7UftRb4u+ellkqTQR97iHnchll/ggXh57g/sewl3r14YKPAAcPYDHe+DN/MnajJ46eDTxsyeLIvL0mYQpbWbJu2HFZEXkf8K4NsALlLVn4uIAPhztIp5/wrA51T1JzaO5aWd89PzIE0Wjo39276wxN3f6++8j8UDuyeleAZhqr5lY3FZGh+iLHbTfmhB3Gwyi7yIzAHwuwC8a6Y/DWC+87MYwH3Ob6vUIT+9Kth0TCzbZTFunn0cXAO1Ikr5ecXUT9L3NEmVr6g+ee0OXAtiABT6hmBj4vXPAHwVE7V1FYD7tcU+ADNE5DILxyIZcCdHvT9p9pFluw0O37Ui92OYSBsO8fu5+0l6N7P62i5suXkhumZ0ZvKuD7MgJs0g00heRFYBGFXVn8pEN7suAN5vyStO22sB+9gAYAMAdHebsw2IfYLi3jaIk2OflrJNxtKEQ4B8RNOG3XSTLIiDJurzLu1YByJH8iLylIg8F/CzCsDXAHw9SwdUdZuq9qpq70UXXZTob5uQn14WaQS+6DBNkKC7fbbVlyS+Olmqb1VVNJtiQWyaqI8q5t4ORI7kVfWmoHYRWQhgLgB3FD8bwE9E5DoAowC8HqaznTarnBjsa1R2TZFVhpIIfFkxeFMf3UlNb7/i3DVM75AJtshBvjpBCGDMeukQCRRwv0iaXlcGcd6rulkQm74redfwrQOpwzWqegTAxe5zETkJoNfJrnkcwBdF5CG0Jlx/oaqTQjU2qKug+7FtAZwFW8eLk7pp8plPkvYZRdaMkbDvWNxwh8nPvWii3lNm1zSPvPLkd6GVPjmCVgrl53M6Dqk5Z9XORUUQnFUlAI5vse/RHhXu8o/k/X7ufsrOWKpKH4h9rIm8qvZ4HiuAjbb2Tewz/+LzAkWq7InNtBQZuosznxEU7ijazz3o/WgqTfs+26T2K17bgTxi9bu/ckPsbIQi5wqycMIX4lEkK8gd9zyjRvBVCHeYFgo2lSTf53ZDtCKTQQDQ29urw8PDZXejFNLEn4sQ2rB+xTl+1vNKcoHJ2te4FHWcLFT1+0TyQUQOqGpv0DaO5CuCzYnGJOQ9Sk9zXt7RN4WnGPg+NxeKfIVI6+melqIyetKMvKtM3eO/FPT2giJPAokjwD39O9sq7hn2noS9D6a7pai7qCzGYWHZRqS9oMiTTLz4xrvjYlXUCLGMieC0cfikzpum9qTGYU1bKEjSQ5GvOKYRX+fUDgwdHI3lXVKUKBaxcCsqxFSXTKC0PLj/5UiRN71H6tvWpPeFmKHI14CuAO/wuLa3RXm4pyGpKF9z5xOh+7Nxkan6RSLKGiHJ59rTvxO3Xt9NL/mGQ5GvAXWqAhRUhCNMOOMKqO3yfEMHRydVY7rt4UOBr40Szrj1a21g2zgsiZc8wz/1hIW8a4DJ3jat7W2euEU4XGzdSdgW+KAi2LbI8y7JZBw2dHAUSweftnKMIFvksCpspNpwJF8DbFUBajeSiK2/KLaftOsYwkJScfbnuleGhVLci1bUOcQlKCTEKmz1hSJfA9IWe86LshZuJSGP/vnDMHHTTE3tXqvjLHMBW588Zk3ggfp5yZNwKPI1IW0VoLwE+ZLzz50UmnHb8yBOXdcqTZDGwVtIPEvf087N3Hp9d6D98fSp52Bu/87SBxPEDhR5MokpAiz+2EzsOX5q0ral82YCAPZvXhZ48YhThCMNh+9aMWny9YJpHaXWew0j7oU16EKZlFkB2VdAKytrT/+NkX/vZtcIgHPOEbz7fus99s5V1Hlx1eKB3RPe57jFYpoCDcoaThajqnXb904Q+qXzZmLH+iUAgCs375pQZclleofg6MBE//aqLV4yEZSqCqRb7JSEtEZv7t8FxeQ7p3YkLuy9dPDp0ItFHbNr/ALv0jShp0EZSYUr6EBL1PccPxUpakHC7xexddv3TtiP9+IRRt4Xi1ffOp14f0lLEKYhagGYrTmbqFTdqgt6EKY7JRt3UHWBIl9jgv75u2Z0TvhHt4Fp1J60bycH+ybdHQDAnuOnsG773lChT2qmlmYuIklKqi1BtzWHkXbOxosp7FPFVF0SH4p8RYkatZpExv0nTZv7HSSYNgQ+rB1AYPw/K+65xBXkIIHLO4vo9XfeR0//zth3M3nCVN1mQpGvILYsgOOk1dmOjU/vkMQXhSoR5oGTJ/5QWBmZQlVL1bVB0VlgVSSTyIvINwCsB/Cm0/Q1Vd3lbNsE4A8AjAH4kqo+meVYVcV2yTHbKwhNWREutk3Fjg6sTC2QRdoDJOlHWX0oS+jrLOp+9m9e1vbZNTZsDf5MVRc5P67AXwVgDYCrAawA8B0R6bBwrEoRVND5xTfexbJ7nkm1v6Cl41mZNaMzs1hM7ygnUc6W2NYtf96Pqf91P6+icNchuD/tJPBAfuGaVQAeUtUzAE6IyAiA6wDszel4pWAq6BxW6DmMPIIcNuKpRwdWppp8tUlWG+Gk8fmqQUEnabExkv+iiBwWke+JyEecti4AXpejV5w2Yom4//S2br2PDqxMLKg2qbJlchBTDDc/Jwf7xheUEVIEkSN5EXkKwKUBmzYDuA/AN9EahH4TwJ8C+EKSDojIBgAbAKC7uzvJn7YdfvEsenQa9zhFC2+cCUubfZoiwNmAm5opAoxsad1x+Ld7+7Vj/ZLAVNK8SFvVijSDSJFX1Zvi7EhEtgP4e+fpKACvJ+pspy1o/9sAbANaK17jHKsq2C7oHLZ0PGwhSpwMgqzhjjgiWQXjsqBShLZj+6b9ndX4mVHedMm4n0uazy+OFz6FvtlksjUQkctU9TXn8R8DWKyqa0TkagDfRysOPwvAPwCYr6qhOX11tDXII7smzdLxvDIIyhTtrEJdtTh8FjGNMxpPUgglS7/88zNBVhakWPK0NfiWiCxCawB6EsAfAoCqPi8ijwB4AcBZABujBL6uZBH0INIuHc8jY6AsceTIMjl+/xrbhVBcgibg3xtTXLl5F4W+omQSeVX9jyHbBgAMZNk/aT/KGH3bOqYpVl8EQZ7yNj3mXUwZVnVeANd0uOKVlE6SYhwm/5uy8U66lkEV6/16afcFSWVCkW8YUfHSMmx/3WPYOrZNgc9y7kmM0fJ+j03mYlEU8dkH2f16C6aQfKHIN4ioeKktT5y0BAlgXscuM65v89hRdwbusYLMxfLom8mbKGxVNO1+y8XGYihSEWzHS20LZdhFxmaao2k9QdBr47zOxNLBpzF0MDAz2ApRYStvf1df24UtNy9EV862wEcHVk4SdGbXVBuO5EkoQcJnKsOXVajjjuqXzgsvTRhEHvYH3gyWKph6ueZiYX23cedEQa8XFHmSGFNd1SLi0a5ImYqPZFndGWQ4F8XpX49h65PHJol8WXMfcTC9R3n1j3a/5UKRL4HLN01c9u5mZmQlTbzUNjaEwiQKXtJMvoaNYtMIvIs/syXN3IfNrJy0K4/jjvKTfn9p91sujMkXjP8fBGjlVl++Kfs/eVS8tC6WtVH//FkEcfHA7sD2MIF349+meLci27xCHmmXeX2mab+/7W73WyYcyReMabGMrUU0UfHSsrNOyrYYyJLRESeDpYjzK/MzzPv7m5W0tiBNhiLfRgR5mxQ9YRg2sVmGeCURZW95vDQ56UBrxJskNFe1u6wqE1R0R532dhb6thP5ooSuaiOKMG+TMjJDosTr3lsWBZprmdrTEEfg/Y6icTJYwkg64s2SDVOFOyebsfgoM0DTW1uRm4zSaCuRL0rowkYUYV7keWLyNgnKDLFJ2hH71iePJWqPS5KMjqyOokUTJKhJiPO5JPn+2lzpGlZqs06fURm01cRrmNDZJGxEMbKlb9I/hK3smjBM3iZ5ep6kqebkTmCawiGvvnU69cg2ySjy5GCfUTxsj45thGRMgmrikvPPnTARGrcPSb6/Nle62i612U601Ui+DKELIm9BD8LkbTIr5xWSSYgjnhrzdUG8/s776OnfiaXzZk4o2lEked2xRQl60AUgqKJWnHBmGd/fOIQV3Wln2krk6yB0eRGUGdI5tcNKoe+6sef4KazbXnxN+bzu2KIuet67l7C7q3tvWRQrnJnXOo+snBjsq9xcWBVoK5EvSuiqOKLwZoaYRmlx4+d3DB3Bg/tfxpgqOkSwdvEc3L16YaL+zNu0K9XfRXHBtI4Jlgsm9hw/VbhjZJgQluVe6SXOvE1Ynrz//OIsaotL3FKb7S7oQbSVyMcROhtUdUThZoYEEXeV5h1DR/DAvpfGn4+pjj9PIthp/y6MDhEcvmvFpAn2MKqUophHX5JMvsYJZybJk9+/eZk159PdX7nBeqnNdqGtRB4IFzqblCnoeY4IH9z/srHdL9ZxUviC/i4taxe3ascHjUjLJuz9t3FnZMIbe4/6DpjCmVnmQWxCQU9HW2XXtANpMlqSMGYo/G5qj8rc8P5d2gtRhwhuvb57XBjjTKSHOVbaIk7mintn5L4P7h3OHUP267NGfQduX74AnVM7rB+XlEvbjeRJNjpEjIKe5u86ZOJMRRIDL9Nro6oklZld4yfJnZEJW4uebnv4EO69ZdF4ODPJp5z3Og+/66h/3su1uyaTySzyIvJfAGwEMAZgp6p+1WnfBOAPnPYvqeqTWY9Fymft4jkTYvJe1m3faxRP09+5IZYokozyTRPsW25eWAnfdy9hd0aucMcRMO/7k0XwveHMuPvJO7smyFba/669fWYM19z5BIU+gEwiLyKfArAKwG+o6hkRudhpvwrAGgBXA5gF4CkRuUJVqxUoJYkF4e7VC40iH2b/645K84o9h53HtCnnVFLg4/L2mbFJ51d2NaYkRViyzhHFtZWOk1XVjmQdyf8nAIOqegYAVPUNp30VgIec9hMiMgLgOgDFJycTI0VPpt29eqH1lEkg+jzOnP0Afzv80gSRjyp4Hkac8IhJxIJWpqbBW7u3CoTlzlcpg6kdyTrxegWA3xaR/SLyf0TkN532LgDeYOMrTtskRGSDiAyLyPCbb76ZsTukLp7xReMdDYYVPI9LlCWAexFYPLB73Kqhp3+n1eLVptq9Yd+BuN+DoYOjRv98/z7yrJFAshM5kheRpwBcGrBps/P3MwFcD+A3ATwiIh9L0gFV3QZgGwD09va2u2FcIEnzg20KelBVqTQ1VquEzYLnUT48ZRH2HYgKoQStMwib08jbY970ffNzwbTgzKB2z6+PHMmr6k2q+vGAn8fQGqE/qi2eBfABgAsBjALwzqjNdtpIQsLc9/LGFMLYsX7JJEGvUsZKUVQhdzwtQXch7t3GbQ8fKsTILy5B3zf/0MM0OV3m/09VyBqTHwLwKQA/EpErAJwL4OcAHgfwfRG5B62J1/kAns14rLakLPe9qLuBKgh6EpGty11GEmzU7k3yHhZt5Ocl7feN7pXZRf57AL4nIs8BeB/AZ1VVATwvIo8AeAHAWQAbmVlTPargl5KWpKNor0iYCp67+7Wd+ROHsIpZQZjusvL8PE1GfmXVSCDxyCTyqvo+gFsN2wYADGTZP8mfKgv6NXc+MSEtLs2ClyAxPDqwMnDy1SUPX50w/AvCTER9VrZ8YoIIM/Ib2dJXWWdKwhWvlSeu+56Lf+FIXWPlfoEH4i94OTnYN56q+N5Ya1GRv2CIV/jnbdoVuCjJpq9OGHEXhGXBPcckdyldMzojjfyqaMTnJen/TxMRTbFEPS96e3t1eHi47G5UjrjZAUErA4HqCX2c88kyqWmyuDVVhgo7VtrRs2k/cczI0oZckrxnrtdPlnMPKnMJVE/o2yG7RkQOqGpv0DaO5GtA3C+kKc0s7orBIiiiVmfSsnNxfXWC8AthlEDHWRBWRAjNvUvJMi9Tl8LZTRP0pFDkSaHkne2Qxqwrq6+O//h1wIb7J6kHFHlSSeJWeLIhUHn76rjkncmU5AIXd7LXJnXN5Ko7FPkGUfeVqF4O37UicPI1DmExeRN5+eq45Jn5kgYbk71hZS6DwnJBlHX+7QSLhjSIOqxENWU1BLVHZdGYFgPt37xskqCbJl3bDQEmFFgJw+u54/54OTHYN2nlqQC43JDRQsqBI/mGUSVBD8Jmrc4wB0YK+mTuvWVRbLvluHceQVk0dbZ7SEueJRyzQpEnhZNE0Ou8KhcoXvDq+H7Vrb9+bBW3zwuKPMkVGyOcOv3DeylrRFvU+5WnONcpVm+jhGOeUORJblR9hFMF6iJkfsLCOaZVpjYYOjg6XoM2bCVukSQtbl80FHmSG2EjnMcPjmb2pakzdRX3OCSZd0lyt+P3uR996zQ2PXoEAEoV+iyL6YqAIk9yI2yEk9aXhhRD1th+HqtMtz55zOhzX6bI21xMlwcUeZIbphGOiXYtxFzFEARQvbsNk599mT73QHGL6dJCkSe5YRrhVIEiMjrijIaDQhC3PXwItz18KHa/qpy+F0WSO4ZZMzoxGiDoJp/7Iom7mK4Ml1i6UFqkHdzukhIkQGHCX8ToMYvzom2WDj4dKFx+TP3yT267xF3wlIWiUx+T1p6tGnm6xNKFsgCKcFesI0EjHP+kq4tbiLmdLpZZQw1lpu8VfUF0hbyKoa04lOUSS5G3BGtJxifIl8bNrmm3i6UpBBGXstL3ylrAtPrartqIelWgyJNSMGXRtNvF8vblCyaFIJJgM30vbpWnqpmtkXAyGZSJyMMicsj5OSkihzzbNonIiIgcE5HlmXtKSANZfW0Xtty8EF0pJw9NaXpJ0/eCqjyp007sYHKDzdslNpPIq+otqrpIVRcB+AGARwFARK4CsAbA1QBWAPiOiHRk7GulSeKuSMrFNNosaxS6+tou7Om/MVW/7l69ELde3z0+cu8QSTXpWpcqT3WmLJdYK+EaEREAnwFwo9O0CsBDqnoGwAkRGQFwHYC9No5XRWy6K7YzRRVermpYIW6/6m7qlZSolbF1WTFdhkusrZj8bwN4XVVfdJ53Adjn2f6K0zYJEdkAYAMAdHd3W+pOOVDQs2PzYlnVRUZZabeYeBzrA66YNhMp8iLyFIBLAzZtVtXHnMdrATyYpgOqug3ANqCVJ59mH6RZ2LhYVtXnxEQZi2S8hFV58lNVO+N2XTEdRaTIq+pNYdtFZAqAmwF80tM8CsA78zPbaSOkEJL6nJQpWkGLZPYcP4V12/cWJvQnBvtiZ9cA5Qs6iY+N8n83ATiqqq942h4HsEZEponIXADzATxr4ViExCKJz0lY+KMIylok42Xo4ChmzeiEAOia0Yl7b1lkFHhSL2zE5NfAF6pR1edF5BEALwA4C2CjqvJeihRGlX1OqkYRoa0kdwlpcVdMk4lkHsmr6udU9bsB7QOqOk9VF6jq/856HEKScPvyBeicOvGfvnNqB25fvqCkHtnDdgpoWGjLBllz8OOcV12ya8qAK15JI6mTz8nSeTONxlUmbMbE87bwtZGDzzmA9FDkSWOpi8/JjvVLSs2uYWir2VDkSdtThZTAPAX98k07cdYzbJ4iwMiWD88tyD+nKaEtQpEnBEBzwwF+gQeAs9pqd4U+79BWkhx8Yh+KPCENxi/wpvY8Q1tJc/CJXSjyhJDcoaCXh43FUIQQQioKRZ6QBjPFEPg2tZPmwXANITG5cvMuvDf2YWR5eofg6MDKEnsUzciWvsjsGtJsRHOuBZmE3t5eHR4eLrsbhEzCL/AudRB60nxE5ICq9gZt40iekBgECXxYe9Mo2wqZpIciTwgJxZYV8h1DR/Dg/pcxpooOEaxdPCdxmUKSHIo8ISQUG1bIdwwdwQP7Xhp/PqY6/pxCny/MriEkBtM7gtNRTO1kIg/ufzlRO7EHRZ6QGBwdWDlJ0DnpGp8xQ4KHqZ3Yg+EaQmLSroKexgrZT4dIoKB3CO+E8oYiT0gNKNP7xYYV8trFcybE5L3tJF+YJ09IxQmqrARUy+QrjlUzs2vyg3nyhNQYG5WV8iSsELpX6O9evZCiXgKZJl5FZJGI7BORQyIyLCLXOe0iIn8hIiMiclhEPmGnu4QQQpKQNbvmWwDuUtVFAL7uPAeATwOY7/xsAHBfxuMQQghJQVaRVwAXOI//JYBXncerANyvLfYBmCEil2U8FiFtiSn/hHkpJA5ZRf42AFtF5GUA3wawyWnvAuBd5fCK0zYJEdnghHqG33zzzYzdIaR5nBjsmyToVZp0JdUmcuJVRJ4CcGnAps0AfgfAH6vqD0TkMwD+GsBNSTqgqtsAbANa2TVJ/paQdqHKgl6FQujETKTIq6pRtEXkfgBfdp7+LYC/ch6PAvAmwM522gghDYSCXl2yhmteBfCvncc3AnjRefw4gN93smyuB/ALVX0t47EIIYQkJGue/HoAfy4iUwC8h1YmDQDsArASwAiAXwH4fMbjEEIISUEmkVfVHwP4ZEC7AtiYZd+EEEKyQxdKQghpMBR5QghpMBR5QghpMBR5QghpMBR5QghpMBR5QghpMBR5QghpMBR5QghpMBR5QghpMBR5QghpMBR5QghpMBR5QghpMBR5QghpMBR5QghpMBR5QghpMBR5QghpMFkrQxFCCMnAuu17sef4qfHnS+fNxI71S6ztnyN5QggpCb/AA8Ce46ewbvtea8egyBNCSEn4BT6qPQ0UeUIIaTCZRF5EfkNE9orIERH5XyJygWfbJhEZEZFjIrI8e1cJIYQkJetI/q8A9KvqQgB/B+B2ABCRqwCsAXA1gBUAviMiHRmPRQghjWLpvJmJ2tOQVeSvAPCPzuPdAH7PebwKwEOqekZVTwAYAXBdxmMRQkij2LF+ySRBt51dkzWF8nm0BH0IwH8AMMdp7wKwz/O6V5y2SYjIBgAbAKC7uztjdwghpF7YFPQgIkfyIvKUiDwX8LMKwBcA/GcROQDgfADvJ+2Aqm5T1V5V7b3ooouSnwEhhBAjkSN5Vb0p4iW/CwAicgWAPqdtFB+O6gFgttNGCCGkQLJm11zs/D4HwB0AvutsehzAGhGZJiJzAcwH8GyWYxFCCElO1onXtSLyzwCOAngVwP8AAFV9HsAjAF4A8ASAjao6lvFYhBBCEiKqWnYfxunt7dXh4eGyu0EIIbVCRA6oam/gtiqJvIi8CeBnFnZ1IYCfW9hPFWnqufG86kdTz62O5/VRVQ3MXKmUyNtCRIZNV7W609Rz43nVj6aeW9POi941hBDSYCjyhBDSYJoq8tvK7kCONPXceF71o6nn1qjzamRMnhBCSIumjuQJIYSAIk8IIY2mcSIvIiedIiaHRKS2K6tE5Hsi8oaIPOdpmykiu0XkRef3R8rsY1oM5/YNERl1PrdDIrKyzD6mQUTmiMiPROQFEXleRL7stNf6cws5ryZ8ZtNF5FkR+alzbnc57XNFZL9T+OhhETm37L6mpXExeRE5CaBXVeu2mGECIvKvAPwSwP2q+nGn7VsATqnqoIj0A/iIqv5Jmf1Mg+HcvgHgl6r67TL7lgURuQzAZar6ExE5H8ABAKsBfA41/txCzuszqP9nJgDOU9VfishUAD8G8GUAXwHwqKo+JCLfBfBTVb2vzL6mpXEj+aagqv8IwF/NdxWAv3Ee/w1a/2i1w3ButUdVX1PVnziP3wHwT2jVUaj15xZyXrVHW/zSeTrV+VEANwL4n0577T4zL00UeQXwQxE54BQkaRKXqOprzuP/B+CSMjuTA18UkcNOOKdWIQ0/ItID4FoA+9Ggz813XkADPjMR6RCRQwDeQKvC3XEAb6nqWeclxqJHdaCJIv9bqvoJAJ8GsNEJDTQObcXZmhRruw/APACLALwG4E9L7U0GRORfAPgBgNtU9W3vtjp/bgHn1YjPTFXHVHURWnUvrgNwZbk9skvjRF5VR53fb6BVXLxJtWVfd+Kjbpz0jZL7Yw1Vfd35Z/sAwHbU9HNz4ro/ALBDVR91mmv/uQWdV1M+MxdVfQvAjwAsATBDRNyiSrUuetQokReR85yJIYjIeWhVrXou/K9qxeMAPus8/iyAx0rsi1VcEXT4d6jh5+ZM4v01gH9S1Xs8m2r9uZnOqyGf2UUiMsN53AlgGVpzDj8C8O+dl9XuM/PSqOwaEfkYWqN3oFXa8PuqOlBil1IjIg8CuAEt29PXAdyJVsH0RwB0o2XJ/BlVrd0EpuHcbkDrtl8BnATwh544di0Qkd8C8H8BHAHwgdP8NbTi17X93ELOay3q/5ldg9bEagdag95HVPW/O1ryEICZAA4CuFVVz5TX0/Q0SuQJIYRMpFHhGkIIIROhyBNCSIOhyBNCSIOhyBNCSIOhyBNCSIOhyBNCSIOhyBNCSIP5//7Ay7xPu1zFAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "### Remove the burn-in samples" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 28, + "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()" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "NLL values\n" + ] + }, + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABHV0lEQVR4nO2dd3wdxbn3f3OOqmXLTa7YxjY2bjSDQomB0EPxpSSBBJJAAgm5IbzpIXATEnJv7r1cyKUm3ECAhBIgkNAJ3TR3bOOKi2Rbsi0X2bJk2eo6Z94/ts3uzrajPWVXz9cff7Rnd3Z2Znf22WeeeeYZxjkHQRAEEX0S+S4AQRAEEQ4k0AmCIGICCXSCIIiYQAKdIAgiJpBAJwiCiAlFubxYVVUVnzhxYi4vSRAEEXmWL1++j3M+witdTgX6xIkTsWzZslxekiAIIvIwxur9pCOTC0EQREwggU4QBBETSKATBEHEBBLoBEEQMYEEOkEQREwggU4QBBETSKATBEHEBBLoEWZZ3X7M27An38UgCKJAyOnEIiJcvvTHRQCAutsvynNJCIIoBEhDJwiCiAkk0AmCIGICCXSCIIiYQAKdIAgiJpBAJwiCiAkk0AmCIGICCXSCIIiYQAKdIAgiJpBAJwiCiAkk0AmCIGICCXSCIIiYQAKdIAgiJpBAJwiCiAkk0AmCIGICCXSCIIiYQAKdIAgiJpBAJ0Ln1hfX4vI/Lsx3MQii30ErFhGh88Ti+nwXgSD6JaShEwRBxAQS6ARBEDGBBDpBEERMIIFOEAQRE3wNijLG6gAcBJAC0Ms5r2aM3QngXwB0A9gM4Juc85YslZMgCILwIIiGfibn/DjOebX6+20AR3HOjwGwCcAtoZeOIAiC8E3GJhfO+Vuc817152IA48IpEkEQBJEJfgU6B/AWY2w5Y+x6yfFrAbweXrEIgiCIoPidWHQq57yBMTYSwNuMsQ2c8w8BgDH2CwC9AP4qO1H9AFwPABMmTAihyARBEIQMXxo657xB/dsI4AUAJwIAY+wbAOYC+CrnnDuc+xDnvJpzXj1ixIhQCk0QBEHY8RTojLEKxtggbRvAeQDWMsbOB3ATgIs55+3ZLSZBEAThhR+TyygALzDGtPRPcc7fYIzVAiiFYoIBgMWc83/NWkkJgiAIVzwFOud8C4BjJfunZKVEBEEQREbQTFGCIIiYQAKdIAgiJpBAJwiCiAkk0AmCIGICCXSCIIiYQAKdIAgiJpBAJwiCiAkk0AmCIGICCXSCIIiYQAKdIAgiJpBAJwiCiAkk0AmCIGICCXSCIIiYQAKdIAgiJpBAJwiCiAkk0AmCIGICCXSCIIiYQAKdIAgiJpBAJwiCiAkk0AmCIGICCfSIwjnPdxEIgigwSKBHlMVb9ue7CARBFBgk0CPKgY6efBeBIIgCgwQ6QRQ4n2xrxrI66pER3hTluwBEppANvb9w2QMLAQB1t1+U55IQhQ5p6BGFxkQJgrBCAj2ikDwnCMIKCXSCIIiYQAI9opDJhSAIKyTQIwonowtBEBZIoBNESKTTHN95Yhk+2LQ330Uh+ikk0CMKmVwKj46eFN5ctwc/e25VvotC9FNIoEcUkueFS0s7zeIl8oOviUWMsToABwGkAPRyzqsZY8MA/A3ARAB1AK7gnDdnp5iEFQrOVXjQEyHyTRAN/UzO+XGc82r1980A3uWcTwXwrvqbIAiCyBN9MblcAuAxdfsxAJf2uTQEEWGo10TkG78CnQN4izG2nDF2vbpvFOd8l7q9G8Co0EtHOEKygyAIK36Dc53KOW9gjI0E8DZjbIN4kHPOGWNSEaN+AK4HgAkTJvSpsIQB+aEXHvREiHzjS0PnnDeofxsBvADgRAB7GGNjAED92+hw7kOc82rOefWIESPCKTVBFDIs3wUg+iueAp0xVsEYG6RtAzgPwFoALwO4Rk12DYCXslVIwg6ZXAoYejZEnvBjchkF4AXGmJb+Kc75G4yxjwE8yxi7DkA9gCuyV0zCCgn0woOeCZFvPAU653wLgGMl+5sAnJ2NQhHekOwoYMjkQuQJmilKEGER4lf21y+txe/n1YSXIdEvIIEeUcjnOd48tqgev3trU76L0a9YWLsPt764Nt/F6BMk0CMKifPCg1xJo81VDy/BE4vr812MPkECnSAIIiaQQI8qpAwSBGGBBHpEoe594UHDGkS+IYEeUUh4EGGytuEAdh/ozHcxiD5CAp0gQiKq39jGg52Ye/98fOGBBfkuCtFHSKBHlKgKj/5A1OYVtXWlAAA7SUOPPCTQIwqZXAqPqM4NSKXT+S4CERIk0CMKDYoSYZEieR4bSKATRD+nlzT02EACPaJks3f/yPytuOLBRWho6cjeRWJIVPtMJM/jg98Vi4gCI5vC4663NqKtO4X1O1tx2JDyLF6JKARSEbX9E3ZIQ48qWXwJtRecXvNgRFUu0qBofCCBThD9HBoUjQ8k0CNKLpTBqLrh5Yuoeh7VNbXluwhESJBAjyjZlLUsclNjCgsWtduntqWiRNQKTlghgU44Ek19M/9ErWOj9SwGlCTzXBKir5BAjyi5MIdETTDlnYjeL+05s8h1LbJDlE2NJNAjSjabHL3XfYPuH5EvSKBHlNwoEdHVVPJBVO9WVMudLSKsoJNAjyK9qTT+74PN+S4GERMMk0t+yxEHXlrZgA827c3b9UmgR5DavYew92AXAKC8uG8DWQc7e/Bx3X7psShrKvmA7lc86Mtj/MEzK3HNo0tDK0tQSKBHEG1i32FDylGU7Jta9e+vfIrL/7gILe3d+j4tR5JP/YOo+s8Tdkig93PeV7uH3b00XbC/optc8lsMIgRIoEeQbGhUshzJhBCMbDyXTXsOhp4n4Q65LRJ5gQax4s+O5vasXyO64ouwQgI96mTxbSTbajCyodjlJHCWWnCaWKQQ5VZPAj2CZLtHSC923wgzFk4qnYMZwVm/ApErSKBHmGzJXc2GuHhLU3YuEFOyIRhzac+lz3h4tHb25OW6JND7ObKXuK07BQB4cvG23BYmi3DO8fqaXXh51U70RCgAeC5WE4rwGGBWCON+LN6cH2WIBHqEoTC3/tmw+yC++9cV+P7Tn2DJFvlEqkIkJyYX3Yae9Uv1G3Lw2KT4FuiMsSRj7BPG2Kvq77MZYysYYysZY/MZY1OyV0xCBr2A/ukS/Oy7U6msXCMb5pE0qc85JxxngPw8tyAa+g8ArBd+/x+Ar3LOjwPwFIBfhliurJBK80j7mMoIqzYxuy02xOcepbrmYrnPCN2OyJCvNuZLoDPGxgG4CMDDwm4OoFLdHgxgZ7hFC5dDXb044t/+iUm3/BO7DnTkuzh9IkoCqT+RFbfFnNrQqcsXFvl6Rf1q6PcAuAmAqC98C8A/GWM7AHwdwO3hFi1cmtuMWCX3vlOTx5KEB71+/onqNzCdL2NsPyaMb2jBauiMsbkAGjnnyy2HfgTgQs75OAB/BnCXw/nXM8aWMcaW7d2bx7CSwnqJfQ1oRUQP8QWLUg/nySX1Wb+G9XZs3H0Q3/jzUtTT4tEZk69JeX409DkALmaM1QF4BsBZjLHXABzLOV+ipvkbgM/KTuacP8Q5r+acV48YMSKMMmeEuP5tUSIezj00ASgIOfTnDuGxVKjre7Z1ZWcAV8Tq5fLcsu14f+NevP3pnqxfO64UrIbOOb+Fcz6Ocz4RwFcAzANwCYDBjLEj1WTnwjxgWnAkhLesOOIauvb1D7MWcZ/mb9LQc3itTEnk8WNNVp6+k69bWJTJSZzzXsbYtwH8gzGWBtAM4NpQSxYy4utRlIyHhg5kx1VuUFkRDnb2hp5v3MmGVhY3r6z+Qr6eWyCBzjl/H8D76vYLAF4Iv0jZQTRPFCeiraFnA7H9TRxegTUNB3DxsWPzV6CQiapYzFa5u3vTSDBFuaF46Gai/A2Nj6rqgdiDnTpqUP4KEgJ6g5O8gak0DzS7UNaz18wvEW7XNsyDotmpWZhmKy0nt6Ju39+Oz935Hl5eFdxj+Ny7P8CXH1qsXitOT7owKFgbelRYub0FtY3OiwGINzjqNnQNWS2+8eel+Mafg69pKF/gIj4vek6DXIXYvNyE7dZ9bahvasdjC+sC51vf1I7l9c2mfTTGrhDlD1xGNvRC5NI/LAAA1N1+kfS4+JBiJKdsfFSzL7S84nSbuMN2qNfIig3d5Zj6t6/hAeL8PuSLQnZbJAoM3eISokolnRofoxc9qkLLVaCrB/tat4jemoKGTC7ZJodua7kkTnUhDAxh7fyEw/7uUvROhVjPFI0L/UHw7WntzPhcWQOMsi3RSi5MbllZ4MLPQTK5FByFHssl8kR16rcMfWafZf+76xtDyt/8NxZEtC7uNnTloB+npu7eNO59p0a6ClWcPtxhEErw3Dy9PP1GoMcRqwk97NjZcRLouahKmC+xYU7xztNPmk93teLudzbh+09/4piGvFzCgzT0LNMftJBMauhmN43TPTPL2ujUy31QVPnrJ2Z6t7rAx95DXYGu0R8J5cNMNvTsYo7lEe0WnBVbrWTQmF70YGTjdrmZU3TTmI98elWpnyQ1PCeQ22JI1DYeQkd39iPUFQIMzCx0wza5hJpbfsnlPIT27hQue2AB1uw4EEJuPrxcfFSoN6WkSbiEvSBRHx7k5RIS59z1AW58aoVtfxblXiyIeq/Fi1w/80+2teDPC7dmfL6fgekgpgE3DT1OM4LDIJRB0RDyyITYCXQAmF9rny0Zp0arB1OyvJth1TCsCSuFRG4GRc2/w1htyE8OfgbDNQ09adHQn/14O3731qZMilZw1Ow5iDve2IADHT35Lgpp6GEi9amO6cQiJ4J+wOTJVbe4NI/8UmimmbA5uuaLK3fig019W6XLz8SiTXsOuQqxVdtb8OLKBgDmhV4A4OO6/fp21BdM+d1bG/HA+5uxaLPdNTMI4YyJkg09NHKxsG4hEravvZbHVQ8vxtz75/c9w9hjv+nXPBo8UJpIc7uzoBafsdvC5/fPq8E/1+wGAIysLDMdi9NyjFoM/0LojeerCLEJziXi1QUthAeebYLWUJZe27d4y37J0WgRtTEUfxqekaa0KOmYqifFcey4wRhYVoTOHrOPY5CVkXpTaTy6YCu+dMJ4DKso8X1erojCc802sdTQvUwu0UedKcqYxXvD2PY7ycjtfY7Vhy8HVfnjB1uyfxEB8fGUFjm/yhxKW0kwZnumVpu6Gws3N+G//rkB971bE7SoOSG0OP7RdUOPp0CXEXcvDit+ZbHbxy9Odyzbz7+zJ4W/L9/hO319UxvOu/sDXP3o0ow/nOJZbh9wzrluO7cOhQS5dHu3YtLY2eJs3sknBaV/RGEJOqKwsOpWJrNCQAEmEyrZaJNtXb1IJhjKip1NBNkg2xPLgoZd2LTnkP6/J8VRUtQ3W7bb5dOcGxq69bwM7kWhj532PZxwbryTskH/0dBDHjDMJ05ui7I0XuT65Zz16zdx9v9+kNuLIvvPPGj+4jKBsiUD/eTnt02n04p3C2P2D7f406stFPp7o0+0KoC+JbktZpn8P+Ls05dGJNPus3XPGvLcZQ96n9q6enHvOzV4fFGdr/R+ljgUhXivn2AsEv7njQ36tqvJBYKG7pBscHlxwQtsTwqo/BRtMcuY/ZAL6MlniTAiL8ZpULQvNVlatx93v7MJv3ppHVrau0PJX3St9SvPrXMBtu1vN465PKs0V8xzTJKOA6gaWIJzZozyV4gI0GeTCw2KFh4HXPx3o47WWJQutH0/ELxRZkN2/+qltQXzUehLOVIpd/NIJvmnM9DQ75vn7F3iHsSLI8GY4hUlHRQNZncr1JWNCklRI5NLyLR0mDWpwnnUuaEv9Q1rgYvHF9Xb/J7zBXfY7ks+jvn7uICnDV1yTtMhee/AvWRKeRIJRQGwa/Lc9zhKob9DYXlnhVFP0tBD5tuPL0NXrxF1MU6DohpWTcm80HPQSsqESt9vVLODiSLXhDa+EFLbMdvQ5ZleePRo3HjmFNzxxWMwrKLE9Xm4aehpzsHAbNP+AaU+DMEGxwvVy6WQXmuyoYfMpj2HsH2/OPhWSI+7b/jycpGex9GT8q8xL6htQs2eg8EKZ+G11bv6dH42CB7nxnv8JbCXi3CCGE9FZMKwCvz089NwxWfGexo5PG3oTFEAbDZ0HlxAF6pC5Gdh7bgTW4GuEM7LF0Vkdbzr7U34zH++45lePLW+qd2WNlA5CuZDGlI5fNhc/Fxp1wFjQe/l9c2+ruMeStflGJQp/omEPR0H13t6/VkQioRxH8iGHiLVhw8FkH0zS3dvGnNun4d3Pt0TfuZ9RNYo759XixbLYLGTcjZh2IBA16ttPITaRrs2XyhBGsNz6XRK4/8C3b1pffp8kduCE8y87XYFr5miXhq6HyXdT88wnxRCU9PuDUVbDIjbV1QLOOT0IoZ1q/e0dqKhpQO/fnldSDn6Q6s7g7N9N7CXS8D9Vs656wOcc9eHtv1hL1ydKWEMElu3ndJ4oXm1nD9rNAaVFfk819mH3Ov6adXLBQzYvLcNew8a64pyxGelIu0e3PbyOtz28rqMV4zqS1vR7iVp6AHoSaVx8n+/i6eXbpMn0L6SMRwI9UvfBJgR+6Ov3c9Cue9h9aL8aF5e90w7PHvCEFuANScUHcXdk8XtWIIZva61Ow+YjjGJAhQmqTTHim3NWLGtGb0BxnCCcqhLiTXT3N6DxxbV+Z4IJuPhj7bg6NvexNZ9bRmdT14uATjU2Ys9rV24Q5gpJ5Lw6PaEZStkIQm9PhdAxSnyoh+sycN6yQtlYYyetHhvgp7tfW6QHqB2XDN1SAOkWXJxSqfhPSjKcN7MUbYCatfJphnltTW78IUHFuILDyzEK6t3Zu06mufQMeMGY+zgckdz36LNTbjovo8c2ybnwLPLtuNgZ69p8pZxnOPrjywJFIwtV0RSoHeq7ohOsZx1k4uk4frlly+uwbPLtrumydcKL7pA8JEmU8LqOjq9VFv2HupbxgFJpdO+puR74WiayuBGMaiTfRyPC9seRW93WRidc666JmofafMIrn8/9MyE/8FOY9zmkLoIRTbwa9678k+LsW5nK+olwlpDkyEpyaSv7lQaH9Xsw0+fW+V4PplcAtDcpjQQpxjQUoEe0AvhycXbcNPfV7um0YWej/xyzbwNjYHSm7R7iC9t32rn9JLJ1n3NJj0pHij2t4ivQFkB0ovC3zrbN9Ny/Ozv7sIlwZj0Iy0+62wJITHfbHbYtLZmvJfuF3NrDYZAtx9zu0/Sj2YOiaRAf3TBVgDAkAHyVVNkI83ZaKzadQpm4E8ohvusQm+Y9KMYvJ7WcyrLimz55oJUmqM4oTT3wKGFxe0QC254l/iwocPd1l6cdH6V05zrM0UB+zNl6j8vdC+XgMOo2bp/VnRlmjGlrhleisP4+Ms0dDciMyjKGEsyxj5hjL2q/maMsf9kjG1ijK1njH0/e8U0k1Rb5mlHVkmPyzR0EyHd7HzHtHBrtOLL/+a63Zh482uuedls6D7S+CFME0Vf6EmlQ1k/MwwvFzGpk4ZuH9NwvsbE4e4upmmL0Lba+7NuOTRNzMrmZcwael9IJDLT0PNNkAUufgBgPYBK9fc3AIwHMJ1znmaMjQy5bI7otjyHRydroNnoAmW7q+qEr4E54ccj87c65iUdB+AAS0jy9F1CA2vvRXtRcv1OpNIcRaoWG9il05fJxX+mYn7MxR3R5IcO53vGJEvLiXT2pJXegGQQX5v6D2TPTCDmmg2TS0t7Nz7YtFeffavVNeNLcUD79ssCp7ndp3w7SvjS0Blj4wBcBOBhYfd3Afw75zwNAJzzYEbbPqDdK6cb62VDDwvtRejscR6QyiZ+NZEgk0aMc8IxuVhfYOawP5vsae3Ews1NrpN43PBlusugPkw1DfhzW/TwQ3fYf9fbm9DQ0oFXhRAMdg2dZVVL72sb8uLHz67CD55ZiT2tin+9EirY/SPnhWFyySyPQje53APgJgDi5+oIAF9mjC1jjL3OGJsqO5Exdr2aZtnevXv7VlqfeLothqyJtHb2orWz7+F6G1s7sbbhgO+GaNWuw3xxZOMQTm37icX1jvk4jS/kUoOp2aN41MwaO7jPeYXSdkwaur+Xv6GlA/9YsUMaJsAtj21Nhh+1sw3dH/ppgWO/ZHcsq/Fgp+l3X73Ptuxr0/Po6k2jsyeF7l5D9LkOivoeFckOngKdMTYXQCPnfLnlUCmATs55NYA/AXhUdj7n/CHOeTXnvHrEiBF9LjAg3CzHrqqHDT0LtLRlLtB7U2lwznHZAwsx9/75WLGtxTV9UKHivlSdPS9Na1OOG/udtJVbX1zrkr+1LLkfd9A+KhcePRqAk82aY+LNr+H/3t8sOSbfNqXJoFyaK6HsXKf8nl8h8X120fLTpo+H3S+LawXJIuYeQfgvpWxA2I/JxakpPr6oTvegu+X5NZh+6xuY8as38Mba3QA88s2TGVbDj4Y+B8DFjLE6AM8AOIsx9iSAHQCeV9O8AOCYrJRQgmFykWNo6PZzrNthkekyYm1dvZj16zfx83+s1pdm22vROPziZBpwC7MrO1c5R9tvsCmDyIvWyRv58ALQBLqb26JWzDvetE9W8zOOEGxQ1LD1iuWzIhsj6k3Z05Z4eLjo+enXMxUmo2eyvL7Z91yCbL97VoHO4L/nI3LkqIHq+QwjBpYCAG46fxpuOn8aUmmOT3e1AvDXuyxYt0XO+S2c83Gc84kAvgJgHuf8awBeBHCmmuxzADZlq5BBkUWPy8YNFnN0imntxYGOHnT1pvHsMkPzykQ4uJXNTSl20pi13WKXPZXB22izoefB1VN3twth9qvTy5zxoGhA9zqxnY2qVITOzDGVvrxvpCYX8MA29HSa44v/txCX/GGBr/TZHhSdNmoQAGCQ6hJ73qy+LafH1ac5flg5bjhjCm44Y4p6gKvHncm322IQLxcrtwP4K2PsRwAOAfhWOEXyRl/E2OGuJWQeGqZGHFI5hIxkmpMf/K5WI0tgfQedNCE/L6vVrqoNLL+2Zrc0jV/qmqyxMHJvY9Q1dB83wskco2+HVirN5CLP06ltiz3BoQNKcMy4IYrJxaFgJg1dvffiPtHLxS/a+Qd9zvrMtmI1flg5AGDxLWejN8VRWV6Ep5duD3wlvefPtbEF486Iz8mf11N+CCTQOefvA3hf3W6B4vmSe4QbLyMfNvRMR8OlAt1nVu4LXPSt8rombe2fB2TehkYs2tyEU44YbtqfW5OL8jdDJxcTXprwmMFlGDO4TB8HGVCStKcVtt28MWTP16o4aB4dTshMLlZFR7ZfRqHGCTeeL8PgAYpIU0wuwS6mpU5zbvPPD2zCyZOK3hcNvWAxYrzINauw7rUoNDO1ocvMGNm2v4kvsWdaPy57HrQIy9B5eSBlA31KuP6hl31E/ZbH3cT1w3Om4sufmQAAuPj38zGswj6bmQvlCeovvXJ7i+m37nPtUH5N2L39o9Nt11fKrWmi/r92ffHjz4Z3k/F8g51nLYqx4pG952KKiulSBQ7PJFklklP/vW6aEfpVOCcLDUnMcuPuzJZqk0V8sxb1vQ2NOOa2N/HUEiVcsHbYTTNzGxT1QtFOtO65eb8MJ813zOAy23n5mIylPfswYrk4a+iqUBG76S7pAWPqf5B7UVlebCuLmxmec45ZYysxddQg6b0XV07KFtkOwWGMkQg7M5hYZJUr4viS+Jz8hVAOePGQiKZAF76kMmQeGubtcO62mMuGDAW6TEO3Dhiu3nEArZ29WL2jxb08DrbeoDZ0QJwA5P0yarutE3fcQjDk0g/danKR2qyF7Yvu+8g0WcwauMwVs1rnea2gbovmNOqUfhcbeirNhR6r3Ya+escBY/auz0eSqW1auXbAk33lr2SasAhgL2yhIARTrrWYQW3o+SKSAt0Lz1guWSBT92rNhv7t0ybhd5cfK03ToQoX6Sx9Z9XM8Zp+bovTJBS3S8leAmW/YMeVzEDNNpoAS/jU0NftbEVTW7DgZrL6MAC1ew7iycX1+MuCrVhev9923M3WKyutLK1b4K40Nz5kUrdFAP9+yazszhQ1bYf74O95ZxN+95biYGcLp+1wKc3N09lKym0TrsQQDa41CKDFZ4NICnSja+TwIkhspdnwhQ1jBpwm0E+aNBwnSNZCBYxV4cVReAC2N74vnjwcHH/7eJs6U9V4OUwaulceloJfctxYWxonoZJNxEEzJ5zsqdZjXs/ZeoWdBzrxyxfX4rZXPtXDMRvPj2VkGhDLYtjQ5WnSnOvvg+FSZ0585MhBga+bafqwP+TaLGDAbPrzsxLU6Xe+hz2thslJ1MA5YOltBRPS+dLiIzkoahNsFuRd63DvcFdvSr5aewAOtPegvkkJsp9MivP4zGXVfI+tphi/SpWf2Zk7Wzrx83+sAaAsVebX5CIz8wwuL8alx43Fpccdhj+8t9liz9fS5q7FG11y/eI+znHY7zQoqttxza5ugLJo+cjKUqxtaLXlwXyWxwmvQFTcpKE79I70InNTucPCaSWtf67ZhZ8+twqpNEdJUQJ/+eaJulLjF9FLLJNZyI2tXRhVWWYqm0xfEp+Tv4lFCr2pNNrUxUcGlCRdwxyHQSQFuhdhLHDhxQsrGnDz82uE/IPnet49H+gBhQYUJ50HDD3GDPRk9lMAyPzVub63olRxq/veUytMaQy3RTF/ewHS3EhvuqYwWcVqM7aWL9v48UN3+8D4GhTVhKGwT9suK06irCgp/SA7aZJ+7o/4WXDT0BNWDd1yPSMuu99rS9pBmuNQdy8qy4rtqR3evQ27D6K9O4UrT5yAp5duw9Z9bYEFuvMsW5d6CJWVjY9wiRHdZEN3KY8xR0b5ffmDi/CJ6sL6l29+BmdMy25Q2kibXJyQ2W7DRluQ9rl/PQVDBhRnZEJobuvBuTNH4dFvVKN64jBh4oc5nVNDcpJPfiP4AcBnJg5TzjFdz5hUwV2+hD/620r8+mUljktCMHMZ57i4CfoqXThoMa2DaHBO1fY0uUguwVRjOZc8SHfB415efXEKFzXfJNAdzF1+74qsZ9ylLgd537waHHPbW9jvMfYgCuBUOo2iBMMNZxyh5hu8VTi9d25mKKfzRRu5NoNWz0+YL+DvY6sk2tHcgdkThuDWuTNxxIiB3if2kWgKdP3G2u/skaMGYuZYNWR7yFq5uQzK3xljKpHwYa870N5jC7Ob5hxHjhqIs6aPQjIhaLSWemk/tZfBs9sPc6Nzkwvay25z6ZNo19arvvBJA55cvE3Nx1wGr+tyznGgowe3vbwOjy2sc04YAnXqyu1JPRa7t0YcNFKndFBUF6RMar/3sn/7xe2joCwQraWT946svaYglos5t89DbyqN19TwvNbIh0q+crNdb5qbBqozuQ0ZhZAQTpGdz7k6PiHsE5+Tr96c8PGbMaYS1506CeOHuS9EEgaRNrmIt3X8sHJUHz4Md3/5OCzZ0mQ7zh1/ZIbuOcG8fYmfWboNNz+/BkMGFGPJv52N0qKknofJ1crBrS5ttCQbsstaNVG399PaRbSe48dt0bimuakbg3DC9QQtZ3n9fvxFFebXfHaiSyn7hvbB0MxLfpCV2brflF79a51dCBhtxPgga8eZq4eKZxnVi7jb0DmSaiwMpzg6DP6EuKyc+w51ozuVdo0frl0uwYDHF9XjpZU7kUwwVJYXoUhQZDK5DW6BzYL6i2vp05zrA85GfsGKJ5pvsuhAZCOaGrrDfu3GhbUeptt5KV2gu6/cvnH3QX3FoJb2HrR1GVq6oj2JAl0u0Z3kuS2KomDHdbfzyfaZhZbh5SJPYyWp33MuCDcXjRjcM/7NtX/5GEf+4nXUZBDlUaS0OIFTp1RhqLoGbTD7tL8XWTqxSN1MMLOvuLUXw7kymzbTNuoWPiDNjdhGTgpD4BmW6l8tOFiaG70fmSuglv7n50/H+UeNxkmThmHb/nasbWhFUn1/lHQ5MrmINnSXj7X5WQrjFH7aj/DxzqZLqJVIaOi/eWUdXl+zG3/99kmKHUryYoh42dCDNJst+9rwzNJtSCQYrpszCSP1EXHluPGyynP9/D0fmn73qAZdm+cF7INWNXsO4jtPLMcW1WSga3geXeMEYxah4TIYqJtzzPtl5h83eWP1LBIH2mTnce7tujhvg7II1o7mDkwd5d+1rq2rFw99uAWpNMc350xEKs0xqaoi0Ivl5XPvhOwaCQapWU5L29zejeP/42184fhxprkIsuKae5qGdu1ULLMNXT6mYV2UwXPw3TLIzDnXJ5VJQ1mou7512mQkEwypNMfznzQoeSSZ1Fznl3SaY1BZEX550YzgJ0Pei1fGRO29GL036yM/Iy/u6i4bNpHQ0EuSCexu7US9JXKfox+6dtykXWbG35fvwJ8+2ooHP9iCtz7do+/Xpuz7MbmIaCufGPZuu0anCbqNew7qwlw8xwlxLNLmGueBrRuuntTRLc6YVMIQ/Mv98/HSygZTepNnkeWDI3sOaW6+pjUEwo7mduGcYE9vydYm3PtuDX7/Xi3e3dCI3hQPvEA0d/jx1qe7rUnt6VWMZ6t89LUqij0pxhj2HuxCmittLWgZGWM+bOhMvx4giVHP7L09r+sCRr5KL0AzuTivwaldIZlgukafZPLFq/2S5hzTRw/S4+fIyul6vjTaqd3kIg5qB+nhpbn/QecwiIRAn3uMMkFF81bwGtSSdS3N5hf/1+4Vlv02CSBBQ7dqxCJfPH4cSpIJ3PElZf2P7c3t2L6/XfctN2voZlORfaDOjKOXiyyxmI+DxixuawK6TRTonOPDmr1Y03AA720wLyFr7RU5rTKvUdt40GRvtcaTv+yBha7ldUM05aTSHD3ptGkRCFl2tntt+gh591Kk+3WTi/JDanIBcKjLeSawH9x7YNw2U9TLE8UpO2sd9TDVgoZ+99s1joP6Yr7a8zA7A7gWS0oqzaX191qDVS+bqZxGD9j2rpnOcc5Yy+OR+VuxflerYkMnDd2M1nCsJgfbA2PmjTDiMIuCRryeGOFN0b6c8x8xqBQDSxXr1lV/WoLT7ngPd6or44ij/NZIhNYcrYNqVkRzh4irt4nzIWlaTaOxmr9l8UDczDZFiYTFhc2cYdOhLv1+BHUJFdN39aTQo2rowYKUyS/q3I609mBcQ4syqXz07Xkyptj39x1S5iJobcSvLV0ccHO2oXPdNKLV/755tbZyBLmuhmFyge6SN792H5rbzcsxWjV6QKk3oKw2pPccfF6/bl8bNqurJXHuPr8gleZ49uPtOOiw5q+Tl4u1vMoHItj92ba/Pec29GgIdAfbn9PtdRr80c8L8GDELplZQ1c0H73L63Qt9cikqgrT/r0Hu+yJLZqKrZxSzdpJ8GQBbthIrV1r7aVaVteMg1291tNspDg33U9r+OGiRAJzplQp5/t4Xuk0xx/eq8Vb63ab0t/2yqdIpbnuWaTkJ6uaXKu0bjt9XBqFFec1jp8wFICikJhNLgb/ddnR+O8vHA0AmK2m15AJgqDjQum06D4pT2P4sgdD+4inOceE4YZLni2UtOSG/+biWbj+9MmK7dvjfRVJpTnO+N37OPt/P0BHd0p5DyVSTOsbztvQiJv+sRp/XlAnz5DbNzWTiyk/YZzCqTnuO6SYzs6dOUovK+fBzFl9JRKDoppA100uHk9ev32m/lRm1xYHeawaujjY5FYmxoDDhpbrv5MJhh7d5CJoAR6DU/aBGktD4Vy3h4q++u+sN5tH/PRWZAM5HFx/BlYPlYFlRWhq68bXHlmilI25ayaptJEXAHz29nm4de5MXFE9HoCqWSbM98ONrU1tuPPNjQCAP1x1vOnYTedPwxePH+cjFwOnpuP0nGsaFY1RjH+uCTzdPVEyGDljTCVmjKnEM0u3BX7tNe2PCRoF5xyLtjTp3lStnT0Yx5S2Z42GqZfD3ozcr6seFz2hTO+GVZ5LrnHJcYfhkuMOAyAoNz4+3GJP7sEPN2PXgU5MHlFhS6d5uaxS48fvbpWHCZZOLOJKqS0mdEdvM0AxIZ57t+IA8bY61tabVoN8kYZuJmk1uWgmCQctSuYGJT6EzXvbHLtgVsQGJOahRLEzNB9HIal134RdyQRDKqUJdLHc6ilqZazdwXTafNxR4xIGcGobjeBF00fbPUWc3iGpdsiNHou1bF/5zAQ8c/3JOHu6ZGqzSQtSfvSk0npeFx0zBl29aaxrOKCnSwt2WT8dqpRDTwoAbjhjih6vwwk3G7p5v/yAdr+mCfdYKz9jynM2PtT2RiG6vsouUSwZ1FW0P/UDru5b29CKq/60BN9+fBm+/fgy7GjuwGA1hvrwgaV6b8BUdi0/Sf6y39rupPCeifdF5uniJtO8etQi4rO9550aNLR0YNxQ+YSdxVua8Pv3FPOSLCSBck1u2+aAbVBUDNEgawP7DnXb7lc6rZwRxkpZfomEhm6MptuanDy9dtRB4D+xuB4JBvzmkqM8r51KK4JF+9pqpIUvr5vPq6adiBpvkjFhUFTU0M1l9dLQZTAA7d0pPDx/Kw6vqsAxhw0GANx/5Wxs29+ODbs3euYhlsV0fcHkYh3ELC1K4OTJw/HKqp36+W7+xb0pruf1y4tmYGHtPpO2JH4wfZlcTOMl7viacOIwEOp0plVrBQzfbG2ugvYBk1WHMeDDTXsx8ebXcJb6UfTbVRfdZtu6FXPXf152FI4dNwQAMGWkMeX8yhMn4BYhBpFyvnKltq5e9KTSnhqlVaHgNg2dW9K7D9zK3lcntOf8g7On4tpTJwEAKsvsYowx84D+xOFyoa8V9canVuhxlcDtvQqThi4pp8wW35vmpmiXuSASAt0aytXpxmovgNuo+Z1fOgZ3vb3JNnDjRCqtvJi9aW77QBhBj9wHTBjMU7+LEky6bJbVT9hLY7JiPXzri2vxheOVbq3zLEl7ppwr0e9kKfVBUYnrG+D+gRLpTaeF2bbK/TGer/JXE4gNLR2mnoZGRWkSYwYr5gRTEDEnLVr9+8zS7Thv5mhUlCbx7vpGDB+oLrYs4BSUTObmBgiD5MK+Il2ga/nIy6OkMX5t2NUqKbtiVzF/XLhtDEcrxxEjBuIo9WPuBQOwpuEAelJcmcyl+vx7ffjE91K8LR/X7ccF936ku+j2ptMokhm6tesH+HBr7W5gaZHe85DXySxEZb0GrewA8Na6PZg2ahC2NrWZPLWMMrorCrIJVal0Wu9F5YpICHSta6ebHDzSy9zltIc0saoCg8uLTXFV7nu3Bh9u2ivNSwsg1AXzV7i28ZBp+r+z5mYX3MmkoaGLX2/ri28zuVi6vDKYpbvw/IoGNW/jOn40ocqyYjxyTTWue2yZqS7GoKhFoOvXF8ui/H3m4+34VBVSrR2KBrliWwvW7VT2Ka5rTK+f7hKq3pDfvrYev31tvbScb/7wdEwbPcg2YO3GmoYD+MvCrRg/dIAeMXPRLWeZ0gQdfJRp6AlBQ09IJIJ1iTNXXBKIj1xWDs+sGfQP+ILafbpA9yIpmMTE+7VpzyEc6urF108+HANV7XnqSOfAVG4rSVnRZIDfxUqM84zcTWdypS13p9K44OjR+LhuPzp70pKPi/s4mZOGrmj6pKGb0O7HexsbMbKyVBcmXvbfW55fgzU7WvDj86YZrkgASouT6Ow1PqnPr9iBQ10pTButNLr9bYb2nuKQDs4NHVCCLjUPURhZ0bpu4jMtSjDdS0Tmhz6/Zh8YgAEl1scj0abFbZcGJ/r7erFTXWfy7BmjTPvrmtr02ZsLNzeZjg1Rp9WL4wpafVZub8H6Xa0YUJJEeUkSZ88YqU/DHz24DMMrSkxufbJwtz8970hMGG4Mfm3Zewj3vFODvQe7THZrwGUlGqH+bV0pfSUoAOjsMZ/kZK7Tytbc1o09BztxxIiBKE4KLpiW56zt83JtFQWwXxu0bkMX7LuZLJjMGMODX6vGtx9fhqV1+zG/Zp80nbVUopeL2P41zfzGs7zHLgBBAQtgcvGaJ2atv9PauGnO9QHTsuKk3tu2mVyE9ullchleUYKmtm4s2bIf3b3eJqwwiYRA1wTqq6t34dXVuzyn9mvsO9SF++bV4sazpgq+sEBJkqGnV5wwBJw2tQp3f/k4AIot8+pHlyrH0oa3hTlQFcf4YUp3n8G9u8hgt61qXiKmLrraSt/d0Ih3NzTidssAls3EJGkpTm3H6qvb3t2L+TX7sG1/uzT9SDVOh8h//XODqQ7ah/W175+KmWMq7WURLvmNz07ELRc6T89OMKYLYu0+i14ZZ0wbaTIhrNregnveqdHDtwbR0AGgO5U2veTiBDI3tKzn3j8fDS0duO7USbh17kyhHkbaMtVVsrQoYTKLSGS/6YcR7lc47PRgmXKq9kGyrp/ql8EDivH/zp6Crz+yFBvV2DmcA9v3t2PM4DIUiROzLNfgHKbGqT0T34tyB9DQ9RhKHnlrRxNMuSeyoGGAsv/M370PQOmVah0pq6nEZEOXmimVvz/7/DRcOvswnH/3h3hZHU/K5aBoJLxcrN1Hr/fV5oZlegBmey1gHuAEYLIfK5HkErbrprgQo8HV5KIn0REn1YgNs6K0CL+5eJZQbjPWqIvWduL2gbO+AE8v3Y7rn1iOBbVN0vRW7dyK+IJMH12pf1xEV07xil4vd0LQYHWBYZp0ZT5fe0aaNihzP3Oju9fcre6xuGE6DYRq5zSrk4YamjvU69tNaFedNAG/vfQo3HDGFCQSzFUgiLfHad1Qe1nUuqhfgNU7WgQTnz8pIiY7beoIfXIToMwoPe2O93DPOzXScokTi2QaerGL3VxWhiCD335NStr8D6ePfG+ao7s3jaMPG4zLZlvcKEWTmMSsJSvXaVOrcNiQciz+t7Mxe8IQ5dwcWtEjItDl++0aq/oX9g+A2FisAl0c4AQUjUrj7U/36DP5zF4uRoMW7aOcc/xj+Q48/NEWPPzRFize0gRrPOxEAvi4rhmAXdBdcNRox/r50WCc2nkyYZ4p2a5O/nn1/52KH5w91Zbe6WX85pyJ+NE5R9ryll6fydPIyy3a0O0mF+v5mkB/fe1u0znW7RvPnCIUx8jj1dU7TXZ562QYk9AV24r6V/ug9VrMf2Ixhw8sxddOPhzj1SX9nOLlWMsmNbk43D4Ghs8dOQIA0CS4zvkVeNZU4mmtqmvvK6t3mtJo98Y8KGqUWYt5lPQZP0dL5cvkotnQveqnHteWfEs52NA3q4Pt5x81GuUlSbR29mJHc4fdvCS6LcrKZbnvFaVFGFCSFIuSEyJhcqkolRfT6fnb7WfcZHJJJIC0sNaEGO8CgB5RUaO0KIGu3rTNNUtQ0PUGXd/Ujp88t8p0/uTSIosN3RCWn5812pRWFFxWAfD+xr34ZFuzYz3dXgirPNUa4MwxlVhWt9+W3imYVVlxUh/o8rqOKKQ8NfSE8cG0Dooq55vTjxmsPKNNunnA/LEFgI9uOhPjhAldoi/3hUeP0QeMAbsbptO9TOtlVP5qYyH6UnxO2pjQi9N7bUJS8fvpZB5wQpvMlBLapG9rh6URSSeUORRHSyq+XwBwoEP5EDhNZHIqgx93Uv1j76GKalfWBPpWIchdSTKhj389sbgegFHvow+rxKrtSk/HyeTiVi6T84Meh540dBMVpUW46fxpmDXWbqfVEF9o6YCIYKawm1zMDfmwIeVY9evzsPW/L8Sfrq7Gs985RbmGKU/Dts6YEvf8+RU79Ibyu8uPxZdOMGYnig9Vs1szBlQNNNuqRWEvaz+Lt9iFrx+sXi4poQGeM9NuXpFNZAGA4gTD+UeNxjfnTJQeFyP7WQeCvcpnFZZFLiaXASVFuOjoMXo4Ylnc9tKihOm+DxlQgvuunI1nv3MKbv/CMab8OrvNq0lxKLP/5t7/kSn2idZstOvpGjrsL7QIc7PLwfwh0DwyHDo7Qlm4aY5DShiczFRDFx+TNtvUKSyC6CyQ5oqr6+JbztbT+bWhB9HQtY+dXyGp9bY3CjH1T5o8HIDynhfrQcKUY8MGGDN9rS7Fb6zdjbP/931co46viRjhsAUlRFD4ckUkBDqgzPZ78OsnmPY5ag6WW5g2CXvVhm6ZWWhtIIPLi8EYw7kzR+FodTDOGkxKe3gzx1Ria1Mbbvr7asFPNomyYm1lAXP53Lx0xG7qPkm8F0UbMpuLrnl0Ka55dCleWb3TUUO0ermk00YsmnFDB5iiEQJw9BsuSiZw2JBy/PpfZkmPm80I4vXdm1pCMLnwtFFmsfxWSooSuu3b/Dy1stjPufjYsThx0jCbBvs/b5onXHHOsabhANY2tGLmmEpcO2cSqgaW2Ewu1mfpJGeKEgzdqTR+++qnRm/R4udsLb9fDPdB8+Q3P1jLKwokbVDP6i5sFeg/e24VXl21CwyK19JotYfrFjRLVgavkqfTHD95dpWvvLXDRUmGasvC08VJhqkjB+K48UN0U5tWb61X+FHNPtP7+Z3PTcbp00Zg+phKHD1OkQdFpt40TPkAMCl8uSISJhcN2+CopQloR22miLQ5rTbyraHEjHa+rjFoYz5HK889X5mNw4dX4N53a0yDNkaUO/+IjeTed2tsx8Xu+JwjqtDa0YMWtYs7rKIUp06pwtNLt9vOG2gxW4k9DAC44cwjsLahFe+sV+yfgxzMKp62cK22lmReGrro1rdws+I2Zxp3kDygkmTCiC8v7K9VtTG3S4r1KE4ybNxtnszDYQyU/vbSozB+2AC8snonOOemj4dhQ9e0avlFv/yZ8fj9e7VYv7tVOjdB/PjIbeiG0BbLyIS6pNLQVTT/Gro5XRDzwHHjhyCV5ujsSWFUZak+w/XJb52EbfvbfPuK+3Vb3NfWhaWqefAEi5C256mQTDBUlhVjxbZm/O9bGzGqsgyvr92NqSMHIpFg+jPW19Y19eiMiUtXnzIRV58yUf9dXrwKC2sN905xToqG/pHIoUSPtEB3+qRbb5+2RqB2zD4oyl1ffunLxM1R3jSBpXkcKEH8g3eAnATfCYcPxfL6ZpNAP2nyMHz79Mm2tFaB/tr3T8XkEQPBuSKsOZQJQuLL+0N1oHPiza8BUAb0RIYOKEZze4/nQrcJByFVNahEklo8T4lH0ptK47t/XQHAWOIMkH9ISouVhU9O+q93cJ06DRyAHu3RzdYvlu2uK47Dc8t3mCaXcW4EINMGYBNM8+gw99QAb3fB8cMG4KRJwxyX3RPPk81GlZtc1DEh9WCKcyT08Qd5ObwylpXfSesfP2wAfnLeNNv+KSMHmsINeBZB19DdJbq22MpdVxyLiVX2gFwyEoxh9oSh+GDTXtwvmM6mj6lEUYLpJjut3mIP+fSpI5zLbPkt6xUmMlDo+kpkTC5AkIEe82+TQGdKt8oazMnrK2qdDZqynKNpI5rGmEgw3S5n1XruUf3dZZeUCa77r5yN3181W9divTSZSVUVmHvMGP33rLFKF/GFT5QVcV5dtUvxr5cU4DcXzzL5VZ88eRhOnDgMS39xDtb+5vO4+NixputYcTK5XKpG1nNCEZZcH4P41qmTcJFQB9l9ufLECbj8hHHY09qFT3caGnYqrcwREEPmujH3mDF4/NoTceeXjsFhQ7RBVK6/7HqQLSiKgKhB72zpwAPv1+KhD7coaVzaUVFSaXeyxyeepWn9i7Y0yRNYEGPeBLWhO+UlYi2vdQWiviLrAYvsbOnA5+58D5f8YQEA2YQ7WZ5M//vdM47Axt9egG8IC5Hf/oWjkWDGfBDt/RXfCa9VrkxjapIFazSFZNhAd2UmTCKloVtflpU7WnDby+sAQDc7AMDhwyvwtZMn4MnF2wCYzStKXBVn84nbtd0GUrWXXhMCScaE8KlmzpymdE0HFNsFDmMMn9x6Lu59twZ/UVerP/3IERhcXoykpQxOvPnD01GUYHh19S7T/rp9ymDs5n2HFLdLiZC8Rmj0APDM9afo28UWO/tr3z8VXZYZlmJ8GyeTggxtYpEm0McNLZdqOyIzxlTi5xdMx3PLd5gCMXX3uscOsaJd5/Lq8RhVWYarH10Kzo1nqU2q0TV0ocq7DnTijjcM+7u7mSeB3nRKKrhk9auRxK/ZvLcNP1W9qFo6usHABJMLR4LZBYsb1mRyDV3bsPwOCS8/7c17D6G+qR3nzRyFCcMG4JQjhvvIU0GUyeL8kuJkAskEBA1dFegJ8Z12iz9jvg/GJCLj/FvnzsR1p07WJyDmgkgJdGtjazrUjRfUxWaLkwkcM36Ivv3bS4/G9NGV+OWLa/WpvEY+7hOLZDBYPgJps5kmaRHoog3divblHzVYPi16aEWJKfCQ7oqWYKYY4k4vgtZwJ1dV4HPTjG5jUZKhOwWkUtzk4pYpA0qKMMCifDDTtv8LMMbQ1NaFpVsVrbS0OOlrYlK5+lFcvaNF39fdm/Y/S9FWDuXvr19eh0Oq6Ubz+NF85bW2c8sF03HNZyeCc2DGr97Q0zhRZOkZimnL1HqcNrUKk6sq8NiierPnlvq3tCiBRWrYheEVpaieONTwcklzaJ2STF3lZOftbu3EG2t3YXVDCwBDMw3LHU/LxinwmaY03HjWFFsQNS9EASsO/BcnmR50T0xnUtJcNHQlZLFRXi0ct+gmW5RMmBb+yAW+BTpjLAlgGYAGzvlcYf99AK7lnPs3mmWIeLOPGz8EL35vjq/0SgB+YyAqkWC2hSu8TS4MW/e14R110sT82n04ceIw/XhSN7loXTjnUe6K0iLc/eVjcfJkZ01DdBs0tA3my+Si8daPTrdoHMq2FtYzU6HnhuklD5D9oNIiLK3br0+4GlJebPbpdXg+5cVJnDa1Clv2Gn7G3am0bx9oKzPGVOKMaSPQ3pVCRUkRTpw4TP9oMKa8uJpPczLBdEHsB0OA2E0W/3bRDJwzcyTmTKlCgjE8tqje3LNU67/g5rNsrq47W5TZquIEOr+1d4qaCSgfl4/UuC7/+uQKfb/mChiayUX9u7u1E2uFmPgatepyc0HutVYPsT2KGrp1sp8m62Xvi1P+4nvY0ZPCWdNHmkIk5IMgGvoPAKwHoDuDM8aqAQwNu1BOiKPmfpQDI4obx6MLtqr7lAfZ0Z3CpztbMV2N1Of1/g8eUIzX1+7WZyYCZg1bF+iiyUW3bSpppo4ciDZV67tstvsKOmLjM7SHYJNOrI1LM5k8sbgewytKfLuVBUEc4NI0onM8wggAwO+/Ohu1e5QXt6QogdkThuozdAHnQb5EguGJ604CADz80Rb89rX16OqDhl41sBR/+eaJ0mNlxUm89ekefSZkeYl/AQOYg7JZOWxIud4mmtVFnM3avPpXcq7WPp5dtt0U4dEP1kl77YLpShOGPz73SJw7cxTSnOOi++brQr5YEmI5ExKMobQogb8u2Ya/LtnmmM5pkQo3xFfgtKlVWFC7DzPUuEPmGc7KtqiVW02MIowBjQe7cOkfFiCV5qjZcwiHD/c3UJtNfAl0xtg4ABcB+E8AP1b3JQHcCeAqAJdlq4Ai5lmIftIbGrrmnnTkqIEYWJrErgOduPC+j/Cnq6t9DYq+9v1TseeAIWAYA6aOMjolWuN4eeVO/bfm+qdNAX7zh6f7KLXChUePwa4DnRg7uFx/6ayDuUHlsdZYxw0tx+jKMk/Xr0xIGBId5SVJvHDDZzG5yrvzNnJQGUYOMpug/JhcZNdWbOjhf6x+f9VsbNiluEQWJZnupueXZIKhvqkdNz71CQDn56dHMZR5u0hOGlpRjGPHD0GDqqnPGlspDa5m5d6vHGcTQprp4No5k/SFmGeqy+QBwCPXVKOuqR2lRQn56lQZkEgwPH/DZ7GzRb5MHKDMhh3tYKKUccoRVahraseJk4xe8OwJQ/G37xhjQqZ5Dup9PXfmKGzf34Fkwt3LRWPl9hacOW0ERgyqwhePdx/4zwV+NfR7ANwEQIxTeiOAlznnu1xXI2HsegDXA8CECRMyK6VK0JF70TbX1Wt0iW4+fwZOOaIK33/6E7yxdjd6U96risgEjsjkqoEoL07ivY2NGFxejLFDyjFr7GCTC1eQGM7jhg6wTd5JJhhqGg9m7MGgDfI8+51TfIU1zQStZJoosi58HARRA3bTljS0F3RHcwcmZGEgavroSkwf7Txb2YuLjx2L1k6lhzZh2AB9EWkrmkljhhDBskgPEGcX8qVFSbzkYX4UGTqgGO3dKX1NTxFN8Tl2/GBdoIu9I6+gbZkya+xg3RsrDH587pH48blHuqa5ono8elMcRUmG09V4OCMHleHmC6b7uIKq0ScY/uzQo8sHngKdMTYXQCPnfDlj7Ax131gAlwM4w+t8zvlDAB4CgOrq6j6Nj5tWw/Eh1ESTR1dPWn9RBg8oxgVHjUZlWRH+sUJx5avqo2vRKUcMx/r/ON+2/4xp4WgxADBiYCkW1DZhQW0TEsw88cEPj117Ip5fsQMjB3lrb5lSqs6OdQodEIRBZcWY95PPIZlgvgT6lJHKR7WlvRufn+UteN79yeccB+KywXmzRuM8S+weGWXFSTz7nVNMcd5vOn8a7nhjo2NcoyC89L1TbcHIZGXQyGW0wFwyY0wl/uNS72UoZWjipzQks1NY+GkdcwBczBi7EEAZFBv6OgBdAGpVwTqAMVbLOZ/inE3fKS9J4oYzjkD9fsWFyQtNs7j60SVoaOnATCEWTHEygY9+fhYOtPeAMfPodKHy4vfmoOmQYl+tKE3aJv94MWXkQNx0vh/tI3OuqB6PweXFoWlyk0f4H2ufM6VK+lF14ogAeXvx4NdP0F0uw+DEScNMv6+oHo/LTxgXimeJm+fFxceOxXsbG3HEiArH2cKEYf7N9yCoFRYk9oOqof9U9HJR9x/y4+VSXV3Nly1b5pUsNHYd6MAdb2zUA+5/9aTDMWdKVc6uTxBRprbxIB5bWI+ffn6a6/qd/ZEV25rx5KJ6HDdhiCkkQLZgjC3nnFd7pouzQCcIgogDfgV6oD4V5/x9AO9L9mfdB50gCIJwp7AMQARBEETGkEAnCIKICSTQCYIgYgIJdIIgiJhAAp0gCCImkEAnCIKICSTQCYIgYkKgiUV9vhhjewHUZ3h6FYB9nqniB9W7/9Ff6071duZwzrln+MecCvS+wBhb5memVNygevc/+mvdqd59h0wuBEEQMYEEOkEQREyIkkB/KN8FyBNU7/5Hf6071buPRMaGThAEQbgTJQ2dIAiCcIEEOkEQREyIhEBnjJ3PGNvIGKtljN2c7/L0FcbYo4yxRsbYWmHfMMbY24yxGvXvUHU/Y4zdp9Z9NWPseOGca9T0NYyxa/JRlyAwxsYzxt5jjH3KGFvHGPuBuj/WdWeMlTHGljLGVqn1/o26fxJjbIlav78xxkrU/aXq71r1+EQhr1vU/RsZY5/PU5UCwRhLMsY+YYy9qv6Ofb0ZY3WMsTWMsZWMsWXqvuy3c855Qf8HkASwGcBkACUAVgGYme9y9bFOpwM4HsBaYd8dAG5Wt28G8D/q9oUAXoeyzPjJAJao+4cB2KL+HapuD8133TzqPQbA8er2IACbAMyMe93V8g9Ut4sBLFHr8yyAr6j7/wjgu+r2DQD+qG5/BcDf1O2ZavsvBTBJfS+S+a6fj/r/GMBTAF5Vf8e+3gDqAFRZ9mW9nee94j5uzCkA3hR+3wLglnyXK4R6TbQI9I0AxqjbYwBsVLcfBHClNR2AKwE8KOw3pYvCfwAvATi3P9UdwAAAKwCcBGV2YJG6X2/nAN4EcIq6XaSmY9a2L6Yr1P8AxgF4F8BZAF5V69Ef6i0T6Flv51EwuRwGYLvwe4e6L26M4pzvUrd3AxilbjvVP9L3Re1Oz4airca+7qrZYSWARgBvQ9EyWzjnvWoSsQ56/dTjBwAMRwTrDeAeADcBSKu/h6N/1JsDeIsxtpwxdr26L+vtPNCaokRu4Jxzxlhs/UkZYwMB/APADznnrYwx/Vhc6845TwE4jjE2BMALAKbnt0TZhzE2F0Aj53y5usB8f+JUznkDY2wkgLcZYxvEg9lq51HQ0BsAjBd+j1P3xY09jLExAKD+bVT3O9U/kveFMVYMRZj/lXP+vLq7X9QdADjnLQDeg2JqGMIY05QqsQ56/dTjgwE0IXr1ngPgYsZYHYBnoJhd7kX86w3OeYP6txHKB/xE5KCdR0GgfwxgqjoyXgJlsOTlPJcpG7wMQBvFvgaKfVnbf7U6En4ygANqt+1NAOcxxoaqo+XnqfsKFqao4o8AWM85v0s4FOu6M8ZGqJo5GGPlUMYN1kMR7F9Sk1nrrd2PLwGYxxUj6ssAvqJ6g0wCMBXA0pxUIgM457dwzsdxzidCeW/ncc6/ipjXmzFWwRgbpG1DaZ9rkYt2nu/BA58DDBdC8YjYDOAX+S5PCPV5GsAuAD1Q7GLXQbEVvgugBsA7AIapaRmAP6h1XwOgWsjnWgC16v9v5rtePup9KhTb4moAK9X/F8a97gCOAfCJWu+1AH6l7p8MRTDVAngOQKm6v0z9Xasenyzk9Qv1fmwEcEG+6xbgHpwBw8sl1vVW67dK/b9Ok1m5aOc09Z8gCCImRMHkQhAEQfiABDpBEERMIIFOEAQRE0igEwRBxAQS6ARBEDGBBDpBEERMIIFOEAQRE/4/U2+hpOuatLAAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Samples\n" + ] + }, + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoLklEQVR4nO2df5Ad1XXnv4eRkKZYWFnFT400HkUIUWBhYU8Qqkl2MUGRrElFCtkYqUTiHykp2ZXXJmzhjCytMWumNGU5xEmljEtKnAqFzI/EZGAzClhEZrNWSaJGliwBkcKoJAMDi3CpMBgLYQ1n/3jdQ8+bvv3z9s/3/VRNzXu33+u+/fq9b5977rnniKqCEEJIPTmv6A4QQgjJDoo8IYTUGIo8IYTUGIo8IYTUGIo8IYTUmClFd8DLxRdfrF1dXUV3gxBCKsWBAwd+qqqX+G0rlch3dXVheHi46G4QQkilEJGfmLbRXUMIITWGIk8IITWGIk8IITWGIk8IITWGIk8IITWmVNE1pLXp6hua1HZyoLeAnhBSH2jJk1LgJ/BB7YSQaFDkCSGkxlDkCSGkxqQSeRHZKiJHReSwiPyjiMzwbNsoIiMickxElqXuKSGEkNikteR3AfiIql4H4N8BbAQAEbkGwGoA1wJYDuBbItKW8liEEEJikkrkVfX7qnrOeboPwGzn8UoAD6vqWVU9AWAEwA1pjkXqjSmKhtE1hKTDZgjl5wA84jzuQEP0XV5x2iYhIusBrAeAzs5Oi90hVYOCToh9QkVeRJ4GcLnPpk2q+rjzmk0AzgHYEbcDqroNwDYA6O7uZlVxQgixSKjIq+otQdtF5DMAfgvAb6iqK9KjAOZ4XjbbaSOEEJIjaaNrlgP4EoDfVtVfeDY9AWC1iEwTkbkA5gN4Ns2xCCGExCetT/6vAEwDsEtEAGCfqv6xqj4vIo8CeAENN84GVR1LeSxSAdZu34s9x0+PP++ZNxM71i0psEeEtDbygYeleLq7u5WVoapLs8C7UOgJyRYROaCq3X7buOKVWMNP4IPaCSHZwyyUZALMBElIvaAlT8ZhJkhC6gdFnlijZ97MWO2EkOyhyBNr7Fi3ZJKgc9KVkGKhT55YpSyCfuXGIZzzBI5NEWBkC+cWSOtBS54AAAYP1mdBcrPAA8A5bbQT0mrQkrfA4MFRbH3qGF598wxmzWjHXcsWYNX1vvnYSsngwVFsfOyIcXvVomuaBT6snZA6Q5FPiSuQZ37ZWNA7+uaZccGsitBvferYeP+9dMxox56+mwvoESHEFhT5lPgJ5JlfjmHrU8cqI/KvvnkmVnsQVR/VEFI3WlrkbSz8sSmQRTFrRjtGffo7a0Z7rP2UZVQzRfxdM1Mkty4QUhpaduLV1sIfkxA2t6/dvhddfUPjf2u37411nCy5a9kCtE+dWJ2xfWob7lq2INZ+TKOaOx45lLaLsRjZ0usr6OcUuHrTzlz7QkjRtKzI2yKKQPol7tpz/HRhQj94cBQ9A7sxt28IPQO7AQBbbl2IjhntEDR88VtuXRjb+g4aveS9anZkSy+mt01W+nfHlEJPWoqWdtfYwBXCID90mRJ3mVwqW25dmGqSdW7fEMoWvPLumH+PTO2E1BGKvAVWXd9RmcnFIJdK3z/8GEf7V8TeZxkFnhDSgCLfYgS5VFxXRlyhz0Lgl973DF489c748/mXXoBdd96UwZEIqTct65M3RdFksfCnTIm7wiJmyuDKaBZ4AHjx1DtYet8zsfbj55MPaiekjrS0JZ/XSs4d65ZEKosXJ6Tz6k07Jwjy9DaJZIHftWzBBJ98XsT5rJsFPqzdxNH+FYk/J0LqQkuLfJ6EJe4KCulsFshm4QKiu1rcuQObYY0Cf5eNADhRcEoECjppdVK5a0TkayJyWEQOicj3RWSW0y4i8pciMuJs/5id7hIgfdTIqus7rLoyTgz0ovldSQWeBUoIsUtaS36rqv5PABCRLwD4CoA/BvBJAPOdv8UA7nf+k5Jg25WR1mKPIu7zL70g1TEIaUVSibyqvuV5egE+GLWvBPCAqiqAfSIyQ0SuUNXX0hyP2KUsroyo1jujawiJT2qfvIj0A/gDAD8D8AmnuQPAy56XveK0TRJ5EVkPYD0AdHZ2pu1OSzC9TXxdM3FdLSZxrVpqYUKIGWkY2wEvEHkawOU+mzap6uOe120EMF1V7xaRfwIwoKo/dLb9C4A/VdXhoGN1d3fr8HDgS2pNHtE1QceKclzTe9PcGKJa8lnefGyfEyF5IiIHVLXbb1uoJa+qt0Q8zg4AOwHcDWAUwBzPttlOGwkgjqgU5WqJEwVkk7wF3m2n0JOqk8pdIyLzVfVF5+lKAEedx08A+LyIPIzGhOvP6I8vB26+96qQ54iCkDqS1ic/ICILALwP4CdoRNYADYt+BYARAL8A8NmUxyEWaE5OVhZODvTGFuw41jdvBqSVSRtd87uGdgWwIc2+iX1MZf6yoLmY9hRppP81kZXo5u2KYWUsUja44rWiJBGTqNWqkopf0ATqOW0If5DQB+2nrJZ30DlXsd4vqR8U+QqStMyeqcyfqWC3X74dk2slCn4l+ZqxVbErDkncRUC0PlWt3i+pHxT5CpK0eLhfcjJTmb+galZRfN5lICxr5dL7nhlfYJXlSMEdQUVJUkeIbSjyFSRp8fAoVaxcolaziltKz70hhPnok+IV67CslW764qxX0s6a0W68aTbfIJklk9iGIl9BTG6XsFzxgN0qVn7ZMKPS7KNPOhpIa4HHTV8cF3ekFDXrZ9LCLVHYPHgED+1/GWOqaBPBmsVzcO+qhdaPQ8pFyxYNqTJRiofnQdoCI66P3qa7p6tvaPyvSNIURM+icMvmwSN4cN9LGHNWuI+p4sF9L2Hz4BHrxyLlgpZ8BYnjdklKz7yZvi6buNWsLrvwfLz+9nu2uhVIUmFPGs2TdMK2CB7a/7KxndZ8vaHIV5Ssi4ebBD7uRGGYwBdtcZuIGkcfVdBNN828GDPkqDK1k/pAdw2ZhEl4m0Wq6FqpXoENczucHOgtNB/9jnVLIo+Csvhc28R/n6Z2Uh9oybcgtlwMfoVHTEyRaHHyUfCrOmVyR7h09Q0VXnTEOwoy3Uiziq5Zs3gOHtz3km87qTcUeUvM7RuaUOe0DPVN/bC9zP/jXR+K5IbwhkumddF4o4ji7CtuJE1WqQ+C+pxV+KTrd2d0TetBkbdAs8ADjRJZc/uGSin0Noki8JddeH6sfbrC6pdQzRtFlIc/P+90w5sHj2QmvPeuWkhRb0Eo8hYweSHKNKVV1ATnZReej/2blk5ou2haG946G54oLY8oorLx4L6X0P3hmbU+R5IvFPkWIK7AB+Wn8baHWbim7YfvWR65T6YoojrHdzOpGbEJRZ744gp0kBh39Q1Zi6ePw3V3PxlpJGAibu6dvF02TGpGbMIQSguYgtCqFJzmJ2JRrG1TaKCbl8X26tO0Al8GotwwoqSFHjw4ip6B3ZjbN4Segd0YPMgKm2QyFHkLnBjonSToZY2u8SOtlbpj3RKcHOiNVcnJ9NqwvkQR+Ljnc9G0tvAXWebkQC9uv7HTuD0sD5E7KT365hkoPkg3TaEnzdBdY4mqCHpU8pioDRPjpfc9MyHsMUqcuyucU84TnHt/8tT3lPMmj68O37O8kBHCvasWovvDMyOnf/aSNN00aT0qL/JViU8vkjiFPrIW96gLsZoFHogW5/7gvpdw76qFvgIvAL7xex/1fV/YZHDPwG5jZE+axWVJI4iSppsmrYcVkReR/wHgGwAuUdWfiogA+As0inn/AsBnVPVHNo7lpZXj07MgSRSOjf3bvrGY9qdIHrFiqr5lY3FZkjxEadJNN8MUxPUmtciLyBwAvwnAu2b6kwDmO3+LAdzv/LdKFeLTy4LNjIlFZ1mMGmcfhzxL+dm4DnGqfAXhpiB2cVMQA6DQ1wQbE69/DuBLmKitKwE8oA32AZghIldYOBZJgTs56v1Lso80221w+J7lmR/DRFp3iK3RzKrrO7Dl1oXomNGeKnd9UApiUg9SWfIishLAqKr+WCZms+sA4P2WvOK0veazj/UA1gNAZ6c52oDYx8/vbYMoMfZJKTrJWBJ3SFbYSDddpxTEfhP1WZd2rAKhlryIPC0iz/n8rQTwZQBfSdMBVd2mqt2q2n3JJZfEem8d4tOLIonA5+2m8RN0t89FuIyKqL6VNXVJQWyaqA8r5t4KhFryqnqLX7uILAQwF4Brxc8G8CMRuQHAKABvDtPZTptVTgz01iq6Js8qQ3EEvigfvKmP7qSmt19RRg1BcelBCGCMemkT8bV6yyySUT6rqqUgNn1Xsq7hWwUSu2tU9QiAS93nInISQLcTXfMEgM+LyMNoTLj+TFUnuWpsUFVBb8Z2CuA02DpelNBNU575OGGfYaSNGAn6jlXN3RH2mTK6pn5kFSe/E43wyRE0Qig/m9FxSMU5p3ZuKgL/qCoBcHyL/RztYe4ukyVf5rqwZegDsY81kVfVLs9jBbDB1r6JfeZfeoGvSBU9sZmUPF13UeYzgtwdeYmp3+dRV+r2fbZJ5Ve8tgJZWH677rwpcjRCWS3PZk40WcmKeAW5o55nmAVfBneHaaFgXYnzfW41REvkO+zu7tbh4eGiu1EISfzPeQhtUL+iHD/tecW5waTta1TyOk4ayvp9ItkgIgdUtdtvGy35kmBzojEOWVvpSc7La31TePKBn3N9ociXiKQ53ZOSV0RPEsu7zFTd/0tBby0o8sSXKALc1TfUUn7PoM8k6HMwjZbCRlFpEocFRRuR1oIiT1Lx4ql3xsUqLwuxiIngpH74uLlqTO1xE4fVbaEgSQ5FvuSYLL72qW0YPDgaKXdJXqKYx8KtMBdTVSKBkvLQ/pdDRT4o1XKcQuykHlDkK0CHT+7wqGlv88rhnoS4onzd3U8G7s/GTabsN4mwlbRxrmtX3xBuv7GTueRrDkW+AlSpCtDi/l3Yv2nphLYg4YwqoLbL8w0eHJ1UjemORw75vjZMOKPWr7WB7Zw4cXLJ0/1TTVjIuwKY0tvmnfY2SvTI62+/N+G5rZGEbYH3K4JtiyxHSaaVtIMHR9EzsNvKMfxyyQdVYSPlhpZ8BbBVBSgNrnWaVQ76LIgjts1FsZtJuo4hyCUVZX9ulssgV4p70wo7h6j4uYRYha26UOQrQNJiz2kwuRzcMMEy+PSDyKJ/zZ9J1DBTU/tlF54/7tpKMxew9alj1gQeKHeaZBIfinxFSFoFKKuVtJddeP4k14zbngVR6rqWaYI0Cq+//d74HEaaviedm7n9xs4JPnmX6VPPw9y+oVyMCZI9FHkyiSkRDLn9m5b63jyaJ11tcfie5ZMmXy+a1lZovdcgot5Y/W6UcZnlE30FNKKy9vTdHPp+N7pGAJx3nuCd9xqfsXeuosqLqxb375rwOXtHUK0AE5TVnKwSVV29aSfeHZv83ZneJjjaPzF/e9kWL5nwC1UFki12ikPSRG/u+/x88u1T22IX9u4Z2B14s6hidE2zwLvUTeiZoIykxiTqzfi9plnE1m7fO0G0eubNxI51S0L3nfXN4tU3z8TeX9wShEkIWwBma84mLFS37ILuh2mkZGMEVRUo8hXG78ffMaN9wg/dBlEFPqxvJwd6sXb7Xuw5fnpC+57jp7F2+95AoY+bTC3JXESckFRbgm5rDiPpnI0Xk9sn71BdYheKfEkJs1pNIuP+SJPGfvsJpg2BD2oHMEn4beCeS1RB9hO4rKOIXn/7PXT1DUUezWRJGUJ1iX0o8iXEVgrgKGF1tn3j09sk9k2hTATlwMmSPcdPF55XpohQ3azJOwqsjKQSeRH5KoB1AN5wmr6sqjudbRsB/CGAMQBfUNWn0hyrrJS95JgpKsLFdlKxo/0rEgtknukB4vSjqD4UJfRVFvVm9m9a2vLRNTbSGvy5qi5y/lyBvwbAagDXAlgO4Fsi0mbhWKXCb/Xni6fewdL7nimmQz7MmtGeWiymtxUTKGdLbKsWP9+Mqf9VP6+8cNchuH+tJPBAdu6alQAeVtWzAE6IyAiAGwDszeh4hWBa3l+mZf82/KlH+1ckmny1Sdo0wnH982WDgk6SYsOS/7yIHBaR74jIh5y2DgDeLEevOG3EElF/9LaG3kf7V8QWVJuUOWWyH6YFZScHetEzb2a+nSEtTaglLyJPA7jcZ9MmAPcD+Boabt+vAfgzAJ+L0wERWQ9gPQB0dnbGeWttiWq15m2dRj1O3sIbdcLSVr+mCHDOZ1AzRYCRLY1r17zd268d65b4hpJmRdKqVqQehIq8qt4SZUcish3APzlPRwF4c6LOdtr89r8NwDagseI1yrHKQpYFneP8+KJEEKR1d0QRyCKiUpoxlSK0uTLVtK9zGj0yyhsuGfW6JLl+UXLhU+jrTaq0BiJyhaq+5jz+EwCLVXW1iFwL4Lto+OFnAfgXAPNVNTCmr4ppDcoSXZNVBEGRop12NWkeq1HjkEZMo1jjcQqhpOlX8/yMXyoLki9ZpjX4uogsQsNdcxLAHwGAqj4vIo8CeAHAOQAbwgS+qpQlXDKLiIGihNG2ZVkGgc+a5vw1tguhuPhNwL87prh6004KfUlJJfKq+vsB2/oB9KfZP2k9ioiCsXVMk68+D/xyytvMMe9iirCq8gK4usMVr6Rw4vjPTflvisY76VoEZaz366XVFyQVCUW+ZoT5S4tI++sew9axbQp8mnOPkxgt68/YlFwsjDyuvV+6X2/BFJItFPkaEeYvtZUTJyl+ApjVsYuMGLF57LCRgXssv+RiWfTNlJsoaFU00/0Wi43FUKQk2PaX5jUB2tU3ZDWFgWk9gd9ro7zORM/Abgwe9I0MtkKY28rb31XXd2DLrQvRkXFa4KP9KyYJOqNryg0teRKIn/CZyvClFeqoVn3PvJm+LpuglaRZpD/wRrCUIamXm1wsqO82Rk4U9GpBkSexMdVVzcMf7YqUqfhImtWdfgnnwjjzyzFsferYJJEvau4jCqbPKKv+Md1vsVDkC+DKjROXvbuRGWlJ4i+1jQ2hMImClySTr0FWbBKBd2mObEky92EzKifpyuOoVn7c7y/T/RYLffI50/wDARqx1VduTP8jD/OXViVlbdiPP40gLu7f5dseJPCu/9vk71akm1fIIuwyq2ua9Pvb6ul+i4SWfM6YFsvYWkQT5i8tOuqk6NWnaSI6okSw5HF+RV7DrL+/aZnbNzShSI6gmgXIbUKRbyH8cpvkPWEYNLFZhHjFEWVvebwkMelAw+KN45or2yirzDQLPNAYZc3tG2ppoW85kc9L6MpmUQTlNikiMiRMvL552yLf5Fqm9iREEfjmjKJRIliCiGvxpomGKcPIKU9fvOmjLckgozBaSuTzErogiyIoF3mWmHKb+EWG2CSpxb71qWOx2qMSJ6KjbPV6w/AT1DhEuS5xvr9c6VoOWmriNUjobBJkUYxs6Z30g7AVXROEKbdJljlPklRzcicwTe6QV988k9iyjWNFnhzoNQq8bevYhkvGJKgmLrvw/AkToVH7EOf7y5Wu5aClLPkihM6PrAXdD1Nuk1kZr5CMQxTx1Iiv8+P1t99DV98QeubNnFC0I0+yGrGFCbrfDcCvolYUd2YR398oCPwNrGLK0JeHlrLkTYJWJqHLiruWLUD71LYJbe1T26wU+q4ae46fxtrt+deUz2rEdvWmnYHbveGLJrr6hsbdmaNvnoHiA3dmlqkbbHJioHeSoBc9F1YGWsqS9wuBy0LoymhReCNDTFZaVP/55sEjeGj/yxhTRZsI1iyeg3tXLYzVn3kbdyZ6XxgXTWubkHLBxJ7jp3PPGBkk8Gn6YiuXu+15myiL2mzT6oLuR0uJfBShs8GJgd7SRdcAH0SG+BF1lebmwSN4cN9L48/HVMefxxHspO8Lok0Eh+9ZPmmCPYgyhShm0Zc4k6+23Zn7Ny0tPPMpaTGRB4KFziZFCnqW1ulD+182tjeLdZQQPr/3JWXN4kbteD+LtGiCPn8bIyMTXt972HegCvM2JD4t5ZNvBZJEtMRhzFD43dQe5gv2vi/pjahNBLff2DkujFEsz6CMlbaIErnijozcz8Ed4WwejF6fNWpuorDvAOdt6knLWfIkHW0iRkFP8r42mShQcRJ4mV4bViWpyOiaZuKMjEwc7V/hWzAmLnc8cgjfvG1R4aui/WjOOto87+WmuyaTSS3yIvLfAWwAMAZgSFW/5LRvBPCHTvsXVPWptMcixbNm8ZwJPnkva7fvNYqn6X2uiyWMOFa+aYJ9y60LSyFYXoJGRu6NLYqAhZV4jEpe7sw4+KWVbv7U3jo7huvufpJC70MqkReRTwBYCeCjqnpWRC512q8BsBrAtQBmAXhaRK5S1XI5SklsQbh31UKjyAel/3Wt0qx8z0HnMW3KeaUU+Ki8dXZs0vlVpRqTjQimqGmlo0RVtSJpLfn/CmBAVc8CgKqectpXAnjYaT8hIiMAbgCQf3AyMZJ3XpN7Vy20HjIJhJ/H2XPv4++HX5og8mEFz4OIMqFsEjG/lalJ8NbuLQNBOeYZRVMsaSderwLw6yKyX0T+j4j8qtPeAcDrbHzFaZuEiKwXkWERGX7jjTdSdodUJWd83nitwaCC51EJSwngdxOwJfAuJh980Hcg6vdg8OAoegZ2Y27fUGgt2yxrJJD0hFryIvI0gMt9Nm1y3j8TwI0AfhXAoyLyK3E6oKrbAGwDgO7u7lZPGOdLc9WisMRZNgXdL3IjSY3VMmGz4Hmceqp5LgwK+g6EuVDiJvLLOse86fvWzEXT2nzb4/5+6kaoJa+qt6jqR3z+HkfDQn9MGzwL4H0AFwMYBeCdUZvttJGY+JWle/HUO1h63zOZH9vkwtixbskkQS9TxEpeRHF3FZ3q14TfKMRNDnfHI4dySeQXFb/vW7PpYZqcLvL3UxbS+uQHAXwCwA9E5CoA5wP4KYAnAHxXRO5DY+J1PoBnUx6rJTGVpUtajzQqYaOBMgh6HAEtcpSR1epOG7V743yGeSfy85L0+1bU76dMpBX57wD4jog8B+A9AJ9WVQXwvIg8CuAFAOcAbGBkTfnIO3eLTeJayF6RMBU8d/drO/InCkEVs/wwjbKyvJ6mla9F1Ugg0Ugl8qr6HoDbDdv6AfSn2T/JnjIL+nV3PzkhLC7Jghc/MQxbPJRFXh0g+oIwE2HXKss8MUErX0e29AZG15Bi4YrXkjP/0gt8h5bNZelcmheOVNVX3izwQPQFLycHescjWd4daywqai4Y4hX+eRt3+oqvzbw6QPoFYWlwzzHOKKVjRnvoytcyJuLzEvf3U0eYu6bk7LrzpklfSFN0gN/KwKJypwex9L5nxif5uvqGfCfBTAtb/BYGNRNUds6PuPl4XOJax/euWojbb+wct9ybc+4E7TOtJZ4kN86evptxYqAXe/pujiTwwAdlLstCnN9PXaElXwGifiFNYWZRVwzmQVC0g60fXtyyc2ncKM3iG+YTj7IgLA8XmjtKSTMvU5XC2a0k6H5Q5EmuZB3tEGU1ajM23ShlnuPwYiP7J6kGFHlSSqJWeLIhUFnn1XHJOpIpzg0u6mSvTaoayVV1KPI1ouorUb0cvme57+RrFExl54KqJGWVV8elbBWSbEz2BpW59HPL+cEKUdnDidcaUYWVqKaoBr/2sCga02Kg/ZuWThL05uiaVkWASZO9JryT4+6fF1Ph7CsNES2kGGjJ14wyCbofu+68yVoukaAMjBT0yXzztkWR0y1HHXn4hUuWNZVDlmRZwjEtFHmSO3EEvcqrcoH8Ba+Kn1fV+tuMreL2WUGRJ5liw8Kp0g/eS1EWbV6fV5biXCVfvY0SjllCkSeZUXYLpwxURciaCXLnmFaZ2mDw4GjpatAmXUyXFxR5khlBFs4TB0dT56WpMlUV9yjEmXeJM9qJm+c+L9LmJMoaijzJjCALJ2leGpIPaX37Wawy3frUMWOe+yJFvsicRFGgyJPMMFk4Jlq1EHMZXRBA+UYbpnz2Rea5B/JbTJcUijzJDJOFUwbyiOiIYg37uSDueOQQ7njkUOR+lTl8L4w4I4ZZM9ox6iPopjz3eRJ1MV0RWWJFSzI5ADRqvA4PDxfdjcS0ei1JP/wEKEj487Aeg/zAeVuvPQO7fYWrGVO/mie3XaIueEpD3qGPzTdEoJHnfsutC0sx8gnDL0ssYEfoReSAqnb7baMlb4k8sitWET8Lp3nS1cUtxNxKN8u0roYiw/fyviG6Ql5G11YUisoSS5G3BGtJRscvL40bXdNqN0uTCyIqRYXvFbWAadX1HZUR9bJAkSeFYIqiabWb5V3LFkxyQcTBZvhe1CpPZUu2RoJJlaBMRB4RkUPO30kROeTZtlFERkTkmIgsS91TQmrIqus7sOXWhehIOHloCtOLG75XhSpPVceUDTbrLLGpRF5Vb1PVRaq6CMD3ADwGACJyDYDVAK4FsBzAt0SkLWVfS02c7IqkWLIqsZeUVdd3YE/fzYn6FaWkYBSqUuWpyhSVJdaKu0ZEBMCnANzsNK0E8LCqngVwQkRGANwAoFzFRi1iM7tiK5NX4eWyuhWi9qvqSb3iErYytiorpovIEmvLJ//rAF5X1Red5x0A9nm2v+K0TUJE1gNYDwCdnZ2WulMMFPT02LxZlnWRUVpazSceJfUBV0ybCRV5EXkawOU+mzap6uPO4zUAHkrSAVXdBmAb0IiTT7IPUi9s3CzLmufERBGLZLwEVXlqpqzpjFt1xXQYoSKvqrcEbReRKQBuBfBxT/MoAO/Mz2ynjZBciJvnpEjR8lsks+f4aazdvjc3oT8x0Bs5ugYoXtBJdGyU/7sFwFFVfcXT9gSA1SIyTUTmApgP4FkLxyIkEnHynAS5P/KgqEUyXgYPjmLWjHYIgI4Z7fjmbYuMAk+qhQ2f/Go0uWpU9XkReRTACwDOAdigqhxLkdwoc56TspGHayvOKCEp7oppMpHUlryqfkZVv+3T3q+q81R1gar+c9rjEBKHu5YtQPvUiT/69qltuGvZgoJ6ZA/bIaBBri0bpI3Bj3JeVYmuKQKueCW1pEp5TnrmzTQmrjJh0yeedQpfGzH4nANIDkWe1Jaq5DnZsW5JodE1dG3VG4o8aXnKEBKYpaBfuXEI5zxm8xQBRrZ8cG5++XPq4toiFHlCANTXHdAs8ABwThvtrtBn7dqKE4NP7EORJ6TGNAu8qT1L11bcGHxiF4o8ISRzKOjFYWMxFCGEkJJCkSekxkwxOL5N7aR+0F1DSESu3rQT74594Fme3iY42r+iwB6FM7KlNzS6htQb0YxrQcahu7tbh4eHi+4GIZNoFniXKgg9qT8ickBVu/220ZInJAJ+Ah/UXjeKToVMkkORJ4QEYisV8ubBI3ho/8sYU0WbCNYsnhO7TCGJD0WeEBKIjVTImweP4MF9L40/H1Mdf06hzxZG1xASgelt/uEopnYykYf2vxyrndiDIk9IBI72r5gk6Jx0jc6YIcDD1E7sQXcNIRFpVUFPkgq5mTYRX0FvE46EsoYiT0gFKDL3i41UyGsWz5ngk/e2k2xhnDwhJcevshJQriRfUVI1M7omOxgnT0iFsVFZKUuCCqF7hf7eVQsp6gWQauJVRBaJyD4ROSQiwyJyg9MuIvKXIjIiIodF5GN2uksIISQOaaNrvg7gHlVdBOArznMA+CSA+c7fegD3pzwOIYSQBKQVeQVwkfP4PwJ41Xm8EsAD2mAfgBkickXKYxHSkpjiTxiXQqKQVuTvALBVRF4G8A0AG532DgDeVQ6vOG2TEJH1jqtn+I033kjZHULqx4mB3kmCXqZJV1JuQideReRpAJf7bNoE4DcA/Imqfk9EPgXgbwDcEqcDqroNwDagEV0T572EtAplFvQyFEInZkJFXlWNoi0iDwD4ovP07wH8tfN4FIA3AHa200YIqSEU9PKS1l3zKoD/7Dy+GcCLzuMnAPyBE2VzI4CfqeprKY9FCCEkJmnj5NcB+AsRmQLgXTQiaQBgJ4AVAEYA/ALAZ1MehxBCSAJSibyq/hDAx33aFcCGNPsmhBCSHmahJISQGkORJ4SQGkORJ4SQGkORJ4SQGkORJ4SQGkORJ4SQGkORJ4SQGkORJ4SQGkORJ4SQGkORJ4SQGkORJ4SQGkORJ4SQGkORJ4SQGkORJ4SQGkORJ4SQGkORJ4SQGpO2MhQhhJAUrN2+F3uOnx5/3jNvJnasW2Jt/7TkCSGkIJoFHgD2HD+Ntdv3WjsGRZ4QQgqiWeDD2pNAkSeEkBqTSuRF5KMisldEjojI/xaRizzbNorIiIgcE5Fl6btKCCEkLmkt+b8G0KeqCwH8I4C7AEBErgGwGsC1AJYD+JaItKU8FiGE1IqeeTNjtSchrchfBeBfnce7APyu83glgIdV9ayqngAwAuCGlMcihJBasWPdkkmCbju6Jm0I5fNoCPoggN8DMMdp7wCwz/O6V5y2SYjIegDrAaCzszNldwghpFrYFHQ/Qi15EXlaRJ7z+VsJ4HMA/puIHABwIYD34nZAVbepareqdl9yySXxz4AQQoiRUEteVW8JeclvAoCIXAWg12kbxQdWPQDMdtoIIYTkSNromkud/+cB2Azg286mJwCsFpFpIjIXwHwAz6Y5FiGEkPiknXhdIyL/DuAogFcB/C0AqOrzAB4F8AKAJwFsUNWxlMcihBASE1HVovswTnd3tw4PDxfdDUIIqRQickBVu323lUnkReQNAD+xsKuLAfzUwn7KSF3PjedVPep6blU8rw+rqm/kSqlE3hYiMmy6q1Wdup4bz6t61PXc6nZezF1DCCE1hiJPCCE1pq4iv63oDmRIXc+N51U96nputTqvWvrkCSGENKirJU8IIQQUeUIIqTW1E3kROekUMTkkIpVdWSUi3xGRUyLynKdtpojsEpEXnf8fKrKPSTGc21dFZNS5bodEZEWRfUyCiMwRkR+IyAsi8ryIfNFpr/R1CzivOlyz6SLyrIj82Dm3e5z2uSKy3yl89IiInF90X5NSO5+8iJwE0K2qVVvMMAER+U8Afg7gAVX9iNP2dQCnVXVARPoAfEhV/7TIfibBcG5fBfBzVf1GkX1Lg4hcAeAKVf2RiFwI4ACAVQA+gwpft4Dz+hSqf80EwAWq+nMRmQrghwC+COBOAI+p6sMi8m0AP1bV+4vsa1JqZ8nXBVX9VwDN1XxXAvg75/HfofFDqxyGc6s8qvqaqv7Iefw2gH9Do45Cpa9bwHlVHm3wc+fpVOdPAdwM4B+c9spdMy91FHkF8H0ROeAUJKkTl6nqa87j/wfgsiI7kwGfF5HDjjunUi6NZkSkC8D1APajRtet6byAGlwzEWkTkUMATqFR4e44gDdV9ZzzEmPRoypQR5H/NVX9GIBPAtjguAZqhzb8bHXytd0PYB6ARQBeA/BnhfYmBSLyHwB8D8AdqvqWd1uVr5vPedXimqnqmKouQqPuxQ0Ari62R3apncir6qjz/xQaxcXrVFv2dcc/6vpJTxXcH2uo6uvOj+19ANtR0evm+HW/B2CHqj7mNFf+uvmdV12umYuqvgngBwCWAJghIm5RpUoXPaqVyIvIBc7EEETkAjSqVj0X/K5K8QSATzuPPw3g8QL7YhVXBB1+BxW8bs4k3t8A+DdVvc+zqdLXzXReNblml4jIDOdxO4ClaMw5/ADAf3FeVrlr5qVW0TUi8itoWO9Ao7Thd1W1v8AuJUZEHgJwExppT18HcDcaBdMfBdCJRkrmT6lq5SYwDed2ExrDfgVwEsAfefzYlUBEfg3A/wVwBMD7TvOX0fBfV/a6BZzXGlT/ml2HxsRqGxpG76Oq+r8cLXkYwEwABwHcrqpni+tpcmol8oQQQiZSK3cNIYSQiVDkCSGkxlDkCSGkxlDkCSGkxlDkCSGkxlDkCSGkxlDkCSGkxvx/Gy+NPMSo82IAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "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 have to write about 5 lines of code to compute the lower and upper bounds." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 25, + "source": [ + "# YOU CODE HERE\n", + "slope = -samps[:,0] / samps[:, 1]\n", + "slope_mean = np.mean(slope)\n", + "slope_std = np.std(slope, axis = 0)\n", + "lower = slope_mean - 2*slope_std# YOUR CODE HERE\n", + "upper = slope_mean + 2*slope_std# YOUR CODE HERE\n", + "\n", + "print(f'The confidence interval is [{lower}, {upper}]')" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "The confidence interval is [0.19559960296585568, 0.520658828359935]\n" + ] + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "### Try tinkering with the $\\sigma$ value in the scripts above\n", + "Then answer the following..." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "#### What happens is sigma is too big?\n", + "\n", + "We will have more number of rejections, hence, slow mixing" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "#### What happens is sigma is too small?\n", + "\n", + "Short movements, thus, slow mixing" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "#### Why did you remove the burn-in samples? \n", + "\n", + "It may take many iterations for the algorithm to search the parameter space for likely values. The starting values can be complete garbage." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "# Plotting the classification line" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 37, + "source": [ + "plt.scatter(X[ind,0], X[ind,1], color='blue')\n", + "plt.scatter(X[~ind,0], X[~ind,1], color='red')\n", + "x = np.linspace(-0.3, 0.3, 100)\n", + "y = slope_mean * x\n", + "plt.plot(x, y)\n", + "plt.show()\n" + ], + "outputs": [ + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAulElEQVR4nO3deXxU1dnA8d8hhH0PCAgkAdlkUyFA3dBWxR1cqFVR0WpR+9rW+tYGBK2iWLUtaqtVwQ0Lvi6AiAFZVRZBBBTDImsI+xq2ACEkmef9Y4YSw0wyk7kz9965z/fzmQ+z3Jl7zmQ4z71nea4REZRSSnlXFbsLoJRSyl4aCJRSyuM0ECillMdpIFBKKY/TQKCUUh5X1e4CVEbjxo0lPT3d7mIopZSrLFu2bJ+INCn7vCsDQXp6OkuXLrW7GEop5SrGmM3BnteuIaWU8jgNBEop5XEaCJRSyuM0ECillMdpIFBKKY/TQKCUUh6ngUAppTxOA4FSSrnA2l35vDB9DbG4dIArF5QppZRXnCj28e+vNvDqlxuoWyOZO36WxpkNalq6Dw0ESinlUMu3HiRzQjZrd+fT/9wzeeK6TqTUqW75fjQQKKWUwxScKGHUrLW8tWATZ9StwVuDMrjs7KYx258GAqWUcpCFG/cxdNIKNucd4/beqQy5uiP1aiTHdJ+WDBYbY64yxqw1xmwwxgwJ8nofY8x3xphiY8yAMq8NMsasD9wGWVEepZRym8PHixg6KZvbxyzGAB8M/hnP3tg15kEALDgjMMYkAa8CVwDbgCXGmCkisrrUZluAu4E/lXlvI+AvQAYgwLLAew9EWy6llHKL2at3M2zyCvbmF3J/nzY8fHl7alZLitv+rega6gVsEJEcAGPMB0B/4L+BQERyA6/5yrz3SmCWiOwPvD4LuAr4PwvKpZRSjrbvSCFPfbaaz37YQcdmdRlzVwbdWjaIezmsCAQtgK2lHm8Dekfx3hbBNjTGDAYGA6SmpkZeSqWUcggR4dPlO3jqs1UcKSzmkSva88AlZ1Gtqj1Lu1wzWCwio4HRABkZGdavqFBKqTjYcbCA4ZNX8sWaPZzbqgEvDOhG+6Z1bS2TFYFgO9Cq1OOWgefCfe+lZd77lQVlUkopR/H5hPe/3cJzn6+hxCc8fl0n7r4gnaQqxu6iWRIIlgDtjDGt8TfstwK3h/neGcCzxpiGgcd9gaEWlEkppRxj076jDJmYzeJN+7mwbQp/vbEbqSm17C7Wf0UdCESk2BjzEP5GPQl4W0RWGWNGAEtFZIoxpifwCdAQuN4Y85SIdBaR/caYp/EHE4ARJweOlVLK7YpLfLy1YBOjZq2jWtUqPH9zV27JaIUx9p8FlGZikcAo1jIyMkQvXq+UcrLVOw6TOTGbFdsPcUWnpjxzQxea1qtha5mMMctEJKPs864ZLFZKKTcoLC7hlS828NpXG2lQK5lXb+/ONV2bOe4soDQNBEopZZHvthzgzxOy2bDnCDed14LHr+tEw9rV7C5WhTQQKKVUlI6dKObvM9bxzsJNNK9Xg3fu6cnPO5xhd7HCpoFAKaWi8PWGfQyZlM3W/QXc+bM0Mq/uSJ3q7mpa3VVapZRyiEMFRTw79Uc+XLqV1o1r89H959OrdSO7i1UpGgiUUipCM1bt4vHJK8k7eoIHLjmLhy9vR43k+CWJs5oGAqWUCtPe/EKenLKKqSt2cnbzerw1qCddW9a3u1hR00CglFIVEBE++X47I7JWc6ywhEev7MDgPm1ITrInSZzVNBAopVQ5th8s4LFJK5i7bi890hry/M3daHtGHbuLZanECGdKKeuMHw/p6VCliv/f8ePtLpEtfD7hvUW59B01lyW5+3ny+k58dP/5CRcEQM8IlFKljR8PgwfDsWP+x5s3+x8DDBxoX7niLGfvETInZrMk9wAXt2vMszd2pVUj5ySJs5rmGlJKnZKe7m/8y0pLg9zceJcm7opKfIyZn8NLs9dTo2oVHr+uEwN6tHR0eohIaK4hpVTFtmyJ7PkEsmrHIf48IZtVOw5zdZdmPNW/M2fUtTdJXLxoIFBKnZKaGvyMIIEvD3u8qIR/fbGe1+fm0LBWNV4b2J2ruza3u1hxpYFAKXXKyJE/HSMAqFXL/3wCWrZ5P3+ekM3GvUcZ0KMlw689mwa1nJ8kzmoaCJRSp5wcEB42zN8dlJrqDwIJNlB8tLCYv81Yy9hFuZxZvyZjf92LS9o3sbtYttHpoyp8Oq3QGwYO9A8M+3z+fxMsCMxbt5e+L85j7KJcBp2fzow/9vF0EAA9I1Dh0mmFyuUOHjvBM1N/ZMKybbRpUpuP7z+fjHR3Jomzmk4fVeHx+LRC5W6fr9jJ45+u4sCxE9zfpw2/v8zdSeIqS6ePquh4eFqhcq89+cd5YvIqpq/aRecz6zH21z3pfKb7k8RZTQOBCo8HpxUq9xIRJizbxtNZqzle7CPzqo7cd3HrhEkSZzX9VlR4Ro70TyMsLYGnFXpSgkwG2Lr/GHe9/S2PTsimQ7O6fP6Hi3nw0rM0CJRDzwhUeDwyrdCzEmAywMkkcS/MWIsBnu7fmYG906hSJTHSQ8SSDhYr5Xbjx0cfoF0+GWDDnnwyJ65g2eYDXNK+CSNv7ELLhombJK6ydLBYqURk1ZG8SycDFJX4GD0vh5dnr6dW9SRG3XION57XImGSxMWLdpp5RYL0/6oyhg37aToI8D8eNiyyzwk16O/gyQArtx+i3ytf87cZa7miU1Nm/fESbuqeOJlC40nPCLwgAfp/VQhWHcm7KMfQ8aISXpq9njHzc2hUuxpv3NmDKzs3s7tYrqZnBF5g1VGjch6rjuQHDoTRo/1jAsb4/x092nEHCt9u2s81L8/n9bkbubl7C2b/8RINAhbQMwIvcGn/rwqDlUfyAwc6ruE/6UhhMS9MX8N7izbTsmFNxt3bm4vaNba7WAlDA4EX6GKwxOWBab1frd3DsE9WsuNQAfdcmM6jV3agVjVtuqyk36YXuKj/V1WCg4/ko3Hg6AmenrqaSd9tp+0ZdZjwwAX0SGtod7ESkgYCL/DAUaNKHCLCtBW7+MuUlRw8VsTvftGWh37RlupVvZckLl50sNgrYp1jXqenKgvsOXyc+/+zjP95/zua16/JlIcu4n/7dtAgEGMaCBJRvBvlk9NTN28GkVPTUzUYqDCJCB8t2cplo+Yyd91ehl7dkU9+ewGdzqxnd9E8QVNMJJqyawbAPx4Qy6mALk9PoOy1df8xhk5awYIN++jVuhHP39yN1o1r212shBQqxYQGgkRjR6NcpYr/TKAsY/xdUUoFUeIT3l2Yy99nrCWpimHI1R25vVeqJomLIc015BV2rBnQ6akqQut35/Pnidl8v+UgP+/QhJE3duXMBjXtLpZnWTJGYIy5yhiz1hizwRgzJMjr1Y0xHwZeX2yMSQ88n26MKTDGLA/cXreiPJ5mR84YvVaBCtOJYh//nLOea/45n9x9R3npV+fy9t09NQjYLOpAYIxJAl4FrgY6AbcZYzqV2exe4ICItAVeBJ4v9dpGETk3cHsg2vJ4nh2NskvSE1hJJ0lFLnvbQfq9soBRs9ZxVZfmzH7kEm7QTKGOYEXXUC9gg4jkABhjPgD6A6tLbdMfeDJwfwLwitG/fmzYtWYgQRc1BaM5/CJTcKKEl2avY8z8HJrUrc6YuzK4olNTu4ulSol6sNgYMwC4SkTuCzy+E+gtIg+V2mZlYJttgccbgd5AHWAVsA44DAwXkfkh9jMYGAyQmpraY3OwPmml4kAnSYXvm5w8hkzMJjfvGLf1asXQa86mXo1ku4vlWU4dLN4JpIpInjGmBzDZGNNZRA6X3VBERgOjwT9rKM7lVOq/NIdfxfKPF/Hc52sYv3gLrRrV5P37enNBW00S51RWBILtQKtSj1sGngu2zTZjTFWgPpAn/tORQgARWRY4U2gP6NxQ5Vg6Sap8X6zZzbBPVrL78HHuu6g1j/Rtr0niHM6KWUNLgHbGmNbGmGrArcCUMttMAQYF7g8AvhARMcY0CQw2Y4xpA7QDciwok1Ixk/CTpCo5Ep53pJA/fPA9v353KXVrVGXigxcw/LpOGgRcIOpAICLFwEPADOBH4CMRWWWMGWGM6RfY7C0gxRizAXgEODnFtA+QbYxZjn8Q+QER2R9tmVQZOsXFUq6fJFXe76ES6UJEhCk/7OCKF+cxbcVO/nBZO7J+dzHnpWqmULfQlcWJzo6UE8q5Kvo9RDgSvuvQcYZPXsHsH/dwTsv6PD+gGx2baX4gp9IUE16lU1xUaRX9HsJMFyIifLBkK89O/ZEin48/9e3APRe2JknTQzhaqECg2UedyMquHJ3iEjtu7HKr6PcQxsr03H1HuX3MYoZOWkGXFvWZ8XAf7ru4jQYBF9NA4DRWp3S2I+WEF7g19XZFv4dyRsJLfMKYeTlc9fI8Vm4/xHM3deX93/QmLUUzhbqeiLju1qNHD0lYaWki/qblp7e0tMp93rhxIrVq/fSzatXyP5/oxo3zf2/G+P+1ss5W/53iJZzfQ5Dvbc3Ow9LvX/MlLTNL7n33W9l5sMCuGqgoAEslSJtqe6NemVtCBwJjgjcwxlT+M2PZIDpVrANgLP5O8RLB7+F4UbGMmrlW2j42VbqPmClTlm8Xn88Xt6Iqa4UKBDpY7DQ6uGuNWH+PHvg7fb/lAJkTs1m3+wg3nHsmT1zfmUa1q9ldLBUFHSx2C6euVnLbwGisB8md+neywLETxTydtZqbXltI/vFi3r47g5duPU+DQALTQOA0TlytZMPAaNRxJ9aD5E78O1lg4YZ9XPXSfN5asInbe6Uy8499+EVHzRSa8IL1Fzn9ltBjBE4U54FRS7r34zVIPm6cSErKqX2kpMRkDCZot76FYz+HCk7IkIk/SFpmllzywheyaOM+awquHAUdLFaVFueBUcviTiUayojeMm6cSHLy6QWtVs3SYBAspt2dPE6KqlkT6Gat2iW9Rs6S1kOy5Nmpq6XgRLFlZVfOEioQ6GCxqlicB0bDXNxquYizcYT6XsDS7ybYbjaRTjrR/U32HSnkqc9W89kPO+jYrC4vDOhGt5YNoi2ucjAdLFaVF+eBUcu798MccBg27KdBAPyPhw0L8bnlDTxbuHI72EelUvnBcBFh8vfbuWLUXKav3MkjV7RnykMXxSUIuG3OgWcEO01w+k27hmwQx7UIlnbvR/BhEfeAherDsnj8JNhuNhFi3xXsd/uBY3LPO99KWmaW3PDqAlm367Bl5ayIl9c2OgU6RqDcxLK4E8GAQ8RjEy4aIygp8cl/FuVK5yemS8fhn8ub83OkuCS+C8Pcuhg7kWggUN4UwWF+pY5YXTBrKGfvEbnl9YWSlpklA8d8I1vyjlpevnC4eTF2oggVCHSwWCW2CAe6x4+HxX8YzyN5w0hlC8dSUqnz8khXrg8oLvHx5oJNvDhrHdWqVuHxazvxy4yWGGNPllAPLMZ2PKdevF6p2Bo5MvhUoBAD3QO//i0D978O+A+Q6uQFFs+Bq4LB6h2HyZyYzYrth+jbqSlP39CFpvVq2Fae8ePhyJHTn0+Qxdiup7OGVGKLZAXw+PHw+uunz10td+qQsxQWl/CPmWvp98oCdh4q4NXbu/PGnT1sDwKDB0Ne3k+fT0lJiMXYCUG7hpT9xo/3N7RbtvjniI60qSumvHUBsV7EYIFlm/1J4jbsOcJN3Vvw+LWdaOiA/EDaJeQcuo5AOVMs8xhFOmm9vDn4lVjEEK8580cLi3nqs1UMeH0hxwqLefeenoy65VxHBAHQi+S5QrARZKffdNZQAonVnMLKTAEKVRZjIp4JFK858/PX7ZULn5sjaZlZ8vjkFZJ/vCjqz7R6yYhOG3UOdPqocqRYzSmsTOsTrPU2RuTBByPe/e9Sxskm0qQEI5tIk9sYZ2njd/DYCXn04+WSlpkll/7tS1mck2fJ544b518GYeWyCF1I5hwaCJQzxepwsbIBxorD4XHj5Ag/bfmOUEtuY5wlc+anr9wpPZ+ZJW2GTpU7//6jpLUptuzovfSSiNK3lJToPteLF8lzolCBQAeLlT1ODhBv3uwfiC39Oyw301uY7ByhDLHvXNK4NC230rvfm1/Ik1NWMXXFTs5uXo+f1+7GU3+oH36SvDCUt8TAhU2FKkMHi5VzlB4gBn8Lc7IFsuoCL3ZeQSzEKGgqWyq1exFh0nfbuOLFucxavZtHr+zAlIcu5N/P1A87SZ4me1PlCnaa4PSbdg25XLxGD0P1R8S6nyJE/fJT0iL+qG0Hjsldby2WtMwsuenfX8v63fn/fS3c3q9I+uhj1TWknAEdI1COYWfSmXiMXFqwj5ISn7y3cJN0evxzOfvxz+WdBTlSUiZJXLjxNJK4GyyPXnKy9uknCg0EyjnsnE9o99lIGDbsyZcBr30taZlZcseboZPEhRtvIo27OrCbuEIFAh0sVvEX8aXALGTX5c/CUFTiY8z8HF6avZ6ayUk8fl0nbu7eotwkceEsytaVveokTTqnnONkS2VHWonU1OCtYqUvf2aNldsPkTkxm1U7DnNV52aMuKEzZ9StOD/QwIEVf20R5t1THqSzhpQ9Bg70H476fP5/wwkCVkx9sXM2URDHi0r424w19H/1a3YfLuS21O5Mf7wHzerXsGx2TyR595Q3adeQcgcru5MckuRuae5+/jwxm5y9RxnQoyUdCzrxhweTbekxU94QqmtIA4FyhwTq6D5aWMwL09fw3jebObN+Tf56U1f6tG+SSFVUDqULypS7JUgKy3nr9tL3xXm8981mBp2fzsw/9qFP+yaAe6qoi9MSjwYCFTk7WoJQg7mVHeSNcx0OHjvBnz7+gbve/pYayVX4+P7zebJfZ2pXPzVfw+oqxkIss4YrGwWbU+r0m64jsJGVC7IimbBe3n4jnfj+4IOnT66PYTrMadk7pMfTs6R15lRJu+5HMVWLgxbTDVk6NaW0u6ELypQlrGoJKtPqBWvwI/2cceNCr7BKSrJ0FdXuwwXywH+WSlpmlvzsL/OkburBCovp9MVcdi4KV9GLaSAArgLWAhuAIUFerw58GHh9MZBe6rWhgefXAleGsz8NBDayqiWwKqBE+jmhti97i+JQ3OfzyUdLtkjXv0yXdsOmyatfrpe09BKnL2gOi54RuFvMAgGQBGwE2gDVgB+ATmW2+S3weuD+rcCHgfudAttXB1oHPiepon1qILCRVS2BVQEl0s8Jtb1FrduWvKNyx5vfSFpmlgx47WvZsCff0uqWxyVplJSNQgUCKwaLewEbRCRHRE4AHwD9y2zTHxgbuD8BuMz41833Bz4QkUIR2RQ4M+hlQZmUVcoOql5zjTULsqwaGY30cyL5/Aim6/h8wrtfb+LKl+bx3eYDjOjfmQ8Hn89ZTepYVpyKDBtG2GmpK0sXpyUmKwJBC2BrqcfbAs8F3UZEioFDQEqY7wXAGDPYGLPUGLN07969FhRbVSjYFJGxY2HQoOhbAqtW+Eb6OSNHQrUwL+oeZiu9YU8+v3xjEU9+tpqe6Y2Y+cgl3HV+OlWqnMoRFI8FzfGaflqZReHK2VwzfVRERotIhohkNGnSxO7ieEOoQ8xp06JvCaw6tKzM5/i7JU+pUuX04BBGK11U4uOVL9ZzzcsL2Lj3CKNuOYd37+lJiwY1LSlmpJw8/VTXHjhcsP6iSG7A+cCMUo+HAkPLbDMDOD9wvyqwDzBlty29XXk3HSOIk0ScIhJqjMMY/9VXwhxlzd56UK56aZ6kZWbJb8cvkz2Hj8el+OVxav+9U8vlRcRwjGAJ0M4Y09oYUw3/YPCUMttMAQYF7g8AvggUagpwqzGmujGmNdAO+NaCMikrOPkQs7JC9ZOIQEEB/Oc/5Z7lHC8q4bnP13DDv78m70ghb9zZg1dv706TutVjV+YwObX/Ph5jFyo6luQaMsZcA7yEfwbR2yIy0hgzAn/0mWKMqQH8BzgP2A/cKiI5gfcOA34NFAMPi8jnFe1Pcw3FiZ3XDYiVUAl9Tionsc/inDyGTFrBpn1H+VVGKx679mzq10yOSTETiYMvAeE5oXINRd01ZMdNu4YiFM3kcqevcIpUsH6KCrq9DheckGGfZEtaZpZc9PwcWbB+rw0Fdy9de+AcxLBrSDlZtMlhEm2KyMn+k6Sk4K+X6fb6cu0ernxxHuMXb+HXF7ZmxsN9uLBt4zgUNDbsGLR12CUgIuOVUe5g0cHpNz0jiECowzGL0ym4TgUjmPuPFMrDH3wvaZlZctk/vpKlufsr/LhoTpxKvz8lJaJx64j2EatB24rq78oTywQc5UZzDXlUOCtpXf7jrrQgrZPP55PPftgu3UfMlLOGTpV/zFgjx4uKy31rSopItWqV/0or6q2y6s8Tqy6aBGwv/RKwT0sDgVeFm1vHTT/uGB1e7jpUIPeNXSJpmVly/b/my+odh0LuvryGO9KvNJw/kRV/nljNBk7A9tIvAadPhwoEOkaQ6IJ10AbjtKufhBKDhPgiwodLtnD5qLnMW7eXx67pyKQHL+Ds5vWCbh9sOmQw4X6l4WxXepvKdlvHajawWy6oE7FEnD4dSrDo4PSbnhFEqPQRdFKSuw/fLD783LzvqNw2epGkZWbJLa8vlE17j1T4nnDz1oU7DBPJGUE03TCx6sJJ2DOCBOzzQruGlIi4/8dt0el6cYlP3pyfIx2Hfy6dn5gu477JlZISX1jvDbe3LdyvOJIxgmgb3Vj0qrn9J1UuV45yh6aBQJ3i5h+3BYefa3cdlv6vLJC0zCy5551vZcfBYxEVIVjDl5x8aqZPZU66wp015NRuazt/Um7+OcebBgKVGKI4/CwsKpGXZq2Tto9NlfNGzJTJ328Tny+8s4BgxQjV+MSqsR43zv09e1ar7M/Bq8FDA4Gyj9X/6yrxecu3HJArX5wraZlZ8tD738m+/NgliYtFn/m4cf6zjki7nSLdh9sax8p81wndlVUBDQTKHjb/rztWWCwjp66W1kOypPfI2TJr1a6Y7zMWVU5JCd7gGePtfv7KnH0l7OB2GEIFAp0+qmLLxtSTizbmcdXL8xg9L4df9Uxl5iN9uLxT05jv92QWi5SUU8/VPP0SBUD4U0Hz8oI/L2JN1g+3ZgitzAzPhJ3uGgUNBCq2bPhfd/h4EUMnreC2Md8A8P5vevPXm7pSr0Z8M4UWFJy6n5d3+nKHGCyJCEuw4OPWxrEyeYy8tDwgbMFOE5x+066hSrCrAzjO5+FzftwlvUfOlvTMLEnrt1qqJBfb0t8dTrUj+WpCdQ2lpERWrlBdQKE+3w3dJZH+tN3aDWYFdIzAw+z85cdp3/vyj8vv/+87ScvMkl5PzJV6rQ/Y+h89nL7rSPq3x407PZ9RtWqRz44pL6A4tnGMwUGMGwfGraCBwGuctJo4hv/rfD6ffLp8u5w3Yqa0fWyqvDhrraS1LrH96NbqMwIRa458Q91ODjo7rnH08uF7DGgg8JJwWwC7VyFFacfBY3Lvu99KWmaW9HtlgazZeVhEnLHoKpz2K9ZtXCQroJ3SBVQ2GOWnpDm7wC6jgcBLEjHjaCklJT4Z/81m6fLEdOkwfJqMmbdRikulhwj3SDvWR8DhfH4syxBuTiSnHGAHC4wlOCCqJxANBF6SwNcgyN13RG59w58k7tY3FknuvtOTxDnhaLz0fuzqbgkVEFNSHNgFJMHLu4kgT7r4IMZuGgi8JAGvSlZc4pPRczdKh+HTpMsT0+X9xZvLTQ9RUQMcj8lMdndv273/SAU7frmNcXIEF1XC4TQQeEksWwAbDnHX7Dws/f41X9Iys+Ted7+VnQcLov7McMcRoqmuE1awOnIAOIRQ39fvUlxUCYfTQOA1CZBvuLCoREbNXCttH5sq3UfMlCnLt1c6SVxZ4TTS0VbXCYPWbuK2Mxg30kCgohfHQ9zvNu+XK0Z9JWmZWfKH//tO8o4UWvr54TQ60VbXCWcEbuOmMxg3ChUINMWECl+M8xCMHw/pZxXT6LLV3PjqQnbvL+btuzN46dbzaFS7miX7OOlkPqC0NDDG/+/o0T/N2xNtdSuT/iAalb2EpZMMHAi5ueDz+f+1Io+SqlhVuwugXCQ11Z8QJ9jzURo/Hv5nxD5qX76Ceg2Pkf99KvsWd2Rnx2ToGPXHBzVwYPkNTbTVPfnZw4b5Pycp6aeJ3Kxs5E7mLTqZOO5k3iKr96MSk54RqPDF6BD3UEERwyZn0+DGxSCwa/zP2D+zK0cPJVuW/bIyR8tWVHfgwFOfU1Lify4WyeXcmj1UOUSw/iKn33SMoJKs6IC1uBN35qpd0mvkLEl9NEsaXLpaTNViywdWoxmEtKK68Rgr0IFpFQ5CjBEY/2vukpGRIUuXLrW7GO5Stu8A/IepZTvG42TfkUKenLKKrOyddGxWlzVju7F5eYPTtktL8/cVRyM9PXgXjxWfHY4qVfzNclnG+PvCrWB3HZU7GGOWiUhG2ee1a8grHNJ3ICJM/n47V4yay8xVu/nfK9rz2e8uYuSfGgTthrnmmugHQO3OtR+P/PfBurGM8QcHtw4cq/jRQOAVdreGwI6DBfz63SU8/OFy0hvXZurvL+J3l7UjOalK0Fk8gwbB2LHRX7jF7guRxGP2UOnvD/zf4cmzkHhd8Ea5WLD+IqffdIygEmyc1F5S4pP3FuVK5yemS8fhn8tb83N+kiQumiKHm9jN7oVK8Zwfr+sXVCjogjKPs6k1zNl7RH75+kJJy8yS28cski15R8N+b0UDoJFUKZYNsdMWQenAsQpFA4GKa4tVVFwir3+1QdoPmyZd/jJdPvx2S8TpISo6sg0nt96DD8Y+1bTdZxtl6RmBCiVUINBZQ8pyq3ccJnNiNiu2H6Jvp6Y8fUMXmtarEfHnlDfRCeCOOyIvm9UTpZw4W8dhE8SUg+isIRVzhcUl/GPmWvq9soCdhwr498DuvHFnj0oFAQidBgJOrZqNlNUTpRwwBn+acNJnKFWanhF4zfjx/pZwyxb/tJmRIy1pIZZtPkDmxGw27DnCTd1b8Pi1nWhocX6gk0IdhYdL5+8rrwp1RqC5hrwkBglpjhYW8/eZa3l3YS7N69XgnXt68vMOZ1hU4OCiPdq2ev5+sG6YWCWWUyoWouoaMsY0MsbMMsasD/zbMMR2gwLbrDfGDCr1/FfGmLXGmOWBW2xbEK+zeFHZ/PV7ufKlebzzdS539E5j5iOXxDwIQOiGPCXl9Pn6ZcVy/r52wyi3inaMYAgwR0TaAXMCj3/CGNMI+AvQG+gF/KVMwBgoIucGbnuiLI8qj0Ud2oeOFfHoxz9w51vfUq1qFT5+4HyevqELdarH5wQz1AKtl18+vVF+8EHrGulQievinTo5EdJNK4cJNpUo3BuwFmgeuN8cWBtkm9uAN0o9fgO4LXD/KyAj0v3q9NFKsmBe4ecrdkrGM7OkzdCp8vznP0rBieKYFbc88Z6775Rpok4ph3InYjF91BhzUEQaBO4b4MDJx6W2+RNQQ0SeCTx+HCgQkb8bY74CUoASYCLwjIQokDFmMDAYIDU1tcfmaEYLvSqKeYV78/1J4qau2Emn5vV4YUA3urSoH+MCO4dTBoWdUg7lTpUeLDbGzAaaBXnpJx3LIiLGmEijykAR2W6MqYs/ENwJvBdsQxEZDYwG/6yhCPej4KdXSglz1pCIMOm77YzIWk3BiRIevbIDg/u0ITnJWzOPnTJN1CnlUImlwkAgIpeHes0Ys9sY01xEdhpjmgPB+vi3A5eWetwSf5cQIrI98G++MeZ9/GMIQQOBskhFl+UqZfvBAh6btIK56/bSI60hz9/cjbZn1IlxAZ0phhdnc2U5VGKJ9rBuCnByFtAg4NMg28wA+hpjGgYGifsCM4wxVY0xjQGMMcnAdcDKKMujLODzCe8tyqXvqLksyd3PU/068/H953s2CIA/HXYkz8dKvK+DrLwh2mkezwEfGWPuBTYDtwAYYzKAB0TkPhHZb4x5GlgSeM+IwHO18QeEZCAJmA2MibI8Kkob9x5hyMRsluQe4OJ2jXn2xq60alTBnEwPmDYtsudjpRK9e0pVSFcWKwCKSnyMnpfDy3PWUzM5ieHXns2AHi3xzwFQ8bjKmFKxpiuLVUgrtx8ic2I2q3Yc5pquzXiyX2fOqFu5/ECJSvvmVSLTQOBhx4tKeHnOekbPy6FhrWq8NrA7V3dtbnexHElTSahEpoHAo5bk7idzYjY5e4/yyx4tGX5tJ+rXSra7WI6lffMqkWkg8JgjhcW8MH0N7y3aTIsGNXnv173o076J3cVyhQhm3irlKhoIPOSrtXsY9slKdhwq4O4L0nn0yg7UjlN+IKWUc3lreahHHTh6gkc+Ws7d7yyhRnIVJjxwPk/26/zfIKBJzJTyNj0cTHDTVuzkiU9XcvBYEb/7RVse+kVbqldN+u/rMbhEgVLKZXQdQYLac/g4T3y6iumrdtGlRT2ev7kbnc88PUmcJjFTyjt0HYFHiAgfL9vGM1mrOV7sI/Oqjvzm4tZUDZEkTpOYKaU0ECSQrfuP8dgnK5i/fh+90hvx3M1dadOk/PxAulBKKaWBIAGUBJLEvTB9LVUMPN2/MwN7p1GlSsXpIXShlFJKA4HLbdiTz58nZPPdloNc2qEJI2/sSosGNcN+vy6UUkppIHCpohIfr3+1kX99sYFa1ZN48VfncMO5LSqVJE4XSinlbRoIXGjFtkM8OuEH1uzK59puzXmqX2ca16lud7GUUi6lgcBFjheV8OLsdYyZl0PjOtUZfWcP+nYOdhVRpZQKnwYCl1ick8eQSSvYtO8ov8poxWPXnk39mpokTikVPQ0EDpd/vIjnp69h3DdbaNWoJuPv682FbRvbXSylVALRQOBgX67Zw2OfrGD34ePce1Fr/rdve2pV0z+ZUspa2qo40P6jJ3g6azWffL+ddmfU4d8PXsB5qQ3tLpZSKkFpIHAQESEreydPTlnFoYIifn9ZO/7n52f9JEmcUkpZTQOBQ+w+fJxhn6xk9o+76dayPuPu683ZzevZXSyllAdoILCZiPDhkq2MnPYjJ4p9DL26I/deFDpJnFJKWU0DgY225B1jyKRsFm7Mo3frRjx/czfSG9e2u1hKKY/RQGCDEp/w7sJc/j5jLUlVDCNv7MJtPVPDShKnlFJW0/6HOFu3O5+bX1vI01mrOf+sFGY90ifsTKHKZnpNT5Wg9IwgTk4U+3jtq4288uV66tZI5uVbz6XfOWdWKkmcsoFe01MlML1UZRz8sPUgmROzWbMrn+vPOZMnr+9EiiaJcxe9pqdKAHqpShsUnPAniXtzfg5n1K3Bm3dlcHmnpnYXS1WGXtNTJTANBDGyaGMeQyZlsznvGLf3TmXI1R2pV0OTxLmWXtNTJTAdLLbY4eNFDJ20gtvGfAPA+7/pzbM3dtUg4HYjR/qv4VmaXtNTJQg9I7DQ7NW7GT55JXvyjzO4Txv+eHl7albT9BAJQa/pqRKYBgIL5B0p5KnPVjPlhx10bFaXN+7swTmtGthdLGU1vaanSlAaCKIgIkz5YQdPfbaa/ONFPHx5O357aVuqVdUeN6WUe2ggqKSdhwoY/slK5qzZwzmtGvDCzd3o0Kyu3cVSSqmIaSCIkM8nfLBkK3+d9iNFPh/Drz2bey5sTZKuDFZKuZQGggjk7jvKkEnZfJOznwvOSuGvN3UlLUWTxCml3E0DQRiKS3y8/fUm/jFzHdWSqvDcTV35Vc9Wmh5CKZUQogoExphGwIdAOpAL3CIiB4JsNx34GbBARK4r9Xxr4AMgBVgG3CkiJ6Ipk9XW7DpM5oRsfth2iMvPbsozN3ShWf0adhdLKaUsE+30liHAHBFpB8wJPA7mb8CdQZ5/HnhRRNoCB4B7oyyPZQqLSxg1ax3X/XMB2w4U8K/bzmPMXT00CCilEk60gaA/MDZwfyxwQ7CNRGQOkF/6OePvV/kFMKGi98fb91sOcP2/FvDPOeu5rltzZj1yCddrplClVIKKdoygqYjsDNzfBUSSUS0FOCgixYHH24AWoTY2xgwGBgOkxii/y7ETxfxj5jre/noTzerV4J27e/LzjmfEZF9KKeUUFQYCY8xsoFmQl4aVfiAiYoyJWU5rERkNjAZ/GmqrP//rDfsYMimbrfsLGBhIEldX8wMppTygwkAgIpeHes0Ys9sY01xEdhpjmgN7Ith3HtDAGFM1cFbQEtgewfstcaigiL9O+5EPlmyldePafDj4Z/RukxLvYiillG2i7RqaAgwCngv8+2m4bwycQXwJDMA/cyii91th5qpdDJ+8kryjJ3jgkrN4+PJ21EjWJHFKKW+JNhA8B3xkjLkX2AzcAmCMyQAeEJH7Ao/nAx2BOsaYbcC9IjIDyAQ+MMY8A3wPvBVlecKyN7+QJz9bxdTsnZzdvB5vDepJ15b147FrpZRynKgCgYjkAZcFeX4pcF+pxxeHeH8O0CuaMkRCRJi8fDtPfbaaY4Ul/Klve+6/5CySkzRJnFLKuzyzsrioxMfg95by5dq9dE9twAsDutH2DE0Sp5RSngkEyUlVaNOkDn3aN+Gu89M1SZxSSgV4JhAAPH5dJ7uLoJRSjqOd40op5XEaCJRSyuM0ECillMdpIFBKKY/TQKCUUh6ngUAppTxOA4FSSnmcBgKllPI4IxKzSwjEjDFmL/4kd5XRGNhnYXHslCh1SZR6gNbFqRKlLtHWI01EmpR90pWBIBrGmKUikmF3OayQKHVJlHqA1sWpEqUusaqHdg0ppZTHaSBQSimP82IgGG13ASyUKHVJlHqA1sWpEqUuMamH58YIlFJK/ZQXzwiUUkqVooFAKaU8LuEDgTGmkTFmljFmfeDfhkG2STPGfGeMWW6MWWWMecCOslYkzLqca4xZFKhHtjHmV3aUtTzh1COw3XRjzEFjTFa8y1gRY8xVxpi1xpgNxpghQV6vboz5MPD6YmNMug3FDEsYdekT+P9RbIwZYEcZwxFGPR4xxqwO/L+YY4xJs6Oc4QijLg8YY1YE2qwFxpjorrolIgl9A14AhgTuDwGeD7JNNaB64H4dIBc40+6yV7Iu7YF2gftnAjuBBnaXPdJ6BF67DLgeyLK7zGXKlQRsBNoEfjs/AJ3KbPNb4PXA/VuBD+0udxR1SQe6Ae8BA+wucxT1+DlQK3D/QZf/TeqVut8PmB7NPhP+jADoD4wN3B8L3FB2AxE5ISKFgYfVce6ZUjh1WSci6wP3dwB7gNNWEtqswnoAiMgcID9OZYpEL2CDiOSIyAngA/x1Kq10HScAlxljnHih7ArrIiK5IpIN+OwoYJjCqceXInIs8PAboGWcyxiucOpyuNTD2kBUs36c2uBZqamI7Azc3wU0DbaRMaaVMSYb2Ir/CHVHvAoYgbDqcpIxphf+I4qNsS5YhCKqhwO1wP87OWlb4Lmg24hIMXAISIlL6SITTl3cINJ63At8HtMSVV5YdTHG/I8xZiP+M+zfR7PDhLh4vTFmNtAsyEvDSj8QETHGBI2cIrIV6GaMOROYbIyZICK7rS9t+ayoS+BzmgP/AQaJSNyP5Kyqh1JWM8bcAWQAl9hdlmiIyKvAq8aY24HhwKDKflZCBAIRuTzUa8aY3caY5iKyM9A47qngs3YYY1YCF+M/pY8rK+pijKkHTAWGicg3MSpquaz8mzjQdqBVqcctA88F22abMaYqUB/Ii0/xIhJOXdwgrHoYYy7HfzBySanuYKeJ9G/yAfBaNDv0QtfQFE5FykHAp2U3MMa0NMbUDNxvCFwErI1bCcMXTl2qAZ8A74lI3ANZmCqsh8MtAdoZY1oHvu9b8deptNJ1HAB8IYGRPYcJpy5uUGE9jDHnAW8A/UTEyQcf4dSlXamH1wLro9qj3SPkcRiBTwHmBL6o2UCjwPMZwJuB+1cA2fhH57OBwXaXO4q63AEUActL3c61u+yR1iPweD6wFyjA3096pd1lL1W2a4B1+MdfhgWeG4G/kQGoAXwMbAC+BdrYXeYo6tIz8P0fxX9Ws8ruMleyHrOB3aX+X0yxu8xR1OVlYFWgHl8CnaPZn6aYUEopj/NC15BSSqlyaCBQSimP00CglFIep4FAKaU8TgOBUkp5nAYCpZTyOA0ESinlcf8PuMAlH5DeQsMAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": null, + "source": [], + "outputs": [], + "metadata": {} + } + ], + "metadata": { + "kernelspec": { + "name": "python3", + "display_name": "Python 3.6.9 64-bit" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.9" + }, + "interpreter": { + "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/homework9/utility.py b/homework9/utility.py new file mode 100644 index 0000000..e22ae20 --- /dev/null +++ b/homework9/utility.py @@ -0,0 +1,140 @@ +import numpy as np +from numpy import sqrt, sum, abs, max, maximum, logspace, exp, log, log10, zeros +from numpy.random import normal, randn, choice +from numpy.linalg import norm +from scipy.signal import convolve2d +from scipy.linalg import orth + +# For unit testing + +def check_gradient(f, grad, x, error_tol = 1e-6): + y = normal(size=x.shape) + y = y/norm(y)*norm(x) + + g = grad(x) + rel_error = np.zeros(10) + for iter in range(10): + y = y/10; + d1 = f(x + y) - f(x) # exact change in function + d2 = np.sum((g * y).ravel()) # approximate change from gradient + rel_error[iter] = (d1-d2)/d1 + #print('d1=%1.5g, d2=%1.5g, error=%1.5g,'%(d1,d2, rel_error[iter] )) + min_error = min(np.abs(rel_error)) + print('Min relative error = %1.5g'%min_error) + did_pass = min_error < error_tol + print(did_pass and "Test passed" or "Test failed") + return did_pass + +def check_adjoint(A,At,dims): + # start with this line - create a random input for A() + x = normal(size=dims)+1j*normal(size=dims) + Ax = A(x) + y = normal(size=Ax.shape)+1j*normal(size=Ax.shape) + Aty = At(y) + # compute the Hermitian inner products + inner1 = np.sum(np.conj(Ax)*y) + inner2 = np.sum(np.conj(x)*Aty) + # report error + rel_error = np.abs(inner1-inner2)/np.maximum(np.abs(inner1),np.abs(inner2)) + if rel_error < 1e-10: + print('Adjoint Test Passed, rel_diff = %s'%rel_error) + return True + else: + print('Adjoint Test Failed, rel_diff = %s'%rel_error) + return False + +# For total-variation + +kernel_h = [[1,-1,0]] +kernel_v = [[1],[-1],[0]] + +# Do not modify ANYTHING in this cell. +def gradh(x): + """Discrete gradient/difference in horizontal direction""" + return convolve2d(x,kernel_h, mode='same', boundary='wrap') +def gradv(x): + """Discrete gradient/difference in vertical direction""" + return convolve2d(x,kernel_v, mode='same', boundary='wrap') +def grad2d(x): + """The full gradient operator: compute both x and y differences and return them all. The x and y + differences are stacked so that rval[0] is a 2D array of x differences, and rval[1] is the y differences.""" + return np.stack([gradh(x),gradv(x)]) + +def gradht(x): + """Adjoint of gradh""" + kernel_ht = [[0,-1,1]] + return convolve2d(x,kernel_ht, mode='same', boundary='wrap') +def gradvt(x): + """Adjoint of gradv""" + kernel_vt = [[0],[-1],[1]] + return convolve2d(x,kernel_vt, mode='same', boundary='wrap') +def divergence2d(x): + "The methods is the adjoint of grad2d." + return gradht(x[0])+gradvt(x[1]) + + +# For logistic regression + +def buildmat(m, n, cond_number): + """Build an mxn matrix with condition number cond.""" + if m <= n: + U = randn(m, m); + U = orth(U); + Vt = randn(n, m); + Vt = orth(Vt).T; + S = 1 / logspace(0, log10(cond_number), num=m); + return (U * S[:, None]).dot(Vt) + else: + return buildmat(n, m, cond_number).T + + +def create_classification_problem(num_data, num_features, cond_number): + """Build a simple classification problem.""" + X = buildmat(num_data, num_features, cond_number) + # The linear dividing line between the classes + w = randn(num_features, 1) + # create labels + prods = X @ w + y = np.sign(prods) + # mess up the labels on 10% of data + flip = choice(range(num_data), int(num_data / 10)) + y[flip] = -y[flip] + # return result + return X, y + + +def logistic_loss(z): + """Return sum(log(1+exp(-z))). Your implementation can NEVER exponentiate a positive number. No for loops.""" + loss = zeros(z.shape) + loss[z >= 0] = log(1 + exp(-z[z >= 0])) + # Make sure we only evaluate exponential on negative numbers + loss[z < 0] = -z[z < 0] + log(1 + exp(z[z < 0])) + return np.sum(loss) + + +def logreg_objective(w, X, y): + """Evaluate the logistic regression loss function on the data and labels, where the rows of D contain + feature vectors, and y is a 1D vector of +1/-1 labels.""" + z = y * (X @ w) + return logistic_loss(z) + + +def logistic_loss_grad(z): + """Gradient of logistic loss""" + grad = zeros(z.shape) + neg = z.ravel() <= 0 + pos = z.ravel() > 0 + grad[neg] = -1 / (1 + exp(z[neg])) + grad[pos] = -exp(-z[pos]) / (1 + exp(-z[pos])) + return grad + + +def logreg_objective_grad(w, X, y): + return X.T @ (y * logistic_loss_grad(y * X @ w)) + +# For gradient descent +def estimate_lipschitz(g, x): + # Your work here + y = x+normal(size=x.shape) + L = norm(g(x)-g(y))/norm(x-y) + return L