pyerrors/examples/03_fit_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('..')\n",
    "import pyerrors as pe\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read data from the pcac example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_obs = {}\n",
    "p_obs['f_P'] = pe.load_object('./data/B1k2_f_P.p')\n",
    "\n",
    "# f_A can be accesed via p_obs['f_A']\n",
    "\n",
    "[o.gamma_method() for o in p_obs['f_P']];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now define a custom fit function, in this case a single exponential. __Here we need to use the autograd wrapped version of numpy__  (imported as anp) to use automatic differentiation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import autograd.numpy as anp\n",
    "def func_exp(a, x):\n",
    "    y = a[1] * anp.exp(-a[0] * x)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit single exponential to f_P. The kwarg `resplot` generates a figure which visualizes the fit with residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit with 2 parameters\n",
      "Method: Levenberg-Marquardt\n",
      "`xtol` termination condition is satisfied.\n",
      "chisquare/d.o.f.: 0.00287692704517733\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAEsCAYAAAA8UOGyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyyUlEQVR4nO3dd3hUZd7/8fd3hoQkBAi9JRi69ABR7CKI/QHrWlGUlbUrrrtr2cdVHn8ra2Mta8GyNLE3RFERRLChAaI0KQKa0KWXQEJy//44Q00hkEzOTPJ5Xde5ppwzM5+cBL5zn3Of+zbnHCIiIhJdAn4HEBERkcOnAi4iIhKFVMBFRESikAq4iIhIFFIBFxERiUIq4CIiIlGomt8BDkf9+vVdamqq3zFEREQqxMyZM393zjUoal1UFfDU1FQyMjL8jiEiIlIhzOzX4taF5RC6mcWZ2fdm9qOZzTOzB4vYprqZvWFmS8xshpmlhiOLiIhIZRSuc+C7gN7Oua5AGnCWmR130DaDgI3OudbAcOBfYcoiIiJS6YSlgDvPttDDmNBy8Jit/YFRoftvA33MzMKRR0REpLIJ2zlwMwsCM4HWwH+cczMO2qQZkAXgnNttZpuBesDvB73PYGAwQPPmzcMVV0REIkxeXh7Z2dns3LnT7yhhFxcXR3JyMjExMaV+TdgKuHMuH0gzsyTgPTPr5JybewTvMwIYAZCenq6ZV0REqojs7Gxq1qxJamoqlfkArXOO9evXk52dTYsWLUr9urBfB+6c2wR8AZx10KoVQAqAmVUDagPrw51HRESiw86dO6lXr16lLt4AZka9evUO+0hDuHqhNwi1vDGzeKAv8PNBm40HrgndvxiY4g4xt+kZTbeWc1IREYlklb1473EkP2e4WuBNgC/M7CfgB2CSc26CmQ01s36hbV4G6pnZEuBO4O5DvenQtNWw4MMwRRYRETlQMBgkLS1t77J8+XJOOOEEAJYvX864ceN8yxaWc+DOuZ+AbkU8f/9+93cClxzO+87fHEeXt66Fy8ZB2zPKHlRERKQE8fHxZGZmHvDcN998A+wr4FdccYUPyaJsLPTbZjSFRh3hjavgly/8jiMiIlVQYmIiAHfffTfTp08nLS2N4cOHV3iOqCrg23YHYcB7UK81vHY5LP/a70giIlKJ5eTk7D18fsEFFxywbtiwYZx88slkZmYyZMiQCs8WVWOhA5BQF67+AEaeA+P+4BX0lGP9TiUiIuE08W5YPad837NxZzh7WImbFHUIPVJEVQt8r8QGcPV4SGwIYy+CFbP8TiQiIlKhoq8FvketJl4R/+85MOYCGDjB+zYlIiKVzyFayn6oWbMmW7f6d3lzdLbA90hKgWvGQ0wCjD4f1h58qbmIiEh4dOnShWAwSNeuXX3pxBa9LfA96raAaz70zomPOs+737C936lERKQS2LZtW7HPxcTEMGXKlIqOtFd0t8D3qN8aBn4EFoSR58Ga+X4nEhERCavKUcAB6rfxingwxmuJr5nndyIREZGwqTwFHPa1xIPVvZZ4eV9yICIiEiEqVwEHqNfK65EeEw+j+sGqn/xOJCIiUu4qXwGH/Yp4AozuB6t+9DuRiIhIuaqcBRygbkuviMcmei3xlZl+JxIRESk30X8ZWUnqtvCK+Mj/8VriV70Lyel+pxIRkXI2fNIinpy8uNDzt/dpw5C+bY/4fYPBIJ07dyYvL49q1apx9dVXM2TIEAKB4tu/y5cv55tvvgn7LGWVtwW+R51Ur4jH14HR/TUBiohIJTSkb1uWDzuXni3q0rNFXZYPO5flw84tU/GGfWOhz5s3j0mTJjFx4kQefPDBEl9TUfOEV/4CDlDnKLh2ItRqSt7oC7j63n+SevdHByzDJy3yO6WIiJRBfoFj445cVmzMYfKCNeQXuHJ9/4YNGzJixAieeeYZnHMsX76ck08+me7du9O9e/e984QfPM1ocduVlTlXvj9gOKWnp7uMjIwjf4Nt62DsBbh1C3mw+l18XnAMD/bvSK92DQkGrPyCiohImS1YsID27Us3smZ+gWPAyzP4bul6ChwkxAZJS0lizKCeZfr/PTExsdBobElJSSxcuJCaNWsSCASIi4tj8eLFXH755WRkZDB16lQee+wxJkyYAMCOHTuK3K40P6+ZzXTOFXnut3KfAz9YYgPyB3zI0n+fxd+3D2ND3o3c+lpuufySRUTEP1MXriUzaxN7Gt07cvPJzNrE1IVr6dO+UVg+My8vj1tuuYXMzEyCwSCLFhV9JLe02x2uqnEIfT9Tf8vl8l1380PB0fw75ln65U/a+0sWEZHoNG/lFnJy8w94Lic3n/krt5Tr5yxdupRgMEjDhg0ZPnw4jRo14scffyQjI4Pc3NwiX1Pa7Q5XlSvg81ZuYX1uLAPz/sqXBV0YFvMSl+VPKPdfsoiIVJyOTWsRHxs84Ln42CAdmtYqt89Yt24dN9xwA7fccgtmxubNm2nSpAmBQIAxY8aQn+99gTh4mtHitiurKlfA9/ySdxHLn/LuZGL+MdwfM4ZzNo6FKOoPICIi+/Rq15C0lCT2nAndcw68V7uGZXrfnJwc0tLS6NixI6effjpnnHEG//jHPwC46aabGDVqFF27duXnn3+mRo0aQOFpRovbrqyqVic2Cnd0qBkLz9Z4mZNzJkPPG+HMf0IJ1/eJiEjFOJxObOD9/372k9PYsSs/KjsoH24ntipXqYIBY8ygnrRumEhyUjz/vjydE+56yyveM56D9/4E+Xl+xxQRkcMUDBh1EmJpVieePu0bRVXxPhJVqxc6hUfrGTTKa9Hf3vtahvRpAJOHQs5G+MMoiC2fwxwiIhJeB//fnnr3R0DZR2KLZFWugA/p27aEX2Y7SKgHE4bA6PPhijcgoW5FxhMRkSNQ8v/tlVNYDqGbWYqZfWFm881snpndXsQ2vcxss5llhpb7w5HlsPUYCJeMglWZ8N+zYfMKvxOJiFRZ0dRPqyyO5OcM1znw3cCfnXMdgOOAm82sQxHbTXfOpYWWoWHKcvg69IOr3vGK9ytnwu+FB8gXEZHwiouLY/369ZW+iDvnWL9+PXFxcYf1urAcQnfOrQJWhe5vNbMFQDNgfjg+LyxanOJNgjL2Iq+IX/k2NOvudyoRkSojOTmZ7Oxs1q1b53eUsIuLiyM5OfmwXhP2c+Bmlgp0A2YUsfp4M/sRWAnc5ZybF+48h6VpGgz6DMacDyPP8zq2tenrdyoRkSohJiaGFi1a+B0jYoX1MjIzSwTeAe5wzh081Nks4CjnXFfgaeD9Yt5jsJllmFmGL9/C6rWCQZO823GXwsxRFZ9BRETkIGEr4GYWg1e8X3XOvXvweufcFufcttD9j4EYM6tfxHYjnHPpzrn0Bg0ahCtuyWo2hms/hlanwYe3wZSHNGqbiIj4Kly90A14GVjgnHuimG0ah7bDzI4NZVkfjjzlonpNuPx16DYApj0K790Au8tnQHoREZHDFa5z4CcCA4A5ZpYZeu5eoDmAc+554GLgRjPbDeQAl7lI72oYjIF+T0PSUfDFQ7B1FVw6BuJq+51MRESqmCo3Fnq5yXwNxt8C9dvClW9B7cPrPSgiInIoGgs9HNIu9y4t25QFL/WF1XP9TiQiIlWICnhZtDoNrvvEu//KWbB4kr95RESkylABL6vGneCPn0PdVBj3B/juefVQFxGRsFMBLw+1m8G1n0Dbs+GTv8FHd2pKUhERCSsV8PJSPREuHQsn3g4Zr8CrF0POJr9TiYhIJaUCXp4CAeg7FPo/C8u/hpf7wvpf/E4lIiKVkAp4OHS7Eq7+ALavg5f6wPKv/E4kIiKVjAp4uKSeCNdPgRoNYPT5MHus34lERKQSUQEPp7otvYlQUk+CD26GT+6B/N1+pxIRkUpABTzc4pO8AV963gjfPQtjL4AdG/xOJSIiUU4FvCIEq8HZw7zObb/NgBGnwuo5fqcSEZEopgJekbpdCddO9A6jv3wGzC00y6qIiEipqIBXtOQeMHgqNO4Mb18Lnz8IBfl+pxIRkSijAu6Hmo3gmgnQYyB89QS8dpkGfRERkcOiAu6XarHwP0/CecP
      "text/plain": [
       "<Figure size 576x355.995 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mass, Matrix element:\n",
      "[Obs[0.2102(63)], Obs[14.24(66)]]\n"
     ]
    }
   ],
   "source": [
    "# Specify fit range for single exponential fit\n",
    "start_se = 8\n",
    "stop_se = 19\n",
    "\n",
    "a = pe.fits.standard_fit(np.arange(start_se, stop_se), p_obs['f_P'][start_se:stop_se], func_exp, resplot=True)\n",
    "[o.gamma_method() for o in a]\n",
    "print('Mass, Matrix element:')\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The covariance of the two fit parameters can be computed in the following way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Covariance:  0.003465486601483565\n",
      "Normalized covariance:  0.8360758153764549\n"
     ]
    }
   ],
   "source": [
    "cov_01 = pe.fits.covariance(a[0], a[1])\n",
    "print('Covariance: ', cov_01)\n",
    "print('Normalized covariance: ', cov_01 / a[0].dvalue / a[1].dvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effective mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the effective mass for comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_eff_f_P = []\n",
    "for i in range(len(p_obs['f_P']) - 1):\n",
    "    m_eff_f_P.append(np.log(p_obs['f_P'][i] / p_obs['f_P'][i+1]))\n",
    "    m_eff_f_P[i].gamma_method()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the corresponding plateau and compare the two results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Effective mass:\n",
      "Obs[0.2114(52)]\n",
      "Fitted mass:\n",
      "Obs[0.2102(63)]\n"
     ]
    }
   ],
   "source": [
    "m_eff_plateau = np.mean(m_eff_f_P[start_se: stop_se]) # Plateau from 8 to 16\n",
    "m_eff_plateau.gamma_method()\n",
    "print('Effective mass:')\n",
    "m_eff_plateau.print(0)\n",
    "print('Fitted mass:')\n",
    "a[0].print(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAEKCAYAAAD6h5dQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsiUlEQVR4nO3de3xU5Z0G8OeXCYFEkgCSQG4akExu3BNAqpaUi0AlUAXLRd3axRKlWg2pldbqUls/UjVL1brdRFmwCqi1LpKCF4pGt1SQBAKYAIEimpBAIpCbCQmZefePmcEQcpkkM7wzJ8/38+FD5syZk2cSmGfOe86cV5RSICIiImPy0R2AiIiI3IdFT0REZGAseiIiIgNj0RMRERkYi56IiMjAWPREREQG5qs7QGuDBw9W0dHRumMQEXmV/Pz8r5VSIbpzkOfxuKKPjo5GXl6e7hiGUVZWhvDwcN0xiMjNRORL3RnIM3Ho3uCys7N1RyAiIo1Y9ERERAbGoiciIjIwjztGT66VkpKiOwIReZn8/PxQX1/flwGMBHcIPZ0VwOfNzc33JCUlVbS1Aove4Fj0RNRVvr6+Lw8dOjQ+JCTknI+PD2c+82BWq1UqKysTTp069TKAuW2tw3dqBpeZmak7AhF5n5EhISE1LHnP5+Pjo0JCQqphG31pe50rmIc0qK2t1R2BiLyPD0vee9h/V+32OYueiIg8jslkSoqLi0uIiYlJnD179vDa2lofAAgICBjX0eO+/vpr0+rVq91y4aB33323f0JCQryvr2/SunXrBrrje7gDi97gwsLCdEcgIuqyvn37Wg8fPlx09OjRwj59+qjMzEynyvvMmTOmtWvXhroj0/Dhw5vWrVt3IjU19Yw7tu8uLHqDS0tL0x2BiKhHbrzxxrpjx471bbmsurraZ/LkyeaEhIR4s9mc8Nprrw0AgIyMjMiSkpK+cXFxCWlpaZEA8Nhjjw0ZOXJkvNlsTkhPT794qdDp06dfl5iYGD9ixIjEZ599drBjectRg3Xr1g2cP39+NADExsY2TZo0qcHHx7uqk2fdG1xOTg5SU1N1xyAiL/WzTfuudcd2n188zqlL9l64cAHvv/9+0M0331zTcnlAQIB169atxwYNGmQtLy/3nTRpUtySJUuqMjMzS+fMmeN/+PDhIgB4++23g44dO9bvwIEDh5RSmD59+oh33323/+zZs+s2bNhwYsiQIZa6ujoZN25cwp133nlu6NChFnc8X51Y9AaXn5/Poicir9PY2OgTFxeXAACTJk2qffDBB79ueb/VapWHHnoocteuXf19fHxQUVHhV1paelmnvffee0GffPJJUEJCQgIA1NfX+xw+fLjf7Nmz637/+98P2bp16wAAOHXqVJ/CwsJ+Q4cO/eYKPL0ryqmiF5FZAJ4DYALwslJqdav7VwC4B0AzgEoA/66U+rLF/UEAigBsVkrd76LsRETkZs7uebua4xh9e/dnZWUNOnPmjO/BgwcP9e3bV0VERIxqaGi4bExdKYWHHnqo/OGHH77kjcLf/va3wI8//jgwLy/vcGBgoHXixImxjseLyMX1GhoaBF6u0wMNImIC8CKA2QASACwWkYRWq+0DkKyUGg3gLQBPt7r/twA+6XlcIiIioLq62jR48OALffv2VTk5OYFlZWV+ABAcHGz55ptvLnbb7Nmza1599dXB1dXVPgDwxRdf9Dl58qRvVVWVKTg42BIYGGjdt29fv/3791/leMzVV199Ye/evf0sFgveeecdrzm7vj3OnFEwEcAxpdRxpVQTgNcBzGu5glLqI6VUvf3mLgCRjvtEJAnAEAAfuCYydUVGRobuCERELnfPPfec3b9//1VmsznhlVdeuXrYsGHnAWDo0KGWpKSkupiYmMS0tLTI2267reb2228/O2HChDiz2Zxw6623XldVVWWaP39+dXNzswwfPjzx4YcfjhgzZszFIfvf/OY3J+fNmzdi/PjxcUOGDLngWP7xxx8HDBkyZPS2bdsGpqenXztixIhEHc+9q0Spjq+JICILAMxSSt1jv30XgEntDcGLyB8BnFJK/U5EfAB8COBOANNh2+u/7HEisgzAMgC45pprkr78ktMqu8qRI0cQGxurOwYRuZmI5Culkl2xrf37958YM2bM152vSZ5i//79g8eMGRPd1n0uPRlPRO4EkAxgin3RcgDblFKlLY95tKaUygaQDQDJycm8GpMLbdq0CatWrdIdg4i6YM32Yjy34+hlyx+cFoP0GWYNicibOVP0JwFEtbgdaV92CRGZDuBRAFOUUo32xZMB3CQiywH0B+AnInVKqZU9i01EZFzpM8xIn2HGwqxPAQBvpE3WnIi8mTNFvwdAjIgMg63gFwFY0nIFERkHIAu2If6L0+Qppe5osc7dsA3ds+SJiIiukE5PxlNKNQO4H8D7AA4BeFMpVSgiT4iIY0q8Z2DbY/+LiBSIyBa3JaYu4WfoiYh6N6eO0SultgHY1mrZ4y2+nu7ENtYDWN+1eNRTSUlJuiMQEZFG3nXBXuoynohHRNS7seiJiMjjeOI0tatWrRpy3XXXJZrN5oTJkyebi4uL/dzxfVyNRU9ERB7HE6epTUpKqi8oKDhUXFxc9IMf/OBcenp6ZOeP0o9Fb3BmMz9zS0TezVOmqU1NTa0NDAy0OjKVl5d7xR49Z68zuCVLlnS+EhFRe95a6pZparFgrVdPU5uVlRUyffr06q4/8SuPRW9wGzduZNkTkdfx5Glq/+u//mvQ/v37A7Kyso645Mm6GYve4IqLi3VHIKJusFgVztU3ob7Rgh2HTiMlNhQmHw0zpjq55+1qnjpN7ebNmwOfffbZsP/7v/874u/v7xWXbOcxeiIiD2OxKty1djeOVdShtKoBD2zah7vW7obF6hW9ckXomKZ2586d/g888MC177zzzrGIiIjmK/l8e4J79EREHib3SAUKSqrg6PX6JgsKSqqQe6QC0+KH6A3nIe65556zs2fPHmE2mxNGjx5d39Y0tVOnTq3OysoqLSws7DdhwoQ4wHZsf8OGDV/Mnz+/Ojs7O2T48OGJw4cPP9/WNLWDBg1qHjNmTL3jjcPDDz8cVV9fb7r99tuvA4Dw8PCmDz/88JiO598VnU5Te6UlJyervLw83TGIiLR5fsdRrNlejJavzgJgxQwzHpgW0+ZjOE1t79bRNLUcuje4/Px83RGIqIsSw4Pg72e6ZJm/nwkJ4UGaEpE3Y9EbXE5Oju4IRNRFKbGhGBs1AI5z7wL8TBgbNQApsW65DgwZHIueiMjDmHwEry6dhBGh/RE5wB8vLB6HV5dO0nPWPXk9noxHROSBTD6CgQF+GBgAnoBHPcI9eoNbvHix7ghERKQR9+gNLjw8vPOViIi66bd/Kwpf+48vwlovX3rjsPLH5iSU6chEl+IevcFlZmbqjkBEBvbYnISyE6tvyR8TGVw3JjK47sTqW/JPrL4lv6cl75im1vHnV7/61VAAeO+99/qPGDEiMS4uLqGurk7S0tIiR4wYkeiYwKYrVq5cObTl7XHjxsX1JLOn4h49ERH1SLPViuqGC6aGJotpS0FZ8PdHD6329enZfmR7l8D985//PGjFihXly5cvPwsAGzduHHzu3LkCX9+u19nzzz8ftnr16lOO2/v27Tvco9Aeinv0RETUbc1WK27/709jvjpb73+6ttHvF3/dP/z2//40ptlqdfn3+s///M/BW7duHfTkk09GzJ07d9jUqVNH1NfXm0aOHJnw0ksvDSwrK/OdOXPmdSNHjowfOXJk/AcffHAVYJvSdsGCBdFmsznBbDYnrF+/fsDy5csjHBPnzJ07dxjw7fS0c+bMGf76668HO77v/Pnzo9etWzewubkZaWlpkY4pb5955pnBrTMeOXLEb9iwYYnz58+Pjo6OHjl37txhmzdvDhw/fnzctddeO/Kjjz4KAICPPvooYOzYsXHx8fEJ48aNi9u/f39fAMjLy+s3atSo+Li4uASz2Zxw8ODBvjU1NT4pKSkjYmNjE2JiYhJfeumlga2/b0e4R29wSUlJuiMQkYFtO3Aq+FB5TX/H5XrPX7D6HCq
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pe.plot_corrs([m_eff_f_P], plateau=[a[0], m_eff_plateau], xrange=[3.5, 19.5], prange=[start_se, stop_se], label=['Effective mass'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting two exponentials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also fit the data with two exponentials where the second term describes the cutoff effects imposed by the boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func_2exp(a, x):\n",
    "    y = a[1] * anp.exp(-a[0] * x) + a[3] * anp.exp(-a[2] * x)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can trigger the computation of $\\chi^2/\\chi^2_\\text{exp}$ with the kwarg `expected_chisquare` which takes into account correlations in the data and non-linearities in the fit function and should give a more reliable measure for goodness of fit than $\\chi^2/\\text{d.o.f.}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit with 4 parameters\n",
      "Method: Levenberg-Marquardt\n",
      "`xtol` termination condition is satisfied.\n",
      "chisquare/d.o.f.: 0.05399877210985092\n",
      "chisquare/expected_chisquare: 0.7915235152326285\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAEsCAYAAAA8UOGyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1AElEQVR4nO3dd3xUVfrH8c+TSUJCAgQCoSMgRYpSjGJFpYldd1UUZe2sva2LZfe36q5r2XXFuqu4VhRlde2KiiCyioVQRIrSBOkJNZAEEibn98cdIIQkhCSTO5P5vl+vec3MbfOcTHlyzj33HHPOISIiItElzu8ARERE5MApgYuIiEQhJXAREZEopAQuIiIShZTARUREopASuIiISBSK9zuAA9G0aVPXvn17v8MQERGpFTNmzFjvnGtW1rqoSuDt27cnKyvL7zBERERqhZktL2+dmtBFRESikBK4iIhIFFICFxERiUJRdQ5cRERiR1FREStXrmT79u1+hxJ2SUlJtGnThoSEhErvowQuIiIRaeXKlTRo0ID27dtjZn6HEzbOOTZs2MDKlSvp0KFDpfdTE7qIiESk7du3k56eXqeTN4CZkZ6efsAtDVFVA7+x23qY8hAkp0GTgyH9YGjcHur4mysiEqvqevLepSrljKoEfn77zTDl/r0X1m8KBx0N7ftDt9OhYSuCxY4pP2Uzb3UuPVo15MSuGQTiYuNDICIiNScQCHDooYfufv7OO+8wfPhwpk2bxrJly5g2bRrDhw/3JbaoSuDHTehE1rdfQ/5G2LgEcn6CFd/B8q9gwfswYRSu7VE8v/VInt7Yl42FCSQnBujdNo2xV/RTEhcRkQOSnJzM7Nmz91o2bdo0AJYtW8a4ceN8S+C+nQM3s65mNrvELdfMbt7vjoEEaNAcDjoGMi+Dc/4FN8+B66bDSXexdt1qrtr8GJPtGv4QP5a0wnVMW7KBm8fPCn+hRESkzktNTQXgjjvu4H//+x+9e/dm9OjRtR6HbzVw59xPQG8AMwsAq4C3q3zAZl3ghFG8UXQ2X056nxHxn3JJ4FN+E/iU8cGT2Nno5poIW0REYkhBQQG9e/cGoEOHDrz99p409eCDD/Lwww/zwQcf+BJbpDShDwSWOOfKHfO1shblbOM7dwjfFR1CK9Zzbfy7DAt8jk3/AgK/hRNGQVKjGghZRERqzYQ7YO0PNXvMFofCKQ9WuElZTeiRIlIuI7sAeK2sFWY20syyzCwrJydnvwd6dFgfjjk4nfqJAdbQlPvjRvK7Fi8S6DUMvn4KnjgcZr4MxcU1XQYREZFa43sN3MwSgTOBO8ta75wbA4wByMzMdPs7XiDOGHtFP6b8lM381bl0D/VCj4s7E468EibcDu/dADNegrOehIxuNVoeEREJg/3UlP3QoEEDtm7d6tvrR0IN/BRgpnNuXU0dMBBnDOzWnBsGdmZgt+Z7ep+36gOXfwLnjIGNS+GZ/jD17xAsqqmXFhGRGHHYYYcRCATo1atXbHViK+FCymk+Dwsz6DUMDh4AE0bB5Ptg/rvwq2dVGxcRkb1s27at3GUJCQlMnjy5tkPazdcauJmlAIOBt2r9xVObwXkvwLBXYetaGHMiTH8O3H5b6UVERHznawJ3zuU559Kdc1t8C6Lb6XDNNDjoWPjwVhh/sTdQjIiISASLhHPg/kvNgIvehCH3wcJP4OnjYMV0v6MSEREplxL4LnFxcMwNcOVEb7S3F06BrBf8jkpERKRMSuClteoDV30OHfrDBzd7l5wV1f3J5EVEJLoogZelfhO46A04/jZv0JcXT4UtK/2OSkREZLdIuIwsMsUFYOD/Qave8PY1MOYkGP46tD7c78hERKSU0RMX8tikRfssv2lgZ24Z3KXKx901nWhRURHx8fH85je/4ZZbbiEurvz6b21NM6oEvj/dzoD0zjDuPHjhNPjVGOh+pt9RiYhICbcM7sItg7sw7JmvARj/26Nr5Lglx0LPzs5m+PDh5Obmcu+995a7T21NM6om9MrIOASunOwNfP+fESx+6y88/tlCJi1YR7BY142LiESCYLFjU34hqzYVhOX3OSMjgzFjxvDkk0/inGPZsmUcf/zx9O3bl759++6eJ7z0NKPlbVddqoFXVmozHm3zDzouH8WZcx6m+c4srt55OUXEc8OATvxuSFe/IxQRiVnBYseI575lcfY2ih3c8NoserdNY+wV/fYMp10DOnbsSDAYJDs7m4yMDCZOnEhSUhKLFi3iwgsvJCsra59pRvPz88vcrrqUwA/AoQc154ZpN/HzzpbcFP8WrS2HW+w2erdN8zs0EZGYNuWnbGav2MyuSnd+YZDZKzYz5adsBnZrHpbXLCoq4vrrr2f27NkEAgEWLlxYre0OlBL4AZi3OpeCwmJGcy7LizN4KOFZXnD38u2ypyFMHxAREdk/7/c5uNeygsIg81fn1mgCX7p0KYFAgIyMDO69916aN2/O999/T3FxMUlJSWXuM3r06Eptd6B0DvwALMreyq4zKm8V9+fKotvoaGs4a+ZlsGGJr7GJiMSyHq0akpwY2GtZcmKA7q0a1thr5OTkcPXVV3P99ddjZmzZsoWWLVsSFxfH2LFjCQa9fyBKTzNa3nbVpQR+AB4d1odjDk6nfmIAA6bH9+WBjL+RnlAIzw2BVTP9DlFEJCad2DWD3m3T2HW6u35igN5t0zixa0a1jltQUEDv3r3p0aMHgwYNYsiQIdx9990AXHvttbz00kv06tWLH3/8kZSUFGDfaUbL2666zEXR7FuZmZmuJk78V0ew2DHlp2zmr86le6uGnNg1g8DGJfDKOZC3AYaNhU4DfY1RRKQuWLBgAd26VX6a52Cx45THppK/I8i9Z/Xwfp9rsANbuJVVXjOb4ZzLLGt71cAPUCDOGNitOTcM7MzAbs29D0fTTnDFRGjSEcadDz+86XeYIiIxJxBnNK6fSOvGyXt+n+swdWKrKQ1awGUfwmvD4b9Xwo6tkHmZ31GJiMSE0iOxtb/jQ6D6I7FFMiXwmpTUCC5+E/7zG28ilB1b4dgb/Y5KRKTO2zUSWyxRE3pNS0iGYa9C97Nh4v/B5/dDFPUzEBGJJNHUT6s6qlJOXxO4maWZ2Ztm9qOZLTCzmhm81m/xiXDu89DnYvjiIfjkLiVxEZEDlJSUxIYNG+p8EnfOsWHDhgO+PtzvJvTHgI+dc+eaWSJQ3+d4ak5cAM54AhIbwDf/9JrTz3jMWy4iIvvVpk0bVq5cSU5Ojt+hhF1SUhJt2rQ5oH18S+Bm1gjoD1wK4JwrBAr9iics4uJg6AOQ1NCriRdug3PGeDV0ERGpUEJCAh06dPA7jIjlZw28A5ADvGBmvYAZwE3OuTwfY6p5ZnDSXZCY6p0TL8yD88dCQs0MpSciIrHJz3Pg8UBf4F/OuT5AHnBH6Y3MbKSZZZlZVlQ3oxx7I5z+KCyaCK8Ng8J8vyMSEZEo5mcCXwmsdM59G3r+Jl5C34tzboxzLtM5l9msWbNaDbDGZV4GZ/8Tln7hDfiyY5vfEYmISJTyLYE759YCK8xs10TaA4H5fsVTa3oPh189C8unwSu/hu25fkckIiJRyO/rwG8AXjWzOUBv4H5/w6klh50H5z4Hq7Jg7DlQsNnviEREJMr4ehmZc242UOYg7XVej3MgkAj/uQRePhNGvAP1m/gdlYiIRAm/a+Cx7ZDT4IJxkP0jvHQG5K33OyIREYkSSuB+6zIEhr8OGxbDi6cRzF3LpAXreHzSIiYtWEewuG6PQCQiIlXj90hsAnDwALjoDdy4Yax7fCB/KfwDywsbkRyakH7sFf3q/LR4IiJyYFQDjxQd+vNoiwdpWLSeF7mHlqwnvzDItCUbuHn8LL+jExGRCKMEHkEC7Y9hROGdNLGtjE/8C20sGwO6ZDTwOzQREYkwSuARZFH2Vma5zlxUeBcNLJ/xiX+hna1lYfZWv0MTEZEIowQeQR4d1odjDk5nSUJnhhf+gWQKeSf5Ph4dmOJ3aCIiEmGUwCNIIM4Ye0U/nriwD6cMGsLCoa+RlhxP4KXTYF3dH6ROREQqz6JpovTMzEyXlZXldxi1K2e
      "text/plain": [
       "<Figure size 576x355.995 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit result:\n",
      "[Obs[0.2146(65)], Obs[15.15(88)], Obs[0.623(60)], Obs[-9.64(74)]]\n"
     ]
    }
   ],
   "source": [
    "# Specify fit range for double exponential fit\n",
    "start_de = 2\n",
    "stop_de = 21\n",
    "\n",
    "a = pe.fits.standard_fit(np.arange(start_de, stop_de), p_obs['f_P'][start_de:stop_de], func_2exp, initial_guess=[0.21, 14.0, 0.6, -10], resplot=True, expected_chisquare=True)\n",
    "[o.gamma_method() for o in a]\n",
    "print('Fit result:')\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting with x-errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first generate pseudo data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(Obs[0.16(35)], Obs[0.15(25)])\n",
      "(Obs[2.21(35)], Obs[0.88(25)])\n",
      "(Obs[3.72(35)], Obs[-1.70(25)])\n",
      "(Obs[6.10(35)], Obs[-1.58(25)])\n",
      "(Obs[7.55(35)], Obs[-0.18(25)])\n"
     ]
    }
   ],
   "source": [
    "ox = []\n",
    "oy = []\n",
    "for i in range(0,10,2):\n",
    "    ox.append(pe.pseudo_Obs(i + 0.35 * np.random.normal(), 0.35, str(i)))\n",
    "    oy.append(pe.pseudo_Obs(np.sin(i) + 0.25 * np.random.normal() - 0.2 * i + 0.17, 0.25, str(i)))\n",
    "\n",
    "[o.gamma_method() for o in ox + oy]\n",
    "[print(o) for o in zip(ox, oy)];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And choose a function to fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(a, x):\n",
    "    y = a[0] + a[1] * x + a[2] * anp.sin(x)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then fit this function to the data and get the fit parameter as Obs with the function `odr_fit` which uses orthogonal distance regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit with 3 parameters\n",
      "Method: ODR\n",
      "Sum of squares convergence\n",
      "Residual variance: 0.03576834451052203\n",
      "Parameter 1 : Obs[0.02(40)]\n",
      "Parameter 2 : Obs[-0.225(75)]\n",
      "Parameter 3 : Obs[1.59(39)]\n"
     ]
    }
   ],
   "source": [
    "beta = pe.fits.odr_fit(ox, oy, func)\n",
    "\n",
    "pe.Obs.e_tag_global = 1 # Makes sure that the different samples with name length 1 are treated as ensembles and not as replica\n",
    "\n",
    "for i, item in enumerate(beta):\n",
    "    item.gamma_method()\n",
    "    print('Parameter', i + 1, ':', item)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the visulization we determine the value of the fit function in a range of x values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnA0lEQVR4nO3dd5hV1b3G8e+PGbrSHFSagggqIs2RMhZUUNEo2DBgiUSUeG1Yco2acnONiSaa2OI1ImPvsaLBgmLDUWBApEhHlKYM0nWAKev+sQ6KMI05ZZ3yfp7nPOecvTd7v3OY+c2etddey5xziIhI+qsTOoCIiCSGCr6ISIZQwRcRyRAq+CIiGUIFX0QkQ2SHDlCVnJwc1759+9AxRERSxrRp09Y451pWtC6pC3779u0pLCwMHUNEJGWY2ZeVrVOTjohIhlDBFxHJECr4IiIZQgVfRCRDqOCLiGQIFXwRkQyhgi8ikiGSuh++xJhzsHI6rJ4H3xVB6VaovyfkdIZOA0OnE5E4U8FPd87Bt4sh50Awg5cvg6J5P92m00k/FvzJD0DnQdB8/8RnFZG4UsFPV87B3HHw7l9g40q49nN/Nn9WPmQ3gCat/POWDVC2zf+b79fCm7+FN26Aw4bCgP+Bpm3Cfh0iEjNqw09H3y6GxwbDc78ADE65wxd3gH27+rP9eo2hThY0agF77uvXNWoBoz+DfpfDnJfhn7lQcC+Ul4f6SkQkhlTw083GlfDAMbDqM1/oL50E3X8OWXVr9u+btoETb4ErpkCH/vD+7bD5m/hmFpGEUJNOumnSGo7/PRxyKjRtW/v9NG8Pw5+GdUt9849z/peJmnhEUpbO8NPB1k2++WbFNP++76XRFfvtzKBFB//6k/vh//rC0knR71dEglDBT3VbNsJjp8Pc16BofvyOc8hpsGcreHIoLP0ofscRkbhRwU9lWzfBE2fBqhlwzqPQ49z4HatZOxjxmv/L4cmhsGxKjf/pnRMWxC+XiNSYCn6qKtkCT57jm3HOftifgcfbHnvDha/CnvvAM+f6Xzg1cPc7C+McTERqQhdtU5bzhffMMdBlcOIOu+e+cN7zsHqu79cvIilDBT8VlW6Fug39mb1Z4o+/V0f/AFg1E/Y51PfpF5GkFpMmHTN7yMxWm9nsStabmd1jZovMbKaZ9YrFcTPSzH/DA/1h8+owxX5H38yBMcfCB7eHzSEiNRKrNvxHgEFVrD8Z6BR5jALuj9FxM8vqufDqVdCwGTRsHjoN7N3FD8Hw3m2w5L3QaUSkGjEp+M65D4C1VWwyBHjMeZ8AzcysVSyOnTG2boJnL4B6e/imnJreORtPZnDqP/xomy9cDN+tCZ1IRKqQqF46bYBlO7xfHlm2CzMbZWaFZlZYVFSUkHAp4dWrYe1iOPshf+drsqjXGIY+DMXrYfyvQ6cRkSokXbdM59wY51yucy63ZcuWoeMkhy0bYc0COPYm6HB06DS72udQGPAH2KerBloTSWKJ6qWzAmi3w/u2kWVSEw2awMXvgCXd7+cfHXlVhYsLFq/54TmvY04iE4nIThJVQcYBv4j01ukLbHDOrUrQsVNXeTlMuhOK10F2PchKgV60c1+DCX8AfJEf+UghACMfKfyh+ItIGLHqlvk08DFwkJktN7ORZnapmV0a2WQ8sARYBDwIXBaL46a9wnx4+4++iKaKldPho7sp+Og9Rj5SSHFJGQDFJWUq+iKBmXMudIZK5ebmusLCwtAxwtiwHO7rA+36wPkvhO9zX1PbvqfgzuGMXD+CYldvl9UN62aRPyJXzTsicWJm05xzuRWtS4E2ggw1/nooL/PdHpOs2N85YUE14+OMqnRNcUkZ5z44udL1owd04poTOkeRTkQqo4KfjOa+CvP/Ayfc7CciSTLXnNC5yqJcsHgNI8d+RLHb9dtLZ/gi4SRxt48M1ron9LsC+qbmpY68jjnkn9OBhlk/bS5UsRcJSwU/GTVtCyf9OTnupq2lvJ7dyL+oLw3r+kHVVOxFwlPBTyar58ETZ8P6r0IniYm8jjnk9/a9b1XsRcJTwU8WzsEbN8DyKVC3Ueg0MZPXzA+xlFdnbuAkIqKCnywWvAFL3oVjb4TGaXQmvP06xMRb/C81EQlGBT8ZlG6FN2/yo04ecXHoNLFVt6F/XvaJ/4UmIsGo4CeDqfmwdgkMujWlL9RWqUlbePcvOssXCUj98JNBz/P9AGkHDgydJH5OvdNP3JJkN5GJZBIV/NCc88W+5/mhk8RX5xNDJxDJeGrSCWn9VzCmP6z8NHSSxNiyAV4dDQsnhE4ikpF0hh/Su7f6vveN03uil9EDOvkXdRvBoon+a+50QthQIhlIZ/ihfDMHPnsa+ozyd9amsR/G3cmqC3lX+B47X34cNpRIBlLBD+Wdm6F+Ezjq2tBJEqvnBdBoL/jortBJRDKOCn4Iy6b6G62OuhoatQidJrHqNYLev/Jf/zefh04jklFiNePVIDObb2aLzOyGCtaPMLMiM5sReaTZ3UW7qXVPOP1+6HNp9dumo96XQO5IqL9n6CQiGSXqi7ZmlgXcB5wALAemmtk459zOp2/POueuiPZ4aSErG3qcGzpFOI1a+IldRCShYnGG3xtY5Jxb4pzbBjwDDInBftOPc/Ds+fDpk6GTJIdlU+DTJ0KnEMkYsSj4bYBlO7xfHlm2s7PMbKaZPW9m7SrbmZmNMrNCMyssKiqKQbwksniin82qtDh0kuQwdSy8fgNs2Rg6iUhGSNRF21eB9s65bsAE4NHKNnTOjXHO5Trnclu2TKP+6c7Be7f6MWV6XhA6TXLocyls2wQz9BePSCLEouCvAHY8Y28bWfYD59y3zrmtkbdjgcNjcNzUsugdWD4VjrkOsuuHTpMc2vSCdn1g8gN+wnYRiatYFPypQCcz62Bm9YBhwLgdNzCzVju8HQxk1mwY28/um+4HPdJ8zJzd1edSWPcFLHwrdBKRtBd1Lx3nXKmZXQG8CWQBDznn5pjZzUChc24ccJWZDQZKgbXAiGiPm3KO/50f9z67XugkyeWQ02Dfw6B4XegkImnPXBKPT56bm+sKCwtDx6iVOycs+HFIAamacxo2WSRGzGyacy63onW60zZO7n5noX+x6B144ybYuilsoGRm5tvwv54VOolIWlPBjyfn4P2/wtxxkN0gdJrk9s7NMHagmnZE4kgFP56+LIBlkyHvqvSdujBWup4JpVvgs2dDJxFJWyr48fTh3/1Y973U775arbpDm1wofEjz3orEiQp+PC1+B/pdDnUbhk6SGnIvgjXz/V9GIhJzKvjx1PN8Pyqk1MyhZ0CDpjD7hdBJRNKSpjiMpyH3hU6QWuo1govegpxOoZOIpCWd4cfD1PzQCVLX3gdDnazQKUTSkgp+rK37koJXHwKgYPGawGFS1KdPQP6JUF4eOolIWlHBj7GC8Y8zctt1AIx8pFBFvzayG/jurEsmhk4iklZU8GOoYNZCRs7qQjF+NMzikjIV/do45DRolAOFD4dOIpJWVPBjpGDxGkY+PfeHYr+din4tZNf3U0DOfx02fRM6jUjaUC+d3XDnhAU/jpFToYovNhaXlHHug5Or3PfoAZ002NqOel4ABffAzGfgyNGh04ikBRX83XDNCZ0rLcoFi9cw8pFCikt2ncijYd0s8kfkktcxJ94R00fLztD/BtivX+gkImlDTTqxsHUzeW4G+Rfm0rDuT8/yVeyjcNyN0K536BQiaUMFPxamPwpPnEVeo+Xkj/ix6KvYx8DqeTD7xdApRNJCTAq+mQ0ys/lmtsjMbqhgfX0zezayfrKZtY/FcZNC6VYouBfaHw2te5DXMYf8EX7uARX7GCi4F8ZdCVs3h04iGeTOCQtCR4iLqAu+mWUB9wEnA12A4WbWZafNRgLrnHMHAncCf432uEnjs6dh0yo4+rofFm0v8ir2MdDrAti2GT5/JXQSySBVd85IXbE4w+8NLHLOLXHObQOeAYbstM0Q4NHI6+eBAWZpMKddWSlMugta94IDjg2dJj216wN7dYJPHw+dRCTlxaLgtwGW7fB+eWRZhds450qBDcBeFe3MzEaZWaG
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_t = np.arange(min(ox).value - 1, max(ox).value + 1, 0.01)\n",
    "y_t = func([o.value for o in beta], x_t)\n",
    "\n",
    "plt.errorbar([e.value for e in ox], [e.value for e in oy], xerr=[e.dvalue for e in ox], yerr=[e.dvalue for e in oy], marker='D', lw=1, ls='none', zorder=10)\n",
    "plt.plot(x_t, y_t, '--')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also take a look at how much the inidividual ensembles contribute to the uncetainty of the fit parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameter 0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADrCAYAAADKbEVrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfwUlEQVR4nO3deZxcVZ338c+p7k4n6YTs+0KF7CE3CRoIq2AIIjYq4rwURjCiPj5IXEYcpeYZ0Rp8fE0NboyjwqgjlIiij+KAFjKyxAACIQmdpLOTQCXpJCQhnb69L1X3PH9UQ9budFVX3XOX3/v16le6a/0mqfr2rXvPPUdprRFCCOGOiOkAQggRJlK6QgjhIildIYRwkZSuEEK4SEpXCCFcJKUrhBAuktIVoaeU+pJSarNSapNS6tdKqYGmM4ngktIVoaaUmgR8AVistZ4PlAE3mE0lgkxKVwgoBwYppcqBwcB+w3lEgEnpilDTWu8DvgPsAQ4Attb6L2ZTiSCT0hWhppQaAXwQmAZMBKqUUjeZTSWCTEpXhN0y4HWt9WGtdRfwCHCx4UwiwKR0RdjtAS5USg1WSingSmBrfx9UKTVcKfU7pdQ2pdRWpdRF/U4qAqHcdAAhTNJar1ZK/Q54BcgANcBPivDQ/w48obX+O6XUAHIH6IRAydSOQhSXUmoYsB44R8sbTJxEdi8IUXzTgMPA/UqpGqXUz5RSVaZDCW+Q0hWi+MqBdwD3aq3PA1qAmNlIwiukdIUovjqgTmu9uvvn35ErYSGkdIUoNq31G8BepdTs7ouuBLYYjCQ8RA6kCc+IxlJDgUnkTlKYdNz344BKoILcR/eKd7dVvLa4o3w6kCU36qADONL99Wb311vfHwbqVty3tMWtv4tSahHwM2AA8Bpwi9b6qFvPL7xLSle4KhpLDSb3UfsCYCEwmWMFO6Svj/Oe1opVCzvLL8/z6Q8Ar3Z/bQU2A5tW3Le0Ls/HEaJgUrqiZKKxVDkwHzifXMleAJxLbiavfimwdHtSD7wIPAc8D6xZcd/SziI9thAnkNIVRRONpcqAy4BrgQuB8yjRSQFFLt2TtQNryZXwc8CqFfctbS3Rc4mQkdIV/RKNpQaQm7/geuADwBg3nrfEpXuyNuBp4LGhjbsf/fivbjnk0vOKAJLSFXnr3i97DbmivRY4y+0MLpfu285bf8+mEQ2vNgC/BX4zd9tWKWCRF5l7QfRJNJaqAD4E3AhcDQwym8gAnd0/ouHV+d0/XQp8b+ucuU8CSeC/527b2mEunPALKV3Rq2gsNRH438BngPGG4xg13H5tJ7khbG8pJ7fFfw1Qv3XO3Ae2TeLHH3p66y4jAYUvSOmK04rGUhcDXwKuQ14nAEzZ+0xvu1FGArf/19VlH/h60toO3FO7vPYpl6IJH5E3k3hbNJZS5A6GfRWZyPtE2jky+kit1dtN7MHU7B6nzgNmANVW0tpEborHB2uX18quBwHIacCC3FCvaCx1C7lTVf8bKdxTDGnet0Whex1f/NvLIicX63zgp8AOK2l93Epa8n4TUrphF42llpGb+/XnwByzabxryr6VA3u7Pqs48PQidX4PV08ld7BtnZW03lP0cMJXZPdCSEVjqVnAd8kN+RK90bpp3MF1C3q7yYtz1XYnoiac4ZEWAf9jJa2ngK/WLq+tKVZE4R9SuiETjaVGAN8AbiM3gYw4g0FthzZFdKbHNc40dP3iysi5eTzkMnJbvb8Cvla7vDbd34zCP6R0Q6J7HoRbgTgwymwaf5m0/7lezyDaP5I1DUNUvvvBFfAx4O+spPUj4K7a5bV2oRmFf8g+3RDo3m+7EfgPpHDzo3XHxAMv9DpqIbks0ufZ0U6jErgdqLWS1pX9eBzhE7KlG2DRWKoSuBv4PLktK5Gnyo6GjeXZjp4OkNFRzo710yO97u/toynAk1bS+iFwR+3y2rYiPKbwINnSDahoLDUfWAN8ASncgk1448Vex9c+fr46WMSnU+R+QdZYSavHohf+JqUbQNFY6vPkCrfXj8XiDLTOTt63am6PV4P9yMWRUqx9Nht4wUpad1lJSz6NBoyUboBEY6mx0VgqBfwA6HVcqTiziq6W2gFdzT3uA98ylQ0dA0q2tHo5cCfwkpW0eix+4T9SugERjaWuIXew7H2mswTFuENrGnu6ToP++VVlU1yI8U7gFStprSj2AyulypRSNUqpPxX7sUXP5KOLz3UPBfs28EVk321RTal7ZkZP19mDqdk7Vrm1rPpA4IfdW7xfrF1emy3S436R3Fpxrs+HHGaypetj0ViqCngU+AekcIuqLNO+ZVB7/cServ/NuyJdbubptgL4o5W0hvb3gZRSk4FqcisWCxdJ6fpUNJYaC/wV2Z1QEmPe3NDjihBZxf6VC9ViN/Mc5xrgb1bSmtrPx7mH3GxyTr8TibxI6fpQ97wJLwKm3viBN3XvU2f3dN3f5qkdTkT1e0XjfrCA1VbSKuj/Xyl1LXBIa72uuLFEX0jp+kw0lroIeAE4x3SWoIo4XbuGtOyfdrrrNHQ+uDSveRZKZTywykpa1xdw30uADyil0sDDwFKl1C+LGU70TErXR6Kx1HXkVqWVU3lLaGT91rqerts3ijX2EOXKisd9MBj4nZW0vprPnbTW/6S1nqy1jgI3AM9orW8qRUBxKildn4jGUiuA3xPGBSFdNnXvU2N7uu6BZRGvHelXwL9ZSet7poOIvpHS9YFoLBUHfoj8f5WccrL7htu7TnsyQnsF2zeeE/HqWX5fspLWt/O9k9b6r1prmVPZRfIm9rhoLHUHuflvhQuG2bt29nRd6nzV44gGj/hHK2klTIcQvZPS9bDuXQryJnLR1Lpnhp/ucg32Hy6OvNPlOIW4w0pa/9d0CNEzKV2PisZSnyA3/61wi3YOjzqy6bS7DzafrdZ3VqjBbkcq0D9bSesrpkOI05PS9aBoLPV+cmcKyVlmLhravHebQp/ynsjNsxDpcdyuR91tJa1bTIcQp5LS9ZhoLLWE3NhJk4PvQ2ly3V9POzKkoYpX6saoqMtxiuGnVtL6oOkQ4kRSuh4SjaVmAn8iN/5SuEnrxnGHTr/i78OXR4o1wYzbyoCHraR1qekg4hgpXY/onkvhCWC06SxhNLj14KaIzg44+fJshLq/LjA2z0IxDAR+ayWt8aaDiBwpXQ+IxlIR4NfIqb3GTNr/7Gn3nz8/T+3SSvn9fTKB3Bav7LLyAL+/mILia8BS0yFCS+v2CW+8dMqoBQ0dv/TGPAvFcDnwr6ZDCCld46Kx1BXIyQ9GVXYc3Vie7ThlGfW60ay1q1SQdvd8xUpaHzIdIuykdA2KxlJjgIeQ/wejJh54ofN0lz+wLDLM7SwuuN9KWj2uiCFKT97shkRjKQU8CPS4OoFwgdbZSfufnXfyxe0VbK2dFplvIlKJDQN+byUtmTjJECldc2LA1aZDhF1FV3PtgK6WkSdf/scl6k0TeVyyALjXdIiwktI1IBpLXQLcZTqHgPEH19gnX6ah4dELI34eJtYXy62k9SnTIcJIStdl0VhqFLkzzmQlZtO01pPrVs46+eLaqNrQWaHC8PH7uzJ+131Suu77HjDZdAgBZdn2LYM66iccf5kGff9VkaihSG4bBnzfdIiwkdJ1UfduhZtN5xA5Yw5vOGW/7dEhrNs3Wvltcpv+uMFKWleZDhEmUrouicZSZcCPkJnDPOPsvU+eUq4PXx4J45LkP7aS1kDTIcJCStc9twILTYcQOZFs586q1jeix1+WiVC3yvL1PAuFmgH8H9MhwkJK1wXRWGo08E3TOcQxI+u37Dv5sufmB2KehULdYSWt2aZDhEFYX2Bu+1dghOkQ4pipe58ad/zPGjoeuiKQJ0P01QDgx6ZDhIGUbolFY6nzgU+aziGOUU62bnjj63OOv2zPGNY0VqlRpjJ5xFIrad1kOkTQSemWUPepvrJ0uscMt3fuOvmyB66KyCeRnO/IKcKlJWVQWrcAF5gOIU40Ze/TJxRs2wC2bD47MFM49tc44DOmQwSZlG6JdA8Ru9N0DnES7RweVb/lhH23f1wSqTcVx6O+YiWtU1bREMU
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Parameter 1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADsCAYAAADXaXXTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAg7UlEQVR4nO3deXyU1aH/8c+ZJQkQCAQUZXMAFwZZXECRqiDaFbF2sa2t1mrtr2mx1fbaOu3trWOt1dafXrcUEPdWq63V1jrW9mpcquwIksAEZF9lERgIBLLMuX9M7EUByUxm5szyfb9e86qBzPN8yyv55uQ85zmPsdYiIiLZ4XEdQESkmKh0RUSySKUrIpJFKl0RkSxS6YqIZJFKV0Qki1S6IiJZpNIVEckila6ISBapdEVEskilKyKSRSpdEZEsUumKiGSRSldEJItUuiIiWeRzHUCkPQKhiA/oAngB78Q9foY2+yzQCrQADZOnTtDm0JLzjDYxF5cCochRwMC21yAgAPQGugM92l7dga4Hvu8bu0rfPCru+dgBfxQHYsAOYHvb/+4AtgCrgBVtr5WTp07Ym6n/P+lkjPkBcDVggVrgSmvtPreppKNUupIVgVCkC3AqcFrbayRwPFCeyvEOUbrJeJdEAS8C5rW9Fk+eOqE1xeOlnTGmL/AGMNRa22iM+SPwgrX2EbfJpKM0vSAZEQhFBgIXAOcCo4ATyZ1rCMe0vQ4s7b3VVTULSRTwTODlyVMnbHWQ7UA+oJMxphnoDGx0nEfSQCNdSYtAKNILmECiaC8gMV2QMR0c6baHBRYA/wD+Cbw5eeqE5gye7yDGmGuBW4BG4J/W2q9l8/ySGSpdSVkgFOkDfBG4hMSo0WTr3Fko3Q9rAGqAp4G/TJ46YXcmT2aM6QH8GfgysBP4E/C0tfb3mTyvZJ5KV5LismgP5KB0D9QIPA/8AXhh8tQJ+9N9AmPMJcCnrLXfbPv468AYa+13030uyS7N6coRBUIRL3AhUAV8gtyZm3WlE4kfOpcAseqqmmeAByZPnTAjjedYC4wxxnQmUfLnk5hvljynka4cVtuo9uq2V3/HcT7A8Uj3cOYD9wJPpmP0a4y5icT0QguJ+eWrrbVpH1VLdql05SCBUOQM4MfAZ8nR34ZytHTftxWYDkyZPHXCetdhJLeodOXfAqHIucB/kVh9kNNyvHTf1wI8Afxi8tQJK1yHkdyg0hUCocgngJ8B57jO0l55UrrvawEeA26ePHXCasdZxDGVbhELhCJnA7cDY1xnSVaele77moGHgV9Onjphnesw4oZKtwgFQpHjgd8An3OdJVV5WrrvawLuJjHyzeh6X8k9Kt0iEghFupKYs70WKHEcp0PyvHTftxn4CfCIdkgrHsW+3rJoBEKRLwPLgB+R54VbQHoDDwFvVFfVjHAdRrJDI90CFwhFjgGmABc7jpJWBTLSPVAriSmHn02eOqHRdRjJHI10C1ggFLkCWEKBFW6B8gI/BOZXV9WMch1GMkcj3QIUCEX6AfcDn3adJVMKcKR7oBbgZuBXk6dOaHEdRtJLI90CEwhFJgJvU8CFWwR8wE3Am9VVNSe6DiPppdItEIFQxBcIRX4D/A2odJ1H0uIMYEF1Vc03XAeR9NH0QgFom054kg8+CaGgFfj0wqH8Frgu2xupS/pppJvn2m7hXUARFW6R+i5QU11Vc4zrINIxKt08FghFvgO8APRynUWy4mwSqxvOch1EUqfSzUOBUMQTCEXuJPErp9d1HsmqPsCr1VU1V7gOIqlR6eaZQCjSmcSzs37gOos4UwI8Ul1Vc4PrIJI8XUjLI4FQpDeJ1QmjXWdxrQgvpB3OncD12rshf2ikmycCoUhf4HVUuPJBPwQera6qycknfMjBVLp5IBCKDABeA7RQXg7lcuC56qqazq6DyJGpdHNcIBQJkCjcwY6jSG77NPDX6qqaUtdB5KOpdHNYIBQZTKJwA46jSH64AHi6uqrG7zqIHJ5KN0e1TSm8CgxwHEXyy4XAE9VVNVpKmKNUujkoEIpUAv8A+rnOInnpi8DD1VU1xnUQOZhKN8cEQpFOwPPAENdZJK9dTuLmGckxKt0cEghFvMBTgG7zlHSoqq6qud51CPkglW5umQpMch1CCsqvq6tq9DWVQ1S6OSIQivwHcLXrHFJwPCQurI10HUQSVLo5IBCKnAf82nUOKVjlJG6e6O06iKh0nWvbgPwptFuYZNYAEjdPlLgOUuxUuikwxnzKGLPUGLPcGBNK9TiBUKSUxI5hR6UvnchhnYl+o3JOpZskY4wXqCZx2+VQ4FJjzNAUD3cPiedgiWTLddVVNRNdhyhmKt3knQEst9autNY2kXg22WeTPUggFPkS8P/SHU6kHR6prqrp4zpEsVLpJq8vsO6Aj9e3/Vm7BUKRY4Ep6QwlkoRewO+rq2r0/e+A/tHdeBA9Jl3cOg/4iesQxUilm7wNQP8DPu7X9mftEghFvk1iPljEtXB1Vc0I1yGKjUo3eXOBE4wxA40xJcBXgOfa88a2rRrvyGQ4kST4gOmaZsgu/WMnyVrbAlxDYhewKPBHa+3idr79QaBLprKJpOAMEl/PkiV6MGWWBEKRK4BHXOcoFHowZVo1AEMnT52w7oifKR2mkW4WBEKRHsDtrnOIHEY5ibXnkgUq3ey4Gd11JrltUnVVzeddhygGKt0MC4Qiw4Aq1zlE2uF27c2QeSrdzLsLbWYj+WEQ8D3XIQqdSjeDAqHIBcD5rnOItFfJ/tjF0SHBHq5zFDKf6wAF7mbXAUTaw9Pa9M5Jy57ceezm2WcDPwJ+6jpTodKSsQwJhCITSTxgUjJAS8bSxLZuGrj6xRWBNX8fa7Dv/+a7BxgcrI9udhmtUGmkmwGBUMQAv3CdQ+SwrI0d++7MBSe+89QYb7zl7A/9bRcS+zJcl/1ghU9zuplxMXCa6xAiB7F2f4/t0dfOefOGeHDp4+O98Zayw3xmVXRIsP9h/k46QCPdzLjRdQCRD7A23mXPppkj6qYd12nftnHteEcp8H0S87uSRprTTbO2h0zWuM5R6DSn234l+3fOG143vWvF7tUnJfnW7UC/YH20MRO5ipVGuul3resAIgDeln3RYP1j+47e9vaoFA9RSWIXvYfTGKvoqXTTKBCKDAQmuc4hxc3EW9YOXvmX9f3Xv3KWAdPBw01GpZtWKt30ugZdnBRXbHxbvw2vLT5+xTNjPTY+IE1HPT06JHhmsD46O03HK3oq3TQJhCLlwDdd55AiZO3eXtvenjO0/rHTfa3723ORLFnXACrdNNGoLH0uBSpch5AiYm1Lt12r/zV25s92j1g8fbyvdX/XDJ3pkuiQoHbJSxONdNPnctcBpHiUNW6bNaJuWu/yPRvPycLpSoGrgVuzcK6CpyVjaRAIRQYAq+n4RQtpp2JdMuZr3rvo5CUPmZ47osOzfOq1wKBgfbQ1y+ctOBrppsfXUOFKBnnizStOfOeP2/psmnGmowgDgHFoDXqHqXTT4zLXAaRA2fi7gTX/eGfg6shYgx3sOM1FqHQ7TKXbQYFQ5BRgqOscUmCs3dV785y3hix78gxvvCkb87btMQltgtNhKt2O+5LrAFJArG3qvnPZzGFLHhxe0rxnvOs4HzIoOiQ4LFgfrXMdJJ+pdDtuousAUgCstZ33vjtjRN20AZ0bt2ZirW26XASodDtApdsBgVCkLzDCdQ7Jb/6mXW8Nr5veqfuulfmwGuMi4FeuQ+QzlW7HfNp1AMlfntb9S4P1v9/Te+tb+bT38hnRIcHeeqpE6lS6HaPSlaSZeOv6QaueWzNg3Utj07AhTbYZ4ELgQddB8pVKN0WBUMQPXOA6h+QRa7f33fiv2hOW/3mMx7b0cx2nAy5CpZsylW7qzgS6uQ4hecDaxp7b62YPXfLoqf7Wxly+SNZeF0SHBEuC9dEm10HykUo3dWNdB5AcZ21r14Z1M4bX3X982f4d413HSaPOwDDgLddB8pFKN3WubseUPFC6b/ucEXXTenZtWJ8rNzak2+modFOi0k3dGNcBJPd4WxrrTl7ycGuv7YvPcJ0lw/JpxUVOUemmIBC
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Parameter 2\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADsCAYAAADXaXXTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAe2ElEQVR4nO3deZhU1YE28PdUdTe9F3QDsm+y1AUKRFBQBPxQo4lbFo1L1Io6Os+ExIzGTHwSYypjMlHzxfmcqGlAxXYZE4P6GXXiEtxQEBFBL3CLfZd96a7eu6vO/HGbiESgq+ree+7y/p6nH4R03/umod8+feqec4SUEkRE5IyQ6gBEREHC0iUichBLl4jIQSxdIiIHsXSJiBzE0iUichBLlzxBCHGBEGKNEGK9EOIO1XmIciX4nC65nRAiDGAtgPMAbAewFMBVUsrVSoMR5YAjXfKC0wGsl1JulFK2AfgjgEsVZyLKCUuXvKA/gG1H/H57558ReQ5Ll6iLhBC3CiFWCSFWCiGeEUIUq85E3sPSJS/YAWDgEb8f0PlnjhFC9AdwC4BJUsqxAMIArnQyA/kDS5e8YCmAEUKIoUKIIphl9xcFOQoAlAghCgCUAvhMQQbyOJYuuZ6UsgPA9wG8BsAA8KyUcpXDGXYA+L8AtgLYCaBOSvm6kxnIH/jIGFEXCCF6AHgOwBUADgH4M4D5UsqnVOYi7+FIl6hrzgWwSUq5V0rZDuB5AGcqzkQexNIl6pqtAKYIIUqFEALAOTCnOoiywtIl6gIp5RIA8wF8DECH+bUzR2ko8iTO6RIROYgjXSIiBxWoDkAEAEhEBICTAPQC0LPz1yPfegIoh/lvtqBOlh18avdTPQFkAKQBNAOog/lkweG3gwD2AtgCYNOsmpkpx/7/EB0DpxfIWYlIGMBwABqA0Z1vGoAozAUHXdIqCzY+svvPw7K8+0EAmwBsBrAewCcAlgNYM6tmZibLaxHlhKVL9kpE+gKY0fl2JsxyLcr3sjmW7rE0AfgUZgEvA7BwVs3MtRZdm+gLWLpkrURkED4v2RkwR7WWa5UFmx7Z/eehdly7004A7wBYAOCNWTUzt9h4LwoQli7lLxEZD+ByAJcBGOXELR0o3aOtgbkg4rlZNTOXOXhf8hmWLuUmEZmAz4t2hNO3b5MFm+Y6W7pH2ozOAgaweFbNTH4RUZexdKnrEpEhAG6Cuf/AySqjKC7dI20F8CiAx2bVzNyuOgy5H0vXA4QQAwE8AfORKglgjpTyAUdubj7KdR7MXb4uhEue7XZR6R6WBvBXAHMBvDKrZmZacR5yKZauBwgh+gLoK6X8WAhRAfMV9q/bejBjIlIJ4HoA3wMw0rb75KhNFmyeu/vPQ1TnOIbPADwE4OFZNTMPKc5CLsPS9SAhxIsAHpRSvmH5xRORwQB+AuBamIsRXMnlpXtYPYDZAP5zVs3MnarDkDuwdD1GCDEEwLsAxkop6y27sPk87c9gztnm/Ryt3dpkeMvc3fMHq87RRa0wp4fum1Uzc73qMKQWS9dDhBDlMJ8d/bWU8nlLLpqIVAO4A8AsACWWXNMBHivdwzoAzAOQmFUzk0f9BJQrXhShExNCFMJ8ROlpSwo3EYkgEfl3mMtib4eHCreTUB0gBwUwf5JYv2jyxXcYUc210zdkH5auB3Rumv0oAENKeX9eF0tEBBKRGwFsAPBzABX5J6RsdGs5oPeoW/8bAGuNqBY3opoXv4FQjli63jAV5gtbM4UQKzrfvpb1VRIRDcDbAB4BUG1pQocJ89E575EyPU6vOfy57wvgcQBLjKh2hrpQ5CRu7egBUsr3kM+P04lIMYA7AfwYHniRrIs8OTqsTG1ZVNG4Y9pRf3wagEVGVHsCwC1a0qhTEI0cwpGu3yUi58I8XuZn8E/hepOUDbGVc463N8V1AFZw1OtvLF2/SkTKkYjMA/AGbNrpSzHPjXR7712+rFtbXe8TvNsQAAuNqHaXEdXCDsQih7F0/SgROQ3ACgDfVRuE/k5m9mjJJyd18b3DAH4J4C0jqg20MRUpwNL1E/PJhH8D8D4Ub0jjAE+NdAdvfWNNONNWluWHTQPwqRHVLrcjE6nB0vWLRKQKwEsA7gVQqDgNHUFk2jcO2/TSmTl+eHcAzxpR7VEjqmVb2uRCLF0/MKcTPoa5Cxi5zKi1f9onIPOdn70BwDIjqp1qRSZSh6XrdYnItwEsBOC1JbF5EZCemF4oaG/8pN+uxadbdLlRABYbUe0mi65HCrB0vcycv/0jgG6qoyjg/sURUsrYqrlWPwtfBGCOEdV+avF1ySFcHOFF5jHmvwfwL6qjKOT6kW5J854PehxaZ9czt782olo1gNu1pOH+b0D0dxzpek0iUgbgRQS7cN1Pyrbxek1/m+9yG4B5RlTj4MlDWLpekoj0gbm1I18wc7nuh9Z9UNq8Z5ADt4oDeM6IasUO3IsswNL1ikSkP4D3AExUHcUl3Du9IGXd2NWPjHXwjpcAeNWIapUO3pNyxNL1AvNUhzfh/wUPvtB31wcritobqxy+7QwAbxtR7UTLjEkxlq7bJSK9ASyACw+HVMyVI12RSe8Yue6PkxXdfgKA94yoFqjHB72GpetmiUhPmIWrqY5CXTN08yubw5kOlfOrIwC8b0S1YQoz0HGwdN3KXNb7NwBOzg16hnDhSDeUbl0zeOvrbtiWsT+Al42oFlEdhP4RS9eNEpEKAK8DGK86CnWdlnyqQUC65WtKAzCfj5O5j1v+gdBhiUgI5iozPqVwfK5aEFDUWrfspL0fu+3v7FwAD6kOQV/E0nWf3wLI/vyz4HHP9IKUmXErZ7v1gM+bjah2m+oQ9DmWrpskIjfAXGVEJ+SegW55447Flaktbn665LdGVLtEdQgysXTdIhGZBuAPqmN4iDtGulI2x/TZbn9SIATgv42oNkF1EGLpukMiMhTA8+DBkZ5TvX/lhyWtB/qqztEFZQBeMqJaP9VBgo6lq5q5gc1fAPRUHcVj1I90ZWb/GONxL40e+8Ms3lLVQYKMpaveA+CzuJ40YMc7qwrSLV7b7+BUAE+qDhFkLF2VEpFvALhRdQzKnsh0bBm+4Xk3LITIxTd5+oQ6LF1VEpF+AOaqjuFhSqcXhm94fmdIZrx8AOj9RlQbqjpEELF0VUhEBIDHAVQrTkI5CHc0rxq4450pqnPkqRxArRHV2AEO4ydcjR8COE91CC9TuffC2NWPpVXd22LTAPxIdYigYek6LRGJAbhHdQwfULI6orhl/4fVB1aPU3Fvm9xtRLUxqkMECUvXSea+Co8imKf3Ws35ka6UHeP0ml6O39de3QDMNqKa+kfwAoKl66wbAZymOgTlprJ+0+Lyxs/8+OLTVAA3qw4RFCxdpyQiPQD8h+oYlCMpG2Kr5kZVx7DRPUZU66M6RBCwdJ3zK3DVmZUc/XG4955ly7q11fttauFI3WEu1CGbsXSdkIicAuCfVcegHMnMbm3N05NUx3DAt42odoHqEH7H0rWb+UzuQwDCqqP4jGMj3cFbX18bzrSVOXU/xX6tOoDfsXTt9x0AZ6oOQbkJpds3DNv0cpD+/k41otpFqkP4GUvXTolIAYBfqo7hT9KRke6odc8cEJBB+ynlF6oD+BlL117XAHD7Btd0DIVtDSv67loSxEf8JhlRjUdG2YSla5dEJAzgp6pjUI6klGNXzQ3ypvJ3qQ7gVyxd+1wJYITqED5m6/RCadPuxT3q1o+28x4uN9mIauerDuFHLF07mMt9f6Y6BuVIyrZxK2sGqo7hApzbtQFL1x6XA9BUh/A520a6PQ6tWVzavJelC5xhRDXuhmcxlq49OMq1mW2NK2Xd2FWP+WkXsXxxtGsxlq7VEpGzAcRUxwgAW3q3785FKwo7GnvYcW2PmmpEtZmqQ/gJS9d6PHvKo0QmvX3Uume9fiKEHe5QHcBPWLpWSkSqAHxLdYyAsHwT82GbXtoakh3c6/gfnWNEtb6qQ/gFS9da14EblDvF0umFcLrVGLTtDa+e7mu3EMxHIMkCLF1rcWrBozSjtkXluWsecLXqAH7B0rVKIjIVQJAfpneaZQVZ1Hroo97
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "for i, item in enumerate(beta):\n",
    "    print('Parameter', i)\n",
    "    item.plot_piechart()\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting with priors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When extracting energy levels and matrix elements from correlation functions one is interested in using as much data is possible in order to decrease the final error estimate and also have better control over systematic effects from higher states. This can in principle be achieved by fitting a tower of exponentials to the data. However, in practice it can be very difficult to fit a function with 6 or more parameters to noisy data. One way around this is to cnostrain the fit parameters with Bayesian priors. The principle idea is that any parameter which is determined by the data is almost independent of the priors while the additional parameters which would let a standard fit collapse are essentially constrained by the priors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first generate fake data as a tower of three exponentials with noise which increases with temporal separation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbkAAAEKCAYAAACPCivzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAceUlEQVR4nO3df5TV9X3n8dd7BrFQdYgyjMrgYmCgoFFYZyUxnh7qlCz+QLqNiZI0m91Qqdv6s+52id0ct9tN4+45WQzGk5QI1XC6amqzLSS0CWcSj8khNUAlCcjCANLlhzAEwyCFiMx97x/3e/HeO3N/zNwf38/93ufjHA9zv/fX+3yT4c3n83l/3h9zdwEAkEQtcQcAAECtkOQAAIlFkgMAJBZJDgCQWCQ5AEBikeQAAIk1Ju4Aipk4caJPnTo17jAAoGFs3br15+7eHnccoQg6yU2dOlVbtmyJOwwAaBhm9k9xxxASpisBAIlFkgMAJBZJDgCQWEGvyQEA4rN169ZJY8aMeUbStQpzUJSStP3cuXO/e8MNN/QP94K6JTkzmyXpIUkTJfW6+1fq9d0AgJEbM2bMM5dffvms9vb2X7S0tATXzT+VStmxY8dmHzly5BlJdw73mooys5mtMbN+M9ued32hme0ysz1mtlyS3H2nu98n6eOSPlzJ9xYzmHL17jyqlb196t15VIOp4P53AYBGcW17e/vJEBOcJLW0tHh7e/uA0iPNYVU6kntW0pclfT1zwcxaJT0taYGkg5I2m9k6d3/dzO6U9B8kra3we4c1mHJ9avWr2nbghM6cHdS4sa2aM2WC1i6dp9YWq8VXAkCStZSb4P70W69fufqHb1yRf33pzVe/+bk7Zh+ufmhpUXwFB2wVjeTc/RVJb+VdvlHSHnff5+5nJb0gaXH0+nXufqukTxb6TDNbZmZbzGzLsWPHRhTPwy++pk17j+v02UG5pNNnB7Vp73E9/OJrI/ocAMDIfO6O2Yf3P3H71us7205d39l2av8Tt2/d/8TtWytNcB/72MemXnrppdd3dXVdM5r312IhcbKkA1mPD0qabGbzzWylmf25pA2F3uzuq9y9292729tHtmm/a9LFyh+vmaQZky4e0ecAAEbuXCqlgTPvth4Z+OXYddsOt51LpSr+zM985jM/X7duXd9o31+3whN3f1nSy7X8jr7+t5U/rnZJu/vfruXXAkDTO5dK6WNf/VHX/3vr9LiUS3/01z95/19seuPUX933ob4xLaMfT916662ndu3aNXa076/FSO6QpClZjzuja2Uzs0VmtmpgYGBEX/zk3XN107TLNH5sq0zS+LGtumnaZXry7rkj+hwAwMhs+OmRtp1vnrwoU+v3y3dTLTvfPHnRhp8eaYszrlqM5DZL6jKzq5VObvdI+sRIPsDd10ta393dfe9I3tfaYlq7dJ5e3tWv1w+f1OwrL9H8mZMoOgGAGvvZoYHx77ybyhk4vfNuqmX74YHxd865cmQjliqqKMmZ2fOS5kuaaGYHJT3u7qvN7H5J35HUKmmNu+8Y4ecukrRo+vTpI46ptcXUM6tDPbM6RvxeAMDofGBy2+kLL2hJ/TIr0V14QUvq2ivbTscZV6XVlUvc/Qp3v8DdO919dXR9g7vPcPdp7v75UXzuendf1tYW6ygXAFCm2667fGDWFZecykyc/coFLalZV1xy6rbrLo9tFCeF2aYFANBgxrS06K/u+1DfVZeOP9NxyYVn/+dHr99XadGJJC1atOjqm2+++dfeeOONCzs6Oq5bsWLFxBHFVdG310gl05UAgHiMaWlR27gLBtvGXTBYrXW49evXv1FRTNUIotpGW3gCAIhHfseTqcu/fYNU+44npQSZ5GppMOV6eVe/dhw+qWuovgSAqvjcHbMPx5nMCgkyydVqupLelgDQXIIsPKlVdeXLu/q17cCJnN6W2w6c0Mu7hj2GCACaXSqVSgU9AojiK9g/LMgkVytf+8E+nT47mHPt9NlBPfODitY1ASCpth87dqwt1EQXnSfXJml7odcEOV1ZK+0XXzjs9YkXj7otGgAk1rlz5373yJEjzxw5ciT4k8ELvSDIJFerNbkn756r46fODlmTo7clAAx1ww039KvAiduNwtyDPPBVktTd3e1btmyp6mdmqivpbQkgicxsq7t3xx1HKIIcydUSvS0BoHmEOMcKAEBVkOQAAIkV5HRlnL0r6YgCAMkRZJKLq3clHVEAIFmYrszy8IuvadPe4zkdUTbtPa6HX3wt7tAAAKNAksvSNeli5Y/XTNKMSRfHEQ4AoEIkuSx9/W8rf9egS9rd/3Yc4QAAKhTkmlxc6IgCAMkSZJKLq7qytcW0duk8OqIAQEI0XVsvAEgy2nrlYk0OAJBYJDkAQGKR5AAAiRVk4UmoaPkFAI2FJFemL353l5763p4h1x+4Zboe/cjMGCICAJTCdGWZ5kyZoPFjW3OujY/20QEAwhRkkjOzRWa2amBgIO5Qzttx+KTOnB3MuXbm7KBeP3wypogAAKUEmeTcfb27L2tra4s7lPNo+QUAjYc1uTLR8gsAGg9JrkzltPyi+hIAwkKSG4HWFlPPrA71zOoY8hwHrgJAeIJck2tEHLgKAOEhyVUJB64CQHhIclVyzZWXaFzePrpxY1s1+8pLYooIAECSq5JtB07odN4+utNnB7XtwIl4AgIAcJ5cNWWqKzlwFUBcOE8uF9WVVVSs+hIAUH9MVwIAEqtuIzkz+y1Jt0u6RNJqd/9uvb4bANCcKhrJmdkaM+s3s+151xea2S4z22NmyyXJ3f/G3e+VdJ+kuyv5XgAAylHpdOWzkhZmXzCzVklPS7pV0mxJS8xsdtZL/kv0fNMZTLl6dx7Vyt4+9e48qsGUj+h5AMDIVDRd6e6vmNnUvMs3Strj7vskycxekLTYzHZKekLS37n7Pxb6TDNbJmmZJF111VWVhBeUUm2/aAsGANVXi8KTyZIOZD0+GF17QNJvSrrLzO4r9GZ3X+Xu3e7e3d7eXoPw4lGq7RdtwQCg+upWeOLuKyWtLOe1ZrZI0qLp06fXNqg6Srf9ejPnTLrstl+lngcAjFwtRnKHJE3JetwZXStbiIemVqpU2y/aggFA9dUiyW2W1GVmV5vZWEn3SFpXg+9pKKXaftEWDACqr6K2Xmb2vKT5kiZKOirpcXdfbWa3SXpSUqukNe7++RF+bma68t6+vr5RxxeaUm2/aAsGoFK09cpF70oASBCSXC7aegEAEivIJGdmi8xs1cDAQNyhAAAaWJBJLonVlQCA+gsyyQEAUA1BnieXxM3gtZapzNxx+KSuoTITACQFmuTcfb2k9d3d3ffGHUtICiWyL353l5763p4hr3/glul69CMzY4gUAMIQZJLDUMUS2ZwpEzR+bGvOZvLxUYPnDEZ6AJoRSa5BFEtkOw6f1Jm8bilnzg7q9cMn1TOrgxMOADStIAtP2EIwVLFEVqrvJSccAGhWQSY5thAMVSyRlep7mT7hIBcnHABoBkxXNohiiezRj8zUw785o2Dfy0yCzH4/JxwAaAYkuQZRKpG1tph6ZnWoZ1bHkPcWS5DDvR4AkoIGzU2CEw6A5kCD5lxBjuTYDF59xUZ6ElsMACRTkEmOzeD1xWZyAEkVZHUl6iuzBy9b/mZyAGhEJDkU3YOXMZhy9e48qpW9ferdeVSDqXDXcgEgI8jpStRXqS0GTGcCaFSM5FByMznTmQAaVZAjOaor66vUHrxSvTEBIFRBJjmqK+uv2BYDOqYAaFRMV6KkUtOZABAqOp6gLHRMARoDHU9yBTldifBU0jGFbioA4kKSQ8WKbTGQxPYDALEhyaFixU4tz/xc6DkAqCUKT1CxYlsMyummAgC1EuRIjn1yjaXUFgO2HwCIS5AjOXdf7+7L2tra4g4FZSi2xYDtBwDixBYCVEWxLQZsPwDqhy0EuUhyAJAgJLlcQU5XAgBQDSQ5AEBiBVldieZSqiMKHVMAjBZJDrEqdSArB7YCqATTlYhVqQNZObAVQCVIcohVqY4odEwBUAmSHGKV6ZaSLbsjSqnnAaCYuiU5M3u/ma02s5fq9Z0I3/y
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m1 = 0.18\n",
    "m2 = 0.5\n",
    "m3 = 0.8\n",
    "\n",
    "A1 = 180\n",
    "A2 = 300\n",
    "A3 = 500\n",
    "\n",
    "px = []\n",
    "py = []\n",
    "for i in range(40):\n",
    "    px.append(i)\n",
    "    val = (A1 * np.exp(-m1 * i) + A2 * np.exp(-m2 * i) + A3 * np.exp(-m3 * i))\n",
    "    err = 0.03 * np.sqrt(i + 1)\n",
    "    tmp = pe.pseudo_Obs(val * (1 + err * np.random.normal()), val * err, 'e1')\n",
    "    py.append(tmp)\n",
    "    \n",
    "[o.gamma_method() for o in py];\n",
    "\n",
    "pe.plot_corrs([py], logscale=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As fit function we choose the sum of three exponentials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func_3exp(a, x):\n",
    "    y = a[1] * anp.exp(-a[0] * x) + a[3] * anp.exp(-a[2] * x) + a[5] * anp.exp(-a[4] * x)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can specify the priors in a string format or alternatively input `Obs` from a previous analysis. It is important to choose the priors wide enough, otherwise they can influence the final result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "priors = ['0.2(4)', '200(500)', \n",
    "          '0.6(1.2)', '300(550)',\n",
    "          '0.9(1.8)', '400(700)']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to chose a sufficiently large value of `Obs.e_tag_global`, as every prior is given an ensemble id."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "pe.Obs.e_tag_global = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fit can then be performed by calling `prior_fit` which in comparison to the standard fit requires the priors as additional input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit with 6 parameters\n",
      "Method: migrad\n",
      "chisquare/d.o.f.: 1.0925587749193326\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAEsCAYAAAA8UOGyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6hUlEQVR4nO3deZxU1Z3//9enqvduoFkaBLoRVGRHQBT3EBEVw4h5jIlLYkyiw/h1iTH5JdFMJm7JxGQyQzSZLBgzMUaNicZRiQQJrnFBGmiRVRYhNCDdsvcC3V11fn/U7aabbqCXqr51q97Px6MeVXXuqVuf27frfu4599x7zTmHiIiIBEvI7wBERESk45TARUREAkgJXEREJICUwEVERAJICVxERCSAlMBFREQCKMPvADqiX79+bujQoX6HISIi0i2WLl36sXOuqK1pgUrgQ4cOpbS01O8wREREuoWZbTnaNHWhi4iIBJASuIiISAApgYuIiARQoI6Bi4hI+qivr6e8vJyDBw/6HUrC5eTkUFxcTGZmZrs/owQuIiJJqby8nB49ejB06FDMzO9wEsY5x65duygvL2fYsGHt/py60EVEJCkdPHiQvn37pnTyBjAz+vbt2+GehkAl8E8P2ed3CCIi0o1SPXk36sxyBiqBTx90wO8QREQkjYTDYSZMmND02Lx5M+eccw4Amzdv5oknnvAttkAdAy/KbvA7BBERSSO5ubmUlZW1KHvrrbeAwwn82muv9SGygLXAi3IifocgIiJprqCgAIA777yTN954gwkTJjBnzpxujyNQLfCCzCgcOgDZPfwORURE0kBtbS0TJkwAYNiwYTz77LNN0x544AF+/OMfM2/ePF9iC1QCB2D/DihSAhcRSSvz74SP3o/vPE8YBzMeOGaVtrrQk0WgutABOLDd7whERER8F8wWuIiIpJfjtJT90KNHDw4c8O/sKLXARUREOmH8+PGEw2FOO+00DWI7nv11IXqqBS4iIt2kqqrqqGWZmZm8/PLL3R1Sk0C1wCsPZsABJXAREZFAJfCKgxmwX13oIiIiHU7gZvYbM6sws5XNyv7TzNaa2Qoze9bMCptNu8vMNpjZOjO7pFn5pV7ZBjO7sz3frRa4iIhITGda4L8FLj2ibCEw1jk3HvgAuAvAzEYDVwNjvM/83MzCZhYG/geYAYwGrvHqHlPFwQyo2gkRXVJVRETSW4cTuHPudWD3EWUvOecas+o7QLH3ehbwB+fcIefch8AG4EzvscE5t8k5Vwf8wat7TJUHw+CiUF3R0bBFRERSSiKOgX8ZmO+9HgxsbTat3Cs7WvkxVRz0Bs2rG11ERNJcXE8jM7N/AxqAx+M4z9nAbIDLJg6KFe7f0Y50LyIi6WLOwg94cNH6VuW3TxvOHdNP7fR8w+Ew48aNo76+noyMDL7whS9wxx13EAodvf27efNm3nrrrYTfpSxuCdzMvgjMBKY555xXvA0oaVat2CvjGOUtOOfmAnMBLj5ngoMqtcBFRKSFO6afyh3TT+WqX70NwFP/enZc5tv8WugVFRVce+217N+/n3vvvfeon+mu24zGpQvdzC4Fvglc7pyraTbpeeBqM8s2s2HAcOBdYAkw3MyGmVkWsYFuzx/ve/bUhSGUqVPJRESklUjUsaemjm17alm0ZieRqDv+hzqgf//+zJ07l5/97Gc459i8eTPnn38+kyZNYtKkSU33CT/yNqNHq9dVHW6Bm9mTwFSgn5mVA3cTG3WeDSw0M4B3nHM3OedWmdkfgdXEutZvcc5FvPncCiwAwsBvnHOrjvfdDoMeJ6gFLiIiLUSijuseWcyGiiqiDm57cjkTSgp57IYphEMWt+856aSTiEQiVFRU0L9/fxYuXEhOTg7r16/nmmuuobS0tNVtRmtqatqs11UdTuDOuWvaKH7kGPW/D3y/jfIXgRc7+v30HAT7yjv8MRERSV2vrqugbOteGhvdNXURyrbu5dV1FUwbNSAh31lfX8+tt95KWVkZ4XCYDz74oEv1OipQ10IHoFcxbF/udxQiIpJEVm3fT21dpEVZbV2E1dv3xzWBb9q0iXA4TP/+/bn33nsZMGAA7733HtFolJycnDY/M2fOnHbV66hAXUoViCXwfeUQjfodiYiIJIkxg3qSmxVuUZabFWb0oJ5x+47Kykpuuukmbr31VsyMffv2MXDgQEKhEI899hiRSGwH4sjbjB6tXlcFMIGXQKQOqiv9jkRERJLE1BH9mVBSSOPh7rysMBNKCpk6on+X5ltbW8uECRMYM2YMF110ERdffDF33303ADfffDOPPvoop512GmvXriU/Px9ofZvRo9XrKjt8xlfymzx5sit9/Hvw5FVw4yIonux3SCIikiBr1qxh1KhR7a4fiTpmPPg6NYci3DtrDFNH9I/rALZEa2t5zWypc67NZBfAFrh3ldZ9W49dT0RE0ko4ZPTOy2Jw71ymjRoQqOTdGcEbxFboXf9FI9FFRMRz5JXYht75F6DrV2JLZsFL4Dm9ILsn7FULXEREYhqvxJZOgteFDodHoouISEoL0jitrujMcgY0gZfAvn/4HYWIiCRQTk4Ou3btSvkk7pxj165dHT4/PHhd6BBrgZe/63cUIiKSQMXFxZSXl1NZmfqnDefk5FBcXNyhzwQzgReWQO0eOFQF2QV+RyMiIgmQmZnJsGHD/A4jaQW3Cx10HFxERNJWQBN447ngSuAiIpKeAprAG1vgGsgmIiLpKZgJvMcJYGG1wEVEJG0FM4GHwtBzsC7mIiIiaSuYCRygcAjsVRe6iIikp+Am8N5DYc9mv6MQERHxRbATeNVHUFfjdyQiIiLdrsMJ3Mx+Y2YVZrayWVkfM1toZuu9595euZnZQ2a2wcxWmNmkZp+53qu/3syu73DkfbyT+/du6fBHRUREgq4zLfDfApceUXYnsMg5NxxY5L0HmAEM9x6zgV9ALOEDdwNTgDOBuxuTfrv1Hhp7Vje6iIikoQ4ncOfc68DuI4pnAY96rx8FrmhW/jsX8w5QaGYDgUuAhc653c65PcBCWu8UHJsSuIiIpLF4HQMf4Jzb4b3+CBjgvR4MND/Xq9wrO1p5++X1hawesPvDTgUsIiISZHEfxOZi932L273fzGy2mZWaWWmLO9KYaSS6iIikrXjdjWynmQ10zu3wusgrvPJtQEmzesVe2TZg6hHlr7Y1Y+fcXGAuwOTJk5t2DOYs/ICR23M4ZcdKpt/5l6b6t08bzh3TT+3yAomIiCSzeCXw54HrgQe85+eald9qZn8gNmBtn5fkFwD/0Wzg2sXAXR35wjumn0o0ehYNb/+Skl7Z3HPFOKaO6E84ZHFZIBERkWTWmdPIngTeBkaYWbmZ3UAscU83s/XARd57gBeBTcAG4GHgZgDn3G7gfmCJ97jPK2u3SNTx29WOLOqp27eD255cznWPLCYSjVvvvYiISNLqcAvcOXfNUSZNa6OuA245ynx+A/ymo9/f6NV1Fby1pydfDsEQq2BJXR/Ktu7l1XUVTBs14PgzEBERCbDAXolt1fb9bKjvB8CJoZ0A1NZFWL19v59hiYiIdIvAJvAxg3qyJ3MAEWeUWGzMXG5WmNGDevocmYiISOIFNoFPHdGfMSX92E4/hlgFeVlhJpQUMnVEf79DExERSbjAJvBwyHjshilUZgxiREYFP71mIo/dMEWj0EVEJC0ENoHPWfgBJ3/7RVYe6k9xdBs3PLqEk7/9InMWfuB3aCIiIgkXr/PAu90d00+NXbDlnX/AXxey+TtToKDI77BERES6RWBb4E36Do8971rvbxwiIiLdKPgJvN8pseePlcBFRCR9BD+B9yqBcLZa4CIiklaCn8BDYehzEuza6HckIiIi3Sb4CRxi3ejqQhcRkTSSGgm873DY8yFE6v2OREREpFukRgLvNxyiDbBni9+RiIiIdIvUSOB9vZHouzb4G4eIiEg3SbEEruPgIiKSHlIjgef1gby+GsgmIiJpIzUSOMQGsimBi4hImkidBN5/JFSuAef8jkRERCThUiiBj4baPVC10+9IREREEi6uCdzM7jCzVWa20syeNLMcMxtmZov
      "text/plain": [
       "<Figure size 576x355.995 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Obs[0.1861(27)], Obs[210(12)], Obs[0.701(60)], Obs[321(433)], Obs[0.711(51)], Obs[435(433)]]\n"
     ]
    }
   ],
   "source": [
    "beta_p = pe.fits.prior_fit(px, py, func_3exp, priors, resplot=True)\n",
    "[o.gamma_method() for o in beta_p]\n",
    "print(beta_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now observe how far the individual fit parameters are constrained by the data or the priors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADxCAYAAABoIWSWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUQklEQVR4nO3defAmRX3H8XfvAXKoIDcGbKKJEA9QKcUjRoVAdBTkMBK1iKh4Jho0SCNgBAUmYgWvKAJKFOJtVKA9gICJGoyLsrJAUlFxVK4gaqDAhV32N/ljZtlld1l+8zzPzLdn5vOqemprF6p+n91tPvSvZ7rblWWJiIh0Y4F1ABGRMVHpioh0SKUrItIhla6ISIdUuiIiHVLpioh0SKUrg+Oc+yfn3CrnXLnWZzvrXCKg0pXh+gVwOnBX/fNLDLOI3EelK73lnPPOuV+tM6OdA7YCtgGeCvwWmAN2MYwqch+nHWnSV845D/ys/unFwI3AK4F7qSYUK4CHAA4oyrLcrfuUIvenma703e0AZVkeAJxPNbNdDKwEtgVWAWX9ayLmVLrSdw/0rdqmwE3AQqqZ7o6dJRLZiEXWAUSmtBWAc+5rVCW7NXAnUACbUE0stgS+ZZJOZB0qXem726iWEZ6/1q9dR/UQbW3bdpZIZCNUutJ3cwBlWTrrICLzoTVdEZEO6ZUxEZEOaaYrItIhla6ISIdUuiIiHdLbC5IMH+Jiqk0MO63z2RHYgmq8LgIWH3z9t3/2jb0222T5HXvvQLXdd/XnDqrtwKs/NwA3FXm2ouPfjsgGqXSlUz7EBcDuwFPqz2OBnanKdVuq3WMPavHcqkth7qHA0+bxr5c+xNtYU8S/BK4FlgBLizy7p+nvQ2RSKl1pjQ9xIbAHawr2ycBeVLPWqVRv5S6Y77u5Dtiu/uy1zj9b6UNcRlXAV9Y/Xlvk2b3TZhTZEJWuzJQPcQ/gRUAG7A1s3soXckA5k0cSi6n+Z/Bk4HX1ry33IV4FfAe4ALiiyLO5WXwxEZWuTMWHuAh4FnAgVdk+pouvO8ccJa6tXWibAc+oP28HbvUhXgh8Bbi0yLO7W/q6MgIqXWnMh/hwqrMOXlT/uHXXGUrnoL3SXdf2wKvrz10+xG9SFfBFRZ79tqMMMhAqXZkXH6ID9gWOAl5MdYKXHedKZ3PcwhbAIfXnXh/i5cDHgS/rDQmZD5WubJQPcUfgVVSzvN83jnOfElyHM90Hsgj40/pzqw/xXOCsIs+ut40lKVPpygb5EPcG3gL8Odaz2g0rKef99kIXtgeOBd7uQ/w6cEaRZ5caZ5IEqXTlPvUSwqHA0VQPkZJVgptjQYo7Kh3wAuAFPsSrgfcDn9a7wLJaioNWDPgQnw9cBXyBxAsXVr+na7688GCeCHwCuN6HeFT93rKMnEp35HyI+/gQvwV8DdjTOM68lVDWrzD0wc7AWcAyH+KB1mHElpYXRqrexHAq1ZsIPeRcgx1pqdgD+KoP8dvA24s8+551IOmeSndkfIi7ACcBR1DdlNtLJSWJPUhr4o+BK3yIXwLeUeTZ/1gHku6odEeiPsHrBKodVg8xjjO1nqzpPphDgYN8iGcDJxR59hvrQNI+remOgA/xSVSHubyTARQuQAlQuiGM30XAG4BrfIgvtA4j7dNMd8Dq2e2JwHEM8O+6HNacYSfgwnqDxd8UeXaHdSBpx6BGrayx1uz2RIZYuI4SBjHTXdeRVG857GcdRNoxuP8Yx26ttdt3MOC/39LhZnS0Y4p2BS72IZ4JHFPk2V3WgWR2Bjtqx8iH+HiqQ7jfyYALF6DElW7Y49dRrfX+yIf4LOswMjtDHrSj4kM8GLiCHm1wmEa5ADeQB2kP5tHAv/kQj7MOIrMx6NnQGNTnJZwIvIt53i82EOWI5gwLgFN9iH8EvEbnOPTbaEbtEPkQNwc+T7XZYUyFS0k51AdpG/MK4Fs+xB2sg8jkxjZoB8OHuCvwXeAw6ywWSqDs7460aewDLPEh7mUdRCaj0u2h+sHKlax/s+14LHCMcKa72i7Ad+p1fOmZsQ7a3vIhvga4jOo68dEqwbV4MWUfbAF8yYd4vHUQaUal2yM+xGOBs6muDR+1krIcydsLG+OA9/gQz/Uhjv3Pojf0F9UTPsSTgNw6RypKB2X/jnZsyyuB83VIej+odHvAh/j3VBsepDbnGMqBN7PyF8BnfIh6DTRxGrSJ8yG+j+o4RllLSckAjnactZcAn1fxpk2lmzAf4inA26xzJMnhyjQvprR2MNVSg/5sEqW/mET5EE+kOrRGNmDO9eqOtK69FDin3q0oiVHpJsiH+DbgZOscKSuBgR94M60jgQ9Zh5D1adAmxof4YuB06xzJG8Z1PW17U/0dkyREpZsQH+ITgPMY2TkKE9PbC/Nxkg/xIOsQsoYGbSJ8iNsCFwBbWmfpBedAD9LmwwHn1SeUSQI0aBNQ3/bwRcAbR+mN0pXo7YV5eyjwVR/i1tZBRKWbig8Bf2IdoncGdjNlyx5DtXlCu9aMadAa8yG+EXiddY6+mQNGfMrYpA4ATrMOMXYatIZ8iM8FPmCdo79UuhM4xof4MusQY6ZBa8SHuD3wOXRl0mRcSclCjd/JnONDfJJ1iLHSoLXzEUZ+Ju5UynJUl6TN2GZUbzRsYh1kjDRoDfgQDwcOtc7RZ6VbMKejHafyOKoLTaVjKt2O1ZcKftg6R9+VTmN3BoKWGbqngdu9M4FtrEP0nd4Wm4lFwLn1e+LSEY3cDtVPjV9snWMI5lyppYXZ2BM4zjrEmKh0O+JD3BGd+jQzczrWcZaO9yE+3jrEWKh0u/Mx4BHWIQbDaWPVDG1CtcygP9QOqHQ74EM8GDjQOseQlE4z3RnbG/hb6xBjoNJtWT17ONU6x9CUpY6/bMEJ9aYdaZFKt31HArtbhxiaOR1g3oYt0RVRrVPptsiHuBnwLuscQ1Q6bQFuyet9iLtahxgyDdx2vRl4pHWIIZrTRLctm6KJQqtUui2pD4wO1jkGq9QW4BYd4UPUklhLVLrtOQ7YyjrEUJXWAYZtIfBu6xBDpdJtgQ/x94C/ts4xZLqTsnWH+hCfYh1iiDRy2/F3wEOsQwzZnG5MbptDrzq2QqU7Y/V23yOscwxd6TTV7cD+PsRnWocYGg3c2Xs91bZKaVGp82668hbrAEOj0p0hH+KmVKUrLdM24M4cXD+jkBlR6c7W4cAO1iHGYM7plbGOLALeYB1iSFS6s/VX1gHGotSDtC4dVX8XJzOg0p2R+tqTva1zjIXuR+vUdsAh1iGGQqU7O6+1DjAmJXp7oWOvtg4wFBq4M+BD3AJ4mXWOMdEdaZ17ng9xN+sQQ6CROxsvAR5mHWJMSrdAY7dbDniVdYgh0MCdjcOsA4xNqfN0LWjTzwyodKdULy3sa51jfFS6Bnb1Ie5lHaLvVLrT2x+ds9C5UmPXygutA/SdBu70dOGkAb0yZiazDtB3Kt0p+BAXoEFoYk6Va+WpPsTtrEP0mUp3Ok+nenFcOlaitxeMLABeYB2izzRwp6OlBTN6kGZI67pTUOlOR6VrRK+Mmdrfh7jYOkRfqXQn5EN8DKDL+4ys0uYISw8Dnm0doq80cCf3HOsAYzansxesaV13Qhq4k9Olfaa0vGDsadYB+kqlOzmVriEdYm5uz/qVSWlIf2gTqB8iPNE6x5iVpdZ0jW0J/IF1iD7SwJ3M4wGdpG/Jaewm4MnWAfpIA3cyWlowpvN0k/Ak6wB9pJE7GV3LY05vLyRApTsBDdzJaKZrbE5jNwUq3Qlo4DZUP0R7gnWOsdMpY0nYxoe4q3WIvlHpNrc7eohmThdTJkOz3YY0cJvbxTqAAE6lm4g9rAP0jQZucztbBxAoHVpeSMOO1gH6RqXb3E7WAUSnjCVEpduQSrc5zXQToEPMk6HSbUgDtzmVbgJUusnYwTpA32jgNqfSTYEepKVCM92GNHCb05puAnQFezK28iHqFcoGNHAbqI+y07dTCdDyQlI0221AA7eZ7YFF1iEESi0vpEQTkQY0cJvZyjqArKbSTYhKtwEN3GZ0A2oqNNNNyebWAfpEA7cZLS0kQpsjkrLQOkC
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADxCAYAAABoIWSWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUc0lEQVR4nO3defBeVX3H8ff5ZQGRIAgEgqUelwpUrCjKoLYjigv1ClgQwcFRQFmqbUVxOSxuIHpVLIKjKIrggkqd2hQ4MLJJp45SowECdemitwiNC3UDxSy/5/aPeyE/QhJ+93nuc793+bxmngmBGfJJcvL5nZx77jkuz3NERKQZM9YBRESGRKUrItIgla6ISINUuiIiDVLpiog0SKUrItIgla70jnPuEufcrHMun/PZ2TqXCKh0pb/uAD4E/K78/rWGWUQeoNKVznLOeefcLzaa0Y6A7YEdgf2AXwEjYHfDqCIPcHojTbrKOeeBH5ffvQa4CzgGWE8xoVgLbA04IMvz/HHNpxR5MM10pet+A5Dn+YuBL1DMbBcB64CdgFkgL/+diDmVrnTd5v6qthXwv8ACipnuro0lEtmChdYBRCa0PYBz7iqKkt0BuBfIgMUUE4ttgRtN0olsRKUrXXc3xTLCX875d9+jeIg2106NJRLZApWudN0IIM9zZx1EZD60pisi0iBtGRMRaZBmuiIiDVLpiog0SKUrItIg7V6Q1vAhzgC7AMs2+uwKPJJivC4EFh3+o29kP3jCXev/PT/SA2soXvldQ3HAzWrgzjmf1VmarG/2ZyOyaSpdaZQP0QF/AuxbfvYAdqMo16UUb5A9rK1m1183k8+u48H7czdn5EP8GRtK+C7g+8AK4NYsTdZW/XmIjEulK1NTzlyfRFGuTy+/fRqwXR3//9zNuM2+BPxgM2yYNT9zo/+21od4K/BtihJeAfwgS5NRHRlFNqbSlVr5EJ8IHAIkFG+FbTutH2vk3HxLd0sWUxTx3DK+x4f4XeCbwJXATVmaaG+l1EKlKxPxIS4AnkVRtAcDezbzI+fkuGkV4RLggPJzGrDah3g5sBy4QcsRMgmVrlTmQ1wCvJiiZF+CwbkGDhi5xjbfLANOLD+/8SFeTVHAV2Vpck9TIaQfVLoybz7EA4DjgcMoDgc35BjlMxZbHh8FHFV+1vgQrwcuBpZrh4TMh0pXtsiHuDNwLPA6il0HLZGTz8yUx92Y2Ypipv8SiiWIzwAXZmlyh2kqaTWVrmySD/FpwBspZnRbGcfZhNyNcG06WWwZcDpwqg/xSuDcLE1utI0kbaTSlQeUe2hfBrwJ+AvbNA9v1M4XKmcoHioe4kNcCZwLXJalyTrbWNIWrRy10jwf4oso9qh+lQ4ULjhyZ7KmW8XTgc8D/+1DPK7ctywDp0EwcD7EZ5YPg75G8fJCR+RN7l6Y1O7ARcCtPsTEOozY0vLCQPkQ9wDOBg63zjKWfEROm5Z052Vv4Eof4teBt2Zp8l3rQNI8le7A+BAfA7ybYkfCvM45aKuRm+lc65aeB6zwIV4GnJalyY+tA0lzVLoD4UNcBJxafoz32NYjp3tT3Tkcxc6Qw3yIHwfOzNLkV8aZpAGdWRST8fkQn0pxoMt76EnhQqdnunMtBk4GbvchzufENOk4zXR7rJzdnkaxf3SRcZxaOXAt3TI2rt2Aq3yInwLenKXJvdaBZDp6NWplgzmz23fTs8K938i16uWIuhwPrPIhPtc6iEyHZro940NcSDG7PYOelu39OrBPd1yPA77uQzwPODVLkz9YB5L69HXQDpIP8U/ZsHbb68IF8p4tL2zMUaz13uxD3M84i9So16N2SHyIhwA3UdzM0HvFukJvZ7pz7Ql804f4VusgUo8hDNre8yGeQXG+6xLjKI0a9XJJd5MWAB/0IV7sQ1xsHUYmo9LtMB/iNj7ELwNn0e09q2PJhzHTnesY4PryuE3pqKEN2t7wIe4OfAM40jqLhZxGb45okz8Hvu1D3Ns6iIxnkKO263yIzwG+w0DWbzcpx+UMs3UBT7HOq8NzOmiog7azfIjHATcAS62zWMqBWRYMbklljiXA5T7EU6yDSDUq3Q7xIb6F4ojAwT9McZD3eJ/ufM0A5/gQL9RZvd2h36iO8CG+E/iQdY7WcHkOvTh7oQ7HAxereLtBv0kd4EN8H8ULD1Ia4SBX6c7xauBzPsROH9c5BCrdlvMhfoDiOEaZI3fFfT3WOVrmaOBSFW+7adC2mA/xTOBt1jnaKIdcw3eTjgQuKS8ZlRbSqG0pH+JpwDusc7RV7oB88A/SNudVwAXWIWTTNGhbyId4MsX9ZbIZ5UxXs7nNO9GHeI51CHkolW7L+BAPBj5snaPtcufQmu7DOsWHGKxDyINp0LaID/HJwKXo9+XhOVw+7Jcj5utsvbnWLvrD3RI+xB2ByxnYSWHjGrlca7rzM0Oxo2EP6yBS0KBtgfK2h68Aj7fO0hU5Ls8ZztmOE3oUsNyHuJ11EFHptsV5wPOsQ3TJCM10K9qTYsarL1TGNGiN+RBPAl5vnaNrcgcDuTmiTi8FzrQOMXQatIZ8iAcA51vn6CLt0x3b6T7Ew61DDJkGrZHy9P/L6P8FklMxcoDWdMfhKN5Y0yHoRlS6dj7GwM/EnYRmuhPZFviiD1Ff8A1o0BrwIb4COMI6R5eVuxc0fsf3FOB06xBDpEHbMB/iUopZrkwgd7lOGZvcaT7Ep1iHGBoN2uZdAOxkHaLrckeu3QsTW0Rx+LmOgmyQBm2DfIivBA6zztEHuYNca7p12Bd4q3WIIdGgbYgPcRfgo9Y5+iInB63p1uVdPsQ9rUMMhQZtcz4B7Ggdoi/K3Qv6a3E9tgYu0h1rzdAvcgN8iIcCL7PO0Scj59A+3Vo9G/g76xBDoNKdsvIhxfutc/RNccqY00y3Xmf6EPWQd8pUutP3GmAv6xC9o90L07AEXYI6dRq0U+RD3BpdnT4VOTjNdKfi9T7Ex1iH6DOV7nT9LfBH1iH6aORcDgs0fuu3NfBO6xB9pkE7JT7E7dFf1aYmn8kBzXSn5Dgf4hOtQ/SVSnd6ArCDdYj+yt1Ia7rTshCduzs1GrRTUK6JafvNFOXO5doyNlVH+RD/zDpEH6l0p+NdwCOsQ/TZSEN32hzwXusQfaSRW7Pydd9XW+fou1xLC0042Ie4n3WIvtHArd9JwFbWIfpu5LS00JA3WQfoG5VujXyIi4G/ts4xBCNmVLrNONyHuMw6RJ+odOt1FLCLdYghyDV0m7IIONE6RJ9o5NbrDdYBhiJ3WtNt0Am6T60+Grg18SHuA+ihQ0NGOVpeaM4y4FDrEH2h0q3PCdYBhmQ0s0Cl26zXWgfoC5VuDXyI2wBHW+cYkpEmuk17kQ9R54jUQKVbjyOA7axDDMkI3RrRsBngGOsQfaDSrcfLrQMMTY7LrTMM0LHWAfpApTuhcmnhQOscQzNy2qdr4PE+xL2tQ3SdSndyL0TnLDQud7oJ2MhLrQN0nQbu5LSVxkCuN9KsJNYBuk6lO4HyymoNQgM68MbMs3yIj7YO0WUauJPZH1hqHWKItKZrZgFwkHWILlPpTuYQ6wBDlat0LWlddwIq3cmodI3kuR6kGTrIh6h90mPSwB1TeXHfXtY5hmpW5+la2gF4tnWIrlLpju+51gGGbOT0RpoxPUAek0p3fPtaBxgynadrTjPdMWnkju8Z1gGGLMdppmtrHx+ilnjGoNIdQ3mgs66nNqTjdM0tAZ5gHaKLVLrj2RtdPmlKa7qtsI91gC5S6Y5H67nGZt0CjV17T7MO0EUauOPReq6xHG0ZawGV7hhUuuPRTNeYDjFvBZXuGFS6FfkQF6OHaOZGejmiDXb1Ie5qHaJrVLrV7QEstg4xdLqCvTU0261IA7c6Xc7XAjrasTWebB2gazRwq9vNOoBAjg68aYll1gG6RgO3OpVuC+i6ntbYxTpA12jgVqfSbYFcuxfaQg/SKlLpVqfSbQHt020NzXQrUulWpzWsFhjpQVpbqHQr0sCtTjPdFsidxm5L7OhDXGgdoks0cCsob//VV/Y20IE3bTED7GwdoktUutUsBfRVvQXyXMs
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADsCAYAAADXaXXTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhnklEQVR4nO3deZwU5b3v8c8zwwz7oigKuLQi2pW4Ytgx7snNHa4nOTHBaEzMYuLJOSceJUc70RONR3RER4254BqXiBqMJuZKn6AsMsCwi7JWIyjDvsvqLN3T/dw/qjFkYJburq6nqvv3fr3mJUJP1Vcy880z1VW/R2mtEUII4Y0S0wGEEKKYSOkKIYSHpHSFEMJDUrpCCOEhKV0hhPCQlK4QQnhISlcUNKXUi0qp9UqpD9MfF5rOJIpbB9MBhPDAf2qt3zAdQgiQla4oEEqpkFIqppR6RSllK6XeUEp1MZ1LiOakdEUhOQeYqLW2gAPAz9K/P04ptVwp9ZhSqqO5eEJI6YrCsklrXZP+9SRgFPBLIAwMBo4H7jSUTQhASlcUluaDRLTWept2NAIvAEMM5BLic1K6opCcppQanv719cBcpVRfAKWUAr4OrDSUTQgAlEwZE4VAKRUCpgJLgIuB1cCNwBTgREABHwK3aK0PmUkphJSuKBDp0p2itT7XdBYhWiOXF4TwAaXUE0opWYEXAXk4QhQErXUtEMhVrlLqS8BxpnMIb8hKVwgPtPTwhlKqFHgYuMN0RuENuaYrhAfS15zXA6O01jVKqedx3uxLACVa68eUUoe01t1M5hT5J5cXhPBO84c3IkAX4DJjiYTnpHSFb4Qi0eOBM4ATgN44T5D1PuKjJ1AOlAKls+meKEGV4qwWG/fHd6+fuuX3nYB9R3xsAz4GPhk7ecpnXv73HEPzHysHA43AOuc2YroopdZprc/yPJnwjFxeEJ4KRaLlwBdw5iScDQw84p/HZ3KsOXRvUKhOh//9UGLfgujmp4e18inbcQr4yI+1wPKxk6c0ZHLuTB1xeWGE1nq+Uuo5wNZaVx3xGrm8UARkpSvyKhSJ9geGpz+GAYOATq1+UpaadCLZxktOTn+MbPb78aoxoz8AaoB5QM3YyVO25yHiGuBfj7ie+2QeziF8Tla6wlWhSPQkYDTwFWAEcEq+ztV8pbu7YcvsGdsmfdmlw6/HKeB5wLSxk6eszeVg8vCGOExWuiJnoUj0QuD/4JTtYJxHbj2XSDW6uYI4I/1xA0DVmNErgTeBP4+dPGW5i+cRRUZWuiIroUh0KM5sg2uAU01kaL7S3XBo9awFu96+zINTrwX+jFPAizw4nyggUrqi3dKXDr4H3ITzZphRzUt37YGl1Uv3TLvU4xibgNeBZ8ZOnvKRx+cWASSlK1oVikTLgArgh8DX8NElqealu3rfvDkr9s65xGCk94CncVbACYM5hI/55htI+EsoEu0O/AS4DehvOE67xJMNpYYjXJ7+2F41ZvRE4Mmxk6fsNpxJ+IysdMU/CEWifYBbcfYX62U2Teuar3QX7oouqj200k87QzTgPHn26NjJU2zTYYQ/SOkKAEKR6JnAL4AfkKf7aN3WvHTn7Hhz2da6dReYzNSCFPAy8Ouxk6dsNB1GmCWlW+TSK9vfAD8mYJebmpfuzK2v2LsaN1smM7WhEXgKGDd28pRdpsMIM6R0i1QoEu0E3I6zO24Pw3Gy0rx0p25+fv3+xK4zTGZqp4PAo0DV2MlTDpoOI7wlpVtkQpGowrnh/wEM3V/rlual+/bGidvqkgf7msyUoV3AOGCi3O1QPGSIeREJRaJDgEU41xcDXbjHktDxrqYzZOhE4HFgadWY0YMNZxEekZVuEQhFol1xVrb/RgH9H23zle7r68cnNdr0bWPZSuIU8H+NnTyl3nAWkUcF8w0oji0UiV4NrAJ+TgH/7621bghw4YIzI3gssLxqzGivn6oTHpKVboEKRaLdcPbeusV0lnw5cqWrtd79eu34E0xncokGngHuGDt5ygHTYYS7CnblU8xCkegwYDkFXLjNaVKF9CO5An4KrKoaM/pqz06q1O+VUsuUUsvTG2fKQPU8kNItMKFI9OfAbJyxhEUjpVN53fnBkFOAqVVjRt9dNWa0F+Myb9NaX6C1Ph/YiPMegHCZlG6BCEWi3UOR6GTgt0CZ6TxeS+lko+kMeVIC/DfwVtWY0T3dOGBL28FrrQ+k/1wBnTl6TzfhAindAhCKRL8ILAa+bTqLKUndVKile9g1wOKqMaPd2nniHGCi1toCDuDM2kAp9QLOXnJh4HcunUscQUo34EKR6PXAQpxvoqKV1IlieLhgILCwaszo77hwrObbwY8C0Fr/AOgH2MAYF84jmpHSDbBQJHoP8AoQtIcCXNeOTSkLRRfg1aoxox+vGjM6l1kZzS8dfP7vWusk8EfgmzkcX7RASjeAQpFoaSgSfQa413QWv0ik4sVSuofdCrxdNWZ0lyw//zSl1PD0r68H5iqlzoLPr+leA8Ryjymak9INmFAk2gV4C7jZcBRfSaQaU6YzGPC/gJlVY0b3zuJzD28HbwPH4WwH/5JSagWwAugL3OdaUvG5QI3yK3ahSPREYArgp0HdvpBINRjZgdgHhgJzH/3O9Vff/tqrmzP4vCat9Xeb/d5IF3OJFshKNyBCkehpQA1SuMcUTzUWa+mC6lxX3uOmmRNumXmW6SiibVK6ARCKRPvjbHo40HQWv4onG4rza1mVr+zY80dnK9VhIDBnwi0zw219ita6Vmvt1q1nIkPF+YUaIKFI9GRgJnCm6Sx+Fk/VB3nYTZY6rO3Y40enKFV++HHdk4H3Jtwys6hvH/Q7KV0fC0WiJwDTgbNNZ/G7eKqhyN6fKNnQsecPe6qSzr2a/cHJwMwJt8wcYCCUaAcpXZ8KRaLH4xTuF01nCYJ4sr7cdAbvqG3lPW4qVSXd+rTwgn7A9Am3zOzvZSrRPlK6PhSKRLsD7wB+3NnWl+Kpxo6mM3hkT3n379aXlPY6pdVXaX1qqPZvj9lhq5c3sUR7Sen6TCgSLQFeA75kOkuQJFINxVC6B8q7j9lZ0uHE1q/va9149trJi8+snfIt4C07bBXD301gSOn6zyNAhekQQRNPNXY2nSHP6su6fX19SYf+rW8xr/XB81Y9s/qUrXOGpX/nUuAlO2wV7y11PlNkbz74WygSvRm4zXSOIGpKNRby/IlEWZevrSotO7P1n350ategDx7b3evAJxc1+5MxwBac7YCEYbLS9YlQJHo5MMF0jqAK4E7A7ZXq0PmyJaUdrVYLV6WSm4Yuvr+u14FPWloJ326HrevzkE9kSErXB0KR6EDgDYpw+LgbtNb1FOjXcmmnoTUdOg0a3tprSpLxtcMX/Lq8a92O09s43NN22JLbDw0ryC/UIAlFop1xBtgcbzhKgOlDphPkQ2n5edVlnUde0tprOiQ+Wz5y/q/6dIrvO6kdh+wGvG6HrU5tvlLkjZSueY8BXzAdIsg0upA2pQSgpGxAdVnXq1vdir1jw6eLRs6/a2BZU30m2/hcADyeUziREyldg0KR6Ddwdn0VOSi0TSlVh/7V5d3+qdXC7XZo89zhC+8ZVJpKZHPXxk/tsHVdlvFEjqR0DQlFoqcCvzedoxAU0qaUqqR3TXm3b3+5tdcc/+nq6sFLHhxZolO53H30jB22ZICSAVK6BqQfgJiEMzxa5KhQNqVUJT0Wlve4cVh654Zj6ru1pvrC5RMuVZDrfbfdca7vyoMTHpPSNeMuoNXVjGi/Jp1oMp0hZ6rz0vIeN12oVMmxp6VpnTxj/ZS51kevtnrZIUMXAlUuHk+0g5Sux0KR6LnAr03nKCTJVMA3pfz7TNxjrzq1rg+vmfT+GRv+NioPZ/+ZHbaG5uG4ogVSuh4KRaIKeAp5EtBVCd0Y4NI9aibuP9J6//krnvyo3/YF+doxRAET7bBVhPOIzZDS9daPkH2oXBfcTSlbnInr0KkdX1r68I4TPl2V72lzg4B/yfM5RJqUrkfSA8kfMp2jEAVzf7TWZ+KqVNOGYYvui/c4uMGrJ8jut8NWex6wEDmS0vXOI8hTZ3m
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADsCAYAAADXaXXTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYxklEQVR4nO3de7RcdX338fcvN4GonCCXJo4IBQoqRYoCchlIRwk8tVZdtV7w0oeiMHWepz4qPh20Xb15mbpa8fIMjlUUW2ltxdYL2GpgN2FIQEBAQJJgJGAGCJcAOefk5Fxm5vv8sSfhGBJyZmbP/u09+/Naa9aZczKz57OyzvnMb3577992ZoaIiMRjnu8AIiJZotIVEYmRSldEJEYqXRGRGKl0RURipNIVEYmRSlcSzTn3qHPOOrcHnHMH+M4k0g+VriTdU8CfAQa0ga/4jSPSH5WuJIJz7i2d0Wyr87XpnHspcDDwt52HrQRO9ZdSpH8qXUmKZZ2vl5mZA1rAWuBpM2t2/i0POB/hRKKi0pWkOKHz9XXOuTuBCeCwWf/ugNs7PxdJLZWuJIUDMLMTzexE4F+AGWDEOfdXncdcBTzkJ55INFS6khR3ATjnPtf5/p3AL4EHgLcT7kh7G/BdH+FEouK0ypgkgXPuLcC3CI9QmEdYsm8EvrfbQ7eY2dKY44lEZoHvACKzmJnN3+1n2nEmQ0XTCyIiMdL0gohIjDTSFRGJkUpXRCRGKl0RkRjp6AVJjGoxGAGOBI7ofH0pMAIsmn2bpPXka0cWcTDzDiM8XbgJjANPAFs7t9n3HwE25yp57cAQ77QjTWJTLQaLgGP51WI9Ytb9A+eynSbtdeeNLGovxr2ii5ffAWwENnRu93W+rs9V8tu62I5IX1S6MjCdkj0F+G1gOXAasH+/223SXn/eyKJWl6W7NwasJ1xcZw2wNlfJb4hguyJ7pNKVyFSLwULgZJ4p2dOByBcdb9LecN7IomZEpbsnTwA3ATcC1+Yq+Z8N6HUkg1S60pdqMTgZeB1hyZ4BLB70azZp33feyKKZAZbu7n5BuObDd4E1uUq+FdPryhBS6UrXqsXgMOA9wB8Bx8X9+h5Kd7atwLXAfwA/yFXy0x4ySIqpdGVOqsVgPvA7hEX7u3g88qVF++fnjiya9lS6sz0O/BPw5Vwlv95zFkkJla48p2oxOIawaP8QSMTqXi3aG88dWTSVgNKdbQ3h9dv+LVfJa6F12SuVrjxLtRgcALwFuBA4y3OcZ0lo6e60jXD0+3e5Sv5B32EkeVS6sku1GBwIfBD4AOFJCYnUov2Lc0cWTSa0dHdqEl7p4lM6BE1mU+kK1WLwAsKi/RCwxHOcfWph9587snBHwkt3pzbwbeATuUr+p77DiH8q3QzrTCP8CXAJ8CLPceYsZaU727XAn+Uq+Tt9BxF/VLoZVC0G8wgP+fo48GLPcbrWwjadO7JwIoWlC+HI96vAx3KV/GO+w0j8VLoZUy0GBeDvgRM9R+lZG3tgxcjC7Skt3Z1Ggb8BPper5Gd8h5H4qHQzoloMjgY+A7zBd5Z+tbEHV4wsHE956e50H/DhXCV/je8gEg+VbgZUi8HFhKPbgZ+iG4e22YMrlgxN6e70A+DiXCXf8B1EBkulO8SqxeBQ4ArCM8iGRhvbvGJk4eiQlS6EUw6X5Cr5L/sOIoOj0h1S1WLwBsIzpA71nSVqQ1y6O/0IuCBXyT/sO4hET6U7ZKrFYDHh3O1FvrMMShtrrBhZuG2ISxfgSaCYq+S/5TuIREulO0SqxeAU4BvAMb6zDFLb7KEVSxY+PeSlu9PXgT/OVfI7fAeRaKh0h0BnBbCPAX9OBq571zZ7eMWShU9lpHQB7gDerLUchoOuBpxy1WKQA+rAX5GBwu1wvgPE7LeA2xrl+mt9B5H+qXRTrFoMXk54ba/TfGeJU9Yat+Ng4IeNcv3DvoNIfzS9kFLVYnAG8H1SsEBN1MxsyzlLFm7N0PTC7v4FeK/W7U0njXRTqFoM3ghcRwYLtyOjg91d3gEEjXL9IN9BpHsq3ZSpFoOLCJcK3M93Fo+yXroApwI3NMr1Zb6DSHdUuilSLQZ/AXwJmO87i2cq3dArgDWNcv1o30Fk7jSnmwKdpRgvBy72nSUJzOzxc5YsfCzDc7q7exQ4V4ukp4NGuglXLQb7AVejwp1Nv7e/6jBgVaNcP8N3ENk3/fImWOeU3h8Bb/adJWE0vfBsI4SHlL3GdxB5birdhOqcZfZNIO87SwKpdPdsMXBto1w/3ncQ2TuVbnJ9niFbkjFCKt29Owj4UaNcP9J3ENkzlW4CVYvBJcD7feeQ1FoKrGyU67/mO4g8m0o3YarF4A+AT/vOkXAa6e7bUYQj3hHfQeRXqXQTpFoMTgf+EZXKvuj/Z25+E/h+o1xf5DuIPEOlmxDVYnAM8D2yfabZXKl05+5Mwv0DkhAq3QSoFoODgf8EXuQ7S0qodLtzcaNcf5/vEBJS6XrWOfnhe4RzcDI3Kt3u/T8dw5sMKl3/vk7G1sONgEq3e4uAbzfK9aW+g2SdStejajF4H/BW3zlSSKXbm2XA1Y1yfaHvIFmm0vWks+PsMt85Ukql27vTgU/6DpFlKl0PqsVgAXAV4Wmb0iWn39t+fahRrp/pO0RW6ZfXj78ETvYdIsU00u3PPODKRrmuN30PVLox65wAcanvHCmn0u3fUcDf+Q6RRSrdGFWLwfOAK9D/e79UutEoNsr1Fb5DZI3++OP158BxvkMMAYcueBKVr2p9hnipdGNSLQavBP6v7xxDwqHWjcqLgb/1HSJLVLox6CxIfgWg4yOj4JzT/EKk3tso10/0HSIrVLrx+N/Aq3yHENmLeWhRnNiodAesWgyeD3zMdw6Rfcg3yvW3+Q6RBSrdwfsT4GDfIUTm4NONcn1/3yGGnUp3gKrF4EDgEt85RObocLSzd+BUuoP1QWCJ7xAiXfjTRrm+zHeIYabSHZBqMTiIsHRF0mR/4CO+Qwwzle7gXAK80HeIYaUDdQfqoka5fojvEMNKpTsA1WJwCOEONJE0OoCIP6U55650zm1yzt3ZuZ0Y5fbTRKU7GGW0bKOkW2kApwd/xMxO7NzujHjbqaHSjVi1GCwF/th3DpE+vZAePq05545wzq13zl3lnFvnnLvaOXfAAPKllko3eh8l3BkhknYfaJTrz+/heccCl5vZy4BR4P2dn3/COXeXc+4y59zzIkuZMirdCHWOy32v7xwiETkI+KMenrfZzNZ07n8DOJNwDenjCBfvPwj400gSppBKN1pvBfbzHSILtOBNbC7u4Tm7H1liZvaIhaaArwGn9B8tnVS60XqP7wAiEXt5o1w/q8vnHO6cO61z/3zgRufcUti1QNybgHuii5guKt2IVIvBrwNn+M4hMgDFLh+/ASg559YRnpH5ReAq59zdwN2Ea5F8PNqI6bHAd4Ah8m70qVeG05sb5fpIrpJ/eo6Pb5rZu3b7WSHiTKmlkW503u07QLbo/S1G+wHv8B1iWKh0I1AtBmcQXl1VZFhdMJcHmdkDZnb8oMOkmUo3GtqBJsPu5Ea5/uu+QwwDlW6fOpdVf6vvHCIxeKPvAMNApdu/3wNGfIcQiYFKNwIq3f5pasED7Ubz4sxGuf4i3yHSTqXbh2oxeAFwnu8cIjGZD7zed4i0U+n25zR0rLNki6YY+qTS7c+ZvgOIxOzcRrmu9UX6oNLtj0pXsmYxcLrvEGmm0u1RtRgsAE71nUPEA60x0geVbu9OIryWlEjWaKTbB5Vu7zS1IFl1WqNc11F7PVLp9i7vO0CW6S/eqwOBV/gOkVYq3d5pXkuyTL//PVLp9qBaDI4FDvGdQ8QjlW6PVLq90XyuZN2rfQdIK5Vub/QuL1l3VKNcn+87RBqpdHujRZol6xYBR/oOkUYq3d5oMWcR+A3fAdJIpdulzspiWt5OBI71HSCNVLrd00cqkZBKtwcq3e5pakEkpNLtgUq3exrpioSO9h0gjVS63VvmO4DoNOCEONh3gDRS6XbvMN8BRBJiv0a5vth3iLRR6XZPpSvyDI12u6TS7Z5KV+QZOnyySyrd7h3qO4A
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADsCAYAAADXaXXTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgV0lEQVR4nO3deZhU1Z3/8fep7uqNplmUfelCUOoKCCqbuzEaZ1LEGGOCSZQsY34Tk0zm+YUklpN1EuNUFrIYo2bMxGhiYhMzidqFxJ1FkUVW4RaIUogogkADsjR015k/bqmIgF1Vt+rc5ft6nn7obrrrftDqT586995zlNYaIYQQlRExHUAIIcJESlcIISpISlcIISpISlcIISpISlcIISpISlcIISpISlcEklLq90qpDUqp5fm3caYzCQFQbTqAEGX0da31faZDCHE4GekKX1NKxZRSGaXUPUopWyl1n1KqwXQuIY5FSlcEwUjgVq21BewGvpj//A+VUiuVUj9XStWaiyfE26R0RRBs0lo/lX//j8C5wA1AHJgA9AauN5RNiHeQ0hVBcOQCIlpr/ap2tAN3AhMN5BLiXaR0RRAMVUqdlX//k8B8pdQAAKWUAi4HnjOU7S1KqZuVUm+YziHMktIVQbAW+JJSygZ6AbcB9yilVgGrgBOBGw3mQyk1Pp9NhJySpR2FnymlYkCr1nq0R7LMBp4FzgBWA9OAduBRnFH481rrRlMZhXlyna4Q7hoJ/IvW+iml1O9wrqQ4BDygtX7Vme0QYSYjXSFckh/pztVaD81/fBGQBBqAC7XWHUqpN2SkG24y0hXCXUeOYibgTC+sz49yG5RS67XWIyqeTHiCnEgTwl1HXklxo9a6v9Y6prWOAfukcMNNSlcIdx3tSgoh3iLTC8ITYsm0AvoAg4DBR/zZC4jiPF+jQPVddHvpYNvigyt3zmkGDubf9gHbgC1HeXttektrZwX+KR1a66uP9ZcynyukdEVFxZLpgcAknLnOk3i7WAcCNV19nHrUwkMq0g6c38Vvyc2YOmUrsAZYiXP97kpg9fSW1v1d/xcIURq5ekGUTSyZrgfOBCbjFO1knJIt2UwaF+5qW9S+cuecrpbuseSAF3AKeBnwJLBwektrR4mPK8RRSekK18SS6V5AAjgbp2DHUKZXUzNpXNjW9szBVTvnnVeGh98DzMW5oeGx6S2tq8pwDBFSUrqiJLFk+gSctQ0+BlyEM+dadjNpXLhj59MHV7c9VY7SPdJrwOPAQ8D901tad1fgmCKgpHRFwWLJ9InAR4ArcYq24ucGZtK4cPvO+YfWtC04t8KHbse51bcFeGB6S+veCh9f+JyUruiSfNFegTOivRDDJ2Fn0rhw2465HZldC88xGOMN4H+Bu4Anpre0yg+TeE9SuuK4Ysn0BODfgY9ToamDrphJ48Kt25/sXLt78dmms+S9BNwO3D69pXWn6TDCu6R0xbvEkulq4KM4ZXvWe3y5ETNpXLhl++O553c/67V8e3EWTf/59JbWF02HEd4jpSvekr/E61rga8BQw3GOayaNC199/VG9fs+yyaazHEMO+DswY3pL69OGswgPkdIVxJLpHsCXcEa2fQ3H6ZKZNC7c/PrDvLhnxSTTWbrgGeCm6S2tD5o4uFLqf4DxgALWAZ/RWssOFoZI6YZYLJmuBb6Ks2ljD8NxCjKTxoUvb5utNryxyk97n80Fvja9pXVxJQ+qlGrSWu/Ov/8zYKvWOlXJDOJtsuBNSMWS6Sk4OxvchM8K9005cn5bEfx8YOGMqVP+PGPqlJjbD66UiimlMkqpe5RStlLqPqVUw2GFq4B63r38pKggKd2QiSXTI2LJdBp4EBhuOk8ptPZd6YLzEv8qIDNj6pSfzpg6xe1900YCt2qtLWA3zs4VKKXuxFn4Jw78yuVjigJI6YZELJnuFkumb8LZFfeDpvO4QZPz8/O3FpgOrJ8xdcqXZkyd4tYvkE1a66fy7/8ROBdAa/1ZnEWFbGCqS8cSRfDzk1Z0USyZvgrIADfg/LAHQk5rP450j9QbuAV4csbUKW688jhy6uCtj7XWncC9OJcDCkNkaccAiyXTI4HfABeYzlIO2n9zusdzPrByxtQp/wHcXMLdbUOVUmdprRfg7FwxXyk1Qmu9Pj+nexnOL2BhiIx0AyqWTE/D2Qo8kIULkNO+nl44mgbgF8CcGVOnFLulz9F2rrhLKbUKZw3hAcD3XcgqiiQj3YCJJdMNwK+BzxiOUnY+PZHWFefhjHq/CfxyektrroDvPdrOFSbXpxBHCNpIIdRiyfRoYAkhKFzw/Ym091IP/AyYNWPqlBNMhxHuCfKTNlRiyfS1wCLAMp2lUgI4vXA0lwLLZkyd8p533mmts1rr0RXIJEoQhidtoMWS6cZYMn0PcAfO6Cg0Aj7SPdwQ4LFfXfvna00HEaULy5M2kGLJ9FhgKc5Z6tAJyUgXgKraM5+NVPe749dfePx/fv2Fx7u8gafwntA8aYMmlkxfDMwHTjadxZSwjHRVVZ/50YYLnA04tb56zHO/udOOWzLP61OheNIGTSyZvhJIA42ms5ikQzHSjdo13T8x/s2PhmVnLerz+spPAs/YcSu0v3D9LARP2mDJnzBrAUL/ElOjg/783V7TdE13parrAHrvWPPksI2z3twTbgRO8fphaUtxmKA/aQMllkxfj3PCTP6/EejrdAE6o92mbIxU9RwMUHtgx6KxK289/4iv6Q08YsetIz8vPEx+eH0ilkz/GJA1UA+TC/BIt6pm1PyqmlPOAIh0tq+dtPjGU9XR/73dgYfsuPWByiYUxZI70jwulkxXAf8NfM50Fq/R5KpMZygHFemxINrtUuf2bZ3bNmnxDxurO9uPN3/fADxgx62PWRnbyO4UousCO1IIgvzODn9BCveognkireqFmqZrxgCg9YFxK255rf7A9kFd+MZa4K923PpYWeOJkgXwSRsM+RFuC/AR01m8KoAn0nbXNH0qolRNI8DwF//+bO+2tYXcYRYF/mzHrQ+VJ55wQ9CetEFyG/Bh0yG8LKcDNb2gqxsusSNVJw4DOHHb8iebNz1azEI1VUCLHbdkkRuPktL1oFgy/X3g86ZzeF2Q5nQj0eFzq2vHTAKo379twZjVd5SyJGc98KAdt2QdBg+S0vWYWDJ9HfBt0zn8QOuATC+obkui3S47D6Cq44A9cfFNY5Wzl1opegGz7bg1tPSAwk3BeNIGRCyZ/iCyaWCXBWOkG3mptmnaCKVUBJ3bMmnxD3pV5Q42uPTgg4B/yC3D3iKl6xH5tXDvxZmTE12gtfb7f6t9Nd2nHlCR+p5ove+M5b/YWdfe1t/lY8SB++y4JZeHeoSUrgfEkum+QCvOhe6iizT+Lt3q+vOWRaoHnILW+pT1M1f03PVCudZCvhD4cZkeWxRIStew/LW4fweaDUfxHT+vMhapHjynum7COQD9ti6ZM3jz3LPKfMj/b8etT5T5GKILfPukDZAfA+X+gQsijV+fv6puRbTxynMAGva++tQo+/cXVujIv7Xj1pgKHUscgz+ftAERS6YvAf7NdA4/0loXu0W5YerV2qZpA5WKVFd17H9uwrOpMyt48Abgb3bc6lXBY4ojSOkaEkumewO/p/RLg0Iphy9L92C08YrtKtLYR+U6N09e9P1+VbmOugpnGA78rsLHFIeR0jXndmCg6RB+paGQbck9oap2wjNV0ebRaP3Gmctm7K09uLuPoSiX23Hr04aOHXpSugbEkulpgCxMUgLts5Guquo7L9pw3vlonYuv/eOapj0bTzEc6Zdy44QZUroVFkumY8gNECXz15xudE1N96smAgx49el5A7c8M9F0IqAHcKcdt2R6q8KkdCsolkxHgLuBJtNZ/M5Hc7qv1zRN66lUdW3jnk3zrHV/KmVNBbddBHzFdIiwkdKtrK8D55kOEQQ+mV7ojHa7bFOkqsfA6kN7V4xf+hMv7mf2X3bcMj3VESpSuhWSv833+6ZzBIUfSreqZsz8qpoRp6tcx0uTF/3nkIju9OJmovXALaZDhImUbuX8FNnB1zVeL10V6bkg2u2SC9B69/h
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADsCAYAAADXaXXTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZSUlEQVR4nO3debxkZX3n8c9z6R26m00EqbCD4rCZAMpSLBWJ8YXjyMRJFI1xGc0hpdFIiOWYyYuELDUYxYlTWBEEjaJxMCMqLmhzBKppkEWgQZsWZC1BNqXv7Xt7uVX1zB+nGpqWhntqOb+zfN+v133d2923T31f3ae//dRzznke571HRESSMWEdQESkSFS6IiIJUumKiCRIpSsikiCVrohIglS6IiIJUulKajnnVjvnfP+j55xbYp1JZFgqXUmzDvAJoN3/8cWGWURGQqUr5pxzb+6PZDduNar9IrAPUAO6RAX8atOgIiOg0pW0cMDDwAJgCjgd6HrvO/1f36H/PSKZptKVNNkRuBnYDOwMzO///G5AD5ixiSUyOipdSZMveO+PAi4HNgE959zfEo1yZ4FfGGYTGQmVrqTJnzrn9gDOBB4AHgHeAjxBNOr9hl00kdGYZx1ApM8DS4DH+l8/ABzBc+dxP+Wc28N7f27i6URGRKUrqeG9X2SdQWTcNL0gIpIgp0XMRUSSo5GuiEiCVLoiIglS6YqIJEh3L0hqNIJwd+CA/sf+wMuAxcBCYNGWjy8vWf/khQuWu92Y2INn12VYDzwFPLnN56eAR4GHS/WyLmCIOZWuJKYRhAuIyvSArT5vXbLL5nKcebi7FuEmgFfGePkN7VrrXmDtth+lenldjOOIDEWlK2PTCMKFwPFApf9xDM+upzCwCe97xJ8aWwwc3v/Ymm/XWmuB67d8lOrlnw2bUWR7dMuYjEwjCOcRFeuWkj2eaEpgpP5j8fQdn1m4fP6OuDgj3TieBFYBK4HvlerlO8f0OlJAKl0ZWCMIHXAkUcH+LlAGlo77db++ePq2xsLlC8dYutu6j2jdhyuIRsLdhF5XckilK7E1gvBlwLuA9xDNxSbqG4umf/zpRcsXJVi6W3sS+DZRAX+nVC9vNsggGabSlTnpTx2cDvx34PVEyy2a+Obi6Vv/ZeHyxUalu7UngX8DLirVy3cbZ5GMUOnKC2oE4QFERftOYC/bNJErF03f/KlFy3dMQelubSXRHm7/t1Qvb7AOI+ml0pXf0L/r4L8Sle2ppGybnG8vmr75gvSV7hbrgC8CnyzVy/dbh5H0UenKMxpBuAtwNnAWsKtxnO367sKZmz6xeNlOKS3dLTrAl4F/LNXLa63DSHqodIVGEC4DPgR8GFhum+bFXbVw5saPL162LOWlu0UP+BrwD6V6ebV1GLGn0i2wRhDuCHwAOIcUj2y39f2FMzeen53S3cIDVwLnlurlH1uHETsq3QJqBOEE0S1ffw/saRwnthULZ1bVFy/bOWOlu0UP+DzwP0r18mPGWcSAVhkrmEYQngrcSnSlPXOFC+DSdV0vrgng3cDP2rXWX7ZrraEfi5ZsUekWRCMID24E4TeAEDjKOM5QJqK36lm3DPg4cFe71nqDdRhJjha8ybn+o7ofBOpESyRm3oTPRelucQjwrXat9V0gKNXLD1kHkvHSSDfHGkG4F/A94AJyUrgAEy5XpbvF64lGve+zDiLjpdLNqUYQngHcCfyedZZRm/CZntN9IUuBf23XWt9v11p7W4eR8VDp5kwjCHdsBOHFwP8DdrPOMw4T+J51hjE7DVjdrrXebB1ERk+lmyONIDwWuJ1o9a/cyvjdC3O1K3B5u9b6fLvWWmIdRkZHpZsDjSDcoRGEf02088FB1nnGbSKPM7rb9yfAqnattZ91EBkNlW7GNYJwX+Aa4DwKcjeKy8ctY3EcCdzSrrV+1zqIDE+lm2GNIDwKuBE40ThKoiaKMb2wrd2Aq9q11tnWQWQ4Kt2MagThKcC1ZPSpsmHk5OGIQewA/HO71vpyu9ZabB1GBqPSzaBGEL6Z6P7bOW1Znjc6aXkr0GrXWi+xDiLx6fzNmEYQngV8lRw97BCXK+o497l+B7iuXWuVrINIPCrdDGkE4d8BF1Lwv7cJV8g53efzCmBlu9Y62DqIzF0hrnZnXSMIdyAqWz0iika629iXaKrhdaV6+Q7rMPLiCj1iyoJGEC4CLkeF+4wCX0jbnpcC17RrreOtg8iLU+mmWH8bnauAM6yzpMlEyjbKTImdgR+0a61C3T6YRSrdlGoE4XzgP4CTrLOkTQEfjpirJUTLRB5hHUS2T6WbXk3gtdYh0mjCa6T7AnYmeojiAOsg8vxUuinUCMKPEW3pIs+jIAveDGNPoqmGwj04kwUq3ZRpBOGZROsoyHbopJ2TA4hGvDtbB5Hn0vmbIo0gLAOXogtFL0h/OHN2BHBlu9ZaZB1EnqXSTYlGEB4CXAEsMI6SejppYzmB6B5vSQmdvynQCMKXAN8hWrhaXoRO2tje1a61zrIOIRGdv8b6Dz98EzjQOktWOK/ngAfwv/XwRDqodO19EXiNdYgsUeMOZD7wNd3RYE+la6gRhO8HtPlgTE69O6i9iIp3vnWQIlPpGmkE4aHA+dY5skiPAQ/lBHTemVLpGug/4nsZoNX/B6DGHdoH27XWqdYhikqla+M84FXWIbLKoQtpQ3LAJe1aa6l1kCJS6Sas/wDEOdY5skxzuiOxH/AJ6xBFpNJNUP/2sM+hP/ehqHFH5r3tWuv3rUMUjf7xJ+tcQFurDEmrjI3UxVqfIVkq3YQ0gvBVwNnWOfJA0wsjtTfwSesQRaLSTUAjCOcRTStoT7oRUOmO3Dvbtdax1iGKQqWbjA+huxVGRqU7cg74lHWIolDpjlkjCJcDH7POkSe6ZWwsjmvXWm+1DlEEKt3x+zDRFioyIhrpjs3/atdaS6xD5J1Kd4waQbgr0dSCjJBKd2x+C91DPnYq3fE6B1hmHSKHdN6Oz1+1a62SdYg808k7Jv2FyT9gnSOPNNIdqyXAR61D5JlKd3w+AuxoHSKPnB6OGLd3a93d8VHpjkEjCPcE/sw6R145nbfjtojoAvDIOOc+75y73zl3e//jqFEeP0t08o7HR9GyjeOkke74Be1aa5cRH/Mc7/1R/Y/bR3zszFDpjlgjCPcG/tQ6R55pTjcRSxngmoRzbj/n3N3Oucucc2ucc19zzuk2tK2odEfvY8BC6xB55nA6b5PxwXattdMAv+/lwIXe+0OBSZ6davsH59xq59wFzrnC/hvRyTtCjSDcHXiPdY4C0HmbjF2B9w7w+x723l/f//pLwIlEU26vAI7pH/cjI0mYQTp5R+utwALrEHmnC2mJOqtda8WdzvHb/th7/6iPbAIuBQq7wI5O3tH6Y+sARaA53UQdDMTdT20f59xx/a/PBFY65/YCcM454E3AXSNLmDEq3RFpBOGWt04yfjpvkxX3wvBaoOqcWwPsAnwGuMw5dydwJ7A78PejjZgdWt91dN5hHaBAVLrJelO71tqtVC8/Ncfv73jv377Nz1VGHSqrdPKOQCMIHfA26xxFoemFxC0gul4hI6DSHY1TgH2sQxSIztvk/clcvsl7/4D3/rBxh8kynbyjoQtoCdLdCyaObtdah1qHyAOdvENqBOFi4M3WOQpG562NN1kHyAOdvMM7g+iRSUmOzlsb/9k6QB7o5B2ephYS5mAH6wwF9ep2rbWHdYisU+kOob8dz2nWOQpI562NCeB06xBZp5N3OCegUZcFnbd23mgdIOt08g7nROsABaXz1s5p7VqrsCuEjYJO3uGodG3o3YWdHYm/FoNsRaU7oEYQLgKOts5RUDpvbWmwMQSdvIM7Bi3jaMNrpGvseOsAWabSHZz+t7ej0rV1bLvW0t/BgFS6g1PpGtFjwOZ2BA63DpFVOnkH0F9VTG+xrDinUZa94178W+T5qHQHcxiws3WIInP+N7aEkWSpdAek0h2Mphak6FS6A1LpDuYE6wBFp1XMzR3YrrUWWYfIIpXuYI6wDiBizAEHWofIIpXuYPa3DiCSAgdbB8gilW5MjSDcHdjJOkfRaXo
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "[o.plot_piechart() for o in beta_p];"
   ]
  },
  {
   "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.6.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}