infer_the_merger_fraction_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import binom\n",
    "from scipy.special import logsumexp\n",
    "import pymc3 as pm\n",
    "import theano.tensor as TT\n",
    "from astropy.coordinates import SkyCoord\n",
    "import astropy.units as u\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "from astropy.io import fits\n",
    "from multiprocessing import Pool\n",
    "#import accuracy_analysis as aa\n",
    "import mergers_w_humans_tools as m_tools\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The tools work with two main inputs - (1) an array of classifications (n x m) where m is the number of objects\n",
    "#being classified and n is the number of classifiers (2) an estimate of the classifiers accuracy, one can \n",
    "#measure this directly using simulated or well-known truth data sets \n",
    "\n",
    "#for the purpose of this notebook, we will use simulated observations, please replace the following arrays with\n",
    "#your own data \n",
    "\n",
    "# obs: measured classifications of UNKNOWN merger state sample, shape = n x m, where n is number of classifiers, m is\n",
    "# number of galaxies, classifications are binary: 1 (for merger) or 0 (for non-merger)\n",
    "\n",
    "# obs_train: measured classifications of KNOWN merger state sample or TRUTH sample, shape = n x m1, where n is number of classifiers,\n",
    "# m1 is number of KNOWN merger state galaxies, classifications are binary: 1 (for merger) or 0 (for non-merger)\n",
    "\n",
    "# true_samp: intrinsic classifications of known merger state sample or truth sample, len(true_samp) = n1  \n",
    "\n",
    "#making a synthetic sample\n",
    "#make_obs(number_of_classifiers,intrinsic_merger_fraction_of_unknown_sample,intrinsic_merger_fraction_of_known_sample\n",
    "#, #_of_galaxies_in_unknown_merger_state_sample,#_of_galaxies_in_known_merger_state_sample)\n",
    "classifications_unknown_sample, true, classifications_known_sample, true_train, truths_known_sample = m_tools.make_obs(10,0.4,0.5,150.,50.)\n",
    "\n",
    "#replace with your data here \n",
    "\n",
    "obs = classifications_unknown_sample\n",
    "obs_train = classifications_known_sample\n",
    "true_samp = truths_known_sample\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "CompoundStep\n",
      ">NUTS: [rM, rI, fM]\n",
      ">BinaryGibbsMetropolis: [true_class]\n",
      "Sampling 4 chains, 0 divergences: 100%|██████████| 6000/6000 [00:40<00:00, 146.94draws/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Median Merger Fraction: 0.39$\\pm$0.04\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEaCAYAAAAboUz3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcYElEQVR4nO3de5hcVZnv8e+PhEuAQIBEDgRCc1EQQQWjMooOgoxyU5lhBEYEEU6OckYZQQUcFFQY41ERxMOjKAiCAzqACKgcUK5yNQnXEPBwCRBASbgm4Rp454+1GiqVru7d3bWru1Z+n+epp6v2Ze137ap+a+21d62tiMDMzMqzwkgHYGZm9XCCNzMrlBO8mVmhnODNzArlBG9mVigneDOzQjnBW8dJWiRpkxHa9o8kfbVNZU3JdRmTX18l6eB2lJ3L+72kA9pV3iC2e5ykBZL+2ult10XSsZLOHuk4Os0Jvg0kzZX0kqSJTdNvlRSSekYotH5JOiPHvajhsXebt7FM0ouI1SPi/nZuJ29rrqTnJS2U9LSk6yV9RtJrn/OI+ExEfLNiWR/sb5mIeCjX5ZU2xL5MAoqIXSLizOGWPcg4NgQOB7aMiP/RhvK2z+/DM5KelHSdpHcOP1Krwgm+fR4A9u19IWlrYNxQC5M0th1BVSjv/+Qk1fv4ZdN6akyQXWCPiBgPbARMB44ATmv3Rtr9/owiGwFPRMTjwy1I0hrAJcDJwNrAZODrwIvDLduq6aZ/3NHuLGD/htcHAD9vXEDSypK+K+khSX/L3QXj8rwdJM2TdEQ+NP6ZpHGSzpT0lKQ5kr4saV5DeetLOl/SfEkPSPp8w7xjJZ0n6WxJzwKfqlqR3Oo+XtJ1wHPAJpIOzDEslHS/pP/VtM5H8xHLs5Luk/RhSccD7wN+mI8OfpiXDUmb5edrSvp5rsODko7u/UKR9ClJf8r77Klcx12q1CEinomIi4C9gQMkbZXLPEPScfn5REmX5Nb+k5KulbSCpLOAKcDFOe4vS+rJcR8k6SHgioZpjcl+U0k35xbrbySt3fj+Nu2zuZI+KOnDwFeAvfP2bmt4Hw7Oz1fI++ZBSY/nfbZmntcbxwH5s7VA0r/38/72uc+VjlguB9bPcZzRx7pr5X02P78nl0jaoMWm3pTfi3Mi4pWIeD4iLouI23NZm0q6QtITOeZfSJrQtH++JOl2SYslnSZpXaWuq4WS/iBpraZ9ME3So5Iek3R4P/tgO6Uji6cl3SZph1bLdrWI8GOYD2Au8EHgHuDNwBjgYVJrKICevNyJwEWk1sx44GLgW3neDsAS4NvAyqTW/3TgamAtYAPgdmBeXn4FYCbwNWAlYBPgfuBDef6xwMvAx/Ky4/qI+wzguD6mXwU8BLwFGAusCOwGbAoI+HtS4t82L/8u4Blg57ytycAWDWUd3FR+AJvl5z8HfpP3Rw/wF+CgPO9TuQ7/M+/TzwKPAurvfehj+kPAZ5vrDHwL+FGu34qkLyP1VVaOLXK8q+X3p3fa2Ia6PgJslZc5Hzi74f2d1yre/H6d3cf7cHB+/mng3vw+rw5cAJzVFNtPclxvI7WS39xiP/W3z5eJs2nddYB/AlbN6/8XcGGLZdcAngDOBHYB1mqav1n+zKwMTAKuAU5s2j83AuuSPlOPA7OAbfI6VwDHNO2Dc/K+3xqY39f+zWU9AexK+rzunF9PGulc0u6HW/Dt1duK3xm4m/TPDqSuDlKi+kJEPBkRC4H/APZpWP9V0gf2xYh4Hvg48B8R8VREzAN+0LDsO0kfyG9ExEuR+rR/0lTeDRFxYUS8msvryxdzK+ZpSQsapp8REbMjYklEvBwRv42I+yK5GriMlBABDgJOj4jL87YeiYi7B9pZSicn9waOioiFETEX+B7wyYbFHoyIn0Tq5z4TWI/0Dz8Yj5K+VJu9nMvbKNfx2sgZoB/HRsTifvbnWRFxZ0QsBr4KfDzXc7g+AZwQEfdHxCLgKGCfpqOHr0dqJd8G3EZK9EupuM9biognIuL8iHguf4aPJ33h97Xss8D2vP7lM1/SRZLWzfPvzZ+ZFyNiPnBCH2WdHBF/i4hHgGuBmyLiloh4Efg1Kdk3+np+f+4AfkZDt2mD/YDfRcTv8uf1cmAGKeEXxQm+vc4C/oXU8vx507xJpFbPzN6EClyap/eaHxEvNLxen3Qk0Kvx+UakQ+mnG8r7Cksnv8blW/luREzIj8aTxEutK2kXSTfmroynSf8MvctvCNxXYVvNJpKOPh5smPYgqYXV67UrOSLiufx09UFuZzLwZB/Tv0NqFV+m1O10ZIWyBtqnjfMfJB0ZTGyx7GCsz7L7aSxLv9+NV708R9/7qco+b0nSqpJ+nLt2niW1uie0+hKLiDkR8amI2IB0ZLM+6UgWSW+QdK6kR3JZZ7Psvvpbw/Pn+3jdXMfm/b9+H2FtBPxz0//O9qQv+6I4wbdRRDxIOtm6K+kQutEC0gfyLQ0Jdc2IaPyANrceHyN1zfTasOH5w8ADDWVNiIjxEdHYChnOUKGvrStpZVJ3w3eBdSNiAvA7UndNbyybDlROHxaQWtEbNUybQsORz3ApXbExGfjTMoGlFuzhEbEJsAdwmKSdeme3KHKgfdr4Hk0h1W8BsJj0Bd8b1xiW/nIfqNxHWXY/LWHphFfFcPf54cDmwLsjYg3g/Xm6Wq+S5KO6M0iJHlIXWQBvzWXtV6WcATTv/0f7WOZh0pFW4//OahExfZjbHnWc4NvvIGDHfIj+moh4lXSY+n1JbwCQNFnSh/op61fAUfnE1mTgXxvm3Qw8q3RSdpykMZK2Uj2XoK1E6vOcDyxROtH5Dw3zTwMOlLRTPlk3WdIWed7fSP3Gy8jdLr8Cjpc0XtJGwGGkltywSFpD0u7AuaS+1zv6WGZ3SZvl7rNngVfyo9+4B7CfpC0lrQp8Azgv1/MvwCqSdpO0InA0aZ/2+hvQo9ZXLJ0DfEHSxpJWJ3Xv/TIilgwmuDbs8/GkhsrTSieQj2m1oKQtJB3eexJW6RLMfUn96r1lLcplTQa+NJi6tPDVfJTxFuBA4Jd9LHM2sIekD+X/m1WUToK3OlnctZzg2yz3U89oMfsIUpfAjfmQ9A+k1lAr3wDmkY4K/gCcR77ELP+j7gG8Pc9fAPwUWHP4tVha7mv9PCkxPEXqhrqoYf7NpH+m75NOtl7N6y3Ek4C98hUXjecQen2O1Lq9n9TK/k/g9GGEe7GkhaRW2r+T+nUPbLHsG0n7dRFwA3BKRFyV530LODofwn9xENs/i9RK/SuwCmm/ERHPAIeQ3qNHSHVuvKrmv/LfJyTN6qPc03PZ15De7xdI+24ohrPPTySdyF1AStSX9rPsQuDdwE2SFufl7yQdBUC6ZHJb0mfmtyx71DsUV5P+x/5I6n68rHmBiHgY+CipS3M+6bPyJQrMh71XDFgXkPRZYJ+I6POkltnySunHhA8AKw72qKZkxX1jlUTSepLem7s9Nie1fH490nGZWXco9dd4pVgJ+DGwMfA0qT/5lJEMyMy6h7tozMwK5S4aM7NCjaoumokTJ0ZPT89Ih2Fm1jVmzpy5ICIm9TVvVCX4np4eZsxodYWhmZk1k/Rgq3nuojEzK5QTvJlZoZzgzcwK5QRvZlYoJ3gzs0I5wZuZFcoJ3sysUE7wZmaFcoI3MyvUqPolqy2/eo78bZ/T507frcORmJXDLXgzs0I5wZuZFcoJ3sysUE7wZmaFcoI3MyuUE7yZWaGc4M3MCuUEb2ZWKP/QybpK4w+i/CMos/65BW9mVigneDOzQjnBm5kVygnezKxQTvBmZoXyVTQ2qrUaRtjMBuYWvJlZoZzgzcwK5QRvZlYoJ3gzs0I5wZuZFcoJ3sysUE7wZmaFcoI3MyuUE7yZWaGc4M3MCuWhCmzEeBgCs3rV2oKX9AVJsyXdKekcSavUuT0zM3tdbQle0mTg88DUiNgKGAPsU9f2zMxsaXX3wY8FxkkaC6wKPFrz9szMLKstwUfEI8B3gYeAx4BnIuKyurZnZmZLq7OLZi3go8DGwPrAapL262O5aZJmSJoxf/78usIxM1vu1NlF80HggYiYHxEvAxcA72leKCJOjYipETF10qRJNYZjZrZ8qTPBPwRsJ2lVSQJ2AubUuD0zM2tQZx/8TcB5wCzgjrytU+vanpmZLa3WHzpFxDHAMXVuw8zM+uahCszMCuUEb2ZWKCd4M7NCOcGbmRXKo0la12ocjXLu9N1GMBKz0ckteDOzQjnBm5kVygnezKxQTvBmZoVygjczK5QTvJlZoXyZpHWUb7Rt1jluwZuZFcoJ3sysUE7wZmaFcoI3MyuUE7yZWaGc4M3MCuUEb2ZWKCd4M7NCOcGbmRXKCd7MrFBO8GZmhXKCNzMrlBO8mVmhnODNzArlBG9mVigneDOzQjnBm5kVygnezKxQTvBmZoVygjczK5QTvJlZoZzgzcwK5QRvZlaoWhO8pAmSzpN0t6Q5kv6uzu2ZmdnrxtZc/knApRGxl6SVgFVr3p4tp3qO/O1rz+dO320EIzEbPWpL8JLWAN4PfAogIl4CXqpre2ZmtrRKXTSSdpc02O6cTYD5wM8k3SLpp5JWG3SEZmY2JFVb8PsAJ0k6H/hZRMypWPa2wOci4iZJJwFHAl9tXEjSNGAawJQpUyoHbt2jsfvEzDqnUqs8IvYDtgHuI7XIb5A0TdL4flabB8yLiJvy6/NICb+57FMjYmpETJ00adIgwzczs1Yqd7tExLPA+cC5wHrAnsAsSZ9rsfxfgYclbZ4n7QTcNbxwzcysqkpdNJI+AhwIbAqcBbwrIh6XtCowBzi5xaqfA36Rr6C5P5dhZmYdULUPfi/g+xFxTePEiHhO0qdbrRQRtwJThx6emZkNVdUumseak7ukbwNExB/bHpWZmQ1b1QS/cx/TdmlnIGZm1l79dtFI+ixwCLCppNsbZo0HrqszMDMzG56B+uD/E/g98C3SNey9FkbEk7VFZWZmwzZQgo+ImCvpfzfPkLS2k7yZ2ehVpQW/OzATCEAN84I0HIGZmY1C/Sb4iNg9/924M+GYmVm7DHSSdZmhBRpFxKz2hmNmZu0yUBfN9/qZF8CObYzFzMzaaKAumg90KhAzM2uvgbpodoyIKyT9Y1/zI+KCesIyM7PhGqiL5u+BK4A9+pgXgBO8mdkoNVAXzTH5r0eBNDPrMlVv2beOpB9ImiVppqSTJK1Td3BmZjZ0VQcbO5d0f9V/Ig0dPB/4ZV1BmZnZ8FUdD37tiPhmw+vjJH2shnjMzKxNqrbgr5S0j6QV8uPjgO+kbGY2ig10meRCXh+D5jDg7DxrBWARcEyt0ZmZ2ZANdBXN+E4FYtYuPUe+fnA5d/puIxiJ2ciq2gePpLWANwKr9E5rvo2fmZmNHpUSvKSDgUOBDYBbge2AG/BYNGZmo1bVk6yHAu8EHszj02xDulTSzMxGqaoJ/oWIeAFA0soRcTeweX1hmZnZcFXtg58naQJwIXC5pKeAR+sKyszMhq9Sgo+IPfPTYyVdCawJXFpbVGZmNmyDuYpmW2B70nXx10XES7VFZWZmw1Z1sLGvAWcC6wATgZ9JOrrOwMzMbHiqtuD3BbZpONE6HZgFHFdXYGZmNjxVr6KZS8MPnICVgfvaHo2ZmbXNQGPRnEzqc38RmC3p8vx6Z+BP9YdnZmZDNVAXzYz8dybw64bpV9USjZmZtc1Ag42d2ftc0krAm/LLeyLi5ToDMzOz4ak6Fs0OpKto5pKGDt5Q0gEebMzMbPSqehXN94B/iIh7ACS9CTgHeEddgZmZ2fBUvYpmxd7kDhARfwFWrCckMzNrh6ot+JmSTgPOyq8/QTrxamZmo1TVFvxngNnA50lDB9+Vpw1I0hhJt0i6ZGghmpnZUAzYgpe0AjAzIrYCThjCNg4F5gBrDGFdMzMbogFb8BHxKnCbpCmDLVzSBsBuwE+HEJuZmQ1D1T749Ui/ZL0ZWNw7MSI+MsB6JwJfBlrevFvSNGAawJQpg/4OMTOzFqom+K8PtmBJuwOPR8TMfB19nyLiVOBUgKlTp8Zgt2OjU8+Rvx3pEMyWewONRbMK6WTqZsAdwGkRsaRi2e8FPiJpV9JAZWtIOjsi9htOwGZmVs1AffBnAlNJyX0X0g+eKomIoyJig4joAfYBrnByNzPrnIG6aLaMiK0B8nXwN9cfkln7NHYVzZ2+2whGYtZ5AyX41wYUi4glkoa0kYi4Co9AaWbWUQMl+LdJejY/FzAuvxYQEeFr283MRqmBhgse06lAzMysvaoOVWBmZl3GCd7MrFBO8GZmhXKCNzMrlBO8mVmhnODNzArlBG9mVigneDOzQjnBm5kVygnezKxQTvBmZoVygjczK5QTvJlZoZzgzcwKVfWm22YD8o22zUYXt+DNzArlFrwtN3x/VlveuAVvZlYoJ3gzs0I5wZuZFcoJ3sysUE7wZmaFcoI3MyuUE7yZWaGc4M3MCuUEb2ZWKCd4M7NCOcGbmRXKCd7MrFBO8GZmhfJokjYsHgPebPRyC97MrFBO8GZmhaotwUvaUNKVkuZImi3p0Lq2ZWZmy6qzD34JcHhEzJI0Hpgp6fKIuKvGbZqZWVZbCz4iHouIWfn5QmAOMLmu7ZmZ2dI60gcvqQfYBripj3nTJM2QNGP+/PmdCMfMbLlQe4KXtDpwPvBvEfFs8/yIODUipkbE1EmTJtUdjpnZcqPWBC9pRVJy/0VEXFDntszMbGl1XkUj4DRgTkScUNd2zMysb3W24N8LfBLYUdKt+bFrjdszM7MGtV0mGRF/AlRX+WZm1j+PRWPLpcYxdOZO320EIzGrj4cqMDMrlBO8mVmhnODNzArlBG9mVigneDOzQjnBm5kVygnezKxQvg7eBq20+7D6mngrlVvwZmaFcoI3MyuUE7yZWaGc4M3MCuUEb2ZWKCd4M7NCOcGbmRXKCd7MrFBO8GZmhfIvWa2S0n69arY8cAvezKxQbsGbNfC4NFYSt+DNzArlBG9mVigneDOzQjnBm5kVygnezKxQTvBmZoVygjczK5QTvJlZofxDJ7MW/KMn63ZO8NaSx58x625O8GYVuDVv3ch98GZmhXKCNzMrlBO8mVmhak3wkj4s6R5J90o6ss5tmZnZ0mo7ySppDPB/gZ2BecCfJV0UEXfVtU0bGl8tMzg+4Wrdos6raN4F3BsR9wNIOhf4KOAEPwo4qbfHYPejvxCsk+pM8JOBhxtezwPe3byQpGnAtPxykaR7hri9icCCIa7brVznLqNvD3qVrq7vELnOg7NRqxl1Jnj1MS2WmRBxKnDqsDcmzYiIqcMtp5u4zuVb3uoLrnM71XmSdR6wYcPrDYBHa9yemZk1qDPB/xl4o6SNJa0E7ANcVOP2zMysQW1dNBGxRNK/Av8PGAOcHhGz69oebejm6UKuc/mWt/qC69w2ilimW9zMzArgX7KamRXKCd7MrFBdleAHGvpAyQ/y/NslbTsScbZThTp/Itf1dknXS3rbSMTZTlWHuJD0TkmvSNqrk/HVoUqdJe0g6VZJsyVd3ekY263CZ3tNSRdLui3X+cCRiLNdJJ0u6XFJd7aY3/78FRFd8SCdqL0P2ARYCbgN2LJpmV2B35Ouwd8OuGmk4+5And8DrJWf77I81LlhuSuA3wF7jXTcHXifJ5B+BT4lv37DSMfdgTp/Bfh2fj4JeBJYaaRjH0ad3w9sC9zZYn7b81c3teBfG/ogIl4Ceoc+aPRR4OeR3AhMkLRepwNtowHrHBHXR8RT+eWNpN8bdLMq7zPA54Dzgcc7GVxNqtT5X4ALIuIhgIjo9npXqXMA4yUJWJ2U4Jd0Nsz2iYhrSHVope35q5sSfF9DH0wewjLdZLD1OYjUAuhmA9ZZ0mRgT+BHHYyrTlXe5zcBa0m6StJMSft3LLp6VKnzD4E3k34geQdwaES82pnwRkTb81c33bKvytAHlYZH6CKV6yPpA6QEv32tEdWvSp1PBI6IiFdS467rVanzWOAdwE7AOOAGSTdGxF/qDq4mVer8IeBWYEdgU+BySddGxLM1xzZS2p6/uinBVxn6oLThESrVR9JbgZ8Cu0TEEx2KrS5V6jwVODcn94nArpKWRMSFHYmw/ap+thdExGJgsaRrgLcB3Zrgq9T5QGB6pA7qeyU9AGwB3NyZEDuu7fmrm7poqgx9cBGwfz4bvR3wTEQ81ulA22jAOkuaAlwAfLKLW3ONBqxzRGwcET0R0QOcBxzSxckdqn22fwO8T9JYSauSRmad0+E426lKnR8iHbEgaV1gc+D+jkbZWW3PX13Tgo8WQx9I+kye/yPSFRW7AvcCz5FaAF2rYp2/BqwDnJJbtEuii0fiq1jnolSpc0TMkXQpcDvwKvDTiOjzcrtuUPF9/iZwhqQ7SN0XR0RE1w4jLOkcYAdgoqR5wDHAilBf/vJQBWZmheqmLhozMxsEJ3gzs0I5wZuZFcoJ3sysUE7wZmaFcoI3MyuUE7yZWaGc4K12kkLSWQ2vx0qaL+mSEY7rlTy+eu+jpw1lTpB0SNO064dbbi5nnKSrJY1pmPY9SXdJOrnFOitJukZS1/yo0drHb7p1wmJgK0njIuJ5YGfgkaor5+FiNZyRBFuU8XxEvL3N250AHAKc0jshIt4zuGhb+jRpyOBXcnybAO+NiC1brRARL0n6I7A38Is2xWFdwi1465TfA7vl5/sC5/TOkLSfpJtzK/rHksZI6pE0R9IpwCxgQ0lflXS3pMslnSPpi4MtY6AgW2z3wjxE72xJ0xqW3T/feee2hiOU6cCmOY7v5OUW5b+HSbozP/6taXs/yeVfJmlci/A+QRqTBkmbA1cDG0m6RdJq/VTrwryuLW9G+i4nfpT/ABYBbyUNDLYKaQjYHYBLSON9XwysmJc9Bdgf6CGNubJdnj41rzcOGA/8f+CLeV6lMvqI65Vc5q3Ar/O0ZdYB1s5/xwF3ksb+eQtwDzCxaZkemu7Yk+v/DtKY5quRbl4xG9gmL78EeHte9lfAfn3EuhLw16ZpxwEHV9j/Y4D5I/058KPzD7fgrSMi4nZSMtuXNKhSr51Iye/Pkm7NrzfJ8x6MdGcbSOPc/yYino+IhaSEPtgymj0fEW/Pjz0bpjev83lJt5HumLUh8EbSGOXnRR78KiL6u1NPb/y/jojFEbGINALo+/K8ByLi1vx8Jmk/NZsIPN00bWvSre6Q9JCkj+TnX5F0Ze9Ckbp0XpI0foAYrTDug7dOugj4Lqn1vk6eJuDMiDiqccF8wnNx46R+yq1aRlWvrSNpB+CDwN9FxHOSriIdhYjB3Yyhv/hfbHj+CulIodnzebuN3gLMlrQhcD2wtaTbSTfHuKVp2ZWBFwYRrxXALXjrpNOBb0TEHQ3T/gjsJekNAJLWlrRRH+v+CdhD0iqSVuf1/vzBlDEUawJP5eS+BelmyL3b/LikdXq3macvJHUhNbsG+JikVXN/+Z7AtVWDiHTf3TGSVsnbGw+8HBHPkY5eLgbWBb4IXEU6EiAvuw6pi+blqtuzMjjBW8dExLyIOKlp2l3A0cBlufV5ObDMjYYj4s+kI4DbSN0bM4BnBlPGEF0KjM3lfpPUTUNEzAaOB67O3Tcn5OlPANflE6nfaYh/FnAG6W5EN5HGc29uZQ/kMl6/JeNWpPMBkBL8TNL9O28g3b91ZsN6H2DpbjFbTng8eOsaklaPiEVKdzS6BpiWE+dyQdI2wGER8cmm6b8kndsYExEv976OfHmnpAuAoyLino4HbSPKffDWTU6VtCWpL/rM5Sm5A0TELZKulDQmnzjtnb53fvpq02uUbod3oZP78skteDOzQrkP3sysUE7wZmaFcoI3MyuUE7yZWaGc4M3MCuUEb2ZWKCd4M7NC/TcHBshLPYA57AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtQAAAEWCAYAAABG5QDSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0uUlEQVR4nO3deZxcZZn3/883SYcEaE0mhMUEAhFkFYQEEhxQEB0TREH9KSCMgDKAgjLuO/KMj7igzwwIGJiAoIZF2UQHxA2GtZE0IPsSg4GwSAyNiZBAmr5+f5xTUKlUd9WpqlNL9/f9etWr6+zXuc+pU1fdfZ/7KCIwMzMzM7PajGp1AGZmZmZmncwJtZmZmZlZHZxQm5mZmZnVwQm1mZmZmVkdnFCbmZmZmdXBCbWZmZmZWR2cUJtZbiSFpK1rXPYvkt4+yLS9JT1Ubl5JX5Y0v7aIM8f4XkmPS/qHpF2bsc0K8QxaZu2o9Dg2eN3XSDoij3WbmZVyQm1ma0mTslVpkvhXST+StGGr4yoWETdGxLaDTDslIo4GkLRlmtSPySmU7wEnRMSGEXFnvSuTNFPSryT1SXpO0v2SvilpYgNizUXJ+dIn6X8kbV7NskMdx3pFxNyIuCCPdReTdHJ6jn2yZPy/p+NPzjsGM2s9J9RmVs67I2JDYDdgd+CrpTPkmKR2kmnAfY1YkaQ3A9cDNwPbRcQEYA7QD+zSiG3kqHC+bAb8FfhBi+NptoeB0trwD6fjM2v0Z8ufVbP8OaE2s0FFxBPANcBO8EoTjuMlPQI8ko77N0mLJD0r6SpJrytZzf6SFkv6m6RTJY1Kl3u9pD9IWp5OWyBpQsmyu6e1tH1pTfm4dNl9JC0tF3NaY/jTdPCG9O9zaQ3qW9M431g0/8ZpDevkMusaJemrkpZIekbSjyW9VtJ6kv4BjAb+JOnPg8RyWtokZIWkXkl7D1bWwHeBH0XEtyLirwAR8VhEfD0irs9QZoVt7yHp1rSm+ylJZ0gam057c7r85unwLul820n6nKTLStb1A0n/NUTspPGuBi4Fdihadj1J35P0WPofj3mSxqfT1jqOaW33ZyXdLenvki4pHPN0+ufTfXlS0tEaokmRpOslFf5TcaSkm9I4+iQ9KmnuYPsh6YuS/ixpZXr+vbfCrt8OrC9px3T5HYHx6fji9R4g6a60rG+RtHPJvn9B0t3A85LGSPpweu4tl/Q1rd20aVRRnMsl/UzSP6XTCv+Z+aikx4A/VIjfzOrkhNrMBpUmXPsDxc0ZDgJmATtIehvwLeCDJLWTS4CLS1bzXmAmSW33gcBHCqtPl30dsD2wOXByybKHAe8EXg+8gTI15RW8Jf07IW2W8b9pfIcXzXMo8LuIWFZm+SPT177AdGBD4IyIeDGtkQXYJSJeP8j2bwfeBPwTcCHw8+IEsUDSBsCewGWl00pnpXKZFbwMfArYKF33fsDHASLiFuBs4II0uf0J8NWIeBD4KTCnkKgrqd08OJ1n6OCk9dN5e4pGf4fk2L0J2BqYApw0xGo+SFIzvxWwM0n5I2kO8Gng7el63lopnhKzgIdIyuO7wLmSNMi8fwb2Bl4L/B/gp5I2q7D+n5DUSkNSW/3j4omSdgPOA44FJpGU/1WS1iua7VDgXcAEkjI7i+QzsFkay5SieT9J8ll8K8n50AecWRLTW0nOk3dWiN3M6hURfvnll1+vvIC/AP8AniNJkM8CxqfTAnhb0bznAt8tGt4QWANsWTT/nKLpHwd+P8h2DwLuLInjuKLh/YE/p+/3AZaWzPv29P3JwE/T91umMYwpmncW8DgwKh1eCHxwkJh+D3y8aHjbdP/GFO3f1hnKto8kAS8dPzVd13ZF476bHoPnSZLdasvs7YPM++/AFUXDXUAvcA/wa0BF064B/i19fwBwf5XnSz/wJPDGdJrS+F9fNP+ewKNDHMfDS8pgXvr+POBbRdO2Hqr8SZrPHJ2+PxJYVDRt/XTZTas8bncBBw4y7WSSHyFbAI+l5foYyY+dnwInp/P9EPhGybIPAW8t2vePFE07CbioJOaXePU8fwDYr2j6ZoVzk1fP++nVnpt++eVXfS/XUJtZOQdFxISImBYRH4+IVUXTHi96/zqSpBuAiPgHsJy1a9KK51+SLlNoanGxpCckrSBJPjYqiaPssvWIiNtIkry3StqOJDG7apDZ19q/9P0YYJNqtiXpM5IeSJsvPEdSy1i6j5Ak2gMkSVEhzs9H0o76inSb1ZZZYdtvUHKD49PpvKcUzxsRa4DzSZrzfD8iomjxC3i1Fv9wKtdOH5TGuh5wAvC/kjYFJpMkgr1pM4fnSJL3dZrXFHm66P0LJD/SIDkWxedD8ftqvLLeiHghfVv2Ztu0qcVdRTHvxCDlXLTOx4BFJOX8SESUxjcN+Exhnel6N2ftc7r0s/XKcBrz8pL1XVG0rgdI/itRfG5mLSMzq5ETajPLqjjxepLkix14penCJOCJonmKe3zYIl0GkqYLAewcEa8hSdxK/wU/2LK1xFqskDD+K3BpJG1/y1lr/9IY+kluvBuSkvbSXyBpwjAxTTj/zrr7SEQ8D9wGvK/Caqsps4IfAg8C26Tzfrl4XklTgK8DPwK+X9L04EpgZ0k7kdRQL6gQV2E/Xo6Iy0kSu72AvwGrgB3TH2gTIuK18WpzmSyeIqnJL6iqJ5GsJE0D/pvkh8Gk9Ljdy+DlXOzHwGcoae6Rehz4ZlE5TIiI9SPioqJ5is/XtfY3bZozqWR9c0vWNy6S+x7Krc/McuSE2szqcSFwlKQ3pQnZKcBtEfGXonk+J2li2h77ROCSdHw3aVOBNLn7XJn1Hy9panqz1ZeLlq3WMpKa3+kl439C0rb7cMonPwUXAZ+StJWSrgNPAS6JiP4qtt1NknwvA8ZIOgl4zRDzfx74SHqj2cYAkqaStCUuXmelMiuedwXwj7Qm/mOFCWnb4fNJmux8lCR5+0Zherx6c+GFwB/T2teKlDgQmAg8EBEDJMnpfxbt0xRJtbTp/RnJubZ92lZ7qHbY9diAJBFdBiDpKNKbcqtwCfAvJLGW+m/gOEmz0nLaQNK7JHUPsq5LgXcruYF0LElb7uKkfh7wzfQHAJImp2VvZi3ghNrMahYRvwe+RnIz3VMkNw8eUjLbL0ja6t4F/A9JEgdJgrAbSa3t/wCXl9nEhcBvgMXp6/9mjO8F4JvAzem/xmen45cCd5AkTjcOsYrzSJLvG4BHgdXAJ6rc/LUkbZEfJmkqspoh/gUfETcBbyO5kfLhouYR1/NqN3TVlFnBZ4EPAStJkrniHyOfJGka8LW0qcdRJMlqcS8kFwBvpIqbEYFfKun1ZAVJeR8REYXuBL9A0hSiJ2168juStuiZRMQ1wOnAden6bk0nvZh1XRW2cz/w/XT9fyUpg5urXHZVRPyupIlUYdpC4N+AM0ia+CwiveFykHXdR3KuXUzy2VoJPMOr+3saSVOl30haSXIj6Kxq4jSzxtPazebMzEYGSecBT0ZE1p5DRgRJW5A0Gdk0Ila0Op5SkrYnaYqxXpX/Meho6X9IniNpwvNoi8MxsxKuoTazEUfSliTtlc+tMOuIpKSv8E8DF7dTMq3kUe9jlTw58jvAL4dzMi3p3ZLWT+9N+B5Jjyx/aW1UZlaOE2ozG1EkfYOkZvNU1/StK03eVgDvILlpsZ0cS9K2+c8kNz5+bOjZO96BJDfGPglsAxwS/reyWVtykw8zMzMzszq4htrMzMzMrA5jWh1AvTbaaKPYcsstWx2GmZmZmQ1zvb29f4uIdR5O1fEJ9ZZbbsnChQtbHYaZmZmZDXOSlpQb7yYfZmZmZmZ1cEJtZmZmZlYHJ9RmZmZmZnVwQm1mZmZmVgcn1GZmZmZmdWhaQi3pPEnPSLp3kOmSdLqkRZLulrRbs2IzMzMzM6tVM2uozwfmDDF9LsmjVbcBjgF+2ISYGqZ3SR9nXreI3iV9TVl/I7aXxzrrWX8t2897Hxoh73JtdTytiKGW4551m83+TLdDTHnI43Nfz/baZZ1Z1t+M873eGLMu34zvsKzD9WrHc68ZMbVjjK3QtH6oI+IGSVsOMcuBwI8jeRZ6j6QJkjaLiKeaE2Htepf0cdj8Hl7qH2DsmFEsOHo2M6ZNzG39Jx2wI//xq/vq2l4e66xn/bWUYd77UNhGz+LlzJ4+qey6qpleul/AkMtUiqeecy3vc7UZMdRy3LMeh2Z/pms531tx7LKqVO5ATeXQzOPWiHXWE3OjzvdGXs+b8Zmtd51Zh1tdRnmssxkx1VuOnXhdG0w7taGeAjxeNLw0HbcOScdIWihp4bJly5oS3FB6Fi/npf4BBgLW9A/Qs3j5OvPU84uudP3X3PtUxe21Yp31rL9n8XJWr8m2/bz3ofBB//5vHuKw+T3rHKveJX18YN4tnHpt+enlYrzsjqVDrjPrPmfdx3qXb4RG7EPxuXLNvU9VPHeyHoe8y6mW9bfDscuq9FiVlvvldyxdZ5+GujZW+kzmUUb1rrOWmIvLoJbrXKPLodHXnWo+s7Wss57hVpdRHutsRkz1lmMnXtcG004JtcqMi3IzRsQ5ETEzImZOnrzO0x+bbvb0SYwdM4rRgq4xo16peSkovaBeeNtjFROxodY/d6fNhtxeLTE3Yp31rH/29EmM68q2/bz3oVKSXzw82IWgNEZBXV8klc61vJdvhEbsQ/G5MnenzSqeO+WOw1AX8bzLqZb1t8Oxy6r0WJWWe8Ba+zRx/bFDXhsrfSbzKKN615k15tIymLj+2MzXuUaXQ6OvO9V8ZmtZZz3DrS6jPNbZjJjqLcdOvK4NRkkLiyZtLGny8auI2KnMtLOB6yPionT4IWCfSk0+Zs6cGe3w6PGh/qV35nWLOPXahwAYLXjz1htx4yN/e2X40/+yLcfvu3Wm9VdqalBLzI1YZz3rr2X7ee5DoQZ6IGBc1+D/il3TP0DXEP+qKo4JqGqZLPtcy3418jjXotH7UG59+592AytW93PaIbuuMw9UPg55l1MjzvdOUKnc4dUmID2Ll691rSy9Npb7TBYvn8d1rHQfam2aUO25Vq4MCuOz7GPe1/N6l2/Gd1jW4Xq127nXrJjaMcY8SeqNiJnrjG+jhPpdwAnA/sAs4PSI2KPSOtsloR5K6ZfASQfsyFevvGfQRM3aRzN+BFjjVfoxVJjHx6H5qmlPXG3yCVQ8zu0gy7lWzblrZq3T8oRa0kXAPsBGwF+BrwNdABExT5KAM0h6AnkBOCoiKmbKnZBQw9C1Zb5YmjVW6X+FqvkvkLWHLNfG4Xqc/f1g1r4GS6ib2cvHoRWmB3B8k8JpuqtPfMtawzOmTfSF0iwns6dPYpRgIDq/XV6eSn/ot4Ms18bhepz9/WDN0o7XgE7VTjclmpk1xIxpE9lu026mThzvf5kPondJHw8+vZKlfatq6nGmHfg4m9VuOFwD2okTajMblrrHdTFlwngnWYPoWbycgbTFXyd3V+XjbFab4XINaBdNa/JhZtZMlxy7Z6tDaGvDtbmEmVXH14DGypxQS9oAWB0RL+cQj5mZNUGhuUSnt5/0Dyez2gyXa0C7qJhQSxoFHAIcBuwOvAisJ2kZcDVwTkQ8kmuUZmbWcN3juuge1+UvUrMRyteAxqmmDfV1wOuBLwGbRsTmEbExsDfQA3xb0uE5xmhVOPjsWzn47FtbHYaZmZnZiFNNk4+3R8Sa0pER8SxwGXCZpK6GR2ZmZmZm1gEq1lCXJtOS/iDpe5IOlfSGcvNY861cvYYnnlvlbm/MzMzMmqyWbvN+B0xIlz08fQKitZD7kjQzMzNrncy9fETEKZKmAZ8FbomIkxoflmVRri9J32BgZpW4hwyzkc3XgMbJXEMt6QDgQ8AA8D5JoxselWVS6EsS3JekDU++6dbMzNpZLQ92ORO4Gfgl0Ov+qFvPfUmamZmZtU4tTT6mSZoKzCBpQ71NRBza+NAsC/claWZmZtYaVSfUkuYCJ5PckPgn4D/dftrMzMzMRrosbajPAj4NzAbOAU6V5JppMzMzMxvRsjT5+GtE3Jy+/52kW4HbAHebZ2ZmZmYjVpYa6r9I+r+SxqbDa4CVOcRkZrYWP7jIzMzaWZaEOoD3AY9LuglYBFwvaZtcIjMzww8uMjOz9ld1k49CTx6SxgE7Abukr/mSpkfE5vmEaGYjmR9cZGZm7a5iQi1JERGF4YhYDSxMX4V5anmEuZlZRYUHFw2EH1xkZmbtqZpE+DpJn5C0RfFISWMlvU3SBcCH8wnPrLH8xL3OU3hw0dSJ41lw9GzXTpuZWduppsnHHOAjwEWStgKeA8YBo4HfkPRHfVdeAZqZ+cFFZmbWziom1GkTj7OAsyR1ARsBqyLiuZxjMzMzMzNre5kePR4Ra4CncorFzMzMzKzj+GZCMzMzM7M6OKE2MzMzM6tDVQm1Eu5n2szMzMysRFVtqCMiJF0JzMg3HKvVJcfu2eoQrEEK3fr5mJqZmXWGLE0+eiTtnlskZmZmZmYdKEtCvS9JUv1nSXdLukfS3XkFZpaHlavX8MRzq+hd0tfqUMzMzGyYyNJt3tzcojBrgt4lfTz49EoGAg6b3+On7pmZmVlDZKmhfgzYGzgiIpYAAWySZWOS5kh6SNIiSV8sM/21kn4p6U+S7pN0VJb1dzI/Ejt/PYuXMxDJ+zX9A/QsXt7agMzMzGxYyJJQnwXsCRyaDq8Ezqx2YUmj0/nnAjsAh0raoWS244H7I2IXYB/g+5LGZojRbFCzp09ilJL3XWNGMXv6pNYGZGZmZsNCliYfsyJiN0l3AkREX8Zkdw9gUUQsBpB0MXAgcH/RPAF0SxKwIfAs0J9hG2aDmjFtIttt2s2K1f2cdsiubu5hZmZmDZGlhnpNWsscAJImAwMZlp8CPF40vDQdV+wMYHvgSeAe4MSIWGcbko6RtFDSwmXLlmUIwUa67nFdTJkw3sm0mZmZNUyWhPp04ApgE0nfBG4CTsmwvMqMi5LhdwJ3Aa8D3gScIek16ywUcU5EzIyImZMnT84QgpmZmZlZY1Xd5CMiFkjqBfZLRx0UEQ9k2NZSoPhpi1NJaqKLHQV8OyICWCTpUWA74I8ZtmNmZmZm1jRV11BL+jSwP7Be+por6aOS3lTlKm4HtpG0Vdr2+hDgqpJ5HiNN2CVtAmwLLK42RjMzMzOzZstyU+LM9PXLdPhdJEnycZJ+HhHfHWrhiOiXdAJwLTAaOC8i7pN0XDp9HvAN4HxJ95A0EflCRPwt0x6ZmZmZmTVRloR6ErBbRPwDQNLXgUuBtwC9wJAJNUBEXA1cXTJuXtH7J4F/yRCT2bCzcvUaVqzup3dJn2+eNDMz6wBZbkrcAnipaHgNMC0iVgEvNjQqsxGq8DTHpX2rOGx+jx+RbmZm1gGy1FBfCPRI+kU6/G7gIkkbsHZf0mZWo3JPc3QtNVxy7J6tDsHMzGxQVSXU6YNWzidprrEXSfvm4yJiYTrLYblEZzbCFJ7mOBB+mqOZmVmnqCqhjoiQdGVEzCBpL21mOfDTHM3MzDpPljbUPZJ2zy2SEW7l6jU88dwqt5k1P83RzMysw2RJqPclSar/LOluSfdIujuvwEYS34hmZmZm1rmy3JQ4N7coRjjfiGZmZmbWubLUUD8G7A0cERFLgAA2ySWqEaZwIxr4RjQzMzOzTpMloT4L2BM4NB1eCZzZ8IhGoMKNaFMnjmfB0bNdO21mZmbWQbI0+ZgVEbtJuhMgIvokjc0prhGne1wX3eO6nEybmZmZdZgsNdRrJI0maeqBpMnAQC5RmZmZmZl1iCwJ9enAFcDGkr4J3ASckktUZmZmZmYdouomHxGxQFIvsB/JkxIPiogHcovMzMzMzKwDZGlDTUQ8CDyYUyxmZmZmZh2n6oRa0kzgK8C0dDmRPJV855xiMzMzMzNre1lqqBcAnwPuwTcjmpmZmZkB2RLqZRFxVW6RmJmZmZl1oCwJ9dclzQd+D7xYGBkRlzc8KjMzMzOzDpEloT4K2A7o4tUmHwE4oTYzMzOzEStLQr1LRLwxt0jMzMzMzDpQlge79EjaIbdIzMzMzMw6UJYa6r2AIyQ9StKG2t3mWce55Ng9Wx2CmZmZDTNZEuo5uUVhZmZmZtahsjx6fEmegZiZmZmZdaIsbajNzMxsGDr47Fs5+OxbWx2GWcdyQm1mZmZmVoeaE2pJ/yzpzEYGY2ZmZmbWabLclIikNwEfAj4IPIof6mJmZmZmI1zFhFrSG4BDgEOB5cAlgCJi35xjMzMzMzNre9XUUD8I3Ai8OyIWAUj6VK5RmZmZmZl1iGraUL8feBq4TtJ/S9qP5KEumUmaI+khSYskfXGQefaRdJek+yT9by3bMTMzMzNrlooJdURcEREHA9sB1wOfAjaR9ENJ/1LthiSNBs4E5gI7AIeWPspc0gTgLOA9EbEj8IFq129mZmZm1gpV9/IREc9HxIKIOACYCtwFlK1lHsQewKKIWBwRLwEXAweWzPMh4PKIeCzd5jMZ1m8t5n5MzcxsuPJ3nA2lpm7zIuLZiDg7It6WYbEpwONFw0vTccXeAEyUdL2kXkkfLrciScdIWihp4bJly7IFb2ZmZmbWQM18sEu5dtdRMjwGmAG8C3gn8LW0l5G1F4o4JyJmRsTMyZMnNz5SMzMzM7MqZeqHuk5Lgc2LhqcCT5aZ528R8TzwvKQbgF2Ah5sTopmZmZlZNs2sob4d2EbSVpLGkvRtfVXJPL8A9pY0RtL6wCzggSbGaGZmZmaWSdU11JLWI+lCb8vi5SLiP6pZPiL6JZ0AXAuMBs6LiPskHZdOnxcRD0j6NXA3MADMj4h7q43RzMzMzKzZsjT5+AXwd6AXeLGWjUXE1cDVJePmlQyfCpxay/rNzMzMzJotS0I9NSLm5BaJmZmZmVkHypJQ3yLpjRFxT27RmBmXHLtnq0MwMzOzDLIk1HsBR0p6lKTJh4CIiJ1ziczMzMzMrANkSajn5haFmZmZtczK1WtYsbqf3iV9zJg2sdXhmHWcLI8eX1LulWdwZmZmlq/eJX08+PRKlvat4rD5PfQu6Wt1SGYdp2JCLWmlpBVlXislrWhGkGZmZpaPnsXLGUifW7ymf4CexctbG5BZB6rY5CMiupsRyEjnG9HMzKwVZk+fxCjBQEDXmFHMnj6p1SGZdZxmPnrczMzM2syMaRPZbtNuVqzu57RDdnUbarMaOKE2MzMb4brHddE9rsvJtFmNqr4p0czMzMzM1lV1DbWkT5cZ/XegNyLualhEZmZmZm3GXQvaULLUUM8EjgOmpK9jgH2A/5b0+caHZp1m5eo1PPHcKne5ZGZmw4q7FrRKsiTUk4DdIuIzEfEZkgR7MvAW4MgcYrMO4ouNmZkNV+5a0CrJklBvAbxUNLwGmBYRq0geRW4jmC82ZmY2XBW6FgR3LWjlZenl40KgR9IvAAHvBi6StAFwfx7BWedwP6ZmZjZcuWtBq6TqhDoiviHpamCvdNSxEbEwfX9YwyOzjuKLjZmZDWfuWtCGUnWTD0nrAdsCGwCvBfaXdFJegVnn6R7XxZQJ432xMTMzsxElS5OPX5B2k4fbTJuZmZmZAdkS6qkRMSe3SMzMzMzMOlCWXj5ukfTG3CIxMzMzM+tAWWqo9wKOlPQoSZMPARERO+cSmZmZmZlZB8iSUM/NLQozMzMzsw6Vpdu8JXkGYmZm2Rx89q0AXHLsni2OxMxsZKvYhlrSTenflZJWlP7NP0QzMzMzs/ZVsYY6IvZK/3bnH46ZmZmZWWfJ8mAXSTpc0tfS4c0l7ZFfaGZmZmZm7S9Lt3lnAXsCH0qH/wGc2fCIzMzMzMw6SJZePmZFxG6S7gSIiD5JY3OKy8zMzMysI2SpoV4jaTQQAJImAwO5RGVmZmZm1iGyJNSnA1cAG0v6JnATcEouUZmZmZmZdYgs/VAvkNQL7EfylMSDIuKB3CIzMzMzM+sAWWqoiYgHI+JM4GLgwawbkzRH0kOSFkn64hDz7S7pZUn/X9ZtmJmZmZk1UzUPdpkt6XpJl0vaVdK9wL3AXyXNqXZDafvrM0keYb4DcKikHQaZ7zvAtdWu28zMzMysVaqpoT6DpK30RcAfgKMjYlPgLcC3MmxrD2BRRCyOiJdIarkPLDPfJ4DLgGcyrNvMzMzMrCWqSajHRMRvIuLnwNMR0QNJ84+M25oCPF40vDQd9wpJU4D3AvOGWpGkYyQtlLRw2bJlGcMwMzMzM2ucahLq4q7xVpVMiwzbUplxpcv/F/CFiHh5qBVFxDkRMTMiZk6ePDlDCGZmZmZmjVVNLx+7SFpBkhCPT9+TDo/LsK2lwOZFw1OBJ0vmmQlcLAlgI2B/Sf0RcWWG7ZiZmZmZNU3FhDoiRjdoW7cD20jaCngCOIRXH2Ne2NZWhfeSzgd+5WTazKy8lavXsGJ1P71L+pgxbWKrwzEzG7EydZtXj4joB04g6b3jAeBnEXGfpOMkHdesOMzMhoPeJX08+PRKlvat4rD5PfQu6Wt1SGZmI1bVD3ZphIi4Gri6ZFzZGxAj4shmxGRm1ol6Fi9nIL0LZU3/AD2Ll7uW2sysRZpWQ21mZo0ze/okRqW3eneNGcXs6ZNaG5CZ2QhWzYNdfpL+PTH/cMzMrBozpk1ku027mTpxPAuOnu3aaTOzFqqmyccMSdOAj0j6MSXd30XEs7lEZmZmQ+oe10X3uC4n02ZmLVZNQj0P+DUwHehl7YQ60vFmZmZmZiNSxSYfEXF6RGwPnBcR0yNiq6KXk2kzMzMzG9Gq7uUjIj4maRdg73TUDRFxdz5hmZmZmZl1hqp7+ZD0SWABsHH6WiDpE3kFZmZmZmbWCbL0Q300MCsingeQ9B3gVuAHeQRmZmZmzXHJsXu2OgSzjpalH2oBLxcNv0xJjx9mZmZmZiNNlhrqHwG3SboiHT4IOLfhEZmZmZmZdZAsNyX+P0nXA3uR1EwfFRF35hWYmZmZmVknyFJDTUTcAdyRUyxmZmZmbcntzG0oWdpQm5mZmVmVDj77Vg4++9ZWh2FN4ITazMzMzKwOVTf5kHQCsCAi+nKMxzqY/x1mZmZmI1GWGupNgdsl/UzSHEnuMs/MzMzMRryqE+qI+CqwDUlXeUcCj0g6RdLrc4rNzMzMzKztZWpDHREBPJ2++oGJwKWSvptDbGZmZmZmbS9LG+pPAkcAfwPmA5+LiDWSRgGPAJ/PJ0QzMzMzs/aVpR/qjYD3RcSS4pERMSDpgMaGZdaZCt0j+QZNMzOzkSNLk4/1SpNpSd8BiIgHGhqVmZmZmVmHyJJQv6PMuLmNCsTMzMzMrBNVbPIh6WPAx4Hpku4umtQN3JxXYGZmZmZmnaCaNtQXAtcA3wK+WDR+ZUQ8m0tUZmZmZmYdomJCHRF/B/4OHJp/OGZmZmZmnaWaJh83RcReklYCUTyJpGvq1+QWnZmZmZlZm6umhnqv9G93/uGYmZmZDQ8rV69hxep+epf0MWPaxFaHYznK9KREMzMzM6usd0kfDz69kqV9qzhsfg+9S/paHZLlqGJCLWmlpBXp39LXimYEaWad7eCzb33loTdmZiNBz+LlDKQNZdf0D9CzeHlrA7JcVdPkw009zMzMzDKYPX0SowQDAV1jRjF7+qRWh2Q5qqaG+qb0b6Gmeq1Xlo1JmiPpIUmLJH2xzPTDJN2dvm6RtEuW9ZuZmZm1gxnTJrLdpt1MnTieBUfPdhvqYa5pNyVKGg2cSfLExaXA7ZKuioj7i2Z7FHhrRPRJmgucA8yqZ7tmzeQbUMzMrKB7XBfd47r8fTACNPOmxD2ARRGxOCJeAi4GDiyeISJuiYhCq/0eYGoT4zOri29AMTMzG5mqeVIiAJLGkTyCfC+S/qhvAn4YEaurXMUU4PGi4aUMXfv8UZInNJaL5RjgGIAtttiiys2b5avcDSiulbA8XXLsnq0OwczMyFZD/WNgR+AHwBnA9sBPMiyvMuOizDgk7UuSUH+h3PSIOCciZkbEzMmTJ2cIwSw/hRtQwDegmJmZjSRV11AD20ZE8U2C10n6U4bllwKbFw1PBZ4snUnSzsB8YG5EuI8Z6xiFG1BWrO7ntEN2de20mZnZCJGlhvpOSbMLA5JmATdnWP52YBtJW0kaCxwCXFU8g6QtgMuBf42IhzOs26wtdI/rYsqE8U6mzczMRpCKNdSS7iFpmtEFfFjSY+mkLYD7B12wRET0SzoBuBYYDZwXEfdJOi6dPg84CZgEnCUJoD8iZmbYHzMzMzPrUIWHgHXaPSLVNPk4oFEbi4irgatLxs0ren80cHSjtmdm7cHdCZqZ2XBWsclHRCwpvIAVwCbAtKKXmdmg3J2gmZkNd1W3oZZ0NHADSZON/5P+PTmfsMxsuCjXnaCZmdlwkuWmxBOB3YElEbEvsCuwLJeozGzYcHeCZmY23GXpNm91RKyWhKT1IuJBSdvmFpmZDQvuTtDMzIa7LAn1UkkTgCuB30rqo0w/0mZmpbrHddE9rsvJtJmZDUtVJ9QR8d707cmSrgNeC/w6l6jMzMzMzDpE1Qm1pHHAx4G9SPqlvolsbbDNzMzMzAbVqd2sZkmIfwzsCPwAOAPYHvhJHkGZmZmZ2cjSyd2sZmlDvW1E7FI0fJ2kPzU6IDMzMzNrD818cmG5blY7pZY6Sw31nZJmFwYkzQJubnxIZmZmZjbSdHI3qxVrqCXdQ9Jmugv4sKTH0klbAPfnGJuZmZlZx2pGre5w0sndrFbT5OOA3KMwMzMzsxGvU7tZrZhQR8SSwntJuwB7p4M3RoTbUHeQZraDMjMzMxspqm5DLelEYAGwcfr6qaRP5BWYmZmZmVknyNLLx0eBWRHxPICk7wC3knSjZ2ZmZmbDTKf2C91sWXr5EPBy0fDL6TjrECtXr+GJ51Z1VL+OZmZm1hqd3C90s2Wpof4RcJukK9Lhg4BzGx6R5aLwoRgIOGx+DwuOnu1fmjlw+3QzMxsuOrlf6GarqoZakoCfA0cBzwJ9wFER8V/5hWaNVO5DYWZmloeDz771lRvhrXN1cr/QzVZVDXVEhKQrI2IGcEfOMVkOCh+KgfCHwprPNfdmZp2nk/uFbrYsbah7JO2eWySWq8KHYurE8W7uYWZmufI9O8NH97gupkwY77yhgiwJ9b4kSfWfJd0t6R5Jd+cVmDWePxRmjeN/aZuV5xvZbCTKclPi3NyiMDMzs2HBN7LZSFQxoZY0DjgO2Bq4Bzg3IvrzDszMrJ25b1az8nzPjtWjU++5qabJxwXATJJkei7w/VwjMjNrc/6XttngfM+OjUTVNPnYISLeCCDpXOCP+YZkZtbe/C9ts6F1j+uie1yXPxc2YlSTUK8pvImI/qRLajOzkWuk/Eu7cNNlp/4L1szq589/dapJqHeRtCJ9L2B8OiySLqpfk1t01lD+UJg1xkjpm9XtxM3MqlOxDXVEjI6I16Sv7ogYU/R+xCbT+592A3t95w9uO2k2Qg33bijdTrxz9C7p48zrFvkYdSDnEsNHln6oLeUvGjMb7sq1E7f207ukjw/Mu4VTr33I30cdxrnE8OKEugb+omkM16o0h8vZalFoJw7Dq534cPs8tMv3UWm5XnLsnm5mWEZxObXLsbPGyPJgF0t1wg1JhQ/r7OmTav6XdCPWMdS6PzDvFgYCxnWNqrprpUbHlOc+toNay7mW7RSX43AvV+iMexLqOQ7t0k68kedSsz4PzdQO30fDsVzzUFpOJx2wY8uPnTVOUxNqSXOA04DRwPyI+HbJdKXT9wdeAI6MiDuaGWM12uWLZjC9S/o4bH4PL/UPMHZMbRe3chdIoGFfbNV2O1b8ZQrUvV+l685zH7PEkdc2qynnarY/1DzlviS+euU9dX25joSEvFSjfwQDdSc5V5/4lpriaJRGJ2rDsbvDVn0fFZ9reZXrcPuhXlpOfS+8VPexa8cyaceYmqFpCbWk0cCZwDuApcDtkq6KiPuLZpsLbJO+ZgE/TP+2nUpfNPWeUPUs37N4OavXDAC1X9xKP/iX3bGUy+9YulYyW5ivlhirqVUp/WHw/t2mZt6vocqxmn2sJfnMopofP5W2OdSXTqVyLrf9QtkUr2+oeUrL8Zp7n8r85ZrnD6dqZS3nStOzJAODnQdZzrdyn5dy/05u1T7WohGJWqXPQy37UG85NLqcqvnh08gYS8+1RtS0louvdBv/8av7cr0u5H0cy51/x++7daZ1VrpWQn2f8Vr2q3TZar5Tsqy/U35YNbOGeg9gUUQsBpB0MXAgUJxQHwj8OCIC6JE0QdJmEfFUE+OsW701xPUuP3v6JMZ1jWJN/0DNF7fSdQjWSmarST6HMmPaRH5+3JuH/FCU/jAIyLRflcqx0j6Wfnk3oua/0j5m3WY1XzpDlXPp9ssd10rznHTAjmuV49ydNuP2vzxb83Gq5YdTvbKWc9bjUCkZKHceQLYfFpU+LxPXH9vSfaxFvdeycvtU/HkAMidt5cqh+D8yldaRx3UkaznUG2Ppudb3wksVr+dZ4it33bnm3qdyvS7UUkZZ/3tS6Xsv62ew9FpZ6Xu52gqces7Par4vslwnWvHDqlbNTKinAI8XDS9l3drncvNMAdZKqCUdAxwDsMUWWzQ80Hr1LF7OS/0DDERtH/x6l58xbeIrF6Raf8GVrgPgsjuWrpV81hNjYRtDLVP6Zfr+3aby/t2mVr1flcqx0j6WfnnXe1yq2ces2yydfs29T60z//H7bj1onOV+VJQuX2mevhdeWud823bT7pqPU9YfTo2QtZyzHodyx2WoH3eFmv8s51ulz0ur97EW9V7Lyu1T8efhzOsWZd6HcuVQUM068riOZC2HemMsd75Wup5nia/cdSfrD/V6Y6imjAqyHMehyinrZ7D0Wlnpe7mac6/e87PS90XW60TW49JKzUyoVWZc1DAPEXEOcA7AzJkz15nearOnT2LsmNoTgnqXh8rJai3ryJJ8NsJgX6bV7lc15TjUPpZupxHHpdz269lm6fSsXzrV/Kiodp7SHyu1HqesP5waIWs5Zz0OlY7LYOdBlvOtms9LK/exVvVcy/L4/NRbDnlcRyrJup+VYmxEpU2l7ZXbRpYf6vXG0IrjmPV8Lb1WwtDfy9XEXO9+Vfq+yHqdqPc7rpmUtK5owoakPYGTI+Kd6fCXACLiW0XznA1cHxEXpcMPAfsM1eRj5syZsXDhwlxjr0Ur21A3y0iMsRX7XG+71kZvr9p5Gr3NvLWyDXWtMQ2HfcxbHp+fdmtDXY12j7Edzp12KKNGf0ZriTnv606nt6GW1BsRM9cZ38SEegzwMLAf8ARwO/ChiLivaJ53ASeQ9PIxCzg9IvYYar3tmlCbmZmZ2fAyWELdtCYfEdEv6QTgWpJu886LiPskHZdOnwdcTZJMLyLpNu+oZsVnZmZmZlaLpvZDHRFXkyTNxePmFb0P4PhmxmRmZmZmVg8/etzMzMzMrA5OqM3MzMzM6uCE2szMzMysDk6ozczMzMzq0LRu8/IiaRmwpEWb3wj4W4u2PVy4DBvD5dgYLsf6uQwbw+XYGC7HxnA5vmpaREwuHdnxCXUrSVpYri9Cq57LsDFcjo3hcqyfy7AxXI6N4XJsDJdjZW7yYWZmZmZWByfUZmZmZmZ1cEJdn3NaHcAw4DJsDJdjY7gc6+cybAyXY2O4HBvD5ViB21CbmZmZmdXBNdRmZmZmZnVwQm1mZmZmVgcn1DWQNEfSQ5IWSfpiq+PpFJI2l3SdpAck3SfpxHT8P0n6raRH0r8TWx1ru5M0WtKdkn6VDrsMM5I0QdKlkh5Mz8k9XY7ZSfpU+nm+V9JFksa5HCuTdJ6kZyTdWzRu0HKT9KX0O+chSe9sTdTtZZAyPDX9TN8t6QpJE4qmuQzLKFeORdM+KykkbVQ0zuVYhhPqjCSNBs4E5gI7AIdK2qG1UXWMfuAzEbE9MBs4Pi27LwK/j4htgN+nwza0E4EHioZdhtmdBvw6IrYDdiEpT5djBpKmAJ8EZkbETsBo4BBcjtU4H5hTMq5suaXXyUOAHdNlzkq/i0a681m3DH8L7BQROwMPA18Cl2EF57NuOSJpc+AdwGNF41yOg3BCnd0ewKKIWBwRLwEXAwe2OKaOEBFPRcQd6fuVJAnMFJLyuyCd7QLgoJYE2CEkTQXeBcwvGu0yzEDSa4C3AOcCRMRLEfEcLsdajAHGSxoDrA88icuxooi4AXi2ZPRg5XYgcHFEvBgRjwKLSL6LRrRyZRgRv4mI/nSwB5iavncZDmKQcxHgP4HPA8W9V7gcB+GEOrspwONFw0vTcZaBpC2BXYHbgE0i4ilIkm5g4xaG1gn+i+QiN1A0zmWYzXRgGfCjtOnMfEkb4HLMJCKeAL5HUoP1FPD3iPgNLsdaDVZu/t6pzUeAa9L3LsMMJL0HeCIi/lQyyeU4CCfU2anMOPc9mIGkDYHLgH+PiBWtjqeTSDoAeCYielsdS4cbA+wG/DAidgWex80SMkvb+B4IbAW8DthA0uGtjWpY8vdORpK+QtLMcEFhVJnZXIZlSFof+ApwUrnJZca5HHFCXYulwOZFw1NJ/sVpVZDURZJML4iIy9PRf5W0WTp9M+CZVsXXAf4ZeI+kv5A0N3qbpJ/iMsxqKbA0Im5Lhy8lSbBdjtm8HXg0IpZFxBrgcuDNuBxrNVi5+XsnA0lHAAcAh8WrD9twGVbv9SQ/kv+UftdMBe6QtCkux0E5oc7udmAbSVtJGkvSOP+qFsfUESSJpM3qAxHx/4omXQUckb4/AvhFs2PrFBHxpYiYGhFbkpx7f4iIw3EZZhIRTwOPS9o2HbUfcD8ux6weA2ZLWj/9fO9Hcm+Ey7E2g5XbVcAhktaTtBWwDfDHFsTX9iTNAb4AvCciXiia5DKsUkTcExEbR8SW6XfNUmC39LrpchzEmFYH0Gkiol/SCcC1JHe0nxcR97U4rE7xz8C/AvdIuisd92Xg28DPJH2U5Av6A60Jr6O5DLP7BLAg/WG8GDiKpJLB5ViliLhN0qXAHST/Xr+T5BHFG+JyHJKki4B9gI0kLQW+ziCf44i4T9LPSH709QPHR8TLLQm8jQxShl8C1gN+m/zGoycijnMZDq5cOUbEueXmdTkOzo8eNzMzMzOrg5t8mJmZmZnVwQm1mZmZmVkdnFCbmZmZmdXBCbWZmZmZWR2cUJuZmZmZ1cEJtZlZC0jaRNKFkhZL6pV0q6T3VlhmS0n35hjTTEmnZ1zmZEmfzSsmM7NO4H6ozcyaLH0IypXABRHxoXTcNOA9rYwrIhYCC1sZg5lZJ3INtZlZ870NeCki5hVGRMSSiPgBvFITfaOkO9LXm0tXMNg8kt4r6XdKbCbpYUmbpvO+qWj5myXtXLLOfST9Kn1/sqTzJF2f1qJ/smi+r0h6SNLvgG2Lxr9e0q/TGvcbJW0naYyk2yXtk87zLUnfbEgpmpm1CddQm5k1344kTxcczDPAOyJitaRtgIuAmdXMExFXSHo/cDwwh+SpZ09Lmg8cCfy7pDcA60XE3RXi3A7YF+gGHpL0Q2Bnksfe70ryHXIH0JvOfw5wXEQ8ImkWcFZEvE3SkcClaVI+B5hVYbtmZh3FCbWZWYtJOhPYi6TWenegCzgjrVF+GXhDmcWGmucTwL0kj12+KB33c+Brkj4HfAQ4v4rQ/iciXgRelPQMsAmwN3BFRLyQxn5V+ndD4M3Az9NHPkPyCOjC44p/AvwS2DMiXqpi22ZmHcMJtZlZ890HvL8wEBHHS9qIV9svfwr4K7ALSdO81WXWMdQ8U4ABYBNJoyJiICJekPRb4EDgg6xb413Oi0XvX+bV74woM+8o4LmIeNMg63oj8BxJUm5mNqy4DbWZWfP9ARgn6WNF49Yvev9a4KmIGAD+FRhdZh1l55E0BvgR8CHgAeDTRcvMB04Hbo+IZ2uM/QbgvZLGS+oG3g0QESuARyV9II1DknZJ378PmAS8BThd0oQat21m1pacUJuZNVlEBHAQ8FZJj0r6I3AB8IV0lrOAIyT1kDTleL7Magab58vAjRFxI0kyfbSk7dPt9gIrSBLuWmO/A7gEuAu4DLixaPJhwEcl/YmkFv7AtOb928BHI+Jh4AzgtFq3b2bWjpRc183MbLiT9DrgemC7tGbbzMwawDXUZmYjgKQPA7cBX3EybWbWWK6hNjMzMzOrg2uozczMzMzq4ITazMzMzKwOTqjNzMzMzOrghNrMzMzMrA5OqM3MzMzM6vD/A8qyInWX8GkUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#5000 draws takes < 5 minutes on a 32GiB, 8 core machine, ubuntu18 OS and python 3.6\n",
    "f_Ms_obs, r_Is_obs, r_Ms_obs, classif_obs, p_Ms_obs = m_tools.analyze_sample(obs, obs_train, true_samp, with_plots='y',draws=5000)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}