-
-[](#contributors-)
-
-
-# CosmoStat Tutorials
+## Importance sampling
---
-CosmoStat is committed to the philosophy of reproducible research, endeavouring
-to provide source code and data for all publications. In this spirit, we have
-additionally put significant effort into providing useful educational
-materials. The aim being to provide other researchers with an in-depth
-understanding of the various tools we use in our work.
-
-We always welcome new tutorial requests, just click [here](https://github.com/CosmoStat/Tutorials/issues/new?assignees=&labels=tutorial+request&template=tutorial-request.md&title=%5BTutorial%5D+Your+idea+for+a+tutorial).
-
-## Contents
-
-1. [Career Development](#Career-Development)
-1. [Cosmology](#Cosmology)
-1. [Data Science](#Data-Science)
-1. [Software Carpentries](#Software-Carpentries)
-
-### Career Development
-
-1. **[Presentation Tips](https://cosmostat.github.io/Tutorials/presentation_tips/presentation_tips/index.html)** | [](https://github.com/CosmoStat/Tutorials/tree/presentations) [](https://www.youtube.com/playlist?list=PLquvp9RIoLGCmfAZTOKJz7qDO0NqCnMAV)
- *Authors:* [@pettorin](https://github.com/pettorin)
- This tutorial provides tips on how to adapt presentations for different goals. The tutorial is mainly meant for scientists, but several tips can be useful for other types of talks.
-
-### Cosmology
-
-1. **[Cosmology with Python wrappers for Einstein-Boltzmann Codes](https://github.com/CosmoStat/Tutorials/tree/boltzmann)** | [](https://github.com/CosmoStat/Tutorials/tree/boltzmann) [](https://www.youtube.com/watch?v=EMpCxUvF3lc)
- *Authors:* [@santiagocasas](https://github.com/santiagocasas)
- This tutorial is comprised of a series of Jupyter notebooks with simple demonstrations and exercises on how to use `CAMB` and `CLASS` using python wrappers. The code is designed for non-experts in the field, therefore it includes relatively simple explanations of cosmological concepts. It intends to show a general overview of the things that are possible with Einstein-Boltzmann codes and python.
-
-2. **[Power Spectrum](http://github.com/CosmoStat/Tutorials/tree/power-spectrum)** | [](http://github.com/CosmoStat/Tutorials/tree/power-spectrum) [](https://colab.research.google.com/github/CosmoStat/Tutorials/blob/power-spectrum/WL_power_spectrum.ipynb)
- *Authors:* [@b-remy](https://github.com/b-remy), [@dlanzieri](https://github.com/dlanzieri), [@EiffL](https://github.com/EiffL)
- This tutorial presents how to compute a power spectrum from a lensing map (in the flat sky approximation) and how to compute the corresponding theory spectrum using a cosmology code.
-
-### Data Science
-
-1. **[Introduction to Python](https://github.com/CosmoStat/Tutorials/tree/python)** | [](https://github.com/CosmoStat/Tutorials/tree/python) [](https://www.youtube.com/watch?v=cbekcGxm70Q&list=PLquvp9RIoLGBFfsz8dqzPtLeaEngHYybM) [](https://mybinder.org/v2/gh/CosmoStat/Tutorials/python)
- *Authors:* [@sfarrens](https://github.com/sfarrens), [@santiagocasas](https://github.com/santiagocasas)
- * [Tutorial 1: Beginner Topics](https://github.com/CosmoStat/Tutorials/tree/python#tutorial-1-beginner-topics)
- The objective of this tutorial is to provide a first look at Python for beginners. The level is aimed at individuals with little or no experience whatsoever with Python. Experienced users are unlikely to benefit from this tutorial.
- * [Tutorial 2: Intermediate and Advanced Topics](https://github.com/CosmoStat/Tutorials/tree/python#tutorial-2-intermediate-and-advanced-topics)
- The objective of this tutorial is to provide a more in-depth look at object-oriented and pythonic coding. The level is aimed at individuals with some experience with Python and good knowledge of basic object types. This tutorial will likely benefit all except the most advanced users.
-
-2. **[Python Optimisation and Memory Management](https://github.com/CosmoStat/Tutorials/tree/profiling)** | [](https://github.com/CosmoStat/Tutorials/tree/profiling)
- *Authors:* [@sfarrens](https://github.com/sfarrens)
- This tutorial introduces some techniques for determinisitic and memory profiling of Python scripts, followed by some tips on how to implement some basic optimisation.
-
-3. **[C++](https://github.com/CosmoStat/Tutorials/tree/CPlusPlus)** | [](https://github.com/CosmoStat/Tutorials/tree/CPlusPlus) [](https://www.youtube.com/watch?v=47Ldg6i2B8Y&list=PLquvp9RIoLGD5yDykupoCueNp2EHEjuZB)
- *Authors:* [@kansal9](https://github.com/kansal9)
- This tutorial aims to help newcomers learn C++ and solve their programming problems. It is assumed that readers are already familiar with C, or at least that they do not have any difficulty reading C code. In other words, those who have experience in C and peo·ple who desire to quickly understand the features of modern C++ in a short period of time are well suited to follow this tutorial.
-
-4. **[Sparsity](https://github.com/CosmoStat/Tutorials/tree/ada)** | [](https://github.com/CosmoStat/Tutorials/tree/ada) [](https://mybinder.org/v2/gh/CosmoStat/Tutorials/ada)
- *Authors:* [@sfarrens](https://github.com/sfarrens)
- This tutorial is comprised of a series of Jupyter notebooks that demonstrate how the tools implemented in sparsity work as well as showing the applicability of these tools to various simple problems.
-
-5. **[Low-Rank](https://github.com/CosmoStat/Tutorials/tree/low-rank)** | [](https://github.com/CosmoStat/Tutorials/tree/low-rank) [](https://mybinder.org/v2/gh/CosmoStat/Tutorials/low-rank)
- *Authors:* [@sfarrens](https://github.com/sfarrens)
- The objective is to provide a beginner level introduction to the concept of low-rank approximation, in particular as a regularisation method for solving linear inverse problems.
-
-6. **[TensorFlow](https://github.com/CosmoStat/Tutorials/tree/tensorflow-tutorial)** | [](https://github.com/CosmoStat/Tutorials/tree/tensorflow-tutorial)
- *Authors:* [@zaccharieramzi](https://github.com/zaccharieramzi), [@EiffL](https://github.com/EiffL)
- * *First Steps with TensorFlow*: [](https://www.youtube.com/watch?v=kawHQpxytLo) [](https://colab.research.google.com/github/CosmoStat/Tutorials/blob/tensorflow-tutorial/TensorFlowFirstSteps.ipynb)
- A short introduction to the basic concepts underpinning TensorFlow, in particular automatic differentation.
- * *MRI reconstruction with TensorFlow*: [](https://www.youtube.com/watch?v=HqNm5fowdTU) [](https://colab.research.google.com/github/CosmoStat/Tutorials/blob/tensorflow-tutorial/MRIReconstructionWithTensorflow.ipynb)
- An illustration of how to use TensorFlow to solve a toy MRI reconstruction problem using classical iterative reconstructions, and a deep learning approach.
-
-7. **[Blind Source Separation](https://github.com/RCarloniGertosio/bss_tutorial)** | [](https://github.com/RCarloniGertosio/bss_tutorial) [](https://mybinder.org/v2/gh/RCarloniGertosio/bss_tutorial/master)
- *Authors:* [@RCarloniGertosio](https://github.com/RCarloniGertosio)
- The goal of this tutorial is to present Blind Source Separation (BSS) problems and the main methods to solve them. This tutorial does not provide in-depth mathematical explanations for every methods; the emphasis is rather on illustrations and intuition.
-
-7. **[Introduction to MCMC and Bayesian inference](https://colab.research.google.com/drive/1EUF7-CwamhAT6wedntStdb4x3dDbihoo?usp=sharing)** | [](https://www.youtube.com/watch?v=EfWEGBHeA3k) [](https://colab.research.google.com/drive/1EUF7-CwamhAT6wedntStdb4x3dDbihoo?usp=sharing)
- *Authors:* [@EiffL](https://github.com/EiffL)
- This tutorial is a practical introduction to Bayesian inference using emcee and
- touching on questions related to measurement errors and covariance.
-
-### Software Carpentries
-
-1. **[Git Tutorial](https://github.com/zaccharieramzi/git-tuto)** | [](https://github.com/zaccharieramzi/git-tuto) [](https://docs.google.com/presentation/d/1vfsG__2-T7xJYGKFs9HfPKmaoMN1Je0V0h7gLyiY1AU/edit?usp=sharing) [](https://www.youtube.com/watch?v=S1A2qSA0TWo)
- *Authors:* [@zaccharieramzi](https://github.com/zaccharieramzi)
- This tutorial will help you practice the basics of the GitHub flow and how to work on open source projects.
-
-2. **[Jekyll Tutorial](https://github.com/sfarrens/jekyll_tutorial)** | [](https://github.com/sfarrens/jekyll_tutorial) [](https://www.youtube.com/watch?v=qRxbbSaVW7M)
- *Authors:* [@sfarrens](https://github.com/sfarrens)
- The objective of this tutorial is to introduce Jekyll and show you how to build a website that you can host on GitHub for free.
-
-3. **[Make and CMake Tutorial](https://github.com/sfarrens/make-tutorial)** | [](https://github.com/sfarrens/make-tutorial) [](https://www.youtube.com/watch?v=K27-uncFZgM)
- *Authors:* [@sfarrens](https://github.com/sfarrens)
- This tutorial is designed to provide a first look at using Make and CMake to build C/C++ projects.
-
-4. **[Introduction to Docker for Data Scientists](https://cosmostat.github.io/Tutorials/docker/docker-introduction)** | [](https://github.com/CosmoStat/Tutorials/tree/docker) [](https://www.youtube.com/watch?v=N7cYFfIWjdw)
- *Authors:* [@EiffL](https://github.com/EiffL)
- This tutorial demonstrates how to create a Docker container to distribute a complete Jupyter notebook environment.
-
-5. **[Bash and Scripting Tutorial for Researchers](https://github.com/CosmoStat/Tutorials/tree/bash_tutorial/README.md)** | [](https://github.com/CosmoStat/Tutorials/tree/bash_tutorial) [](https://www.youtube.com/watch?v=NU8dYd3RM-E&t=37s) [](https://mybinder.org/v2/gh/CosmoStat/Tutorials/bash_tutorial?filepath=intro.ipynb)
+> Author: [Martin Kilbinger](http://www.cosmostat.org/people/kilbinger)
+> Year: 2022
+> Email: [martin.killbinger@cea.fr](mailto:martin.killbinger@cea.fr)
- *Authors:* [@fadinammour](https://github.com/fadinammour), [@JulienNGirard](https://github.com/JulienNGirard)(screen segment)
- This tutorial gives the key concepts of bash and its scripts to help researchers manage and process their data with more ease.
+This is a very brief tutorial about how importance sampling for Bayesian parameter inference.
+Specifically, it describes the case where an existing sample (e.g. a Monte-Carlo Markov chain) can be importance-sampled
+to obtain samples under a joint posterior distribution.
-## Contributors ✨
+For a more general tutorial of Bayesian paremeter inference see the Introduction to MCMC and Bayesian inference linked from CosmoStat Tutorial site.
-Thanks goes to these wonderful people ([emoji key](https://allcontributors.org/docs/en/emoji-key)):
+This tutorial contains a single Jupyter notebook with description and executable code.
-
-
-
-Python (require >=3.5)
+* NumPy (recommend >=1.16.2)
+* Emcee (recommend >=3.1.2)
+* Matplotlib (recommend >=3.0.3)
+* Jupyter (recommend >=1.0.0)
-This project follows the [all-contributors](https://github.com/all-contributors/all-contributors) specification. Contributions of any kind welcome!
+All of the packages listed above can easily be installed using `pip`.
diff --git a/docs/docker b/docs/docker
deleted file mode 160000
index d8f223d..0000000
--- a/docs/docker
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit d8f223d4381c2603b39681bdcd044050bda63cb5
diff --git a/docs/presentation_tips b/docs/presentation_tips
deleted file mode 160000
index 2700adc..0000000
--- a/docs/presentation_tips
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit 2700adc9882714e371ede3a22ba6d6033996482c
diff --git a/docs/tensorflow-introduction b/docs/tensorflow-introduction
deleted file mode 160000
index 1592869..0000000
--- a/docs/tensorflow-introduction
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit 15928695d2b6a1bbc4cab31529bd9399ae4f949d
diff --git a/example_MC.jpeg b/example_MC.jpeg
new file mode 100644
index 0000000..5bade99
Binary files /dev/null and b/example_MC.jpeg differ
diff --git a/importance_sampling.ipynb b/importance_sampling.ipynb
new file mode 100644
index 0000000..5f9ae56
--- /dev/null
+++ b/importance_sampling.ipynb
@@ -0,0 +1,591 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "3814beb8",
+ "metadata": {},
+ "source": [
+ "# Importance sampling\n",
+ "\n",
+ "Illustrated with a specific example of importance-sampling an existing MCM chain."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d87e396b-8f4b-48ba-a88c-a71f0f338d51",
+ "metadata": {},
+ "source": [
+ "## Monte-Carlo integration"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c6c62e88-a680-45bc-b319-e0625c5f2a24",
+ "metadata": {},
+ "source": [
+ "An integral\n",
+ "\n",
+ "$$F = \\int_a^b {\\rm d} x \\, f(x) $$\n",
+ "\n",
+ "can be approximated by evaluating $f$ at random points, drawn from a uniform distribution\n",
+ "between $a$ and $b$:\n",
+ "\n",
+ "$$\\int_a^b {\\rm d} x \\, f(x) \\approx \\frac{b-a} N \\sum_{i=1}^N f(x_i),$$\n",
+ "\n",
+ "with sample points written as\n",
+ "\n",
+ "$$x_i ~ \\sim U_{[a; b]}.$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e9ac7929-7982-4437-be1d-0a7912a69795",
+ "metadata": {},
+ "source": [
+ "Proof:\n",
+ " \n",
+ "$$E \\left[ \\frac{b-a} N \\sum_{i=1}^N f(x_i) \\right]\n",
+ "= \\frac{b-a} N \\sum_{i=1}^N E \\left[ f(x_i) \\right]\n",
+ "= \\frac{b-a} N \\left( N \\frac 1 {b-a} \\int_a^b {\\rm d} x \\, f(x) \\right)\n",
+ "= F\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "578fec40-49d4-419f-bef9-2b286af2c591",
+ "metadata": {},
+ "source": [
+ "This is illustrated in the following example.\n",
+ "\n",
+ "The goal is to estimate the grey area under the curve.\n",
+ "\n",
+ "We pick uniformly drawn random points $x_i$\n",
+ "and compute the area of the rectangle between $a = -2$ and $b = 5$, which is $(b - a) f(x_i)$.\n",
+ "\n",
+ "An (unbiased) estimate of the grey area is the average of all rectangles, which is the above equation.\n",
+ "\n",
+ "![](\"example_MC.jpeg\")
\n",
+ "(From https://towardsdatascience.com/the-basics-of-monte-carlo-integration-5fe16b40482d)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "91158d0f-0890-4208-8011-93cb6902d192",
+ "metadata": {},
+ "source": [
+ "## Drawing from arbitrary distributions"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "cce38eff-88a6-4ce1-9a83-0b6adf0466bf",
+ "metadata": {},
+ "source": [
+ "$$F = \\int_a^b {\\rm d} x f(x) = \\int {\\rm d} y \\, p(y) f(y)$$\n",
+ "\n",
+ "with variable change\n",
+ "\n",
+ "$${\\rm d} y \\, p(y) = \\frac{{\\rm d} x}{b-a}$$,\n",
+ "\n",
+ "where $p$ is a normalised probability distribution $\\int {\\rm d} y \\, p(y) = 1$.\n",
+ "\n",
+ "This is approximated by\n",
+ "\n",
+ "$$F \\approx \\frac 1 N \\sum_{i=1}^N f(y_i)$$\n",
+ "\n",
+ "with random points $y_i$ drawn from the distribution $p$, $y \\sim p$.\n",
+ "\n",
+ "Note that the distribution $p$ is not evaluated in the sum. It only serves as a kind of\n",
+ "weight. The density of sample points $y_i$ is proportional to this weight. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fc4c50ff-8b7b-4785-b885-c4969d258f27",
+ "metadata": {},
+ "source": [
+ "Proof:\n",
+ "\n",
+ "$$E \\left[ \\frac 1 N \\sum_{i=1}^N f(y_i) \\right] = \\frac 1 N \\sum_{i=1}^N E \\left[ f(y_i) \\right]\n",
+ "= \\frac 1 N \\sum_{i=1}^N \\int {\\rm d} y \\, p(y) f(y)\n",
+ "= \\frac 1 N N F = F\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d6631af1-dbaa-4cf3-aada-0aeb0110e703",
+ "metadata": {},
+ "source": [
+ "## Monte-Carlo Markov chain"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b74e2ebc",
+ "metadata": {},
+ "source": [
+ "A Monte-Carlo Markov chain of points $\\{x_i \\}_{i=1}^N$ represents, when converged, a sample under a posterior distribution $p(x) \\propto \\pi(x) L(x)$, with $\\pi$ being the prior and $L$ the likelihood function. We write as above\n",
+ "\n",
+ "$$x \\sim p.$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "75142390",
+ "metadata": {},
+ "source": [
+ "This allows us to estimate (high-dimensional) integrals of the form\n",
+ "\n",
+ "$$\\int {\\rm d} x \\, p(x) f(x) = \\int {\\rm d} x \\, \\pi(x) L(x) f(x).$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "84ec49e5",
+ "metadata": {},
+ "source": [
+ "The (unbiased) estimator is\n",
+ "\n",
+ "$$ \\int {\\rm d} x \\, p(x) f(x) \\approx \\frac 1 N \\sum_{i=1}^N f(x_i),$$\n",
+ "\n",
+ "which is just the normalised sum over the function $f$ evaluated at all chain points."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f4e53939-134f-483b-b7c9-47d5ff468468",
+ "metadata": {},
+ "source": [
+ "For example, for the mean of a parameter we set $f(x) = x$. The mean estimated by MCMC is simply\n",
+ "\n",
+ "$$\\bar x \\approx \\frac 1 N \\sum_{i=1}^N x_i.$$\n",
+ "\n",
+ "Another example: To plot the 2D $68\\%$ credible region (contour), pixelise yout prior range in 2D. (Ignoring further dimensions is equivalent of marginalising over those dimensions.) Then sort the pixels by density (= number of sample points in each pixel), and set $f = 1$ in the first $68\\%$ of pixels, and $0$ for the remaining pixels, and plot $f$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "06c67f01-269f-4ffb-ac4f-79195bef87fe",
+ "metadata": {},
+ "source": [
+ "## Importance sampling"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5ab283bf",
+ "metadata": {},
+ "source": [
+ "Importance sampling is based on the following identity\n",
+ "\n",
+ "$$\\int {\\rm d} x \\, p(x) f(x) = \\int {\\rm d} x \\, q(x)\\frac{p(x)}{q( x)} f(x),$$\n",
+ "\n",
+ "which holds if the support of the proposal distribution $q$ is not smaller than the support of $p$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "42545e20-0175-4c5d-962c-79ae7a3af804",
+ "metadata": {},
+ "source": [
+ "This integral can be estimated with samples under $q$ instead of $p$,\n",
+ "\n",
+ "$$x \\sim q,$$\n",
+ "\n",
+ "if the samples are multiplied by the importance weights\n",
+ "\n",
+ "$$w_i = \\frac{p(x_i)}{q(x_i)},$$\n",
+ "\n",
+ "since\n",
+ "\n",
+ "$$ \\int {\\rm d} x \\, p(x) f(x) = \\int {\\rm d} x \\, q(x)\\frac{p(x)}{q( x)} f(x) \\approx \\frac 1 N \\sum_{i=1}^N w_i f(x_i).$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1c4385ad-ebf0-43bf-a937-09a7f149b7d9",
+ "metadata": {},
+ "source": [
+ "Importance sampling can be useful if sampling under $p$ is very time-consuming of difficult, but evaluating $p$ at points $x_i$ is not. One then chooses a simple proposal function $q$, e.g. a (multi-variate) normal.\n",
+ "\n",
+ "Importance sampling permits to separate sampling (fast; under $q$) and evaluating the original posterior (potentiall slow, but can now be done in parallel using the already sampled points $x_i$)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ded5de23-66f3-48c5-831b-6e6f1b425cd3",
+ "metadata": {},
+ "source": [
+ "### Importance-sampling an existing MCM chain"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "01e3f214",
+ "metadata": {},
+ "source": [
+ "A recurring case is the following. We have already run an MCMC and obtained samples $x \\sim p$ under the posterior $p(x) = \\pi(x) L(x)$.\n",
+ "\n",
+ "Next, we would like to get samples under a new, combined, posterior $p_{\\rm comb}(x) = \\pi(x) L(x) L_2(x)$. For example, $L$ is the CMB likelihood, $L_2$ is the one for weak lensing. We want to get the combined weak-lensing + CMB constraints.\n",
+ "\n",
+ "Now, instead of re-running a new chain sampling from $p_{\\rm comb}$, we can see that the integral\n",
+ "\n",
+ "$$\\int {\\rm d} x p_{\\rm comb} f(x) = \\int {\\rm d} x \\pi(x) L(x) L_2(x) f(x)$$\n",
+ "\n",
+ "has the same form to the modified one above, if we identify\n",
+ "\n",
+ "$$q(x) = \\pi(x) L(x) \\quad \\mbox{and} \\quad L_2(x) = \\frac{p(x)}{q(x)}.$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "678f6d9b",
+ "metadata": {},
+ "source": [
+ "So we can use the existing chain points $x_i$ sampled under $p=q$ and multiply each data point $x_i$ with the \"importance weight\" $L_2(x_i)$, and can thus obtain estimators of integrals with the combined posterior $p_{\\rm comb}$!"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "825765c7",
+ "metadata": {},
+ "source": [
+ "Let's test this."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "d23712d6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%matplotlib inline\n",
+ "\n",
+ "import matplotlib.pylab as plt\n",
+ "import numpy as np\n",
+ "import emcee"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "02225b65-2c7d-4fc3-a752-ff8159be3520",
+ "metadata": {},
+ "source": [
+ "First, we sample a 1-D distribution with emcee to see whether this simple case works."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "26860f9e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Define simple log of a 1D normal distribution with mean μ and width σ\n",
+ "def log_prob(x, μ, σ):\n",
+ " diff = x - μ\n",
+ " return -0.5 * (diff / σ)**2"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "4e7a619b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Input values\n",
+ "μ = 0\n",
+ "σ = 3"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "45f3f7d9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Set x array for plotting\n",
+ "n = 5\n",
+ "x = np.linspace(-σ*n, σ*n, 1000)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "5da74549",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Set up MCMC\n",
+ "ndim = 1\n",
+ "nwalkers = 5\n",
+ "\n",
+ "# Random starting points\n",
+ "p0 = np.random.rand(nwalkers, ndim)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "4dcffc21-a30c-48d5-b014-db396b6a900f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Initialise emcee samples\n",
+ "sampler = emcee.EnsembleSampler(nwalkers, ndim, log_prob, args=[μ, σ])\n",
+ "n_burn = 100\n",
+ "n_chain = 5000"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "37e1b0a8-ae4d-4e0b-a768-e896f7f27c64",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Run burn-in phase and remove points from chain\n",
+ "state = sampler.run_mcmc(p0, n_burn)\n",
+ "sampler.reset()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "6803b569-df99-42da-a7ec-b554ecf75453",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Run rest of the chain\n",
+ "sampler.run_mcmc(state, n_chain);"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "15bbb805",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0eUlEQVR4nO3deXxU5fX48c/JzpawJKxhU8IeQQiLsogLihtYhQpqq62tv7ZqW21t6bfWUmu/31qrtFa0tdZKVVy+fqulihsiCgo1AZQdEjBCWEISMIQl+/n9cScwhEkyCXNzZ5Lzfr3yysxdMucmmZw893me84iqYowxxtQW5XUAxhhjwpMlCGOMMQFZgjDGGBOQJQhjjDEBWYIwxhgTUIzXAYRKcnKy9uvXz+swjDEmoqxZs6ZQVVMC7WsxCaJfv35kZWV5HYYxxkQUEfmirn12i8kYY0xAliCMMcYEZAnCGGNMQC2mD8IY07JVVFSQl5dHaWmp16FEpISEBFJTU4mNjQ36HEsQxpiIkJeXR4cOHejXrx8i4nU4EUVVKSoqIi8vj/79+wd9nt1iMsZEhNLSUrp06WLJoQlEhC5dujS69WUJwhgTMSw5NF1TvneWIIwxxgRkCcIYY8LAa6+9xubNm088v++++1i6dGmjv05ubi7Dhw8PSUzWSW1MKM1Ph+JdkNQH7trgdTQmgrz22mtcddVVDB06FID777/f44isBWFMaBXvgnnFzmfT4lxzzTWMHj2aYcOG8eSTTwLQvn17fv7znzNixAjGjx9Pfn4+AP/+978ZN24c5557Lpdccgn5+flUV1eTlpZGQUEBANXV1QwYMIAPPviAxYsXc8899zBy5Eh27NjBLbfcwiuvvAJAZmYm559/PiNGjGDs2LGUlJSQm5vLpEmTGDVqFKNGjeLjjz8O+fVaC8IYt1mrIuR+9e9NbN57OKRfc2jPRH559bB6j3n66afp3Lkzx48fZ8yYMVx33XUcPXqU8ePH85vf/Iaf/OQn/PWvf+Xee+9l4sSJrF69GhHhqaee4ne/+x0PP/wwN910E88//zw//OEPWbp0KSNGjOCCCy5g+vTpXHXVVcycOfOU1ywvL+f666/npZdeYsyYMRw+fJg2bdrQtWtX3n33XRISEsjOzmbOnDkhr0fnagtCRKaJyDYRyRGRuQH2TxaRtSJSKSIzA+xPFJE8EXnMzTiNcZW1KlqMRx999ERLYffu3WRnZxMXF8dVV10FwOjRo8nNzQWceRuXXXYZ6enpPPTQQ2zatAmAb37zm/zjH/8AnITzjW98o97X3LZtGz169GDMmDEAJCYmEhMTQ0VFBd/+9rdJT09n1qxZp/RfhIprLQgRiQYWAFOBPCBTRBarqv9V7AJuAX5cx5f5NfChWzEaYyJTQ//pu2H58uUsXbqUVatW0bZtW6ZMmUJpaSmxsbEnhpBGR0dTWVkJwJ133sndd9/N9OnTWb58OfPmzQOgd+/edOvWjWXLlvHJJ5/w/PPPNyme+fPn061bNz777DOqq6tJSEgIyXX6c7MFMRbIUdWdqloOvAjM8D9AVXNVdT1QXftkERkNdAPecTFGY4wJSnFxMZ06daJt27Zs3bqV1atXN3h8r169AFi4cOEp+771rW9x0003MWvWLKKjowHo0KEDJSUlp32dQYMGsW/fPjIzMwEoKSmhsrKS4uJievToQVRUFM8++yxVVVWhuMxTuJkgegG7/Z7n+bY1SESigIepu2VRc9xtIpIlIlk1nT7GGOOGadOmUVlZyZAhQ5g7dy7jx4+v9/h58+Yxa9YsRo8eTXJy8in7pk+fzpEjR065vTR79mweeughzj33XHbs2HFie1xcHC+99BJ33nknI0aMYOrUqZSWlvK9732PhQsXMmLECLZu3Uq7du1Ce8Hg1Ohw4wOYCTzl9/xrwGN1HPsMMNPv+R3AT3yPb6nrPP+P0aNHqzGe+2XiqZ/r2mYabfPmzV6HEDKZmZk6ceLEZn/dQN9DIEvr+Lvq5iimPUBvv+epvm3BOA+YJCLfA9oDcSJyRFVP6+g2xphI8tvf/pYnnniiyX0PzcnNBJEJpIlIf5zEMBu4IZgTVfXGmscicguQYcnBGNMSzJ07l7lzI+PPmWt9EKpaiXOr6G1gC/Cyqm4SkftFZDqAiIwRkTxgFvAXEdnkVjzGGGMax9WJcqq6BFhSa9t9fo8zcW491fc1nsHpozDGGNOMbCa1MW5I6gPzkk4+NiYCWYIwxg1WUsO0AJYgjDkTVmfJOzXf+1AJo5/hlClT+P3vf09GRoancViCMOZM1NRZqrmd5HO4tILnV+/izY372J5fQlx0FOuB19fv5YrhPYiKqrW6lyWaxqv53odKrZ+hsXLfxoTcyuxCLnn4Ax58aytx0VHcMLYv00f2BOCOReuY89fV7P3y+KknWUG/sHf06FGuvPJKRowYwfDhw3nppZe4//77GTNmDMOHD+e2226rmejLlClTuOuuu8jIyGDIkCFkZmZy7bXXkpaWxr333gs4C/sMHjyYG2+8kSFDhjBz5kyOHTt22uu+8847nHfeeYwaNYpZs2Zx5MgRwBkuO3ToUM455xx+/ON6i040mSUIY0Ls60//h8Q2sbx2+wRe+e753Hf1UB64Jh2AB69L55G9X6PnH7pT8XDzF5wzTffWW2/Rs2dPPvvsMzZu3Mi0adO44447yMzMZOPGjRw/fpzXX3/9xPFxcXFkZWXxne98hxkzZrBgwQI2btzIM888Q1FREeBUav3e977Hli1bSExM5PHHHz/lNQsLC3nggQdYunQpa9euJSMjg0ceeYSioiJeffVVNm3axPr1608knVCzBGFMiLy1cR8A55+dzOI7JjCyd8fTjrl+TB96SQFjY/+P2JK801sSJmylp6fz7rvv8tOf/pQVK1aQlJTE+++/z7hx40hPT2fZsmUnSnqDU2+p5rxhw4bRo0cP4uPjOeuss9i92ylT17t3byZMmADATTfdxMqVK095zdWrV7N582YmTJjAyJEjWbhwIV988QVJSUkkJCRw66238s9//pO2bdu6cs3WB2FMiNz98mdMi4Knbs4gITa63mOfvXUce55IodcfujsbbChs2Bs4cCBr165lyZIl3HvvvVx88cUsWLCArKwsevfuzbx58ygtLT1xfHx8PABRUVEnHtc8rykJXlMmvEbt56rK1KlTeeGFF06L55NPPuG9997jlVde4bHHHmPZsmUhu9YTsYb8KxrTypRWOGWW28c7/281lBwABnXvwLY5H9GvdBG/zlhlHdMRYO/evbRt25abbrqJe+65h7Vr1wKQnJzMkSNHTiwP2hi7du1i1apVACxatIiJEyeesn/8+PF89NFH5OTkAE4/yPbt2zly5AjFxcVcccUVzJ8/n88+++wMry4wa0EYc4bmL93Oz4BHvjoSGlF/7aLB3bj5vL78beXnXDy4K+cPSG74JHOS/2TEUH29emzYsIF77rmHqKgoYmNjeeKJJ3jttdcYPnw43bt3P7HiW2MMGjSIBQsW8M1vfpOhQ4fy3e9+95T9KSkpPPPMM8yZM4eysjIAHnjgATp06MCMGTMoLS1FVXnkkUca/drBkJpe90iXkZGhoV6P1ZgGzUvirLJF7Iy/4eRw17qGXtbs8zumtKKKS+d/SFxMFG/+YBKxv+4U2qGbLciWLVsYMmSI12GETG5uLldddRUbN25sttcM9D0UkTWqGnDChd1iMqaJav656tg2rslfIyE2mvuuGkrOgSMs/Dg3RJEZExqWIIxpomVbDwDww0vSTt7uaEJn88VDujJlUAp/XJod6hBNGOvXr1+zth6awvogjGmC6mrlwbe2cjEwZ2wfOK/pncwiwk+nDebyP66A0K8736Ko6mkjfUxwmtKdYC0IY5rg3S35bM93ZrTGRgf5NqqnlTGkRyJXntMDgINHy0MWZ0uSkJBAUVFRk/7QtXaqSlFREQkJjfsPxFoQxjSSqrLg/Rz6dmkLRxtxYgNDWX94cRpsh79/9Dk/unTQmQXZAqWmppKXl0dBQYHXoUSkhIQEUlPrXX7nNJYgjGmklTmFrM8r5n+uTa+1HNaZSevWAYDnVn/B7RcOCGo+RWsSGxtL//79vQ6jVbFbTMY00pMf7qRbYjzXjurlytc/dKyCf67d48rXNqYxLEEY0wg7C46wIruQm8b1JT7Gnf/w03sl8beVO6mutnvtxluuJggRmSYi20QkR0TmBtg/WUTWikiliMz02z5SRFaJyCYRWS8i17sZpzHBenb1F8RGC7PHulc76daJ/dlRcJQVOYWuvYYxwXAtQYhINLAAuBwYCswRkaG1DtsF3AIsqrX9GPB1VR0GTAP+ICId3YrVmGAcK6/klTV5XD68Bykd4hs+oYkuT+9O53ZxvPiJrQ1hvOVmC2IskKOqO1W1HHgRmOF/gKrmqup6oLrW9u2qmu17vBc4AKS4GKsxDfrXp3spKa3k6+f1dfV14mOiuW5UL97dnE9BSZmrr2VMfdxMEL2A3X7P83zbGkVExgJxwI4A+24TkSwRybKhb8Zt/5u1m7Su7Rndt5PrrzV7bB8qq5VX1uS5/lrG1CWsO6lFpAfwLPANVa2uvV9Vn1TVDFXNSEmxBoZxz86CI6zd9SUzR6c2y0zes587j9yEG/ju8lHo/OGuv54xgbiZIPYAvf2ep/q2BUVEEoE3gJ+r6uoQx2ZMo/xz7R6iBL5yrjtDW09TvIt/Tt9Ev9JFSPHuho83xgVuJohMIE1E+otIHDAbWBzMib7jXwX+oaqNX4XDmBCqrlb+uTaPSWkpdE1svmJJlw7rTkJsWDfyTQvn2m+fqlYCdwBvA1uAl1V1k4jcLyLTAURkjIjkAbOAv4hIzYKuXwUmA7eIyKe+j5FuxWpMfVbtLGJvcSnXjfaVKZif7lvToWnVW4PVPj6GqUOdJUkrqk67w2qM61wttaGqS6hVjEBV7/N7nIlz66n2ec8Bz7kZmzHBenXdHjrEx3Dp0G7OhuJdzbaozzUje8I2+HB7ARcP6dYsr2lMDWu/GlOP8spq3t6037nd89hI11sNtU0e6Ay+eO3Tvc32msbUsARhTD1W5hRQUlrJVef0ONlyaKAq6xmpVRI8NjqKQ7Hd+dO2KVTbaCbTzKyaqzH1eH39PhITYpgwILl5XjBA8sm5cRWz/ryKXG5onhiM8bEWhDF1KKus4t1N+SyP/T5xD3Rq1ltL/kb36eRqaQ9j6mItCGPqsDK7kJKySjrL/mbrlA4kKkq4bFg3+BSOl1fRJs7WiTDNw1oQxtThjfX7SGoT63UYAFw+3FmO9IPtVlLGNB9LEMYEUFZZxbub808ObfXY2P6dAXhr4z6PIzGtiSUIYwJYvfMgJWWVTBve3etQAGc0E8B7Ww5QVlnlcTSmtbAEYUwASzfn0yY2uvlGLwWppKySj3OKvA7DtBKWIIypRVVZuiWfSWnJJMSGV4dwh/gY3tq43+swTCthCcKYWjbtPcy+4lKmhkn/g7/Jg1L4wcZrncl089O9Dse0cJYgjKnl3c35iMBFg7t6HcqpkvqwYPuFVKvy2a1fODO7jXGRJQhjalm6JZ/RfTrRpX2YTU67awMH7ylgcvmjvLf1gNfRmFbAEoQxfvZ+eZxNew+H5e0lgM7t4hjVpxPLtuZ7HYppBSxBGONn6RbnD+8lYZogAC4a0pWNew57HYZpBSxBGOPn8qWXkptwA2cv6BW2ncAXDw7f5GVaFksQxviUVlSRUpXPr0Z/fLL2UjOv/xCMgd3ak9qpjddhmFbAivUZ45OVe4iJwKQ03+Q4N9d9OAMiwsWDu7JnbTK95iU5G5P6hG28JnJZC8IYnxU5TiG8cf27eBxJwy4a0o0JZY/y/uxsp7VjQ16NCyxBGOOzMrsQgHbx4d+wHte/MwmxUVbd1bjK1QQhItNEZJuI5IjI3AD7J4vIWhGpFJGZtfbdLCLZvo+b3YzTmKIjZWzaGzkjgxJioxl/VhdLEMZVriUIEYkGFgCXA0OBOSIytNZhu4BbgEW1zu0M/BIYB4wFfikindyK1ZiPdkReAbzJaSl8XniU3QePeR2KaaHcbEGMBXJUdaeqlgMvAjP8D1DVXFVdD1TXOvcy4F1VPaiqh4B3gWkuxmpauZXZBWGzOFCwLhiUAtgiQsY9biaIXsBuv+d5vm0hO1dEbhORLBHJKiiwN4lpGlVlZXYh558d/p3T/s5Kbkevjm340BKEcUlEd1Kr6pOqmqGqGSkpKV6HYyLUzsKj7C0uZWJaeK390BARYfLAFD6OwNtjJjK4mSD2AL39nqf6trl9rjGNUjN6adKAyPsn44KByRwpq/Q6DNNCuTmeLxNIE5H+OH/cZwM3BHnu28B/+3VMXwr8LPQhmlZvfjo3F+/i0jYp9OiS43U0jXb+gGSio+TkhvnpzpwImzhnQsC1FoSqVgJ34Pyx3wK8rKqbROR+EZkOICJjRCQPmAX8RUQ2+c49CPwaJ8lkAvf7thkTWsW7GK4v00Mj8z5+YkIs5/bueHJD8S6bOGdCxtUZQaq6BFhSa9t9fo8zcW4fBTr3aeBpN+MzBnBu0SR4HUXTTR6YAvlw8Gg5nb0OxrQoEd1JbUwoRAlUJ/YOy8J8wbhgYAp5mkznh1IiMn4TvsK/poAxLktP7UjU7Ru9DqPJhvdKIiP6CS4a3I2HvzrC63BMC2ItCNNqHS6tAGDSgMga3lpbdJQwMS2FD7MLUFWvwzEtiCUI02qt9s0fiLT5D4FMSkumoKSM7flHvA7FtCCWIEyrtTLHmf8wqk/kl/ma4GsF1VyTMaEQVIIQkX+KyJUiYgnFtBg1E+TiYiL/17pXxzacldyOjyxBmBAK9p3xOM4kt2wR+a2IDHIxJmNcl3foGDsLj3odRkhNGJDM6p1FVFTVrn1pTNMElSBUdamq3giMAnKBpSLysYh8Q0QiqwSmMZxsPbQkEwYkc6y8ik93f+l1KKaFCLptLSJdcNZu+BawDvgjTsJ415XIjHHRipxCuiXGex1GSJ13VheiBLvNZEIm2D6IV4EVQFvgalWdrqovqeqdQHs3AzQm1KqrlY9zCk907LYUSW1jSe+VZAnChEywLYi/qupQVf0fVd0HICLxAKqa4Vp0xrhg097DHDpWwaQWMLy1tgkDklm360uvwzAtRLAJ4oEA21aFMhBjmsuKHKcwX0trQQBMHJBMZbVyvG0vp3TI/HSvQzIRrN5SGyLSHWcltzYici5QU1c4Eed2kzHhr6YEts/3gOvadKVrh2zvYnLJqL6diI+J4qEhr3Df1UOdJGFMEzVUi+kynI7pVOARv+0lwH+5FJMxoVVTAhs4Xl7FiF+9w/bY2R4H5Y6E2GjG9OvMxzusH8KcuXoThKouBBaKyHWq+n/NFJMxrvkk9yDlVdXQggdnTxiQzINvbeVASSldvQ7GRLSGbjHdpKrPAf1E5O7a+1X1kQCnGRO2VmYXEBcd+TOn6zNxQDIPAqt2FDHD62BMRGvondLO97k90CHAhzERZUV2IRn9Ir/2Un2G9kykY9vYFjkZ0DSvhm4x/cX3+VfNE44x7ikoKWPr/hJ+Mm2Qs0p6CxUdJZx/dhebD2HOWEO3mB6tb7+qfj+04Rjjnpo/mJMGpMByb2Nx2/lnJ7Nkw/6IXkrVeK+hUUxrmiUKY5rBiuxCOrWNZVjPRK9Dcd3EFjjHwzS/YEYxNZmITMOp2RQNPKWqv621Px74BzAaKAKuV9VcXwHAp3BqPcUA/1DV/zmTWEzrpqqszCng/AHJREVJwydEuL5d2tKrYxso9ToSE8kausX0B1X9oYj8GzhtLUNVnV7PudHAAmAqkAdkishiVd3sd9itwCFVHSAis4EHgeuBWUC8qqaLSFtgs4i8oKq5jbw+YwDIOXCE/MNlEb+8aLBExGlFbISqaiW6FSRFE3oN3WJ61vf590342mOBHFXdCSAiLwIzAP8EMQOY53v8CvCYiAhOMmonIjFAG6AcONyEGIwBnNtL0DKWFw3WhDQnQWzcU8yI3h29DsdEoHqHuarqGt/nD3BqLx0CDgKrfNvq0wvY7fc8z7ct4DGqWgkUA11wksVRYB+wC/i9qh6s/QIicpuIZIlIVkFBQQPhmNZsZU4h/ZPbkdqp9VSIOf/sLoAtQ2qaLthy31cCO4BHgceAHBG53MW4xgJVQE+gP/AjETmr9kGq+qSqZqhqRkpKiovhmEi3emfRqR23SX2cOkVJfbwLymXJ7Z31Lmy4q2mqhm4x1XgYuFBVcwBE5GzgDeDNes7ZA/T2e57K6aPPa47J891OSsLprL4BeEtVK4ADIvIRkAHsDDJeY05xrLzq1NtLd23wLphmlpV7iOPlVbSJi/Y6FBNhgq05UFKTHHx24hTsq08mkCYi/UUkDpgNLK51zGLgZt/jmcAyVVWc20oXAYhIO2A8sDXIWI05TXSUcJ7vlktrU15VTWbuQaeqrZUAN43Q0Cima30Ps0RkCfAyTgfyLJwEUCdVrRSRO4C3cYa5Pq2qm0TkfiBLVRcDfwOeFZEcnL6NmhKbC4C/i8gmnBLjf1fV9U26QmOAEalJJCa04Ap99YiLjuKjnEIm11S1tRLgJkgN3WK62u9xPnCB73EBzuiieqnqEmBJrW33+T0uxUk2tc87Emi7MY1VfKyCJGBiWivto0rqw/bi2fAJLbq/xbijoYly32iuQIxxw6qdhUyDFrm8aFDu2sBjy7L5/TvbWfPjS2idN9lMUwU7iilBRG4XkcdF5OmaD7eDM+ZM1cx/GNmK5wHULK368Y4ijyMxkSbYTupnge44K8x9gDMiqaFOamM8VzMHILaFrwFRn3NSO9IhIcbKf5tGC/ZdM0BVfwEc9dVnuhIY515Yxpy53QeP8UXRMa/D8FxN+e+VOYU4gwSNCU6wCaLC9/lLERmOM1/BVjM0YW2F/cd8wsS0FPZ8eZzcomMnJwnacFfTgGAnyj0pIp2AX+DMXWjve2xM2FqZU0CPpAQo8zoS79XMIl+ZU0j/mkmCNtzVNCCoFoSqPqWqh1T1A1U9S1W71qw2Z0w4qqpWPsopsnURfPr5yn+vzLaaZSZ4wY5i6iIifxKRtSKyRkT+ICI2Ys6ErY17iik+XtGqqrfWp6b898c7iqiqtn4IE5xg+yBeBA4A1+GUxCgEXnIrKGPOVM3oJWtBnDQxLZmS0krW533pdSgmQgSbIHqo6q9V9XPfxwNANzcDM+ZMrMguYFjPRLr4Kpqak+W/rbqrCVawCeIdEZktIlG+j6/i1FgyJuwcLatkzReHmNRay2vUoUv7eIb1TLT1IUzQ6k0QIlIiIoeBbwOLcFZ2K8e55XSb++EZ03irdxZRUaVMtv6H00wckMyaLw5xrLzS61BMBGhoRbkOqpro+xylqjG+jyhVTWyuII1pjBXZhSTERjG6XyevQwk7E9OSqahSPvn8tAUajTlNsPMgEJHpwGTf0+Wq+ro7IRlzZj7MLmD8WV2Ij7EFcmob068zcTFRrMwuZIrXwZiwF+ww198CPwA2+z5+ICL/42ZgxjRF3qFj7Cw4av0PdUiIjWZMv07WD2GCEmwn9RXAVFV9WlWfBqbh1GMyJqzUlNew/oe6TRiQzNb9JVQl9raSG6ZejSlx2dHvsc3RN2FpRbZTXmNA1/ZehxK2auaGvH7R284Kc8W7PI7IhKtgE8R/A+tE5BkRWQisAX7jXljGNF5VtbIyu5BJacmIiNfhhK1hPZPo2DbWihmaBjXYSS0iUUA1MB4Y49v8U1Xd72ZgxjTW+rwvOVxaebL/YX6689+xLbV5ihPlv7Od8t+WSk1dGmxBqGo18BNV3aeqi30fQSUHEZkmIttEJEdE5gbYHy8iL/n2/0dE+vntO0dEVonIJhHZICIJjbkw0/p8uL0QkZMrqFG8y7mFUlO91JwwOS2F/YdL2Z5/xOtQTBgL9hbTUhH5sYj0FpHONR/1nSAi0cAC4HJgKDBHRIbWOuxW4JCqDgDmAw/6zo0BngO+o6rDgCmcXJPCmIBWZBeQ3iuJzu3ivA4l7F0wyGllLd92wONITDgLNkFcD3wPZ7nRLL+P+owFclR1p6rWzL6eUeuYGcBC3+NXgIvFuXl8KbBeVT8DUNUiVa0KMlbTCh0urWDd7i+ZlJbs3Fqal2S3lurRI6kNg7t3YPk2K/9t6hZsghiK0xr4DPgU+BMwrIFzegG7/Z7n+bYFPEZVK4FioAswEFARedtXYvwnQcZpWqlVvjLWk9NS7NZSkC4YlELWFzaj2tQt2ASxEBgCPIqTHIZy8j9/N8QAE4EbfZ+/IiIX1z5IRG4TkSwRySoosP+EWrMPtxfQLi6ac/tYeY1gXTAwhYoqWxvC1C3YUhvDVdW//+B9EdncwDl7gN5+z1N92wIdk+frd0gCinBaGx+qaiGAiCwBRgHv+Z+sqk8CTwJkZGTYb3orpap8sL2A887uQlxMY6b2tG4ZfTvTLs7KkZi6BftuWisi42ueiMg4Gu6DyATSRKS/iMQBs3HWs/a3GLjZ93gmsExVFaeUeLqItPUljgtwSnwYc5odBUfIO3ScCwd39TqUiBIXE3VixJfztjPmVMG2IEYDH4tIzZTLPsA2EdkAqKqeU/sEVa0UkTtw/thHA0+r6iYRuR/IUtXFwN+AZ0UkBziIk0RQ1UMi8ghOklFgiaq+0fTLNC3Zsq3OSJwpgyxBNNaUQV1hp5NkB3Tt4HU4JswEmyCmNeWLq+oSYEmtbff5PS4FZtVx7nM4Q12Nqdf7WwsY1K0DvTq28TqUiHPBoBR4E5ZvK7AEYU4T1C0mVf2ivg+3gzSmLiWlFWTmHrTbS01Uk1RtuKsJxHr0TERbmV1IZbVy4SAr730mPvn8IEfLbJU5cypLECaivb/tAB0SYhjd14a3nonyqmpW7SjyOgwTZixBmIilqry/rYDJA1OIibZf5TPRNi6aD7b73WaqmY1ua0W0avauMhFr097DFJSUcaGNXjpj55/dhfe3HUDnD3cSA9haESb4NamNCTfvnxjemnKytDdYDaYmuHhIN5ZuOYAk7HYSgzFYgjAR7P1tBxiRmkRy+/iT9ZdMk1xso8BMAHaLyUSkg0fLWbf7S5scFyJdExMY0buj12GYMGMJwkSk97bkowpTh3bzOpQW41Lf9/LA4VKPIzHhwhKEiUjvbM6nZ1ICw3omeh1K5EvqA/OSuP2D0eRpMku32CJCxmF9ECbiHC+vYkV2Addn9MZZX8qcEd+6GarKnIfeJ21LPjeMs45+Yy0IE4FWZBewNOoOfrVugo3TDyER4ZIh3ViZU8ixcptVbSxBmAj0zuZ8UqWQil8csnH6ITZ1aDfKK6v5cHuh16GYMGAJwkSUyqpq3tuSD0CszZ4OuTH9OpOYEMNS3/fYtG72DjMRZc0Xhzh0rOLkBl8Hq02OC43Y6CguHNyVZVsPUFVtiwi1dtZJbSLKO5vzifNvOfg6WE3oXDKkG//6dC9rvjjEWK+DMZ6yFoSJGKrKu5vzmTCgi9ehtGgXDu5KXEwUSzbs8zoU4zFLECZibNlXwq6Dx5g6tLvXobRo7eNjmDIwhTc3WoJo7SxBmIjxxoa9RAlcNsxmT7vtynN6kH+4zOswjMcsQZiIoKq8sX4f55+dTJf28V6H0+Jd5LvNZFo3V38DRGSaiGwTkRwRmRtgf7yIvOTb/x8R6Vdrfx8ROSIiP3YzThP+Nu09TG7RMa48p4fXobQKHRJimZzmLONabaOZWi3XEoSIRAMLgMuBocAcERla67BbgUOqOgCYDzxYa/8jwJtuxWgixxsb9hEdJVw2zPofmsuV5zjf63W7D3kcifGKmy2IsUCOqu5U1XLgRWBGrWNmAAt9j18BLhZfcR0RuQb4HNjkYowmnPmWvdT5w323l7rQuV2c11G1GhcPcfp63li/3+NIjFfcTBC9gN1+z/N82wIeo6qVQDHQRUTaAz8FflXfC4jIbSKSJSJZBQUF9R1qIpFvESAp3s2ug8e4ym4vNavEhFgA3ty4z24ztVLh2gs1D5ivqkfqO0hVn1TVDFXNSElJaZ7IjCdi7PaSZ/YVl9ptplbKzZnUe4Defs9TfdsCHZMnIjFAElAEjANmisjvgI5AtYiUqupjLsZrwlB1tRIFTBiQTMe2dnvJC/ExUby2bi+j+3b2OhTTzNxsQWQCaSLSX0TigNnA4lrHLAZu9j2eCSxTxyRV7aeq/YA/AP9tyaF1ysw9CMA15/b0OJLWa+rQbvx7/V7KK6u9DsU0M9cShK9P4Q7gbWAL8LKqbhKR+0Vkuu+wv+H0OeQAdwOnDYU1rdur65xG52XDup/otLbifM3r2lG9+PJYBR9st36+1sbVYn2qugRYUmvbfX6PS4FZDXyNea4EZyLCG+v38VuBtnExJzqtTfOalJZCl3ZxvLouz9YAb2WsmqsJayVllZDgdRStW2x0FO9F30HH7P3O8JGkPlZFt5UI11FMxgDQIykBTeptt5W84ltvo21cDP1KF/HCFRtsFb9WxFoQJiwVlJSRAlxzbi9k2kavw2m9fC2FWFXOeuQDXl23hzn+izRZS6JFsxaECUuLP9sLwLXn1p5babwgIlw3KpVPPj9I7tdWO31B1pJo8SxBmLCjqryc6UzCT+vWweNoTI3rRqUSHSW8mLm74YNNi2AJwoSdtbsOsS2/xOswTC3dkxK4aHBXXlmz2+ZEtBKWIEzYWfSf3bSLi/Y6DBPADWP7UHiknKVb8r0OxTQDSxAmrBQfq+D19XuZYX0PYWnywBR6JiXwwifW/9AaWIIwYeXVdXmUVVZzw1gb0hqOoqOE68f0YUV2odehmGZgCcKEDVXlhU92MyI1ieG9krwOx9Thq2NSiRKvozDNwRKECRurdx5kW34JN4yz1kM465HUhosGOyU3SiuqPI7GuMkShAkbf1v5OZ3bxTFjpPU/hLtvTugHwL8+rV3B37QkliBMWMgtPMp7W/O5cVwfEmJtBFO4O+/sLuRLV65/Ix2dP9zrcIxLLEGYsPDMx7nERAlfG9/3ZFlvq70UtkSED65cRr/SRUixTZxrqawWk/Fc8fEKXs7azdXn9KRrYoKV9Y4Q00f05HdvbYVK34b56SfLb1idphbBEoTx3EuZu3hHbid1S+HJctIm7CXERnPjuL7wEeQcKGGAf2KfZ6PQWgJLEMZTpRVVPPnh59wmhdZqiEBfP89JEI+/v4NHvA7GhJz1QRhPVT48jKzK6yhtZyOXIlGX9vEA/MtXfde0LJYgjGfKKqtoX7qPWd3fJP7Hm7wOx5yBaJs51yJZgjCe+b81zhj6Oy9KQ8T+wESy2WN6A7Dny+MeR2JCydUEISLTRGSbiOSIyNwA++NF5CXf/v+ISD/f9qkiskZENvg+X+RmnKb5lVVW8fjyHAAmpSV7HI05U9+54GwAHn8/x+NITCi5liBEJBpYAFwODAXmiMjQWofdChxS1QHAfOBB3/ZC4GpVTQduBp51K07jjedX7yLvkPPfprUeIlxSH3r+oTuHYrvzYuZudhYc8ToiEyJutiDGAjmqulNVy4EXgRm1jpkBLPQ9fgW4WEREVdepak2v1yagjYjEuxiraUaHSyv407JsJg6wlkOLcNcGmFdM5ffXkxATxUNvb/M6IhMibiaIXoD/FMs837aAx6hqJVAMdKl1zHXAWlUtq/0CInKbiGSJSFZBQUHIAjfu+vPyHRw6VsHcywd7HYoJoZQO8dw2+Wze3Ljf61BMiIR1J7WIDMO57fT/Au1X1SdVNUNVM1JSUpo3ONMk+4tLefqjz5kxsqeV9G6BvjWpP8m+oa+q6nE05ky5mSD2AL39nqf6tgU8RkRigCSgyPc8FXgV+Lqq7nAxTtOMYv90DlujZ/Pwnpu8DsW4oF18DHdPHQjAYpsbEfHcTBCZQJqI9BeROGA2sLjWMYtxOqEBZgLLVFVFpCPwBjBXVT9yMUbTjD7KKaRLZT7zJ2QSU5LndTjGJdf7hrw+8MYWDpdWeByNOROuJQhfn8IdwNvAFuBlVd0kIveLyHTfYX8DuohIDnA3UDMU9g5gAHCfiHzq++jqVqzGfWWVVfziXxsB+O6Us516S1axtUWqmTRXdKSMh63DOqK5WotJVZcAS2ptu8/vcSkwK8B5DwAPuBmbaV5/Xr6TnQVHIcEp8maVPlu+r43vy7fXzoB1hVbdNUKFdSe1aRk27S3mT8uyuXpET69DMc3oR5cNIlUKuTRp8cky4CaiWIIwriqrrOJHL39Gp3Zx/HrGMK/DMc0oMSEWgO35NnEuUlmCMK565N3tPHXoG2RWXEfH36VYn0Mr9LXxfQH4OKfQ40hMY1mCMK55b0s+f/lgJ6k1az3MK7b70K2F3yCEn13hTIi8++XPKDxy2nxXE8ZswSDjit0HjzH4pQnkJhSgSb2xakutjN8/Am19nw8dK+fORet49taxxETb/6aRwH5KJuSOlVfy3efX0IsCvrhzL3LXRq9DMmHgN19JZ9XOIh58a6vXoZggWYIwIVVVrXz/hXVs3nsYgL5d2nkckQkLSX2Y+e9hrOtwN39d8TkvfmKjmiKBJQgTMqrKr1/fzNItB/jl1TZiyfjxVXztVLGfCwam8PPXNvLelnyvozINsARhQkJVeejtbTzzcS6fJt7Nze+MsBFLJqDHbxzFsJ6J3L5oLat3FsH8dKdDe36616GZWixBmDOmqsxfms3jy3cwZ2wfOpbvtxFLpk7t4mN4+pYx9O7Ullv+/okziW5esU2mC0OWIMwZqapWfvXvzTz6XjZfzUjlN9cM9zokE858w1+Tf9+Vt7idfr4+KrvdFJ5smKtpsuPlVfzofz9lyYb9fGtif/7riiFERdmAVlMPv1Zl9LwkXrhnPDwE3/5HFjvjndaoLUEbPqwFYZrk88KjfOXxj3hz437uvXII9141lKg/nmMVWk3wkvrQ6aEUqpN6M3VoNwDm/t8GSiuqPA7M1LAWhGkUVWXxZ3u599WNREcLf79lDFMG+Sqx19xLNiYYvtZEFPBEtcL98FLWbtbtPsQfZ5/LkB6JJ4+dn+78fllV2GZlCcIEbX9xKfe+tpGlW/L5pO0P6Vp9AF7wO8BaDqaJoqIEkvqQyw1QDHl/TmbBlLf51qT+xMdEn/znY54tU9ucLEGYBh0vr+Lpjz7nieU7qKiq5udXDKHrsgPWWjCh5WsZFB0pI/X3XXno7W1c9+HldNcDVq7FI5YgTJ2OlVfyv1l5PLF8B/sPl3LJkG7ce+UQ+iW3g2VeR2daqi7t4wFY+M2xdF90gH6lizi3a0e+v+0AF3ocW2tjCcKc5ouio7yctZvnVu+i+HgFo/t24o+zRzLurC6n3gs2xkUXDEwB4L+/ks6C93P4xt8zyU2AZ1flMuPcXiQ+Mcr6JVxmCcKc/KPvE63JPF7+KJcN7c63J/dndN/OpyYGu7Vk3OZXLvyGcX2YOTqVV9flkf9mV7729gh4G4piuvHxV7Zw9atDvI62xRJV9TqGkMjIyNCsrCyvw4gYVY8MJ/rwbgD2S1fGH/8DAAO7teed4uns/eF+enZsc/KEeUmWGIznVJX1ecW8mLmbdzbtp+hoObkJNzCn59uM7d+Zcf07M7JPR9rG2f++wRKRNaqaEWifq99FEZkG/BGIBp5S1d/W2h8P/AMYDRQB16tqrm/fz4BbgSrg+6r6tpuxRhz///rraGKXV1aTf7iUvV8eZ0fBUXYUHCHnwBG27S9hdflu+pUuIi4mirH9OnPvoBQuHNyVs5Lbwa9wkkPt1zDGYyLCiN4dGdG7Iw9cM5w1XxyChVBSVsGslZeT+tGpq9Ydju/B8iuWkdqpDamd2pDSPt4m4jWCawlCRKKBBcBUIA/IFJHFqrrZ77BbgUOqOkBEZgMPAteLyFBgNjAM6AksFZGBqhp2M2hUFVVQ32OoeQyKs+/ksSe3BToe3/6KKqW8qpryymoqfJ/L/B6XV1ZzSfEuFl2+gcOlFXzn/VH84rWNlJRWUHy8gvzDZeQfLqXoaPmJ114Z/31ukJNvnpKEHvzr1gkM6ZFIXEwd8yVtXoMJY9FRwtj+nQF4/c5JMK+QZbO3sz6vmK37Sti6/zDPHfkW018dQp4mM7bsUVbGf59UKeRAVFfm9llEp7ZxdG4XS8e2cXz9P1fRoXQfx9r05KOrl9M2Lpo2cdHEx0QRExVFTLQQEyXEREc5n6Pk5PZoQRBEQHASmfPZiTVSk5KbLYixQI6q7gQQkReBGYB/gpgBzPM9fgV4TJzv5AzgRVUtAz4XkRzf11sV6iCLjpQx6Xfvn/bHG98f9tp/0Gv+3nt1Z67mFzxPk/mvV51Ww3cS4PX1e+mQEEtimxj+cfhWkqvyIeHkeZUdUqm+68sTpTA6ACPqehG/+7/GhD2/39eLBnfjosHdTuw6Vr6N7YeO03vhOHLlBr6M684PzlrOvdnX83TuJad8mTxNJr1sESv1+0x9eeBp+yaWPXrGoQZKIIKz0f95Y/PJOalJvHjbeWcc32nxutUHISIzgWmq+i3f868B41T1Dr9jNvqOyfM93wGMw0kaq1X1Od/2vwFvquortV7jNuA239NBwLYzCDkZaAmrqreU6wC7lnDVUq6lpVwHnNm19FXVlEA7IronR1WfBJ4MxdcSkay6OmoiSUu5DrBrCVct5VpaynWAe9fiZrG+PUBvv+epvm0BjxGRGCAJp7M6mHONMca4yM0EkQmkiUh/EYnD6XReXOuYxcDNvsczgWXq3PNaDMwWkXgR6Q+kAZ+4GKsxxphaXLvFpKqVInIH8DbOMNenVXWTiNwPZKnqYuBvwLO+TuiDOEkE33Ev43RoVwK3N8MIppDcqgoDLeU6wK4lXLWUa2kp1wEuXUuLmShnjDEmtGzBIGOMMQFZgjDGGBNQq04QIjJLRDaJSLWIZPht7ycix0XkU9/Hn72MMxh1XYtv389EJEdEtonIZV7F2BQiMk9E9vj9LK7wOqbGEJFpvu97jojM9TqeMyEiuSKywfdziKjCZyLytIgc8M29qtnWWUTeFZFs3+dOXsYYrDquxZX3SatOEMBG4FrgwwD7dqjqSN/Hd5o5rqYIeC21ypZMAx73lUGJJPP9fhZLvA4mWH7lZi4HhgJzfD+PSHah7+cQafMHnsH5/fc3F3hPVdOA93zPI8EznH4t4ML7pFUnCFXdoqpnMvs6bNRzLSfKlqjq50BN2RLjvhPlZlS1HKgpN2Oamap+iDNS0t8MYKHv8ULgmuaMqanquBZXtOoE0YD+IrJORD4QkUleB3MGegG7/Z7n+bZFkjtEZL2vaR0RtwF8WsL33p8C74jIGl+Zm0jXTVX3+R7vB7rVd3AECPn7pMUnCBFZKiIbA3zU95/cPqCPqp4L3A0sEpHE5om4bk28lrDXwHU9AZwNjMT5uTzsZayt3ERVHYVzy+x2EZnsdUCh4pugG8lj/l15n0R0LaZgqOolDR912jllQJnv8RpfEcGBgKcdc025FiKgbEmw1yUifwVedzmcUAr7731jqOoe3+cDIvIqzi20QP13kSJfRHqo6j4R6QEc8DqgplLV/JrHoXyftPgWRFOISEpNR66InIVT6mOnt1E1WUSXLfG9cWt8BaczPlIEU24mIohIOxHpUPMYuJTI+lkE4l/q52bgXx7Gckbcep+0+BZEfUTkK8CfgBTgDRH5VFUvAyYD94tIBVANfEdVm6VTqKnquhaPypaE0u9EZCRO8z8X+H+eRtMIdZWb8TispuoGvCrOQgUxwCJVfcvbkIInIi8AU4BkEckDfgn8FnhZRG4FvgC+6l2EwavjWqa48T6xUhvGGGMCsltMxhhjArIEYYwxJiBLEMYYYwKyBGGMMSYgSxDGGGMCsgRhjDEmIEsQxhhjArIEYYyLROR9EZnqe/yAiPzJ65iMCVarnkltTDP4Jc6s/K7AucB0j+MxJmg2k9oYl4nIB0B7YIqqlngdjzHBsltMxrhIRNKBHkC5JQcTaSxBGOMSX4XN53FWLjsiIoGWiTQmbFmCMMYFItIW+CfwI1XdAvwapz/CmIhhfRDGGGMCshaEMcaYgCxBGGOMCcgShDHGmIAsQRhjjAnIEoQxxpiALEEYY4wJyBKEMcaYgP4/pVeMSJMGD+sAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ "
