From 8d1a93b2165a8b2bba26016b6ea1db1743b4a1a9 Mon Sep 17 00:00:00 2001 From: rcw5890 Date: Thu, 4 Feb 2021 13:53:00 +0000 Subject: [PATCH] Add broken MH notebook --- .../notebooks/broken_MCMC_MH_1DGaussian.ipynb | 264 ++++++++++++++++++ 1 file changed, 264 insertions(+) create mode 100644 change_detection/notebooks/broken_MCMC_MH_1DGaussian.ipynb diff --git a/change_detection/notebooks/broken_MCMC_MH_1DGaussian.ipynb b/change_detection/notebooks/broken_MCMC_MH_1DGaussian.ipynb new file mode 100644 index 0000000..65d804d --- /dev/null +++ b/change_detection/notebooks/broken_MCMC_MH_1DGaussian.ipynb @@ -0,0 +1,264 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Detection of Broken MCMC Samplers - MH/1D Gaussian\n", + "In this notebook, we deliberately break a standard Metropolis Hastings sampler by causing it to incorrectly accept proposals with a certain probability. The sampler is used on the 1D Gaussian distribution." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Dependencies\n", + "This notebook uses hildensia's implementation of the ''Bayesian Online Changepoint Detection'' algorithm by Ryan Adams and David MacKay, available at the following link\n", + "\n", + "https://github.com/hildensia/bayesian_changepoint_detection" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from functools import partial\n", + "import random\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import pints\n", + "import pints.toy\n", + "import pints.plot\n", + "import bayesian_changepoint_detection.online_changepoint_detection as oncd" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "class BrokenMH(pints.MetropolisRandomWalkMCMC):\n", + " \"\"\"Broken version of Metropolis Hastings.\n", + "\n", + " At each MH step, with probability given by error_freq, it will always\n", + " accept the proposal.\n", + " \"\"\"\n", + " def __init__(self, x0, sigma0=None):\n", + " super().__init__(x0, sigma0)\n", + " self.error_freq = 0.0\n", + "\n", + " def set_error_freq(self, error_freq):\n", + " self.error_freq = error_freq\n", + "\n", + " def tell(self, fx):\n", + " if self.error_freq == 0.0 or random.random() > self.error_freq:\n", + " # Run MH step correctly\n", + " return super().tell(fx)\n", + " else:\n", + " # Always accept it even if it is bad\n", + " self._acceptance = ((self._iterations * self._acceptance + 1) /\n", + " (self._iterations + 1))\n", + " self._iterations += 1\n", + " self._current = self._proposed\n", + " self._current_log_pdf = fx\n", + " self._proposed = None\n", + " return self._current" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "class NormalDist(pints.toy.GaussianLogPDF):\n", + " \"\"\"Same as Pints, except it doesn't complain about KLD for 1d distribution.\n", + "\n", + " # todo: check if there is a problem in pints\n", + " \"\"\"\n", + " def kl_divergence(self, samples):\n", + " \"\"\"\n", + " Calculates the Kullback-Leibler divergence between a given list of\n", + " samples and the distribution underlying this LogPDF.\n", + "\n", + " The returned value is (near) zero for perfect sampling, and then\n", + " increases as the error gets larger.\n", + "\n", + " See: https://en.wikipedia.org/wiki/Kullback-Leibler_divergence\n", + " \"\"\"\n", + " m0 = np.mean(samples, axis=0)\n", + " m1 = self._mean\n", + " s0 = np.cov(samples.T)\n", + " s1 = self._sigma\n", + " cov_inv = np.linalg.inv(s1)\n", + "\n", + " if s0.ndim < 2:\n", + " s0 = np.array([[s0]])\n", + "\n", + " dkl1 = np.trace(cov_inv.dot(s0))\n", + " dkl2 = np.dot((m1 - m0).T, cov_inv).dot(m1 - m0)\n", + " dkl3 = np.log(np.linalg.det(s1) / np.linalg.det(s0))\n", + " return 0.5 * (dkl1 + dkl2 + dkl3 - self._n_parameters)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "def get_klds(num_runs, start_breaking, error_freq):\n", + " \"\"\"Run MCMC multiple times, break at some point, and get the KL divs.\n", + " \n", + " It uses the BrokenMH sampler and the 1D Gaussian distribution.\n", + "\n", + " Parameters\n", + " ----------\n", + " num_runs : int\n", + " Total number of runs\n", + " start_breaking : int\n", + " Which run to break the MCMC algorithm\n", + " error_freq : float\n", + " Error probability per MH step once the algorithm is broken\n", + "\n", + " Returns\n", + " -------\n", + " list\n", + " List of kl divergences from samples to posterior for each run\n", + " \"\"\"\n", + " posterior = NormalDist(np.array([1.0]), np.array([0.1**2]))\n", + " x0 = [np.array([1.0])]\n", + "\n", + " klds = []\n", + " for run in range(num_runs):\n", + " mcmc = pints.MCMCController(posterior, 1, x0, method=BrokenMH)\n", + " mcmc.set_max_iterations(1000)\n", + " if run >= start_breaking:\n", + " for s in mcmc.samplers():\n", + " s.set_error_freq(error_freq)\n", + "\n", + " mcmc.set_log_to_screen(False)\n", + " chains = mcmc.run()\n", + " kld = posterior.kl_divergence(chains[0])\n", + " klds.append(kld)\n", + "\n", + " return klds" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "def run_changepoint(data):\n", + " \"\"\"Run changepoint detection using the library.\n", + " \"\"\"\n", + " ## Set hyperparameters\n", + " # timescale of hazard function\n", + " lam = 250\n", + "\n", + " # T distribution parameters\n", + " # df=2*self.alpha,\n", + " # loc=self.mu,\n", + " # scale=np.sqrt(self.beta * (self.kappa+1) / (self.alpha * self.kappa))\n", + " alpha = 0.1\n", + " beta = 0.01\n", + " kappa = 1.0\n", + " mu = 0.0\n", + "\n", + " R, maxes = oncd.online_changepoint_detection(\n", + " data,\n", + " partial(oncd.constant_hazard, lam),\n", + " oncd.StudentT(alpha, beta, kappa, mu))\n", + "\n", + " return R" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We run the following experiment to test broken MCMC sampling and changepoint detection.\n", + "\n", + "The MCMC sampler is run 60 times. From the 30th run onwards, the sampler is deliberately broken by insisting that at each accept/reject step, the proposal will be automatically accepted with probability 0.15. \n", + "\n", + "Next, we attempt to detect any changepoints in the series of KL divergences obtained." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXibZ5Xw/++RZFleJC+xEzt7mqRJ0zZJm7SltBQoy7QFWpYCLcsAw0xZ2nlnBpiZ9je8MDDvXC/LOwydoQNToKwDZYdQWkpb2lLokqVJ2uxJs3uP992WdH5/PI9s2ZZt2ZYsyT6f69Jl63keyfdtyzq6t3OLqmKMMWb+8mS6AMYYYzLLAoExxsxzFgiMMWaes0BgjDHznAUCY4yZ53yZLsBUVVRU6MqVKzNdDGOMySm7du06p6qVic7lXCBYuXIlO3fuzHQxjDEmp4jIqfHOWdeQMcbMcxYIjDFmnrNAYIwxGbL7dCvv+vqz9IcjGS1HWgOBiFwnIodF5JiI3Jng/IdF5EUR2SMifxSRDeksjzHGZJNnjjfz9EvNNLT3Z7QcaQsEIuIF7gGuBzYAtyZ4o/+Bql6sqpuBLwBfSld5jDEm2zR3DQDQ0TeY0XKks0VwOXBMVY+r6gBwP3BT/AWq2hF3twiwDHjGmHmjpdsNBL2ZDQTpnD66BDgTd/8scMXoi0TkduBjgB+4NtETichtwG0Ay5cvT3lBjTEmE851OV1Cc7lFkBRVvUdVVwP/CHxynGvuVdWtqrq1sjLheghjjMk5sa6h9gy3CNIZCGqAZXH3l7rHxnM/8OY0lscYY7LKcNdQOKPlSGcg2AGsFZFVIuIHbgG2xV8gImvj7r4BOJrG8hhjTNZQVZq7s6NrKG1jBKoaFpE7gIcBL3Cfqu4Xkc8CO1V1G3CHiLwWGARagfelqzzGGJNNOvvDDEac+TFzebAYVX0QeHDUsU/Fff836fz5xhiTrWLjAwAdfXO3a8gYY8w4mruGF5HN5cFiY4wx42h2B4pDAV/Gu4YsEBhjTAbEuoZWVRZnfLDYAoExxmRArGto5YLCOT191Bhj5qVf7anhV3smWjbldA0FAz4qivOtRWCMMXPJnjNtfOzHe/nqEy9NeF1z9wAVxfmUFOTRMxBhMBKdpRKOZYHAGGNSpLs/zN/ev5tIVKlr75vw2pbufsqL/IQCziz+TA4YWyAwxpgU+ZcHDnCqpYdr1y+kvXeQnoHx+/6buwZYUOQnVJAHZHYtgQUCY4xJgd/uq+f+HWf4yCtX86ZN1QDUT9AqONc1wIJiP6GAGwisRWCMMbmroaOPO3/+AhcvKeFvX3s+VaECYPxAEI0qrT0DLCjKj2sRWCAwxpicFI0qn/jJXvoHo3z5ls34fR6qSwIA1I4TCNp7B4lElQXFfkpigSCDU0gtEBhjzAx86+mTPHX0HJ984wWsriwGoMoNBPXtvQkfE8s6Wl7kJ1TgDBZnMs2EBQJj5qHnjjfz9q89TX84kumi5DRV5Z7Hj/GKtRW86/Lh3RMDeV7Ki/zjzhyKrSquKM4fHiOwriFjzGx65ngzO062cqYl8SdWk5z6jj5augd43YZFiMiIc1WhwLhjBLE8Q+VFfgr9XrwescFiY8zsaup0uibqxum6MMk5UNsBwAXVoTHnqksC47cI3ECwoNiPiDiJ56xFYIyZTY2xQNA28aInM7GDdU4gWF8VHHOuqiRAfcd4XUPuGEGhH4CSgjwbLDbGzK5Yi6DWWgQzcrCuk+XlhQTdfv541SUBWroH6BscOw7T3DVAaWEePq/zFhwqyLMWgTFmdjVZiyAlDtZ1cEH12NYAQFXJ+GsJWrqdVcUxoUCezRoyxsweVbUWQQr0DIQ50dzNhuqShOcXu1NIE40TnOvqZ0FR/tD9UEFmN6exQGDMPNPeO8iAm+lyohQIZmKH6jtRZYIWgbuWoGNssG3udtJLxIQCeXM315CIXCcih0XkmIjcmeD8x0TkgIi8ICKPiciKdJbHGDPcLVRRPP48dzO52EBxohlDMBwIEv2OW0YFAmeweA62CETEC9wDXA9sAG4VkQ2jLtsNbFXVjcBPgS+kqzzGGEcsEGxaWkpXfzjjm6LkqgO1HQQDPpaWFSQ8X+j3UVKQN2YcJuLmGSof0TWUR384mnBgeTaks0VwOXBMVY+r6gBwP3BT/AWq+riq9rh3nwWWprE8xhiGp45uXFoK2IDxdDkDxaExC8niJVpL0NozgKrTIosZ2pMgQ0E5nYFgCXAm7v5Z99h4Pgg8lOiEiNwmIjtFZGdTU1MKi2jM/BNrEWxc5gxy2oDx1EWjyqH6TjaM0y0UU10SGDNGEEsvUR4/ayjDieeyYrBYRN4DbAW+mOi8qt6rqltVdWtlZeXsFs6YOaaxs49AnofzFzmDnNYimLrTLT30DETGHSiOqSopGDMgH1tMNmLWUIbzDfnS+Nw1wLK4+0vdYyOIyGuBfwJeqar9aSyPMQanRVAZzGdRMB+PWJqJ6TjgDhSPN3U0prokwLmuAfrDEfJ9XmA4vcSIrqGCzG5Ok84WwQ5grYisEhE/cAuwLf4CEbkE+G/gRlVtTGNZjDGuxs5+FgYD+LweFgYD1FqLYMoO1nXg9QhrFxVPeF1s5lBD+/Bn3KH0EkXxs4ZiYwRzrGtIVcPAHcDDwEHgx6q6X0Q+KyI3upd9ESgGfiIie0Rk2zhPZ4xJkabOfiqLnW6J6tKAtQim4WBdB+dVFBHI8054XfXQFNLh33Fz9wAegdLCkesIIHMtgqS7hkSkEkBVkx6tVdUHgQdHHftU3PevTfa5jDGp0djZz5WrFwCwuKRgqJvDJO9AbQeXrSqf9LrqoUVlw62u5u4Bygr9eD3Ds41iXUOZSjMxYYtAHP8sIueAw8AREWkSkU9N9DhjzEi7TrXwpUeOZLoY9IcjtPcODrcISgLUtvWiqhkuWerVtvVyrLEr5c/b1jNAbXvfuAvJ4sXyDcVPIW3u6h+xmAwg3+fB7/Vk7fTRvwOuAi5T1XJVLQOuAK4Skb9Le+mMmSN+uuss//HYUdp7Mrt4KzZ1tDIY6xoqoD8cpTXD5UqHf962n/d/a3vKg9zBuk5g/BXF8YrzfQQDvhEzh5yEc/kjrhMRJwNplk4ffS9wq6qeiB1Q1ePAe4A/T2fBjJlLzrY6fcQH66feDdMfjnDzV5/mj0fPzbgcsUCwMOS8EcUSo9W2zb1xglPNPZxt7eVUc8/kF0/BwaEZQ5MHAhhudcU0dw1QPqpFAG7iuSxtEeSp6phXnztOMDYBtzEmoRr3jeDgNPrjD9d3svNUK7/eWzvjcgy1CIqdAFBdOrbrYi5Q1aHf+dMvNaf0uQ/UdVBRnD/UqppMVUnBiDGCc139VBQlCASBzOUbmiwQDEzznDHGpapDnwinEwj2u9sh7jrdOuOyNI7TIphrM4c6esN09TvdLH96aeYtqXgT7UGQSHVoOM3EQDhKR194RJ6hmFAGE89NFgg2iUhHglsncPFsFNCYXOfsUuWkfZ7ODJ39te0AHGvsmvEYQ1NnPyIMbYpSUZxPnlfmXIvgbJvTHVRSkMczLzUTjaZmnGAwEuVoQxcbFifXLQTOWoJzXf0MhKO09gzvVTyas29xFo4RqKpXVUMJbkFVta4hY5IQ66JYVVHEkYYuwu5eAMnaX9tBod+Zr/78mZm1Cho7+1lQ5B/aItHjERaFAtTNsTGCGndM5qbNi2npHuBQfWdKnvd4UzcDkWjS4wMAi0sDqDqpPWJ5hioSBIJMpqKe9oIyETmdyoIYM1fF3pRee8FCBsJRjp/rTvqxkahyqK6TGzctxusRnj81s0DQ1NlPRfHIbonFJQXUzrEWQawr7u1bnCw3T0/QPdQ3GOHRAw1JPe+BOqd1lsyMoZj4LSubu2OrisfpGuobzMhU3pmsLB4/96oxZkisRfDaCxYBUxsnOHGum97BCFtXlrO+KsiuGQeCPhaGAiOOzcXVxTVtveT7PFy0JMR5FUUTDhj/95PH+cvv7uRA7eR/l4N1nfh9Hs6rKEq6LLFFZbXtwy2CxF1DeQxGdKgbcTbNJBDMvRUoxqRBTVsvhX4vl64ow+/1JPWGExMbH7hwcYgtK8rYe6Ztyl1L8eLTS8RUuxkyU9WPng1q2npZUlqAiPDyNQt47ngzgwl+b+FIlPt3OJ0bu5PodjtY18G6RcGhrrVkDG1Z2d47lHBuQaJZQ26+oUysLp5sZfHHxrl9HCdHkDFmEjWtzptSntfD2kXFUxowPlDXgd/rYc3CYrasKKN7IMLhhun1d6sqTV39QzOGYqpLAgxGlHPdcyf5b01rL0vcncOuWl1B90CEF862jbnuicNNQwPle8+MPR9PVTlQO7UZQwDBfB9Ffi917X00d/Xj88hQbqF4mUxFPVlYC45zKwbuTm/RjJkbatqG35QuqA4NrUxNxoHaDs6vKibP6+HS5WUA0x4naOsZZDCiCVoE7hTSOZSFNNYiAHjZeQsQgaePje0e+sH201QG83nF2gr2nmmf8Dlr2/to7h7gwsUTp54eTUSoKgk4YwRdA5QX+fF4xvasl2QwFfVkgeA+Vf1MohuwazYKaEyuq417U7qgOsS5rv6hhV0TUVX213ZwoZvzfmlZAZXB/AnHCZ493sxzxxP3hzeOSi8Rs3hoUdncGCfoG4xwrmtgqF5lRX42VIfGrCeoaevlicONvHPrMrasKONIY+fQ2oNEtp9wfq+XrZw82dxoi0sLnBZB98CI9NPxhvYkyMIWwSMisnL0QRH5ANYiMGZSPQNhWnsGh96UYt0KyQwY13f00dI9wIVLnBkqIsKW5WU8fzpxF8ZAOModP9jNP/1yX8LzQ+klgolbBLmyL8Fks2piM4ZiwRfgqjUVPH+qjd6B4c3hf7T9NArccvkyNi0rRRVePDt+q2D7iVZCAR/rqqbWNQRQFXIG5Ju7x87aihnatzgD+YYmCwQfA34nImtjB0TkLvf4K9NZMGPmgtjU0aVu11Bs/nky4wSxQeX4OetbVpRxuqWHxs6xb9oP76/nXFc/LzV1Jfxk29TlPGZ0i6C8yE++z5MTLYLdp1u56nO/52e7zo57TWyWVqw7DuDlqxcwEImy81QL4AwS/2jnGV55fiVLywrZtLQUgL0JxhFitp9o5rKV5SPSRyeruiRAY2c/jR39udcicPcT+AjwkIhcJCJfBt4EXKOq4/8ljDEAnB316bS00M/ikkBSLYL9tR2IjJyzfumK2DjB2Des7z17Cp9Hxv1k29gRSy8xcvqoiDiJ0bJ8LcFTR5t49zeeo7a9jyeOjL8tSiz4xrcILl9Vjs8jQ9NIHzvUSENHP++6fDngBMPl5YXjDhg7AbY7qT0IEqkqKUDVCVKJpo7C8GBxJjLUTjoHSlUfAz4APAGcB1yrqjNPemLMPFCb4NOpM2CcTCBoZ9WCIoryh/ePumhJCL/Xw+5ReYeONHSy/UQLH3zFKiDxJ9umzn4K8rwU+cfuqlVdUpDVq4t/80Idf/HtHSwvL2TLijIOT5DFtaatF48MT9sEKPT7uGR5KU8fc8YJfvDcaapCAa5dv3Doms3LSscNBDtOOC2Jy6cZCKrjyjJe15Df56Egz5t9LQIR6RSRDpxdxkLAa4DGuOPGmAnUtPbi8wgLg8NvBBdUh3ipqZu+wcgEj3RaBKNz2uT7vFy0JDRmwPj7z57C7/PwoWtWs6y8IOFUycZOZ+qoyNiuDWdRWXa2CH7w3Gnu+OHzbFpayo8+dCVXrCp3Uj2EE6+nqGntZVEoQN6ouf4vX13BizXt7K9t5w9Hm3jHZctGrAfYtKyU2vY+GjvG/h62n2yhIM/LRVOcMRQTH5TG6xoCNxV1to0RuDmFYrmF/KpaFHc/+TXWxsxTNW29VJUERvQrX1AdIhJVjjaMv3tWe+8gZ1t7EyY327KijBdq2ukPO4Gkqz/Mz5+v4Y0bqykv8rNpaWnCqZCJFpPFLC4poKGjb0aL1dLhnseP8f/94kVedX4l3/vgFZQU5LGuKkg4qhw/l/j3Fz91NN5VayqIKnziJy8gwC2XLRtxfvMy501+T4JWwfYTLVy6ohS/b3prcBeXDJcn0WKymFAgL/taBMaYmYktJosXe3OfqHsoNlCcaM76lhVlDISjQ+mpf7m7hq7+MO952QoANi0tpaatd8wU1cbOvjGLyWKqSwNEdXiKaTZ4+tg5vvjwYW7ctJh7/3wrBW6XVmzWzuFxEsnFr9uIt3lZKQV5Xg7WdfDqdQuHZnLFXLi4BK9HxnSrdfQNcqCuY1rTRmNCBT4K3I3uxxsjcK6bg4FARK4TkcMickxE7kxw/hoReV5EwiJyczrLYkwm1CZ4U1pRXkih3zvhzKFYaolEWS7jF5apKt9/9hQXLg5xyTJn5ssm9+vo7qHJWgSQXRvU3P3YURYG8/nCzRtHdPOcV1GMzyMJM4pGokp9e1/CFoHf5xnq43/XFcvHnA/keVlfFRzTmtp1shXV6Y8PwPCAPDBmm8p4oYBv3BQTE61xmKm0BQIR8QL3ANcDG4BbRWTDqMtOA+8HfpCuchiTKYORKPUdY9+UPB5hXVVw4hZBXQcLg4l3wVoYCrC0rIDnT7ey61Qrh+o7ee/LVgz1/V+0JIRHRqZM6BuM0NEXHjNjKKa6NLs2qNl+ooXnTrTwoVeuJpA3cnDb7/OwurI4YYugoaOPcFQTtggA3nnZMl61rpJXrVuY8PzmZaXsPds2Iu/S9pMt5HmFS5aVzaBGw+MEE7UISsbZt7i9Z5BNn/kd33v21IzKMJ50tgguB46p6nFVHQDuB26Kv0BVT6rqC0B2dUwakwL17X1ElYSfTi+oDnGgrmPcxVEHaju4cILNT7asKGPXqVa++8wpggEfN25ePHSu0O/j/EVB9sZNIR3eonKcrqFYiyBLFpX95++PsqDIPzS9c7R1VcGEgaAmwWKyeDdcXM23P3D5uGsBNi0rpbMvzInm4VTh20+0sHFp6VDX1HRVlQTwez0Ux80CG228rqGdp1qIRJU1lelJ8ZbOQLAEOBN3/6x7bMpE5DYR2SkiO5uaxp8/bEw2SbSwKeaC6hCdfeGha+L1DUY42tg1YU6bLSvKaOjo5zcv1nHzlqUU+ke+uWxa6nyyjQWapq7E6SViQgEnMVptFrQI9pxp46mj5/jLV5w37pvvuqogNW29dI5600y0hmAqNrvdanvc1du9brK6mYwPxNx6+XI+8WfnJ5y1FRPbt3j0B4ShVsny0hmXI5GcGCxW1XtVdauqbq2srMx0cYxJSqJUBzEbhlJNjP1Ue6Shk0hUJ2wRxMYJIlEdGiSOt3FZCW09g5xpccoQW0w2XiAQEapLC7KiRfCV3x+ltDCP9145tl4x6xY5v78jozKxThR8k7G6spgiv3dowHj3mVYGI8oVMxgfiLlsZTm3XbN6wmtCBT6iCt0DI6cW7zjRwsVLSsZ0k6VKOgNBDRA/P2upe8yYeSH26XT07BSAdVUhRBLPHBpKLTFBIFhfFaTI7+XlqxewOkF3QSxlwh73DS3WIhidZyhedUnmN6jZX9vOowcb+YurVk3YhTI8c2jkFNKatl7KCvPGtJCS5fUIFy8tGRpf2X6iBRHYsnJm4wPJGkpFHTdg3DcY4cWa9mmvak5GOgPBDmCtiKwSET9wC7AtjT/PmKxS09ZLRbE/4ae44nwfK8oLEwaC/bUdBPN9LCsrHPe5fV4P973/Mj7/to0Jz6+rCpLv8/CC+4bW1NGHR2DBOGMEQFakmbjn8WME83287+UrJ7xuSWkBRX7vmBXG8fsQTNfmZWUcqOugPxxhx8kWNlSHEu4fkA6xVNTxM4d2n25jMKJcnoLuqfFML2wmQVXDInIH8DDgxUlpvV9EPgvsVNVtInIZ8AugDHiTiHxGVS9MV5mMmU3jLWyKiQ0Yj7a/tp0LFocS5qyPd8V5C8Y9l+f1cOHi0FAXR1NXP+VF+RMmTKsuKeBcVz8D4ei0F06NNhCOcrCug+7+MF39YboHwnT3R/D7PFy7fuGIdAtHGzp5aF89t79qzdAb4ng8HuH8quCYKaQ1bb2srkx+G8lENi8rYTCivHC2nV2nWrnlssQD1ukQSrAnwY6TTqtk64ocDAQwlLTuwVHHPhX3/Q6cLiNj5pyatl7WT5Cy+ILqEA/tq+eft+1nQ3WIC6pDrF5YxKH6Tt6xddm4j0vWpmWl3L/9DOFIlMaO/gm7hQAWlwZQdaZgLisfvzWSrK7+MO/++rMjZi/F84iT9uGNG6u57qIq7nn8GAV5Xv7i6lVJPf/6qiAP7atHVRERVJWa1l6uWTuzccTYOozvP3uKvsFoSsYHkjW8S9nwFNIdJ1tYtyhISWH6WiVpDQTGzFeqSm1bL9eOM18d4IaLq3jicCM/2nGGXjfvkAioMuFAcbI2LS3lW386ydHGLpq6+scdKI6JjTV87qFD/L+3b5rRdMmBcJQPf28X+2o7+JebLmTNwiDF+T6K8r0U5/s41zXAb16s5YEX6rjz5y/yyV/uI6LKX73ivAlz8cRbtyjID7efoamzn4WhAG09g/QORlhcmnitRLKqQgEWBvN54IU6gLT2zY8W27c41iIIR6I8f6qVt16a3s/LFgiMSYPm7gH6BqMT9levWRjk5x+9ikhUOd3Sw8G6Dg7VdVDX3sdrL1g04zLEPtnuPdNGY0f/0Eyb8WxZUcY/XreeLzx8iJPN3dz751unNQ0zGlU+/pO9/PHYOb5480benqB1szAUYMPiEJ94/Tr21XTwwAu1HKjr4LZrzkv655zvtrYO1XeyMBQYmjG0dIZjBCLC5mWl/O5AA6sri8bNFpoOo/ctPlDXQfdAJO3ByAKBMWkw0dTR0bweYVVFEasqirjh4uqUlWHlgkJCAR97zrRxLokWgYjwkVetZl1VMX/zwz3c9JU/8rX3bGHrFAYpVZXPPnCAX++t5c7r1ycMAqN/5sVLS7h46dSzeq6vclpNh+s7ueb8Ss4OrSGYebfWJjcQzCStxHQE3V3KYoPF22Ppr9M4UAw5so7AmFwztLBphp9OZ0JE2LSslCePNBGO6qRjBDHXrl/EL25/OcX5Pm79+rP8aMdpBpPMSnrP48f49tMn+atXrOJDU/h0Px3lRX4qg/kcdtcSzHQNQbzYwq0rVo0/IJ8OPnflcSzNxI6TLSwrLxiRxjotPzetz27MPDVZqoPZsmlpKU8ddTZjqQwm/2ayZmGQX91+NXf88Hn+8Wcv8o8/exGfRyjweynI81Lgd/r6Q4E8SgryCBX4iEThZ8+f5a2XLOGu6y+YcAVtqqyPSzVR09pLQZ6XshQMql553gK+9p5LU9JFN1WhgI+OPmd18c6TrbxyXfoX0VogMCYNzrb2UuT3TjoNMt02xnW5TNY1NFpJYR7fev9l/GJ3DfXtffQORugdjNA3GKFnIEJXX5iOvkFOnOumvXeQjr5B3rixms/fvHHSqa+pcv6iIN9/9hSRqFLT1sPi0kBKApCIcN1Fqeumm4pQgZNm4qWmbpq7B9LeLQQWCIxJi1j66dn4VDyRWO4cmHhV8Xh8Xs+k/fyZtK4qSH84yqnmbmrb+lgywSK8XBHbnGbHSWd8YDZmLdkYgTFpMNlistmyMBSgyk09PdUWQS6IrdM40tCZNb/zmQq5qah3nGihotjPeRUzWyCXDAsExqRBTVtvwhxDmbBpWQlFfi9FE+TuyVVrFwYRgd1n2mjpHpjx1NFsECpwNqfZfrKFrSvKZ6VVOfdeGcZkWHd/mLaewYzOGIr319eu5fUbqjJdjLQo8HtZUV7I44cagcwPzqdCKJBHfUcfkajygauSW2U9UxYIjEmxqawhmA0XLSnhoiVTn6efK9ZVBXl4fwOQ2em6qRIqyCPi7pA2GwPFYF1DxqTc2RStcDXJWVc1nI4jW7rjZiLkLior8nu5oHri1eCpYoHAmBSbaB8Ck3qx1Blej7BoDgyIx6YcX7qiDJ93dt6iLRDMgr1n2vjyo0dGbIht5q6atl58HmHhFBZwmemLbVJTFQrM2htnOsVSUc9WtxDYGEHatfUMcNv3dtLQ0U9FcX7CbQXnqv9+8iV2nGzlP27dPO0do3KNqnKgtoPq0sCEuf9N6qxcUIjf55kT4wMAKxYU4hFmZUVxTO6HzyymqvzTL/fR3DXAhuoQn3voEPUZ3gFqtjz0Yh3/96FDPHqwgb/+wW7CSeaqyWXRqPLJX+7jySNN3LRpSaaLM2/4vB6uv6iKa9ZWZLooKbG+KsTeT7+ejUvTs1F9IhYI0mjb3lp+80Idf/e68/nqey4lHI3yyV/uQ3XmXUSnm3to7xmc/MIZau8dnHJ5jzR08vGf7GXzslI+9cYNPHaokU9t25+SemercCTKJ36yl/957jQffuVqPv768zNdpHnl7lsu4Y5r12a6GCkTnKWtMWMsEKRJbVsvn/zlPrasKOPDr1zNigVFfOx15/PowQYefLF+2s8bjSr3/uElrv23J3jjV57iTEtPCks9TFX5xlPH2fIvj/DWrz7N7tOtST2uvXeQD31vF4V+H197zxb+4upVfORVq/nBc6f5rydemnG5atp6aejIrlZVfzjCHT/Yzc931/D3f7aOO69fn/HUEsZMxfzouJ1l0ajyiZ/sJRJVvvSOTUN9xX9x1Sp+vbeOT2/bx1VrFlBamNxOTDGNnX18/Md7eeroOV69rpLnT7fx9q89w//81RVDu0ulQnvvIP/w0708vL+Bq9dUcLihk7f819O85ZIl/MN166guSdwXG40qH/vRHs609PDD2142lDr371+/jrq2Xr748GEWlwZ4yyVT322puz/Mv/3uCN9++gRRhcUlAS5ZXsYly0u5ZHkpFy8pTWqf3WhUU5oQrXcgwoe+v4s/HGni02/aMGsLgIxJJQsEafCtp0/y9EvNfO6tF7NiwXCeEJ/Xw+fedjE3fuVP/OtvDvLFt29K+jmfPNLEx3+8h86+MP/nzRfx7iuWc6i+k/d+8zne+d/P8L0PXsEF1TPf3vDFs+189Ae7qGvr45NvuIAPXr2K7oEIX33iGF9/6gS/3VfPh1+5mlsuX8bCYP6IT753P3aUxw418tmbLuSyuBkPHo/whZs30dDRzz/89ETdys8AAB+9SURBVAUWBgNctSb5/tzHDzfyyV/so6atl3dfsZzzKovZc6aN3adb+c2LznaCVaEAt796Ne+4bBn5vpFbLKoqTx5p4suPHuVAXQev27CIm7cs5RVrKqY1y6S+vY8dJ1vYcbKFPxxp4nRLD19420becVn2JmczZiKSzn5bEbkOuBvwAt9Q1c+NOp8PfBfYAjQD71TVkxM959atW3Xnzp3pKXACp5t7+PULtTxyoIGFwXxec8FCXr1uIQtDY6cGhiNRnj/dxnu++RzXrK3g63++NWEXwed/e4ivPvES3//gFVw9wQBXNKq8WNPOz58/y3eeOcX5i4r5z1svHZouB3CssYt3f+NZ+gajfO+Dl09rgCkSVRo7+/jd/gb+9TcHWVDs5yvvupQtK8pGXHempYfP/fYQv3H3ci3O97GqoojzKosoK/Tz7adPcvOWpXzx5o0J693eO8g7vvYMx891sXVFOVevreAVayu4cHFJwhk2zV39fPaBA/xqTy1rFhbzubdePGa3rMbOPnadbOW+P51gx8lWFpcEuP3aNbx9yzLyvMITbgDYe6aNJaUFXL2mgt8dqKe1Z5CFwXzecukS3nBxNUX5PqJRJarO7yMcjdLSPUBTZz/nupyvDZ19vHC2jTMtzjqBQr+XLSvKeN+VK3nthtnPW2/MVIjILlXdmvBcugKBiHiBI8DrgLPADuBWVT0Qd81HgY2q+mERuQV4i6q+c6LnnW4gaOzo41RLD119YTr7w3T1henqH6R/MEpJYR5lhX4WFPkpL/aT7/Py+KFGtu2tZc+ZNsDZuq6po49ad9bPxUtKePX6hRT5vRyu7+RQfSfHmroYCEdZUOTnt397zbjZHvsGI1x/91P0D0Z4y6VLqC4pYElpAdWlAUoK8th+ooUnDzfx5JEmmrsHEIF3Xb6c//3GDQTyxm4ofrq5h3d941naega58/r15Ps8DEaUwUiUwUiUgUiUgfDwbTASpas/Qm1bLzVtvdS19zIYcV4Hrzy/kn9/5+YJNxB/8Ww7z59u5XhTF8fPdXO8qZuatl62rijj+395RcIyDv0dOvv4xlMn+MORJg65G4qUFuaxdUUZkajS2Rd2b4Oc6x5AVfnoq9bw0VevHvNJP56q8sdj5/jSI0fYfdp5068o9rP3bDtLSgu449o1vO3Spfh9HgbCUX5/qJGf7jrL44cbh5bzT6TI76UimM8FVSEuW1XOZSvL2FAdmhPz1s38kKlAcCXwz6r6Z+79uwBU9f/GXfOwe80zIuID6oFKnaBQ0w0EX3vyJT730KEpPWZDdYgbNy/mDRdXs6y8EFXlcEMnjx1s5PeHGnn+dCuqsCiUz7qqEOurgqxbFOTqtRUsStBiiLfrVCsf//EezrT2JnwjKivM45rzK3n1uoW8Ym0FCybZQLuuvZd3f+M5jjd1j3tNnlfI83rw+zwU5nmpLnUC0JIy5+vKBUW8fPWCafWh9w1GyPd5pjRI2tTZz5+OneOpo+d44WwbgTxn16tgwEcwkEdpYR7vvGwZ50+y6Xq8WDfQ3Y8dpb1nkL+65ryhADBeGZ493kxUFY8IXo/gEfB6PJQX5VFZHKAi6J836yDM3JWpQHAzcJ2q/qV7/73AFap6R9w1+9xrzrr3X3KvOTfquW4DbgNYvnz5llOnTk25PKebezjd0kNRvpdgwEdxfh5F+V7yfV7aewdp6R4YunX2DbJ1ZTlrFk48ANvWM4AqlE3w6XkysW6Z2rY+6tp7ae4aYOPSEjYuLZ3ygqT+cISa1t6hN/s8r2f4zd/rmbVdo4wx2WeiQJATH3NU9V7gXnBaBNN5juULClm+IPHuRZXB/Glt2jHVWT+JeD1CdUmBOxOnbNLrJ5Lv83JeCmcPGWPmh3R2cNYA8dMolrrHEl7jdg2V4AwaG2OMmSXpDAQ7gLUiskpE/MAtwLZR12wD3ud+fzPw+4nGB4wxxqReuqeP3gB8GWf66H2q+q8i8llgp6puE5EA8D3gEqAFuEVVj0/ynE3A1AcJHBXAuUmvyh1zqT5zqS5g9clmc6kukHx9Vqhqwkx2aQ0E2UZEdo43WJKL5lJ95lJdwOqTzeZSXSA19bFJ0MYYM89ZIDDGmHluvgWCezNdgBSbS/WZS3UBq082m0t1gRTUZ16NERhjjBlrvrUIjDHGjGKBwBhj5rl5EwhE5DoROSwix0TkzkyXZ6pE5D4RaXTzM8WOlYvIIyJy1P06sxwVs0RElonI4yJyQET2i8jfuMdztT4BEdkuInvd+nzGPb5KRJ5zX3M/chdW5gQR8YrIbhF5wL2fy3U5KSIvisgeEdnpHsvV11qpiPxURA6JyEERuTIVdZkXgcBNiX0PcD2wAbhVRDZktlRT9m3gulHH7gQeU9W1wGPu/VwQBj6uqhuAlwG3u3+PXK1PP3Ctqm4CNgPXicjLgM8D/66qa4BW4IMZLONU/Q1wMO5+LtcF4NWqujluvn2uvtbuBn6rquuBTTh/o5nXRVXn/A24Eng47v5dwF2ZLtc06rES2Bd3/zBQ7X5fDRzOdBmnWa9f4exbkfP1AQqB54ErcFZ7+tzjI16D2XzDyQv2GHAt8AAguVoXt7wngYpRx3LutYaTi+0E7iSfVNZlXrQIgCXAmbj7Z91juW6Rqta539cDObdNloisxEkx8hw5XB+3K2UP0Ag8ArwEtKlq2L0kl15zXwb+AYi69xeQu3UBUOB3IrLLTWkPuflaWwU0Ad9yu+2+ISJFpKAu8yUQzHnqfBzIqbnAIlIM/Az4W1XtiD+Xa/VR1Yiqbsb5NH05sD7DRZoWEXkj0KiquzJdlhS6WlUvxekavl1Erok/mUOvNR9wKfBVVb0E6GZUN9B06zJfAkEyKbFzUYOIVAO4XxszXJ6kiUgeThD4H1X9uXs4Z+sTo6ptwOM43Selbnp1yJ3X3FXAjSJyErgfp3vobnKzLgCoao37tRH4BU6gzsXX2lngrKo+597/KU5gmHFd5ksgSCYldi6KT+P9Ppy+9qwnzn6W3wQOquqX4k7lan0qRaTU/b4AZ7zjIE5AuNm9LCfqo6p3qepSVV2J83/ye1V9NzlYFwARKRKRYOx74PXAPnLwtaaq9cAZEVnnHnoNcIBU1CXTAyCzONByA3AEp+/2nzJdnmmU/4dAHTCI88nggzh9t48BR4FHgfJMlzPJulyN03x9Adjj3m7I4fpsBHa79dkHfMo9fh6wHTgG/ATIz3RZp1ivVwEP5HJd3HLvdW/7Y//7Ofxa2wzsdF9rv8TZ1nDGdbEUE8YYM8/Nl64hY4wx47BAYIwx85wFAmOMmed8k1+SXSoqKnTlypWZLoYxxuSUXbt2ndNx9ixOWyAQkfuA2OKUixKcF5z5yTcAPcD7VfX5yZ535cqV7Ny5M9XFNcaYOU1ETo13Lp1dQ99mbJK0eNcDa93bbcBX01gWY4wx40hbIFDVPwAtE1xyE/BddTyLs3KxOl3lOd3cw+/212PTZY0xZqRMDhYnnQhORG4TkZ0isrOpqWlaP+zBfXXc9r1d9AxEpvV4Y4yZq3Ji1pCq3quqW1V1a2VlwrGOSQUDznBIZ194kiuNMWZ+yWQgmNVEcMFAHgCdfYPp+hHGGJOTMhkItgF/Lo6XAe06nFM75YL5boug31oExhgTL53TR3+Ik7SqQkTOAp8G8gBU9WvAgzhTR4/hTB/9QLrKAtY1ZIwx40lbIFDVWyc5r8Dt6fr5o1nXkDHGJJYTg8WpUOy2CLqsRWCMMSPMm0BgXUPGGJPYvAkExX4fItY1ZIwxo82bQODxCMV+n80aMsaYUZIaLBaRrcArgMVAL852fI+oamsay5ZyxQGfdQ0ZY8woE7YIROQDIvI8cBdQABwGGnH2nH1URL4jIsvTX8zUCAZ81jVkjDGjTNYiKASuUtXeRCdFZDNO9tDTqS5YOgQDedYiMMaYUSZsEajqParaKyJXjT4nIlep6h5VfSx9xUut4nwfXTZGYIwxIyQ7WPyfSR7LakEbIzDGmDEm7BoSkSuBlwOVIvKxuFMhwJvOgqWD0zVkYwTGGBNvsjECP1DsXheMO94B3JyuQqWLtQiMMWasCQOBqj4JPCki31bVcfe7zBXBfB/94SgD4Sh+37xZQmGMMRNKNulcvojcC6yMf4yqXpuOQqXLcJqJQRYU52e4NMYYkx2SDQQ/Ab4GfAPI2b0eYxlIu/rDFgiMMcaVbCAIq+pX01qSWVBsieeMMWaMZDvKfy0iHxWRahEpj93SWrI0iHUNddjMIWOMGZJsi+B97te/jzumwHmpLU56hWJdQ9YiMLPsD0eaePRgA5+96aJMF8WYMZJqEajqqgS3nAoC4KwsBusaMrPvoX11fPeZU4Qj0UwXxZgxJltQdq2q/l5E3provKr+PD3FSo/4WUPGzKaGjn4A2nttxprJPpN1Db0S+D3wpgTnFMixQBDbt9haBGZ21bf3AdDaM2CBwGSdyRaUfdr9+oHZKU56+X0e8n0eSzxnZl1jpxMIWrqtNWqyT1JjBCJSIiJfEpGd7u3fRKQk3YVLh2DAR4e1CMwsGghHOdc1ADgtAmOyTbLTR+8DOoF3uLcO4FvpKlQ6WeI5M9uauvqHvm/ttkBgsk+y00dXq+rb4u5/RkT2pKNA6RYM2J4EZnY1dPQNfd/aYx9CTPZJtkXQKyJXx+64G9Uk3LUs2xXnWwZSM7sa2uMDgbUITPZJtkXwEeA77riAAC0MLzLLKcGAj3PnujNdDDOPxFoEBXleWqxryGShpAKBqu4BNolIyL3fkdZSpVEwkGcri82squ/oJ88rrKwoos1aBCYLJTtraIGI/AfwBPC4iNwtIgvSWrI0sa4hM9saO/pYGAxQXpRnLQKTlZIdI7gfaALehrMzWRPwo3QVKp1CAR9dA2GiUc10Ucw8Ud/RR1VJgLJCvw0Wm6yUbCCoVtV/UdUT7u3/AIvSWbB0CQbyUIWuAWsVmNnR0NHHolA+5UV+Gyw2WSnZQPA7EblFRDzu7R3Aw+ksWLrE8g3ZOIGZLQ0d/SwKBSgt9NPeO0jEWqMmyyQbCP4K+AEwAPTjdBV9SEQ6RSSnBo5tcxozm7r6w3T1h1kUClBe6LRG23ute8hkl2RnDQXTXZDZMpx4zv4ZTfrFpo5WhQKIOMdaugcoL/JnsFTGjJTsrCERkfeIyP927y8TkcvTW7T0GEpFbauLzSyIBYKFoXzKCp03fxsnMNkm2a6h/wKuBN7l3u8C7klLidIsaJvTmFkUCwSLQoGhVoDlGzLZJtmVxVeo6qUishtAVVtFJCfbttY1ZGZTbEOaRaHA0GIyaxGYbJNsi2BQRLw4m9EgIpXApHvuich1InJYRI6JyJ0Jzr9fRJpEZI97+8splX4abNaQmU317X0U5/sozvcNtwhsLYHJMsm2CP4D+AWwUET+FWdR2ScneoAbOO4BXgecBXaIyDZVPTDq0h+p6h1TK/b0Ffq9eMS6hszsaOx01hCAk2so3+exriGTdZKdNfQ/IrILeA1O0rk3q+rBSR52OXBMVY8DiMj9wE3A6EAwq0TETTNhn8pM+tW397EoFACc115Zod/STJisk+ysoXKgEfghznqCBhHJm+RhS4AzcffPusdGe5uIvCAiPxWRZeP8/Ntiu6M1NTUlU+QJOZvTWIvApF9DRz9VbiAAKCuyNBMm+yQ7RvA8Tn6hI8BR9/uTIvK8iGyZwc//NbBSVTcCjwDfSXSRqt6rqltVdWtlZeUMfpwjGPDZ9FGTdtGo0tjZx8L4QFCYZ4PFJuskGwgeAW5Q1QpVXQBcDzwAfBRnamkiNUD8J/yl7rEhqtqsqrF9/L4BzCSoJC0YsK4hk36tPQMMRpQqd4wA3BaBdQ2ZLJNsIHiZqg7lFlLV3wFXquqzQP44j9kBrBWRVe5U01uAbfEXiEh13N0bgcnGHVLCuobMbKiPW0MQU15oiedM9kl21lCdiPwjTo4hgHfijBN4GWcaqaqGReQOnOR0XuA+Vd0vIp8FdqrqNuB/iciNQBhn17P3T78qyQsGfLzUZIHApFdjbA1ByciuoTY38ZzXI5kqmjEjJBsI3gV8Gvile/9P7jEv8I7xHqSqDwIPjjr2qbjv7wLumkJ5U8I2pzGzIVGLoKzIjyp09A5SZvmGTJZIdvroOeCvxzl9LHXFmR1O19AgqoqIfSoz6TGUZygYN0bg5htq6RmwQGCyRlKBQETOBz4BrIx/jKpem55ipVcw4GMwovSHowTyvJkujpmjGjr6qCj2k+cdHoori883NPMJcMakRLJdQz8BvoYzsyeSvuLMjlDcngQWCEy6NHT0szAYGHGsvNDSTJjsk2wgCKvqV9Naklk0vDnNIJXB8SY9GTMz9e3OXsXxSguddZg2hdRkk2Snj/5aRD4qItUiUh67pbVkaRTMj2UgtQFjkz7xeYZihhPPWSAw2SPZFsH73K9/H3dMgfNSW5zZMZSB1FYXmzQZjEQ51zUwYsYQOEkP/T4PLRYITBZJdtbQqnQXZDbFdw0Zkw6NncP7EMRzEs/lWdeQySrJtggQkYuADcDQK1tVv5uOQqVbyN2cpsO6hkyaxO9VPFpZoSWeM9kl2emjnwZehRMIHsTJNfRHICcDgW1OY9KtoX14r+LRygot35DJLskOFt+MsxdBvap+ANgElKStVGlWZPsWmzSbqEVQXmT5hkx2STYQ9KpqFAiLSAhnb4KEewfkgjyvh4I8r40RmLSp7+gnzytDK4njlRXlWdeQySrJjhHsFJFS4OvALqALeCZtpZoFwYDPZg2ZtGns6GNhMIAnQWK5skI/bT0DRKOa8Lwxsy3ZWUMfdb/9moj8Fgip6gvpK1b6OXsSWCAw6VHfMXYNQUxZoZ+oQkffIKUJWgzGzLZku4YQkSUi8nJgOVAqItekr1jpVxzIo8O6hkyaNHSMXVUcE1tUZnsXm2yR7Kyhz+PsQXCA4VxDCvwhTeVKu5C1CEwaNXT084q1ibPKDaWZsAFjkyWSHSN4M7AublvJnBcM+Khzp/gZk0pd/WG6+sNjFpPFDKWZ6LYWqckOyXYNHQfy0lmQ2eZsTmP/iCb1hqaOlow/RgBYmgmTNSZsEYjIf+J0AfUAe0TkMWCoVaCq/yu9xUsf27fYpEssECwKJm4RxPYkaLNAYLLEZF1DO92vuxi18XyuCwZ89AxEbO9Yk3JDgWCcweIivxe/10OLdQ2ZLDFZIHgQqFTVA/EHReRCnEVlOSvo5hvq6gtTUjiner1MhjV0JE44FyMilFriOZNFJhsj+E+gIsHxcuDu1Bdn9gTdNBM2hdSkWn17H8X5Porzx/+cZWkmTDaZLBCsUdUxU0RV9SlgY3qKNDtsTwKTLrVtvQmTzcVzMpBaIDDZYbJAEJzgXE73p8S6hmzA2KRS32CEp19qZsvysgmvs3xDJptMFgiOicgNow+KyPU4U0pzlm1OY9Lh0YMNdPWHefMlSya8zlJRm2wy2WDx3wEPiMg7cGYOAWwFrgTemM6CpVswYKmoTer9ak8tC4P5vOy8BRNeF+sassRzJhtM2CJQ1SPAxcCTwEr39iSw0T2Xs4YCgY0RmBRp6xngicON3Lhp8aRTksuKnMRz9kHEZIPJFpSJm1biW5NcoykvWZoF82NjBNY1ZFLjwRfrGYwoN22euFsIoLzIef219AzY9GWTcZONETwuIn8tIsvjD4qIX0SuFZHvAO9LX/HSJ5DnwecR+0RmUuZXe2o4r7KIi5aEJr02ln7aZg6ZbDBZILgOJ9voD0WkTkQOiMgJ4ChwK/BlVf12msuYFiLibE5jgcCkQG1bL8+daOHNm5cgMnmff3ksENiAsckCE3YNqWof8F/Af4lIHs7isl5VbZuNwqWbk2/IuobMzG3bWwvATZsXJ3X9UOI5CwQmC0w2RhAAPgysAV4A7lPVOfMR2slAOmeqYzLol7tr2LyslBULipK6vswdI2iztQQmC0zWNfQdnOmiLwI3AP+W9hLNomDAZ7OGzIwdru/kUH0nb06yNQDOh5A8r1gqapMVJltHsEFVLwYQkW8C29NfpNkTDORR09ab6WKYHPerPTV4PcIbNiYfCJzEc7aozGSHyVoEQ+3WudQlFONsYG9NczN90ajyqz21XLWmgsrgxPmFRiu3fEMmS0zWItgkIh3u9wIUuPcFUFWdfJ5cFgsGfJZ0zszI86dbqWnr5eOvP3/Kjy0ryrPtKk1WmGzWkHe2CpIJQXcDe1VNasqfMfFUlR/vPEMgz8PrL6ya8uPLCv0ca+xKQ8mMmZpk9yyeFhG5TkQOi8gxEbkzwfl8EfmRe/45EVmZzvKMFgzkEYkqvYOR2fyxJse1dA9w7x9e4tp/e5If7zzLmzYunnDvgfGU2Z4EJktM/dWbJBHxAvcArwPOAjtEZNuo3c4+CLSq6hoRuQX4PPDOdJVptNg/b2dfmEJ/2n4VJseFI1Hq2vt4qamLX+yu4aEX6xmIRNm6ooy/vnYNb9hYPa3nLSt0UlFbi9RkWjrf/S4HjqnqcQARuR+4CYgPBDcB/+x+/1PgK7OZuyiWeO6+P56grMiPACIgCPZ/mTvGe7UoY0+oQlQh6j4oGlUiqoQjymA0ymBYCUej9A1GqG3r40xrDzWtvYSjzvXBgI93XbGcWy9fzrqqibbrmFxZoZ9IVPnyo0fx+zwjXnvxL7/Ya1FI/kWZzOt3Kv9l2fT/oOr8bWN/y9j3w+eH7ywozufWy5cneBYTL52BYAlwJu7+WeCK8a5R1bCItAMLgHPxF4nIbcBtAMuXp+6PurqymDyv8N9/yOmtFUwKeAR8Xg9+rwefV/B7PVSXFrBxaSlv3FjNsrJClpcXcsnyMgr8qRk6u6A6hNcj3P3Y0ZQ8nxlrfVXQAkEScqI/RFXvBe4F2Lp1a8paCxctKWH/Z64jHI26PweUkZ8oJixXqgpipiTRh9PxulYSHY2liPaI4BHnsZOljU6Hq9ZUcPhfriOiw59oY592Y4aOT+F5p9KgTqZLKhuTC0vsbxfXevfE1SWbWjC5IJ2BoAZYFnd/qXss0TVnRcQHlADNaSzTGH6fB396x8yNGZfP68mNT2NmTkvnO+AOYK2IrBIRP3ALsG3UNdsYTmN9M/D7XNzbwBhjclnaPoy4ff53AA8DXpyEdftF5LPATlXdBnwT+J6IHANacIKFMcaYWSS59gFcRJqAU9N8eAWjBqJz3Fyqz1yqC1h9stlcqgskX58VqlqZ6ETOBYKZEJGdqro10+VIlblUn7lUF7D6ZLO5VBdITX1slNQYY+Y5CwTGGDPPzbdAcG+mC5Bic6k+c6kuYPXJZnOpLpCC+syrMQJjjDFjzbcWgTHGmFEsEBhjzDw3bwLBZHsjZDsRuU9EGkVkX9yxchF5RESOul/LMlnGZInIMhF5XEQOiMh+Efkb93iu1icgIttFZK9bn8+4x1e5+2wcc/fd8Ge6rMkSEa+I7BaRB9z7uVyXkyLyoojsEZGd7rFcfa2VishPReSQiBwUkStTUZd5EQji9ka4HtgA3CoiGzJbqin7NnDdqGN3Ao+p6lrgMfd+LggDH1fVDcDLgNvdv0eu1qcfuFZVNwGbgetE5GU4+2v8u6quAVpx9t/IFX8DHIy7n8t1AXi1qm6Om2+fq6+1u4Hfqup6YBPO32jmdVHVOX8DrgQejrt/F3BXpss1jXqsBPbF3T8MVLvfVwOHM13GadbrVzgbGOV8fYBC4HmclOvnAJ97fMRrMJtvOAkiHwOuBR7ASeKak3Vxy3sSqBh1LOdeazhJOU/gTvJJZV3mRYuAxHsjLMlQWVJpkarWud/XA4syWZjpcLcnvQR4jhyuj9uVsgdoBB4BXgLaVDXsXpJLr7kvA/8ARN37C8jduoCTxft3IrLL3dsEcvO1tgpoAr7ldtt9Q0SKSEFd5ksgmPPU+TiQU3OBRaQY+Bnwt6raEX8u1+qjqhFV3YzzafpyYH2GizQtIvJGoFFVd2W6LCl0tapeitM1fLuIXBN/Modeaz7gUuCrqnoJ0M2obqDp1mW+BIJk9kbIRQ0iUg3gfm3McHmSJiJ5OEHgf1T15+7hnK1PjKq2AY/jdJ+UuvtsQO685q4CbhSRk8D9ON1Dd5ObdQFAVWvcr43AL3ACdS6+1s4CZ1X1Off+T3ECw4zrMl8CQTJ7I+Si+P0c3ofT1571xNkW65vAQVX9UtypXK1PpYiUut8X4Ix3HMQJCDe7l+VEfVT1LlVdqqorcf5Pfq+q7yYH6wIgIkUiEox9D7we2EcOvtZUtR44IyLr3EOvwdkDfuZ1yfQAyCwOtNwAHMHpu/2nTJdnGuX/IVAHDOJ8MvggTt/tY8BR4FGgPNPlTLIuV+M0X18A9ri3G3K4PhuB3W599gGfco+fB2wHjgE/AfIzXdYp1utVwAO5XBe33Hvd2/7Y/34Ov9Y2Azvd19ovgbJU1MVSTBhjzDw3X7qGjDHGjMMCgTHGzHMWCIwxZp6zQGCMMfOcBQJjjJnnfJNfYowBEJEI8CLO/80J4L3qLCAzJqdZi8CY5PWqk8HyIqAFuD3TBTImFSwQGDM9z+AmXhORJ0Rkq/t9hZueARF5v4j8XER+6+aK/0LmimvM+CwQGDNF7v4WryG5NCWbgXcCFwPvFJFlk1xvzKyzQGBM8grcVNOxVL+PJPGYx1S1XVX7cPLCrEhnAY2ZDgsExiSvV51U0ytwNmuJjRGEGf5fCox6TH/c9xFsgobJQhYIjJkiVe0B/hfwcTc180lgi3v65vEeZ0y2skBgzDSoaizb6K3A/wM+IiK7gYqMFsyYabDso8YYM89Zi8AYY+Y5CwTGGDPPWSAwxph5zgKBMcbMcxYIjDFmnrNAYIwx85wFAmOMmef+fxScOHGQBa5YAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "random.seed(1000)\n", + "np.random.seed(1000)\n", + "\n", + "klds = get_klds(60, 30, 0.15)\n", + "R = run_changepoint(klds)\n", + "\n", + "runs = np.arange(len(klds))\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot(2, 1, 1)\n", + "ax.plot(runs, klds)\n", + "ax.set_ylabel('KLD')\n", + "\n", + "ax = fig.add_subplot(2, 1, 2, sharex=ax)\n", + "window_size = 10 # How many time points to get before evaluating\n", + " # probability\n", + "ax.plot(runs[:-window_size], R[window_size, window_size:-1])\n", + "\n", + "ax.set_xlabel('Run')\n", + "ax.set_ylabel('P(Changepoint)')\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}