{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Approximate Bayesian Computation: Synthetic likelihood (parametric approximation)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Florent Leclercq
\n",
"Imperial Centre for Inference and Cosmology, Imperial College London
\n",
"florent.leclercq@polytechnique.org"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"init_cell": true
},
"outputs": [],
"source": [
"import scipy.stats as ss\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from cycler import cycler\n",
"np.random.seed(123457)\n",
"%matplotlib inline\n",
"plt.rcParams.update({'lines.linewidth': 2})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this notebook we try to infer the unknown mean $\\mu$ of a unit-variance Gaussian distribution.\n",
"\n",
"The data consist of $N_\\mathrm{samp}$ realizations. For inference, we run $N$ simulations per value of $\\mu$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def simulator(mu, sigma, Nsamp=10, random_state=None):\n",
" mu, sigma = np.atleast_1d(mu, sigma)\n",
" return ss.norm.rvs(mu[:, None], sigma[:, None], size=(1, Nsamp), random_state=random_state)\n",
"\n",
"def summary_statistic(sim):\n",
" return np.mean(sim, axis=1)\n",
"\n",
"def neg_log_likelihood(mu, Phi_0, Nsamp):\n",
" return 1/2.*np.log(2*np.pi/Nsamp) + Nsamp/2.*(Phi_0 - mu)**2\n",
"\n",
"def sample(mu, sigma, N):\n",
" sims=np.zeros(N)\n",
" for j_ in range(N):\n",
" sim=simulator(mu, sigma)\n",
" sims[j_]+=summary_statistic(sim)\n",
" return np.mean(sims)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"data: d_0=[[1.72558743 1.96316607 3.90707517 1.95831498 1.57465709 0.93464301\n",
" 0.23803412 1.5998699 2.56876862 2.08744352]]\n",
"summary statistic of the data: Phi_0=1.855755990636889\n"
]
}
],
"source": [
"# Set the number of samples and number of simulations per mu\n",
"Nsamp=10\n",
"N=2\n",
"\n",
"# Set the generating parameters that we will try to infer\n",
"mean0 = 2\n",
"sigma0 = 1\n",
"\n",
"# Generate some data\n",
"d_0 = simulator(mean0, sigma0)\n",
"print(\"data: d_0=\"+str(d_0))\n",
"Phi_0 = summary_statistic(d_0)\n",
"print(\"summary statistic of the data: Phi_0=\"+str(Phi_0[0]))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# Set the prior\n",
"mu = ss.uniform(-2,6)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generate simulations"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Generate some data points\n",
"N_sims=100\n",
"mu_sims=mu.rvs(N_sims)\n",
"sims=np.zeros(N_sims)\n",
"for i_ in range(N_sims):\n",
" sims[i_]=sample(mu_sims[i_], sigma0, N)\n",
"S1_sims=neg_log_likelihood(sims,Phi_0,Nsamp)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAF5CAYAAADXmP0gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt4pXV56P3vzQSdcLAigh0YTn3FbmCAQSLJKOWVpLrRilorlt3tLuxXi7vtbm2swEzUjJPLvLXa2hns6Z3Wylhba4si9LC1mkiVVxOb4RCsSFERnMKWEYQiCZXF3PuPrIyZIYeVZK31rMP3c13rWut51nO4syaT3Pkd7l9kJpIkSaqvQ4oOQJIkqR2ZhEmSJBXAJEySJKkAJmGSJEkFMAmTJEkqgEmYJElSAUzCJEmSCmASJkmSVACTMEmSpAKYhEmSJBWgo+gAKvHc5z43Tz755KLDkCRJWtLu3bu/l5nHLHVcUyRhJ598MhMTE0WHIUmStKSIuLeS4+yOlCRJKoBJmCRJUgFMwiRJkgrQFGPCJElqdk8++SR79uzhiSeeKDoUVcnatWtZv349hx566IrONwmTJKkO9uzZw5FHHsnJJ59MRBQdjlYpM3nooYfYs2cPp5xyyoquYXekJEl18MQTT3D00UebgLWIiODoo49eVctmWyVhpVKJgYEBenp6GBgYoFQqFR2SJKmNmIC1ltX+e7ZVd+Tg4CA7duxgamqKyclJIoLh4eGiw5IkSW2orVrCRkdHmZqaAmB6epqRkZGCI5IkSe2qrZKw3t5eOjs7Aejs7KSvr6/giCRJai5HHHHE/tcvfvGLl33+I488wh/+4R8esG8l11mN+WIoQlslYUNDQ/T399Pd3U1/fz/btm0rOiRJkgqRmezbt29V1/jSl7607HPmS4BWcp3VMAkrQEdHB8PDw4yNjTE8PExHR1sNiZMkNZFaTCb79re/zWmnncav/Mqv8MIXvpDvfOc7fPSjH+W8885j48aNvOUtb+Gpp54C4LWvfS3nnnsuZ5xxBjt37pz3ekcccQR//Md/zMaNG9m4cSOnnHIKF1544f7357vG5s2b+eY3v8nGjRu58sor919n1gc+8AE2bNjAhg0b2L59+wFx/9Iv/RJnnHEGL3/5y5menubxxx/nZ37mZzj77LPZsGEDH//4x58WY6UxHOylL30pd911FwAPPfQQGzZsWNZnXZHMrNkD6Af+Bfgq8DFgLXAKMA7cDXwceMZS1zn33HNTkqRm9rWvfW1Zx2/ZsiUPO+ywBLKzszMHBgZWHcM999yTEZFf/vKX98f0qle9Kn/4wx9mZuYv//Iv565duzIz86GHHsrMzKmpqTzjjDPye9/7XmZmHn744fuvN/f1D3/4wzz//PPzxhtv3L9vvmvcc889ecYZZxwQ1+x1JiYmcsOGDfmDH/wgH3vssTz99NPzlltuyXvuuSfXrFmTt956a2ZmXnLJJfnnf/7ned111+Wb3/zm/dd55JFHnvY1VxrDwY4//vh86qmnMjNzdHQ0L7300nmPm+/fFZjICvKkmrWERcTxwK8DXZm5AVgDXAr8NvB7mXkq8H3gTbWKQZKkZlWryWQnnXQSPT09AIyMjLB7925e9KIXsXHjRkZGRvjWt74FwDXXXMPZZ59NT08P3/nOd7j77rsXve5b3/pWent7ufjii/fvW+41br75Zn72Z3+Www8/nCOOOILXve51fPGLXwTglFNOYePGjQCce+65fPvb3+bMM8/kc5/7HFdffTVf/OIX+bEf+7GnXXO5MQDce++9HH/88RxyyEyaNDk5yVlnnbXkectV6+7IDqAzIjqAw4AHgF7guvL7u4DX1jgGSZKaTq0mkx1++OH7X2cml112Gbfddhu33XYbd911F+9+97u56aab+NznPseXv/xlbr/9ds4555xFi5Jee+213HvvvWzdunX/vuVeYzaehTzzmc/c/3rNmjWUSiVe8IIXsHv3bs4880y2bNnC0NDQAeesJAaA22677YCka/fu3c2VhGXmvwG/A9zHTPL1KLAbeCQzZzu29wDH1yoGSZKaVT0mk/X19XHdddfx4IMPAvDwww9z77338uijj3LUUUdx2GGH8fWvf52xsbEFr7F7925+53d+h49+9KP7W46ABa9x5JFH8thjj817rQsuuIBPfepTTE1N8fjjj3P99dfzUz/1Uwve+/777+ewww7jjW98I29/+9u55ZZbDnh/JTEA3H777fuTtbvvvpsbbriBM888c8HjV6qW3ZFHAa9hZgzYccDhwCvmOXTetDciroiIiYiY2Lt3b63ClCSpIdVjMtnpp5/Oe97zHl7+8pdz1lln8bKXvYwHHniAiy66iFKpxFlnncW73vWu/d2X8/n93/99Hn74YS688EI2btzIm9/8ZoAFr3H00Ufzkpe8hA0bNjxtUPwLX/hCLr/8cs477zy6u7t585vfzDnnnLPgve+44479kwqGh4d55zvfecD7K4kBZlrC9u3bx9lnn83Q0BCnnXYau3btquxDXYZYrOlvVReOuAS4KDPfVN7+RWATcAnw45lZiohNwLsz8z8vdq2urq6cmJioSZySJNXDnXfeyWmnnVZ0GKrA85//fG699VaOPPLIJY+d7981InZnZtdS59ZyTNh9QE9EHBYziyv1AV8DPg+8vnzMZcANNYxBkiSpYo899hiHHHJIRQnYatVyTNg4MwPwbwHuKN9rJ3A18LaI+AZwNPChWsUgSZK0HEceeST/+q//Wpd71bRaaWZuBbYetPtbwHm1vK8kSVKja6uK+ZIkSY3CJEySJKkAJmGSJNVJrSoSqBir/fc0CZMkqQ7Wrl3LQw89ZCLWIjKThx56iLVr1674GjUdmC9JkmasX7+ePXv2YAHy1rF27VrWr1+/4vNNwiRJqoNDDz2UU045pegw1EDsjpQkSSqASZgkSVIBTMIkSZIKYBImSZJUAJMwSZKkApiESZIkFcAkTJIkqQBtl4SVSiUGBgbo6elhYGCAUqlUdEiSJKkNtV2x1sHBQXbs2MHU1BSTk5NEBMPDw0WHJUmS2kzbtYSNjo4yNTUFwPT0NCMjIwVHJEmS2lHbJWG9vb10dnYC0NnZSV9fX8ERSZKkdtR23ZFDQ0NEBCMjI/T19bFt27aiQ5IkSW0oMrPoGJbU1dWVExMTRYchSZK0pIjYnZldSx3Xdt2RkiRJjcAkTJIkqQBtmYRZK0ySJBWt7Qbmg7XCJElS8dqyJcxaYZIkqWhtmYRZK0ySJBWtLbsjrRUmSZKKZp0wSZKkKiq8TlhE/GRE3Dbn8e8R8RsR8ZyI+GxE3F1+PqpWMUiSJDWqmiVhmXlXZm7MzI3AucAUcD2wGRjJzFOBkfK2JElSW6nXwPw+4JuZeS/wGmBXef8u4LV1ikGSJKlh1CsJuxT4WPn18zLzAYDy87F1ikGSJKlh1DwJi4hnAK8G/maZ510RERMRMbF3797aBLcIq+pLkqRaqkeJilcAt2Tmd8vb342IdZn5QESsAx6c76TM3AnshJnZkXWI8wBW1ZckSbVUj+7I/8KPuiIBbgQuK7++DLihDjFUZG7r17XXXmtVfUmSVDM1bQmLiMOAlwFvmbP7vcBfR8SbgPuAS2oZw3LMbf3q6Oigo6ODUqlkVX1JklR1NU3CMnMKOPqgfQ8xM1uy4cxdU7JUKnHcccdxwgknWFVfkiRVXVsuW7SQ3t5eJicnmZ6eprOzk8svv9xxYJIkqSZMwuZwTUlJklQvrh0pSZJURYWvHSlJkqSFmYRJkiQVwCRMkiSpACZhkiRJBTAJkyRJKoBJmCRJUgFMwiRJkgpgEiZJklQAkzBJkqQCmIRJkiQVwCRMkiSpAC2VhJVKJQYGBujp6WFgYIBSqVR0SJIkSfPqKDqAahocHGTHjh1MTU0xOTlJRDA8PFx0WJIkSU/TUi1ho6OjTE1NATA9Pc3IyEjBEUmSJM2vpZKw3t5eOjs7Aejs7KSvr6/giCRJkubXUt2RQ0NDRAQjIyP09fWxbdu2okOSJEmaV2Rm0TEsqaurKycmJooOQ5IkaUkRsTszu5Y6rqW6IyvlLEpJklS0luqOrJSzKCVJUtHasiXMWZSSJKlobZmEOYtSkiQVrS27I51FKUmSiubsSEmSpCpydqQkSVIDq2kSFhHPjojrIuLrEXFnRGyKiOdExGcj4u7y81G1jEGSJKkR1bolbAfw6cz8T8DZwJ3AZmAkM08FRsrbkiRJbaVmSVhEPAu4APgQQGb+MDMfAV4D7Coftgt4ba1ikCRJalS1bAn7CWAv8OGIuDUi/jQiDgeel5kPAJSfj61hDJIkSQ2plklYB/BC4I8y8xzgcZbR9RgRV0TERERM7N27t1YxSpIkFaKWSdgeYE9mjpe3r2MmKftuRKwDKD8/ON/JmbkzM7sys+uYY46pYZiSJEn1V7MkLDP/N/CdiPjJ8q4+4GvAjcBl5X2XATfUKgZJkqRGVevZkb8G/EVETAIbgf8XeC/wsoi4G3hZebutlUolBgYG6OnpYWBggFKpVHRIkiSpxmq6bFFm3gbMVzHWxRrnGBwcZMeOHUxNTTE5OUlEMDw8XHRYkiSphhZMwiLidYudmJmfrH447Wl0dJSpqSkApqenGRkZMQmTJKnFLdYSdnH5+VjgxcBoeftC4CbAJKxKent7mZycZHp6ms7OTvr6bCiUJKnVLZiEZeZ/B4iIvwNOn63tVZ7R+Af1Ca89DA0NERGMjIzQ19fHtm3big5JkiTVWGTm4gdEfDUzN8zZPgSYnLuv1rq6unJiYqJet5MkSVqxiNidmfONiT9AJQPzb4qIzwAfAxK4FPj8KuOTJElqa0smYZn5PyPiZ5lZBxJgZ2ZeX9uwJEmSWlulJSq+BJSYaQn7Su3CkSRJag9LFmuNiDcwk3i9HngDMB4Rr691YJIkSa2skpawdwAvyswHASLiGOBzzKwFKUmSpBWoZNmiQ2YTsLKHKjxPkiRJC6ikJezTc2ZHAvw88A+1C0mSJKn1LdmilZlXAv8fcBZwNjOzI6+udWDtwIW7JUlqX5V2K/7/zNQGGym/VhXMLtw9Pj7O+9//fk466SSTMUmS2oSzIws0d+HuUqnE/fffz/bt29m6dWvBkUmSpFqrpCVsdnbkZZn5i8B5wLtqG9bKNFv3Xm9vL52dnQfsm56eZmRkpKCIJElSvVQyML9pZkfOdu9NTU0xOTlJRDA8PFx0WAuaXbj7wx/+MHv37qVUKtHZ2UlfX1/RoUmSpBqrJJn6dER8JiIuj4jLgb+nQWdHzu3ea4YWpY6ODoaHh7nvvvu46qqr6O7upr+/n23bthUdmiRJqrFK1o68MiJ+DngJEDTw2pG9vb1MTk4yPT3dVC1Ks8lYI7faSZKk6qpo7cjM/ATwiRrHsmqz3XsjIyP09fXZoiRJkhpWZObiB0S8Dvht4FhmWsICyMx8Vu3Dm9HV1ZUTExP1up0kSdKKRcTuzOxa6rhKWsLeB1ycmXeuPixJkiRBZQPzv2sCJkmSVF0LtoSVuyEBJiLi48CngP+YfT8zP1nj2CRJklrWYt2RF895PQW8fM52AiZhkiRJK7RgEpaZ/72egUiSJLWTxbojr8rM90XEB5lp+TpAZv56TSOTJElqYYt1R84Oxrc2hCRJUpUt1h35t+XnXfULR5IkqT0s1h35t8zTDTkrM1+91MUj4tvAY8BTQCkzuyLiOcDHgZOBbwNvyMzvLytqSZKkJrdYd+TvVOkeF2bm9+ZsbwZGMvO9EbG5vH11le4lSZLUFBbrjvyn2dcR0QmcmJl3VeGerwFeWn69C7gJkzBJktRmlqyYHxEXA7cBny5vb4yIGyu8fgL/GBG7I+KK8r7nZeYDAOXnYxe47xURMRERE3v37q3wdpIkSc2hkmWL3g2cBzwCkJm3MTOeqxIvycwXAq8AfjUiLqg0sMzcmZldmdl1zDHHVHqa5lEqlRgYGKCnp4eBgQFKpVLRIUmS1PYqWcC7lJmPRsSyL56Z95efH4yI65lJ5r4bEesy84GIWAc8uOwLa1kGBwfZsWMHU1NTTE5OEhEMDw8XHZYkSW2tkpawr0bELwBrIuLUcvHWLy11UkQcHhFHzr5mZtmjrwI3ApeVD7sMuGFFkatio6OjTE1NATA9Pc3IyEjBEUmSpEqSsF8DzmBm8e6/BP4d+I0KznsecHNE3A58Bfj7zPw08F7gZRFxN/Cy8rZqqLe3l87OTgA6Ozvp6+srOCJJklRJd+SxmfkO4B2zOyLiRcA/L3ZSZn4LOHue/Q8BZgF1NDQ0REQwMjJCX18f27ZtKzokSZLaXmQuWI915oCIW4CLM/PfytsXAH+QmWfWIT4Aurq6cmLC1ZMkSVLji4jdmdm11HGVdEe+BfhURPx4RLwSuAZ45WoDlCRJamdLdkdm5j9HxK8D/wg8AbwsMy3cJUmStArLWTvyMOBR4EMRUdHakZIkSZpfPdaOlCRJ0kEqWjtSkiRJ1bVYd+TNmXl+RDzGgd2SAWRmPqvm0UmSJLWoxVrCzi8/H1m/cCRJktrDYi1hz1nsxMx8uPrhSJIktYfFBubvZqYbcr6VuxP4iZpEJEmS1AYW6448pZ6BSJIktZNKKuZLkiSpykzCmlypVGJgYICenh4GBgYolUpFhyRJkiqw5LJFamyDg4Ps2LGDqakpJicniQiGh4eLDkuSJC1hyZawiOiJiCPnbB8ZEd21DUuVGh0dZWpqCoDp6WlGRkYKjkiSJFWiku7IPwJ+MGf78fI+NYDe3l46OzsB6OzspK+vr+CIJElSJSrpjozM3F8xPzP3RYTdmA1iaGiIiGBkZIS+vj62bdtWdEiSJKkCMSe/mv+AiE8CN/Gj1q9fAS7MzNfWNrQf6erqyomJiXrdTpIkacUiYndmdi11XCXdkf8DeDHwb8AeoBu4YnXhSZIktbcluxUz80Hg0jrEIkmS1DYWWzvyqsx8X0R8kJllig6Qmb9e08gkSZJa2GItYXeWnx2MJUmSVGWLrR35t+WXU5n5N3Pfi4hLahqVJElSi6tkYP6WCvdJkiSpQouNCXsF8Erg+Ii4Zs5bzwJcoFCSJGkVFhsTdj8z48FeDeyes/8xoL+WQUmSJLW6xcaE3Q7cHhF/WT7uxMy8q26RSZIktbBKxoRdBNwGfBogIjZGxI2V3iAi1kTErRHxd+XtUyJiPCLujoiPR8QzVhS5JElSE6skCXs3cB7wCEBm3gacvIx7vJUflbsA+G3g9zLzVOD7wJuWcS1JkqSWUEkSVsrMR1dy8YhYD/wM8Kfl7QB6gevKh+wC6rYGpSRJUqOoJAn7akT8ArAmIk4tV9D/UoXX3w5cBewrbx8NPJKZs7Mr9wDHLydgSZKkVlBJEvZrwBnAfwAfA/4d+I2lToqIVwEPZubcmZUxz6FPWxKpfP4VETERERN79+6tIExJkqTmUckC3lPAO8qP5XgJ8OqIeCWwlpn6YtuBZ0dER7k1bD0zpTDmu+9OYCdAV1fXvImaJElSs1qyJSwiXhAROyPiHyNidPax1HmZuSUz12fmycClwGhm/lfg88Dry4ddBtywivibVqlUYmBggJ6eHgYGBiiVrH8rSVI7WbIlDPgb4I+ZGVz/VBXueTXwVxHxHuBW4ENVuGbTGRwcZMeOHUxNTTE5OUlEMDw8XHRYLa1UKjE4OMjo6Ci9vb0MDQ3R0VHJfwFJkqqvkt9Apcz8o9XcJDNvAm4qv/4WMyUv2tro6ChTU1MATE9PMzIyYhJWYwcnvplJRJiUSZIKsdjakc8pv/zbiPgV4HpmBucDkJkP1zi2ltbb28vk5CTT09N0dnbS19dXdEgt7+DEd9euXTz66KO2RkqSCrHYn/27mZm5ODuj8co57yXwE7UKqh0MDQ0REYyMjNDX18e2bduKDqnmiu4OPDjxjQhbIyVJhVls7chTACJibWY+Mfe9iFhb68BaXUdHB8PDw3X5pV908jOr6HFwBye+Tz31FNdcc42tkZKkQlTym/hLwAsr2KcGVXTyM6vocXAHJ76lUok1a9a0VWukJKlxLFiiIiJ+PCLOBToj4pyIeGH58VLgsLpFqFWbL/kpQm9vL52dnQAN0fI0m5SNjY0xPDzsoHxJUl0t9lvnPwOXM1NQ9QNz9j8GDNQwJlXZciYB1LLrsh3HwUmStJDIXLwYfUT8XGZ+ok7xzKurqysnJiaKDKGplUoltm7dekDys1BiNTAwsL/rsrOzk/7+fgerS5K0DBGxOzO7ljyugiTsmcDPASczp+UsM4dWGWPFTMLqp6enh/Hx8f3b3d3djI2NFRiRJEnNpdIkrJIFvG8AXgOUgMfnPNSCGm3cliRJraqSwT7rM/OimkeihlCLcVuNUiJDkqRGUkl35E7gg5l5R31Cejq7I5ub48wkSe2kmt2R5wO7I+KuiJiMiDsiYnL1IapdNEqJDEmSGkklSdgrgFOBlwMXA68qP0sVOXic2SGHHEJPTw8DAwOUSqWCo5MkqRiVDMz5NeDPMvNrtQ5GradUKrFv3z6OOuoonv3sZ3PSSSdx2223MT097aLZkqS2VklL2NeBP4mI8Yj4HxHxY7UOSq1jcHCQD37wg9x///088sgj3HfffUxPTwN2TUqS2tuSSVhm/mlmvgT4RWZqhU1GxF9GxIW1Dk7NoVQqMTAwMG8X48HjwTLTEhiSJFFZdyQRsQb4T+XH94DbgbdFxFsy89IaxqcmsNgC4QcvmXTZZZdxyCGHuHSRJKntVVKi4gPMDMQfBT6UmV+Z895dmfmTtQ3REhWNbrEq+8tZMkmSpFZQaYmKSn4bfhV4Z2ZOzfPeecuOTC1nsQXCOzo6GB4edvC9JEkHqWRg/l1AAETEGyPiAxFxEkBmPlrL4NQchoaG6O/vp7u7m/7+/lV1MS42vmw5x0iS1OgqaQn7I+DsiDgbuAr4EPAR4P+uZWBqHtVs7VpsfNlyjpEkqdFV0hJWypmBY68BdmTmDuDI2oaldlVJdX0r8EuSWkElSdhjEbEFeCPw9+WZkofWNiy1q4Or689XwqKSYyRJanSVdEf+PPALwJsy839HxInA+2sbltrV0NAQEbFoCYtKjpEkqdEtWaLigIMjXpWZf1fDeOZliQpJktQsKi1RUUl35FxDK4xHakjOtJQkFWW5VTOjJlGoLZVKJQYHBxkdHaW3t5ehoaG6F3J1pqUkqSjLbQl7S6UHRsTaiPhKRNweEf8SEdvK+08pLwZ+d0R8PCKescwY1CJmE6Dx8XG2b9/O1q1b6x6DMy0lSUVZMgmLiNfNPoD15dd9EXHsEqf+B9CbmWcDG4GLIqIH+G3g9zLzVOD7wJtW+TWoSTVCAuRMS0lSUSrp+3kTsAn4fHn7pcAY8IKIGMrMP5/vpHJtsR+UNw8tPxLoZWa2JcAu4N3MFIRVm1lsuaN6caalJKkolSRh+4DTMvO7ABHxPGaSpm7gC8C8SVj52DXAbuD5wB8A3wQeyczZ0c97gONXHL2aWiMkQK5tKUkqypIlKiLijsw8c852AHdk5oaIuDUzz1nyJhHPBq4HBoEPZ+bzy/tPAP5h7vXnnHMFcAXAiSeeeO699967jC9LkiSpGJWWqKikJeyLEfF3wN+Ut18PfCEiDgceqSSYzHwkIm4CeoBnR0RHuTVsPXD/AufsBHbCTJ2wSu4jSZLULCqZHfmrwIeZGVx/DjPjuH41Mx/PzAsXOikijim3gBERncBPA3cyM7bs9eXDLgNuWHn4kiRJzWnJJKw8wP5mYBT4HPCFrKzM/jrg8xExCfwz8Nlytf2rgbdFxDeAo4EPrTR4tT6LqUqSWtWS3ZER8QZm1oq8iZlirR+MiCsz87rFzsvMSWZazg7e/y3gvBVFq7ZjMVVJUquqpDvyHcCLMvOyzPxFZhKod9U2LGlGI9QSkySpFipJwg7JzAfnbD9U4XnSqh1cTPXCCy+0e1KS1BIqmR356Yj4DPCx8vbPA/9Qu5DUila6TuTBtcSeeuoprrnmGrsnJUlNb8k6YQAR8XPAS5gZE/aFzLy+1oHN1dXVlRMTE/W8papsYGBg/9iuzs5O+vv7V5Q89fT0MD4+vn+7u7ubsbGxaoYqSWpBK20MWIlq1gkjMz8BfGLVUaltzTe2ayVJWCMsdSRJaj6NONFrwSQsIh5jZq3Hp73FTOWKZ9UsKrWcaiVPjbDUkSSp+VSrMaCaFkzCMvPIegai1lat5Mm1HiVJK9GIPSkVjQkrmmPCJEnSapRKJbZu3XpAY0DRY8JMwiRJkqqo0iTMel+SJEll9VwurzbtcJIkSU3o4FmUmUlE1KS0RSVrR843S/JRYAL4zfJakJIkSU3v4FmUu3bt4tFHH61JaYtKuiM/AFwJHA+sB94O/AnwV8CfVSUKSZKkBnDwcnkRUbM1jCtpT7soM7vnbO+MiLHMHIqIgapFIkmSVLCFlsurRWmLSpKwfRHxBuC68vbr57zX+FMrJUmSKnRwPcpSqcSaNWtqUiR8yRIVEfETwA5gEzNJ1xjQD/wbcG5m3ly1aBZgiQpJktQsqlaiIjO/lZkXZ+ZzM/OY8utvZOZ0PRIwSatTz+nWkqTKLZmERcQLImIkIr5a3j4rIt5Z+9Ck5tVIic/sdOvx8XG2b9/O1q1bC4tFkvQjlcyO/BNgC/AkQGZOApfWMiip2TVS4jPforWSpOJVkoQdlplfOWif/RnSIhop8Tl4unUjLForSapsduT3IuL/ojwTMiJeDzxQ06ikJtfb28vk5GRNpjQv18HTras5s0eStHKVzo7cCbwY+D5wD/DGzPx2zaMrc3akKlEqlRgcHKzJ0hIriWXr1q0HJD5FxSJJqq9KZ0cumYTNueDhwCGZ+dhqg1sukzBVYmBgYP96X52dnfT391dtaQlJkipVtRIVEfHMiPgF4K1Af0QMRsRgNYKUqqmRxmFJkrSUSgbm3wC8hpnB+I/PeUgNxQHokqRmUskglfWZeVHNI5FWyQHokqRmUsnA/J3ABzPzjvqE9HSOCVM9NdIAf0lS86l0TFglv1nOBy6PiHuA/wACyMw8a4kATgA+Avw4sA/YmZk7IuI5wMeBk4FvA2/IzO9XEIdUF7OFVqemppgZROVcAAAQRUlEQVScnCQiHOAvSaq6SsaEvQI4FXg5cDHwqvLzUkrAb2bmaUAP8KsRcTqwGRjJzFOBkfK21DAc4C9JqodKFvC+d75HBec9kJm3lF8/BtwJHM/MIP9d5cN2Aa9defhS9TXCAP9GWntSklQbdRnoEhEnA+cA48DzMvMBmEnUIuLYesQgVaoRBvjbJSpJ82ulcbsVF2td8Q0ijgD+CRjOzE9GxCOZ+ew5738/M4+a57wrgCsATjzxxHPvvXfJxjepZfT09DA+Pr5/u7u7m7GxsQIjkqTG0AyFuatWrHWVQRwKfAL4i8z8ZHn3dyNiXfn9dcCD852bmTszsyszu4455phahik1nFp0idrFKakVtNK43Zq130VEAB8C7szMD8x560bgMuC95ecbahWD1Kyq3SVaKpW44IILGB8fZ9++fXZxSmpavb29TE5OMj093fSFuWvWHRkR5wNfBO5gpkQFwAAz48L+GjgRuA+4JDMfXuxa1gmTVmdgYID3vve9zP3/bhenpGZUKpXYunXrAX+kNtqYsGrWCVuRzLyZmZpi82netFVaQiMOGh0dHT0gAYuIpv7rUVL76ujoYHh4mOHh4Yb8ebsczROp1CQacWZjb28vt99+O0888QQRwaZNm1zWSVLTa8Sft8tR04H5UjtayaDRWg+aHxoa4m1vexvd3d1s2bKFf/qnf2qqvxYlaT7NPkjfJEyqspXMbJz9a258fJzt27ezdevWBY9dScI223w/NjbG8PCwCZikljD3521EcMghhzTVzG+TMKnKhoaG6O/vp7u7m/7+/oq6/Zbz19xyEjZJamVDQ0Ns3LiRiCAzufXWW5vqZ6JJmFRlK2l1Wk7rWbM3v0tStXR0dLBv3779E4+eeOKJpvqZaJ+E1ACWUxeslWrkSNJqNfPPxJovW1QN1gmTfqQZauRIUr004s/ESuuEmYRJkiRVUUOsHSlJkqT5mYRJkiQVwCRMkiSpACZhkiRJBTAJk9pErZdGktTa/BlSfc5rl9pAqVTiggsuYGxsjMzk9ttvb7qFbiUVq9kXy25EtoRJbWBwcJDx8fGmrSotqXiu1lF9JmFSGxgdHWXfvn37tyOiqapKSyrecpZXm0+1uzNXc71G6Vq1O1JqA3OX9YgINm3aVNHC4pI0aznLq82n2t2Zq7leo3StWjFfagONuKyHpPbS09PD+Pj4/u3u7m7GxsaWfZ1SqcTg4CDXXHMNjz/++IquV61YFmLFfEn7dXR0MDw8zNjYGMPDwyZgkuputd2Zs2ZbseYmYMu9XrViWS1/EkuSpJpbbXfmrLkTBAAOP/xw3vrWty7retWKZbXsjpQkSU1jYGCA7du3Mz09TWdnJ/39/Q1XKqPS7khbwiRJUtNolFasarAlTJIkqYocmC9JktTATMIkSZIKYBImSZJUAJMwqY00ylIdkqQaJmER8WcR8WBEfHXOvudExGcj4u7y81G1ur+kp5stcjg+Ps727dvZunVr0SFJqiH/8GpstWwJuxa46KB9m4GRzDwVGClvS6qTuUUOp6enGRkZKTgiSbXkH16NrWZJWGZ+AXj4oN2vAXaVX+8CXlur+0t6ukZZqkNS7cxt/br22mv9w6uB1btY6/My8wGAzHwgIo6t8/2lttZKRQ4lzW+29Wtqaoo1a9bQ0dFBqVSq6x9es4tsj46O0tvby9DQkGvWzqNhP5GIuAK4AuDEE08sOBqpNcwu5N1oS3xIqp65ww6eeuopjjvuOE444YS6/uE1NxGcnJwkIvy5M496z478bkSsAyg/P7jQgZm5MzO7MrPrmGOOqVuAkiQ1s4OHHVx++eWMjY0xPDxct9Yox59Wpt4tYTcClwHvLT/fUOf7S5LU0hph2EFvby+Tk5P7F9l2/On8arZ2ZER8DHgp8Fzgu8BW4FPAXwMnAvcBl2TmwYP3n8a1IyVJah6lUomtW7cekAi205iwSteOdAFvSZKkKnIBb0mSpAZmEiZJklQAkzBJkqQCmIRJqrpK1qtzTTtJ7a59pipIqptKCjVazFFSu7MlTFLVVVKo0WKOktqdSZikBa20y7CShcJdTFxSu7M7UtKCVtplWEnF7kao6i01MhfBbn0Wa5W0oJ6eHsbHx/dvd3d3MzY2VmBEUvsYGBjY/0dQZ2cn/f39jptsEhZrlbRqdhlKxXHcZOuzXVPSguwylIrjItitz+5ISZKqpJrjuNp9EexmZnekJElVVMls4dnJLOPj42zfvp2tW7eu+H4dHR0MDw9z8803k5mcf/75FjZuMSZhkhqWVfXVSCpJsFYyjmup7/NqJnZqLCZhkmpqNYmUv3zUSCpJsFYymWWx7/NSqcS1117rAP0WZRImqabe+c538v73v5/x8XHe97738a53vavic50dpkZSSYI1NDREf38/3d3d9Pf3VzSZZbHv88HBQfbu3bt/u6OjwwH6LcQkTFJNfeQjH9nf+vXUU0/xkY98pOJzLZGhRlJJgjU7jmtsbIzh4eGKBtIv9n0+Ojp6QOvxscceu/++dtc3P6dZSGpYlshQI5lNsKpdMHWx7/ODy1Rcfvnl+xO7la5oocZhiQpJNbV582Z+93d/l1KpREdHB29/+9v5rd/6raLDkubVaEsFLVamwhUtGlelJSpsCZNUU+95z3tYs2aNrVlqCo3WurRY65vFXJufLWGSWk49WjMarcVE1dFMrUsWc21ctoRJalqrTXDq0ZrRaC0mqo5mal2q1Rg11Y+zIyU1nNXWB6tHaYtGKZ/hDLnqWkmJCWmlTMIkNZzVJjj1KG3RKOUzLGhbPbMtsCMjI/T29tq9p5rzu0tSw1ltl1A9Sls0SvmM+RLWandPtcv4N7uYVW+t979IUtNbbYJTj7EyHR0dbNu2jcxkZGSEzCwkOanHGKZ2SU7mS2i3bdvWFgmoiuF3kqSG0ygDjpdqAWqE5KQeLXL1aG1rBPMltI3wb6zW5ZgwSVrAUuOtFhq7Nt9g+YUG0K92YP1KlslZruWMf6v2RIF6TjyYb1B+o0zAUIvKzLo/gIuAu4BvAJuXOv7cc89NSaq37u7uBPY/uru7D3h/y5Yt2dnZmUB2dnbmpk2bsru7Ozdt2nTA/oGBgdyyZUsedthhB+ybvcZ8+xvJk08+mQMDA9nd3Z0DAwP55JNPLnhstb+eoj+fg/+NG/HfR40HmMhK8qFKDqrmA1gDfBP4CeAZwO3A6YudYxImqQhL/QKem5zMTbwi4mnJ20IJ3VKJ3pNPPplbtmzJ7u7u3LJly6IJUBEOju+8885b9OtZrqU+n1pbTgIqzWrkJGwT8Jk521uALYudYxImqQjL+QV8cLIwm4jNbQmbL6Gbb//cxGbTpk1VawlaKKFbLNF78skn8+qrr85169blunXr8uqrrz7g/SuvvPKApHPdunXzfj2LXWOx+9sSpWbUyEnY64E/nbP934DfX+wckzBJje7gZOHFL37xAcnbQgndfPvndsHN16q2UldffXV2dHQkkB0dHbl58+b9sS+U6G3ZsmX/ObPnzX3/8MMPPyA+4Glf+1LXWOz+tkSpGTVyEnbJPEnYB+c57gpgApg48cQTa/QxSVJ1VDNZWKpVbaXWrVt3wHWPO+64ee83N9E7+L2D3z/00EMXfb+SaxTd5ShVW6VJWBGzI/cAJ8zZXg/cf/BBmbkzM7sys+uYY46pW3CStBLVnKU4dzbi2rVr2bRpU02X0Vls9mNvb+8BX0tHR8cB75977rkHXCsinjZ7sre3lzVr1ix4jUZZfUCqu0oytWo+mKlN9i3gFH40MP+Mxc6xO1JSO6lVF9xC3ZGL3e/JJ5/MzZs353HHHZfHHXdcbt68+YD3p6enc9OmTXnooYfmEUcckVddddXT4l3qGnY5qtVQYUtYzBxbXxHxSmA7MzMl/ywzF61819XVlRMTE3WJTZJaValUYuvWrQcUdrX6u1R9EbE7M7uWPK6IJGy5TMIkSVKzqDQJs2K+JElSAUzCJEmSCmASJkmSVACTMEmSpAKYhEmSJBXAJEySJKkAJmGSJEkFMAmTJEkqgEmYJElSAUzCJEmSCmASJkmSVACTMEmSpAKYhEmSJBUgMrPoGJYUEXuBe5c47LnA9+oQTrvzc64PP+f68HOuDz/n+vBzro9KPueTMvOYpS7UFElYJSJiIjO7io6j1fk514efc334OdeHn3N9+DnXRzU/Z7sjJUmSCmASJkmSVIBWSsJ2Fh1Am/Bzrg8/5/rwc64PP+f68HOuj6p9zi0zJkySJKmZtFJLmCRJUtNoqSQsIt4fEV+PiMmIuD4inl10TK0oIi6JiH+JiH0R4UycKoqIiyLiroj4RkRsLjqeVhURfxYRD0bEV4uOpVVFxAkR8fmIuLP88+KtRcfUiiJibUR8JSJuL3/O24qOqZVFxJqIuDUi/q4a12upJAz4LLAhM88C/hXYUnA8reqrwOuALxQdSCuJiDXAHwCvAE4H/ktEnF5sVC3rWuCiooNocSXgNzPzNKAH+FW/n2viP4DezDwb2AhcFBE9BcfUyt4K3Fmti7VUEpaZ/5iZpfLmGLC+yHhaVWbemZl3FR1HCzoP+EZmfiszfwj8FfCagmNqSZn5BeDhouNoZZn5QGbeUn79GDO/uI4vNqrWkzN+UN48tPxwsHcNRMR64GeAP63WNVsqCTvI/wP8r6KDkJbheOA7c7b34C8ttYCIOBk4BxgvNpLWVO4iuw14EPhsZvo518Z24CpgX7Uu2FGtC9VLRHwO+PF53npHZt5QPuYdzDSF/0U9Y2sllXzOqrqYZ59/0aqpRcQRwCeA38jMfy86nlaUmU8BG8vjoK+PiA2Z6XjHKoqIVwEPZubuiHhpta7bdElYZv70Yu9HxGXAq4C+tP7Gii31Oasm9gAnzNleD9xfUCzSqkXEocwkYH+RmZ8sOp5Wl5mPRMRNzIx3NAmrrpcAr46IVwJrgWdFxEcz842ruWhLdUdGxEXA1cCrM3Oq6HikZfpn4NSIOCUingFcCtxYcEzSikREAB8C7szMDxQdT6uKiGNmKwFERCfw08DXi42q9WTmlsxcn5knM/OzeXS1CRi0WBIG/D5wJPDZiLgtIv646IBaUUT8bETsATYBfx8Rnyk6plZQnlTyP4HPMDOI+a8z81+Kjao1RcTHgC8DPxkReyLiTUXH1IJeAvw3oLf88/i2ciuCqmsd8PmImGTmD7nPZmZVyieo9qyYL0mSVIBWawmTJElqCiZhkiRJBTAJkyRJKoBJmCRJUgFMwiRJkgpgEiZJklQAkzBJkqQCmIRJajsRcVNE/GT59dER4RIvkurOJExSO3o+cHf59VnAHQXGIqlNmYRJaisRcRLwb5m5r7zrLGCywJAktSmTMEntZiMHJl3nYhImqQAmYZLazdnAWoCIOBV4DXZHSiqASZikdrMROCQibgcGgTuBy4oNSVI7iswsOgZJqpuI+AZwTmY+VnQsktqbLWGS2kZEHAnsMwGT1AhsCZMkSSqALWGSJEkFMAmTJEkqgEmYJElSAUzCJEmSCmASJkmSVACTMEmSpAKYhEmSJBXAJEySJKkA/wcvs2s7rMQIugAAAABJRU5ErkJggg==\n",
"text/plain": [
"