Test_LikelihoodBernoulli.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:072ba15e3fefdd444e83774b0946732bad16b1c3f57064a4c6ce2c4d3e737634"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy as sp\n",
      "import scipy.special\n",
      "import scipy.linalg\n",
      "import hips.distributions.polya_gamma as pg\n",
      "%matplotlib inline\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "if 0:\n",
      "    # Generate test data\n",
      "    K = 10\n",
      "    trials = 1\n",
      "    dN = np.random.rand(K, trials) > .9\n",
      "    #dN = np.bitwise_xor.accumulate(dN, axis=0)\n",
      "    dN = np.sum(dN, axis=1)/5\n",
      "    sigma2e = 0.1\n",
      "else:\n",
      "    # Load thaldata\n",
      "    import pandas as pd\n",
      "    dN = pd.read_csv('thaldata.csv', header=None).values[0]#[500:540]#800]\n",
      "    trials = 50\n",
      "    sigma2e = 0.12"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Likelihood of a state-space model with Bernoulli observations\n",
      "\n",
      "The log-likelihood is\n",
      "\\begin{align*}\n",
      "\\log p(N_{0,K}, x; \\theta) &= \\log p(N_{0,K}|x) + \\log p(x; \\theta) \\\\\n",
      "&= \\sum\\limits_{k=1}^{K} \\log p(\\Delta N_{k}|x_k) + \\sum\\limits_{k=1}^{K} \\log p(x_k|x_{k-1}; \\sigma_{\\epsilon}^2) \\\\\n",
      "&= \\sum\\limits_{k=1}^{K} \\left[\\Delta N_{k} \\log\\left(\\frac{e^{x_k}}{1+e^{x_k}}\\right) + (1-\\Delta N_{k}) \\log\\left(\\frac{1}{1+e^{x_k}}\\right) \\right] + \\pi(x; \\theta) \\\\\n",
      "&= \\sum\\limits_{k=1}^{K} \\left[ \\Delta N_{k} \\, x_k - \\log(1+e^{x_k}) \\right] + \\pi(x; \\theta)\n",
      "\\end{align*}\n",
      "where\n",
      "$$\n",
      "\\pi(x; \\theta) = - \\frac{K}{2} \\log(2\\pi\\sigma_{\\epsilon}^2) - \\frac{1}{2\\sigma_{\\epsilon}^2} \\sum\\limits_{k=1}^{K} (x_k - x_{k-1})^2\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Marginal likelihood using Polya-Gamma $PG(1,0)$ latent variables\n",
      "\n",
      "Using the Polya-Gamma latent variables $w_k \\sim PG(1,0)$, we can rewrite $p(\\Delta N_{k}|x_k)$ as (see http://dx.doi.org/10.1080/01621459.2013.829001)\n",
      "\\begin{align}\n",
      "p(\\Delta N_{k}|x_k) &= \\frac{e^{x_k\\,\\Delta N_k}}{1+e^{x_k}} = \\frac12 e^{w_k \\, y_k \\, x_k} \\int\\limits_{0}^{\\infty} e^{-w_k \\, x_k^2 / 2} \\, p(w_k) \\, dw_k \\\\\n",
      "&= \\int\\limits_{0}^{\\infty} \\frac{e^{w_k \\, y_k^2 / 2}}{\\sqrt{2 \\pi^{-1} w_k}} \\left[ \\frac{1}{\\sqrt{2\\pi w_k^{-1}}} e^{-\\frac{1}{2 \\, w_k^{-1}} \\left(y_k - x_k\\right)^2} \\right] p(w_k) \\, dw_k\n",
      "\\end{align}\n",
      "where $y_k = \\frac{\\Delta N_k - 1/2}{w_k}$.\n",
      "\n",
      "So, the data likelihood can be written as\n",
      "\\begin{align}\n",
      "p(N_{0,K}, x; \\theta) &= \\int\\limits_{w \\in \\mathbb{R}_{+}^K} c(w,y) \\, p(y,x; w^{-1}, \\sigma_{\\epsilon}^2) \\,  p(w) \\, dw\n",
      "\\end{align}\n",
      "where\n",
      "\\begin{align}\n",
      "c(w,y) &= \\frac{e^{\\sum_{k=1}^K w_k \\, y_k^2 / 2}}{\\prod_{k=1}^K\\sqrt{2 \\pi^{-1} w_k}}, \\\\\n",
      "p(y,x; w^{-1}, \\sigma_{\\epsilon}^2) &= \\prod_{k=1}^K \\left[\n",
      "\\frac{1}{\\sqrt{2\\pi w_k^{-1}}} e^{-\\frac{1}{2 \\, w_k^{-1}} \\left(y_k - x_k\\right)^2} \\,\n",
      "\\frac{1}{\\sqrt{2\\pi \\sigma_{\\epsilon}^2}} e^{-\\frac{1}{2 \\, \\sigma_{\\epsilon}^2} \\left(x_k - x_{k-1}\\right)^2}\n",
      "\\right]\n",
      "\\end{align}\n",
      "The term $p(y,x; w^{-1}, \\sigma_{\\epsilon}^2)$ is the data likelihood of a Gaussian state-space model with observations $y_k$, measurement noise with variance $w_k^{-1}$ and process noise $\\sigma_{\\epsilon}^2$.\n",
      "\n",
      "We can now compute the marginal likelihood as\n",
      "\\begin{align}\n",
      "p(N_{0,K}; \\theta) &= \\int\\limits_{x \\in \\mathbb{R}^K} \\; \\int\\limits_{w \\in \\mathbb{R}_{+}^K} c(w,y) \\, p(y,x; w^{-1}, \\sigma_{\\epsilon}^2) \\, p(w) \\, dw \\, dx \\\\\n",
      "&= \\int\\limits_{w \\in \\mathbb{R}_{+}^K} c(w,y) \\left[ \\int\\limits_{x \\in \\mathbb{R}^K} p(y,x; w^{-1}, \\sigma_{\\epsilon}^2) \\, dx \\right] p(w) \\, dw \\\\\n",
      "&= \\int\\limits_{w \\in \\mathbb{R}_{+}^K} c(w,y) \\, p(y; w^{-1}, \\sigma_{\\epsilon}^2) \\, dw\n",
      "\\end{align}\n",
      "where\n",
      "\\begin{align}\n",
      "J &= \\nabla_x^2 [-\\log p(y,x; w^{-1}, \\sigma_{\\epsilon}^2)]\\\\\n",
      "x_{\\mathrm{MAP}} &= 0_K - J^{-1} \\nabla_x [-\\log p(y,x; w^{-1}, \\sigma_{\\epsilon}^2)] \\big|_{x = 0_K}\\\\\n",
      "p(y; w^{-1}, \\sigma_{\\epsilon}^2) &= (2\\pi)^{K/2} \\, p(y, x_{\\mathrm{MAP}}; w^{-1}, \\sigma_{\\epsilon}^2) \\, |J|^{-1/2}\n",
      "\\end{align}\n",
      "(see http://dx.doi.org/10.1007/s10827-009-0150-x, then do some algebra).\n",
      "\n",
      "Finally, we can estimate the marginal likelihood for the Bernoulli observations by drawing $M$ vectors $w$ whose $K$ elements are sampled from $PG(1,0)$, then\n",
      "$$\n",
      "p(N_{0,K}; \\theta) \\simeq \\frac{1}{M} \\sum_{m=1}^{M} c(w^{(m)},y^{(m)}) \\, p(y^{(m)}; (w^{(m)})^{-1}, \\sigma_{\\epsilon}^2).\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Do_Kalman_Likelihood(y, sigma2obs, sigma2e):\n",
      "    \"\"\"MAP solution, inverse covariance matrix, and marginal loglikelihood of state-space model\n",
      "\n",
      "    :param y: Observations (K,)\n",
      "    :param sigma2obs: Variance of observation noise (can be scalar or vector)\n",
      "    :param sigma2e: Variance of process noise\n",
      "    :return: x_map, L, marginal_loglikelihood, joint_loglikelihood\n",
      "    \"\"\"\n",
      "    # Build diagonals of information matrix\n",
      "    sigma2obs *= np.ones(len(y))\n",
      "    D = 1. / sigma2obs + 2. / sigma2e\n",
      "    D[-1] = 1. / sigma2obs[-1] + 1. / sigma2e\n",
      "    B = -np.ones(len(D)) / sigma2e\n",
      "    B[-1] = 0.\n",
      "    \n",
      "    # Solve, assuming x_init=0 for simplicity\n",
      "    L = sp.linalg.cholesky_banded((D, B), lower=True)\n",
      "    x_map = sp.linalg.cho_solve_banded([L, True], y / sigma2obs)\n",
      "\n",
      "    # Compute joint and marginal probabilities\n",
      "    joint_loglikelihood = -.5 * ((np.sum(np.diff(x_map)**2) + x_map[0]**2) / sigma2e +\n",
      "                                 np.sum((y - x_map)**2 / sigma2obs) +\n",
      "                                 (len(y) * np.log(2*np.pi*sigma2e * 2*np.pi) + np.sum(np.log(sigma2obs))))\n",
      "    marginal_loglikelihood = len(y)/2. * np.log(2*np.pi) + joint_loglikelihood - np.sum(np.log(L[0]))\n",
      "    return x_map, L, marginal_loglikelihood, joint_loglikelihood"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Marginal_Likelihood_PG(dN, sigma2e, n_samples, ret_accum=False, trials=1.):\n",
      "    K = len(dN)\n",
      "    dN = dN.astype(float)\n",
      "    marginal_loglikelihood_sample = np.zeros(n_samples)\n",
      "    b, c = trials * np.ones(K), np.zeros(K)\n",
      "    log_c0 = K * (.5 * np.log(2*np.pi) - trials * np.log(2)) + np.sum(np.log(sp.special.binom(trials, dN)))\n",
      "    for i in range(n_samples):\n",
      "        w = pg.polya_gamma(a=b, c=c) # For b=1 c=0, same as J* distribution (http://dx.doi.org/10.1016/j.spl.2009.07.028)\n",
      "        y = (dN - trials/2.) / w\n",
      "        marginal_loglikelihood_gaus = Do_Kalman_Likelihood(y, 1. / w, sigma2e)[2]\n",
      "        log_c = log_c0 + .5 * (np.sum(w * y**2) - np.sum(np.log(w)))\n",
      "        marginal_loglikelihood_sample[i] = log_c + marginal_loglikelihood_gaus\n",
      "\n",
      "    offset = max(marginal_loglikelihood_sample)\n",
      "    if ret_accum:\n",
      "        accum_marginal_loglikelihood = offset + np.log(np.cumsum(np.exp(marginal_loglikelihood_sample - offset))\n",
      "                                                       / np.arange(1, n_samples+1))\n",
      "        return accum_marginal_loglikelihood\n",
      "    marginal_loglikelihood = offset + np.log(np.mean(np.exp(marginal_loglikelihood_sample - offset)))\n",
      "    return marginal_loglikelihood"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Simple test of Do_Kalman_Likelihood only (K=1, Gaussian)\n",
      "import scipy.stats as stats\n",
      "y = 3.\n",
      "sigma2obs = 10.\n",
      "x_map, L, marginal_loglikelihood, joint_loglikelihood = Do_Kalman_Likelihood(y * np.ones(1), sigma2obs, sigma2e)\n",
      "j = stats.norm.pdf(x_map, 0, np.sqrt(sigma2e)) * stats.norm.pdf(x_map, y, np.sqrt(sigma2obs))\n",
      "m = np.sqrt(2*np.pi) * j / L[0,0]\n",
      "assert(abs(np.log(j) - joint_loglikelihood) < 1e-9)\n",
      "assert(abs(np.log(m) - marginal_loglikelihood) < 1e-9)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Simple test of Marginal_Likelihood_PG (K=1)\n",
      "n_samples = 10000\n",
      "#state = np.random.get_state()\n",
      "#np.random.set_state(state)\n",
      "accum_marginal_loglikelihood = Marginal_Likelihood_PG(np.ones(1), sigma2e, n_samples, ret_accum=True)\n",
      "# normpdf(x,m,s2) := 1/sqrt(2*%pi*s2) * exp(-.5*(x-m)^2/s2);\n",
      "# bernoulli(y, x) := exp(x*y) / (1 + exp(x));\n",
      "# quad_qag(normpdf(x, 0, .5) * bernoulli(1, x), x, -25, 25, 3);\n",
      "#   [0.5000000000000001, 3.188337035858471E-10, 217, 0]\n",
      "# We can actually prove it's 0.5 because = integral( (.5+odd(x)) * gauspdf(x;0,s2e) ) = .5 + 0\n",
      "print(np.log(.5), accum_marginal_loglikelihood[-1])\n",
      "plt.plot(accum_marginal_loglikelihood)\n",
      "plt.plot([0,n_samples-1], [np.log(.5),np.log(.5)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(-0.69314718055994529, -0.69314244856156493)\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "[<matplotlib.lines.Line2D at 0x7f1e3689b810>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEACAYAAAB78OvLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH9VJREFUeJzt3XuUFOWd//H3Z5hIFOSqgFEQFTFEEwkocLJGRxCDJ66Q\naFyTddVo1N3EbE4uRvwdz5E97u6R7IkmMdH8NBEx+WmiMQkYQ0AcJu5uQkQUUUG8RUSNg3gBFQFh\nvr8/nmqnGXu4dQ3ddH9e5/Tp6qqnqp6qqalP11OXVkRgZmaWh4ZKV8DMzGqHQ8XMzHLjUDEzs9w4\nVMzMLDcOFTMzy41DxczMclN2qEjqK2mepBWS5krq3Um5wdnwZZIekzQk6z9B0mJJD0u6X9KhHcY7\nXVKbpFHl1tXMzLpWHkcqU4H5EXEE0Axc3km5W4HpEfERYAywOut/PfD5iPg4cDtwRWEEST2BfwUW\n5lBPMzPrYnmEymRgZtY9E5jSsYCkEUC3iGgGiIj1EbEhG9wGFI5uegMvFY16FXA1sDGHepqZWRfL\nI1QGREQrQES8DAwoUWY4sFbSXVlT13RJyoZdCMyR9DxwNilEkPRx4KCImJNDHc3MbDfYoVCRdK+k\npUWvR7P300oUL/Xcl0bgOOAbwLHAYcB52bCvA5MiYghwC3BtFjjXAN8srsaO1NXMzCqncUcKRcTE\nzoZJapU0MCJaJQ2i/VxJsReAJRGxMhvnt8BYSXcDR0fEg1m5XwJzgJ7AUUBLFjCDgFmSTouIhzrM\n3w8vMzPbBRGR+5f1PJq/ZtN+1HEuMKtEmUVAH0n9s8/jgceB14FekoZl/U8GlkfEmxGxf0QcGhGH\nkE7U/33HQCmICL8iuPLKKyteh2p5eV14XXhdbPvVVXboSGU7pgN3SDofWAmcCSBpNHBxRFwUEW2S\nvgU0Z6dSFgM/iYgtki4Efi1pCylkzi8xj8DNX2ZmVa/sUImI14CTSvRfDFxU9Pk+4OgS5WZR+uim\nuMz4cutpZmZdz3fU15CmpqZKV6FqeF2087po53XR9dSVbWu7g6TY05fBzGx3k0RU6Yl6MzMzwKFi\nZmY5cqiYmVluHCpmZpYbh4qZmeXGoWJmZrlxqJiZWW4cKmZmlhuHipmZ5cahYmZmuXGomJlZbhwq\nZmaWG4eKmZnlxqFiZma5caiYmVluHCpmZpYbh4qZmeXGoWJmZrlxqJiZWW4cKmZmlhuHipmZ5cah\nYmZmuXGomJlZbhwqZmaWG4eKmZnlxqFiZma5KStUJPWVNE/SCklzJfXupNzgbPgySY9JGpL1nyBp\nsaSHJd0v6dCicc6U9LikRyX9vJx6mpnZ7qGI2PWRpenAqxHxHUmXAX0jYmqJcguAqyKiWdI+QFtE\nbJC0Avj7iHhS0r8Ax0bE+ZIOB34BnBgR6yTtFxFrOqlDlLMMZmb1SBIRobynW27z12RgZtY9E5jS\nsYCkEUC3iGgGiIj1EbEhG9wGFI5uegMvZd1fAn4UEeuycUoGipmZVZfGMscfEBGtABHxsqQBJcoM\nB9ZKugsYCswHpmaHFxcCcyStB9YB44rGQdL/kILv3yJibpl1NTOzLrbdUJF0LzCwuBcQwBUlipdq\nh2oEjgNGAquAO4DzgBnA14FJEfGgpG8B15KCphEYBhwPDAHul3RU4cjFzMyq03ZDJSImdjZMUquk\ngRHRKmkQsLpEsReAJRGxMhvnt8BYSXcDR0fEg1m5XwJzisZZGBFtwHOSngQOBxaXqse0adPe625q\naqKpqWl7i2VmVldaWlpoaWnp8vnkcaL+tYiY3tmJekkNpDA4KSJelXQz8ABwE/A34BMR8bSkC0hH\nLZ+T9Cng8xFxnqT9svFHRsTrJergE/VmZjupq07Ulxsq/UjNWYOBlcCZEfGGpNHAxRFxUVZuAnBN\nNtpi4KKI2CxpMnAVsAV4HTg/Ip7LxvkuMAnYDPx7RNzZSR0cKmZmO6kqQ6UaOFTMzHZetV5SbGZm\n9h6HipmZ5cahYmZmuamJUPEpFTOz6uBQMTOz3DhUzMwsNw4VMzPLjUPFzMxy41AxM7Pc1ESomJlZ\ndaiJUPGRiplZdXComJlZbhwqZmaWm5oIFTMzqw4OFTMzy41DxczMcuNQMTOz3NREqPhEvZlZdaiJ\nUDEzs+rgUDEzs9w4VMzMLDcOFTMzy01NhIpP1JuZVYeaCBUzM6sODhUzM8uNQ8XMzHLjUDEzs9w4\nVMzMLDc1ESq++svMrDqUHSqS+kqaJ2mFpLmSendSbnA2fJmkxyQNyfpPkLRY0sOS7pd0aFH5ZkkP\nSVoi6ZRy62pmZl0rjyOVqcD8iDgCaAYu76TcrcD0iPgIMAZYnfW/Hvh8RHwcuB24Iut/BfDLiBgF\nfD4rZ2ZmVSyPUJkMzMy6ZwJTOhaQNALoFhHNABGxPiI2ZIPbgMLRTW/gpaw7gF5Zdx/gxRzqamZm\nXUhR5gkJSa9FRL/OPmf9JgNfAjYBQ4H5wNSICEnHAb8F1gPrgHER8ZakQcA8oC+wD3BSRDxcYv7x\n5ptBz55lLYaZWV2RREQo7+k27uDM7wUGFvciHUlcUaJ4qZRqBI4DRgKrgDuA84AZwNeBSRHxoKRv\nAtcCF5KavGZExLWSxgE/B44sVb//+I9pdO+eupuammhqatqRxTIzqxstLS20tLR0+XzyOFJZDjRF\nRGt2dLEgIkZ0KDMWuDoiTsw+nw2MBf4NWBgRw7L+g4E5EXGUpMeAT0XEi9mwZ4CxEbGmw7Rj3bpg\n333LWgwzs7rSVUcqeZxTmU066gA4F5hVoswioI+k/tnn8cDjwOtAL0nDsv4nA8uz7pXASfDeOZnu\nHQPFzMyqSx5HKv1IzVmDSUFwZkS8IWk0cHFEXJSVmwBck422GLgoIjZn51uuAraQQub8iHguC5Kb\ngJ6kk/mXRsR9JebvIxUzs53UVUcqZYdKpTlUzMx2XjU3f1XcHp6LZmY1oyZCxczMqkNNhMqqVZWu\ngZmZQY2Eypw5la6BmZlBjYRKQ00shZnZnq8mdscOFTOz6lATu2PlflGcmZntipoIlW7dKl0DMzOD\nGgkVN3+ZmVWHmtgdO1TMzKpDTeyOHSpmZtWhJnbHDhUzs+pQE7tjh4qZWXWoid2xQ8XMrDrUxO7Y\nlxSbmVWHmggVH6mYmVWHmtgd77NPpWtgZmZQI6FywAGVroGZmUGNhIp/+dHMrDrURKiYmVl1cKiY\nmVluaiJU3PxlZlYdaiJUzMysOjhUzMwsNzURKm7+MjOrDjURKmZmVh0cKmZmlhuHipmZ5aasUJHU\nV9I8SSskzZXUu5Nyg7PhyyQ9JmlI1n+8pMWSlkqaIamhaJwfSHpK0hJJI7dVD59TMTOrDuUeqUwF\n5kfEEUAzcHkn5W4FpkfER4AxwGpJAm4BzoyIjwErgXMBJJ0CHBYRhwMXAz8us55mZrYblBsqk4GZ\nWfdMYErHApJGAN0iohkgItZHxAagP7AxIp7Jis4HTi+a7q1Z+b8AvSUNLLOuZmbWxcoNlQER0QoQ\nES8DA0qUGQ6slXRX1tQ1XZIiYg3QKGlUVu4MYHDWfSCwqmgaL2b9SnLzl5lZdWjcXgFJ9wLFRwkC\nAriiRPFSu/dG4DhgJCko7gDOA2YAZwHfk7QXMA/YshN1f8/NN0/jvvtSd1NTE01NTbsyGTOzmtXS\n0kJLS0uXz0dRxtd8ScuBpoholTQIWBARIzqUGQtcHREnZp/PBsZGxFc7lJsIXBARZ0n6cTatX2bD\nngBOKBwVdRgv7r8/+OQnd3kxzMzqjiQiQnlPt9zmr9mkow5IJ9lnlSizCOgjqX/2eTywDEDS/tl7\nd+Ay2k/IzwbOyYaNA94oFSgFbv4yM6sO5YbKdGCipBXABOBqAEmjJd0IEBFtwLeAZkmPZOPdlL1f\nKmkZsASYFREt2Ti/B/4q6Wng/wJfLrOeZma2G5TV/FUNJMUf/xgcf3yla2Jmtueo1uavqrCH56KZ\nWc2oiVAxM7Pq4FAxM7Pc1ESouPnLzKw61ESomJlZdXComJlZbmrikuKPfjRYurTSNTEz23N01SXF\nNREqED6vYma2E3yfipmZVT2HipmZ5cahYmZmuXGomJlZbhwqZmaWG4eKmZnlxqFiZma5caiYmVlu\nHCpmZpYbh4qZmeXGoWJmZrlxqJiZWW4cKmZmlhuHipmZ5cahYmZmuXGomJlZbmoiVM4+u9I1MDMz\nqJFQ8a8+mplVB4eKmZnlxqFiZma5KStUJPWVNE/SCklzJfXupNzgbPgySY9JGpL1Hy9psaSlkmZI\nasj6f0HSI9nrfyR9dFv1cKiYmVUHRRl7ZEnTgVcj4juSLgP6RsTUEuUWAFdFRLOkfYA2YCOwEjgx\nIp6RNA1YGREzJI0DlkfEWkmTgGkRMa6TOgSEg8XMbCdIIiKU93TLbf6aDMzMumcCUzoWkDQC6BYR\nzQARsT4iNgD9gY0R8UxWdD5wRlZmYUSszfovBA4ss55mZrYblBsqAyKiFSAiXgYGlCgzHFgr6a6s\nqWu6JEXEGqBR0qis3BnAQSXG/xIwp8x6mpnZbtC4vQKS7gUGFvcCAriiRPFSjVCNwHHASGAVcAdw\nHjADOAv4nqS9gHnAlg7zPhH4Yjb+Nkxj2rTU1dTURFNT07aLm5nVmZaWFlpaWrp8PuWeU1kONEVE\nq6RBwIKIGNGhzFjg6og4Mft8NjA2Ir7aodxE4IKIOCv7/DHgLmBSURNZqTr4nIqZ2U6q1nMqs0lH\nHQDnArNKlFkE9JHUP/s8HlgGIGn/7L07cBnw4+zzEFKg/NO2AsXMzKpLuUcq/UjNWYNJV3KdGRFv\nSBoNXBwRF2XlJgDXZKMtBi6KiM2SvgOcSmpSuz4irsvK3wR8NpumgHcjYkwndfCRipnZTuqqI5Wy\nQqUaOFTMzHZetTZ/mZmZvcehYmZmuXGomJlZbmoiVA44oNI1MDMzqJFQ8Ul6M7Pq4FAxM7PcOFTM\nzCw3DhUzM8tNzdz82NYGyv02HjOz2uSbH7fjppsqXQMzM6uZULn++krXYPs2b650DczMulbNhEpD\nlS/JO+/ABz6Qmui6usUxAp59FjZt6tr5mJl1VOW74h1XqVB54AH4wQ+27rdmDfzxj6n7m9+Eu++G\nffZpH97QACNHwu2379y8fvrT9lBasyZ1S3DLLXDssfClL8HXvpamf9hh0L07PPYYrF273UmbmeUj\nIvboFxAQMWZM7HZbtkSkXXzEXXel96lTI/r0Sd3nnts+HCKuvXbrzxBx223bn8+UKe8fb3uvk0/e\n+vNnPtNe5+XLIzZt6tJVY2ZVLu3+898n18zVX+PGwZ//3DXzKFxZVnx12aOPwvPPw6mn7vh0Cqv6\n2WfhmWfg5JPT540b4cUXYejQredx881wwQXtn/v2haeegs99Lr0//TS8+iq88gocfjj89a9p2ief\nnI5SIqClBf70J7iixI8/n3sujB8P55yz48tQ7O23YfFiOP74XRu/lIi0Dh58MB2ZbdkCI0bAEUfA\nBz8IxxwDvXrlNz+zeuXfU+lEIVQ+8Qn43//tmnkMGgSjR8Pvfw/f/W7aqV14YRrW1JS6//EfYdq0\n9PrhD+GSS9LwtrYUdp/4xPunu2gRjCn502Nbu+SSFCS9eqVms11xwQUppAAWLoQTTkhhBvCzn8F/\n/icsX771ODfcAMOGwbhx8NxzadlvuQV69IDhw+Hhh9vLHngg3HMPHH301tPYsAFeeAGuvDLNc9Mm\n+NWvUvPgIYekIJwxA+68MwVUa+uOL9PkyamZ77jjUuAOG5bqUe2Xlre1Vfc5wI0boVu39Cp8oWpr\ng/Xr05eVl16COXPgb39LZR94AHr2hCOPTH+Pl16Cfv3SF58lS9K2M3Qo7LsvrFuX/hcOPzxNb+zY\nNM4rr6T/3w99KD3Lb999U38JVq+Gd99N23+vXmm8zZvbX717p3pB+lKyYkXazo48Mi3DzohIy9Ta\nCo2N8Prrad5Dh6YvOM9kv0M7fHiabzl/x3ffhbfeStvuunVpvvvvn4a1taXzsJs3p2Hr1qVlj0jL\n1r9/Gnf9ethvvzS8Wzd44400fNOmNP23307jbd4Mffqkv8vAgWl9OVQ6UQgVyO8E+H//NzQ3p2/3\nEekEe2f+8Af41Kfe37+wwy5s7J356lfTN/Kjj047+2Jf+Qpcd13X7CTffTf9s196KczKfgR6zJj0\nT/LJT8J//Vfp8Zqa0tHPkUemcGhrS93FDj8cpkxJ/5i33tp5HZqa0sb/pz+lz9dck849bdmSAuPA\nA7cuHwH33Qfz56ed2tKlaf0W1nXB0KFpZwXw6U+nMH3ggfQPCOnobPXqVPePfzztKB55JP0TF4wY\nkcJq/Ph0dLR8earvO+/Am2+m81T77pvObW3cmM5pzZ2b/tEPOAB+8pO0Uxo5Mi3Pz36W6n/YYfC7\n36V5TJ2a+vXvn3bKxxyT6rNsGQwZkv4+t9ySjoi/8IU0zgMPwG23wahR8NGPpukvW5aW5dFHU+Cv\nXZuOWBsbU30aGuCoo9Ln1lYYMCCV/8Mf0t9uzRpYuTL1O+CAdPRb2Obffvv9f7cBA+BjH4ODD047\nqhNOSOM/9VSqw6ZNcOihqX49eqRpHHgg7L13Ws7Zs1P5N95I/2uQpnXooSlcNm2CVava/x577dV+\n0UmvXmkH+sEPpuVavz69H3lkmnfhnGPfvmnZe/aECRPSci1alLaXhoa0HtauTcNbW1O/YcPgtddS\nC0C/fmk6gwal9fbkk2n+hxySluPxx9O8+vVLr/32S/2HD091evvtNM/C36KhIW1DPXqk7aa1tT0s\nIS3Phg0pqArB0daW/g59+qRx/va3NH4hcHr0SGXffLM9cPr3T+ur8Nq4Mb2kFCyvvZbm27MnvP66\nQ6UkScG0StfCzGwPMw2HSil5HKm8/nr6Zj5rVvrW2vGbb2Hara0waRI89FD1N7FUUrU37+yKd96B\nJ55I3z732y99e92ypb2J6Omn07fUcm3alL6dvvtumkdDQ/qG2dhY/rTNirn5qxN5hErHgPjwh9M/\ndaH9dN48mDixjEqamVUZP6ali7z55vv7vfxyah/+4Q9TUDlQzMx2TN2FyoMPpiOTRx5JJzhLXZ56\nySWpzFe+svvrZ2a2J6upllopXW2xrfsYnngivY8c2X5ZMKQrUd54I13SV3z3u5mZ7biaO1IpXK7Z\nmeLr1gtPNj711HQp38EHO1DMzMpRc6Gy116dD4tI94V0dPfdXVcfM7N6UnOhsq0bFZctSzd2HXts\ne789/OI3M7OqUnOhsi2f+Ux6v+229P7971euLmZmtagmTtSPGZMeXwHpkQmdeeqp9H7IIek+lJ19\nLpCZmW1bWUcqkvpKmidphaS5knp3Um5wNnyZpMckDcn6j5e0WNJSSTMkNXQY71hJ70r67LbqUfx0\n4o7NWX/5C1x9dXoWz6RJ6dlK3bqlO5R9V7yZWb7KuqNe0nTg1Yj4jqTLgL4RMbVEuQXAVRHRLGkf\noA3YCKwEToyIZyRNA56PiJuzcRqAe4F3gJsj4ted1CEi4r2AaG6GE08sHr51+d/9Lj1k0MysnlXr\nHfWTgZlZ90xgSscCkkYA3SKiGSAi1kfEBqA/sDEisoehMB84vWjUrwK/AlbvTIW2bNn68xFHbP15\nzZqdmZqZme2MckNlQES0AkTEy8CAEmWGA2sl3ZU1dU1XOrxYAzRKGpWVOwM4CEDSgcCUiLgB2Kkk\nbWtL7y++CNdem35bodiZZ+7M1MzMbGdsN1Qk3Zud8yi8Hs3eTytRvFRbWiNwHPAN4FjgMOC8bNhZ\nwPckLQTWAYXjjGuBy4qrsb16Pvhgei8cqVx3HXzjG+8vt/fe25uSmZntqu1e/RURnT5OUVKrpIER\n0SppEKWbql4AlkTEymyc3wJjgRkR8Rfg+Kz/RNJRDcAxwC8kCdgPOEXSuxExu1Q9pk2bBqQf2Vmy\npIlTTmni6ae3t2RmZvWjpaWFlpaWLp9PHifqX4uI6Z2dqM9OuC8GToqIVyXdDCyKiBsk7R8Rr0jq\nDtwD/HtEtHQYfwZw9/ZO1AOcdlr62dzJk7c+Qb9gQTp5P3Ro+vlaM7N6V60n6qcDEyWtACYAVwNI\nGi3pRoCIaAO+BTRLeiQbL3vqFpdKWgYsAWZ1DJTMDqde4QeTOir8Pvxzz+3olMzMbFeUdfNjRLwG\nnFSi/2LgoqLP9wFHlyj3beDb25nH+Ttan4aG91/9NX16+/PAfv7zHZ2SmZntipq4o76gW7f3h8qX\nv5ze/YwvM7OuV1PP/uoYKmecAT17Vq4+Zmb1pqZCpaEhnVMpnFc5f4cbzszMLA81FSqFI5V589Ln\nU06pbH3MzOpNTYbK889XuiZmZvWp5kKlrS39LPDw4dsvb2Zm+aqpUClcUvzWW3DkkZWujZlZ/amp\nULnpJvjnf06h4qu+zMx2v5oKlYK333aomJlVQk2GyltvQY8ela6FmVn9qdlQ8ZGKmdnuV1Oh8tOf\npqcRX3VV+m16MzPbvWoqVPr3bz9CmTOnsnUxM6tHNRUqxfbfv9I1MDOrPzUVKhLcc0/qvvXWytbF\nzKwe1VyoFB4mOXp0ZetiZlaPaipUNm1q795778rVw8ysXtVUqNx5Z3t39+6Vq4eZWb2qqVAp/n36\nxpr6TUszsz1DTYVKQ9HSSJWrh5lZvaqpUDn88ErXwMysvikiKl2HskiKwjKsWgVDhqT+e/himZl1\nKUlERO5tOjV1pPKhD6X3G2+sbD3MzOpVTYVK4TyKr/wyM6uMmgyV5csrWw8zs3pVk6Hy1luVrYeZ\nWb2qqVAp8El6M7PKqMlQOeaYStfAzKw+lRUqkvpKmidphaS5knp3Um5wNnyZpMckDcn6j5e0WNJS\nSTMkNRSN0yTp4az8gp2p1+mnl7NUZma2q8o9UpkKzI+II4Bm4PJOyt0KTI+IjwBjgNWSBNwCnBkR\nHwNWAucBZOH0I+DUiDgK+NyOVmjKlPr9ffqWlpZKV6FqeF2087po53XR9coNlcnAzKx7JjClYwFJ\nI4BuEdEMEBHrI2ID0B/YGBHPZEXnA5/Nur8A3BURL2bjrNnRCv3mN1s/rqWe+B+mnddFO6+Ldl4X\nXa/c3e+AiGgFiIiXgQElygwH1kq6K2vqmq50G/waoFHSqKzcGcDgonH6SVogaZGkfyqznmZmthts\n91m+ku4FBhb3AgK4okTxUtddNQLHASOBVcAdpGauGcBZwPck7QXMA7YUjTMKGA/0AP4s6c8R8fT2\nF8nMzComInb5BSwHBmbdg4DlJcqMBRYUfT4buK5EuYnAL7Luy4Ari4b9BDi9kzqEX3755ZdfO/8q\nZ//f2avcXx2ZTTrqmA6cC8wqUWYR0EdS/4h4lXT0sQhA0v4R8Yqk7qQg+fdsnFnAdZK6Ad1JwXRN\nqQp0xQPRzMxs15R7TmU6MFHSCmACcDWApNGSbgSIiDbgW0CzpEey8W7K3i+VtAxYAsyKiJZsnCeA\nucBSYCFwY0QsK7OuZmbWxfb4R9+bmVn12KMvvpU0SdITkp6UdFml65M3SQdJapb0uKRHJf1r1r/T\nm04l/UDSU5KWSBpZ1P/cbD2tkHROJZYnD5IaJD0kaXb2eaikhdmy3S6pMeu/l6RfZOviz4UbbrNh\nl2f9l0s6uVLLUg5JvSXdmS3D45LG1ut2Ienr2U3SSyX9v+xvXxfbhaSfSmqVtLSoX27bgaRR2Xp9\nUtL3dqhSXXGiZne8SIH4NHAw8AFSE9qHK12vnJdxEDAy6+4JrAA+TGp2/Ha0X9RwddZ9CnBP1j0W\nWJh19wWeAXoDfQrdlV6+XVwnXwd+DszOPv8S+FzWfQNwcdb9L8D1Wfc/0H4RyEeAh0lXGA7NtiFV\nerl2YT3cAnwx627M/rZ1t10AHwKeBfYq2h7OrZftgvYra5cW9cttOwD+Ahybdf8e+NR261TplVLG\nyhwHzCn6PBW4rNL16uJl/i1wEvAEJa66A34M/ENR+eWky8HPAm4o6n9Dcbk95QUcBNwLNNEeKq8A\nDR23CeAPwNisuxuwutR2AswplNtTXkAv4JkS/etuu8hCZWW2Y2wkXTw0EVhdL9sF6Yt1cajksh1k\n4y4r6r9Vuc5ee3Lz14Gk+14KXsj61SRJQ0nfSBaSNpjim04L9xF1tk469n+RPXNdXQtcSrocEkn9\ngdcjXQwCW28D7y1zRGwh3YDbj9pYF4cAa5Sel/eQpBsl7UMdbhcR8RLwXeB5Uv3XAg8Bb9ThdlHQ\n8ab0Xd0ODszKdCy/TXtyqNQNST2BXwFfi4i3yHaqRTq72qJmLreW9GmgNSKWsPVy7egy1sy6oP3m\n4B9FxCjgbdI37XrcLvqQHhd1MOmopQcwaWcm0RX1qjK7dTvYk0PlRWBI0eeDsn41JTvB+CvgZxFR\nuA+oVdLAbPgg0qE+pOUfXDR6YZ3Uwrr6O+A0Sc8Ct5Pud/o+0FvtT7cuXq731kV2v1OviHiNztfR\nnuQFYFVEPJh9vosUMvW4XZwEPBsRr2VHHr8hbSt96nC7KMhrO9ildbInh8oiYJikg5Ue83IWqT21\n1txMatf8flG/wk2nZO+zivqfAyBpHKkJoJV0z8/E7IqhvqQ257ldX/X8RMT/iYghEXEo6W/dHBFn\nAwtof4p18Q24s7PPZMObi/qflV0FdAgwDHhgdyxDXrK/6SpJw7NeE4DHqcPtgtTsNU7SByWJ9nVR\nT9uF2PqoI5ftIGs6WytpTLZuz6H0De5bq/RJpjJPUE0iXRH1FDC10vXpguX7O9Lz0JaQrkx5KFvm\nfqSnOq8gPTOtT9E4PyRdufIIMKqo/3nZenoSOKfSy1bmejmB9hP1h5CuUHmSdMXPB7L+3UnPmXuK\ndB5qaNH4l2fraDlwcqWXZxfXwdGkL1ZLgF+Trtypy+0CuDL7Wy4lPS39A/WyXQC3AS8BG0kB+0XS\nRQu5bAfAaODRbNj3d6ROvvnRzMxysyc3f5mZWZVxqJiZWW4cKmZmlhuHipmZ5cahYmZmuXGomJlZ\nbhwqZmaWG4eKmZnl5v8DaMEi8FQdhU8AAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f1e36d4cc10>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Test Marginal_Likelihood_PG\n",
      "n_samples = 12500\n",
      "%time accum_marginal_loglikelihood = Marginal_Likelihood_PG(dN, sigma2e, n_samples, ret_accum=True, trials=trials)\n",
      "plt.plot(accum_marginal_loglikelihood)\n",
      "print(accum_marginal_loglikelihood[-1])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 11min 20s, sys: 25.9 s, total: 11min 46s\n",
        "Wall time: 11min 46s\n",
        "-25640.8372177\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEACAYAAACODmB7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHMdJREFUeJzt3X+QF/Wd5/HnS3DUxICiESog4g80QDSIWWOVRr+hIkH3\nFnN7pZJUThHXeOGuLpeyuEBInbOp7GZhL1kTc8S7OneDXhRRsglsAAFlkrtVcBHBH0x01KiokSTy\nw40ijMz7/ujPMB2cYWB6+jvf+c7rUdU1n/70p7s/356Z72v6093fUURgZmbW247q6w6YmVl9csCY\nmVkpHDBmZlYKB4yZmZXCAWNmZqVwwJiZWSkKBYykBZKaJW2WtFTSkFR/mqR3JG1K08LcOkdL+p+S\nnpW0VdK/TfUNkhZLapH0qKTRuXXmpvpmSVOK9NnMzKqj6BnMamBCREwEWoC5uWXPR8SkNM3K1c8D\ntkfEORExHvhFqr8R2BERY4HbgAUAksYD1wDjgCuAhZJUsN9mZlayQgETEWsjoi3NrgdG5RZ3FQIz\ngW/ntrEjFa8CFqXyA8DkVJ4GLI6I9yLiJbIgu7BIv83MrHy9eQ1mJrAyNz9G0uOS1km6BEDS0LTs\nW2nZfZI+nOpGAtsAImI/sFvSsHx98lqqMzOzGtZtwEhaI+nJ3PRU+vpnuTbzgNaIuCdVvQ6MjogL\ngFuAeyQdDwwmO8v5f2nZeuC/d7XrAq/LzMz62ODuGkTE5YdaLmkGcCUdQ1pERCuwM5U3SXoBODuV\n346If0xN7yc784HszORU4HVJg4AhEbFDUnt9u1GpbWd98QermZn1QET0+h/1Re8imwrMBqZFxN5c\n/cmSjkrlM4CzgBfT4uWSPp3KnwG2pvIy4PpUvhp4OFc/Pd1ldnra1mNd9Ski+u1066239nkf3P++\n78dA7H9/7ns99L8s3Z7BdON2oAFYk27sWh/ZHWOXAt+UtA9oA26OiF1pnTnA3ZL+DvgdcEOqvzPV\ntwBvAtMBImKrpCVkQdQKzIoyj4iZmfWKQgET2S3FndX/BPhJF8teAS7rpH4v2e3Ina3zbXJ3npmZ\nWe3zk/w1pFKp9HUXCnH/+1Z/7n9/7jv0//6XRfU02iTJo2dmZkdIElFrF/nNzMy64oAxM7NSOGDM\nzKwUDhgzMyuFA8bMzErhgDEzs1IUfZLfzOyIrF8Pl18ObW0ddfn/8NRV+XDb9eb6Ze3z85+H+fOp\new4YM6uq7dvhU5+CJUs66vKPr3VVrma7svf1wQ8yIDhgzKyq2trgmGPg+OP7uidWNl+DMbOq2r8f\nBg3q615YNThgzKyqHDADhwPGzKpq/344yu88A4K/zWZWVT6DGTgcMGZWVQ6YgcMBY2ZV1dbmIbKB\nwt9mM6sqn8EMHH4Opp974QV4/fWOJ4Sl7qfDbVfGNqvRrr3OapMDZuBwwPRzU6fC3r1w6qnZfET3\nUz2360pnH9dxuF+9Tu+u8+ab8Od/jg0ADph+7rjjYOlSOO+8vu5J7cmHUFd1XX09krbVWqee+jRm\nDDYAOGD6uUP91T7QebjMrG/5In8d8BuomdWiQgEjaYGkZkmbJS2VNCTVnybpHUmb0rQw1R8v6YlU\n94Sk30n6blrWIGmxpBZJj0oandvP3FTfLGlKkT7XG5/BmFmtKnoGsxqYEBETgRZgbm7Z8xExKU2z\nACLiDxFxfqo7H3gZWJra3wjsiIixwG3AAgBJ44FrgHHAFcBCyX+zm5nVukIBExFrI6L93watB0bl\nFh8yBCSdDXw4Iv45VV0FLErlB4DJqTwNWBwR70XES2RBdmGRftcbx62Z1aLevAYzE1iZmx8j6XFJ\n6yRd0kn7a4H7cvMjgW0AEbEf2C1pWL4+eS3VGR4iM7Pa1e1dZJLWAMPzVUAA8yJieWozD2iNiHtS\nm9eB0RGxU9Ik4KeSxkfEH3LbmQ588VC7PoLXcUBjY+OBcqVSoVKp9GQz/YrPYMzsSDQ1NdHU1FT6\nfhQF/wSWNAO4CZgcEXu7aLMOuCUiNqX584AlEfHRXJuVQGNEbJA0CPhNRJwiaQ4QETE/tVsF3BoR\nGzrZTxR9Pf3N+PFw//0wYUJf98TM+itJRESv/6la9C6yqcBsYFo+XCSdLOmoVD4DOAt4Mbfq54F7\nD9rccuD6VL4aeDiVlwHT011mp6dtPVak3/XGZzBmVouKPmh5O9AArEk3dq1Pd4xdCnxT0j6gDbg5\nInbl1rsauPKgbd0J3C2pBXiTbAiNiNgqaQmwFWgFZg2405RD8JEws1pVeIislgzEIbJx47KPihk/\nvq97Ymb9VVlDZP6omE789rfwzjvZJ74eztTX/9vCQ2RmVoscMJ0YPhxGjMjCY//+7ic4vBCSsq+H\nUz7cti+/7I8+N7Pa5IDpxHHHZf9n5QMfOLz2bW3dh1BE1q79a5Fyvm7wYBg7ttzjYWbWEw6YThzp\nZZz2M4qjjy6nP2Zm/ZE/TbkTEb6uYWZWlAPGzMxK4YDphM9gzMyKc8B0wQFjZlaMA6YTA+xZTTOz\nUjhguuAzGDOzYhwwnfAZjJlZcQPqOZgIeOWVrNzQ0DEdc0z2DEv+rMVnMGZmxQyogHnhheyp99Gj\nYd8+2Ls3+7pvH7S2ZiHT0JA9Hd/Xny9mZtbfDaiA2bs3+9ThZ555/7KILGT27cvOXvz5XmZmxQyo\ngDkUqWPIzMzMivNAkJmZlWJABYzvDjMzq54BFTDgu8PMzKplwAWMmZlVhwPGzMxK4YAxM7NSDKiA\n8UV+M7PqGVABA77Ib2ZWLYUCRtICSc2SNktaKmlIqj9N0juSNqVpYW6dz0t6Mq2zQtKwVH+ipNWS\nnpX0oKShuXW+L6klrTOxSJ/NzKw6ip7BrAYmRMREoAWYm1v2fERMStMsAEmDgNuAy9I6TwH/KbWf\nA6yNiHOAh9u3JekK4MyIGAvcDNxRsM9mZlYFhQImItZGRFuaXQ+Myi3ubDCqve5DkgQMAV5LdVcB\ni1J5UZpvr78r7W8DMFTS8J71tydrmZlZT/TmNZiZwMrc/BhJj0taJ+kSgIh4D5hFdubyKjAO+PvU\n/pSI2J7avQG0h8hIYFtuu6+luh7xNRgzs+ro9sMuJa2h480esrOQAOZFxPLUZh7QGhH3pDavA6Mj\nYqekScBPJY0H3gW+DHw8Il6SdDvZUNhf8/4znh6dbzQ2Nh4oVyoVKpVKTzZjZla3mpqaaGpqKn0/\nioLjRpJmADcBkyNibxdt1gG3kJ0xfTsiLk/1nwK+FhH/RlIzUImI7ZJGAOsiYpykO1L5vrTOr8iu\n4WzvZD9xqNfz5JPwxS9mX83MLCOJiOj18Z2id5FNBWYD0/LhIulkSUel8hnAWcCLZMNb4yWdlJpe\nDjSn8jJgRirPAH6Wq78ubesiYFdn4XI4fA3GzKx6iv4/mNuBBmBNds2e9emOsUuBb0raB7QBN0fE\nLmCXpL8E/m9a9jIdoTIfWCJpZqq/BiAiVki6UtLzwNvADQX7bGZmVVB4iKyWdDdEtmULXHdd9tXM\nzDI1OURmZmbWFQeMmZmVYkAFTB2NBpqZ1bwBFTDgBy3NzKplwAWMmZlVhwPGzMxKUfQ5mJqzfz9c\ncAGcdBKMGAEnnwzHH59NO3f2de/MzAaOuguY1tbsOZfVq2H7dvj97+Htt+Gtt2DPHviLv+jrHpqZ\nDQx196Dlnj3BCSfAu+/2dW/MzPoHP2hpZmb9St0FTB2dkJmZ9Wt1FzDgZ13MzGpBXQaMmZn1PQeM\nmZmVou4CxtdgzMxqQ90FDPgajJlZLajLgDEzs77ngDEzs1I4YMzMrBR1FzC+yG9mVhvqLmDAF/nN\nzGpBXQaMmZn1PQeMmZmVolDASFogqVnSZklLJQ1J9adJekfSpjQtzK1zraQtkp6S9O1cfYOkxZJa\nJD0qaXRu2dxU3yxpyqH65GswZma1oegZzGpgQkRMBFqAubllz0fEpDTNApA0DFgAfDoizgVGSPp0\nan8jsCMixgK3pXZIGg9cA4wDrgAWSoe+yuJrMGZmfa9QwETE2ohoS7PrgVG5xZ29zZ8BPBcRO9L8\nQ8C/S+WrgEWp/AAwOZWnAYsj4r2IeIksyC4s0m8zMytfb16DmQmszM2PkfS4pHWSLkl1zwPnSBot\naTDwOeDUtGwksA0gIvYDu9MZz4H65LVUZ2ZmNWxwdw0krQGG56uAAOZFxPLUZh7QGhH3pDavA6Mj\nYqekScBPJY2PiF2SvgwsAfYDjwBndrXrnrygv/qrRvbtg8ZGqFQqVCqVnmzGzKxuNTU10dTUVPp+\nFAWvikuaAdwETI6IvV20WQfcEhGbDqq/CTgzIuZIWgXcGhEbJA0CfhMRp0iaA0REzE/rHGjXyX7i\nrbeCj3wE/vVfC70sM7MBQxIR0etXr4veRTYVmA1My4eLpJMlHZXKZwBnAS+m+Q+nrycCs4D/nVZb\nBlyfylcDD+fqp6e7zE5P23qsSL/NzKx83Q6RdeN2oAFYk27sWp/uGLsU+KakfUAbcHNE7ErrfE/S\nx8mG2f4yIp5P9XcCd0tqAd4EpgNExFZJS4CtQCswKw5x2vXLX8KePQVflZmZFVZ4iKyWSIrRo4Md\nOzxEZmZ2uMoaIqu7gNm9OxgypK97YmbWf9TkNZha5IcszcxqgwPGzMxK4YAxM7NS1F3AmJlZbai7\ngPEZjJlZbXDAmJlZKRwwZmZWCgeMmZmVou4CxszMakPdBYzPYMzMaoMDxszMSuGAMTOzUtRdwJiZ\nWW2ou4DxGYyZWW1wwJiZWSkcMGZmVgoHjJmZlaLuAsbMzGqDA8bMzErhgDEzs1I4YMzMrBQOGDMz\nK0WhgJG0QFKzpM2Slkoaklt2nqRHJD0taYukhlQ/SdKTkp6TdFuu/YmSVkt6VtKDkobmln1fUkva\nz8QifTYzs+ooegazGpgQEROBFuDrAJIGAXcDX4qIjwEVoDWt80Pgxog4Gzhb0mdT/RxgbUScAzwM\nzE3bugI4MyLGAjcDdxTss5mZVUGhgImItRHRlmbXAyNTeQqwJSKeTu12RkRIGgF8KCL+JbW7C/hc\nKl8FLErlRWm+vf6utJ0NwFBJw4v028zMyteb12BmAitS+WwASaskbZQ0O9WPBF7NrfMqHaE0PCK2\nA0TEG8Dw3Drbcuu8llvHzMxq1ODuGkhaQ8ebPYCAAOZFxPLUZh7QGhH35rZ7MfAJ4F3gIUkbgbeO\noG9xBG0PaGxsPFCuVCpUKpWebMbMrG41NTXR1NRU+n4U0aP38Y4NSDOAm4DJEbE31V0LTI2IG9L8\nN4A9wI+BdRExLtVPBy6LiC9LagYqEbE9DaWti4hxku5I5fvSOr9K62zvpC9R9PWYmQ00koiIXv+g\nraJ3kU0FZgPT2sMleRA4V9KxkgYDlwHPpKGv3ZIulCTgOuBnaZ1lwIxUnnFQ/XVpfxcBuzoLFzMz\nqy2FzmAktQANwJupan1EzErLvkB2V1kb8POIaL8r7ALgR8CxwIqI+EqqHwYsAU4FXgauiYhdadkP\ngKnA28ANEbGpi/74DMbM7AiVdQZTeIisljhgzMyOXE0OkZmZmXXFAWNmZqVwwJiZWSkcMGZmVgoH\njJmZlcIBY2ZmpXDAmJlZKRwwZmZWCgeMmZmVwgFjZmalcMCYmVkpHDBmZlYKB4yZmZXCAWNmZqVw\nwJiZWSkcMGZmVgoHjJmZlcIBY2ZmpXDAmJlZKRwwZmZWCgeMmZmVwgFjZmalcMCYmVkpCgWMpAWS\nmiVtlrRU0pDcsvMkPSLpaUlbJDWk+m9JekXSWwdtq0HSYkktkh6VNDq3bG6qb5Y0pUifzcysOoqe\nwawGJkTERKAF+DqApEHA3cCXIuJjQAVoTessA/6kk23dCOyIiLHAbcCCtK3xwDXAOOAKYKEkFey3\nmZmVrFDARMTaiGhLs+uBkak8BdgSEU+ndjsjIlL5sYjY3snmrgIWpfIDwORUngYsjoj3IuIlsiC7\nsEi/zcysfL15DWYmsCKVzwaQtErSRkmzD2P9kcA2gIjYD+yWNCxfn7xGR5CZmVmNGtxdA0lrgOH5\nKiCAeRGxPLWZB7RGxL257V4MfAJ4F3hI0saIWHcEfevRMFhjY+OBcqVSoVKp9GQzZmZ1q6mpiaam\nptL3ozRy1fMNSDOAm4DJEbE31V0LTI2IG9L8N4A9EfGd3HpvRUT+poCVQGNEbEjXcH4TEadImgNE\nRMxP7VYBt0bEhk76EkVfj5nZQCOJiOj1a9tF7yKbCswGprWHS/IgcK6kYyUNBi4Dth68+kHzy4Hr\nU/lq4OFUXgZMT3eZnQ6cBTxWpN9mZla+otdgbgeOB9ZI2iRpIUBE7AK+C2wENgGPR8RKAEnzJW0D\njku3K/+3tK07gZMltQD/BZiTtrUVWEIWUCuAWT5NMTOrfYWHyGqJh8jMzI5cTQ6RmZmZdcUBY2Zm\npXDAmJlZKRwwZmZWCgeMmZmVwgFjZmalcMCYmVkpHDBmZlYKB4yZmZXCAWNmZqVwwJiZWSkcMGZm\nVgoHjJmZlcIBY2ZmpXDAmJlZKRwwZmZWCgeMmZmVwgFjZmalcMCYmVkpHDBmZlYKB4yZmZXCAWNm\nZqVwwJiZWSkKBYykBZKaJW2WtFTSkNyy8yQ9IulpSVskNUg6TtI/pXWekvTXufYNkhZLapH0qKTR\nuWVzU32zpClF+mxmZtVR9AxmNTAhIiYCLcDXASQNAu4GvhQRHwMqQGta528jYhxwPnCJpM+m+huB\nHRExFrgNWJC2NR64BhgHXAEslKSC/TYzs5IVCpiIWBsRbWl2PTAylacAWyLi6dRuZ2T2RMQvUt17\nwCZgVFrnKmBRKj8ATE7lacDiiHgvIl4iC7ILi/TbzMzK15vXYGYCK1L5bABJqyRtlDT74MaSTgD+\nDFibqkYC2wAiYj+wW9KwfH3yGh1BZmZmNWpwdw0krQGG56uAAOZFxPLUZh7QGhH35rZ7MfAJ4F3g\nIUkbI2Jdaj8IuAe4LSJe7mrXPXg9NDY2HihXKhUqlUpPNmNmVreamppoamoqfT+KiGIbkGYANwGT\nI2JvqrsWmBoRN6T5bwB7IuI7af5O4K2I+GpuOyuBxojYkALoNxFxiqQ5QETE/NRuFXBrRGzopC9R\n9PWYmQ00koiIXr+2XfQusqnAbGBae7gkDwLnSjpW0mDgMmBrWudbwJB8uCTLgetT+Wrg4VReBkxP\nd5mdDpwFPFak32ZmVr5CZzCSWoAG4M1UtT4iZqVlXyC7q6wN+HlEzJXUfj2lGdhHNtT2g4j4e0nH\nkN15dn7a3vR0UR9Jc8nuMmsFvhIRq7voj89gzMyOUFlnMIWHyGqJA8bM7MjV5BCZmZlZVxwwZmZW\nCgeMmZmVwgFjZmalcMCYmVkpHDBmZlYKB4yZmZXCAWNmZqVwwJiZWSkcMGZmVgoHjJmZlcIBY2Zm\npXDAmJlZKRwwZmZWCgeMmZmVwgFjZmalcMCYmVkpHDBmZlYKB4yZmZXCAWNmZqVwwJiZWSkcMGZm\nVopCASNpgaRmSZslLZU0JLfsPEmPSHpa0hZJDal+paQnJD0laaEkpfoTJa2W9KykByUNzW3r+5Ja\n0n4mFumzmZlVR9EzmNXAhIiYCLQAXweQNAi4G/hSRHwMqACtaZ2rI+L8iDgXOAW4OtXPAdZGxDnA\nw8DctK0rgDMjYixwM3BHwT7XrKampr7uQiHuf9/qz/3vz32H/t//shQKmIhYGxFtaXY9MDKVpwBb\nIuLp1G5nREQq/wFA0tFAAxBpnauARam8KM2319+V1t0ADJU0vEi/a1V//yF1//tWf+5/f+479P/+\nl6U3r8HMBFak8tkAklZJ2ihpdr6hpFXAG8BbwAOp+pSI2A4QEW8A7SEyEtiWW/01OoLMzMxq1ODu\nGkhaQ8ebPYDIzjrmRcTy1GYe0BoR9+a2ezHwCeBd4CFJGyNiHUBETE3XZH4MTAYe6mTX0UmdmZn1\nFxFRaAJmAP8MHJOruxb4h9z8N4BbOln33wPfT+VmYHgqjwCaU/kO4NrcOr9qb9fJ9sKTJ0+ePB35\nVDQLOpu6PYM5FElTgdnApRGxN7foQWC2pGOB94DLgO9K+iDwoYh4Q9Jg4E+BX6Z1lpGF1fz09We5\n+v8I3CfpImBX+1DawSJCRV6PmZn1HqW//Hu2stRCdqH+zVS1PiJmpWVfILurrA34eUTMlXQK8E9p\nnaOAdcBXI6JN0jBgCXAq8DJwTUTsStv6ATAVeBu4ISI29bjTZmZWFYUCxszMrCt18yS/pKmSfiXp\nOUlf6+v+AEgaJelhSc+kB0v/c6o/4odKJV2fXtuzkq6r8us4StImScvS/BhJ61N/7k3DnUhqkLQ4\n9f9RSaNz25ib6pslTali34dKuj/t9xlJn+xPx1/SV9PDyk9K+nE6xjV7/CXdKWm7pCdzdb12vCVN\nSsfiOUm3Van/h3qgvNPj2tX7UVffu7L6nlt2i6T2kaL2uvKPfRkXdqo9kQXl88BpwNHAZuCjNdCv\nEcDEVD4eeBb4KNl1pv+a6r8G/E0qX0E2nAjwSbIhR4ATgReAocAJ7eUqvo6vAv8HWJbm7yN7YBbg\nh8DNqfxlYGEqXwssTuXxwBNkdxeOSd8rVanvPyIbViXtf2h/Of7AR4AXgYbccb++lo8/cAkwEXgy\nV9drxxvYAPxJKq8APluF/n8GOCqV/wb49qGOK4d4P+rqe1dW31P9KGAV8GtgWDWPfem/4NWYgIuA\nlbn5OcDX+rpfnfTzp+mH9cCdcBz6jrlmslvEpwM/zNX/MN+u5D6PAtaQfRpDe8D8LvcLd+DYpx/i\nT6byIOC3nX0/gJXt7Uru+xDghU7q+8XxJwuYl9Mv/WCyG14uB35by8ef7I01/wbdK8c7rbs1V/9H\n7crq/0HLPgfcfajjyiHejzr53VlVdt+B+4Fz+eOAqcqxr5chsoMfxnyVGnsYU9IYsr8u1pP9sh3O\nQ6Xtr6MvHzb9O7I7BQNA0knAzuj4BIf8sT7Qz4jYD+xOp+R91f/Tgd9L+oc0xPe/JH2AfnL8I+J1\n4DvAK2mfu4FNZHdS9ofj3+5wH6Lu7niPTG0Obl9N+QfKD9XP972uLn53PlJmZyVNA7ZFxFMHLarK\nsa+XgKlpko4n+8SCr0T2UTkH31nR1Z0WfXrbtaQ/BbZHxOaD+nK4/err28YHA5OA/xERk8juQpxD\n/zn+J5B9VNJpZG9EHyS7m/KwN1FGv3pBTR7v7uj9D5T3aDO91Z9udyQdR3Yn762H07yMPtRLwLwG\njM7Nj0p1fS5dxHuA7LS6/dme7UqfpyZpBNmQB2R9PjW3evvr6KvXdzEwTdKLwL1kn7rwPbLPg2v/\n2cn35UD/lX3g6ZCI2EHXr6tsr5L99bYxzS8lC5z+cvw/A7wYETvSGck/kn1PTugnx79dbx3vPnsd\nkmYAVwJfyFUfUf8j4k26/t6V4Uyya0NbJP067W+TssdFqnPsyxqHreZENt7cflGtgeyi2ri+7lfq\n213Adw+qm0/HmOwcOi56XknHhbeL6PzCW3v5hCq/jsv444v810bHGO1/SOVZdFxkns77LzI3kA1b\nVfMi/y+As1P51nTs+8XxBy4EngKOJfsL80ekh45r+fiTvak9VcbPO9kQ84XpeKwAplah/1OBZ4CT\nDmrX6XGl8/ej/EX+933vyur7Qct+DZxYzWNf+i94tab0Q/As2b8NmNPX/Ul9uhjYn37AniAbP58K\nDAPWpv6uzr9ZAT9IP5xbgEm5+hnptT0HXNcHryUfMKeT3VHyXPqFOTrVH0P2sGxL+mEck1t/bnpd\nzcCUKvb748C/pO/BT9IvTr85/mSh2Aw8SfYp40fX8vEH7gFeB/aSXTu6Ib1R9crxBi4gC90W4HtV\n6n8L2c0Wm9K0sLvjShfvR11978rq+0HLXyRd5K/WsfeDlmZmVop6uQZjZmY1xgFjZmalcMCYmVkp\nHDBmZlYKB4yZmZXCAWNmZqVwwJiZWSkcMGZmVor/D8VOdkajo9+ZAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f1e33b18fd0>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Laplace approximation at the MAP"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Do_Kalman_Likelihood_Bernoulli_LaplaceMAP(dN, sigma2e, tol=1e-8, trials=1.):\n",
      "    \"\"\"MAP solution, inverse covariance matrix, and marginal loglikelihood of state-space model\n",
      "    computed using Laplace approximation around MAP state.\n",
      "\n",
      "    :param dN: Observations (K,)\n",
      "    :param sigma2e: Variance of process noise\n",
      "    :param tol: Convergence criterion on the gradient of the log-likelihood\n",
      "    :param trials: Number of trials for binomial observations (1 for Bernoulli)\n",
      "    :return: x_map, U, marginal_loglikelihood, joint_loglikelihood\n",
      "    \"\"\"\n",
      "    x = np.zeros(dN.shape)\n",
      "    dN = dN.astype(float)\n",
      "    while True:\n",
      "        # Build gradient of joint\n",
      "        d2x = np.convolve(x, [-1, 2, -1])[1:-1]\n",
      "        d2x[-1] -= x[-1]\n",
      "        G = -dN + trials * (1. / (1. + np.exp(-x))) + d2x / sigma2e\n",
      "        # Build Hessian of joint\n",
      "        D = trials / (np.exp(x) + 2. + np.exp(-x)) + 2. / sigma2e\n",
      "        D[-1] -= 1. / sigma2e\n",
      "        B = -np.ones(len(D)) / sigma2e\n",
      "        B[0] = 0.\n",
      "        U = sp.linalg.cholesky_banded((B, D), lower=False)\n",
      "        # Check convergence\n",
      "        if np.dot(G, G) < tol:\n",
      "            x_map = x\n",
      "            break\n",
      "        # Update estimate of map\n",
      "        x -= sp.linalg.cho_solve_banded([U, False], G)\n",
      "\n",
      "    # Compute joint and marginal probabilities\n",
      "    joint_loglikelihood = (np.sum(np.log(sp.special.binom(trials, dN)) + dN * x_map - trials * np.log(1 + np.exp(x_map))) -\n",
      "                           .5 * ((np.sum(np.diff(x_map)**2) + x_map[0]**2) / sigma2e + len(dN) * np.log(2*np.pi*sigma2e)))\n",
      "    marginal_loglikelihood = len(dN)/2. * np.log(2*np.pi) + joint_loglikelihood - np.sum(np.log(U[-1]))\n",
      "    return x_map, U, marginal_loglikelihood, joint_loglikelihood"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Test Marginal_Likelihood using the Laplace approximation around the MAP\n",
      "%time x_map_l, U_l, marginal_loglikelihood_l, joint_loglikelihood_l = Do_Kalman_Likelihood_Bernoulli_LaplaceMAP(dN, sigma2e, trials=trials)\n",
      "print(marginal_loglikelihood_l)\n",
      "plt.plot(x_map_l)\n",
      "plt.plot(dN)\n",
      "plt.show()\n",
      "plt.plot(np.exp(x_map_l)/(1+np.exp(x_map_l)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 8 ms, sys: 0 ns, total: 8 ms\n",
        "Wall time: 9.81 ms\n",
        "-3102.74986369\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEACAYAAAC+gnFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYFMX5x7+1u9wg98qxwCKnogYVUUTNemCIxqgY45Go\nGKN4J9EYNR6AV9RfvBWJZ4wa8QZEI0cCKKAcIi73fbMsywLLwrLAztTvj3dqpqen7+npmdl5P8+z\nz8z0Tk1Vd1e/b71HVQkpJRiGYZjcJS/dDWAYhmHSCysChmGYHIcVAcMwTI7DioBhGCbHYUXAMAyT\n47AiYBiGyXF8UQRCiDeEEOVCiFLNsZFCiC1CiIWRv6F+1MUwDMP4i18WwVsAfmZw/Bkp5YmRv698\nqothGIbxEV8UgZRyFoDdBv8Sfvw+wzAMkzpSHSO4VQixSAjxuhCiZYrrYhiGYTyQSkUwBkAPKWV/\nANsBPJPCuhiGYRiPFKTqh6WUFZqPrwH43Oh7Qghe7IhhGMYDUkpf3O9+WgQCmpiAEKKD5n/DACwx\nKyilrLd/I0eOTHsb+Pz4/HLt3HLh/PzEF4tACPFvACUA2gohNgEYCeAsIUR/AGEAGwCM8KMuhmEY\nxl98UQRSyqsMDr/lx28zDMMwqYVnFqeYkpKSdDchpfD5ZS/1+dyA+n9+fiL89jW5boAQMt1tYBiG\nyTaEEJAZGCxmGIZhshBWBAzDMDkOK4KAuG7CdbjykyvT3QyGYZgEOEYQEA0eaYC6cB3kyPp/rgzD\npB6OEWQhgtffYxgmQ2FFEBBCsCJgGCYzYUUQEGwRMAyTqbAiYBiGyXFYEQQEu4YYhslUWBEEBLuG\nGIbJVFgRMAzD5DisCBiGYXIcVgQBwTEChmEyFVYEAcExAoZhMhVWBAHBFgHDMJkKKwKGYZgchxVB\nQLBriGGYTIUVQUCwa4hhmEyFFQHDMEyOw4ogIA6HDqe7CQzDMIawIgiIA3UH0t0EhmEYQ1gRMAzD\n5DisCBiGYXIcVgQMwzA5DisChmGYHIcVAcMwTI7DioBhGCbHYUXAMAyT47AiYBiGyXFYETAMw+Q4\nrAgYhmFyHFYEDMMwOQ4rAoZhmByHFQHDMEyO44siEEK8IYQoF0KUao61FkJMEUKsFEJMFkK09KMu\nhmEYxl/8sgjeAvAz3bF7AUyTUvYB8D8A9/lUF8MwDOMjvigCKeUsALt1hy8C8Hbk/dsALvajLoZh\nGMZfUhkjKJRSlgOAlHI7gMIU1lWvKd9Xjtq62nQ3I2Vs3bsVdeE6V2UqaypRfbA6RS3KDLbs3YJQ\nOJTuZtiyqWpTupvAJElBgHVJs3+MGjUq+r6kpAQlJSUBNCd76PB0B1x/wvV4/Zevp7spKaHo2SI8\nde5TuHvw3Y7LdHqmE04tOhUzh89MYcvSS5dnu2DM+WNw88k3p7sppny7+Vuc9uZpkCNNH2/GJ2bM\nmIEZM2ak5LdTqQjKhRBHSinLhRAdAOww+6JWETDGlO0rS3cTUsrOmp2uvn8odAgb9mxITWMyiMoD\nlelugiX7Du1LdxNyBv0gefTo0b79tp+uIRH5U0wEMDzy/loAE3ysi2EgZf0fhebCOTLpx6/00X8D\nmAOgtxBikxDiOgBPABgihFgJ4JzIZ8YjLBASkebexnpDLpwjk358cQ1JKa8y+de5fvw+wxiRC8ox\n089RCGH/JSbj4ZnFWQI/cImEZTjdTUg5bBEwQcCKgMlackFIZrpFwNQPWBEwWUsuCMlMV3YCbKnW\nB1gRBEwuuDO84MX1Vb6/HJU17tMrq2qrAp0EtffgXmzcs9FTWbfK7lDoEFbuXOmprmxh3e51mLtl\nbrqbUa9gRRAwny3/LN1NqFec9+55rsv85tPfoNtz3VLQGmOGjx+O4ueLPZV1axE8991z6PtyX091\nZQs9XuiBU984Nd3NqFewIgiYg6GD6W5CvWLHftN5iqYEPUlr14Fdnsu6tQjq+7IbTGpgRRAw+SI/\n3U2oV2RDnCCZjK9MjxEw9QNWBAGTn8eKwE+8CMqgA5x5wvtjlumKjtOa6wesCAImGaHAJJLpghJI\nTvGwRcAEAUulgGFF4C/ZICiTcg1lgaJLF3xt/IOlUsBwjMBfvAiDoN0Z9TnXvj6fWy7BiiBgOEbg\nL/XeInB5ftlwPZjMgxVBAOw/tD/63u0EKDXiXbpjqes0xHlb52F15WrX9b1b+i42V212VQ4Adh/Y\njX8s+IerUbpK/1y0fZHr+gBvFoHXCVdl1WVYs2uN63LJuAO9Xpegmb91vuPvlu8rx6rKVfhm4zdJ\n1clKzz9YEQTA32b9Lfr+2vHXuipbvr8cALCxaiOun3i9q7KnvH6K6zKrd63G1Z9djQGvDXBVDgBG\nTBqBm764CYt3LHZc5tIPLwUATF472XV9gDdh4HUewbnvnIteL/ZyXS4Z98nUdVMDq8sLytoZ+PpA\nx2UuGncR+rzUB2f+88yk5lgw/sGKIAAOhQ5F37s29TUj3gOHD7iuu6KmwlN9WivGKTWHa1yX8XJO\nWoIMGHrdjSvImEQ2jJL3H471LY4xZAasCAIgnemDbtc2Um6MZISXm4B4sjGTIAWfV6XDwi4e7fXg\njKrMgBVBhqPt7F4eGrcPSzKCWQllN7+RbBZVNgiDIC2CbFM62dbe+gorgizCi9BzO2JWD2YyD6ib\n4Giy8yqCXM3Vq/VRn4Wdl3PTKkaemZwZsCIIgHSuNePVNeQFpajc/EYupNPyJMJ44lxDPOs6I+Ae\nmuEk6/pwqwj8UFquYgTJuoayQBhwsJjJdFgRBIBfox4vD7lX14mXurxYBMmOlrMiRlCPXUNeYHdQ\n5sGKIMPRCrrDocOOy32+8nMAQCgcclWfUhxuUyW379uOb7d8C4DmPDglWddQ9aFqVOx3niK7uNz5\nHAd9uS17twBwd02raqvw2Qr3mxF53YFLKR0pJSaunOiq7I/bf8R1E64LdBKb2wHH1LWxeRXZMAjI\nFlgRBICXzVMU2gdl9S7ns4R/Oe6XAICt1Vvd1efx4RoxaQT2HtwLABjyzhDH5fwYLd/+n9sdf/f4\nscd7qkNbbs7mOY7LaScTusHrDlyqv+w7tA8XjbvIVdn+/+iPfy76J074xwmOy3gZ3Sdzz73sSMfY\nw4ogiwhiBOTVx6xtmxt3lB+B1KD94m7qc2uR+UW2xAp4VJ8ZsCIIAL8mzQSRKun1wfQq0Nlf7C9a\n11B9J1uUXTbAiiDDSTZY7BavyiaXBLobIZsuayWThWQu9ZVsgRVBACSVNRS0ReBRgGgtAjeCMhsz\natxco3SNzLPFIshkhZVLsCIIAL8mlAUSI/DBNeTmfOv76LC+CzpPM4t9Uv7ZouyyAVYEWUQQFoHX\nOuIUgYsHPSstAjeuoXRZBFmigFiYZwasCDIc7YMSxMMd9Ho69d0iSBeZLGD9uufZouyyAVYEacBN\nSqHeNfT6wteTXsPfite+fy36/tPlnzoup7UIDoedT3zTKpB5W+c5LqdFu99DEKzbvc7xd99a9FbS\n9bkR6o998xiVcSkk1RwQt2iF+qxNs5yV0dzzrs919bSPBeMvrAgCQO9ucSNI9MHiGz6/ATM2zHBc\nvqS4xPF3AeC5uc9F36vdw5zgdZSnLTfojUGefiOZCXteuHPKnY6/W3WwKun6tBu5pIoPlnyQ9G+c\n8dYZrsvsO7QPU9ZOSbpuJjlYEQRAMmZ6sumjDfMbeq7bDenctCVoN0jQq4l6idu4vSaB7pmgq8vz\nelgZ7P7KNlgRBIBfvkwvHT/Tg7F+CKCgfcVBX1NPiiCga+LHtQhyTwnGGFYEAZDMQxn0PIKgyUaL\nIOgAdzIWgdNrE6Ry09eVrmU4mBgFqa5ACLEBQBWAMIDDUsqBqa4z0/DbNeRqQlOGZ1Zko0WQFa6h\nDL/vWoJcKp0xJuWKAKQASqSUuwOoKyPJBYsgnds4Bm4RBOwa8jJijloEkBnvHszkfp0rBDG0EQHV\nk7Ek09GDnlnsFc/B4iycRxC0RRCS7hVBUMLVywDAr2Ax4x9B9GgJYLIQYr4Q4oYA6st4vK5V4+WB\nyfTRoC8WQT13DdWF61yXiboRHSpor9fQDyHOWUPpJ4gePVhKOQDA+QBuFUKcHkCdGcW7pe/GfV5d\n6XyDmRU7V0Tfe0kFnbx2Mro/3x21dbWuywbBB0tj+eteBcKCbQs8lVuza42ncmX7yjwJITf3Xcv7\ni993XabLs10AxPcfK7xMfNtUtQnPfPuM63J6V5dTi0ftEMf4T8oVgZSyLPJaAeAzAAnB4lGjRkX/\nZsyYkeompR3tpC07tDtc3XHKHZ7q27BnA1ZVrvJU1inZGLh7ce6Lnst6UaxehCYA3Pvfez2VA4An\nZz/p6Htudl1TvLHwDUxYOcF1Of0kO6cDgH/9+C/XddUnZsyYEScr/SSlwWIhRFMAeVLKfUKIZgDO\nAzBa/z2/TyrTceNa0PpTlUXgZTSaajM6G830ZJSXJzddFsZDgsDptdS7EbNx8JEMJSUlKCkpiX4e\nPTpBlHom1VlDRwL4TAghI3W9J6XM+fnkXv3ifqWhMsnDAU7v6Pux02sZdGwml0ipIpBSrgfQP5V1\n5BKZLHyyUdEko1i9ZPJkeuA+XbAiSD98ZdOAGxdB3M5fPs1HYIhkrme2zIbNBneU02upPxfu0/7B\niiANeN24JZNdQ7n2UHqxCPgeGMMWQfrhK5sG3AgE7SgoqYlpqQ4W55prKAMtAqN7kInuKH07WRGk\nH76yWYTTtYaMBFyqBXUmxy9SQSbGCNJx372gb6fTa6lXBJl4btkKK4IUY7R7VvOGzR2X1+7a5TQn\n3OgB8Tr6nb5+uu13QuEQxq8YH3fsnqn3eKovlWzftz3u85gFY/DnKX+2Lbd+9/qEYx8u/dB1/WMW\njMEPZT+4LgcAwz4YZvsdr4Jxc9XmhGN7avfYlvOiDI24e+rdaPdUOxysO2j5vc9WfBb3eeLKib7U\nz7AiSDkV+ysSjl3S95LA2+FVSDz89cO23zHaQeupOU95qi+VGCm1p7992rac0Zadd025y1Mbxi4Y\n66mcXggaYWSVObFCvlj9RcKxDXs22Jbzc8vUygOV2Fmz0/I7+p35bvvyNt/qz3VYEdRDDF0EKYwR\nZItbKBNcCU6yeLzeK6NyXutzck+9+uz9ug9+WSQMK4K0kI7ArdeHz0lbsyVbJRPa6WSE7vle+Vgu\nG5R7NrQxW2BFkAbSkUbIFkFm4GSEHvRqnF77S7rve7rrr0+wIkgDqe7AfloETgRXtpjomeAacoJX\nge7nTl9Ofivd15MVgX+wIkgD6XBR+Dla1JOJOfWZihO/up8C3Ykryqg+J0I+3fsIZIKrr77AiiAN\npMU1lMI6s8YiyADBkdIYgY/K3olyT/eIPFv6XTbAiiAN1LdZvmwR+EvQm7kblXOyK1q6FWu6FVF9\nIusVgZQSvV/sHf3c56U+EKNFwiYWszfNhhgtsK16m6d6Hv36Uby96G3868d/4ZGZjyTVZtWBT3r1\nJOw9uNd1+VcWvGL5f6MH9Iy3znA0SUjPzI0zbb+jn6il+Hzl55blZm+anXDsqdne5h9cP+F62+98\nvPxjw+NidGpm/BrtqPXqwlfx9zl/tyzndFcxPUaC0cnOY0b95cHpD3qqL1tYsG0BrvzkSiwsW4gr\nPr7C9HuPff0Y/rnon6b/v3Pynbb93IoL37/Q8/32k6xXBCEZwupdsS0A1U5cN39xc9z31HaRXrc1\nfHD6gxg5YyQenP4gHprxkMfWEmoEtrBsITZVbXJd/svVX1r+3+wBNZoh6wc1h2ui76/rf130/ROz\nn7As997i9wAAPdv0jB67Z5q3GclvLnrT9jteZ6J6FXjfbfku4dih0CHcPfVuy3KT10z2VJ+fI/Rv\nNn1jX1/Aaa53DboLQ44a4qmsno+WfoRxS8bhk2WfxG2XqueB6Q9YKsVnv3sWL8x7wXM7Jq2ahKlr\np3ou7xdZrwhU509Yv0TnrlDZL8k8LF4m0Bhl3WjbkIr1Z7xsdp4MWkHZp22f6Hu7a50v8gEAp3fN\n7G2svfqivSqQTHANOcGrS9CsX9hlqDXIa4Czis/yVKfbunKNrFcE6qGRkHEdzOxhSsZ/np+X77ms\nFm3bUtEhgw6ieT0fpVgzcYVMLV4Va9CKIOgsHr9dQ3btCMuwbyuQ+tnnkv2tTFBKqd6qMuUooSel\nhPZ+6DupullBWwR2GTypEIJBB2+1isfNRjrZogj8HvnaEfiEsoAtCa+EZdg3oennktbJtikT+n/O\nWQTJjGI8KQIfVwJ1SjotgjhFYHOe6ruZvs580K4hMwFrdz29CuagLRCvSEj/LIIMGIVnEpn9BDog\nqgikjHsQ9A9FNEaQxCjGL4sg1a4hM1dGqjp/3Pm4GN1ELYIMfyg9WwQ+zwewE7ypnDRoWC4dFoFP\no2f1O370PXYNZQBai8CqX6bNNWSy3INqR31wDZlaBPXFNZQpFoHN9Qw8yOyzArETiFL6bxGkey5E\nplBvLIKwDDvq0EkFi4X7YLHZgl6pNKvNBFeqOr1XC0cF392OiIJ+eL0qVr9dLrYWQZa7hpwEi/0a\nPfsaLOYYgT8c+fcjIUaLuL8mjzVB/sP5KHi4AGK0wJS1UwAAt3xxC8bMH4PjXzkepeWl6D+2f/Q3\n9BObxGiBfmP60fvIzXplwSu49YtbcfbbZ2P6+ukY8s4QTFs3LVrm/cXv48pProQYLdDs8WYAgDv+\ncwcA4MfyH6N5/2K0wCmvnxJX3/3/vT/a/vW716PnCz2xdvfa6P9VvvwdX92BgkfIGLN6mErLS03/\nV1tXa/o/M8Flt9mI2YQzMVrgdxN+Z1pOew6dW3S2rCPudyMPgH6UJ0YLw12zFN2e62Z4XN0nu/rc\nsvOA8YYpV35ypWU5q3t777R7Tf9XWVNpeLzbc91Qvq/ctNzQd4caHr/ty9vw8ryXTcuZ9RerPgYA\n75S+k3BMjBZ4Y+EbGL9iPC4ed7FhObO5M8MnDMdZbxunh1bWVOKFeS9AQKBD8w7R4wvLFlq2Ucv0\n9dNR/FwxerzQw1B4f7DkA8PJZWpi4OcrP0+QUwoxWqD3i71xKHQIYrTA/f+9H+/8+A6uHX9t9Lsr\ndq7A+e+dH/2dSasmAQB27N+Bwv8rRIe/d4j77Yr9FdH3/cf2xzcbv4EYLRxNnnRLRiiCHft3JByr\nratFWIajo9vPltMOTa8seAUvznsRi3csxjcbv8HmvSQw9h7ci9WVqxN+Z1nFsrjPMzbMwJgFYzB9\nw3RMWjUJ09ZNw4QVE6L/f3PRmxi3ZByA2ESpF+e9aNhu7TaSAPD4rMej70vLS7F291p8v+17AMDU\nq6di2S3xbQGshcWsTbMAAKcWnZrwv+qD1abl6sJ1KGxWiH337Us4boXVrGurGarqHB47+zEMOzq2\npaLTYLGRgFb31Qj1v6uPvxprbl8TPW52nxR3nGKtKMzo176f4XHVT8xQAvb6E67HhCsmxP3PatvR\n3m1jM+XX/yE2CXD7vu0o32+uCBbvWAwAqPlrTdzxl+e/jOfnPm/eThMLsqq2yrSMllW3rYr7/LdZ\nf8N7i9/DhJUTTEoYM23dtIRdyBRq97KCvAKsu2Nd9PjMDfYz3xWfr/ocG6s2Yt3udYYxgn/++E/L\nyWXjllrf79W7VkdlxuOzHserC1+NW+Fg7pa5+M+a/+D9Je8DAD5a9hEAYP2e9aioqUi4t9pB5I/l\nP0YnRDqZPOmWjFAEXtGbwn5ky/jtdlBtLGxWiAb5DRL+b6UIlKBsXNDY9HeNCMkQ2jVth2YNmyUc\nt8Kr/1WdQ3GrYm/zCAzKOGlLpxad4kaHdni1CLy6QNT1bpTfCC0atvBUX3GrYtf1NmnQJOGY1fU0\nswicupqKjihKLJsi911+Xr7h+TlBew38DNBq+5WXvmJ2b4KM9WW1ItBfdKub4FQIpMrvaXqzLYSz\nmesEsG5nKBxCQV5iHoBdx/IqKNXv6h/+ZILFTpdPDiLjwqtQU9dFCOGqX6UiG8fqOplZisk8C6mK\ngXmJ0ym0fSpVKcuOdvRD/GoIZm0JcoWArFEEcamhDpeVcPObVseSQf2e2c12YhEYdX6rDheSIcMy\ndhaBV6Eal7mlwXaJiUiw2OjaOF23380DHXQwVXu93fxGKoSolWJNNrlA32/cKj4nqHuXzOx+bTtT\nFSy2Om91DvrnxawtrAgcor/olqNrh0IuaIvAkSIw6PxW5erCdYZl7BRlsq4ht1i5hpxu6eimzUGn\nV6r6BFxaBClwq3hxDTlts5Eg8/s5Ur+XzEg+Va4hLY52dpPxCsHsnA6FDvnXMBvqlyKwEHRpcw35\nYBF4cQ15sgiS9KH76RpyumeumzYHvWaQ1jXkRplYjio9KglLRWDSL5J5Fny3rH2YdxOnCFK01pCT\n5WP0z4vZvTkcPuxXE22pV4rADyGeSRaBGrW4FeohGfJkEfjtGrLDyiJwOrJyM0L0qgg8xwhS4Bry\nKmBTaRH4Xdbq97LJIjBzlepdRKaKIBSvCFI59yhrFIFWu5otF+FEqNqRjFluOLJNkUVgJdTNgsWp\n8jmaWQR2WFoEDgSeW9dQuiwCwN21sTp3r22x6v9mgwqnWXj6NkkpU5Z9l4wiSNXELe21tTpvfZDY\nLu6h76+pVAQi3VOshRASo9LahCj9O/THou2LfPmthvkN43x8FXdXoF3Tdgm7YfVu2xsrb1sZd+y0\nN07DiJNGYPiE4QCAewbfk5B3fma3M/H1xq+x9969OOKJI3DbybfhxfMpj/7Rrx/FR8s+wo83/Wi4\n+9agokGYc/2chON2O3UZtVVb7q2L3sLw/sPjfqdnm544udPJmLFhBsb9ahzO7HYmAOCTZZ/gVx/9\nCgDwl9P+gqfmxO9MNvt3s3H31LtRW1eLAR0H4NWFr+KSvpfgf+v/h6qDlN9+y4Bb8PIFLye0+77T\n78Pj5zyOGz+/Ea0at8JTQ57ChBUTcPEHxhOcAOCpc5/C3YNpw5iZG2biqk+vwrbqbSjIK7BUIHJk\n7Pl5es7T+HrT15i0ahJeueAVjJg0AgDwygWvYGDngTjp1ZPiyj557pO4Z9o9mPv7uRjYeWD0+MMz\nH8bIGSPRqnEr7L5nt+F9eennL+HWgbdi36F9aPG3FpAjZfR72vdmbdb+/4pjr7CdEwEA3Vt1x4Y9\nGxAeGUazx5tFc+YPPnAQjR5tZFqXYuyCsQkbRhnRuKAx5v1+Hh6a8RDGrxiPkuKS6PyCN3/5Jq47\n4TrD85t/w3yc/NrJ6Ne+H5bcsgRA7J543ZRIcUGvC/DF6i+S+g3fGQVIKX3RblljEQSBX0oASAz0\ntGvaDgDQsXnHuONqRzUt3275FhNXxTru5f0uB4A4Ifz1xq8BxGYDvzT/pej/QuEQhvagWaaHHzyc\nMNv32y3fWra97sE6vDfsvYTjRm21Y82uNXh/yfso21cWtzXlV2u+ir5vVNAIdQ8mjn7mbJ6DhWUL\n8erCVwEAn634LKoE1HcAYMEN8bvOvbbwteir2tZTzUwHgJsH3Iw/nfqnuDLaLTlnbZoVnVynVQLP\nnPcMAOCq464CADRt0DTuN95d/C4mrpyIsAxHZ92eVXwWRpw0Aid2PDHh2oyZPwYAMH/r/IRzB4BP\nf/0pAODRsx5NKPvhsg8BxE/6Orv72Rh/+XgAdN9X3Gq8BaJ+ZFlzuAY3D7gZhx6wDk6u37M+OopV\nSmDCFRPQML8hpl873bIsALz5Q/xEqJLiEsPv1dbVorS8FONX0LloJ5lpd7NrlB+vfNR1XFqxNHpM\n3ZNkyTgl4DOsCALGi3mnTE+94DEjJEPR7xbkFbj2K+fn5fu2rpLTcnrz2GmwGIChG8yKBnkNEso4\ncR2q7zTMbxj3akXTBk09r6Wk3CCNCoxH2wntg4i77/oJhQq9WzEUDqFxQWPDCY92qGvgpb94cdVo\nXUMtG7eM+1/Qq6HWJ1gRBIwnRaCmwzt8cELh+GBxUO4/P32YTn7La5A62UChup5O5q14yXvXn5dT\nIatPEjArp/f9h2XYtTLV1+FFCHtZCUCrCPTxgqBX3a1PsCIImGSEpeH+xwYPoH5CmZc60z26crOS\nrFGw0gq7dFXTZZJ15awEmfo9T5aVrn4nAVIpZULasJkS0gvMunCd5xm7bhSdm+QOM9ymwaY7Bpot\npFwRCCGGCiFWCCFWCSHuSXV9mY4XAasUgNN8e/2EMi91enmA/Nzy0JVFENDDrs9WswoiK6FkJbjM\nrpf+3J0Iv7AMJ9x3M+Gub7fZBEQnRC0CD/cgGTcpkHhdgpyJW99IqSIQQuQBeAnAzwD0A3ClEKJv\nKuvMdPy2CIwmnejTR4MU6n6RateQH8rDyhXhZw6+kzkAIRlCSMbfdzN3j37kbLYkiROSWfIhWYsg\nwUJj15BnUm0RDASwWkq5UUp5GMA4ABeluM6MJilFYGARGE1D1z/YgVkEPo7MvSzepbCLARgucudi\nHZpojMDCNeTkPpvVo59FayVsVRtC4ZBn11AyMQJVzlOMwIPgtpoUZnQ/MmEbyGwg1YqgMwDtwvJb\nIsdylmTWbzHq1IaKQBcsDipGYFdPnLC1mYTjxiKwihEYrTufdLBYs1iY2X7NTtwUSnDp26M/H6PR\nuqpLbxFo77uZokmwCMLOyhnhxjWUkK3kYW0wqxVE2SLwTtbvWZwNHH/k8dH3F/a+MLoxhaLBI4lp\ne58s/yT6vn3T9gCAJgWJ67BrJymp3wmFQ3h32LvR4+f3Oh+T10xGRU2FZZ1atPnadm1V9G1n7fW7\n77/34cHpDwKIF5T9O/RP+O75/z7f8rcA4JTOtEOcfk+CygOV0XZWH6pGg0caxNV3Usf4iV0AMGnV\npGgZIyHevVX36AY1g4oG4bPln6FJgyZo+EjDqNDSllNzUrTXpHXj1thduzv6We18deuXt+IPX/0h\nelwJtE4tOgEAjm5/dEJ7Zm6ciQaPNIgK4FZPtAIAtGwUS6k0SwcteiZ+/4Bvt3yLa35yDQCga8uu\nGFQ0yHKDFm0fUCmcRnsS6L+rvT4/7fZTDOw80HSHsd98+hvD44XNCqPvT+l8StxuZw/NeCihXo4b\nOCOlM4tukPqAAAAgAElEQVSFEKcCGCWlHBr5fC8AKaV8UvMd+dBDsRtYUlKCM396JvJEXmz1RiEg\nIOJGYhK072++yE8YWeSL/OgexnkiLzriUb+pRhJSyjj3glrfXn0vX+RDbTSfJ/LilteNjshkKPp9\nAYH8vHyEwiFIyLj9BPQjYOWntuqoqi1xG8LLxEClkc9bKwS0U9rtRk1mm+fYlSvIK0g4R7v6hBCG\nLglVn/a89O/zRX6iayBy3c36dJ7IM3SX2N0H1Xf0hMIh2yVCjO696q/afqBHfx/U9QzLsOn5GV1P\nq3LaDdz19y8Upn4dkiHTdhrVF92bwsRtZ7airnr+rO6f/j6oa6nkgdX1VMcL8gqi9QiImCyIfK4L\n10XljfYaKRmijmvPz0guqPaqfqy+q66pVm6pdmj311C/q65ZKBzCjBkzMHNmbOLjIw8/4tvM4lQr\ngnwAKwGcA6AMwDwAV0opl2u+IznFi2EYxh0R5eaLIkipa0hKGRJC3AZgCige8YZWCTAMwzDpJyMW\nnUt3GxiGYbINPy0CnlnMMAyT47AiYBiGyXFYETAMw+Q4rAgYhmFyHFYEDMMwOQ4rAoZhmByHFQHD\nMEyOw4ogIDZuBNavT3crGIZhEuEJZQHRrRuwaROQA6fKMEwA+DmhjBVBQKj1vHLgVBmGCQCeWZzF\n7NwZXF3V1cDBg8HVxzBMdsKKIGCWLXNfRgjgiSfcl+vVC7jySvflamrcl2H8RwhgwoR0t4LJBVgR\nBMCBA0CjRvT+5pvdla2tpddt29zXW14OrFzprsy+fUCzZt4C22PHkvBya4X06we8+KK7MuEwsHt3\ndiit8nLgxx+9lZ0719+2pAK2OrMfVgQBUFYGdOgAXHMNcM457spWVtLriy8CCxa4r1spIKco19V7\n77mv69tv6dWN0vrf/8hKuuMOd3UNGQK0aQOccYa7cgDw9NPA737nvtxnnwEvv+w+zvPrXwP9Ezdh\nc0TY5S6j5eXAP//pra4vvwTGj3dXZtkyoHFj4JJL3Ne3ebP9d/Ts30+Djddfd1+WMYcVQQBUVACF\nhcBpp8VG+E5RigAA3n3X/HtmuBUk1dX0+uCD7spVVgI7dtD77793Xu4D8x0RLfnf/+h14UL3VsGf\n/wy89Zb7OocNA267DZg+3XmZw4eBr792X5dSNk8+af09Pf/4B3Ddde7rA4ALLnAv0FV/dqtAxo4F\nunYlC9QNzZvTqxcXK2MOK4IAqKwE2rYFWrcmd4YbtMHlVq2clwtFdog8fNhdfdoH083I94EHgK++\novduHtKOHem1SeJ2zJZcfHHs/XffuSvrZe/6PXti79242555JvbezfVUCtktSlAG5TJzO9BQqP65\ndq238s8+660cYwwrggDYuTOmCD7+2N0oaNeu2Hs3Qn3hQiAvj4SyG8FcXQ106ULv58xxVw6gALXW\nirFDSqB3b6BpU+dlALKwFI895rzcjh3eUnhHjIi9dyO81L1u2TL+XtqhvYZu+stdd9HrvHnOy6g6\nGjak926E+4EDsfdurmtVFb327x+z7pyg+iYA1PG+9L7BiiAA1q8HiouBBpH9yLdudV62qgo44YTY\ne6cMHBh7oPv1c16uuhoYMICU1vPPOy+nrJWhQ4EXXnBerqaGYidVVe4Eifa7WmFkxy23OP+ultat\nY+83bnRerlEj4K9/BY4+GigtdV7uF7+Ivb/pJuflFHfdBVx0kfPvb98OFBUBLVoAe/c6L6eNBzlV\nWFIC998f++x04CAlxRV++IE+Z0OiQLbAiiAA1qwBevYEjjmGPitfuhOqqoAzzwSuvTbePZEq9u0j\n98LYsTHz3QkNGlAQ9uGH3dVXU0NKpK4OOPFE5+UOHCB//ZgxQLt2zstNmxZ77+b8mjalFN5XXyWr\nzmms55NPqH29etHMcqdorTgvgfuFC4GJE527mCZOJIuzuhp46inn9VxxBXDKKfT+iCOcldm/P/6z\nU4G+PLLbec+elHzhNr7AmMOKIAC2baPRVmEhjdLcuE6qqkhQ/vrX7iajCeFupKyorqZRYZs27uIZ\nVVXk/mjZEigoAA4dclZu/35KVwWARYuc13fgAGUMHXusO5fLVVcBjz5K792MtKuq6Jqcey59dppe\nu2gRudgKC92nAJ9+OgVwAefWUpcu8YHbFSuclbvrrlgWz9/+5swqUBanWyvrv/+l1wsvpFen/ayq\niizd5s2pz+gVCuMdVgQBsHdvbLTUvLm7QKASsEVFztPtDh8mRdCoEY0MATL9naAUQfPm7h40pbCE\ncBcUr6mh0Xa3bs7rAmLXtKgIWLXKebnqaqBTJ3rvJrd/zx46v+7dSRhVVDgv+/DDZBG4qe+EE4Dn\nngMmTaJ74dRds2cPUFIS+7xmjfM677sv9t5JQLysjBTcNddQ9tewYc7qefxxelWKxGlf2b075qJb\nuzY75lhkC6wIAkCrCFq0cGfSKkVw5JHOhc+BAyRchSCBcvTRzkfNW7fGFIHTdobDwKefxs7RjTWh\nFMHMmfGBQDt276Z6OnWi66LmMNjx7rsxa8XNRKjycroHAAVi337bWbkmTei8+vYFtmxxVkZK8oOr\n69m2rTMrsq6OrmeLFrFjTvqMlJRYoHXrORk47NoFtG8fa6PTPjZsGPCnPwFnnUWfH34Y+OIL+3Ja\nRQCQm47xB1YEAeCHRdC2LT0IXbsCGzZYlzlwID4ds2VL54Hml1+mkV6zZuSndhIwVgpDuXiaNHGe\nqbR3L10TpeicZkYpodCoEaWSOnErqd9u2pRG26Wlztxt+/YBs2fHYhF33OFsol4oRMqmWTOgc2fn\nSQKqf/TsSa+NGzuzeqqqqJ/l5QGzZtG8FSeK4MAByhgqKADuvJPOs6zMvpyKJwHurMADB0hZ3Xln\nTKHOmGFfTil/gPqpmyQIxhpWBAGQrEVwxBH0kNbVkXvILjXQSBG4CTSPGBF7wCdNsv/+/v0kCE47\njT4vWgRceqmzutRku8aNqZ1OrZ5du2Kjw5ISYMkS+zKzZtHrb38bm8msRrRWKAF+9NH0OmCAM8Wq\n4h9CkAurstKZgN2zh76v5jt07w6sXu2snMreGjwYuPxyZ/ddLSsCUMB/xAhnFoFWEbixAmtrqX8q\nNyIQb8WYobUICgvJSmP8gRVBiqmro46vHjS3FsHevSQgAbIGAPtZnMlYBE2bkuBRIy+V8mrF/v2x\n7wOUzeN0XsCOHbE5AeXlzsz9cDgWkwCovU5SOvPzKQArhLtJZXv2ACefHCvTsaP7EXPjxrS8iFJG\nVmgD6ACV005Ms2qndtLhEUc4iy3s3x9rJ+Dt/NxaBI0b03uVueVkToBW+bduHUwWXa6QEYqgPq/R\nv28fjXaUEHGbp61cQ0DM/WEnxIwUwZVX2j+o4XAsvpCfD4wa5WxykV5w/e53FBy1o66OroVWiYwe\nbV9u716qr6CAPjsVClrBBcR84napoNp7AJCgdOLmqa6Or699e2fCUsVNFIcP27sDAe+KQH9dOnRw\nbxG0aEHtduLaUxYBEFuDycn907qGWrViReAnGaEI3C6DYEeXLpkzBV3rFgJoFq0TN4ZCK4Rat6Zg\np51wrqmJVwTqYV261LqcGqnlRXrFeec5e9j0iqB1a8qQsVv8bNcueqBVfWPHAiedZF+fPmjoVCgo\npaxQ6ynZWVh6AXvUURT4tQuO6gVsmzbOAqp6RXDZZbHjVuivyxFHuHNhKZxeT225vDznLkitRVBc\nDLz/vjPFoz0/vxTBwYPOU2zrMxmhCNwuxGbHli0U3POT//6XJs+4tV70iuD44ymQ6vR3jEajdg+N\nGtVrP6vfssJIIHgVJAAwZYp1Ob2A/fnPneXa6wWeU7eEXjArbrjBXTubNKEsILulJvSKp1s3Z75+\n/fXs2ZMGN3bxExVPUrhRkNrr4jSOpS/n9D7U1sYUAUAuTyfLmaRCEbz2Wiz2k8tkhCJIxQxBs1Hz\nnj3uZngqzj2XcqXdLp27d2+8MGjdmgRJaam9kD14kHyo2ofGSZBM7xpS7gG7FEQ/RoZATBgVFVmX\n0yu5Ll2o7XbKQOsicNNOI0Xwr38Bv/ylu3YCJNTt4hL6+vr0cbZOkd4iAJxZE/v3x/c1pwFc/f1z\nGsfSn9+RRzpzmen7Z1ERlbM7P61CbtmS+rXXRe8UvJcCkRGKQHW6mhp367gYoUba+uUDqqtpTf8R\nI+ghXrPGXScYPJhe//Mfd4td6UdpAAVg+/enTm1l7isBpI0JqDRSK/QPmsrPdyJItALBqamvL5eX\nR+mZdtdXL2CFoMCo3VrzeotAZebYtdVIEXToYK9Y9RYBQH3Ibnaxvr6uXZ0NQvTXE6DzdeKK0rvo\nnLiijCwCL4qgd29g3Tr7cnqLQPVPu/unve8FBdTOZK0C9WyZKZQxY2iZEEV9jWdmlCL461/JZ+jH\nb+k7yIcfUsrgxx/T51693MURGjUi19BNNznLpFEYjSa1I16rdYeMyqqH26pD6hXBo4/SEgJ2CkQv\ngJo0IYVqt1SFXgABzgKVRgL2Jz+JLUFghl4RCEEuN7tF3cwUgZ2rzeg+HHWUvSLQB4vVyNdujSMz\ni8Du/unPr2VLuqd29envnxps2A149PU5cVsCif1TCLrvdhayvr+0a+duhrcRqo+aKcw5c+JnMF9/\nPfCznyVXpxsWLXK3NplXMkYRdO0a25/1pptoESwnVFTEC0U1utObqCogqf2um+ydXbvid5k6/XTj\nHcPmzQPmz4991qY5Kr75htbu79XL2i2mTR1VNG5MGT1WlkRNTfyDXVBA19etRSAEuTPsgsxGI1gn\ngUojAXv11RRDsUog0KYRKn7yE/tJZXrBDJAisEuVNFJYrVvbn58+RtC4MZWzE5b6+6fqc6LIteen\nArhO7ru2XNOmdP/shKxRtpGTtFO9RQDQ9bU6PzVrWp+F5Wb9LSNUebNz3bEjXpZMm0axr1CIZMnW\nrald8+iEE2Iz2lNJRiiC+fPJ965S5P7xj9huXOPHA1OnmpctLKQOrwR8eTlw3HF0g7RCX91w7TE7\nf7+2c1RUxCuC2bNj1oWWU06htWgURsLu9NNpVNG6tbUJblQWsPcX79+fOKJ0Kkj0AujMM+13ETNT\nBHaK1uj8iovpflrFM/QWAUAru9plfxhZLm3b0j2wcmMZKXMnis5I8ZSX2+8CZnT/QiHgxhutLUGj\n8+ve3VlQW99OJ6NtI0Xw8svuY1gA0KOHdSBd9ZU8jcTywyJQgzmz3ykvj1cESoEVFFBbiorISjDi\nhReAhx6i91bPXjgc29QJIKV33XX+Z1NakRGK4J57Yu/VgmDKZ3f55ZTGeOhQ7JiUNItV62JR7osd\nO8hslzK+Q27bFlsW4G9/o1RFqwdk6VJSMsuXk9Dduze2P66aAKXvPEqYKJ/nAw8A77xjvkyynS/W\nKL4A2CuCZEaU+nL9+sUsNatyekHiVREA9uscGSmCwkL70aFRfXl5VNbK/N6zJ7GcClZasWNH4szl\nP/0pNjHQDCPX0OWX06uVNWF0H447zt5SMrrv7dvbC1mjiWiAvTVvpOg6dLC/B/p73qaN+xjBggXA\nrbfS+1CIrs3PfkbX1ehZ3LEjtkZUZaWxstIPlI47jvrvww8DjzxCkwjbtCE30513JpZfs4Yy5rZu\npZjE7bdT6vXcuSQjmjWzd5cmS0YoAgA4/3x6VUFZ1SlUh2nUKLZs8IoVtMjYv/8dK68eyhdfpAf7\n4MFYxwRIEahJTjffTBt/WCkCpUSOOYa2QvzJT6gtdXWxti1eHF/mz3+m17IyUlyPPUZzBsxW1nSi\nCLxYBEaCxEng0Egg/PSn5IpyWy4ZRWB3XYwUgRPBZVZfhw7WmUoVFYkC3cls7bKy+D4IkLWYZ/PU\nGV3PoUNpyW2rczSLgdx8s70loS/nxO2it0DUgMduvkNlZeLgqLDQ+tyM3HN27iQjTj6ZhO3SpcDZ\nZ9N5dutGilY/6AqH6Roo78LTT9NxNSDs2DE26FODwLo6euYXL45Zter7gwdTXHLBAnK7qhiMWiF2\n0SJSUmPH0ufvvydLqWtXylpUiTSpCFhnjCIYN44utLIIvvkGGDky3oRcsoQm16gNXiZNog7Uqxc9\nlLW1tHiVELGUwF/9il63biXL44or6CHu2JE6stmIYssWWlMGoIlHqsPn55OVsHQptUdrvtXU0MzY\nHj3ilxw2GwHa5WsbjUSBYF1DnTrZp3Oqpau1JGsR2CkCbfoo4EwRGMVcADpHK9+2duVRhZPzM1IE\nTu6DkSIH7Od1GLmG1KQ57Y5gerxaBHoFcvTRtLyIVbnaWhok6fuLXX1misCNRaDNSPzmG+Drr0kx\na5W89rmqrIyt4ltTQ+1+9NHYM19WRu7lI44g186TT8bKq+03jRYnvPpqen31VXKLT59On7Wxxbw8\n4I9/pH6uJlmWldGxHj3I0vCTlCkCIcRIIcQWIcTCyN9Qs+/++c+x1Qjvu4+CMd27k2lVVka+tksu\noSCi1i8/c2ZMsFdVxUbqffoAn31G71Xq17ZtwKBBNIsRoAt9+uk0KlA7HwH0kC5fTmlwQ4fSUrnr\n1sVuHkCTiY45hkYDb70FfP45Ha+ooOyV3r3j/YpWisBK4O3YYRwosvONGrmGnGad6Mu1akUPgFVA\nTD9XAvCeNQR4swiKiijGZOVXNVM8Voqgtja2qJ4WJ66hsjIakWtxch+MFLmq00rwGbmGmjalZ8Qq\nA8jo/vXsaZwMocXIkjjmGOu+uXMn9V/9Mil2y3b4oQg2bKCYV2FhTKF26xbr8yeeGL8Pw5dfUj9W\nA68lSyhOeNNNJA9696bvHX88rYB7770xK+Dxx8ndc/Ag/cZzzwHDh5MgV66mUIgU0d//TrJm9GiS\nKePG0SAYoPv3zjvkMfnwQ1oNeP36mBzzi1RbBM9IKU+M/H1l9qX/+7/Y+yOPBIYMoYuhuP12Wu9e\n3bzKythMxKOOogfk73+nUXr//nSxteb3oUOkCPSjs8suo3Xf1UUHgD/8gTrz6NEkXGpqqLMddVRi\nuwcPpnkJv/wltU2NADt3jvkS77qLPhthJ/CMRqKAfe67mWtozx5rgWA0MhTCfhGy6upEs7pZMxKi\nXgSzndvFSBG0bk0PrNWa/15cbUoZ6wWXndUiJfmdjSwCJ/MBjFbjtBN8RsFpgNwKVsLZ6LoMGUK+\nbSs3hFF9diP7igrjmNmxx9IoXe0eZ1Subdv4Y24VgVJCHTuSrBg2jAT2iSeS6/eYY+KXT9+wgYS3\numcrVtAg8NpraSLi8cfT91TMEIgpgpqamLv76KNJrrz1FsmnffuoHdp4z8030+v+/eSmUkFmNdho\n1y5+wULt4NUPUq0IXKzxGM/AgXRzhmrsiMaNSai2aQOceiod69uXOvFHH5FvrX//mBL4wx/oddUq\ncjHpheN559HrRx/FJsJoO/6558bW4z/uuMQ2ahdWmzgx9uAfcwxlOnXoQArKbJE4O2Gyc6fxMsl2\nue9GiqBRIxr5Wi1eZqQIAHtFYDSizMuz31DFqyIwSh8FyIq0Oj+z+qzcNWbKWJn8Zlty7t1LbkS9\noHRiERi52gD762JWzsn901+XY44hIWtWTkrj63nkkWSNmykDsz7dvj0tXW62L4GRUu3UyT61WV93\nu3bUxtmzKTWzcWMKFi9aRPEwFZStqaFnul8/umdqCXHtbPm//IVikmqw06ABnfeZZ9LnSy6huUt/\n/GOsjNpDYfDgeAtIbUmq5BpA3pE//YneN21KrqO77qLPyc630pNqRXCrEGKREOJ1IYTBI2jNokU0\nk1cxbVpsJyMhSCCcd16s80+cGD9yf+45uvClpcajcu1Mzx49aEKb1qfXowelg0ppPNJSiuCaa2gk\nUVFBnaxPH7I07DZVt4sR7NqVOAoCnMUIjAR6ly7WKbNm5Tp1sh5pG1kEAJ2/VcDRzDVk5QsPh0lw\nGZUrLjaf5HXwIJXV568D1tdz+3bzPG4rRa5dXltLq1Yxd5MZZgK9dWtrxWpWrndv642CzLKp+vQx\n3+qypoYEn94HrlKszzmHBkF6zCwCgAZuZpMXjdxsJ59MAdXHHjMuY1R327aUJLJuXeLvdesWs7Rv\nu422eR00iPrHypXUx7Xne/LJ9L0xY0hWtWpF17lHD5IZnTvTQPKKK2JlBgwgF9OIEbF055deokHD\njh3xLp/HHyeXkyoH0DX94QdnW4m6ISlFIISYKoQo1fwtjrxeCGAMgB5Syv4AtgMwXVF91KhR0b8Z\nFlsVDR5MVoKidWtSCNqZfr/4RXyZtm1JEaggtJ4uXUirA5QdpPz9TjjjDFpy+fzzKYWsSxfqKEVF\nJAD1Ixg9dq4hs5Gvl6whwLsiOOUUCqyZYWQRAPaKwEwoWI18q6pIABtlMhUXm1sEKvBuZJ3ZWQR6\ngaGwun9qwx09eXmkCJQ1aoSZQLda3ygcNo4RANQfy8vN3TxmllLPnuaKwCq+89JLlDVz992J/1ej\nciOsrCWzwHtREaVpK1eKFdu3071UiRwnnxz/f21/LS+nQWTXrnTdZ80ylyFFRSSXOnaktthNAOvX\njzwMc+eSFavSWdu3Nx5QATRXYdq0GRg1ahTGjx+Fxx8fZX/CLihIprCUcojDr74GwFTEjho1Kplm\n4Mor6W/TpsTAbLt2dOPNbiJAWr22NtZxL7zQflkFgG74G2+Qibd+fSxDqUcPejUSBFqcKAJ9dgxg\nn3liFmy0UwRGE6cAOh+rDVX0K6wqrPzFBw/SNTdzDZm1c+dOYysJoAd25kzj/1VWmu9GZnU9rSwC\nK4vOaA6B4oYbrNfkMVME3bvHp0xr2b8/NutcT+PGtBWlkQvIzMUDkMVrNsnLTBEAsXRJI6VrlIqr\nsIqfGLmGAFL8xxxjPVBRlJWR+6dZM+NnRJsyqwKyanb9hx/GD0KNUH192DD7tigXk12qrZZzzinB\nOeeURD+PdrJ5h0NSmTWkHUcNA+BiFX5vGGXnFBdTapedT03NXxg2jGYzT57svF7ldjrxRHpV5mOf\nPtblrBSBlNYWgZUiMAqmAs4UgZFAsPL1W41ErSwCZQ0YCQsr19D69caBe4AUgZlFYDUStRJAVhaB\nlWvIzCIAYgMGM6wUgZnry6yMwiyutHs3KQr9TF/A2iLYvdtcEfTrR+6Qvn0T/7dpUyz3Xo9SyEaW\ni5FrCCDF9/HH9qmuoRBNjFT7QBsNlFQ23h//SMFYtTx1nz5Uv9l8IMXvfkfpunpLw4ziYvvJhUGR\nlEVgw1NCiP4AwgA2ABiRwrpMOeMM2ubvN7+x/p4S2lddZT/hx4j16+NHLOvWOXMNmY0o9+0jH6zR\nA2qXeWKlCL780ryc2bwFK0Wwfz+10WgkajUpyWpkaOUaWrMmZnHpOfbYxEl+Cq8uifJySjM2wqtF\nYDU3IxxOXFNH0bUrCaTDhxMXPjSzyhSdO5Plqh+cbN5sLph79TJXBFu2mGfD5eeTf/ussxL/t2wZ\n8PvfG5dr0oSePf1+GuGwtUI+8kj7hdm++IKUgVHSh6JxY7q2zz9PcqNhQzqurpnZAERhlvFkxty5\nxvMM0kHKFIGU8ppU/bYbLr6YOpLd9o5t2tCNUaN6t+gtju7d7cvYjSitBGV1NZ2XXmkdOkRuFyNB\n0revdZaFmWvIamRvNRItLDQPalmNmK1m+q5ZExvVGdVXW2uc327lUkrGNWR2/8rLzRWWimUY3T+l\nWI0GIw0a0LXZvDlRKJll4yiUItCzeLG5IlAWgZSJz8/69dZ93GjOQ2UlZfCdcop5OTXI0SqCsjJ6\nPs2EZtu2sdToAhOJdtFFlJZp9n/FE0/Qc6DSOQFS3IMHk8vYT+xcx0GSMTOLU4nTjcoHDrTvKH5i\nNfK1yq7Iz6cHw2gUpKwBo3Pu3p0ErNlcAqs8ezOT3cyKAEjAmK2/b6XolAtk6tTESVtr15orAiHM\n3V9WFkHz5qQ8jVJBrUaiVhO8VGDSrL7WrY0zsexG9mZbnVpZIICxIqisjJ8YpadlS1Kqn38OvPde\n/P/sFMERR8QGK4qlS2lErkbaRhhNliwttd5FLD/fPqOqSxea+WvHPfeQNaNVjkJQjKw+72SWE4og\nU7EaaduN8Mz8xWZuIYBGlIWFxqPtUIhcEkaj+wYN6LhRW60ErFWWi5UiaNOGBMh55yVuUmPlGlJ1\nGsUJrNophLFVIKW1C8QqGG5lSQDmgVgrywWg3Hcjq84sXVXRsyfFvrR7EyhftlpDx4jTT6fR9G9/\nG2/92CmCggIKymrLbN1qfi0VRv165kygpMS6XGEhXXOjZIFDh0ih29Wdy7AiSCPNmsWCrXqsBCVA\nroFVqxKPmwWYFWaj9L17aaRqFh/p1cu4PqMFxBRdu5orAqsRrBAxF5XWIgiFKPZipQh69jQXsFbz\nOowUQVkZXRMr15eZb3rzZutMNbONbaxcZoB5iqydIvj5z8n1qZZcmTOH6l+0yHpRQe1cAK1VYKcI\nALqH2liW0ex+PT16JMYlSktJAVrRuTNl+ZSUJCrnLVvoXgRp7WcbrAjSiBDmo0qzLAnFGWcYrzdi\n5yIoLjZOXbQbifbpY+zvtwvCNm1qvE+AnaJ7+mmKaSihLgSl/rVubTzXQWEW4LQ7P6MA/LJl5m4T\nwFwRHDhA988quNijh/Hqt3bXxcwStLvvXbrQUiqqz0ybRrnpdimRAwfS4mgPPURLJACkkDdvts+i\n0e9YZpYCqqVXr/jlEz74gBSYdi8QI44+mlI8gcTJc+vX27c112FFkGbMFMG6ddaC5OqrKcVVP99h\n61bz4B9AmTVGrgWzpRQUP/lJ/JZ9CjuXy6WXAq+8kvg/O4F32WU0ev33v2MZX7Nnx1ZiNKNbN2OL\nZ+tW6xG6kUUwZIhx5ovCTBHMmUMKxGoE2qcPBSa1s04B++vSty8FePXxmo0b7YXd739P7qH9+0kp\nXHut9fcBuoc33EBr7iiX4tat1EajWdpa9NlRVnETxfnn01r869dT/77iClLQVn0aoH6xfj1ZNwsX\nxiV4r9EAAA1wSURBVP9v1Sprhc6wIkg7ZopgzRprRdCiBQkFlY6ngs5btsSvh6LnuOOMUyztFMFp\npxlvcGLncjn/fOPlDex86ABNFBo6lJSBcg3oZ47rMdsgfvNm65xtfQqpCgJrFyTUY6YIzj2XXDFW\nnHMOvX7wQbxr0G5k3707BXGnTYs/bpVNpSgqomswbRqNoNXELyd060Z967bbyJJxkhWnVwRGs4ON\n6jn7bFo64r77qI9/8YV9wscll9DilM8+SwMGLcuX28/pyXVYEaQZI2ESDtN6InbmcJ8+NJnmjTdo\najtgrwh69jR2DdkpgqIi4ywXOxeWWVDUKn9dy0cf0evll1Mwe4TNbBSjrKF9+yit1Mo11L59vBtj\n6VJa38XKf15YmLgKbE0NjZTtcspVlsupp8YLrq1b7YXlJZfEW2cHDpDys4qdKPr2JeXjNgMmLw8Y\nNYq2ovz+e7Is7fBiEQDkjlq+nM7rww9jq3ha0bw5LVc/ZEj8XiA1NbRy8aBB9r+Ry7AiSDMdOyZm\n8bz1Fgk0q5EhQJ37iSdoRK5GwWvXWs+iVsFivWvBboTeoQNZLvrU002brF0SxcVUThv0PXSIjlm5\nahTNm9Po/C9/MZ5cp6d9e8pU0U7d37SJFJnVqLJnz5jPftcuWjXypz+1rqtDB/qu1j23dCkpaP2E\nLyPatKFZuFqf/8aN9qPtHj3iV8t94gn6HSfXR8WW7LJwjBg5kizKsWOdzbfp1CmWsiolnZsT5X/c\ncbSxy5Yt7l063bvTQED109deI4vLau4Cw4og7RgF/159FbjlFvuyV1xBQuHJJymQKyXlmKulbo1o\n0YLyuPWBUaN1mrQ0aJC4hrqTcgUF5FZ68MGYcN62jYSo0ywOs8XijMjLI+WjDcQuW2Y/AtZaLl99\nRX5m7aQiIwoKSAmq+7d2La1gq5YrcYL+/jvx9R91VEwRbNpEGzipJdft+P3vaWOV3/7WeRu1XHop\nnaeTyVXdu8fuQ3k5TQizymhTDBlCy1EPHux+5m2jRmSpKet19WpSYHbbreY6rAjSTPfuia6a8nLr\n1SkV3bqRb/g3v6GRV2kpjaCtXCAACW69+8QuOA0kuofq6kgx2I3sb7yRzPbx4+nzmjXOfMxe6d+f\nXGuffkr+5dJS++wYrSJYvZr8005cLdr1eB5+mPzvZksoGKFVBIcP2wf7AbpPSsCuWEHxBrslVBQd\nOlDw3olANuKBB2KK3I5+/chCkpKuqXb/DivatydF8Npr3tqovT7aNYMYc1gRpBl9PvnBgyRcnZjQ\nikaNyPS95x5nJn+XLon5/VaLuSm0IzyAhFZhob0b5LLLaM9cNQ9h6VJnPmavDBhAS4pfeikFl2fP\ntl8IrFs3ChZXVtIGIGr9dzu08yuqq8n/bpfzri+v0nJXraJ22GXjdOlCdVVUkD/cygL0m/x8+xiG\non17sj7VIMVNO3/6U3fPgBbtUiqsCJzBiiDNdOtGAdfaWvq8YgWNRK2m4Rtx442Ubnfjjfbf1c8J\nqK0loWIVZAYSt/Jbvtx4hUkjtAvCLV6cWuF10UW0laBi/nzKRLEiP59iLjNnUiB24EBndQ0YEAvc\n/vijewV33HE0Wq6pofJq+0O7tp56KqWpzp5NrrdM5cQTgXnzgG+/DS5gO2gQXZtNm8hq9apQcglW\nBGmmYUMSzKWl9Lm01Jkw0PPrX9Mo0YlFcOyx8evVrFhB1oCdH1WvCBYvtl7NUYsapUtJe9OmUnj1\n6kWunVdeoaWHP/7YehKa4pJLaJJVcbG9UlSccQadz4YNdP2dKkZF48akFBcsoDX1nQY1S0poX+9v\nvrEPaqeToUPp+k+ZYq+M/eL000mhT5lC18ZpfCmX4UnXGcCAAZSSN3Ag+ba9KALAeMVRI449lnaR\nUixcaD9RS5XTpuaVllpPuNLSowdltYwbRy4Yr+folPvvd19m+HAanV95pfMyXbuSC+Saa2jfWS9L\nmA8dStdl0iQSXk5QCuvyy53569PFb39L9+K881IbF9LSowcNbG65JZZ+zNggpUzrHzUht3nzTSl/\n9St6f/zxUs6Zk9r6amqkbN5cyupq+nzLLVI+84x9uVBIynbtpNy0iT537y7lkiXO633+eSkBKUeN\nct/mTGbyZCkHDZJyzRpv5TdtkrJhQyl/8Qt35Wpr6Z5kOjU1Uh4+HGydGzZI+f77UobDwdYbJBHZ\n6YscFtJsI9OAEELIdLch3ezcSaOYyZNp/4QtW1K/QFZJCXDvvTQa7duXFhRzYhVcdhmVOfNM+tu2\nzbnpLSVZPiecwOl8eqqraV0mvi6MU4QQkFL64vhi11AG0K4dmfinn06Tg4JYJfGCC8gd0bkzCSGn\nmS5XX01pkj/8QLu5ufG/CuE8GyfXsNpmkmFSDVsEGcLhwxR8PeGEYIJblZUUpC4qoolp997rrFwo\nRFk5S5ZQ8DeT/dMMU5/x0yJgRZDDzJlDwvz2250ticAwTObAioBhGCbH8VMR8DwChmGYHIcVAcMw\nTI7DioBhGCbHYUXAMAyT47AiYBiGyXFYETAMw+Q4rAgYhmFyHFYEDMMwOQ4rAoZhmByHFQHDMEyO\nw4qAYRgmx2FFwDAMk+OwImAYhslxWBEwDMPkOKwIGIZhcpykFIEQ4ldCiCVCiJAQ4kTd/+4TQqwW\nQiwXQpyXXDMZhmGYVJGsRbAYwCUAZmoPCiGOBvBrAEcD+DmAMUIEsQFj5jFjxox0NyGl8PllL/X5\n3ID6f35+kpQikFKulFKuBqAX8hcBGCelrJNSbgCwGsDAZOrKVup7Z+Tzy17q87kB9f/8/CRVMYLO\nADZrPm+NHGMYhmEyjAK7LwghpgI4UnsIgARwv5Ty81Q1jGEYhgkGXzavF0JMB3CXlHJh5PO9AKSU\n8snI568AjJRSzjUoyzvXMwzDeMCvzettLQIXaBs0EcB7QohnQS6hngDmGRXy60QYhmEYbySbPnqx\nEGIzgFMBTBJC/AcApJTLAHwIYBmALwHcIv0wPRiGYRjf8cU1xDAMw2QvaZ1ZLIQYKoRYIYRYJYS4\nJ51t8YoQYoMQ4kchxA9CiHmRY62FEFOEECuFEJOFEC01338hMtFukRCif/pabowQ4g0hRLkQolRz\nzPX5CCGujdzXlUKIa4I+DzNMzm+kEGKLEGJh5G+o5n+GEyMzte8KIYqEEP8TQiwVQiwWQtwROZ71\n99Dg3G6PHK8X908I0UgIMTciSxYLIUZGjhcLIb6LtPV9IURB5HhDIcS4yPl9K4ToqvktdxN6pZRp\n+QMpoTUAugFoAGARgL7pak8S57EOQGvdsScB/CXy/h4AT0Te/xzAF5H3pwD4Lt3tNzif0wH0B1Dq\n9XwAtAawFkBLAK3U+3Sfm8X5jQRwp8F3jwbwAyiWVhzpryKT+y6ADgD6R943B7ASQN/6cA8tzq0+\n3b+mkdd8AN9F7skHAC6LHH8FwIjI+5sBjIm8vxw0dwsAjjE6b6t602kRDASwWkq5UUp5GMA40ES0\nbEN1LC0XAXg78v5txM7rIgD/AgBJGVQthRBHIoOQUs4CsFt32O35/AzAFClllZRyD4ApAIYiAzA5\nPyBxUiRgPjEyY/uulHK7lHJR5P0+AMsBFKEe3EOTc1Pzk+rL/auJvG0EEuQSwFkAPokcfxvAxZH3\n2nv6MYCzI+9/CZcTetOpCPSTzrYgOyedSQCThRDzhRC/jxw7UkpZDlDnRWweRrZOtCt0eD7qHmbj\ned4acY28rnGbmJ1HVvRdIUQxyPr5Ds77ZFbcQ825qZT0enH/hBB5QogfAGwHMBVkie2RUoYjX9G2\nNXoeUsoQgCohRBt4uHe8+mjyDJZSDgBwPqgzngFSDlrqW0Te7HyyNRV4DIAeUsr+oAfw6TS3J2mE\nEM1Bo8Q/REbPTvtkxt9Dg3OrN/dPShmWUp4AsuIGglxfTvF879KpCLYC6Kr5XBQ5llVIKcsirxUA\nxoNuXrly+QghOgDYEfn6VgBdNMWz5Zzdnk9W3VspZYWMOFcBvIaYGZ2V5xcJJn4M4B0p5YTI4Xpx\nD43Orb7dPwCQUu4FMAPAIACthBBKVmvbGj0/IUQ+gCOklLvgQc6kUxHMB9BTCNFNCNEQwBWgiWhZ\ngxCiaWR0AiFEMwDngVZknQhgeORrwwGoh3EigGsi3z8VZPKVB9hkpwgkThAcHnk/HPbnMxnAECFE\nSyFEawBDIscyhbjziwhGxTAASyLvJwK4IpKd0R2xiZGZ3nffBLBMSvm85lh9uYcJ51Zf7p8Qop1y\nawkhmoCu+TIA0wFcFvnatYi/d9dG3l8G4H+a40bnbU6aI+RDQZH/1QDuTWdbPLa/Oyjj4AeQArg3\ncrwNgGmRc5sCoJWmzEugKP6PAE5M9zkYnNO/AWwDcBDAJgDXgTJIXJ0PSNisBrAKwDXpPi+b8/sX\ngNLIvRwP8qer798XOb/lAM7L9L4LYDCAkKZfLoy01XWfzLR7aHFu9eL+ATguck6LIudzf+R4d1As\nZBUog6hB5Hgj0MTd1aA4ULHdeZv98YQyhmGYHIeDxQzDMDkOKwKGYZgchxUBwzBMjsOKgGEYJsdh\nRcAwDJPjsCJgGIbJcVgRMAzD5DisCBiGYXKc/weMFAUIAPQrhwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f1e2d41ad10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "[<matplotlib.lines.Line2D at 0x7f1e2d21ef10>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEACAYAAABRQBpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYFcXV/79ngIFBQUBQkB1MEHBBg+Iuvi5g1OCamOgr\n8Y3RX+LyxjcxxhiDRPNGTdRojIl5o4ka4564xKjEhWwqQgQ1gLKIIKsg+w4z9fvj3PL27am+XVXd\n07fvnfN5nnnu1jVVvZ1vnXOqqkkpBUEQBKF1U1fpBgiCIAiVR8RAEARBEDEQBEEQRAwEQRAEiBgI\ngiAIEDEQBEEQYCkGRDSWiN4lojlEdJXh9yuIaCYRzSCivxBR38Bv4wvl3iOi89NsvCAIgpAOFDfP\ngIjqAMwBcByApQCmAjhHKfVuYJtjAExRSm0lov8HYLRS6hwi6gpgGoCDABCAfwE4SCm1rkX2RhAE\nQfDCxjM4BMBcpdRCpdQOAA8DGBfcQCn1V6XU1sLH1wH0LrwfA2CSUmqdUmotgEkAxqbTdEEQBCEt\nbMSgN4APA58Xo2jsTXwFwHMRZZfElBUEQRAqQNs0/xkRnQfgMwCOSfP/CoIgCC2LjRgsAdAv8LlP\n4bsSiOh4AFcDOLoQTtJlR4fKvmIoKwskCYIgeKCUorT+Udk/AG0AzAPQH0A9gBkAhoa2ObCwzeDQ\n910BzAewW+B9F0MdqpaZMGFCpZvQosj+VTe1vH+1vG9KKVWwnbF23OYv1jNQSjUS0aXg5G8dgHuU\nUrOJaCKAqUqpPwG4GcAuAB4jIgKwUCl1mlJqDRFdDx5RpABMVJxIFgRBEHKEVc5AKfU8gCGh7yYE\n3p9QpuxvAfzWr3m1xbPPAqecAsiq4YIg5A2ZgZwBo0ePBgC8/XZl29FS6P2rVWT/qpda3re0iZ10\nlkkjiFQe2tHS3HgjcPXV4hkIgpAORJRaAlk8gwypk6MtCEJOEfOUIZTOADBBEITUETHIEPEMBEHI\nK2KeMkQ8A0EQ8oqIQYaIZyAIQl4R85QhIgaCIOQVMU8ZImEiQRDyiohBhohnIAhCXhHzlCEiBoIg\n5BUxTxkiYiAIQl4R85QhkjMQBCGviBhkiHgGgiDkFTFPGSKegSAIeUXEIEPEMxAEIa+IecoQ8QwE\nQcgrIgYZIp6BIAh5RcxThohnIAhCXhExyBDxDARByCtinjJEPANBEPKKiEGGiGcgCEJeEfOUISIG\ngiDkFTFPGbJ0aaVbIAiCYIaUUpVuA4hI5aEdLY3OGbSCXRUEIQOICEqpVLKR4hkIgiAIIgaCIAiC\niIEgCIIAEQNBEAQBIgYVQRLIgiDkDRGDCtDUVOkWCIIglCJiUAFEDARByBsiBhWgsbHSLRAEQShF\nxKACiGcgCELeEDGoACIGgiDkDRGDCiBiIAhC3hAxqAAiBoIg5A0RgwogYiAIQt4QMagAIgaCIOQN\nEYMKIGIgCELeEDGoADt3VroFgiAIpYgYVIDVqyvdAkEQhFJEDCrA2rWVboEgCEIpIgYVYNOmSrdA\nEAShFBGDCiBiIAhC3hAxqAA7dlS6BYIgCKWIGFQAEQNBEPKGlRgQ0VgiepeI5hDRVYbfjyKifxHR\nDiI6I/RbIxG9SUTTiejJtBpezYgYCIKQN9rGbUBEdQDuBHAcgKUAphLRU0qpdwObLQQwHsC3DP9i\nk1LqoDQaWyuIGAiCkDdixQDAIQDmKqUWAgARPQxgHIBPxEAptajwm+npvpRCO2sKEQNBEPKGTZio\nN4APA58XF76zpT0RvUFErxLROKfW1SgyA1kQhLxh4xkkpb9SahkRDQTwMhG9rZRaEN7ouuuu++T9\n6NGjMXr06AyaVhmUyX8SBEGIYfLkyZg8eXKL/G8bMVgCoF/gc5/Cd1YopZYVXhcQ0WQABwIoKwa1\njoiBIAg+hDvKEydOTO1/24SJpgLYm4j6E1E9gHMAPF1m+09yBETUpVAGRNQdwOEAZiVob00gYiAI\nQt6IFQOlVCOASwFMAjATwMNKqdlENJGITgEAIhpJRB8COAvAL4nonULxoQCmEdF0AC8B+FFoFFKr\nRMRAEIS8YZUzUEo9D2BI6LsJgffTAPQ1lHsNwP4J21hzZCkGmzYBK1cCAwZkV6cgCNWHzECuAFmK\nwRVXAAMHZlefIAjViYhBBchSDD76KLu6BEGoXkQMKoCPGFx1FbDPPu7lGhvdywj54brrgBtuqHQr\nhNaAiEGG9C1kVXzE4MUXgffecy8nYlDdTJzIf4LQ0mQx6UwoMGoU0K6dnxj4hpaamvzKCflBRp8J\nWSCeQYY0NQFt2mQrBiQrQwmCYIGIQYYkEQPfHn6bNn7lhPwgnoGQBSIGGVIJz6BOzrAgCBaIqciQ\nJGIwf3767REEQdCIGGRIEjHYvDn99gjVgYSJhCwQMciQJGIgCILQkogYZIiIgeCDXC9CFogYZEgl\nxEAMiSAINogYZIh4BoIg5BURgwxpagI2bgRmz650SwRBEEoRMciQpiZg+nTg8ccr3ZJ4vvIVYMqU\nSrdCEISskLWJMqQS6wT5hqTuvZfXURo1Kt32CIKQT8QzyJCkYtC/fzrtsGXjxmzrEwShcogYZEhS\nMch60bkdO9zLbN8O3Hln+m0RBKFlETHIkGoKE/mWfecd4LLL3Ms+9xzw1FPu9QmCkA4iBhlSbZ5B\nEiG5/3637c84AzjtNP/6BEFIhohBhlTbg2Z8xEAL1uLF6bZFEISWRcQgQ6rNM/Bpr25jtQmfILR2\nZGhphjQ18fNs//nP7OpMEupJIgauz16WJ7IJQmXJjWfQGpZoaGoC2rf3f/pYNeUM2ko3QxCqityI\nQWtAr03kG0KpBjHQbWzXLt22VDsffQQceGClWyEI0YgYZEhTE/eYq8ULylIMtmxxr6uamDULmDGj\n0q0QhGhyIwbVYiCT0Bo8A414BqX4hgYFIStyIwatgWrzDHxES+9bVjmDRx8FDjggm7qSIGIg5B0R\ngwyphGeQ9QxkXaa+3r9eW3buBH7yE+Dtt1u+rqTU1fid9otfANu2VboVQhJyc4lWS285CUk9g2oI\nE2mhy6In/OCDwNSpLV9PGtS6Z/D1r2d7Lnbu5HWwhPTIjRi0BpJ6BlmTRAyyEPdqSjrXuhgA2V7X\n550HDByYXX2tARGDDGlNnkEWYlBNE9VEDNLljTeApUuzq681kBsxaC1hoqw9g6xnIOsyLmWD7n5W\n18Hy5cBhh2VTF1D7OQMAuOmm7M6f6wx3IZ5WcInmh2rzDF56yb2Mj2ewfn3z8jYkOR4zZgCvv+5f\n3pVqEgPf6/P557Mz0iIG6VNFl2j1k9QzmD0buPXWdNuUNj5iEDSULuWSiEHWeRvd1rx7wM88k0y4\nstq/1auzqac1kRsxyPtNkgZpzDP49a/Ta08cHTu6l9E9Nl8xyHtyvbEReOst93L6eGSxf0oB3br5\nlZ03L1ndWZ2/aho8UC3kRgxaA2nkDGbPTq89cfTq5V4maQI5K8/At32PPAKMGOFfXxbhjaYmYM2a\nZMuJ+NIaOnW1iohBhmgxyPKGSVKXj2HwEYOgOGaVM/AV5LVr/cpl7RlkVVeYvHt2QjS5EYPW0KNQ\nisNE1XLDZCUGwW19y7niW/aSS5LVl8W513X4eCHV5Bnsumt2dbUGciMGrYFqG1qad88gyXG8917/\nsoD7M57vuINfswgTJRGepGLgWuf99wOLFiWrU0gHEYMMSSOBnOUQxbx7BknE4Kmn/MrpZxL85jdu\n5e65h1+z9AwqIQau1/b48e4j5CZN4tdqmnRYDeRGDFpDmCgNzyDL41RNnsGLL9qXS8JBB/Gr73nI\n0jOoRJjI59p2XdTwd7/j19ZgM7IkN2LQGlAq+wRyEnySpZXyDBYutC8H+Bu9pLF/8QyK6LkCrmIg\nHkHLIGKQIdW2UN3y5e5lknoGvuVcDYTv8xaSGFog/55BUnw8O9cHIelzLaKQLrkRg2rpLSchjaGl\nWR8nV6PnYyyD++QbJnI9Lscd57Z9uB7f85B3zyBpTspnsqGrGFTT0h7VhBzWDKm20UQAsHmz2/Y+\nM5B9jXqSfTvmGL9y1RAmqpacgW6n64quraHjWAmsxICIxhLRu0Q0h4iuMvx+FBH9i4h2ENEZod/G\nF8q9R0Tnp9XwakTPM6imi9lVDJLmDFyMSdeuxfdZhQySzrDOagZy8DVLfDoBrudOLwciYaJ0iY2c\nElEdgDsBHAdgKYCpRPSUUurdwGYLAYwH8K1Q2a4Avg/gIAAE4F+FsutSan/VoFQxgVwtOQMA2LTJ\nbXu9bxs3upcB3IxJ795+5ZJQ657Bj3+crG6fMJ+rUd+xw217wQ4bz+AQAHOVUguVUjsAPAxgXHAD\npdQipdS/AYRvyTEAJiml1iml1gKYBGCsqZJq6i37oPevrq669tXXM7juOvsyvp5BsNzOnfblklAN\nCeQkbXz//WR1+3gGe+3lVkdW57q1YSMGvQF8GPi8uPCdDeGySxzK1hRKsRAQVVfOwOXGW70aOPdc\n9zrSyBk0NLjX60OtewZJ8fEM2rd3q2P33flVwkTp4jnALn1uuOG6T0YVjB49GqNHj65ga9KnqYnF\noNo8A5eb+4MP+LVLF+DII+3LpeEZZIXkDMrjclx8BhsAwNlnA//8p1uZWmHy5MmYPHlyi/xvGzFY\nAqBf4HOfwnc2LAEwOlT2FdOG11xzndf6+dVCUAxcDUIaRk8pv56US1v1qBDXxfjS8AyyMnw+Q0u3\nbSu+r5ZVS/X16lPOdVvX61u3a12ryzw27yhPnDgxtf9tc7qnAtibiPoTUT2AcwA8XWb7oMl5AcAJ\nRLRbIZl8QuG7VkdTExtjH88gjecYZzEuXk/kSiIGvp5BVl6Cz4qgwX3K0jNIUpfvI0F9cgau507C\nQy1DrBgopRoBXApO/s4E8LBSajYRTSSiUwCAiEYS0YcAzgLwSyJ6p1B2DYDrAUwDMAXAxEIiudWh\ncwZ1df4TuXzrDb66koVnkMYS1nn2DJK0c9gw9/rSCBO5Lu8RrttlWxGDfGCVM1BKPQ9gSOi7CYH3\n0wD0jSj7WwC/ja/DpiXVSzBM5HqTpuEZ6AlvvuVt0P/fdfhsGp5BnsXAd//C9dkawSQJ5DZtuJzv\nch1ZeAZCyyAzkDNCi4HPaKKksV/AvyfrY1BcBS8NzyDrMJGv2LkezyRhKZ/r5rOf5VefjoNrnb5i\nIOLRMogYZESSnEHWYaKkS0on8QyyWvPHt55ly9zrSyIGPsngJJ6BLtO9u3vZYN0udWV5PwjR5EYM\nal3tk+QMsk4g+3oGQTFwKZfG0NKsDMQbb7jXFzwWvp6Bj/j4PEFMt8933kYWnoGIQcuQGzGodZLk\nDNLwDLJcRTQrzyC4revQ6yzX7Q/un+vsWR8x0G275hq3uoCiGPhec1nkDNLwJIXm5GbSWa0TDBNl\n6Rlo4+NrvHw8A5+cQZcu/OcrWo89Zl8uXNaHLI5nsB4fcb3oIre6gORikLVn4DswQmhObjyDWlf4\nSg0t3b69WL8tlfAMBg/m0ERWS1gn/R8+x8VnwmESzyDJY0urxTPYssWtrBBNbsSg1kkaJtplF796\n9QqPrsarQwfg5JP9PIPNm92Nl89IqyRikHT+hY/Ra2hwDxMl8Qx0R8AF36RuuG6Xum6+2a2OYNvu\nuMOtrBCNiEFGJJ1n4Dvu29czIPKP/c+d617OZ6RVJcXAx+h17JjN0FK9Tz6rezY2+l2jGp/rZdYs\n4KGH/OrYsMG+nFCe3IhBrYeJkuQMksRFtWfgamR91lFKMiooa88g6YQnH6PXr59/AvmRR9zL+IpB\nu3bZhokA4MYb/coJ6ZEbMah1kg4t9X3uq48YaOHKaiZxNXoGrsdz4EAWA98E8muvuZephBhkdb1o\nar0TmSUiBhmRNGfgKwY+o0N8hcv35vb1DJLgE4tfvrz43vW41NWxuGYxtDSJZ9DUVBnPwLXceee5\nlxPKkxsxqPWTGuz9NjYCLzis3arL+uA7tDTYVpdypve29RHl2zN4OrBWr68YZDG09Pvf51efx0Pq\ndYmy9gxcw5GDBxffC+mQGzGodbZvB+rri0b9gQfsyyYJEzU2shHKMoEcfm9Tn48nEtync86xLwcA\nehl4l+MSPAeuyy7U1bGR9fUMXIylFq285wyC++Qjrq71CeURMciIrVt5uKYWA5+eug8+N7dvArmS\nnoGvWLq0M3gOXM9fmzbZeQYaXzFIsgR5VjkDfb6HDrUvJ5QnN2JQ6wq/ZYv/ei9JPYMdO4CXX7Yv\no43z448D48e7lTO9jyMoPlmvWupSrlev4nsf49W2bTZDSzVJcga+PXzfnIEL+niedBLQs6ff/xCa\nkxsxqHW2bSt98HeWnkFjI/ClL9mX0cZ561a3unxjwEHPIIseadT/iKNLF+Cww9zrSyOBnJUY+HgG\nwXp8PQOX61sfzywHHLQGciMGte4ZhB8Y4jPu3xWl/JKIvuKTdJ5BEs8gq9Ev+jy4GOdgzsA3TJS1\nZ+ArBlmNJvJdDl6IRsQgI7RB8MF3aOk55/gdV51AfvFF4JBD7Mvpm/v22/1zBr6eQVaTx7RIuno+\nOmfgKs66fSNGuJUDshtNlIZn4IJ4Bi1DbsSg1tGjejRZhIn++lf3MkCxp96pk/vNfeKJwOmn81r6\n8+a51Zd3zyDoobmKQV0dL0fhGnpTCjj3XKCv8aGy5cnKMwiKjsvxPPfc4nuX63vnThYs1wEHQnly\nIwa1flKTiEGSMJEPWnxcwxrBHhtgb/jS8Az+8Af7clH/I46gh+aTM9h1V2DjRrf26dFgtZgz8F1X\nSLdRwkTpImKQEWExcCGNGL5rOZ9x8dro6cXx6uvd6qvEze0rBj7HpaGBV3S15fXXWVAbGvyeWpbV\nPAPdU+/f3//8uXoGbdpImChtciMGtU5YDFwu/iRDS30Irk3kY/QGDSp+dqkvq4XqZswordsW38lO\nOl/kOs/g4ov59W9/A378Y/tymp07eRTZe+/Zl2lq8vMM+vQB9t47u8lqbdqIZ5A2uRGDWj+pSXMG\nPmKwcqV7GaCYQPYNEwE8GcjWMCTNGUycCAwZYl/uySeb/w8btGjdeivnU1zKaTHwicX7hlOeeYaX\nhn7mGfsyOgTjclx27Mg2bBPMGYhnkB4iBhmRRpjo/vuB3r3TbZeJJGEi7fG4zF5O4hmMH8+9X5d2\nJs3dXHSRm0ju3MmhF9cZ3VoMfJ9loXEVPJ8wUbt2yYyzi6es7yVJIKdLbsSg1gn37n3WCjr66GzC\nRb5homA4y6UXHFy11HcNJZd2Bo+h70xil/p0z9k1TOQrBmefDRx7bPGzawjGJ0yUpWdw//2ce5Ew\nUbrkRgxq/aSGPYPOnd3LJr34bcsGPQPfMJGPZ+CzUJ0OZ7kY5yRrDPmIge45u4qBTsS7ikGHDsVV\nPV3xHVqaddhm+nQJE6VNQgc0PVqTGFx+ObDnnu5lk178evKTzXY+RjYoBj6eAQCsXm1fn29uw9dD\n0/unRUvXH4fuOft6Bu3a2ZfR7QyO5Ko1zwAohqVq3W5kSW48g1onKAadOrkvZ5CGZ2Br2H173Ek9\ng0mTSicitVQ7g3kXnwSyFgDbGb7BMJGLkfX1DJKKgW/O4LXXgFNPtS/Xti1w8sn22wcZOjTZs5qF\n5uRGDGpd4YNi4POcAC0GSS5+W+Ose+q+Q0sBf8/ABd+cQZcuPFO6Vy+/NZQ0d9xhVy7Yc867Z+Az\ntHTlSmDVKjevDgD23Rc49FB+75JA7tSJl1oRzyBdRAwyIigGrkZBj1N3DRONGQPcdlvxs63BTCNM\n5CJcSVZl9W1nhw7uQyjDgwDef9+unG+YSHsG+ti4HE+fFXKVKorBtm327Tz9dGDOHPvtNdoLcUWv\n8yQJ5HTJjRjUOmHPIIswUV0d8OlPl/4fG4Kege+kLBfBS7rcRpJlM5KIga3B9BUDva1uo21YKsmM\ndSJg2jTOa7U0OlntSlo5NKGU3IhBrSt8MHnramSDYuAz2kbj6hm43mxBI+Syj0mNl+98CNfj6SsG\nOinr86QzoHhvaE8hjqYmv2WlFyzgbT/6yK19voSXdXcpJ/MM0kfEICOCS1j7hIl8ekJhI+vSU/dZ\nLz7Yw8/KM/ARA995DeHj6eIZ+D72EvDzDHzE4PTT+dV2Tamk6JCUK7IcRcuQGzGodSoRJgobWZfR\nRD4PqPdNICf1DPRxcamvro57wi5LNfgez2AC2Ses4eMZBJ9DYXvNaHHLSgzS8AwkTJQeuRGDWlf4\ncALZd2EuFxEJh4laehKYbwI5qWcAuOUNguJzxRX29YXDRLb7lzRM9J3v8KuLGLRvXxy22dJicPLJ\nwNNPu5UBSnMGPstRiGeQLjLpLCOCYuAT7mnTxi8c4pMzSMMz8Jln4EpYDPR4d5tyPuITFgMXz8An\nTPStb3H8/owzeCVYlzCRHgDgwsKF/OoqBhs2uC3cp7GdBBlEqdIwn3gG6ZEbz6DWSeoZ1NWxoXN5\nlGHYePksHAfYtzXpPIM99rDbPlguLAau7XQhXG74cLtyvqOJ6uuLo8Hq6908Ay3mgHtHS4uB7bW2\ndSsP1b3pJrd6fMJEwWHWkkBOl9yIQa2f1LAY+Iwm0qEQl3Hjvp6BLufS+0rqGTQ0FOt3bafLxLMk\nYqDru+wy+xVkfcNEwXa6dATCnoHrvaW9qx/8wG57fX1+/vP8gBtbfBLIWlgBCROljYhBRiQJEwUT\nZi5GTxvLqVPdYurhUUFZeQY+k6t0fb45AxfCq7La1qfDRD5hN91OH8/AVwx0ufnz7bb3Gfo8fTqw\ndKmfZ+B7HwnlyY0Y1DrB+KhvAhlw6yFqMRg5kkMarvMMdFt9RgX5lNPGztbIBo+LS5gojZyBqxj4\neAbBdvp4BrrssmX2dQaxfYa1jxh8+CG/ajGwFeck4VahPLkRg1o/qeF5Bj6eAcBGwScc8tZbxRsw\njrBn4DMqyGeegTY+LS0GaeQMXDwRHSbq0AHYtMmvviSegW05jetCfD6j3bQIiGeQH0QMMiKNMBHA\nN4+rZ6D51rfsyiXp4SeZZ3DSSfzZRwyyyhkkCRPttZfb7N5wzqClw0RHHln6nA3bx236eAY6L5FU\nDGrdbmRJbsSg1kni3gZDTD5hIs0ZZ7iXyyKBrD2D3/3ObXnv4GiULHIG4eU2XMNESUaD1de7h4kW\nLLCvCwD22Qe45ZbiPg4caFdOX9suHQAtBh07urVRwkQtR27EoNZPatLRRD6x42DY5rzzgGHD3Mtl\nkUBOYmSzzBkExcelnbpcklVgfcJEf/87f7Y9f8GROgCw99525Xw8A72q6i672G0frguQMFHaiBhk\nRFphItecgY+RTSNM5Ls2kauRzTJnEDSWPqKVxDPwSSAHP9u2Uz++EvDLGdjW1aEDcMAB7h6ahIla\njtyIQa2TxnIUgH/OwMV4+SaQfZej8BWtJDkDnzBRWAxc1ybK2jPQ2F5r+olsuqxtW33EQF+bum0u\nS3voc/7008CVV9qVE+LJjRjUusKnMekM8A8TZTVENMk8A8B9lE5QJH1XSfUJo/iMJvLxDPR5yLNn\n4LMQXxpiMG+eXRnBDisxIKKxRPQuEc0hoqsMv9cT0cNENJeIXiOifoXv+xPRZiJ6s/B3V1QdtS4G\nwRs0aZjIxyi4egZpzDPwCUu59LiDx2XrVvtx8b7LdCQNE/kutQ0k8wxcPRjXcj5DS8Ni4NMBENIl\nVgyIqA7AnQDGABgO4ItEtE9os68AWK2U+hSAnwK4OfDbPKXUQYW/r6fU7qojzTCRT2zcN0zkIlzh\n2L9vOZ/RRLNnA1/9ql25sLF07QG7tjPYcyZyE8kkQ0s1ru30zRm4nnOiYh0+YjBqlF0ZwQ4bz+AQ\nAHOVUguVUjsAPAxgXGibcQDuK7x/HMBxgd+sorO17hkkCRP5Di31FYNwD99nNFESz8BnNBEAzJrl\nXh9gb2SDYrBoEXDPPXblwkNgbc/funWly1H4holcxUBfm1nkDPS8hoULecSbbV0AL4x39NF29Qnx\n2IhBbwDBuauLC98Zt1FKNQJYS0TdCr8NIKJ/EdErRHRkVCWtSQx8wkRJR5X4JmZdbvA33ywakKw8\ng6AYnHiiXbmZM4EVK4BTTuHPPmJw6632s4mD5VxGgz32GHDXXcVyWXkGevuWHE2kxWDYsOLM+Acf\ntK8LcA+7CeVpqecZ6H7XMgD9lFJriOggAE8S0TCl1MZwgTvvvA7du/P70aNHY/To0S3UtMqQVpjI\nVwxcbtRt24rjwF3K/elPwJQpwPXXZz+a6JZb7JfbuPdeflWKZwX7hIluuQX45jfty/mMBgOAQw/l\nVx/P4JlngFNPdc8Z6OOfhWcAuOUAWrsYTJ48GZMnT26R/20jBksA9At87lP4LshiAH0BLCWiNgA6\nK6VWF37bDgBKqTeJaD6ATwN4M1zJJZdch6FDHVtfRaQ1msg1Z+BjZLdt43HgPm1du7ZYLst5Bh06\n2CeQg7gkZoNi8JnP2IcogmEiF8/g8MOBc88tlnP1DI44oli/DUk9A50QDs98NxEepGCL7zDrWiHc\nUZ44cWJq/9vmNEwFsHdhZFA9gHMAhB9y9wyA8YX3ZwN4GQCIqHshAQ0iGgRgbwDvmyppTWGirIaW\n+oaJ/v1vYGPBd3MNaem2+c5AdhG7SoqBiyFKUk7vn49n0LUr8NJL7mKgn9PgIwa210u4A2BLa/cM\nWpJYMSjkAC4FMAnATAAPK6VmE9FEIipEXXEPgO5ENBfANwAUntqKowG8TURvAngUwMVKqbVp70Q1\nEEwCZzW01LfHfcMNwJ//zO+3bOF1513xXTiuQ4fis3jjSEMMXHrcDz1U9HxcevhBo+5SLuhR+A4t\n9Zkcd++9wC9+YVdOqdJr27YTEPZabfGdjS/EY5UzUEo9D2BI6LsJgffbAHzeUO4PAP5gV4fNVtVL\neAlr15yBDtu0awfcfTfQowevMlkOX88gyMqVHG6wbe/hh7vXt3Nn8VGLHTqwANmW08ayocFeDE44\nAfiv/+L39fX25dauBX75S+C733X3YJJ6Br4JZJ8hsA0NQK9edu3URt11XkpaYSIRg/SQGcgZkdbQ\n0rZtgRfUo3k5AAAgAElEQVReAK64wq6c7zMUfBgyBPj1r/m9ixHSyyAAbmLg6xnU1QFduvD7vn15\nmKgtWpR9w0RJPAOf8KCLwQyHs2zKhUd0+YiBhInyQW7EoNZXH0xzaKn+Lo40PAMXfG/UnTuL+9W+\nvXusGnATgxdeKNa3yy5u4aXf/IZffcNEWXsGrufBNRwZnrWcRAxsPITWnkBuSUQMMiKtoaW6Z2or\nBj6jifbYA7joIvv2BdupDYNSxTHycfj0SHV9vjkDva5N+/bA4sV2ZXr3Bvr1K7bTRbSCnoFPOV/P\nwFW0KuEZBL3XlmijYIeIQUakNZpor7341TZJ5+MZHH00cNxx8duFCfYs//QnYPXq8tsHy2UtBrqO\nFSuAq5qtthXfziRGNgvPIDg6qyVFK80wkY9nIGKQHiIGGZE0TKTL6oeB+ISJfFaGdCFYbto0+3KV\nEAOd6NZ12RzPYDjL16h362Yvkr6ege/oM59RTyYxsDmWJs/AZllx35VqhXhEDDIirTCRNrKunoHr\nWkFJxcCFJGLgk3iur8cnExxdFmVLo529e9uHpXw9g/DoJR9xtRW78Dm37XSkMZrIZdisEI+IQUYE\nDWySMFGnTvxqc6Mmmdnr+1hI3U49dNOGNMIovkNSfcXAN0zUpw+wJDx/P4KgUW/f3n7+hU8iOFxf\nEs/AVQyC37nUl8WgiNaEiEFGBA2sa5goKCQ6TGQTakhjoToXgsbk7rvt/0caYSLbMEp4kpTGptcd\nHALrK1pduxYnrrmUcwmD+S5/4TPqKYkYJO1wSJgoXUQMMiJpmCg84sLGMPiOJvINEwWNia7PZj+T\njMN37SXqY6mPi361EYM05gv4lnOZVOfraWWZQDY9etSmnISJWg4Rg4xIazSRHvliKwa+CeSkvTYi\n+5s12ONu08Z+aWifXmLQ4IXbUI6mpuaT/7Jcm8g1DOazZIPPWkhJxCB8jbkmkF06HEI8IgYZkdZo\nop497cv5JpDTGE0E2BuioLF87DHg6qvt6guO7rEVnvADcWw9Ay0i4QX1bD2fpGsTNTQACxbYx9V9\nPQPXhLXvaCKTKLuKgfbwat12ZIWIQUak5Rm4kORJZ7q+FSuA3Xf3a6dt/D8oBi6L4oU9Ch8jZJtA\nNs209TF8vp5BQwOweTPw4ovlyzQ1lcbjW3qSW/ic25YLizLgLgaAJJHTRMQgI4IXcceO9i5/uKwL\nSRLIupzr8hBBg2lr+HbsKC5U17GjXV1A83H/tnMFgsfyttv4Na4XHBYDXaer2PkaZz3z/OOP7csE\n2+jqwdiukhquzzbRbbqm9cOtXMpJEjk9ciMGtX5CwxPHNjZ71ptdWQA46ii73npYDObMATZscKvP\n1iiEe6SAfUhk27aiGJji+VEEjaxOysd1KsLGa8gQYOTI+GGbQS9E0769neELj9JxWRpaH8+GhuL3\ncXX5eDAmj8I2qR68Nn3F4J13gN12cy8nSeT0yI0Y1LpnELyxd93VPkmqywZvgHvuKa66WY7wlP+p\nU4FLLnGrz6WHGO7p2Rq+7duLj9kcMCB+e03Y8Nl4P6bwxK67xotz0AvR2M4mTjJKR58/2zWpTPtn\nI8rh+nzDRC5iEPZgbENuEiZqGUQMMiJ4EdsYn3DZYI/b9oYLJ5ABYM0at/p0m+NuuCRiEPQMJk2K\n314T7q3bGBRT4nLyZH5OQTlMYaLdd7cTA595FGGjrs9HXHjRtH824bqomcSu59322gzvn++8BgkT\npUfNisHSpex65oVwmGjTJvshceEbwHY2ajhMBLjP8gTsvAOTGDQ12YnPli3Fnm/37sVn6cZh8gzi\nDK2p53zhhcD++7vVBfAEMlsxcJ0RHDUEVk86tKlLY1OnKbFuWy6NMJHt0NkdO0o9NAkTpUfNisGZ\nZ8bf4FkSHh/t8oStrMUgHJayFYOw8VqyBBgxIr6+deuKYa82bfjPdkRK2DD4eAb9+sWfC5MY2A4E\n8FkryGTUv/CF+BE3UZ6Bj0jahIrSEgPbY2kSA/EM0qFmxWDz5nT/X1LCF79LEtn3hkviGQTDUvX1\n8eJjMia2bNxY2uO1fQ5y2EDbGr2wsdSeWjlMCeSOHeOvM31duw71NBl1m06Aaf98PANdzmaUVdA4\n2/bww/U1NNiLges5F+yoWTHIW28hfPF36wasWmVXNnzD6R5buWOmjX54zfgsw0S2BIeWAvYGJa0E\n8j/+Afz0p/Hlwglkm95suD7fUTqAvRj4JJCjPAPXIbe2y2aYRiFt3x5vB1rCM9BLY7T2mcw1KwZ5\nO7HhHk337vFjxjXbt5feAETxN2p47RfdM/UJE7Vv3/JiENw/m566LpdGAtlGlE09bj0RzKW+Hj2A\nlSvj6zO108YjDB9LwC+BDNiFiUxiYOsZhGeC2wzVDe+f7Wi3cmiBTfp/qp2aFYOWYNMmXhvIpycS\n7lna3jSA+QaPMwymeDpgP8szHCZqKTFQqrnh69TJbj6EaR99esC33Qb06hVfzidnEC7XubNdeNDX\nMwh7WYC9Z+ATJgoLskuYKLx/NveEKSzl8lAjE/qYukwErUVEDBy46irg5puBhx92Lxu+aVzEYPv2\n5jd4nGEI1xdMXsfhGyZymTAWbGdwfDtgP/TWJ0xkamfnzvEzn6M8gzgPJmzUXXrOPjmDsBcJtLxn\n4NPJMe3fmjXAokXly4U7Ri73URQiBoyIgQMffMCvPj1g001j26OJ8gzKGYYoz8DGYIdXlPSJHQMs\nnp/+dHy58L7ZeAbao3BdC8nU47Y5F6YEss3kQd/lGkzttA0TmTwDnwSyb87ANt9juo/inkcdvhfE\nM0iPmhWDlsgZ6F6nbeI3XC5oYNPwDMrdBGHjFU4kx7XX1TMwicGZZxafzFaunWEx6N6dF8iLqy/s\nUdiM0PKNxZtEy3bmcrA+3dv28WCSeAY+Imk7miitMBEQP49CPIOWo2bFoByLFvHYcldOOYVfXcUg\nKsTw5JPAE0/YxYHDN3icYQgbL31T2wxr9Jln4DtO3bRve+8NzJoVX59pRnBcUt63x22qz1YMwolS\nm6GzvjkDU8chiWfgmkB2mXsR3r+DDwZOO618OZMYJPUM9LUtYpATshSDefOADz90L6fj28uXu5Uz\n9Sq1GJx1Fq81FFdv+AaPM2CmkTaA3fh93wSyT5jBJAbDhwPvv1++nOmY2oqBr2eQRphI12ezrIRv\nmCh8XHyHlvp4BklyBvvt5955cHnoTxT6nogaGfb228CPfpSsjmogN2IQdJvXrWu5en71q2IIYtMm\ntzWCtm0DDjoI+Pe/3eo0xZv1KpRAfEjL5Pq7egZHHglccYXfonO2k858R6OYwi82iVmftYKihl7q\n31zq8/EMALverEnsfD0DmwSyqT5bMU8rZ2BTX3j/0gwTRYnBz34GfPe7yeqoBnIjBvoi2LCBlyZI\nGvM3lV+zBrj4YuBLX+LPRxzBxt2WrVuBgQOB6dOBr3zFvpzpRvvDH4rvd921fHmTZ+AzmujMM+2X\nsahkmKihAfjzn+Pr8wkTmeoD4nvdJkG2yVGYjLOvJ9KSQ0tNHZaWHk1kCoPFXWdbtpR2pNJMIEeJ\nQfi8L10KPPRQsjrzSG7EQF9A+mZ2WeLZlvAIlbfeAubOtS+/bRtPGtq8Gbj3XvtyphsteOHF3XAm\nQ2Qzz8DHmADNw0S2k87S8gwOPri4cF0UpmNiIwYm4wzEH8+tW5u3ycYziBKDOINpumZsDJ/v0FJT\nh8UnTGSzRAcQPZTVVQzS9AyibE74uDzwQLFDmQWNjdlMiMutGNisdlkOnYMI5iLWr+fXHj38/qcW\nA83MmXblTIZ5yhTg+OOBr33NzvXX6/1rXMNEgP+zCXxHE9mU27at+b7tuisnWssZ2rBRAFrWM9i2\nzU8MTD31WbOAl14qX843TOQ7tDTKM8gyZ2ATjmxJMbD1DLp149esVjm47LLm90hLkCsxuO8+4Nln\n+fOECRzS8UWf4ODFpT0DX6HZurV4IQAcKrr55vhyJgO0117AX/4S3+vW7r3ryJKoJ3PZeAZhA11f\n7xfDtwkzmIwsETB4MPDee9HlTD317t3jR3pFiUFcHH/rVrNoxR0Xk2cwbJjdGjxZTjqL6jy01HIU\nvjmDlgwT2XoGceKRNn/8I7+2tPjkSgy+/GUWAQD4zW842WvDG29wyCeIvkCCN4/2DHxXOdy2rfRC\nnDLFPEnm5ptLx79v3BidF4i7AXzDGr49S6WaG8yhQ+O9IN/RKCYjCwCf+lT5EUUmz6B79/h1f7L2\nDEzn7/TT4wdJmM6f76Sz3XaLr88kPj7LUSTJGXTpEj8AYOvW9D0Dfc1EGXe9f9oY607lhx8CixcD\nL7wA/PCHydpQjpNO4teWCJ0HyY0YlLvIly7lAx7FqFG8bv6cOaX/L/zMANOsVhf3y9QbNfGLX5R+\nXruWH4RiIs41NoWIAD/PwHa0Rrt2pTmDYcNKj62JcjmDcj0ak5EF4nvdJjHo3Dl+5rLJWALxi86Z\nwlnamyg3gczUU//nP4GJE8u3MyrnM2eOe302i+P5dh5MOYN584AFC8qXM4lBz57ARx+VLxd8EBKQ\njmewYAFPjow6//r61cdCX2NDhwJ9+/Kw0+99z1x2+XLg1FPt2hG+T+6/H1i2jO0H4D6/yZXciIFe\n6sHET38KjB0b/bu+ud94g1+VYkPSrVvzMFHw2cH77cfblutZfOELxRtp1Sq7B9GPGsWvuu41a6Kf\nWRwXJgo+EjJczjVnYHNzm4xer17Aiy+WL2cyXvpBNeU8sSjPIG7lUpMw2yQvTcYSiH+esam+urp4\nETF5BpddVr6NgFnM6+q4jTfc4FbfnnvGz+g21WezeqxpNBEAnHxy+XImMbBZhqQlcgYLFgD77hu9\nr7pN+jXsDerr19Tpeftt4NVX+f0tt5QXu7q60vts/HjuWGqvfMmS8vuRlNyIwXPPFd+Hh3vWRbTy\nK1/hUJI2tPqi3LiRT1DnzqUXyvr1QO/e/P43v+ETFddrevTR4qzI5cu59/I//wOMHh1dRo9Q+ugj\n7gGedlq0GNh4Bj5hIt+Ys0kM9tuv2JYootaaibtZozyDOENk8gxsH1JjEoO4fIPpuADxoSKTJ3LA\nATxEuRwmMddPjdtjj/L1hcv17Ru/AJypPlsxCIeJgGJvNgrTefcRA59JZ9u3AwMGFEV8wQKe6Bgl\n6rqTsHEjG/zHHy9tgw4JByey9uvHOdD33mO7s2MH8K1vATfdxMcrKmw3dSrw8ss8JwjgY7JpEw+D\nf+21ln2QT27EIMiAAcX3jY3FA/DNb/IIHM2993KSuXNn/qxPyurV3NObM4ddOZ2s27CBY9FAcc2c\ncmKgy736Kl9Ay5dzL+uWW4CRI4vbBXsE69cDb77JN/uKFcBjj/H3vp6Bb5gozdE9dXV8nMol3k31\nAXycda7GhK9nsHx5c6No4xlEiUHcSKT164vXWZA4MTCJue0opPDxJAK+/vXyyWdTff36AU89Vb7T\nk8QzCK9NBER34DSbNzcX87hrBTB7Bq5hojPOABYu5PDzzJlFz+Ddd/kYh4+vvu43bOBtly0DDjmk\n+HufPvwazHF9+CEPhpkzh4/Rr3/N3996K9u0WbN4MtuyZfy9tiGrVgEnnFB82FLbtnytDB8OXHkl\n8IMfuO2rC7kSg9//nl+HDCl+N2dO8Sa99VYektfYWLyZzjqLY5RHHFG8KAYM4JOh4/TavVq/Hthn\nH36/55782qNHdI8w2Et5+WWuUxugUaPYiO22W6mR/Ne/+HXYMD7h+sIt5xn4hInibtSoceONjeVj\nzlE94G7dyhvLqHbGxfGj8jC77FLesK9e3Txkpw1DOWMZJT5du5YXu7VrzecwLrexcaOfBxMlrnGC\nZ/JEdAeoXN7N1zMIi0hdHfDZz5Z26Exs2dJ82fA4z2DnTj63SReqmzGDX5ctYxFYs4YF8x//4O9n\nzy7dXnsGGzawLdlrL+Dcc4u/6/yINuyaDRuKI+K+/vXS3w4/HLj88mJ+Uf+P1atLr9+1a/ka6tuX\nPy9ezMIRl8PzIVdi0L07K/b48cXvhg1rboSuvrrYs58xg0M/w4c3vyh0rK1fPz6gGzbwQZ05k8UD\nYDGIWmtoxQpg0CDgi1/kGcPDhxcv/LPOYsPSr19psmzNGg4LDR3KPQ99Y0bNbYjr4Zt6UEC88Yrq\nWXbuXL73FSUGcUs9bNxoXqE0rre3aZN5pco4Q7RhQ/P66urYwJQzKBs2sICHiTue69aZy8X18pcu\nZeMRRC/mVk60TMYZiD8uppyIFrH//M/ocqbrJdzRMbF5c3Ojfu218UNSTeXixECPJAqO1PNJIHfu\nzJ3BYOgseA2GR86tXs29/w0b2D4dc0yxIzJoEHcUBw7kDuiwYcVyzz/PRtt0P+mFMq+/Hpg/vzjj\nfunS4nV2/PHsCTQ2AhdeyN/17s15hSFD0p8FnRsx+MY3OPTSrx/3ZJ57Dvj2t/m3+fOB22/nk7j/\n/sXQC8BewcCBzXsIDQ2c+NTJ3KVL+WR27swnTF9QQ4dyuEn36DVr1/LJ7duXt5k0ydzbGTGCw0jP\nP8+fly3jvEK3bnwjaXfZd2jpqlUskmHiEp5R4ZC4clHGOc4ziBo+G+cZmIwC4CcGAF8/5eLjUeGe\nLl38PYNyYrBkSTFPpdGJ53I9WpNxBvw8AwB45hngM5+JLmfKwQwfDrzzTnQZwHy9xHU4APN53203\nPpZRcXFTG308g+XLed+C14m+dkeNKl2V4OOPuXc/fDhfc3rF43335fOoRwrttx8PM589u7Rzt2pV\nMef23HNs4x54oBgWOvVUXqX3sst4kMxLL7E9uO8+7vRqevbk4as7dwInnsjfTZnitt9x5EYMbrut\nGNapq+MDo4drzZoFjBnDvbPhw3nk0be/Xbx5tcs5d27xRtHqri+U5cv5Ag0bkLFjgb/9rTQHMGcO\nt+W119gIDhrEPQLtqgU55hi+CE46iW+MZctYhLTR1YbiwAPN+x33XIJVq8xeRVxP1teoa8EMUw2e\nAcCCXW5YY5QYJAkTuYoBEL+uUXg8vSZODKLCYP37l7/OTMdTG8xyI2A2bzaLQdy8BpPHqztG4ZCK\nZt265ufOxnsJsm0bH/fBg3nf6uvZvgwdyr+fdlqpGMydy7ahT5/ik9j69eOHNjU28pwRgI+VJpib\nGT262OaxYzlBfOCB3NFs06b0OJ1wAv/PLVvYi/uP/+BwkqZz51IBu/12+/22ITdiYKJTJ76IgaKr\nrZ9V+6MfFW/OTp04ln/XXXzAd9+9OFrj0EP5dfnyYgI4SLC3pC96PQzv2mv5xOhk0b77Nm/jwQcX\nBeeUU4qewac+xQmpDRv4gtIjQcJ06VJ+5MXKlWYxiOvhR/XUBw1ibyoKk2Dq+uJEJEvPIMqoDxxY\nfphy1P75ikGcUY8Sg65dy5+/qOPZqVN5Qxslkt26lR8tZSrXvj0brvCEziCbNjU/f3vswfdBuXau\nWGEeFXXIIcD//Z+5zOLFxWStpndvrst2dq6ut1s3vk5GjmQh6NqV/8+RR3InUDN7Nt/LPXuy/Qg/\nC2X//flV24Z27fj/aq+8c2f+TdsxoGibunYtFdp9920+GOaoo4rnZffdi3nVliDXYgDwgV28uHhA\nzjqLZyrr8MvUqcDdd/PJAoC//730grnzTuCCC/giWLiw9KQAfND1yd9zT3bNgomgCRP4Ymhs5HrD\n6GR3//58U+gE00EHcT5j9eryT/uKG8USJQZxxqucGJTrOUd5Bnvu2TxBZlNfnGewcaOfZ2BKIAN8\no5WbuRwlIjqsF8XatdE5g3LtjBKDvfbi0GUUUZ5Wv358HUcRtX+9enE7o/YxSkRGjGgeQg1i8ux0\nmGrwYOCRR5qX2byZhSLcMQM4nBU1l8ckBp06ce7l+9+PbmOQFSvYVnTqxKHd4L3VsyeLw/z5LC7v\nvAP8139xh7JnT77+Fy0qtSFdu/K2Rx3Fw0GPPpo7gXvswds++CDwk5+Uhts6dmQ788ADvE/6OBx9\nNK+I/OijxW1/+MOiYGjv5emnOYf5wAN2+2xL7sUAKL2ZDjuM5whoRo5kly041Ovgg4vv27Xjk7dw\nYfQFeOihxaUvbryxVH31xRc1VK59e85B/PjHPG/hrbfYZezSpWhgyi2MZyMGUTmDuBi+6ebWxyKK\nqJ7zyJGcNylXn0kMunYt387Fi83GMs7I6uHDYQ44oDhaJKqcqYdfztP6+GM+/6bjUi5MtH07n3/T\nNdehQ+kw6TBRnkGc5xMl5nV1HLuOygFEid3++zcfXRPk44/N5+HBB/m3c85p/tuiRRxyNd1TPXrw\n8TSFwkxiAHBY6YYbSkMqUSxaxNfbN7/Jn8PXXocOfPw+/phnigMcAu7dm8suXGh+SmLfvjzasUcP\n4KtfZYHq25c7m+3aNb92/vd/OYrRtWuxo9qhA3fWzj67uF2bNsXRdgccAEybxnmG008Hzjsvfn9d\nqAoxsGHQII4HTp0K/Pznpb/16sWZ9969o436V7/KXgXAY7L/+79ZJEyGKszMmfysgI8+YldSJ5q1\nESg35nr33cu771E5g86duUcS5Yqb4qsAt62cGEQZk0MPLfaYTKxZYzYmQ4ZwTymK99/ncxcmzjNY\nscJ8XPbdN/qRmUo179lpdNjGtH+zZ3OvLDiKRVNODJYt4x6iaTLetdeay2iiRLlPH+4gROUN1q41\nnz+ADfvbb5t/+/BDs6EdPry8uOqJmGG++MXi+3CC9/33o4eeEkUPAohq409+wqNyXnstflLlnDl8\nvXXowMdRj+cP0qsX1//oozz081Of4g7mCy/wvWy6zjX6Wo5aniKM7rjq8Hc5iMoPAkhKzYgBwO7p\nyJHNR1McdRS75HFT5I88kl/HjOGL5LXXzAbARF0dXzz3318sc/318Y+03GWXYtLIxNKl5guFiI1a\n1MgZncgOM3w4u6JRwxpXrTL39HbZhW+gqN7zihXR9UUZ56YmFibTbNyePTnEYmrnmjUsFKb6evbk\nnphpHPbKlZy0NBnZDh24nMmwz55dOmQwSDkxiAoRAezhdu4c3RGIMrJt23L4JWqcuWkoqyZKDJTi\nvJZJlA88kOsyXZ9r1/LoJdP1QsS5KdNSGH/5S/mHSg0YYPZ+5s2LFpHvfY9H5ZQLEc6fz08s++xn\n+XP37uYRd0OGsNF95RWeoAbweWxsNHckglx6KYdwopLgYb75TfZqTO3IGisxIKKxRPQuEc0hombr\ndBJRPRE9TERzieg1IuoX+O3qwveziejENBtvy9ChbEDuuCN+29mzgd/+1q+eceNKx3KPGcMxx3IQ\ncQ/XNGJDKb6Ao5YuGDQoutet3eEw/fpxzybq0Z1RPWeAv49KPi9ebDbO++zDhsY0XHDZMm6LKYHc\npQv31v/4x+ZD6P74R+C448xCTcSjMP72t+a/LVhQfjJUr168H2F06M9E587RAvnuu6UTKMPtPO44\n9kJNRIXPABYS0wSyzZv5OjeJCMBiMG1a6blYvpw7UFu2mM97fT2Hl669luPywX39xz/YaEZ1mAYP\n5v8ZzDXt3Mm97YsuMpcB+HoPP3Rq9Wr2+oPh4DCf/jQPA40aqn3llcCxx/KonXJMmMBh33feKSa5\niTh8/LOflS/bq1dxhJENhx8OXHON/fYtilKq7B9YMOYB6A+gHYAZAPYJbfM1AHcV3n8BwMOF98MA\nTAfQFsCAwv8hQx2qlnnllVfK/j5mjFJPPdX8+3ffVWqvvZRqajKXu+supQ49VKl160q/37JFqY4d\nldqwwVzuG99Q6qqrzL8NG6bUm29Gl/v2t5t//8QTr6iuXaPbOXSoUs880/z7p55S6oQTzGWUUurS\nS5UClKqrU2r6dKWef16piROVOuQQpZ58Mrrcgw8qdeKJzb//+c+VuuCC6HJnnKHU739f+t38+Up1\n7/6Kmj7dXGbKFKUOOKD5942NSp1yilK/+lV0fX/8o1KjRin11lul32/apFRDg1Jbt0bXOXCgUtu3\nl37/0ktKHXxwdH0bNvDxPPlkLvvmm/x5zz1fUfPmRZebMYO3A7i9O3fy91/4Al+D5bjggtJtpk5V\navjw8mUeeUSpww/n67ipSanHH+e6LrywfLkf/pCPOaDU3Xfzd/re275dqe7dlfrgg/L/o9oo2M5Y\nO27zZyMGhwJ4LvD5OwCuCm3zPIBRhfdtAHxk2hbAc3q7UPmWO1o5YMKECWV//9GPmhupDz5Q6phj\nlPrOd6LLbd+u1NixSn3uc0qNHs03zNatbICOOy663AcfKNWtm1KLF5d+/957SvXo0dzIaBYt4nJh\nw3HqqRPUqadG1/fKK3wjPvGEUitXFr8fP16pG2+MLrd5s1LTpin1pS/xldq1K7fv+OOj26jLdeum\n1L//zcZk61be/rDD2NBEcffdSp10EpfZuVOpe+5RqlcvpQ47bEJkma1buU0zZ/LnpialXn2VBf4z\nn1Fq/fro+rZtU2rECN63m28ufn/ffeXPX1OTUqeeyn9agJ98kg317bdHl1NKqY0blTr2WDaae+yh\n1G9/G399KsXnbeNGFuJ27ZS65ho+px9/XL7cQw8pddRRLI5KKXXJJdEdEc327UqdeSZf01/+Mnds\nRo9WaunS+Dbq+6BHD6UWLOB9W7VKqXHj+JzUGlmLwZkAfhX4fB6AO0LbvANgr8DnuQC6AfgZgC8F\nvv81gDMMdbTc0coBcTfbxx8r1aePUhddxAbphhvYI7j2WjYY5Vi/XqnvfY/LnXaaUrvswmWjerLF\nNnGP9le/YkPy858rtc8+Sv34x+XL3XQT9/TvvJN70RMnKtXQMEG9+mr5ck8+ycZxzz2V+trX+OYc\nPFip1avLl9MsW6bU8uV22yrFBrVLF6V692bj1b0717ljR3SZzZu5Zz1ihFJDhih1xBHcK407f7/8\nJYvGBRfwPg4apNSVV5avK8gHHyjVs6dSn/88C+Tuu6vY47ltG4vbyJFsbPv04XOne+3lWLOGr7GX\nXj/3oM4AAAYmSURBVOLPNmKg2bGDPbQLL1Rq8mS77Y88ko/LsceyR7NihV25735Xqcsu4/a6cvvt\nfB8MHz5B9e6t1OWXN/ega4E0xcAw4T0VLNOuAsAJuGnTeMjsG29wvPyJJ4oT5srRqRMnqgGOw65e\nzXF40wiWIBMmcDz/uee4TI8evLzu5z5XvtyVV3Js9tlneeRR797A+edzHLsc48bx34wZHM8/7DAe\nJ11uDkaQqDh4FOefz5MA9dDe9evNI1GCNDTw0iKvv14cjFBXF78kw8UX8zIG06bx8y+OPz7++Afp\n358Tu88+y/Hu666LX+itvh7461/5b+dOnglvmrFsoksX/zh127acCxszxn77V17hc751K4+lj1qa\nJVwuydPDLr+cZ/9edx3nOqImfQoB4tQCHCZ6PvDZFCb6JPyD8mGiT8JJofJK/uRP/uRP/tz/svQM\npgLYm4j6A1gG4BwAXwxt8wyA8QCmADgbwMuF758G8CAR3QagN4C9AbwRrkApJZ6EIAhCBYkVA6VU\nIxFdCmASeGTRPUqp2UQ0EcBUpdSfANwD4AEimgvgY7BgQCk1i4geBTALwA4AXy/EuQRBEIQcQWKb\nBUEQhIrPQI6b0FYtENEHRPQWEU0nojcK33UloklE9B4RvUBEuwW2v6MwGW8GEeUuvUVE9xDRCiJ6\nO/Cd8/4Q0fjCuX2PiM7Pej9MROzbBCJaTERvFv7GBn4zTpzM67VLRH2I6GUimklE7xDR5YXva+X8\nhffvssL3VX8Oiag9EU0p2JF3iGhC4fsBRPR6oZ0PEVHbwvfpTfhNK/ng8weLCW3V8gfgfQBdQ9/d\nBODbhfdXAbix8P4kAM8W3o8C8Hql22/YnyMBjADwtu/+AOgKYD6A3QB00e9zum8TAPyPYduhMEyc\nzPO1C6AngBGF97sCeA/APjV0/qL2rybOIYCOhdc2AF4vnJNHAJxd+P4XAC4uvE804Tf4V2nP4BAA\nc5VSC5VSOwA8DGBchdvki764gowDcF/h/X0o7ts4APcDgFJqCoDdiMiwtmXlUEr9A0B4wWPX/RkD\nYJJSap1Sai047zQWFSZi3wDzkOhx4Btsp1LqA/AcmkOQ42tXKbVcKTWj8H4jgNkA+qB2zp9p//Ti\nHVV/DpVSehnC9mBjrgAcC+CJwvf3ATit8D54Th8H8B+F95+DeZ8jqbQY9AbwYeDzYhRParWhALxA\nRFOJqPDEUuyplFoB8AUMQBv88H4vQXXs9x6W+6PPY7Xt5yWFMMmvAyGUqH2oimuXiAaAvaDXYX89\nVs35C+yfXsGq6s8hEdUR0XQAywH8BeyRrVVK6WUbg+38ZB+UUo0A1hFRN3icu0qLQS1xhFJqJIDP\ngi/Io8ACEaTWsvVR+1ONQ4XvAjBYKTUCfBPeUuH2JIaIdgX3Fv+70IO2vR6r4vwZ9q8mzqFSqkkp\ndSDYmzsEHAKzxfvcVVoMlgAIPiqiT+G7qkMptazwuhLAk+CTuEKHf4ioJwC9NukSAMEnKlfLfrvu\nT9WcX6XUSlUItgL4PxRd6qrct0KC8XEADyil9NqoNXP+TPtXa+dQKbUewGQAhwHoQkTaXgfb+cm+\nEVEbAJ2VUqvhYWMqLQafTGgjonrw/ISnK9wmZ4ioY6GXAiLaBcCJ4PWangbw5cJmXwagb8qnAZxf\n2P5QsAsYWvU9FxBKexqu+/MCgBOIaDci6grghMJ3eaBk3wrGUXMGAL3I99MAzimM2hiI4sTJvF+7\n9wKYpZQKPja9ls5fs/2rhXNIRN11eIuIGsDHfBaAV8ATegGe4Bs8d+ML78MTfk37HE0OMudjwaMB\n5gL4TqXb47kPA8EjEaaDReA7he+7AXixsH+TAHQJlLkTnOF/C8BBld4Hwz79HsBSANsALAJwAXh0\nidP+gI3OXABzAJxf6f0qs2/3A3i7cB6fBMfX9fZXF/ZtNoAT837tAjgCQGPgmnyz0Fbn6zGn5y9q\n/6r+HALYr7A/Mwr7ck3h+4HgvMgc8MiidoXv2wN4tND+1wEMiNvnqD+ZdCYIgiBUPEwkCIIg5AAR\nA0EQBEHEQBAEQRAxEARBECBiIAiCIEDEQBAEQYCIgSAIggARA0EQBAHA/wcMKTYH53BmWAAAAABJ\nRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f1e2d2acb90>"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "s2e = np.arange(.01, 1., .01)\n",
      "marginal_loglikelihood_l = np.zeros_like(s2e)\n",
      "joint_loglikelihood_l = np.zeros_like(s2e)\n",
      "for i in range(len(s2e)):\n",
      "    marginal_loglikelihood_l[i], joint_loglikelihood_l[i] = Do_Kalman_Likelihood_Bernoulli_LaplaceMAP(dN, s2e[i], trials=trials)[2:]\n",
      "plt.plot(s2e, marginal_loglikelihood_l)\n",
      "#plt.plot(s2e, joint_loglikelihood_l)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "[<matplotlib.lines.Line2D at 0x7f1e2d4d4710>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEACAYAAACtVTGuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucndO9x/HPNyJF3InEZSIhF3GtS0lPtcatoqeIulQV\nOYhKaKpUXSqVCTlEeqp6O0obFRVUqUsqSByGKpEgokjEJYm7BgkiyO13/ljPyG7MZWf27Nmz93zf\nr9d+zbPXfi6/eV4z85u11rPWUkRgZmbWkA6lDsDMzNo2JwozM2uUE4WZmTXKicLMzBrlRGFmZo1y\nojAzs0YVlCgkXSRphqTpku6R1C0r7yvpEUmfSDprlWMGSJolabakc3PKe0iakpXfKKljIbGZmVnL\nUCHjKCStGxGLsu1hwPYRMVRSF6A7MBBYEBGXZ/t0AGYD+wNvANOAYyJilqQ/A7dExF8kXQk8FRFX\nFfLNmZlZ4QqqUdQliUxnYEVWPj8ingCWrXLInsALETEvIpYCNwGHZZ/tB9yabY8DDi8kNjMzaxkF\nN+9IGgWcACwE9m1i9y2BV3PevwbsKWkTUs1jRU75FoXGZmZmhWuyRiFpsqSnc17/zL4eAhARwyOi\nOzAeGFZALCrgWDMzK5ImaxQRcWCe57oBmAjUNLLP66S+izpbAa9HxLuSNpTUIatVbJXtWy9JnqDK\nzKwZImK1/ykvqOlJUq+IeDF7OxCYWd9uOdvTgF6StgbeBI7JXgD3A0cBfwYGAXc0du3WnswwAh5/\nHG6+Ob0AttsOtt0WevaETp1g2TJYvhwWLIC5c9Nr9mxYe23Yc8/02m8/2H13WGONlomrpqaGmpqa\nljlZmfO9WMn3YiXfi5Wk5jXcFNpHMVpSH1In9jxgSBZMV+BxYD1ghaQzSE9ELZL0fWASqdlrbETM\nys51HnCTpIuB6cDYAmNrMdOmwWmnwfvvw7e/DX/7G+y4I+RzzyNgzhyYOhUefRROPBHefhsOPBC+\n+U045BBYf/3ifw9mZs1VUKKIiCMbKH8bqGrgs3uAvvWUzwH2KiSelrZgAVxwAdx2G4wZA8cdl19y\nyCXBNtuk1zFZ3enVV+Hee+Gmm1ICqq6Go4+Gww+HddZp8W/DzKwgHpndgFmzYOed0x/6556D449f\n/STRkKoqGDwYJkyAV16BI4+E8eNhyy3h5JPh739PNZF8VFdXt0xQFcD3YiXfi5V8LwpX0IC7UpEU\nxYx75kw44AC49FI44YSiXeZz3nwTrr8e/vjHlCiGDEnX32ij1ovBzCqXpGZ1ZjtRrKJUSSJXRKpV\n/O53MHEiHHEEnHFGquGYmTVXcxOFm55yzJ2bksTo0aVLEpCauL72NbjhhvTU1DbbwIABsP/+qblq\nxYqmz2Fm1lJco8hEwEEHpcdXzzuvRU/dIpYsgb/8BS6/HBYvhrPPTp3rX/hCqSMzs3LhpqcC/elP\n6Y/wtGnQsQ3PWxsBtbXpKawZM+CHP0x9GX7E1sya4kRRgPnzYaed4K670mC4cvH003DZZTBpEgwd\nCj/4AWy6aamjMrO2yn0UBTjrLPjud8srSUDq3B4/Pg3ke+st6NsXzjknDegzM2sp7T5RTJoEDz8M\nF11U6kiar1cvuPpqeOop+Phj6NcPzjwzPW5rZlaodp8oamrgZz+Dzp1LHUnhqqrg17+GZ55J73fY\nwQnDzArXrhPFjBlpZPTAgaWOpGVtsQX84hfw7LOp83uHHeBHP4J//avUkZlZOWrXieKqq+CUU9r2\nU06F2HxzuOKKVMNYsiQ1SZ1/Prz7bqkjM7Ny0m4TxYcfpkn5Bg8udSTFt8UWqUlq+nR4773U6T1y\nJHzwQakjM7Ny0G4TxY03wj77pIn42ovu3VMtaupUePnl1Al+2WVpAJ+ZWUPaZaKISPMoDRlS6khK\nY5ttYNw4ePDBtBhTr17w29+m5ikzs1W1y0QxbRosXJgWD2rP+vVL04JMmJAWY+rbF667Lq3SZ2ZW\np6BEIekiSTMkTZd0j6RuWXlfSY9I+kTSWascMzfnmKk55RtJmiTpeUn3StqgkNgac9VVcOqp0KFd\npsnP2313uPvuVMu4+mrYZRe4/fb818Qws8pW0BQektaNiEXZ9jDScqdDJXUBupPW0V4QEZfnHPMy\nsHtELFjlXJcB70bEGEnnAhtFRL3T8xUyhceSJdClS5qVtWvXZp2iokWkqUwuuCCt9X3ppbDvvqWO\nysxaQkmm8KhLEpnOpLWziYj5EfEEsKyew9TAdQ8DxmXb40hJpsVNnZra5J0k6ieltbynT09zRw0e\nnKY4f/LJUkdmZqVScOOLpFGSXgGOBS7M45AA7pU0TdIpOeWbZWttExFvAZsVGlt9amv9H3I+OnSA\nY49NCzkdeij853+mNb9ffLHUkZlZa2tyqJmkyUDu/98i/bG/ICImRMRwYHjWXDQMqGnilF+JiDez\n5qnJkmZGxMP17Ndo21JNzcrLVFdX570u7gMPpGktLD+dOsFpp6WFnK64Avr3h6OPhgsvhG7dSh2d\nmTWmtraW2trags/TYtOMS6oCJkbETjllI4APc/soVjnms88lzQSqI+LtrFP8gYjo18Bxzeqj+PTT\nNA33a6/BBkXrKq9s77yT+i2uvTZNbf7jH/tempWLkvRRSOqV83YgMLO+3XL2X0fSutl2Z+DrQDaF\nHXcC/5VtDwLuKCS2+jz2WHok1H/Ymm/TTeHnP099GK+/Dr17p/effFLqyMysWAp96ukWoA+pE3se\nMCRrVuoKPA6sl322CNge6ALcRmpW6giMj4jR2bk2Bm4GqrJzHR0RCxu4brNqFCNHwkcfpdXhrGU8\n+yz85CdpivORI+H442GNNUodlZnVxyvc5aG6Gs49Fw4+uOVjau/+8Y90bxcuTE1T3/xmeoLKzNoO\nJ4omfPxxGj/x5puw3npFCqydi0gjvM8/HzbcMM0j9ZWvlDoqM6vjpVCbMGUK7Lijk0QxSXDIIWmd\nj1NOSY/XHnroyoWUzKw8tZtE8cADHj/RWtZYAwYNguefT/d8v/3gxBPTIlFmVn6cKKxo1lorjVl5\n4YU0nfuuu6aV9rxwkll5aReJYvHi9Din28tLY4MNYNSo1AS1eHGapfa//zs9gWZmbV+7SBSPPQY7\n7wydO5c6kvZt883hyivh0Ufh6afTGIwrr4SlS0sdmZk1pl0kiueeS4nC2obeveHPf07rYNx2WxoE\nedNNsGJFqSMzs/q0i0Qxezb06VPqKGxVu+8Okyal1QZ//nP40pfS+zJ8YtusojlRWMkdcECa/v38\n82HYMNh//9RcaGZtQ7tIFC+84ETR1klw5JFpSpBjjoEjjkivmfXNHmZmrariE8WSJWm22J49Sx2J\n5aNjR/je91Jy798f9tkHTjrJYzDMSqniE8XLL0P37rDmmqWOxFbH2munKcxnz4YttkhjMM48E+bP\nL3VkZu1PxSeK2bPTUzZWnjbcMI3BePZZWLYMttsOamrggw9KHZlZ+9EuEoX7J8pft27w61/D44+n\nWmLv3nD55V4Hw6w1OFFYWenZE667Dv7v/+Chh1LC+MMfUm3DzIrDicLK0o47wu23wy23wI03wvbb\np0F8HrRn1vIKXQr1IkkzJE2XdE+21jWSjs3KZ0h6WNLOOccMkDRL0mxJ5+aU95A0JSu/UVLHQmKr\n40RR2fbaK9Uu/vd/U1PUbrulNTE8aM+s5RS6FOq6EbEo2x4GbB8RQyX1B2ZGxPuSBgA1EdFfUgdg\nNrA/8AYwDTgmImZJ+jNwS0T8RdKVwFMRcVUD181r4aIPP4SuXWHRIuhQ8XUni4A774QLLoD114dL\nLkmrGppZUpKFi+qSRKYzaX1sImJKRLyflU8Btsy29wReiIh5EbEUuAk4LPtsP+DWbHsccHghsUF6\nFr93byeJ9kKCww5LCyedfjoMHgwHHphGfZtZ8xX8J1TSKEmvAMcCF9azy2Dg7mx7S+DVnM9eA7aU\ntAmwICJW5JRvUWhsbnZqn9ZYA7773TSq+6ij0gjvww5LM9aa2eprsh9A0mSga24REMAFETEhIoYD\nw7P+hmFATc6x+wInAnvnEctqVYdqaj67DNXV1VTX08bgRNG+rblmGuV9wglp4sGvfz01RY0cmdbE\nMKt0tbW11NbWFnyegvoo/u1EUhUwMSJ2yt7vTGpKGhARL2Vl/Un9FQOy9+cBERGXSZoPdI2IFdl+\nIyLi4AaulVcfxXHHpaaHQYNa4ju0crdoURqLcfnl8M1vwoUXemoXa19K0kchqVfO24HAzKy8OylJ\nHF+XJDLTgF6StpbUCTgGuCP77H7gqGx7UE55s3lUtuVad900Q+0LL0BVFeyxBwwdmuYCM7OGFfrU\n0y1AH1In9jxgSES8Ken3wLeyMgFLI2LP7JgBwC9JSWpsRIzOynuSOrc3AqYDx2Ud3vVdt8kaRQRs\nvHH6o7Dpps3+Fq2CvfMOjBkDY8fC8cenJNK1a9PHmZWr5tYoWqzpqTXlkyjmz0/t0O++m56GMWvI\nW2/B6NFpxPfgwXDOOf7nwipTSZqe2rK6jmwnCWtKt25wxRXpqagPP0z/YAwfDgsWlDoys7ah4hOF\nWb622gquvBKeeCLVMnr3Tk9Ivf9+08eaVTInCrNV9OiRJhp87DGYMwd69UpTnXtqc2uvKjZRzJ3r\nRx+tMNtuC9deC//4B8yalRLGpZemx2zN2pOKTRTvvANdupQ6CqsEffrA9dfDgw+mfoxtt4XLLnPC\nsPajohOFn1yxltSvX5rS/P77Yfr0lDDGjIGPPip1ZGbFVdGJYpNNSh2FVaIddoCbbkrTmz/+eEoY\nP/uZE4ZVropNFO++6xqFFdeOO8LNN8N9961MGGPGuEnKKk9FJorFi9PI7HXWKXUk1h7suGNaXS83\nYYwencZkmFWCikwUdf0THmxnramuhnH//fDUUylhXHKJH6u18lexicL9E1YqdX0YDz4Izz67chyG\nB+5ZuarIROH+CWsL+vWD8ePh739PA0B79UojvT01iJWbikwUfjTW2pK+fdOEg48+CvPmpYQxfHj6\nOTUrBxWbKNz0ZG1Nr15wzTWpw7tuduNzzoG33y51ZGaNq8hE4aYna8t69oSrrkod3h9/nJqozjjD\nCyhZ21WRicJNT1YOqqrS0qzPPpvW9955ZxgyJE1EaNaWFLoU6kWSZkiaLukeSd2y8mOz8hmSHs7W\nz647Zm7OMVNzyjeSNEnS85LulbRBc+NyorBysvnm8D//A88/n1Zl3GOPtM77rFmljswsKbRGMSYi\ndomIXYG7gBFZ+cvA1yJiF2AUcHXOMSuA6ojYtW551Mx5wH0R0Ze0fvb5zQ3q3XfdR2Hlp0uXNO7i\nxRdTf8bXvgZHHw0zZpQ6MmvvCkoUEZE7WUFnUhIgIqZERN1T41OALXP2UwPXPQwYl22PAwY2Ny7X\nKKycbbQR/PSn8PLLsNdecPDBcMghMGVKqSOz9qrgPgpJoyS9AhwLXFjPLoOBu3PeB3CvpGmSTskp\n3ywi3gaIiLeAzZobkxOFVYJ114Uf/SgljG98A445BvbfP01GWIZL3VsZUzTxEydpMtA1t4j0x/6C\niJiQs9+5wNoRUZNTti/wG2DviFiQlW0eEW9K6gJMBr4fEQ9Lei8iNs459t2IqLcBSVKMGDHis/fV\n1dVUV1d/9n7ttVPzk+d6skqydGma5vzSS2G99eAnP4FDD4UOFflIirWE2tpaamtrP3s/cuRIImK1\nJzdqMlHkfSKpCpgYETtl73cGbgUGRMRLDRwzAvgwIi6XNJPUd/F21in+QET0a+C4aCjuxYtT/8TH\nH7fAN2XWBq1YAbfdlvozPvkEzj8/1TY6dix1ZNbWSWpWoij0qadeOW8HAjOz8u6kJHF8bpKQtI6k\ndbPtzsDXgWeyj+8E/ivbHgTc0ZyY3Oxkla5DBzjiiDRw7xe/gLFjoXdv+O1v/Q+SFUdBNQpJtwB9\nSJ3Y84AhWbPS74FvZWUClkbEnpJ6AreRmq46AuMjYnR2ro2Bm4Gq7LijI2JhA9dtsEbx5JNw8slp\nBTKz9mLKlNQk9dhjafDe0KGw4YaljsramubWKFqs6ak1NZYoJk9O6xnfd18rB2XWBjzzTPr5nzgR\nTjkFfvhD6Nat1FFZW1GSpqe2yE1P1p7tuCP86U/wxBNpadbtt4dTT01jM8yay4nCrAL16JGmB5k1\nCzbbDL785TR474knSh2ZlaOKSxSeENBspc02g4svTmMx+veHgQPhwANTE20ZtjpbiVRcovAU42af\nt956cNZZ8NJLcNxxqe9ijz3SWt/LlpU6OmvrKjJRuEZhVr9OndKEg//8J9TUwG9+A336pEdrFy8u\ndXTWVjlRmLVDHTqk+aP+/ne4/vrUFNWjR0oe8+eXOjprayouUbiPwmz1/Md/wO23p6Txxhtp5b3T\nTvOTUrZSxSUK91GYNU/fvnD11fDcc2kG2y9/GY480rPWWgUOuPOEgGYtY9Ei+OMf4fLLYaut4Oyz\nU3OVJyEsXx6ZTeqM23RTd8qZtaRly+Cvf4Wf/Qzefz89PXXCCf5nrBx5ZDZudjIrho4d02C9qVPh\n979P04P06AEXXghvv13q6Kw1VFyicEe2WXFIsM8+cOed8NBDKUlst12aU+q550odnRVTRSUKP/Fk\n1jq22w6uugpmz4aqKthvv7QKn1ffq0wVlShcozBrXV26pCaouXPhW9+CH/wAdt0VrrsOliwpdXTW\nUiouUbiPwqz1rbUWDB6cpjkfPTrNYNujB4walX4vrbwVnCgkXSRphqTpku7JljFF0qE55VMlfSXn\nmEGSZkt6XtIJOeW7SXo6++yK1Y3FNQqz0pJgwIA00vvee2HOnLT63qmnwsyZpY7OmqslahRjImKX\niNgVuAsYkZXfl1N+MvAHAEkbARcCXwL2AkZI2iA75krg5IjoA/SRdNDqBOI+CrO2Y6ed0jKts2bB\n5pvDvvvCwQenBOJ+jPJScKKIiEU5bzuTlkUlInJHM6xbVw4cBEyKiPezpU4nAQOymsh6ETEt2+86\n0jrceXONwqzt6do1zSE1d256zPbHP04LLF19tcc8lYsW6aOQNErSK8CxpNpCXflASTOBCcBJWfGW\nwKs5h7+elW0JvJZT/lpWljf3UZi1XWutBSeeCDNmpEWV/vY32Hpr+MlP4PXXSx2dNSavRCFpctZ3\nUPf6Z/b1EICIGB4R3YHxwLC64yLi9ojoR6oZjCrGN5DLTU9mbZ+UHqe980545BH48MPUTPWd73he\nqbaqYz47RcSBeZ7vBmAiULPK8Q9L2kbSxqQaRHXOx1sBD2TlVauUN/h/Rk3NyktUV1dTXV3tpiez\nMtO7d6pdjBoF11wDxx6bHrk944w0IWGnTqWOsLzV1tZSW1tb8HkKnutJUq+IeDHbHgZ8NSKOlrRt\nRLyUle8G3BERVVln9uPAbqQazePA7hGxUNIU4AfANFLH+K8i4p56rlnvXE/rrJOanzwHjVl5Wr48\nNUn98pfw/PMwdCh873tpSVcrXCnnehqdNUM9BRwAnJGVHyHpGUlPAr8GjgaIiAXAxaQE8RgwMuvU\nBjgdGAvMBl6oL0k0ZMUK+OSTNHusmZWnNdaAww6D+++Hu+9OHeB9+6a+jenTSx1d+1Uxs8d+/DFs\nvHH6amaV45134A9/SMu19ugBw4bB4YfDmmuWOrLy0+6nGV+4EHr2hAULShSUmRXVsmVpJb5f/Qpe\nfjk1S51yipulVke7n2b8k0/S43dmVpk6dkwd3A89lPox5sxJzVKDBsHjj5c6uspWUYniC18odRRm\n1hq++MXUHPXii7DDDimB9O8P48fDp5+WOrrKU1GJwjUKs/Zlk03gnHPgpZfg/PPh2mvTIL6f/hRe\ne63Jwy1PThRmVvbqnpaaPBkeeCD1We68Mxx1FNTWem6pQlVMovj0UycKM4N+/dIgvrlzoboaTjst\njfy+8kpYtKipo60+FZMoXKMws1zrrw+nnw7PPpuelJo8Gbp3T4/Xesrz1VNRicKd2Wa2qrq5pf76\n1zQh4YYbpinP998fbr0Vli4tdYRtX0UlCtcozKwxVVVw8cXwyitpRb5f/jIN4hs5Et54o9TRtV1O\nFGbW7nTqlGarfeghuOceeOuttEbGkUem6UPc+f3vnCjMrF2r6+ieOzc1UZ1xRuoQv+IKz/RQp2IS\nhZ96MrNCrL9+ekLq6afh97+HqVPTtEAnnQTTpjV9fCWrmEThzmwzawkSfPWrcMMNMHs29OmTlnDd\nY480Gvyjj0odYeurqEThGoWZtaTNNoPzzksjv0eNggkTVj5i+8wzpY6u9ThRmJk1oUMHGDAA7rgD\nnnoqLWlw0EGw997wpz+lvz+VzInCzGw1VFWlx2nnzoWzz04TEVZVwVlnwaxZpY6uOApKFJIukjRD\n0nRJ90jqlpUfmlM+VdJXco5ZLunJ7LPbc8p7SJoiabakGyXltZ53HXdmm1lrWnNNGDgwPV772GOp\nj7S6Or1uuKGyZrEtaOEiSetGxKJsexiwfUQMlbRORCzOyncCbo6Iftn7DyJi/XrO9Wfgloj4i6Qr\ngaci4qoGrvu5hYtOPx223z59NTMrhSVLUvPUVVelUeAnnJDW/O7bt9SRJSVZuKguSWQ6Ayuy8sU5\n5evWlWcaCnI/4NZsexxw+OrE4qYnMyu1Tp3SjLX33QdTpqRaxz77rKxllGtfRsF9FJJGSXoFOBa4\nMKd8oKSZwATgpJxDvpA1Rz0i6bBs302ABRFRl1BeA7ZYnTicKMysLdl2Wxg9Ok0XMmwYjBuX+jLO\nPBOee67U0a2eJhOFpMmSns55/TP7eghARAyPiO7AeGBY3XERcXvW3DQQGJVzyq0jYk/gu8AVknrW\nXaqQb8SJwszaok6d4Igj4N57U1/G2munCQn33jslj8WLmz5HqTXZYRwRB+Z5rhuAiUDNKsc/LGkb\nSRtHxHsR8WZWPkdSLbBrRPxV0gaSOmS1iq2A1xu7WE3NystUV1fz6afVThRm1qZtsw1cckl6auqu\nu9II8LPOgmOOgVNOSUu8tqTa2lpqa2sLPk+hndm9IuLFbHsY8NWIOFrSthHxUla+G3BHRFRJ2hBY\nHBFLJG0KPAIcGhGzss7sv0bEn7PO7BkR8bsGrvu5zuwDDkgDYw44oNnfjplZq3v1VbjmGhg7Frp2\nTbPafuc7aUqRllaSzmxgdNYM9RRwAHBGVn6EpGckPQn8Gjg6K+8HPC5pOvB/wCURUffk8XnAWZJm\nAxsDY1cnEDc9mVk5qqqCESNgzhy46CKYNCmN/j7pJHjkkbYxk21BNYpSqa9Gscce8Lvfpa9mZuXs\n7bdT/8XYsWlU+ODB6VHbLl0KO2+pahRthmsUZlYpunaFc85JI72vvjrNaNu7d1ov4+67Yfny1o2n\nYhKFR2abWaWpm8l23DiYNw8OPDA1U229NQwfDi+/3DpxVEyi8DTjZlbJNtgATj01rZMxcSIsWgR7\n7ZUWWxo/Hj7+uHjXrpg+ii5d0iCWQtvwzMzKxaefpilDrrkmLa707W+nTvDdd0+1kVU1t4+iYhLF\neuulxdHXW69EQZmZldCrr8K116ZX585w4olw3HH//s9zu08Ua66ZRjiuuWaJgjIzawNWrICHHkq1\njPnzU+d3nXadKJYvTwli+fL6q1tmZu1RxL//TWzXj8fWPfHkJGFmtlJL/U2siEThMRRmZsXjRGFm\nZo1yojAzs0ZVRKL49FMPtjMzK5aKSBSuUZiZFY8ThZmZNcqJwszMGuVEYWZmjSo4UUi6SNIMSdMl\n3SOp2yqff0nSUknfyikbJGm2pOclnZBTvlu2Yt5sSVfkG4M7s83MiqclahRjImKXiNgVuAsYUfeB\npA7AaODenLKNgAuBLwF7ASMkbZB9fCVwckT0AfpIOiifAFyjMDMrnoITRUQsynnbGViR834YcAvw\nr5yyg4BJEfF+RCwEJgEDsprIehExLdvvOmBgPjE4UZiZFU/HljiJpFHACcBCYN+sbAtgYETsK2nP\nnN23BF7Nef96VrYl8FpO+WtZWZOcKMzMiievRCFpMtA1twgI4IKImBARw4Hhks4l1SJqgCuAc1s2\n3JVqamo+237nnWrWWqu6WJcyMytLtbW11NbWFnyeFp1mXFIVcFdE7CypbjVXAZsCHwHfA9YBqiNi\nSHbM74AHgAeBByKiX1Z+DLBPRAyt5zr/Ns34ZZfBe++lr2ZmVr+STTMuqVfO24HALICI2CZ79ST1\nU5wWEXeSOrYPlLRB1rF9IHBvRLwFvC9pT0kiNWXdkU8MbnoyMyueluijGC2pD6kTex4wpJ59Pvv3\nPyIWSLoYeDwrH5l1agOcDlwLrAVMjIh78gngk09g/fWb/w2YmVnDKmKFuzPPhO7d01czM6tfu17h\nzk1PZmbFUxGJwiOzzcyKpyIShWsUZmbF40RhZmaNcqIwM7NGOVGYmVmjKiJRuDPbzKx4KiJRuEZh\nZlY8ThRmZtYoJwozM2uUE4WZmTWqIhKFO7PNzIqnIhKFaxRmZsXjRGFmZo0q+0SxbFn62rFFVv82\nM7NVlX2icG3CzKy4CkoUki6SNEPSdEn3SOq2yudfkrRU0rdyypZLejI75vac8h6SpkiaLelGSXnV\nEdyRbWZWXIXWKMZExC4RsStwFzCi7gNJHYDRpDWyc30UEbtFxK4RMTCn/DLg5xHRB1gInJxPAK5R\nmJkVV0GJIiIW5bztTFo3u84w4BbgX6sc1tAyfPsBt2bb44DD84nBicLMrLgK7qOQNErSK8CxwIVZ\n2RbAwIi4ks8nhi9ImirpEUmHZftvAiyIiLpE8xqwRT7Xd6IwMyuuJvsBJE0GuuYWAQFcEBETImI4\nMFzSuaRaRA1wBXDuKsfU2Toi3pTUE7hf0tPABzRc06hXTU0NAG+8AUuWVAPVq3O4mVnFq62tpba2\ntuDzKCIKjwaQVAXcFRE7S3q5rhjYFPgI+F5E3LnKMX8EJkTEXyX9C+gWESsk9QdGRMTBDVwr6uJ+\n5BE4++z01czMGiaJiFitf8qh8KeeeuW8HQjMAoiIbbJXT1I/xWkRcaekDSV1yo7dFPgK8Fx2/APA\nUdn2IOBTCyYUAAAIiUlEQVSOfGJw05OZWXEVOkxttKQ+pE7secCQevbJrbL0A66StJyUpC6JiFnZ\nZ+cBN0m6GJgOjM0nACcKM7PiKihRRMSReexzUs72o8DODew3B9hrdWNwojAzK66KGJntAXdmZsVT\n9oni009dozAzK6ayTxRuejIzKy4nCjMza5QThZmZNaoiEoU7s83MiqfsE4U7s83MiqvsE4WbnszM\nisuJwszMGuVEYWZmjaqIROHObDOz4qmIROEahZlZ8ZR9ovBTT2ZmxVX2icI1CjOz4nKiMDOzRlVE\nonBntplZ8RS6FOpFkmZImi7pHkndsvJ9JC2U9GT2Gp5zzABJsyTNlnRuTnkPSVOy8hsl5bWokmsU\nZmbFVWiNYkxE7BIRuwJ3ASNyPnsoInbLXqMAJHUAfgMcBOwAfEfSdtn+lwE/j4g+wELg5HwCcGe2\nmVlxFZQoImJRztvOpLWz66ieQ/YEXoiIeRGxFLgJOCz7bD/g1mx7HHB4PjG4RmFmVlwF91FIGiXp\nFeBY4MKcj/pnTVJ3Sdo+K9sSeDVnn9eALSVtAiyIiBU55Vvkc30nCjOz4mqyH0DSZKBrbhEQwAUR\nMSEihgPDs/6GYUAN8ASwdUQslnQwcDvQp6lLrU7gNTU1AHz4ITz2WDUHH1y9OoebmVW82tpaamtr\nCz6PIqLwaABJVcDEiNipns/mALuTkkVNRAzIys8DIiIukzQf6BoRKyT1B0ZExMENXCsigghYYw1Y\nsgQ65tX1bWbWfkkiIlbrn3Io/KmnXjlvBwIzs/KuOfvsSUpI7wHTgF6StpbUCTgGuCPb9X7gqGx7\nUE55g5YtS4nCScLMrHgK/RM7WlIfUif2PGBIVn6kpKHAUuBj4NsAEbFc0veBSaQkNTYiZmXHnAfc\nJOliYDowtqmLu3/CzKz4WqzpqTXVNT3Nnw/bbw/z55c6IjOztq8kTU+l5lHZZmbFV/aJwk1PZmbF\nVdaJwqOyzcyKr6wThWsUZmbF50RhZmaNKvtE4c5sM7PiKvtE4RqFmVlxlXWicGe2mVnxlfWAu4UL\nYcEC6Nmz1BGZmbV9zR1wV9aJwszM8tcuR2abmVnxOVGYmVmjnCjMzKxRThRmZtYoJwozM2tUoSvc\nXSRphqTpku6R1C0r30fSQklPZq/hOcfMzTlmak75RpImSXpe0r2SNigkNjMzaxmF1ijGRMQuEbEr\ncBcwIuezhyJit+w1Kqd8BVAdEbtGxJ455ecB90VEX9KyqOcXGFu70BILp1cK34uVfC9W8r0oXEGJ\nIiIW5bztTEoCdRp6VlcNXPcwYFy2PY60Brc1wb8EK/lerOR7sZLvReEK7qOQNErSK8CxwIU5H/XP\nmpfukrR9TnkA90qaJumUnPLNIuJtgIh4C9is0NjMzKxwTSYKSZMlPZ3z+mf29RCAiBgeEd2B8cCw\n7LAngK2zJqnfALfnnPIrEbEH8A3gdEl7N3BpD702M2sDWmwKD0lVwMSI2Kmez+YAu0fEe6uUjwA+\njIjLJc0k9V28nXWKPxAR/Rq4lpOImVkzNGcKj46FXFBSr4h4MXs7EJiZlXeta0aStCcpIb0naR2g\nQ0QsktQZ+DowMjv+TuC/gMuAQcAdDV23Od+omZk1T0GJAhgtqQ+pE3seMCQrP1LSUGAp8DHw7ay8\nK3BbViPoCIyPiEnZZ5cBN0s6KTvX0QXGZmZmLaAsZ481M7PW06ZHZksaIGmWpNmSzq3n806SbpL0\ngqRHJXUvRZzFlsd9OFPSs5Keyh4+qCpFnK2hqXuRs98RklZI2q0142tN+dwLSUdnPxv/lHR9a8fY\nWvL4HamSdH82APgpSQeXIs7WIGmspLclPd3IPr/K/m4+JemLTZ40Itrki5TEXgS2BtYEngK2W2Wf\nocD/ZtvfBm4qddwlug/7AGtl20Mq8T7key+y/dYFHgQeAXYrddwl/LnoRXoCcf3s/aaljruE9+Iq\n4NRsux8wp9RxF/F+7A18EXi6gc8PBu7KtvcCpjR1zrZco9gTeCEi5kXEUuAm0qC8XLmD9G4B9m/F\n+FpLk/chIh6MiE+yt1OALVs5xtaSz88EwMXAaODT1gyuleVzL04BfhsRHwBExDutHGNryederADW\nz7Y3BF5vxfhaVUQ8DCxoZJfDgOuyfR8DNpDUtbFztuVEsSXwas771/j8H8DP9omI5cBCSRu3Tnit\nJp/7kOtk4O6iRlQ6Td4LSbsCW0VEpd6DOvn8XPQB+kp6WNIjkg5qtehaVz73YiRwvKRXgb+xcsxX\ne7Tq/XqdJv65LPSpp7amXT82K+k4YHdSU1S7I0nA5aTHqz8rLlE4bUFHUvPT14DuwEOSdqyrYbQz\n3wH+GBG/kNQfuB7YocQxlY22XKN4nfTDXWcrPl9dfA2oApC0Bqkt9j0qSz73AUkHkCZSPCSrflei\npu7FeqRf/tpskGd/4I4K7dDO9/fjzohYERFzgdlA79YJr1Xlcy9OBm4GiIgpwFqSNm2d8Nqc18n+\nbmbq/ZuSqy0nimlAL0lbS+oEHEMalJdrAiv/ezyKNOtspWnyPmTNLb8DDo2Id0sQY2tp9F5ExAcR\nsVlEbBMRPUn9NYdExJMlireY8vn9uB3YFyD7o9gbeLlVo2wd+dyLecABAJL6AV+o4D4bSDXphmrT\ndwInAGS1q4WRDZBuSJtteoqI5ZK+D0wiJbSxETFT0khgWkT8DRgL/EnSC8C7pB+QipLnfRhDmr33\nL1nzy7yIqLjZd/O8F/92CBXa9JTPvYiIeyV9XdKzwDLg7IhorJOzLOX5c3E28HtJZ5I6tgc1fMby\nJukGoBrYJJuwdQTQCYiIuDoiJkr6hqQXgY+AE5s8Z/aIlJmZWb3actOTmZm1AU4UZmbWKCcKMzNr\nlBOFmZk1yonCzMwa5URhZmaNcqIwM7NGOVGYmVmj/h+YVAuKYjeSwwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f1e2d1c27d0>"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Conditional_on_one_axis(dN, sigma2e, x_map, axis, values, trials=1.):\n",
      "    \"\"\"Compute the joint probability as a function of one of the coordinates of the state 'x'\n",
      "    keeping all other coordinates at the MAP solution. Returns both the exact value and the Laplace approximation.\n",
      "\n",
      "    :param dN: Observations (K,)\n",
      "    :param sigma2e: Variance of process noise\n",
      "    :param x_map: MAP solution (e.g. computed by Do_Kalman_Likelihood_Bernoulli_LaplaceMAP)\n",
      "    :param axis: Along which axis the section of the joint should be computed (e.g. 3 for x_3)\n",
      "    :param values: Values of x[axis] for which the joint should be computed\n",
      "    :param trials: Number of trials for binomial observations (1 for Bernoulli)\n",
      "    :return: joint_loglikelihood, joint_loglikelihood_map\n",
      "    \"\"\"\n",
      "    binom = np.sum(np.log(sp.special.binom(trials, dN)))\n",
      "    # Build Hessian of joint\n",
      "    D = trials / (np.exp(x_map) + 2 + np.exp(-x_map)) + 2. / sigma2e\n",
      "    D[-1] -= 1. / sigma2e\n",
      "    B = -np.ones(len(D)) / sigma2e\n",
      "    B[-1] = 0.\n",
      "    L = sp.linalg.cholesky_banded((D, B), lower=True)\n",
      "    joint_loglikelihood_map = (binom + np.sum(dN * x_map - trials * np.log(1 + np.exp(x_map)))\n",
      "                           -.5 * ((np.sum(np.diff(x_map)**2) + x_map[0]**2) / sigma2e + len(dN) * np.log(2*np.pi*sigma2e)))\n",
      "    joint_loglikelihood_map = joint_loglikelihood_map - .5 * (values - x_map[axis])**2 * (L[0][axis]**2 + L[1][axis]**2)\n",
      "    x = x_map.copy()\n",
      "    joint_loglikelihood = np.zeros(values.shape)\n",
      "    for i in range(len(values)):\n",
      "        x[axis] = values[i]\n",
      "        # Compute joint and marginal probabilities\n",
      "        joint_loglikelihood[i] = (binom + np.sum(dN * x - trials * np.log(1 + np.exp(x))) -\n",
      "                                  .5 * ((np.sum(np.diff(x)**2) + x[0]**2) / sigma2e + len(dN) * np.log(2*np.pi*sigma2e)))\n",
      "    return joint_loglikelihood, joint_loglikelihood_map"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "values = np.arange(-6., 2., .01)\n",
      "axis = 12\n",
      "joint_loglikelihood, joint_loglikelihood_map = Conditional_on_one_axis(dN, sigma2e, x_map_l, axis, values, trials=trials)\n",
      "plt.plot(values, joint_loglikelihood)\n",
      "plt.plot(values, joint_loglikelihood_map)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 46,
       "text": [
        "[<matplotlib.lines.Line2D at 0x7f1e2c8bee10>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEACAYAAAC3adEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeUVeX1//H3nkbvVUCaCqKgI01QlBFQitQoIkYFAaMx\nJuarxqCYiEnUSH4aC2JFBKWLQRBFJDLYQIqgiAgoFoqANEEQmLJ/f5yDzhAGZrh35s4Mn9dad3Hu\nc855zmYWuuepx9wdERGRQ+JiHYCIiBQuSgwiIpKNEoOIiGSjxCAiItkoMYiISDZKDCIikk1EicHM\nRpjZKjNbbmbTzKx8WN7JzJaY2cdmttjMLspyT3Mz+8TM1pjZI1nKK5nZHDNbbWZvmlmFSGITEZHj\nE2mLYQ5wprsnA2uBO8Py74Hu7n42MBB4Mcs9TwKD3b0R0MjMOoflQ4G57t4YeDtLXSIiUoAiSgzu\nPtfdM8OvC4E6YfnH7r45PF4JlDSzRDOrCZRz98XhPeOA3uFxL2BseDw2S7mIiBSgaI4xDALeOLzQ\nzC4HPnL3NKA2sCHL6Q1hGUANd98CECaV6lGMTUREcinhWBeY2VtAjaxFgAPD3H1meM0wIM3dJxx2\n75nAA8DFxxGb9uoQEYmBYyYGdz/q/9TNbCDQDehwWHkd4BXgGnf/OizeCJyc5bI6YRnAZjOr4e5b\nwi6nrUd5ppKGiMhxcHc71jWRzkrqAvwJ6OnuB7KUVwBeA/7s7guzBLQZ+MHMWpuZAdcCr4anZxAM\nVAMMyFJ+RO5e6D/33HNPzGMoLnEWhRgVp+Is7J/cinSM4XGgLPCWmX1kZqPC8puBU4C/mtmy8FzV\n8NzvgNHAGmCtu88Oyx8ELjaz1UBH4J8RxiYiIsfhmF1JR+Pup+VQfh9wXw7nlgLNjlC+A+gUSTwi\nIhI5rXzORykpKbEOIVeKQpxFIUZQnNGmOGPD8tLvVFiYmRfFuEVEYsnM8PwefBYRkeJHiUFERLJR\nYhARkWyUGEREJBslBhERyUaJQUREslFiEBGRbJQYREQkGyUGERHJRolBRESyUWIQEZFslBhERCQb\nJQYREclGiUFERLJRYhARkWyUGEREJBslBhERyUaJoQhJz8hk3/400jMyYx2KiBRjCbEOoLjJyHA2\nbNvNV5t38M3W7WzcsYNNu7azZfd2tu/bwc79O9iX/iP7M/dyIHMvB30vB20v6baXjLi9ZMbvwy0d\nt3SISwfLCP6MTw8ekBkPcRngBpkJwcfjscwE8AQsM5H4jDLEZ5Yh0cuQSBlKWFlKxpWhZHwZyiSW\no3KpylQrW4WTylfhpEqVqVu1CvVrVOGUk6pQqVzJ2P4ARSTm9M7nPPhu+14+WPUVq9ZvZN333/Ht\nrk1s/nETO9I2sds3sT9xExmlNkNGCRIOViEpvQolqUzZuCqUS6hMpZJVqFyqEuVLlqV8ybJUKF2G\niqXLULFMGSqXLUPlcmWoXK40pZISKZmUQGJCPCUTEyiZlEBSYjwJ8b808NIzMtl/MP3nz4G0dA6m\nZ7B3/0G2797Ljh/3sjP87Nq3lx9++pHd+/eye/8edvy0g50HtrMnfTv7fDv743aQlrCdzBLbITOR\nxP21KJNZiwrxtahWsha1y9WiXuVaNKpZm3NOrU2LU0+mRGJigf/8RSQyuX3nsxJDFgfTMvlg5Xo+\nXPMlKzZ8xZfb17Fx31fs8HXsK/EVnribpJ/qUTbzZCol1KJ6qZOoXb4W9avUotFJtTjj5Fo0a1CT\n8qWL5m/dmZnOhu/3sOLr71i1fhNfbN3Etzs2senHjWw7sIkfMjexL3E9maU2k7C/FuUzGlI9qQH1\nyzekcY0GtGzYkPZnncLJVarG+q8iIkegxHAUO344yFsfreX91av4ZNPnrNuziu99FfvLrCY+rSLl\n0k+lekJDTi7fgMbVGnB2vYa0adyAM+vVJD5OwzK79hxkwapvWfrFV3y6cR1f7lzHpjCB7i+9ljiS\nqHCwCXVKNuGMak1o1aAJnc5uwln1TsbsmP8mRSSfKDEAmZmw+LMtvLZkOR98tZzPdy1ja9xy0st+\nTYmf6lGVJjQo14RmJ53OeY2acHHy6dSoWL4A/gbFV1qas3Dld8xf+TlLvl7F6h2r2HRwFXtKroKk\nPVQ4eAYNSifTsnYyl5x1Dp2Tm1GuZJlYhy1yQjjhEkNGBsxbuoFXFi3kw/Uf8dW+5ewquRxL3E+l\nA8mcWjaZ1iefwyVnn02HsxpTukSJGEV/YnKHT9bsYtbiT3l37XI+27Gc71hGWoVVlDpQj5MTk2le\nO5luyS3o2aI1FUoqQYtEW7FPDF98u5fJ85cy9/OFfLprIdtKfkhc4kFqpp9L0yotaXdKMj1aJXNW\nvbrqvijEtnyfxoyFq5i7YjkfbVrGN+lLSKu6jHLp9WlUug3tT2lD37ZtaFW/CfFx8bEOV6RIK/aJ\ngWGlqXCgKaeXbUP7U86l3/ltOKd+AyWBIs4dvvo2janvfMKczxayYueHbCuxECu3hdrWkrZ1zqd/\n2/Zc0qQtpRNLxzpckSKlQBKDmY0AegAHgC+B69x9t5m1Ap7Jcum97j49vKcL8AjB4rrR7v5gWF4f\nmARUBpYC17h7eg7P9X0Hf6JUYtGc/SN5c/AgvPX+Nia//yHvfv0+38anQvWPqcE5nFerPf3bptDl\njPMok6SxCpGjKajE0Al4290zzeyfgLv7nWZWEjgYltcEPgZOCm9bA3QENgGLgSvd/XMzmwy87O5T\nzexJYLm7P53Dc2OyjkEKh7Q0eHfhXsa/+wGpX83nG+bjNZdxUnwzOjbswMB2nWlXry2J8VprIZJV\ngXclmVlv4DJ3v+aw8gbAB0BtoDVwj7t3Dc8NJUgmD5rZ90CNMJm0AYa7e5ccnqXEID87cADmzt/H\nuLcXMu/bueyo+CZxVb/kzDIp9D2nM1e26kzDSg1jHaZIzMUiMcwAJrn7hPB7a+B5oC5Bt9CrZnYZ\n0NndfxNeczVBsrgXWODujcLyOsDr7n5WDs9SYpAcbdoEL7+xlYmL3uKj3W+SWX8O5UuWo1P9S7im\nTTc6ndKRkgnqhpQTT9QSg5m9BdTIWgQ4MMzdZ4bXDAOau/tlR7i/MTAOuIBgPCKnxLDQ3U8Ly5UY\nJCrS0+H9953Rr33CrNVvsrf2a3j1j2ldrSMD2/ak5+mXUq1MtViHKVIgcpsYjrmJnrtffIwHDQS6\nAR1yuH+1mf0INAU2ErQgDqkDbHT37WZW0czi3D3zUPnRnjt8+PCfj1NSUkhJSTnWX0VOQAkJ0L69\n0b792bifzcqVdzBh+jYmzn2d3344g982vIVG5c/i16168qszetK4auNYhywSNampqaSmpub5vkgH\nn7sADwEXuvv2LOX1gfXunmFm9YD3gbOAH4DVBIPP3wGLyD74/Iq7Tw4Hnz9296dyeK5aDBKxDRtg\nyiv7ee6/81iXOIP4M2ZQtWwFBrS8giubXcEZ1c6IdYgiUVVQs5LWAknAoaSw0N1vCruIhgIHgUyC\n6aqHup26AI/yy3TVf4blDQimq1YClgFXu3taDs9VYpCo+vZbmDQ5k+fnfMj6clOIa/oy1cpX4NoW\nfenX9AqaVGsS6xBFIlbsF7gVxbilaPjiC5g4KZMxby1ke80p+Bkvc1LFivw6+Qr6N+3PaVVOi3WI\nIsdFiUEkQu6wbBmMeSGTF+ctoGybKeypO5lGNeozMPka+jXtR9XS2mJcig4lBpEoOnAAZs2C519I\nJ3X9W1TvNI6tFV6n06kXcc1Z19C9UXdKJGhjRinclBhE8smWLfDSS/Dk87vZ33AaZdq+yPfxH9P3\njMu57pzrOLf2udqzSwolJQaRfOYO8+bBU0/BmwvW07jvS2yu9TwVypRkyDlDuObsa6hcqnKswxT5\nmRKDSAHavBmefx6efsYp3WQ+lTo+y2fps7i00aUMOWcIKfVT1IqQmFNiEImBjAyYPRsefRQ+WbuD\nFte9xFeVniWN/Qw5ZwjXnXMd1ctUj3WYcoJSYhCJsU8/DRLE1Jed9lctIr71M8zb/Ao9GvXg5tY3\n07p261iHKCcYJQaRQuL774NxiFGj4PRzdnBa3zHM3f0EVUpX4eZWN9OvaT9t6icFQolBpJA5cAAm\nTYJ//QuSSmTQ+eY3+ChxJMu++4jB5wzmxpY3Uq9ivViHKcWYEoNIIZWZCTNnwv33ww8/wMDb1rCp\n9ijGf/oiHRt05La2t3FunXNjHaYUQ0oMIoXcoemu998Pa9fCzbftwZNHM+qjR6hTvg63n3c7PRr1\nID4uPtahSjGhxCBShCxeHCSIDz+E2+9Ip3r7V3hsyf9j5/6d3NrmVgYkD6B0YulYhylFnBKDSBG0\nfDkMHw5LlsCddzqnd36Px5c8xAfrP+CGFjfwh3P/oBcLyXHLbWKIK4hgRCR3kpNh+vTgM2uWMajj\nBVy6ZzrzrnmPrXu30nhkY/5v9v+xcfdR32MlEhElBpFCqGVLeP31YBbTlCnQ47xGnL/zaZb/ZgVm\nRrMnm3HDzBtYt3NdrEOVYkiJQaQQa9sW3noLXngBnn4aeqTUpjMPs+b3a6hWphqtn23NNf+5hs++\n/yzWoUoxojEGkSLCPehiGjoU6tWDESOgwek/8MTiJ3j0w0e5oO4F3H3h3STXTI51qFJIaYxBpJgx\ngz59gq02+vSBrl3h97+pwK/r3sW6P6yjXd12dBvfjcumXManWz+NdbhShCkxiBQxiYnw29/CmjXQ\noAE0bw5/+0sZBp3xR774wxe0rdOWjuM60n9af1ZvWx3rcKUIUmIQKaLKlYN77w1aEN9/D6efDlPG\nl+bWNrfzxe+/oFn1ZrQb044B0wfw5Y4vYx2uFCEaYxApJhYvhj/8Idhy4/HHoXVr+GH/D/x74b8Z\nuWgkfU7vw90X3q39mE5gGmMQOcG0agXvvw+/+10wBnHddfDTrgoMTxn+8yym5s8059Y3b2X7vu2x\nDlcKMSUGkWIkLg6uvRY+/xyqV4emTeHhh6F8YmXu73g/K29ayf70/TQe2ZgH3n2AfWn7Yh2yFELq\nShIpxtasCVoQ338frIM4N9y0dc32NQx7exgL1i/g3pR7GZA8gIS4hNgGK/lOeyWJCBCsf5g4EW67\nLehiuv9+qFgxOPfhhg+5Y+4dbNu3jQc6PkCPRj30bupiTIlBRLLZuRPuvBNmzAi6l/r1C9ZGuDuz\n1s5i6NyhVCpViYcueUivHS2mlBhE5IgWLIAbboCTTgpeOdqgQVCekZnBC8tf4C/z/kKnhp14oOMD\n1C5fO7bBSlRpVpKIHFHbtrB0KXToEMxkevzxYIprfFw8g5sPZvXNq6lTvg5nPXUWf5//d35K+ynW\nIUsBU4tB5AS2ejUMGgTx8TB6NJx22i/nvtr5FX96608s3rSYEZ1GcMWZV2j8oYgrkBaDmY0ws1Vm\nttzMpplZ+cPO1zWzPWZ2a5ayLmb2uZmtMbM/Zymvb2YLw/KJZqYpEiL5rHFjeOcd+NWvgpbEww9D\nRkZwrkGlBrx8xcuM6z2Of77/Ty4YcwFLNi2JbcBSICLtSpoDnOnuycBa4M7Dzj8EvH7oi5nFASOB\nzsCZQH8zOz08/SDwkLs3AnYBgyOMTURyIT4e/vjH4LWiM2ZAu3awatUv59vXb8+S65cwMHkgPSb2\nYNCrg9i6d2vsApZ8F1FicPe57p4Zfl0I1Dl0zsx6AeuAlVluaQ2sdfdv3D0NmAT0Cs91AKaFx2OB\nPpHEJiJ5c8op8PbbcM01cMEF8O9/B2MPEIw/DGk+hNU3r6ZiyYo0HdWUJxc/SUZmRmyDlnwRzcHn\nQcAbAGZWBrgDuBfI2p9VG1if5fsGoLaZVQF2ZkkyG4BaUYxNRHIhLg5uuiloPUydCpdcAhs2/HK+\nfInyPNz5Yf577X+Z+OlEzn3uXBZtXBS7gCVfHLMf38zeAmpkLQIcGObuM8NrhgFp7j4hvGY48G93\n35eHwao8jWoNHz785+OUlBRSUlLycruIHMUppwRjDw8+GGzr/eij0L//L+eb1WjG/IHzeemTl+g1\nqRc9GvXggY4PUKV0ldgFLf8jNTWV1NTUPN8X8awkMxsIXA90cPcDYdk7/NKtVAnIAP4KfAQMd/cu\n4XVDAXf3B83se6CGu2eaWRvgHnfvmsMzNStJpIAsWQJXXx0kiCeegEqVsp/ftX8Xf533VyavnMw/\nLvoHg5sPJs40E74wKpAFbmbWhWCA+UJ3P+J2jWZ2D7DH3R82s3hgNdAR+A5YBFzp7p+b2WTgFXef\nbGZPAh+7+1M51KnEIFKA9u2DO+4IBqdfeCFYA3G45ZuXc9Osm8jwDJ7p/gxn1zy7wOOUoyuoBW6P\nA2WBt8zsIzMbdbSL3T0DuJlgNtNKYJK7fx6eHgrcamZrgMrA6AhjE5EoKV0aRo6EZ54JWg9/+Quk\np2e/JrlmMu8Neo/rm1/PxS9ezNC5Q7V7axGlBW4ikiebNwdbe+/fDxMmQJ06R7jmx838cfYfWbxp\nMU93f5pODTsVfKDyP7Qlhojki5o1YfZs6NoVWraEmTOPcE3Zmky6fBKPdXmMwTMGM2D6ALbt21bw\nwcpxUWIQkTyLiwt2ap02DW6+Gf7v/+DAgf+97tJGl7LyppVUKVWFpqOa8tInL6HWfuGnriQRiciO\nHTB4MKxfD5MnB1Ndj2TJpiVcP/N6qpWuxtPdn6ZBpQYFG6ioK0lECkblyvDKK8G4Q9u28NprR76u\nZa2WLBqyiE4NO9Hq2VY8segJMn9e0yqFiVoMIhI1CxbAFVfAwIEwfHiwD9ORfL7tcwa9Ooik+CRG\n9xzNKZVzaGZIVKnFICIFrm3bYEHce+9Bt26wLYfx5tOrns67171Lj0Y9OPe5c3n8w8fVeihE1GIQ\nkahLT4dhw4Ixh6lTgxcC5WT1ttUMmjGIhLgERvcczamVTy24QE8wajGISMwkJAT7LP3733DppcHC\nuJw0rtqYdwa+Q+/GvWnzXBseXfioWg8xphaDiOSrNWugTx+48MJgM76kpJyvXbt9Lde9eh1xFsfY\n3mM1cynK1GIQkUKhUaNgUHrjRrj4Yvj++5yvPa3KacwfOJ9ejXvR+rnWjFk2RuseYkAtBhEpEJmZ\n8Ne/wvjxMH06nH2MPfZWbFnB1f+5moaVGvJM92eoVqZawQRajKnFICKFSlwc/OMf8MAD0KlTsGr6\naJrVaMaiIYtoVLkRZz91Nq+tyWGBhESdWgwiUuCWLg3GHQYNCloRccf4FfWdb95hwPQBXNLwEh7q\n/BBlk8oWTKDFjFoMIlJotWgBixbBnDnBgrh9x9id+8J6F/LxjR9zMPMgyU8ls2D9goIJ9ASlxCAi\nMVGzJsybB6VKwUUXwZYtR7++fInyjOk1hhEXj6DP5D7cm3ovGZkZBRPsCUZdSSISU+5w770wblyw\nz9IZZxz7nk17NnHNf64hPTOdl/q8xMkVTs7/QIsBdSWJSJFgFuyrNHw4pKTAf/977HtqlavFnKvn\n0PXUrrR8tiX/WfWffI7yxKIWg4gUGqmp0K9fMHNp0KDc3bNww0KumnYVnU/pzMOdH6ZUYql8jbEo\nU4tBRIqclBSYPx/uuy/YaykzFztjtKnThmU3LGPXgV20erYVn279NN/jLO6UGESkUDn9dFi4MBiY\nvvZaOHjw2PdUKFmBCb+awO3n3c5FYy/iycVPasV0BNSVJCKF0k8/wZVXwv79wWK4srlcurB622r6\nT+tP/Yr1Gd1zNJVKVcrfQIsQdSWJSJFWqlSQEOrWDaazbt2au/saV23MgsELqFO+Di2eacHSTUvz\nN9BiSIlBRAqthIRgy+6uXaFdO1i3Lnf3lUgowWNdH2PExSPoOr4roxaPUtdSHqgrSUSKhFGjgkHp\nWbMgOTn3963dvpa+U/vSpFoTnun+DOVKlMu/IAs5dSWJSLFy003B+xwuuQTefjv3951W5TQWDF5A\nuaRytHy2JSu2rMi/IIsJJQYRKTIuvxymTAkGpadOzf19pRJL8UyPZ7j7grvpMK4DLyx/Id9iLA7U\nlSQiRc7y5dCtW9C1dN11ebt35daV9J3alzZ12jCy20hKJ5bOnyALIXUliUixlZwcrHO45x54/PG8\n3Xtm9TNZdP0iDmYcpM1zbfhyx5f5E2QRFlFiMLMRZrbKzJab2TQzKx+W1zOzfWb2UfgZleWe5mb2\niZmtMbNHspRXMrM5ZrbazN40swqRxCYixVvjxvDOO8G4w/335+3eskllebHPi9zY8kbOe/483lj7\nRv4EWURF2mKYA5zp7snAWuDOLOe+cPfm4eemLOVPAoPdvRHQyMw6h+VDgbnu3hh4+7C6RET+R/36\n8O67MGECDB0a7NSaW2bGTa1u4pUrXmHIzCH8ff7fyfRc7MFxAogoMbj7XPeff5ILgTpZTv9PP5aZ\n1QTKufvisGgc0Ds87gWMDY/HZikXEcnRSScFm+/NnQs335y7/ZWyOr/u+Sy+fjGzv5xNn8l9+GH/\nD/kSZ1ESzTGGQUDW9lh9M1tqZvPMrF1YVhvYkOWaDWEZQA133wLg7puB6lGMTUSKsapVgymsK1YE\ng9Hp6Xm7v1a5WswbMI+Ty59Mq2dbsXLryvwJtIhIONYFZvYWUCNrEeDAMHefGV4zDEhz9wnhNZuA\nuu6+08yaA9PNLBev38jmqI3C4cOH/3yckpJCSkpKHqsXkeKkfHmYPTt4l3T//kH3UmJi7u9Pik9i\nZLeRjF0+lpSxKYzqNoq+Z/bNv4ALQGpqKqmpqXm+L+LpqmY2ELge6ODuB3K4Zh5wG0HCmOfuTcLy\nK4H27v5bM1sFpLj7lrDL6efrjlCfpquKyBEdOBCsd0hMhEmTICkp73Us3bSUy6ZcRr8z+3Ffx/tI\niDvm79BFQoFMVzWzLsCfgJ5Zk4KZVTWzuPC4IXAqsC7sIvrBzFqbmQHXAq+Gt80ABobHA7KUi4jk\nWokS8PLLkJEBV1yRu227D9eiVguW/GYJS79bSrfx3dj5087oB1qIRdRiMLO1QBKwPSxa6O43mdmv\ngL8BB4FM4K/u/np4TwvgBaAk8Lq73xKWVwamACcD3wBXuPuuHJ6rFoOIHNXBg8Hb4DIyglXSJUrk\nvY70zHRun3M7r699nRn9Z3B61dOjH2gBym2LQSufRaTYSksLts84cCDYwvt4kgPA6I9Gc+d/72Rc\nn3F0ObVLdIMsQEoMIiIEyaF//+DFP9OmQcmSx1fPe9++R9+pfbm97e3c2vZWgt7wokWJQUQklJYG\nv/417NkD//nP8SeHb3/4lp4Te5JcM5mnuj9FyYTjrChGtFeSiEgoMTGYvlqxIvTuHbwu9HjUrVCX\n9we9z48Hf+SisRex+cfN0Q20kFBiEJETQkICvPhisN7heGcrAZRJKsOUvlPockoXWj/bmo+++yi6\ngRYC6koSkRNKWhr07RskikmTgj+P17TPpnHjrBt58tInufyMy6MXZD7RGIOISA4OHAi6lCpXhnHj\nID7++Ota9t0yek7qyc2tbuaO8+8o1IPSSgwiIkfx009w6aXQoAE8+yzERdCxvmH3BnpM7EGLk1rw\n5KVPkhifh704CpAGn0VEjqJUKZg5E9asCXZljeR3zTrl6/Dude+yZe8Wuozvwq79R1ybW2QoMYjI\nCatMGZg1C5YuhVtvjSw5lE0qy/R+02larSltR7dl3c510Qu0gCkxiMgJ7dCurKmpMGxYZHXFx8Xz\naNdH+V2r33H+8+ezYP2CqMRY0DTGICICbNsG7dvDgAFwxx2R1zdrzSwGvjqQkV1H0q9pv8grjAIN\nPouI5NHGjdCuHdx1F1x/feT1fbz5Y3pM7MENLW7grgvuivmMJSUGEZHjsHZt0HJ49NFgvUOkNu3Z\nRPcJ3WlZqyWjLh0V03c7KDGIiByn5cvhkkvgpZeCPyO158Ae+k7tS0JcApMvn0yZpDKRV3ocNF1V\nROQ4JSfDK6/A1VfDgiiMH5crUY6Z/WdStXRVLhp7EVv3bo280nykxCAicgTt2sHYscEK6U8/jby+\nxPhExvQaQ+dTOnPe6PP4YscXkVeaT5QYRERy0LVrMNbQpQusi8KyBDPj7x3+zp/O+xMXjLmARRsX\nRV5pPigeb7gWEcknV14JO3fCxRfDe+/BSSdFXucNLW+gVrlaXDrhUsb0GkP3Rt0jrzSKNPgsIpIL\n990XvDv6nXeCRXHR8OGGD+k9uTd/S/kb17eIwvzYY9CsJBGRKHIP9lRavRpefx2SkqJT79rta+k6\nvitXNbuKe1Puzde1DkoMIiJRlpERrG0oVSp46U8kO7JmtXXvVrqN70bLWi15otsTxMdFsA/4USgx\niIjkg59+CsYb2raFf/0revXuPrCb3pN6U7V0VV7s8yIlEkpEr/KQ1jGIiOSDUqVgxoxgV9ZHHole\nveVLlOf1X79OemY63Sd258eDP0av8jxSYhARyaPKlYMdWR96CKZMiV69JRNKMqXvFOpVqEfHcR3Z\ntm9b9CrPAyUGEZHjULdu0Gq4+eZgy+5oSYhL4Nkez9KhfgcuGHMB639YH73Kc0mJQUTkOJ11Fkye\nDP36wYoV0avXzHig0wMMOWcI7ca04/Ntn0ev8lxQYhARicBFF8Fjj0G3brBhQ3Trvu2827g35V4u\nGnsRizcujm7lR6FZSSIiUfCvf8GECcECuHLlolv3jNUzGDJjCBMvm0jHhh2Pu54CmZVkZiPMbJWZ\nLTezaWZWPsu5s8zsAzP71Mw+NrOksLy5mX1iZmvM7JEs11cyszlmttrM3jSzCpHEJiJSkG6/HVq3\nDrbQSE+Pbt09G/fk5Stepv+0/kz/fHp0Kz+CSLuS5gBnunsysBa4C8DM4oEXgd+4e1MgBUgL73kS\nGOzujYBGZtY5LB8KzHX3xsDbwJ0RxiYiUmDMYOTIICn88Y/BSulourDehbzx6ze48bUbmbhiYnQr\nP0xEicHd57p7Zvh1IVA7PL4E+NjdPw2v2+nubmY1gXLufqizbBzQOzzuBYwNj8dmKRcRKRISE4Pp\nq/PnB+MO0daiVgvmXjuX29+6nec+ei76DwhFc3fVQcChNNYIwMxmA1WBye7+L4LEkXV4ZgO/JJMa\n7r4FwN1GnCu9AAAMfklEQVQ3m1n1KMYmIlIgKlQIprG2bQv160OvXtGtv2n1pqQOSKXTi53Ye3Av\nt7S5JboPIBeJwczeAmpkLQIcGObuM8NrhgFp7n4oMSQA5wMtgf3Af81sCbA7D7EdtSE2fPjwn49T\nUlJISUnJQ9UiIvmnbl149dVgplKdOtCiRXTrP63Kabwz8B06juvI3rS93HXBXUe8LjU1ldTjWGQR\n8awkMxsIXA90cPcDYVk/oIu7Xxd+vxv4CRgPzHP3JmH5lUB7d/+tma0CUtx9S9jl9PN1R3imZiWJ\nSKE3fXqwAO6DD4JkEW3f7fmOTi92olfjXtzX4b5j7sxaULOSugB/AnoeSgqhN4FmZlbSzBKA9sBK\nd98M/GBmrS34G1wLvBreMwMYGB4PyFIuIlIk9e4Nt90G3bvD7rz0l+TSSeVOYv7A+cz+YjZ/nP1H\nMn8e8o1MRC0GM1sLJAHbw6KF7n5TeO4qgllKmcAsd78zLG8BvACUBF5391vC8srAFOBk4BvgCnff\nlcNz1WIQkSLh0HscvvwSXnsNEvLhvZm79u+i2/hunFHtDJ7u/nSO23Zr220RkUIiPR0uvRSaNInu\njqxZ/XjwR3pN6kX1MtUZ13scifGJ/3ONtt0WESkkEhKCPZVmz4Znn82fZ5RNKsusq2ax+8Bu+r3c\nj4MZB4+7LrUYREQKyNq10K5dsNahffv8ecaB9AP0e7kfmZ7J1L5Ts73wRy0GEZFC5rTTYPz4YDfW\ndevy5xklEkowte9UkuKT6DO5D/vT9+e5DiUGEZEC1KkT/OUv0LNn/sxUAkiMT2TiZRMpX6I8PSf2\nZF/avjzdr64kEZEC5g433RRs0z19OsQfeRJRxNIz07nu1evYuHsjM/vPpGyJspqVJCJSWKWlQefO\n0LIljBiRf8/JyMxgyMwh7E/fz6TLJykxiIgUZtu3w7nnBl1LAwbk33MyPZPNP26mdvnaSgwiIoXd\nqlXBDKXp0+G88/L3WZqVJCJSBDRpAmPHwuWXR//VoMdLiUFEJMa6doVbboFf/Qr25312adSpK0lE\npBBwD9Y3lC0Lo0cHb4SLNnUliYgUIWbw/POwZAmMGhXjWIrib95qMYhIcfXll8Eg9NSpcOGF0a1b\nLQYRkSLolFNg3Di48kpYvz42MSgxiIgUMp07B4PRl10Wm8FodSWJiBRC7kGroXTpYOwhGoPR6koS\nESnCDg1GL10KTzxRwM8uir95q8UgIieKdeugbdvoDEarxSAiUgw0bAgvvhiscSiowWglBhGRQu6S\nS4LB6CuugIPH/8bOXFNXkohIEeAOffrAySfD448fXx3qShIRKUbM4IUXYPZsmDAhn59VFH/zVotB\nRE5Un3wCHTvCvHnQtGne7lWLQUSkGDrrLHjooWDxW369M1otBhGRIujGG2HbtmAaa24Xv6nFICJS\njD36KHzzDTz8cPTrVotBRKSI+uab4J3RU6bkbvGbWgwiIsVcvXrBa0H794fvvotevRElBjMbYWar\nzGy5mU0zs/Jh+VVmtszMPgr/zDCzs8JzLczsEzNbY2aPZKmrkpnNMbPVZvammVWI7K8mIlL8de4M\nN9wQrIxOS4tOnZG2GOYAZ7p7MrAWuBPA3Se4+znu3hy4Bljn7p+E94wCBrt7I6CRmXUOy4cCc929\nMfD2obpEROTo7r47eCXonVH6v2ZEicHd57p7Zvh1IVDnCJf1ByYBmFlNoJy7Lw7PjQN6h8e9gLHh\n8dgs5SIichRxccF+SlOnwowZUagv8ip+Ngh44wjl/YCJ4XFtYEOWcxvCMoAa7r4FwN03A9WjGJuI\nSLFWpQpMmgTXXw9ffx1ZXcdMDGb2VjgmcOizIvyzR5ZrhgFp7j7hsHtbA3vd/bPjiE3TjkRE8qBt\nW7jjjmC8IZLN9hKOdYG7X3y082Y2EOgGdDjC6Sv5pbUAsBE4Ocv3OmEZwGYzq+HuW8Iup61He+7w\n4cN/Pk5JSSElJeVol4uInBBuvRXeeQf+/Gfo1SuV1NTUPNcR0ToGM+sCPARc6O7bDztnwHqgnbt/\nnaV8IfAHYDEwC3jM3Web2YPADnd/0Mz+DFRy96E5PFfrGEREcrBzJzRvHix+69Pnl/LcrmOINDGs\nBZKAQ0lhobvfFJ5rDzzg7ucddk8L4AWgJPC6u98SllcGphC0KL4BrnD3XTk8V4lBROQoFi2C7t1h\n4cLgZT9QQIkhVpQYRESO7ZFHYPx4eO89KFFCiUFE5ITnHuzCWqcOPPZY7hPDMQefRUSkaDKD558P\nxhtys5fSIdorSUSkGKtYMdhkb9my3N+jriQRkROEdlcVEZHjosQgIiLZKDGIiEg2SgwiIpKNEoOI\niGSjxCAiItkoMYiISDZKDCIiko0Sg4iIZKPEICIi2SgxiIhINkoMIiKSjRKDiIhko8QgIiLZKDGI\niEg2SgwiIpKNEoOIiGSjxCAiItkoMYiISDZKDCIiko0Sg4iIZKPEICIi2SgxiIhINkoMIiKSTUSJ\nwcxGmNkqM1tuZtPMrHxYnmBmL5jZJ2a20syGZrmni5l9bmZrzOzPWcrrm9nCsHyimSVEEpuIiByf\nSFsMc4Az3T0ZWAvcGZb3BZLc/SygJXCDmdU1szhgJNAZOBPob2anh/c8CDzk7o2AXcDgCGOLudTU\n1FiHkCtFIc6iECMozmhTnLERUWJw97nunhl+XQjUOXQKKGNm8UBp4ACwG2gNrHX3b9w9DZgE9Arv\n6QBMC4/HAn0iia0wKCr/WIpCnEUhRlCc0aY4YyOaYwyDgDfC45eBfcB3wNfA/3P3XUBtYH2WezYA\ntc2sCrAzS5LZANSKYmwiIpJLx+zHN7O3gBpZiwhaBMPcfWZ4zTAgzd0nhNe0BtKBmkAV4F0zm3us\nR+UxdhERyQfm7pFVYDYQuB7o4O4HwrKRwAJ3Hx9+H03QmtgADHf3LmH5UMDd/UEz+x6o4e6ZZtYG\nuMfdu+bwzMiCFhE5Qbn7MX8Jj2jmj5l1Af4EXHgoKYS+JRgzGG9mZYA2wMPA58CpZlaPoJvpyvAD\n8DbBoPVkYADwak7Pzc1fTEREjk9ELQYzWwskAdvDooXuflOYDMYAZ4Tlz7v7w+E9XYBHCcY3Rrv7\nP8PyBgSD0ZWAZcDV4QC1iIgUoIi7kkREpHgp0iufzez34QK7FWb2z1jHczgzu8fMNpjZR+GnS6xj\nOhozu83MMs2scqxjORIz+5uZfWxmy8xstpnVjHVMR5LTws/CxswuN7NPzSzDzJrHOp6scloIW9iY\n2Wgz22Jmn8Q6lpyYWR0zeztcbLzCzP5wrHuKbGIwsxSgB9DM3ZsB/y+2EeXoYXdvHn5mxzqYnJhZ\nHeBi4JtYx3IUI9z9bHc/B5gF3BPrgHKQ08LPwmYFwXqh+bEOJKtjLIQtbMYQxFmYpQO3uvuZQFvg\nd8f6eRbZxAD8Fvinu6cDuPu2GMeTk6IyUP5vgokEhZa7/5jlaxkgM6drY+koCz8LFXdf7e5rKXz/\nRo+2ELZQcff3gJ2xjuNo3H2zuy8Pj38EVhGsKctRUU4MjYALw/2V5plZy1gHlIPfhV0Kz5lZhVgH\ncyRm1hNY7+4rYh3LsZjZP8zsW+Aq4K+xjicXsi78lNw54kLYGMVSrJhZfSAZ+PBo1xXqjeqOsrju\nboLYK7l7GzNrBUwBGhaiGIcBo4C/ubub2T8IpuzGZA+oY/ws7yLoRsp6LiaOtaDS3e8G7g77nX8P\nDC/4KI974WeBy02ccmIws7IEu1Lccljr+38U6sTg7hfndM7MbgReCa9bHA6aVnH37Tndkx+OFuNh\nngVi9h9iTnGaWVOgPvCxmRlBt8dSM2vt7lsLMEQgTz/PCcDrxCgxHCvOcOFnN4L1PDGTh59nYbIR\nqJvle52wTI5TuFv1y8CL7p7jGrFDinJX0nTC/+jMrBGQWNBJ4VgOmzXzK+DTWMWSE3f/1N1runtD\nd29A0Gw/JxZJ4VjM7NQsX3sT9JUWOlkWfvY8bOFnYVaYxhkWEy6ENbMkgkWwM2Ic09EYhevndyTP\nA5+5+6O5ubjIrmMws0SCv2wywe6tt7l7YZtdMY4gvkyCzQRvcPctMQ3qGMxsHdDS3XfEOpbDmdnL\nBGNLmQSzp2509+9iG9X/ymnhZwxDOiIz6w08DlQl2Op+eU7b0BS0nBbCFjZmNgFIIdgTbgvBVj5j\nYhrUYczsfOAdglloHn7uOtosySKbGEREJH8U5a4kERHJB0oMIiKSjRKDiIhko8QgIiLZKDGIiEg2\nSgwiIpKNEoOIiGSjxCAiItn8f/wYin/pwVVmAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f1e2c8485d0>"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Expectation-maximization to optimize $\\sigma_{\\epsilon}^2$ using the Laplace approximation at the MAP"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def cov_from_chol_precision(U):\n",
      "    \"\"\"Given the Cholesky factorization (U) of the posterior precision matrix (J), with U^t * U = J,\n",
      "    return the tridiagonal part of the covariance matrix.\n",
      "\n",
      "    :param U: Cholesky factorization (U) of J, given as [0, A; D] where A is the upper diagonal and D the main diagonal\n",
      "    :return: Cov_tri: Tridiagonal part of the covariance matrix returned as [0, C_i,i+1; C_ii; C_i+1,i, 0]\n",
      "    \"\"\"\n",
      "    assert(U.shape[0] == 2 and U[0,0] == 0)\n",
      "    A, D = U # Unpack matrix into first (above) diagonal and diagonal\n",
      "    Cov_tri = np.zeros_like(U)\n",
      "    C, V = Cov_tri # Obtain _views_ into the first diagonal and diagonal\n",
      "    # Compute last element of diagonal\n",
      "    V[-1] = 1. / (D[-1] ** 2)\n",
      "    # Recursively compute other elements of main diagonal and first diagonal\n",
      "    for i in range(len(D)-1, 0, -1):\n",
      "        iD = 1. / D[i-1]\n",
      "        iDA = iD * A[i]\n",
      "        N = -iDA * V[i]\n",
      "        C[i] = N\n",
      "        V[i-1] = iD ** 2 - N * iDA\n",
      "    return Cov_tri"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def EM_fit_sigma2e_Laplace(dN, sigma2e_init, tol=1e-6, trials=1.):\n",
      "    \"\"\"Optimize sigma2e using the EM algorithm using the Laplace approximation at the MAP.\n",
      "\n",
      "    :param dN: Observations (K,)\n",
      "    :param sigma2e_init: Initial estimate of sigma2e\n",
      "    :param tol: Convergence criterion for the EM\n",
      "    :param trials: Number of trials for binomial observations (1 for Bernoulli)\n",
      "    :return: x_map, U, marginal_loglikelihood, joint_loglikelihood\n",
      "    \"\"\"\n",
      "    sigma2e_old = sigma2e_init\n",
      "    while True:\n",
      "        x_map, U, marginal_loglikelihood, _ = Do_Kalman_Likelihood_Bernoulli_LaplaceMAP(dN, sigma2e_old, trials=trials)\n",
      "        #print(sigma2e, marginal_loglikelihood)\n",
      "        Cov_tri = cov_from_chol_precision(U)\n",
      "        sigma2e = (np.sum(Cov_tri[1]) + np.dot(x_map, x_map) # E[x_k^2]\n",
      "                   + np.sum(Cov_tri[1,:-1]) + np.dot(x_map[:-1], x_map[:-1]) # E[x_{k-1}^2]\n",
      "                   - 2 * np.sum(Cov_tri[0]) - 2 * np.dot(x_map[1:], x_map[:-1])) / len(dN) # E[x_{k-1} * x_k]\n",
      "        if (abs(sigma2e - sigma2e_old) < tol): break\n",
      "        sigma2e_old = sigma2e\n",
      "    return sigma2e"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 86
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%time sigma2e_opt = EM_fit_sigma2e_Laplace(dN, .0019, trials=trials)\n",
      "print(sigma2e_opt)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 1.11 s, sys: 4 ms, total: 1.12 s\n",
        "Wall time: 1.11 s\n",
        "0.104697603827\n"
       ]
      }
     ],
     "prompt_number": 87
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}