{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model selection using refnx and dynesty\n",
    "\n",
    "refnx + dynesty can be used to obtain the Bayesian evidence, which allows you to perform model selection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import dynesty\n",
    "\n",
    "from refnx.analysis import Objective, Model, Parameter, Parameters\n",
    "from refnx.dataset import Data1D\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss(x, p):\n",
    "    A, loc, sd = p\n",
    "    y = A * np.exp(-((x - loc) / sd)**2)\n",
    "\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll synthesise some experimental data from two Gaussians with a linear background. We'll also add on some noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(3, 7, 250)\n",
    "rng = np.random.default_rng(0)\n",
    "\n",
    "y = 4 + 10 * x + gauss(x, [200, 5, 0.5]) + gauss(x, [60, 5.8, 0.2])\n",
    "dy = np.sqrt(y)\n",
    "y += dy * rng.normal(size=np.size(y))\n",
    "\n",
    "data = Data1D((x, y, dy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8WklEQVR4nO29e5Rkd3Xf+93nUa/uru6e7p6enpdGEkIwQkiALjDA9cUmMQZjhFdAIxIbhUWiwRd87bWSlQtxEnSTy1pOYjuOb1bEyMAFO34MfhBkG4Mx8Q0mGmQkoUiakawH0jx6eqanX/V+nMe+f/zO79Sp6uru6ke992etXlN9zqmqX5+p+p59vr/925uYGYIgCMJgYXR7AIIgCMLeI+IuCIIwgIi4C4IgDCAi7oIgCAOIiLsgCMIAYnV7AAAwPT3Nx44d6/YwBEEQ+orHH398iZlnmu3rCXE/duwYHnvssW4PQxAEoa8gogsb7RNbRhAEYQARcRcEQRhARNwFQRAGEBF3QRCEAUTEXRAEYQARcRcEQRhARNwFQRAGEBF3QRCEAUTEXRh6Tp4+i5Onz3Z7GIKwp4i4C0KAiLwwSIi4C4IgDCAi7oIA4PxCFucXsjt+vkT9Qq8h4i4MHSLEwjAg4i4IgjCAiLgLQ83J02dbtmMk4hf6CRF3QRCEAUTEXRAEYQARcReEBsR+EQYBEXdhKDm/kN1SwEXkhX6mJ3qoCkI3qTgeDINgmwaYGUTU8nM3En+9/cypE3syRkHYLiLuwlDj+4yqxyCPwQw8cXENdx6ZgGm0LvA62+b4XFoifaFnEFtGGGoKVRcAwAAqrg/XZ5Qcr7uDEoQ9QMRdGGoKFSXk0Ti9vANxLzseLiwX92hUgrB7RNyFoeNqpoxc2QUzI19xQQCSMRNJW30ddiLurs/IlZ09Hqkg7Bzx3IWh4+KKirArro98xYVpUOCxE+IMlBx/W6/HrPx61+c2jFYQdoZE7sLQkYqZAICVQhUV14cRmTxN2GbTyN31fLx0PY+n5zPrJk050HSvQdxbSbcUhHYhkbswNGihjVkGClUPVzJlAIAVEfekbSBXdsKUSC3Q//6Dd2ApX0XMWh8P+YG6uz5vO5VSENqFRO7C0BGNtJO2WZf2mLBN+AxUvXprJhv46czrrZdowN4YvQtCt9hS3InoCBH9FRGdJ6JzRPQLwfYHiGieiJ4Mft4bec6niehFIvpbInp3O/8AQdguXkSgp0ZidfsStrJsyg2+e7akxL2ZdvuR1xPfXegVWrFlXAD/hJmfIKIxAI8T0beCff+BmX8lejARHQdwL4DbABwE8JdE9GpmluRhoSfwfQYBMA3C1GgMmUiWS8xU8Y7j+VguVOEEEfxmkXt00wuLeRyfS+Natlwn+oLQabYUd2ZeALAQPM4R0bMADm3ylLsB/D4zVwC8TEQvAngzAJlZEnoCjxmmQXjTDZPr9lmmsmhcn7GUK6PqBuJeUoudooH5Ur6CZ+YzcIOLBUNF8c9dzSFfcZv684LQKbb16SOiYwDeAODRYNMniegpIvoiEelvyiEAlyJPu4wmFwMiup+IHiOix65fv779kQvCDvF91K9aiqD9d9djOD7DZxWtZyPR/bkrKmPm/JUsClV1Q2pEnpevuCBSdwiAFCATukPL4k5EowD+CMAvMnMWwIMAbgZwJ1Rk/6vbeWNmfoiZ72Lmu2ZmZrbzVEHYFR5znbYfn0vj+FwaAGAQwSSC5/twA0um4vqh5w7UovelfAUAELcMxIMoXU+oTiRtmVwVukpL4k5ENpSw/w4z/zEAMPM1ZvaY2Qfwm1DWCwDMAzgSefrhYJsg9ATac2/k+FwaZ06dgGUSqq4finjJ8ZAtu+FxzIzzC1n8+l8+DwCwTYJOuOHg93TCBgOhrSMInaaVbBkC8AUAzzLzr0W2z0UO+2kAzwSPHwZwLxHFiehGALcA+Ju9G7Ig7ByfGQwAm+SimwZhLRKpX1gu4k+fuhJeELToO17tInHbwfFQ4OOWiZG4yrrRhckEodO0ki3zdgA/C+BpInoy2PbPAXyYiO6EClZeAXAKAJj5HBF9BcB5qEybT0imjNBtdCPsV+8fBQAcSCfq9kfrrlsG1U2c+sxwfUbCNlFyvDBjxvF82KYRLloyDYLvMeKWgVRMfbUKFReTqfp0S0HoBK1ky3wXzaefvr7Jcz4L4LO7GJcgtAUt2uYmi0gts/6G1vMZnsewTULZqY/c7cgLWYYBx/MQtw2YBoFofb68IHQKKT8gDAXMjGLVw/VgElRntzTrlBQtR2AS4AVFwRK2AcOgMH/d8XzETCNcFKWfF7eUJUNonhcvCJ1AEnGFocD1GZ7PWCsqL93YxHO3ItG4To2suB7+zmtnMTsWD9MjlS0TOTZ4nIjkt0vcLnQLEXdhKNCVHiuu+nezNnqWUfta6ON8BtJJO6wo6bPKaY9aOLXIvbZNInehW4i4C0OB9r4dT4ntZi1So7ZMtBxwOmEjFVdOpudzkPZY+wrpx3plKqF5LRpB6AQi7sJQUHbrE7Z+9Z47m/rtAPBLP/laACpbMnoNSCet0HLRBcKitszceAKpmAkiwplTJzCasOEH9s31XGUP/xpB2BqZUBUGGr3sv9KQtTIS2/ijP5Gyg2NM3HZwHN9/ZUXZMgkbRGrBkheKe33kHrV7iFRRseV8FRdWiljMlrG/IQVTENqFiLswsOjc9uNz6XXdlZKBd96M8aTKS9feu0EqQyadVKJvmUa48vQ/3vsG3HpgLHw/vco1+jydTbOYq4i4Cx1DbBlhKCg3lAFIbSLukykt4ioK1757OqFioZhJoV0zm45v+DpGELnrSVWdhikInUDEXRh4XM+H53OdZZK0N4vclbhry0Xb6gcnkgAAIsJI3MRtB9OY2GT1KUFF7npS9d/8yXmpDil0DLFlhIGnFFgyJgEeVERtbJYKaRq4aXoEY0GkbpkG3jiXxpF9qfAYIsJofPOvj47co4ueBKFTSOQuDDzFoOa6adY89K2YGYuHLfeA+onTVnn/nQdhm4aqHw8Rd6GziLgLA0djc4xi1YNpUGivbLaAaSvOnDoR1n7fiphloOr5kchdlQoWa0boBCLuwsBTqLhhWqNJ1FLkrtmOmDcSt0x4PostI3QFEXdhoGFmFB0vLMFrmwSjQ5/6WMOCJ706VhA6gYi7MNCoIl8Im2ckYmZYtbHdxMz61nsSuQudRLJlhIFGC6vOa3/VzGjb3quxnEHcro/cXZ/DFn0nT5/dsPyBIOwFErkLA0uh4qLi+hiJm2Feu2nQriZUt4OO3KMlgMWYETqFiLswEDRmyADAlUwZBODV+8fCVnidRHvuuYoblgGWCsBCpxBxFwaS8wtZrBWrMAwKRXavaDWDRnv7zMDRYAGUL+oudAgRd2FgYVaVGdvBmVMntvTMo0073nTDJOJWreCY5LsL7UYmVIWBhFk109iL6GWnE5/RO4ZkzMSBdAIXVorr6twIQjsQcRcGirLj4en5TGh/aK9dC3Sr0fJeZLJEI/eEbWJmLI4LK0U4ng/T6Ew6pjC8iLgLA8WVtVJYSwbYvJ1eu4lG7gnLhGkQLIPClEhBaCfiuQsDhdfQtLQbWTKa6GIpnfNuGhTUeO/WqIRhQcRdGChKDR2XeidyV491823XlyJiQnsRcRcGilLV2/qgDhEV93iwiMowCES1VauC0C7EcxcGAtf3UXW5rp0e0XpbppNL/usnVGuPLYPgeOK7C+1FxF0YCBbWyljMVeoid4OormF1p2mcUNXYpgHH8yR6F9rKlrYMER0hor8iovNEdI6IfiHYvo+IvkVELwT/TgbbiYh+g4heJKKniOiN7f4jhOEkWnKg6vlwfUah4ob7p0ZiXS3OFTOjtkztsUEAAXClBLDQRlrx3F0A/4SZjwN4K4BPENFxAJ8C8G1mvgXAt4PfAeA9AG4Jfu4H8OCej1oQGtBZMiuFKgBVBfLIZLKbQ6q3ZSwzLFtARLBMlRLpS/QutIktxZ2ZF5j5ieBxDsCzAA4BuBvAl4PDvgzgA8HjuwH8Fiu+B2CCiOb2euCCEEVHwcuBuBukGlt3EyIKo/do5A7UsmZykTsNQdhLtvXpJ6JjAN4A4FEAs8y8EOy6CmA2eHwIwKXI0y4H2wShbWj/ejlf6fJI6tF3FI0NQnT5gUzJ6fiYhOGgZXEnolEAfwTgF5k5G93Hatp/W/eXRHQ/ET1GRI9dv359O08VhHW4vsqS0ZF7r6CTdRJBKmTUmjEJyIq4C22iJXEnIhtK2H+Hmf842HxN2y3Bv4vB9nkARyJPPxxsq4OZH2Lmu5j5rpmZmZ2OXxAA1CLk5UIVSVs1w+4FdDPueJOyw6ZpoFD1sNpjFyRhMGglW4YAfAHAs8z8a5FdDwO4L3h8H4CvRbZ/JMiaeSuATMS+EYQ9p+x40POSy/kKkrHeKcrVGLlHMYN9LyzmOzgiYVhoJc/97QB+FsDTRPRksO2fA/hlAF8hoo8BuADgnmDf1wG8F8CLAIoAPrqXAxaERj780PfCxyuFKiZTsS6Opp4wcrfXx1F6gdWn//gpTI/GpaeqsKdsKe7M/F2otNxmvKvJ8QzgE7sclyBsia6LHl0MVKx6mBvvnaoaYeQemVA9c+oETp4+i3NXMgAQNvAQhL2kd74FgrANLi4X8diFVeQrLjy/Xhx7yZbRkbttro+PiFSz7mbi3qwnrCBsByk/IPQll1eLAJTf3tjVKGH1krgHK1I3KD0cMw1UPInchb1HxF3oS/TiH9X4on5fMmbitz/2li6Maj1EFEbvzYhF+qoKwl4itozQN0StCl1DxvU49NyNTTJTuoWK2tdv1/nu0abZYsUIe4mIu9CXaHH3fD/03KdH4wCAZA+Ju2UQLGPjr1nMNKTGjNAWxJYR+pI6WwZqOf9owsJirtJT4n54MrWu9V8UXRa46vlISNNsYQ8RcRf6kqgtw6wi5GuZMoD6xhjdJtZkZWqz/RXX7yk7Seh/eudbIAjboFBRTTlcX3nulkEwAtM90UOpkJtx5tQJ/Od/oNodyKSqsNeIuAt9ST6I3ItVF7myA8skmMHMZS/ZMlsxPabmCdxIrr7vM65ly/ClDZ+wC8SWEfqSfFmJu9a/uGWGWTP9JO6jMfUVjPryP7i0hleWi7C7XI9e6G9E3IW+4vxCFidPnw29ai2JMcsAO8qq6aUVqlthGLSuhML1nKpJ7/pi1Qg7R0IDoW/wfAYHoXq+oYNR3DJqnnsPrVBtBcugushdtwqUHqvCbpDIXegLTp4+iycuroalBrQto4lbBkrV/ppQ1ZgGwfUY2bKDquuH3aRcyX0XdoGIu9AX+MyqZnsgeIWKi7GEhVwg8nHLhGGox/3kuQO1yP1qpox8xQ27SUnkLuwGsWWEvqASpAr6DJy7ksG1XAWHJpLhfrtPs2UAhJ6746mfrz2pGpeJ5y7sBhF3oS+oBJOlgBJ4z2ccnlTiTlQrnwsAyVh/fax15O4E1SGL1VoOvyDsFLFlhL6g7NbngQMII3dddVGLeyrWOx/rVrorqci99vfpuxSxZYTd0F8hjjC0RCN3L8iYOTSpxV1tn0jaODaVwmsOjHV8fLvBNAz4HE4nhJQdD7c/8E2pFCnsCBF3oS8ou37Y+EKnDU6PxjGRtMOI3TAIs+nEho0xepV//L/e2HQ7A2HqpyBsFxF3oS+oOD6MoPGFjnBH4hZuPTCGOw5P9HVz6XTSXrftQDrRhZEIg4SIu9Dz+D6j7HogQl1LvdF473jruyGdWC/uN+8fAVBbgbuXSFOQ4UDEXeh5lgoVMKuJUysi7iODIu6RyF3/eTfPjALAuhaCgtAqIu5Cz5MtqcVJRyaTuO1gOtw+MJF7svZ36BIKNXHfnbpLlD68iLgLPY+uI2MaVDdZOhLvr8VKG6FtGYOAg+PKa79lNhD3Nr5vVPjlIjB4DEboIww0ubIDoOa3J20DpmGEPVP7nfHAlrFNAzOjccQtE6/aL7aMsDtE3IWeRxcJ0+JumQaOz6UHpt55Kih0ZpkEyzSwbySGiWQMgOqtei1b3vP3PL+QBQAcn0tvcaTQrwzGt0MYaHRxMLPP8tdbhYKJYtuofR1jloH9Y3GAa/XdBWE7SOQu9Dy5Sn3kvhH9nOs+GrfWTRDfOD2C1WIVVa9WmkD74rv5W1cKVfjMYdmGKHvx+kJvIOIu9DyNtkyUQRGhWzcomUAAqh7D9XxYe2BDMTNeWMwDAMYS9V9/sWoGiy0/LUT0RSJaJKJnItseIKJ5Inoy+HlvZN+niehFIvpbInp3uwYuDCbNsjZyZQdGUPnxzKkTQyU+Ojvojn/9F3uSzVJxt1dGWLJo+pdWQoEvAfiJJtv/AzPfGfx8HQCI6DiAewHcFjznPxPRYOSrCV0jX3G3tGQGFf1nN2bNtCq682ulsG0fUCsnDKgGKMWq2+xpwgCwpbgz83cArLT4encD+H1mrjDzywBeBPDmXYxPEJAruzi6LzUwFsx20JG7v8OcyKuZMpbytQnZQqT3bLHi4en5LD70uUfwylKhro+r0P/sxsT7JBE9Fdg2k8G2QwAuRY65HGxbBxHdT0SPEdFj169f38UwhEHj/EK2LirNVVyMNqm/MgzQBpF7KzAzPJ/r6sKXIqWT9dZMycG1XCVsFiIMBjsV9wcB3AzgTgALAH51uy/AzA8x813MfNfMzMwOhyEMA/myg7EBKTWwFWdOnai7Q9Fm1E7KEJQcDwzAiTQCidoymtWiWiQmkftgsaNvDDNf04+J6DcB/Gnw6zyAI5FDDwfbBGFHnDx9Fs8u5PDOW2sBwDDZM0QEQmtlCBrTGNcC0Y5G7tpjT9gGDCIUqx7WisqT91ldRHZqAQm9xY4idyKai/z60wB0Js3DAO4lojgR3QjgFgB/s7shCsNI1JrxfB6YImHb4cypE3j6gXcjFTND4V2M2CeN9lUjmVIg7j7jns89gpOnz6IURO6qNr46zomIf6Hq4fELq5IhMwBs+Y0hot8D8E4A00R0GcBnALyTiO6ECiheAXAKAJj5HBF9BcB5AC6ATzDz+vtAQdiAapNUPdf3MTaknjugVquWHA+XV0u4kinDmUyG/WM3Q4s7oATeNqnOlqFA4H1Wi6jyFRfM6kuto/er2TLSieG7sA4CW/6vMfOHm2z+wibHfxbAZ3czKGE4eeFaDj+4tIakbYQLdpRNAIwOuMBsZjXZpuqxeiWjasz4zHA8P2wUvhHalgGUNWObCCP3V+8fxUtLhaCzFWM0bqLievBZ3SlVHB+ez7iwXMTcuHSF6kektozQM+iVk9omOL+QDVdNDsuEajOmR2OwDMJsOg6CEuqLK0UUnc1virORyP2F6+rcas9d143X1kzCNpGwTcQtJQklxwurccpEa38i4i70DFfWSgCUhaCzQ7SsNC6VHybGEjaSMRPHpkYQtw24PqPi+mDefMVp1JbR51NfEHQRNi3yCVutNdRiX3a8sEmKTLD2J8P7jRG6TjS74+Tps7iwXAj3aY9Y68qg2zKtYhkGXM8PJ1XzlY1XmK6VaitT9XksVVUvWiJVQ8b1GdcCX30eteyckuMhK5F7XyORu9AzVFwfCdsA0fqJ1WHMlmmGbRIcj0PrShdVA1S0/cx8Jrxo/sFjl2vlC4JjilUPKdsMV75aBuHQRBJf+fjbwpo9hkFB5K7EXbS9PxFxF3qGqusjbpmIBROIXsSeGeZsmSiWYcDx/DCajkbuLy8VUI748J7PsE0DJlHNlql6SMasdYulohgElBw/tHW2mrgVehMRd6FnqLg+4pYB21RRZdX1oVOwD0jGBgDVrcn1o3npbniX4/oMj4FzV1T07vp+0N2J4HqM56/lkC07YeenjTCIVKbMShEA4Inn3peIuAs9ge8zXJ8RswwQEexAxBzXx2jcaimvexjQFz5AWSrMwEJGTUTraF5LsesxLMOAaRAYqszAS4v5TcX9+FwaRybVuX4xyF6SyL0/EXEXegLdbSgW5LfHgpQ8hkoFFBQ//2O3hI91pstCkP+uxV2LseszLIPC3HYAuLBcRHKLyF33pv1hkD7pMUtd9z5ExF3oCXRK32d/+nYcn0vDCPqKAsDUiIi7Zl/kXOgg/pe++jROnj4b2ic60PZ8hmkQbpweCY8tOd6Gkbv24e3gwqojd50mKfQXkoIgdJRo+mOmVKv2qH3jqP2SsNXE6l60l+tnohOfT1xcDR+/di6NJy6uoer68JnDdEf1mMPIfWYsjmu5MgoVFcGnYpt/7XXkng0ycYpVD8wcZtgI/cFwf2uErrGYLeO5qzlcDxpJaHGfHY+HxxDR0HZg2oh9KRW5G6RE2DQI1Uj2DKAidv27FYTshFr54CcurGIzTKotZtKvx9i6UJnQW0jkLnSFpbxaYKNFver5sAxC3BILYDMmA1smnJswDRW5R8Td59r51VG6mqQ2UPX80KsHmte00cdGV782Tqo2lhcWeg+J3IWuoHOo9UrLqueHk6jD1gR7O+gKjXZk4lmljCrxjQWR+oWVIpK2WVfRUWfamC3cDNkNVpikQ/YfIu5CV8gES+OrQSK74/rrBEVYj04T1XZLLIjGtQ1z4/RIKOJz44k6n1xPlBoteOf6NXSQ70sHvr5Dvk1CV4hG7idPn0Wx6oVWQ5Tjc2m59W/gyGQKB9JqUVfMUuUI9MIm0yDELQM3TKUw1ZBCqs+vscU8xplTJ/C+1x+se07J8eBKj9W+QsRd6CiuzyhVvbDWeKnq4dyVDBhKqDRizWzMzFgc6aSNM6dOhHnvuuyAaRCICAfSiboI/fhcGn//LUfVMS1E7jNjamJbR/vzayWUHL+uvIHQ24i4Cx1lIVPCuYVsGLkzahULm0XuwubosgylqoqqNxPu/WPq2K0id6Am7vr/RIv61Wx554MVOopkywgdxXG5rm4JUJus01Gi0Dpz42pdQCkQ32bCrW2tvzyv+tq3MqE6MxqIe/B/opNlrucqyJUdlB1P0lR7HBF3oSPUml2rCPP5q7lwn06z+7V77qx7jnjtW6Mj96gtsxE3zowAqAn2ZhydSgEAkpHVqQQl8ueuZPHc1dxQN1DpB+R/R+goOkp/ealQy/QIosLZtFR+3C7phAWD1FxGzDQ2zYS5eWYUdx6ZCFvpbcarZ8dw+6F0XQaTaahibheWC6i4PmxHJlh7GbkPFjqKGylqdVMQSXo+gyA1ZHYCEUWsE66biG5Ws70VYdekYladh6/vCv79N/8WAFBxZXK1lxFxFzpKdJn8q/aPho9t02hpok9Yj570bIcHHr0RUCUPCLmg5ozjMT704CN7/p7C3iC2jNBRouI+PVqrI3NoUiyZnRJrWJy0l3MVRASDlNeu7hIoLEAG1Eo1SzmC3kPEXegoUXFPJ23cOjuKS6ulME1P2JpGAW1n5K5f1/cYREDcMuvEvdLQ61ZEvncQW0boGMxc12x5ImljIhWTlLpdoiP3dp1HgwgE4La59DrPviKLmnoWEXeh7fg+46nLGVzNVuq2jyel6fVe8C9+8jiA1vLXd4JpqCqR0clb/V6NkbvQO4i4C23nSqaEkuOFq1J1UaqJlIj7XqBz3dtpy8QsA2dOncBnfkpdSCZHYohbBhZzFanx3qOI5y60nY996TEAtVWUSduE47kSue8Rc+OtlxXYjheuj33fb/x1uO3QhFrcFLdMVCwfxYq7naEKHWTLyJ2IvkhEi0T0TGTbPiL6FhG9EPw7GWwnIvoNInqRiJ4ioje2c/BCf6BXT+rGHBMpGzHLwI3TI90c1sCwbySGmGm0rdfpSNzCSNAO8dh0CuNJGxMpGwnbhMdqLkVonU41G2/FlvkSgJ9o2PYpAN9m5lsAfDv4HQDeA+CW4Od+AA/uzTCFfkR/iBsrCY4lbLzhyASmIqmQws4hItxxZByzY+0/n3HLxGsOjGE0boX9b1eLDh67sIpCxcXl1RLmV0ttH8dO6JSo9gpbijszfwfASsPmuwF8OXj8ZQAfiGz/LVZ8D8AEEc3t0ViFPqXUIO6NlQulZvvuMYg63sBa15a5uFKE5zMKFRfLhQqWC9WOjkNozk4991lmXggeXwUwGzw+BOBS5LjLwbYFNEBE90NF9zh69OgOhyH0A+WGGiTmFj08hf4gbhkg1DJmKq4fWm/M3PGLjVDPrrNlWBlu2zbdmPkhZr6Lme+amZnZ7TCEHsXzWTVljnzPJa99MCCiuv/LYtWDz2o1q27GInSPnYr7NW23BP8uBtvnARyJHHc42CYMKcWqyqbQFR8JgGj74GAaaoFT3DKwVqoJ+s9+4dHuDarHYeawMXw72am4PwzgvuDxfQC+Ftn+kSBr5q0AMhH7RhhCrqyVYRqE6aCfp24DJwwGtkm488gERuP1Dm8vLm4qVl2sFbszHxCdzF0pVPHkpTVky+29u9nScyei3wPwTgDTRHQZwGcA/DKArxDRxwBcAHBPcPjXAbwXwIsAigA+2oYxCz3EZrVE8mUXayUHhyeTYZqeWDL9x2bzInrV6rqyBB0U91br2VxZKyPXA3n5ZdeHz6qrVTrRvrUeW4o7M394g13vanIsA/jEbgclDAY6Mpkdi8Mggm0qj1YmUfub6P+fFtbG7k7VHqz17vkMz9tdTv5eFEbTxfMypS5H7oKwUxyPYRBgBVULbdOQyL1NdPuC+Zmfug0f/dL3kbANEGjLapHdqB7p+QyPVQ/fdn4Ot/rbXE/EXehzXN+HZdQiuqP7UhC7fTA5OKEadcctsy49spfQLR5zZQcTqe51/XKDPsLZNou7FA4T2obrMaxIqcLxpN1Wj1HoHgcnVDZUzDIQs4zeFPfADsmW1vvuraxePXn6LM4vZHf0vq8sFVCqenXjaHe6qETuwq6YXyttmNro+D4ssWGGgrGEjamRGPalbJQcD57PWCtWOxIhFyourqyV4Hh+XUPvRvwgcm93lopGXyx+5UN34FquEo7NFc9d6AeW89UNrRbXYyTi5jqfVRgcor6y7onLQbrh3//N7+Hrv/AjbXtv/XlaKzlYKTq4li3j8GRqw+Mtw4DjeZvaIVv55cwcrsJtFV2OQd/NTI3GUFwptV3cxZYRdoXj+Sg7XtPKgK5fb8sIw0EySHuNlp0oVFysRnLMF7Pl0KbYLXpB0GZi6fkc1jjaTeReqnr4waW1TSthVlwPfqTl2HK+Em4HaraQRO5Cz+J6fniLuVZ0MDkSCyOf//KP3gLPZ/yDt9zQzSEKXUDnvEcLxs2vlZArK1GruB5eXi7iQHpv+uY63sZeukavlG48bjt3k1XXh86kLDkeUrH18snMeGY+i9l0Aocn1STzcl5d1HIVF/d87pHw4iKRu9CzrEYmhH6mYbm5jtImR7qXlSB0ByKCQagr9ex6DNdn5MoOFoN2i56/N5OurUTu0abezY6ruj5eup7ftCzA1Wx509cAgHzFhRu5SwBqtgyzuoPQQb9E7kJP0CxPORoNNWZHfOxL3wcA7ItMqHU7F1voHAZRnS2jRfPDD30P/+qnbgNQm1jcDJ2hsllZaP3a2bKzLhLXz8lHVqY2s2UyJQdL+SoMog2byESft1Z0MDeeXHfMakEdU3X9cIW2tmWA+rsZSYUUehYnstqvEvnQnl/I4sXFPABgckRSH4cRHbnrtD8n+LfseFjIqGYeXoO4b7eZxivLBSxmyxFbZrPIPWrLrD9OXyAWc5W6oCWK5zH0DNJGUfdKUU+eeljMVTC/VsKFlWJtHME8Qzph4YfXC21NMpDIXdgx0VvYq9ly3QdVf233iS0zNEQj6x/9lf8PLy8VcHm1iAPjiVDIy66Pqxllb7QSuW/GUr6KpG1G8tc3Efeo515eL966LLXPwEqh+eu4PoNI2Ssb5aivBhaM4zGWCypi/8HFVYzFLeSClE0AODqVwvkr28+Z3w4SuQs7RkdMcctA4/dU+4r7urgSUOgeelJ1fq2ElUhnporj4Yv/42UA6yP3rYhG9oWKC8/nOrulMZo+v5CNHF+7s/yr5xbDx66v5gKqrh8WQLvWEKgAgB8cZwR5vxtdSKJ/q/7zlvJVvO7QeN3ffGQyFdS+b1//WYnchZYoVNywRozmaqYEApCKmVgrOfCZcWFZtVzTH9luLvMWukcs+Kxcy5YxFq9Zc2XHD9NjG4uNbYfo5KamWUSu0baMbVLdReWlxTwYanI3FtQ+ama5aL/9QDqBy2slrJWalw5e3aCk8KHJJCyDwruVo/tUPv52L3DbQSJ3oSVeWMxjflV5h6Wqh1eWC/AZsExC0jbBrFK+FnMVVFwfzAwzKAcrDB/6//1atoKlwJ5I2AbKrodqYOflK25dPngjVdfHheUi8mUXi7lK3b5rmXpxN0hF7rmyg2euZPDMfKbONtQRftwy6uygQtVFoeKi6jJiloGkbcJn4NyVTF30rm2YmKUuABvZMisb9I+dGo2FPWcB4Egg7u4uq1RuhnzzhC3Rq/K0DfP4hVVcy1bg+gzbNEJfXU8cecEtrCxgGl5MQ6VDXs2U8W/+5DwAYCxuw/E4zKJhRpj7ni07ePzCap04PnMlg6vZMhjrbZDGyP3wZArLhSqeXcih6vqouD7Kjh+mY+pJ0phlhoW7ChUXTpCiWfV83HPXEaRiagFW9Jpz8vRZfOJ3nwAA/Mv3HcdE0t5wQnW1WEU6IuJ69fbUSAxTkfmnQ0GhNXeP0kGbIbaMsCE6cvn8fXeBUfsg6mwHQJXxTcVMGKREfSJlY63ogLl26ykMJzHLwGKuHEbQM2MxXM9X4PmMmGmg6vn4yBcfxQ+XCjg0kQzy4GvWylIQrROUaOqUSAD4tW89X/deN82M4OxLy8HnlMMZ/aV8FSdPn8Xl4K4zFrFlGn312XQCyVDc6yNqN/gbxlM2xpN2XUvBaJrwSqGKufEkitU8fGZMjcRxPV/B1Eg8sCgL4esAwMtLKmOmHWnCIu7CluhoSt9CLkRuiW1Ttc2zTVUJcP9YHFXXB0EyZYadmGngO88voeJ6MAgYjVvhxT8VM1Et+aHQak/8er6C2x/4Jo7PpfGBNxwCoO4CXJ/BzGGA4bg+TIOQsAyUHR+zY4lwrUVUl5fyFRyaUNk6IzETlqEm/x3PR7lhbcZsOoE/+rm349Z/8efwfa6bkNVWzkTSxnjKRra0PqceAB55cRkgZf8wlBUFKFtGtZhU47txSuXSt9GVEVtG2JpQ3P2auOvbTTtsxEG4aXoEE0lb3ZJLNcihJ2ap6JxZfU6ICIcnkiBC6D/rz5TO/45GzHrxj26sUXZ8PHU5Az+wUWKmgdGEhWTMDCPhuvc3VdOQfMXDcqGKqueHr5UpOXUraAHgwHgCMcvASNyC43Fd/Rgd2EymYvjh9QKeuLDa9G92fB+2QTi6L4Ub9qUwEvSW1XexIzETbzo6gcmRGBKWsemcw24RcR9yTp4+i9sf+Oa6KCRauzoq7r7PuJopIWmbiFsG9o/FAagl5zNjcWl+LQBQFsUH33QEjuuDgXD+ZSRu4a4bJjGeVGKsI/diELkzIxTVpXy1LlBwfYbjMRZzlTB18ei+FF57YKze5w5+YpYBIhW9M6uLhB1M9F7PVVB2PBAAPTU0m1af5RumUmDUr7p2fR9EQDpp12W9NKJ6GBhIJ22kk8rCufPwOG6aURUziSjMOhuJW0E5gvYIvIj7ELHdFYDMjIrj1U1yfehzj+DRl1fCpgy68bUgNHIgHQdD5YhH66wbRKHYX1wpgplRdv2wmqTWzeVCFXYwMRvlwnIhiNwJBinx1xcLg2oNQ4gIk6kYVgpVcPC+seB9r2XLqDg+RuMWpkaVqM8E/47GLdgmwfEYjudjpaCywAxSPYAts7m4+xskEsSDv+vMqRPhnAEAjCYsMBBmD+01Iu7Chrg+48nLGTxxcbVuW8X18VOvPxh+UKMf2sYPsDC8zAZVHxm1MsAa3X7R51r2zNRoLNwGKFvmdYfG8cwD765r+vKpP3oKjsd1abbpUNxV+u2dRyZwfC6N6dEYXL/WN1Xn3y9mKyi7HhK2gYMTCdyyf7RuHYe+GD13NYcXFvNwPA7HYBkURtznF7J4ej6D1WI1zIW3jdZkdTSwbPKV9jQTF3EXQqKRfbQpwf94cTk8puqqSbAD43tTrlUYXGYjn5HZdAJnTp0Is0K0VutUxKmRGP7rJ94OQFkbhYqL5XwVU6MxEFGYkQUgzFSJW7ULhhb32XS8LrgYT9rhhYUZoS2jKkAyEraJuGWum/zXqZzFSM15HZGbgXhrS6ni+nh5qRDe4baaAqzTLvObLL7aDZItM+QsZMooVlxkSk54a3vy9Fk8fSUb3h5fjBQ+0h/2uU3EXao/CgBwMKiaaJsUliPQ6LLAPqv9r9o/iunROOJB/9VzV7IwDcKbjk0CAG6eGcXzizkUq164gEhHvgDCz25U8AEVyd88M4JnrmTBrMoHWAbh0ZdXAGBTW1GnayZjJgoVL7zb0BH8U/OZMIr3oXx8tb+1mPkPPv42vOfXvxOmX+41Iu59wFatv1rF8znM19WsFavwWN1+3nF4PPyw6+O096i/iFrcGyN3EXShkQPjCbxq/yiuRtZFRNHBQ9Ryec2BMZxfyIYR/XQQUccsQ/nrRPBYWSQ6zRAAbj80jkMTiab9A0biFpK2gZuDSU3bNPDkpTUAqHsNzZlTJ3Dy9Fmcu5LB7YcmcT1fQaFSDMepxVjXqAeU9fR8UAk1ZtGGrSUbvyf6jqMdiLj3CTrndjciupAp4Vq2AmYOs1oqQb6w7zOuZcu4YWpE5RMHHmPSNuF4LhK2iWLVCyvdzY0nRdCFLZkaieFaw2pS/bk59qk/A1BLdQRUJG2ZqtcpgHCyEwCOz6Xx8lIBi7kKxhIWiKjuM9jYPzUqsJZpIBlT/Xzv++Lf4L8/fx3A+kg/CjVM1mq7ZTRu4X+5YRKXVkt1K2V1lcfo5HE3vyMi7kNExVVt8TIlB6d++/HQV7dNwsRIDNfzVRyeTKEapK+ZBiFhm8iWa+Lus0oz28yWEYRWuGFfClez5XXps1HLWk+yanSknU5sHPFuJag65ZFQf2FpfA0ddcctQ31HIlG2YRCmRmN14v7sQhaE+juRbiITqj3CRmmK0Xzz6Iq5neA0rDB1PFW9Ueeoez4jW3bCuhmvOTCGj/9vNwNQ9oz+IozELUmBFHaNXjTUiPbjAWBqJF63byToW7obO0Nn8ej3iE70NoNIBTmNFU5Hgkle/b147mo2XKzVC4i4d5FGQd9IvKuur6LpFhc7bHSh0D66rg2js2EMAhLBl8zxGJmyAyJ1i3zTjFombRlGGJFEq9sJwlZs1iLv+FwaTz/w7nX7dd30xsg9nbTxxqMT+JOff0fLlkdjeu7+QNz3jcRbeo2N0nt1Fo++myg7PmyrN4Qd2KUtQ0SvAMgB8AC4zHwXEe0DcAbAMQCvALiHmZuv1RW2pFh1w5Vyjre+oFGz3qbRAktRouUDAOClJTUB9KpIjq/r+ag4Pszgy3XL7BgANaFlGYQKRNyF9nJ8Lo3FXBmvLBfDKDv6+bbN3cWks8Gq6maTqduFiHD8wBi+/8oqGMDbbp7Gb37krpaf305Pfi++pT/KzEuR3z8F4NvM/MtE9Kng9/9zD95nKNFCrFPEdOf4ZkQtnGbo6nxXM2Uw17qw/84/fivSCRs3ffrPwlV5+0Zi4Qfv+FwaIzETK8FkajQFTRA2Y7vipY+/53OPYDxYvr8X7xXddyxogJ3aIgVxq7FHfXk7SJvU5Th6gXbYMncD+HLw+MsAPtCG9+go21223+prnruSAbOq17KYq8Dn+kp0a8UqlvNV2CaFnvfltRI++OAj2x5P2fHC1LPfPnsBj19cU4WUiMLJKSv4gLpBSVbNWMKCEWTOmIRdR06CsBVEtGkmy2549ewYbj+U3taFI0ozj94OZoH3j/VOosFuQzAG8BdExABOM/NDAGaZeSHYfxXA7C7fY2CpuD6eu5rD4xdX8fKSqvMcvVV86bpqAWYZqqxuzCSUHB+5srNp+7qy4+Gl6/nw95Onz9YVQYrWoo5F3s82CKUgj91ussru6L4UjuxLSQqksCd083OUim1P+rYaqwp4POxP907kvltxfwczzxPRfgDfIqLnojuZmQPhXwcR3Q/gfgA4evToLofRHjaKjrezqGijY72gwh2zF+YBExCszlMtvu656wiA2uSSrgedKbmYSMXUStL5TOghMjOcoGPSasGpy2fXk6m6nrQmunLwzTfuw1+/oBy2ZtE5EaF3posEoXfQwdDMaO+I+67ur5l5Pvh3EcBXAbwZwDUimgOA4N/FDZ77EDPfxcx3zczM7GYYG9IOO2WnNI5FN9J1fcbCWs1XZ67Vjr6wrKL54wdVRsFtB8cxFrfqWnxVXB8vLxdRqLg4v5ALGxB4zFgrOqEPrydTU0EK43hQdz16a7pvJB4e9+8+dEdbzoMg7Iat0hY7jR6PDoZ6KXLfsbgT0QgRjenHAH4cwDMAHgZwX3DYfQC+tttB7obd5oYDwKWVIp6/ltvyuI0uJucXsusmOpfztTK6z15V+yxT5fdWPJX2+PJyEfFg2bVmPGmj5HhBT9NaJ5v//vx15CsuYmZtWXa0JoyeTNXNA8aTFlIxs66MQDTtbCYyMdRrXyhB6DV0vn4vFdTbjS0zC+CrwW2/BeB3mfkbRPR9AF8hoo8BuADgnt0Ps7NogS5VVUnQtaBrCzPj3oe+h2fmMyDCjssBfPDBR7BWckKL5NmFXNiuLha0DXtqPgPzSnZdwaXxpI1LqyXkyg6cSE3pv35BLaeOli3V4u4z4/Kqym0fjVtYzFUwnrSx2tDBPdrAd3pUWuQJ3aPfgonp0TgStjkYE6rM/EMA6+7dmXkZwLt2M6id0OhtMzMurhTXFcpqlWLVxdPzWRybSqFU9cAAVosOKq6HQtVrOuHYOBZmxg8urcEyVKaLjrJzQd3nQxNJXF4t4YVrOdx6YAyjcQvnrmRgkFpgRET48JuP4P/+wO3h31Z2PLzmX34DJcdHtuyAoGa1v/O88sqjK5//7Teew1jcQiFSL3p6NIaxhIWvffId6+4ydNlT02ieqdBvXzhB2Ii9/iw3Wpy9wMDmtL3/P30XC5kySo4fNt/dDteCfPKLK0XdSB0XV4q4mlHbXW/r9lj5igvHY5QcH8Wqh8cvrOKezz0SeuZaTF2fMT2qVsvddnAcMUs18fV8xrGgka4mEbS3K1ZdFCouLFNNcs6vlcLo/7aD47BNQr7s4oXFWtaMzrrZqHTAdDAZtNmFSxCE/mCgVqNEKyfqms+Eeu8Z2LwMp84/L4YNe2vHPTOfwWKuHEbL+chF4+Tps/jBxVXcMDUSivZS4KvrcrkAUKh4yJTcIDquXVujHrdl1Gpd39Ag7oAqOaotlWOTKcyvlVBxfSRsM1zq/brPfDNMeYxZBqquv66JQGP0oscteeyC0P/09bc4OoGpPXHNWtGBQWqSMtpNpRUcj+FzvQcNAJ/9s2fhcy0X/YXFfN37Vz0Oq8RVXNVx3TJU/YmRuIqWr+XKKDleIOCEyaBru46ao9E7AXjt3Ni68UVblo3EzfAi8d7XzYWCrbeNJ23Eguh+q/ZfWtxjIu6C0PcMxLf4Q597BE9dzqDkqCyTtWIVuYobCqjrM5byGy/bB9SF4vYHvqnSBj0fCdsIy9ombRO2SSg5HszAP7cMgutxmJOuXz9XdvE/L6/hyUsZeD6HVokRVLpbyldBBNx2MB146GpOYKZh2bJtGnjTDZPralQDwD/98VvDx0nbDG2WG6Zrx+oL0OxYPLBiDByaTG56DnS2jNgygtD/9L24M6v65AzlUT89n8Hfe/ARAMCr9o+FnvXPfP5RADX7RPc7bMRnhsfA9EgcX/vkO2AGkbeeYBxLWLjt4DhumEqpzu6sxnA9ktpYcXzELQOvPzRel72iS4POjMZD60MLabPslI1qTd8yqzrKmKSXaavXujFi4UyNxDE3nsBEqlZaYDxpb5rWmIpZOJCOr+snKQjC5vRiunBfi3uu7OCp+QwW1pQPTqQW9awVHVgGYTRuhsWBSo6yZkpVZZ8s5SvIlR08fy0Xlr6NNoWeGo3BNAivOTCGI/uSoYDqOix6+b/j+VguVFF1/TDXlaEsjmRMeeD659bZMdVtPZILq0W+Ma+8WVVHjW4XNjOmJmH133jrgZqFk4yZOLovparWBe/fCjdMjWBsk0YIgiD0B30t7gYRfB/Kggma8PoMLBeqmEjZICL88f/+NpiRmimLOeWJV10fq0UHq0UHd/+n7+LclQwqwdJ926xllIzGLcStWs3m8aSag7YCa8bxVMplKmaG3jYA7IvUftFX9ZG4hTsOTyAe8cz1BWE7y5ZH4hZunR3FXNCAeDxp447D47gpEP0ovRhRCILQfvo6W2YkbuG2g2lcyZSQL7sgQpjJoltiERFStonlfAUfevCR0D6pBA0wAFVIq+T4oUf+hiMT6wRxelR51w9/8h0gIpw8fValLDoeHI9x03QKEykbz1/LYa3ohBOozYi+tp441Q0Emh3TjGjhsM3SGwVBGE76WtwBJY7RXPAnL63B9bluQUEyZiJXcbGYr8DzGaahJlm1VXMtW4bPKsMkZjVvk5WwTRyaSIb7oh3Sb5oZDTNYbpwegeczvvLxt7U0/v2jcYzFrbYvgJDoXRCGi74X90Zsk/D6Q+P4g5+rievMWByLuQouLhcDzzuJHy4V6nLZDWqeJRLtANMMIqpLTbRNA9sJog2Dwnove8lOxVwuAoIwGPS1uDcKb7QzSpTRuBUu/98/lqhbPPSTt8/hz55eQDyI2Dfr97idsWx2zF4jgiwIQiN9Le6aVsQtbhmIWwZmRmNhWVsA+PHbZnE1U8LzkWX6rb7HRheTTiCCLgjCZgyEuLeCaRBeHTR7NoyauB/Zl4JlGruO2Df6fbvPFwRB2AsGTty3irD1/lt+6etwPMbhLVZtCoIg9CN9nee+G2KWAYN6qy2WIAjCXjFwkftGNEb0KdsKVrWSWCOCIAwcQyPujRwLasMIgiAMIkMr7sYGRbkEQRAGgaH13AVBEAYZEXdBEIQBZGhtGZlEFQRhkJHIXRAEYQARcRcEQRhARNwFQRAGEBF3QRCEAUTEXRAEYQARcRcEQRhARNwFQRAGEBF3QRCEAUTEXRAEYQAh5u7XRiSi6wAu7PDp0wCW9nA4e0Wvjgvo3bHJuLaHjGt7DOK4bmDmmWY7ekLcdwMRPcbMd3V7HI306riA3h2bjGt7yLi2x7CNS2wZQRCEAUTEXRAEYQAZBHF/qNsD2IBeHRfQu2OTcW0PGdf2GKpx9b3nLgiCIKxnECJ3QRAEoQERd0EQhAGkL8SdiBJE9DdE9D+J6BwR/V9NjokT0RkiepGIHiWiYz0yrn9IRNeJ6Mng5x+1e1yR9zaJ6AdE9KdN9nX8fLU4rq6cLyJ6hYieDt7zsSb7iYh+IzhfTxHRG3tkXO8kokzkfP2rTowreO8JIvpDInqOiJ4lohMN+zt+zloYU1fOFxHdGnnPJ4koS0S/2HDMnp6vfmmzVwHwY8ycJyIbwHeJ6M+Z+XuRYz4GYJWZX0VE9wL4twBO9sC4AOAMM3+yzWNpxi8AeBZAusm+bpyvVsYFdO98/Sgzb7SY5D0Abgl+3gLgweDfbo8LAP6amd/XobFE+Y8AvsHMHySiGIBUw/5unLOtxgR04Xwx898CuBNQwQ2AeQBfbThsT89XX0TurMgHv9rBT+NM8N0Avhw8/kMA7yIi6oFxdQUiOgzgJwF8foNDOn6+WhxXr3I3gN8K/s+/B2CCiOa6PahuQUTjAH4EwBcAgJmrzLzWcFhHz1mLY+oF3gXgJWZuXJW/p+erL8QdCG/lnwSwCOBbzPxowyGHAFwCAGZ2AWQATPXAuADg7wW3WX9IREfaPaaAXwfwzwD4G+zvyvlqYVxAd84XA/gLInqciO5vsj88XwGXg23dHhcAnAiswT8nots6MCYAuBHAdQD/b2CxfZ6IRhqO6fQ5a2VMQHfOV5R7Afxek+17er76RtyZ2WPmOwEcBvBmInpdl4cEoKVx/QmAY8z8egDfQi1abhtE9D4Ai8z8eLvfazu0OK6On6+AdzDzG6FujT9BRD/Soffdiq3G9QRUfZE7APw/AP5rh8ZlAXgjgAeZ+Q0ACgA+1aH33ohWxtSt8wUACKyi9wP4g3a/V9+Iuya4zforAD/RsGsewBEAICILwDiA5W6Pi5mXmbkS/Pp5AG/qwHDeDuD9RPQKgN8H8GNE9F8ajunG+dpyXF06X2Dm+eDfRSgv9M0Nh4TnK+BwsK2r42LmrLYGmfnrAGwimm73uKCiysuRO9U/hBLWKJ0+Z1uOqYvnS/MeAE8w87Um+/b0fPWFuBPRDBFNBI+TAP4ugOcaDnsYwH3B4w8C+G/c5hVarYyrwTN7P9REYlth5k8z82FmPgZ1C/jfmPlnGg7r+PlqZVzdOF9ENEJEY/oxgB8H8EzDYQ8D+EiQ0fBWABlmXuj2uIjogJ4rIaI3Q32n2x7UMPNVAJeI6NZg07sAnG84rKPnrJUxdet8RfgwmlsywB6fr37JlpkD8OVgltkA8BVm/lMi+tcAHmPmh6EmUX6biF4EsAIlHr0wrv+DiN4PwA3G9Q87MK6m9MD5amVc3ThfswC+GnznLQC/y8zfIKKPAwAzfw7A1wG8F8CLAIoAPtoj4/oggJ8jIhdACcC97b5IR/h5AL8TWA0/BPDRHjhnW42pa+cruED/XQCnItvadr6k/IAgCMIA0he2jCAIgrA9RNwFQRAGEBF3QRCEAUTEXRAEYQARcRcEQRhARNwFQRAGEBF3QRCEAeT/B0aHSY4HspQfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is our model that we want to fit.\n",
    "# It will have a linear background and a number of Gaussian peaks\n",
    "# The parameters for the background and each of the Gaussian peaks\n",
    "# will be held in separate entries in `p`.\n",
    "def n_gauss(x, p):\n",
    "    y = np.zeros_like(x)\n",
    "    \n",
    "    # background parameters\n",
    "    a, b = np.array(p[0])\n",
    "    y += a + b*x\n",
    "    \n",
    "    for i in range(1, len(p)):\n",
    "        g_pars = p[i]\n",
    "        A, loc, sd = np.array(g_pars)\n",
    "        y += gauss(x, [A, loc, sd])\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "5907it [00:05, 1075.69it/s, +500 | bound: 14 | nc: 1 | ncall: 28160 | eff(%): 22.752 | loglstar:   -inf < -4233.325 <    inf | logz: -4244.268 +/-  0.198 | dlogz:  0.001 >  0.509]\n",
      "12557it [00:18, 696.48it/s, +500 | bound: 48 | nc: 1 | ncall: 53297 | eff(%): 24.499 | loglstar:   -inf < -1048.585 <    inf | logz: -1072.989 +/-  0.313 | dlogz:  0.001 >  0.509]\n",
      "17403it [11:47, 24.60it/s, +500 | bound: 2186 | nc: 1 | ncall: 1701762 | eff(%):  1.052 | loglstar:   -inf < -913.031 <    inf | logz: -947.230 +/-  0.372 | dlogz:  0.001 >  0.509]\n",
      "19222it [02:40, 119.67it/s, +500 | bound: 261 | nc: 1 | ncall: 461894 | eff(%):  4.270 | loglstar:   -inf < -911.813 <    inf | logz: -949.640 +/-  0.392 | dlogz:  0.001 >  0.509]\n"
     ]
    }
   ],
   "source": [
    "# the overall parameter set\n",
    "pars = Parameters(name=\"overall_parameters\")\n",
    "\n",
    "# parameters for the background\n",
    "bkg_pars = Parameters(name='bkg')  \n",
    "intercept = Parameter(1, name='intercept', bounds=(0, 200), vary=True)\n",
    "gradient = Parameter(1, name='gradient', bounds=(-20, 250), vary=True)\n",
    "bkg_pars.extend([intercept, gradient])\n",
    "\n",
    "pars.append(bkg_pars)\n",
    "\n",
    "# now go through and add in gaussian peaks and calculate the log-evidence\n",
    "model = Model(pars, n_gauss)\n",
    "logz = []\n",
    "for i in range(4):\n",
    "    if i:\n",
    "        A = Parameter(5, name=f\"A{i}\", bounds=(40, 250), vary=True)\n",
    "        loc = Parameter(5, name=f\"loc{i}\", bounds=(3, 7), vary=True)\n",
    "        sd = Parameter(5, name=f\"sd{i}\", bounds=(0.1, 2), vary=True)\n",
    "        g_pars = Parameters(data=[A, loc, sd], name=f\"gauss{i}\")\n",
    "        \n",
    "        pars.append(g_pars)\n",
    "    \n",
    "    objective = Objective(model, data)\n",
    "    nested_sampler = dynesty.NestedSampler(\n",
    "        objective.logl,\n",
    "        objective.prior_transform,\n",
    "        ndim=len(pars.varying_parameters())\n",
    "    )\n",
    "    nested_sampler.run_nested()\n",
    "    logz.append(nested_sampler.results.logz[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-4244.267669942061, -1072.989448001302, -947.2295781758878, -949.6400447012933]\n"
     ]
    }
   ],
   "source": [
    "print(logz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbe3bd87130>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD4CAYAAAAD6PrjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAghklEQVR4nO3de3RU9bn/8fcj4X4LNwEJGC4BFLxUUuqlKipCtK1Ya1vbrorWX6lVT7W1hfa4Vv392tO1Sjlq66Vaqp5qjy32aBXaihgQLz0VNCBeEAIBBBIQIuEeIAl5fn/MThhhQobMJHsun9das7Ln+/3umefLhP3s2zwxd0dERATgpLADEBGR1KGkICIijZQURESkkZKCiIg0UlIQEZFGOWEHkKi+fft6fn5+2GGIiKSVZcuWfezu/Y5uT/ukkJ+fT0lJSdhhiIikFTPbGKtdp49ERKSRkoKIiDRSUhARkUZKCiIi0khJQUREGikpiIhIIyUFERFplPbfUxCR5DhYe5id1TXs3F/LruoadlbXsrO6hl3VNdTU1R8ZaHZk8dgmDIs1NPbY6CdxvpY1MzbGSx7zXrHjPnZsU/0NHScSS3NjaS6WGOtPPL0/PTq1J5kSSgpm9mXg/wKnAePdvSSq7yfATcBh4HvuviBoLwJ+A7QDHnX3XwbtQ4E5QB9gGfBNd69JJD6RbFRf7+w9WBfZwFfXsKu6lqr9R5ajf+6sbkgANRysrW/yNRs2QvrzK6ll4Q8uTq2kALwPXAP8LrrRzE4HrgPGAKcAC81sZND9EHA5UA68ZWbz3P0DYCZwn7vPMbNHiCSUhxOMTyStHao7fGQDHmMP/shG/cjGfld1DfVNbLxPMsjt0oHcLu3p1aUDg3I7MeaUHvTq0p7cLh3o1aXDkeWukTG5XdrTMafdceNs+GNd0UnDY/V/oi167LHrE2OsR71CS9+LOF8r+g+QxXqt6PVjLcaaX1zvFePfoKnXGtiz87GDE5RQUnD3VRDzEHAKMMfdDwEbzKwMGB/0lbn7+mC9OcAUM1sFXAp8PRjzBJEjECUFyQjuzp6Ddcdu1I/Z0H/yZ3XN4SZfs1P7k4KNdmRDftqAHo0b+4afvbpGNvC9gw1+9045nHRSE+dXEhDrdMtRI5L+ntI6WuuawiBgSdTz8qANYPNR7Z8hcspol7vXxRh/DDObBkwDGDJkSJJCFolPTV09uw4EG+79n9xbbzgVU7X/yPKu6lp2HajlcBO772bQs/ORjXn/Hp0YNaD7J/faY+zBd2p//L13kZZoNimY2UJgQIyuu9x9bvJDap67zwZmAxQWFuosp7SIu7PvUN2x59f3xz4t0/Bz36G6Jl+zQ85J9I7aUx81oHvjnnz0Xn10W4/O7WnXCnvvIi3RbFJw94kteN0KYHDU87ygjSbadwC5ZpYTHC1EjxdpVu3h+sbz6Uefc99ZXcOu/cdeYN19oIbaw03vU/TolEOvrpENeZ9uHRhxcrcjp2Wi9uBzu7SnV9dIW+f27WLeUSOSLlrr9NE84E9mdi+RC80FwJtETiwWBHcaVRC5GP11d3czWwxcS+QOpKlAKEchEi53p7rmcOw7ZPY3fYF178Hj7L23O+kT59qH9+vWeK495umZLu3p2bk9Oe30NR7JPonekvpF4AGgH/APM1vh7pPdfaWZ/QX4AKgDbnX3w8E6twELiNyS+ri7rwxebgYwx8z+A3gbeCyR2CT9vF+xm6mPv8mO/U3fidy9Yw65XY+cihnat+uRjXrX9kedqon87NJBe+8i8TJP8xuPCwsLXX9kJ/25O994dCmrP9rLtIuGxdyDz+3SnvbaexdJCjNb5u6FR7frG82SEl5f+zH/WreDu79wOjdeMDTscESylna7JHT19c6vFqwmr1dnvv4Z3WIsEiYlBQndP97byvsVe/jB5SOb/easiLQuJQUJVe3heu55qZTRA7oz5ewmv68oIm1ESUFC9fRbm/lwRzU/mjxKX+ASSQFKChKa6po6frNoLZ/O78Wlo08OOxwRQUlBQvRf//shlXsPMaNotL5HIJIilBQkFLuqa3jk1XVMPO1kCvN7hx2OiASUFCQUD7+yjn2H6vjR5NFhhyIiUZQUpM1t3X2AP/zrQ774qUGMGtA97HBEJIqSgrS5XxevxR2+P3Fk84NFpE0pKUibKtu+j/9ZtplvnDuEwb27hB2OiBxFSUHa1H8uKKVLhxxuu2RE2KGISAxKCtJmVmzexYsrP+LbFw6jT7eOYYcjIjEoKUibcHdmzl9Nn64duOlCVUEVSVVKCtImXlv7MW+s38G/XTqCbh1VsV0kVSkpSKurr3d+9WKkNPbXVBpbJKUpKUir+/t7W1m5ZQ93TlJpbJFUp6QgreoTpbHPUmlskVSnpCCtas5bm9m4o5rpRaM4SaWxRVKekoK0muqaOu5ftJbx+b25ZJRKY4ukg4SSgpl92cxWmlm9mRVGteeb2QEzWxE8HonqG2dm75lZmZndb0HNZDPrbWbFZrY2+NkrkdgkfI2lsa8YpdLYImki0SOF94FrgNdi9K1z97ODx81R7Q8D3wYKgkdR0P5jYJG7FwCLgueSpnbur+GRV9Yx8bT+jDtVpbFF0kVCScHdV7l7abzjzWwg0MPdl7i7A08CVwfdU4AnguUnotolDT386jr21dQxvWhU2KGIyAlozWsKQ83sbTN71cwuDNoGAeVRY8qDNoD+7r41WP4I6N/UC5vZNDMrMbOSysrKpAcuidmyK1Ia+5pP5TGyv0pji6STZr9aamYLgQExuu5y97lNrLYVGOLuO8xsHPC8mY2JNyh3dzPz4/TPBmYDFBYWNjlOwvHrhWvA4fuXF4QdioicoGaTgrtPPNEXdfdDwKFgeZmZrQNGAhVAXtTQvKANYJuZDXT3rcFppu0n+r4SvrLte3lmWTk3nD+UvF4qjS2Sblrl9JGZ9TOzdsHyMCIXlNcHp4f2mNm5wV1H1wMNRxvzgKnB8tSodkkjsxpKY1+q0tgi6SjRW1K/aGblwHnAP8xsQdB1EfCuma0AngFudveqoO8W4FGgDFgHzA/afwlcbmZrgYnBc0kjb2/ayYKV25h20TB6d+0Qdjgi0gIJlat09+eA52K0Pws828Q6JcDYGO07gMsSiUfC4+7MfHE1fbt14KbPqjS2SLrSN5olKV5dU8mS9VX826UFdFVpbJG0paQgCYuUxi5lcO/OfG28SmOLpDMlBUnY397dwgdb93Dn5aPokKNfKZF0pv/BkpCaunrueWkNowd056qzTgk7HBFJkJKCJOTptzaxqaqaGUWjVRpbJAMoKUiL7T9Ux28WlTF+aG8mjOoXdjgikgRKCtJi//W/G/h43yFmFI1WaWyRDKGkIC2yc38Nv3t1PZef3p9xp+pPX4hkCiUFaZHfvlLG/po6fjRZpbFFMomSgpywil0HeOKNjVxzjkpji2QaJQU5Yb8ubiiNPTLsUEQkyZQU5ISs3baXZ5eX883zTmVQbuewwxGRJFNSkBMya0EpXTvkcOslKo0tkomUFCRuyzft5KUPVBpbJJMpKUhc3J2Z81fTt1tHvqXS2CIZS0lB4vLKmkqWbqjie5eNUGlskQympCDNaiiNPaR3F677tEpji2QyJQVp1t/e3cKqrXu4c9JIlcYWyXD6Hy7H1VAa+7SBPfjCmSqNLZLplBTkuOYEpbGnF41SaWyRLKCkIE3af6iO+xet5TNDezNhpEpji2SDhJKCmc0ys9Vm9q6ZPWdmuVF9PzGzMjMrNbPJUe1FQVuZmf04qn2omS0N2p82M90IH7LH/7mBj/fVMOMKlcYWyRaJHikUA2Pd/UxgDfATADM7HbgOGAMUAb81s3Zm1g54CLgCOB34WjAWYCZwn7uPAHYCNyUYmySgan8Nv3ttPZNO7885Q1QaWyRbJJQU3P0ld68Lni4B8oLlKcAcdz/k7huAMmB88Chz9/XuXgPMAaZYZDf0UuCZYP0ngKsTiU0S89vFZVSrNLZI1knmNYVvAfOD5UHA5qi+8qCtqfY+wK6oBNPQHpOZTTOzEjMrqaysTFL40qBi1wGefGMjXzonjwKVxhbJKs1+NdXMFgIDYnTd5e5zgzF3AXXAU8kNLzZ3nw3MBigsLPS2eM9scl/xGjC4Q6WxRbJOs0nB3Scer9/MbgA+D1zm7g0b6ApgcNSwvKCNJtp3ALlmlhMcLUSPlza0Ztte/rq8nG9dMFSlsUWyUKJ3HxUB04Gr3L06qmsecJ2ZdTSzoUAB8CbwFlAQ3GnUgcjF6HlBMlkMXBusPxWYm0hs0jIqjS2S3RK9pvAg0B0oNrMVZvYIgLuvBP4CfAC8CNzq7oeDo4DbgAXAKuAvwViAGcAPzKyMyDWGxxKMTU7Qso07Kf5gG9+5eBi9VBpbJCslVO4yuH20qb5fAL+I0f4C8EKM9vVE7k6SELg7M19UaWyRbKdvNAsAr5RW8uaGKm6/bARdOqg0tki2UlIQ6usjRwlDenfhqyqNLZLVlBSEee9sYfVHe1UaW0SUFLJdTV099xSXcrpKY4sISgpZ789vbmJz1QGVxhYRQEkhq+0/VMcDL6/l3GG9uVilsUUEJYWs9lhDaewilcYWkQglhSy1Y98hZr+2nslj+vMplcYWkYCSQpb67SvrVBpbRI6hpJCFyndW88c3NnLtuDxGnKzS2CJyhJJCFrqveG2kNPZElcYWkU9SUsgypR/t5a9vl3PD+fmcotLYInIUJYUsM2tBKd065PDdi4eHHYqIpCAlhSyybGMVC1dt4+YJw1UaW0RiUlLIEu7OzPml9OvekRsvyA87HBFJUUoKWWJx6Xbe/LCK711WoNLYItIkJYUsUF/v/OrFUk7t04XrPj24+RVEJGspKWSBue9UBKWxR9G+nT5yEWmathAZrqaunnteWsOYU3rw+TMGhh2OiKQ4JYUM96elGynfeYDpRaNVGltEmqWkkMH2HarjgZfLOG9YHy4q6Bt2OCKSBhJKCmY2y8xWm9m7ZvacmeUG7flmdsDMVgSPR6LWGWdm75lZmZndb0HNZjPrbWbFZrY2+KnSnQl67PUN7Nhfw4wrVBpbROKT6JFCMTDW3c8E1gA/iepb5+5nB4+bo9ofBr4NFASPoqD9x8Aidy8AFgXPpYUipbHXUTRmAGcPzg07HBFJEwklBXd/yd3rgqdLgLzjjTezgUAPd1/i7g48CVwddE8BngiWn4hqlxZ4aPE6DtQe5ocqjS0iJyCZ1xS+BcyPej7UzN42s1fN7MKgbRBQHjWmPGgD6O/uW4Plj4D+Tb2RmU0zsxIzK6msrExS+JmjfGc1/71kI18eN5gRJ3cLOxwRSSPNfrXVzBYCA2J03eXuc4MxdwF1wFNB31ZgiLvvMLNxwPNmNibeoNzdzcyP0z8bmA1QWFjY5LhsdW/xGszgjssLwg5FRNJMs0nB3Scer9/MbgA+D1wWnBLC3Q8Bh4LlZWa2DhgJVPDJU0x5QRvANjMb6O5bg9NM209wLkKkNPZzb1cw7cJhDOyp0tgicmISvfuoCJgOXOXu1VHt/cysXbA8jMgF5fXB6aE9ZnZucNfR9cDcYLV5wNRgeWpUu5yAWQtW061jDt+doNLYInLiEq2M9iDQESgObnlcEtxpdBHwMzOrBeqBm929KljnFuAPQGci1yAarkP8EviLmd0EbAS+kmBsWafkwyoWrtrOjyaPIreLSmOLyIlLKCm4+4gm2p8Fnm2irwQYG6N9B3BZIvFkM3dn5ourVRpbRBKibzRniJdXb+etD3dyu0pji0gClBQywOGgNHZ+ny58VaWxRSQBSgoZYO6KCkq3qTS2iCROW5A0d6juMPcWr2HsoB58TqWxRSRBSgpp7k9LN0VKY09WaWwRSZySQhrbd6iOB18u4/zhfbhQpbFFJAmUFNLYo6+vj5TGLlJpbBFJDiWFNPXxvkP8/rX1XDF2AGepNLaIJImSQpp6aHEZB+vqVRpbRJJKSSENba6q5qklm/jyuDyG91NpbBFJHiWFNHRfQ2nsiSPDDkVEMoySQppZ/dEenltRwQ0X5DOgZ6ewwxGRDKOkkGZmvVhK9445fPdilcYWkeRTUkgjb31YxaLV27l5wnCVxhaRVqGkkCbcnZnzV3Ny947ceP7QsMMRkQylpJAmFq3aTsnGndw+sYDOHdqFHY6IZCglhTRwuN6ZtaCUoX278pVClcYWkdajpJAGnn+7oTT2SJXGFpFWpS1MimsojX3GoJ5cOValsUWkdSkppLinlmyiYtcBpheNUmlsEWl1SgopbO/BWh5cXMYFI/pwYUG/sMMRkSyQcFIws5+b2btmtsLMXjKzU4J2M7P7zaws6D8nap2pZrY2eEyNah9nZu8F69xvWV4P+tHXN1C1v4bpk0eHHYqIZIlkHCnMcvcz3f1s4O/AT4P2K4CC4DENeBjAzHoDdwOfAcYDd5tZr2Cdh4FvR61XlIT40tLH+w7x6OvrufIMlcYWkbaTcFJw9z1RT7sCHixPAZ70iCVArpkNBCYDxe5e5e47gWKgKOjr4e5L3N2BJ4GrE40vXT34cqQ09p2TVBpbRNpOTjJexMx+AVwP7AYuCZoHAZujhpUHbcdrL4/RHuv9phE5+mDIkCGJTyDFbK6q5qmlG/lKoUpji0jbiutIwcwWmtn7MR5TANz9LncfDDwF3NaaAQfvN9vdC929sF+/zLsAe2/xGk4y4/bLVBpbRNpWXEcK7j4xztd7CniByDWDCiD667d5QVsFMOGo9leC9rwY47PKqq17eH5FBd+5aLhKY4tIm0vG3UcFUU+nAKuD5XnA9cFdSOcCu919K7AAmGRmvYILzJOABUHfHjM7N7jr6HpgbqLxpZtZC1QaW0TCk4xrCr80s1FAPbARuDlofwG4EigDqoEbAdy9ysx+DrwVjPuZu1cFy7cAfwA6A/ODR9Z4c0MVL6/ezoyi0fTs0j7scEQkCyWcFNz9S020O3BrE32PA4/HaC8BxiYaUzpyd2a+uJr+PTpyw/n5YYcjIllK32hOEQtXbWfZxp3cftlIlcYWkdAoKaSASGns1Qzr25WvFOY1v4KISCtRUkgBz71dwZpt+7hz0ihyVBpbREKkLVDIDtYe5r6G0thnDAg7HBHJckoKIXtqaaQ09oyi0WR5/T8RSQFKCiHae7CWhxaX8dkRfflsQd+wwxERUVII0+8bSmMXqeidiKQGJYWQVO6NlMb+3BkDOTMvN+xwREQAJYXQPPjyWg7V1XPnJBW9E5HUoaQQgk07qvnTm5v4SuFghqk0toikECWFENxbXEq7k4w7JhY0P1hEpA0pKbSxD7bsYe47W7jxgqH076HS2CKSWpQU2tisBavp3jGHmy9SaWwRST1KCm1o6fodLC6t5JZLRqg0toikJCWFNhJdGnvqeflhhyMiEpOSQhsp/mAbyzft4o6JKo0tIqlLSaENREpjlzKsb1e+PE6lsUUkdSkptIG/Li9n7fZ9/HCySmOLSGrTFqqVNZTGPjOvJ1eMVWlsEUltSgqt7L+XbGTL7oMqjS0iaUFJoRXtCUpjX1jQlwtGqDS2iKS+hJKCmf3czN41sxVm9pKZnRK0TzCz3UH7CjP7adQ6RWZWamZlZvbjqPahZrY0aH/azDokElsqePS19eysrmX65NFhhyIiEpdEjxRmufuZ7n428Hfgp1F9r7v72cHjZwBm1g54CLgCOB34mpmdHoyfCdzn7iOAncBNCcYWqsq9h3j0nxv43JkDOSOvZ9jhiIjEJaGk4O57op52BbyZVcYDZe6+3t1rgDnAFIucbL8UeCYY9wRwdSKxhe2BoDT2DyfpD+iISPpI+JqCmf3CzDYD3+CTRwrnmdk7ZjbfzMYEbYOAzVFjyoO2PsAud687qr2p95xmZiVmVlJZWZnoFJJu045q/rR0E1/99GCG9u0adjgiInFrNimY2UIzez/GYwqAu9/l7oOBp4DbgtWWA6e6+1nAA8DzyQza3We7e6G7F/br1y+ZL50U9xSXktPOuP0ylcYWkfSS09wAd58Y52s9BbwA3B19WsndXzCz35pZX6ACGBy1Tl7QtgPINbOc4GihoT3trNyym7krtnDLhOEqjS0iaSfRu4+id4WnAKuD9gHBdQLMbHzwPjuAt4CC4E6jDsB1wDx3d2AxcG3wWlOBuYnEFpZZC0rp2bk937lYpbFFJP00e6TQjF+a2SigHtgI3By0Xwt818zqgAPAdcGGv87MbgMWAO2Ax919ZbDODGCOmf0H8DbwWIKxtbkl63fwSmklP7liND07qzS2iKQfi2yr01dhYaGXlJSEHQbuzjUP/4utuw7yyo8m0Km9KqGKSOoys2XuXnh0u77RnCQvfbCNtzft4o6JBUoIIpK2lBSSoLE0dr+uXKvS2CKSxpQUkuDZ5eWUbd/HjyapNLaIpDdtwRJ0sPYwvy5ew1l5PSlSaWwRSXNKCglSaWwRySRKCgnYc7CWB4PS2OerNLaIZAAlhQT8/rX17KquZUaRSmOLSGZQUmih7XsP8ujrG/j8mQMZO0ilsUUkMygptNADi8qoPazS2CKSWZQUWmDjjv38+c1Iaex8lcYWkQyipNAC97y0hvbtTlJpbBHJOEoKJ2jllt3Me2cL3/psPierNLaIZBglhRP0qxcjpbGnXaTS2CKSeZQUTsAb63bw6ppKbr1kuEpji0hGUlKIk7sz88XVDOzZievPyw87HBGRVqGkEKcFK7exYrNKY4tIZlNSiEPd4Xr+86VShvfrypfOUWlsEclcSgpx+Ovyikhp7MkqjS0imU1buGYcrD3MfQvXcNbgXCaPUWlsEclsSgrN+OMbG9m6+yAzikapNLaIZDwlhePYc7CWh14p46KR/Th/uEpji0jmS1pSMLM7zczNrG/w3MzsfjMrM7N3zeycqLFTzWxt8Jga1T7OzN4L1rnfQt41n/1qpDT29Mkqeici2SEpScHMBgOTgE1RzVcABcFjGvBwMLY3cDfwGWA8cLeZ9QrWeRj4dtR6RcmIryW27znIY//cwBfOOkWlsUUkayTrSOE+YDrgUW1TgCc9YgmQa2YDgclAsbtXuftOoBgoCvp6uPsSd3fgSeDqJMV3wu5/eS21h+u58/KRYYUgItLmEk4KZjYFqHD3d47qGgRsjnpeHrQdr708Rnus95xmZiVmVlJZWZngDI714cf7mfPmZq4br9LYIpJdcuIZZGYLgVj3Y94F/DuRU0dtxt1nA7MBCgsLvZnhJ+ye4khp7O+pNLaIZJm4koK7T4zVbmZnAEOBd4JrwnnAcjMbD1QAg6OG5wVtFcCEo9pfCdrzYoxvU+9X7OZv72zhtktGcHJ3lcYWkeyS0Okjd3/P3U9293x3zydyyuccd/8ImAdcH9yFdC6w2923AguASWbWK7jAPAlYEPTtMbNzg7uOrgfmJhJfS/xqQSm5Xdoz7eJhbf3WIiKhi+tIoYVeAK4EyoBq4EYAd68ys58DbwXjfubuVcHyLcAfgM7A/ODRZv617mNeW1PJXVeeRo9OKo0tItknqUkhOFpoWHbg1ibGPQ48HqO9BBibzJjiFSmNXcrAnp345nmnhhGCiEjo9I3mwIKVH/HO5l18f+JIlcYWkaylpECkNPasBaWMOLkb15wT8y5YEZGsoKQAPLu8nHWV+/nhJJXGFpHslvVbwIO1h7mveC1nD85l8pj+YYcjIhKqrE8KT77xIR/tOciMotEqjS0iWS+rk8LuA7U8tHgdF4/sx3nD+4QdjohI6LI6Kcx+bR27D9QyvUilsUVEIIuTQkNp7KvOOoUxp6g0togIZHFS+M2itdQddu6cpNLYIiINsjYpDO7dhWkXDePUPiqNLSLSoDVrH6W0my8eHnYIIiIpJ2uPFERE5FhKCiIi0khJQUREGikpiIhIIyUFERFppKQgIiKNlBRERKSRkoKIiDSyyJ9STl9mVglsbOHqfYGPkxhOmDJlLpkyD9BcUlWmzCXReZzq7v2Obkz7pJAIMytx98Kw40iGTJlLpswDNJdUlSlzaa156PSRiIg0UlIQEZFG2Z4UZocdQBJlylwyZR6guaSqTJlLq8wjq68piIjIJ2X7kYKIiERRUhARkUZZkRTMrMjMSs2szMx+HKO/o5k9HfQvNbP8EMJsVhzzuMHMKs1sRfD4P2HEGQ8ze9zMtpvZ+030m5ndH8z1XTM7p61jjEcc85hgZrujPpOftnWM8TKzwWa22Mw+MLOVZnZ7jDEp/7nEOY+0+FzMrJOZvWlm7wRz+X8xxiR3++XuGf0A2gHrgGFAB+Ad4PSjxtwCPBIsXwc8HXbcLZzHDcCDYcca53wuAs4B3m+i/0pgPmDAucDSsGNu4TwmAH8PO8445zIQOCdY7g6sifE7lvKfS5zzSIvPJfh37hYstweWAuceNSap269sOFIYD5S5+3p3rwHmAFOOGjMFeCJYfga4zMysDWOMRzzzSBvu/hpQdZwhU4AnPWIJkGtmA9smuvjFMY+04e5b3X15sLwXWAUMOmpYyn8ucc4jLQT/zvuCp+2Dx9F3ByV1+5UNSWEQsDnqeTnH/oI0jnH3OmA30KdNootfPPMA+FJwWP+MmQ1um9BaRbzzTQfnBYf/881sTNjBxCM4BfEpInum0dLqcznOPCBNPhcza2dmK4DtQLG7N/mZJGP7lQ1JIZv8Dch39zOBYo7sPUh4lhOpMXMW8ADwfLjhNM/MugHPAne4+56w42mpZuaRNp+Lux9297OBPGC8mY1tzffLhqRQAUTvMecFbTHHmFkO0BPY0SbRxa/Zebj7Dnc/FDx9FBjXRrG1hng+t5Tn7nsaDv/d/QWgvZn1DTmsJplZeyIb0qfc/a8xhqTF59LcPNLtcwFw913AYqDoqK6kbr+yISm8BRSY2VAz60DkQsy8o8bMA6YGy9cCL3tw1SaFNDuPo87tXkXkXGq6mgdcH9ztci6w2923hh3UiTKzAQ3nd81sPJH/c6m2wwFE7iwCHgNWufu9TQxL+c8lnnmky+diZv3MLDdY7gxcDqw+alhSt185LV0xXbh7nZndBiwgcgfP4+6+0sx+BpS4+zwiv0B/NLMyIhcNrwsv4tjinMf3zOwqoI7IPG4ILeBmmNmfidwB0tfMyoG7iVxEw90fAV4gcqdLGVAN3BhOpMcXxzyuBb5rZnXAAeC6FNzhaHAB8E3gveAcNsC/A0MgrT6XeOaRLp/LQOAJM2tHJHH9xd3/3prbL5W5EBGRRtlw+khEROKkpCAiIo2UFEREpJGSgoiINFJSEBGRRkoKIiLSSElBREQa/X/ZvcSc1L8rPgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(logz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The log-evidence points to the use of 2 Gaussians to fit the data. There is a sufficient increase in evidence over 1 Gaussian. However, 3 Gaussians is not justified, the logz term does not increase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}