pyerrors/examples/04_fit_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyerrors as pe\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('./base_style.mplstyle')\n",
    "plt.rc('text', usetex=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read data from the pcac example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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": 4,
   "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": 5,
   "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": "iVBORw0KGgoAAAANSUhEUgAAAt0AAAHECAYAAADlBpY8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAABSLElEQVR4nO3deXzdVZ3/8ddJui9JuligG226sAgIXUCUGcC2uAsDXVBUFKV1RxSJzMwPl1lqcRtwRFodUUeU0iIi6gAtiIAK0hYoO21TaJsCpUuatnTP+f3xvUnTLO1NmuRur+fj8X0k93u+328+tze9eefkfM8JMUYkSZIkdZyiTBcgSZIk5TtDtyRJktTBDN2SJElSBzN0S5IkSR3M0C1JkiR1MEO3JEmS1MEM3ZIkSVIHM3RLkiRJHaxLpgvIpBBCAAYD2zJdiyRJknJWX2B9PMSqkwUdukkC97pMFyFJkqScNxSoaqmx0EP3NoC1a9dSUlKS6VokSZKUY2pqahg2bBgcZuRE1oTuEEIZMD31cBRQBlTEGKsPc97MBg/LYozXtfZrl5SUGLolSZLUYbLpRso5wJIY47wYY0Vq34JDnZAK3GWpc+YBlSGEOR1dqCRJktQa2RS6y4HJDR6vavS4ORXAwroHMcaFwMyWD5ckSZI6X9YML4kxTmm0axSwuKXjU8NRymOMlY2aykII42KMy9q5REmSJKlNsqmnu14Ioa7Xe9YhDitvYX/1IdokSZKkTpc1Pd11UuO0ZwGzmunFbqh/C/s3t9QWQugOdG+wq2+bipQkSZJaIet6ulM3RY4HKkIIV7fz5a8BtjbYnKNbkiRJHS7rQncDc4A5qaEmzdncwv7+h2ibDZQ22IYeUYWSJElSGrIidIcQykIIC1I3R9apG1rS0gwmlXXnNtpf1uDcg8QYd8cYa+o2XP5dkiRJnSArQjcHpgtsOBa7LPWxpQBdnWprMn7bmUskSZKUTbIidKdC8rxGN07OAJbFGBdDMqNJo9UnIRmCMrXuQaq9AkmSJCmLZNPsJbMbrSZZBkxq8HgySaCeV7cjxjgvhHB13cqUwIAGq1mm5Z//oVubC5YkSZLSEWKMma4hY0IIJfFrJVt3nvtNep59RabLkSRJUo6pqamhtLQUoDR1z2CzsmJ4SSZ956+76fmna+GJX2e6FEmSpIJSWVnJrFmzCCHQr18/Kioq6rdp06YRQqCi4sAghrr9uSibhpdkxFcW7eYLMz9Gtzs/A937wAnvz3RJkiRJBaG8vJy5c+eyZMkSJkyYwJw5cw5qr66u5vLLL69/PHHiRAYMGNDkOvPmzWPmzMa3/mWXgu/pBtg1eTaceAEsvAxW3Z/pciRJkgpK//7NLzReVlbGxIkT6x9PnTqVq69uunbiokWLOqy29mLoBigqhn+aC+XnwK2XwJpHM12RJElSwaqsrKS6uhqAyZNbWrIl6QmvqKigsrLZGaaziqG7TpduMP0XMHgc3DINXlme6YokSZIK0uLFi9m8OVlgfNy4cQAsW7aM8ePHM378+IOOq6yspLKysn4seF1YzzYFP3sJsHXr1q2UlJQkO3fVwC8+ANVr4bK7YeCYjNYoSZKU76ZMmUJlZSVTp06lsrKShQsXsmrVKsrLyw86bvHixcyaNYtVq1YdtK+iooKlS5d2dtlA+rOXFPyNlE30KIEP/wZufjf84vwkeJcNz3RVkiRJTe15Aza+mOkqEgPHQrdebT598uTJ9TdSXnfddc0e0ziE5xJDd3N69YeP/BZuflcqeN8DfQZluipJkqSDbXwR5p2d6SoSM/8Mg09tl0sdahx3rjJ0t6TkGPjonfDTd8EvLoCP/T4J45IkSdli4Ngk7GaDgWPb7VJ147jborKyMit7xA3dh9JvRBK8b34P/O8F8NHfQc+yDBclSZKU0q1Xu/Uu54tly5ZlZeh29pLDedNxSfCuXgO/vDC50VKSJEntqm62kkNpbmaS8vLy+ikDs7WXGwzd6Tn6pGSM98aVyXSCu7dnuiJJkqScV7cM/OLFi+tnJlm8eHGzxy5btqx+Tu5Zs2bV7y8vL2fmzJlUVFSwePHiIxqa0pGcMrDxlIGHsm5JMr578KnwoduO6A5dSZIk5b50pwy0p7s1hk6ASxZA1TK49UOwd1emK5IkSVIOMHS31rFnwoduhTV/g9s+Cvv2ZLoiSZIkZTlDd1uM/Ee4+FdQ+SdY+HHYvzfTFUmSJCmLGbrbavQkmP6/8OI98JuZsH9fpiuSJElSljJ0H4nj3gXTboZn74Q7Pwu1+zNdkSRJkrKQoftInfB+uOjH8NRt8LvPG7wlSZLUhCtStoeTLoLaWrhjJsQI5/83FBVnuipJkiRlCUN3ezllWvLxjplAhPN/aPCWJEkSYOhuX6dMgxDgN5cnPd4X3GjwliRJkqG73Z08NQnet18OsRb+6SaDtyRJUoEzdHeEky4CAtz+SSDCBTexYcc+Nmzb3eIpg/p2Z1BJj04rUZIkFa4NNbvMJZ3M0N1RTrow6fFe+AmIkV+VfJn/un91i4dfMWkMV04Z24kFSpKkQnXLo2u4/r4VLbabS9qfobsjvfmfgAALL+NTY/cx+bPfg6IurNywnS/Of4L/mnEqowf1AZLfKCVJkjrDJWcMZ8qJRwGYSzqJobujvfkCCIEeCy/jpC6B1Wdfz8MrNwLw8MqNvGVYGSMH9s5sjZIkqaAMKunRZPjI6EF9OGlIaafWUVlZyZw5c5g3bx5lZWXMnDmzvq26uprx48cftC+XGbo7w4nnw9SbqV3wcZ5dvp7f7fss0IXfLFvHb5atY85FpzBtwrBMVylJktSpysvLmTt3LpWVlZSXlzNnzpyD2q+77jqmTZvGggULWnXdefPmZV1Yd0XKTrJ60CQ+vefzTClawn93uZ5u7KU2Qm2EituX89LGHZkuUZIkFZjVG3fws7++BMDP/voSq7Msj1x99dVUV1dz3XXXteq8RYsWdVBFbWfo7iS3LVnL4ng6s/Z+ibOLlvPjrt+lB8ldwyEE5i9Zm+EKJUlSIbltyVomffcBfrNsHQC/WbaOSd99gAVZlkmmTZvG7Nmz0zq2urqaiooKKisrO7iq1jN0d5J1W3YSY+RPtafx8b1fYWLRC9zc9dv0YhcxRtZt2ZnpEiVJUoFYvXEHX719ef1f3YGs/Qv89OnTqa6uZtmyZUAyDnzhwoUsXLiQiooKFi9eXH/s4sWLqayspLKykoqKCioqKqiurj7seZ3B0N1JhvbrSQgBgL/WnsRH91RwUtFq/rfbbPqGnQzt1zPDFUqSpEJx25K19bmksWz7C3xZWRkAS5YsAWDWrFlUVlYydepU5syZw6xZs+qD9dSpU5k1a1b9+PA5c+bUn3+o8zqDobuTTJ8wjBhj/eMl8Xgu2fPPjArr+WWXf+eDb3YGE0mS1Dnq/gLfnGz9C3xdQJ47d+5BN0mWl5en1Wvd1vPai6G7k4wc2Js5F51CUYCi1C+WTzOKD+39V8b02Mrwu6bD9g2ZLVKSJBWEhn+BbyyEkFV/ga8L2+PGjQOSsLx582bmzZvHwoUL2bx5M5s3bz7sddp6XntxysBONG3CMCaO6M9//2klC5eu48JxQ/ncuefQI06Cn38Abn43fPR3UDok06VKkqQ8Nn3CMOb+eVWzbTFGZmTRVMZ1w0omTJgAJNMIPvbYY/XTCM6fP/+Q59dNR9ja89qbPd2dZEPNLp6u2sr23fs4a/RAAM4aPZDtu/fx9J6j2TTtt7BvdxK8t7yc2WIlSVJea+4v8HWfz7noFEZk0cJ9Dcdm181O0nDe7rqe8LobLRtbtmxZm85rb/Z0d5JbHl3D9fetOGjfF+c/Uf/5FZPGcOXH/3hwj/fA0Z1cpSRJKhTN/wV+dFYF7oqKCiCZrxuoHw5SXV1df4Nk3b7KykrGjRtHeXl5/ZSBdb3c6ZzX0UJLg+gLQQihBNi6detWSkpKOvRrbajZxYZtu1tsH9S3e7Ica80r8IvzYecW+Mhv4OiTO7QuSZJU2J6u2sr7fvAwv//8WRldBr68vJypU6cCSTjevHkzEydOrA/cda677jpWrVr
      "text/plain": [
       "<Figure size 800x494.438 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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Covariance:  0.003465486601483565\n",
      "Normalized covariance:  0.8360758153764554\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": 7,
   "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": 8,
   "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": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwkAAAGLCAYAAACMfN52AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAABJhElEQVR4nO3de3RV533n/8+jC5KQdHQkgQAbHFmASR2HJgLilTTNz4lF09qZaRtzGTedlWTSoDWraZ2sTpH9R6fpL7MWEW3SuLdZIs3EvaRdBuxp0jidFFynzcWTBBSXZrWJAUFBNiCQdHQkobue+WPvLc45nKuQzt7a5/1aC6R9Pd9z0d7P9zw3Y60VAAAAAHjK/A4AAAAAQLCQJAAAAABIQpIAAAAAIAlJAgAAAIAkJAkAAAAAkpAkAAAAAEhCkgAAAAAgSYXfAdwpY4yRdJekUb9jAQAAK0q9pNctk0YBt1nxSYKcBKHf7yAAAMCKtFHSa34HAQRNGJKEUUm6fPmyIpGI37HAZ1evXtUXv/hFffjDH9b69ev9DgcAEFDnzp3Tjh07JFoiAGmFIUmQJEUiEZIEaHx8XFVVVaqvr+fzAADIqK6uzu8QgECj4zJCpaamRu3t7aqpqfE7FABAgFVXV/sdAhBoZqX31THGRCSNjIyM8M0xAADISzweV0NDgyQ1WGvjfscDBA01CQiVmZkZDQwMaGZmxu9QAAABxn0CyI4kAaFy48YN/cmf/Ilu3LjhdygAgAAbHBz0OwQg0EgSAAAAACQhSQAAAACQhCQBAAAAQBKSBIROeXm53yEAAACsaAyBCgAASg5DoALZUZMAAAAAIAlJAkLl+vXr6unp0fXr1/0OBQAQYAyVDWRHkoBQmZ2d1ZUrVzQ7O+t3KACAAOM+AWRHkgAAAICiMcZEjTFtfseB7EgSAAAAsGjGmDZjTI8xxhpjho0x3e6/HmPMMWNMR8K+7ZJelHTCv4gLZ4zpMMYc8zuOYqrwOwAAAACsXNbaPkmdbu1An7W2y9tmjIlKOm2M6bbWHrHW9hpjPiqp4AK3MeaAtfbIkgWe32O2S9ovKSqppGo/qElAqESjUe3du1fRaNTvUAAAAeYOf4plZq2NSeqR1J2wOrbI0+2+03gKZa3tdZOeFVXzsRRIEhAqNTU1etOb3qSamhq/QwEABBj3iaKKSYq6tQoFc/swdGuJv8lP7RvhLS82zrApuLmRMeZAwmLUWns4j2MOur/uUko1VJp9T1hri54pIhzGxsb0L//yL3rzm9+suro6v8MBAATU2NiY3yGUks1yyn+xTDu4hfV2d3GXpBPW2pPucoecBKHNTRYk6ZB3Pnfd9919Yl6TpGzndJsReYnHZnefpyQdlNQpqajNmoKooCTBTRAWEgNjzB63jVm2Qn93Stu0Y8aYY9bavWn23SPngwAsyujoqL7+9a+rtbWVJAEAkFEQk4TWJ1/oUH5Nai5d/PSjf5xy7K9K+qikr+c49sTFTz/qFb7V+uQL1ZJ+J9P2O+UW1A9IejjHrj1yCvGHJR03xpw3xuyw1sastceNMTFJballTmPMCUndCYX/Y8aYPnc52zl73eSixzuXtbbLTR6gwmsSupTw4XXftM+762/jVtd0GGOiCdnjITkdWNrcji6J+5ZUhxAAAIAENXI6yOYylGZdnaTv5nF8unZW0RzbC7HT+1JZUrO77t5stQiuTiU/rz45Xxwfz3SAW6DfmVDjIDl9BzolnczjnOlex1xxloy8kwSvEJ9YsHdFjTHt1treDIe2uf+87X0J6xPPtU9O1U5ixxYAAIBSMaH8CqnpqkHG8jx2Is26WI7thTi1mBGIrLV9bn+AfXIK703uv2x2ShpyW6J4onKaHi32nHAVUpOQ6Vv+mJKTgAVu1tiY4TyJtQjtkk4VEAsAAECouM18FtXUJ7X5UQHHTSpDi5Bicvuv7vKaoxtj9ufYf6Fcaq1NW9tQ6Dld0XxjDrtCkoRMmddQlm3pdEo6mVIjsdNaeySf3uTGmCpJVQmr6gt4bIRcVVWVtm3bpqqqqtw7AwBK1qpVq/wOAS63/NdtrTUJq6PutkytVbwvmG/7EttNIIbyOGdMt5dhqWlwFXUyNbfGoEPSjoR1ewqslnpK0m8vdWwIh6amJj3++ON+hwEAWGID8UkNjE5l3N5SX6WWSHXe52tqoiy4TPJ5YaPpjknpw+qdx2ut0uf+7iUBfW7n4+OJk6y5CUe7e0yucw4lxuIemxpbIc8rVIy1Nr8dnQL+6ZSMTMaYYUkfzVTVk7LvCUl7E4asiiqhw4m7PJz6GCnnSFeT0D8yMqJIJJLXc0F4zc3NaXJyUtXV1SovL/c7HADAEvn9E6/q6RfPZtz+xMNb9Ynd9+V9vuHhYS9RaLDWxu88wtLlFto75QwfGpPTx/RQus7KCUOPdkg67I1W5DYN2iyn43FMTlLQI+lYQgLg9Vs9n/gFs7t+0D1uyCuT5nlOr5N1n7tPp5wkostaezJhxuU97vojcsrDoR8itZAkISppWFJj4ptujLGSdmTpuOzt1yOn2iexL8IBOW/eoLtqs5xhsrrkZIj5JB4RSSMkCZCkK1euqKenR52dndqwYYPf4QAAlkhiTcK5gTF9/NlX9Ln9b9GWFme460JrEl599VVt27ZNIkkA0sq7uZG1NmaM6ZNT3RJL2ZYrQTighATBzTijqVmYN5ZuPhO0AQCA0tESqb4tCdjSUqcH7m7wKSIg3MoK3L9bTnWLpIXCf+JEaW0pMzJ7E6RF5cyS1+Eudyl5+FNPtMB4AAAAACyxgjouuyMQHUycJCNl5rsOOQlAYueRYxnO1Zm47J7TG6LqmKSelMkxAAAAABRBwaMbZWsK5DYfOpKwHJOUsRNytmMBAAAA+KOoQ6ACy23dunV66qmnVFlZ6XcoAIAAa2lp8TsEINBIEhAqZWVlTKQGAMiprKzQbplAaeEvBKEyODiov/iLv9Dg4GDunQEAJWtoaMjvEIBAI0lAqExPT+v8+fOanp72OxQAQIBxnwCyI0kAAAAAkIQkAQAAAEASkgQAAAAASUgSECqRSESPPPKIIpGI36EAAAKsvr7e7xBKljEmaoxp8zsOZMcQqAiV2tpave1tb/M7DABAwNXW1vodQmi4Bf4uSQckxXRrctyopCZJPdbak+6+7ZI+727bXORQF8UYc9D9dZekPmttl5/xFAtJAkJlYmJCZ8+e1datW1VTU+N3OACAgJqYmPA7hNCw1vZJ6nSThaRCtDEmKum0MabbWnvEWttrjPmopGOFPo4x5oC19kjuPZeOG3fi8zlmjDlmrd1bzDj8QHMjhEosFtPzzz+vWCzmdygAgAAbGRnxO4SSYK2NSeqR1J2wOrbI0+2+03gK4SY4He5PzyFJe0qhuRRJAgAAWDEu3BjXM9+5KEl65jsXdeHGuL8BIR8xSdGUwnbe3D4M3ZKWtGCe2jfCW06Jsy3lcfsS1ocazY0AAMCKcPTUZT353JmF5ed7+/V8b7+6H9uuvTs3+RgZctgspxlSLNMObmG93V3cJemE149BUofcwrqbLEjSIe987rrvu/vEvCZJ2c7p9o3wEg+vb8RTkg5K6pR0xD1/Y0qoXnLQp5AjSQAAAIF34ca4nnzujObtrXXe713PndGu1ia1rqEzctC4BfUDkh7OsWuPnEL8YUnHjTHnjTE7rLUxa+1xY0xMUltqp2FjzAlJ3QmF/2PGmD53Ods5e93kosc7l7W2y00esumUdNLthxFqJAkIlcrKSm3cuFGVlZV+hwIAWEJHT12WMUay9rZtxhg9e+qyun72jXmfL5D3iU82dCi/dveX9MmRP0459lclfVTS13Mce0KfHDm5sPTJhmpJv5Nxe+F2GmMOyBm9qNldd2+2WgRXp6ShhOU+OTUIxzMd4BbodybUOEjSCfdcJ/M4Z+I2T8Y43cfrkLQj0z5hQpKAUFmzZo1+5Vd+xe8wAABLrH94QjZNgiBJ1lr1Dxc2WlFzc3PunYqvRk7hOpd0hds6Sd/N4/h0Q/9Fc2wvxKnFjEBkre1z+wPsk/P8mtx/2eyUNGSM2ZOwLiqn6dF
      "text/plain": [
       "<Figure size 640x395.55 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": 10,
   "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": 11,
   "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.7915235152330492\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt0AAAHECAYAAADlBpY8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAABSLklEQVR4nO3deXxU9b3/8dc3gYQ1hIDgAgoBFZeqBdTa9tZawS7W1ltBrF3sJtTe9trFSu29XW9bivXebr+q0EW7F6GLbW1dsIt6W70CtlYtLgQUooIIIawBku/vj5mE7ExCJmeW1/PxmMdk5syZfDiOJ+988znfb4gxIkmSJCl7SpIuQJIkSSp0hm5JkiQpywzdkiRJUpYZuiVJkqQsM3RLkiRJWWboliRJkrLM0C1JkiRlmaFbkiRJyrIBSReQpBBCAI4EtiddiyRJkvLWcODZ2M2qk0UdukkF7g1JFyFJkqS8Nw6o7WpjsYfu7QDr16+noqIi6VokSZKUZ+rr6xk/fjwcpHOi2EM3ABUVFYZuSZIkZY0XUkqSJElZZuiWJEmSsszQLUmSJGWZoVuSJEnKMkO3JEmSlGWGbkmSJCnLDN2SJElSlhm6JUmSpCwzdEuSJElZVvShe3JVCWH7c9DQ7cqdkiRJUq+FGGPSNSQmhFARP1OxreWJQSNgxNEwahIccSoceRocNR0GuUS8JEmSOqqvr2fEiBEAI2KM9V29ruhD96uOKd32+1t/zpCSfbBtPdQ9Ay88Ac/9HfZuh5IBMO4MmPQaOP51MPZkCAGAtZt3csuK9WzYuptxIwdz8fTxTBw9NOF/lSRJkvqLoTsDIYQKYNu2bduoqGg3mt3UBC8+BevugTV/hJo/p0L46OPZcdyFLGk4iy/8ZRcBiNByf+W5x3LpGUczpmJQv/97JEmS8klNTQ0LFy5k8eLFVFZWMnfu3Dbbli1bxtVXX83ChQsBmD9/PjU1NSxdujSpkjsoitAdQlgKzOpk06oY47QM9u86dLfXuC8Vvh9ZRsMjv2Zg4x7+2HQaP2g8j3uaXkJs1R5/2VnH8Lk3n9yzf4wkSVKRmjZtGtOnT2fRokVtnq+rq+Pyyy9vCdnLli2jpqaGq6++us3rFi9e3Caw96dMQ/eA/ispK7YAM9P3zeYAS/r8O5UOhOPOg+PO4ythLttXLOGdpXfw/bKF1DQdzuLGN/LzxlfRGPL9kEqSJPWvqqqqTp+vrKzk9NNPb3k8a1ZnY61w1113JRa6M5XvCfGuGOPy5gchhErgxRjjqmx+02e2w11Nr2ZJ49lMDU/y3gG/40sDvsuHBvySG/e/iRfqLs7mt5ckSSpoNTU1VFVVUVlZyYwZM7p8XV1dHQsWLKCmpqYfq+udvJ4yMMa4rN1T18QYr832992ycy+prpzAqngc/7bvw7x270JWNB3P5wbczOeevgweXprqC5ckSVKPLF++nC1bUo0MU6dOBWDVqlVMmzaNadOmtXldTU0NNTU1zJ8/n/nz51NXV5dEyQeV1z3drYUQpgLVnQTx1q8pB8pbPTUc2JBRT3crK9ZtYfaiv9LZoZtUUsuyyXcx8pk74cip8NovwTFnZfzekiRJxWbmzJnU1NQwa9aslgso16xZQ3V1dZvXLV++nHnz5rFmzZo2z82fP5+VK1f2d9lA8fR0t3ZNjHH2wV4DfOZQv9H0CVVce9EpzP/5w4QQiDG23L//La9n5PS5sO4+uOM/4KbXwUlvgdctgOGHH+q3liRJOmDvLtj8RNJVpIw+DsqG9Hr3GTNmtMxScu21nTcutA/h+aQgQncIIdP/AguA/2n1eDiwoTffc/b08Zw+oYolrebpnjN9PBOa5+me8Eq4/I/w8BK48z/h/50BMz4N094DJXnd1SNJknLF5idg8dlJV5Ey98+phQX7QHd93PmqIEI3qWkDD9pBH2NsABqaH4f0Ije9NWH0UOa/bkrXLygpgdPeCse9FpZ/Bm77GPz9Z3DB12HsSYf0vSVJkhh9XCrs5oLRx/XZWzX3cfdGTU1NTo6IF0rongMsOuirkjKkCt70TTj1rfCbD8Ois+GcT8IrroSS0qSrkyRJ+apsSJ+NLheKVatW5WToLpQ+h2raztWdm455Obz/Xjjr3+Duz8NNr4cX1xx8P0mSpALXPFtJdzqbmaS6urplysBcHeWGwgndlUBdwjVkZkA5zPwcvPv3sGMj3PgvsOJ7dDoViiRJUgGrqalh3rx5LF++vGVmkuXLl3f62lWrVrUsAz9v3ryW56urq5k7dy7z589n+fLlh9Sakk0FMWVgCGErcG5PF8Xp0TLw2dCwI3WR5cqbYMob4c3/DwaP7P86JEmS1CuZThlYECPdMcaR2V6FMivKh8EFX4NLfgLr7oVFr4INycwxKUmSpOwpiNCd96acD/PuhaFj4HvnwV+/ZbuJJElSATF054qRx6T6vM98P9zxSfjZpbC7LumqJEmS1AcM3blkQBm89ovw1p/B0/8L334NvPB40lVJkiTpEBm6c9Hxr0+tZllaBt8+F1bflnRFkiRJOgSG7lw1ahK8bzlMOifVavLHBdDUlHRVkiRJ6gVDdy4rHwYX/wBe8yn480JY8jbY0+VMNJIkScpRhu5cFwK86iq49BZY97/wvddC3TNJVyVJkqQeMHTni+POg/fdBXt3pvq8a53PW5IkKV8YuvPJYcfD++6GkRPgpvPhsVuTrkiSJEkZKIhl4Hsr8WXge2vfbvjVB+DRX8CMz8Errky1oUiSJGVgU/0eNm1v6HL7mOHljKkY1I8V5a9Ml4Ef0H8lqc8MHAwXfTc1w8nyz8CLT7H2ZV/gloeeZ8PW3YwbOZiLp49n4uihSVcqSZJy0I8feIav3/1kl9uvPPdYPjLzuH6sqPA50p2PI92t1N//A4be8RHubTyJD+67kp0MIgCR1P8wl55xtL+pSpKkNlqPdD+1aQcfXvI3vjbnNCaPGQY40t0TjnQXies2TmVNw8dZNPCr/GjgF3jP3qvZQuoXiK8tf5KtO/fyuTefnHCVkiQpl4ypGNQhVE8eM4yTjxrRr3XU1NSwcOFCFi9eTGVlJXPnzm3ZVldXx7Rp09o8l88M3XkuAH+NL2HO3k9xc9lClpV9lnfu+wQb4hhKbPOWJEk5rLq6mkWLFlFTU0N1dTULFy5ss/3aa69l9uzZLF26tEfvu3jx4pwL685ekuee27aHCDwaJ3LR3s8RiPyi7LOcGNa1bJckSerM2s07ufkv6wC4+S/rWLt5Z7IFtXP11VdTV1fHtdde26P97rrrrixV1HuG7jy3Zedemtvyn4ljmbX3szwfR7Kk7L84MzzKlp17ky1QkiTlpFtWrOfc//4Tv1i1AYBfrNrAuf/9J5auWJ9wZW3Nnj2bBQsWZPTauro65s+fT01NTZar6jlDd577xOuntJkt8EVGcMneT/FQ02RuHriQa6c8lVxxkiQpJ63dvJNP/PxhmiI0pQfvmr+e//OHWZdDI94XX3wxdXV1rFq1Ckj1gS9btoxly5Yxf/58li9f3vLa5cuXU1NTQ01NDfPnz2f+/PnU1dUddL/+YOjOc9MnVHHtRadQEqC0JFASoKFkMJfv/zjPj3st1X/+d1hxU9JlSpKkHHLLivWELtb4CCGwJIdGuysrKwFYsWIFAPPmzaOmpoZZs2axcOFC5s2b1xKsZ82axbx581r6wxcuXNiyf3f79QcvpCwAs6eP5/QJVSxZsb5lnu4508dzTNUFcPt8+O2HoWE7vOLfky5VkiTlgA1bd9PVtNExRjZs3d3PFR1cc0BetGgRVVVVLc9XV1ezfPlyZs2a1e3+vd2vrxi6C8SE0UOZ/7opHTe8/loor4C7PgUN9XDOf7h6pSRJRW7cyMGpke5OgncIgXEjBydQVeeaw/bUqVOBVFiuqanhlltuoaqqii1btrBly5aDvk9v9+srhu5CFwKc+ykYVAF3fRr21MPrvgwldhZJklSsLp4+nkV/XtPpthgjc6aP7+eKutbcVjJ9+nQgNY3ggw8+2DKN4JIlS7rdv3k6wp7u19cM3cXiFVdC+XD47Udh7w644BtQ6n9+SZK
      "text/plain": [
       "<Figure size 800x494.438 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": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(Obs[-0.08(35)], Obs[0.19(25)])\n",
      "(Obs[1.85(35)], Obs[0.34(25)])\n",
      "(Obs[4.01(35)], Obs[-1.39(25)])\n",
      "(Obs[6.10(35)], Obs[-1.30(25)])\n",
      "(Obs[8.08(35)], Obs[-0.37(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": 13,
   "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": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit with 3 parameters\n",
      "Method: ODR\n",
      "Sum of squares convergence\n",
      "Residual variance: 0.5988333933914471\n",
      "Parameter 1 : Obs[-0.01(29)]\n",
      "Parameter 2 : Obs[-0.165(55)]\n",
      "Parameter 3 : Obs[0.89(23)]\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": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk8AAAFyCAYAAADsyz6AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9mUlEQVR4nO3deXxU1cHG8d9NgEQgCyBBloIGgoDsqLiDQtx3UVTq1ipqbRH1rdSli+1rLbaiYmuV12pbRVSw1rpUDbjgCiiCKMoWNeBCEMjCkgDJef84Ewwhy0wyM+fOnef7+eSTzMyd5DGS5Jlzzz3HM8YgIiIiIuFJcR1AREREJJGoPImIiIhEQOVJREREJAIqTyIiIiIRUHkSERERiYDKk4iIiEgEVJ5EREREItDKdYBo8zzPA7oB5a6ziIiISELJAL42TSyCGbjyhC1O61yHEBERkYTUA/iqsQOCWJ7KAdauXUtmZqbrLCIiIpIAysrK+MEPfgBhnLkKYnkCIDMzU+VJREREok4TxkVEREQioPIkIiIiEgGVJxEREZEIqDyJiIiIREDlSURERCQCKk8iIiIiEVB5EhEREYmAypOIiIhIBFSeRERERCKg8iQiIiISAZUnERERkQgEdm87kUZ99gKsXQgVJVBRBimpMHg85OVDRSlUlkNWD9cpRUTEh1SeJPgqy+GTf8OK/8Lp90G7TrD8WSh6F9p2grQMqK62pQng0+fg2Wugw/6QdwIMOR+6DQPPc/lfISIiPuEZY1xniCrP8zKB0tLSUjIzM13HEZc2roF3/wwfPQU7t8EPRsKpd0NOf1uWUho4a71tE3zxFnw+Hz79D2xZD8MugjP+HN/8IiISN2VlZWRlZQFkGWPKGjtW5UmCa+a58M1SGHEpDL+4eafhqnZB4evQtiN0Hw5F78GWYuh/mkaiREQCJJLypNN2EhxVO+Hdv0DXwdD7ODhtOuzTAVqnN/9zpraCvLHf3175Mrw1Dfrkw8l3QsfclucWEZGEoqvtJBjWL4cZx8K82+Cbj+x9mV1bVpzqM+ZXcP7jsOEzeOBoWPpEdD+/iIj4nkaeJLEZAx88Ai/dBB0OgCtetZO7Y8XzoN8pcMAx8OLP4ZVboe+JsE82xWUVFJdXNvkpcjLSyMmMcqkTEZG4UXmSxLZzG7x1Dwy9EE74PbTeJz5fNy0DznoAytfDPtmwfTMz5xdy71vfNvnUa8fkcV1+39hnFBGRmNCEcUlMW7+zc5wyu8L2EltgXHr6CorXLKX4xL9Cpz4ArC7ewuQnl3DP+KH0yWm/+1CNPImI+I8mjEuwfbcKHjsHOh8IE2a7L04A+b8lZ8O55Lx4Nlz0jL0yL6RPTnsGds9yGE5ERKJJE8YlsaxfDo+cDK3S4ZS7XKf5XmZXuOR52DcP/nmGXdJAREQCSeVJEsfXS+Dvp0D7LnDZi5Dd03WiPe2TbUedug6Brz90nUZERGLE96ftPM/LNsaUuM4hPvDdSruu0g/n2PWb/CgtAy5+1u6V91Wp6zQiIhIDcSlPnudNrHUz2xhzZxPHjwUKat0uBPKNMYUxiih+tn2zLUuDz4OB59hi4mc1+Va8BLS3k9vRnCcRkaCI+Wm7UHHKNsbMMMbMAAo9z5vaxNOygRGht97GmN4qTkmqdJ1djPK9B+xtvxen2nqMsO9fvhl2bnebRUREoiYec56mAHNqbhhj5gATGz58t0JjzGKVpiS2fTM8ejbg2b3kEk27zvb9ps/h2Wvsgp4iIpLwYlqePM/LBnLrKUDZnucNr+cpIlbVTph9KWxZbydhZ3V3naj5jr0ZPn4aFj3kOomIiERBrOc8NbRraknoscWNPPc8z/M2hT4+xBgzpb6DPM9LA9Jq3ZURaUjxoQUPwBdv2eK0bx/XaVqm97HQ4W+Qd7zrJCIiEgWxLk8dG7h/UyOPARQSOm0H4HleR8/zHjTGXFnPsTcBv25ZTPGdQ66A/QbZPeSCYNA4+37zl5Daxq4LJSIiCcmX6zyF5jrVHpWaC0wMnQas6w7spUw1bz1in1Bi5ou34dtl0Dodcke7ThNd1dXw+Hh4+nKo2uU6jYiINFOsy9OmBu7v2Mhje6k1Z2qv04DGmEpjTFnNG1AeeUzxhbJvYPYl8NodrpPERkqKXRW96B2Y3+hqHSIi4mOxLk+FsHvieG3ZNY/V5Xletud5mz3Py619X4zyiV9U7YKnfwwpreC0e12niZ39j4RRv4D5f4R177tOIyIizRDT8hRaGbyQeuY31TktV9f7da7Qyw3jOZLIXvtfux/cuIehfWfXaWLr6Bug61B4frKWLxARSUDxmPM0FRhXcyO0aOaUWrdza69AHipcBezpptrPkYDZtgkW/xPG/Ap6HeE6TeyltoJzHoJxj4DnuU4jIiIRivn2LMaYGZ7n3Viz0jjQqc6yA2OxxWhGrefc6XnejaGbvYGC0OrkEkRtO8JP3oO2+7pOEj+detv3O7dDyVro3NdtHhERCVtc9rZrbC+7UCnaqxg1tf+dBIAx8M50GHYRtM9xncaNF26AL9+Bn7wLrfdxnUZERMIQl/Ik7hSXVVBcXtnkcTkZaeRkpschUS0fPAIFv4LO/aFvki4geeRkWDYb5v8JxvzSdRoREQmDylPAzVxQxL3zVjV53LVj8rguP46njr5bDS/fAiMuTd7iBPZ03VHXw5t32YU0c/q7TiQiIk1QeQq4CSN7kj+gy+7bq4u3MPnJJdwzfih9ctrvvj8nI62+p8dG1S54ZiJk7AfH3x6/r+tXR10HH8+B56+Dy/6rSeQiIj6n8hRwOZnp9Z6O65PTnoHdsxwkAr58G779GC57EdLaN3180LVOh9P/DBWlKk4iIglA5UniL3cUTF4GGV2aPjZZ9DrcvjcGqqvscgYiIuJLvtzbTgKqugo+nGlP26k47c0YeOpimPcb10lERKQRKk8SP4segmd/At8sdZ3EnzwPug6G9x6wE+pFRMSXVJ4kPkqKYO5tcMjl0GOE6zRxMWth0R7vw3L4TyGjK7xyS4xSiYhIS6k8SewZA89Nhn06wJhfu04TluKyCu4uWElxWUWznj993ipmLrClaeaCIqaHsVwEYBfKPP53sPIlipcWtCiDiIjEhsqTxN6aefbt1LshPdN1mrAUl1dy77xVYS0wWtf0eauYVrByj/umFawMv0ANOAP6nkhx8bfNziAiIrGj8iSx13sMXPpCUiyGWV9xqhF2gfI8uOAJGHB6lNOJiEg0qDxJbG1cY8vA/ke5ThJzjRWnGhEVqMot9uNdGnkSEfETLSYjsfPtMnhwFJz3T+h/qus0zbK6eEtYx81aWLR7jlNTphWsZH1ZBRcc2rPxr73uW/vBJ89AryvD+twiIhJ7Kk8SG9XV8MIN0KkP5CXu6brJTy6JyeeduSD8ssWHj8Fx50O6oxXhRURkDypPEhtLH4e1C+CS56FVG9dpmq3uHoANiWTkCeyeg02OPIX2IWRXJbw9Hcb8MuzPLyIisaPyJNG3bRMU/AoGnQsHHO06TYuEuwfg7WcNoktmepNzngCuz+/LpDF54YcYNA7euxsOvwbadgz/eSIiEhMqTxJ9rdJgxGVw6BWuk8RVTSFqrEBFXJwAhl4Aw49UcRIR8QmVJ4kuY6BNu6Q9xdRYgWpWcQJIy4TuR9rvbdXOhD4NKiISBFqqQKLHGHjqInj/YddJWiwnI41rx+SRk5EW8XMnjcnj+vy+e9zXnOK0RwZj4LGz4dXfRpxHRESiS+VJomdVAXz6HLTLcZ2kxXIy07kuvy85menNev6kMXlMGGknhE8Y2bNZI057ZPA86DYcFv0NtmxoViYREYkOlSeJjqqd8PLNsP/R0O8U12l8oeZquqauqgvb4deAlwLv3hedzyciIs2i8iTRsehvsHE1nHiHHSWR6GvbEQ6dCAsfgq0bXacREUlaKk8SHesWwfCLYb9BrpME2+E/hXadoHi56yQiIklLV9tJdIz7G+za4TpF8LXrBJOWQEqq6yQiIklLI0/SMt+thhUv2avBdAl9fKSkQvl6+Hy+6yQiIklJ5UlaZu6v4cWf2wnjEj/z74Q5P4K
      "text/plain": [
       "<Figure size 640x395.55 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": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameter 0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAFbCAYAAADlb5X5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvvUlEQVR4nO3deZxcVZ338c+t6s6+L5ANqGwSIGUWVsNqABVaRhBQB1nUQXCeuD8qhYgW6miDj4ADasS1Z9wQdFQsRNmGJawh6aSTkJCtspA96a5snfRS5/mjqqHJRi9V9bu37vf9etWrqU7q1ldMim+fc+45nnMOEREREUsR6wAiIiIiKiQiIiJiToVEREREzKmQiIiIiDkVEhERETGnQiIiIiLmVEhERETEnAqJiIiImFMhEREREXMV1gFERKx5nndD/h8HAUOB7zrnGswCiYSQComIhJrneV8B7msrIJ7nDQJuB240jCUSOpqyEZGwu7D9aEj+n8eZpREJKRUSEQm7IflREhExpEIiImF3E3C753mPep43yPM8TdeIGFAhEZFQc849BlwIXADUAy8751bZphIJHxUSEQk1z/PGAdOBwcB9wAPt7roRkRLxnHPWGUREzHie94Bz7sp2zy8AHgXGa6REpHQ0QiIioeV53nTgLaUjP4VzB7kpHBEpERUSEZGDreSAoiIixaVCIiKh5ZybB0zPb4bW3sn5kRIRKRGtIRGRUMuXkZvzT7ejreNFTKiQiIiIiDlN2YiIiIg5FRIRERExp9N+ReQgsURqADAI6A8MyH/t3+55byBK7oeayJcbemeBLNDa7rGf3M6nO8itzdiRf2RmzZ6puWIReQutIREJmVgi1R+YBIwFRgOj2n1te/TtzDW/1NDLeXheB397K28WlbaysgZY3u6xetbsmc2dySAiwaZCIlKmYonUaOAkcuXj+PzXSeQKR0F1spB0RCsHl5TlwJJZs2euKeD7iIhPqJCIlIFYIjUEOBU4Lf/1VGBEqd6/CIXkSDYBL+YfLwAvz5o9c3eJ3ltEikSFRCSAYonUROB84BxyJWS8ZZ4SF5IDtQKLyZWTtsdSrVMRCRYVEpEAiCVSI8gVkLbHsbaJ3sq4kBxKPfA48AjwyKzZM183zlNQnud9BWjIP93hnHvQMI5IQaiQiPhQLJHygBnAZcD7yK0F8S0fFpIDLQL+DjwEPDdr9sxW4zxd5nneo8CNzrlV+cMBX3HO+fnfvUiHqJCI+EQskaoE3g18EPgAJVwD0l0BKCTtbQNSwF+Bf8yaPXOPcZ4O8zzvBmC8c+6mdt+bnj+TRyTQVEhEDOVLyEXAlcD7ye39ETgBKyTt7QP+AvwK+Oes2TOztnGOzPO8euBKHfwn5UiFRMRALJGaAnwcuAoYbhyn2wJcSNp7Hfg1UDNr9sxXrcMcKH8IYD258jok/+23jJaIBJkKiUiJxBKp4cBHgY8BU2zTFFaZFJL2XiI3avL7WbNn1htnAXJTM8Ar5NaP3Jf/3gX551eahhMpABUSkSKLJVIzgVnAJUClcZyiKMNC0mY/ubUmvyC33sTsAzNfPh4FTm6/ZsTzPEdupGSVVTaRQlAhESmCWCLVB7gG+DQw2ThO0ZVxIWlvCfB94NezZs9sKvWbe543DlgJDHbONbT7fj3wSd36K0GnQiJSQLFEaiTwGeBG3pznL3shKSRtNgL/CcyeNXtmQynfOD8acqgRkitVSCToVEhECiCWSE0AbiG3SLWHcZySC1khabML+Blw16zZM9eV4g09z3sFuKn9XTaaspFyoUIi0g35LdxvJVdEosZxzIS0kLRpAe4Hvjdr9swFxXyj/DqSK51zN+af3wBcqEWtUg5USES6IJZIvQP4GiEvIm1CXkjaewhIzJo9c0mx3qBtc7S257rtV8qFColIJ+SnZr4B/CsqIm9QIXmLVuCXwNdnzZ650TqMSFCokIh0QCyRGgx8ndztu2V56253qJAc0l7gTuCOWbNn7rIOI+J3KiQiRxBLpCqA/0NuVCQ0d810lgrJEW0Bvgn8ZNbsmS3WYUT8KmIdQMSvYonUJeROif0BKiPSdUcB9wKLf/ipJy63DiPiVxohETlALJGaRO4/IOdbZwkKjZB0yrPAp2bNnrnYOoiIn6iQiOTFEqkewM35R0/jOIGiQtJpzcD3gG/Nmj1zn3UYET/QlI0IEEukzgTmA0lURqT4KoGvAnU//NQTGokTASqsA4hYiiVSA4DbyW31rp/wpdQm9Nu9PvnqpBOuBr5wwtJXG6wDiVjRCImEViyRej/wKvApVEbEgstunVb7g8nAx4AltZNPuNg4kYgZjZBI6ORP4r2T3KiIiJmJKx5cUdmy9135pyNrzo/cek1N/IPA5+uuq9ttmU2k1DRCIqESS6SmA/NQGRFjvfZtf/GY159qKyOsG8acx6ZHzgD+DVgQr4mfaZdOpPR0l42EQiyRigBfBr6FdlotON1l00nOZWa88LXGXvsbRgC0emz65OeivXb39ga1+11Z4A7g63XX1TVbxBQpJY2QSNmLJVJjgMeBalRGxAeOW/vPhW1lBODHVZF1B5QRyH0+J4Cn4jXxUaXMJ2JBhUTKWiyRuoDc7bznGUcRAaCyeXftuNV/Pavt+aqjefbpeOTUI7zkXcC8eE387OKnE7GjQiJlK5ZI3QQ8AgyzziICgHP7ps+/a5CXv6ur1eP1266KvrMDrzwaeCJeE/9ccQOK2NFdNlJ2YolUP+BXgM4NEV8ZsfnFF/ru3XQegAN396WRzY29vNEdfHkFcHe8Jn4a8Mm66+r2FiuniAWNkEhZiSVSxwMvoTIiPhNp3b900rLfvDFV89ponnlxUmR6Fy51FfB8vCY+rnDpROypkEjZyJ/O+xJwgnUWkbdwrnXqgntdxGUrAFoirPn2R6Ind+OK7wTmxmviFxUmoIg9FRIpC7FE6tPAn4EBxlFEDjJkx5JnBu1cdQKAg+wdV0Qy+3t4fbt52cHA3+I18Vu6n1DEntaQSKDFEikP+H/AF62ziByKl21ZE1/809PbntfFvGdqx0fOLdDlI8C34zXx0cCn666ryxbouiIlpxESCaxYItULuB+VEfGxyUt+viOabe4N0Bxl1e1XRk5/u9d0wb8Dv4/XxHsU4dqd4nneo9YZJJhUSCSQYonUUOAx4ErrLCKH03/XmmeGb1s4DcBB63c+HGlsrvB6FentrgQejtfE+xfp+m/L87wrgAus3l+CTYVEAieWSB0HPAforA/xL5fdPHXBPW/sMfLKBO/ZxcdFTiryu54PPBmviQ8v8vscxPO8QYDu/JEuUyGRQIklUhOAp4F3WGcROZLjX7t/dWVL40CA/RW8ducHI+96u9cUyMnAnHhNPFai92vzIeC+Er+nlBEVEgmMWCI1CXgKONY6i8iR9G7c+sLojc+eAeCg+ZtXRV1L1Cvl+o6J5ErJ5FK8med504G5pXgvKV8qJBIIsURqMrkyokPGxN+ca5g+/643pi6eO8Gbs3y0d7xBklHA0/Ga+IwSvNcpzrl5JXgfKWMqJOJ7sURqGvC/wFHGUUTeVmzN3xf1bMocBdBYyav3/EvkrLd7TRENBh6J18SPdHhft3ied4VzTlM10m0qJOJrsUTqVOBxYKh1FpG3U9m0a/64dOosAAf7v3F1tEc24lnv99SfXCmJF/rC+YWsDYW+roSTCon4ViyRmgL8k9xPeSL+5lzj9Pl3DWl7+uQ7vefTI7zxlpHaGQL8M14Tn1jg634IuNDzvK94nvcV4HaA/PMrCvxeUuY855x1BpGD5A/JexpN0wTClxp6OQ/Ps85hadSGOU9Neu235wLs6UndJ74QPcl5nt9+6FsLnFV3Xd26Ylzc87xxwErnXKj/LEjX+O0vi0jbPiOPoTIiARFt2bfk+Nd+1zZVs/fWa6L9fVhGIHeH2uPxmvjRRbr+oCJdV0LAj39hJMRiidRwctM0Y6yziHSIcy1TF9wT9XBRgEdO9l5eP9yLGac6konAo/GaeEGnQj3Pu4E3p2we8DxPO7ZKp6iQiG/EEql+wMNo0zMJkKHb6+YM3JU+HmBnb2p/eWH
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Parameter 1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAFdCAYAAAAzNnbkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAyIElEQVR4nO3deXiU1d0+8PvMZGbYA2EN6xBACCgom8puFK2mLv1VrfZ1qVYxGKvVqsyrVedVq9G6azRq3attRVxaR22hkV2QVUESZQsYIAESGAhknTm/P2aCEVmyzMz3We7Pdc1FZpI8z+0FJvec5zznKK01iIiIiCQ5pAMQERERsZAQERGROBYSIiIiEsdCQkREROJYSIiIiEgcCwkRERGJYyEhIiIicSwkREREJI6FhIiIiMSxkBAREZE4FhIiIiISx0JCRERE4lhIiIiISBwLCREREYljISEiIiJxLCREREQkjoWEiIiIxLGQEBERkTgWEiIiIhLHQkJERETiWEiIiIhIHAsJERERiWMhISIiInFJ0gGIyBi8voACkAwgJfro3ODjhs/bAnBGH0kAnL/Z5ynuGnb0AhAGEIo+DgDYCyB42J/1H+8BUJydl1GRgP88IjI4pbWWzkBECeL1BdoBGBR9DGzw8QAAXREpGU12Q9CzrIN2jGlmrHIARQC2NHjUP9+UnZcRbOZxichEWEiILMjrCyQDGAtgFIDB+KF8dI/H+VpYSI7newBfHfbYkJ2XEY7T+YhIAAsJkcl5fQEngJMAnArgtOifQwCoRGWIcyE5kgMA1iJSTpYBmJedl7E+gecnohhjISEyGa8v0AnAFACnI1I+RiEyr0OMQCE5ku0A5gOYh0hBKRDOQ0RNwEJCZAJeX+AkAJkAzgMwDs2c6xEvBikkhytFpKDMBfBJdl5GkWgaIjomFhIiA/L6Am0BnIlIATkPQB/ZRMdm0EJyuK8AfAjgw+y8jNWyUYjocCwkRAbh9QW6A7gEwM8RuSTjEQ3UBCYpJA1tAfARIgVlfnZeRkg2DhGxkBAJ8voCHgAXArgawDkw2KWYxjJhIWmoHMD7AF7PzstYJB2GyK5YSIgEeH2BcYiUkEsBdJRN03ImLyQNfQfgdQBvZudlbBPOQmQrLCRECeL1BfoBuCr6GCgcJ6YsVEjqhQD8B8BrAP6ZnZdRLZyHyPJYSIjiKLoc+3kAbgWQgQSuDZJIFiwkDZUD+CuAp7PzMjZJhyGyKhYSojiIzg25AsBtAIYKx4k7ixeSemEAHwB4PDsv4wvpMERWw0JCFENeXyAFwI0AbkKclmk3IpsUkoa+APA4gA+4hD1RbLCQEMWA1xcYgMhlmWsAtBGOk3A2LCT1NgJ4CsBr2XkZB4SzEJkaCwlRC3h9geEA7gXwCwAO4ThibFxI6pUB+DOAZ7PzMg5KhyEyIxYSombw+gJ9ATwI4H9g4yJSj4XkkFIAOQDysvMyqqTDEJkJCwlRE0TniNyFyBwR06ykGm8sJD+xDYAfkUs5XAWWqBFYSIgawesLtAJwCwAfLLCQWayxkBxVAYC7svMyPpQOQmR0LCREx+D1BZyIrKj6fwB6C8cxLBaS41oI4KbsvIyvpIMQGRULCdFReH2BsxC5g2KYcBTDYyFplBCA5wDcm52XsU86DJHR2H4yHtHhvL5Ad68v8A6A2WAZodhxInLZrzA3K//X0mGIjIYjJERR0WXebwDwMDhPpEk4QtIsnwPIzs7LKJAOQmQELCREALy+QDqAVwCcLp3FjFhImq0WkcuCfq5fQnbHSzZka15fIMnrC9wNYBVYRijxXADuALA6Nyv/VOkwRJI4QkK25fUFTgHwKoCThaOYHkdIYiKEyOXC+7PzMmqlwxAlGkdIyHa8voDy+gI+AF+CZYSMwwngjwCW5Gblp0uHIUo0jpCQrXh9gW4A3gRwjnQWK+EIScxVIbIi8FPZeRn8IU22wBESsg2vL3AGgNVgGSHjawXgCQD/zc3K7yMdhigRWEjI8ry+gMPrC/gBzAGQKhyHqCnOALAqNyt/qnQQonjjJRuyNK8vkArgbUR+sFOc8JJN3IUB3APgYV7CIaviCAlZltcXOAfAV2AZIfNzAPgTgPdzs/I7SIchigcWErKk6CWaTwF0FY5CFEsXAViWm5XPLQ3IclhIyFK8voDH6wu8DeA+AEo6D1EcnABgaW5W/qXSQYhiiXNIyDK8vkBnAB8CmCAcxXY4h0RMDoC7OK+ErIAjJGQJXl9gEIAlYBkhe/EBeCs3K98tHYSopVhIyPS8vsBERMrIQOksRAL+B8AnnOxKZsdCQqbm9QWuQGR9kRTpLESCzgQwPzcrv6d0EKLmYiEh04reSfMWAA5XEwEjENkHZ6h0EKLmYCEh04lujpeHyJ00RPSDPgAW5WblT5IOQtRULCRkKl5fQAF4CcAN0lmIDKojgP/kZuVnSgchagoWEjINry/gAPAXANdJZyEyOA+AWSwlZCYsJGQK0TLyKoBrpbMQmQRLCZkKCwkZXrSMvA7gauEoRGZTX0rOkw5CdDwsJGRoXl/ACeBNAFdKZyEyKQ8im/KxlJChsZCQYUXLyFuILPxERM3HUkKGx0JChhS9m+ZNAJdLZyGyiPpSMlU6CNGRsJCQUT0C4NfSIYgspn5OySnSQYgOx0JChuP1BW4EcId0DiKLag8gkJuV3y8RJ1NKTYs+7lRKPaKU6piI85L5sJCQoXh9gQsAPCOdg8jiUgF8lpuVH9c9oJRSdwJ4V2v9ktb6UQAPIzL6SfQTLCRkGF5fYCyAvwFwSmchsoEhAP6Zm5XfKo7nmKq13lv/JPpxWhzPRybGQkKG4PUF0gD8C0Ab6SxENjIewF9zs/Lj9bsgJTpKQnRcLCQkzusLdAbwKYBu0lmIbOiXAJ6M07FnAHhEKTVbKdVRKfUIuA8VHQULCYny+gKtAHwE4ATpLEQ2dnNuVv70WB9Uaz0HwFQAZwHYA2CZ1npTrM9D1sBCQtJeRmTYmIhkPZWblX9aLA+olEoDMBJAJ0R26Z6plJoWy3OQdbCQkBivL3ADgCukcxARAMAN4L3crPxYXjp9RGv9qNZ6r9b6BkRGS16MFhWiH2EhIRFeX+AUAE9L5yCiH+kF4O+5WfktvtNNKTUSwI8uz0Qv4TyKyCUcoh9hIaGE8/oCyQBmIrJqJBEZyxmIrBcSLxtxWFEhAlhISMZrAAZIhyCio7ojNyv/ly05gNZ6JYCRR1iZdVR0pIToR5TWWjoD2YjXF7gVwBPSOSi2bgh6lnXQjjHSOSim9gMYm52XUdjcA0TLyP9Gn5YB6Azg4YaLpRHVYyGhhPH6AqcBmA/AJZ2FYouFxLLWAhiTnZdRJR2ErI+XbCghooufvQuWESIzORFAjnQIsgcWEkqU1wH0kQ5BRE12c25W/lTpEGR9LCQUd15f4CoAP5fOQUTNogC8kZuV31k6CFkbCwnFldcX6I747ZNBRImRCuA56RBkbSwkFG+5AFKkQxBRi12Wm5V/sXQIsi4WEoobry/wS0R2EiUia3g+Nyu/q3QIsiYWEooLry+QgsjoCBFZR1fw/2uKExYSipenAHSXDkFEMXdJblb+2dIhyHpYSCjmvL7AuQCulM5BRHHzXG5Wvls6BFkLCwnFlNcXaA/gRekcRBRXgwDcIR2CrIWFhGLtfnABNCI7uDs3K7+fdAiyjiTpAJQYSqlpDZ521Fo/GutzeH2BgQCyY31cIjKk1gCeBnCRcA6yCI6Q2EC0jHTUWr+ktX4JwCal1CNxOFUOuFcNkZ1cmJuVf550CLIG7vZrA0qpjQCmaq03NXhtj9a6U6zO4fUFxgFYFKvjkblwt19b2wRgGHcEppbiCInFKaU6AkhrWEaiOiqlRsbwVI/H8FhEZB5pAG6WDkHmx0JifWlHeX3vMT7XJF5f4FcATovFsYjIlGbkZuV3kA5B5sZCYn1H20em/BifazSvL+AG8HBLj0NEppYC4HbpEGRuLCTUUr8D0F86BBGJu5X73FBLsJBYX/lRXk85xucaJbpfzd0tOQYRWUY7AP8rHYLMi4XE+jYBhya3NtSx/nMt8L8AYna
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Parameter 2\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAFjCAYAAAANRhA5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxEUlEQVR4nO3deZzcVZ3v//eprdPZqpNOCIQ9YekGKmEnrIkI43ad6yjg1RHrijPQGEBABL06/mrGUQevjsvYbDo6jTMqwvwARUXRlrAEZBW+JNUJJAFCFrKR7k6601ud+0d1MISk00tVfb5V9Xo+HvVoujpd9eaRTtW7zznfc5z3XgAAAJYi1gEAAAAoJAAAwByFBAAAmKOQAAAAcxQSAABgjkICAADMUUgAAIA5CgkAADBHIQEAAOZi1gEAhJNz7tJdPq3z3n/dLAyAiscICYC3GSwjdd7727z3t0la6Zy70ToXgMrlOMsGwO6ccyskne+9X7nLfW9476cYxgJQwSgkAN7COVcn6Q3vvdvtfi/pJO/9MybBAFQ0pmwA7G7WXu7fOsTXAGBMKCQAdjd1L/dvGeJrADAmFBIAAGCOQgJgd1v2cv/UIb4GAGNCIQGwu5XSm4tbd1W382sAUGgUEgBv4b3fqnzxeNt6Ea6wAVAsFBIAe3KjpAt2fjK4UdoNdnEAVDr2IQGwR86565W/1LdOUr33nkICoGgoJAAAwBxTNgAAwByFBAAAmKOQAAAAcxQSAABgjkICAADMUUgAAIA5CgkAADBHIQEAAOZi1gEAhNsNX7hu4h0D75giaZKkCZJqJY2TlJAUlxR57/Z4z7F9sbikgV1uPcrv9PrG4K1j4S3nshPjPgxu0y8N7pAr6WuD5wsBFY2dWoFqlEk6STMkHSTp4F1uBw3e9lP+cL26P+dmP/aB3i+fPdTD/W1nzcMzByJD/hlJOUnt+ktB2Xlbo/xhfjtvqxbecu6OUf6flbXB7fpv21lABk9cvtF7f5llLqAUGCEBKlkmGZd0lKRjBm/HDn48QlLNcB5isrqiBUoTkTRl8DYU39zUuk7SKv2lpCyX9KykZQtvOTdXoDxhdL73/us7P/Heb3XOzbIMBJQKhQSoFJlkjaQTJM2TdJqkuZKO1Bj/nU9w3cMqLgXkJM0cvJ2529e2NTe1/lnS05KeGvxYSSVlqnPu+l1LCVAtKCRAucokD5V0hvIFZJ6k45Vf11FQteopdSEZykRJZw3edtq1pCyS1LrwlnPbDbIVwg2SHnDOnS/pQkmfl8R0DaoCa0iAcpFJTpb0TknnS/orSbNL8bQ9Prby6J7bh5w2GOYaklIZUH705IHB22MLbzm3zzbS8DnnzlM+tyRd6L2/yzIPUCqMkABhlUlGJJ0q6V3KF5BTZfBvNqaBCaV+zjGKKj9ldZqkLyo/gvKQBgvKwlvOXWIZbiiD60VOVH6dzY2S7nTOXea9v802GVB8jJAAYZJJRiXNl3SBpL+RtL9tIMl7dR/e85Paof5MyEZI9mWVpJ9LumPhLec+ax1mV865O733F+7y+c7Rktne+5V2yYDio5AA1jLJmKRzlS8hH5A03TTPHhyx4/a+fsXie/t6mRWSXS3XX8rJC5ZBnHMnSvqw9/6G3e6/UdIKRklQ6ZiyAaxkkidJ+qSki5TfACu0ktresVnJUGccpaOUn9b5YnNT61JJdyhfTpbZxnqLFcpf+gxUNEZIgFLKJKdI+pjyRWSucZphe2fP/31lhT/w0L19vYxHSPbmCUk3S/pZKTdpc849oPxC1q273HcrG6OhGjBCAhRbflfUc5UvIX+j/LbrZWWqOrtWWIcorVMHb99sbmr9oaSbF95ybilGKS6U9HnnnCRtVn7k7IYhvwOoEBQSoFgyyVpJH5d0taQG2zBjU+86dqg6B1OnSrpO0meam1p/K+kmSb8q1kZsgyMjFBBUJQoJUGiZ5AxJCyVdLmmacZqCmObae6wzGHOS3j14e7m5qfU2SbcuvOXcLbaxgMpBIQEKJZM8VtK1kv5WwzwnplxMc+391hlC5DBJX5X0+eam1mZJ31x4y7mbbCMB5S9iHQAoe5nkCcok75X0gqRLVGFlRJLq1TFgnSGEJkn6nKRVzU2tX29uag3d5dpAOaGQAKOVSR4/WESekfTX1nGKqd51WEcIs4mSPqv8VM43mptaZ1gHAsoRhQQYqUyyQZnknaqCIrLTFHVaRygH4yV9RvkRk281N7UeYB0IKCcUEmC4MsmZyiR/qPzUzAXKL3SsCkm3nfVmw1er/JVVK5qbWv+5ual1onEeoCxQSIB9ySTHKZP8gvLbjH9C+cPbqsokdVFIRq5W0hckLW9uav1kc1Mrr7fAEPgHAgwlk/yQpKWS/llSuZ16WzATXE/FLdQtoQMk/UDSM81NredYhwHCit96gD3JJOdI+rakdxgnCYVx6i273WVDaK6kRc1NrT+VdN3CW85dax0ICBNGSIBdZZITlUl+R/kFq5SRQXH1j7fOUEE+ImlZc1Pr9c1NrXs9QRmoNhQSYKdM8q+UX7B6lapwnchQohqYZJ2hwkyUdKOkJ5ubWsvmkEWgmJiyATLJqZK+pfy5M9gDl98EDIU3V/lS8hVJX1l4y7nsiIuqxQgJqlsmeYHyi1YpI0NwTrHx6t5unaNCxSVlJD3R3NQ6xzgLYIZCguqUSU5VJnmXpDslsbPmMNRpO9u1FtcJyo+WfLG5qZXRa1QdCgmqTyY5X9Jzkj5kHaWcTHGdXdYZqkBC0pclPd7c1HqsdRiglCgkqB6ZZEyZ5D9LapV0kHWccjOVQlJKJym/b8mnrYMApUIhQXXIJA+T9JDyO2fycz8K09TeY52hyiQkfbu5qfUOtp9HNeCFGZUvk7xI0p8lnW6cpKxNcx291hmq1EXKry05xjoIUEwUElSuTDKqTPKbku6QlLSOU+6mufYB6wxVrEH5q3A+ah0EKBYKCSpTfm+R+yVdax2lUtS7DgqJrQmS/qu5qbW5uak1YR0GKDQKCSpP/hyapySdZx2lkkxRp3UE5H1K0sPNTa2HWAcBColCgsqSSV4oabGkw62jVJoprpPt9MPjVElPNze1si4KFYNCgsqQSbrBS3p/rvzQNgpssrooJOEyTdIfmptaP2AdBCgECgnKXyYZl9Si/CW9KJKJrpt1C+FTK+m/m5tar7QOAowVhQTlLZOcKOmXki62jlLpatVTY50BexSR9N3mptZvNDe1OuswwGhRSDBqzrnrB293OuduLHmATHKGpAclvavkz12FatQ33joDhvQZST9rbmqlOKIsUUgwKs65G733Xx+8XShplnPuzpIFyCSPUH7x6kkle84qF9cAhST8LpL0QHNT6xTrIMBIUUgwYs65OknnDX7c6WuSLnDOzSp6gEzyROXLSPGfC2+KKDfZOgOG5WxJi5ubWmdaBwFGgkKC0ZqltxaClbvcXzyZ5CmS/iBpelGfB2/jnCZENdBvnQPD0iDpQUoJygmFBCPmvd/qvZ/ivX9ml7t3FpGVe/qegsgkT5X0gKS6oj0HhjRZ29kdrXwcKemPlBKUCwoJCuUySb/33henkGSS85QvI5xJY6jObaOQlJejlC8lB1gHAfaFQoIxc86dqPw27RcW5QkyydMl/VYSaxiMTdG2LusMGDFKCcoChQSFcKOkk7z3Wwv+yJnkGaKMhMY017HDOgNG5WjlS8n+1kGAvaGQYEycc7dKuqxIZWSupF9LmlTwx8ao1Lv2HusMGLWjlV/oSilBKFFIMGrOuUsl3bhz3Yhzbtbg9M3YZZKzJN0v1oyEyjS1c5VNeTta+fNv2KcEoUMhwag45y5Q/mqXWc658wY/v0GFuMomk9xP+WkafpMLmWmug0JS/o6RdHdzUytnEyFUKCQYscEN0e5Ufu3IA4O3OyVdOuapm0xykqTfSDpiTI+Doqh3HdYRUBjzJf0HZ98gTCgkGLHBfUjcnm5jeuBMMiHpbkmFmfZBwU0RV/1WkI9I+qp1CGAnCgnCIZN0km6X9E7rKNi7Orctap0BBfW55qbWy6xDABKFBOHxj5I+bB0CQ5uk7ph1BhRcc3N
      "text/plain": [
       "<Figure size 640x395.55 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": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAq4AAAGLCAYAAAAGdhAeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAysklEQVR4nO3dX2xc53nn8d9DZy0pjIcUQ9C1RQq0bDWyamgDibaBIC1iUAoSZC8W0T8gCorsFivBUGC7CGDKV91e2TTQje1YgKXdLoqufSFZ6sUCQYBKlOo2WMBeSgiEQBJC/YNJ2SjBlciJGUppxGcv5gx1OJwz4hnOzJlz5vsBBHLOOeS8xyNSP7/zvM9r7i4AAACg2bUlPQAAAABgOQiuAAAASAWCKwAAAFKB4AoAAIBUILgCAAAgFQiuAAAASAWCKwAAAFLhS0kPoFbMzCQ9Lum3SY8FAACkyiOSPnOa2ze9zARXFULrRNKDAAAAqdQr6WbSg0BlWQquv5Wk8fFx5XK5pMcCAABSIJ/Pq6+vT+Id21TIUnCVJOVyOYIrAABABrE4CwAAAKlAcAUAAEAqEFwBAACQCgRXAAAApALBFQAAAKnQVF0FzGyXpGlJGyTJ3Y8mOiAAAAA0jaYJrmbWKek1d98WPHZJBFcAAABIqlNwNbPtkg64++4y5/aHHna6+5uS5O7TZjYYXLNV0ol6jA0AAADpVNPgGgTOvZI6FbzdX3J+v0Jh1cx2mdmwuw9JC+F1l6QDkpaEXgAAALQuc/faf9NC+Fx42z90/KqkHe5+LXTstruvLbmuU9L10uMPeM6cpJmZmZma75x1fWpWx0fHNXF7Tr1r12jPQJ+e6G6v6XMAAIDGy+fz6ujokKQOd88nPR5U1rAa1yCMbgiH1kCnmW119/PFA8HM6y0z2+7upxs1xnKOj47r0MkLMjO5u8xMRz66quGdW7R7oC/JoQEAAEiS+g/9/DFJj1W45PMbb3zv80aNp14auThrSelAYFrSBjMbkLTN3Q8Ex7sklYbchhq9cUtDJy/IXVJxZjr4+OrJC9rQ3a5t/V3JDRAAAKDggKS/qnD+ryX918YMpX4a2cc1KuHdCs4dl3TKzLab2bCk/1Jmdrah3vjFZUVVUrhLr//icmMHBAAAUN4RSduCP/uCY/tCx44kNK6aapp2WO4+rfudBB5YHmBmqyStCh16pNZjal/1kExSuexqVjgPAACQtKAM4HNJ6j/08+Lhyzfe+N75yC+qo0odplaikcH1VsTxrgrnKnlNlafEV2z27r2yoVUqzLjO3r1Xz6cHAABIlQd1mFqpRpYKXJMWFmmFdaq6WtbXJXWE/vSuYGxlHfruJpmVP2cmvfbdTbV+SgAAgKr1H/r5Rkk/Dh7+OHjcMO5+Pmhzeqoe379hwTUoBbimMrWu4Y4CMb7fXXfPF/9I+u3KR7nYQH+X3ty5RW0mPdRmiz6+uXMLC7MAAEDT6D/08/8k6bKkPw8O/bmky/2Hfv6jxAZVY/UqFYhKdMOSdkkqbkCwX9LQSp7IzA5KOqg6hfDdA316tr9Lx0J9XPcO9KmfPq4AAKBJBDOr/0OL81BxMc7f9h/6+S9vvPG9K40fWW3Va+esXSq0uDoi6Zy7H5Ukdz9qZq8Wd9CS9NXirlnVcvfDkg4XNyBY0Q1E6O9u19B3KAsAAABN6z9Lmlf5iTyX9BcqrA9KtZoG1+At//OqMIta3O4VAAAANdOvyu8+9zdmGPXVyMVZAAAAqI8bKsy4VjqfeqkPrmZ20MwuSvok6bEAAAAk5H8qOteZpL9t4FjqJvXB1d0Pu/tmSc8lPRYAAIAk3Hjje2Mq1LHOSyo2mr8XPP6LLCzMkjIQXAEAACDdeON7fyfpa5L+Pjj095K+FhxvtLr0DCW4AgAAZEQws/pu8PDdRs+0mtlWMxtWYaH+VjM7EnSTqolGbvlaF/Xu4woAAIDlWU6HqZVIfdijxhUAAKA1pH7GFQAAoNX1H/r5Y5IeCx4Wd03a1H/o58VLPr/xxvc+b/jAaozgCgAAkH4HJP1VybEPQp//taT/2rDR1AnBtQ6uT83q+Oi4Jm7PqXftGu0Z6NMT3e1JDwsAAGTXEUn/u8L51M+2SpK5e9JjqAkzy0mamZmZUS6XS2wcx0fHdejkBZmZ3H3h4/DOLdo90JfYuAAAwFL5fF4dHR2S1OHu+aTHg8pSvzirmXbOGr1xS0MnL2jepXvzvujjqycv6NyNW0kPEQAAILVSH1ybqavAG7+4rKgJbHfp9V9cbuyAAAAAMiT1wbWZtK96SBZxzqxwHgAAANUhuNbQ7N17iqoYdi+cBwAAQHUIrjV06LubZBFTrmbSa9/dVP4kAAAAHojgWkMD/V16c+cWtZn0UJst+vjmzi3a1t+V9BABAABSK/V9XM3soKSDapIQvnugT8/2d+lYqI/r3oE+9dPHFQAAYEXo4woAAFoWfVzTpSlmKQEAAIAHIbgCAAAgFQiuAAAASAWCKwAAAFKB4AoAAIBUSH1wNbODZnZR0idJjwUAAAD1k/rg6u6H3X2zpOeSHgsAAADqJ/XBFQAAAK2B4AoAAIBUILgCAAAgFb6U9ABQcH1qVsdHxzVxe069a9doz0CfnuhuT3pYAAAATYPgmrDJ/B198PGneufMmEySSzJJ7310VS8PbtQPnluvntzqhEcJAACQPEoFEvbu2St6e2RM7tK8a9HHt06P6fDZK0kPEQAAoCkQXBNmktqs/Lmo4wAAAK0o9cE17RsQ3Prdv63oPAAAQKtIfXBN+wYEXV/+dys6DwAA0CpSH1zTzlWoaS0n6jgAAEAroqtAwn78wlPqan9Yb48UugoUubTQVQAAAACSuWdjWs/McpJmZmZmlMvlkh5ObDemZnUs1Md170Cf+iv0caXvKwAAK5fP59XR0SFJHe6eT3o8qIzgmkLHR8d16OQFmZncfeHj8M4t2j3Ql/TwAABIDYJrulDjmjKjN25p6OQFzbt0b94XfXz15AWdu3Er6SECAADUBcE1Zd74xWVFTZK7S6//4nJjBwQAANAgBNeU6Wp/WFZhw4Ku9ocbOyAAAIAGIbimzGMdq1VpQ63HOlY3bCwAAACNRHBNGfq+AgCAVkUf15Sh7ysAAGhVtMNKqbh9XwEAwFK0w0oXgisAAGhZBNd0SX2Nq5kdNLOLkj5JeiwAAACon9QHV3c/7O6bJT2X9FgAAABQPyzOahHXp2Z1PFQTu2egT09QEwsAAFKE4NoCjo+O69DJCzIzubvMTEc+uqrhnVu0e6Av6eEBAAAsS+pLBVDZ6I1bGjp5QfMu3Zv3RR9fPXlB527cSnqIAAAAy0Jwzbg3fnFZUY0j3KXXf3G5sQMCAACoEsE147raH5ZF7BHbZoXzAAAAaUBwzbgne76itojkamZ6sucrDR4RAABAdQiuGTe4qUfzEbUC8+7avqmnwSMCAACoDsE14/5lbKpijes/j001dkAAAABVoh1Wxu17fr12bH5Un03P6R8v/qsm83fUk1utb29+VI93rlHPI6uSHiIAAMCymEdNx6WMmeUkzczMzCiXyyU9HAAAkAL5fF4dHR2S1OHu+aTHg8ooFQAAAEAqEFwBAACQCk1V42pm+4NPt0k65e4nkhwPAAAAmkfTBFcz2yrpVjGsmpmb2Vp3n052ZAAAAGgGdQmuZrZd0gF3313m3P7Qw053fzP4fIOkHZKKs6zXgmPn6zFGPNj1qVkdHx3XxO059a5doz0DfXqiuz3pYQEAgBZV064CwazpXkmdkgbcfVvJ+f0KhVUz2yXpWXcfCh53uvu0mXVKuu7ua2M8N10FamQyf0cffPyp3jkzJpPk0sLHlwc36gfPrVdPbnWygwQAoAboKpAuNV2c5e7ngxB6KuKSId2fUVVQFrA/9Hg6+PS/S1oyW4vGePfsFb09MiZ3ad616ONbp8d0+OyVpIcIAABaUMO6CgSzqBvc/VrJqc5gprZ43auSjrj76UaNDYuZpDYrfy7qOAAAQL01sh3Whojj08VzQenAeXc/bWZbw4EWjXPrd/+2ovMAAAD10MiuAl0Rx29J6gpC6oeSps1
      "text/plain": [
       "<Figure size 640x395.55 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": 18,
   "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": 19,
   "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": 20,
   "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": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit with 6 parameters\n",
      "Method: migrad\n",
      "chisquare/d.o.f.: 1.100354109100944\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAAHECAYAAAAUDc2xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAABXI0lEQVR4nO3deXxc9X3v//dXkiVvksbygjcZWV6wwSEg26QkNE6CTbmhSXrBS0ia4NggB0qbzVj43l+z9D5SV4akTdIAdmKgIaQxNrlN0tIkGJJwaQLFFgSy4MgeG1sCY2yt3mXp+/vjnJFH0sxoljNzZnk9H495zMz5nnPmo6PR0We+8znfr7HWCgAAAEBiivwOAAAAAMhFJNIAAABAEkikAQAAgCSQSAMAAABJIJEGAAAAkkAiDQAAACSBRBoAAABIAok0AAAAkIQSvwPwkjHGSJoqqdvvWAAAAJCzyiW9boeZuTCvEmk5SXSL30EAAAAg502X1BprhXxLpLsl6fDhw6qoqPA7FgAAAOSYrq4uVVdXS3FUOORbIi1JqqioIJEGAABAWnGxIQAAAJAEEmkAAAAgCSTSAAAAQBJIpAEAAIAkkEgDAAAASSCRBgAAAJJAIg0AAAAkgUQaAAAASAKJNAAAAJCEvEuk3zYp734kAAAAZKG8yzp/sXqM3yEAAACgAORdIl01yki9PX6HAQAAgDyXd4m0JJkzHX6HAAAAUHCCwaDWrVsnY4zGjRunhoaG/tuKFStkjFFDQ0P/+qHluarE7wDSwZxukzTL7zAAAAAKSm1trbZs2aLdu3dr0aJFamxsHNDe0dGh2267rf/54sWLNX78+CH72bp1q+rr69Meb6ryOJEGAACAH6qqqiIuDwQCWrx4cf/z5cuXR1zvySefzIlEOj9LO063+x0CAAAAXMFgUB0dHZKkpUuXRl2vo6NDDQ0NCgaDGYosNXmXSJ/vszKnj/sdBgAAAFy7du1SW5tTMVBXVydJampq0sKFC7Vw4cIB6wWDQQWDwf7a6lACno3yrrTj+CmrSko7AAAAfLVr167+3uWdO3dq//79A9rr6urU2NiodevW9S9bvny5AoGAgsHgkPrqbJR3ifSxU1YBSjsAAECuOndKOvZHv6NwTJgrlY5OatOlS5f2J8ObN2+OuE5tbW3SoWWDvEuk3zplNYceaQAAkKuO/VHausTvKBz1v5SmXpHybmLVReeyhBNpY8xSSeustTEH/TPGPGmtXTZoWfjllwFr7eZE2uNx7JRVEYk0AADIVRPmOglsNpgw15PdhOqikxEMBrO25zruRNoYUydplaSApJg/jTFmuaSlg5bVKyw5NsYsN8Y0Wmsb4mmP17FTluHvAABA7iod7UkvcL5oamrK2kQ67lE7rLVNblL7ZKz1jDEBRU60GyTtDNvfTkn1CbTHhUQaAADAf6FROmKJNCJHbW1t//B32dwbLaVn+LuVkraGLwgl19bawYMCBowxdcO1J/LiJNIAAAD+CE0RvmvXLu3atav/cSRNTU39o3qEj9xRW1ur+vp6NTQ0aNeuXSmVhaSbsdYmtoFTtrHRWrswQlvoJw1KarfWmrDle0LPw9Zvl3Sbu37Udrd3Op7YKj7ytpLOR28cLf2vN5K+yhQAAACFqaurS5WVlZJUaa3tirWu1z3Si6y1TRGWR54nUmpz24Zrj9uxU+4HA3qlAQAAkEaeDX9njFlurd06/JreMcaUSSoLW1Ten0ifPCZVTs9kOAAAACggnvRIuzXOHTFWidY9XOW2DdcezUZJnWG3lv5E+hTThAMAACB9vOqRXilpVliN9CxJMsZskFP/vMt9HrDWdoRtF3Dbg8O0R7NJ0lfDnpcfO2VbJEmnKO0AAABA+niSSA8u6TDG1EqqD59QxRgTlNPD3DFo26Z42qO87llJZ8NeQ6d6JFsyUoYeaQAAAKRRMqUd8Vz8F4iwrFHS8tATdwKWhgTa42ZHVVHaAQAAgLSKO5F2x3tulJPc1hljtgya0ju0Xr2cpFjGmB3ulOL9vdbGmHq35GNWeI/1cO2JcBLpY8lsCgAAAMQl4XGks5kxpkJSZ8+DN6hkTJW08jt+hwQAAIAc4uc40lnBjqziYkMAAACkVV4m0n2jx1MjDQAAgLTybEKWbGJHjXMmZAEAACgAR7vO6Gj32ajtk8rLNKliZAYjKgz5mUiPnuBcbNjXJxXlZac7AABAv0efP6SvPdUctf1T187RZ5bNzWBEhSF/E2nbJ51uk8ZM8DscAACAtProO2Zo2aUXSZL2HT2hT29/Sf+06grNnjRWktMjDe/laSI90Xlw4iiJNAAAyHuTKkYOKd2YPWmsFkyrzGgcwWBQjY2N2rp1qwKBgOrrL4yU3NHRoYULFw5YluvyMpHuGz3eeXDyqKRLfY0FAACgUNTW1mrLli0KBoOqra1VY2PjgPbNmzdrxYoV2rFjR0L73bp1a1Ym4HlZQHyhR/otfwMBAADIoAPHTurhXx2UJD38q4M6cOykvwENsmHDBnV0dGjz5sTm3HvyySfTFFFq8jKRVukYacRot0caAAAg/z22+7Cu/cov9IOmFknSD5padO1XfqEduw/7HNlAK1as0KZNm+Jat6OjQw0NDQoGg2mOKjl5WdohSRoz0amRBgAAyHMHjp3U3Y+/rL6wCatDjxsef1mLa6pUM2GMP8ENsnLlSq1bt05NTU2qq6tTMBhUU1OTJOmFF17QsmXLtHTpUknSrl27FAwGFQwG1dDQIEnauHGjAoFAzO0yJX8T6bGTGEsaAAAUhMd2H5YxRrJ2SJsxRtt3H1bD9fN8iGyoQCAgSdq9e7fq6uq0bt06LVu2TBs2bNDy5cs1a9Ys7dmzR4FAQMuXL+9PmgfXW8faLlPys7RDksZMorQDAAAUhJb207IRkmhJstaqpf10hiMaXkdHhyRpy5YtAy4krK2t1a5du4bdPtntvJTHPdITpddf8jsKAACAtJs+blTMHunp40b5EFVkoQS6rq5OkpMAB4NBPfbYY6qqqlJbW5va2tqG3U+y23kpfxPpMROlk4zaAQAA8t/KRdXa8sv9EdustVq1qDrDEUW3e/duSdKiRYskOUPivfDCC/1D4m3fvj3m9qGh9RLdLh3yvLTjrYifzAAAAPLJzAlj1HjT5SoyUpFxloUeN950edZcaChJjY2NamxsVCAQ6B+VI3xc6VCPdehCwsGampqS2i4d8jeRHjtR6j0nnenwOxIAAIC0W7GoWk9/7j26sW66JOnGuul6+nPv0Yos6o0OjbyxYcMGSeovxQglweHLQkPehUo4Qstqa2vj2i4T8ri0Y5Jzf+ItadQ4f2MBAADIgJoJY7T6nTXauadFq99Zk/Ge6NDoGrt27VJtbW1/4tzR0aG2tjYtXrx4wOQqodkPGxoatGzZMgUCAe3YsUPr1q3TihUr+tepr69XQ0ODZs2a1X+B4XDbZYKJdoVnLjLGVEjq7OzsVMW5N6V/XiSt/g+p5hq/QwMAAMiI37Z26s+/8az+/a+v0YJplX6Hk3O6urpUWVkpSZXW2q5Y6+Zxj7Q7TTgXHAIAgDx3tOuMjnaflSTtO3piwL0kTSov06SKkb7Els/yN5EeWSkVlzqlHQAAAHns0ecP6WtPNQ9Y9untL/U//tS1c/SZZXMzHFX+y99E2hh3CDwmZQEAAPnto++YoWWXXhS1fVJ5WQajKRz5m0hLTiJ9gkQaAADkt0kVIynd8EH+Dn8nSWMnUSMNAACAtMjvRHrMJHqkAQAAkBYJl3YYY5ZKWmetHTJInzFmg/twsaSgtbZhUHt92NOAtXZzIu0JGztROvhMSrsAAAAAIom7R9oYU2eMaZS0QlJthPZGa+1m97ZCUq0xZkdYe72c5HirtXarpKC7v7jakzJ2stMjnUdjZQMAACA7xJ1IW2ub3B7mJwe3GWMCkpa69yGbJC03xoSS7gZJO8P2t1NSeA/0cO2JK58snT/DNOEAAADwnJc10rUa2FMdmui81k2wa621gyc/D7g93THbk46ofLJz330k6V0AAAAAkXiSSFtrO6y146y1TWGLQ0l1UBFKQVwdGpq
      "text/plain": [
       "<Figure size 800x494.438 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Obs[0.1795(32)], Obs[186(14)], Obs[0.578(70)], Obs[597(170)], Obs[1.42(83)], Obs[239(173)]]\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": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAFnCAYAAACW11IvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAr8UlEQVR4nO3deZwcZ2Hm8eftOSTNrRvbkl2YGxIOcwUCcQBDnBQBQhxgMSSEw7AhgUASKEOykHCVIRBzJRiWxQFMDpsAC8URjhgMBAg4a8DY2Ngu67Is6xiNZnqmr3r3j+qR5pZG091vVdfv+/nMR9PdNdIjkHueea8y1loBAAC4VHIdAAAAgEICAACco5AAAADnKCQAAMA5CgkAAHCOQgIAAJyjkAAAAOcoJAAywRgz5joDAHcoJACcMcZcYIyxxhgr6Ygx5jZjzLkrXP/VDsYD0EG9rgMAKLQxSY9sfj5urb19uQuNMRdJumDO4zFJz2k+vE/z93q9tXa8DTkBtBkjJABcu91ae/1JysiYpIUjJ5dJ+qG19sPW2tc3n7u6ef3rjDHXGGNuMcbcZYx5XVuSA2gZCgmAPHiOpA8veO5czRkxkXSbpAuMMZdJGpG0WdJ/SFov6VnN5wFklOHmegBarTm9cmnz4TskbVI6pbLZWvv6Oa+PSZpdFzIo6X7W2l9Z8Hud1/z0dklHrLVmmT/zCkn3V1pGNkl6iqTzJP2NpAdJmrDWjrbi7weg9VhDAqDlrLXXGGOkdArl+HSMMeYyY8wV1tqXz3k9stZ+vvn6l40xuyUdaP5W75D0TEk/kXTmcn9ecyHsBZKeLembkkattbc3y0y1edmIMeY8a+31Lf7rAmgBpmwAtMu40vUhc9eGvEPSJc0CMfv65+e8/npJOyR9SunoxjZr7Qutte+UdD9JMsY0jDEzxph/N8b8yBhzQNKPJb3cWnuDpCcvyNHf/HVCi9ehAMgICgmAdhqf+6C5A2ZcadlY6vUbmp/eX1Is6ZY5L/9R89eSpHWS3matfaSkj0makvTZ5uLVTQsybJX0NUkHl3gNQEYwZQPAlUFJDzPGnDtnSmdszuslSU+ds4bkPs1fK5J6rbXfbD5+h6TZXTSXSXrhnN/j3krXlPyepB+1/G8AoGUoJADaaWzug2bhGJN0vaRfllRdMKXz0Oav10l61JztvLPrRC5ROhqy3hhztaSXWWvHjTHjkkYlGUkPnPP7vUDSTc1rNkk63LK/GYCWYsoGQDudu2DU41JJH26WkClJGxa8/mFJ35e0XwvKjNK1JVK6SLWkdBHrUlMw/yVJxpiPSbpCUqP5/JjSnToAMogREgDtdLvSs0HGla4bOTR31KP5+vuNMdNKR0d+aq29yBhzgZplpjm6cYlObCPeJKlH0vuaO2nGdKK83GKt/Zwx5h5JVypdPzI7uiJ22ADZRSEB0E7j1tprmp9/bZnXX7jE89L8MjMm6R+UTvXMnlvyyeZhZ+dLKisdOXls84yTb0m6WNJk82uvlPSmNf5dALQRhQRAO42t4fVFZaY5cnKo+XV/pXRh7Kika621fnO05Oolfq9zJX3yVEMD6DwKCYCWaxaH1yuddrlswTTNSV9vGlvmt5+RVLHWPmfhC81txUue5Aog2zg6HkCmzCkrF0h654KdNhdI+melJ7H2OYoIoA0oJAByoVlG/l7NE1slfX/hfW8A5BeFBAAAOMc5JAAAwDkKCQAAcI5CAgAAnKOQAAAA5ziHBMCKvCAykjZL2qD0PWPex2PUo/dosFeSVXrfmErzY1rp/WqO7QifmLjIDiA/KCRAQXlBtF7SvSSdscSvcz/fphXeK4Zkrld6n5rl2Kte9Cc/2D99x1lK77Z7SNKBvsHfvrmn/34VSfsk7ZZ0p6Rdr/zQk2tr/bsByB8KCdDFvCAqKT234zxJj5D0MElnKy0bYx2KYSqN8iald+vdcfzJ0uDPJT1gwbXJB1/xjX2SftFfOfrDJ/znG+5697NL/++65K23/Dz8nT0dygvAAQoJ0CW8IOqT9BClxWNuARlymUuSqsnM4MLnTGl44xKXltQsLiVbH5D0mP1jpf+qHO5/tBdERyT9WNINcz5ujEN/po3RAXQIhQTIoeZ0y8N1onicp7SMrHMYa1lLFJJEZnDzSl8zMLW/bKXpfWMDdR2WJG1Uemff8+dcVveC6L+V3nzva5K+E4d+pXXJAXQKhQTICS+Izpb0dEm+pCcpXWSaedZaW0sqwwuePmxMactKXzdyLO6p9GmXalvLK1zWK+nRzY9LJZW9IPq20nLyVUk3xKHPcdRADlBIgIzygqhH0uOUFpCnS/olt4lO2zFJI/Of6j0iacVCMnr0jtFDIzqk8s7VHE8wIOlpzQ9JuscLoq+rWVDi0N+1it8LQAdRSIAM8YJoo6QLlRaQCyVtcpto7azspBYWEtN/7GRfNzy556yb7md+kpS9sTX88VslPa/5IS+IblVaTr4i6ctM7wDZQSEBHPOC6MFKC8jTJT1eUo/bRK1lbWPRlIsxG6ZP8kXj/bVjm2850/TXZ3ae1cI492t+/E9JR7wg+hdJ/xiH/vda+GcAOA0UEsABL4g2SXqBpD9Uuji1azVsY1H5MKWh+kpf09OY2SNp7I4zes3ksdFtbYq2UdIrJL3CC6KfS/pHSZ+IQ5/txYADFBKgQ5pngjxV0oslPVMZ3RHTanVbWzQtYkojS1163IaZQ+NWqsdbBhKddHKnJR4g6e2S3uoF0TckXSnp3+LQX3kkB0DLUEiANvOC6ExJlygtIjsdx+m4elKtLnzO9Iyu+N4zNLm3Ue/R7kZ9a6fPGClJuqD5MeEF0dVKp3Su63AOoHAoJECbeEF0vqQ/lvQsFfi/tVpSWTQ9Y0qjK25ZHpmI1x0e0t2a3uHyBqAjkl4i6SVeEN2udErnijj073aYCehahX2TBNrBC6JBSS+U9Erld5tuS1WTmUXngJjS6MJzSeYZmbhj88/PMHvt1L1H25dsVc6V9NeSAi+IPiLpnXHo73WcCegqFBKgBbwgGpX0Z5JeJSkr30QzoZJML1FIljw2PmVtMli+a+etZ+ru+syOrE1xbVD6//ErvCD6mKQwDv3YbSSgO1BIgDVojoi8StJfKN21gQUqjbJZ8FRdZmDZ81WMTfb1JPUdd9yrz05WRre2Od7p6pf0cqXTOVdJensc+rc4zgTkGoUEOA1eEK1TepZFIGm74ziZVmmUF7zPmEPGmGX/N+urTe630lm3bRuUdrc73Zr1SvoDSS/0guhfJb0tDv2fOs4E5BKFBFgFL4h6le6W+Suld6XFScw0yn3zn+kZ1wolbmD6wFRS0t4ZbV60OyfDSkpPg32uF0Sfk/SWOPSvd5wJyBUKCXAKmmeIPF/SmyXdx22afKkk5f55T5h1K54sMnzsTnN0QPs1vap72GSFUbqr6lleEH1J0t9wCixwaigkwAq8IDKSni3pbyQ92HGcXKo0yvO2+JrSwIpni4xOxMN7tpiJpOzlfXHwb0r6TS+IPiXpL+LQ3+c6EJBlefwJBOgIL4ieKumHkq4RZeS0VRvT8wuJGWqsdP3wsTu333qGTH16Z9Z22Jyu50v6uRdEgRdE/Se9GigoCgmwgBdEW7wg+qSkf5d0nus8eVdJZobmPjY9Kxwbb215/czhM+J79dopm9kdNqdjSNI7JN3oBZHvOgyQRRQSYA4viJ4v6WeSLnadpVvUkpl5DcSURvuWu7aUVHcbydy6fXDR2SVd4r6SvuAF0Re8IDrHdRggS1hDAkjygmiHpA9J4qfXFrLWlq3swNznTGlsYLnr11XGDydG9xzt27LitE4X8CX9uhdEb5J0eRz63f73BU6KERIUmhdExguiP1I6KkIZaTk7ufAZUxpZ9tj4oam7KpPrtdfxPWw6ZVDS30r6gRdETA2i8IrwHz2wJC+IHiDpm5I+KGnFe6vg9CRKphY+Z0rDm5e7fmTijr67NmnC5n+HzWqcp7SUvLt58i9QSBQSFI4XRL1eEL1B0g2Snug6TzdLbKO84KmqKW0YW+76kYl44y/OMEkX7bA5VT2SXivpJ14QPcZ1GMA
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAFmCAYAAABdi4GKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAu4ElEQVR4nO3deXhcZ2Hv8d87I8uWV9mWJSfOcmwndkggDiHsUJYE2pu5hRIoO5RCCZSwtCzJgVtuafu0TCgEaFkCBS5QdsJ9CpfD0qQUQgKBBDsbieMkzmRxnFiOrHVG0izn/nGObFmakUbSzLznzPl+nkdPNDNHo1+eOPJP7/ue9zW+7wsAAMCmlO0AAAAAFBIAAGAdhQQAAFhHIQEAANZRSAAAgHUUEgAAYB2FBAAAWEchARA5xphu2xkAtBaFBEAkGGMuMMb4xhhf0hFjzL3GmG22cwFojQ7bAQAg1C3pSeHng77v76/3C40xV/u+/4KmpALQEhQSAFGy3/f9wYV8gTHmZZIumOeabkkvDx9uV1B+Llvo9wLQPEzZAIitsGjUM61zuaSbfN//vO/7l4XPfbdpwQAsGCMkAJoqHMF4f/jww5I2KBih2DitHEx5uTFmIPz8yVVen+nlkj6voHDMZZuCUZTd4eN7JV08I+el4acZSX2+758xz3sCaCAKCYCm8n3/KmOMFIxI7J5aG2KMudwY8znf998SXrpf0lZJU49vNMZcK+mHqlJejDHnSrqpzgwz15dsl3TNtPe6XNK3Jb1CkpF0sjHmJknf8X3/I3X/ywJYNKZsALTCoIL1IdMXqn5Y0sVTd9L4vr/b9/1/Dp8/V9I3JT1bwQiIjDGfm/Ge5ykY+fh5+PrLjDEXG2MuDQtGVeH3u0Bh8QmnfS5QUIjulXR/+PEkSaNzvReAxqGQAGiVwekPwgWlgwrKx8zr9vu+/4vw8TbNKC/GmJeF60GuknRFeN3u8LmPhNfMLDAyxlysYKTmLTPK0bbw4zJJv5Y0GT6/TzOmdgA0B4UEgHXGmG5jzJFp+44MTt8cbXp5CZ8fnPblQ+E11UZf3hbuZzJmjLlXQfm5T9J3jDH/PfXevu+vVzBCsk3SIUmd4fvsl9QdTg8BaCIKCYBW6Z7+ICwW3Tq20PSmKqMW8n1/t473ckkvCKdmLpX05vD9Lg0X0E4vMIcUjHqsDN/vmZK+pWBq5rnGmDtmfr/QJknXhHkGVd+dPACWgEWtAFplmzGme9reH++X9Plpi1yvnnZtd/j6ZeFr3eFzu2eUFhljHpaUmWPxaVHSuKRhSY9XUID+QtLNkh5njPmG7/uvVnD3jxQsrF0r6U/DxwPTXgPQJIyQAGiV/ZIuCLeIv1TSY9PusJHv+x8Jn3+5ghGJ66aVjOPKywyrZz4xY/TlBAVTMHkdW1ibDi8dl/SqGVvUv1bSnWyaBrQWIyQAWmUwXIQqTbvldrqwlFwg6XxJB8LPz1VQXmbtSRIuUp2asvkPSZ/yff8aTSswxph7FKwzKUoaCb/0FZLK4cdY+D2mys7nJL1p2rfZoGCUBEATUUgAtEr3Aq6dt7xIku/7nzfG7FdQIr4mHd3gbGaBmVAwUrI2vI23V8EoySoFUzlScIuxJI2G77MtzNytY2UFQJMwZQOgqcJRjssUrCGpd0+P7gV+m0Hf96/yff8a3/c/UmU9yWYFUzuFsKgcknSLpC8rWC/SK2m5pIOSnhN+/8sUFpEqC2sBNBgjJACaKpxCqTnKMd3M8lLH1vFTuud5/aCkFZIcY8yXFdxF8w1JGxWUkk9Pu/aN4T8vVnB3Tr0ZACyB8X3fdgYAWJRpBeYCSR+psc7kAkmfCR+mw+sHFW7INnM0JZzyGVTt83YANAGFBEBbCwvJV8OHB33ff5LNPACqY8oGQNsKy8jlCha0SlKXxTgA5sCiVgBtK1zk+iRJL1CwjqWbw/KAaGLKBgAAWMcICQAAsI41JABqclxvnYI9PDYr2K9jmYKfG0c/jJT6pdamJfmSKpImJRXCj7yCjcZGFeyWevik7LMnWv3vASD6mLIBEsZxvWWS+hQs9Nw87WPm482qYxFoWir/Iigk9Rr7/ZHrr7598LoTJT0i6YCkhyQ9KOl+BZuRHXjPt3/IDycgQRghAdqY43qbJT1xxsc2ScZirFX50vAJkp5S64Jlq/7kxk+/9Wddku6WdJek2yXdJmnvJVc+f7I1MQG0EoUEaAOO6xlJ2zW7fPTZzFXLWHl41Vyvpzp6HQW7qT5+xkvFT7/1Z3dK2qPgJN8bJe2+5MrnMw0ExBxTNkAMhSMfL5B0noLisUvBGo+WW8SUjX700BfuHyk+dmqNl0dXrH/36gW83aSk3aff/d0fnHzg53f+8ixz7cXfu4PTeYGYYYQEiAHH9VIKpjgulJRRUEJsTrssyXh5dH3tVzselrRjAW/XKelp64/c1Sdp6/efavZ+/P3fPlzwV98g6WpJv8xlM4Wl5AXQfBQSIKIc11sv6Y8UlJA/ktRjN1Fj+L4/UaxM1BzNMalVCx/d8CuHV+cPbq0YHXqoZ8V4YWD1syQ9S9J7JU04rnedgo3R/lPSnlw2w9AwEDEUEiBCHNfbpWAE5EJJT1NwGFxb8eUP6NhW7rOY1PoFrwdZMX7kXkk9j3Tr3lT+lPyMl5dLOj/8+LCknON635D09Vw2c8dCvxeA5qCQABaFi1GfL+kVCkrIFruJmq/il4Y0VyHp6F3who0bjuwdl6TfnW6KlaGz51wwK8mR9AFJH3Bc72ZJX5f0zVw2c2Ch3xdA41BIAAsc19sq6Q2S/kxSrcWdbalYmRyb6/VUum8hC1olSb39uzdK0vWPS/WMD5914gK+9Jzw43LH9X6hoJxclctmhhaaAcDSUEiAFnFcb5Wkl0n6c0l/oBgvSl2KiUphfK7XTXrTpgW9oe+Pdw/evaMiPZbbtHx8YmhV9yJipSQ9L/z4tON6noJy4uWyGW4pBlqAQgI0meN6Z0p6m6TXydKtuVEyXh4rzfHypEmtW8gIhzpK+X0pv3z2o9262xROHV1iPClYc3JR+DHouN73JH02l838rgHvDaAGCgnQBI7rdUh6iaRLJD3HcpxIyZeG53g1fcAYs3Uh77duaP+AJO3ebib8oV1rlhRutm5Jb5L0Jsf1/kvS5bls5uoGfw8AopAADeW4Xo+kt0u6WHMs3EyyfGm49p1DpuuwpAUVkt7+Pask6bozUxvGR888ZWnp5nS+pPMd19st6XJJ38tlM+Umfj8gUSgkQAOEe4a8V9I7JS14UWaSjJWGltd6zaS7F7aBme/7PY/ddrovDe7vXTExMbJy3ZIDzu9cSd+WdK/jeh+V9OVcNjPnuhgA86OQAEvguN46SX8dfiR+fUg98qXhlbVeS6U3LWjDslSleO+yUv60/rX6rcadkaWnW5Dtkj4r6UOO631S0me4OwdYPAoJsAiO662W9C5J75E0xzbomClfGqk5ipFK9823h8hxVo8deFjSaXu2m0Jl6JxGrx+pV5+kf5LkOq73OUkfz2UzBy1lAWKLQgIsgON6KxWsEXmf2mQr91YrlEc31HrNdPRuXMh79Ry+NS1J15+ZWjc+dsb2pWZborUK/ly8y3G9r0j637ls5hHLmYDYoJAAdXBcb4Wkv5R0mYLfiLEIvu/ny36x1pRN2aTWL2in2t7+m0/xpaF9m7uKk/d22RohmalT0pslvcJxvX+Q9MlcNlO0nAmIvAVv0QwkieN6nY7rXSLpXklXiDKyJL4qR2q/mjpoTLqz/jerPLqycOjkgTXap8LWKK7dWCvpnyXd5rjeH9oOA0QdhQSowXG950q6XdKnJC1osy5UVw7OsanOLO9fyHt1FQ7fJ0k3bzNjlaFdrbi7ZrF2SvqJ43o/cFzP9rQSEFlM2QAzOK63QdJHFWzxjgYqViZqnmNjUusWdJfMhoE7JiTpV48zawpjZ5yx1Gwt8MeSXui43hWS/jGXzcx5pg+QNIyQANM4rvdKSXeKMtIUE+VCzXNhUumeBd3y29e/p9eXRveesKpU1IoF3Z1j0XJJ75d0l+N6r7YdBogSRkgASY7rnSrpM5I
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAFaCAYAAAAuM0ZcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzrElEQVR4nO3deXiU5aH+8e+bjbBJAMENdQRxBxFwQQnWiHaZ2Pb00KY9dm+Ptub0tD3+qmP3emodwNj22Gqk1iNWW7cea3VsC+64CwFBFlFw2MISlkDINtvz+2MmNoYAWWbmeWfe+3Ndc8Esmbm9rhjuPO+zOMYYRERERGwqsB1ARERERIVERERErFMhEREREetUSERERMQ6FRIRERGxToVERERErFMhEREREetUSERE+sFxnDLbGUTygQqJiEgvOY4z03Ec4ziOAfY4jrPOcZyxtnOJ5LIi2wFERHJQGTAl9fdGY8x6i1lE8oIKiYhI36w3xjTaDiGSL3TJRkRERKzTCImISIrjOLOAG1J3bwZGkLw8M9IYc32Xl3/GcZzdqb+f283zItILjk77FRH5p1QpeRgY1zE3xHGc2UCZMebq1P3JAMaYutT9q4ApHc+LSO+pkIiIdOI4zkzgTmPMuE6PlQF76FRSunzNWGAdMFzzSkT6RnNIREQO1Nj5TqpkNAKTu3txp5Kipb8ifaRCIiLSC47jlDmOs6fzviPaHE2k/1RIREQOVNb5TqpwlAF1qYcWd7l0Mxb+OadERHpPhURE5EBju4x63ADMM8Z07D2ysMvrbwC0ykakH7TsV0TkQOuBmY7jNJKcN7Kr87JeY8wcx3GuS90dByw0xszLfkyR/KFVNiIinaRW2cw2xkw57ItFJG10yUZE5EBltgOIeI0KiYhISmp05HqSc0hm284j4iW6ZCMiB/AFQgXAEcCwg9wGA4Ukf6kpvJ7SxisoGQbEgHinWyuwC9jZ6bZrTLA8ltX/IBFxPRUSEY/xBUJFwBjghG5ux6eeGwY4PX3PuQx8fhrFF/fw5QbYywdLyk5gB8ndTtcAq8cEyxt6+vkikvu0ykYkT/kCoULgZGBCl9tY7F6udUjO0Sgjma9bmwOLdgFvkywoa4DVqT/fGxMsj2c8pYhklUZIRPJAqnxMAqYD55AsHmcApdn4/F6OkPRXO/AusARYBCwaEyx/O0ufLSIZokIikoN8gdBA4AKSBaQcmAYMsZUny4WkOzuAF1O3RcBSjaKI5BYVEpEckJr3UQ58FJhBcrOuYquhOnFBIelqP/AKqREU4NUxwfI2u5FE5FBUSERcyhcIjQA+BlQCH8bFe2O4sJB01Uxyu/fHgCfGBMt3Ws4jIl2okIi4iC8QOh24InWbRnJprevlQCHpLE5y9OQvwCNjguUb7MYREVAhEbHOFwgdB/wb8AWSk1FzTo4Vkq5eBx4CHhoTLN9kO4yIV6mQiFjgC4SGAP9KsoRcQo7vmpzjhaSDAV4D7gXuGxMsb7KcR8RTVEhEsiS1++nlJEvIJ4FBVgOlUZ4Uks6agPuBO8YEy5fbDiPiBSokIhnmC4SGAV8D/gM4yXKcjMjDQtLZy8DtJOebtNsOI5KvVEhEMiQ1QfVbwBdJnv2St/K8kHTYCdwN3DkmWL7edhiRfKNCIpJGqcsyfpJFZCa9OA8ml3mkkHQwwD9Ijpo8MSZYrh+iImmgQiKSBr5AqBj4Esmj6w96Pku+8lgh6Ww58DPgURUTkf7J6Zn9Irb5AqEBvkDoGpJnq/wOD5YRj5sI/Bmo2xxY9EnLWURymkZIRPrAFwiVAF8Hvg8cZzmOdR4eIemqDvjpmGD547aDiOQaFRKRXkidKfMl4EfAiZbjuIYKyQEWkywmIdtBRHKFLtmI9JAvEKoEVgJ3oTIihzYVeGJzYNFrmwOLPmo7jEgu0AiJyGGklu/+kuQBd9INjZAc1kLgP8YEy9faDiLiViokIgfhC4SGAz8FrgGK7KZxNxWSHokAc4GbxgTLW22HEXEbXbIR6cIXCBWmVs68A/wnKiOSHiXAD4BVmwOLPm47jIjbqJCIdOILhC4AlgK/BUZajiP5yQc8tjmw6PHNgUV5eZSASF/oko0I4AuEBgI3Ad9GRb3XdMmmz1qBm4E5OidHvE4/eMXzfIHQxSR33Pwu+n9CsmsgcCOwYnNg0eW2w4jYpBES8SxfIDQEmAN8A4+cOZMpGiFJm7uA74wJljfbDiKSbfptUDzJFwhdDrwFfBOVEXGPrwPLNgcWXWA7iEi2afWAeIovEBoA3AL8h+0sIgdxMvDi5sCiXwA3jgmWx2wHEskGjZCIZ/gCoZOBV1AZEfcrBH5Ut+upe2qqKsfaDiOSDSok4gm+QOhzJA8+O8d2FpGe2N2+7cV39i25ElhaU1VZZTuPSKbpko3ktdRy3tuAr9nOItJTcRN775mtf+woz0cAD9RUVV4GfOvaB5/QLq+SlzRCInnLFwidAbyOyojkEGNM5LmtD7TFTXRwl6e+BrxUU1V5go1cIpmmQiJ5yRcI/SvwBnCW7SwivbG+6c1XdrZvOf0gT58DvFFTVTk9m5lEskGFRPKOLxD6IfAwMMh2FpHeaI01L1m86x8zDvOy0cAzNVWVV2Ujk0i2aA6J5A1fIFQK/B74N9tZRHrLGLNrYf09Y+jZvjjFwJ01VZVnA9++9sEntDRYcp5GSCQv+AKho4BnURmRHPXGzr+ta43vP6qXX3YNsLCmqvLITGQSySYVEsl5vkDobJKTV7W7peSkhrbNL7y3f8V5ffzyD5GcV3JmGiOJZJ0KieQ0XyBUCbwIaOWB5KR4Ivruc9se6GsZ6eADFtVUVU5LQyQRK1RIJGf5AqHPA48CQ2xnEekLY0zb01v/aBImXpqGtxsOPFVTVfnhNLyXSNapkEhO8gVC1wD3oonZksPe2bfk9T2RbePT+JaDgMe1s6vkIhUSyTm+QOj7wG/RKb2Sw1pi+95Yuvvpwy3x7Yti4I81VZXfyMB7i2SMConkFF8gNBu4yXYOkf4wJrFjwZb5J2XwIwqAO2qqKn+Ywc8QSSsVEskJvkCowBcI3QFcZzuLSH8YY8wrDY9vak+0ZGOp7n/XVFXemoXPEek3FRJxPV8gVADMBzQELTlve9uGFzY1r5mSxY/8bk1V5S+z+HkifaJCIq7mC4QcYB7wedtZRPorloi8vWjbwzaW5n6npqryvy18rkiPqZCI292KTuuVPGCMaXmq/r7iBIkSSxF+WFNVGbD02SKHpUIiruULhG4EvmM7h0g6rN776pK90YaxlmPcXFNV+Z+WM4h0S4VEXMkXCP0/4Ee2c4ikw/5o46sr9rxQbjtHyq9qqio16iiuo0IiruMLhK4G5trOIZIOCZPYurB+/qm2c3TiAPNqqip1EKW4igqJuIovEPo34HbbOUTSwRiTeGnHo9sjibbhtrN0UQDMr6mq9NsOItJBhURcwxcIlQP/i74vJU/Ut767qL7l3Um2cxxEEfBATVXlRNtBREA/+MUlfIHQWJIH5dlagSCSVtFE+8qXtj96ke0chzGE5Nk3R9sOIqJCItb5AqFhwBPASNtZRNLBGNO0sP7eIQaTC4c/ngA8VlNVmY4Th0X6TIVErPIFQkXAQ8DptrOIpMuKPYvebIruPtF2jl44j+ScEh1YKdaokIht/wNcbjuESLrsi+x6efXeV6bbztEHnwF+ZjuEeJcKiVjjC4S+BXzTdg6RdEmY+Oantv7hTNs5+uFHNVWVV9oOId6kQiJW+AKhCkAHfkneMMbEX9j+yJ5oon2Y7Sz99PuaqsqptkOI96iQSNb5AqFjgT8BhbaziKTLpuY1L25vDU+wnSMNBgAP1lRV5nqxkhyjQiJZlZrE+gAw2nYWkXSJxFuXv9rweC7OGzmYscBdtkOIt6iQSLbdBLjlTA+RfjPG7F1QP3+kweTbiN+smqrKa2yHEO9QIZGs8QVCHwG+ZzuHSDot2/3MyubY3uNs58iQW2uqKifZDiHeoEIiWeELhI4B7iV5sJdIXtjTvv3FtfsWX2g7RwYNAB6qqaocajuI5D8VEsk4XyBUAPwBGGU7i0i6xE0
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAFdCAYAAAAzNnbkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA27UlEQVR4nO3deXzU9Z0/8Ndnrpwk4RYEGbkREQUUUDmM2NZOvI944NXdqm1qt7vZ6rj9dbfddbejbbrbbtlG7N2KjdraasYLoUJAlPu+jwFCOAI5SMg98/n9MRMdYoAcM/P+fuf7ej4eecBMJt/vyzZMXvl8P5/PV2mtQURERCTJJh2AiIiIiIWEiIiIxLGQEBERkTgWEiIiIhLHQkJERETiWEiIiIhIHAsJERERiWMhISIiInEsJERERCSOhYSIiIjEsZAQEcWZUipHOgOR0bGQEBHFgVJqnlJKK6U0gGql1D6l1EjpXERG5ZAOQESUpHIATI38vUZrvV8wC5HhsZAQEcXPfq11jXQIIjPgJRsiIiISxxESIqJuUErdDeDZyMMfAOiH8OWZ/lrrZzq8/F6lVFXk71d38nkiilBaa+kMRESmEiklrwEY1T43RCn1PIAcrfUTkcdTAEBrvT7y+HEAU9s/T0RnYyEhIuompdQ8AC9qrUdFPZcDoBpRJaXD14wEsA9AX84rIfo8ziEhIuqZmugHkZJRA2BKZy+OKilc+kvUCRYSIqIYU0rlKKWqo/cd4eZoROfHQkJE1DM50Q8ihSMHwPrIU2s7XLoZCXw2p4SIzsZCQkTUMyM7jHo8C2Ch1rp975HFHV7/LACusiE6By77JSLqmf0A5imlahCeN3Iqelmv1voFpdTTkYejACzWWi9MfEwic+AqGyLqlNvrVwD6AhgEYGDkz0EA+gNIQfgXGnvkT8dU2E//BBk5ANo6+WgGUAngOIBj7X8O881qSeB/UsxEVtk8r7WeesEXE1GXcISEyILcXr8DwFgAEwGMwWdlI7p8DEA33iNcUJsBXNGdHOXeshp0KClRfx4DsAvA/mG+WUb8zSlHOgBRMmEhIUpibq/fhvDlgokALo/6cywAl2C0djmRj3HneU1dubdsC4BNADZG/twyzDerId7hOhMZHXkG4Tkkz3P3VaLY4CUboiTh9vpdAGZEPiYhXD7GA0hLxPlnwrH5h0jv1ghJL4QA7MVnBWUTgE3DfLPKE3R+IooxFhIik3J7/XYA0wDkRj6uBZAulSfBheRcTgFYCeBdAO8O8806IJyHiLqIhYTIJCKTTCfjswIyC0CWaKgoBikkHe1BpJwA+FDqMg8RXRgLCZGBub3+IQBuB3AjgLkIr3AxJIMWkmjNAMrw2ejJNuE8RBSFhYTIYNxefz8AdwG4H8AcmGQDQxMUko7KAbwH4B2EC8oZ4TxElsZCQmQAbq8/A8BtCJeQLwJwyibqPhMWkmj1AN4A8HsAS4b5ZoWE8xBZDgsJkZDIqpibES4ht0BwQmosmLyQRDsCYBGA3w3zzdoqHYbIKlhIiBLM7fXPBvAwwpdlcmTTxE4SFZJo6wC8BGDRMN+sOukwRMmMhYQoAdxevxNAPoB/RPi+J0knSQtJu3oArwBYOMw3a610GKJkxEJCFEdur78vgCcBfAPAUOE4cZXkhSTaegA/R/iSjinvxUNkRCwkRHHg9vpHIzwa8giADOE4CWGhQtLuMIDnAfximG9Ws3QYIrNjISGKIbfXPwfAPwHIg0mW68aKBQtJuwoAPwTw4jDfrEbpMERmxUJC1EuRHVTzAfwzAMvejt7ChaTdcQBFAP6Pe5oQdR8LCVEvuL3+eQBeAHCVdBZpLCSfOgngvwH8L1fmEHUdCwlRD7i9/kkID9N/UTqLUbCQfE41gJ8A+Mkw36wa4SxEhsdCQtQNbq//YgDPIbyPiKXmiFwIC8k51SI8YvIC55gQnRsLCVEXuL3+LABeAN8CkCabxphYSC7oAICnhvlm+aWDEBkRCwnReUQ2NHsSwL8CGCAcx9BYSLrsrwC+Ocw365B0ECIj4ZAz0Tm4vf5bAGwD8FOwjFDs3AZgR7m3zFvuLTPdTRSJ4oWFhKgDt9ff3+31LwLwJoAx0nkoKaUD+AGATeXeshukwxAZAQsJURS3138nwqMi90tnIUuYAGBpubfs5XJv2UXSYYgkcQ4JEQC31z8AwAIA90pnMSvOIem1WgDfRXhjtaB0GKJE4wgJWZ7b678HwHawjJCsbITnK60p95ZNlA5DlGgsJGRZbq9/kNvrfw3AqwAGSuchirgK4VLy99JBiBKJhYQsye3134fwXJG7pbMQdSINwEvl3rJF5d6yPtJhiBKBhYQsxe31Z0dGRV4Bl/KS8d0PYH25t8zy90qi5MdCQpbh9vqvBLAeHBUhcxkNYFW5t+wp6SBE8cRCQpbg9vq/AmAVgJHSWYh6IAXAT8u9ZX8u95blSIchigcWEkpqbq8/1e31/xLALwGkSuch6qU7AGwo95ZNlw5CFGssJJS03F7/JQA+AvAV6SxEMeQGUFbuLft2ubdMSYchihUWEkpKbq9/FoC1CC+hJEo2TgAvAHiLq3AoWbCQUNJxe/1PAFgC7i1Cyc8DYEW5t2yYdBCi3mIhoaTh9vqdbq//5wCKEf4NksgKrgDwcbm37ErpIES9wUJCScHt9WcA8AN4UjoLkYCLEZ5XcrN0EKKeYiEh03N7/TkAFgO4STgKkaRMAG8Fnl7CLefJlFhIyNTcXv9gAB8CmCkchUic1rpp+fHX/7EoP+/70lmIuouFhEwrsqy3DMBk6SxE0rTWbZ+c9G+vbDp8GYB/LcrPKy7Kz+N7PJkGv1nJlNxe/1gAKwCMkc5CZATbalZ+fLB+29VRTz0BoKQoP48TvMkUWEjIdNxe/2SER0aGS2chMoKD9ds/3Faz8vpOPnU3gFdZSsgMWEjIVNxe/0yE54wMEo5CZAgnm44s/7jyrbnnecntAF5jKSGjYyEh03B7/fMQXk2TIxyFyBDOtNZ8suToHzobGenoNgCvF+XnueKdiainWEjIFNxe/y0ASgFkSGchMoKWYNPmt4/8YjK6/j5+K4A/sZSQUbGQkOG5vf7rAbyK8C3YiSwvGGrb5y9/cXhIB7t7B+s8AH8uys/jvyUyHBYSMjS31z8RwJsAuvvGS5SUQjp07J0jv0htCTX17eEhPABe5pJgMhp+Q5Jhub3+4QDeBdDTN16ipKK1rl1S8YfTZ9pqL+7loe4C8LNYZCKKFRYSMiS3198X4TLCu5gSAdBat6w88cb+qpajY2N0yK8V5ed9N0bHIuo1FhIyHLfXnwbgLQCXSWchMgKttd5YtXTtkYY9V8X40P9elJ/31Rgfk6hHWEjIUNxevx3AKwCuk85CZBT76jYs33167bVxOvzPi/LzbovTsYm6jIWEjGYBwnsmEBGAY40Hlq07tXhOHE9hB/DHovy8ruxnQhQ3LCRkGG6v/98Qvv8GEQGobTm5ctmxV2cn4FSpCC8HHpGAcxF1ioWEDMHt9X8VwPekcxAZRVPwzPr3j/z6agAqQaccCOAvRfl56Qk6H9FZWEhInNvrn47wpRoiAtAWatnlP7xwTAihRO+qeiWAXyf4nEQAWEhImNvr7wegBABv/EUEIKSD5W+Xv5TTplv6CEW4tyg/71+Ezk0WxkJCYtxevwLwOwC8bk0EQGtd9f6R37Q2BusHC0f5j6L8vDzhDGQxLCQk6RmEt7EmsjytdeOyY68eqW09eal0FoR/NrxclJ83XjoIWQcLCYlwe/2zATwnnYPICLTWwTUn39lyvCkwSTpLlCwArxbl5/E+UpQQLCSUcG6vfxDCm5/ZpbMQGcGO2o8/OlC/5RrpHJ2YBKBIOgRZAwsJJZTb67cBWARgqHQWIiM4fGbXh1uql8+SznEeXy/Kz7tdOgQlPxYSSrR/A3CjdAgiI6hqPlr20Ym/zJXO0QW/LMrP440uKa5YSChh3F7/TQD+n3QOIiNoaDu9+oOK38fr/jSx1g/hSa68zEpxw0JCCeH2+gcDeBn8niNCa6h529vlL12uoc30A342+AsFxRF/OFCi/BT
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAFkCAYAAAAQQyCBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvUUlEQVR4nO3deZxcVZ338e+pysISpNlBo5YJOrhjwN0MY4trNc4zirTOwoyyObaPOrZKOTPquIwUzrT7aIOCjss80wKC0CViQAUXCAlhR9ZQIYqQhXT2dHfde54/blXSaTpJL7fqd2/dz/v16ldSne66X14J1d8659xznPdeAAAAlnLWAQAAACgkAADAHIUEAACYo5AAAABzFBIAAGCOQgIAAMxRSAAAgDkKCQBIcs51WGcAsoxCAiCznHMnO+e8c85L2uCce8g5t8A6F5BFs6wDAIChDkkn1H8/5L1faZgFyDQKCYCsW+m9H7IOAWQdUzYAAMAcIyQA2o5z7lRJH68/PE/SoYqmZw7z3p877stPc849Uf/9Syf4cwAt4DjtF0A7qpeSSyQtbKwNcc6dL6nDe39O/fEiSfLer6g/PlvSCY0/B9A6FBIAbck5d7KkC7z3C8d8rkPSBo0pKeO+Z4GkhyQdwroSoLVYQwKgnQ2NfVAvGUOSFk30xWNKCrf+Ai1GIQGQSc65DufchrH7jrA5GmCHQgKgnXWMfVAvHB2SVtQ/tXzc1M0CadeaEgCtQyEB0M4WjBv1+LikC733jb1Hloz7+o9L4i4bwAC3/QJoZyslneycG1K0bmT92Nt6vfdfcM59rP5woaQl3vsLWx8TAHfZAGhL9btszvfen7DPLwZgjikbAO2swzoAgMmhkABoO/XRkXMVrSE53zoPgH1jygYAAJhjhAQAAJijkAAAAHPc9gtgN4VS5VBJx0h6av3XYyQdqOj1Ypak2WN+P+tqHbT1ILkDJIWSAkmjkjZL2lj/2DTu142SHp9fXjzawv8sAAnHGhIgI+pF42navWg8ddyvR0vabyrP+ysdtHqW3NOnGCeU9Lik1ZIeqf+6eszjh+aXF6+f4nMCSDEKCdCGCqXKUyWdMO7jmGZca5qFZDIel3RX/ePO+q93zy8v3tKEawEwRiEBUq5QqjxNTy4fR7fq+k0sJBPxklZJul3SUkk3Slo2v7x4a4uuD6BJKCRAitSnXV4jo/IxkRYXkokEku6Q9BtJN0i6YX558RrDPACmgUICJFyhVDlW0l9KequkV0vK2ybaXQIKyUTulXSNpIqk6+eXF48Y5wGwDxQSIGEKpUpO0isVFZC3SjrONtHeJbSQjLVF0rWKyslP55cXP2qcB8AEKCRAAhRKlQMlvUFRASlKOsI20eSloJCMd6uicnLJ/PLiO6zDAIhQSAAj9TthTlFUQjo1xdttkyKFhWSsOyX9UNL/zC8vXm0dBsgyCgnQQoVSZa6kt0s6S9JJkpxtoplLeSFp8JJ+raicXDK/vHiDcR4gcygkQAsUSpXnKiohp0s6zDhOrNqkkIw1Iumnki6U9LP55cW8SAItQCEBmqQ+GvIOSecoulW3LbVhIRnrQUnfkPSd+eXFQ8ZZgLZGIQFiVihVjpb0PkVF5EjjOE3X5oWkYaui6Zyvzy8vvtM6DNCOKCRATAqlygmSPiTpNElzbNO0TkYKyVg3SPqqpMvnlxeH1mGAdkEhAWaoUKr8laReRZuWZU4GC0nDfZLOl/QDTi4GZo5CAkxToVR5raIfSC+1zmIpw4Wk4ZHbn/jlJ+/dePP/9g4MDluHAdKKQgJMUaFUOV5SWdIbjaMkQtYLSejDxy6rfvGQUMFaSedJ+nbvwCBb1QNTRCEBJqlQqjxL0uckvUttsH9IXLJeSB7YdMv1K9Zfe9KYT62W9FlJF/UODLLGBJgkCgmwD4VS5QhJn1B010xmFqtOVpYLiffh2stWfWle4Gv7T/DHd0j6cO/A4HWtzgWkEYUE2INCqTJP0WLVXkkHGcdJrCwXkoc23Xb98vXXnLSPL7tK0kd6Bwbvb0UmIK0oJMA4hVJltqLRkE8oA/uIzFRWC4n3fv1lq760X+BHD5zEl48q2mDt070Dg2xLD0wgZx0ASJL6Lby/l/Q1UUawF9Utd901yTIiSbMlfVDSg33dXR/o6+6a1cRoQCoxQgJIKpQqh0n6L0nd1lnSJosjJN77DT9e9eVZNT8y3am8WyWd0TsweGucuYA0Y4QEmVcoVf6PpLtFGcEkrd567x0zKCOS9BJJN/d1d5X7urv2iysXkGaMkCCzCqXKoYqmZv7aOkuaZW2ExHu/8fJHvqLRcPjgmJ7yfkln9g4M/jqm5wNSiRESZFKhVDlF0l2ijGCK/rjtgdtiLCOS9BxJ1/d1d32jr7uLu7mQWYyQIFMKpUqHpK9IOt04StvI0giJ937zFY98tTYS7jikSZdYLems3oHBa5r0/EBiMUKCzCiUKm9WNCpCGcG0/Gn7Q7c0sYxI0tMlXd3X3fWffd1ds5t4HSBxGCFB2yuUKk+R9CVJ77HO0o6yMkLivd/6k0e+vn043HZ4iy65TNI7ewcGV7boeoApRkjQ1gqlSqeiURHKCGbk8R3V5S0sI1J0ivStfd1dp7XwmoAZRkjQtgqlyoclfUFS3jpLO8vCCIn3fvuVq7+xeUewxWqzvG9J+mDvwOB2o+sDTccICdpOoVSZUyhVLpbUJ8oIYrB2x+qbDcuIJJ0laVlfd9fzDDMATUUhQVsplCpHSvqFpHdbZ0F78N4PL107+BzrHJKeL2lpX3fXW62DAM1AIUHbKJQqxytaCPhq4yhoI+uHH126Ldh8jHWOunmSLu/r7jrXOggQN9aQoC0USpW3SfqepMkedoaYtPMaEu/9aOUPF6zZWtv4NOssE/iepLN7BwaHrYMAcWCEBKlWKFVcoVT5lKRLRRlBzDaMPHZTQsuIFO2n88u+7q6jrIMAcaCQILUKpcoBkgYk/ZskZ5sG7cZ7X7txzVUF6xz78EpFh/S92DoIMFMUEqRSoVSZL+nXkt5hnQXtaePI2qVbahvSMBX1DEm/7evuKloHAWaCQoLUKZQqr1C0eHWRdRa0J+99cOPaK5M6VTORAyVd0dfd9S7rIMB0UUiQKoVS5U2SfinpaOssaF+bRtcv3TS6vmCdY4pmSfpBX3fXe62DANNBIUFqFEqVUyRdIWk/4yhoY9778Ka1V6V1oWhO0jf7urtK1kGAqaKQIBUKpcrbJV0maa51FrS3LbUNS4dG1iy0zjFD5/V1d5WtQwBTwT4kSLxCqfIuRXsuzLLOgidrp31IvPd+yaPfe3DDyGPPts4Sk35JPb0Dg6F1EGBfGCFBohVKlb+X9ANRRtACW2sbl7ZRGZGk90r6Xl93F6/1SDz+kSKxCqXK6ZK+I/6dokVuWjvYYZ2hCf5G0gV93V3s1YNE44UeiVQoVU6TdLHY8Awtsq226eb1w388zjpHk5wp6YvWIYC9oZAgcQqlylsVTdPkrbMgO5aurbT70QMf6uvu+ox1CGBPKCRIlEKp8gZJP5I02zoLsmN7bcsta3Y88nzrHC3wib7uro9ahwAmQiFBYhRKlZMU7TPCrb1oqaXrKnOsM7TQF/q6u/7ROgQwHoUEiVAoVV4iaVDS/tZZkC07gq23Pr69+kLrHC32X33dXX9rHQIYi0ICc4VS5UhJP5E0zzoLsmfZup9lceG0k3RxX3fXX1gHARooJDBVKFVmK9qBtS021kK6DAfbb39024PHW+cwMlvSZX3dXe207wpSjEICa1+X9BrrEMim5euuCawzGDtU0mBfd9ch1kEACgnMFEqV90k62zoHsmkk2HHnH7bdt8g6RwI8R9Ilfd1d7IYMUxQSmKjfUfNl6xzIrhXrlwxbZ0iQ1ykarQTMUEjQcoVS5ZmSLhV7jcDIaDh8z6qt95xonSNhzunr7vqQdQhkF4UELVUoVQ5UdEfN4dZZkF23rv/FFusMCdXX1931ZusQyCYKCVrtu5JebB0C2VULR+57eMsdL7POkVA5ST/o6+56hnUQZA+FBC1TKFX+VdKp1jmQbbc/8asN1hkS7lB
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAFbCAYAAADlb5X5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA12klEQVR4nO3deXyU5b028OueLQkhC7tggBFQAUGQzQWCNOI6aWvdotXaqmjtSdvTc3Jax/a8fdue89bBNj3dbCO2avdGbevRjCtSIKKCyL6vAwkBErJA9mXmfv+YiQwxCVlm5vc8M9f388mHPMNknsuPyeTivu/nuZXWGkRERESSLNIBiIiIiFhIiIiISBwLCREREYljISEiIiJxLCREREQkjoWEiIiIxLGQEBERkTgWEiIiIhLHQkJERETiWEiIiIhIHAsJERERiWMhISKKIaVUpnQGIiNiISEiijKl1FKllFZKaQC1SqmDSqlJ0rmIjMQmHYCIKAFkApgb+rxOa31IMAuRIbGQEBHFxiGtdZ10CCKj4pQNERERieMICRHRACml7gDweOjwCQDDEZyeGaG1fqzL0+9SStWEPp/fzd8TJTSltZbOQERkWqFS8iKAyZ1rQ5RSywFkaq2/HDqeAwBa602h40cAzO38eyJiISEiGhSl1FIAT2utJ4c9lgmgFmElpcvXTAJwEMAwrishCuIaEiKiwasLPwiVjDoAc7p7clhJ4aW/RCEsJEREUaSUylRK1Ybfd4Q3RyP6JBYSIqLByww/CBWOTACbQg9t7DJ1Mwk4u6aEiFhIiIgiYVKXUY/HAazQWnfee+TtLs9/HACvsiEKw8t+iajPnG6vDcAwAA4AVgTfQ2wAbDNhxa+RagHQ0c1HO4C6LE+2XyR49B0CsFQpVYfgupHq8Mt6tdZPKqW+FTqcDOBtrfWK2MckMi5eZUOU4JxubzqA0aGPMWGfdz0eg2AZUd29znxYd/wPUmf0cqoAgCoAJwAc7+XP41me7MZB/4fFSOgqm+Va67nnfTIR9YgjJEQJwun2jgAwM/QxI/TnZQDSYxTBgmCpGQNgVm9PLHeXNiBYUCoA7EJwLcZmANuzPNmtUc45EJnSAYjMjiMkRHHG6famAJiOs+Wjs4CMjeZ5+zBCEgkdCBaUzThbUrZkebLro3zeboVGRx4DsBTAk7z7KtHAsZAQmZzT7R0N4DoEfykuAjAFAgvWY1RIuqMBHMC5JWVjlie7ptevIiJDYSEhMhmn25sKYDGCBWQpgiMg3a7riCXBQtKdAICNAEoAeAFszvJk882OyMBYSIgMzun2WgEswNkCchWCV7kYisEKSVfHAbyGYEFZmeXJbhDOQ0RdsJAQGVDoypfbAdwKYAlit/B0wAxeSMK1AViD4MhJSZYn+6BwHiICCwmRYTjdXjuAmwHcB+DTAJJlE/WPiQpJV3sRKicA1sbxvVKIDI2FhEiY0+29BsEScheAEcJxBszEhSRcBYA/AHg2y5O9TzoMUSJhISES4HR7L0GwhHwewTt3ml6cFJJw6wA8C+AFrjkhij4WEqIYcbq9GQDuB/AFAPOF40RcHBaSTo0AigH8KsuT/ZF0GKJ4xUJCFGVOtzcLwL8BeBhAmnCcqInjQhJuPYBfIjhq0iYdhiiesJAQRYnT7Z0B4JsA7gFgF44TdQlSSDpVAvgNgKIsT3aZdBiieMBCQhRhTrd3CYBvIXjFTMJIsELSyQ/gTwC+l+XJPiwdhsjMWEiIIsDp9loA3IbgiMgC4TgiErSQdGoH8ByA/8ryZJdLhyEyIxYSokFwur3JAB4A8O8I7iGTsBK8kHRqBVAE4IksT/ZJ6TBEZsJCQjQATrdXIXjFzA8BjBOOYwgsJOdoRHDx65Pc5I+ob2K+IyiR2Tnd3oUANgB4Hiwj1L1UAI8BOFzuLv1eubvU8Lf+J5LGERKiPnK6vU4ATwK4UziKIXGEpFc1AH4E4OdZnuwm6TBERsRCQnQeTrc3DcC3AXwDJttfJpZYSPrkJIDvA3g6y5MdkA5DZCQsJEQ9CF058wCA/wZwgXAcw2Mh6Zf3ACzL8mTvlg5CZBRcQ0LUjdC9RD5C8OZXLCMUadcA2FLuLv1uubs07m+aR9QXHCEhCuN0e0cAeApAnnQWs+EIyYDtQHC0ZL10ECJJHCEhCnG6vS4EfzmwjFAszQDwXrm79Kfl7tJU6TBEUlhIKOE53d6hTrd3BYAScHqGZFgA/CuAHeXu0hukwxBJYCGhhBa6p8hWBHfiJZLmBPBmubv09+Xu0hHSYYhiiYWEEpLT7XU43d7lANYCmCSdh6iLLwDYVe4uvVs6CFGssJBQwnG6vZcD+BDBHXn5M0BGNRrAX8rdpX8pd5cOkQ5DFG18M6aE4XR7LU639zEEy8jl0nmI+uhuAB+Uu0sTevNGin8sJJQQnG5vFoA1ADwAHMJxiPprJoCN5e7SXOkgRNHCQkJxz+n2ZgPYCGCRdBaiQcgA8Eq5u/QH5e5SvndT3OE3NcU1p9v7LwDeATBGOgtRBCgA/wfAq+Xu0mHSYYgiiYWE4pLT7U1yur2/QfCuq7w1N8WbWxCcwuFaKIobLCQUd5xu7wUAVgN4SDgKUTRNAvB+ubv0XukgRJHAQkJxJXRJ73oAV0lnIYqBIQD+WO4u/Vm5u9QmHYZoMFhIKG443d6bAbwLYIJ0FqIY+zqAVeXu0tHSQYgGioWE4oLT7f0qgFcBpElnIRKSDWBdubvUKR2EaCBYSMjUnG6vcrq9PwPwCwBW6TxEwqYgWEoukw5C1F8sJGRaTrfXAuC3CA5XE1HQOABry92lXEdFpsJCQqYUKiPPAXhAOguRAQ0HsLLcXXqDdBCivmIhIdNxur1WAL8HcL90FiIDS9Vav/K3h9y83TyZAgsJmYrT7bUB+BMA3nuB6Dz2nt6w3tew4x+Febl50lmIzoeFhEzD6fbaAfwFAN9cic7j4Jkta7bWrl4MwAbgT4V5uZ+XzkTUGxYSMoVQGXkBwB3SWYiM7mjD7tUbq9+8NuwhK4DfF+blfkEqE9H5sJCQ4TndXgeAvwG4VTgKkeEdbzq0+v2qV5Z081dWAM8X5uV+KbaJiPqGhYQMzen2JgH4B4BPS2chMrpTLcfWrj354pJenmIB8NvCvFxenUaGw0JChhWapnkZwZ1NiagXda2V775z/I/ZfXiqBcCKwrzcG6Odiag/WEjIyFYAuEk6BJHR1bfXvP9WxfNXA1B9/BIbgBcK83JnRDEWUb+wkJAhOd3ebwP4knQOIqNr6jiz4fXy387T0P3dOiEdgLcwL/eCaOQi6i8WEjIcp9t7F4D/ls5BZHQt/sZN3vIVl2sE7AN8iQkAXi3Myx0SyVxEA8FCQobidHuvAvA79H3omSghtflbtpWUPX1pQPuTB/lS8xC8Twl/H5AofgOSYTjd3osA/C+Awb7BEsW19kDb7pLyool+3Z4aoZe8FcCTEXotogFhISFDcLq9GQC8AEZLZyEyso5A+/6SsqIx7YHWjAi/dEFhXu6jEX5Noj5jISFxof1pXgIwTToLkZH5td/nLV+R0RZoHh6lU/yClwOTFBYSMoJfA1gqHYLIyALaX/5a+TOOFn9DNEcReTkwiWEhIVFOt/ebAJZJ5yAyMq0DJ9849qy/qeP0uBicjpcDkwgWEhLjdHtvBLBcOgeRkWmtq9+q+F1DfXvNxBietvNy4JQYnpMSHAsJiXC6vaPBy3uJeqW1Pv3O8T9V1bVVThY4/TzwHwwUQywkJOVZAGOkQxAZlda6Yc2J4rLq1mNTBWN8tTAv9wbB81MCYSGhmHO6vV8F4JLOQWRUWuvmdyv/fuBkyxHpxaUKwHOFebnRuqqH6GMsJBRTTrd3BoAfSecgMiqtddv6qpKdFU0HZktnCRkH4GnpEBT/WEgoZpxubzKAP4N3YiXqltba/1H1W5uONO6aJ52lizsK83Lvlw5B8Y2FhGJpOYCZ0iGIjEhrHdhWu+aDg/VbrpLO0oNfFOblOqVDUPxiIaGYcLq9NwP4unQOIqPaffqDdXtOr18onaMX6QB+z034KFr4jUVRF7rE93npHERGdeD
      "text/plain": [
       "<Figure size 640x395.55 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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}