logistic_regression_assignment_MontaserMohammedalamen.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistic regression tutorial "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Do **not** change this cell, and do **not** import\n",
    "## any other modules anywhere in the notebook.\n",
    "import numpy as np\n",
    "import numpy.random as rn\n",
    "from scipy import optimize, stats\n",
    "import scipy.linalg as linalg\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "In this tutorial we're going to cover the basics behind logistic regression. For simplicity we will only consider the binary classification case, in which target variables are $y \\in \\{0,1\\}$. \n",
    "\n",
    "In logistic regression, the probability of a data point $\\boldsymbol x$ being of class 1 is given by\n",
    "\n",
    "$$p(y = 1 | \\boldsymbol x, \\boldsymbol\\theta) = \\sigma (\\boldsymbol x^\\top \\boldsymbol\\theta) ~ ,$$\n",
    "\n",
    "where $\\sigma(z) = 1/(1+\\exp(-z))$ is the _sigmoid_ function.\n",
    "\n",
    "Combining this with a Bernoulli likelihood and summing over all datapoints $\\{\\boldsymbol x_i, y_i\\}_{i=1}^N$ we end up with a negative log-likelihood function that looks like this:\n",
    "\n",
    "$$-\\log p(\\boldsymbol y|\\boldsymbol X, \\boldsymbol\\theta) = -\\sum_i\\left(y_i \\log \\sigma(\\boldsymbol x_i^\\top \\boldsymbol\\theta) + (1 - y_i) \\log ( 1 - \\sigma(\\boldsymbol x_i^\\top \\boldsymbol\\theta))\\right)$$\n",
    "\n",
    "You will see this expression in many other classification problems, especially in deep learning, where it's known as the _cross-entropy loss_.\n",
    "\n",
    "Your goal in this tutorial is to learn how to perform inference over the parameters $\\boldsymbol\\theta$ in logistic regression, including point estimates $\\boldsymbol\\theta_{\\mathrm{ML}}$ and $\\boldsymbol\\theta_{\\mathrm{MAP}}$ and approximations to the posterior $p(\\boldsymbol\\theta | \\boldsymbol X, \\boldsymbol y)$.\n",
    "\n",
    "Let's do it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Maximum likelihood estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start easy. First, let's generate a toy 1D binary dataset with two paramaters:\n",
    "\n",
    "* A **jitter** parameter that controls how noisy the data are; and\n",
    "* An **offset** parameter that controls the separation between the two classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Data generation parameters\n",
    "N = 50\n",
    "D = 2\n",
    "jitter = 0.7\n",
    "offset = 1.2\n",
    "np.random.seed(42)\n",
    "# Generate the data\n",
    "x = np.vstack([rn.normal(0, jitter, (N//2,1)), rn.normal(offset, jitter, (N//2,1))])\n",
    "y = np.vstack([np.zeros((N//2, 1)), np.ones((N//2, 1))])\n",
    "x_test = np.linspace(-2, offset + 2).reshape(-1,1)\n",
    "\n",
    "# Make the augmented data matrix by adding a column of ones\n",
    "x = np.hstack([np.ones((N,1)), x])\n",
    "x_test = np.hstack([np.ones((N,1)), x_test])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now on to the regression. First, let's code up the logistic log-likelihood as a separate function. This will come in handy.\n",
    "\n",
    "**Task 1**\n",
    "\n",
    "* Write a function to calculate the log-likelihood of a dataset given a value of $\\boldsymbol\\theta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## EDIT THIS FUNCTION\n",
    "from numpy.linalg import inv\n",
    "from numpy import linalg as LA\n",
    "import scipy\n",
    "def sigmoid (a):\n",
    "        return (1+np.exp(-a))**(-1)\n",
    "def log_likelihood(X, y, theta):\n",
    "    # X: N x D matrix of training inputs\n",
    "    # y: N x 1 vector of training targets/observations\n",
    "    # theta: parameters (D x 1)\n",
    "    # returns: log likelihood, scalar\n",
    "    \n",
    "    L = - np.inf ## <-- EDIT THIS LINE\n",
    "    sigma=sigmoid(X @ theta)\n",
    "    L=(y.T @ np.log(sigma) + (1-y).T @ np.log(1- sigma))\n",
    "    return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-35.23792467])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_likelihood(x,y,np.array([1,2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it's time to optimize it to fit the maximum likelihood parameter,\n",
    "\n",
    "$$\\boldsymbol\\theta_{\\mathrm{ML}} = \\mathrm{arg}_\\theta \\max p(\\boldsymbol y | \\boldsymbol X, \\boldsymbol\\theta)$$\n",
    "\n",
    "For linear regression, the likelihood function had a closed-form minimum, which made our lives easy. Alas, that is not the case for logistic regression. We will have to resort to _numerical optimization_.\n",
    "\n",
    "In the lectures you saw how to derive the gradient and all that jazz. For this tutorial you can do it that way, or any other way you want. The optimization is convex, so this should be easy peasy.\n",
    "\n",
    "**Task 2**\n",
    "\n",
    "* Write a function to optimize the log-likelihood function you've written above an obtain $\\boldsymbol\\theta_{\\mathrm{ML}}$. Use any optimizer of your choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## EDIT THIS FUNCTION\n",
    "def max_lik_estimate(X, y):\n",
    "    # X: N x D matrix of training inputs\n",
    "    # y: N x 1 vector of training targets/observations\n",
    "    # returns: maximum likelihood parameters (D x 1)\n",
    "    \n",
    "    N, D = X.shape\n",
    "    theta_ml = np.zeros((D,1)) ## <-- EDIT THIS LINE\n",
    "    fun = lambda w: -log_likelihood(X,y,w)\n",
    "    theta_ml=scipy.optimize.minimize((fun),np.random.normal(size=(D,1)),method='BFGS',options={'gtol': 1e-10, 'disp': False})\n",
    "    return theta_ml.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.13633653  2.57591922]\n"
     ]
    }
   ],
   "source": [
    "theat2=max_lik_estimate(x,y)\n",
    "print (theat2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 3**\n",
    "\n",
    "* Write a predict function to evaluate your estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## EDIT THIS FUNCTION\n",
    "def predict(X, theta):\n",
    "    # Xtest: K x D matrix of test inputs\n",
    "    # theta: D x 1 vector of parameters\n",
    "    # returns: prediction of f(Xtest); K x 1 vector\n",
    "    \n",
    "    prediction = sigmoid(X @ theta)\n",
    "    \n",
    "    return prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this we're in a good position to fit a logistic regression to our toy dataset and start visualising the results. Have a go.\n",
    "\n",
    "1. Use the function you wrote above to estimate $\\boldsymbol\\theta_{\\mathrm{ML}}$ on the toy dataset.\n",
    "2. Visualize the results, including:\n",
    "    1. The data $x$ and target labels $y$. \n",
    "    2. The labels predicted by the model.\n",
    "    3. The probability assigned by the model, $\\sigma(x\\theta)$ as a function of $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ac,curecy= 78.0 %\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f4975f46208>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYXHWd7/H3t5Ze0h2SkDQkZKFBcCTIEuhBEQRUVECv\nyKID47ghZnT0yjjOc0VRvHIHxeVBLuDViciIGyq4RQ0jQVaVAB0MhBAhMSjEBNIkpDudXqvqe/+o\nE1LprkpO6pzU6bI+L57z1KnfOVW/DyfV9a2zm7sjIiKNJ5V0ABERSYYKgIhIg1IBEBFpUCoAIiIN\nSgVARKRBqQCIiDQoFQARkQalAiAi0qBUAEREGlQm6QC7M2PGDO/s7Ew6hohI3Vi+fPnz7t4RZt4J\nXQA6Ozvp7u5OOoaISN0ws7+EnVebgEREGpQKgIhIg1IBEBFpUCoAIiINSgVARKRBRS4AZjbXzO4y\ns9VmtsrMLikzz2lm1mtmK4Lh8qj9iiTBRx+nsO1qCtuuwUfXJB1HJJI4DgPNAR9z94fNbDKw3MyW\nuvvjY+a7z93fHEN/IokobLsatn8LGAEM334j3v5hUu0LE04mUp3IawDuvtHdHw7GtwGrgdlR31dk\nIvHRJ4Iv/yGgAOSL4/3X4bmnE80mUq1Y9wGYWSewAHigzOQTzewRM7vNzI6Ms1+Rfc2HlgKj5abA\n8J21jiMSi9jOBDazduDHwL+6e9+YyQ8DB7t7v5mdBfwMOLzC+ywEFgLMmzcvrngikZhlcKzMlBQT\n/IR6kYpiWQMwsyzFL//vuftPxk539z537w/GlwBZM5tR7r3cfZG7d7l7V0dHqMtZiOx7LWcC6TIT\nHFpeX+s0IrGI4yggA74JrHb3qyvMMzOYDzM7Ieh3c9S+RWrFMgfD5EuBZqAFaC2O73cFlj4w2XAi\nVYpj3fUk4J3ASjNbEbR9EpgH4O5fB84HPmhmOWAQuMDdPYa+RWom1fYOvOV1MHwXYNB8OpYuuyIr\nUhciFwB3/y2U3ThaOs/1wPVR+xJJmqVnwqQLk44hEgudCSwi0qBUAEREGpQKgIhIg1IBEBFpUCoA\nIiINSgVARKRBqQCIiDQoFQARkQalAiAi0qBUAEREGpQKgIhIg1IBEBFpUCoAIiINSgVARKRBqQCI\niDQoFQARkQalAiAi0qDiuCfwXDO7y8xWm9kqM7ukzDxmZtea2Voze9TMjovar4iIRBPHPYFzwMfc\n/WEzmwwsN7Ol7v54yTxnAocHwyuArwWPIiKSkMhrAO6+0d0fDsa3AauB2WNmOxv4thctA6aa2ayo\nfYskwX0Q96GkY4hEFus+ADPrBBYAD4yZNBt4puT5esYXiR3vsdDMus2su6enJ854IpEUcn9h8Pnz\nGXj2SAaePZLBzf9IIb8h6VgiVYutAJhZO/Bj4F/dvW/s5DIv8XLv4+6L3L3L3bs6OjriiicSiRcG\nGNx8DoXR5UAeyFEYWcbQ8+fiPpJ0PJGqxFIAzCxL8cv/e+7+kzKzrAfmljyfA+ink9SN3NCvoDAI\nFEpa87j3kR+6I6lYIpHEcRSQAd8EVrv71RVmWwy8Kzga6JVAr7tvjNq3SK0Ucn8GBsZP8GEK+adr\nHUckFnEcBXQS8E5gpZmtCNo+CcwDcPevA0uAs4C1FP+K3htDvyI1k87OJ2dt4Nt3nWDNpLJHJBNK\nJKLIBcDdf0v5bfyl8zjwoah9iSQl3fIGbNuX8fwzwGjQ2oSlDybd9Ooko4lUTWcCi4RglqV1xk/J\nTHo72FSwaWQmvYPW6bdgpj8jqU9xbAISaQiWmkrzlM/RPOVzSUcRiYUKgEhI7gVyQ0sZGfwJkKJp\n0tvJNJ9G8TgIkfqjAiASgrsz8MIljA7fDl48Gmh0+E6aWs9l0tTPJ5xOpDraeCkSQn50+S5f/gD4\nACMDt5IfXZ1cMJEIVABEQhgdugf3AQo4OQrkKFDAcXKMDt+bdDyRqmgTkEgY1k4ew0vOBM7jGAWK\nF8EVqT9aAxAJIZU5FCc/rt3JYZnDE0gkEp0KgEgIo6OPUP58xwy50eW1jiMSCxUAkRBSqXYgO36C\nZSleCFek/qgAiITQ3Ho2lf5cmlvPqm0YkZioAIiEkE7PYvK068Bai7/4rR2zdqZMu4FUav+k44lU\nRUcBiYTU0noWzc2nMTLyOyBFU/OrMGtNOpZI1VQARPaCpSbR3PL6pGOIxEKbgEREGpQKgIhIg4rr\nnsA3mtkmM3uswvTTzKzXzFYEw+Vx9CsiItWLax/At4DrgW/vZp773P3NMfUnIiIRxVIA3P1eM+uM\n471EJip3p39oKVu334qRZmrb22lr0f0ApH7V8iigE83sEWAD8O/uvqqGfYtE4u78dctH6BtYQoEh\nAPoGb2da29uZtb/uByD1qVY7gR8GDnb3Y4DrgJ9VmtHMFppZt5l19/T01CieyO4Njixn68Bi8gzh\ngAN5hti8/fsMjfwx6XgiValJAXD3PnfvD8aXAFkzm1Fh3kXu3uXuXR0dHbWIJ7JHW/t/HFwN1HYZ\nnDy9AxV/z4hMaDUpAGY204INpWZ2QtDv5lr0LRKH4fzTlaflKk8Tmchi2QdgZjcDpwEzzGw98BmC\nSye6+9eB84EPmlkOGAQucHePo2+RWmjJzqd/+J6y01qz82ucRiQecR0FdOEepl9P8TBRkbq0f/sF\nbO6/AWd0l3Yjy9S2cxNKJRKNzgQWCaEpcyhus3Fnl8HSh5PNHJR0PJGqqACIhLB18Df05jawrdDE\nkKcZ9DR9hSa2jqylb6g76XgiVVEBEAlhfe8NwS3gUwx7lhHP4qQoYGzoXZR0PJGqqACIhDCcr3xO\nylD+2RomEYmPCoBICK2ZwypOm9R0ZA2TiMRHBUAkhObsoRRP/hrLaM501jiNSDxUAETCsCbyFC8B\nUeQ4UCAFlk4slkgUKgAiIUxrOQ6nneFChuFCmuFChqFChgKtTGk+Nul4IlXRPYFFQpjWcjKD3s6Q\nw85NQU4b+zOl+fgEk4lUT2sAIiFsHv4Do15g1/0AxlBhkK26GqjUKRUAkRB6Bh4k74PAzrOAAQqe\no2fwoQSTiVRPm4BEQmhKT8Vppj9v5IPfTRnytKdTNKWmJJxOpDpaAxAJYeak0+jLp4Mv/+K9AHKk\n6c0bM9tOSTqeSFVUAERCeKr/HvzFG8HsYDjG09t+l1QskUhUAERCeGb7sorTnt5+fw2TiMRHBUAk\nhLQ1U+4WRu6QsdbaBxKJgQqASAgtmdnFm8GXFAH34pnBrZk5ScUSiSSWAmBmN5rZJjN7rMJ0M7Nr\nzWytmT1qZsfF0a9I7bSwJdfG9nwz2/LN9OWb6c83sSXXBqls0uFEqhLXGsC3gDN2M/1M4PBgWAh8\nLaZ+RWoi5wUGvYkthTZeyLexNd/GC4V2Br2JvO5uLXUqlgLg7vcCW3Yzy9nAt71oGTDVzGbF0bdI\nLWweeb7kKKCdg2NsGX4+2XAiVarVPoDZwDMlz9cHbSJ14YXRzRWnbR6tfLMYkYmsVgWg3IXUy644\nm9lCM+s2s+6eHv1hycTQnt6fSvcDaM9MrXUckVjUqgCsB+aWPJ8DbCg3o7svcvcud+/q6OioSTiR\nPTmp45yK006ecV4Nk4jEp1YFYDHwruBooFcCve6+sUZ9i0TW2X4UR04+me25Jp4bbue54Xa257Mc\nN+UNHNh6SNLxRKoSy8XgzOxm4DRghpmtBz4DZAHc/evAEuAsYC0wALw3jn5FasXd2TjSTG9uGqM+\nCkBhtI3nRpsSTiZSvVgKgLtfuIfpDnwojr5EkvBk/5M83vf4i1/+ACM+SvcL3bxh+xs4uO3gBNOJ\nVEdnAouEsKp3FcOF4XHt+UKex/seTyCRSHQqACIhtGfaydr4M34zqQxtmbYEEolEpwIgEsIrpr8C\nK3sYKHRN66pxGpF4qACIhDA5PZmB3FQKzotD3o3RfAeTMpOSjidSFd0SUiSEX238PZuGh4H9aErl\nARgppIFt3NfzB17dsSDRfCLV0BqASAh3b3o4GDNGChlGChl2nBl8+7O6KbzUJ60BiITQni3e9KXg\nMJIv/tk0pXOkDCZrE5DUKa0BiITwzoPPZHA0w5bBSfSPNNE/0sSWwUkMjaZ5Z+furoQuMnFpDUAk\nhJb0JIZyrUBhl/ah3CSyqeZkQolEpDUAkRCWblxZ9iDQtKW489lVNc8jEgcVAJEQcoU8hTJ3hfdg\nmkg9UgEQCeHUA48gm0qPazfg1Qe8rPaBRGKgAiASwkv3m8U/HHwiLekshpHCaE5lueglpzG3bXrS\n8USqop3AIiFd8rIzOH3my1n67EpSGGccdAwv3U+3tpb6pQIgshcOmzyT3sE8hnFI+wFJxxGJRAVA\nJKS716/jw/f8/MXnaUvx9decw4mz5iWYSqR6sewDMLMzzOwJM1trZpeWmf4eM+sxsxXBcHEc/YrU\nyqaBfj5w10/pHx15cegdGeJ9v7mV3uGhpOOJVCVyATCzNPBV4ExgPnChmc0vM+sP3f3YYLghar8i\ntbT4qdU44w8DBfjvvzxR4zQi8YhjDeAEYK27r3P3EeAHwNkxvK/IhNE3MsRIPh8c+G/FwWG0kKd3\nRGsAUp/iKACzgWdKnq8P2sY6z8weNbNbzWxuDP2K1MzJB3WSKWSgPw2DqeLQn8ZyKU6a1Zl0PJGq\nxFEAyp0hP3Zd+RdAp7sfDdwB3FTxzcwWmlm3mXX39PTEEE8kupnN7eS2g435b3Q7HNg6Oel4IlWJ\nowCsB0p/0c8BNpTO4O6b3X3HHbW/ARxf6c3cfZG7d7l7V0dHRwzxRKK7+oHfV5x27YOVp4lMZHEU\ngIeAw83sEDNrAi4AFpfOYGalZ8u8BVgdQ78iNfPXbX3j12sD6/v6ahtGJCaRC4C754APA7+m+MX+\nI3dfZWZXmNlbgtk+YmarzOwR4CPAe6L2K1JLr5xdYbeVwynzOmuaRSQusZwI5u5LgCVj2i4vGf8E\n8Ik4+hJJQkdTG+QAg1Su2FbIAgWY3qQ7gkl90sXgREK4c+060oNGdpuRHiwO2b7i4+1Prk06nkhV\ndCkIkRC2DQyRytu4Q95Sedg+NJJIJpGotAYgEsLocOWbvowM52qYRCQ+KgAiIcxoawMHy0FqqDhY\nDnCY3qZ9AFKfVABEQrjoxONJD0N6CNK54o7g9BCkR+D9J/190vFEqqICIBLCfk3NZAs79wFYMGTz\nRmtKu9KkPqkAiIRw3x+fwgvjzwRLm/HbJ/5c+0AiMVABEAkhm06XPRE4ly/QnNEagNQnfXJFQmhJ\nZV68FIQFBwR5GtygNa0/I6lP+uSKhNC7fYjUMGSG2XlNIIN8M7zQP5hkNJGqqQCIhNCayZAZGnPt\ncy8eCTS5pTmpWCKRaB+ASAiPPPnXitOWP/5MxWkiE5kKgEgIQyO5snc+MmBodLTWcURioQIgEsKb\nT5xfcdrZJ728hklE4qMCIBLCW151JHMPmDqu/aVzOjj1mJckkEgkOhUAkRDMjMsvfB375dOkgktB\nTPUsn/rH1yUdTaRqKgAiIWztG+B/XfUzRvpGadpWoGlbgaGtw3z0/9zKwKAuBy31KZYCYGZnmNkT\nZrbWzC4tM73ZzH4YTH/AzDrj6FekVu743RMUylwKolBw7n5gTQKJRKKLXADMLA18FTgTmA9caGZj\n95i9D3jB3Q8DvgJ8IWq/IrW0pXc7wyPBdf/diwMwMprnhd6BBJOJVC+ONYATgLXuvs7dR4AfAGeP\nmeds4KZg/FbgdWZW7qg6kQlpwfy5tGTTpAbzpAcKpAcKpAbzZCzFsfPnJB1PpCpxFIDZQOmZMOuD\ntrLzuHsO6AWmx9C3SE0cefhM6M+RKuy8FHSqAOmhPIcf3JF0PJGqxFEAyv2SH7uxNMw8xRnNFppZ\nt5l19/T0RA4nEod7fvskmUx6XHsmleL+B/+UQCKR6OIoAOuBuSXP5wAbKs1jZhlgCrCl3Ju5+yJ3\n73L3ro4O/bKSieGvG7YyNDT+jN/hkRwbnt2aQCKR6OIoAA8Bh5vZIWbWBFwALB4zz2Lg3cH4+cCd\n7l52DUBkIjr0kBmU22vlBefQTv1QkfoU+Wqg7p4zsw8DvwbSwI3uvsrMrgC63X0x8E3gO2a2luIv\n/wui9itSS1lSxY2W7rxYCYLfMNmUTqeR+hTL5aDdfQmwZEzb5SXjQ8Db4uhLJAlr1jxLajBHIW14\npviFb7kCaYw1Tz7Lccd1JhtQpAq6H4BICB0HTCY9nCc9dsNlCjoO2C+RTCJRad1VJIStz/fvvCUk\nJYe1FWB0cDihVCLRqACIhLD4loeAXY9n3jH+nRvvq3kekTioAIiEkMsVyp7MApAbydc0i0hcVABE\nQjjlNUeUP3MROPMtC2qaRSQuKgAiIXzoY2eQSad2KQIONDWl+aeLTkkqlkgkKgAiITQ1Zfj+4kv4\nu/kHgRVPBTjqmHncvPijpNP6M5L6pMNARUKatn87137jvTz60FOkUyle3tWJLmor9UwFQCSkxd/7\nPV/73C93Hg5q8NErz+P1bz0+2WAiVdK6q0gIz6zbxNeu/OUu17B1h6s/+WO29PQlF0wkAhUAkRC+\n8qkfF0fcoVAoDsG1gL7yqZ8kmEyketoEJBLCpo29kA+++EulUjy7vuyVzUUmPK0BiIRwyGEHjv/y\nBygUOOLouePbReqACoBICB0zJ7+4yWcX7kw/YHLtA4nEQAVAJIRn/rSp/I1NzVi3emPN84jEQQVA\nJIRD/m4W5SsAHDZ/dm3DiMREBUAkhLf/82tIpWzXzUDupNPGW997cnLBRCKIVADMbH8zW2pma4LH\naRXmy5vZimAYe79gkQlvxsypfOr6d5FtzpBKG6m00dyS5YobLmLylElJxxOpikW5N7uZfRHY4u5X\nmdmlwDR3/3iZ+frdvX1v37+rq8u7u7urzicSt+GhUVY+8CcsZRx1wktoataR1DKxmNlyd+8KM2/U\nT+/ZwGnB+E3A3cC4AiDyt6K5JUvXqS9LOoZILKLuAzjQ3TcCBI8HVJivxcy6zWyZmb11d29oZguD\nebt7enoixhMRkUr2uAZgZncAM8tMumwv+pnn7hvM7FDgTjNb6e5/Kjejuy8CFkFxE9Be9CEiInth\njwXA3U+vNM3MnjOzWe6+0cxmAZsqvMeG4HGdmd0NLADKFgAREamNqJuAFgPvDsbfDfx87AxmNs3M\nmoPxGcBJwOMR+xVJxOhIjkd++wQrf/8kudFc0nFEIom6E/gq4Edm9j7gaeBtAGbWBXzA3S8GjgD+\n08wKFAvOVe6uAiB1p/vOVXz+4kXsOHIunUlz+U0f5KhXvTThZCLViXQY6L6mw0BlotjyXC/v7bqM\n4cGRXdpb2pr57qNfoF3nAsgEsTeHgepMYJEQ7vnpQ7iXuRqoO7/9xcO1DyQSAxUAkRC2bd3OyPD4\nbf650Tz9WwcSSCQSnQqASAgLTjmCdCY9rt2BBaceUftAIjFQARAJoX3KJAq5/LiLwXkuT9t+LckF\nE4lABUAkhPuX/KH4c98dCjuHVMpYdtsjSccTqYoKgEgIqXQ6uBw0xSIQrAmkUilSKf0ZSX3SJ1ck\nhFeffTzpdPk/l1e9eUGN04jEQwVAJITZLzmQi/73eTS1ZGlqydLc2kRTS5YPffkdzDio7G0wRCY8\nXcxcJKS3/vPpnPSm47j/thWkUsaJb1rA9JlTk44lUjUVAJG90DFnf97y/tcmHUMkFtoEJCLSoFQA\nREQalAqAiEiDUgEQEWlQKgAiIg1KRwGJ7IVCocC6R57GUsahR8/DzJKOJFK1SGsAZvY2M1tlZoXg\nLmCV5jvDzJ4ws7VmdmmUPkWSsuLuVZwz42I+dOJl/MsrPsm5B76f1Q+sSTqWSNWibgJ6DDgXuLfS\nDGaWBr4KnAnMBy40s/kR+xWpqRc2beXjZ3yOwf4hvOB4wdm+dYB/e81nGdg2mHQ8kapEKgDuvtrd\nn9jDbCcAa919nbuPAD8Azo7Sr0itffc/foIXxt8+NZ8r8KMv/yKBRCLR1WIn8GzgmZLn64M2kbqx\n7tGnK05bu+LPtQsiEqM9FgAzu8PMHiszhP0VX24vWcU70ZvZQjPrNrPunp6ekF2I7FuHHdtZcdrh\nxx1auyAiMdpjAXD309395WWGn4fsYz0wt+T5HGDDbvpb5O5d7t7V0dERsguRfeufPn0eqTKXg05n\n07z9Y29KIJFIdLXYBPQQcLiZHWJmTcAFwOIa9CsSmynTJ/OlpZfRNmUSljIsZew3vZ1r7v0sre2t\nSccTqUqk8wDM7BzgOqAD+JWZrXD3N5rZQcAN7n6Wu+fM7MPAr4E0cKO7r4qcXKTGjjr5CH6y6Rv8\nedV6Uilj3hGzdR6A1DVzr7g5PnFdXV3e3d2ddAwRkbphZsvdveJ5WaV0KQgRkQalAiCyl/L5Avl8\nIekYIpHpWkAiIW16+nmu+cAilt+xEjM44cwFXPK19zN9lu4JLPVJawAiIQwPjvA/X3UZy+9YSSFf\nIJ8r8OBtf+CSkz9NbjSXdDyRqqgAiIRw763LGNw2SKFk008+V6Bv8zYe+NUfEkwmUj0VAJEQ1j+5\ngcH+4XHtI4OjrF9T8bxGkQlNBUAkhENePo/W9pZx7U2tWTqPnJdAIpHoVABEQjjpnBOY0rEf6Wz6\nxbZMU4YD5s6g643HJJhMpHoqACIhZJsyXHf/lZz2thNpasnS3NrE695xMl+597Oky1wjSKQe6JMr\nEtKah9ex7JfdZLIp0hnj9z9/kKdWVr5MtMhEpwIgEsILz23ls+d9ie29Awz0DTLQN8i2zf186s2f\nZ3vv9qTjiVRFBUAkhLtu/l3ZO4K5O/feuiyBRCLRqQCIhNC3ZRsjw6Pj2nMjOfpf0BqA1CcVAJEQ\njjv9aFomNY9rT2fTLDj9qAQSiUSnAiASwlGvPoLjX380LW07i0BLWzOnnH8ihx17SILJRKqni8GJ\nhGBmfPqWj3HvLcu4/aa7SKVSnHHRazn53FckHU2karohjIjI35Ca3RDGzN5mZqvMrGBmFTs0sz+b\n2UozW2Fm+kYXEZkAom4Cegw4F/jPEPO+xt2fj9ifiIjEJFIBcPfVgG6MLSJSh2p1FJADt5vZcjNb\nuLsZzWyhmXWbWXdPT0+N4omINJ49rgGY2R3AzDKTLnP3n4fs5yR332BmBwBLzeyP7n5vuRndfRGw\nCIo7gUO+v4iI7KU9FgB3Pz1qJ+6+IXjcZGY/BU4AyhYAERGpjX1+HoCZtQEpd98WjL8BuCLMa5cv\nX/68mf1lnwaMZgZQDzu26yUn1E9W5YxfvWSd6DkPDjtjpPMAzOwc4DqgA9gKrHD3N5rZQcAN7n6W\nmR0K/DR4SQb4vrtfWXWnE4iZdYc93jZJ9ZIT6iercsavXrLWS84woh4F9FN2frmXtm8AzgrG1wG6\nZZKIyASjawGJiDQoFYBoFiUdIKR6yQn1k1U541cvWesl5x5N6GsBiYjIvqM1ABGRBtXQBcDMzjCz\nJ8xsrZldWmb6wWb2GzN71MzuNrM5JdO+YGaPBcM/lLR/L3jPx8zsRjPLBu2nmVlvcEG8FWZ2ecI5\nv2VmT5XkOTZoNzO7NujrUTM7LmzOfZj1vpKcG8zsZ0F7lGV6o5ltMrPHKkyvuBzM7N1mtiYY3l3S\nfnxw0cO1wWstaN/fzJYG8y81s2lJ5TSzSWb2KzP7oxUv5HhVyfzvMbOekuV5cdic+yJr0H538Hna\nkemAoL3ZzH4YvNcDZtaZVE4zm1ySb4WZPW9m1wTTIi3Tfc7dG3IA0sCfgEOBJuARYP6YeW4B3h2M\nvxb4TjD+JmApxaOo2oBuYL9g2lmABcPNwAeD9tOAX06gnN8Czi/T31nAbUH+VwIPJJ11zOt/DLwr\nyjINXnsKcBzwWIXpZZcDsD+wLnicFoxPC6Y9CJwYvOY24Myg/YvApcH4pcAXksoJTKJ4YUaCf6P7\nSnK+B7g+wt/UvlimdwNdZd7rX4CvB+MXAD9MMueY1y8HToljme7roZHXAE4A1rr7OncfAX4AnD1m\nnvnAb4Lxu0qmzwfucfecu2+n+EV3BoC7L/EAxS+EOUSzT3LuxtnAt4P/hWXAVDObNRGymtlkikXj\nZyHzVOTFS5Fs2c0slZbDG4Gl7r7F3V+gWLTOCKbt5+73B//23wbeWvJeNwXjN5W01zynuw+4+13B\ne48ADxP9M7pPsu6hu9Jleivwuh1rXEnmNLPDgQMoFtYJr5ELwGzgmZLn64O2Uo8A5wXj5wCTzWx6\n0H5msDo9A3gNMLf0hVbc9PNO4L9Lmk80s0fM7DYzO3IC5LwyWMX9ipntuNdhmP6SyLpj/t+4e19J\nWzXLNIxK/y+7a19fph3gQHffCBA8HpBgzheZ2VTgf7CzIAOcF3wmbjWzscs/qaz/FWw++XTJl/yL\nr3H3HNALTE84J8CFFNdGSo+u2ZfLNJJGLgDlfi2MPSTq34FTzewPwKnAX4Gcu98OLAF+T3Ezz/1A\nbsxr/x9wr7vv+CXwMHCwux9D8ezpsL9i91XOTwAvA/6e4irtx/eiv1pn3eHCYNoO1S7TMCr9v+xt\n+75WVR4zy1Bcltd68WRNgF8Ane5+NHAHO39hJ5n1He5+FPDqYHjnHt4rDlH+jS9g18/ovl6mkTRy\nAVjPrr8w5wAbSmdw9w3ufq67LwAuC9p6g8cr3f1Yd389xQ/Gmh2vM7PPULw8xr+VvFefu/cH40uA\nbPBLN5Gc7r4xWMUdBv6L4uabUP3VOitAsJZwAvCrkveqdplG+X/ZXfucMu0Az+3YjBY8boopYzU5\nd1gErHH3a3Y0uPvm4PMA8A3g+BhzVpXV3f8aPG4Dvk+Zz2lQzKaw+806+zRnkOMYIOPuy3e01WCZ\nRrO7HQR/ywPFnY3rgEPYucPyyDHzzKB4ITuAK4ErgvE0MD0YP5rindEywfOLKf6KbR3zXjPZed7F\nCcDTO54nlHNW8GjANcBVwfM3sesOsAeTXqZB2weAm+JYpiWv76TyjsCyy4Hi2tJTFHcCTgvG9w+m\nPRTMu2Mn8FlB+5fYdSfwF/fysxp3zv+guDM9Nea9ZpWMnwMsq+LvKraswedpRjBPluK2/g8Ezz/E\nrjuBf5R18hkVAAABE0lEQVTkMg2mXwV8Nu5lui+HxAMk+j9f3Nv/JMUjVy4L2q4A3hKMn0/xV+iT\nwA1Ac9DeAjweDMuAY0veMxe834pguDxo/zCwiuKX4jLgVQnnvBNYSfGL9rtAe9BuwFeDvlZS5giM\nWmcNpt9NcSdmaVuUZXozsBEYpfjL7n0Ui8wH9rQcgIuAtcHw3pL2rmB5/gm4np3FaTrF7exrgsf9\nk8pJ8VerA6tLPqMXB9M+X7I87wJetpf/9nFnbaN4RM2jQa7/C6RLPi+3BPM/CBya5L99MG3d2GUW\ndZnu60FnAouINKhG3gcgItLQVABERBqUCoCISINSARARaVAqACIiDUoFQESkQakAiIg0KBUAEZEG\n9f8BERicR2aFA+4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4977fd2898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## ADD CODE HERE\n",
    "# Fit and plot the logistic regression\n",
    "## ADD CODE HERE\n",
    "# Fit and plot the logistic regression\n",
    "p=predict(x,theat2)\n",
    "print ('ac,curecy= {} %'.format(100*np.sum(np.byte(p >= 0.5) == y.reshape(N,))/N))\n",
    "plt.scatter(x[:,0],x[:,1],c=p)\n",
    "# plt.plot(y,'o',label='logistic')\n",
    "# plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There you go! That should be a nice and easy fit. There are a few things we can start playing with at this point:\n",
    "\n",
    "* Evaluate the performance of your model: plot the decision boundary, likelihood and accuracy on held-out test sets, etc.\n",
    "* Write a gradient-based and a non-gradient-based optimizer. Do they arrive at the same result? Which one takes longer? Which one evaluates the likelihood function more times?\n",
    "\n",
    "(Warning: if the plot looks odd and you get several warnings, it may be that the data is linearly separable and the sigmoid is saturating, leading to `np.log(0)` numerical problems. Add more noise and retry.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian logistic regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MAP estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's move to Bayesian inference on the parameters $\\boldsymbol\\theta$. Let's put a prior on them. Because that's what we do. We put priors on things.\n",
    "\n",
    "More specifically, let's use a Gaussian prior parametrized by a mean $\\boldsymbol m$ and a variance $\\boldsymbol S$:\n",
    "\n",
    "$$\\boldsymbol\\theta \\sim  \\mathcal{N}(\\boldsymbol m, \\boldsymbol S)$$\n",
    "\n",
    "Given that $\\boldsymbol\\theta_{\\mathrm{ML}}$ had no analytical solution, it should really come as no surprise that $\\boldsymbol\\theta_{\\mathrm{MAP}}$ doesn't either. That should be no problem for a machine learning expert like you:\n",
    "\n",
    "**Task 4**\n",
    "\n",
    "1. Write down the equation for the full unnormalized posterior $p(\\boldsymbol\\theta | \\boldsymbol X, \\boldsymbol y) \\propto p(\\boldsymbol y | \\boldsymbol\\theta, \\boldsymbol X) p(\\boldsymbol\\theta)$.\n",
    "2. Write a separate function for it, as we did with the log-likelihood above.\n",
    "3. Optimize it to find $\\boldsymbol\\theta_{\\mathrm{MAP}}$ and use it to make predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## EDIT THIS FUNCTION\n",
    "def postior_log(X, y,m,S, theta):\n",
    "    N, D = X.shape\n",
    "    theta=theta.reshape(-1,1)\n",
    "    log_theta=(theta-m).T @ inv(S) @ (theta-m)\n",
    "#     print('log_theta.shape',log_theta.shape)\n",
    "    theta_map = np.zeros((D,1)) ## <-- EDIT THIS LINE\n",
    "    pl= log_theta + log_likelihood(X,y,theta)\n",
    "    return pl\n",
    "def map_estimate(X, y, m, S):\n",
    "    # X: N x D matrix of training inputs\n",
    "    # y: N x 1 vector of training targets/observations\n",
    "    # m: D x 1 prior mean of parameters\n",
    "    # S: D x D prior covariance of parameters\n",
    "    # returns: maximum a posteriori parameters (D x 1)\n",
    "    \n",
    "    N, D = X.shape\n",
    "    theta_map = np.zeros((D,1)) ## <-- EDIT THIS LINE\n",
    "    fun = lambda w: -postior_log(X,y,m,S,w)\n",
    "    theta_map=scipy.optimize.minimize(fun,np.random.normal(size=(D,1)),method='TNC',options={'gtol': 1e-100, 'disp': True})\n",
    "    return theta_map.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you can perform a similar model evaluation as you did before. How does your prior influence the MAP estimate and the model's performance?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "th_map [-1.99136229  4.17167365]\n",
      "accurecy= 98.0 %\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f4975e007f0>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGzFJREFUeJzt3X1wJHd95/H3d0bSSvvkfZC0u/audzF+wjZ+wLIxcODg\nmGMNuXMMJAeVws4FysXdpepSl6sKKe64InXckaQqIRy5Sm3ubANHQgIECGBjL9gGUsEG2fHDGgd7\n/byrXa/2SdLoebq/90e37JE0o/XOtLrV05+Xa0qj/rV6Pu6VvvrpO7/pMXdHRETaRynrACIikiwV\ndhGRNqPCLiLSZlTYRUTajAq7iEibUWEXEWkzKuwiIm1GhV1EpM2osIuItJmOLB60t7fXd+3alcVD\ni4jk1kMPPXTU3ftOtV8mhX3Xrl0MDg5m8dAiIrllZi+8lv3UihERaTMq7CIibUaFXUSkzaiwi4i0\nGRV2EZE2k8mqGJGVxmf34eO3QfASdF2Drb4FK/dmHUukKSrsUnjh5F4Y+V1gBghh9kl84m+h91tY\neWvW8UROm1oxUmjuAYx9EpgCwnjrDPgoXvl8hslEmqfCLsUWDIFP1BuA6R+nHkckCSrsUmyldeDV\nBmMb0s0ikhAVdik0K22AjgvqD/a8P90wIglRYZdCcw8hOFhnxGD2qdTziCRBhV2KLTgIPllnwGFG\nPXbJp5YLu5l1m9lPzexRM3vCzD6VRDCRVJTWAkGDsTNSjSKSlCRm7NPAde5+GXA5sNvMrknguCLL\nzkob1WOXttPyC5Tc3YFK/GlnfPNWjyuShqjHPlRnxGD2mdTziCQhkR67mZXN7BHgCLDX3R+ss8+t\nZjZoZoPDw8NJPKxI65bssf8w9TgiSUiksLt74O6XA9uBq83skjr77HH3AXcf6Os75Ts7iaRjyR77\n+lSjiCQl0VUx7n4SuB/YneRxRZZL1GM/r/5gzwfSDSOSkCRWxfSZ2Yb4fg9wPfDPrR5XJA1Rj/1w\nnRGD2edSzyOShCRm7NuA+8zsMeBnRD327yRwXJHlFxwEH68z4DBzX+pxRJKQxKqYx4ArEsgikr7S\nGl69quPCMfXYJZ/0ylMpNCttgo5z6w92ax275JMKuxRa1GM/UmfEoPpC6nlEkqDCLsUWDIFX6gw4\nzNybehyRJKiwS7GVVtO4x74u1SgiSVFhl0JTj13akQq7FFrUY693iQuD6kup5xFJggq7FFswBD5W\nZ0A9dskvFXYptiV77GtSjSKSFBV2KbSox/76+oPdulaM5JMKuxSauzfusQcHUs8jkgQVdim24GDj\ndezTP0g9jkgSVNil2NRjlzakwi6FFvXYz6k3At3vSz2PSBJU2KXQoh77sfqDda/TLrLyqbBLsQUH\nG69jn96behyRJKiwS7FZDw177LY61SgiSVFhl0Kz8mboeF29EehRj13ySYVdCi3qsZ+sP1j3Ou0i\nK58KuxRbOAQ+UmfAYfqe1OOIJEGFXQqum8Y99p5Uk4gkRYVdCi3qse+qN6Ieu+SWCrsUWtRjH60/\nWPcaMiIrnwq7FFs4BF7vyVOtY5f8UmGXgluqx96dahKRpKiwS6FZeTOUd9UbgZ6b0o4jkggVdik0\nd4ewUY+9wTVkRFY4FXYptnAI/ESdAa1jl/xSYZeCWwV4/SFblWoSkaS0XNjNbIeZ3WdmT5rZE2b2\nH5MIJpIGK/dCeWe9EfXYJbc6EjhGFfhdd3/YzNYBD5nZXnf/eQLHFllW7t7gsr1A2OAaMiIrXMsz\ndnc/5O4Px/fHgCeBs1o9rkgqwkMNCrjD1PdSjyOShER77Ga2C7gCeDDJ44osn04ar2NXj13yKbHC\nbmZrga8Dv+Pui9aPmdmtZjZoZoPDw3qptqwMVu5r3GPv/tXU84gkIZHCbmadREX9y+7+d/X2cfc9\n7j7g7gN9fX1JPKxIy6Ie+3iDwQbr20VWuCRWxRjwf4En3f1PWo8kkqLwEIT1XojkMHVX6nFEkpDE\njP1twIeB68zskfj2ngSOK5KCThqvY+9KNYlIUlpe7uju/wBYAllEUmflPry8A4JnFo6oxy65pVee\nSqFFPfbJBoOVdMOIJESFXYotPKweu7QdFXYpuA4a99g7U00ikhQVdim0aB37jkWl3THovjGTTCKt\nSuJaMSK55e6EPonhi4t7o/XtIiucZuxSbOEhCIfxRf+F+NR3s04n0hQVdim4pXrs+oNW8kmFXQrN\nyv1Q3l63x27qsUtOaUoihebuhOE05gt77Ib7VEapRFqjGbsUW3gYD48QRl31ue46IQHB1HeyTifS\nFBV2KbgytT322lm76Q9aySkVdik09dilHWlKIoXnPgOL1rEbzkw2gURapBm7FJoHh/DwCA4LbiHh\n1LezDSfSJBV2KbgSDdex68dDckrfuVJoVt4C5bPq9thLuh675JR67FJ4gc9SWrCO3ePtInmkGbsU\nWhgcJgwOUyUkiNewB9EqdqpT38o6nkhTNGMXiXmdlTEieaQZuxRaqbwVK59Zt8fe2aMeu+STZuxS\neLMElFncY68SZBVJpCWasUuhBcEhgmCIGUKqOAHOLCGzBExPfiPreCJN0YxdJBYuasioxy75pBm7\nFFq5vI1Sgx77qp6bMskk0irN2KXwpsKQzgXr2EOMWW/0ilSRlU0zdim0anCYajDEhDsz7sy6M+XO\nlAeMT/xd1vFEmqIZuxTcq7Py6hJjInmiGbsUWkd5G6Xytro99jWr35dJJpFWacYuhTfuXXQtqOwB\nxpR3ZhNIpEWJzNjN7DYzO2Jm+5I4nkhaZqqHmawOcdI7GPcyk16i4mVGvcTxCa1jl3xKqhVzB7A7\noWOJpGhuqm7MUGKKMrOUiN5BKcwymEjTEins7v4j4HgSxxJJU1fHNjrLW1m4stEpsXG1rhUj+aQn\nT6XwTvpaHAjj4h46VL3EhK/ONJdIs1Ir7GZ2q5kNmtng8PBwWg8rsqSp6jCV2ZcZDtcx5t1Uwi5G\nwh6OhT0crHwn63giTUmtsLv7HncfcPeBvr6+tB5WZElOABYtb5z0LirezTSdqMcueaZWjBRaT8dW\nusr9dXvs29bckE0okRYltdzxr4GfABeY2QEz+0gSxxVJw3C1jxCr6bEbM2GZsXBjtsFEmpTIC5Tc\n/UNJHEckbePVYwzPHMW9j3WlKTosYMo7GQ9XMTpyL1f13px1RJHTpleeSqGFHgBGSImRcP4qmLDO\n1WNE8kA9dim0tR19rC5vqjNS4vx170w9j0gSVNil0MyMwzNbCNwIPHrHpMCNqaDMaLA543QizVEr\nRgrt5MxJDkyexH0rZ3RO0GkBE0EXY0E3laM/4YZtevWp5I8KuxRaSIBhVClxfHbtvLHAg4xSibRG\nrRgptI2dm1jXeUadEWNg45tTzyOSBBV2KTQz48hEH0H46jr2wI3poIPKrNaxSz6pFSOFdmx6jP2V\nCqH3s6Frgs5SwES1i9HZbqZmH+Hmc67POqLIaVNhl0KrekAJmPUSx6bn99hnQ/XYJZ/UipFC6191\nBhu71i7abhjXbXljBolEWqfCLoVmZkxXNhCGEMYXcwxDqFZLlGbrPakqsvKpsEuhHZmssP/4GMPD\n66lUuhkf72JkdDVHjq7la/ufyDqeSFPUY5dCmw0CzIwwLDE+0T1vbCZQj13ySTN2KbQz16xn86rF\nb4FXwnjvzgszSCTSOhV2KTQz44yJ1RDCK2+YFIJXYUNV73kq+aTCLoV2dHycF18epeNIF6VKGRsv\nUR7poHykk28+9mTW8USaoh67FNpMNaBkhgVGuTL/x2FaPXbJKc3YpdC2rV/Hxp6eRdtLwO4Lzk0/\nkEgCVNil0MyMTZNdi3rsVKHPFxd8kTxQYZdCOzY6zsGhk/QcgY4KlMehawRWHYFvP6Aeu+STeuxS\naNPVaB27hdBVmT82NTubTSiRFmnGLoW2beM6Nq7pXrTdgOsvOy/9QCIJUGGXQjMzNs90Qejg8QXZ\nQ8cC2FJWj13ySYVdCu3EyARDL5xg9bDTWXE6JpxVo07PkZDv3atrxUg+qccuhTY9U6VUMkpVWFUB\n8FfGpmaqmeUSaYVm7FJoW3rXsX5tnR67wTuu0jp2yScVdik0M2OzdUT99ZoeO4Gzpc4Ll0TyQIVd\nCu3EiXEOPnOUVcdnKE8ElKYCOsardB2b4Z67Hss6nkhT1GOXQpuZjnrsNgudE/OvDTM1pXXskk+J\nzNjNbLeZ/cLM9pvZx5M4pkga+resZ22DHvu/ePv5GSQSaV3Lhd3MysCfAzcAFwEfMrOLWj2uSBrM\njE1dC3rsHvXY+9erxy75lMSM/Wpgv7s/6+4zwFeAGxM4rsiyO3l8nINPvUzHyBQ2VcWmq5QmZimf\nnOSebzycdTyRpiTRYz8LeKnm8wPAmxM4rsiym5qcoVQuYbMBHZPze+pTE9MZpRJpTRIzdquzzRft\nZHarmQ2a2eDw8HACDyvSui1nbmDNuvo99rdep46i5FMShf0AsKPm8+3A0MKd3H2Puw+4+0BfX18C\nDyvSOjNjw+rORT12D52+TWuyDSfSpCQK+8+A88zsdWbWBXwQ+PsEjiuy7E4eq3DgF0NQmYSZWZit\nwtQMjE2w96sPZh1PpCkt99jdvWpmvw3cDZSB29xdV0+SXJiemKFUKmEewPT8HvtERT12yadEXqDk\n7ncCdyZxLJE09W/fyJr1PUwveDGSmXHNuy7OKJVIa3RJASk0M2Pjhu7FPXYP2XLmhmzDiTRJhV0K\nbeRYhRefHILpGagGEARRn31ymr1/9Y9ZxxNpiq4VI4U2NT5NqWww41Cdf/31ibGpjFKJtEYzdim0\nvu0bWb1u8aUDrGRcs/vSDBKJtE6FXQqtVCqxcfOauK8e9djdHQ9Dtm7fmHE6keaosEuhjRwb48Un\nDuAzsxAEeBBANcCnZrjnS/+QdTyRpqjHLoU2OTZFqVyCmWr05GmN8dGJjFKJtEYzdim0/rM301Pv\neuwl4827L8sgkUjrVNil0EqlEpu3rgcW9tids87Zkm04kSapsEuhjRwd44WfH8CrQXzxrxDCEK9W\nufuLP8w6nkhT1GOXQpusTFEqxfObMJw3Vjk5nkEikdZpxi6FtnSP/fIMEom0ToVdCq1UKtE7d02Y\nBddj337e1uyCibRAhV0KbfTYGM/veym6RszchcDcIQi45wvqsUs+qbBLoU2MTlLuiH8M3KM+ezxz\nHzuhHrvkkwq7FFr/zl6619TvsV/1bq1jl3xSYZdCK5VK9MXXhFm4jn3nG87MMppI01TYpdBGj1d4\n/vEXo2vEzF0IzB0PAu6+4/6s44k0RevYpdDGRyaia8XA/HdRIir6InmkGbsU2padvXSvXrVou3rs\nkmcq7FJopVKJ/rP7gMU99l0X78gymkjTVNil0EaPV3jusRfw8NWijoOHzvduvz/bcCJNUmGXQhsf\nmaDUUY4+8fgWGz02lkkmkVapsEuh9Z/dS/fqrkXbzYwr36X3PJV8UmGXQiuXS/Tv7Fu03d15/WW7\n0g8kkgAVdim00eMVnt/3Yt2xu++4L+U0IslQYZdCGz85Tqlcrjs2clQ9dsknFXYptP6dfazqqd9j\nf9P16rFLPqmwS6GVyyW27uqNPqm9Hrs7579pV2a5RFrRUmE3s18zsyfMLDSzgaRCiaRl7ESF5x5/\nATx+W7xXinvIXbfdm10wkRa0OmPfB7wP+FECWURSVzkxTnluHfsrC9mj4j5yZDSrWCItaekiYO7+\nJET9SJE86t/ZS+eqTqbGp+dtNzOuuP6NGaUSaY167FJo5XKZM1+/+L1N3Z3zB87NIJFI605Z2M3s\n+2a2r87txtN5IDO71cwGzWxweHi4+cQiCRo7UeHZx56vO3b37eqxSz6dshXj7tcn8UDuvgfYAzAw\nMOCn2F0kFXM99tnp6qKxky+PZJBIpHVqxUih9Z8d9dgXMoPLr7skg0QirWt1ueNNZnYAeAvwXTO7\nO5lYIukod5Q589x6PXa48Gr12CWfWirs7v4Nd9/u7qvcfYu7vzupYCJpqJwc59lHn687pvc8lbxS\nK0YKbex4pWYd+3wnDp9MOY1IMlTYpdD6z+6lo2vxGgIzuPyd6rFLPqmwS6GVO8psP3/bou3ucOE1\n52WQSKR1KuxSaJWT4zzzyPN1x+7R9dglp1TYpdCW6rEfO6Qeu+STCrsUWqMeOwaXXXtR+oFEEqDC\nLoUW9djPXDzgcPHbLkw/kEgCVNil0MZHxnnmkefqjqnHLnmlwi6FNnqscY/96NCJlNOIJEOFXQqt\n/+xeyp2LC7uZcek73pBBIpHWqbBLoZU7ypx9wVmLtrs7l7xdhV3ySYVdCm18ZJz9Ddex359qFpGk\nqLBLoUU99vo/BkcPHk85jUgyVNil0Pp2bG7YY3/j27XcUfJJhV0KraOzgx0NeuyXXXtxBolEWqfC\nLoUWrWN/vu7Y3bdrHbvkkwq7FNpSPfZh9dglp1TYpdD6dmymVF78Y2BmXPK2CzJIJNI6FXYptI7O\nDnZetH3RdnfXG21IbqmwS6GNj06w/5+erzum9zyVvFJhl0IbGR5t3GM/cCzlNCLJUGGXQus/u7dh\nj/3it6rHLvmkwi6F1tHZwc6Ldyza7u5ccZ167JJPKuxSaFGPvcH12L94f7phRBKiwi6FNjI8SrlO\nKwbgyIvqsUs+qbBLoS21jv0N15yXQSKR1qmwS6F1dnWyq0GP/cp3XZZBIpHWqbBLoU2MTfJ0gx77\nXvXYJadU2KXQTh4Zadhjf/n54ZTTiCSjpcJuZn9sZv9sZo+Z2TfMbENSwUTS0LdjM1ayRdvVY5c8\na3XGvhe4xN0vBZ4Cfr/1SCLp6ezq5HWXnL1ou7tz5b+8PINEIq1rqbC7+z3uXo0/fQBYfDUlkRVs\nYmyy4Tr273/phymnEUlGkj323wLuajRoZrea2aCZDQ4Pq3cpK8PJIyOUOha/NR7AYfXYJac6TrWD\nmX0f2Fpn6BPu/q14n08AVeDLjY7j7nuAPQADAwPeVFqRhPXt2EzJ6vfYL7zq3AwSibTulIXd3a9f\natzMbgF+Bfhld1fBllzp7Ork5k/9Orf/l68wPTENgBmsWt3FzZ/69YzTiTTnlIV9KWa2G/g94Fp3\nn0gmkki63v87v0LvmZv48qe/ztGh47zhzefxkf/xG3WfVBXJA2tlkm1m+4FVwNxFNR5w94+d6usG\nBgZ8cHCw6ccVESkiM3vI3QdOtV9LM3Z3VxNSRGSF0StPRUTajAq7iEibUWEXEWkzKuwiIm1GhV1E\npM20tNyx6Qc1GwZeSP2BX7te4GjWIV6DvOSE/GTNS07IT9a85ISVn3Wnu/edaqdMCvtKZ2aDr2Wt\naNbykhPykzUvOSE/WfOSE/KVdSlqxYiItBkVdhGRNqPCXt+erAO8RnnJCfnJmpeckJ+seckJ+cra\nkHrsIiJtRjN2EZE205aF3cx2m9kvzGy/mX28zvhOM/tB/Cbc95vZ9pqxPzSzffHt39Rs/3J8zH1m\ndpuZdcbbf8nMRszskfj2yRWQ9Q4ze64m0+XxdjOzz8WP9ZiZvSnjnD+uyThkZt+Mtzd9TuN/myNm\ntq/BeMNzYGa3mNnT8e2Wmu1Xmtnj8dd8zix6Zw4z22Rme+P995rZxqxymtlqM/uuRW8u/4SZfaZm\n/980s+Ga8/nR15pzObLG2++Pv5/mMvXH21eZ2d/Ex3rQzHZlldPM1tXke8TMjprZZ+Oxls7psnP3\ntroBZeAZ4BygC3gUuGjBPl8FbonvXwd8Kb7/XqI36O4A1gCDwPp47D2Axbe/Bv5dvP2XgO+ssKx3\nAB+o83jvIXr7QgOuAR7MMueCr/86cHMC5/QdwJuAfQ3G654DYBPwbPxxY3x/Yzz2U+At8dfcBdwQ\nb/8j4OPx/Y8Df5hVTmA18M54ny7gxzU5fxP4fAs/U8txTu8HBuoc698DfxHf/yDwN1nmXPD1DwHv\nSOKcLvetHWfsVwP73f1Zd58BvgLcuGCfi4AfxPfvqxm/CPihu1fdfZyogO0GcPc7PUb0g57EG3cv\nS9Yl3Ah8Mf7feADYYGbbss5pZuuIfhl88zVkWZK7/wg4vsQujc7Bu4G97n7c3U8Q/TLaHY+td/ef\nxP/2XwR+teZYX4jvf6Fme+o53X3C3e+Ljz0DPExCby6fdNZTPFztOf0a8MtzfyFlmdPMzgP6iX5h\nrnjtWNjPAl6q+fxAvK3Wo8D74/s3AevMbHO8/Yb4z9pe4J3AjtovtKgF82HgezWb32Jmj5rZXWZ2\n8QrJ+un4z80/NbNVp/F4aeec2/8H7j5as63Zc3oqjf5fltp+oM52gC3ufggg/tifYc5XmNkG4F/x\n6i9agPfH3w9fM7OF5z+rrLfHbYz/WlO8X/kad68CI8DmjHMCfIjor4fa1SbLeU5b0o6Fvd5v94VL\nf/4zcK2Z/RNwLXAQqLr7PcCdwD8StVt+QvQm3bX+N/Ajd5/7zf0w0ct8LwP+F6c361yurL8PXAhc\nRfTn5e+dxuOlmXPOh+KxOa2c01Np9P9yutuXW1N5zKyD6Fx+zt2fjTd/G9jl7pcC3+fVGXGWWX/D\n3d8IvD2+ffgUx0pCK//GH2T+9+hyn9OWtGNhP8D8GeF2YKh2B3cfcvf3ufsVwCfibSPxx0+7++Xu\n/i6if/Cn577OzP4b0Af8p5pjjbp7Jb5/J9AZz0wzy+ruh+I/N6eB24laKa/p8dLMCRDP6q8Gvltz\nrFbOabP/L0tt315nO8DLc62s+OORhDI2k3POHuBpd//s3AZ3PxZ/LwD8JXBlgjmbyuruB+OPY8Bf\nUed7NP4ldQZLt1eWNWec4zKgw90fmtuWwjltzVIN+DzeiJ6kexZ4Ha8+0Xfxgn16gVJ8/9PAH8T3\ny8Dm+P6lwD6if1CAjxLNOnsWHGsrr74e4GrgxbnPM8y6Lf5owGeBz8Sfv5f5Tx79NMuc8baPAV9I\n6pzGX7OLxk+g1T0HRH/ZPEf05NnG+P6meOxn8b5zT56+J97+x8x/8vSPTvN7Nemc/53oSejSgmNt\nq7l/E9F7E5/uz1ViWePvp954n06iXvrH4s//A/OfPP3bLM9pPP4Z4FNJn9PlvGUeYFn+p6Jnv58i\nWsnxiXjbHwD/Or7/AaJZ41PA/wFWxdu7gZ/HtweAy2uOWY2P90h8+2S8/beBJ4iK3QPAW1dA1nuB\nx4mK6P8D1sbbDfjz+LEep86qhDRzxuP3Ez35V7ut6XNK9OfyIWCWaCb2EaJfHh871TkAfgvYH9/+\nbc32gfhcPgN8nld/6Wwm6mM/HX/clFVOolmmA0/WfI9+NB77nzXn8z7gwtP8Hk066xqiFSaPxbn+\nDCjXfL98Nd7/p8A5Wf7bx2PPLjxnrZ7T5b7placiIm2mHXvsIiKFpsIuItJmVNhFRNqMCruISJtR\nYRcRaTMq7CIibUaFXUSkzaiwi4i0mf8PPzOvhiRFYy0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4975e9e8d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## ADD CODE HERE\n",
    "# Fit and plot the MAP logistic regression estimate\n",
    "# np.random.seed(50)\n",
    "m = np.zeros((D, 1))\n",
    "S = 5*np.eye(D)\n",
    "th_map=map_estimate(x, y,m,S)\n",
    "print('th_map',th_map)\n",
    "p=predict(x_test,th_map)\n",
    "print ('accurecy= {} %'.format(100*np.sum(np.byte(p >= 0.5) == y.reshape(N,))/N))\n",
    "#print ('p',p)\n",
    "plt.scatter(x_test[:,0],x_test[:,1],c=p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Laplace approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we have hinted above, in logistic regression the posterior distribution over $\\boldsymbol\\theta$ doesn't have an analytical solution. This is the first example in the course of _approximate Bayesian inference_: The exact posterior is analytically intractable so that we have to approximate it using one of various techniques. The one we'll use in this part of the tutorial is called the **Laplace approximation**.\n",
    "\n",
    "In brief, **the Laplace approximation is a Gaussian centered at the peak of the pdf of interest with the same curvature**. Let's make this a bit more rigorous below.\n",
    "\n",
    "Let's say we have a probability distribution $p(\\boldsymbol z)$ we want to approximate. The distribution $p(\\boldsymbol z)$ is of the form\n",
    "\n",
    "$$p(\\boldsymbol z) = \\frac{1}{Z} \\tilde{p}(\\boldsymbol z) ~ ,$$\n",
    "\n",
    "where $\\tilde{p}(\\boldsymbol z)$ is an unnormalized distribution that we can evaluate easily, but $Z$ is unknown. Formally, the Laplace approximation results from a second-order Taylor expansion of $\\log \\tilde{p}(\\boldsymbol z)$ around $\\boldsymbol  z_0$:\n",
    "\n",
    "$$\\log \\tilde{p}(\\boldsymbol z) \\approx \\log \\tilde{p}(\\boldsymbol z_0) + \\frac{d}{d\\boldsymbol z}\\log \\tilde{p}(\\boldsymbol z)\\Big|_{\\boldsymbol z=\\boldsymbol z_0}(\\boldsymbol z -\\boldsymbol  z_0) + \\frac{1}{2}(\\boldsymbol z-\\boldsymbol z_0)^\\top\\frac{d^2}{d\\boldsymbol z^2} \\log \\tilde{p}(\\boldsymbol z)\\Big|_{\\boldsymbol z=\\boldsymbol z_0}(\\boldsymbol z-\\boldsymbol z_0)$$\n",
    "\n",
    "Now let's evaluate this expression at the mode of $p(\\boldsymbol z)$ &ndash; which is the same as the mode of $\\tilde{p}(\\boldsymbol z)$. We define the mode $\\boldsymbol z^*$ such that\n",
    "\n",
    "$$\\frac{d}{d\\boldsymbol z} \\tilde{p}(\\boldsymbol z) \\Big|_{\\boldsymbol z = \\boldsymbol z^*} = \\boldsymbol 0 ~ .$$\n",
    "\n",
    "At this point, the $\\mathcal{O}(\\boldsymbol z)$ term of the expansion vanishes and we are left with\n",
    "\n",
    "$$\\log \\tilde{p}(\\boldsymbol z) \\approx \\log \\tilde{p}(\\boldsymbol z^*) - \\frac{1}{2}(\\boldsymbol z-\\boldsymbol z^*)^\\top\\boldsymbol A(\\boldsymbol z-\\boldsymbol z^*)$$\n",
    "\n",
    "Or, equivalently,\n",
    "\n",
    "$$\\tilde{p}(\\boldsymbol z) \\approx \\tilde{p}(\\boldsymbol z^*) \\exp\\big(-\\tfrac{1}{2}(\\boldsymbol z - \\boldsymbol z^*)^\\top\\boldsymbol A(\\boldsymbol z - \\boldsymbol z^*)\\big) ~ ,$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\\boldsymbol A = - \\frac{d^2}{d\\boldsymbol z^2} \\log \\tilde{p}(\\boldsymbol z)\\Big|_{\\boldsymbol z=\\boldsymbol z^*} ~ .$$\n",
    "\n",
    "And now this distribution we know how to normalize, because it's one of those Gaussians we know and love. By inspection, we can identify the mean and the covariance, and write down the Laplace approximation of $p(\\boldsymbol z)$ as\n",
    "\n",
    "$$q(\\boldsymbol z) = \\mathcal{N}(\\boldsymbol z | \\boldsymbol z^*, \\boldsymbol A^{-1})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, let's use the unnormalized distribution $\\tilde{p}(z) = x e^{-x/2}$. When normalized properly, this is in fact the $\\chi^2$ distribution with $k=4$ degrees of freedom. Have a go yourself:\n",
    "\n",
    "1. Plot $p(z)$.\n",
    "2. Take the first derivative of $\\tilde{p}(z)$ (or the first derivative of its log), and find its maximum $z^*$ analytically.\n",
    "3. In the same plot, draw a vertical line at $z = z^*$ to verify you got the right answer.\n",
    "4. Take the second derivative of $\\log \\tilde{p}(z)$ and evaluate it at $z^*$.\n",
    "5. Plot the corresponding Gaussian $q(z)$ and verify the approximation looks reasonable.\n",
    "\n",
    "**Task 5**\n",
    "\n",
    "* Write a function that evaluates the Laplace approximation $q(z)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## EDIT THIS FUNCTION\n",
    "def laplace_q(z):\n",
    "    # z: double array of size (T,)\n",
    "    # returns: array with Laplace approximation q evaluated\n",
    "    #          at all points in z\n",
    "    return stats.multivariate_normal.pdf(z,2,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f1660898128>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VGX2wPHvmUkjEGpCTeihg5SAIIJKEyyAiNIs2FB3\nXcWuu66rrmv5WdeGooKI0kRUQBRRREVpoRNqgACBAKGFkJD+/v64EzfGhEzIzNwkcz7PM8/M3Pve\ne89d2Tm5bxVjDEoppZTD7gCUUkqVD5oQlFJKAZoQlFJKuWhCUEopBWhCUEop5aIJQSmlFKAJQSml\nlIsmBKWUUoAmBKWUUi4BdgdQGuHh4aZp06Z2h6GUUhXK2rVrjxljIkoqV6ESQtOmTYmNjbU7DKWU\nqlBEZJ875bTKSCmlFKAJQSmllIsmBKWUUoAmBKWUUi6aEJRSSgGaEJRSSrloQlBKKQVUsHEI6tzO\npKVx/PA+Uo/u5+zxRHJPHQSBWtE9adbhIoJCw+wOUSlVjmlCqMCyc/NY9dM3VFn1Os0ytlFbUqlW\nVMFdkPO1g4TAppyq3YngJt2J7NSPsKh2vg5ZKVWOaUKogJJTM1m2ZD6NN7/BxWYTJ6nBnoh+7A5r\niLNGQ4JrR1ItPIpa9ZuQmZnBvs3LObt3FdWObaL5kSXUOPolrIE9Yd2pe8XjVGvTD0Tsvi2llM3E\nGGN3DG6LiYkx/jx1xZaDKSxb/AXdEibTyxFHiqMmRzvdSYvB9+IIKfLZ4E/OZmazfesGklZ+RvfD\ns4iQFI5W70CtQY8Q2O5qcGizklKVjYisNcbElFhOE0L5l5Obx7uL19Jp5f30dWwmNaAOmRfeQ/gl\nd0FQ6Hmfd/uBo6yY9yb9js+iieMoqWEtqNbvQeSC0eBwevAOlFJ20oRQSSSlnOWJT5fxwOFHae08\nRM5lTxLS6w4IrOKxa/y0PYkV8z9g2Jk5tHXsJ6N+N0JGTobwlh67hlLKPpoQKoGl24/w3OxlTMp7\nhubOZJxjZ0DLAV65Vm6e4fPYA2z85n0eMR9SzZmLc+DT0GOCViMpVcG5mxDc+n+6iAwWkR0iEi8i\njxWxv6+IrBORHBEZWWD7ZSKyocArQ0SGu/Z9JCJ7C+zrXJobrMyycvJ4duFW/v7RYj6Sp2kReALn\njXO9lgwAnA7h+h6NufOex7mj6pv8kt0Gvn0UPh4KJxO8dl2lVPlR4hOCiDiBncBAIBFYA4wxxmwt\nUKYpUB14CJhvjJlbxHlqA/FApDEmXUQ+AhYWVbY4/vCEkJyaye0fx3I8cRfzw16kFqeRGz6Hxj19\nFkNKejZ3TY8lcv88ng35lCCnIIOehW7jtTeSUhWQJ58QegDxxpg9xpgsYBYwrGABY0yCMWYTkHeO\n84wEvjHGpLtxTb90Ii2LGz5YRcbhXSyp+QK1HWnITfN9mgwAaoQGMu22C8m74Ab6pT3HroDWsHAi\nzLkJstJ8GotSynfcSQiNgAMFvie6tpXWaGBmoW3/EZFNIvKaiASfxzkrjZT0bG78cBUZx/exoNp/\nqCJZcPNCiOxmSzxBAQ5evq4Towf25vIT9zM97A7M9oUwdQicPmRLTEop73InIRRVR1CqlmgRaQB0\nBBYX2Pw40AboDtQGHi3m2AkiEisiscnJyaW5bIWRmpHNTVNXs+tIKvOiPiMo54yVDBp0sjUuEeFv\n/aN5fXRX/n2iP8+GPYk5vhsmXwYH19kam1LK89xJCIlAVIHvkUBp/0S8HvjCGJOdv8EYk2QsmcBU\nrKqpPzHGTDbGxBhjYiIiSlwjusJJy8zhlqlriDuYwucXJ1In6Sfo/yTUKz/TSgzr3Ig3x3ZhSnIr\nnqr7OsYZCFOvgLgv7Q5NKeVB7iSENUC0iDQTkSCsqp/5pbzOGApVF7meGhARAYYDW0p5zgovIzuX\n26fFsm7/Sd69JoqOm56DyB5WV89y5vL29fnXVe2YFh/KS1GTMA06wWc3w08vQQXquqyUKl6JCcEY\nkwPcg1Xdsw2YY4yJE5FnRGQogIh0F5FE4DrgPRGJyz/e1QMpCvip0Kk/FZHNwGYgHHi27LdTcWTm\n5DJh+lpW7j3Oq9d3ZkDCK1aD7bC3yu0o4fG9m3FHn2a8E3uaqc3/C51GwY/Pwhd3QW6O3eEppcrI\nrcntjDGLgEWFtj1Z4PMarKqkoo5NoIhGaGNMv9IEWpkYY/jHF1v4eWcy/3dtJ4aHrIe4L6DfExDR\n2u7wzunxIW05dCqDZxbvod6Yp7iydgtY9hxkp8O1H0JAkN0hKqXOkw5BtcGsNQeYuzaRe/tHc32H\navD1A1CvI/SeaHdoJXI4hFeuv4CYJrW4/7ONrGl6B1z+PGybD7NvgOwMu0NUSp0nTQg+tinxFP/6\nKo6+rSK4r380fPcEpB2zqoqcgXaH55aQQCfv3xRDZK0q3D4tlt0tb4IrX4Vdi2HmKB2roFQFpQnB\nh06mZXH3J+uICAvm9VGdce79EdZ/Ar3vhYYVa+aOWlWD+Gh8DwKdwh3TYknrdDMMnwR7f4ZPRkJm\nqt0hKqVKSROCj+TlGSbO3kByaibvjOtK7YAsWHAf1GkJlxQ5BKPca1wnlDfHdGXv8TSeXhAHncfC\ntR/AgVXw8XA4e9LuEJVSpaAJwUfeXBrPTzuTefLqdlwQVRNWvAWn9sPQtzw6lbWv9WpRh79c2oI5\nsYks3HQIOlwL138MSRth2lBNCkpVIJoQfOCnncm8/sNORnRpxLgLG0P2WVg9GVoNgSa97A6vzCYO\naEXnqJo8Pm8ziSfToe1VMGYmHN0Gn14HmWfsDlEp5QZNCF6WeDKd+2atp3W9MP5zTUdEBDbOhPTj\ncNE9dofnEYFOB2+M7oIxMHHWBnJy8yB6IIycYk1xMWuM9j5SqgLQhOBFeXmG+2dvICfXMOmGblQJ\nckJeHqx4Gxp2gSa97Q7RYxrXCeXZ4R2I3XeSt36Mtza2GwrD37Eamj8bD7nZ5zyHUspemhC8aNqK\nBNYknORfV7ejWXhVa+POb+F4PPS6p9KtLTC8SyNGdGnEGz/sYk3CCWvjBaPhipdh5zcwbwLk5dob\npFKqWJoQvGTf8TT+79sdXNo6gpHdCgziXvEW1IiCdsPtC86Lnh7WnshaoUyctYGUdNcTQY87YMDT\nEDfP6lmlcx8pVS5pQvCCvDzDY59vxukQnstvNwBIXAv7foWefwGnW7OGVDhhIYG8MaYLR05n8Pcv\nN/9vx8UToe/DsH46LP67JgWlyiFNCF4wY/V+Vuw5zj+ubEvDmgW6lK54E4JrQNcb7QvOBzpH1WTi\ngGi+3pTE4rjD/9tx2T/gwrth5Tvwyyv2BaiUKpImBA9LPJnO84u2cXHLcEZ3L7CMxMl9sPUriBkP\nwWG2xecrd17Sgjb1w/jXV3GkZriqjkTg8ueg4/Ww9N+w/lN7g1RK/YEmBA8yxvD4vM0Y4PkRBaqK\nAFZOAnFAjztti8+XAp0Onh/RkSOpGbzy3c7/7XA4YNjb0PxSmP832PW9XSEqpQrRhOBBn8Um8suu\nYzw+pA1RtUP/t+PsSVj3MXQYCTXOZznqiqlL41rc2LMJ01YksPHAqf/tCAiC66dbq8LNuUmX41Sq\nnNCE4CGHUzL499dbubBZbcZd2OSPO9d+BNlplWYgWmk8dHlrIqoF8/i8zdaAtXwh1WHcXKhaB2Zc\nDyf22BekUgrQhOAxT3y5mezcPF68thMOR4GqopwsWPWeVUVSv6Nd4dmmekggTw9tz9ak00z9NeGP\nO8Pqww3zrLEJ00fAmWRbYlRKWTQheMAP247w/bajPDCwFU3zB6Dl2/I5pCZBr7/ZE1w5MLhDfQa0\nrcurS3Zy4ET6H3eGR8PYOZB6GGbovEdK2cmthCAig0Vkh4jEi8hjRezvKyLrRCRHREYW2pcrIhtc\nr/kFtjcTkVUisktEZotIhVx7MSM7l2cWbqVl3Wrc0rvZnwus/QjCW0HL/j6PrbwQEZ4e1gERePKr\nLZjCYxCiusN1U60ZUj+/XUczK2WTEhOCiDiBt4EhQDtgjIi0K1RsPzAemFHEKc4aYzq7XkMLbH8R\neM0YEw2cBG47j/ht9+Hyvew7ns5TV7cn0Fnof87Th+DASuh4XaWbpqK0GtWswgMDW/HjjmQWbT78\n5wKth8CQ/7OmuFj8d98HqJRy6wmhBxBvjNljjMkCZgHDChYwxiQYYzYBeUWdoDCx+mP2A+a6Nk0D\nKtxcDodOneWtpfEMbl+fi6PD/1xg2wLrvd2wP+/zQ+MvakqHRtV5akEcpzOKmOiuxx3Q86+w6l1Y\n+a7vA1TKz7mTEBoBBwp8T3Rtc1eIiMSKyEoRyf/RrwOcMsbknOc5y4XnFm0jzxj+cWXbogvEfQkR\nbSGitW8DK6cCnA6ev6YTx85k8vbS+KILDfo3tL4SFj8OO77xbYBK+Tl3EkJRdR2lmYimsTEmBhgL\nvC4iLUpzThGZ4EooscnJ5acXyordx1m4KYm7L23xxzEH+VIPw/4V0L7CPfh4VcfIGozsGsmUX/ey\n73janws4nHDt+1C/E8y9FQ5t8H2QSvkpdxJCIlBgDgYigUPuXsAYc8j1vgdYBnQBjgE1RSR/hrdi\nz2mMmWyMiTHGxERERLh7Wa/Kyc3j6QVxNKpZhbsuaVF0oW0LAKPVRUV4+PLW1kjmRduLLhBUFcbO\nhtA6MGMUpCT6NkCl/JQ7CWENEO3qFRQEjAbml3AMACJSS0SCXZ/Dgd7AVmN1M/kRyO+RdDPwVWmD\nt8snK/ex/XAq/7yqLSGBzqILbf0KwltD3WKqk/xY3eoh/OXSFnwbd5gVu48XXSisvtUdNTvdSgoZ\np30bpFJ+qMSE4KrnvwdYDGwD5hhj4kTkGREZCiAi3UUkEbgOeE9E4lyHtwViRWQjVgJ4wRiz1bXv\nUeABEYnHalP40JM35i3Hz2Ty6pKdXNwynMvb1y+60Jmj1jTXWl1UrNv7NKdRzSr8e+FWcvOKqYGs\n1w6u+8ham3nurZCbU3Q5pZRHuDUpvzFmEbCo0LYnC3xeg1XtU/i434Aih+e6qpB6lCbY8uClxTtI\nz8rlqaHt/jh5XUHb5oPJ0+qicwgJdPLYkDb8beZ65q49wKjujYsu2LI/XPkyLLwfvnsChrzg20CV\n8iM6UrkUthxMYXbsAcZf1JSWdc8xhfXWr6BONNQtPFxDFXRVpwZ0a1KLlxbv/N8U2UWJudVaVGjV\nJFjzge8CVMrPaEIohRe/3U6NKoH8rX908YXOJEPCcqu6yM8Ho5VERPjnVe04diaTd5btPnfhQc9C\n9OWw6BHYvdQ3ASrlZzQhuOnnncn8susY91zWkhpVAosvuH2BVheVQueomozo0ogPf9n753mOCnI4\nYeSHENEG5oyH5B0+i1Epf6EJwQ15eYbnv9lOZK0q3NirybkLb/0KareAeh18E1wl8PDg1jgc8MI3\nxXRDzRccBmNnWespzLge0orpoaSUOi+aENzw5YaDbEs6zcOXtyY4oJhupmD9QO39RauLSqlBDWs8\nx9ebk1i998S5C9dsDKNnwukkmH0D5GT6Jkil/IAmhBJkZOfy8uIddGxUg6s7NTx34e0LwORCO+1u\nWlp39m1BverBvPDNtj/PhlpYVHcY/g7s/w0WTISSyiul3KIJoQTTfkvgUEoGj1/R5o8L3xRl61dQ\nq5lfLoRTVlWCnNzXvxXr9p9iydYjJR/QcSRc8hhsnAG//tf7ASrlBzQhnMOp9Cze/jGeS1tHcFGL\nImYzLSj9BOz5SauLyuD6mEiah1flpcU7ih+sVtClj0H7EfD9U7D9a6/Hp1RlpwnhHN7+MZ7UzBwe\nG9Km5MLbF2p1URkFOB08dHlrdh09w7x1bsxfJGJVHTXsAp/fAYc3ez9IpSoxTQjFOHAinWm/7WNk\n10ja1K9e8gFbv4KaTaDBBd4PrhIb0qE+nSJr8NqSnWRku7FyWmAVGDMTQmrAjNGQ6kZ1k1KqSJoQ\nivHKdzsQgQcGtSq5cHaG1buozZVaXVRGIsKjg9twKCWDT1buc++gsPpWd9SzJ2DWWOu/h1Kq1DQh\nFGHLwRS+3HCIWy9uRoMaVUo+4GAs5GZCs77eD84P9G4ZTp/ocN7+Mb7oldWK0uACuOY967/F/Hu0\n55FS50ETQhFe+W4HNaoEFr/WQWF7fwFxQONe3g3Mjzw6uA0n07N5/+c97h/Ubij0fxI2fwY/v+y9\n4JSqpDQhFBKbcIIfdyRz5yXNzz1FRUEJy60VvqrU9G5wfqRDoxpc1akBH/yyl6OppagCuvgB6DQa\nfnwW4r7wXoBKVUKaEAowxvDS4h2EVwtm/EVN3TsoOwMS10DTi70amz96aFBrsnPzePOHYtZfLooI\nXP1fiLoQvrgbDq7zXoBKVTKaEApYHn+MVXtPcM9lLQgNcmupCEhcbbUfNO3j3eD8UNPwqozqHsXM\n1fuLXn+5OIEhMOpTqBYBM8dAykHvBalUJaIJwcUYw8uLd9CoZhXGXFjMYi1FSVhutR800fYDb7iv\nfzSBTgevLtlZugOrRcCY2ZCVBjNHW+9KqXPShOCyZOsRNiamcF//6HNPYFdYwnKrh0tIDe8F58fq\nVg/h5ouaMn/jIXYeSS3dwfXawcgpcGQLzJsAeXneCVKpSsKthCAig0Vkh4jEi8hjRezvKyLrRCRH\nREYW2N5ZRFaISJyIbBKRUQX2fSQie0Vkg+vV2TO3VHq5eYZXvttJ8/CqjOjayP0Ds89q+4EP3Nm3\nOVWDAnittE8JAK0GweXPWSPJlz7j+eCUqkRKTAgi4gTeBoYA7YAxIlJ4bcj9wHhgRqHt6cBNxpj2\nwGDgdREp2BXnYWNMZ9drw3neQ5kt3HSIHUdSmTiwFQHOUjw0HVgNuVnQVMcfeFOtqkHcdnEzvtly\nmC0HU0p/ggvvgm63wPLXYEPhf6JKqXzu/Pr1AOKNMXuMMVnALOAPy4EZYxKMMZuAvELbdxpjdrk+\nHwKOAhEeidxDsnPzeG3JTtrUD+Oqjg1Kd3B++0Hjnt4JTv3utj7NqFElsPRtCWD1PLriJWh2Ccy/\nFxJ+9XyASlUC7iSERsCBAt8TXdtKRUR6AEFAwcVz/+OqSnpNRIJLe05PmLs2kYTj6Tw0qHXJ01sX\nlvALNOgMIW7MdaTKpHpIIHde0pyl24+ybv/J0p/AGQjXT4NaTWH2ODhewhrOSvkhdxJCUb+SpZoX\nQEQaANOBW4wx+U8RjwNtgO5AbeDRYo6dICKxIhKbnJxcmsuWKCM7lzd+2EXnqJr0b1u3dAdnpUNi\nLDTT7qa+cnOvptSpGsSr353HUwJAlVowdjYg1hKcZ88jsShVibmTEBKBqALfI4FD7l5ARKoDXwNP\nGGNW5m83xiQZSyYwFatq6k+MMZONMTHGmJiICM/WNs1YtZ+klAweubw1UtpJ6RJXQ162jj/woarB\nAdx9aQuWxx9jxe7zXE+5TgsY/Smc2g+zb4ScLM8GqVQF5k5CWANEi0gzEQkCRgPz3Tm5q/wXwMfG\nmM8K7WvgehdgOLClNIGX1dmsXN5ZtptezetwUcsSFr8pyt5fQJzafuBjN/RsQr3qwby6ZEfJS20W\np8lFMPQtq8rv6/t1IjylXEpMCMaYHOAeYDGwDZhjjIkTkWdEZCiAiHQXkUTgOuA9EYlzHX490BcY\nX0T30k9FZDOwGQgHnvXonZVg+soEjp3JdG9666IkLLcWZgkO82xg6pxCAp3c0y+aNQkn+WXXsfM/\n0QWj4JJHYf0n8OvrngtQqQrMrfkZjDGLgEWFtj1Z4PMarKqkwsd9AnxSzDn7lSpSD0rLzOHdn/bQ\nJzqc7k1rl/4EWWlwcC30+qvng1MlGhUTxbvLdvPKdzvoEx1e+uq+fJc+bjUuf/+UtRZ2e13tTvk3\nvxyp/NFvCZxIy+KBgef5dHBgldV+oA3KtggKcHBf/2g2Jqbw/baj538iERj2NkT2gC/uhMS1ngtS\nqQrI7xJCakY2k3/eQ782denSuNb5nSRhudV+EKXtB3YZ0bURTeuE8sp3O8jLK0MbQGCItQRntXow\ncxScTPBYjEpVNH6XEKYsTyDlbDb3DzjPpwOwEkKjrhBczXOBqVIJcDqYOKAV2w+nsjjucNlOVjUc\nxs2F3Gz49Drtjqr8ll8lhJT0bD5YvodB7erRMfI8J6PLPGO1H2h3U9tdfUFDWkRU5bXvd5JblqcE\ngIhWMHqG9YQw6wbIyfRIjEpVJH6VED5cvofUjBzuP9+2A3C1H+TohHblgNMhTBzQip1HzvD15qSy\nn7Bpbxj2DuxbDl/puszK//hNQjiZlsWUXxO4smMD2jYow1QTCcvBEWCtyKVsd2XHBrSqV43XPfGU\nANDpOuj3BGyeAz/+p+znU6oC8ZuEMPmXPaRl5XDfgOiynSjhF2jUTdsPygmHQ7h/QCv2JKcxf6OH\nVkbr8xB0uRF+fgnWTffMOZWqAPwiIRw7k8lHvyYw9IKGtKpXhoFkWWnWGr1NensuOFVml7evT9sG\n1fnv97vIyfXAIjgicNVr0PwyWDgRdi8t+zmVqgD8IiG899NuMnNyubd/GZ8ODm8GkwtRRU67pGxi\nPSVEk3A8nS/We+gpwRkI138M4a2tOY+SNnrmvEqVY36REAKdDkb3aEyLiDJW8+T/KDSwbXE3VYyB\n7erRoVF13li6i2xPPCWANa35DXMhpKbVHfXkPs+cV6lyyi8SwiOD2/Cf4R3KfqKkjVC1LoTVL/u5\nlEeJCA8MbMWBE2f5fG2i505cvaGVFHIy4JNrIf2E586tVDnjFwkBOP/5bgpK2ggNLrDqmFW5c1nr\nulwQVZM3l8aTleOhpwSAum1hzCxryuwZo6y1tJWqhPwmIZRZ9lk4us1KCKpcyn9KOHjqLHNiD5R8\nQGk0uQhGTIbENTD3NsjL9ez5lSoHNCG468hWq0FZE0K51jc6nG5NavHW0ngysj38o91+OAx5EXZ8\nDYse0oFrqtLRhOCupA3WuyaEck1EeHBgKw6fzmDm6v2ev8CFd0LviRA7BX5+2fPnV8pGmhDclbTR\nWpO3ZmO7I1El6NWiDj2b1+btH3dzNssLVTv9/wWdRsGPz0LsVM+fXymbaEJwlzYoVxgiwoODWnPs\nTCbTVyZ4/gIOh7WOQsuB8PUDEPel56+hlA00IbgjJwuObtXqogqke9Pa9IkO592f9nAmM8fzF8gf\nuBbZHebdAXuWef4aSvmYWwlBRAaLyA4RiReRx4rY31dE1olIjoiMLLTvZhHZ5XrdXGB7NxHZ7Drn\nG+KRfqFekrwNcrM0IVQwDw5qzYm0LKb9luCdCwSFwtjZUKclzBpnTYuuVAVWYkIQESfwNjAEaAeM\nEZF2hYrtB8YDMwodWxv4F3Ah0AP4l4jkL1M2CZgARLteg8/7LrxNRyhXSJ2jatK/TV0m/7yH0xnZ\n3rlIlVpwwzwIrQ2fjITknd65jlI+4M4TQg8g3hizxxiTBcwChhUsYIxJMMZsAgqPBrocWGKMOWGM\nOQksAQaLSAOgujFmhTHGAB8D5XeF86SNEBRmLcSuKpT7B7Yi5Ww2H/6y13sXqd4AbvwSHE6Yfg2k\neHCktFI+5E5CaAQUHOWT6NrmjuKObeT6XOI5RWSCiMSKSGxycrKbl/Ww/AZlhza5VDQdGtVgcPv6\nTFm+l1PpWd67UJ0W1pNC5mkrKaQd8961lPISd37hiqrbd3dETnHHun1OY8xkY0yMMSYmIiLCzct6\nUG4OHN6i7QcV2P0DW3EmK4fJP+/x7oUadPrfFBfTh8PZU969nlIe5k5CSASiCnyPBA65ef7ijk10\nfT6fc/rW8V2Qc1YTQgXWun4YV3dqyEe/JXD8jJfXSm7aG0Z9Cke3w6cjrTW4laog3EkIa4BoEWkm\nIkHAaGC+m+dfDAwSkVquxuRBwGJjTBKQKiI9Xb2LbgK+Oo/4ve+QjlCuDO4bEE1Gdi7v/rTb+xeL\nHgDXTbUWU5o5WifDUxVGiQnBGJMD3IP1474NmGOMiRORZ0RkKICIdBeRROA64D0RiXMdewL4N1ZS\nWQM849oGcDfwARAP7Aa+8eideUrSRggMhfAyLq6jbNUiohrXdInk4xX7OHI6w/sXbHs1XPOutQb3\n7ButsSxKlXNiKtAEXTExMSY2Nta3F50yxJrU7rbvfHtd5XH7j6fT75VljO4RxbPDO/rmoms/ggX3\nQduhMHIqOAN8c12lChCRtcaYmJLKabeZc8nLg8ObtLqokmhcJ5TRPaKYtfoA+4+n++ai3cbD5c/D\ntvnw1V+tf1NKlVOaEM7lxB7IOqMJoRK5t180AU7hte99OICs11/gsidg0yxYcK8mBVVuaUI4F53y\nutKpWz2E8Rc148sNB9l++LTvLtz3Iej7CKyfDgvv06SgyiVNCOeStAGcQRDRxu5IlAfddUlzqgUH\n8PJiHz4liMBlf4c+D8G6j+Hr+zUpqHJHE8K5JG2Eeu2tmS1VpVEzNIg7+zbn+21HWLf/pO8uLAL9\nnoCLH7Aamxc9qKuuqXJFE0JxjHFNWaET2lVGt/RuRni1IF76dgc+7WknAv2f/N+qa19rUlDlhyaE\n4pzaBxkp2n5QSVUNDuCvl7VkxZ7jLI/38bxDIjDgKbjoXoj9EBY9rElBlQuaEIrz+5TXmhAqq7EX\nNqZRzSq8tNjHTwlgJYWBz8BFf4M178Oih7RNQdlOE0JxkjaCIwDqFl76QVUWwQFOJg6IZlNiCt9u\nOez7AERg4L+tJ4U1H8CCv0GeF9aAVspNmhCKc2gDRLSFwBC7I1FedE2XRrSIqMrL3+0gN8+Gapv8\nJ4VLHoX1n8C8CZDrpcV8lCqBJoSi/N6grNVFlV2A08FDg1qzOzmNz9fZtLBNfpfU/v+CLXPhs/E6\n95GyhSaEopw+BOnHoKH2MPIHgzvU54Komrz63U7OZtlYZdPnARj8AmxfCLPHQbYPJuFTqgBNCEU5\nvNl6r++jCdCUrUSEf1zRlsOnM5jyqxeX2nRHz7vhqtdg13cw43rISrM3HuVXNCEU5ZhrBGtEa3vj\nUD7To1ltBrarx6Rluznm7UV0ShJzKwyfBAm/WMtxnvXh4Dnl1zQhFOXYTqhaF6rUsjsS5UOPDm7D\n2exc3vhB3d0PAAAdaElEQVRhl92hQOexMHKKtcjO1Csh1YZeUMrvaEIoyrFduiCOH2pZtxqju0cx\nY9V+9iSXg6Uv218D4+bAyQT4cJA1+65SXqQJoSjHdmpC8FMTB7QiOMDBi99utzsUS4t+cPMCyEyF\nDy+HpE12R6QqMU0IhaUdh7MnILyV3ZEoG0SEBXPnJS1YHHeENQknSj7AFyK7wa3fWpMsfnQlJPxq\nd0SqknIrIYjIYBHZISLxIvJYEfuDRWS2a/8qEWnq2j5ORDYUeOWJSGfXvmWuc+bvq+vJGztvx131\nx5oQ/NbtfZpRNyyY5xZt8/2UFsWJaG0t4xpWHz4ZATvK5xLkqmIrMSGIiBN4GxgCtAPGiEjh+Rxu\nA04aY1oCrwEvAhhjPjXGdDbGdAZuBBKMMRsKHDcuf78x5qgH7qfs8nsYaZWR3woNCuDBQa1Yv/8U\nizaXo8bcGpFwy7fWdCqzxlqzpSrlQe48IfQA4o0xe4wxWcAsYFihMsOAaa7Pc4H+IiKFyowBZpYl\nWJ84thMCQqBGlN2RKBuN7BZF63ph/N/i7WTllKNJ56rWsdoUWg6AhffD90/ppHjKY9xJCI2AAwW+\nJ7q2FVnGGJMDpAB1CpUZxZ8TwlRXddE/i0ggAIjIBBGJFZHY5ORkN8Ito2O7oHYLcDi9fy1Vbjkd\nwmNXtGHf8XQ+WbnP7nD+KLgajJ4J3W6B5a/BvDsgx+axE6pScCchFPVDXbhi9ZxlRORCIN0Ys6XA\n/nHGmI5AH9frxqIuboyZbIyJMcbEREREuBFuGWkPI+VyaasI+kSH8/r3OzmRVs7mFnIGWCOaBzxl\nzX80fYQOYFNl5k5CSAQK1p9EAoeKKyMiAUANoGAXjdEUejowxhx0vacCM7CqpuyVk2n1+dYGZYU1\npcWTV7UjLSuXl7/bYXc4fyYCF98P134Iiautbqkny9nTjKpQ3EkIa4BoEWkmIkFYP+7zC5WZD9zs\n+jwSWGpc3TNExAFch9X2gGtbgIiEuz4HAlcBW7Dbib1g8jQhqN9F1wvjpl5NmLl6P1sOptgdTtE6\njoQbv4Azh+GDAZAYa3dEqoIqMSG42gTuARYD24A5xpg4EXlGRIa6in0I1BGReOABoGDX1L5AojGm\n4DDLYGCxiGwCNgAHgffLfDdlpT2MVBEmDmhFrdAgnpofV366oRbW9GK4bQkEhcLUK2DTHLsjUhVQ\ngDuFjDGLgEWFtj1Z4HMG1lNAUccuA3oW2pYGdCtlrN6XnxDqtLQ3DlWu1KgSyCOXt+axeZuZv/EQ\nwzoX7lNRTkS0htuXwpybrIbmo9ug3z/BoeNPlXv0X0pBx3ZB9UZWLw6lCrguJoqOjWrw/KLtpGXm\n2B1O8arWsaqPut4My1+F2TdAZjmYl0lVCJoQCtIeRqoYTofw1NB2HD6dwTvL4u0O59wCguDq/8Lg\nF2HnNzDlcji13+6oVAWgCSGfMa5ZTrVBWRWtW5PajOjSiPd/3su+4+V84RoR6HkXjJsLpw7A5Mt0\nDiRVIk0I+c4cgaxUTQjqnB4d0oZAp/DvhdvsDsU9LfvDHT9AlZow7WpY8Y71x49SRdCEkE97GCk3\n1Ksewj39ovl+2xGW7Sgf02+VKDwa7lgKrYfA4sfh89t1aU5VJE0I+X5PCPqEoM7t1oub0rROKM8s\n2EpGdq7d4bgnpAZcPx36Pwlx86zxCsd32x2VKmc0IeQ7tguCqkFYA7sjUeVccICTZ4Z1YM+xNN5Z\nVoF+VB0O6PMg3PC5tSTn5Mt0Gm31B5oQ8h3baY0/KHqOPaX+oG+rCIZ3bsikZfHsOpJqdzil06If\n3PkT1G4GM0fD909DbjnuSqt8RhNCPu1hpErpiavaUTU4gMfnbSYvr4I11NZsDLcuhq43WeMVpl0F\nKYl2R6VspgkBICsdUg5oQlClEl4tmH9c0ZbYfSeZuaYC9vMPDIGhb8KID+DwZnj3Ytjxrd1RKRtp\nQgA47hpopD2MVCmN7BZJr+Z1eOGb7Rw9nWF3OOen03Vw58/WimwzR8Hif0BOOZvuW/mEJgTQHkbq\nvIkIz43oSGZOHk8v2Gp3OOevTgu47XvoMQFWvGWNbj6ZYHdUysc0IYDVfiAOqN3c7khUBdQsvCr3\n9Y/m681JfL/1iN3hnL/AELjiJat76ondMOli2DBDB7L5EU0IYD0h1Gxs/R9CqfNwR5/mtK4XxpNf\nbeFMeZ78zh3thsJdy6FBJ/jybmv21PQTJR+nKjxNCKA9jFSZBQU4eG5ER5JOZ/BKeVxdrbRqNoab\nF8CAp62xCu/0gvjv7Y5KeZkmhLw8q1FZE4Iqo25NanFjzyZ89FsCaxIqwV/UDidcPNGa9qJKTfjk\nWlj0CGSftTsy5SWaEE4nQs5Z7WGkPOKRwW2IrFWFB+ZsqPhVR/kadIIJy+DCu2H1e/BeXziw2u6o\nlBdoQtAeRsqDqgUH8Nr1nTl48izPLqzAvY4KC6wCQ16wFt/JSocPB1ndU7PS7Y5MeZBbCUFEBovI\nDhGJF5HHitgfLCKzXftXiUhT1/amInJWRDa4Xu8WOKabiGx2HfOGiE1zRhzbZb1rQlAeEtO0Nnde\n0oJZaw6wpCL3OipKi37wlxXQbbzVPfXd3rDvN7ujUh5SYkIQESfwNjAEaAeMEZF2hYrdBpw0xrQE\nXgNeLLBvtzGms+t1V4Htk4AJQLTrNfj8b6MMju2EkJoQWseWy6vK6f4BrWjboDqPz9vEsTOZdofj\nWSHV4erX4ab5kJcLU4fAood1qc5KwJ0nhB5AvDFmjzEmC5gFDCtUZhgwzfV5LtD/XH/xi0gDoLox\nZoUxxgAfA8NLHb0n5Pcw0kntlAcFBTh4fVRnTp/N4fF5mzGVsS9/80vg7t/gwrtg9fswqRfs0p5I\nFZk7CaERcKDA90TXtiLLGGNygBQg/0/uZiKyXkR+EpE+BcoXnEmrqHMCICITRCRWRGKTk5PdCLeU\ntMup8pLW9cN4ZHBrlmw9wmdrK+nEccHVYMiLcMs34AyGT6+Fz8bD6SS7I1PnwZ2EUNSfzoX/3Cmu\nTBLQ2BjTBXgAmCEi1d08p7XRmMnGmBhjTExERIQb4ZZCRgqcOaw9jJTX3Nq7GT2b1+bp+XEcOFGJ\nG2Cb9IK7f4XLnoDti+Ct7rDqPatKSVUY7iSERCCqwPdI4FBxZUQkAKgBnDDGZBpjjgMYY9YCu4FW\nrvKRJZzT+47lT2qnTwjKOxwO4eXrLsAhwgNzNpBb0abJLo2AYLjkYavROao7fPMIvN8PDq23OzLl\nJncSwhogWkSaiUgQMBqYX6jMfOBm1+eRwFJjjBGRCFejNCLSHKvxeI8xJglIFZGerraGm4CvPHA/\npaNdTpUPRNYK5amh7VmTcJI3fthldzjeV6cF3DAPRk6B1CQrKSy8H9KO2x2ZKkGJCcHVJnAPsBjY\nBswxxsSJyDMiMtRV7EOgjojEY1UN5XdN7QtsEpGNWI3Ndxlj8odw3g18AMRjPTn4fi2/4/EgTqjV\nxOeXVv5lRNdGjOjaiDeW7mLZjqN2h+N9ItDhWrhnDXS/A9ZOgze7WNVIujpbuSUVqfdDTEyMiY2N\n9dwJP78DDqyCiZs8d06linE2K5dr3vmVw6czWPi3i4msFWp3SL5zdBt8+xjsWQYRbWHw89DiMruj\n8hsistYYE1NSOf8eqXxqvzWJl1I+UCXIyaQbupGba/jrp+vIzPGjBte6beHGL2H0DGuqmOnDYdY4\nOLHH7shUAX6eEPZBTa0uUr7TLLwqL113ARsTU3h24Ta7w/EtEWhzJfxlFfR/Enb/CG/1sCbMSztm\nd3QKf04IOZlWg5c+ISgfG9yhPhP6Nmf6yn18uf6g3eH4XmAI9HkQ7l0HXW6ANR/AfzvDzy9BVprd\n0fk1/00IKa6BQpoQlA0eubw1PZrW5vF5m9lxONXucOwRVt+aAuMvK61Rz0ufhTe6wtqPtOHZJv6b\nEE7ts941ISgbBDgdvDW2C1WDA7j7k7WkZmTbHZJ9IlrB6E/h1sVWj78F98E7PWHzXGu9EuUzfpwQ\n9lvvmhCUTepWD+GtsV3YdyKd+2ZtICfXz3/8Gve0ksKoT8EZCJ/fBpMugrgvNTH4iH8nBEcAhDWw\nOxLlx3o2r8PTQ9uzdPtRnloQVzknwSsNEWh7Fdz1qzWwzeTCZzdbi/Js/xr8/X8fL/PvhFC9ETgD\n7I5E+bkbejbhzr7N+WTlfib/rN0wAXA4rIFtf1kJI96H7HSYNRYmXwJb5+sTg5f4d0LQ6iJVTjw6\nuA1XdmrA899sZ+Em30/rVW45nNDpevjrahj2DmSmwpwb4Z0LYcMMyPXjthcv8O+EoFNWqHLC4RBe\nue4CYprU4oE5G4lNOFHyQf7EGQBdxsE9sVZVkjMIvrwb3ugCqyZD9lm7I6wU/DMh/D4GQROCKj9C\nAp28f1MMjWpW4faPY9mTrCuQ/YnDaVUl3bUcxn4G1RvCNw/D6x1h2Ys6wK2M/DMhnHKt96NVRqqc\nqVU1iI9u6Y5DhPFT13C8si2/6Ski0GqQ1Stp/CJo2AWWPQevtoP5f7PmTlKl5qcJQccgqPKrSZ2q\nfHBzDEdOZzB+6hpSzmo9ebFEoGlvGPcZ/HUNdB4Lm+ZY4ximj4D477VnUin4aULQMQiqfOvauBaT\nbujK9sOnuenDVZoU3BHRyhr5fP9W6PcEHNkCn1xrrd62chKcPWV3hOWe/yYEHYOgyrl+beoxaVw3\ntiad5qYpqzntz6OZS6NqHej7MEzcDNe8B1VqWlNvv9rWqk5K2mh3hOWW/yaEGpFWA5VS5diAdvV4\nZ1w3th5K4cYPNSmUSkAwXDAabv8eJvwEHUfCps+sQW4fDID1n+pkeoX4b0LQ6iJVQQxsV4+3x3Zl\n66EUbtKkcH4adoahb8KD22DwC3D2JHz1F3i5Ncy/FxJjta0BNxOCiAwWkR0iEi8ijxWxP1hEZrv2\nrxKRpq7tA0VkrYhsdr33K3DMMtc5N7hedT11UyXShKAqmEHt6/P22K5sOZjCzVNW+/dkeGVRpRb0\nvNsaz3DLN9D2atj8GXzQH97pBSve9uuuqyUmBBFxAm8DQ4B2wBgRaVeo2G3ASWNMS+A14EXX9mPA\n1caYjsDNwPRCx40zxnR2vXyz0Gx2Bpw5rGMQVIUzqH193h7Xlc2JVvWRdkktAxFochFcMwke3AFX\n/xeCqsLiv8MrreHT663ZVrPS7Y7Up9x5QugBxBtj9hhjsoBZwLBCZYYB01yf5wL9RUSMMeuNMfnj\n8OOAEBEJ9kTg503XQVAV2OXt6/POuK5sSzrNiEm/6eA1TwipDt3Gwx0/wN0roNdf4fBma7bVl6Ph\ni7th91LIq/xLnrqTEBoBBwp8T3RtK7KMMSYHSAHqFCpzLbDeGFPwz5qpruqif4qIlCry8/X7GAR9\nQlAV06D29Zk5oSdnMnIYMek3Vu/VaS48pl47GPgM3B8HNy+E9tfA9oUw/Rqrl9KihyHh10qbHNxJ\nCEX9UBdufTlnGRFpj1WNdGeB/eNcVUl9XK8bi7y4yAQRiRWR2OTkZDfCLYEOSlOVQNfGtZj3l4uo\nHRrEDR+s4qsNfrgUpzc5HNCsDwx7Cx7aBdd/DFE9YN3H8NEV1ojoSpgc3EkIiUBUge+RQOHpGH8v\nIyIBQA3ghOt7JPAFcJMxZnf+AcaYg673VGAGVtXUnxhjJhtjYowxMREREe7c07md2g+OQGv5PqUq\nsCZ1qjLvLxfROaom983awNs/xut6Ct4QGALthsGoT+Dh3XDthxDV/Y/JYeH9sOt7a560CsydxQDW\nANEi0gw4CIwGxhYqMx+r0XgFMBJYaowxIlIT+Bp43Bjza35hV9KoaYw5JiKBwFXA92W+G3foGARV\nidQMDWL67T14ZO4mXlq8g/3H03lmeHuCA/Tft1cEV7PGM3QcCZlnYOe3sPUr2DgbYqdAUBhED4DW\nV0L0QGtQXAVSYkIwxuSIyD3AYsAJTDHGxInIM0CsMWY+8CEwXUTisZ4MRrsOvwdoCfxTRP7p2jYI\nSAMWu5KBEysZvO/B+yqedjlVlUxwgJPXR3WmSe1Q3lgaT1xSCm+O6Uqz8Kp2h1a5FUwO2Rmw9ydr\nVbcd30DcF9ZsCI17QfQgKzlEtLF6N5VjUpEeMWNiYkxsbGzZTvJyK+s/0LC3PBOUUuXIkq1HeOiz\njeTk5vHciI4M61y4/4fyurw8OBhrJYddS+BonLW9RhS0HGD9/jTrayUUHxGRtcaYmBLL+VVCyD4L\n/6kPlz0BlzzsucCUKkcOnTrLfbPWsybhJNfHRPLU0PaEBulSsbZJSbRmXd21BPYsg6wz1tNDZA9o\ncRk0vxQadvXqcr6aEIpybBe8FQPXTIYLRnkuMKXKmZzcPP77wy7e+jGeFhHVeGtsF9rUr253WCon\nC/avsMY17FnmmmjPQHB1aHoxNLvEeq/bzurp5CHuJgT/+rNBu5wqPxHgdPDgoNb0bF6HibM3MOyt\nX7l/YCtuu7gZgU7/nMKsXAgIguaXWC+A9BOw92fY86OVIHYssrZXqQVNeluvpr2hXgefdITxs4Sg\n6yAo/9K7ZTiL7u3D37/YzAvfbOfL9Qf5zzUd6Naktt2hKYDQ2tB+uPUC6zcq4VdIWA77lluD4gBC\nalgrw9Xv4NVw/CshnNznGoOg6yAo/xERFsz7N8XwXdxhnpofx7WTVjCmR2MeG9yGGqGBdoenCqrZ\nGDo3hs5jrO8piVaC2Pcr1Gnh9cv7V0I4tR9qRnm0bk6pimJQ+/r0bhnOa0t2MvW3BJZsPcwTV7Zj\nWOeG+GrmGFVKNSKt9k4ftXn61y+jjkFQfq5qcABPXNWO+ff0plGtUCbO3sDId1fw227/nfJZ/Y8m\nBKX8UPuGNZh390U8P6IjB0+eZez7qxgzeSVr9+lEef7MfxJC9llIO6oJQSkXp0MY06Mxyx6+lCev\naseuo6lcO2kF46euZnNiit3hKRv4T0I45ZrBW6e9VuoPQgKd3HpxM35+5DIeHdyGDQdOcfVby7n1\nozX8vDOZvLyKM1ZJlY3/NCprl1Olzik0KIC7L23BDT0bM2V5AtNXJnDTlNU0D6/Kjb2aMLJbJGEh\n2iupMvOjJwQdlKaUO8JCArlvQDS/PtaP10d1pkZoIE8v2ErP537gn19uYcfhVLtDVF7iX08IjkCo\npusgKOWO4AAnw7s0YniXRmxKPMXHK/YxO/YA01fuo039MK6+oCFXdWpAkzo6q2pl4T9zGX12CyRt\ngHvXezYopfzI8TOZLNh4iIWbkojddxKATpE1uLpTQ67s1ICGNavYHKEqik5uV9j7/awJpG760rNB\nKeWnDp46y9ebDrFgYxKbD1q9ktrUD6NPdDh9oiPo0aw2IYG6UE95oJPbFXZqP7S+wu4olKo0GtWs\nwoS+LZjQtwUJx9L4Nu4wv+xKZtpv+3j/l70EBTi4sFlt+kSHE9O0Nu0bVteV3Mo5/0gIWemQlqwN\nykp5SdPwqtx1SQvuuqQFZ7NyWbX3OL/sOsYvu5J5btF2AIKcDto2rE6XqJp0aVyTLlG1iKpdRafN\nKEf8IyGk6BgEpXylSpCTS1vX5dLWdQE4ejqDdftPsf7ASdbvP8XsNQf46LcEAMJCAmhVL4xW9aq5\n3sOIrleNiGrBmihs4FZCEJHBwH+x1j/+wBjzQqH9wcDHQDfgODDKGJPg2vc4cBuQC9xrjFnszjk9\nSscgKGWbutVDGNyhPoM7WD38cnLz2HEklQ0HTrE9KZUdR1L5ZsthZq4+8Psx1UMCiKodSlStUBrX\nCSWqVhWiaocSWSuUetWDqRYcoAnDC0pMCCLiBN4GBgKJwBoRmW+M2Vqg2G3ASWNMSxEZDbwIjBKR\ndsBooD3QEPheRFq5jinpnJ6jYxCUKjcCnA7aN6xB+4Y1ft9mjCH5TCa7jpxhx+FU9h5L48DJdHYe\nTWXpjqNk5eT94RyhQU7qVQ8hIiyYetVDqBsWTO2qQdQMDaRWaJD1qmp9rh4SSEigQxOIG9x5QugB\nxBtj9gCIyCxgGFDwx3sY8JTr81zgLbH+1x8GzDLGZAJ7RSTedT7cOKfnnNoPziCoVs8rp1dKlY2I\nUDcshLphIfRuGf6HfXl5hqOpmRw4mU7iyXSOns7kaGomR05ncDQ1k82Jpziamkl6Vm6x5w9wCNVC\nAqgWHEBYSCBhIQFUDXISGhRASKCT0CAnVYKcVAm03oMDHAQFOAgOsD7nfw9yOghwOgh0CoFOh+sl\nBDgcOJ1CgENwOgSnCE6n690hOH5/p1wnJncSQiPgQIHvicCFxZUxxuSISApQx7V9ZaFjG7k+l3RO\nzzm1H2roOghKVUQOh1C/Rgj1a4TQvWnxK71lZOdyKj2bk+lZnEzL4mR6NifSs0jNyOZMRg5nMnM4\nk5HD6YwcUjOyOZ6WxYGTZzmblcvZ7FzSs3LIyM4r9vweux8Bh1hJQn7/zO/fRf6XOATXu8Bnd/ai\nabh3BwG6kxCKSmeFBy8UV6a47UX9Mhc5IEJEJgATABo3Ps8qn/qdtEFZqUouJNBJ/RpO6tcIOe9z\n5OUZMnJyycrJIysnj8ycPDJzcsnItj7n5OaRnWvIzs1zvazPuXmG3DxDTp4hNy/P9e56GUNeniE3\nj98/5xmDAevd4NqW//2P+wxWlVpocPlYUzkRiCrwPRI4VEyZRBEJAGoAJ0o4tqRzAmCMmQxMBmtg\nmhvx/lmfB87rMKWUf3E4hNCgAEKD7I7EHu7UoawBokWkmYgEYTUSzy9UZj5ws+vzSGCpsYZAzwdG\ni0iwiDQDooHVbp5TKaWUD5X4hOBqE7gHWIzVRXSKMSZORJ4BYo0x84EPgemuRuMTWD/wuMrNwWos\nzgH+aozJBSjqnJ6/PaWUUu7yn7mMlFLKT7k7l5F2u1FKKQVoQlBKKeWiCUEppRSgCUEppZSLJgSl\nlFJABetlJCLJwL7zPDwcOObBcCoCvWf/oPdc+ZX1fpsYYyJKKlShEkJZiEisO92uKhO9Z/+g91z5\n+ep+tcpIKaUUoAlBKaWUiz8lhMl2B2ADvWf/oPdc+fnkfv2mDUEppdS5+dMTglJKqXPwi4QgIoNF\nZIeIxIvIY3bH400iEiUiP4rINhGJE5H77I7JV0TEKSLrRWSh3bH4gojUFJG5IrLd9d+7l90xeZuI\n3O/6d71FRGaKyPmvhlNOicgUETkqIlsKbKstIktEZJfrvZY3rl3pE4KIOIG3gSFAO2CMiLSzNyqv\nygEeNMa0BXoCf63k91vQfcA2u4Pwof8C3xpj2gAXUMnvXUQaAfcCMcaYDlhT54+2Nyqv+AgYXGjb\nY8APxpho4AfXd4+r9AkB6AHEG2P2GGOygFnAMJtj8hpjTJIxZp3rcyrWj0Sjcx9V8YlIJHAl8IHd\nsfiCiFQH+mKtRYIxJssYc8reqHwiAKjiWpkxlGJWWqzIjDE/Y60rU9AwYJrr8zRguDeu7Q8JoRFw\noMD3RPzgBxJARJoCXYBV9kbiE68DjwDeXyW9fGgOJANTXdVkH4iId1dgt5kx5iDwMrAfSAJSjDHf\n2RuVz9QzxiSB9UcfUNcbF/GHhCBFbKv0XatEpBrwOTDRGHPa7ni8SUSuAo4aY9baHYsPBQBdgUnG\nmC5AGl6qRigvXPXmw4BmQEOgqojcYG9UlYs/JIREIKrA90gq4WNmQSISiJUMPjXGzLM7Hh/oDQwV\nkQSsKsF+IvKJvSF5XSKQaIzJf/qbi5UgKrMBwF5jTLIxJhuYB1xkc0y+ckREGgC43o964yL+kBDW\nANEi0kxEgrAaoebbHJPXiIhg1StvM8a8anc8vmCMedwYE2mMaYr133epMaZS/+VojDkMHBCR1q5N\n/bHWLq/M9gM9RSTU9e+8P5W8Ib2A+cDNrs83A1954yIB3jhpeWKMyRGRe4DFWL0Sphhj4mwOy5t6\nAzcCm0Vkg2vb340xi2yMSXnH34BPXX/o7AFusTkerzLGrBKRucA6rN5066mEI5ZFZCZwKRAuIonA\nv4AXgDkichtWYrzOK9fWkcpKKaXAP6qMlFJKuUETglJKKUATglJKKRdNCEoppQBNCEoppVw0ISil\nlAI0ISillHLRhKCUUgqA/wceEqBx6O32OQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1660898ef0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## ADD CODE HERE\n",
    "# Find the Laplace approximation of x*exp(-x/2) with pen and paper and then plot it.\n",
    "z = np.linspace(0,10)\n",
    "p = stats.chi2.pdf(z, 4)\n",
    "c=laplace_q(z)\n",
    "c = c/c.max() * p.max()\n",
    "plt.plot(z,c)\n",
    "plt.plot(z, p)\n",
    "# plt.xlim(1.8, 2.2)\n",
    "# plt.ylim(.125, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayesian logistic regression (for real this time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have obtained the mode (peak) of the posterior through the MAP estimate above, it's time to go all the way and calculate the posterior over $\\boldsymbol\\theta$. However, as we mentioned above the posterior doesn't have an analytical form, so we'll use &ndash; you guessed it &ndash; the Laplace transform.\n",
    "\n",
    "**Task 6**\n",
    "\n",
    "* Write a function, based on your previous code, that will calculate the Laplace approximation $q(\\boldsymbol\\theta)$ of the true posterior $p(\\boldsymbol\\theta | \\boldsymbol X, \\boldsymbol y)$ and return the mean and variance of $q$.\n",
    "\n",
    "To visualize the behavior and the diversity of $q$, draw a number $j = 1, ..., J$ of samples $\\boldsymbol\\theta_j \\sim q(\\boldsymbol\\theta)$. For each sample, plot its predicted class probabilities $\\sigma(x \\boldsymbol\\theta_j)$.\n",
    "\n",
    "_Hint_: the extension of the Laplace approximation to multivariate distributions is straightforward, and in this case the variance of the Gaussian is the Hessian of the negative log likelihood $\\boldsymbol A = - \\nabla_\\theta \\nabla_\\theta \\log p(\\boldsymbol\\theta | \\boldsymbol X, \\boldsymbol y)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## EDIT THIS FUNCTION\n",
    "\n",
    "def Hessian(X,Y, theta):\n",
    "    return X.T @ (np.identity(X.shape[0])*sigmoid((X @ theta) - Y)) @ X\n",
    "\n",
    "def get_posterior(X, y, m, S):\n",
    "    # X: N x D matrix of training inputs\n",
    "    # y: N x 1 vector of training targets/observations\n",
    "    # m: D x 1 prior mean of parameters\n",
    "    # S: D x D prior covariance of parameters\n",
    "    # returns: maximum a posteriori parameters (D x 1)\n",
    "    #          covariance of Laplace approximation (D x D)\n",
    "\n",
    "    mu_post = map_estimate(X, y, m, S)#np.zeros((D, 1)) ## <-- EDIT THESE LINES\n",
    "    S_post  = np.linalg.inv(Hessian(X,y,mu_post) + np.linalg.inv(S))#np.eye(D)\n",
    "\n",
    "    return mu_post, S_post"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.25621289 -0.1440992 ]\n",
      " [-0.1440992   0.10139411]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/ipykernel/__main__.py:15: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VUX+x/H33PTeG2l0kE4IvSuICIJiQ1esK/ZesKA/\ndUVX3dVVQRQsYC+ALE3pEHrvhBrSe+/JLfP7IyEkJNHgkv59PY8PyTlzzp2JySeTuXNmlNYaIYQQ\nLYuhsSsghBDi8pNwF0KIFkjCXQghWiAJdyGEaIEk3IUQogWScBdCiBZIwl0IIVogCXchhGiBJNyF\nEKIFsm6sF/b29tZt27ZtrJcXQohmad++felaa58/K9do4d62bVv27t3bWC8vhBDNklIqpi7lZFhG\nCCFaIAl3IYRogSTchRCiBWq0MfeaGI1G4uPjKS4ubuyq1Ct7e3uCgoKwsbFp7KoIIVqoJhXu8fHx\nuLi40LZtW5RSjV2deqG1JiMjg/j4eNq1a9fY1RFCtFBNalimuLgYLy+vFhvsAEopvLy8WvxfJ0KI\nxtWkwh1o0cF+XmtooxCicf1puCulvlRKpSqljtZyXimlPlJKnVFKHVZKhV3+agohhLgUdRlzXwDM\nBr6u5fx4oFP5fwOBueX/tgivvfYazs7OPPvsszWeX7p0KZ07d6Zbt24NXDNx3tIDCby3+iSJ2UW0\ncXdgdFcfNp5Iq/L5ikNJZBcZAfBwtOH/ruvO9X0Da7zfzKVH+G5nLOd3F3aytWLWDT1rLH/xaz83\nrkut973U8pfr3uePJ2QXYaUUZq3xcLQhv9iI0VL1HoEXff3cHW0oMZopLC/o7mDDa5O6szcmkx92\nxWHWGiulGNTeg+iMolrreqltaWkao/1/Gu5a6wilVNs/KDIZ+FqX7bS9UynlrpQK0FonXaY6NmlL\nly5l4sSJEu6NZOmBBF5ccoQioxmAhOwivt0ZW3H+4s8BsgqNPLfoEEC1H7CZS49UK19QauaZX6qX\nr+m1X1xypMb7Xmr5y3XvvTGZLN6XUHHcrHXF16AmF3+9Li6XXWTkyZ8OVjlm1pptZzOr3KNyXS+1\nLS1NY7X/csyWCQTiKn0eX36s3sO9vn4bzpo1i4ULF+Lr60twcDD9+vVj/vz5zJs3j9LSUjp27Mg3\n33zDwYMHWbZsGZs3b+bNN99k8eLFbNiwoVo5R0fHy9BaUZP3Vp+s+KG5FEaz5r3VJ6t9v/ywK67G\n8mZL9fI1vXaR0Vzjfc+Xf63zVwDMOHbnH5b/K/euqfz53nVDq1zXS21LU6S1xlRSjKk4H2NxQfm/\nhRQWFZBfVEBhcTFFJUYKS0opMpooMpkpNlsoNVuITMlnhDNopbAYrLAYFEVZmvdW2zb5cK/p3cEa\nv5uUUtOB6QAhISH/04vW12/Dffv28eOPP3Lw4EFMJhNhYWH069ePKVOmcP/99wMwc+ZMvvjiCx57\n7DEmTZrExIkTuemmmwBwd3evsZyoH4nZRZf12j8KwovL1/baf3Tcs3N6ncr/lXvXpDGC/bzzdbrU\nOl9upaWlFBUVUZiXQ1FWMoVZyeSmppCRX0BGsYlMsyZbW5GjbMi1tqPAypYiK1uKbWwptralxNqW\nEmsbjFbWmKysMRmsMFpZYbKyQSsPwKMsBe3L/7tYDVE3bu8mEjPqt/2XI9zjgeBKnwcBiTUV1FrP\nA+YBhIeH/0/fdfXVG9iyZQs33HBDRW970qRJABw9epSZM2eSnZ1Nfn4+48aNq/H6upYTl0cbdwcS\n/mJItHF3qHbs/Jh0XcrX9to13ff8cStL9R+5msr/lXvXVP6P2lPfztf1UttSV0ajkZycHLKzs8lO\nTSA3PprctFRSC4qIV7ak2DmTYe9Cvp0Dhbb2FNraUWRrT4GtH8WeweBV+6w1O2MpdsYSHIyl2BtL\n8SzMx85sxMFswt5iwlabsTObsbGYsbKU/2s2YTCbUCYjymgEYymUllCcX4i1qQRrS9l5K7OJg/Y9\naRPyv7X/z1yOcF8GPKqU+pGyN1JzGmK8vT57AzVNVbz77rtZunQpvXv3ZsGCBWzatKnGa+taTlwe\nz43rUuUvuLqysVI8N65LteO3DQyuNuYOYGWoXr6m13awsarxvufLq7NVv7dqK/9X7l1T+Rv7BVYZ\nc28olet6qW2pzGg0kpGRQWpqKinR0WTGnSMrO4tUiyLRyZ0sR1eyHF3IcXAiz6EteZ27YbSu+uS3\ntdmEW1EBHiUFBOVn4GMx4mdQ+Nrb4uXohJeLG55urni4uuBmZwslhRRmZZKfkUNeRjp5WenkZaRT\nkJlBQU42Rbm5aG2pXlmlsHd2wcHZBXuXsn/tvN1JKNBsjs6nQNlQYmtLqcGOHCd/Xq5D+/8Xfxru\nSqkfgFGAt1IqHvg/wAZAa/0psAq4FjgDFAL31FdlK6uv3sCIESO4++67eeGFFzCZTCxfvpwHHniA\nvLw8AgICMBqNfPfddwQGlv114OLiQl5eXsX1tZUT9eP8X2mXa7bMm9f3BKjTbJmaXvuP3ve5vm8g\nm86WjVkq+MPyf+XetZUPD/Vs1NkydW1LXl4e8fHxxJ06RUpcHOlZmeSYzOQ6OJHq4kGqizvZvp3J\nbOdKvoNTxXW2JiPBhdn0MuURYiwmxNGNUN8AQn29CHK0x8PaqkqHrbggn8yEeLJTkshOTiLn9GHO\npiSTk5pMQXZWta+to5s7Ll7euPj44t+pC05u7ji6e1T86+jqjoOLC3ZOThgMVjX+//G56P3Blxtg\ntozSjfQnW3h4uL54PffIyEiuuOKKOl1/8Zg7lPUG3p5S85S1S1H5DdWQkBDCwsJwcnLi3XffxcfH\nh4EDB5KXl8eCBQvYtm0b999/P3Z2dixatIg1a9bUWO5il9JW0XJsWnQNAKNu+r2Ra9K4TCYTCQkJ\nxJ4+TdypUySmp5NvsWBSBlJcPciycybNxZ0EL9+KILczG+lsyqeLgw1d/Pzp4u1DFyd7guxtsarh\nr+3S4iLSos+RkRBLRlwsGQlxZMTFkJ91YWYPSuHs6YW7nz/ufgG4+wXg6u2Di5cPzl7eOHt6Yd3E\n1oBSSu3TWof/abnmGu7QvOfOSri3Tq013LXWZGZmcurQIU4dOkRcdjYmpdBAqTaQ5uhKtF8Ap9u0\npcTGDoBOuoC+rg6E+QYQ5uHGFU4O2BhqHic3FheTEn2W1KgzJEedISXqDJmJ8VCeb9Z2dngFBuMV\nFIJXUAiegcF4BLTBzccPa1vbhvoyXBZ1DfcmtXDYpbq+b2CzCXMhWhuLxULUyZMc2bKFqMREzg9e\nOubnYzLYcDYkiKMhnUhy8QWggyriDg9HRgWGMNDdGVfrmoc4AApzc0iIPEZ85FHiIo+SHhNdMQ7u\n5OGJX/uOdB0yAt92HfAODsXV2wdlaHKrrdSrZh3uQoimRWtNYnQ0+9esITIhgUKDAWujEb+0dNyd\nbYlsG8Cy/leT6OiLozYywsHCk218GeXrRaiDXa33LSksIObwAWKPHiY+8igZ8WVvelvb2tGmcxcG\nTrkF/w6d8GvXEWdPr4ZqbpPW5MJda93iF9ZqrKEwIepLdkoKe5Yv51hMDNlWVhjMZtpkZNDL1ZbI\ndq78OHAMx1w6YaUtjLIr5pV2bRjn542jVc29aa01WUkJRO3bTdSBvSScOIbFbMbG3oHArt24Ytgo\ngrr1xL9DR6ysm9aYeFPRpMLd3t6ejIyMFr3s7/n13O3ta3raQYjmJW7fPrasWMlpswltMOCTm8sI\nDw98urrxk40LszyHU2TlQB+Vx5vBrkwODsbHtuYw1lqTcvY0J7Zv5uze3WSnlM2o9g4Opd/EG2jf\nN5w2na/AYFX7cI24oEmFe1BQEPHx8aSlpTV2VerV+Z2YhGiOLKWlHF/6X3bs3UuCowPWJhPdLJoh\nVw4j1zmVOQnZLHUfBEpxg0MRD3fvxBWuTrXeLzMxgRPbNnFi22aykhKxsrYmpGcf+k24nvZh/XH1\n8W3A1rUcTSrcbWxsZHciIZooU14+e7/4nN2xcWS6umBvbcUgZxeGTr2VM4k7mZEYzzqrMBzdS7jX\npYTpPcMJqmUcvSg/j+ObNxC5dRMpUadBKUK696T/5JvoNGAI9k7ODdy6lqdJhbsQoumxFBVxasEC\nNkRGkurpiauDA2M7d2bAjTeSnRXD/+3+ne9dwvF0zud510Lu6dkfj1qGXlKjozjw+wpObN2EyViK\nX/uOjJx2H12HjJA3Qi8zCXchRI0sJSUkff89G7dt42ybNti6uTGubxgDr5uIxWxkwbZfedcYRKFT\nHx62TeXp/qNxriHUzSYTp3dv5+DqFSScOI61nR3dRlxJn3ET8AmVv9Tri4S7EKIKbTaT+csidixb\nxqG2oZgCAwnr0IExN9+Mg4MDO07t4aWzKUTad2WEJYZZvcPp5Fv9mRpjSTGH1qxi38ql5Gdl4ubn\nz8hp99Fj1FjsnWXYpb5JuAshKhSfPMnxN2ex2cuTnM6dCPXyYsKtt+Lr60t+STGPrfuVX6zaEWiw\n5wuPLK7tPanazDZjaQmH1/7G7v8uojAnm5AevRk7/THa9enX6h4kakwS7kIILMXFpM2ew+6ICA70\n7oWDvT23TplC165dUUpxLC2R+w8cJ9o6hCeMR3hi+CQcHd2q3MNYWsKRdb+z+7+LKMjOIqRHLwY/\n/SJBXbs3UqtaNwl3IVq5gu3biX5zFtsDAkgI60uH0FCm3HILTk5OaK35OvIwryQW464Vi90SGdxv\nWpXrLRYzRzesZfui7ynIyiS4W08mPjGDoG49GqlFAiTchWi1zHl5pLw5izPbtrFz+DCK7e25euxY\nBg0ahMFgIM9k5rndO1la4sSo/Eg+7tMTn9Crqtwj8dQJNnz1KSlRZ2jTpRsTHnuW4O69GqlFojIJ\ndyFaoeKTJ4l7/HEOubhybMxVuLu787ebb67Yf+BobgHT9x0kWtvzYuZvPHb13RhcAyquL8zJZssP\nCzm6cS3OHp5c+/hzdB0yosU+Wd4cSbgL0cpkL11Kwhv/YPfAAcT4+dGzZ08mTJhQsSTG2pR07j92\nDvfSQhYbNzD4+ufBuuxhJIvZzKG1q9j287cYi4sJv24Kg2+ciq2DbALf1Ei4C9FKWEpKSJn1FmlL\nlrBzwrUkOjgwduxYhgwZUtHj/jUxjcdOxNIt/yzfeefhM3gmlJ/LiI9l1ex/k3ruLCE9+3DlPQ/g\nFRj8Ry8pGpGEuxCtQGl8AglPPEH2mTNsv+Vm0rXm+smT6dOnT0WZBXEpvHg6kUE5h/k6EFwG3A+A\ntlg48PtyIr5fgK29AxOffIHOg4bKEEwTJ+EuRAtXsHMn8U88Sb6dHVtvv408o5HbbrmFzp07A2Wr\nMX4Unczb0SmMzdzBvDYGHAY8CEBeZjqr535IzOEDtA/rz9UPPI6Tu0djNkfUkYS7EC1Y7tq1JD79\nDAVXdGVDv36Yteauu+4iOLhsOEVrzRtnEpgbn86NKWv4T6AtNkOeBODkji2smz8Hk8nImL8/Qq8x\n10hvvRmRcBeihcpe8itJM2eSN6A/6zt3xs7GhrvuuANf37IldM1a8/zJOL5LyuSehCXMCnTEMPxJ\njMXFrPt8Dse3bMS/Y2fGP/IMnm1kO8vmRsJdiBYoc+FCUt7+J6UjR7ChbVucnJy46667cHMre6pU\na83MU/F8l5TJUzELeT7IDTVqBjmpKfz3X2+SFhvNoBtvY9CUW7GylphojuT/mhAtiNaa9I8/Jv2T\nuahrrmF9cBDWWjNt2rSKYAf4NC6NrxIzeDDuR2b4O8BV/0dc5FGWv/82FrOZKS+8Rrs+/RqxJeJ/\nJeEuRAuhLRZSZr1F1nffYXfTjazy96ckL4977rkHD48Lb4IuS83m9bOJTEzbzKt28eirv+bQmlVs\nXDgPd78AJj/3igzDtAAS7kK0AFprkma+Qs6SJbjefTe/+3iTlZjIHXfcgb+/f0W5Xdn5PHY8mgF5\nkcxO/QHLXSvZMH82RzasoX1Yf6597FnsHGvfEk80HxLuQrQAae+/T86SJXg+9CAbfH2JO3mSm2++\nucq2lWcKi7n7yDkCS9JYEPkalpu+45d3/0niqUgG3nALQ275GwaDbD7dUki4C9HMZSxYQMb8z3Gb\neiu72rXj5P79jB8/nu7dLyy1m1Zq5PZDURiMBXx/4HGsR8zkp9mfk5uWysQnZ9Bl8PBGbIGoDxLu\nQjRjOctXkPrPd3AZO5bTo0axf9Mmhg8fzsCBAyvKFJjNTDt8jrTiYhYfeBLX9hP48bsNlBQUcOPL\nbxB0hSzN2xLJtihCNFP5W7eR+NJLOPbvT9EjD7Nx0yZ69erFlVdeWaXcjJPxHMorZG7kGwRZe/HT\n2iRMpaXc8upbEuwtmIS7EM1Q0ZEjxD/+OHbt2+P27jv8umwZPj4+TJw4scpTpL8kZ7IoJYtnkn+l\nR2o0Px9wxMrGlqmvv4Nf+46N2AJR3+oU7kqpa5RSJ5VSZ5RSL9RwPkQptVEpdUApdVgpde3lr6oQ\nAqA0Opq46Q9g7eFB4Kef8uvvv2M0Grn55puxtbWtKBdVWMKMU/EMMiZw/f7vWXymLY7unkx94x08\n2wQ1YgtEQ/jTcFdKWQFzgPFAN+A2pVS3i4rNBH7WWvcFpgKfXO6KCiHAnJND7P3TAQj+fD7bIo8T\nExPDhAkTKpYVACi1WHjwWDS2FhPPbnqN5fFd8QwMZerr7+Dq7Vvb7UULUpc3VAcAZ7TWUQBKqR+B\nycDxSmU04Fr+sRuQeDkrKYQoe0gp8fkZGJOTCV24kDizmYiICPr27Vtl6V6AWVFJHM4v4t0jc9gT\n5Ytfx87c+NIbMoe9FanLsEwgEFfp8/jyY5W9BtyhlIoHVgGP1XQjpdR0pdRepdTetLS0v1BdIVqv\njHnzyN+8Gb8ZMzB17MCSJUvw9fVl/PjxVcqtz8jls7g0bs49Stb2RLwCApjy4usS7K1MXcK9pjU+\n9UWf3wYs0FoHAdcC3yilqt1baz1Pax2utQ738fG59NoK0Urlb9tG2ocf4TpxIm63TWXx4sU1jrOn\nlBh5PDKWTqqE9r/8iIurIze99i/snZwbsfaiMdQl3OOByntpBVF92OU+4GcArfUOwB7wvhwVFKK1\nMyYmkvjMs9h17EDAG6+zadMmYmJiuO6666jcSbJozaORMRSYTFz1yxycDSZufuMDHN3cG7H2orHU\nJdz3AJ2UUu2UUraUvWG67KIyscBVAEqpKygLdxl3EeJ/ZCktJf7Jp9BGI4EffURiZiZbt26lb9++\n9OrVq0rZObGpbMnKZ+y2ZfhnJ3Hzow/g4h/SSDUXje1Pw11rbQIeBVYDkZTNijmmlHpDKTWpvNgz\nwP1KqUPAD8DdWuuLh26EEJco5e23KT58mIC338IqOJhly5bh7OzMuHHjqpSLKizhvXNJdIs/Re9j\n27jp6ra4D7ihkWotmoI6LT+gtV5F2RullY+9Wunj48DQy1s1IVq3nP/+l+wffsTzvntxvfpqNm/e\nTGpqKrfddhv29vYV5bTWPB8ZgyotZczmJdzUJRHvW35sxJqLpkDWlhGiCSqJiiLp/17DsX9/fJ96\nirS0NCIiIujevTtdunSpUnZxcgZbcwsZu3M1f/PYgf/NH4ODbGLd2km4C9HEaJOJxBdexGBnR5t/\n/wttMLB8+XJsbW2rTXvMMpp46dg5AtKTeSp7EaH9h0C3SbXcWbQmEu5CNDEZn39B8eHDBL7/b2x8\nfdm9ezexsbFcf/31ODtXndL4dMQu8pQ9L8ZtIMw9HsYubqRai6ZGFg4TogkpPnGCtDlzcBl/Da7X\nXktOTg7r1q2jffv29O7du0rZ/+47yG8GJ65KPMldBQsh/F7wlsXARBnpuQvRROjSUhJnvICVmxv+\nr76K1pqVK1eitea6666rstpjUnwsM2PScLe3Z7ZehcHOCUbOaMTai6ZGeu5CNBFpn3xCycmTBLzx\nOtYeHhw7doxTp05x5ZVXVtngujA3hxlLV5Lm4cM7XiW4Ry2HYU+Bkzw3KC6QcBeiCSg6dIiMefNx\nu+EGXK68kqKiIlatWkWbNm2q7KpkNpn46pOP2dAlnKvsFZMP/gNcg2DQQ41Ye9EUSbgL0cgsxcUk\nvvAi1n5++L30IgCbN2+msLCQ6667DoPhwo/plh+/5tuAzthaW/Mv+7OQdBCuegVsHBqr+qKJknAX\nopGlffAfSs+do82sN7FycSEjI4Pdu3cTFhZGQEBARbnTu7ez6MBhzoV0ZkY7PwI2vQr+vaDnLY1Y\ne9FUSbgL0YgK9+0j8+uv8bj9dpyGDAFg7dq1WFtbM3r06Ipy2clJrJr7IdtGTiLU3oZ74hdBThxc\n/Q8wyI+xqE6+K4RoJNpoJPm117EO8Mf32WcAOHfuHCdOnGDYsGG4uLgAYCwtYdkHb3OsQy+SXDx5\nKcgV2y3/gk5XQ/tRjdcA0aRJuAvRSDK/+ZaS06fxf/llDI6OWCwWVq9ejZubG4MHD64ot/Grz0iM\nj2PHsGvp6+LIpKNzoDQPxr7RiLUXTZ2EuxCNwJiURNrs2TiPHo3LVVcBcOjQIZKTkxkzZgw2NjYA\nHN20jiMb1pB66wOkacWrAbaovV9An7+B7xWN2QTRxEm4C9EIUt7+J1gs+L38EgAlJSWsX7+eoKAg\nevToAUBazDnWfzEXjz79WebehnHergw+PBe0BUY815jVF82AhLsQDU1D3po1eD/4ILZBQQBs27aN\n/Px8xo0bh1KK0qJCln/wNnZOTpy47g6KLBZm+tnBvgXQayp4hDZuG0STJ+EuRAMzaLBt1w7Pe+8B\nICcnh+3bt9OjRw+Cg8t2tNywYB7Zycn0euQ5vk/P428BXnQ6OBfMpTD86casvmgmJNyFaECqfH8y\n/1dfwVC+sfW6desAGDNmDACndm7l2KZ1DLzhZr5UztgaDDzrawd7voQeN4FXh0apu2heJNyFaCCl\n0dEoDVqBU/lsmISEBI4cOcLgwYNxd3cnLyOdtfNm49+xM9ZjJ7MiLYeHg33x3f8pGAth+DON3ArR\nXEi4C9EAtNYkvzmr7OMLizuyYcMGHB0dGTZsGNpi4fdPPsBsMjH+kaeZFZ2Cr601D/nYwK55ZZtw\n+HZtpBaI5kbCXYgGkLd6DQVbt1YJ9piYGM6ePcvQoUOxs7Nj38qlxB49xKi77me/nSu7cwp4tq0/\nTns/L5vXLjNkxCWQ9dyFqGeWkhJS33sPuy5d0KoYKOvJb9iwAWdnZ/r3709qdBRbfviajv0H0WP0\nWK7df4Zge1umetjAzk+gy7Xg37ORWyKaE+m5C1HPsr75BmNCAn4znq84du7cOWJiYhg+fDgKzcqP\n3sPB1ZWx0x9jQ1Y+B/MKeSLUD9t9X0BxtvTaxSWTcBeiHpkyM0n/9DOcR46sWBhMUzbW7urqSlhY\nGBHffkVmQhzXPPQkDi6u/PtcMkH2NtziaQc7ZkPHMRAY1rgNEc2OhLsQ9Sh99mwsRUX4Pn+h523E\nRHx8PCNGjCDh2GEOrl5B2PhJtO0dxobMPA7kFfJkqD+2+xdAYQaMeL72FxCiFhLuQtSTkrNnyfrp\nZzxuvRW7DmVz0zVQZCjG3d2dbl06s3rex3gGBjP89rvRWvOv8712L0fY/hG0GwEhA//4hYSogYS7\nEPUk9d33MDg44P3oIxXHjJgwYWHkyJFs+W4BBZmZXPPwk1jb2lb02p8I9cP26M+QnyLz2sVfJuEu\nRD3I37aN/M2b8X7oQaw9PQGwWCwUGYqxwoCrNnF04xrCJ00hoGMXtNb8OzqZQDsbbvV1hx1zymbH\ntBvZyC0RzZWEuxCXmTabSX3nXWyCgvC4446K48ePH8eEBQeLHes+n4NnYDBDbrodgI2ZeezPLeTJ\ntn7YntsI6Sdh8GOgVG0vI8QfknAX4jLL+fVXSk6dwveZpzHY2QFlvfZNmzZhhQE7bagyHKO15l/n\ne+3+nrD9Y3BpA91vaOSWiOasTuGulLpGKXVSKXVGKfVCLWVuUUodV0odU0p9f3mrKUTzYM4vIPXD\nD3Ho0weXa66pOH7kyBHS09NxtNiitKViOAYu6rWnHIVzm2HgdLC2baxmiBbgT59QVUpZAXOAsUA8\nsEcptUxrfbxSmU7Ai8BQrXWWUsq3viosRFOW8cXnmNPS8fv4Y1T5kIrFYiEiIgJfHx9sLQqNYmj5\ncEyVsXZ/T1j6Mtg4Qb+7G7EVoiWoS899AHBGax2ltS4FfgQmX1TmfmCO1joLQGudenmrKUTTZ0pL\nI3PBQlyvHY9Dnz4VxyMjI8nIyMCttBAFoKywLl/ud3NWHvtyy2fI5CfD0UUQNg0cPBqnEaLFqEu4\nBwJxlT6PLz9WWWegs1Jqm1Jqp1LqGoRoZdLnzkUbjfg88UTFMa01ERERuDo7kbRzM1oZqrxJ+lFM\nKv62Ntwa4Am755VtoTfwwcaovmhh6hLuNb1dry/63BroBIwCbgM+V0q5V7uRUtOVUnuVUnvT0tIu\nta5CNFmlsbFk/fwL7jffhG3ohS3wTp8+TUpKCoakWLzaBIGyqji3L6eA7dn5PBjsg52xEPZ+CV0n\ngme7xmiCaGHqEu7xQHClz4OAxBrK/FdrbdRanwNOUhb2VWit52mtw7XW4T4+Pn+1zkI0OWkffYyy\ntsb7oYcqjp3vtdtZW2GKj2bsA49Vuebj2BTcra2Y1sYLDn4HxTkw5LGLby3EX1KXcN8DdFJKtVNK\n2QJTgWUXlVkKjAZQSnlTNkwTdTkrKkRTVRwZSe6KFXjeeSc2vhfmEkRHRxMfHw/xUfQZO56grt0r\nzp0sKOb39FzuDfLGyUDZsr5BAyB4QCO0QLREfxruWmsT8CiwGogEftZaH1NKvaGUmlRebDWQoZQ6\nDmwEntNaZ9RXpYVoSlL/8x8Mbm54/f2+KscjNm/GymLGQ2mG335XlXOzY1NwMBi4L9AHTqyErGgY\n8mgD1lq0dHXarENrvQpYddGxVyt9rIGny/8TotUo3LOHgs0R+D77DFaurhXH4+LiOBcdjV1aImPu\newg7R6eKc2kGd35NyeKeQG+8bK3LlvV1Dy0bbxfiMpEnVIX4i7TWpL7/Ada+vnj87W9Vzq1fuwZl\nNnFFp46LszgDAAAgAElEQVR0DK+6quNK++EAPBjsC/F7IW4XDHoYDFYIcblIuAvxF+Vv3ETRgQN4\nP/IIBgeHiuNJiYlEx8bhkJfJ2HurTmvMVU5ssO/PjX6eBNrbws65YOcKff928e2F+J9IuAvxF2iz\nmbQPPsA2NBT3KVXXgFmxZBGYTVx17USc3Ks+jPS7/RCMWPNIiC/kJsHxpdB3Gti5NGT1RSsg4S7E\nX5C7YgUlp0/j8+QTKBubiuOxZ8+QkJaBl5UibOz4Ktfkm8ysth9CeOlxOjvZl81rt5hhwN8buvqi\nFZBwF+ISWUpLSfvoY+y7dcNl3Lgq55b+8D1ozZQ776pYW+a8bxIzKDA4MLl4E5hKysK98zXg2b4B\nay9aCwl3IS5R9i+/YExIwOepp1CGCz9Ch7ZsItNoItTXm8D2HatcU2Kx8FlcGt2NZ+hoioejS6Aw\nHQY+0NDVF62EhLsQl8BSVET6p5/iEN4Pp2FDK46XFhWyevkyQDH5tupvji5KziK51Mjkok1lB3bN\nBZ+u0H5UQ1RbtEIS7kJcgqzvf8Cclo7vE09UGXbZ9MPXFNo70blDOzy9vKpcY9GaT2JT6eXsQE/j\nmbKDSYfKeu2y05KoJxLuQtSROT+fjPnzcRo6FMf+/SuOp0SdYe/+A2CwYuz4CdWuW52ew9miEh4O\n8b2wCp+9G/S6tWEqLlolCXch6ihz4ULM2dn4PHlhSV+L2czv8z7G6OlLp44duXhBPK01s2NTCbG3\nZaJPpYVSw+4EWyeEqC8S7kLUgTk7m8yvFuA85iocevasOH7g9+UkZuehDVaMGDmy2nW7cwrYl1vI\ng8E+WBsqDcH0v78hqi1aMQl3Ieog44svsBQU4PP44xXHctNT2frTd1gCQmjbti3BwcHVrpsTm4qn\njRVTA7zAWHThhEdotbJCXE4S7kL8CVNaGpnffIvrhAnYd+4MlA23rP/yU0qc3TBqGDZsWLXrThYU\nsyYjl3sDfXC0MsDRxQ1dddGKSbgL8SfSP5tXtn3eYxeW5D29eztn9+1GBXfA39+fDh06VLtubmwq\nDgbFPYHeoDXs+rQhqy1aOQl3If6AMTGR7J9+wn3KDRXb55UUFrDxq89wat+FgpIShg0bVu1p1KSS\nUhanZHFbgFfZsr6xOyD5SGM0QbRSEu5C/IH0uXMBqmyft/XHr8nPzka3aYuHhwdXXHFFtevmx6Vj\n1poHgstnz+ycCw4e1coJUV8k3IWoRWl0NNlLfsV96lRs2rQBIPHUCQ6uWUW7kWNIzchg6NChWFlV\nXYc912Tm68R0Jvm6E+pgB9mxcGIFhN1V08sIUS8k3IWoRdrHs1G2tnhPL5u2aDaZWDd/Ns4enuQ5\nuOLk5ETv3r2rXfd1Qjr5ZgsPh5Tvp7rnc0BBf1n9UTQcCXchalB88iS5K1fiOW0a1uUPJu1buZS0\n2Gh633gb56KjGTx4MDaVlvuFsgXC5senMcLDmV4ujlBaAPsWwhUTwb36VEkh6ouEuxA1SPvwIwwu\nLnjddy8AOanJ7Fj0Ax37DyImIwc7OzvCw8OrXbc4JYuUUhOPhPiVHTj8MxRnw8CHqpUVoj5JuAtx\nkaKDB8nfsAGv++7Dys0NrTXrPv8EZTDQe/KtREZGMmDAAOzt7atcZ9GaubGp9HB2YISHc/n0x8/A\nvxeEDGqk1ojWSsJdiIuk/udDrLy88Jx2BwAnt0cQfWg/w6ZO4+CxY1hbWzNoUPWwXp2ew+nCEh4N\n8S2bGhm1CdIiYdBDsvqjaHAS7kJUUrBjB4U7d+L9wHQMTk4U5+ezceF8/Dt0ot2gYRw+fJh+/frh\n5FR10S+tNR/FpBJaeYGwXZ+Bozd0n9IILRGtnYS7EOW01qR+8B+sAwJwv7VsOd6I77+iKC+XsdMf\nY8eOnQAMGTKk2rXbsvM5kFfIIyG+ZQuEZUbBqd8h/F6wsa9WXoj6JuEuRLn8jRspPnwYn0cexmBn\nR3zkUY6sX03YtZNx9PZl//799OrVCzc3t2rXzo5JxcfWmlv8PcsO7J4PBquycBeiEUi4CwFoi4W0\nD/6DbWgobtdfj6m0lDWffYyrjx9Db/4bu3btwmQy1bhA2OG8QjZl5TE9yAd7KwOU5MGBb6H7DeAa\n0AitEULCXQgAcleuouT0abwffwxlbc3OJT+RlZTA2OmPYgZ2795Nt27d8Pb2rnbt7NhUXKwM3BVY\nfu7g91CSCwMfbNhGCFGJhLto9bTRSNrHH2PXpQuu48eTFnOOPcsW0X3kVbTt1Zc9e/ZQUr5A2MWi\nCktYkZrN3YHeuFpbgcVS9kZqYDgEVZ8HL0RDkXAXrV724iUYY2PxefIJNJrVn36EvbMLI6fdh9Fo\nZOfOnXTo0IE25evLVPZJbCo2BsX9QeULhJ1ZB5lnpdcuGp2Eu2jVLIWFpM2ZjUNYGM6jRrF/1TJS\nok4z+u7pOLi4cuDAAQoKChg+fHi1a5NLjPycnMmt/p742pUvQ7BrLjj7Q7fJDdwSIaqqU7grpa5R\nSp1USp1RSr3wB+VuUkpppZT8PSqahYwFCzCnpeP77LPkpKaw7advad9vAF0GD8dsNrNt2zaCg4MJ\nDa2+Ld68uDRMWvPI+QXCko/C2Q0w4O9gbdvALRGiqj8Nd6WUFTAHGA90A25TSnWroZwL8Diw63JX\nUoj6YMrIIPPzL3AZOwaHvn1YO382BisDY+57GKUUhw8fJicnp8bNOLKNJhYmpjP5/LK+ADtmg40T\nhN/XCK0Roqq69NwHAGe01lFa61LgR6Cmvzn/AbwLFF/G+glRb9I/mYulpASfp57meMQGYo8cZPjt\n9+Di5Y3ZbCYiIoKAgAA6l++bWtnChAwKzBYeDS1fICwnAY78AmHTwNGzgVsiRHV1CfdAIK7S5/Hl\nxyoopfoCwVrrFX90I6XUdKXUXqXU3rS0tEuurBCXS2lMDFk//YT7jTdi8nRn08L5BHbtRu8x1wBw\n+PBhsrKyGDVqVLVee4HZzLz4NK70dKG7s0PZwd2fgbaUrSMjRBNQl3CvacUjXXFSKQPwAfDMn91I\naz1Pax2utQ73KV8jW4jGkPbhhygbG7weeZj1X8zFWFLM2OmPoQyGP+21f52QQYbRxJPne+3FubD3\nq7I3UT3aNmxDhKhFXcI9Hqi8y0AQkFjpcxegB7BJKRUNDAKWyZuqoqkqOnKE3FW/4Xn3XUSdOcHp\n3dsZcssdeAWWfZsfOXKErKwsRo4cWa3XXmi2MCc2leEezgxwdy47eOCbsoeWhjzW0E0RolZ1Cfc9\nQCelVDullC0wFVh2/qTWOkdr7a21bqu1bgvsBCZprffWS42F+B9orUn917+x8vDA/qYbWf/lpwR0\n6kL4dTcAYDab2bx5M/7+/nTp0qXa9V8npJNuNPFMW/+yA2Zj2ebXoUMhsF9DNkWIP/Sn4a61NgGP\nAquBSOBnrfUxpdQbSqlJ9V1BIS6ngi1bKNy1C6+HHmLD9wswlZRwzcNPYTCUbXJ9vtde01h7odnC\nnLhUhrk7M+h8r/3YUsiJk167aHKs61JIa70KWHXRsVdrKTvqf6+WEJefNptJ/de/sQkJISXYn7Mr\nf2LkHffi2SYIoGKsvbZe+7eJ6aSVmpjXvbzXrjVs/wi8O0OncQ3ZFCH+lDyhKlqNnGXLKTl1Cqfp\nf2fjN1/QpvMVhE24MKv36NGjZGZm1jjWXmS2MDs2lSHuzgw+32s/FwHJh2Hwo2CQHyXRtMh3pGgV\nLAUFpH3wAXY9e7Dz7HFMRiPjHnqyYjjm/Fi7n58fXbt2rXb9t4kZpJaaeKat34WD2z8GJx/odWtD\nNUOIOpNwF61C+mfzMKWmknPdeKL272HY1DvxbHPhcY3zvfaaxtqLzRZmx6Yw2N2JoR4uZQdTjsOZ\ntTDgAdlpSTRJEu6ixSuNiyPzq6+wmTCebRtXE9i1O2Hjr6s4f36s3c/Pr+ax9qQMUkorzZAB2DEH\nrB2gvyw1IJomCXfR4qW88w7axoYDzjaYzSbGPfQEqtIY+ZEjR8jIyGDkyJEYLho7LzZbmB2TyiA3\nJ4aeH2vPjoPDP0HfO2SpAdFkSbiLFq1g+3by160ndfwYYk8cY9S0v+Phf2FddqPRyMaNGwkICKhx\nrP27pAySS408287/wnDN1g/K/h36REM0QYi/RMJdtFjaZCL5rbcobBvCgeiTdAgfRK/ytWPO27Nn\nDzk5OYwdO7bmXntsKgMr99pz4sueSO17B7gHI0RTJeEuWqysH36k6GwUB0N9sXdx5eoHHqvyZmlR\nURERERF06NCB9u3bV7t+QUI6SSVGnm17Ua9daxj+dEM1Q4i/RMJdtEimrCzSPv6Y02Hdyc7OYvzD\nT+Po6lalzNatWykuLmbs2LHVrs82mvhPTAqjPFwY7lk+QyYnAfZ/DX1uB/eQhmiGEH+ZhLtokdI+\n+ogkg4UoUxHh100htFefKudzcnLYuXMnvXr1wt/fv9r1H8WkkmMyM7NDwIWD2/5Ttqzv8D9dAFWI\nRifhLlqc4hMnSF68mKPtA/Ft24FhU6dVK7Nx40YARo8eXe1cfHEpXySkcZO/Bz1cHMsO5ibBvoVl\nvXaP6lvuCdHUSLiLFkVbLCT/400Ot/XHbDBw7ePPYmVtU6VMSkoKhw4dYsCAAXh4eFS7xzvnkgCY\n0e7iXrtZeu2i2ZBwFy1K9uLFHIs9Q7qDLaPvnl6xRntl69atw9bWluHDh1c7dyy/iEXJWdwX6EOQ\nffkm17lJZZtx9J4qm3GIZkPCXbQYxtRUTnz4ASfbeNOx/2B6XlV9pcbo6GhOnz7N8OHDcXR0rHb+\nzbOJuFlb8Xio74WD2z4EiwmGP1uf1RfispJwFy1GzBtvsM/XFVcv77KnUC9aI0Zrzdq1a3F1dWXg\nwIHVro/IzGNjZh5PhPrhblO+GnZeMuz7CnrfBp7tGqIZQlwWEu6iRchZu5atUZGY7GyZ9Pyr2Ds5\nVytz/PhxEhISGD16NDY2VcfhLVrzj7OJBNnbcE+g94UT2z4q221phIy1i+ZFwl00e+b8fDb95z0y\nXBwYe/+j+Lat/kBSaWkpa9aswdfXl969e1c7/2tKFkfyi3ihXQD2VuU/FnnJsPfLsiV9PavfU4im\nTMJdNHv7X3uFM862dOs7gO5XXl1jmYiICHJycpgwYUKNywy8fS6JHs4OTPGrNHtmwz/KxtpHPlef\n1ReiXki4i2Ytcd1atkefxNPBibFPv1BjmdTUVLZv306fPn0IDa0+R31+fBrxxUZe6dAGw/lx+qRD\ncOA7GPiA9NpFsyThLpqtkrw8Vn76IUopbnj9HaxtbauV0VqzatUqbG1ta1xmIKaohPejkxnv7cbI\n88sMaA2rXy5bzneE9NpF8yThLpolrTWrXnqGXINm7ORbcA9tW2O5I0eOEB0dzZgxY3Bycqp2j5dO\nJaCU4s1OF3Zl4sRKiN4Co14EB/d6bIUQ9UfCXTRLO+Z/QlRqIj3cfLjijrtqLFNUVMTq1asJDAwk\nLCys2vmVaTmsz8zl+bb+BJ5/YMlUCmtfAZ+u0O+e+myCEPVKwl00Oyc2b2DH+t8ILDJy5Tvv11pu\nw4YNFBYW1vgmar7JzMzTCXR3tufvQT4XTuyeB5lRcPUssLKuryYIUe/ku1c0K4mnTvD73A/wKChm\n/AuvYeNZ8zZ3iYmJ7NmzhwEDBtCmTZtq5985l0RKqZEve7TF2lD+JmpBBmx+FzqOgU5j6rMZQtQ7\n6bmLZiMnNYWls17BrriUMcPG4DZsWI3lLBYLK1aswMnJiSuvvLLa+cN5hXwRn86dbbwIc6s0Dr/p\nbSjNL+u1C9HMSbiLZqG4IJ8ls2ZiKihgqI0Lwc/U/sTovn37SExMZNy4cdjb21c5Z9aa50/G42Vr\nzUvtK636mHqi7IGl8HvAt/peqkI0NzIsI5o8s8nE8vffJis5iQFJWXT9fi7qouUDzsvKymLt2rW0\na9eOnj17Vju/MCGdg3mFzO0WiptNpW//NTPB1hlGvVRfzRCiQUnPXTRpWmvWfzmX2KOH6BmbSo9n\nZ2Bbw4NIAGazmSVLlqCUYvLkydUWDkspMfJ2VBIjPJy53rfSFMdTa+DM2rInUZ286rM5QjQYCXfR\npG3/5TuOrF9Nh9Rsug0Zgdv1k2stu23bNuLi4pgwYQLu7lXnp2utmXEqjlKt+Wfn4AvBX5wDK54s\nm/o4YHp9NkWIBlWncFdKXaOUOqmUOqOUqvaMt1LqaaXUcaXUYaXUeqWU7EMm/me7lv7CzsU/ElJs\nprvBHv/X/q9ab/y8hIQENm3aRI8ePWoejknM4Pf0XF5sF0B7R7sLJ1a/DHlJMPkTsLardp0QzdWf\nhrtSygqYA4wHugG3KaW6XVTsABCute4FLALevdwVFa3L/t+WsfWHhYRYO9D9bDyB772HlYtLjWVL\nSkpYvHgxzs7OTJgwodovgMj8Il47k8BoTxemB1ea0356HRz4BoY+AUH96rM5QjS4uvTcBwBntNZR\nWutS4Eegyt/GWuuNWuvC8k93AkGXt5qiNTm8fjUbF8wj2M2LbvuOEvDKKziG9a21/OrVq8nMzGTK\nlCk4ODhUOVdktvDQ8Ricraz46IqQCwuDFefA8sfLhmNGvVifzRGiUdQl3AOBuEqfx5cfq819wG//\nS6VE6xW5ZSNr588myD+Qblt24zVtGh633lJ7+chI9u/fz9ChQ2nbtm2182+cTeREQTEfXRGCj22l\nGTarXypbr/16GY4RLVNdpkLWNMipayyo1B1AODCylvPTgekAISEhdayiaC1O79rOb598QJugULqv\n2YLrkKH4zXi+1vJ5eXksW7YMf39/Ro8eXe386vQcvkpI54FgH670cq30QmvhwLcw7GkIlOEY0TLV\npeceD1TeQj4ISLy4kFJqDPAyMElrXVLTjbTW87TW4VrrcB8fn5qKiFbqzJ6drPjwXfyC29J7+34c\n2rQh8IP3UdY19z8sFgtLly7FaDRy4403Yn1RuaSSUp46EUtPZ4eqDysVZcOyx8HnChhV8/rvQrQE\ndQn3PUAnpVQ7pZQtMBVYVrmAUqov8BllwZ56+aspWrKjG9ey7P238AlpS7/Ic1iZzATN/QQrV9da\nr1m/fj1nz55l3LhxXNxRMGvNY8djKTJr5nYPxa7yomGrX4b8FBmOES3enw7LaK1NSqlHgdWAFfCl\n1vqYUuoNYK/WehnwHuAM/FI+UyFWaz2pHustWog9yxYT8d1XhPbsQ1hcGiVR5wiePw+7du1qvebg\nwYNs27aN8PBwwsPDq52fE5vK1ux83u8aTEfHSssPnFoDB7+F4c9AYPUlgIVoSeq0/IDWehWw6qJj\nr1b6WJbQE5dEa03Ed1+xd/kSugweTp/sQvI2R+D3ykychgyp9bqYmBiWLVtGu3btGD9+fLVpj7+l\nZfN2VBLX+7pzm3+lFSMzo2DJ/eDbHUbOqK9mCdFkyNoyosFZzGbWfPYRxzavp/fVE+iWnEnOjz/j\ned+9eNx+e63XZWVl8dNPP+Hh4cEtt9yClZVVlfP7cwt4+HgMfV0deb9ryIXgL8mHH/9W9vHUb2U4\nRrQKEu6iQRlLS1jxwT+J2r+HwTfdTvtz8WR98y0ed07D99lna30Ctbi4mO+//x6LxcLtt99ebT57\nTFEJdx4+h4+tDQt7tsPRqnyc3WKBXx+AtBNwxxLZ7Fq0GhLuosHkZaaz/N9vk3T2FFfd+xBtTpwh\n48uvcL9tKn4vvlhrsJvNZhYtWkRGRgbTpk3Dy6vq4l7ZRhN3HI7CpDXf9WpfdT57xHtwYgWMews6\nVJ8uKURLJeEuGkR85FGWf/BPjCUlTHr6Rdz3HiR97qe433wT/q+8UmuwA6xdu5YzZ84wceJE2l30\nRmupxcK9R6OJKSrlpz4d6ORU6Q3UyBWw6S3ofRsMeri+miZEkyThLuqV1pqDa1ayaeF83Hz9uOXV\nt9C/rSbto49xmzwZ/9dfRxlqn5EbERHBzp07GTRoULWZMVprnj4Rx/bsfOZcEcJgd+cLJ1Mjy4Zj\n2oTBxP/AH/zyEKIlknAX9cZYWsL6zz/h2Ob1tA/rzzWPPE3BDz+S9u/3cZ0wgYC3Zv1hsG/atIlN\nmzbRq1cvrr766mrn34tOZlFKFjPa+XNj5ZkxhZnww21g6wRTvwMb+2rXCtHSSbiLepGbnsqyf79F\nStQZBt90G4Mm30zKrLfI/uknXK8dT5t3/om6aLbLeVprNm7cSEREBH369GHSpEkYKv0S0FrzfnQK\n70enMNXfkydD/S5cXFoIP98JuQlw90pwrb45thCtgYS7uOxObI9g/RdzsZhNTH7uFdp17U78w49Q\nsHUrXvffj89TT9baY9das379erZu3UpYWBgTJ06sFuyvn03k07g0bvX35F9dKm28UVoIP9wK0Vth\nyjwIHtAQzRWiSZJwF5dNUV4u67+Yy8kdW/Dv2JnxjzyDC4qY226nJCoK/3+8gcfNN9d6vdaaNWvW\nsGPHDsLDw7n22murBLtZa2acjOfbpAzuC/TmH50CLyzhW1oA398KMdvghs+gV+0rSQrRGki4i8vi\n3IG9rP7sI4pycxh66zQGTL6JkuORnHv4IXRRMSHz5/3hk6daa37//Xd27drFgAEDqj19arRoHouM\nYWlqNk+G+jGjnX/Vh5S+vwVid8AN86BX7b9AhGgtJNzF/6S0qJBN33zBkfWr8Q4O5YYZ/4dfuw7k\nrllD4vMzsPb0JPjLL7Hr1Kn2e5SWsmLFCg4fPsygQYMYN25clWAvNluYfiyaNRm5zGwfwKOVx9hL\n8uG7myFuJ0yZDz1vqs/mCtFsSLiLv0Rrzdm9u9i4cD656an0n3QjQ265A4PJRNLrr5P9w4/Y9+pF\n8CdzsPb2rvU+55cUSE5OZvTo0YwYMaJKsOcYTdx3NJpt2fn8s3MQdwdWuldJXnmw74YbP4ceN9Zn\nk4VoViTcxSXLiI9l48L5xBw+gFdQCLe+9k+Cunan6NgxEp97ntKoKDzvuQefp57EYGtb632ioqL4\n5ZdfKpYU6Ny5c5Xzh/IKuf9oNIklpXx8RQg3VZ7umJ9atl5Mwj646QvofkN9NVeIZknCXdRZcUE+\nO375ngOrV2Dr4MDou6fTe2zZm54Zn39O6ocfYe3hQciXX/zp+PqOHTtYu3Yt3t7eTJ06tcqSAlpr\nFiZm8OrpBHxsrflv3070c3O6cIO43WXTHYuy4eYF0E1WlxbiYhLu4k9ZzGaOblrL1h++pig/j15X\njmPo1Gk4urphTEoifsYLFO7ejcvVV+P/+mtYe3jUeq/S0lKWLVvG0aNH6datG5MnT8bO7sIqjfkm\nM8+ejGNpajZXebrycbcQPG3Kv021hj2fw+8vglsg/H0t+Pes7+YL0SxJuItamU0mIrduYtevP5Gd\nnERg126MvvsB/Np1QBuNZC5cSNrsOWA2EzBrFm5TbvjDNWLOnTvH8uXLycrKYsyYMQwdOrRK+cj8\nIv5+NJpzRSW81D6AR0N8K011LIQVT8Lhn6DTOJjyGTjU/ktEiNZOwl1UYzYZOR6xkV1LfyYnJRmf\ntu2Z9MxLdOw/GKUU+Vu2kPL2PymNisJp2DD8X30F2z/Y8LywsJC1a9dy4MABPDw8uPPOO6ssAGa0\naL5MSOPtqCRcra1Y1KcjQzwqrROTGQU/TYOUYzD6ZRj+LPzBsgVCCAl3UYnJaOTYpnXs/u8v5Kal\n4te+E6Ofv5/2YQNQSlESdY6Ud/5JweYIbENDCfp0Ls4jR9baW9dac+zYMX777TcKCwsZOnQoI0eO\nxLbSm6w7svN58VQ8JwqKucrTlf9cEXxhyV6LBQ58DWteLVv462+LoJNs+iVEXUi4C3JSkzm07neO\nblhDUV4uAR27MOa+h2nbpx9KKUyZmWR8No/M777DYG+P7/PP43nH31B/MBMmOzublStXcvr0adq0\nacMdd9xBQEBAxfmUEiP/OJvIopQsguxt+KpHW67xdrvwiyLxIKx8BhL2QsgQuGEueLSt56+EEC2H\nhHsrpS0Wog/t5+CalUQd2ItC0SF8IH3GTSCkR2+UUpTGx5P55VdkL1mCLinB/aYb8XniiT+ct56X\nl8fWrVvZu3cvBoOBcePGMWDAgIot8UzlQzDvnkum1KJ5KtSPx0L9LuycVJQNG2eVvXHq6FW+lMCt\nsmSvEJdIwr2VyU5O4sT2CI5uXENOagqObu4MmnIrPa8ch6u3DwDFJ0+SMf9zcn/7DQwG3CZdh9d9\n92HXvvYt6s6H+r59+zCbzfTp04eRI0fi7u4OQInFws/JmXwSm8q5olJGe7owq1MQ7R3LZ8poDYd/\nhjUzoTAdwu+DK2eCg3u9f02EaIkk3FuBvMx0Tu3Yyoltm0k+exqAoG49GH773XTsPwgraxu0yUR+\nRASZ335LQcQWDI6OeN75/+3de3Bc1X3A8e9v3yvrLVvCD8l2/LYxtrEHY2LM2xhKA3Sg4ZGGTNLS\nTJKZJNPxJEwyLemkMzBuM+mUzDBuoaElEDDFrhNoDQy4SSkQhGQwtiwsP/Sy9bBWq13t8+7eX/9Y\nWcgPyStL8lqr85nR7L1nz73nd658f3N87r26X6X8ka/ivuKK4fcdDvPuu+9SW1s7mNSvv/56yssz\nDxz1p9I8d6KH7a1ddCZTrCry89zK+WyuKM5MwaRTcGAnvPtz6PwUZq+Fh3fArNWX5NgYRr4yyT1P\n9XV1cLS+ls/e+1/aDh0AVSrnL2DTV77Okg0bKZ5eCUCiqYlTO3cS2v0bUt3dOMvLmfG971L24IM4\nS0rOu2/btjl+/Dh1dXU0NDRg2zarVq1i06ZNg0m9O2nxTNsp/rX9FH2pNJvKCnlqWRUbywozST0Z\nhfrn4b1/gmALTF8C9zydmYIxd8IYxpiZ5J4nUpZFe8MBju37kGP1HxE40QZA+aw5XHffQyy5bhPl\ns2Zn6vb2EnjhBfp27iK+fz84nRTecAMl995D0Q03DHuhNBwOs2/fPurq6ujt7cXn87F27VrWr19P\nRbl8qFcAAA39SURBVEUFCdvmt11BXukM8FZPiLTCnTNK+E5NFWuKCzI7iZyC2mfhg6ch2gPV62HL\nk7B4i0nqhjGOTHKfpNKpFF3Hj9DecIDWhk9p/fQTrEQcp8vFnOUruerWO5i/Zi1lMzMJPdnUxKnt\n/0z/3r3E9u0D28a7eDGVP/wBJXfdNexF0lgsxuHDhzl48CCNjY2oKnPnzuWmm25i2bJluFwuakNR\nnmxsZXdXkGAqTaXHxV/MmcFDMysyL6xORuCTHbD/ZTjyNtipzINIG78PczdcysNmGFOGSe6TRCIa\nofNoE20NB2g/dIAThw+RSiQAKK2ayfJNNzN/zTpqVlyF2+cjHQ4Tq6+n85lf0r93L1Z7OwC+5cuZ\n/s1vUnTrLXiXLTvvPeqBQIDGxkYaGxtpbm5GVSksLOS6665jzZo1FJSV83/Bfl461slbPSGa40n8\nDuHOGaXcV1XG9WVFuDQFR96B/Tvg0GtgRaB4Dmz4Nqx6CCqXXtLjZxhTjUnul6FYOETnsSN0HTsy\n8NlEsONk5ksRZsydz8qbNjN76QpmL11OYVk5VmcXsbqP6Hn974nW1ZFobATbRnw+pm3YQMWjj1J4\n4w24q6rOaS8cDtPS0kJLSwtHjx6lu7sbgMrKSjZu3MjixYuJllawN9jPM219vHegnYSt+B3ChtJC\nvj+viruml1AYaoYjL8HRvXD89xDvA19p5uUZK/8UajaYqRfDuERMcs8RVSXaFyTQ3kpPexs9bS0E\n2lsJtLfS3xsYrFdSWUXl/AVceeNtVM5fwMxFS3CGwsQbDhGv/4TgizvoOHgQqy0zxy5+P/7Vq5j+\nrW9RsPZq/GvW4PD5BveXTqfp6emhra1tMKEHApn2XC4X1dXVLLl6LdGZ1XymTp4JRdjXHCTQ1APA\nogIvX5s1nZvKCrlWu/GdrIM/7IWj/wN9rQNB18DyuzPz6AtvA9fwDzsZhjExTHKfQFYiTn+gh77u\nLvo6Owh2nqSvq4O+zk76ujpIRCODdd0+PxVzqpl71Roq5tRQWTOPMq8fR1c3yePHSe4/RGLX67Qc\nOkS6r+/z7ebW4FuxgrKHH6Zg3Vp8S5cibjeqSiQS4fiJE3R2dtLZ2UlHRwfd3d2k02kAPAUFFMxb\nQOGaawgXl9PscPHf0QRN0QQc60aAxdN8bKkoZq0zyvXRRmq6aqFhH5z8GBKhTBC+Epi/KTOH/oUb\nofwL5qEjw8gxk9wvQjIeIxoMEukLEg0FM8vBXvp7e+jvOUV/oIdwoId4f/iM7ZwuF8WVV1BaWcXM\nhYsoLiik2OWhKKV4+vpInezAqjuItfsNom1tRFOpwW0dRUV45s+naPNmvMuW4lu6DJk3j4idJhQK\n0dnbS6Cjg96GBgKBAIFAgGQyieVwEvYVYJVVYFfWkFh4FWFfAR0OFy3JNJYqpIBAlBqPg2Uui/u9\n3VwdPczqnjqKTh3M/OGudHKgE16oWgEr78/ciz5zdWbd4bx0vwDDMC4oq+QuIluAfwScwL+o6hNn\nfe8F/g1YC/QAX1bV4+Mb6vhR28ZKxEnG4yRjURKRCIlohEQ0OvAZIRGJEO8PEevvJx4OEesPE+8P\nEwuHBi9kns1fWERhYREF/gKmz6qhwOnEZ4M/YeGPxHAHerEPHif1+4+w+/szsQAD418C3kICJdMp\nXjiPiltuxpoxA6usjOS0aXzSEebDppM403Hch0/ibD6J7XaScHuIuz1EPV5iXj+pwlIS82YRWeQl\n5HDRz5kjaK+mKeoLUhE5xeZYK8siR1gdaeCa6EGKNTpYz8LFcbuKDx2zaOYOGqwquqYt5d7bb+Hu\ntfPO2/9d9e1s29NIezCGU4S0KmUFblShL2Yxq9TP1tuXcM+a2WP9FQLw4137efGDVtKqOEV4cH01\nP71n5Zjj2FXfzuO7DxCMWYNlpX43j39pxYjbbNvTyIlgbMz93FXfzk9+c4DeqJVV26Pd93jFaVze\nLpjcRcQJ/AK4DWgDPhSR3ap6cEi1bwC9qrpQRB4AngS+PBEBn2ptpqPpM1KWRdpKkkomSVkWqWSC\ntGVhJRKkkgmsRBwrHseKx0glEiTjscx6Io6VTF6wHRHB43ThdbrwiAMPQoEqbtuJS7y4khaueBJJ\nJHDGEjgB29FHyukk6faQ8Hjo8nhIejzEi0tIlhSTWFhJ4qoCEj4/ca+HuMtLrwpBW0k7najTQcrp\nxHK6sJwuki43ybQLK+okWTaDxIZZJFxudJiLkh7boiodpsrqZHH0FDPi3cyOtVMdP0lNvIPq+Emm\nW0EEiKubE1rBSa2gg3Ke080Dy9Np1pkctytIc9ZoPAQf7DqEOtznJIRd9e089up+YlZmyietCjCY\noADagzEee3U/wJgTyo937ef591sG19OqPP9+C8e6+6lr6bvoOHbVt7N1x8dYtp5RHoxZbN3x8bDb\nDO37WPq5q76dra98jJX+vP2R2h7tvscrTuPyl83I/RqgSVWPAojIr4G7gaHJ/W7g8YHlV4CnRERU\n9cwzZBxsf/kl3imdDSIoAgKKEyhARcAH+DOj1cz3QiYIQSXzgwj2wLIOKVcRbHEMLqcdjsEy2/H5\nd/bpsqGfjtOfY5uecNkpfOkEhakoRekoFakAJckwJal+itIRyqwQpakQZVaI8lSIackIPiuOJ5nE\npQ6WzZuTmQP3FUN5CRQsgGnXwrTp/OWrzRxKeAloMWH8wOjnxWNWmm17Gs9JBtv2NA4mjYvZfrRe\n/KD1vOXvHgmctzzbOLbtaTwnsZ9m2Zp134fbv83Ip8S2PY1nJPYLtT0ao4nTmPyySe6zgaFnUhuw\nfrg6qpoSkT6gAjg1tJKIPAo8ClAzwssdRhKqKOfwzIUMpmxVBlO5kikfKJOBMlAcqsiQT1HFgeJQ\ne7DcgeImjRNFsHGp4lQbp23jTNuZZRSXpnFj4yaNS+2BZRsvabyk8UkaLzY+sfGKzTSHTYHY+B3g\nFfA6BJ+AzyHsrG1D0jZ2SsGGlLqwcJFQN//w8LXgquArz7XSpR5aqCDKbCLqI4aXKN4zRtcCHHvk\nj4Y9dm+EXrtAasnOiWAsq7LRbD9a6XEYN1xMP0azzfnKw6mK89TMrv2xHrfRxGlMftkk9/MN784+\ns7Kpg6puB7YDrFu37qLOzie+9uc8ceFqk8aj9W/THj735Jpd6ocVNwNwrNhDexYn4KxS/wW/z2Y/\nF9POaPZ9oTizcXoufSwuph+j2eZ8de9+4IULxjRc+2M9bqOJ05j8snmipA2oHrI+BzgxXB0RcQEl\nQHb/P57itt6+BL/7zKkcv9vJ1tuXnFHH7Rh5CuXsbbJt62xOh4zY1nDtZLPvbOPMxoPrq89b/sUF\n5WOKY6Rj7XZI1n2/2H5uvX0Jbue57Q/X9mj3PV5xGpe/bEbuHwKLRGQ+0A48ADx0Vp3dwCPAe8B9\nwNsTMd+ej07PdY50B8Pp5bPv4HAI2JoZ5Wdz18PZbZUWuIlbaWKWDUBZgZu/+eMVZ9Qp8bsRgWB0\n5LtMhu77Utwt89N7VgKM+90ywx3rke5YyeZ3mK3T20zE3TLjGadx+ZNscrCI3An8nMytkM+q6t+J\nyN8Ctaq6W0R8wL8Da8iM2B84fQF2OOvWrdPa2toxd8AwDGMqEZGPVHXdhepldZ+7qr4OvH5W2V8P\nWY4D9482SMMwDGNimL/iZBiGkYdMcjcMw8hDJrkbhmHkIZPcDcMw8pBJ7oZhGHnIJHfDMIw8ZJK7\nYRhGHsrqIaYJaVikG2i+yM2nc9YfJctjU6WvU6WfMHX6OlX6CZe2r3NVdcaFKuUsuY+FiNRm84RW\nPpgqfZ0q/YSp09ep0k+4PPtqpmUMwzDykEnuhmEYeWiyJvftuQ7gEpoqfZ0q/YSp09ep0k+4DPs6\nKefcDcMwjJFN1pG7YRiGMYJJm9xFZJuIHBKRT0Rkp4iU5jqmiSIi94vIARGxReSyuiI/HkRki4g0\nikiTiPww1/FMFBF5VkS6ROTTXMcykUSkWkTeEZGGgX+33811TBNFRHwi8gcR+Xigrz/JdUynTdrk\nDrwJXKmqVwGfAY/lOJ6J9CnwJ8Dvch3IeBMRJ/AL4A5gOfCgiCzPbVQT5pfAllwHcQmkgL9S1WXA\ntcC38/h3mgBuVtVVwGpgi4hcm+OYgEmc3FX1DVVNDay+T+bdrnlJVRtUtTHXcUyQa4AmVT2qqkng\n18DdOY5pQqjq75gC7xZW1ZOqWjewHAYagLx8l59m9A+sugd+LosLmZM2uZ/l68B/5ToI46LMBlqH\nrLeRp4lgKhKReWRev/lBbiOZOCLiFJF9QBfwpqpeFn3N6jV7uSIibwFXnOerH6nqfw7U+RGZ/wb+\n6lLGNt6y6WuekvOUXRYjH2NsRKQQ+A/ge6oaynU8E0VV08Dqget+O0XkSlXN+XWVyzq5q+qtI30v\nIo8AdwG36CS/p/NCfc1jbUD1kPU5wIkcxWKMExFxk0nsv1LVV3Mdz6WgqkER2UvmukrOk/uknZYR\nkS3AD4AvqWo01/EYF+1DYJGIzBcRD/AAsDvHMRljICICPAM0qOrPch3PRBKRGafv1BMRP3ArcCi3\nUWVM2uQOPAUUAW+KyD4ReTrXAU0UEblXRNqADcBrIrIn1zGNl4GL4t8B9pC58Payqh7IbVQTQ0Re\nBN4DlohIm4h8I9cxTZAvAn8G3Dxwbu4TkTtzHdQEmQm8IyKfkBmovKmqv81xTIB5QtUwDCMvTeaR\nu2EYhjEMk9wNwzDykEnuhmEYecgkd8MwjDxkkrthGEYeMsndMAwjD5nkbhiGkYdMcjcMw8hD/w+v\nUBPUN3m18AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f166087b2e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## ADD CODE HERE\n",
    "# Calculate the Laplace approximation of the posterior for theta, \n",
    "# draw a few samples and plot the corresponding likelihood functions \n",
    "# for each one.\n",
    "m = np.zeros((D, 1))\n",
    "S = 5*np.eye(D)\n",
    "nb_samples = 5\n",
    "\n",
    "theta_map, S_post = get_posterior(x, y, m, S)\n",
    "print(S_post)\n",
    "# print ('theta_map',theta_map)\n",
    "# print ('S_post',S_post)\n",
    "plt.scatter(x[:,1], y,label='data')\n",
    "for i in range(nb_samples):\n",
    "    t=np.random.multivariate_normal(theta_map.ravel(),S_post)\n",
    "#     t=stats.multivariate_normal.pdf(x,theta_map.ravel(),S_post)\n",
    "    plt.plot(x_test,predict(x_test,t))#,label='prediction {}'.format(i)) ## <--EDIT THIS LINE\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Comparing posterior approximations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Laplace approximation is part of a family of methods known as _deterministic approximate inference_. In addition, there's another set of methods known as _stochastic approximate inference_ which, as you can guess includes most of the sampling techniques you have studied.\n",
    "\n",
    "You must be an expert in sampling by now. Let's actually go and check whether this Laplace approximation we just made is legit.\n",
    "\n",
    "* What sampling methods do you know to sample from an unnormalized distribution?\n",
    "\n",
    "For example, let's try the Metropolis algorithm.\n",
    "\n",
    "1. Write a proposal function to move in $\\boldsymbol\\theta$-space.\n",
    "2. Write a function to accept or reject new proposals based on the Metropolis criterion.\n",
    "3. Write a loop and run the Markov chain for a few thousand iterations.\n",
    "4. Check that the sampling worked: did the Markov chain mix properly? What's the acceptance rate? How does it depend on the proposal function?\n",
    "\n",
    "**Task 7**\n",
    "\n",
    "* Write a function to sample from the true posterior $p(\\boldsymbol\\theta | \\boldsymbol X, \\boldsymbol y)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's plot the results and see if both inference methods arrive at roughly the same posterior.\n",
    "\n",
    "In the same axis, plot\n",
    "\n",
    "* The histogram pdf of the MCMC samples (you may want to look at the `density` option in `plt.hist`); and\n",
    "* The Laplace posterior.\n",
    "\n",
    "Make one plot for the intercept ($\\theta_0$) and one for the slope ($\\theta_1$). What do they look like? Do they match? What kinds of posteriors do you think the Laplace approximation will be good or bad at approximating?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## EDIT THIS FUNCTION\n",
    "def metropolis_algorithim(x_intial,S,nb_iter,prop_fn):\n",
    "    np.random.seed\n",
    "    sample = [x_intial]\n",
    "    aceptence_rate = 0\n",
    "    for i in range(nb_iter):\n",
    "        q_cand = np.random.multivariate_normal(x_intial, S)#.reshape(-1,1)\n",
    "#         q_cand=q_cand.reshape(-1,1)\n",
    "        q_back = np.random.multivariate_normal(q_cand, S)#.reshape(-1,1)\n",
    "        u = np.random.uniform()\n",
    "#         print('q_cand.shape',q_cand.shape)\n",
    "#         print('q_back.shape',q_back.shape)\n",
    "#         p_cand=10*stats.multivariate_normal.pdf(q_cand,m.reshape(-1,1),S)\n",
    "#         p_back=10*stats.multivariate_normal.pdf(x_intial,m.reshape(-1,1),S)\n",
    "        p_cand = prop_fn(q_cand)\n",
    "        p_back = prop_fn(x_intial)\n",
    "        if np.all((p_cand*q_back) >= u*(q_cand*p_back)):\n",
    "            sample.append(q_cand)\n",
    "            x_intial = q_cand\n",
    "            aceptence_rate += 1\n",
    "        else:\n",
    "            sample.append(x_intial)\n",
    "    return np.array(sample), aceptence_rate/nb_iter\n",
    "\n",
    "def prop_fn(x,m,S,k=10):  \n",
    "    return 10*stats.multivariate_normal.pdf(x,m,S)\n",
    "\n",
    "def posterior_sample(X, y, m, S, nb_iter):\n",
    "    # X: N x D matrix of training inputs\n",
    "    # y: N x 1 vector of training targets/observations\n",
    "    # m: D x 1 prior mean of parameters\n",
    "    # S: D x D prior covariance of parameters\n",
    "    # returns: nb_iter x D matrix of posterior samples\n",
    "    N,D = X.shape\n",
    "    \n",
    "    mu_post = map_estimate(X, y, m.reshape(-1,1), S)#np.zeros((D, 1)) ## <-- EDIT THESE LINES\n",
    "    S_post  = np.linalg.inv(Hessian(X,y,mu_post))#np.eye(D)\n",
    "    \n",
    "    fun = lambda x: prop_fn(x ,mu_post ,S_post)\n",
    "    \n",
    "    D = X.shape[1]\n",
    "    samples_MCMC, ac_rate = metropolis_algorithim(mu_post + 2,0.01*np.eye(D),nb_iter,fun) #np.zeros((nb_iter, D))\n",
    "#     print('acceptence rate: ', ac_rate*100)\n",
    "    \n",
    "    return samples_MCMC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "acceptence rate:  25.28\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([8.000e+00, 5.000e+00, 5.000e+00, 7.000e+00, 1.100e+01, 2.700e+01,\n",
       "        9.000e+00, 3.960e+02, 3.161e+03, 6.372e+03]),\n",
       " array([0.00540387, 0.00564921, 0.00589455, 0.00613989, 0.00638523,\n",
       "        0.00663057, 0.00687591, 0.00712125, 0.00736659, 0.00761193,\n",
       "        0.00785727]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEa5JREFUeJzt3X+sX3ddx/Hny5WBItKO3S1Nu9BNGmT+wViabYSEKDXd\nDwydCUtKjDRkpn84jSYaLfLHdEgCJgqSyExlxULUbU5xDSyMpkCMfzB2x8ZgG7OXMem1Y73SbQpE\nyODtH9/P1e+6e/v93u7e7+3u5/lIvjnnvM/nnO9553R99Zzv+X6XqkKS1J+fWO0DkCStDgNAkjpl\nAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1Kl1q30Ap3LuuefWli1bVvswJOlF5b777vvP\nqpoaNe6MDoAtW7YwPT292ochSS8qSf59nHHeApKkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAk\nqVMGgCR1ygCQpE6d0d8ElqTVtGXvp1ftvR9//1tX/D28ApCkThkAktQpA0CSOmUASFKnDABJ6pQB\nIEmdMgAkqVNjBUCS9UnuSPL1JI8keWOSc5IcSnKkTTe0sUny4SQzSR5McunQfna38UeS7F6ppiRJ\no417BfAXwGeq6ueA1wOPAHuBw1W1FTjclgGuBra21x7gZoAk5wA3ApcDlwE3zoeGJGnyRgZAkp8B\n3gzcAlBVP6yqp4GdwIE27ABwbZvfCXy8Br4IrE+yEbgSOFRVJ6rqKeAQcNWydiNJGts4VwAXAXPA\nx5Lcn+SjSV4OnF9VTwC06Xlt/Cbg6ND2s622WF2StArGCYB1wKXAzVX1BuB7/P/tnoVkgVqdov7c\njZM9SaaTTM/NzY1xeJKk0zFOAMwCs1V1T1u+g0EgPNlu7dCmx4fGXzC0/Wbg2Cnqz1FV+6pqW1Vt\nm5qaWkovkqQlGBkAVfVt4GiS17bSduBh4CAw/yTPbuDONn8QeGd7GugK4Jl2i+huYEeSDe3D3x2t\nJklaBeP+HPRvAX+b5GzgMeBdDMLj9iTXA98Crmtj7wKuAWaA77exVNWJJO8F7m3jbqqqE8vShSRp\nycYKgKp6ANi2wKrtC4wt4IZF9rMf2L+UA5QkrQy/CSxJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6\nZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMG\ngCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASerUWAGQ5PEkX03yQJLpVjsnyaEkR9p0Q6snyYeTzCR5\nMMmlQ/vZ3cYfSbJ7ZVqSJI1jKVcAv1hVl1TVtra8FzhcVVuBw20Z4Gpga3vtAW6GQWAANwKXA5cB\nN86HhiRp8l7ILaCdwIE2fwC4dqj+8Rr4IrA+yUbgSuBQVZ2oqqeAQ8BVL+D9JUkvwLgBUMBnk9yX\nZE+rnV9VTwC06Xmtvgk4OrTtbKstVn+OJHuSTCeZnpubG78TSdKSrBtz3Juq6liS84BDSb5+irFZ\noFanqD+3ULUP2Aewbdu2562XJC2Psa4AqupYmx4HPsngHv6T7dYObXq8DZ8FLhjafDNw7BR1SdIq\nGBkASV6e5BXz88AO4GvAQWD+SZ7dwJ1t/iDwzvY00BXAM+0W0d3AjiQb2oe/O1pNkrQKxrkFdD7w\nySTz4/+uqj6T5F7g9iTXA98Crmvj7wKuAWaA7wPvAqiqE0neC9zbxt1UVSeWrRNJ0pKMDICqegx4\n/QL17wDbF6gXcMMi+9oP7F/6YUqSlpvfBJakThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMG\ngCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBI\nUqcMAEnqlAEgSZ0yACSpU2MHQJKzktyf5FNt+cIk9yQ5kuS2JGe3+kvb8kxbv2VoH+9u9UeTXLnc\nzUiSxreUK4DfBh4ZWv4A8MGq2go8BVzf6tcDT1XVa4APtnEkuRjYBfw8cBXwkSRnvbDDlySdrrEC\nIMlm4K3AR9tygLcAd7QhB4Br2/zOtkxbv72N3wncWlU/qKpvAjPAZcvRhCRp6ca9AvgQ8PvAj9vy\nq4Cnq+rZtjwLbGrzm4CjAG39M238/9UX2EaSNGEjAyDJLwPHq+q+4fICQ2vEulNtM/x+e5JMJ5me\nm5sbdXiSpNM0zhXAm4C3JXkcuJXBrZ8PAeuTrGtjNgPH2vwscAFAW/9K4MRwfYFt/k9V7auqbVW1\nbWpqaskNSZLGMzIAqurdVbW5qrYw+BD3c1X1q8Dngbe3YbuBO9v8wbZMW/+5qqpW39WeEroQ2Ap8\nadk6kSQtybrRQxb1B8CtSf4EuB+4pdVvAT6RZIbBv/x3AVTVQ0luBx4GngVuqKofvYD3lyS9AEsK\ngKr6AvCFNv8YCzzFU1X/A1y3yPbvA9631IOUJC0/vwksSZ0yACSpUwaAJHXKAJCkThkAktQpA0CS\nOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlT\nBoAkdcoAkKROGQCS1CkDQJI6tW61D0CSRtmy99OrfQhrklcAktSpkQGQ5GVJvpTkK0keSvLHrX5h\nknuSHElyW5KzW/2lbXmmrd8ytK93t/qjSa5cqaYkSaONcwXwA+AtVfV64BLgqiRXAB8APlhVW4Gn\ngOvb+OuBp6rqNcAH2ziSXAzsAn4euAr4SJKzlrMZSdL4RgZADXy3Lb6kvQp4C3BHqx8Arm3zO9sy\nbf32JGn1W6vqB1X1TWAGuGxZupAkLdlYnwEkOSvJA8Bx4BDwDeDpqnq2DZkFNrX5TcBRgLb+GeBV\nw/UFthl+rz1JppNMz83NLb0jSdJYxgqAqvpRVV0CbGbwr/bXLTSsTbPIusXqJ7/XvqraVlXbpqam\nxjk8SdJpWNJTQFX1NPAF4ApgfZL5x0g3A8fa/CxwAUBb/0rgxHB9gW0kSRM2zlNAU0nWt/mfBH4J\neAT4PPD2Nmw3cGebP9iWaes/V1XV6rvaU0IXAluBLy1XI5KkpRnni2AbgQPtiZ2fAG6vqk8leRi4\nNcmfAPcDt7TxtwCfSDLD4F/+uwCq6qEktwMPA88CN1TVj5a3HUnSuEYGQFU9CLxhgfpjLPAUT1X9\nD3DdIvt6H/C+pR+mJGm5+U1gSeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnq\nlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4Z\nAJLUKQNAkjo1MgCSXJDk80keSfJQkt9u9XOSHEpypE03tHqSfDjJTJIHk1w6tK/dbfyRJLtXri1J\n0ijjXAE8C/xuVb0OuAK4IcnFwF7gcFVtBQ63ZYCrga3ttQe4GQaBAdwIXA5cBtw4HxqSpMkbGQBV\n9URVfbnN/zfwCLAJ2AkcaMMOANe2+Z3Ax2vgi8D6JBuBK4FDVXWiqp4CDgFXLWs3kqSxLekzgCRb\ngDcA9wDnV9UTMAgJ4Lw2bBNwdGiz2VZbrC5JWgVjB0CSnwb+EfidqvqvUw1doFanqJ/8PnuSTCeZ\nnpubG/fwJElLNFYAJHkJg7/8/7aq/qmVn2y3dmjT460+C1wwtPlm4Ngp6s9RVfuqaltVbZuamlpK\nL5KkJRjnKaAAtwCPVNWfD606CMw/ybMbuHOo/s72NNAVwDPtFtHdwI4kG9qHvztaTZK0CtaNMeZN\nwK8BX03yQKv9IfB+4PYk1wPfAq5r6+4CrgFmgO8D7wKoqhNJ3gvc28bdVFUnlqULSdKSjQyAqvpX\nFr5/D7B9gfEF3LDIvvYD+5dygJKkleE3gSWpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQB\nIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS\n1CkDQJI6ZQBIUqcMAEnqlAEgSZ0aGQBJ9ic5nuRrQ7VzkhxKcqRNN7R6knw4yUySB5NcOrTN7jb+\nSJLdK9OOJGlc41wB/A1w1Um1vcDhqtoKHG7LAFcDW9trD3AzDAIDuBG4HLgMuHE+NCRJq2NkAFTV\nvwAnTirvBA60+QPAtUP1j9fAF4H1STYCVwKHqupEVT0FHOL5oSJJmqDT/Qzg/Kp6AqBNz2v1TcDR\noXGzrbZY/XmS7EkynWR6bm7uNA9PkjTKcn8InAVqdYr684tV+6pqW1Vtm5qaWtaDkyT9v9MNgCfb\nrR3a9HirzwIXDI3bDBw7RV2StEpONwAOAvNP8uwG7hyqv7M9DXQF8Ey7RXQ3sCPJhvbh745WkySt\nknWjBiT5e+AXgHOTzDJ4muf9wO1Jrge+BVzXht8FXAPMAN8H3gVQVSeSvBe4t427qapO/mBZkjRB\nIwOgqt6xyKrtC4wt4IZF9rMf2L+ko5MkrRi/CSxJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBI\nUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOjXy56Alad6WvZ9e7UPQMvIKQJI6ZQBI\nUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktSpiQdAkquSPJpkJsneSb+/JGlgogGQ5CzgL4Gr\ngYuBdyS5eJLHIEkamPQVwGXATFU9VlU/BG4Fdk74GCRJTP63gDYBR4eWZ4HLJ3wMWkN6/G2ax9//\n1tU+BK0Rkw6ALFCr5wxI9gB72uJ3kzy64kc1cC7wnxN6rzNJr33Di7T3fOAF7+JF2fcyeFH1/QLP\n86vHGTTpAJgFLhha3gwcGx5QVfuAfZM8KIAk01W1bdLvu9p67Rv67d2+NW/SnwHcC2xNcmGSs4Fd\nwMEJH4MkiQlfAVTVs0l+E7gbOAvYX1UPTfIYJEkDE/8fwlTVXcBdk37fMUz8ttMZote+od/e7VsA\npKpGj5IkrTn+FIQkdWpNBcCon5lI8tIkt7X19yTZMrTu3a3+aJIrh+qPJ/lqkgeSTA/V/yjJf7T6\nA0muWen+FrNCfa9PckeSryd5JMkbW/2cJIeSHGnTDZPocSET7nvNnu8krx3q64Ek/5Xkd9q6NXu+\nR/R9xpzvFVVVa+LF4EPlbwAXAWcDXwEuPmnMbwB/1eZ3Abe1+Yvb+JcCF7b9nNXWPQ6cu8D7/RHw\ne2u47wPAr7f5s4H1bf5Pgb1tfi/wgU76XtPn+6T9fxt4dQ/n+xR9nxHne6Vfa+kKYJyfmdjJ4D9w\ngDuA7UnS6rdW1Q+q6pvATNvfi8Gy953kZ4A3A7cAVNUPq+rpBfZ1ALh2hfoaZdJ9nylW+s/5duAb\nVfXvC+xrTZ3vk7Y9ue8urKUAWOhnJjYtNqaqngWeAV41YtsCPpvkvgy+pTzsN5M8mGT/Kl4ar0Tf\nFwFzwMeS3J/ko0le3sacX1VPtH09AZy3vO2MbdJ9w9o938N2AX8/tLyWz/ewk/uGM+N8r6i1FAAj\nf2biFGNOte2bqupSBr9gekOSN7f6zcDPApcATwB/tuQjXh4r0fc64FLg5qp6A/A9Bpf/Z5JJ972W\nz/dgo8GXM98G/MNpH93KmXTfZ8r5XlFrKQBG/szE8Jgk64BXAidOtW1VzU+PA5+kXTpW1ZNV9aOq\n+jHw16zeLaOV6HsWmK2qe1r9DgZ/MQI8mWRj29dG4PiydbI0E+17jZ/veVcDX66qJ4dqa/l8z3te\n32fQ+V5RaykAxvmZiYPA7jb/duBzNfjE5yCwqz1FcCGwFfhSkpcneQVAuxWwA/haW944tN9fma+v\ngmXvu6q+DRxN8tq2zXbg4QX2tRu4cyWaGsNE+17L53tou3fw/Nsga/Z8D233vL7PoPO9slb7U+jl\nfAHXAP/G4FP+97TaTcDb2vzLGFzmzTD4A3DR0Lbvads9ClzdahcxeHrgK8BD8/ts6z4BfBV4kMEf\nsI1rpe9WvwSYbv39M7Ch1V8FHAaOtOk5nfS91s/3TwHfAV550nut9fO9WN9nzPleyZffBJakTq2l\nW0CSpCUwACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6tT/AqO0djreD1tHAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1660849240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## ADD CODE HERE\n",
    "# Plot a histogram of the MCMC posterior samples and the\n",
    "# analytical expression for the Laplace posterior. If\n",
    "# everything went well, the peaks should coincide and\n",
    "# their widths should be comparable.\n",
    "x_int=th_map#.reshape(-1,1)\n",
    "# var=200*S_post\n",
    "var=20*np.eye(D)\n",
    "\n",
    "nb_iter = 10000\n",
    "samples = posterior_sample(x, y, x_int, var, nb_iter)\n",
    "y1=stats.multivariate_normal.pdf(samples,np.zeros(D),var).reshape(-1,1)\n",
    "plt.hist(y1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et violà! Now you're an expert in logistic regression. (Wait, I think that's a big violin. I meant to say: et voilà!)\n",
    "\n",
    "Now we can visualize the posterior we can play around with the data and the inference parameters:\n",
    "\n",
    "* Play around with the data generation process. What happens as you increase/decrease $N$ and the jitter parameter?\n",
    "* What does the joint posterior look like? Make a visualization of the MCMC and Laplace approximations in the $(\\theta_0, \\theta_1)$ plane.\n",
    "* What happens if the model is misspecified? Take out the intercept term in the model (i.e., remove the column of ones in $\\boldsymbol X$), but set the `offset` in the data generation process to non-zero. What happens to the posterior and its Laplace approximation?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE31JREFUeJzt3X+wpmV93/H3Zw/7AwQ0K6I7u4vQuE5caeqPLWidJMwA\nGXBSNjOxEfpDyGh2ppamnXYy0uqQlvxj06lpWmnNNjqgsSKSJm7tMihEm05HkK3iDySEDWpYQRcR\nULuwsHu+/eM82OPhOXvOfT/3nh97vV8zZ85zP891ne91n3v3c65zP/d17lQVkqS2rFnuAUiSlp7h\nL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWrQSUN8kSQfAn4JOFhV5455PcDvAW8G\nDgFXVdUXj/U112V9beAFQwxPkprxQx7/XlW9ZKF2g4Q/cAPwfuDD87x+KbBt9HE+8J9Hn+e1gRdw\nfi4caHiS1Ibb65ZvLabdIKd9qurPgO8fo8lO4MM1407gRUk2DVFbktTdUp3z3ww8NGv7wOg5SdIy\nGOq0z0Iy5rnn/TnRJLuAXQAbOOV4j0mSmrVUM/8DwNZZ21uAh+c2qqrdVbWjqnasZf0SDU2S2rNU\n4b8HeFtmvAF4sqoeWaLakqQ5hrrU82PABcAZSQ4AvwWsBaiqDwB7mbnMcz8zl3r+2hB1JUn9DBL+\nVXXFAq8X8I+GqCVJmtxSveErSSeGNVPd+9R0v1rH8Ta7/nkHSWqQ4S9JDTL8JalBhr8kNcjwl6QG\nGf6S1CDDX5Ia5HX+ko6PjPt7jotwHK9tH8T00aWr1WdNwSKH58xfkhpk+EtSgwx/SWqQ4S9JDTL8\nJalBhr8kNcjwl6QGGf6S1CAXeUk6Pvou1uqzOGylLwzrq+9NYBZhkJl/kkuS3J9kf5Jrxrx+VpLP\nJvlSkq8kefMQdSVJ/Uwc/kmmgOuBS4HtwBVJts9p9h7g5qp6LXA58J8mrStJ6m+Imf95wP6qerCq\nngFuAnbOaVPA6aPHLwQeHqCuJKmnIc75bwYemrV9ADh/Tpt/BXw6yT8GXgBcNEBdSVJPQ8z8x707\nM/fdlyuAG6pqC/Bm4CNJnlc7ya4k+5Lse5bDAwxNkjTOEOF/ANg6a3sLzz+t83bgZoCq+jywAThj\n7heqqt1VtaOqdqxl/QBDkySNM0T43w1sS3JOknXMvKG7Z06bvwIuBEjyKmbC/9EBakuSepg4/Kvq\nCHA1cBtwHzNX9dyb5Lokl42a/XPg15N8GfgYcFXViXphriStfIMs8qqqvcDeOc9dO+vx14E3DVFL\n0jJw4dXyOI7fQ/+8gyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBK/tOXn0WlvTh\nYhTp2Pw/csJx5i9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoNW9nX+XXktsrT6+f94STjz\nl6QGGf6S1KBBwj/JJUnuT7I/yTXztPnVJF9Pcm+S/zpEXUlSPxOf808yBVwPXAwcAO5Osmd00/bn\n2mwD/gXwpqp6PMmZk9aVJPU3xMz/PGB/VT1YVc8ANwE757T5deD6qnocoKoODlBXktTTEOG/GXho\n1vaB0XOzvRJ4ZZL/neTOJJeM+0JJdiXZl2TfsxweYGiSpHGGuNRz3N9dnnut1knANuACYAvwv5Kc\nW1VP/ESnqt3AboDTs9HrvSTpOBli5n8A2Dprewvw8Jg2n6yqZ6vqG8D9zPwwkCQtgyFm/ncD25Kc\nA3wbuBz4u3Pa/AlwBXBDkjOYOQ304IJfOV1/Nk13bD8qc9Lazn3qyLO9armARdJKMPHMv6qOAFcD\ntwH3ATdX1b1Jrkty2ajZbcBjSb4OfBb4zap6bNLakqR+Uit0Jnp6Ntb5U7/YrVM585fUttvrlv9T\nVTsWaucKX0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktSglX0nr+mjnZpn7bpeZTLV/WdgHe35c3Pc\nH8NYSMfvgyQtxJm/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNWtnX+afbRfF1dOmuh1+z\nYX2vftNPPTXwSCSpO2f+ktQgw1+SGmT4S1KDBgn/JJckuT/J/iTXHKPdW5JUkgVvMSZJOn4mDv8k\nU8D1wKXAduCKJNvHtDsN+A3grklrSpImM8TM/zxgf1U9WFXPADcBO8e0+23gd4CnB6gpSZrAEOG/\nGXho1vaB0XM/luS1wNaq+tQA9SRJExriOv9xF+PXj19M1gC/C1y14BdKdgG7ADZwygBDkySNM0T4\nHwC2ztreAjw8a/s04Fzgc5lZtPUyYE+Sy6pq3+wvVFW7gd0Ap2djkY6/mPS86Ul6LNjKKSf3qjW1\nvnuto48/3quWJM1niNM+dwPbkpyTZB1wObDnuRer6smqOqOqzq6qs4E7gecFvyRp6Uwc/lV1BLga\nuA24D7i5qu5Ncl2Syyb9+pKk4Q3yt32qai+wd85z187T9oIhakqS+nOFryQ1yPCXpAYZ/pLUIMNf\nkhpk+EtSg1b2nbxqulPzNaed1qtMprr/DKzNZ/ar9Z3HOvdZc0r31c7TTx/u3Ke3novrJC0fZ/6S\n1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBq3wRV61cJtZRncK615my6bOfQ6d1W9B\n2cnrun/Lp76zrnOfOvho5z4AdXgJF4dJWjbO/CWpQYa/JDXI8JekBg0S/kkuSXJ/kv1Jrhnz+j9L\n8vUkX0lyR5KXD1FXktTPxOGfZAq4HrgU2A5ckWT7nGZfAnZU1c8CtwC/M2ldSVJ/Q8z8zwP2V9WD\nVfUMcBOwc3aDqvpsVR0abd4JbBmgriSppyHCfzPw0KztA6Pn5vN24NYB6kqSehriOv9xF9ePvUA/\nyd8HdgC/MM/ru4BdABvocQOTQ4cWbjSu7lT39QGPvHGqV63pdS/o3OcVH+12UxuANU/+oHMfgKNe\n5y81YYjwPwBsnbW9BXh4bqMkFwHvBn6hqsYmTFXtBnYDnJ6N3VZ4SZIWbYjTPncD25Kck2QdcDmw\nZ3aDJK8Ffh+4rKoODlBTkjSBicO/qo4AVwO3AfcBN1fVvUmuS3LZqNm/BU4FPpHkniR75vlykqQl\nMMjf9qmqvcDeOc9dO+vxRUPUkSQNwxW+ktQgw1+SGmT4S1KDDH9JatDKvplLRzXdc2nA2u4Ltp49\n89lepd74qr/s3Ofzp72ic59XfrDf385bc983O/eZ/uEPe9WStHyc+UtSgwx/SWqQ4S9JDTL8JalB\nhr8kNcjwl6QGGf6S1CDDX5IadEIt8mL6aL9uPRZ5bXxpvztlnbb26c59Nmzs3uepTSd37gNw6p93\nv6sZa/rd1azv8eosPfapr1q6exBl/frOfWop79TW9/u+hN/Dljnzl6QGGf6S1CDDX5IaZPhLUoMM\nf0lq0CDhn+SSJPcn2Z/kmjGvr0/y8dHrdyU5e4i6kqR+Jg7/JFPA9cClwHbgiiTb5zR7O/B4Vb0C\n+F3g30xaV5LU3xAz//OA/VX1YFU9A9wE7JzTZidw4+jxLcCFyVJefC1Jmm2IRV6bgYdmbR8Azp+v\nTVUdSfIk8GLge7MbJdkF7ALYwCkDDG1xTnryqe59PvHiXrVe9a4vdO7z+2/6fOc+f/2ud3buA3DK\nt8/q3CdffaBXrepxM7Ss7f5PNuvWdS8ETP/fQ9071RItXAPqmWeWrFYvLtZa0YaY+Y+bwc896otp\nQ1XtrqodVbVjLd1XL0qSFmeI8D8AbJ21vQV4eL42SU4CXgh8f4DakqQehgj/u4FtSc5Jsg64HNgz\np80e4MrR47cAf1rl74SStFwmPuc/Ood/NXAbMAV8qKruTXIdsK+q9gAfBD6SZD8zM/7LJ60rSepv\nkL/qWVV7gb1znrt21uOngb8zRC1J0uRc4StJDTL8JalBJ9bNXPrq8d7zo+dN9yp163fO7dzn975w\nYec+P/crX+vcB+DgHVsXbjRH7+vNe3zf63D36+iX9AYmS2hNj5u5TD/d/cZAOjE585ekBhn+ktQg\nw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yEVeQA53v6vIOX98pFetbz3W/WYpbO6+iOpvvfAv\nu9cB3vfWV3fu89PTP9Or1vSX7+vVr6s1p/S7MdD0oR43c1lCLtjSJJz5S1KDDH9JapDhL0kNMvwl\nqUGGvyQ1yPCXpAZNFP5JNib5TJIHRp9/akyb1yT5fJJ7k3wlyVsnqSlJmtykM/9rgDuqahtwx2h7\nrkPA26rq1cAlwL9P8qIJ60qSJjDpIq+dwAWjxzcCnwPeNbtBVf3FrMcPJzkIvAR4YsLagzn67Uc6\n99nQ4y5UAGc9/bxfjhb0jb99cuc+H33PL3XuA3DS9nTus+aJH/WqdfgXd3Tus/bT+zr3WemLtaTl\nMOnM/6VV9QjA6POZx2qc5DxgHdBv+akkaRALzvyT3A68bMxL7+5SKMkm4CPAlVU19ga4SXYBuwA2\n0G9JviRpYQuGf1VdNN9rSb6bZFNVPTIK94PztDsd+B/Ae6rqzmPU2g3sBjg9G/udV5EkLWjS0z57\ngCtHj68EPjm3QZJ1wB8DH66qT0xYT5I0gEnD/73AxUkeAC4ebZNkR5I/GLX5VeDngauS3DP6eM2E\ndSVJE5joap+qegy4cMzz+4B3jB7/IfCHk9SRJA3LFb6S1CBv5gJMbd7Uuc+Rb3yrV62jf+2Mzn1+\n+uM/6Nxn6uDjnfsAHHxr9/Ed+asDvWqt/+6jnfuMvUxMUmfO/CWpQYa/JDXI8JekBhn+ktQgw1+S\nGmT4S1KDDH9JapDhL0kNcpEXcOTBby5ZrQ0Pfq9zn3q8+31vjjzxZOc+AFNfPbtzn7xme69a01+6\nt1c/SZNz5i9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMmCv8kG5N8JskDo88/dYy2pyf5\ndpL3T1JTkjS5SRd5XQPcUVXvTXLNaPtd87T9beB/Tlhv1et7B7Cl8vL33dO5z/ShQ71qHb70b3bu\ns/7Wu3vVkvSTJj3tsxO4cfT4RuCXxzVK8nrgpcCnJ6wnSRrApOH/0qp6BGD0+cy5DZKsAf4d8JsT\n1pIkDWTB0z5JbgdeNualdy+yxjuBvVX1UJKFau0CdgFs4JRFfnlJUlcLhn9VXTTfa0m+m2RTVT2S\nZBNwcEyzNwI/l+SdwKnAuiQ/qqprxtTaDewGOD0ba7E7IUnqZtI3fPcAVwLvHX3+5NwGVfX3nnuc\n5Cpgx7jglyQtnUnP+b8XuDjJA8DFo22S7EjyB5MOTpJ0fEw086+qx4ALxzy/D3jHmOdvAG6YpKYk\naXKu8JWkBnknr6W2Zqp7n+mjw49jvlI9Fmyt2bChV61T9n2zc5+l+05IJzZn/pLUIMNfkhpk+EtS\ngwx/SWqQ4S9JDTL8JalBhr8kNcjr/JfaEl6zv1Smn366V796/c907pNHH+1VS9JPcuYvSQ0y/CWp\nQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JalCqarnHMFaSR4FvLfc4OjgD+N5yD2JAJ9L+uC8r\nk/tyfLy8ql6yUKMVG/6rTZJ9VbVjuccxlBNpf9yXlcl9WV6e9pGkBhn+ktQgw384u5d7AAM7kfbH\nfVmZ3Jdl5Dl/SWqQM39JapDh31GSDyU5mORr87yeJP8hyf4kX0nyuqUe42ItYl8uSPJkkntGH9cu\n9RgXK8nWJJ9Ncl+Se5P8kzFtVsWxWeS+rIpjk2RDki8k+fJoX/71mDbrk3x8dFzuSnL20o90YYvc\nl6uSPDrruLxjOca6KFXlR4cP4OeB1wFfm+f1NwO3AgHeANy13GOeYF8uAD613ONc5L5sAl43enwa\n8BfA9tV4bBa5L6vi2Iy+16eOHq8F7gLeMKfNO4EPjB5fDnx8ucc9wb5cBbx/uce6mA9n/h1V1Z8B\n3z9Gk53Ah2vGncCLkmxamtF1s4h9WTWq6pGq+uLo8Q+B+4DNc5qtimOzyH1ZFUbf6x+NNteOPua+\n0bgTuHH0+BbgwiRZoiEu2iL3ZdUw/Ie3GXho1vYBVul/3JE3jn7NvTXJq5d7MIsxOm3wWmZmZrOt\numNzjH2BVXJskkwluQc4CHymquY9LlV1BHgSePHSjnJxFrEvAL8yOq14S5KtSzzERTP8hzduxrJa\nZwdfZGap+N8A/iPwJ8s8ngUlORX4I+CfVtUP5r48psuKPTYL7MuqOTZVdbSqXgNsAc5Lcu6cJqvm\nuCxiX/47cHZV/SxwO///N5oVx/Af3gFg9k/7LcDDyzSWiVTVD577Nbeq9gJrk5yxzMOaV5K1zITl\nR6vqv41psmqOzUL7stqODUBVPQF8Drhkzks/Pi5JTgJeyAo/HTnfvlTVY1V1eLT5X4DXL/HQFs3w\nH94e4G2jK0veADxZVY8s96D6SPKy5869JjmPmX8vjy3vqMYbjfODwH1V9b55mq2KY7OYfVktxybJ\nS5K8aPT4ZOAi4M/nNNsDXDl6/BbgT2v07ulKsph9mfMe0mXMvF+zIp203ANYbZJ8jJkrLc5IcgD4\nLWbe+KGqPgDsZeaqkv3AIeDXlmekC1vEvrwF+IdJjgBPAZevxP+UI28C/gHw1dE5WYB/CZwFq+7Y\nLGZfVsux2QTcmGSKmR9QN1fVp5JcB+yrqj3M/KD7SJL9zMz4L1++4R7TYvblN5JcBhxhZl+uWrbR\nLsAVvpLUIE/7SFKDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhr0/wC6TkUuwz011AAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f166075a668>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist2d(samples[:,1],samples[:,0], bins=25);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}