Skip to content

Instantly share code, notes, and snippets.

@AshNguyen

AshNguyen/Q5 ABC.ipynb

Created March 5, 2020 06:36
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save AshNguyen/0f2157692b6988d92b5165f3b8ebd45c to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save AshNguyen/0f2157692b6988d92b5165f3b8ebd45c to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
Q5 ABC.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"from functools import partial\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import elfi\n",
"import scipy.stats\n",
"\n",
"from elfi.examples import gauss\n",
"\n",
"# prior hyperparameters\n",
"alpha0 = 1\n",
"beta0 = 3\n",
"mu0 = 0\n",
"lamb0 = .9\n",
"\n",
"fn_simulator = partial(gauss.gauss, n_obs=100)\n",
"\n",
"y_obs = np.array([[-4.951515 , -0.39704242, 2.93966133, -2.61614689, -3.91260062,\n",
" -2.92831672, -1.72750631, -2.6125205 , -3.75656916, -5.10181966,\n",
" -0.8333176 , -3.31976079, -5.87175536, 0.77110933, -4.63390199,\n",
" -1.82295875, -0.07399099, -0.2759834 , 2.29914645, -0.60204064,\n",
" -2.77424661, -1.75384278, -4.40727997, -1.88203221, 0.16178288,\n",
" 1.19626051, -6.55802152, -0.10185907, -2.48141801, 5.63406919,\n",
" 3.10279276, -3.48939553, -4.33681171, -6.51492843, -0.49942249,\n",
" -1.97526955, 0.01324536, -1.61270259, -6.78153185, -3.94747663,\n",
" -5.46071129, -1.015685 , -2.55320709, -2.46808139, -2.54721672,\n",
" 4.33884202, 0.60772748, 1.87865953, 4.69863091, -5.40394326,\n",
" 0.06438283, -1.68361822, 6.24149311, -1.72479885, -4.38586099,\n",
" 0.07351717, -1.71502666, 1.81006243, -2.12492885, -3.38041527,\n",
" -2.88373656, -8.16202789, -2.84201145, -2.82528038, -1.10812402,\n",
" -0.27522718, -0.50822677, -2.37870664, -5.15688687, 1.54306991,\n",
" 0.61995686, 2.74833618, 0.2785206 , -3.30246663, -2.02362492,\n",
" -2.63885949, -0.74579404, 0.31596226, -0.8354919 , -3.79755068,\n",
" -2.38694435, 0.55486644, 0.24803359, -0.90377691, 2.34959998,\n",
" -3.36090201, -0.79326858, 0.94455038, -2.65756609, -4.63980848,\n",
" -3.49274672, 4.89780132, -7.96107034, -0.1507432 , -2.47656375,\n",
" 4.75936734, -4.70119197, -1.9313409 , -6.30239165, -3.38322332]])\n",
"\n",
"# These are written as custom classes to allow for inheritance and so var --> sd\n",
"class CustomPrior_x(elfi.Distribution):\n",
" def rvs(sigma, mu0, lamb0, size=1, random_state=None):\n",
" x = scipy.stats.norm.rvs(loc=mu0, scale=sigma / np.sqrt(lamb0), \n",
" size=size, random_state=random_state)\n",
" return x\n",
" def logpdf(x, sigma, mu0, lamb0):\n",
" return np.log(scipy.stats.norm(loc=mu0, scale=sigma / np.sqrt(lamb0)).pdf(x))\n",
" \n",
"class CustomPrior_sigma(elfi.Distribution):\n",
" def rvs(alpha0, beta0, size=1, random_state=None):\n",
" var = scipy.stats.invgamma.rvs(a=alpha0, scale=beta0, size=size, random_state=random_state)\n",
" return np.sqrt(var)\n",
" def logpdf(x, alpha0, beta0):\n",
" var = scipy.stats.invgamma(a=alpha0, scale=beta0).pdf(x)\n",
" return np.log(var)\n",
" \n",
"m = elfi.new_model()\n",
"\n",
"priors = []\n",
"priors.append(elfi.Prior(CustomPrior_sigma, alpha0, beta0, model=m, name='sigma'))\n",
"priors.append(elfi.Prior(CustomPrior_x, m['sigma'], mu0, lamb0 , model=m, name='x'))\n",
"\n",
"elfi.Simulator(fn_simulator, priors[1], priors[0], observed=y_obs, name='gauss')\n",
"\n",
"sumstats = []\n",
"sumstats.append(elfi.Summary(gauss.ss_mean, m['gauss'], name='ss_mean'))\n",
"sumstats.append(elfi.Summary(gauss.ss_var, m['gauss'], name='ss_var'))\n",
"\n",
"d = elfi.Distance('euclidean', *sumstats, name='d')\n",
"\n",
"rej = elfi.Rejection(d, batch_size=10000, seed=251)\n",
"rej2 = elfi.SMC(d, batch_size=10000, seed=251)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"n = y_obs.shape[1]\n",
"xbar = np.mean(y_obs)\n",
"mu = (lamb0*mu0 + n*xbar)/(lamb0 + n)\n",
"lamb = lamb0 + n\n",
"alpha = alpha0 + n/2\n",
"beta = beta0 + (1/2)*np.sum((y_obs-xbar)**2) + (n*lamb0/(n+lamb0))*(xbar-mu0)**2/2\n",
"pos_mean_x = mu\n",
"pos_mean_sigma = beta/(alpha+1+1/2)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Progress: |██████████████████████████████████████████████████| 100.0% Complete\n",
"Progress: |██████████████████████████████████████████████████| 100.0% Complete\n",
"Progress: |██████████████████████████████████████████████████| 100.0% Complete\n"
]
}
],
"source": [
"result_re1 = rej.sample(10000, threshold=.2)\n",
"result_re2 = rej.sample(10000, threshold=.5)\n",
"result_re3 = rej.sample(10000, threshold=.7)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Progress: |██████████████████████████████████████████████████| 100.0% Complete\n"
]
}
],
"source": [
"schedule = [0.7, 0.2, 0.05]\n",
"result_smc = rej2.sample(10000, schedule)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Progress: |██████████████████████████████████████████████████| 100.0% Complete\n"
]
}
],
"source": [
"schedule = [0.7, 0.5, 0.2]\n",
"result_smc2 = rej2.sample(10000, schedule)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAs8AAAHVCAYAAAAO1xbXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X90V9WB9/v3JkSQEhEIaY1BoMYfQGGwUSm32IswUkdHEaUN2iXMPLjs8zC2A1zXUO+4OqxRB3vn1vFHcR6ZMles2qRSZ7BAqUpR+4DaMUrrjwzKI1ggVJFfkyhJ+LHvH/kSCUnIIfnmm1/v11pZOWeffc53n9Nv4eNmn71DjBFJkiRJLevV0Q2QJEmSugrDsyRJkpSQ4VmSJElKyPAsSZIkJWR4liRJkhIyPEuSJEkJGZ4lSZKkhAzPkiRJUkKGZ0mSJCmh3h3dgJPJzc2Nw4cP7+hmSFLa7Nm8hzPO3EX26dl8+j4c7XWQnWefzgWDL+DTTzcD0K/fBR3cSknqecrKyj6OMQ5pqV6nDs/Dhw/ntdde6+hmSFLaPDrpUa7883v5wrgvsGkGVPXfxJ1/N44X/uIF3nhjEgAXXfRCh7ZRknqiEMIHSeo5bEOSJElKyPAsSZIkJWR4liRJkhLq1GOeJUmSlD6HDh1ix44dVFdXd3RTOkzfvn0pKCggOzu7VecbniVJknqIHTt2kJOTw/DhwwkhdHRzMi7GyJ49e9ixYwcjRoxo1TUctiFJktRDVFdXM3jw4B4ZnAFCCAwePLhNPe+GZ0mSpB6kpwbnY9p6/4ZnSZIkKSHDsyRJkjq1mpoaiouLKSwsZPz48Wzbtq1Rne3bt3P55ZczcuRIRo8ezQMPPNAubTE8S5IkqVNbtmwZAwcOZMuWLcyfP5+FCxc2qtO7d29++MMfUl5eziuvvMKSJUt455130t4WZ9uQJEnqgebNg02b0nvNcePg/vtPXufxxx/nwQcfpLa2lvHjx/Pwww+TlZV10nNWrlzJokWLAJgxYwa33XYbMcYG45fPOusszjrrLABycnIYOXIkO3fuZNSoUW26pxMZniVJkpQR5eXllJaWsmHDBrKzs5k7dy5PPPEEq1evZvPmzY3qL1iwgFmzZrFz506GDh0K1PUwDxgwgD179pCbm9vk52zbto033niD8ePHp/0eDM+SJEk9UEs9xO1h3bp1lJWVcckllwBw8OBB8vLyKC0tPel5McZGZc3NmlFVVcUNN9zA/fffzxlnnNH2Rp/A8CxJkqSMiDEye/ZsFi9e3KC8uLj4pD3PBQUFbN++nYKCAg4fPsyBAwcYNGhQo/qHDh3ihhtu4Fvf+hbXX399u9yD4VmSJEkZMWXKFKZNm8b8+fPJy8tj7969VFZWttjzfO2117J8+XImTJjAihUrmDx5cqOe5xgjc+bMYeTIkSxYsKDd7sHZNiRJkpQRo0aN4u6772bq1KmMHTuWK664gl27drV43pw5c9izZw+FhYXcd9993HvvvQBUVFRw1VVXAbBhwwZ+8pOf8Otf/5px48Yxbtw41qxZk/Z7sOdZkiRJGVNcXExxcfEpndO3b1+eeuqpRuX5+fn1AXnixIlNjo1ON3ueJUmSpIQMz5IkSVJCLQ7bCCH8K/DnwEcxxi+lyv4RuAaoBf438Jcxxv2pY3cAc4AjwHdjjL9KlV8JPABkAT+OMd6b/tuRpMyrKKtIXLe2spajRyI1lbUcPeLIOUnqapL0PD8KXHlC2XPAl2KMY4F3gTsAQgijgJnA6NQ5D4cQskIIWcAS4M+AUcCNqbqSJElSl9FieI4xvgTsPaHs2Rjj4dTuK0BBansaUBJjrIkxbgW2AJemfrbEGN+PMdYCJam6kiRJUpeRjjHP/w34ZWr7bGD7ccd2pMqaK28khHBrCOG1EMJru3fvTkPzJEmSpPRoU3gOIfwtcBh44lhRE9XiScobF8a4NMZ4cYzx4iFDhrSleZIkSeoGampqKC4uprCwkPHjx7Nt27Ym6w0fPpwxY8Ywbtw4Lr744nZpS6vfVgkhzKbuRcIp8bNJ9XYAQ4+rVgAce5OmuXJJkiSpWcuWLWPgwIFs2bKFkpISFi5c2OyqhOvXryc3N7fd2tKq8JyaOWMh8H/GGD897tAzwJMhhPuAfOA84LfU9TyfF0IYAeyk7qXCm9rScEmSJLXevLXz2PTHTWm95rgvjOP+K+8/aZ3HH3+cBx98kNraWsaPH8/DDz9MVlbWSc9ZuXIlixYtAmDGjBncdtttxBgbLdGdCS0O2wgh/BR4GbgghLAjhDAH+BGQAzwXQtgUQvifADHGt4GfAe8Aa4G/ijEeSb1ceBvwK6Ac+FmqriRJknqI8vJySktL2bBhA5s2bSIrK4snnniC4uLi+iW1j/957LHHANi5cydDh9YNYujduzcDBgxgz549ja4fQmDq1KkUFRWxdOnSdrmHFnueY4w3NlG87CT17wHuaaJ8DZD+BcYlSZJ0ylrqIW4P69ato6ysjEsuuQSAgwcPkpeX1+wQjGOaWna7qV7nDRs2kJ+fz0cffcQVV1zBhRdeyNe+9rX0ND7FGfolSZKUETFGZs+ezeLFixuUFxcXs3nz5kb1FyxYwKxZsygoKGD79u0UFBRw+PBhDhw4wKBBgxrVz8/PByAvL4/p06fz29/+1vAsSd1JjHWrDlaUVVD7aS2QfMXC/KL89myaJKXdlClTmDZtGvPnzycvL4+9e/dSWVnZYs/ztddey/Lly5kwYQIrVqxg8uTJjXqeP/nkE44ePUpOTg6ffPIJzz77LN///vfTfg+GZ0mSJGXEqFGjuPvuu5k6dSpHjx4lOzubJUuWMGzYsJOeN2fOHG6++WYKCwsZNGgQJSUlAFRUVHDLLbewZs0aPvzwQ6ZPnw7A4cOHuemmm7jyyhMXyW47w7MkSZIypri4mOLi4lM6p2/fvjz11FONyvPz81mzpu6Vui9+8Yv87ne/S0sbTyYdKwxKkiRJPYLhWZIkSUrI8CxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEDM+SJEnq1GpqaiguLqawsJDx48ezbdu2RnU2b97MuHHj6n/OOOMM7r8//UuQO8+zJEmSOrVly5YxcOBAtmzZQklJCQsXLmy0KuEFF1zApk2bADhy5Ahnn312/aIp6WR4liRJ6oHee28eVVWb0nrN/v3Hcd55J+/tffzxx3nwwQepra1l/PjxPPzww2RlZZ30nJUrV7Jo0SIAZsyYwW233UaMsdES3cesW7eOc889t8WVC1vD8CxJkqSMKC8vp7S0lA0bNpCdnc3cuXN54oknWL16NZs3b25Uf8GCBcyaNYudO3cydOhQAHr37s2AAQPYs2cPubm5TX5OSUkJN954Y7vcg+FZkiSpB2qph7g9rFu3jrKyMi655BIADh48SF5eXqMhGCeKMTYqa67Xuba2lmeeeYbFixe3vcFNMDxLkiQpI2KMzJ49u1GwLS4uPmnPc0FBAdu3b6egoIDDhw9z4MABBg0a1ORn/PKXv+TLX/4yn//859vlHgzPkiRJyogpU6Ywbdo05s+fT15eHnv37qWysrLFnudrr72W5cuXM2HCBFasWMHkyZOb7Xn+6U9/2m5DNsCp6iRJkpQho0aN4u6772bq1KmMHTuWK664gl27drV43pw5c9izZw+FhYXcd9993HvvvQBUVFRw1VVX1df79NNPee6557j++uvb7R7seZYkSVLGFBcXU1xcfErn9O3bl6eeeqpReX5+PmvWrKnf79evH3v27GlzG0/GnmdJkiQpIcOzJEmSlJDhWZIkSUrI8CxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkpQx99xzD6NHj2bs2LGMGzeOV199lUmTJnHOOec0WIb7uuuuo3///vX77777LldddRWFhYWMHDmSb37zm3z44YcZb7/zPEuSJCkjXn75ZVatWsXrr79Onz59+Pjjj6mtrQXgzDPPZMOGDUycOJH9+/c3WDylurqaq6++mvvuu49rrrkGgPXr17N79+52W4a7OYZnSZKkHmjtvLX8cdMf03rNL4z7Alfef2Wzx3ft2kVubi59+vQBIDc3t/7YzJkzKSkpYeLEiTz99NNcf/31vP322wA8+eSTTJgwoT44A1x++eVpbXtSDtuQJElSRkydOpXt27dz/vnnM3fuXF588cX6Y1OmTOGll17iyJEjlJSUNFiF8K233qKoqKgjmtyIPc+SJEk90Ml6iNtL//79KSsr4ze/+Q3r16+nuLiYe++9F4CsrCwmTpxIaWkpBw8eZPjw4RlvXxKGZ0nqoirKKlp9bn5RfhpbIknJZWVlMWnSJCZNmsSYMWNYvnx5/bGZM2cyffp0Fi1a1OCc0aNHN+il7kgO25AkSVJGbN68mffee69+f9OmTQwbNqx+/7LLLuOOO+7gxhtvbHDeTTfdxMaNG1m9enV92dq1a3nzzTfbv9EnMDxLkiQpI6qqqpg9ezajRo1i7NixvPPOOw16mUMI3H777Q1eJAQ4/fTTWbVqFQ899BDnnXceo0aN4tFHHyUvLy/Dd+CwDUmSJGVIUVERGzdubFT+wgsvNFm/qqqqfvvCCy9k7dq17dW0xOx5liRJkhKy51mSaNvLd5KknsOeZ0mSJCkhw7MkSZKUkMM2JCkDsqt/D0A4WgXxSOr3GUAkHK0iu/r3dWXH1T3mUN+xmW6uJKkZ9jxLkiRJCRmeJUmSlDH33HMPo0ePZuzYsYwbN45XX32VSZMmcc455xBjrK933XXX0b9///r9d999l6uuuorCwkJGjhzJN7/5TT788EOee+45ioqKGDNmDEVFRfz6179u1/Y7bEOSJEkZ8fLLL7Nq1Spef/11+vTpw8cff0xtbS0AZ555Jhs2bGDixIns37+fXbt21Z9XXV3N1VdfzX333cc111wDwPr169m9eze5ubn84he/ID8/n7feeouvf/3r7Ny5s93uwfAsSZLUE5XNg32b0nvNgeOg6P5mD+/atYvc3Fz69OkD0GAlwZkzZ1JSUsLEiRN5+umnuf7663n77bcBePLJJ5kwYUJ9cAa4/PLLG11/9OjRVFdXU1NTU/8Z6eawDUmSJGXE1KlT2b59O+effz5z587lxRdfrD82ZcoUXnrpJY4cOUJJSQnFxcX1x9566y2KiopavP7Pf/5zLrroonYLzmDPsyRJUs90kh7i9tK/f3/Kysr4zW9+w/r16ykuLubee+8FICsri4kTJ1JaWsrBgwcZPnz4KV377bffZuHChTz77LPt0PLPGJ4lSZKUMVlZWUyaNIlJkyYxZswYli9fXn9s5syZTJ8+nUWLFjU4Z/To0Q16qU+0Y8cOpk+fzmOPPca5557bXk0HHLYhSZKkDNm8eTPvvfde/f6mTZsYNmxY/f5ll13GHXfcwY033tjgvJtuuomNGzeyevXq+rK1a9fy5ptvsn//fq6++moWL17MV7/61Xa/B8OzJEmSMqKqqorZs2czatQoxo4dyzvvvNOglzmEwO23397gRUKA008/nVWrVvHQQw9x3nnnMWrUKB599FHy8vL40Y9+xJYtW7jrrrsYN24c48aN46OPPmq3e3DYhiRJkjKiqKiIjRs3Nip/4YUXmqxfVVVVv33hhReydu3aRnXuvPNO7rzzzrS1sSX2PEuSJEkJGZ4lSZKkhAzPkiRJUkKGZ0mSJCkhXxiUpE4uu/r3rTrvUN+xaW6JJMmeZ0mSJCkhe54lSZKUEXv27GHKlCkA/PGPfyQrK4shQ4YA8Nvf/pbTTjutI5uXiOFZkiRJGTF48GA2bdoEwKJFi+jfvz+33357gzoxRmKM9OrVOQdIGJ4lSZJ6onnzIBVk02bcOLj//lM+bcuWLVx33XVMnDiRV199lX//93/nT/7kT9i/fz8AJSUlPP/88/z4xz/mww8/5H/8j//BH/7wB3r16sWDDz7IV77ylfTex0l0zkgvSZKkHuWdd95hzpw5vPHGG5x99tnN1vvud7/L3/zN3/Daa6/xs5/9jFtuuSWDrbTnWZIkqWdqRQ9xezr33HO55JJLWqz3/PPPs3nz5vr9ffv2cfDgQU4//fT2bF69FsNzCOFfgT8HPooxfilVNggoBYYD24Bvxhj3hRAC8ABwFfAp8BcxxtdT58wGji08fneMcXl6b0WS2l9rp42TJJ3c5z73ufrtXr16EWOs36+urq7fjjF26MuFSYZtPApceULZ94B1McbzgHWpfYA/A85L/dwK/DPUh+2/A8YDlwJ/F0IY2NbGS5Ikqfvp1asXAwcO5L333uPo0aP827/9W/2xP/3TP2XJkiX1+5vSPW67pba1VCHG+BKw94TiacCxnuPlwHXHlT8W67wCnBlCOAv4OvBcjHFvjHEf8ByNA7kkSZIEwA9+8AOuvPJKpkyZQkFBQX35kiVL2LBhA2PHjmXUqFH8y7/8S0bb1doxz5+PMe4CiDHuCiHkpcrPBrYfV29Hqqy58kZCCLdS12vNOeec08rmSZIkqTNbtGhR/XZhYWGjHuTi4mKKi4sbnTdkyBBWrFjR3s1rVrpn2whNlMWTlDcujHFpjPHiGOPFxybNliRJkjqD1obnD1PDMUj9/ihVvgMYely9AqDiJOWSJElSl9Ha8PwMMDu1PRtYeVz5rFDnK8CB1PCOXwFTQwgDUy8KTk2VSZIkSV1GkqnqfgpMAnJDCDuomzXjXuBnIYQ5wB+Ab6Sqr6Fumrot1E1V95cAMca9IYS7gP9I1fv7GOOJLyFKkiRJnVqL4TnGeGMzh6Y0UTcCf9XMdf4V+NdTap0kSZLUibg8tyRJkpSQy3NLkiT1UGVl6b1eUVF6r3dMTU0Ns2bNoqysjMGDB1NaWsrw4cMb1Nm+fTuzZs3ij3/8I7169eLWW2/lr//6r9PeFnueJUmS1KktW7aMgQMHsmXLFubPn8/ChQsb1enduzc//OEPKS8v55VXXmHJkiW88847aW+LPc+S1ITy3eVNlucc3npK1xkxcEQ6miNJ3cbjjz/Ogw8+SG1tLePHj+fhhx8mKyvrpOesXLmyflGVGTNmcNtttxFjJITPlhI566yzOOusswDIyclh5MiR7Ny5k1GjRqW1/YZnSZIkZUR5eTmlpaVs2LCB7Oxs5s6dyxNPPMHq1avZvHlzo/oLFixg1qxZ7Ny5k6FD65YM6d27NwMGDGDPnj3k5uY2+Tnbtm3jjTfeYPz48Wm/B8OzJEmSMmLdunWUlZVxySWXAHDw4EHy8vIoLS096Xl1E7o1dHyv8/Gqqqq44YYbuP/++znjjDPa3ugTGJ4lSZKUETFGZs+ezeLFixuUFxcXn7TnuaCggO3bt1NQUMDhw4c5cOAAgwYNalT/0KFD3HDDDXzrW9/i+uuvb5d7MDxLUhewdd+pjbUGqOyd3WT5yCEj29ocSWqVKVOmMG3aNObPn09eXh579+6lsrKyxZ7na6+9luXLlzNhwgRWrFjB5MmTG/U8xxiZM2cOI0eOZMGCBe12D4ZnSZKkHqq9ppZrzqhRo7j77ruZOnUqR48eJTs7myVLljBs2LCTnjdnzhxuvvlmCgsLGTRoECUlJQBUVFRwyy23sGbNGjZs2MBPfvITxowZw7hx4wD4h3/4B6666qq03oPhWZIkSRlTXFxMcXHxKZ3Tt29fnnrqqUbl+fn5rFmzBoCJEyc2OTY63QzPknqmvQ1XBsiu3t1g/1SnpJMk9QwukiJJkiQlZM+zJKVJxc4mCvc33K2uhqNHP/stSepa7HmWJEmSEjI8S5IkSQkZniVJkqSEHPMsSZLUQ5VVlLVc6RQU5bfPxNE1NTXMmjWLsrIyBg8eTGlpKcOHD29Ub/jw4eTk5JCVlUXv3r157bXX0t4Ww7OkbqOirCJx3ROnppMkdV7Lli1j4MCBbNmyhZKSEhYuXNjsqoTr168nNze33drisA1J6mDV1bB1a93vY9sn/lTsbPgjSV3V448/zqWXXsq4ceP49re/zZEjR1o8Z+XKlcyePRuAGTNmsG7duowsiNIUe54lqR1VfNJwsZWaI9VEjlJ7pJrI54hEao5UU/HJVvoPqgZg7ycu0CKpeyovL6e0tJQNGzaQnZ3N3LlzeeKJJ1i9ejWbN29uVH/BggXMmjWLnTt3MnToUAB69+7NgAED2LNnT6Me5hACU6dOJYTAt7/9bW699da034PhWZIkSRmxbt06ysrKuOSSSwA4ePAgeXl5zQ7BOKapXuYQQqOyDRs2kJ+fz0cffcQVV1zBhRdeyNe+9rX0ND7F8CxJkqSMiDEye/ZsFi9e3KC8uLj4pD3PBQUFbN++nYKCAg4fPsyBAwcYNGhQo/r5+fkA5OXlMX36dH77298aniVJktQ1TZkyhWnTpjF//nzy8vLYu3cvlZWVLfY8X3vttSxfvpwJEyawYsUKJk+e3Kjn+ZNPPuHo0aPk5OTwySef8Oyzz/L9738/7fdgeJYkSeqh2mtqueaMGjWKu+++m6lTp3L06FGys7NZsmQJw4YNO+l5c+bM4eabb6awsJBBgwZRUlICQEVFBbfccgtr1qzhww8/ZPr06QAcPnyYm266iSuvvDLt92B4liRJUsYUFxdTXFx8Suf07duXp556qlF5fn4+a9asAeCLX/wiv/vd79LSxpNxqjpJkiQpIcOzJEmSlJDDNiSpm8o5/G6T5dnVh2DvrpOfPCiz4yAlqaswPEtSD7N131Y+3v3xSet8Wt3ydTL9opEkdQaGZ0ndRvnu8sR1cw67ip8k6dQZniXpBBU7O7oFkqTOyvAsSZLUQ1VWlqX1ejk57TOcq6amhlmzZlFWVsbgwYMpLS1l+PDhDeps3ry5wRR477//Pn//93/PvHnz0toWw7MkSZI6tWXLljFw4EC2bNlCSUkJCxcubLQq4QUXXMCmTZsAOHLkCGeffXb9oinp5FR1kiRJypjHH3+cSy+9lHHjxvHtb3+bI0eOtHjOypUrmT17NgAzZsxg3bp1xBibrb9u3TrOPffcFlcubA17niVJkpQR5eXllJaWsmHDBrKzs5k7dy5PPPEEq1evZvPmzY3qL1iwgFmzZrFz506GDh0KQO/evRkwYAB79uwhNze3yc8pKSnhxhtvbJd7MDxLkiQpI9atW0dZWRmXXHIJAAcPHiQvL6/REIwTNdXLHEJosm5tbS3PPPMMixcvbnuDm2B4liRJUkbEGJk9e3ajYFtcXHzSnueCggK2b99OQUEBhw8f5sCBAwwaNKjJz/jlL3/Jl7/8ZT7/+c+3yz0YniWpgwVqyKICqAFIbZ/c0dq2feaRmtNOfjy74X5W35Ft+0BJAqZMmcK0adOYP38+eXl57N27l8rKyhZ7nq+99lqWL1/OhAkTWLFiBZMnT2625/mnP/1puw3ZAMOzJElSj9VeU8s1Z9SoUdx9991MnTqVo0ePkp2dzZIlS1p8sW/OnDncfPPNFBYWMmjQIEpKSgCoqKjglltuYc2aNQB8+umnPPfcczzyyCPtdg+GZ0mSJGVMcXFxg/mYk+jbty9PPfVUo/L8/Pz64AzQr18/9uzZ0+Y2noxT1UmSJEkJGZ4lSZKkhAzPkiRJUkKGZ0mSJCkhXxiUpB6oasvJ57qrzKppsN+rz/767TPHntkubZKkrsCeZ0mSJCkhe54ldWtHa7Y1XX6k+YVIslq45hHyW98gSepEKspaXpTpVOQXtfzn4z333MOTTz5JVlYWvXr14pFHHmHhwoW8//77fPDBB/WLn1x33XU8//zzVFVVAfDuu+8yb9483n33XbKzsxkzZgwPPfRQu60k2BzDsySdoiQrADYnkAUcpW41wZiuJklSl/Dyyy+zatUqXn/9dfr06cPHH39MbW3dMLIzzzyTDRs2MHHiRPbv38+uXbvqz6uurubqq6/mvvvu45prrgFg/fr17N69O+Ph2WEbkiRJyohdu3aRm5tLnz59AMjNzSU/v663eubMmfUrBz799NNcf/319ec9+eSTTJgwoT44A1x++eV86UtfymDr6xieJXVb27ZCRUXTP7t3N/8jSWofU6dOZfv27Zx//vnMnTuXF198sf7YlClTeOmllzhy5AglJSUNViF86623KCrK7FLizTE8S5IkKSP69+9PWVkZS5cuZciQIRQXF/Poo48CkJWVxcSJEyktLeXgwYMMHz68Q9vaHMc8S5IkKWOysrKYNGkSkyZNYsyYMSxfvrz+2MyZM5k+fTqLFi1qcM7o0aMb9FJ3JHueJakLOtmwE4ekSOqsNm/ezHvvvVe/v2nTJoYNG1a/f9lll3HHHXdw4403NjjvpptuYuPGjaxevbq+bO3atbz55pvt3+gT2PMsSZLUQyWZWi6dqqqq+M53vsP+/fvp3bs3hYWFLF26lBkzZgAQQuD2229vdN7pp5/OqlWrmDdvHvPmzSM7O5uxY8fywAMPZLT9YHiWJElShhQVFbFx48ZG5S+88EKT9Y/N8Qxw4YUXsnbt2vZqWmIO25AkSZISMjxLkiRJCRmeJUmSpIQc8yxJatHRmm3120eq61YGq6xs/fVycjrHYgeSdKra1PMcQpgfQng7hPBWCOGnIYS+IYQRIYRXQwjvhRBKQwinper2Se1vSR0fno4bkCRJkjKl1eE5hHA28F3g4hjjl4AsYCbwA+CfYoznAfuAOalT5gD7YoyFwD+l6kmSJEldRluHbfQGTg8hHAL6AbuAycBNqePLgUXAPwPTUtsAK4AfhRBCjDG2sQ2SJElqjb1l6b3eoJaHZN1zzz08+eSTZGVl0atXLx555BEWLlzI+++/zwcffEAIAYDrrruO559/vn66unfffZd58+bx7rvvkp2dzZgxY3jooYf4/e9/z/e+9z1qa2s57bTT+Md//EcmT56c3vs6TqvDc4xxZwjh/wX+ABwEngXKgP0xxsOpajuAs1PbZwPbU+ceDiEcAAYDHx9/3RDCrcCtAOecc05rmydJkqRO5uWXX2bVqlW8/vrr9OnTh48//pja2loAzjzzTDZs2MDEiRPZv38/u3btqj+vurqaq6++mvvuu49rrrkGgPXr17N7925yc3P5xS9+QX5+Pm+99RZf//rX2blzZ7vdQ1uGbQykrjd5BJAPfA74syaqHuueoXU5AAAYUklEQVRZDic59llBjEtjjBfHGC8eMmRIa5snSZKkTmbXrl3k5ubSp0/di8e5ubnk59etcjhz5kxKSkoAePrpp7n++uvrz3vyySeZMGFCfXAGuPzyy/nSl77ERRddVH+N0aNHU11dTU1NTbvdQ1teGPxTYGuMcXeM8RDwNPB/AGeGEI71aBcAFantHcBQgNTxAcDeNny+JEmSupCpU6eyfft2zj//fObOncuLL75Yf2zKlCm89NJLHDlyhJKSEoqLi+uPvfXWWxQVtTwk5Oc//zkXXXRRfThvD20Jz38AvhJC6BfqBqdMAd4B1gMzUnVmAytT28+k9kkd/7XjnSVJknqO/v37U1ZWxtKlSxkyZAjFxcU8+uijAGRlZTFx4kRKS0s5ePAgw4cPP6Vrv/322yxcuJBHHnkk/Q0/TlvGPL8aQlgBvA4cBt4AlgKrgZIQwt2psmWpU5YBPwkhbKGux3lmWxouSZKkricrK4tJkyYxadIkxowZw/Lly+uPzZw5k+nTp7No0aIG54wePbpBL/WJduzYwfTp03nsscc499xz26vpQBvneY4x/l2M8cIY45dijDfHGGtijO/HGC+NMRbGGL8RY6xJ1a1O7Remjr+fnluQJElSV7B582bee++9+v1NmzYxbNiw+v3LLruMO+64gxtvvLHBeTfddBMbN25k9erV9WVr167lzTffZP/+/Vx99dUsXryYr371q+1+D64wKEmS1FMlmFounaqqqvjOd77D/v376d27N4WFhSxdupQZM+pG/IYQuP322xudd/rpp7Nq1SrmzZvHvHnzyM7OZuzYsTzwwAP86Ec/YsuWLdx1113cddddADz77LPk5eW1yz0YniV1XSfMT5pz+N0G+4N7QVb9O8uSpI5WVFTExo0bG5W/8MILTdY/NsczwIUXXsjatWsb1bnzzju5884709bGlrRp2IYkSZLUkxieJUmSpIQMz5IkST1IT58puK33b3iWJEnqIfr27cuePXt6bICOMbJnzx769u3b6mv4wqAk9SC7dyert/foZ9upVW8ldQMFBQXs2LGD3Un/MOiG+vbtS0FBQavPNzxLkiT1ENnZ2YwYMaKjm9GlOWxDkiRJSsjwLEmSJCVkeJYkSZISMjxLkiRJCfnCoCSpVcp3l7f63KzKz7aL8ovS0BpJygx7niVJkqSEDM+SJElSQoZnSZIkKSHDsyRJkpSQ4VmSJElKyNk2JHU6FWUVieplV+9u55ZIktSQPc+SJElSQoZnSZIkKSHDsyRJkpSQ4VmSJElKyPAsSZIkJWR4liRJkhJyqjpJnUplZRmffrI/Ud0+tVsb7B+tbTjFXVbaWiVJUh3DsySpkUG9ttVv5xxpeKxf1WnNnvdp/xHt1CJJ6hwctiFJkiQlZM+zpA5XVlFWv32kupwD+2oSnZdzJNlKhGqb3Scs5Phf/Zqvuy/7s+0RJ+mEPlJdXr9dWdnKhh0nJ6eo7ReRpATseZYkSZISMjxLkiRJCTlsQ1KbVFaWtVypBcf/E74kSZ2ZPc+SJElSQoZnSZIkKSHDsyRJkpSQ4VmSJElKyPAsSZIkJWR4liRJkhIyPEuSJEkJGZ4lSZKkhAzPkiRJUkKGZ0mSJCkhw7MkSZKUkOFZkiRJSsjwLEmSJCVkeJYkSZIS6t3RDZCkllRUNF0+yP/8lyRlmOFZknRKaj6obfZY9dGa+u0DVQ2PDRjTp72aJEkZY7+NJEmSlJDhWZIkSUrIYRuSpA5Vvru8zdfIqqz7XZRf1OZrSdLJ2PMsSZIkJWR4liRJkhIyPEuSJEkJGZ4lSZKkhAzPkiRJUkKGZ0mSJCkhp6qT1GplFWUcqW77NGOSJHUV9jxLkiRJCbUpPIcQzgwhrAgh/GcIoTyEMCGEMCiE8FwI4b3U74GpuiGE8GAIYUsI4fchhC+n5xYkSZKkzGhrz/MDwNoY44XAnwDlwPeAdTHG84B1qX2APwPOS/3cCvxzGz9bkiRJyqhWh+cQwhnA14BlADHG2hjjfmAasDxVbTlwXWp7GvBYrPMKcGYI4axWt1ySJEnKsLb0PH8R2A38fyGEN0IIPw4hfA74fIxxF0Dqd16q/tnA9uPO35EqkyRJkrqEtsy20Rv4MvCdGOOrIYQH+GyIRlNCE2WxUaUQbqVuWAfnnHNOG5onqavoV7W1wf7RI7UN9gf5arMkqZNoy19JO4AdMcZXU/srqAvTHx4bjpH6/dFx9Yced34BUHHiRWOMS2OMF8cYLx4yZEgbmidJkiSlV6vDc4zxj8D2EMIFqaIpwDvAM8DsVNlsYGVq+xlgVmrWja8AB44N75AkSZK6grYukvId4IkQwmnA+8BfUhfIfxZCmAP8AfhGqu4a4CpgC/Bpqq4kSZLUZbQpPMcYNwEXN3FoShN1I/BXbfk8SZIkqSP5Go4kSZKUkOFZkiRJSsjwLEmSJCVkeJYkSZISMjxLkiRJCbV1qjpJkjrckepyACor03O9nJyi9FxIUrdjz7MkSZKUkOFZkiRJSsjwLEmSJCXkmGdJUtoM6rWtfjvnSMNj/apOa/a8T/uPaKcWSVJ62fMsSZIkJWTPs9TDVFaWpe1ax2Y4kCSpp7DnWZIkSUrI8CxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEDM+SJElSQoZnSZIkKSHDsyRJkpSQi6RIyqitWxuXDTzUcL9md2baIknSqbLnWZIkSUrInmdJ7eLAmzVNlldXNC6r6VXbzq2RJCk97HmWJEmSErLnWZLULnafMHb9v/o1X3df9mfbI0a0T3skKR3seZYkSZISMjxLkiRJCRmeJUmSpIQc8ywpbfpVfTaJ89EjTc+gMcj/ZJckdWH+NSZJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEDM+SJElSQoZnSZIkKSHDsyRJkpSQ4VmSJElKyEVSJEndRvnu8rRc59KcorRcR1L3Y8+zJEmSlJA9z5IknaCysizt18yxN1vqFux5liRJkhKy51nqYdI1JlSSpJ7InmdJkiQpIcOzJEmSlJDhWZIkSUrI8CxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEnOdZ6sTaY5UzSZLUevY8S5IkSQkZniVJkqSEDM+SJElSQo55liR1uIGHttZv96s6tXM/7T8iza2RpObZ8yxJkiQlZM+zJCkjaj6oTVSv6tPGZf0LT0tzaySpdex5liRJkhIyPEuSJEkJGZ4lSZKkhAzPkiRJUkJtDs8hhKwQwhshhFWp/REhhFdDCO+FEEpDCKelyvuk9rekjg9v62dLkiRJmZSOnue/BsqP2/8B8E8xxvOAfcCcVPkcYF+MsRD4p1Q9SZIkqctoU3gOIRQAVwM/Tu0HYDKwIlVlOXBdantaap/U8Smp+pIkSVKX0Nae5/uBvwGOpvYHA/tjjIdT+zuAs1PbZwPbAVLHD6TqNxBCuDWE8FoI4bXdu3e3sXmSJElS+rQ6PIcQ/hz4KMZYdnxxE1VjgmOfFcS4NMZ4cYzx4iFDhrS2eZIkSVLatWWFwa8C14YQrgL6AmdQ1xN9Zgihd6p3uQCoSNXfAQwFdoQQegMDgL1t+HxJkiQpo1odnmOMdwB3AIQQJgG3xxi/FUJ4CpgBlACzgZWpU55J7b+cOv7rGGOjnmdJXcPWrY3LBh76bLvGUVeSpG6oLT3PzVkIlIQQ7gbeAJalypcBPwkhbKGux3lmO3y2JKmLa+p1l//q13z9fdkwYkT7tUeSjpeW8BxjfAF4IbX9PnBpE3WqgW+k4/MkSZKkjuAKg5IkSVJChmdJkiQpIcOzJEmSlJDhWZIkSUrI8CxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEDM+SJElSQoZnSZIkKSHDsyRJkpRQ745ugCRJPUFlZVlar5eTU5TW60lKxp5nSZIkKSF7niVJOkH57vK0XGfkkJFpuY6kzsOeZ0mSJCkhe54lNdCvamuiegMPtXNDJEnqhOx5liRJkhIyPEuSJEkJOWxDUrOqttQ2e6xmdwYbIklSJ2HPsyRJkpSQ4VmSJElKyPAsSZIkJWR4liRJkhIyPEuSJEkJOduG1EWka7lgqSuq+aD5mV+qj9ZwoKr5cweM6dMOLZLUU9nzLEmSJCVkz7OURpWVZR3dBEmS1I7seZYkSZISsudZktTlVVQ0f2xv/+aPjRiR/rZI6t4Mz1IPt3Vrw/2Bhz7bdgluSZIactiGJEmSlJDhWZIkSUrI8CxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEDM+SJElSQoZnSZIkKSHDsyRJkpSQ4VmSJElKyPAsSZIkJWR4liRJkhIyPEuSJEkJGZ4lSZKkhAzPkiRJUkKGZ0mSJCkhw7MkSZKUkOFZkiRJSqh3RzdA6u7Kd5d3dBMkSVKaGJ4lSV3aoF7bTnq8z6HTmizflz2iHVojqbtz2IYkSZKUkOFZkiRJSshhG5IktZN0vvMwcsjItF1LUuvZ8yxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEfGFQ6qb6VW1NVG/goXZuiNSJbU32f5NGRjhFtNRj2fMsSZIkJdTqnucQwlDgMeALwFFgaYzxgRDCIKAUGA5sA74ZY9wXQgjAA8BVwKfAX8QYX29b8yVJ6pkqK8vSer2cnKK0Xk/qrtrS83wY+L9ijCOBrwB/FUIYBXwPWBdjPA9Yl9oH+DPgvNTPrcA/t+GzJUmSpIxrdXiOMe461nMcY6wEyoGzgWnA8lS15cB1qe1pwGOxzivAmSGEs1rdckmSJCnD0jLmOYQwHLgIeBX4fIxxF9QFbCAvVe1sYPtxp+1IlZ14rVtDCK+FEF7bvXt3OponSZIkpUWbw3MIoT/wc2BejPG/Tla1ibLYqCDGpTHGi2OMFw8ZMqStzZMkSZLSpk3hOYSQTV1wfiLG+HSq+MNjwzFSvz9Kle8Ahh53egFQ0ZbPlyRJkjKp1eE5NXvGMqA8xnjfcYeeAWantmcDK48rnxXqfAU4cGx4hyRJktQVtGWRlK8CNwNvhhA2pcr+b+Be4GchhDnAH4BvpI6toW6aui3UTVX3l234bEmS2mTgodatkLIv2xVSpJ6s1eE5xvi/aHocM8CUJupH4K9a+3mSJElSR3N5bqkbOX6pYZfdliQp/QzP6tHSvUKXJEnq3tIyz7MkSZLUExieJUmSpIQctiF1czUf1HZ0EyRJ6jbseZYkSZISMjxLkiRJCRmeJUmSpIQMz5IkSVJChmdJkiQpIcOzJEmSlJDhWZIkSUrIeZ4lSd1aW+Y67zPstDS2RFJ3YHiWmlC+u7yjmyBJkjohh21IkiRJCRmeJUmSpIQctiFJ0inauvXUzxkxIv3tkJR5hmepE2rNX8yS1BaVlWVpv2ZOTlHaryl1NMOz1MkNPGSSlpS+F5lHDhmZlutIPZVjniVJkqSEDM+SJElSQoZnSZIkKSHHPKvLaI+XWSRJkk6FPc+SJElSQoZnSZIkKSHDsyRJkpSQY56lLqDmg9qOboIkScKeZ0mSJCkxw7MkSZKUkOFZkiRJSsjwLEmSJCVkeJYkSZIScrYNSZKa0dRMN9VHaxKd2/f8PulujqROwPAsSZLaRWVlWVqvl5NTlNbrSa1heJYyoF/V1lOqP/BQOzVEkiS1ieFZ7SbdPQ6SJEkdzfCsbqV8d3lHN0FSNzeo17ZE9focOq3B/v53YV/2iFP+vBGnfoqkdmR4ltrR1tRoDYdhSJLUPThVnSRJkpSQPc+SJPUg6RreNnLIyLRcR+pq7HmWJEmSEjI8S5IkSQkZniVJkqSEHPMsZUhTy/xKkqSuxfAsSVKGDDx0aquNAvSrgk/7O9mz1FkYniVJ6sQqKmBf9qmf5+IqUvtwzLMkSZKUkD3PUkJbT/1fWyVJUjdjeFankK5J+yVJktqT4Vk6Ba152UdSz9SWGXb6DDstjS3pPiory9J6vZycorReTz2D4VmSJJ2ydP6LoUt9qyvxhUFJkiQpIXueVS/d/xwmSZLU3Rie1Wq+5CdJknoaw3MXZS9x67VlyjmX2JbUVbTmz7qetrCKLyCqNRzzLEmSJCVkeJYkSZISctiGeqTWztdck+Z2SFISmZpjvl/VZ9uf9u9hYzikhAzPkiR1Mh21wEpFxWfb+7KTn9fWsdLpegHd+aKVCYZndQqtebGltT0xA1t1liSpszOEKxMyHp5DCFcCDwBZwI9jjPdmug0dwdkxOhdnzZCkkzuVDgqHe9Rpj7/rncGj88loeA4hZAFLgCuAHcB/hBCeiTG+k8l2JNEZw253/S/qflVbGXioo1shST1bWzoVKjh+qEjy0J2f/9l2Tw7d6loy3fN8KbAlxvg+QAihBJgGdLrw3FFB9cM39jd77MC+k7+uNmBMn0Sfke7FTQ68WdeunCPb2L07rZeWJJ2ijviXtdZ+5tYPjh+jnSx0Hx+426q5wN6pFgFLtSVdHV/p7snuib3tIcaYuQ8LYQZwZYzxltT+zcD4GONtx9W5Fbg1tXsBsDljDWy9XODjjm5ED+Wz7zg++47l8+84PvuO5fPvON392Q+LMQ5pqVKme55DE2UN0nuMcSmwNDPNSY8Qwmsxxos7uh09kc++4/jsO5bPv+P47DuWz7/j+OzrZHqRlB3A0OP2C4CKZupKkiRJnUqmw/N/AOeFEEaEEE4DZgLPZLgNkiRJUqtkdNhGjPFwCOE24FfUTVX3rzHGtzPZhnbSpYaZdDM++47js+9YPv+O47PvWD7/juOzJ8MvDEqSJEldWaaHbUiSJEldluFZkiRJSsjw3AohhLtCCL8PIWwKITwbQmhyyvYQwuwQwnupn9mZbmd3FEL4xxDCf6ae/7+FEM5spt62EMKbqf+NXst0O7ujU3j2V4YQNocQtoQQvpfpdnZXIYRvhBDeDiEcDSE0O1WU3/30O4Vn73e/HYQQBoUQnkv9XfpcCGFgM/WOpL73m0IITkbQBi19l0MIfUIIpanjr4YQhme+lR3HMc+tEEI4I8b4X6nt7wKjYoz//YQ6g4DXgIupm8u6DCiKMe7LdHu7kxDCVODXqZdPfwAQY1zYRL1twMUxxu48mXtGJXn2IYQs4F3gCuqmpvwP4MYYY6dbRbSrCSGMBI4CjwC3xxibDMZ+99MvybP3u99+Qgj/D7A3xnhvKsgNbObP/aoYY//Mt7B7SfJdDiHMBcbGGP97CGEmMD3GWNwhDe4A9jy3wrHgnPI5TljoJeXrwHMxxr2pwPwccGUm2tedxRifjTEeTu2+Qt1c4cqAhM/+UmBLjPH9GGMtUAJMy1Qbu7MYY3mMsSusuNrtJHz2fvfbzzRgeWp7OXBdB7alJ0jyXT7+f5MVwJQQQlML4XVLhudWCiHcE0LYDnwL+H4TVc4Gth+3vyNVpvT5b8AvmzkWgWdDCGWpJd+VXs09e7/3Hc/vfsfwu99+Ph9j3AWQ+p3XTL2+IYTXQgivhBAM2K2X5LtcXyfVqXIAGJyR1nUCmV6eu8sIITwPfKGJQ38bY1wZY/xb4G9DCHcAtwF/d+IlmjjXMTIJtPTsU3X+FjgMPNHMZb4aY6wIIeQBz4UQ/jPG+FL7tLj7SMOz93vfBkmefwJ+91shDc/e734bnOz5n8Jlzkl9978I/DqE8GaM8X+np4U9SpLvco/+vhuemxFj/NOEVZ8EVtM4PO8AJh23XwC80OaG9QAtPfvUy5d/DkyJzQzajzFWpH5/FEL4N+r+GcoA0YI0PPsdwNDj9guAivS1sHs7hT93TnYNv/utkIZn73e/DU72/EMIH4YQzoox7gohnAV81Mw1jn333w8hvABcBBieT12S7/KxOjtCCL2BAcDezDSv4zlsoxVCCOcdt3st8J9NVPsVMDWEMDD1ZvDUVJnaIIRwJbAQuDbG+GkzdT4XQsg5tk3ds38rc63snpI8e+peLDkvhDAihHAaMBPwrfcM8bvfofzut59ngGMzVs0GGv1LQOrv2j6p7Vzgq4Ava7ZOku/y8f+bzKDuZfIe0/NseG6de0MIb4UQfk/dX05/DRBCuDiE8GOAGONe4C7qvoT/Afx9qkxt8yMgh7p/jt4UQvifACGE/BDCmlSdzwP/K4TwO+C3wOoY49qOaW630uKzT419u426/1AsB34WY3y7oxrcnYQQpocQdgATgNUhhF+lyv3ut7Mkz97vfru6F7gihPAedTNA3AsN/84FRgKvpb7764F7nemkdZr7LocQ/j6EcG2q2jJgcAhhC7AA6FFTMzpVnSRJkpSQPc+SJElSQoZnSZIkKSHDsyRJkpSQ4VmSJElKyPAsSZIkJWR4liRJkhIyPEuSJEkJ/f/KwTwk5I7uRgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 864x576 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(12,8))\n",
"plt.hist(result_re1.samples['x'], color='b', bins=30, alpha=0.2, label='e=0.2')\n",
"plt.hist(result_re2.samples['x'], color='g', bins=30, alpha=0.2, label='e=0.5')\n",
"plt.hist(result_re3.samples['x'], color='y', bins=30, alpha=0.2, label='e=0.7')\n",
"plt.hist(result_smc.samples['x'], color='purple', bins=30, alpha=0.2, label='SMC')\n",
"plt.hist(result_smc2.samples['x'], color='orange', bins=30, alpha=0.2, label='SMC2')\n",
"plt.axvline(np.mean(result_re1.samples['x']), color='b', alpha=1, label='e=0.2')\n",
"plt.axvline(np.mean(result_re2.samples['x']), color='g', alpha=1, label='e=0.5')\n",
"plt.axvline(np.mean(result_re3.samples['x']), color='y', alpha=1, label='e=0.7')\n",
"plt.axvline(np.mean(result_smc.samples['x']), color='purple', alpha=1, label='SMC')\n",
"plt.axvline(np.mean(result_smc2.samples['x']), color='orange', alpha=1, label='SMC2')\n",
"plt.axvline(pos_mean_x, color='r', label='True')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAs8AAAHVCAYAAAAO1xbXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X90ldWB//v3JomgBQUM6YChYkVrQqFoVMq3dL4glVodRYRpUEeYubjsKrVT4HqH+v26ZritDs7qtw5qtaO3zIi/mlRqRwuUqSLWXrR2iDL1R74gV7DhRxVBaKgkQNj3j3MI+UkekpOThLxfa2Wd8+xnP8/Zz7MOyYd99tk7xBiRJEmS1LY+Xd0ASZIkqacwPEuSJEkJGZ4lSZKkhAzPkiRJUkKGZ0mSJCkhw7MkSZKUkOFZkiRJSsjwLEmSJCVkeJYkSZISyu3qBhxPfn5+HDFiRFc3Q+p9Nm5MPX7mMx0+1ccfp861aXdqe+zwVs75x40cOnCIP+4dypmfObPDrytJ0omoqKj4MMY4pK163To8jxgxgvXr13d1M6TeZ+LE1OOLL3b4VK+/njrXpGWp7fVLWjnn8xP5w4Y/sHrFt/nrF/+6w68rSdKJCCG8l6SewzYkSZKkhAzPkiRJUkKGZ0mSJCmhbj3mWZIkSZlz6NAhtm3bRk1NTVc3pcv069ePwsJC8vLy2nW84VmSJKmX2LZtGwMGDGDEiBGEELq6OVkXY2T37t1s27aNc845p13ncNiGJElSL1FTU8OZZ57ZK4MzQAiBM888s0M974ZnSZKkXqS3BuejOnr9hmdJkiQpIcOzJEmSurXa2lpKS0sZOXIk48aNY+vWrc3qVFVVMWnSJIqKihg1ahT33ntvp7TF8CxJkqRubenSpQwaNIjNmzczf/58Fi5c2KxObm4u3//+96msrOQ3v/kNDzzwAG+//XbG2+JsG5IkSb3QvHmwYUNmzzl2LCxZcvw6jz/+OPfddx8HDx5k3LhxPPjgg+Tk5Bz3mGeeeYZFixYBMGPGDG699VZijI3GLw8dOpShQ4cCMGDAAIqKiti+fTvFxcUduqamDM+SJEnKisrKSsrLy1m3bh15eXnMnTuXJ554gpUrV7Jx48Zm9RcsWMCsWbPYvn07w4cPB1I9zGeccQa7d+8mPz+/xdfZunUrr7/+OuPGjcv4NRieJUmSeqG2eog7w5o1a6ioqOCSSy4B4MCBAxQUFFBeXn7c42KMzcpamzVj//79TJ8+nSVLlnD66ad3vNFNGJ4lSZKUFTFGZs+ezeLFixuVl5aWHrfnubCwkKqqKgoLCzl8+DD79u1j8ODBzeofOnSI6dOnc+ONN3Ldddd1yjUYniVJkpQVkydPZurUqcyfP5+CggL27NlDdXV1mz3P11xzDcuWLWP8+PEsX76cyy67rFnPc4yROXPmUFRUxIIFCzrtGpxtQ5IkSVlRXFzMnXfeyZQpUxgzZgyXX345O3fubPO4OXPmsHv3bkaOHMk999zD3XffDcCOHTu48sorAVi3bh2PPfYYL7zwAmPHjmXs2LGsWrUq49dgz7MkSZKyprS0lNLS0hM6pl+/fjz11FPNyocNG1YfkCdMmNDi2OhMs+dZkiRJSsjwLEmSJCXksA2pF6rYUXHc/ecfrAZgUxv1jioZVtLhNkmS1BPY8yxJkiQlZHiWJEmSEjI8S5IkSQm1GZ5DCP8aQvgghPBmg7LvhRD+dwjhdyGEn4UQBjbYd3sIYXMIYWMI4csNyq9Il20OIXw785ciSZKkk1FtbS2lpaWMHDmScePGsXXr1mZ1qqqqmDRpEkVFRYwaNYp77723U9qSpOf5EeCKJmXPAZ+NMY4BNgG3A4QQioGZwKj0MQ+GEHJCCDnAA8BXgGLg+nRdSZIk6biWLl3KoEGD2Lx5M/Pnz2fhwoXN6uTm5vL973+fyspKfvOb3/DAAw/w9ttvZ7wtbc62EWN8KYQwoknZLxts/gaYkX4+FSiLMdYCW0IIm4FL0/s2xxjfBQghlKXrZv6KJEmS1KZ5q+ex4Q8bMnrOsX82liVXLDlunccff5z77ruPgwcPMm7cOB588EFycnKOe8wzzzzDokWLAJgxYwa33norMcZGS3QPHTqUoUOHAjBgwACKiorYvn07xcWZ7a/NxJjn/wP4Rfr5WUBVg33b0mWtlTcTQrglhLA+hLB+165dGWieJEmSuoPKykrKy8tZt24dGzZsICcnhyeeeILS0tL6JbUb/jz66KMAbN++neHDhwOpHuYzzjiD3bt3t/o6W7du5fXXX2fcuHEZv4YOzfMcQvifwGHgiaNFLVSLtBzSW1w/Mcb4MPAwwMUXX9z5ayxKorKy8Xbhn1oub9XOxpslTvssSd1eWz3EnWHNmjVUVFRwySWXAHDgwAEKCgooLy8/7nEtLbvdsNe5of379zN9+nSWLFnC6aef3vFGN9Hu8BxCmA38BTA5HruibcDwBtUKgR3p562VS+qh8kil67q6xuXV1cee19WlNiKnpfe1vPDKqXXVLZZLkk4eMUZmz57N4sWLG5WXlpaycePGZvUXLFjArFmzKCwspKqqisLCQg4fPsy+ffsYPHhws/qHDh1i+vTp3HjjjVx33XWdcg3tCs8hhCuAhcB/jzF+3GDXs8CTIYR7gGHAecBvSfVInxdCOAfYTupLhTd0pOGSJEnqWSZPnszUqVOZP38+BQUF7Nmzh+rq6jZ7nq+55hqWLVvG+PHjWb58OZdddlmznucYI3PmzKGoqIgFCxZ02jUkmarux8ArwGdCCNtCCHOAHwADgOdCCBtCCP+SbvRbwE9IfRFwNfCNGGNdjPEwcCvwH0Al8JN0XUmSJPUSxcXF3HnnnUyZMoUxY8Zw+eWXs3PnzjaPmzNnDrt372bkyJHcc8893H333QDs2LGDK6+8EoB169bx2GOP8cILL9SPmV61alXGryHJbBvXt1C89Dj17wLuaqF8FZD5K5CU2NEhE3U1jQcz5zWpF/g4XZ500LMkScmUlpZSWlp6Qsf069ePp556qln5sGHD6gPyhAkTWhwbnWmuMChJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEDM+SJElSQoZnSZIkKSHDsyRJkrq12tpaSktLGTlyJOPGjWPr1q3N6lRVVTFp0iSKiooYNWoU9957b6e0pd3Lc0tSayobTA9dU5N6PFLXfF9DF/wJDtellvauaLKCd0lJ5tsoSeo5li5dyqBBg9i8eTNlZWUsXLiw2aqEubm5fP/73+eiiy6iurqakpISLr/8coqLizPaFsOzpA7bUt36YiqD+qQWXDkS+x637oi6jzkS6/j4cHXmGyhJauadd+axf/+GjJ6zf/+xnHfekuPWefzxx7nvvvs4ePAg48aN48EHHyQnJ+e4xzzzzDMsWrQIgBkzZnDrrbcSY2y0RPfQoUMZOnQoAAMGDKCoqIjt27cbniVJktQzVVZWUl5ezrp168jLy2Pu3Lk88cQTrFy5ko0bNzarv2DBAmbNmsX27dsZPnw4kOphPuOMM9i9ezf5+fktvs7WrVt5/fXXGTduXMavwfAsSZLUC7XVQ9wZ1qxZQ0VFBZdccgkABw4coKCgoNkQjKZaWna7Ya9zQ/v372f69OksWbKE008/veONbsLwLEmSpKyIMTJ79mwWL17cqLy0tPS4Pc+FhYVUVVVRWFjI4cOH2bdvH4MHD25W/9ChQ0yfPp0bb7yR6667rlOuwfAsSZKkrJg8eTJTp05l/vz5FBQUsGfPHqqrq9vseb7mmmtYtmwZ48ePZ/ny5Vx22WXNep5jjMyZM4eioiIWLFjQadfgVHWSJEnKiuLiYu68806mTJnCmDFjuPzyy9m5c2ebx82ZM4fdu3czcuRI7rnnHu6++24AduzYwZVXXgnAunXreOyxx3jhhRcYO3YsY8eOZdWqVRm/BnueJUmSlDWlpaWUlpae0DH9+vXjqaeealY+bNiw+oA8YcKEFsdGZ5o9z5IkSVJChmdJkiQpIcOzJEmSlJDhWZIkSUrI8CxJkiQlZHiWJEmSEjI8S5IkKWvuuusuRo0axZgxYxg7diyvvvoqEydO5FOf+lSjqeauvfZa+vfvX7+9adMmrrzySkaOHElRURFf/epXef/997Pefud5liRJUla88sorrFixgtdee42+ffvy4YcfcvDgQQAGDhzIunXrmDBhAnv37m20eEpNTQ1XXXUV99xzD1dffTUAa9euZdeuXXzyk5/M6jUYniVJknqh1fNW84cNf8joOf9s7J9xxZIrWt2/c+dO8vPz6du3LwD5+fn1+2bOnElZWRkTJkzg6aef5rrrruOtt94C4Mknn2T8+PH1wRlg0qRJGW17Ug7bkCRJUlZMmTKFqqoqzj//fObOncuvfvWr+n2TJ0/mpZdeoq6ujrKyskarEL755puUlJR0RZObsedZkiSpFzpeD3Fn6d+/PxUVFfz6179m7dq1lJaWcvfddwOQk5PDhAkTKC8v58CBA4wYMSLr7UvC8CxJkqSsycnJYeLEiUycOJHRo0ezbNmy+n0zZ85k2rRpLFq0qNExo0aNatRL3ZUctiFJkqSs2LhxI++880799oYNGzj77LPrt7/4xS9y++23c/311zc67oYbbuDll19m5cqV9WWrV6/mjTfe6PxGN2F4liRJUlbs37+f2bNnU1xczJgxY3j77bcb9TKHELjtttsafZEQ4NRTT2XFihXcf//9nHfeeRQXF/PII49QUFCQ5Stw2IbUI1VUtO+4urrU45bqzLVFkqSkSkpKePnll5uVv/jiiy3W379/f/3zCy64gNWrV3dW0xKz51mSJElKyPAsSZIkJWR4liRJkhIyPEuSJEkJGZ4lSZKkhAzPkiRJUkKGZ0mSJGXNXXfdxahRoxgzZgxjx47l1VdfZeLEiXzqU58ixlhf79prr6V///7125s2beLKK69k5MiRFBUV8dWvfpX333+f5557jpKSEkaPHk1JSQkvvPBCp7bfeZ4lSZKUFa+88gorVqzgtddeo2/fvnz44YccPHgQgIEDB7Ju3TomTJjA3r172blzZ/1xNTU1XHXVVdxzzz1cffXVAKxdu5Zdu3aRn5/Pz3/+c4YNG8abb77Jl7/8ZbZv395p12B4lnqAih2NV0Wp3Nu+8+RRmYHWSJJOChXz4KMNmT3noLFQsqTV3Tt37iQ/P5++ffsCNFpJcObMmZSVlTFhwgSefvpprrvuOt566y0AnnzyScaPH18fnAEmTZrU7PyjRo2ipqaG2tra+tfINIdtSJIkKSumTJlCVVUV559/PnPnzuVXv/pV/b7Jkyfz0ksvUVdXR1lZGaWlpfX73nzzTUpKSto8/09/+lMuvPDCTgvOYM+zJElS73ScHuLO0r9/fyoqKvj1r3/N2rVrKS0t5e677wYgJyeHCRMmUF5ezoEDBxgxYsQJnfutt95i4cKF/PKXv+yElh9jeJYkSVLW5OTkMHHiRCZOnMjo0aNZtmxZ/b6ZM2cybdo0Fi1a1OiYUaNGNeqlbmrbtm1MmzaNRx99lHPPPbezmg4YniV1Q5V7G4/xZkfHzlcyrO2P+iRJnW/jxo306dOH8847D4ANGzZw9tln8+abbwLwxS9+kdtvv53rr7++0XE33HADixcvZuXKlVx11VUArF69mrPOOovhw4dz1VVXsXjxYr7whS90+jU45lmSJElZsX//fmbPnk1xcTFjxozh7bffbtTLHELgtttua/RFQoBTTz2VFStWcP/993PeeedRXFzMI488QkFBAT/4wQ/YvHkz3/3udxk7dixjx47lgw8+6LRrsOdZkiRJWVFSUsLLL7/crPzFF19ssf7+/fvrn19wwQWsXr26WZ077riDO+64I2NtbIs9z5IkSVJChmdJkiQpIcOzJEmSlJDhWZIkSUrI8CxJkiQlZHiWJEmSEnKqOkmSJGXF7t27mTx5MgB/+MMfyMnJYciQIQD89re/5ZRTTunK5iVieJYkSVJWnHnmmWzYsAGARYsW0b9/f2677bZGdWKMxBjp06d7DpAwPEuSJPVG8+ZBOshmzNixsGTJCR+2efNmrr32WiZMmMCrr77Kv//7v/O5z32OvXv3AlBWVsbzzz/Pj370I95//32+/vWv8/vf/54+ffpw33338fnPfz6z13Ec3TPSS5IkqVd5++23mTNnDq+//jpnnXVWq/X+9m//lr/7u79j/fr1/OQnP+Hmm2/OYivteZYkSeqd2tFD3JnOPfdcLrnkkjbrPf/882zcuLF++6OPPuLAgQOceuqpndm8eoZnSZIkdblPfOIT9c/79OlDjLF+u6ampv55jLFLv1zosA1JkiR1K3369GHQoEG88847HDlyhJ/97Gf1+770pS/xwAMP1G9vyPS47TbY8yypSw06tAWA3FhDIJAXP2bwocpGdU6rbvnYjwcUdXbzJEld5J/+6Z+44oor+NSnPkVxcTG1tbUAPPDAA3z961/n3/7t3zh8+DCTJk1qFKY7W5vhOYTwr8BfAB/EGD+bLhsMlAMjgK3AV2OMH4UQAnAvcCXwMfDXMcbX0sfMBu5In/bOGOOyzF6KJEmSeopFixbVPx85cmSzHuTS0lJKS0ubHTdkyBCWL1/e2c1rVZJhG48AVzQp+zawJsZ4HrAmvQ3wFeC89M8twA+hPmz/AzAOuBT4hxDCoI42XpIkScqmNnueY4wvhRBGNCmeCkxMP18GvAgsTJc/GlMjvH8TQhgYQhiarvtcjHEPQAjhOVKB/McdvgJJJ5EjBGrIZUuj0rralmvX5bV9xpx+Du2QJGVOe78w+MkY406A9GNBuvwsoKpBvW3pstbKmwkh3BJCWB9CWL9r1652Nk+SJEnKvEzPthFaKIvHKW9eGOPDMcaLY4wXH13rXJIkSeoO2hue308PxyD9+EG6fBswvEG9QmDHccolSZKkHqO9U9U9C8wG7k4/PtOg/NYQQhmpLwfuizHuDCH8B/CPDb4kOAW4vf3NlnqH6uoKAOpqGk/dlmCoryRJ6gRJpqr7Makv/OWHELaRmjXjbuAnIYQ5wO+Bv0xXX0VqmrrNpKaq+xuAGOOeEMJ3gf9M1/vO0S8PSpIkqWtUVGT2fCUlmT3fUbW1tcyaNYuKigrOPPNMysvLGTFiRKM6VVVVzJo1iz/84Q/06dOHW265hW9961sZb0uS2Taub2XX5BbqRuAbrZznX4F/PaHWSZIkqddbunQpgwYNYvPmzZSVlbFw4ULKy8sb1cnNzeX73/8+F110EdXV1ZSUlHD55ZdTXFyc0ba4wqCkbm9HK9+Q+CjB+JVDADsbl3VWz4gkqW2PP/449913HwcPHmTcuHE8+OCD5OTkHPeYZ555pn5RlRkzZnDrrbcSYyS1Pl/K0KFDGTp0KAADBgygqKiI7du3G54lSZLUM1VWVlJeXs66devIy8tj7ty5PPHEE6xcuZKNGzc2q79gwQJmzZrF9u3bGT48NfdEbm4uZ5xxBrt37yY/P7/F19m6dSuvv/4648aNy/g1GJ4lSZKUFWvWrKGiooJLLrkEgAMHDlBQUNBsCEZTqZHBjTXsdW5o//79TJ8+nSVLlnD66ad3vNFNGJ4lSZKUFTFGZs+ezeLFixuVl5aWHrfnubCwkKqqKgoLCzl8+DD79u1j8ODBzeofOnSI6dOnc+ONN3Ldddd1yjUYniVJkpQVkydPZurUqcyfP5+CggL27NlDdXV1mz3P11xzDcuWLWP8+PEsX76cyy67rFnPc4yROXPmUFRUxIIFCzrtGgzPkiRJvVS2v0BdXFzMnXfeyZQpUzhy5Ah5eXk88MADnH322cc9bs6cOdx0002MHDmSwYMHU1ZWBsCOHTu4+eabWbVqFevWreOxxx5j9OjRjB07FoB//Md/5Morr8zoNRieJUmSlDWlpaWUlpae0DH9+vXjqaeealY+bNgwVq1aBcCECRNaHBudae1dnluSJEnqdQzPkiRJUkKGZ0mSJCkhw7MkSZKUkOFZkiRJSsjwLEmSJCXkVHWSulztewc5cuoROBI4UhOpfe9gsgNHdm67JOlkV7GjIqPnKxnWORNH19bWMmvWLCoqKjjzzDMpLy9nxIgRzeqNGDGCAQMGkJOTQ25uLuvXr894WwzPkjps0KEtre7LPaUGgMCRNutKktSSpUuXMmjQIDZv3kxZWRkLFy5sdVXCtWvXkp+f32ltcdiGJEmSsubxxx/n0ksvZezYsXzta1+jrq6uzWOeeeYZZs+eDcCMGTNYs2ZNVhZEaYk9z5KAxj3CubGmWZkkSR1VWVlJeXk569atIy8vj7lz5/LEE0+wcuVKNm7c2Kz+ggULmDVrFtu3b2f48OEA5ObmcsYZZ7B79+5mPcwhBKZMmUIIga997WvccsstGb8Gw7MkSZKyYs2aNVRUVHDJJZcAcODAAQoKClodgnFUS73MIYRmZevWrWPYsGF88MEHXH755VxwwQX8+Z//eWYan2Z4lrKgop3fxzj6SdaW6sy1RZKkrhJjZPbs2SxevLhReWlp6XF7ngsLC6mqqqKwsJDDhw+zb98+Bg8e3Kz+sGHDACgoKGDatGn89re/NTxLkiSpZ5o8eTJTp05l/vz5FBQUsGfPHqqrq9vseb7mmmtYtmwZ48ePZ/ny5Vx22WXNep7/9Kc/ceTIEQYMGMCf/vQnfvnLX/L3f//3Gb8Gw7MkSVIv1VlTy7WmuLiYO++8kylTpnDkyBHy8vJ44IEHOPvss4973Jw5c7jpppsYOXIkgwcPpqysDIAdO3Zw8803s2rVKt5//32mTZsGwOHDh7nhhhu44oorMn4NhmdJkiRlTWlpKaWlpSd0TL9+/XjqqaealQ8bNoxVq1YB8OlPf5r/+q//ykgbj8ep6iRJkqSEDM+SJElSQoZnSZIkKSHDsyRJkpSQXxiUlFG17x1stN1n+BEA4pGW90uS1JPY8yxJkiQlZM+zJElSL1Vd3c4lcFsxYEDnzBtdW1vLrFmzqKio4Mwzz6S8vJwRI0Y0qrNx48ZGU+C9++67fOc732HevHkZbYvhWZIkSd3a0qVLGTRoEJs3b6asrIyFCxc2W5XwM5/5DBs2bACgrq6Os846q37RlExy2IYkSZKy5vHHH+fSSy9l7NixfO1rX6Ourq7NY5555hlmz54NwIwZM1izZg0xxlbrr1mzhnPPPbfNlQvbw55nSZIkZUVlZSXl5eWsW7eOvLw85s6dyxNPPMHKlSvZuHFjs/oLFixg1qxZbN++neHDhwOQm5vLGWecwe7du8nPz2/xdcrKyrj++us75RoMz5IkScqKNWvWUFFRwSWXXALAgQMHKCgoaDYEo6mWeplDCC3WPXjwIM8++yyLFy/ueINbYHiW1KIjNUe6/bRygw5tabPOYeC0mqalnfOFFknS8cUYmT17drNgW1paetye58LCQqqqqigsLOTw4cPs27ePwYMHt/gav/jFL7jooov45Cc/2SnXYHiWJElSVkyePJmpU6cyf/58CgoK2LNnD9XV1W32PF9zzTUsW7aM8ePHs3z5ci677LJWe55//OMfd9qQDTA8S5Ik9VqdNbVca4qLi7nzzjuZMmUKR44cIS8vjwceeKDNL/bNmTOHm266iZEjRzJ48GDKysoA2LFjBzfffDOrVq0C4OOPP+a5557joYce6rRrMDxLkiQpa0pLSxvNx5xEv379eOqpp5qVDxs2rD44A5x22mns3r27w208HqeqkyRJkhIyPEuSJEkJOWxDUq9TkdnVaClx8g5J6jXseZYkSZISMjxLkiRJCRmeJUmSpIQc8yxlSHV16wNp6+qy2BA1s6W6stH2nmYrDp64ooEOdJbU8+2o2JHR8w0rGdZmnbvuuosnn3ySnJwc+vTpw0MPPcTChQt59913ee+99+oXP7n22mt5/vnn2b9/PwCbNm1i3rx5bNq0iby8PEaPHs3999/faSsJtsbwLEmSpKx45ZVXWLFiBa+99hp9+/blww8/5ODBgwAMHDiQdevWMWHCBPbu3cvOnTvrj6upqeGqq67innvu4eqrrwZg7dq17Nq1y/AsnQwqdzXu6dxS3UUNkSSpG9m5cyf5+fn07dsXgPz8/Pp9M2fOpKysjAkTJvD0009z3XXX8dZbbwHw5JNPMn78+PrgDDBp0qTsNj7NMc+SJEnKiilTplBVVcX555/P3Llz+dWvflW/b/Lkybz00kvU1dVRVlbWaBXCN998k5JuMi+o4VmSJElZ0b9/fyoqKnj44YcZMmQIpaWlPPLIIwDk5OQwYcIEysvLOXDgACNGjOjStrbGYRvSSWbQoS1d3QRJklqVk5PDxIkTmThxIqNHj2bZsmX1+2bOnMm0adNYtGhRo2NGjRrVqJe6K9nzLEmSpKzYuHEj77zzTv32hg0bOPvss+u3v/jFL3L77bdz/fXXNzruhhtu4OWXX2blypX1ZatXr+aNN97o/EY3Yc+zpJNaLs174vM6cL5DFHXgaEnqXpJMLZdJ+/fv55vf/CZ79+4lNzeXkSNH8vDDDzNjxgwAQgjcdtttzY479dRTWbFiBfPmzWPevHnk5eUxZswY7r333qy2HwzPknqw2vcOtu/AkZlthyQpmZKSEl5++eVm5S+++GKL9Y/O8QxwwQUXsHr16s5qWmIO25AkSZISMjxLkiRJCRmeJUmSpIQMz5IkSVJChmdJkiQpIcOzJEmSlFCHpqoLIcwHbgYi8AbwN8BQoAwYDLwG3BRjPBhC6As8CpQAu4HSGOPWjry+JEmSOmBPRWbPN7ikzSp33XUXTz75JDk5OfTp04eHHnqIhQsX8u677/Lee+8RQgDg2muv5fnnn6+frm7Tpk3MmzePTZs2kZeXx+jRo7n//vv53e9+x7e//W0OHjzIKaecwve+9z0uu+yyzF5XA+0OzyGEs4C/BYpjjAdCCD8BZgJXAv8cYywLIfwLMAf4YfrxoxjjyBDCTOCfgNIOX4Gk42rPXMhHao50QkskSb3dK6+8wooVK3jttdfo27cvH374IQcPpv5ODRw4kHXr1jFhwgT27t3Lzp0764+rqanhqquu4p577uHqq68GYO3atezatYv8/Hx+/vOfM2zYMN58802+/OVXCGprAAAa6ElEQVQvs3379k67ho4O28gFTg0h5AKnATuBy4Dl6f3LgGvTz6emt0nvnxyO/tdCkiRJJ72dO3eSn59P3759AcjPz2fYsNQqhzNnzqSsrAyAp59+muuuu67+uCeffJLx48fXB2eASZMm8dnPfpYLL7yw/hyjRo2ipqaG2traTruGdofnGON24H8BvycVmvcBFcDeGOPhdLVtwFnp52cBVeljD6frn9ne15ckSVLPMmXKFKqqqjj//POZO3cuv/rVr+r3TZ48mZdeeom6ujrKysooLT02QOHNN9+kpKTtISE//elPufDCC+vDeWdod3gOIQwi1Zt8DjAM+ATwlRaqxqOHHGdfw/PeEkJYH0JYv2vXrvY2T5IkSd1M//79qaio4OGHH2bIkCGUlpbyyCOPAJCTk8OECRMoLy/nwIEDjBgx4oTO/dZbb7Fw4UIeeuihzDe8gY58YfBLwJYY4y6AEMLTwH8DBoYQctO9y4XAjnT9bcBwYFt6mMcZwJ6mJ40xPgw8DHDxxRc3C9eSJEnquXJycpg4cSITJ05k9OjRLFu2rH7fzJkzmTZtGosWLWp0zKhRoxr1Uje1bds2pk2bxqOPPsq5557bWU0HOjbm+ffA50MIp6XHLk8G3gbWAjPSdWYDz6SfP5veJr3/hRij4ViSJKmX2LhxI++880799oYNGzj77LPrt7/4xS9y++23c/311zc67oYbbuDll19m5cqV9WWrV6/mjTfeYO/evVx11VUsXryYL3zhC51+De3ueY4xvhpCWE5qOrrDwOukeoxXAmUhhDvTZUvThywFHgshbCbV4zyzIw2XJElSByWYWi6T9u/fzze/+U327t1Lbm4uI0eO5OGHH2bGjFS/awiB2267rdlxp556KitWrGDevHnMmzePvLw8xowZw7333ssPfvADNm/ezHe/+12++93vAvDLX/6SgoKCTrmGDs3zHGP8B+AfmhS/C1zaQt0a4C878nqSJEnquUpKSnj55Zeblb/44ost1j86xzPABRdcwOrVq5vVueOOO7jjjjsy1sa2uMKgJEmSlJDhWZIkSUrI8CxJktSL9Pb5Gjp6/YZnSZKkXqJfv37s3r271wboGCO7d++mX79+7T5Hh74wKEmSpJ6jsLCQbdu20ZsXouvXrx+FhYXtPt7wLEmS1Evk5eVxzjnndHUzejSHbUiSJEkJ2fMsqdcZdGhLu4/9IK8ogy2RJPU09jxLkiRJCRmeJUmSpIQMz5IkSVJChmdJkiQpIcOzJEmSlJDhWZIkSUrI8CxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEDM+SJElSQoZnSZIkKaHcrm6AJPV0FRWZPV9JSWbPJ0nKHHueJUmSpIQMz5IkSVJChmdJkiQpIcc8S93QoENburoJkiSpBfY8S5IkSQkZniVJkqSEDM+SJElSQo55lqQTkEclAHV1mTlfTo6TOktST2J4llrQnkUvGoapLdWZa4skSeo+HLYhSZIkJWTPs3qt6urWu5cz9ZG8JEk6udjzLEmSJCVkz7PUQ9S+d7CrmyBJUq9neJakdthSXZmR8xxKPxYNdNYNSeoJHLYhSZIkJWTPs6Rep71DYPqefQqDDm1p9+t+lHdOu4+VJHUP9jxLkiRJCRmeJUmSpIQMz5IkSVJChmdJkiQpIcOzJEmSlJDhWZIkSUrI8CxJkiQlZHiWJEmSEnKRFAmo3NV4qeUt1V3UEEmS1K3Z8yxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEDM+SJElSQoZnSZIkKSHDsyRJkpSQ4VmSJElKyEVSJCmh2vcOnvAxfc8+pRNaIknqKvY8S5IkSQkZniVJkqSEDM+SJElSQoZnSZIkKaEOhecQwsAQwvIQwv8OIVSGEMaHEAaHEJ4LIbyTfhyUrhtCCPeFEDaHEH4XQrgoM5cgSZIkZUdHe57vBVbHGC8APgdUAt8G1sQYzwPWpLcBvgKcl/65BfhhB19bkiRJyqp2h+cQwunAnwNLAWKMB2OMe4GpwLJ0tWXAtennU4FHY8pvgIEhhKHtbrkkSZKUZR3pef40sAv4txDC6yGEH4UQPgF8Msa4EyD9WJCufxZQ1eD4bemyRkIIt4QQ1ocQ1u/atasDzZMkSZIyqyPhORe4CPhhjPFC4E8cG6LRktBCWWxWEOPDMcaLY4wXDxkypAPNkyRJkjKrI+F5G7Atxvhqens5qTD9/tHhGOnHDxrUH97g+EJgRwdeX5IkScqqdofnGOMfgKoQwmfSRZOBt4FngdnpstnAM+nnzwKz0rNufB7Yd3R4hyRJktQT5Hbw+G8CT4QQTgHeBf6GVCD/SQhhDvB74C/TdVcBVwKbgY/TdSVJkqQeo0PhOca4Abi4hV2TW6gbgW905PUkSZKkruQKg5IkSVJChmdJkiQpIcOzJEmSlJDhWZIkSUrI8CxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEOrrCoKR2qNlUe9z9tX0OZqklkiTpRNjzLEmSJCVkeJYkSZISMjxLkiRJCTnmWZK6UB6VANTVHSurrm7/+QYMKOlgiyRJx2PPsyRJkpSQ4VmSJElKyPAsSZIkJeSYZ0nqRLXvHZuzu+ZI6/N77zttb/3z92vgkxcO7NR2SZLax55nSZIkKSHDsyRJkpSQ4VmSJElKyPAsSZIkJWR4liRJkhIyPEuSJEkJGZ4lSZKkhJznWeokgw5taXVfbZ+Dre6TJEndl+FZkrJkcJ+tre4bwLH/UPU9CH32Dmi0/8jAos5qliTpBDhsQ5IkSUrI8CxJkiQlZHiWJEmSEnLMsyR1M9u3w75TGpd93C/ZsTk5zctKSjreJklSiuFZkrqBXTU7Gm3/8ePG6fmjQ8nOc7Ra0UATsyR1BodtSJIkSQkZniVJkqSEDM+SJElSQoZnSZIkKSHDsyRJkpSQ4VmSJElKyPAsSZIkJWR4liRJkhIyPEuSJEkJGZ4lSZKkhFyeWz1GdXVFVzdBkiT1cvY8S5IkSQkZniVJkqSEDM+SJElSQo551kmhsrJjx2+pzkw7JEnSyc2eZ0mSJCkhe56lDqjZVNvqvto+B7PYEkmSlA32PEuSJEkJGZ4lSZKkhAzPkiRJUkKGZ0mSJCkhw7MkSZKUkOFZkiRJSsjwLEmSJCVkeJYkSZISMjxLkiRJCRmeJUmSpIQMz5IkSVJCuR09QQghB1gPbI8x/kUI4RygDBgMvAbcFGM8GELoCzwKlAC7gdIY49aOvr56p8pdlY22t1R3UUOkLBl0aEuieofTj6fVpB4/7lfSOQ2SpF6qw+EZ+BZQCZye3v4n4J9jjGUhhH8B5gA/TD9+FGMcGUKYma5XmoHXl6STTu17B9t1XM7ZGW6IJKmRDg3bCCEUAlcBP0pvB+AyYHm6yjLg2vTzqelt0vsnp+tLkiRJPUJHxzwvAf4OOJLePhPYG2M8+snhNuCs9POzgCqA9P596fqNhBBuCSGsDyGs37VrVwebJ0mSJGVOu8NzCOEvgA9ijBUNi1uoGhPsO1YQ48MxxotjjBcPGTKkvc2TJEmSMq4jY56/AFwTQrgS6EdqzPMSYGAIITfdu1wI7EjX3wYMB7aFEHKBM4A9HXh9SZIkKavaHZ5jjLcDtwOEECYCt8UYbwwhPAXMIDXjxmzgmfQhz6a3X0nvfyHG2KznWZLUcVuqUzPS7KnhWBdGB5QMc9YOSYLMzLbR1EKgLIRwJ/A6sDRdvhR4LISwmVSP88xOeG0po9qaHqy2T/tmRJAkST1TRsJzjPFF4MX083eBS1uoUwP8ZSZeT5IkSeoKndHzLEnqIrk0/rQkD6ir6dg5c/oVdewEknQScXluSZIkKSHDsyRJkpSQ4VmSJElKyPAsSZIkJWR4liRJkhIyPEuSJEkJGZ4lSZKkhAzPkiRJUkKGZ0mSJCkhw7MkSZKUkOFZkiRJSsjwLEmSJCVkeJYkSZISMjxLkiRJCRmeJUmSpIQMz5IkSVJChmdJkiQpodyuboAkqXNt2dKx4w8B7Dy2XVLSsfNJUk9mz7MkSZKUkOFZkiRJSsjwLEmSJCVkeJYkSZIS8guDknQSG3So/d8W/CjvnAy2RJJODvY8S5IkSQkZniVJkqSEHLYhSSeR2vcOnvAxfc8+pRNaIkknJ3ueJUmSpIQMz5IkSVJChmdJkiQpIcc8S0DNptoWy2v7nPj4UUmSdPIyPKtLVFSc+DF1dceeb6nOXFskSZKSctiGJEmSlJDhWZIkSUrI8CxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEnKpOktSmyr0N5pfc0fHzlQwr6fhJJKkLGJ7V6Sp2NJ/UuXLviZ8nj8oMtEaSJKn9HLYhSZIkJWR4liRJkhIyPEuSJEkJOeZZJ71Bh7a0Wae2z8EstESSJPV0hmd1murq1BcF62qaf9EvL9uNkdRuTb+sW1fTsfPl9Cvq2AkkqQs5bEOSJElKyJ5nSerlat9redhSzZHaFsv37U89njG6b2c1SZK6LXueJUmSpIQMz5IkSVJChmdJkiQpIcOzJEmSlJDhWZIkSUrI8CxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEDM+SJElSQi7PLUk6ITt2pB739G/f8YcAdh7bLinpaIskKXvseZYkSZISand4DiEMDyGsDSFUhhDeCiF8K10+OITwXAjhnfTjoHR5CCHcF0LYHEL4XQjhokxdhCRJkpQNHel5Pgz8nzHGIuDzwDdCCMXAt4E1McbzgDXpbYCvAOelf24BftiB15YkSZKyrt3hOca4M8b4Wvp5NVAJnAVMBZalqy0Drk0/nwo8GlN+AwwMIQxtd8slSZKkLMvImOcQwgjgQuBV4JMxxp2QCthAQbraWUBVg8O2pcuanuuWEML6EML6Xbt2ZaJ5kiRJUkZ0ODyHEPoDPwXmxRj/eLyqLZTFZgUxPhxjvDjGePGQIUM62jxJkiQpYzoUnkMIeaSC8xMxxqfTxe8fHY6RfvwgXb4NGN7g8EJgR0deX5IkScqmjsy2EYClQGWM8Z4Gu54FZqefzwaeaVA+Kz3rxueBfUeHd0iSJEk9QUcWSfkCcBPwRghhQ7rsfwB3Az8JIcwBfg/8ZXrfKuBKYDPwMfA3HXhtSZIkKevaHZ5jjP8vLY9jBpjcQv0IfKO9rydJkiR1NVcYlCRJkhLqyLANSdJJbHCfrcfd3/fQKa3u+yjvnAy3RpK6B8OzTio1m2qbldX2OdgFLZEkSScjh21IkiRJCRmeJUmSpIQMz5IkSVJChmdJkiQpIcOzJEmSlJCzbUiSsiqPSurqjm1XV3fsfAMGlHTsBJJ0Aux5liRJkhIyPEuSJEkJOWxDkpR1W6orj23s6ti5cqqhZJhDNyRlh+FZktQute+1vnpnzZHmq30e1e/8vp3RHEnKCsOzmqnYUZGR89TVVLZdSZIkqQcxPKvHGHRoS5t1avu03hMmSZLUUYZnJVLZjk7kvMw3Q5IkqUs524YkSZKUkD3PqlddnRrr3NJYZXuRJUmS7HmWJEmSErPnWZKUcYP7bG11X99DpzTa3rup8f6P8s45odc6BLDz2HaJUz5L6kT2PEuSJEkJGZ4lSZKkhAzPkiRJUkKOeVa3VbOp8fK+LoAiSZK6mj3PkiRJUkKGZ0mSJCkhh21IkrKq9r3jD8GqOVLbYnm/8/t2RnMk6YQYniVJPVoeldTVHduuru7Y+QYMcKJoSa1z2IYkSZKUkD3PkqQeb0t15bGNXR07V041lAyz91lSywzPyqpBh7YkruvUdJIkqbtx2IYkSZKUkOFZkiRJSshhG5KkbmVwn60tlvc9dEqbx36Udw5bko8Oa9EhgJ3Htksc/iypAXueJUmSpIQMz5IkSVJChmdJkiQpIcc8S5J6hLaW9YbmS3u7pLekTLPnWZIkSUrInuceqrq6oqubkFjNpmM9QS58IkmSejLD80mqsrLtOq3ZUp25dkiSJJ1MHLYhSZIkJWTPsyRJTVTuPTY0rnJNx89XNPDYSisuuiL1bIbnk0DlruZjNBx6IUntk0cHxr214BBFGT2fpK5leJYknTSaLu2dZElvSC3rLUlJGJ51wgYd2nJC9Z1hQ1JXSTI3NDg/tKTk/MKgJEmSlJDhWZIkSUrI8CxJkiQl5JhnSVKv194vGkLbXzbMo5K6umPb1RmYDWnAAOe7k7qKPc+SJElSQvY8S5LUybZUN5g7elfHz5dTDSXD7H2WuoLhWYnVbEpN5eTUc5IkqbcyPEuSlEVbTmyq/FZUUrct9awoAwsYOoZaSs7w3AUqdlR0+Bx1NZldPlaSdEzSxVUATmNju19nz5ERQPsWZakfCuIwECmrDM/dROUJZuG8zmmGJKmHyURP9iGAnannJWZo6bgMz1lSXX2st7mlXuNsh+ETXWIbHOssSZJkeO6lTuQjSUlS5zg6v/SJzCt9VFvzS0vqHIbnBByjLEmSJOiC8BxCuAK4F8gBfhRjvDvbbUjieMMsMvNN6cxoz/ALgNoMt0OSlF2nbW7fFxU/HvmZZmV5VLJ5b+pvXd27HWpWvaIhGZgGJM3ZQNSdZDU8hxBygAeAy4FtwH+GEJ6NMb6dzXacTBx+IUk9X3f6XZ6pDqItW1Jh/JxMjC7ZlflPby/99F8BUNHxD5cb8QuXJ79s9zxfCmyOMb4LEEIoA6YCvTo8t7f3GOxBliSdmKMLXp2ofuf3bdffq72bUo/dbYz2li2PZ/R8Iwf+VUbP99t3j7UvE/+hadi+TAX8o23M1H+4jraxu/8HJMQYs/diIcwArogx3pzevgkYF2O8tUGdW4Bb0pufgQ5MoNm6fODDTjhvb+S9zBzvZeZ4LzPHe5k53svM8V5mjvfymLNjjEPaqpTtnufQQlmj9B5jfBh4uFMbEcL6GOPFnfkavYX3MnO8l5njvcwc72XmeC8zx3uZOd7LE9cny6+3DRjeYLsQ2JHlNkiSJEntku3w/J/AeSGEc0IIpwAzgWez3AZJkiSpXbI6bCPGeDiEcCvwH6SmqvvXGONb2WxDWqcOC+llvJeZ473MHO9l5ngvM8d7mTney8zxXp6grH5hUJIkSerJsj1sQ5IkSeqxDM+SJElSQidNeA4hDA8hrA0hVIYQ3gohfKuFOjeGEH6X/nk5hPC5Bvu2hhDeCCFsCCGsz27ru5eE93JiCGFf+n5tCCH8fYN9V4QQNoYQNocQvp3d1ncvCe/l/9XgPr4ZQqgLIQxO7/N92UAIoV8I4bchhP9K38//u4U6fUMI5en336shhBEN9t2eLt8YQvhyNtvenSS8jwtCCG+nf1+uCSGc3WBfXYP3bK//0nfC+/nXIYRdDe7bzQ32zQ4hvJP+mZ3d1ncvCe/lPze4j5tCCHsb7PO92UQIISeE8HoIYUUL+/x92R4xxpPiBxgKXJR+PgDYBBQ3qfPfgEHp518BXm2wbyuQ39XX0R1+Et7LicCKFo7NAf4/4NPAKcB/NT22N/0kuZdN6l8NvNBg2/dl4/sTgP7p53nAq8Dnm9SZC/xL+vlMoDz9vDj9fuwLnJN+n+Z09TV14/s4CTgt/fzrR+9jent/V19Dd/pJeD//GvhBC8cOBt5NPw5KPx/U1dfUne9lk/rfJDX5wNFt35vN79EC4MlW/mb7+7IdPydNz3OMcWeM8bX082qgEjirSZ2XY4wfpTd/Q2qeaTWR5F4eR/0S7DHGg8DRJdh7pXbcy+uBH2ejbT1RTNmf3sxL/zT91vNUYFn6+XJgcgghpMvLYoy1McYtwGZS79deJ8l9jDGujTF+nN709+VxJHxftubLwHMxxj3pv0/PAVd0QjN7hHbcS39nHkcIoRC4CvhRK1X8fdkOJ014bij9scOFpP7H2po5wC8abEfglyGEipBaIly0eS/Hpz9a+0UIYVS67CygqkGdbSQP3ie1tt6XIYTTSP3R/GmDYt+XTaQ/gtwAfEAqdDS9n/XvwRjjYWAfcCa+NxtJcB8bavr7sl8IYX0I4TchhGs7taE9RML7OT09DGZ5COHogmG+L5tI+t5MDyU6B3ihQbHvzcaWAH8HHGllv78v2+GkC88hhP6kwse8GOMfW6kzidQfg4UNir8QY7yI1HCOb4QQ/rzTG9vNtXEvXyO1BvzngPuBfz96WAun6vXzISZ5X5IasrEuxrinQZnvyyZijHUxxrGkekIvDSF8tkmV1t6DvjcbSHAfAQgh/BVwMfC9BsWfiqnlfG8AloQQzu30BndzCe7nz4ERMcYxwPMc6+3zfdlE0vcmqWEGy2OMdQ3KfG+mhRD+AvggxlhxvGotlPn7sg0nVXgOIeSRCihPxBifbqXOGFIfX0yNMe4+Wh5j3JF+/AD4Gb3844m27mWM8Y9HP1qLMa4C8kII+bgEezNJ3pdpM2ny8aPvy9bFGPcCL9L8I+7692AIIRc4A9iD780WHec+EkL4EvA/gWtijLUNjjn6vnw3feyF2WhrT9Da/Ywx7m5wD/8foCT93PdlK4733kw73u9M35vwBeCaEMJWUkMoLwshPN6kjr8v2+GkCc/pMTpLgcoY4z2t1PkU8DRwU4xxU4PyT4QQBhx9DkwB3uz8VndPCe/ln6XrEUK4lNR7aTcuwd5IknuZrncG8N+BZxqU+b5sIoQwJIQwMP38VOBLwP9uUu1Z4OiMBTNIfQEzpstnpr9dfg5wHvDb7LS8e0lyH0MIFwIPkQrOHzQoHxRC6Jt+nk/qD/Tb2Wp7d5Twfg5tsHkNqe8/QGrF3Snp+zqI1L/z/+j8VndPCf+NE0L4DKkvWL7SoMz3ZgMxxttjjIUxxhGk/ha/EGP8qybV/H3ZDlldnruTfQG4CXgjPVYK4H8AnwKIMf4L8PekxvI8mM59h9Mf73wS+Fm6LBd4Msa4OrvN71aS3MsZwNdDCIeBA8DM9D+47rIEe3eR5F4CTAN+GWP8U4NjfV82NxRYFkLIIfUftp/EGFeEEL4DrI8xPkvqPyuPhRA2k+pBmQkQY3wrhPATUn9MDwPfaPJxb2+S5D5+D+gPPJV+D/4+xngNUAQ8FEI4kj727hhjrw0oaUnu59+GEK4h9d7bQ2r2DWKMe0II3yXV8QDwnSZDt3qbJPcSUl8ULEv/3TnK92YC/r7sOJfnliRJkhI6aYZtSJIkSZ3N8CxJkiQlZHiWJEmSEjI8S5IkSQkZniVJkqSEDM+SJElSQoZnSZIkKaH/HyukTecdNpgtAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 864x576 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(12,8))\n",
"plt.hist(result_re1.samples['sigma'], color='b', bins=30, alpha=0.2, label='e=0.2')\n",
"plt.hist(result_re2.samples['sigma'], color='g', bins=30, alpha=0.2, label='e=0.5')\n",
"plt.hist(result_re3.samples['sigma'], color='y', bins=30, alpha=0.2, label='e=0.7')\n",
"plt.hist(result_smc.samples['sigma'], color='purple', bins=30, alpha=0.2, label='SMC')\n",
"plt.hist(result_smc2.samples['sigma'], color='orange', bins=30, alpha=0.2, label='SMC2')\n",
"plt.axvline(np.mean(result_re1.samples['sigma']), color='b', alpha=1, label='e=0.2')\n",
"plt.axvline(np.mean(result_re2.samples['sigma']), color='g', alpha=1, label='e=0.2')\n",
"plt.axvline(np.mean(result_re3.samples['sigma']), color='y', alpha=1, label='e=0.2')\n",
"plt.axvline(np.mean(result_smc.samples['sigma']), color='purple', alpha=1, label='SMC')\n",
"plt.axvline(np.mean(result_smc2.samples['sigma']), color='orange', alpha=1, label='SMC2')\n",
"plt.axvline(np.sqrt(pos_mean_sigma), color='r', label='True')\n",
"plt.legend()\n",
"plt.show()"
]
}
],
"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.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment