docs/notebooks/Analysis_results_showcase.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis Results\n",
    "\n",
    "3ML stores the results of a fit in a container we call an \"Analysis Result\" (AR). The structure of this object is designed to be useable in a *live* sense within an *active* analysis (python script, ipython interactive shell, jupyter notebook) as well as storable as a FITS file for saving results for later.\n",
    "\n",
    "The structure is nearly the same between MLE and Bayesian analyses in order to make a seamless functionality between all analyses.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:59.058491Z",
     "start_time": "2020-04-13T07:50:50.363Z"
    }
   },
   "outputs": [],
   "source": [
    "from threeML import *\n",
    "\n",
    "from threeML.analysis_results import *\n",
    "\n",
    "from threeML.io.progress_bar import progress_bar\n",
    "\n",
    "from jupyterthemes import jtplot\n",
    "\n",
    "%matplotlib inline\n",
    "jtplot.style(context=\"talk\", fscale=1, ticks=True, grid=False)\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.style.use(\"mike\")\n",
    "\n",
    "\n",
    "import astropy.units as u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at what we can do with an AR. First, we will simulate some data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:40.643650Z",
     "start_time": "2020-04-13T07:50:40.387692Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using Gaussian statistic (equivalent to chi^2) with the provided errors.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoEAAAIqCAYAAABbiJCVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dfZBV9X348Q/sWlgF42BbH0qINqJoYoWJM1oftpQl6DqKUKdWgTa4bqdJBWIsNWoyKtpSDCaY1JlExGVrWrfAGFpbNZHCVKC2xYB0HepDoIhMTU1NUQhUZB9+f9hcXZCfe1fuPefc7+s1k5l7bnaXz3ruw3u/59x7B/X29vYGAABJGZz1AAAAVJ8IBABIkAgEAEiQCAQASJAIBABIkAgEAEhQzUZgY2Nj1iMAAORWzUYgAACHJwIBABIkAgEAElT1CJw1a1Zs2bKltP3UU0/F1VdfHc3NzfGVr3wl9uzZU+2RAACSU7UI3LBhQyxcuDA6OztL1/3oRz+KP//zP4+vfOUrsXz58jhw4EB0dHRUayQAgGTVV+sf6uzsjJ6enmhoaChd98QTT8SVV14Z55xzTkRE3HrrrbF79+5qjQQAkKyqRWBra2tERGzcuLF03UsvvRRjxoyJz33uc/GTn/wkLrzwwrjpppuqNRIAQLIyfWHIW2+9Ff/8z/8c8+bNi0ceeST++7//Ox566KEsRwIASELVVgIPZ8qUKXHKKadERMQ111wTS5YsKev729raor29/cgPBgBQwzKNwBNPPDG6urpK2z09PfELv/ALZf2MlpaWaGlpOeR6nxgCAHB4mR4Obmpqir/5m7+J//iP/4hdu3bFX//1X0dTU1OWIwEAJCHTlcDm5ubYtWtXzJ07N3p6euKSSy6JqVOnZjkSAEASBvX29vZmPUQlNDY2xtq1a7MeAwAgl3xsHABAgkQgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIECVdE1oiq4JPh8dyAcRCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEI8BF0TWiKrglNWY8BUDYRCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkKCqR+CsWbNiy5Yth1z/4IMPxty5c6s9DgBAkuqr9Q9t2LAhnn766ejs7Dzk/3v55Zdj2bJlMXbs2GqNAwCQtKqtBHZ2dkZPT080NDT0ub6rqysWLlwYV1xxRbVGAQBIXtVWAltbWyMiYuPGjX2u/+53vxuf+cxn4hOf+ETs3LmzWuMAACStahH4QbZt2xb/+I//GIsXL441a9YM6Ge0tbVFe3v7kR0MAKDGZRaBPT09cc8998SXvvSlGDJkyIB/TktLS7S0tBxyfWNj40cZDwCgpmUWgXv37o0XX3wx5syZ0+f6yy67LJ544omMpgIASENmETh8+PBYu3ZtafvJJ5+M1atXx7333pvVSAAAyfBm0QAACar6SuDy5cs/8Prm5uZobm6u8jQAAGmyEggAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQgAkCARCACQIBEIAJAgEQhQBb1dXaXL3Q8u6bMNkAURCFAFPUvbS5d7Ozqip739sF8LUA0iEKAaXnjhoO0Xs5kD4P+IQIBqOPPMg7bHZDMHwP8RgQBVMPi6maXLg6ZdG4Nnzjzs1wJUQ33WAwCkYFD9ew+3da2tGU4C8C4rgQAACRKBAAAJEoEAAAkSgQAACRKBAAAJEoEAAAkSgQAACRKBAAAJEoEAAAkSgQAACRKBAAAJEoEAAAkSgQAACRKBAAAJEoEAAAkSgQAACRKBAAAJEoEAAAkSgQAD1NvVVbrc/eCSPtsAeScCAQaoZ2l76XJvR0f0tLcf9msB8kYEAgzUCy8ctP1iNnMADIAIBBioM888aHtMNnMADIAIBBigwdfNLF0eNO3aGDxz5mG/FiBvqh6Bs2bNii1btkRExNtvvx3z58+PSy+9NKZMmRLf/OY3o8uJ1UBBDKqvL12ua23tsw2Qd1WLwA0bNsTChQujs7OzdN1f/dVfxSuvvBJLly6NxYsXx5YtW2LlypXVGgkAIFlVi8DOzs7o6emJhoaG0nU//OEPY9q0aXHSSSfFL//yL8fkyZNj48aN1RoJACBZVTt20draGhHRJ/JuvvnmOOGEE0rbL730UvzSL/1StUYCAEhWpiewnHrqqRHx7rmBixcvjnXr1sW3v/3tsn5GW1tbtHtvLgCAsmR+FvOmTZviz/7sz+L000+PJUuWxC/+4i+W9f0tLS3R0tJyyPWNjY1HakQAgJqTaQT+8Ic/jDvuuCNuueWWuPjii7McBQAgKZm+T2BbW1vMnj1bAAIAVFmmK4Hbtm2L+fPnx/z580vXjR07Nr71rW9lOBUAQO2regQuX768dPkHP/hBtf95AADCx8YBACRJBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJKg+6wEAUlG/ZnXWIwCUWAkEAEiQCAQASJAIBABIkAgEAEiQCAQASJAIBABIkAgEAEiQCAQASJAIBABIkAgEAEiQCAQASFDVI3DWrFmxZcuWiIjo7e2N+++/P5qbm2Pq1KmxcuXKao8DAJCk+mr9Qxs2bIinn346Ojs7S9c9/vjj8eyzz8bSpUtjz549cdNNN8VZZ50VZ5xxRrXGAgBIUtVWAjs7O6OnpycaGhpK161atSqmT58eJ554YowePTqamppizZo11RoJACBZVVsJbG1tjYiIjRs3lq7btm1bjB49urR96qmnxrPPPlutkQAAklW1CPwge/fujWOOOaa03dDQEPv27SvrZ7S1tUV7e/sRngwAoLZlGoHDhg2L/fv3l7b3798fw4cPL+tntLS0REtLyyHXNzY2fuT5AABqVaZvETNq1KjYtm1bafuVV16J0047LcOJAADSkGkETpw4MTo6OmLXrl3x/PPPx6pVq2L8+PFZjgQAkIRMDwdPnjw5tm/fHtdcc00MGzYsZs+eHR//+MezHAkAIAmDent7e7MeohIaGxtj7dq1WY8B1LiuCU0REVG/ZnXGkwCUx8fGAQAkSAQCACQo03MCAYrOYWCgqKwEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQCACRIBAIAJEgEAgAkSAQC5FDXhKbomtCU9RhADROBAAAJEoEAAAkSgQAACRKBAAAJEoEAAAkSgQAACepXBC5cuDA2bdoUvb29FRmio6Mjrrrqqpg0aVLMmTMnduzYUZF/BwCAd/UrAl977bWYO3duTJkyJb7+9a/H5s2bj1gQbtmyJZYtWxb33XdfPPbYYzFmzJj4xje+cUR+NgAAH6y+P1+0aNGi+NnPfhbPPPNMrF+/Pm655ZZoaGiICRMmxKRJk+KMM84Y8AB1dXUREdHd3R0RET09PXHssccO+OcBAPDh+hWBERHDhg2LSZMmxaRJk+Kdd96JFStWxMMPPxyPPvpofOITn4gpU6bE5ZdfHkcddVRZA4wZMyYuuuii+L3f+72IiBgyZEi0tbWV91sAAFCWfkdgb29v/Nu//Vv80z/9U6xfvz5ef/31+MxnPhO/+Zu/GYMGDYrly5fHv/7rv8aCBQvKGmDjxo2xbt26WLx4cYwcOTIeeOCBuOuuu+LBBx8s+5cBAKB/+hWBf/qnfxr/8i//Ej/72c9i3LhxMX369GhsbOxz2PZXf/VX44Ybbih7gHXr1sWkSZNizJgxERHx+c9/Ppqbm2P37t39Oizc1tYW7e3tZf+7AAAp61cEvvHGG/H7v//78Ru/8RvxsY997AO/5ld+5Vdi0aJFZQ8wZMiQ6OnpKW3X1dVFXV1dDBkypF/f39LSEi0tLYdc39jYWPYsAACp6NergxctWhSTJ08+bABGvHvO4Nlnn132ABdccEE89dRT8e///u+xb9++WLJkSZx33nn9jkAAAMrX73MCK+Wcc86Jz3/+83H33XfHrl27YuzYsXHzzTdnPRYAQE3LPAIjIpqbm6O5uTnrMQAAkuFj4wAAEiQCAQASJAIBABIkAgEAEiQCAQASJAIBABIkAgEAEiQCAQASJAIBABIkAgEAEiQCAQASJAIBABIkAgEAEiQCAQASJAIBABIkAgFyprerq3S5+8ElfbYBjhQRCJAzPUvbS5d7Ozqip739sF8LMFAiECBvXnjhoO0Xs5kDqGkiECBvzjzzoO0x2cwB1DQRCJAzg6+bWbo8aNq1MXjmzMN+LcBA1Wc9AAB9Dap/76G5rrU1w0mAWmYlEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhDgA3RNaIquCU1ZjwFQMSIQACBBIhAAIEG5iMDXXnst5syZE5MmTYrf/d3fjeeeey7rkQAAalouIvCOO+6Ic845J773ve/FtGnT4u677856JIBM1a9ZHfVrVmc9BlDD6rMe4KWXXoo333wzrrvuuhg8eHA0NzfHyJEjo7u7O+rq6rIeDwCgJmW+Evjyyy/HySefHLfffns0NzfH9ddfH4MHDxaAAAAVlHkEvvXWW7F58+a44IIL4tFHH43m5ub46le/Gvv37896NACAmpX54eCIiFNPPTUuu+yyiIi46qqr4uGHH46dO3fGaaed9qHf29bWFu3t7RWeEACgtmQegSeeeGJ0dXX1ua63tzeGDBnSr+9vaWmJlpaWQ65vbGw8IvMBANSizA8Hn3feefE///M/8bd/+7exb9++WLFiRZxwwgkxcuTIrEcDAKhZmUfg8OHD4957743HH388fuu3fiueeeaZuOuuu2LQoEFZjwYAULMyPxwcEXHWWWfF4sWLsx4DACAZma8EAgBQfSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhDgIL1dXaXL3Q8u6bMNUCtEIMBBepa2ly73dnRET3v7Yb8WoKhEIMDBXnjhoO0Xs5kDoIJEIMDBzjzzoO0x2cwBUEEiEOAgg6+bWbo8aNq1MXjmzMN+LUBR1Wc9AEDeDKp/76GxrrU1w0kAKsdKIABAgkQgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIABAgnIVgdu3b4+mpqb46U9/mvUoAAA1LTcR2N3dHV/72tfiwIEDWY8CAFDzchOBy5Yti0996lNZjwFQOF0TmqJrQlPWYwAFk4sIfPXVV+P73/9+tLa2Zj0KAEAS6rMeoKenJxYsWBBz5syJoUOHlv39bW1t0d7efuQHAwCoYZlH4KOPPhqjRo2Kc889d0Df39LSEi0tLYdc39jY+FFHAxJWv2Z11iMAVFTmh4Ofe+65eOKJJ6KxsbEUblOnTo1169ZlPBkAQO3KfCVw/vz5fbYbGxtj5cqVcfzxx2c0EQBA7ct8JRAAgOrLfCXwYGvXrs16BACAmmclEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhCgwHq7ukqXux9c0mcb4P9HBAIUWM/S9tLl3o6O6PFZ6kA/iUCAInvhhYO2X8xmDqBwRCBAkZ155kHbY7KZAygcEQhQYIOvm1m6PGjatTF45szDfi3A++XuY+MA6L9B9e89jNe1tmY4CVA0VgIBABIkAgEAPqKuCU3RNaEp6zHKIgIptCLe6QAgD0QgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIEDB1a9ZHfVrVmc9BiSrt6urdLn7wSV9tvNMBFJYRb3TAXDkZfkJUj1L20uXezs6oqe9/bBfG5GfT7sSgRRWuXc6AKiIF144aPvFbOYokwikuAp6pwOgxpx55kHbY7KZo0wikOIq6J0OqD15ObxHNgZfN7N0edC0a2PwzJmH/do8EYEUVlHvdGTHEzVuA1TCoPr60uW61tY+23lWjCnhAxx8pwMA+s9KIABAgkQgAECCRCBwCOdNkQduh1BZIpBkeEIBgPeIQACg0Ir0CVJ5mlUEAvCRWGXvvyL9tyrSrEX6BKk8zSoC4SMo0oMkQM0q0idI5WhWEQgAH0GeDu8lKwefIFW/ZnXUr1n94V+Yg1l/TgQmolZXrPp9pwOSV6lYy9Phvaxk/RxTpE+QytOsPjGE3Pn5A4m4g+xU6n6Y5f37kFgbdIQ+bShHh/dSVaRPkMrTrFYCAUhDpWItR4f3oBwiEPhIsj4M1F/O26LcWOvvbTtPh/egHCIQ6KNWY8l5W8VSidthpWLt4MN779+GPBOBJKFWw6YSajaWnLdVKJW4HYo16EsEkoSihU0lDrH2+2fWaiw5b6tYavV2CDkiAkmDJ5T+q9FYct5WZVRslb0Gb4dFOiJRpFkrqSjnPA+UCCQNFXhCycuD5JF+kCpaLPX393cosDIqtcpetNthfxTpiEReZq31CMuaCKwSN+RsVeIJJS8PkkeaWKIsFVplr8nbYZGOSFRwVs+H+SECE5CXFassVeQJpUgP6AXjSaJAavCwbcUU6b9VGbPm5TmmSJ8glZdZReBHUJQnqiKtWOXlwaRfivSADmUo535Yzip7oe7fFVCkQ9zlzFqk5xj6EoEpKNCKVZEeTIr0gA7lKOd+WM4qe5Hu35WQl0Pc/VnAKGvWAj3H0FcuIvDJJ5+Mq6++Oi655JKYPXt2vPrqq1mPVFuKtGJVoAeTvDygZyn1lZ2aVan7YYHu3+XetvNyeC8TRXqOoY/MI3Dnzp2xaNGi+OIXvxgrV66Mc889N+bNm5f1WDWlUCtWHkwKJfWVnZpVqfthDu7f/Y01t+3+K9RzDH1kHoGbNm2KcePGxYUXXhhHH310TJ8+PbZu3Rq7d+/OerSaUaQVKw8mBVOglR36r1L3w0r93Iqcn52D23ZRzjsv0nMMfWUegRdffHHceOONpe2XX345hg4dGsccc0yGU9FfR/pByoNJZQ6x5uUNfYvypJa6St0PC3X/zsGqJVRa5vfAESNGlC6vW7cuFi5cGF/4wheirq6uX9/f1tYW7Zbpj6ifP0kne35Lxg45DDXo3SfMvP3MiHdXdro7OiLCyi21xW2bg/94HnzdzHz/4TIAufhtdu/eHQsWLIitW7fGbbfdFueff36/v7elpSVaWloOub6xsfFIjkhCMo/gMg5D9ftBqsxDW/393Q9e2am2FB6kU5blH6JZ37bJ/v5dqT+e8yTzw8Fvv/123HDDDTFixIh4+OGHywrAovAKynzI+tV7/T4UWsZhqH6fvF6jh7acvA/FUs7zYeb37xycF1ppmUfgqlWr4vjjj4+5c+fG0KFDsx6nIjK/IVMoZZ08388HqSK94KasP5oGsMLpNAfITlnPh1lHWI3+8fx+mR832bp1a2zatOmQw7fLli2Lk046KaOpPlxZy9Rl3pAzPxxJpso6DHXmmRGbN79v+4MfpIp0aKusQzD9/P2pLI9V/Vek/1YVmbWc58OM798pnBea+Urgl770pVi7du0h/8tzAEaU+ddMAn9NpCrrlaUirfD1WxlPEjX5+0MFZX56UhnPh1nfvwv1avYByjwCC8sTFTlQkw9SZTxJ1OTvDwPQ33OOsz49qZznQ/fvyvNfdKDKWKbOw6G4Ih2CIG0pHIKBzFToPLuivKNAhFOu3s9K4ABZ3aucrA+xki1//VdWf1eMKnU/PNI/N/PDmxWSlzd4p7aJwAHyRFW7D75UhtsLlZD14c1KqdTvZQGD90uvXPj/KudVzym8kSYfrr+rOm4vVEQF30Yk0yMSFfq98nA4lvywEkgfhXoPpwqwWlVBNXh7IQdq9fBmrf5e5IoIpK9y38Opz3bxH6QqdQhGXEZN3l7IXq0e3qzV34t8EYEF1u+PIStHgd7DqSIqtFpVblxW4qT8rF9wU5O3FzJXq+dn1+rvRb6IwJzJesUo+fdwqtRqlUOhtXl7AWpa1n88V5pH4ZzJ+uT51E8arth71FXw449q+QGqP1L//aGWuX9XlpXAKun3XxNWjDJVqdUqh0KBSsn6CFK5an11rUhEYN44eb4mORRaHk8S0H+1+l6JlVC0YK40EfgRVOKJqlZXjCryIhYAHEEqg2DuSwTmjBUjoFKsgtSoMo8gJb3SLpj7EIEAibAKUptq9QhSRTjlqg8RCJAKqyA1yRGk/hPMfbmlAKSigm9VlKVaPbRZq79XllJ/G7SDWQksqLyc25P0uSWUze0lW1ZBgPezElhQlXxTaU/SUJusggDvJwKLyrk9FSOCAUiBCCyqGj23p5aJSwDyxDmBBeXcHgDgoxCBOdSfk+eL9JYAeXkRCwDwHhFIxXmDWgDIn/wuH1E7vIgFoKKcc8xAiEAqz4tYAMgJwfweh4OpOC9iAYD8sRJIxXmDWsgPqyDAz1kJBABIkAgEAEiQCAQASJBzAgvMuT0AwEBZCQQASJAIBABIkAgEAEiQCAQASJAXhlAVXsQCAPliJRAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEEiEAAgQSIQACBBIhAAIEG5iMANGzbEtddeG5MmTYrbb7899u7dm/VIAAA1LfMI3LNnT8ybNy/+4A/+IJYvXx69vb2xZMmSrMcCAKhpmUfg+vXr4/TTT4/x48fHcccdFzNmzIg1a9ZkPRYAQE3LPAK3bdsWo0ePLm2fcsopsWvXrti3b1+GUwEA1Lb6rAfYt29fnHDCCaXtIUOGRF1dXezbty+OPvroD/3+tra2aG9vr+CEAAC1J/MIHDZsWLzzzjul7QMHDkR3d3cMHz68X9/f0tISLS0th1zf2Nh4xGYEAKg1mR8OHjVqVGzbtq20vWPHjhg5cmQMGTIkw6kAAGpb5hF40UUXRWdnZ2zcuDHefPPNeOCBB2LixIlZjwUAUNMyPxx83HHHxa233hr33HNPvPXWWzF+/PiYMWNG1mMBANS0Qb29vb1ZD1EJzgkEAIhYu3btB15fsxFYLY2NjYf9j0t+2E/5Zx8Vg/2Uf/ZRMeRhP2V+TiAAANUnAgEAEiQCAQASJAIBABJUd+edd96Z9RBF1tvbG+PGjct6DD6E/ZR/9lEx2E/5Zx8VQx72k1cHAwAkyOFgAIAEiUAAgASJQACABIlAAIAEiUAAgASJQACABIlAAIAEicAB2rBhQ1x77bUxadKkuP3222Pv3r1Zj8T7zJo1K7Zs2RIR774h5/333x/Nzc0xderUWLlyZcbT8eSTT8bVV18dl1xyScyePTteffVV+ylnvve978WUKVPi0ksvjS9/+cvxxhtv2Ec5tn379mhqaoqf/vSnERHxyCOPxBVXXBGXX355PPTQQxlPl7Y//MM/jMbGxtL/7r777tzcl+oz+VcLbs+ePTFv3rz44z/+4xg7dmx8/etfjyVLlsQXv/jFrEdL3oYNG+Lpp5+Ozs7O0nWPP/54PPvss7F06dLYs2dP3HTTTXHWWWfFGWeckeGk6dq5c2csWrQo7rjjjhg3blysWLEi5s2bF1OnTrWfcmL79u3x0EMPxYIFC+KUU06Je++9N77zne/E2LFj7aMc6u7ujq997Wtx4MCBiIjYuHFjrFixIu67775oaGiIP/qjP4rRo0dHY2NjxpOm6fXXX49Vq1bFkCFDStf9/d//fS7uS1YCB2D9+vVx+umnx/jx4+O4446LGTNmxJo1a7Iei4jo7OyMnp6eaGhoKF23atWqmD59epx44okxevToaGpqsr8ytGnTphg3blxceOGFcfTRR8f06dNj69at9lOObN68OS644II4++yzY/jw4XHJJZfE9u3b7aOcWrZsWXzqU58qbT/11FMxZcqU+OQnPxknn3xyXHnllbF69eoMJ0zX/v37o66urk8ARuTneUkEDsC2bdti9OjRpe1TTjkldu3aFfv27ctwKiIiWltb48tf/nJ87GMfK1138P469dRT4z//8z+zGI+IuPjii+PGG28sbb/88ssxdOjQ+NGPfmQ/5cQVV1wRN998c/T29sabb74ZTz31VHz60592X8qhV199Nb7//e9Ha2tr6bpt27bF6aefXtq2n7LzX//1X9Hd3R3XX399XHrppXHLLbfEG2+8kZv7kggcgH379sUxxxxT2h4yZEjU1dWJwJzau3dvn/3V0NBgX2VoxIgRcT67zccAAARZSURBVNJJJ0VExLp16+KWW26JL3zhC/G///u/9lNO1NfXx1FHHRWPPfZYTJ48OdavXx8TJkxwX8qZnp6eWLBgQcyZMyeGDh1auv7g/XT00UfbTxnZvXt3jBw5Mm677bZYvnx5DB8+PP7kT/4kN/cl5wQOwLBhw+Kdd94pbR84cCC6u7tj+PDhGU7F4QwbNiz2799f2t6/f799lbHdu3fHggULYuvWrXHbbbfF+eefHw899JD9lDNXXnllfPazn43HH388vvrVr0ZDQ4N9lCOPPvpojBo1Ks4999w+1w8fPrzPfnr77bftp4ycffbZ8c1vfrO0fcMNN8SVV14Zxx57bC7uS1YCB2DUqFGxbdu20vaOHTti5MiRhxzzJx8O3l+vvPJKnHbaaRlOlLa33347brjhhhgxYkQ8/PDDcf7550eE/ZQnf/mXfxmrVq2KiHdXkaZMmRK7d++OESNG2Ec58txzz8UTTzxRetVpRMTUqVNjx44d9lNOPPPMM7Fhw4bSdldXVxx11FExcuTIXOwjETgAF110UXR2dsbGjRvjzTffjAceeCAmTpyY9VgcxsSJE6OjoyN27doVzz//fKxatSrGjx+f9VjJWrVqVRx//PExd+7cPoew7Kf8OPbYY+O73/1u7Ny5M/bt2xePPPJIfPKTn4yrrrrKPsqR+fPnx9q1a0v/i4hYuXJl3HnnnbFy5cr48Y9/HK+88kqsWLHCc1RG9u7dG9/4xjdi+/btsXv37vjOd74T48ePj0mTJuXivuRw8AAcd9xxceutt8Y999wTb731VowfPz5mzJiR9VgcxuTJk2P79u1xzTXXxLBhw2L27Nnx8Y9/POuxkrV169bYtGnTIW9X0dHRYT/lxGWXXRY7duyI2bNnx/79++PTn/503HXXXXHSSSfZRwXw67/+6/Hiiy/G9ddfH/X19TF9+vQYN25c1mMl6bOf/Wzs3Lkzbrzxxujp6YnzzjsvbrzxxmhoaMjFfWlQb29vb9X/VQAAMuVwMABAgkQgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIABAgkQgAECCRCAAQIJEIMAR8hd/8RcxceLEeO211yIiYufOnTFx4sR48sknM54M4FAiEOAImTZtWpx88snxrW99KyIi7rvvvjj77LOjubk548kADuWzgwGOoM7Ozpg9e3Zcfvnl8YMf/CDa29tj5MiRWY8FcAgrgQBH0K/92q/F5ZdfHn/3d38Xn/vc5wQgkFsiEOAIe/311yMi4ic/+UnGkwAcnggEOIL+4R/+ITZv3hyzZ8+Oxx57LJ5//vmsRwL4QCIQ4AjZs2dP3H///fE7v/M78du//dtx4YUXxr333htdXV1ZjwZwCBEIcIR8+9vfjrq6upgxY0ZERMyZMyd+/OMfxyOPPJLxZACH8upgAIAEWQkEAEiQCAQASJAIBABIkAgEAEiQCAQASJAIBABIkAgEAEiQCAQASJAIBABI0P8DZyEHiojjEdkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 748.8x655.2 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gen_function = Line(a=0.0, b=2) + Gaussian(F=30.0, mu=25.0, sigma=1)\n",
    "\n",
    "# Generate a dataset using the line and a gaussian.\n",
    "# constant 20% error\n",
    "\n",
    "x = np.linspace(0, 50, 50)\n",
    "\n",
    "xy = XYLike.from_function(\n",
    "    \"sim_data\", function=gen_function, x=x, yerr=0.2 * gen_function(x)\n",
    ")\n",
    "\n",
    "xy.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MLE Results\n",
    "\n",
    "\n",
    "First we will demonstrate how AR's work for an MLE analysis on our synthetic data. As we will see, most of the functionality exists in the Bayesian AR's as well. \n",
    "\n",
    "Let's do a simple likelihood maximization of our data and model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:40.816499Z",
     "start_time": "2020-04-13T07:50:40.646361Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best fit values:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>result</th>\n",
       "      <th>unit</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>parameter</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>(9 +/- 4) x 10^-3</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>1.76 +/- 0.11</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>(3.2 +/- 0.4) x 10</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>(2.521 +/- 0.012) x 10</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>1.09 +/- 0.09</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                      result              unit\n",
       "parameter                                                                     \n",
       "fake.spectrum.main.composite.a_1           (9 +/- 4) x 10^-3  1 / (cm2 keV2 s)\n",
       "fake.spectrum.main.composite.b_1               1.76 +/- 0.11   1 / (cm2 keV s)\n",
       "fake.spectrum.main.composite.F_2          (3.2 +/- 0.4) x 10       1 / (cm2 s)\n",
       "fake.spectrum.main.composite.mu_2     (2.521 +/- 0.012) x 10               keV\n",
       "fake.spectrum.main.composite.sigma_2           1.09 +/- 0.09               keV"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Correlation matrix:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table id=\"table4718698384\">\n",
       "<tr><td>1.00</td><td>-0.85</td><td>0.00</td><td>-0.02</td><td>-0.00</td></tr>\n",
       "<tr><td>-0.85</td><td>1.00</td><td>-0.05</td><td>0.02</td><td>-0.08</td></tr>\n",
       "<tr><td>0.00</td><td>-0.05</td><td>1.00</td><td>-0.24</td><td>-0.15</td></tr>\n",
       "<tr><td>-0.02</td><td>0.02</td><td>-0.24</td><td>1.00</td><td>-0.13</td></tr>\n",
       "<tr><td>-0.00</td><td>-0.08</td><td>-0.15</td><td>-0.13</td><td>1.00</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       " 1.00 -0.85  0.00 -0.02 -0.00\n",
       "-0.85  1.00 -0.05  0.02 -0.08\n",
       " 0.00 -0.05  1.00 -0.24 -0.15\n",
       "-0.02  0.02 -0.24  1.00 -0.13\n",
       "-0.00 -0.08 -0.15 -0.13  1.00"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of -log(likelihood) at the minimum:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>-log(likelihood)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sim_data</th>\n",
       "      <td>24.81213</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>24.81213</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          -log(likelihood)\n",
       "sim_data          24.81213\n",
       "total             24.81213"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of statistical measures:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>statistical measures</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AIC</th>\n",
       "      <td>60.987896</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BIC</th>\n",
       "      <td>69.184375</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     statistical measures\n",
       "AIC             60.987896\n",
       "BIC             69.184375"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fitfun = Line() + Gaussian()\n",
    "\n",
    "fitfun.a_1.bounds = (-10, 10.0)\n",
    "fitfun.b_1.bounds = (-100, 100.0)\n",
    "fitfun.F_2 = 25.0\n",
    "fitfun.F_2.bounds = (1e-3, 200.0)\n",
    "fitfun.mu_2 = 25.0\n",
    "fitfun.mu_2.bounds = (0.0, 100.0)\n",
    "fitfun.sigma_2.bounds = (1e-3, 10.0)\n",
    "\n",
    "model = Model(PointSource(\"fake\", 0.0, 0.0, fitfun))\n",
    "\n",
    "data = DataList(xy)\n",
    "\n",
    "jl = JointLikelihood(model, DataList(xy))\n",
    "_ = jl.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get our errors as always, but the results cannot be propagated (error propagation assumes Gaussian errors, i.e., symmetric errors)\n",
    "In this case though errors are pretty symmetric, so we are likely in the case\n",
    "where the MLE is actually normally distributed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:40.893996Z",
     "start_time": "2020-04-13T07:50:40.818569Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>result</th>\n",
       "      <th>unit</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>parameter</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>(9 +/- 4) x 10^-3</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>1.76 +/- 0.11</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>(3.2 +/- 0.4) x 10</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>(2.521 -0.011 +0.012) x 10</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>1.09 +/- 0.09</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                          result  \\\n",
       "parameter                                                          \n",
       "fake.spectrum.main.composite.a_1               (9 +/- 4) x 10^-3   \n",
       "fake.spectrum.main.composite.b_1                   1.76 +/- 0.11   \n",
       "fake.spectrum.main.composite.F_2              (3.2 +/- 0.4) x 10   \n",
       "fake.spectrum.main.composite.mu_2     (2.521 -0.011 +0.012) x 10   \n",
       "fake.spectrum.main.composite.sigma_2               1.09 +/- 0.09   \n",
       "\n",
       "                                                  unit  \n",
       "parameter                                               \n",
       "fake.spectrum.main.composite.a_1      1 / (cm2 keV2 s)  \n",
       "fake.spectrum.main.composite.b_1       1 / (cm2 keV s)  \n",
       "fake.spectrum.main.composite.F_2           1 / (cm2 s)  \n",
       "fake.spectrum.main.composite.mu_2                  keV  \n",
       "fake.spectrum.main.composite.sigma_2               keV  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "jl.get_errors();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to get the AnalysisResults object that is created after a fit is performed. The AR object is a member of the JointLikelihood object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:40.897883Z",
     "start_time": "2020-04-13T07:50:40.895729Z"
    }
   },
   "outputs": [],
   "source": [
    "ar = jl.results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can display the results of the analysis. Note, when a fit is performed, the post display is actaully from the internal AR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:40.942770Z",
     "start_time": "2020-04-13T07:50:40.901268Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best fit values:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>result</th>\n",
       "      <th>unit</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>parameter</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>(9 +/- 4) x 10^-3</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>1.76 +/- 0.11</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>(3.2 +/- 0.4) x 10</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>(2.521 +/- 0.012) x 10</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>1.09 +/- 0.09</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                      result              unit\n",
       "parameter                                                                     \n",
       "fake.spectrum.main.composite.a_1           (9 +/- 4) x 10^-3  1 / (cm2 keV2 s)\n",
       "fake.spectrum.main.composite.b_1               1.76 +/- 0.11   1 / (cm2 keV s)\n",
       "fake.spectrum.main.composite.F_2          (3.2 +/- 0.4) x 10       1 / (cm2 s)\n",
       "fake.spectrum.main.composite.mu_2     (2.521 +/- 0.012) x 10               keV\n",
       "fake.spectrum.main.composite.sigma_2           1.09 +/- 0.09               keV"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Correlation matrix:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table id=\"table4720008208\">\n",
       "<tr><td>1.00</td><td>-0.85</td><td>0.00</td><td>-0.02</td><td>-0.00</td></tr>\n",
       "<tr><td>-0.85</td><td>1.00</td><td>-0.05</td><td>0.02</td><td>-0.08</td></tr>\n",
       "<tr><td>0.00</td><td>-0.05</td><td>1.00</td><td>-0.24</td><td>-0.15</td></tr>\n",
       "<tr><td>-0.02</td><td>0.02</td><td>-0.24</td><td>1.00</td><td>-0.13</td></tr>\n",
       "<tr><td>-0.00</td><td>-0.08</td><td>-0.15</td><td>-0.13</td><td>1.00</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       " 1.00 -0.85  0.00 -0.02 -0.00\n",
       "-0.85  1.00 -0.05  0.02 -0.08\n",
       " 0.00 -0.05  1.00 -0.24 -0.15\n",
       "-0.02  0.02 -0.24  1.00 -0.13\n",
       "-0.00 -0.08 -0.15 -0.13  1.00"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of -log(likelihood) at the minimum:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>-log(likelihood)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sim_data</th>\n",
       "      <td>24.81213</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>24.81213</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          -log(likelihood)\n",
       "sim_data          24.81213\n",
       "total             24.81213"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of statistical measures:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>statistical measures</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AIC</th>\n",
       "      <td>60.987896</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BIC</th>\n",
       "      <td>69.184375</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     statistical measures\n",
       "AIC             60.987896\n",
       "BIC             69.184375"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ar.display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, the equal tail intervals are displayed. We can instead display highest posterior densities (equal in the MLE case)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:40.986221Z",
     "start_time": "2020-04-13T07:50:40.945384Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best fit values:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>result</th>\n",
       "      <th>unit</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>parameter</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>(9 +/- 4) x 10^-3</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>1.76 +/- 0.11</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>(3.2 +/- 0.4) x 10</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>(2.521 +/- 0.012) x 10</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>1.09 +/- 0.09</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                      result              unit\n",
       "parameter                                                                     \n",
       "fake.spectrum.main.composite.a_1           (9 +/- 4) x 10^-3  1 / (cm2 keV2 s)\n",
       "fake.spectrum.main.composite.b_1               1.76 +/- 0.11   1 / (cm2 keV s)\n",
       "fake.spectrum.main.composite.F_2          (3.2 +/- 0.4) x 10       1 / (cm2 s)\n",
       "fake.spectrum.main.composite.mu_2     (2.521 +/- 0.012) x 10               keV\n",
       "fake.spectrum.main.composite.sigma_2           1.09 +/- 0.09               keV"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Correlation matrix:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table id=\"table4719530576\">\n",
       "<tr><td>1.00</td><td>-0.85</td><td>0.00</td><td>-0.02</td><td>-0.00</td></tr>\n",
       "<tr><td>-0.85</td><td>1.00</td><td>-0.05</td><td>0.02</td><td>-0.08</td></tr>\n",
       "<tr><td>0.00</td><td>-0.05</td><td>1.00</td><td>-0.24</td><td>-0.15</td></tr>\n",
       "<tr><td>-0.02</td><td>0.02</td><td>-0.24</td><td>1.00</td><td>-0.13</td></tr>\n",
       "<tr><td>-0.00</td><td>-0.08</td><td>-0.15</td><td>-0.13</td><td>1.00</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       " 1.00 -0.85  0.00 -0.02 -0.00\n",
       "-0.85  1.00 -0.05  0.02 -0.08\n",
       " 0.00 -0.05  1.00 -0.24 -0.15\n",
       "-0.02  0.02 -0.24  1.00 -0.13\n",
       "-0.00 -0.08 -0.15 -0.13  1.00"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of -log(likelihood) at the minimum:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>-log(likelihood)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sim_data</th>\n",
       "      <td>24.81213</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>24.81213</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          -log(likelihood)\n",
       "sim_data          24.81213\n",
       "total             24.81213"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of statistical measures:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>statistical measures</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AIC</th>\n",
       "      <td>60.987896</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BIC</th>\n",
       "      <td>69.184375</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     statistical measures\n",
       "AIC             60.987896\n",
       "BIC             69.184375"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ar.display(\"hpd\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The AR stores several properties from the analysis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:40.991872Z",
     "start_time": "2020-04-13T07:50:40.988368Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'MLE'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ar.analysis_type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:40.997560Z",
     "start_time": "2020-04-13T07:50:40.993754Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.47801527e-05, -3.69124668e-04,  2.79048679e-05,\n",
       "        -1.05067384e-05, -1.38912356e-06],\n",
       "       [-3.69124668e-04,  1.28732318e-02, -2.17478045e-02,\n",
       "         3.18981090e-04, -7.90107253e-04],\n",
       "       [ 2.79048679e-05, -2.17478045e-02,  1.55285025e+01,\n",
       "        -1.09924522e-01, -5.24718986e-02],\n",
       "       [-1.05067384e-05,  3.18981090e-04, -1.09924522e-01,\n",
       "         1.34279848e-02, -1.35822852e-03],\n",
       "       [-1.38912356e-06, -7.90107253e-04, -5.24718986e-02,\n",
       "        -1.35822852e-03,  7.61121218e-03]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ar.covariance_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.166278Z",
     "start_time": "2020-04-13T07:50:40.999086Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "WARNING UserWarning: get_point_source_flux() has been replaced by get_flux()\n",
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0aadddfba8454e5fb0c8891f433462ea",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HTML(value='Propagating errors : '), HTML(value=''), FloatProgress(value=0.0)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>flux</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake: total</th>\n",
       "      <td>(2.00 -0.21 +0.22) x 10^-5 erg / (cm2 s)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                 flux\n",
       "fake: total  (2.00 -0.21 +0.22) x 10^-5 erg / (cm2 s)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>flux</th>\n",
       "      <th>low bound</th>\n",
       "      <th>hi bound</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake: total</th>\n",
       "      <td>2.0042122454660008e-05 erg / (cm2 s)</td>\n",
       "      <td>1.7945578355340123e-05 erg / (cm2 s)</td>\n",
       "      <td>2.2218485618241006e-05 erg / (cm2 s)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                             flux  \\\n",
       "fake: total  2.0042122454660008e-05 erg / (cm2 s)   \n",
       "\n",
       "                                        low bound  \\\n",
       "fake: total  1.7945578355340123e-05 erg / (cm2 s)   \n",
       "\n",
       "                                         hi bound  \n",
       "fake: total  2.2218485618241006e-05 erg / (cm2 s)  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ar.get_point_source_flux(1*u.keV, .1*u.MeV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.194790Z",
     "start_time": "2020-04-13T07:50:53.168134Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "Model summary:<br><br><div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>N</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Point sources</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Extended sources</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Particle sources</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div><br><br>Free parameters (5):<br><br><div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>value</th>\n",
       "      <th>min_value</th>\n",
       "      <th>max_value</th>\n",
       "      <th>unit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>0.00877149</td>\n",
       "      <td>-10</td>\n",
       "      <td>10</td>\n",
       "      <td>s-1 cm-2 keV-2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>1.76234</td>\n",
       "      <td>-100</td>\n",
       "      <td>100</td>\n",
       "      <td>keV-1 s-1 cm-2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>32.0749</td>\n",
       "      <td>0.001</td>\n",
       "      <td>200</td>\n",
       "      <td>s-1 cm-2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>25.2136</td>\n",
       "      <td>0</td>\n",
       "      <td>100</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>1.09261</td>\n",
       "      <td>0.001</td>\n",
       "      <td>10</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div><br><br>Fixed parameters (2):<br>(abridged. Use complete=True to see all fixed parameters)<br><br><br>Linked parameters (0):<br><br>(none)<br><br>Independent variables:<br><br>(none)<br>"
      ],
      "text/plain": [
       "Model summary:\n",
       "==============\n",
       "\n",
       "                  N\n",
       "Point sources     1\n",
       "Extended sources  0\n",
       "Particle sources  0\n",
       "\n",
       "Free parameters (5):\n",
       "--------------------\n",
       "\n",
       "                                           value min_value max_value  \\\n",
       "fake.spectrum.main.composite.a_1      0.00877149       -10        10   \n",
       "fake.spectrum.main.composite.b_1         1.76234      -100       100   \n",
       "fake.spectrum.main.composite.F_2         32.0749     0.001       200   \n",
       "fake.spectrum.main.composite.mu_2        25.2136         0       100   \n",
       "fake.spectrum.main.composite.sigma_2     1.09261     0.001        10   \n",
       "\n",
       "                                                unit  \n",
       "fake.spectrum.main.composite.a_1      s-1 cm-2 keV-2  \n",
       "fake.spectrum.main.composite.b_1      keV-1 s-1 cm-2  \n",
       "fake.spectrum.main.composite.F_2            s-1 cm-2  \n",
       "fake.spectrum.main.composite.mu_2                keV  \n",
       "fake.spectrum.main.composite.sigma_2             keV  \n",
       "\n",
       "Fixed parameters (2):\n",
       "(abridged. Use complete=True to see all fixed parameters)\n",
       "\n",
       "\n",
       "Linked parameters (0):\n",
       "----------------------\n",
       "\n",
       "(none)\n",
       "\n",
       "Independent variables:\n",
       "----------------------\n",
       "\n",
       "(none)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ar.optimized_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving results to disk\n",
    "\n",
    "The beauty of the analysis result is that all of this information can be written to disk and restored at a later time. The statistical parameters, best-fit model, etc. can all be recovered.\n",
    "\n",
    "AR's are stored as a structured FITS file. We write the AR like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.263004Z",
     "start_time": "2020-04-13T07:50:53.196553Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: VerifyWarning: Card is too long, comment will be truncated. [astropy.io.fits.card]\n",
      "WARNING:astropy:VerifyWarning: Card is too long, comment will be truncated.\n"
     ]
    }
   ],
   "source": [
    "ar.write_to(\"test_mle.fits\", overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The FITS file can be examines with any normal FITS reader."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.266896Z",
     "start_time": "2020-04-13T07:50:53.264652Z"
    }
   },
   "outputs": [],
   "source": [
    "import astropy.io.fits as fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.273905Z",
     "start_time": "2020-04-13T07:50:53.268658Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filename: test_mle.fits\n",
      "No.    Name      Ver    Type      Cards   Dimensions   Format\n",
      "  0  PRIMARY       1 PrimaryHDU       6   ()      \n",
      "  1  ANALYSIS_RESULTS    1 BinTableHDU     36   5R x 8C   [36A, D, D, D, D, 16A, 5D, D]   \n"
     ]
    }
   ],
   "source": [
    "ar_fits = fits.open('test_mle.fits')\n",
    "ar_fits.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, to easily pull the results back into the 3ML framework, we use the ${\\tt load\\_analysis\\_results}$ function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.351831Z",
     "start_time": "2020-04-13T07:50:53.275494Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "WARNING YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ar_reloaded = load_analysis_results(\"test_mle.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.359571Z",
     "start_time": "2020-04-13T07:50:53.353442Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>-log(likelihood)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sim_data</th>\n",
       "      <td>24.81213</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>24.81213</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          -log(likelihood)\n",
       "sim_data          24.81213\n",
       "total             24.81213"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ar_reloaded.get_statistic_frame()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can get a DataFrame with the saved results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.388715Z",
     "start_time": "2020-04-13T07:50:53.361105Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>value</th>\n",
       "      <th>negative_error</th>\n",
       "      <th>positive_error</th>\n",
       "      <th>error</th>\n",
       "      <th>unit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>0.008771</td>\n",
       "      <td>-0.003853</td>\n",
       "      <td>0.003748</td>\n",
       "      <td>0.003801</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>1.762337</td>\n",
       "      <td>-0.113744</td>\n",
       "      <td>0.113748</td>\n",
       "      <td>0.113746</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>32.074949</td>\n",
       "      <td>-3.894620</td>\n",
       "      <td>3.847490</td>\n",
       "      <td>3.871055</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>25.213581</td>\n",
       "      <td>-0.112303</td>\n",
       "      <td>0.113055</td>\n",
       "      <td>0.112679</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>1.092612</td>\n",
       "      <td>-0.087786</td>\n",
       "      <td>0.084940</td>\n",
       "      <td>0.086363</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                          value  negative_error  \\\n",
       "fake.spectrum.main.composite.a_1       0.008771       -0.003853   \n",
       "fake.spectrum.main.composite.b_1       1.762337       -0.113744   \n",
       "fake.spectrum.main.composite.F_2      32.074949       -3.894620   \n",
       "fake.spectrum.main.composite.mu_2     25.213581       -0.112303   \n",
       "fake.spectrum.main.composite.sigma_2   1.092612       -0.087786   \n",
       "\n",
       "                                      positive_error     error  \\\n",
       "fake.spectrum.main.composite.a_1            0.003748  0.003801   \n",
       "fake.spectrum.main.composite.b_1            0.113748  0.113746   \n",
       "fake.spectrum.main.composite.F_2            3.847490  3.871055   \n",
       "fake.spectrum.main.composite.mu_2           0.113055  0.112679   \n",
       "fake.spectrum.main.composite.sigma_2        0.084940  0.086363   \n",
       "\n",
       "                                                  unit  \n",
       "fake.spectrum.main.composite.a_1      1 / (cm2 keV2 s)  \n",
       "fake.spectrum.main.composite.b_1       1 / (cm2 keV s)  \n",
       "fake.spectrum.main.composite.F_2           1 / (cm2 s)  \n",
       "fake.spectrum.main.composite.mu_2                  keV  \n",
       "fake.spectrum.main.composite.sigma_2               keV  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ar_reloaded.get_data_frame()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis Result Sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When doing time-resolved analysis or analysing a several objects, we can save several AR's is a set. This is achieved with the analysis result set. We can pass an array of AR's to the set and even set up descriptions for the different entries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.496786Z",
     "start_time": "2020-04-13T07:50:53.394131Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: VerifyWarning: Card is too long, comment will be truncated. [astropy.io.fits.card]\n",
      "WARNING:astropy:VerifyWarning: Card is too long, comment will be truncated.\n"
     ]
    }
   ],
   "source": [
    "from threeML.analysis_results import AnalysisResultsSet\n",
    "\n",
    "analysis_set = AnalysisResultsSet([ar, ar_reloaded])\n",
    "\n",
    "# index as time bins\n",
    "analysis_set.set_bins(\"testing\", [-1, 1], [3, 5], unit=\"s\")\n",
    "\n",
    "# write to disk\n",
    "analysis_set.write_to(\"analysis_set_test.fits\", overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.641211Z",
     "start_time": "2020-04-13T07:50:53.499639Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "WARNING YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.\n",
      "\n",
      "\n",
      "WARNING YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "analysis_set = load_analysis_results(\"analysis_set_test.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.684275Z",
     "start_time": "2020-04-13T07:50:53.643120Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best fit values:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>result</th>\n",
       "      <th>unit</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>parameter</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>(9 +/- 4) x 10^-3</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>1.76 +/- 0.11</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>(3.2 +/- 0.4) x 10</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>(2.521 +/- 0.012) x 10</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>1.09 +/- 0.09</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                      result              unit\n",
       "parameter                                                                     \n",
       "fake.spectrum.main.composite.a_1           (9 +/- 4) x 10^-3  1 / (cm2 keV2 s)\n",
       "fake.spectrum.main.composite.b_1               1.76 +/- 0.11   1 / (cm2 keV s)\n",
       "fake.spectrum.main.composite.F_2          (3.2 +/- 0.4) x 10       1 / (cm2 s)\n",
       "fake.spectrum.main.composite.mu_2     (2.521 +/- 0.012) x 10               keV\n",
       "fake.spectrum.main.composite.sigma_2           1.09 +/- 0.09               keV"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Correlation matrix:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table id=\"table4718696848\">\n",
       "<tr><td>1.00</td><td>-0.85</td><td>0.00</td><td>-0.02</td><td>-0.00</td></tr>\n",
       "<tr><td>-0.85</td><td>1.00</td><td>-0.05</td><td>0.02</td><td>-0.08</td></tr>\n",
       "<tr><td>0.00</td><td>-0.05</td><td>1.00</td><td>-0.24</td><td>-0.15</td></tr>\n",
       "<tr><td>-0.02</td><td>0.02</td><td>-0.24</td><td>1.00</td><td>-0.13</td></tr>\n",
       "<tr><td>-0.00</td><td>-0.08</td><td>-0.15</td><td>-0.13</td><td>1.00</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       " 1.00 -0.85  0.00 -0.02 -0.00\n",
       "-0.85  1.00 -0.05  0.02 -0.08\n",
       " 0.00 -0.05  1.00 -0.24 -0.15\n",
       "-0.02  0.02 -0.24  1.00 -0.13\n",
       "-0.00 -0.08 -0.15 -0.13  1.00"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of -log(likelihood) at the minimum:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>-log(likelihood)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sim_data</th>\n",
       "      <td>24.81213</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>24.81213</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          -log(likelihood)\n",
       "sim_data          24.81213\n",
       "total             24.81213"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of statistical measures:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>statistical measures</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AIC</th>\n",
       "      <td>60.987896</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BIC</th>\n",
       "      <td>69.184375</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     statistical measures\n",
       "AIC             60.987896\n",
       "BIC             69.184375"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "analysis_set[0].display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error propagation\n",
    "In 3ML, we propagate errors for MLE reults via sampling of the covariance matrix *instead* of Taylor exanding around the maximum of the likelihood and computing a jacobain. Thus, we can achieve non-linear error propagation.\n",
    "\n",
    "You can use the results for propagating errors non-linearly for analytical functions:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.697642Z",
     "start_time": "2020-04-13T07:50:53.686151Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Propagating a+b, with a and b respectively:\n",
      "equal-tail: (9 +/- 4) x 10^-3, hpd: (9 +/- 4) x 10^-3\n",
      "equal-tail: 1.76 +/- 0.11, hpd: 1.76 -0.12 +0.11\n",
      "\n",
      "This is the result (with errors):\n",
      "equal-tail: 1.77 +/- 0.11, hpd: 1.77 -0.12 +0.10\n",
      "(1.6607512376061622, 1.8806969631644799)\n"
     ]
    }
   ],
   "source": [
    "p1 = ar.get_variates(\"fake.spectrum.main.composite.a_1\")\n",
    "p2 = ar.get_variates(\"fake.spectrum.main.composite.b_1\")\n",
    "\n",
    "print(\"Propagating a+b, with a and b respectively:\")\n",
    "print(p1)\n",
    "print(p2)\n",
    "\n",
    "print(\"\\nThis is the result (with errors):\")\n",
    "res = p1 + p2\n",
    "print(res)\n",
    "\n",
    "print(res.equal_tail_interval())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The propagation accounts for covariances. For example this\n",
    "has error of zero (of course) since there is perfect covariance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.704012Z",
     "start_time": "2020-04-13T07:50:53.699147Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "This is 50 * a/a:\n",
      "equal-tail: (5.0 +/- 0) x 10, hpd: (5.0 +/- 0) x 10\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "WARNING RuntimeWarning: divide by zero encountered in log10\n",
      "\n",
      "\n",
      "WARNING RuntimeWarning: divide by zero encountered in log10\n",
      "\n",
      "\n",
      "WARNING RuntimeWarning: divide by zero encountered in log10\n",
      "\n",
      "\n",
      "WARNING RuntimeWarning: divide by zero encountered in log10\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"\\nThis is 50 * a/a:\")\n",
    "print(50 * p1/p1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-01-04T15:57:15.657758Z",
     "start_time": "2018-01-04T15:57:15.648995Z"
    }
   },
   "source": [
    "You can use arbitrary (np) functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.711828Z",
     "start_time": "2020-04-13T07:50:53.705707Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "This is arcsinh(a + 5*b) / np.log10(b) (why not?)\n",
      "equal-tail: (1.17 -0.09 +0.12) x 10, hpd: (1.17 -0.11 +0.10) x 10\n"
     ]
    }
   ],
   "source": [
    "print(\"\\nThis is arcsinh(a + 5*b) / np.log10(b) (why not?)\")\n",
    "print(np.arcsinh(p1 + 5 * p2) / np.log10(p2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-01-04T15:41:58.869808Z",
     "start_time": "2018-01-04T15:41:58.835016Z"
    }
   },
   "source": [
    "Errors can become asymmetric. For example, the ratio of two gaussians is\n",
    "asymmetric notoriously:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.718090Z",
     "start_time": "2020-04-13T07:50:53.713718Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Ratio a/b:\n",
      "equal-tail: (2.0 -0.7 +1.6) x 10^2, hpd: (2.0 -1.0 +0.7) x 10^2\n"
     ]
    }
   ],
   "source": [
    "print(\"\\nRatio a/b:\")\n",
    "print(p2 / p1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-01-04T15:41:58.869808Z",
     "start_time": "2018-01-04T15:41:58.835016Z"
    }
   },
   "source": [
    "You can always use it with arbitrary functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.724738Z",
     "start_time": "2020-04-13T07:50:53.719704Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Propagating using a custom function:\n",
      "equal-tail: 1.77 +/- 0.11, hpd: 1.77 +/- 0.11\n"
     ]
    }
   ],
   "source": [
    "def my_function(x, a, b):\n",
    "\n",
    "    return b * x ** a\n",
    "\n",
    "\n",
    "print(\"\\nPropagating using a custom function:\")\n",
    "print(my_function(2.3, p1, p2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is an example of an error propagation to get the plot of the model with its errors\n",
    "(which are propagated without assuming linearity on parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:53.735893Z",
     "start_time": "2020-04-13T07:50:53.726379Z"
    }
   },
   "outputs": [],
   "source": [
    "def go(fitfun, ar, model):\n",
    "\n",
    "    fig, ax = plt.subplots()\n",
    "\n",
    "    # Gather the parameter variates\n",
    "\n",
    "    arguments = {}\n",
    "\n",
    "    for par in fitfun.parameters.values():\n",
    "\n",
    "        if par.free:\n",
    "\n",
    "            this_name = par.name\n",
    "\n",
    "            this_variate = ar.get_variates(par.path)\n",
    "\n",
    "            # Do not use more than 1000 values (would make computation too slow for nothing)\n",
    "\n",
    "            if len(this_variate) > 1000:\n",
    "\n",
    "                this_variate = np.random.choice(this_variate, size=1000)\n",
    "\n",
    "            arguments[this_name] = this_variate\n",
    "\n",
    "    # Prepare the error propagator function\n",
    "\n",
    "    pp = ar.propagate(\n",
    "        ar.optimized_model.fake.spectrum.main.shape.evaluate_at, **arguments\n",
    "    )\n",
    "\n",
    "    # You can just use it as:\n",
    "\n",
    "    print(pp(5.0))\n",
    "\n",
    "    # Make the plot\n",
    "\n",
    "    energies = np.linspace(0, 50, 100)\n",
    "\n",
    "    low_curve = np.zeros_like(energies)\n",
    "    middle_curve = np.zeros_like(energies)\n",
    "    hi_curve = np.zeros_like(energies)\n",
    "\n",
    "    free_parameters = model.free_parameters\n",
    "\n",
    "    with progress_bar(len(energies), title=\"Propagating errors\") as p:\n",
    "\n",
    "        with use_astromodels_memoization(False):\n",
    "\n",
    "            for i, e in enumerate(energies):\n",
    "\n",
    "                this_flux = pp(e)\n",
    "\n",
    "                low_bound, hi_bound = this_flux.equal_tail_interval()\n",
    "\n",
    "                low_curve[i], middle_curve[i], hi_curve[i] = (\n",
    "                    low_bound,\n",
    "                    this_flux.median,\n",
    "                    hi_bound,\n",
    "                )\n",
    "\n",
    "                p.increase()\n",
    "\n",
    "    ax.plot(energies, middle_curve, \"--\", color=\"black\")\n",
    "    ax.fill_between(energies, low_curve, hi_curve, alpha=0.5, color=\"blue\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:59.011748Z",
     "start_time": "2020-04-13T07:50:53.737722Z"
    },
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "equal-tail: 1.81 -0.11 +0.12, hpd: 1.81 -0.09 +0.14\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "37c1e337f0a84225bd8e9cc748c045ea",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HTML(value='Propagating errors : '), HTML(value=''), FloatProgress(value=0.0)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAAIYCAYAAAA7EQvxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeZhcZZ33/8+pqq6ll+pOd2clOwQCGDZZXCAgIBpQRJ1xQ1FadPg5wOM44/LMODKjj+AMDA7qD5QlxIUHWYQBNIAoIgjKEoiBAAGykH3rTnfXXqfqnOePk9p6SXqpvd6v68rVXdWnzrnDpc2H733f39uwbdsWAAAAaoKr0gMAAADA2BHeAAAAagjhDQAAoIYQ3gAAAGoI4Q0AAKCGEN4AAABqSNWFt6VLl1Z6CAAAAFWr6sIbAAAARkd4AwAAqCGENwAAgBpCeAMAAKghhDcAAIAaQngDAACoIYQ3AACAGkJ4AwAAqCGENwAAgBpCeAMAAKghhDcAAIAaQngDAACoIYQ3AACAGkJ4AwAAqCGENwAAgBpCeAMAAKghhDcAAIAaQngDAACoIYQ3AACAGkJ4AwAAqCGENwAAgBpCeAOAUViW8wcAqgnhDQBGsX69tG1bpUcBAIUIbwAwiu3bpd27Kz0KACjkqfQAAKBa7dtX6REAwHCENwAYRW+v5PVWehQAUIhpUwAYRV+fFA5XehQAUIjwBgAjSCalgQHCG4DqQ3gDgBFs2eJ8jUQqOw4AGIrwBgAj2L5d8niovAGoPoQ3ABhBb6/zNRqV0unKjgUA8hHeAGAEfX3OV9OUBgcrOxYAyEd4A4ARZMJbU5O0a1dlxwIA+QhvADCEaTo7TSXJ7Zb6+ys7HgDIR3gDgCG2bi08kJ5NCwCqCeENAIbYts2ZLs0gvAGoJoQ3ABgis94tg15vAKoJ4Q0Ahhga3qi8AagmhDcAGILKG4BqRngDgDymKe3bV/heJCLZdmXGAwBDEd4AIM/27cNPVIjFpESiMuMBgKEIbwCQZ+tWyestfM+2h1fjAKBSCG8AkGfoejfJCXO7d5d/LAAwEsIbAOQZKbwZBuebAqgehDcAyDNSeJNoFwKgehDeAGC/VGr0tW20CwFQLQhvALDfzp1OgBsJlTcA1YLwBgD7bdlSeKZpPipvAKoF4Q0A9uvtdTYnjITKG4BqQXgDgP1G26wgOZU3yyrfWABgNIQ3ANjvQOHNNJk6BVAdCG8AIOdIrAOFN7ebRr0AqgPhDQAk7do1+k5TyQlvvb3lGw8AjIbwBgCSNm8efaep5GxkYNMCgGpAeAMASQMDo+80zWDNG4BqQHgDAEmJxMGvofIGoBoQ3gBAzm7Sg6HyBqAaEN4AQFIyefBrqLwBqAaENwDQ2CpvhDcA1YDwBgAaW+UtHh9byAOAUhpzeLvsssu0du3aYe/ffPPN+qd/+qeiDgoAym0s4S2Vkvr7Sz8WADiQg4a3Z599Vtdcc43WrFkz7Gevv/667rzzzpIMDADKaSzhzeuVdu4s/VgA4EAOGt7WrFkjy7IUCAQK3k+lUrrmmmv0wQ9+sGSDA4ByGct0qMslDQ6WfiwAcCCeg11wySWXSJJWrVpV8P7Pf/5zvf3tb9e8efO0ZcuW0owOAMpkrGvZ2LQAoNImtGFh/fr1evzxx3XxxRcXezwAUBFjmTaV6PUGoPIOWnkbyrIs/cd//If+4R/+QT6fb8IPXr58uVasWDHhzwNAsaRSzh/PGH4jEt4AVNq4w1skEtFrr72mK664ouD9c889VytXrhzzfXp6etTT0zPs/aVLl453SAAwKcmkZFljuzYUKu1YAOBgxh3e2tra9MQTT2RfP/TQQ/r973+va6+9tqgDA4BySSSkdHps11J5A1BpNOkF0PDCYcntHtu1kYhk26UdDwAcyJgrb3fdddeI7y9btkzLli0r2oAAoNzGE94SCSkalVpaSjsmABgNlTcADS8adXq4jYXLJe3dW9rxAMCBEN4ANDzTlAxjbNd6PIQ3AJVFeAPQ8MZz2LxhsOMUQGUR3gA0vLE26M2IxUozDgAYC8IbgIY33vA23usBoJgIbwAa3njD2HimWQGg2AhvABreeMNYKlWacQDAWBDeADQ8pk0B1BLCG4CGN97KG+ENQCUR3gA0PCpvAGrJuA+mB4B6M1Ll7dprpysS2a2ZM0/QF7+4quBnhDcAlUTlDUDDGxrGbNtSJLJbkhQKbR92PRsWAFQS4Q1Awxsa3vr738p+H48PyLbtA14PAOVEeAPQ8IZOm+7d+1r2+1Qqplist+DnhDcAlUR4A9DQbHt4GMsPb5I0MLC54HUqJaXTpR4ZAIyM8AagoY0UxA4W3iyL6huAyiG8AWhoyeTw8Nbb64S3zs7DJBHeAFQXWoUAaGiJhDN1mu9tb/ukpkw5THPnvluDg1s1e/Y7Cn7ucknhsNTeXsaBAsB+hDcADS0UcsJYvhNPvFQnnnjpqJ9xu6VIpMQDA4BRMG0KoKFFIpJnnP8Z63JJ0WhpxgMAB0N4A9DQolHJMHKvd+1ao3XrHlAotF2WldIbbzyk1at/WvAZw2DNG4DKIbwBaGimWRje1qz5hX75yw/pxReXS5LuuOODuv/+zymVShR8jvAGoFIIbwAa2mgNeru7F8vl8igYPESSFAptK7iO8AagUghvABraaA16u7sXS5La2+dKGt4uZKTD7AGgHAhvABpafghLpRLat2+DDMOlzs5FkkYPb1TeAFQK4Q1AQ0vkLWXr63tTtp3WlCkL5fH4JEnBIJU3ANWF8AagoeWHsKFTphKVNwDVh/AGoKHlh7DBwa2SpK6u4eEtGt1b8DkqbwAqhRMWADS0/BD2jnf8L51wwiVKp3OJbuHCs/SNbwzI5wsWfI7KG4BKIbwBaGhDQ5jX2yKpJfva4/HL4/Ef9HMAUC5MmwJoaBOd/mTaFEClEN4ANLRMCBsc3Kbvf3+O7rnn48OuefDBv9OPfnSEdu7867DPAUC5Ed4ANLRMq5C9e1/V4OBWhUI7hl0TCm1Vb+/r6u/flH2PaVMAlUJ4A9DQMhW0kdqEZIzU643KG4BKIbwBaFi2naugHSi8jdTrLZWSLKv0YwSAoQhvABpWfgAbS3gbHMyFt3SaqVMAlUF4A9CwEgknhEnjr7wR3gBUCuENQMOKx52p00RiUKHQNrndPrW3zxt23UjhzTCkWKxsQwWALJr0AmhYoZDkdkvptK33vvdaJRKDcrncw64LBg/R8cdfoo6O+bJtW4ZhyO12Pj99egUGDqChEd4ANKxIxAlvTU3tete7/nHU61wuj84//+aC99xu5/MAUG5MmwJoWPG45Jrgb0HDYM0bgMqg8gagYWXC16uv3ivJ0IIF75Hf3zHitf39b2n37pfV1bVIXV2HS8o1+AWAcqLyBqBhZRrtPv74v+muuz6iffs2jHrtc8/doDvu+IDWrr172OcBoJwIbwAaVqbyFov1SZKam6eOeu1IO04JbwAqgfAGoGFlwlciMSBJ8vmCo147UqNe1rwBqATCG4CGlUxKlpVWMhmWZMjnaxv1WipvAKoF4Q1Aw0omnQa9kuTztckwRv+V2NHhNO8dGNgs27YlEd4AVAbhDUDDcsJbZsq0/YDX+nzt8nrblEyGFY/3Zz8PAOVGeAPQsEwzV3nz+w8c3gzDyFv3tkUS4Q1AZdDnDUDDMk1p+vRj9M1vJmWa0YNef+GFDykQmCKvtzX7eQAoN8IbgIaVqZy53U1yuw9ceZOk9vY5Ba8JbwAqgWlTAA1rsuGLaVMAlUB4A9CwkknppZfu0C23vEPPPPPDg16/Zcuf9bOfnaVHHnEOsafyBqASmDYF0LCSSam/f6O2bXtG8+efcdDrLcvUxo2PyTRj2c8DQLlReQPQkGzbCV/x+NhahUhSa+sMSVI0uleSU3nb3/INAMqG8AagIZmmZFm5Pm8HaxUi5Y7PyrQXSaeZOgVQfmMOb5dddpnWrl0rSYrH47rqqqv0/ve/XxdccIGuv/56pVKpkg0SAIotkXDC11ib9DrXFIY3J/yVbowAMJKDhrdnn31W11xzjdasWZN97/bbb9emTZt022236aabbtLatWt13333lXSgAFBM8Xjm69grbx5PQIbhVioVUzptyrZz9wGAcjloeFuzZo0sy1IgEMi+9/zzz+tTn/qUZs6cqWnTpun888/XqlWrSjpQACimUEhyu8dXeTMMI1t9SyZDcrulcLikwwSAYQ662/SSSy6RpIJw9rWvfU3Tp0/Pvl63bp2mTp1aguEBQGmEw054W7ToPE2ZslDB4Oxh16RSzoaEpqbce2972yeUTjvbTN1uKRIp14gBwDGhViELFiyQ5Kx9u+mmm/Tkk0/qxhtvHNc9li9frhUrVkzk8QAwaYmE5HJJp532z6NeEwgMbwdy3nk3FLyOxUoxOgAY3YT7vL3wwgu6+uqrdfjhh+uWW25Rd3f3uD7f09Ojnp6eYe8vXbp0okMCgDEbS4+2qVOd6dVQaPRr2G0KoNwmFN6ef/55XXnllfrGN76h0047rdhjAoCSc3q0Wdq27Vn5/R3q7l487JqpUyXDKAxvodAORSK71d4+R4FAJ+ENQNlNqM/b8uXLdfnllxPcANSsZNJp+XHrre/UzTefPOzntu2Et/Yh+xgeeeQf9JOfHKc333wkex8AKKcJVd7Wr1+vq666SldddVX2veOOO04/+MEPijYwACgl0zxwm5BkUlq4cPiU6dBeb1TeAJTbmMPbXXfdlf3+kUceKclgAKBcnMpbJrx1DPu5zyfNmCFt2jT0/cLwRuUNQLlxPBaAhpRfeRupx9u0ac5u1O5up2VIRn6fN+dr6ccKAPkIbwAaUiJx4HNNMxvoDznEOQYrg2lTAJVGeAPQkA5Wecv0HW9udvq9ZRDeAFQa4Q1AQzLN0Y/GSqWk2fsPXDAMqa0t9zPWvAGotAk36QWAWpZMOkddzZ79jmHhzTCk+fNzr9vbpf5+5/v5889QT89Tam2dKYnKG4DyI7wBaEimKQUCnQoEOof9rKur8DzTYDD3fXNzt5qbcyfKEN4AlBvTpgAa0oGmO4ee9je0UW++RKI44wGAsSK8AWhIpik988wP9eCDX9SOHS8W/CyzWSEjGHROXJCctW4PPXSFHn74y9n7AEA5Ed4ANBzbdipv69c/rBdeuFmh0LbszyzLac6bb+bMwpD27LM/1Isv3iopc0ZqOUYNAA7CG4CGk0w6IS0ed3Yh5G9YSKelQw8tvH7aNGcTgyR5va377xGWZaWVShU28QWAUiO8AWg4iYQT0kY62zQYlFpbC6/3+aSWFud7w3DJ63V6hyST4WwVDwDKhfAGoOHEYs7Xkfq8DV3vlpG/aSG/15tlsWkBQHkR3gA0nHBY8nhGrryNFt7y24XkhzeXy7kfAJQL4Q1Aw4lEJMNIZw+Xz0yD2razvm0ko4U3j4fwBqC8aNILoOHE45JlJTR37mmyrJRcLrckZ+3a0M0KGfnTplOnHi3btuRyeWQYuWlYACgHwhuAhpNMSk1Nzbr44icK3m9udk5XGEl7u1OZMwzpQx+6teBn9HoDUE5MmwJoOKOFrfyWIEPNnj36xgR2mwIoJ8IbgIaTTErpdFLJZER2Xofd0TYrSFJnZ+F5p5aVVjrtpDYqbwDKifAGoOGYprR+/W919dWtuuOOD2bfP1B4c7ulNmdfg/74x2/rO9/x6Mknr5ZE5Q1AeRHeADScZDLXJsTncxKZaUpz5x74c5kdp01NzZKc3aaZ+wFAuRDeADQc0xzeoNfjcda1HUhmx2l+q5DM/QCgXAhvABpOYeXNSWSdnc7U6IEMDW/JJOENQPkR3gA0nGQyV3nLnK4QCBz8c5lp06GVN6ZNAZQT4Q1AwzHN4ZU3n+/gn+vsdA60Z9oUQCUR3gA0nPw1b5nK21jC2yGHSKkUlTcAlcUJCwAaTjIpnXzy5Tr00HM0Z867JI0tvAWDznXt7XN13nk/VlvbzOz9AKBcCG8AGo5pSnPmvFNz5rwz+57ff/DPGUZm3VuHTjzx7wruBwDlwrQpgIYzUtjyesf22fwD6g90PwAoFSpvABqKbTvTnM8//0MZhqHjjrtYXm/LmKZNpdyO09Wrf6podK9OPvnvJY2hbAcARUJ4A9BQTFOyLOkPf/imEolBLVlyoaTxh7ff//4bCod3asmSTykQmKl0+uB94gCgGJg2BdBQnEPpLSUSIUnOztF0WmppGdvnM9OmXq9zrFYiMSjLkhKJUowWAIYjvAFoKImE9gc3W15vq1wut1KpkdeyjWTqVKd6l98uxLLYcQqgfAhvABpKOCyZZmGDXklqaxvb52fNctbN5Yc3l0uKRIo+VAAYEeENQEOJRHLhLdOg1+Ua2/FYzmeca/PDm9vthEIAKAfCG4CGEotJyWRh5a2pSfKMcftWptcblTcAlUJ4A9BQnN2mSfn9HQoEOiU5pyYYxtjv4fc74c3j8SudTsow6PUGoHxoFQKgoZimtGDBmfr61/fJtm1JYzsaK5/PJy1b9gOde+6PCu4LAOVA5Q1AQ8nfFWrsL7eNN7x5vZJhFP76pFUIgHIhvAFoKCNVyCZSeRvLfQGgFAhvABqKaUpPPPFd/eAHh+qFF26VNLHw9uabD+uGG47WQw9dkb0vAJQD4Q1AQzFNKRTapn37NiiVikuaWHhLp5Pas+cV9fdvlESTXgDlQ3gD0FBMU0okCvu8eb3ju4fPV9gqJHNfACgHwhuAhpJMSvF4YZ+3iWxYILwBqBTCG4CGkkwOr7z5/eO7h98veb2F4Y1pUwDlQngD0FBMc/KVt9ZWyeMhvAGoDMIbgIYydM2bc8j8+O7R3i653U54i8cHZNs206YAyoYTFgA0FNOUTjrpMoXDO9Tc3C3LktraxneP5mapqcmnk066TD5fm2zbkmm6SzNgABiC8AagoZim9O53fzX7OpFwKmnj4fNJbrehc8/9YcF9AaAcmDYF0FBGWpvW0jK+e3g8UlNT4XuENwDlQngD0DBsW4pEwlq37gFt3fqMJMnlGv9uU8mpvu3a9ZLWr/+tYrF9SiYlyyrygAFgBIQ3AA0jnZZ6ezfql7/8kB54oEeS07PNM4EFJF6vtHLl3+sXv3ifdu1aI8ui+gagPAhvABpGMinFYpNrE5Ix9JQFy6JdCIDyILwBaBiJRK7H20Qb9GYMDW+2LcXjRRkmABwQ4Q1AwwiHJdMsrLyN91zTDCe8OfdIJAbldjv3B4BSI7wBaBiRyPDwNtFp06Hnm7rdzv0BoNQIbwAaRjQ6/FzTYq15c7mYNgVQHmMOb5dddpnWrl0rSbJtWz/60Y+0bNkyffjDH9Z9991XsgECQLHkH41V7A0LmfsDQKkddIP8s88+qz/+8Y9as2ZN9r3f/OY3eu6553TbbbcpFArpK1/5io466igdccQRJR0sAEyGaUqnn/6vOumkL8nrdTrzTmbadMmST2nRomVqbu7O3h8ASu2g4W3NmjWyLEuBQCD73qOPPqoLL7xQM2bM0IwZM3TWWWfpscceI7wBqGrJpNTU1Kz29ubse5OpvDU3d6m5uavg/gBQagcNb5dccokkadWqVdn31q9fr0WLFmVfL1iwQM8991wJhgcAxTNSZWwylbex3B8Aim1CB9NHIhG15B0GGAgEFI1Gx3WP5cuXa8WKFRN5PABMiGlKDz/8ZQ0MvKWzzvqeuruPmHCft0BA6uvbpMce+4ZaWqZr2bLrqbwBKIsJhbfW1lYlEons60Qioba2tnHdo6enRz09PcPeX7p06USGBAAHZZrSxo2Pafful3T66VfKtifepDcYlBKJuNauvVOdnYskXU/lDUBZTKhVyNy5c7V+/frs602bNumwww4r2qAAoBSG7jZNp6Vx/ndnVjAoeTzsNgVQfhMKb2effbbuuOMO7du3Ty+99JIeffRRnXHGGUUeGgAUVzJZeDxWKiW1t0/sXoGA1Nwc3H/fkCTCG4DymNC06fnnn6+NGzfqE5/4hFpbW3X55Zdrzpw5xR4bABRVImFlq2Q+X1CplNTcfJAPjcI5YaFFkiHTjMqyUjLNCf1KBYBxGfNvmrvuuiv7vdvt1le+8hV95StfKcmgAKAUwuGwJFtNTS1yuTxyuye+5s3lknw+Qz5fUInEgBKJkJLJKUUdLwCMhOOxADSMUKjwaCyv1wlhE+X3F56ywLQpgHKgxg+gYaRSLh155Efk9zsVspF6tY1HU5O0YMF7FI32yjBchDcAZUF4A9AwmpsP0cc+9qvs64k26M3//AUX/DT7mj5vAMqBaVMADWNouJroereMoeGPyhuAciC8AWgY4fCg+vs35e04ndz9vF4plYorHN6lZDIs05QsqwgDBYADILwBaAjptLR27b26/voFWrnyMkmTD29+v7Ry5eX6r/+aoZde+r+yLKpvAEqP8AagITgNep1mupkdosWovOXvNrUs1r0BKD3CG4CGkEzmTkLwep0zsYqxYSE/vNm2FI9P7p4AcDCENwANIRLJhTefrzjhbWjlze2WwuHJ3RMADobwBqAhRCKSaZa28uZ2O88BgFIivAFoCPnhrViVt6EnLLhcTJsCKD3CG4CGMHTNm21LgcDk7tncLDU1OUEw036EDQsASo0TFgA0BNOUzjjj33XCCZdo5swTlEpJ7e2Tu2d7u9TdfYI+9rF7FQzOzj4HAEqJ8AagIZimNH36Ek2fvkSSM70ZDE7unq2tUkvLNB155IcLngMApcS0KYCGMDRUud2Tnzb1+yXXkN+ihDcApUblDUBDSCalxx//N1lWWu9+91fV1BSc9IYFj0dyuUw9+eQ1SqXies97vs2aNwAlR+UNQEMwTekvf/lvPfnk/5FlpeT1SoYxuXsahuTzufTYY/+iJ574P7Jti8obgJIjvAFoCImEXbDbdLJVt4xAwC2vt1WSnT2cHgBKifAGoCFEIjHZtiW32ye3u6lo4W1oo17CG4BSI7wBaAihUHEb9GYQ3gCUG+ENQEMIh4t7NFaGz5e7ZyIRUiJRnPsCwGgIbwAaQia8Zapkfn9x7tvUROUNQHkR3gA0hGTSVmfnIk2ZskBScStvbW2z1N4+VxJ93gCUHn3eADSEGTNO0OWXv559Xazw5vdLH/7wz7KvCW8ASo3KG4CGkEoVvi5WePN6C18T3gCUGuENQEMYevJBMadN85mmZFnFuTcAjITwBqDuWZb01FP/rauvbtMf/vAtScMrZhPl9UrPPXeDrrlmqh577F9lWVTfAJQW4Q1A3UsmpXh8QMlkWLZty7Kk5ubi3NvvlyzLUjS6V7FYnyxreJUPAIqJ8Aag7jnhLdekN5WS2tuLc+/WVqmpqW3/cwZl21IsVpx7A8BICG8A6l4sJplmrs+bZUnBYHHu3d4ueTyZPm8hud1SOFycewPASAhvAOpeOJwLb15vm9xuKRAozr2bm3NHbiUSg4Q3ACVHeANQ9yIRKZnMTZt6PM7JCMXg80l+v1N5SyZDcrnEEVkASorwBqDuxeO58Ob1tsnnkwyjOPd2u6XW1tzZphK7TQGUFicsAKh7pimdfPJlOuywZerqWlS0Hm8ZXV0zdfrp/6a2tlnZ5wFAqRDeANQ905SOOupvsq+LHd7a2zt0xhlXFjwPAEqFaVMAdW9omCp2eBvplAUAKBXCG4C6Z5rS6tUrtHbtXbKsdEnC2xtvPKQ1a25XKhWnSS+AkmLaFEDdi8fTuv/+iyVJ3/pW8cOb1ys98MDnFQ7v0Pz5p8s0Zxf3AQCQh8obgLoXCjmN17zeVhmGqySVt1yvtxDTpgBKivAGoO6FQrk2IVJp1rz5fLleb0ybAiglwhuAuhcO5xr0Ss40ZzH5fLlgmEgMUnkDUFKENwB1r9SVt6amXOUtkaDyBqC0CG8A6l4mvJWy8pZ/vimVNwClRHgDUPcikVzlzbaLX3nz+yWvN7fmjfAGoJRoFQKg7h155AX6l3+Jy7JMpdNSa2tx7x8ISGeeebXOOusqeb2thDcAJUV4A1D3TFPyeHySfIrHpWCwuPcPBiW3O5idjiW8ASglpk0B1L38MGUYTqWsmIJBybIKn5f/GgCKifAGoK5ZlvT00z/Srbe+S2vW3C6XqzRr3jZvfkw//emZevzxf5NliR2nAEqG8Aagrpmm1Nv7hrZu/bOi0T1qapI8RV4w4vVKptmvTZv+oF27/qp0mvAGoHQIbwDqWjLp9F6TnN2mXq8zdVpMhiG1tmb6vA3KMKRYrLjPAIAMwhuAuhaPO4FKcnqxFbvHW0Zra+5sU7dbCodL8xwAILwBqGvhsGSahZW3UsiEt2QyJJdLikRK8xwAILwBqGuRSC68lbLyFgzmpk1dLimRKM1zAIDwBqCuJRJONUwqbeWtvT03bSrR6w1A6dCkF0BdM03piCMu0PTpx6qlZVrJwltbW6uOPPKj8vnaZdu2TLPIuyIAYD/CG4C6ZprSmWd+J/u6VOEtEHDrYx+7p+C5AFAKkwpvd9xxh+655x6FQiEtXrxY//iP/6h58+YVa2wAMGlD+62VKrwNbfxLeANQKhNe87Z27Vrdeeed+u///m898MADWrx4sa677rpijg0AJi0WS2vbtme1d+86ScU/XSHD65UGB7dq1641Ms0oTXoBlMyEw5vb7ZYkpdNpSZJlWdndVgBQLfr6+nTLLafo1lvfKal0lTe/X/rlLy/Qj398rHbvXkt4A1AyE542Xbx4sU499VRddNFFkiSfz6fly5cXbWAAUAyhUK5NiPO1NM/xeiWfL9cuhGlTAKUy4fC2atUqPfnkk7rppps0e/Zs/eQnP9G3v/1t3XzzzWP6/PLly7VixYqJPh4AxiQX3pxgVarKm9frtCKRnNYkhDcApTLh8Pbkk0/qnHPO0eLFiyVJl156qZYtW6bBwcExTZ/29PSop6dn2PtLly6d6JAAYJhwONfjzbZLV3lrbs6FNypvAEppwmvefEN+A7rdbrnd7mHvA0AlRSLOIaM+X5vSaam1tTTPCQalpqbMtGmINW8ASmbC4e1d73qXfvvb3+qVV15RNBrVLbfcolNOOYXwBqCq5FfeUiknZJVCayvTpgDKY8LTpscee6wuvfRSfec739G+fft03HHH6Wtf+1oxx7NIyhIAACAASURBVAYAkxaJ5DYsGIYUCJTmOT6f5PezYQFA6U2qSe+yZcu0bNmyYo0FAIruqKMu0Be+cIx8vna53aVb8+bzSUuWXKiFC09XR8d8whuAkuF4LAB1zePp1KxZnZIk25Y8JfqtZxjStGnzlU7Pl8QJCwBKZ8Jr3gCgFuSHKK/XCVmlkt+GxDQlyyrdswA0LsIbgLpl29Jzz92mBx/8O23e/FTJerxlhELr9fDDX9bTT18ryxp+rioAFAPhDUDdSqWkDRt+rxdeuEn79q0veXhLJvfomWeu1yuv3K10mvAGoDQIbwDqViLh9FyTnDYepe5k1NqaadIbkmFIsVhpnwegMRHeANSteNzpuSY5rUJKXXnLnC6TSAzK7ZbC4dI+D0BjIrwBqFsDA5Jp5ipvpQ9vuSa9brcUiZT2eQAaE+ENQN0KhcpbeWtvd87ecqZqbSUSpX0egMZEeANQt/KnTctReQsEPGpqapZkyzQj9HoDUBI06QVQt5JJadq0JWpu7pbf317yDQterzRjxnFKpRJKpRIyzdbSPhBAQyK8AahbyaT06U8/nH1d6sqbzyf19DxV8HwAKDamTQHUraFrzkod3oben2lTAKVAeANQt+JxS6lUXLZtSypP5c22Jdu2ZFlpKm8ASoLwBqBuvfXWm/rudwO64YajZNsq+Zq3tjbp7rs/rm9/263XX/81lTcAJUF4A1C3BgednaYeT0DptNRa4v0DwaBkGM5S4mQyRHgDUBKENwB1KxQalOT0eEulnHBVSoGA5PPlTlkgvAEoBcIbgLoVCuV6vBmG1Nxc2uf5fE5QlJxGvYQ3AKVAeANQtyKR3OkKbnfp17x5PJLfn6u8sWEBQCkQ3gDUrXA4V3lzu51wVUqGIbW05M43pfIGoBQIbwDqkm1L0Wjh0ViGUfrnNjcT3gCUFicsAKhLpiktWHCuzj+/W1OnHlXyKdOMRYverXPPvUHTph3NtCmAkiC8AahL8bjU1XW0Zs06WpLU1FSe586de4Samo6QxAkLAEqDaVMAdSkSKXxd6tMVRnpOKiVZVnmeC6BxUHkDUJcGBqT16+9XOLxVixadqzlzFpTluanUoF544W653U066qiLlExKfn9ZHg2gQRDeANSlSER68cWf6M03H1JHx4PyessV3vr14IOXKBicrSOOuEiJBOENQHExbQqgLsXjzo5PydltWq4NC+3tuSa9LpcUDpfnuQAaB+ENQF1KJp0AJTlNesu15i0T3pLJkNxuWwMD5XkugMZBeANQl5LJylTeAgGPPJ6AbNtSOh0dtnECACaL8AagLiUSlam85Z9vmkyGlEiU57kAGgfhDUBdGlp5K2erEJ8vd74p4Q1AsRHeANSlSCQlt9srl8ujpqbmsoW3QCBzlqpXphklvAEoOlqFAKhLqZRH3/jGgGzblmSUbc1bMCh99rPPyu93fr0S3gAUG5U3AHUpc66oYRhKp6XW1vI8t61Nsu3cfxfH4+V5LoDGQXgDUJfyK17ptBOqysHvl9zukccBAMVAeANQlzZsWKUf/nCR7rvvIklSc3N5nuvzSc8885+68cYlWrPmdsIbgKIjvAGoSwMDverre1Oh0Ha5XCrbmjeXS4rH92j37pcVCm0jvAEoOsIbgLpjWVI0muvx5vFITU3le35LS65VCGveABQb4Q1A3UkkpFissMebYZTv+c3NufNNqbwBKDbCG4C6E4vlTlcoZ4PejJaW3AkLpulUAgGgWAhvAOpOKCSZZvmPxspoa8tNm5om7UIAFBfhDUDdGRiQksl+SZLf31H28Nbamqu82bZTCQSAYuGEBQB1Jx6X5s9fKssydcghp5Q9vB1yyEKdeOL/p2nT3iaPR+rvl7q6yjsGAPWL8Aag7iQS0uGHf0CHH/4BSeVrE5Ixb97hOu+8GyRJtu1M4wJAsTBtCqDuZI7Gyih35S0/LBoG06YAiovKG4C6k0hImzY9LsNwa9asE+X1Bsr6fLc7rW3bVsk0o5o//wzahQAoKipvAOpOMin9z/98TitWLFU4vKPslbemJku33HKKfvazs2XbNuENQFER3gDUnURCisf3SZL8/ikV6PPWJI/HL9tOK5WKEd4AFBXhDUDdiUZTSiQGJRny+9vLvmHB63WaA0tOs2D6vAEoJsIbgLqzb1+ux5thuCrQ5y0/vA1SeQNQVIQ3AHVnYMCZMg0Epsi2y7/bNBiUvF7nlIVkMjRs9ysATAbhDUDd6e/vkyQFAp1Kp6W2tvI+v6WlcNqUyhuAYiK8Aag7g4O5zQqVCG8+n+T35843Zc0bgGKizxuAujN79pn68pc3y7YtSVJzc3mf39Qkvf/9P5B0nYLB2VTeABQV4Q1AXUmlJNv2qr19jiTJNMt/PJZhSNOnL1Qq5bxOJJxjsgyjvOMAUJ+YNgVQV+JxKZ3OvfZ4nEpYueUHxlTKCZEAUAyENwB1JRKRXnzxx7rzzo/ozTcfVlNTZSpe69c/rF/96lNavXqF0mmx7g1A0UwqvG3fvl1XXHGFzjnnHH3mM5/Riy++WKxxAcCEhELSjh3P6rXX7tPg4NaytwnJ6O/foJdfvkNbt/5FhuGMCwCKYVLh7corr9Sxxx6re++9V5/61Kf0ne98p1jjAoAJCYWkRMJpFVKJo7EyOjq6JEmxWJ88HmlgoDLjAFB/JrxhYd26derv79fFF18sl8ulZcuWafbs2Uqn03K73cUcIwCMWf65poFAZ9k3K2TkwluvXC5nOhcAimHClbfXX39ds2bN0re+9S0tW7ZMn//85+VyuQhuACoqHpdisdwJC5WrvHVKkqLRXkmiXQiAoplw5W1gYECrV6/W17/+df3zP/+zVq5cqW9+85v65S9/Kd8Y/lN3+fLlWrFixUQfDwAjSiZzlbdKTpt2duamTSXCG4DimVSftwULFujcc8+VJH30ox/Vz372M23ZskWHHXbYQT/b09Ojnp6eYe8vXbp0MkMC0OCSycLKW6WmTbu7c9OmEuENQPFMOLzNmDFDqUwHyv1s2x5T1Q0ASiUWs3TYYe9XIjEgr7etYpW3YLBFc+eeJr+/XZaVUjxOT3QAxTHh3yannHKKrrvuOt1///1673vfq1//+teaPn26Zs+eXczxAcC4pFIufexj92RfVyq8+f2GLr74iexrKm8AimXCGxba2tp07bXX6je/+Y0+8pGP6Omnn9a3v/1tGZz/AqCChjbDrdRkwNDQSHgDUCyTquMfddRRuummm4o1FgCYtMHBiPr796i5uVteb2sFK2+SaSaVSPTJ5wsqkWiuzEAA1B2OxwJQV1555Y+6/voFuvvuv5VUmXNNJamtTbrnno/rv/5rpt544yGOxwJQNIQ3AHVlcNBpzREIdMq2KzdtGgxKPp/T6y0W62XaFEDREN4A1JVQKNfjzbKk1tbKjCMQkJqbc73eCG8AioXwBqBu2LYUDufCWypVufDm8+XCWzTaK9OULKsyYwFQXwhvAOpGMilFo5lp0ymSnLVnleB2S62tuWnTVGr4TlgAmAjCG4C6kX+uqd8/RR6PM31ZKW1tuWlTy5JiscqNBUD9ILwBqBuDg1I8ntuw4PM5FbBK6ejIHZHl8Uj9/ZUbC4D6wXktAOpGKCSdeeZVOumkv9OsWSdVtOomSXPnvk1/+7d3KxicI7fbGR8ATBbhDUDdiEalmTOXaObMJZKcRrmVNGVKl4466m+yr5k2BVAMTJsCqBtD23FUOrwN7TFHuxAAxUB4A1A34nHpkUe+oscf/zdZVlrNFT6RyuuV/vKX6/Xoo19TMhkhvAEoCsIbgLoRCsX0l798X3/609UyDFfF17w54e06Pf30NYpEdhPeABQF4Q1A3di3L9cmxDCMih2NleH1OrteJWfHKX3eABQD4Q1A3ejrc8JbJjBVw5q3QIAjsgAUF+ENQN3o7y88XaEapk3zj8hKJis7HgD1gfAGoG709xceSl8NGxb8/ty0KZU3AMVAeANQNwYGctOmqZTU2VnZ8Thr3nLTpqx5A1AMhDcAdcSrzs5Fam+fK9uW2tsrO5rWVqm1ddb+ExZ8VN4AFAUnLACoG0cf/UktXPhJSVIqVflp07Y26bjjvqRTTvmSJKdJr21LhlHZcQGobVTeANSN/MpWpQ+ll4ZX/kzTCZUAMBmENwB1IZ2WEgk7+7rSO00lp/Ln2T+/Ydu20mnONwUweYQ3AHUhkZDuvPMD+t732rVhw+8q3uNNklwuKRbbpGuvna4bbjhaLpcUClV6VABqHeENQF2IxTKNcAfl8QSqovImSR0dbYpEdisc3imPRxoYqPSIANQ6whuAutDfLyUSmVYhU6omvHV2dkgyFI/3S0orHK70iADUOsIbgLoQiUjxeK5Jb7WEt9ZWt/z+Dkm24vF+TlkAMGmENwB1IRq1FYvljseq9KH0GYFA7ogsTlkAUAyENwB1ob8/IstKyeMJyOPxV8WGBcnZcRoIOEc9RKOENwCTR3gDUBf6+nLr3ZyvlRxNTiDAEVkAiosTFgDUiaDOO+9GGYZLti21tFR6PI7mZun443u0cOF71d29mMobgEkjvAGoC4bRrhNPvFSSlExKHR0VHtB+ra3S4sV/I9f+eQ7CG4DJYtoUQF2IRnPfW1b1hLepU51jsTIIbwAmi8obgLrwxhsv65VXntGMGcdp2rS3V/xQ+oxgUBocfEs7dz6ltrZZmjLljEoPCUCNo/IGoC6sWfOIHnzwEq1Z8wv5fLkzRSvN75d27nxa9957oZ5//kYqbwAmjfAGoOZZlhQKZXabdlZNmxBJMgypo8PZbUqrEADFQHgDUPMiESkSyTXorabwJkkdHU6ft1isV6bphE0AmCjCG4Cat3dv7lzTajoaK6O7O9fnLZUSvd4ATArhDUDN6+uT4vFM5a2zasNbNNory5JisQoPCEBNI7wBqHn5h9IHAtVXeevsbJPL5ZFpRiQl1N9f6REBqGWENwA1LxKRTNNp9Ob3V9+at+ZmQ4FApwzDJdPcp1Co0iMCUMuqZDM9AExcJCJ96UsvyzRjcru98vkqPaJCzc3S5Ze/Ia+3VYbhYtoUwKQQ3gDUvMzpCk1NznxptU2bBgKS1xuUYTivaRcCYDKYNgVQ8yKR3Pe2rao5XSGjo6OwPQiVNwCTQXgDUPO2b9+iG244Wnfd9TdKparnXNOM7m7pqaf+Szfd9Ha9/PKdCocrPSIAtYxpUwA1r7d3r/bseUUul0eWJU2ZUukRFWptlaLRndqx4wX1928ivAGYFCpvAGpaMin19+d6vLnd1Tdt6vVKbW2ZRr29hDcAk0J4A1DTBgelaDR3uoLXKzU1VXhQIwgGnSOyotFeWoUAmBTCG4CatmuXlErlwlu19XjLmDLFqbzF432KRiXTrPCAANQswhuAmjY4mDvXtBoPpc/IhLdotFfptLRvX4UHBKBmEd4A1LRoVIrFqvdQ+ozubmfaNBbrldcr7dxZ4QEBqFnsNgVQ0yIRafbsU/T2t1+qQw45qWorbzNnztKxx35WU6YslGE4FUMAmAjCG4CaFotJixdfoMWLL5BUfacrZMyY0a0LLliRfc2mBQATxbQpgJo2tO1GtVbehrYvIbwBmCgqbwBqWjQqbdnyZ7lcHk2fvkSBQHWmt+Zmae/e1xUO79CsWScqFGqp9JAA1CgqbwBqWiQi3Xnnh3XLLScrGt1btZW3KVOku+/+mH760zPU27uORr0AJozwBqBmWZY0OBhXJLJLLpdHLS0zq+50hYyuLikQyJyy0Ed4AzBhRQlvGzdu1FlnnaXe3t5i3A4AxiQalfr6tkqSgsHZsiy32tsrPKhRtLRIzc25Xm/RqHO0FwCM16TDWzqd1n/+53/KpF04gDLr7ZUGB9+SJLW3z1U6XX2H0me43VJra67Xm2U54weA8Zp0eLvzzjt19NFHF2MsADAuvb1SJLJZkhPe3G6nwlWtgsHctCmNegFM1KTC2+bNm/Xwww/rkksuKdZ4AGDMIhFpYCAT3uZV7aH0GR0duWlTw6BdCICJmXB4syxL3/ve93TFFVfIX63buwDUtXA4P7zNld8vGUaFB3UAHR25aVOJ8AZgYibc5+1Xv/qV5s6dqxNPPHFCn1++fLlWrFgx0ccDgKJR6dxzf6RTT/2GAoHOqm0TkvHOd35QTU1/VWvrTEnDGwwDwFhMOLy9+OKL+tOf/qSVK1dm3/vwhz+s7373uzrttNMO+vmenh719PQMe3/p0qUTHRKABhONSk1NAXV1LZJUvUdjZUyb1qXp07uyr6m8AZiICYe3q666quD10qVLdd9996mrq2uUTwBAcUUiha+rPbwN7UFH5Q3ARNCkF0DN2r17j1asOF0rV14mqXrPNc0IBKSHH/6yfvGL98k0Y1TeAExI0c42feKJJ4p1KwAYk23bNuqtt55QIjEoqfrDW2ur9MYbv1Ff35vq79+ojo6jFI9X/7gBVBcqbwBqUjIp7d2b22kqVf+0aVeX1N6+UJLU17detk2jXgDjR3gDUJNCIWnfvlyPN6n6K1gdHVJHx6GSpH371qupiUa9AMaP8AagJu3aJYXDuaOxbHv4hoBqEwhIXV1OeOvrW0+jXgATQngDUJMGB6VQKDdtmk6rag+lz3C5pBkznPDW379BEjtOAYwf4Q1ATRp6NFY1H0qfb+bMXOVNovIGYPyKttsUAMopGpUWLjxHbW2zNGXKArlc1X0ofcacOQs1b97pmjr1KElU3gCMH+ENQE2KRqWzz7664D2vt0KDGYeOjhZ97nOPZ19TeQMwXkybAqhJI52uUM2H0meMdMqCbVdmLABqE+ENQE3asWOXtm9fpVhsn6TqbxOS0dwspVJx7dnzigYGNisel2KxSo8KQC0hvAGoSS+88IBuvvlE/fa3X5FU/TtNM/x+6U9/+g/dcMPReu65G2UY0p49lR4VgFpCeANQc2xb2rPH2WkaDDqnK3R0VHJEY9fcLHV25hr1ejzS7t0VHhSAmkJ4A1BzotHc6QodHc7pCrXQJkSSgkGpvT0X3mjUC2C8CG8Aak5fnzQ4mDtdwbKcc0NrQXe31NaW6/Vm2zbtQgCMC+ENQM3Zu1caHMydrmCa0ty5FR7UGLW1Sa2tU+X1tiqRGFAs1kd4AzAuhDcANWdwMK3Bwa2SpGBwjpqbnVBUC3w+yes1NGVKbuqUaVMA40F4A1BztmzZKcsy1dw8VU1NAXV01EaPN8kZZyCQ27TQ17eeyhuAceGEBQA1x+ebocsvf1OxWJ+k2tlpmtHWJr3nPd/RGWf8uzo7D1Mo5OygrZUACqCyCG8Aak487t5fuXKqV7UW3jo7lT3bVJKSSWcHbS2czQqg8pg2BVBzotHC17XSJiRj6M5Yw5B27arMWADUHsIbgJpi29KDD16ve+/9tDZvfkqplDRtWqVHNT5Tp0qxWEgPPvhF3XPPx+XxODtoAWAsCG8Aakpfn7Ru3e/10ku3KxLZLcuS5syp9KjGZ+FCyTACWr36Nq1de5dSqRg7TgGMGeENQE154w0pFMr1eGttdXZv1pK2NikY9Ki93Tkdor9/I+ENwJgR3gDUlJ07pYGB3OkKtbZZQXLWuHV10S4EwMQQ3gDUlLfeGlQ83i+PJ6Dm5u6aDG+SE95yjXo3UHkDMGaENwA1w7alN9/cIsmpuhmGUXM7TTMKw9t67dolpVIVHhSAmkB4A1Azdu6Udu3KTZnadu31eMuYMUPq6MiFt1TKWc8HAAdDeANQM9avl1pagjriiPM1d+5pSqWkWbMqPaqJWbBA6uw8UgsWnKVZs06WxyNt2lTpUQGoBZywAKBm7N4tzZt3qubNO1WSM804c2aFBzVBfr+0cOERuuii32Xf27atggMCUDOovAGoGTt3Fr5ua5O83sqMpRiGnrSwfbuUTldmLABqB+ENQE2wLGnPHmnv3tcUieyWbds1u94to6tLSiYj2rXrJUUie5RMShs2VHpUAKod4Q1ATdiyxTnA/Re/eJ+uvXa69u3bULM7TTO6u6WVK7+kH//4GK1bd788HmddHwAcCOENQE146y3J7U5pcHCbJEPB4Oyar7wdcogUDOYa9UrO1CkAHAjhDUBN2LlTGhzcKttOq7V1htxuX82Ht9mzpc7OhZKcdiGSs2nBsio5KgDVjvAGoCbs3i299daTkqSZM49XMinNnVvhQU1SU5M0d26u15skxeO0DAFwYIQ3AFUvlXI2K2zc6LTVWLDgbHk80tSpFR5YERx+eG7a1LZtNTWx7g3AgRHeAFS9TZukVMrWhg1OeFu48Gy1t0tud2XHVQzz50+V19uqRGJAsVifJPq9ATgwwhuAqrd5sxQOb1IotEMtLdM1bdrban6naUZ3t6Hp04+VJG3b9qwkaetW5xxXABgJJywAqHq7d0tTpizQ1762V31962UYRs1vVsiYP18688zrFAx2qrPzMElSLOa0Rqn1NX0ASoPwBqDqZU5WCAQ6dcghnZJq90D6oaZNk+bPP7ngvaYm55B6whuAkTBtCqCqJZPS3r227Lx5RNuWOjsrOKgicrkK/y6pVEKSM3UKACMhvAGoauvXS9u2/UXXXXeIfv/7f5akumgTkq+722mDcuONS/TrX39RktOsl3VvAEbCtCmAqrZtm7Rp0+8UDu9QIjEoSfL5VDcbFiTnjNOWlqnavftlDQ5uVTqdVCjk1fbtzikMAJCPyhuAqrZzZ66/28KFZ0ty1ru56ui3V3e31NW1WFOnHq14vF8bN/5BXq+0bl2lRwagGtXRrz8A9WjLlrC2bPmzDMOl+fPPkFQ/mxUyFi6UTFM68siPSpJeffVXkjjnFMDICG8AqlY4LP31r0/IskzNmnWS/H4ntdVbeOvokAIB6aijnPD22mv/I8tK0+8NwIgIbwCq1u9+56x3k3JTpqmUdNhhlRxV8RmGs+5t2rQl6uw8TNHoHm3e/KTCYenllys9OgDVhvAGoCr19kqrVw9f79bWJh11VCVHVhrd3ZJhGNmp09de+x95PNKvfy0NDFR4cECDSqeleFwaHHTOV37rLSmRqPSo2G0KoEo98ojz9QMf+Ik2bPidZs9+pyRp8eL62qyQMX++9NJL0vHHf16zZ79Thx56jiTnXxR33il94QtOhQ6AZFlOyyDTdL5Go07IisWc91Ip52s67Xyf+ZN5nflc5vtUKvd66Gdt2/lecn73fOIT0jHHVPbvT3gDUHXeekt67TXnF+WcOe/UnDlOcDNN6YQTKjy4EjnxROkPf5Dc7kXq6lpU8LMtW6Tf/lZ63/sqNDhgHEYKVrGY8yeZdIKQaRaGrExYyn//QD9Lp53nWJYTrlwuye12vhbzP3I8+1NSU5PztVrWoBLeAFSdRx8dubo2Y4Y0Z075x1MOLpfzX/N/+UvuXz62bcswDLnd0tNPSwsWSIcfXtlxonbYdmHVKZl0qlOZClV+pSmdLqxKWVbu69BKVH41K/9z+VUs2y5tsHK7nT+NivAGoKqsXStt2iS53bbuvfdCTZ9+rN7xjv8lt9uvxYsrPbrSOv106dlnpWh0n1au/JJ2716rSy/9qwzDkGFI990nXXaZ1NJS6ZFiojJTcPnhJxOqEonc1/wglT+ll6k+5VetRvqaqVRlAlQmTGVCVCZIlWIq3kOyKDn+EQOoGrYtPfaY88u/t/cNvfzyHdqw4VG9+91flWVJ73hHpUdYWs3N0pFHSmvXtmvTpj8qHN6hHTte0KxZb5fkTD/dead08cWsfyumTKDKhKLR1lANrS4Nff9Af/IrV5lAlfkq5QJVMUNVo1en6hnhDUDVePZZafdu5184f/3rzyRJCxactb9Br7PTtN69613SmjUuHXnkR/Tcc/+/HnzwC7root8pEOiUYUgbN0p33SWdemp9Hp2VCVJDp/Gi0Vxlaujap/xK09BpvaEhK/+6/J9l1k9ZljOOoRWqYsncE9UjnU4qlYrv/5NQKhVXOu18bW6eqo6OeZKkSGS3+vs3SzqxsgMW4Q1AlUinpT/+0fmX5rp1D+rJJ6+SZOiEEy6RZUlHH13pEZbH7NnSvHnS0qXf1Pr1v9XOnS/q5z8/R5/5zKMKBKbI7ZZefVVas8ZZA3j44dLJJ5f2rNfM1J5pOuEpU5VKJEYPUEOrTfmLzkdasJ6/fiozzZepTkmlqUxlEKgqy7JSMs3Y/vAUywtScXk8fk2f7mztTKdNvfTS7fuvda7L//6YYz6j2bNPkSS9+up9ev75G/OCWCL71TBcuuKKN7PPv/HGJertfX3EsZ1yypf1/vd/X5K0d+9rWr36p5JuLe0/kDEgvAGouHjcWc8VDku9va/o3nsvlGTrzDO/q4ULz5bLJb397ZUeZfmccoq0efMMffazf9CKFadrx45V+sUv3qfPfOa32VMmmpqcXnhPPy396U/SrFlSa6sz5ezxOD/PrD3KBKL8XXpDQ9fQdVhD2yuUoyqVuTfKyzSjikZ7s9WmwgpUUosWLcteu2bN7QqHdxSEoUxAmjdvqZYs+ZQkac+eV/Tww/9rxGpWKpXQ5z73uLq7nUWs99/fozVrfj7i2GbPfqc+//mnJUm2ben++y8e9e8xc+bbs+EtFNqmDRseHfE6wyj8H5nX2yqvt1Uej19ut08ej18ej09ut0/BYK683do6QzNmHH+wf5xlQXgDUDHptPTEE04ASSalRKJPd9xxvpLJkI4++uM69dT/LcmpLmW26jeCJUuc0yXc7kOyAW779uf03HM36rTT/nfBtYbhBJ5du5w/xWYYLEAvFtu2skHG5fLI53PWAZhmVHv3vlZQccoPUUcccb6am7skSa+9dr+2bv1LQRBKpxNKp5Nqb5+nc865RpJkWWnddtupo4anZct+qOOPd4LQ6tUrtHLl3484ZsNw6V//NSVjf0r/85+v1c6dq0f9O2bCWzIZ1oYNvxv1OtOMZr9vamrZ/yewPzjl/nR3H5m9zu326thjL5Lb7d9/bSD7CGI8YAAAFoRJREFU1ePxZYObJB1++AfV2blIHo8vL5T5suEss5Nbkr74xVWjjjNfV9fh6uysju3ek/q/5EMPPaTbbrtNAwMDOvzww/XVr35Vc+fOLdbYANQp25aef96ZJh0czJ8Kc6ur63D5fEF96EPLZRiGUqnKN8QsN8OQjjvO+efT3j5Hn/3sH7Rq1U069dSvV3poNSU/LBWGl7hcriZNm3b0/utsvfrqvUqnk/v/JPavg3KC0aGHvi+7aWTTpj9q9erb8q4t/MzFF/8pGwpuv/1cbd/+fPaZlmVmx3bCCV/QBz94kyRpz55XdfPNo6+j+sIXnsuGtzfffEirVv1kxOtmzDgu+71huLR16zOSRm5MlkrFst/7fO1qa5u1PwTlh6dM0LGy1aolSy7U/Pln5gWh3Ndp096WvWdX1xH69KcfGbGa5fH41dzcnb32Ax+4UR/4wI2j/v1zfydDF1zw04NeJ0kdHfOya9Xq0YTD25YtW/T9739fV155pY4//njdfffd+vd//3fdemvl54IBVA/blvr7nca7vb3O99u35zYm5K818vvb9clPPqh4fJ+ampolSe3tzg7MRnPqqdJTTznVyY6OeTrrrO9mf/bSS/9X69bdr8MOW6ZDD32f2tpmVnCko0ulEorH+4etOcqsUXI2ozhB59VX79Pg4NaCNUqZ0HPIISfp2GMvkiTt27dBv/71pcMqSZkpvk9+8sFs0Fq58jI9//zIoWD69GN06aV/leSEgnvu+bhsOz3itT5fMHvPffs26K9/HT1AWJYpt9srSYrH+xWN7in4eS7weAvuP2PGcXnhKVctamoKZKfKJemIIz6k9va52ftkrnO7vQWByDAM9fQ8VRCYMt9nqlUZxxxzoY455sJR/0753vWufxrTdX5/e/aUEBTfhMPbCy+8oOOPP17vfve7JUkXXnihli9frsHBQQWDwaINEEDp5TfizG/mGYs5X/ObeR6sHUIyWfgns7i9qakwqGXWNsVi+/T66w/qbW/7hNxur1wud8G/hI44ojHbYni9zvTp6hFmqF555R699tp9Wrv2LklOxWX27HfK52vXzJnH6+ijPyZJSiQG9cYbK2XbliwrLdu2ZNvpbLXoyCM/ml3Ts27dA9q48Q/DKknpdEKtrTN13nk3ZJ9/yy2nKJkMF1SeMuHs7LP/QyeffJkk6eWX7zjgGqV/+ZeYPB6/JGc6bsuWp0e8LpH4dDa8/b/27jy6qTpvA/iTrUm60kaESsu0lNKy2jLKTulAS0GlizguU5z3KDi+LFUWFWR8EdRx6oLLGc958QjoQV8QUDujIEpGkUXfFzwUrLiAhdJWBARPIdDaNE3y/nF7b29ubkpBaHLJ8znnniS3t+0NP5I8/S3f29raHHAuEwC4XI3SfaPRKgsuFp8AY7P5Dn8NGHBb2/eYoddHwGCIkMJRYmL7ZT169x6DwsLVMBjajxHvGwxmn/lUd975z7af2R6wdCr/mW22dNx//76Az0kuPX2yzxy0johXJqGrzyWHt7Fjx2LYsGHS40OHDsFisSCK1SOJLklH1dCdTmETV+opK51f6FY+KV2tnIL8MjPy1X2Xqyq62ax8rl78+OP/Yu/eV/HNNxvQ2tqMysrXcPfddunDXDwH2dtM2Bk7Fjh4UAjRcvn5zyE1dTyqq7egpmYbTpzYL81DGjjwDim8nTv3E959966AP79HjyFSeDt6dDt2735J9biEhL4+j0+d+hYtLedVj3W55MNxsYiKutZveE0IUVZ4PK3SsZmZt6Jnz2y/+Ukmk9Vn3lNc3O9QWvqRTw+V8D1CeIqK6i4dW1CwHAUFywM+f7nbbnu7U8fZbP6XLwskKuraTh1HdLEuObwlJCRI93fu3InnnnsOM2fOhIFLhUjDxB4oedkDp1P48BQDlDIMXWhT1pPqqJyCWoC60tXQxQnvXfHSPXbsSxw5YseBA+vw888HpP19+uRh2LAyv+B2xx1COYxwZbMBM2cCa9YAp061t39CQhqGDZuDYcPmoLW1GbW1O3D69PdwOs9JK/gAYS7TwIG3Q6fTQ6czQK83QKfTSz1L0dHtw60ZGVMQG9tLCkHyniX5sB0A3Hvv59DpDLIepwgYDCYpbIn6978V/fvf2qnnOmrUgk4dFxERhb59eZFXCm86r/fSL7PqcDhQXl6O6upqzJ8/HyMuovz56tWr8cYbb6h+bceOHZd6SqRx8t4nsQ6UvPdJ3Cded09ZzkDt+nvKENVR8FLrgRLDzZUMUMHk8biluUUWSxz0euFvuoaGIzh//oTPHCVxi4y8BhkZhQCEnpZPPnkULlcTXK5GtLScb9sa4XSexeTJr6BPnwkAgO3bn8Bnnz0OQOiVyMq6B0OHzvDr2RGDW0ZGF/5DhLDWVuDtt4VeOP59TNROXhNQXhtQvK/TtU/XEN+/xU0cXZDfV+4zGv2/npcn1GIMpksOb83Nzbjvvvtw/fXXY86cObBYLBf+pk7IyclheOtC4n905fCafJiuudk38Ii9U2KAUvYyqd3vqFcqUIDqyt6nriRfAedbJ6kFsbG9YLUKvdpnzhzFyZNVPqUAxPsej9unp2LbtsfhcNT5rYBrbXUiM7MYw4c/AEDo+dqwYapPEJOvgJs161t07y4MUb333jR8/fX/qD6H5ORRuPfezwEI4e3ppyMDPt+bb/5v3HDDfwIAamq24eDBfyE5eTQyM4t8Jm2LdDohuF3t1zG9WF4vsHWrsIiBBWXpYomf9MpLc3k8/qEG8L2V/wEr9tLL35fF92Z56BH3iQFIfowyJIn75N8rbmq/C/DfbzAI80TF+oYREcKt/Hcpf56WXfKwqd1uh81mw0MPdW7lCalTznNyudoniYsTxdXmOakV1FQLRWLIUnss3oovYDEsyf9akd92RWi6HC8qr9frUxLA7XYhLi5Z+vrPPx9oWwEnhqH2ydkJCX2RlCT0IDscx7B//+s+q+SE44RjJ0z4O+LihNI4O3f+HT/8sEkRnlxwu1vQo8cQ/OlPmwAIl2F56imz/0m3KSxcLdVeOnRoE7ZsKQvw72T0CW8HD/4TJ09WqR4rH0bzej1wOOoVR+hgMgkTu+VzkGy2fujVa7hf7SVhsnd7l5jRaMHEicthMkVKxS7btxif5fqpqX9AauofAj5/vR744x8Z3NTodEBBAdC9O7Bpk/DewPprl4+y10Z+3dFA4UYtQATar3aMvAdIfqv8nWrnoJyTqvwd4n0xOBkMvqHGbBaCTkSE//nKf8fVEHSuRpf80q+urkZlZSVycnJ89q9fvx6JiaG5bF1OfHEqh9zE0CS/9EugauSdne+k7KFS/ixlN29X9jTJ/4q5HDwetxScxMuduN0t6NYtRZoLc+rUdzhzpka6rIn88iZWqw1Dh04HIASNDz74i1+ZAfF29OiFyMwsAgDs2/c6Pv54ntSb5UuHJUvc0iqvioq7AxaZHDr0L1J4O3/+OLZt+6+Az3XkyAVSeGtoOBJwpZx8ArVebwKga5tLZIGyVpJYtBMA4uP7ID39ZkU5gPZj5UUmc3KWwOk86zNXSZyHJJ4jIKxKfPDBozAaLVIo0+tNqivgxo1bgnHjlgR8/iKdToeRI+df8LiOiB+St90WnmVBLsbQocD11wPV1cCRI0LZlZ9+Et6zxCGey0FtOEq+X+yxAfzfq5ThQnyfkYcUZdBQCwzKXh+14S15QFH2+Ch7gdQCiTxQiaEmIqJ9k4cf5XPS+igAadclh7d58+Zh3rx5l/NcutQ77wBffeXb6xTojeNKupzzV5qbz/oVohRrIMXH95FWlf3yyyHU1u6UhaEWn2CUn/+c9GH+738/itOnv/W7lpzb7UT//lORl1cOADhxYj9ee22YzxCc3IwZe9Cr140AgN27X+6wyKQY3gAd9u1bjUBFJh2OH2WPvHA6z0qPlKvbvF43dDrhv3uPHkNgMkX6TMwWl/v36tW+tDEmphfGjFmsWozSYDAjLq69R2n06IeRlfUfquFJrFcGCEFHHiQ7kp5+E9LTb7rgcQAwYMDUTh1nNJpDpnBlaytgtQIpKUDfvkIguUyzL656BoMwH1CcE+jxAHV1Qh09sUSL0yncuttKl11oGEwZhEym9iEoZe+NfGOwIep6Ydvp/lsu+SIMywmBJyIiRvogbmiogdPpgLLqdmurEzEx10nhpbHxFL76ao10jMfjku67XL9i9OhHkJCQBgDYvfsf+O67dxSFKIXbhIS+uOeendJ5PfusLWCRycmT/yHVXqqr+xwffDAj4PMbP/5vMBqFob3a2s/w44//p3rc+fPHpft6vbEtuOlUK3SLk+AB4NprByEtrUB2WZP2S5zIhzd1Oh0KC1dCrzeqlBqwIj6+j3TsoEF3ITOzpG34LwI6XeB+/s5W6I6JSfQpjNoRm62fX92oQDoT3K4mYs+z2SwU3I2LE1ZR9usnhA8Oyfx2er0QglNSgn0mRNQVwja8rVw5GydONMLjaW0LTy54PC60tjrx+9/fL/VkHDz4Pj7+eL7f8J5o8eImaTiwouJu1Nd/rvr7hgyZhpIS4cK7TU2nYLcHnis4ZEipFN4aGo6gtlZ9AYfZ7FsMOSqqOzwet2qFbnm9oWuuyUBW1j0+gUhekFIeLsaP/xucTodU+Vt+HTmLJV46rnv3AfjrX5sDFqGUE0scdEZ29r2dOs5ksvqUKAhnHa2+Avx7YNQmCF9oWEptErDB0D6HRpw4bDYDUVFA795AbCyDGhHR5RC24W379jVoblYvMtmnT7503+VqQkPDYb9j9HoTjEYz3O4WKTTYbOloaTnnM2wmHGdBz57Z0vdGRnbHiBHz/I4zGCJgMlkRH58mHXvjjbOQmVmsenFd+XAcACxYcBydkZw8CsnJozp1bGrq+E4dp9Pppd66q5k8DClXxMrn2ShDj3K5eaAJzMqvi98nXp1AvipLGajkXxeHtOQrr5Srr+RDXmHWGUhEpGm/qc7bldBVpULKytagvt4Dvd4Ivd7UNjQnFJm85poMdOuWAkC4vExj489SRXBxVZ5ez2JLl4vaCq9Aq7wCLStXmyCtnIisnNsjf6zWkyRf5i5+vxiEzGZhfpbV6jsvSAxEREREV0rY9ryNG/dnHDhw4ePM5li/4UmtUfYUeTzqq8PE+/JwJO+lkffqqE10DhSC5Jvaz1EuW5dPlFb7GewlIiKicBa24e1KUQtKgWoEKYfK5CFJHnJMpvZ98vvK3qRA4cloFHqIxN4ii6Xj3iuGIyIiotAVtuEtKgqIj/efe9TRFigsyb8mTtKWb/J5RsqNQYmIiIguRtiGt8mThY2IiIhISzi1moiIiEhDGN6IiIiINIThjYiIiEhDGN6IiIiINIThjYiIiEhDGN6IiIiINIThjYiIiEhDGN6IiIiINIThjYiIiEhDGN6IiIiINIThjYiIiEhDGN6IiIiINIThjYiIiEhDGN6IiIiINIThjYiIiEhDGN6IiIiINIThjYiIiEhDGN6IiIiINIThjYiIiEhDGN6IiIiINIThjYiIiEhDjME+ATU5OTnBPgUiIiKioNqxY4fqfp3X6/V28bmEhJycnID/KBQ62E6hj22kDWyn0Mc20oZQaCcOmxIRERFpCMMbERERkYYwvBERERFpCMMbERERkYYYli5dujTYJxEMXq8X2dnZwT4NugC2U+hjG2kD2yn0sY20IRTaKWxXmxIRERFpEYdNiYiIiDSE4Y2IiIhIQxjeiIiIiDSE4Y2IiIhIQxjeiIiIiDSE4Y2IiIhIQxjeiIiIiDQk7MLbnj17cNddd2HixIlYsmQJGhsbg31KJDNnzhx88803AIRCiK+88gomT56MkpISVFRUBPnsaMuWLbj99ttRUFCAsrIy1NXVsZ1CzHvvvYfi4mJMmjQJCxcuxOnTp9lGIaympgYTJkzAL7/8AgBYu3YtpkyZgltuuQWrVq0K8tmFt1mzZiEnJ0fannzyyZB5LRmD8luD5Ny5c1i2bBkefvhhZGVlYfny5Vi5ciUefPDBYJ9a2NuzZw+2b9+Oqqoqad/mzZvx5Zdf4vXXX8e5c+cwf/58DBgwABkZGUE80/BVX1+PF198EY8//jiys7OxceNGLFu2DCUlJWynEFFTU4NVq1ahvLwcKSkpeP7557FixQpkZWWxjUKQ2+3Gs88+C5fLBQDYu3cvNm7ciJdeeglWqxULFixAeno6cnJygnym4enkyZOw2+0wm83Svk2bNoXEaymset527dqFfv36ITc3F926dcO0adPw6aefBvu0CEBVVRU8Hg+sVqu0z263o7S0FD179kR6ejomTJjA9gqiyspKZGdnY/To0YiMjERpaSmqq6vZTiFk//79GDVqFAYPHoyYmBgUFBSgpqaGbRSi1q9fj4EDB0qPt27diuLiYqSlpeG6665DUVERPvnkkyCeYfhyOp0wGAw+wQ0Inc+lsApvhw8fRnp6uvQ4JSUFDQ0NaGpqCuJZEQDMmDEDCxcuRFxcnLRP2V6pqak4duxYME6PAIwdOxZz586VHh86dAgWiwU//PAD2ylETJkyBY888gi8Xi/OnDmDrVu3YtCgQXwthaC6ujp89NFHmDFjhrTv8OHD6Nevn/SY7RQ8J06cgNvtxvTp0zFp0iQsWrQIp0+fDpnXUliFt6amJkRFRUmPzWYzDAYDw1uIamxs9Gkvq9XKtgqihIQEJCYmAgB27tyJRYsWYebMmfj111/ZTiHCaDTCZDLh/fffR2FhIXbt2oXx48fztRRiPB4PysvL8cADD8BisUj7le0UGRnJdgoSh8OBpKQkLF68GBs2bEBMTAyeeuqpkHkthdWct+joaLS0tEiPXS4X3G43YmJignhWFEh0dDScTqf02Ol0sq2CzOFwoLy8HNXV1Vi8eDFGjBiBVatWsZ1CTFFREfLz87F582Y89thjsFqtbKMQ8u6776J379644YYbfPbHxMT4tFNzczPbKUgGDx6Ml19+WXo8e/ZsFBUVITY2NiReS2HV89a7d28cPnxYelxbW4ukpCS/MW0KDcr2Onr0KPr27RvEMwpvzc3NmD17NhISErBmzRqMGDECANsplLz11luw2+0AhF6b4uJiOBwOJCQksI1CyL59+/Dhhx9KqxgBoKSkBLW1tWynEPHFF19gz5490uPW1laYTCYkJSWFRBuFVXgbM2YMqqqqsHfvXpw5cwavvvoq8vLygn1aFEBeXh7WrVuHhoYGfP3117Db7cjNzQ32aYUtu90Om82Ghx56yGeoh+0UOmJjY/Hmm2+ivr4eTU1NWLt2LdLS0jB16lS2UQh5+umnsWPHDmkDgIqKCixduhQVFRU4fvw4jh49io0bN/IzKkgaGxvxwgsvoKamBg6HAytWrEBubi4mTpwYEq+lsBo27datGx599FE888wzOHv2LHJzczFt2rRgnxYFUFhYiJqaGtx5552Ijo5GWVkZkpOTg31aYau6uhqVlZV+ZQvWrVvHdgoRN910E2pra1FWVgan04lBgwbhiSeeQGJiIttIA0aOHInvv/8e06dPh9FoRGlpKbKzs4N9WmEpPz8f9fX1mDt3LjweD4YPH465c+fCarWGxGtJ5/V6vV3+W4mIiIjokoTVsCkRERGR1jG8EREREWkIwxsRERGRhjC8EREREWkIwxsRERGRhjC8EREREWkIwxsRERGRhjC8EREREWkIwxsRERGRhvw/1FlaGjw4EwIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 748.8x655.2 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "go(fitfun, ar, model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian Analysis Results\n",
    "Analysis Results work exactly the same under Bayesian analysis. \n",
    "\n",
    "Let's run the analysis first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:59.042928Z",
     "start_time": "2020-04-13T07:50:59.013647Z"
    }
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "sample() takes from 1 to 2 positional arguments but 4 were given",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-28-880e2c1bc7df>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0mbs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mBayesianAnalysis\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m \u001b[0msamples\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m20\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m: sample() takes from 1 to 2 positional arguments but 4 were given"
     ]
    }
   ],
   "source": [
    "\n",
    "for parameter in ar.optimized_model:\n",
    "    \n",
    "    model[parameter.path].value = parameter.value\n",
    "\n",
    "model.fake.spectrum.main.composite.a_1.set_uninformative_prior(Uniform_prior)\n",
    "model.fake.spectrum.main.composite.b_1.set_uninformative_prior(Uniform_prior)\n",
    "model.fake.spectrum.main.composite.F_2.set_uninformative_prior(Log_uniform_prior)\n",
    "model.fake.spectrum.main.composite.mu_2.set_uninformative_prior(Uniform_prior)\n",
    "model.fake.spectrum.main.composite.sigma_2.set_uninformative_prior(Log_uniform_prior)\n",
    "\n",
    "bs = BayesianAnalysis(model, data)\n",
    "\n",
    "samples = bs.sample(20, 100, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we grab the results from the BayesianAnalysis object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:59.044981Z",
     "start_time": "2020-04-13T07:50:38.210Z"
    }
   },
   "outputs": [],
   "source": [
    "ar2 = bs.results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write and read the results to/from a file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:59.046278Z",
     "start_time": "2020-04-13T07:50:38.211Z"
    }
   },
   "outputs": [],
   "source": [
    "ar2.write_to(\"test_bayes.fits\", overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:59.047586Z",
     "start_time": "2020-04-13T07:50:38.212Z"
    }
   },
   "outputs": [],
   "source": [
    "ar2_reloaded = load_analysis_results(\"test_bayes.fits\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The AR holds the posterior samples from the analysis. We can see the saved and live reults are the same:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:59.049379Z",
     "start_time": "2020-04-13T07:50:38.214Z"
    }
   },
   "outputs": [],
   "source": [
    "np.allclose(ar2_reloaded.samples, ar2.samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOTE:** *MLE AR's store samples as well. These are the samples from the covariance matrix*\n",
    "\n",
    "We can examine the marginal distributions of the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:59.050824Z",
     "start_time": "2020-04-13T07:50:38.215Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#ar2.corner_plot();\n",
    "\n",
    "# with chain consumer (pretty!)\n",
    "ar2.corner_plot_cc();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can return pandas DataFrames with equal tail or HPD results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:59.052167Z",
     "start_time": "2020-04-13T07:50:38.217Z"
    }
   },
   "outputs": [],
   "source": [
    "ar2.get_data_frame(\"equal tail\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:59.053654Z",
     "start_time": "2020-04-13T07:50:38.218Z"
    }
   },
   "outputs": [],
   "source": [
    "ar2.get_data_frame(\"hpd\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Error propagation operates the same way. Internally, the process is the same as the MLE results, however, the samples are those of the posterior rather than the (assumed) covariance matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:59.055079Z",
     "start_time": "2020-04-13T07:50:38.219Z"
    }
   },
   "outputs": [],
   "source": [
    "p1 = ar2.get_variates(\"fake.spectrum.main.composite.a_1\")\n",
    "p2 = ar2.get_variates(\"fake.spectrum.main.composite.b_1\")\n",
    "\n",
    "print(p1)\n",
    "print(p2)\n",
    "\n",
    "res = p1 + p2\n",
    "\n",
    "print(res)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To demonstrate how the two objects (MLE and Bayes) are the same, we see that our plotting function written for the MLE result works on our Bayesian results seamlessly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-13T07:50:59.056230Z",
     "start_time": "2020-04-13T07:50:38.221Z"
    }
   },
   "outputs": [],
   "source": [
    "go(fitfun, ar2, model)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "3ml",
   "language": "python",
   "name": "3ml"
  },
  "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.7.7"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "4px",
    "width": "254px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {
     "0119c6b30bbe48aea72993e9e48bcba1": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "0437363051864c95a9a3cc831a860c9c": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "0aadddfba8454e5fb0c8891f433462ea": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "VBoxModel",
      "state": {
       "children": [
        "IPY_MODEL_880c9b8b89e04e03bfcaa896655a717f",
        "IPY_MODEL_a5979c69237b45af97267840721cba6c",
        "IPY_MODEL_0f64b5edd6184d3a8100b0b054fefa7f"
       ],
       "layout": "IPY_MODEL_a7ea2993152b4ef1ac80cae5576957ef"
      }
     },
     "0cab54e234f14fdfadaddb8d295aff0c": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "0f64b5edd6184d3a8100b0b054fefa7f": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "FloatProgressModel",
      "state": {
       "layout": "IPY_MODEL_837d1311f21b44c08534012c5eb98d6b",
       "style": "IPY_MODEL_583928f28b6e4e8f89976755ade02680",
       "value": 100
      }
     },
     "110c7bb8733640909d12c14bf6b2b625": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "VBoxModel",
      "state": {
       "children": [
        "IPY_MODEL_e283d5eda344431ab8095e620660d3c5",
        "IPY_MODEL_1aeee3fb9f8b4d2dad67f7ce750f0262",
        "IPY_MODEL_ee888ec72312445c9a8967f617ac911b"
       ],
       "layout": "IPY_MODEL_db1e7d75bfe04240a68b9cd96d01ce9c"
      }
     },
     "155ca44ebef3459bbf3d74dd08db26fd": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "FloatProgressModel",
      "state": {
       "layout": "IPY_MODEL_c7f7739d0d374b88ab67995eb6b29d69",
       "style": "IPY_MODEL_2e83e2ff02ef4103814f04df26aa6da4",
       "value": 100
      }
     },
     "1aeee3fb9f8b4d2dad67f7ce750f0262": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "HTMLModel",
      "state": {
       "layout": "IPY_MODEL_0437363051864c95a9a3cc831a860c9c",
       "style": "IPY_MODEL_37d26d39fa994a448805b70e3f9d317f",
       "value": "100 / 100  in 5.2 s (0:00:00 remaining)"
      }
     },
     "244668962fac45c79cfba5c9d2a3d35b": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "2b8804f15b4f4e84a41a8f193cc8aeea": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "ProgressStyleModel",
      "state": {
       "description_width": ""
      }
     },
     "2e83e2ff02ef4103814f04df26aa6da4": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "ProgressStyleModel",
      "state": {
       "description_width": ""
      }
     },
     "3231ce3ccea84054879a7f53df0a8ec7": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "DescriptionStyleModel",
      "state": {
       "description_width": ""
      }
     },
     "327300af818f4dddb8e073c0dc9e01ea": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "336c71001b41470b8d40c176f6fabc9f": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "36acdc32b54c4412aeb1e9579f0cbe59": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "HTMLModel",
      "state": {
       "layout": "IPY_MODEL_8d9d9e130a0c453f92523b8f742e446c",
       "style": "IPY_MODEL_3231ce3ccea84054879a7f53df0a8ec7",
       "value": "100 / 100  in 5.0 s (0:00:00 remaining)"
      }
     },
     "37c1e337f0a84225bd8e9cc748c045ea": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "VBoxModel",
      "state": {
       "children": [
        "IPY_MODEL_a5af1349ab8340d480aac501ea988396",
        "IPY_MODEL_36acdc32b54c4412aeb1e9579f0cbe59",
        "IPY_MODEL_155ca44ebef3459bbf3d74dd08db26fd"
       ],
       "layout": "IPY_MODEL_0cab54e234f14fdfadaddb8d295aff0c"
      }
     },
     "37d26d39fa994a448805b70e3f9d317f": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "DescriptionStyleModel",
      "state": {
       "description_width": ""
      }
     },
     "4588e71b5c564599a9ac23418fdd7be2": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "FloatProgressModel",
      "state": {
       "layout": "IPY_MODEL_d1415612239744a89313e78709a16d0a",
       "style": "IPY_MODEL_2b8804f15b4f4e84a41a8f193cc8aeea",
       "value": 100
      }
     },
     "4b847daef0ce4d6eb895fc66c72b3a28": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "DescriptionStyleModel",
      "state": {
       "description_width": ""
      }
     },
     "4e0b4aa6462949d3996b74e4aed1bf6c": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "VBoxModel",
      "state": {
       "children": [
        "IPY_MODEL_fc5a23c641ab4fe285c8658e0ab1acf8",
        "IPY_MODEL_ca6d275337514992bdf3bd0440efc422",
        "IPY_MODEL_4588e71b5c564599a9ac23418fdd7be2"
       ],
       "layout": "IPY_MODEL_0119c6b30bbe48aea72993e9e48bcba1"
      }
     },
     "583928f28b6e4e8f89976755ade02680": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "ProgressStyleModel",
      "state": {
       "description_width": ""
      }
     },
     "66aaddcf75054201b443650b70690dd9": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "795a03f7ce6e4fee9f4a100e2ea4cd2a": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "8303bb4d4b72469eb2f5f692091bdbc3": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "DescriptionStyleModel",
      "state": {
       "description_width": ""
      }
     },
     "837d1311f21b44c08534012c5eb98d6b": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "83b686826ee14f258eeb96f1fb3fe8ac": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "8775fb0f32a94948bba5b72b1be4040f": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "DescriptionStyleModel",
      "state": {
       "description_width": ""
      }
     },
     "880c9b8b89e04e03bfcaa896655a717f": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "HTMLModel",
      "state": {
       "layout": "IPY_MODEL_244668962fac45c79cfba5c9d2a3d35b",
       "style": "IPY_MODEL_8303bb4d4b72469eb2f5f692091bdbc3",
       "value": "Propagating errors : "
      }
     },
     "8b3b514603a64d6e89f9f7fdeb2f887c": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "DescriptionStyleModel",
      "state": {
       "description_width": ""
      }
     },
     "8d9d9e130a0c453f92523b8f742e446c": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "a4e5b8c4bed242bd85de5a7f3f764b77": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "a5979c69237b45af97267840721cba6c": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "HTMLModel",
      "state": {
       "layout": "IPY_MODEL_327300af818f4dddb8e073c0dc9e01ea",
       "style": "IPY_MODEL_b46a911954c049c782df69088fd23647",
       "value": "1 / 1  in 12.1 s (0:00:00 remaining)"
      }
     },
     "a5af1349ab8340d480aac501ea988396": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "HTMLModel",
      "state": {
       "layout": "IPY_MODEL_66aaddcf75054201b443650b70690dd9",
       "style": "IPY_MODEL_c532fa77c017468fb4b8491a37c5d645",
       "value": "Propagating errors : "
      }
     },
     "a7ea2993152b4ef1ac80cae5576957ef": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "b46a911954c049c782df69088fd23647": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "DescriptionStyleModel",
      "state": {
       "description_width": ""
      }
     },
     "c0a9438511ab49a4bab63cb317ca6903": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "ProgressStyleModel",
      "state": {
       "description_width": ""
      }
     },
     "c532fa77c017468fb4b8491a37c5d645": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "DescriptionStyleModel",
      "state": {
       "description_width": ""
      }
     },
     "c7f7739d0d374b88ab67995eb6b29d69": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "ca6d275337514992bdf3bd0440efc422": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "HTMLModel",
      "state": {
       "layout": "IPY_MODEL_83b686826ee14f258eeb96f1fb3fe8ac",
       "style": "IPY_MODEL_8b3b514603a64d6e89f9f7fdeb2f887c",
       "value": "1 / 1  in 14.0 s (0:00:00 remaining)"
      }
     },
     "d1415612239744a89313e78709a16d0a": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "db1e7d75bfe04240a68b9cd96d01ce9c": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "e283d5eda344431ab8095e620660d3c5": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "HTMLModel",
      "state": {
       "layout": "IPY_MODEL_a4e5b8c4bed242bd85de5a7f3f764b77",
       "style": "IPY_MODEL_4b847daef0ce4d6eb895fc66c72b3a28",
       "value": "Propagating errors : "
      }
     },
     "ee888ec72312445c9a8967f617ac911b": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "FloatProgressModel",
      "state": {
       "layout": "IPY_MODEL_336c71001b41470b8d40c176f6fabc9f",
       "style": "IPY_MODEL_c0a9438511ab49a4bab63cb317ca6903",
       "value": 100
      }
     },
     "fc5a23c641ab4fe285c8658e0ab1acf8": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "HTMLModel",
      "state": {
       "layout": "IPY_MODEL_795a03f7ce6e4fee9f4a100e2ea4cd2a",
       "style": "IPY_MODEL_8775fb0f32a94948bba5b72b1be4040f",
       "value": "Propagating errors : "
      }
     }
    },
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}