model selection with lmfit

lmfit_model_selection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model selection using lmfit and emcee "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`lmfit` can be used to obtain the posterior probability distribution of parameters, given a set of experimental data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import lmfit\n",
    "import matplotlib.pyplot as plt\n",
    "import triangle\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gauss(x, a_max, loc, sd):\n",
    "    return a_max * np.exp(-((x - loc) / sd)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.linspace(3, 7, 250)\n",
    "np.random.seed(0)\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 * np.random.randn(np.size(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Container object of 3 artists>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl4VdX1/t8VQmROCJAQCPOMMgpI0WIqCk4FtUWxTnVs\na+2gtXVshTpUbbWTVavUVv22VutI/YmKSgRFQBAQCCggQwYykAQIhJn1++O9p3fInZLcOevzPPe5\n5+wz3HVPbt6zztprry2qCsMwDCN1SYu3AYZhGEZ0MaE3DMNIcUzoDcMwUhwTesMwjBTHhN4wDCPF\nMaE3DMNIcYIKvYi0EZFlIrJaRIpE5Deu9mwRWSAiX4rIuyKS5XHMHSKySUQ2isjUaH8BwzAMIzgS\nKo9eRNqpar2IpAP4CMCtAKYD2KWqD4vIbQA6q+rtIjIcwL8AjAfQE8B7AAar6vGofgvDMAwjICFD\nN6pa71rMANAKQC0o9M+62p8FcIFreQaAF1T1iKpuA7AZwIRIGmwYhmE0jpBCLyJpIrIaQAWAhaq6\nHkCuqla4dqkAkOta7gGgxOPwEtCzNwzDMOJEeqgdXGGX0SKSCeAdEfmGz3YVkWDxH6uxYBiGEUdC\nCr2Dqu4Rkf8H4GQAFSLSXVXLRSQPQKVrt1IAvTwOy3e1eRHixmAYhmEEQFWlsceEyrrp6mTUiEhb\nAGcBWAVgHoCrXLtdBeB11/I8ALNEJENE+gEYBGB5AGMT/nXPPffE3YZUsTMZbDQ7zc5EfzWVUB59\nHoBnRSQNvCk8r6rvi8gqAC+JyLUAtgG42CXeRSLyEoAiAEcB3KjNsc4wDMNoNkGFXlXXAhjrp70G\nwJkBjnkAwAMRsc4wDMNoNjYyNggFBQXxNiEsksHOZLARMDsjjdmZGIQcMBWVDxWxiI5hGEYjERFo\npDtjDcMwjOTHhN4wDCPFMaE3DMNIcUzoDcMwUhwTesMwjBTHhN4wDCPFMaE3jABMnw7U1sbbCsNo\nPib0hhGAhQuB7dvjbYVhNB8TesPww/HjwL59QEVF6H0NI9ExoTcMP+zfz/fy8vjaYRiRwITeMPxQ\nV8d38+iNVMCE3jD8YEJvpBIm9EaLZ+1a4PbbvdscobfQjZEKmNAbLZ7PPgMWLfJuq6sDRMyjN1ID\nE3qjxVNaCuzd691WVwfk55tHb6QGJvRGi6esDNizx7utrg4YNMg8eiM1MKE3Wjylpf6Fvm9fYPdu\n4MiRuJhlGBHDhN5o8ZSWUtiPH3e31dUBWVlAly5AVVX8bDOMSGBCb7R4ysr47mTaOMsdOwK5uRa+\nMZIfE3qjRXPsGFBZSUH3DN+Y0BuphAm90aKprASysxmi8cy82bePQt+9u2XeGMmPCb3RoqiqAg4d\ncq+XlgI9egCZmf49+r59gQ0bYm6mYUQUE3qjRXHFFcBTT7nXS0uBnj2BTp38C/355wNvvBF7Ow0j\nkpjQGy2GI0eAjz4C3nrL3VZW5vboPUM3jtCPG8cwjnn1RjJjQm+0GFatAnJyKPYHDrDN8egDhW7S\n0oALLgBeey0+NhtGJDChN1oMixYB550HjB4NfPgh2xyPvlMn/x49wCkFPZ8CDCPZCCr0ItJLRBaK\nyHoRWSciP3a1zxaREhFZ5Xqd43HMHSKySUQ2isjUaH8Bw/Dk8ceBrVv9b/vwQ2DyZODss4H589m2\naxfQrVtgjx4ABg4ESkq4vHUrUFwcPfsNIxqE8uiPALhZVU8EMBHAD0VkGAAF8KiqjnG95gOAiAwH\ncAmA4QDOBvC4iNhTgxEzXnyRIRpfVBmymTwZGD8eKCpi+65dQNeu/oW+QwcuO7n0qsAf/wg8+WT0\nv4dhRJKgIqyq5aq62rW8D8AGAD1dm8XPITMAvKCqR1R1G4DNACZEzlzDCE59PVBb616/4w7g00+Z\nL5+eTtHu1s1d1mDXLubQe4Zujh5lCmb79lxv356x+n37gJ07rSSCkXykh7ujiPQFMAbAUgCnAviR\niFwJYAWAn6nqbgA9XNsdSuC+MRhG1KmvZyEygGL92GOMwasCffqwvVs3Cj8AVFd7e/S//jUHULVv\nz3r0Do5XX1YGHD4c2+9kGM0lrLCKiHQA8DKAn7g8+ycA9AMwGsBOAI8EOVyba6RhhIun0C9cSC/8\nq6+A7du9hX7XLnrue/YAnTu7hf7NN4GHHnLH5x0cod+5032TMIxkIaRHLyKtAbwC4P9U9XUAUNVK\nj+1zAfzXtVoKoJfH4fmutgbMnj37f8sFBQUoKChonOWG4QdPoZ83D5gwgR2oPXu6hf6EE4C2bdme\nlQW0asXQza5dwJYtFP3MTO/z5uS4hT7Nep2MGFFYWIjCwsJmnyeo0IuIAPgbgCJV/YNHe56q7nSt\nXghgrWt5HoB/icijYMhmEIDl/s7tKfSGESkcoVel0P/2t8BvfgP06sXsGYecHA6C6tKF65mZwPr1\nwIknAt/8JrBggfd5c3OBzZuBgwctRm/EDl8neM6cOU06TyiP/lQAlwP4XEScXIY7AVwqIqPBsMxW\nAN8DAFUtEpGXABQBOArgRlW10I0RE1TdQl9dDezfD5xzDnDDDaxZM2WKe99u3Zh507Ur1zMzefy4\nccCNN3J/T3JzgTVrgH79GAY6fBjIyIjVNzOM5hEq6+YjVU1T1dGeqZSqeqWqjlTVUap6gapWeBzz\ngKoOVNWhqvpO9L+CYZDDhzl5yO7dDLN0787QTOvWwIoV7tAN4PboHaHv1Inv48dzrtgbbvA+d04O\n0zZ79uQxu3a5t+3fDzz6aHS/m2E0B4s2GilDfT3fHaHPzeV6//6MrXsKfbdu3qGbNm3ooY8b5//c\nubnAF18AeXneWTsAsG4dcOut7glMDCPRMKE3Uob6eqZE+gp9v34c/NS5s3tfX48eAH71K2DUKP/n\nzsnh00JeHperqpiff+wYR82qAq+8Er3vZhjNwYTeSGpqaymyAIU+J4dtvh59nz7eefHdujH10lPo\n77qLGTn+cM7l6dGffTawfDkLo/XvD/znP5H/foYRCUzojaTm3HPdGTKO0B88SPHNyWG7I/SedOvG\nd0+hD4Yj9D168NiPPwZqapijX1ICXHUVsHYtQ0SGkWiY0BtJy/HjFNdly7heX88RrVlZwJdfusV5\n5kymWHri3AScGH0osrJYQsEJ3cybx3ZH6AcMAMaOZbzeMBINE3ojadmxgxkvK1Zwvb4eaNeOovzF\nF26hz84GRo70PraxHn1aGrN4evbksaWlPKcj9Pn5zNW3ypZGImJCbyQtRUUcBLVyJdc9hX7LFrfQ\n+8Px6MMVegAoLASGDHEfe9llJvRGcmBCbyQtRUWM0R84AJSXewv90aPBhd4R+HBDNwDDMyL06Nu0\nAS68kKNlnVmq8vNN6I3ExITeSFqKiliy4OST6dV7Cj3g9rz9kZEB/PWvDOs0lpNOAu67j528VVUs\ngNamjXn0RuJiQm8kLUVFwPDhFPoVK7yFvkMHLgfjhhuaVqAsKwv42c9YDK1PH3rygAm9kbiEXY/e\nMBKFVas4LWBRETBsGMMnCxdS2J1XsLBNJBkwgNk4gFvoVb1z9g0j3pjQG0nH/Pmczq9/f8bYc3I4\ngMnx6DMyYif0/fszzRNwlzbes4di7zkS1zDiiYVujKSjvh64/npg9WquO1MDeoZuYiX0M2awrDFA\nL75XL+C114DRo2Pz+YYRDubRG0nH/v3enaieHn1+PssR+5YZjhbTpnmv9+oF3Hsvs4AChXBUGXqy\nuXaMWGEevZF07N/vnrgbcHv0+/fTox86FDjvvPjY1qsX69UfP860T39UVQFnnukO+RhGtDGhN5IO\nJ0Tj0K4ds2cqK0Nn2kSb3r1Z7Cwnh7Vw/FFezqqXdXWxtc1ouZjQG0mHr0cPUFi3bYu/0N94I/DU\nU+yIra31v49T+CzQjcAwIo0JvZF0+BP6bt0YMom30HftylGywYS+vJzv1dWxs8to2ZjQG0mHE4v3\nJCeH9eV9bwDxIjs7sMduHr0Ra0zojaTDKUfsiVONMt4evUM4Hr0JvRErTOiNpCNQjB5IHKEP5tGX\nl3NwlQm9EStM6I2kw1/oJpk8+p07WYzNYvRGrDChN5IOf6GbRPbo6+qA005zz21bXk6hN4/eiBUm\n9EbC8/bbwOuvu9cDZd0AiSP0nh59WRnnmN2+nes7d7Lqpgm9EStM6I2EZ8kS4P33uXzkCEeUZmR4\n79OtG8sNnHBC7O3zh6fQ79rF9xUreJM6fJjF0Cx0Y8QKE3oj4amvd4umE5/3rSGTk+O/PV54hm4c\noV+5kmGbvDxW3TSP3ogVJvRGwuMp9P7i8wAHKT35ZGztCoavR5+bS4++vJyTjAfLyjGMSGNCbyQ8\nBw64RdFffB7gbE+XXx5bu4Lh69FPnUqPvqyMHr0JvRFLTOiNhMdf6CbRycriBCTHj1PoTzqJdt90\nE3DOOW6P38nEMYxoElToRaSXiCwUkfUisk5EfuxqzxaRBSLypYi8KyJZHsfcISKbRGSjiEyN9hcw\nUp/6erf3Gyh0k2ikp1PY9+6l0HftCvz0pyx4dv317Exu25bbDSPahJp45AiAm1V1tYh0ALBSRBYA\nuBrAAlV9WERuA3A7gNtFZDiASwAMB9ATwHsiMlhVrfK20WQOHHB7v4FCN4lIdjbtdoT+mmsabq+p\ncU9BaBjRIqhHr6rlqrratbwPwAZQwKcDeNa127MALnAtzwDwgqoeUdVtADYDmBAFu40WRH09cPQo\nRT5ZQjeAOzzjCL0v2dmWYmnEhrBj9CLSF8AYAMsA5KpqhWtTBQBnhs4eAEo8DisBbwyG0WTq6/le\nU5NcHn1+PgdJmdAb8SasOWNdYZtXAPxEVevEI1lZVVVEgnUp+d02e/bs/y0XFBSgwCbQNAJw4AAH\nQtXWJk+MHmAH7Lp1FHN/Qt+5M7B7d+ztMpKHwsJCFBYWNvs8IYVeRFqDIv+8qjoD0StEpLuqlotI\nHoBKV3spgF4eh+e72hrgKfSGEYz6eqBHDwp9MoVuRowAXn6ZtW6yshpuz8xkZo5hBMLXCZ4zZ06T\nzhMq60YA/A1Akar+wWPTPABXuZavAvC6R/ssEckQkX4ABgFY3iTLDMNFfT0HRCVb6Oakk4DFi+m5\np/n5TzOhN2JFKI/+VACXA/hcRFa52u4A8CCAl0TkWgDbAFwMAKpaJCIvASgCcBTAjaqWKWw0jwMH\nKPTJFroZOpShmUGD/G83oTdiRVChV9WPENjrPzPAMQ8AeKCZdhkGAA44OniQoRvHo/cX705EMjKA\nwYNZ18YfWVnApk2xtclomdjIWCOhOXiQHbFOTnoyxegBxukD3ZgyM92dsRUV/vcxjEhgQm8kNAcO\nUNg9hT5ZQjcAMHKku1a+L07o5tAhoF8/4Nix2NpmtBzCSq80jHhRX0+h79yZoZtDh5JL6G+80T0O\nwBdH6KuqeEOrrwc6doytfUbLwITeSFiczte2bd0evWpyhW6ysvynVgLeQg8wDdOE3ogGFroxEpKN\nG4HJk92hm86dWeJ3587k8uiD4Qh9pWsUyr598bXHSF1M6I2EZNcuCrsTusnLY4bKsGHA6NHxti4y\nOKWMHY/ehN6IFha6MRKSujqGavbuZeimVy92xKan0C+2Uydvj76uLr72GKmLefRGQlJXx3h8aak7\nJp9KIg8ArVsDbdoA27Zx3Tx6I1qY0BsJQW0tMHeue93xbouLk6vztbFkZgKbN3PZPHojWpjQGwnB\nmjXAHzyqKTkzL23fztBNquIIfXa2efRG9DChNxKCXbu8p9VzvNsdO1Lbo8/KArZu5YApX49+zRp3\n/N4wmoMJvZEQVFd7C11LCt0cPQr079/Qo//1r4E33oiPXUZqYUJvJASO0Du1TuvqKII7dqR+6Aag\nR+8r9KWlFrc3IoMJvZEQVFez1svBg1yvq6OXe+hQ6nv0rVpx2kFfUS8rM6E3IoMJvZEQOHOnOsK2\ndy+9XCD1hb5rV+bUe3r0x49zFLAJvREJTOiNhMBX6Ovq3EKf6qGbbt1Y48ZT1HftYuzes4PaMJqK\nCb2REPgT+v79uZzKHn1WFoW+Qwdvj77UNdOyefRGJDChNxKC6mqGMBwPtqUIfffuLO/g69GXlXGe\nWRN6IxKY0BsJQXU10LdvywvdXHQR8NRTDT36sjL/ufWG0RRM6I24c/w4p9Tr08e7M7ZHD867msoe\nfVoap0rs2LFh6GbIEBN6IzKY0BtxZ/duCl1WFoXt2DHWoW/fnuGcVBZ6hw4dGoZuhg41oTcigwm9\nEXd27QK6dHHHqffto/ClpQHnnMMc81TH16N3hN6yboxIYEJvxIXaWuDVV7lcXe0t9J5T6s2dy0lH\nUp2MDIawDh/melmZhW6MyGFCb8SF5cuBu+/mcjChbymIuL//ihUs/TBwIEcGHz0ab+uMZMeE3ogL\ntbWccEPVLfSdOjFUsXdvyxN6gOGqv/wFmDEDmDYNyM1tmI1jGE3BhN6ICzU17HCtrDSP3qFjR+Cf\n/wQefJDvrVo1zK83jKZgQm/Ehdpavm/b5l/oO3WKq3lxoUMHToA+bZq7rVMnE3qj+ZjQG3Ghpobv\nW7cCGzcCAwaYR9+xIzB+PJCT491mQm80l5BCLyLPiEiFiKz1aJstIiUissr1Osdj2x0isklENorI\n1GgZbiQ3tbWcPm/bNmDlSgpcSxf6Dh2Ac8/1buvY0VIsjeaTHsY+fwfwZwDPebQpgEdV9VHPHUVk\nOIBLAAwH0BPAeyIyWFWPR8heI0WoqQHGjgU+/ZQDpgYM4NypdXUttzP25z931/dxMI/eiAQhPXpV\nXQyg1s8m8dM2A8ALqnpEVbcB2AxgQrMsNFKS2loK/fz5wMknc3CUI2olJcw4aWmcemrDMQONFfpD\nh4AFCyJrl5H8NCdG/yMRWSMifxORLFdbDwAlHvuUgJ69YXjhePQHDgDjxrHNCVO88w4wZUp87UsU\nAnXGDhniLu3syerVwFVXRd8uI7kIJ3TjjycA/Nq1fC+ARwBcG2Bf9dc4e/bs/y0XFBSgoKCgiaYY\nyUhtLTBmDJcdoW/fnh5pfT0wcmT8bEsk/Hn0hw8DX37J1NQuXby31dZyZqo9e9zz0frjvPOAe+7h\n3+C554BrrwWWLmXn+KWXRv57GE2jsLAQhYWFzT5Pk4ReVSudZRGZC+C/rtVSAL08ds13tTXAU+iN\nlkdNDdCzJ0d/nnIK20TYIXn22Vw2/At9eTnfa/0EVJ1spg0bgIkTA593+3Zg8WKWXbjuOnYCz53L\n403oEwdfJ3jOnDlNOk+TQjci4hlJvBCAk5EzD8AsEckQkX4ABgFY3iTLjJTl4EFWqGzXjqmVvXu7\nt3XsyEJmBvGXdbNzJ9937264vyP+GzYEP+/evewIX7qU6++8A7z7LmvsGKlHSI9eRF4AcDqAriJS\nDOAeAAUiMhoMy2wF8D0AUNUiEXkJQBGAowBuVFW/oRuj5VJbC3TuTK+9VSvvbXffbULviT+P3hHj\nQEKfkRGe0K9YwRIUBQXAH//ImL89SaUmIYVeVf09yD0TZP8HADzQHKOM1Kamhjn0/vj+92NrS6Lj\nKfRVVXwKCiX048YFF3pVd93/998HXnmFYn/55cBLL3G7CPD66/TyH3884l/LiDE2MtaIOY5Hb4Sm\nXz9g3TqK7xVXAE88QaFv1cp/jL62Fpg0KbjQ798PtGnDtFYAmDyZnbIzZrCPxMnm+eILvozkx4Te\niDnBPHrDmzFjKPJvvsn8+M8/Z4x+wIDgHn1JCftCHH74Q6C4mMt79zJtc9w4dtiK0HO/6CJO3+g8\nMezc6e74NZIbE3oj5phHHz4iwGWXMf3xxBMp9GVlwPDhgYU+J4cjbDdudLfPnw8UFXHZEfrrrgPu\nuINtXbty0Fpengl9KmJCb8QMVY7+XLTIPPrGcNlljM/ffz9DKTt2BBb6mhreREeP5uAph+pqoKKC\ny47QDx3Kv4cnPXq4s3p27uT5nFmvjOTFhN6IGdu2AZ98AjzzjHn0jWHQIIZuzj0X6NWL8ffhwwPH\n6Dt3Zshn1Sq2HTlCcfcVen/4evQAB2YZyY0JvREzPvuMIzLPOss7d94IzXnnsQN2xAggPZ3iHyh0\n4yv0zg0hHKH39egHDrTwTSpgQm/EjJUrmenx9ttWj6WpjBhBrzs7u6HQHznCDtiOHSn0a9Zw5KuT\nRROu0JeVMf1S1YQ+VTChN2KGI/RpaXwZjWfkSIpxVlbD0E1tLdtFWAMnM5O1axoj9E7oprycy3l5\nJvSpQFOLmhlG2Hz8MUVn5UpWrDSazvnnsxM1M5OFy5zBTYB7MheH0aMZvsnIYGzfEexgUzU6oZud\nOyny3bu7bxBG8mJ+lRF1HnqIg3jS0ykkRtPJyGBHbOvWQNu2wL597m2+aaujRzMds7qax4Qbutm9\nG1i71i305tEnPyb0RtTZsgW4+WbgmmuslkokycryjtM7qZUOAwYAX31FoR8yhNuPHQsu9BkZrDX0\n+OMU+txcE/pUwITeiCrHjzNOfOutwANWASmi+MbpfT36vn2Z0lpdzUFUWVnArl3BhR4AZs7k4Crz\n6FMHE3ojquzcySyQljgHbLTp3Nnt0d93H/tAPIW+Xz/eZGtq2DnrxNtDCf0553ASGIvRpw4m9EZU\nOHqUHuaWLQwhGJHHCd1UVlLof/97b6Hv2ZMefFkZhd4Jw4QS+nbtgDvvBMaPp9CXlXnXxPc3UMtI\nbEzojajw1lvMEPnqK9ZdMSKPI/QLFtALf+EFDqxyaNUKyM9n5k12NoU+HI8eoNAPH879rr6aov/a\na8C99/I8/uarNRIXS680okJZGcsdjBhhHn20yMlhOYSyMmDqVGDWrIb79O0LfPCB26MPV+gdRIA/\n/xn4z3/4rspc/uXLbYKYZMI8eiMqVFZSFJ5/3oQ+Wtx0E/D008B//wtMm+Z/n379+N6lC8tOfPll\n44TeYeZM3jAWLuRnOVMQGsmBCb0RFaqqOAy/vt6EPlr07QvcdRfQrVvg8FjfvnzPzgYuuAB49VUO\ntGqs0HsycaIJfbJhQm9EhcpK4Lvf5aO/CX30+OlPg4tuv37ACSewg7VvX9a0F2FbUznlFGDZMqbO\nGsmBCb0RFSorKSqrVzNzw4gOIsFLPvfty7CNM1Dtiiua580D7Bvo0sWmGUwmrDPWiApVVRSEESPi\nbUnLZswYZso4XHwx8+qby/jxzNsfNqz55zKij3n0RkRxZiOqrGTs2Igv7dqx9IRDp07AL37R/PMO\nGMBRt0ZyYEJvRIxjx1gUq66Og2q6do23RUa06N2bUxoayYEJvRExdu/mQJqPP6bnmG6BwZTFhD65\nMKE3IoYzWnLhQsbnjdTFhD65MKE3IobTybdwocXnU51evSj0qvG2xAgHE3ojYtTUMFyzcqV59KlO\np06c/MSzwNnGjZHp6DUijwm90SRU2fnqSU0NMG4cB9KYR5/6+IZv1q8H3n8/fvYYgQkp9CLyjIhU\niMhaj7ZsEVkgIl+KyLsikuWx7Q4R2SQiG0VkarQMN+LLc88BZ57p/eheU8M5Ydu1M4++JeAIvTNC\ntro6Mjn6RuQJx6P/O4CzfdpuB7BAVQcDeN+1DhEZDuASAMNdxzwuIvbUkIK8/TaH3j/3nLutupqe\n/PDh5tG3BHr3ZimE/v2ZUltdbbXqE5WQIqyqiwH4/vmmA3jWtfwsgAtcyzMAvKCqR1R1G4DNACZE\nxlQjUVAFFi0C/vEP4Lbb3B5dTQ2LZ33nOxw5aaQ2vXsDv/sdsH07JzSpqWHBNN+QnhF/mupt56qq\nM8FYBYBc13IPACUe+5UA6NnEzzASlC1bWDvl4ouBtm2BzZvZ7gj9zTez8JWR2vTuzb//0KEcCe2k\n13pOWG4kBs0e0qKqKiLBkqz8bps9e/b/lgsKClBQUNBcU4wYsWgRMHkyxX7sWOCzz4DBg91Cb7QM\npk8HBg4E7r+fE5o4Ql9by6JnRvMpLCxEYWFhs8/TVKGvEJHuqlouInkAKl3tpQB6eeyX72prgKfQ\nG8nFokXA6adz2RH6WbNM6Fsa7dszRJeb6+3RW4ds5PB1gufMmdOk8zQ1dDMPwFWu5asAvO7RPktE\nMkSkH4BBAJY38TOMBGXrVnrwAKsjrlrF5epqE/qWSE6OW+h79vTfIXv8uNWvjyfhpFe+AGAJgCEi\nUiwiVwN4EMBZIvIlgDNc61DVIgAvASgCMB/Ajao2di7V2LOHE1MDbo9elZ6cPbK3PDyFfuBA/x79\n7NnAY4/F3DTDRcjQjapeGmDTmQH2fwDAA80xykhs9uwBMjO53L07Zyvato1zkWZlBT3USEFychjO\nq61l+WJ/Hv2mTcChQ7G3zSBWX9BoNJ5CDwAnnwzMnw907Ai0ahU/u4z4kJtLIW/blsv+hL6srPkz\nWxlNxwYzGY1ClZ675z/tRRcBjz9u8fmWSk4OpxXs0oW/AX+hm7IyG0wVT0zojUaxfz9DNa1bu9u+\n9S3m1pvQt0xycoCDByn0nTt7C/qBA3QOSkstvz6emNAbjWL3bu+wDUDvfsYM64htqWRnA2lpfPf0\n6AsLgW98g6G+AwdM6OOJCb3RKHzj8w4//jEwbVrs7THiT1oaaxv5evTvvgt8/jlQUtKwpHEiMXs2\nQ4+pjHXGGo0ikNBPmsSX0TLJzXXH6B1BLyykJ//xx8CQIayHEy6LF3O07be/HRVzvVi/nr/rVMY8\neqNRBBJ6o2WTk+P26GtqgH376M2ffjrw3nusaLp7d/gzUr3yCvDnP0fXZofi4sbdhBrLxo3Af/4T\nvfOHgwm9EZJ//pNliVVN6A3/5OTQm3dCNx9/zLTbk0/mZCT9+jF8s39/eOdbt45lsA8ciK7dAIV+\n587onX/uXO9y3uFy003Az3/OjKXmYkJvhOSRR4Crr2Ys04Te8MeddwKXXML6N0eOAA8/zI7YE0+k\n8PfowZtAuB2y69YBeXnAJ59E1+7Dhyny0fTo33qraaGhf/wDWL0aeOqp5ttgQm+EpKIC+NWvgE8/\nNaE3/HPiiRwlLQJMmQKcdhpw660M2QAU+qys8Dpkq6o4inbWLMb5o0lZGdChQ/Q8+q1bGbppbMbR\noUO8Cc1gAwy7AAAXXklEQVScySeO5mJCbwTl+HH+440fzwkmTOiNUMyfD8yZQwEdNoxtPXu6PfpX\nXw0eklm3DjjpJD4RRFvoi4v5WYcOAfX1kT//W2+xpHdjPXqnQGCvXib0RgzYvZuP44MGuYXe6tkY\n4ZKZCUyYwOkGHY/+2muBf/878DGO0J9yirtgXrTYsYMTqHTvzifXSPPBB5ygx1Po33iDN8JgVFez\nc9uE3ogJFRVMncvMBNLT+ShqHr3RGJYt428oK4ulEnbvBv7+d/f2ffvcy3v3AmvXUuizsoCMDGDX\nLm5T5dNCJCkudgt9NMI327ezwmtdHZ+Oi4qAK68EXnwx+HG+Qt/cm50JvRGUigpmVABAnz78JzSh\nN5pC587sXJ0wgYK/aRPbTzsNWLKEseysLODpp4HRo7mtXz9WRgVYZuPcc4GPPgr9WRs3Ar/9LZf3\n7mUHsT+KiymmeXmR65AtK2N/FsDBYn368Km4rg647z7g9tvpMB0+HPgcjtB36sR+j+bm+ZvQG0Gp\nrKQ3BvAHW1xsQm80jawsCv3IkexofeklTiS+YQOwYAFH0l5zDfPwncF3ffu6hX79elbIvPfe0J/1\n6qvAL3/JUNH55wN//av//Ryh797dv9D/4Q+NH9H7zDMU9MOH+V1ycvjdd+/mTWDiRD5FOHMt+8MR\nepHIhG9M6I2gOKEbgEIPmNAbTaNzZ4rp0KH06teuZbGzw4cZy/7gA2bsdO5MgQPo0W/dyuWiIuD6\n6+mtn3ce8K9/Bf6sZcvoRd9yC0fZOrOg+eIZo/cN3Rw5wrTR1au9255/Hli4MPBnL11KES8r45NC\nq1b8n9mzh45Tt27MRioqch9z9Kj3ORyhB0zojShSWcn4om/oBjChN5qG04k/bBjFfsMGhmNOPhlY\nuZIZNmec4X2Mr0c/ZgwnORk3LrDQq1LoH3iAuegzZ3KUri9ffcXQyqBB/kM3a9YwO6ikhOd86ilO\noXn//Tx3oM9eupQ3px07gPx893d3hD4nx1voH36Y6acvvUQ7Dx5sKPQlJYGvaziY0Bt+mT0buO02\n79BN7958N6E3moKv0G/axFj9yJEU8Px892/NoW9fb49++HA6HOef39ADP34c+MUv6PEDwHXXAd/7\nHvD73/PYY8e897/nHuBHP+KEOf48emewVmkpj//VrzhKfNmywKN2N2/mk0RmJuP0jtBnZlK89+yh\ngDtC/957tO8vf+HAxClTgCeeMI/eiAFHjwIvv8wOMn+hG5spyGgKnTsDbdrQYWjfniGM99/nPLPn\nnMOXL05n7LFjFHAnLz8vr6Ewb9rEDthZs5ia2aoV8OSTzOHPyaEH77B2LfsFbrmF6z16NBTTJUt4\nAyopoYCPG8e+g8xMYMQIlnm48krvEMzSpYzBDxgAfPiht9Bv2cKbXatWFPolS4DLL+eTycyZvIHc\ndx+fXDznX87PN6E3osCHH9KLOHCAecxO6KZfP/5QPScdMYxw6dOHg6Cc6SaHDgXeeYdCf/vtwEMP\n+T9m+3Z69d260fsG6HxUVfEGMHcuc++XL2cRta1b2QfgyYgR3uGbu+/mE6tzvhEjKMR799KbLi2l\nEM+cyeWvvuJYAIczz+TI3+efZ9/CsWPAo49y3RH6RYt4kwH4f7Npk/t/acgQ3kB+9jNeE4ehQ3lD\nM4/eaDYPPcQfdCBefJFe0aRJ/DE6Hn1ODr0Nw2gK/ftzpKjDsGFMORwwgDXt0/yoUYcOfL31Fsss\nOLRuzSeEqioK/dy59IinT+eApKuv9j7PiBH04gGGZFatAn7wA/f2E06gx75kCXDzzbR1/36KcEmJ\nf6Ffs4azq61Zw5vIb3/LENA3v8nvtGePt0fvKfTt2jG0c+ut3nb6E/qRI+lwNWdydatH38I4fJix\nxsmTga99reH2fftYInb1asY8583zjpv26BE7W43UxgnDDBgQfL9+/fibffpp73YnfLNlC99zcoBL\nLwVOPbXhOUaOdHfePvggcNddDCN5cvrpfMKYN4/iXV/PczpCf9ZZ7n0nTQJef51ifMstvEGcdx5v\nOJ7fyVfoJ050n+PkkxvamZPDp4OtW91C3707b1QLFgS/TsEwj76FsW4dxb6khP8c3/mO9/bnnwcK\nCvi4OGkS85bbt4+LqUaKM3Qo0LVr6JIazz/POP3Mmd7teXnszK2vp4e/Zg1j6v446ywOtFq6lOmW\nl1/ecJ/Jk9kROm4cbRs7liJbXc3P8fTo09M5febIkXzKXbzYe+IdX6HPymL4xfHoAyHCzz50yHsO\n5m9/u3k17U3oWxgrV/K9pISPmy+84E4rUwUee4yZCAA7tP70J3dOs2FEklNOCTyQyZPBg/3fDPLy\n2CE6YABw4YUsm9Cunf9zZGdzMNb06aw94895mTiR3vSsWe629HSK85YtfLLwpVMnPvHOm+ct9AMH\nsi8iL4/rTqZaKKEHKPSdOnn3hX3rW8B//xv62ECY0LcQqqo4Mm/FCmY9lJS409Y++IDvCxdS1E8/\nnesZGUxRM4xo0KYNcNFFTT8+L4+e9IABTKO8++7g+99yC0OT11/vf3v79pw79uKLvdt79qSYB3qy\nHTWKT76DB7vbunblE4Yj1s6Nqlu30N9r6FB32MbThkB2h4PF6FsIP/6xe8DG9OnuCZt79WIu73e+\nw6nbbrrJPHgjOcjL41PpWWfRgx44MPT+ZWXBQ0X+xDQ/350p5I9Rozhi1rcz2bPzuLEeva/QA0yi\nePjh0Mf7w4S+BXD0KOuItG/PuPz99wO//jW3XXMNa3Ns3850sOefj6+thhEueXkMN4bqzPWkKSW2\ne/akxx6I664L/WTSGKGfOtUd8okUFrppAXzyCfORn3ySccTBg+nRb9sGTJtGD/7EE4Ebb2Qqm2Ek\nA44YNkbom8I3vgGcfXbg7fn57JQNRmNCN23bNhwH0FxEm1HoWES2AdgL4BiAI6o6QUSyAbwIoA+A\nbQAuVtXdPsdpcz7XCM2WLRTt3FwORklP56g7VWbddOzIDp9169jWpg3zkg0jWfjqK4r8li3eGTGJ\nyP79/H90Zo5qKiICVW10cLW5Hr0CKFDVMarq3INuB7BAVQcDeN+13uJ4+unwZ7yPBvfdx5g7wMEm\n553HZREODsnKon25ufSMTOSNZCMvj6EQpwZTItOuHQdoxWt2tuZ69FsBjFPVao+2jQBOV9UKEekO\noFBVh/ocFzePvr6ePe/hxMqaQ34+q9F5plzFkilTmDXzr38xbFNTQ6/eYexYljjYsCE+9hlGJDh6\n1Pt3nerE06N/T0RWiIjTX52rqs7sixUAcv0fGh+efZYV7qKJKqc/27Ejup8TjOJiplJ+8gkHgPj+\nM+TnszKgYSQzLUnkm0NzL9OpqrpTRLoBWODy5v+HqqqIJFQwvrSUKYbRpL6eI9tiKfT33MPyBN/7\nHm80xcXs1Pn3v/0/VTgj9gzDSH2aJfSqutP1XiUirwGYAKBCRLqrarmI5AHwK6uzZ88GwPzTKVMK\ncMYZBc0xJWx27nRPNuygyhGi7doxlh2sOqMjounpgeu+OOePxOzt4fLRRyzUdPHFfJxt25Zzcb74\nIvDaaw33v/BCy5c3jESnsLAQhYWFzT+RqjbpBaAdgI6u5fYAPgYwFcDDAG5ztd8O4EE/x6rDJZeo\nfutbqkePakw45xzVvn2923buVG3fXnXwYNUnnwx+/Jtvct+MDNUdO/zvs2KFKqA6fbp3++LFqu+8\nE56djzyiumSJanW16p13qh47Fnz/vn1VzzhD9Re/UF25UnXkSNU5c2hHdXV4n2kYRmLj0s5G63Vz\nYvS5ABaLyGoAywC8qarvAngQwFki8iWAM1zrAW4yHH6/bRur0zWVvXs5iW84+PPo169nJblbb2UV\numC8+y6HWp9xRuB5KJ0UKt/QzcsvM5TiyZEjHLDka8999wEXXMD83YcfZrw9EEePcsTfI4+w8FFx\nMTMRxo/nBAfNSecyDCP5abLQq+pWVR3tep2kqr9xtdeo6pmqOlhVp6pPDr0nX3zBEMNjj4VXgvNn\nP6Og+bJgAXDHHQ2nCvNHeTmzbg4edLc5U5RNmhRa6D/4gCI/ciRrWfhj1y5W0fMN3WzZ0lD8V60C\nrr3Wu/2uu/j6/e9p0y23sGiSJ/fdB8yfz+XiYqZJjhzJejbLl7O0wbRpLLtqGEbLJq4jYz/6CPj6\n18Ob/La2ljO4/O53DbctXEjh3rzZ/7EVFaySd+wYRTgnh163w/r1HBk6bBi3B+qsraignWPHsr5F\nIKGvrma9iv37vXPpv/qqodAvX853J45eUsLr8sMfsv7MH/5Az/6NN+ixjx/PzKHHHgO++13e+LZu\n5YCRtDTWmH/xRV7TtDTrdDUMI85Cv3gxhT43lwJ75EjgfZcupRj/4x/eIg1w9vjevf3P9P7llxxO\nfNNNFMTsbNaYrqpy7+MIfVoaS6cuWkSx9E31X7iQNavT0yn0/j4P4Hfp2tV7CrDjx/n5xcXu9EtV\nzjIzYwYn+wBYdvW007wnRZgwgTcfZ/Lgq68G/u//mGFzyy28gTglVE89lU8OyTCIxDCM2BA3oT9y\nhILqCGdODuPnDz3EWLYvS5bQs73oInqzDpWV9IIvu8y/8L7yCqf2GjWKYYy8PIqwE6dXpdAPH871\nSZPoKc+a1fAJ4bXX3LPMOHM++hv9Wl3tFnrHgy8vZ9mBDh1o8xlnsIDY8uXAnXfS9ooKfk/fdMhW\nrRiqefllzo5TWsqpzH76U456Xb/eLfTOsb16Bbz0hmG0MOI23OChh5jOONQ1ZjY/n8L57ruswfLt\nb3vvv2QJY/QDB9Ljvflmerhbt3J97Fh6uatXM4zjTNlVXs56GAcOAG+/TaHPzKTQb9nCjsy0NPdI\n2YsvprB+/jlDKIMGsX3VKt6Y/vY3rqenM9Tz4YcceOTcKACeu0sXetVvv83ywN27045Dhzgn5Pr1\nvAbFxbT9ggv4tPLxx8Af/9jwenmWT3WKOWVn09t/7jn3MRMm0DYTesMwHOLm0d9zDzBnjjuX2xH6\njRuB999n4S2A2Tg/+hFDHBMnsvJiQQGXFy5k6d0bbuCciqtWcYqwe+91f05FBUNDo0ezI7V7d3rb\nVVXAFVdw1OiJJ7rtGDqUHaGTJ1N0HW6/HfjlL72rO44ZwyeMggLviXud0M2oUfS4r7uOnnv//ixH\n8Oab7Dg9dIh2paezXvyf/sSSBP7mkgzEhReyvIHj0bdr5775GIZhAHEU+lGjvGs45+fTy929m2K+\nZAlF/4knWCu9Tx93muBdd9GrnTeP26dPp6dfWcmCXR995I73V1RQ3EeP5ohVJ3RTUUGv/emngZ/8\npKF9p53G8wAMzyxa1HC2pd/9juGmE0/0zopxZnD/6U/5Hc49lxk0AwbQy3/9dd5g5sxxP7mMHcvv\nMGpUw0mLgzFjBt89pzmbNMkGQxmG4SZuoZslS7xnZMnPpwAOGUJh/OtfKaK33cYOx7o6976jRtHr\n96RVK8bpv/99CvKnn1Lwysvp0TvZJ3l5FMFXX2W75/yQnpx0Eo+tqmJ4aOhQFgnzxKlEd/31vGE4\nkxc7Hr3DZZcxf37AAG7bvJmhp8su8z7fQw81vmxCfj5vQpGeqMAwjNQhbh69r9eanw8sW0ZBnTmT\nYZxJkxi2SUtzz9ASjLlz6SlPmeKeB9Xx6DMz6fXm5bH4/8cfB54xHuCN42tfo1e/bh2FPxAXXcS+\ngU8+4brj0TtMm0bhHzjQnQ0zdmzD80yc2HC+ynD4+tfNgzcMIzAJM8NUfj47RocOZbx98WLggQcY\nimksU6bQ4z90iIOjnFrrd91F8e7alX0AwYQe4MwyH3xAoR8xIvB+bdoAf/87wzBffMFMHs/Z6Fu3\n5k1s4kQKfVpa6BlpDMMwIkVCCT3gzsJpDqeeytBNeTm9dydEdO21nP/RCav486o9Oessjrpduza4\n0AMshnbrrcDpp/P8vh62M6Bp+HB27HreCAzDMKJJwlRzdipBRkLoMzMptkuWMGzjizNvYyiPftQo\nZrSUlAQP3TjcfDO9+3ffDbxPhw6cnNswDCNWJIxHn5HBgUNDhkTmfCNG0BvP9TPtSU4OP8vfTcCT\ntDSGgdLTwy8l8IMf+C8LbBiGES8SRugBerpNicn7Y8QI4L33/It5enr4XvXUqTyXdXYahpGsJEzo\nJtKMGOGu6tgcLruM5QoMwzCSlYTy6COJk9USKjwTiowMDtYyDMNIVlJW6AcPZlpjcz16wzCMZCdl\nhb51a4ZvrB67YRgtHVHfouux+FARjcXnOlUkrSPVMIxUQESgqo1WtJQWesMwjFSiqUKfsqEbwzAM\ng5jQG4ZhpDgm9IZhGCmOCb1hGEaKY0JvGIaR4pjQG4ZhpDgm9IZhGCmOCb1hGEaKY0JvGIaR4kRF\n6EXkbBHZKCKbROS2aHyGYRiGER4RF3oRaQXgMQBnAxgO4FIRGRbpz4kFhYWF8TYhLJLBzmSwETA7\nI43ZmRhEw6OfAGCzqm5T1SMA/g1gRhQ+J+okyx8/GexMBhsBszPSmJ2JQTSEvieAYo/1ElebYRiG\nEQeiIfRWltIwDCOBiHiZYhGZCGC2qp7tWr8DwHFVfchjH7sZGIZhNIGEqEcvIukAvgAwBUAZgOUA\nLlXVDRH9IMMwDCMs0iN9QlU9KiI3AXgHQCsAfzORNwzDiB9xmWHKMAzDiB1RGxkrIm1EZJmIrBaR\nIhH5TYD9/uQaWLVGRMZEy56m2igiBSKyR0RWuV53x9JGH1tauWz4b4DtcbuWPnYEtDNRrqeIbBOR\nz102LA+wT9yvZyg7E+h6ZonIyyKywfW/NNHPPolwPYPaGe/rKSJDPD57lcuWH/vZr3HXUlWj9gLQ\nzvWeDmApgNN8tp8L4C3X8ikAlkbTnibaWABgXqztCmDrLQD+6c+eRLiWYdqZENcTwFYA2UG2J8T1\nDMPORLmezwK4xrWcDiAzQa9nKDsT4nq6bEkDsBNAr+Zey6jWulHVetdiBhivr/HZZTp44aGqywBk\niUhuNG3yJQwbAaDRvdyRRkTywT/wXPi3J+7XEgjLTgRpjzXB7EiI6+ki1PWK6/UUkUwAX1fVZwD2\n06nqHp/d4n49w7QTSJzf55kAtqhqsU97o69lVIVeRNJEZDWACgALVbXIZxd/g6vyo2mTL2HYqAAm\nuR6R3hKR4bG0z4PfA/g5gOMBtsf9WroIZWeiXE8F8J6IrBCR6/1sT5TrGcrORLie/QBUicjfReQz\nEXlaRNr57JMI1zMcOxPhejrMAvAvP+2NvpbR9uiPq+polxGTRaTAz26+d8+Y9g6HYeNn4KPTKAB/\nBvB6LO0DABE5H0Clqq5CcG8jrtcyTDvjfj1dnKqqYwCcA+CHIvJ1P/vE9Xq6CGVnIlzPdABjATyu\nqmMB7Adwu5/94n09w7EzEa4nRCQDwDcB/CfQLj7rQa9lTMoUux6P/h+AcT6bSgH08ljPd7XFnEA2\nqmqdE95R1fkAWotIdozNmwRguohsBfACgDNE5DmffRLhWoa0M0GuJ1R1p+u9CsBrYI0mTxLheoa0\nM0GuZwmAElX91LX+MiioniTC9QxpZ4JcT4A39pWuv7svjb6W0cy66SoiWa7ltgDOArDKZ7d5AK50\n7TMRwG5VrYiWTU2xUURyRURcyxPAlFR/cfyooap3qmovVe0HPs59oKpX+uwW12sZrp2JcD1FpJ2I\ndHQttwcwFcBan93ifj3DsTMRrqeqlgMoFpHBrqYzAaz32S3u1zMcOxPherq4FHSW/NHoaxnxAVMe\n5AF4VkTSwBvK86r6voh8DwBU9a+q+paInCsim8HHqKujaE+TbATwbQA/EJGjAOpBAYs3CgAJdi39\n0cBOJMb1zAXwmuv/OR3AP1X13QS8niHtRGJcTwD4EYB/ukIOWwBck4DXM6SdSIDr6bqpnwngeo+2\nZl1LGzBlGIaR4thUgoZhGCmOCb1hGEaKY0JvGIaR4pjQG4ZhpDgm9IZhGCmOCb1hGEaKY0JvGIaR\n4pjQG4ZhpDj/H0alzV+vYglmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x104e23710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def residual(p, just_generative=False):\n",
    "    v = p.valuesdict()\n",
    "    generative = v['a'] + v['b'] * x\n",
    "    M = 0\n",
    "    while 'a_max%d' % M in v:\n",
    "        generative += gauss(x, v['a_max%d'%M], v['loc%d'%M], v['sd%d'%M])\n",
    "        M += 1\n",
    "\n",
    "    if just_generative:\n",
    "        return generative\n",
    "    return (generative - y) / dy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create a Parameter set for the initial guesses\n",
    "def initial_peak_params(M):\n",
    "    p = lmfit.Parameters()\n",
    "    a = np.mean(y)\n",
    "    b = 1\n",
    "    a_max = np.max(y)\n",
    "    loc = np.mean(x)\n",
    "    sd = (np.max(x) - np.min(x)) * 0.5\n",
    "\n",
    "    p.add_many(('a', np.mean(y), True, 0, 10), ('b', b, True, 1, 15))\n",
    "\n",
    "    for i in range(M):\n",
    "        p.add_many(('a_max%d'%i, 0.5 * a_max, True, 10, a_max),\n",
    "                   ('loc%d'%i, loc, True, np.min(x), np.max(x)),\n",
    "                   ('sd%d'%i, sd, True, 0.1, np.max(x) - np.min(x)))\n",
    "    return p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving with `minimize` gives the Maximum Likelihood solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Variables]]\n",
      "    a:        1.51224223 (init= 10)\n",
      "    b:        10.5609028 (init= 1)\n",
      "    a_max0:   193.129507 (init= 141.2683)\n",
      "    loc0:     5.03878771 (init= 5)\n",
      "    sd0:      0.57125988 (init= 2)\n",
      "[[Correlations]] (unreported correlations are <  0.500)\n"
     ]
    }
   ],
   "source": [
    "p1 = initial_peak_params(1)\n",
    "mi1 = lmfit.minimize(residual, p1, method='differential_evolution')\n",
    "\n",
    "lmfit.printfuncs.report_fit(mi1.params, min_correl=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, this doesn't give a probability distribution for the parameters.  Furthermore, we wish to deal with the data uncertainty.  This is called marginalisation of a nuisance parameter. `emcee` requires a function that returns the log-posterior probability. The log-posterior probability is a sum of the log-prior probability and log-likelihood functions. The log-prior probability is assumed to be zero if all the parameters are within their bounds and `-np.inf` if any of the parameters are outside their bounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This is the log-likelihood probability for the sampling.\n",
    "def lnprob(p):\n",
    "    resid = residual(p, just_generative=True)\n",
    "    return -0.5 * np.sum(((resid - y) / dy)**2 + np.log(2 * np.pi * dy**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 310 (-4396.6022688087269, 0.49210220942131855)\n",
      "1 310 (-1045.339451412781, 6.7159454900820492)\n",
      "2 310 (-937.58116589602696, 9.1279318938875349)\n",
      "3 310 (-950.16931281126676, 10.295013691219992)\n"
     ]
    }
   ],
   "source": [
    "# Now lets work out the log-evidence for different numbers of peaks\n",
    "burn = 300\n",
    "thin = 10\n",
    "ntemps = 15\n",
    "workers=4\n",
    "log_evidence = []\n",
    "mini = []\n",
    "res = []\n",
    "\n",
    "# set up the Minimizers\n",
    "for i in range(4):\n",
    "    p0 = initial_peak_params(i)\n",
    "    mini.append(lmfit.Minimizer(lnprob, p0))\n",
    "    out = mini[i].minimize(method='differential_evolution')\n",
    "    res.append(out)\n",
    "\n",
    "total_steps = 310\n",
    "# burn in the samplers\n",
    "for i in range(4):\n",
    "    # do the sampling\n",
    "    out = mini[i].emcee(steps=310, ntemps=ntemps, workers=workers, reuse_sampler=False,\n",
    "                       params=res[i].params)\n",
    "    # get the evidence\n",
    "    print(i, total_steps, mini[i].sampler.thermodynamic_integration_log_evidence())\n",
    "    log_evidence.append(mini[i].sampler.thermodynamic_integration_log_evidence()[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 410 (-4397.2560033105292, 0.44918471141681948)\n",
      "1 410 (-1046.4683593227637, 7.4815351753807136)\n",
      "2 410 (-939.94542406160633, 9.4959066821616034)\n",
      "3 410 (-952.93565515210992, 11.469802405346286)\n",
      "0 510 (-4397.6793454040953, 0.40973253505489993)\n",
      "1 510 (-1047.3160213239353, 8.0491423392570596)\n",
      "2 510 (-941.66433605066152, 9.6985630480652389)\n",
      "3 510 (-955.08217059421099, 12.582124296094207)\n",
      "0 610 (-4397.9596953022838, 0.37919697193956381)\n",
      "1 610 (-1048.0376025386586, 8.4897199923516382)\n",
      "2 610 (-942.92595681693899, 9.8236085882018642)\n",
      "3 610 (-956.65961145532879, 13.516930818580931)\n",
      "0 710 (-4398.1945761284132, 0.35281168432993582)\n",
      "1 710 (-1048.6882834472437, 8.8220957020687365)\n",
      "2 710 (-943.98784364111691, 9.9501969794896468)\n",
      "3 710 (-957.83756672312495, 14.135558317186451)\n",
      "0 810 (-4398.3967256181213, 0.33012047874217387)\n",
      "1 810 (-1049.3368576955197, 9.0329401783053527)\n",
      "2 810 (-944.90031810750895, 10.074995996845018)\n",
      "3 810 (-958.82377961133204, 14.647967017922838)\n",
      "0 910 (-4398.5283877147976, 0.31542166970575636)\n",
      "1 910 (-1049.9134746623272, 9.1783057759282656)\n",
      "2 910 (-945.70917595371918, 10.176225190078185)\n",
      "3 910 (-959.67454304128751, 15.120676241524848)\n",
      "0 1010 (-4398.6317817432455, 0.30344262582366355)\n",
      "1 1010 (-1050.4668892836764, 9.2678993448282654)\n",
      "2 1010 (-946.45561122310528, 10.265865299859797)\n",
      "3 1010 (-960.50625989708055, 15.576003930653883)\n",
      "0 1110 (-4398.7365426230754, 0.29208846922301746)\n",
      "1 1110 (-1051.006721293078, 9.3212622555188318)\n",
      "2 1110 (-947.15669172052594, 10.337691691247301)\n",
      "3 1110 (-961.2669296071399, 15.994385843707164)\n",
      "0 1210 (-4398.843301324262, 0.27995088339775975)\n",
      "1 1210 (-1051.4665881258704, 9.372763350239893)\n",
      "2 1210 (-947.70042367552696, 10.285647497025138)\n",
      "3 1210 (-961.97475088305794, 16.398541104237665)\n",
      "0 1310 (-4398.9267439337, 0.27063312804421003)\n",
      "1 1310 (-1051.8695938759097, 9.4062550576870763)\n",
      "2 1310 (-948.12290229289226, 10.178115630579669)\n",
      "3 1310 (-962.67299781865472, 16.777003427183331)\n",
      "0 1410 (-4399.0036260521774, 0.26202137877953646)\n",
      "1 1410 (-1052.2767479562497, 9.4241563923483227)\n",
      "2 1410 (-948.54271923107467, 10.125763551047612)\n",
      "3 1410 (-963.29661823518734, 17.10479241831024)\n",
      "0 1510 (-4399.0765006333149, 0.25400269263082009)\n",
      "1 1510 (-1052.6521524042705, 9.4278280532398639)\n",
      "2 1510 (-948.92348959849403, 10.083195522663914)\n",
      "3 1510 (-963.88823570703426, 17.394351664440819)\n",
      "0 1610 (-4399.141515858485, 0.2468372792336595)\n",
      "1 1610 (-1052.9923541331409, 9.423718326379003)\n",
      "2 1610 (-949.2442946241639, 10.042476609640858)\n",
      "3 1610 (-964.45730090202721, 17.673802606047616)\n",
      "0 1710 (-4399.195774085063, 0.24094253244129504)\n",
      "1 1710 (-1053.318129304917, 9.4141151790629465)\n",
      "2 1710 (-949.54754436608403, 10.005417921886419)\n",
      "3 1710 (-964.98900777555787, 17.92882120033164)\n",
      "0 1810 (-4399.2417890272663, 0.23586072698708449)\n",
      "1 1810 (-1053.623259202602, 9.4004350987040652)\n",
      "2 1810 (-949.87576321555957, 10.029504080201264)\n",
      "3 1810 (-965.4925715460397, 18.151758665935176)\n",
      "0 1910 (-4399.2771360427232, 0.23213925552045112)\n",
      "1 1910 (-1053.9053202340979, 9.3831121179487127)\n",
      "2 1910 (-950.19561548823754, 10.041374424887522)\n",
      "3 1910 (-965.96998725706817, 18.358786651316223)\n",
      "0 2010 (-4399.300406624694, 0.22997833960562275)\n",
      "1 2010 (-1054.1744818659249, 9.3637019802408759)\n",
      "2 2010 (-950.48270060717118, 10.037358319010536)\n",
      "3 2010 (-966.44783189467739, 18.523598510719808)\n",
      "0 2110 (-4399.3188087653571, 0.22827471152231738)\n",
      "1 2110 (-1054.4271458009343, 9.3410102671191453)\n",
      "2 2110 (-950.78014284322364, 10.042700923908683)\n",
      "3 2110 (-966.87840898703723, 18.689471808336521)\n",
      "0 2210 (-4399.3374998254585, 0.22640657826832467)\n",
      "1 2210 (-1054.6693862868051, 9.3147736119333331)\n",
      "2 2210 (-951.08379255239402, 10.063762785431663)\n",
      "3 2210 (-967.31783185924769, 18.854763150536314)\n",
      "0 2310 (-4399.3559540748456, 0.22449021963348059)\n",
      "1 2310 (-1054.9038628174678, 9.2846819740852879)\n",
      "2 2310 (-951.38424226248958, 10.093324075223904)\n",
      "3 2310 (-967.76300886987246, 19.020488500822694)\n"
     ]
    }
   ],
   "source": [
    "# Now lets do the collection run\n",
    "# set up the Minimizer\n",
    "for j in range(20):\n",
    "    total_steps += 100\n",
    "    for i in range(4):\n",
    "        # do the sampling\n",
    "        res = mini[i].emcee(burn=burn, steps=100, thin=thin, ntemps=ntemps,\n",
    "                         workers=workers, reuse_sampler=True)\n",
    "        # get the evidence\n",
    "        print(i, total_steps, mini[i].sampler.thermodynamic_integration_log_evidence())\n",
    "        log_evidence.append(mini[i].sampler.thermodynamic_integration_log_evidence()[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x104ed5320>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEPCAYAAACOU4kjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucFNWZ//HPN6IIikZW4w28T1xQI4iKuZnJmri4+Qma\nYIQYFyMmhgneQmIEkzhJ1HhHjIqJskHI6pKsBmElKBImbrIiXgBRZIPXCKIbUUERlYHn90fVQDsO\nMz0z3VPdPd/369WvqTp1uuopS+aZc07VKUUEZmZmhfSRrAMwM7PK4+RiZmYF5+RiZmYF5+RiZmYF\n5+RiZmYF5+RiZmYFV3LJRVKtpBWSFqafE3K2jZW0XNIyScfnlA+QtCTdNiGbyM3MrEHJJRcggOsi\non/6+QOApL7AqUBfYBBwsySl35kIjIyIKqBK0qAsAjczs0QpJhcANVE2BLgzIjZExAvAM8BASXsC\nPSJiQVpvCnBSx4RpZmZNKdXkco6kxZImSfpoWrYXsCKnzgpg7ybKV6blZmaWkUySi6Q56RhJ489g\nki6u/YF+wCrg2ixiNDOztuuSxUEj4ov51JN0GzAzXV0J9M7Z3IukxbIyXc4tX9nEvjyJmplZG0RE\nU0MVzSq5brF0DKXBycCSdHkGMEzSdpL2B6qABRHxCrBW0sB0gP90YHpT+46Iiv1ccsklmcfg8/P5\ndbZz6wzn11aZtFxacKWkfiR3jT0PnA0QEUsl/RZYCtQDNbHlzGuAyUA3YFZEzO7wqM3MbLOSSy4R\n8a/NbLscuLyJ8seAw4oZl5mZ5a/kusWsbaqrq7MOoah8fuWrks8NKv/82krt6VMrJ5Kis5yrmVmh\nSCIqYUDfzMzKn5OLmZkVnJOLmZkVnJOLmZkVnJOLmZkVnJOLmZkVnJOLmZkVnJOLmZkVnJOLmZkV\nnJOLmZkVnJOLmZkVnJOLmZkVnJOLmZkVnJOLmZkVnJOLmZkVnJOLmZkVnJOLmZkVXCbJRdIpkp6S\ntFHSEY22jZW0XNIyScfnlA+QtCTdNiGnvKukaWn5fEn7duS5mJnZh3XJ6LhLgJOBX+YWSuoLnAr0\nBfYGHpBUlb6feCIwMiIWSJolaVBEzAZGAqsjokrSqcCVwLCOPBmzziAC1q2Dt95KPmvXbllubv2d\nd0CCj3yk+c8227Rcpy11y3Xf+dZVq19A3DEySS4RsQySdzM3MgS4MyI2AC9IegYYKOlFoEdELEjr\nTQFOAmYDg4FL0vK7gBuLHL5Z2aivh7ffzj8RNLe+bh1svz306JF8dtppy3Lj9X333bLevXsSy6ZN\nzX82bmy5Tkt16+vbt99CxdHeuq2pD8VNcm2VVctla/YC5uesryBpwWxIlxusTMtJf74EEBH1ktZI\n6hkRr3dAvGYFFQHvvde6X/zN1XnvPdhxx5aTwU47Qa9eW9/eo0eyny6l9hvDiChuQhwwoG1xFe1/\nFUlzgD2a2DQuImYW67jNqa2t3bxcXV1NdXV1FmFYhWlrd9HW1qHlZNCjB+y+Oxx0UPMtie7dS7fb\nxApDKmzSr6uro66urt37UTKckQ1J84AxEfF4un4RQERcka7PJunyehGYFxF90vLhwLERMSqtUxsR\n8yV1AVZFxG5NHCuyPFcrLYXuLurateVkkO96165Z/9cx20ISEdHqP1FKoZGbG/QM4A5J15F0d1UB\nCyIiJK2VNBBYAJwO3JDznREk3WlDgbkdFrl1mHy6i1qTHHK7i5r7Zb/TTrD33s0nA3cXmX1YJi0X\nSSeTJIddgTXAwog4Id02DjgTqAfOi4j70vIBwGSgGzArIs5Ny7sCU4H+wGpgWES80MQx3XIpQxFw\nyikwY0ay3HicoK0the7d2zdYadZZtLXlkmm3WEdycilP//mf8JOfwPz5Hj8wy4KTSwucXMrPunXQ\npw/85jdw7LFZR2PWObU1ubhjwErWpZfCZz/rxGJWjtxysZL0v/8Ln/40PPEE7LVX1tGYdV5uuVjF\niIBzz4WxY51YzMqVk4uVnN//HlasSBKMmZUnd4tZSXnnnWQQf/Jk+Pzns47GzNwtZhXh8svhk590\nYjErd265WMlYvjxJLIsXJ0/Fm1n23HKxshYB550HP/iBE4tZJfCMSFYSZsyA55+H6dOzjsTMCsHd\nYpa59euhb1+47TY47risozGzXO4Ws7J1xRVw1FFOLGaVxC0Xy9Szz8LRR8OiRdC7d9bRmFljbrlY\nWTrvPPj+951YzCqNB/QtMzNnJrcf33131pGYWaE5uVgm1q9PWi233ALbbZd1NGZWaO4Ws0xcdRX0\n7w/HH591JGZWDB7Qtw73/PPJ3WGPPw777JN1NGbWHA/oW9k4/3z47nedWMwqWSbJRdIpkp6StFHS\nETnl+0laL2lh+rk5Z9sASUskLZc0Iae8q6Rpafl8Sft29PlY/mbNgqefhjFjso7EzIopq5bLEuBk\n4MEmtj0TEf3TT01O+URgZERUAVWSBqXlI4HVafl44MpiBm5t9+67yTtabrgBunbNOhozK6ZMkktE\nLIuIv+ZbX9KeQI+IWJAWTQFOSpcHA7eny3cBfs67RF1zDRx2GAwa1HJdMytvpXgr8v6SFgJrgB9G\nxJ+BvYEVOXVWpmWkP18CiIh6SWsk9YyI1zsyaGveCy/A9dfDo49mHYmZdYSiJRdJc4A9mtg0LiJm\nbuVrLwO9I+KNdCxmuqRDChVTbW3t5uXq6mqqq6sLtWtrwQUXJM+17Ldf1pGYWXPq6uqoq6tr934y\nvRVZ0jxgTEQ83tx2YBXwx4jok5YPB46NiFGSZgO1ETFfUhdgVUTs1sS+fCtyRmbPhtGj4cknYfvt\ns47GzFqjnG9F3hy0pF0lbZMuHwBUAc9FxCpgraSBkgScDtyTfm0GMCJdHgrM7bDIrUXvvQfnnAMT\nJjixmHUmWd2KfLKkl4BjgHsl/SHd9DlgcTrm8jvg7Ih4M91WA9wGLCe5o2x2Wj4J+AdJy4HzgYs6\n6jysZddeC336wJe+lHUkZtaR/IS+Fc3f/pZM8fLII3DAAVlHY2ZtUc7dYlahvvvd5LkWJxazzqcU\nb0W2CjBnDixcCFOnZh2JmWXBLRcruPffTwbxr78eunXLOhozy4KTixXc+PFw0EFw4olZR2JmWfGA\nvhXUihXQrx88/DAceGDW0ZhZe3lA30rCmDFQU+PEYtbZeUDfCuaBB2DBAvj1r7OOxMyy5paLFUTD\nIP748dC9e9bRmFnWnFysICZMSCalHDIk60jMrBR4QN/abeVKOPxweOghqKrKOhozK6S2Dug7uVi7\nDR+ePIV/2WVZR2Jmhebk0gInl+KYNw/OOAOWLoUddsg6GjMrNN+KbB1uw4bkPS3jxzuxmNkHOblY\nm/3iF9CrF5x8ctaRmFmpcbeYtcmqVXDYYfA//wMf/3jW0ZhZsXjMpQVOLoX19a9D797w859nHYmZ\nFVNbk4uf0LdWe/DB5PP001lHYmalymMu1iobNsB3vpO8vtiD+Ga2NU4u1io33QS77w5Dh2YdiZmV\nsrySi6Tukg4u1EElXS3paUmLJd0taeecbWMlLZe0TNLxOeUDJC1Jt03IKe8qaVpaPl/SvoWK0z7o\nlVfg0kuTu8TU6h5YM+tMWkwukgYDC4H70vX+kma087j3A4dExOHAX4Gx6b77AqcCfYFBwM3S5l9j\nE4GREVEFVEkalJaPBFan5eOBK9sZm23FhRfCmWdCnz5ZR2JmpS6flkstMBB4AyAiFgIHtOegETEn\nIjalqw8DvdLlIcCdEbEhIl4AngEGStoT6BERC9J6U4CT0uXBwO3p8l3Ace2JzZr25z8nT+P/6EdZ\nR2Jm5SCf5LIhIt5sVLapyZptcyYwK13eC1iRs20FsHcT5SvTctKfLwFERD2wRlLPAsbX6dXXJ4P4\n11wDPXpkHY2ZlYN8bkV+StJpQBdJVcC5wP+09CVJc4A9mtg0LiJmpnUuBt6PiDtaEXOb1dbWbl6u\nrq6murq6Iw5b9iZOhF13ha9+NetIzKzY6urqqKura/d+WnyIUtIOwMVAw+D6fcDPIuLddh1YOgP4\nJnBcw74kXQQQEVek67OBS4AXgXkR0SctHw4cGxGj0jq1ETFfUhdgVUTs1sTx/BBlG7z6Khx6KPzp\nT9C3b9bRmFlHK9rElRGxLiLGRcSR6efiAiSWQcD3gSGN9jUDGCZpO0n7A1XAgoh4BVgraWA6wH86\ncE/Od0aky0OBue2JzT7oootgxAgnFjNrnRa7xSQ9AAxtGHdJxzPujIh/bsdxfwFsB8xJbwZ7KCJq\nImKppN8CS4F6oCanuVEDTAa6AbMiYnZaPgmYKmk5sBoY1o64LMdDD8H998OyZVlHYmblJp9usUUR\n0a+lslLnbrHW2bgRjjwSvv99+NrXso7GzLJSzPe5bMx9MFHSfhT2bjErQbfcAjvvnLxl0systfJp\nuQwCfgU8mBYdC3wrp1uqLLjlkr+//z0ZY5k3LxnMN7POq6hT7kvaDTgGCGB+RLzW+hCz5eSSv5Ej\nYaedkjdMmlnnVuwp97cDXk/r900P9mAL37EyNH8+/OEPnk7fzNonn7vFriSZ72spsDFnk5NLhdm4\nMXkS/6qrkvEWM7O2yqflcjJwcES8V+xgLFu33pq8o+W007KOxMzKXT7J5VmSbjEnlwr22mvw4x/D\nAw94On0za798kst6YJGkuWxJMBER5xYvLOto48Yltx1/4hNZR2JmlSCf5DIj/TTcaqWcZasAjzwC\nM2d6EN/MCiffW5G7A/tERNlOBOJbkZu2aRMcc0wykD9iRMv1zaxzKdoT+jlvopydrhfiTZRWIiZN\ngm23hdNPzzoSM6sk+Tyh/zjwTyRT3vdPy56MiLJ6dtstlw9bvTp5Ev+++6BfWc0UZ2YdpZhzixX7\nTZSWkYsvhlNOcWIxs8Ir2psorbQ9+ihMn+5BfDMrjnxaLucAh5DchnwnsBY4v5hBWXFt2gSjR8Pl\nl8Muu2QdjZlVorzuFqsEHnPZYtIkuO02+Mtf4CP5/HlhZp1WwWdFljQzZzVInm/ZvB4Rg1t7sCw5\nuSRefz0ZxJ81C444IutozKzUFSO5VKeLJwN7AL8hSTDDgVcjoqy6xpxcEt/5DkTAzTdnHYmZlYOi\nvc9F0mMRMaClslLn5AILF8IJJ8DSpdCzZ9bRmFk5KOatyN0lHZhzoAOA7q09UC5JV0t6WtJiSXdL\n2jkt30/SekkL08/NOd8ZIGmJpOWSJuSUd5U0LS2fn/tKZtti06ak1XLppU4sZlZ8+SSXC4B5kv4k\n6U/APNp/t9j9wCERcTjwV2BszrZnIqJ/+qnJKZ8IjIyIKqAqff0ywEhgdVo+HriynbFVpClTkve1\nnHlm1pGYWWfQ4nMuETFb0seBfyQZ2F/W3ne7RMScnNWHga80V1/SnkCPiFiQFk0BTiKZkmYwcEla\nfhdwY3tiq0RvvgljxyaTU/ruMDPrCFv9VSPpuPTnV4B/AQ4EDgK+JOnLBYzhTGBWzvr+aZdYnaTP\npGV7Ayty6qxMyxq2vQQQEfXAGknu+Mnxox/B4MFw5JFZR2JmnUVzLZdjgbnAiTQ9xf7dze1Y0hyS\nu8waGxcRM9M6FwPvR8Qd6baXgd4R8YakI4Dpkg5p4RzyVltbu3m5urqa6urqQu26ZC1aBNOm+Ul8\nM8tPXV0ddXV17d5PPneLdUlbBAUl6Qzgm8BxEfHuVurMA8YAq4A/RkSftHw4cGxEjJI0G6iNiPmS\nugCrImK3JvbV6e4Wi4DPfjaZ8fjss7OOxszKUTHvFntO0q8kHScV5gW46WD894EhuYlF0q6StkmX\nDwCqgOciYhWwVtLANIbTgXvSr80AGt5EMpSktWXA1Knw7rtw1llZR2JmnU0+LZcdgP8HDAOOAGYC\n0yLiv9t8UGk5sB3welr0UETUpOM7PwE2kMy8/OOIuDf9zgBgMtANmNXwmmVJXYGpQH9gNTAsIl5o\n4pidquWyZg306ZNMTnn00VlHY2blqmgPUTY6yC7ADcDXImKb1h4sS50tuZx/PqxbB7femnUkZlbO\n2ppc8plyv2EqmFOBQcAjwFdbeyDrOEuWwB13JE/im5llIZ9usReARcA0YGZEvN0BcRVcZ2m5RMDn\nPgfDh8OoUVlHY2blrpgtl8MjYk0bYrIM3HFH0h32rW9lHYmZdWb53C22h6S5kp4CkHS4pB8WOS5r\ng7Vr4cIL4cYbYZuyGhEzs0qTT3K5FRgHvJ+uP0Ey7b6VmJ/8BP75n+GTn8w6EjPr7PLpFuseEQ83\nPOISESFpQ3HDstZ68slkcsqnnso6EjOz/Fouf5d0UMOKpKEkT8xbiYiA0aPhkkvgYx/LOhozs/xa\nLqOBXwEHS3oZeB44rahRWav8x38kMx9/+9tZR2Jmlsj7IUpJO6b13ypuSMVRqbciv/VW8iT+tGnw\n6U9nHY2ZVZpizi0GQPp8y52tPYAV109/Cl/4ghOLmZWWvJ7Qz7F3y1WsoyxdCpMnJ4P5ZmalpLXv\nJVxUlCis1SLgnHOSF4HtvnvW0ZiZfVCrkktEfKNYgVjr/O538NprUFOTdSRmZh+Wz9xiS0jeRJk7\noLOGZALLSyNidfHCK5xKGtB/++1kEP+OO5KXgZmZFUsx5xabDdQDd5AkmGFAd+BVkvernNjag1r7\nXHopVFc7sZhZ6cqn5bIwIvo3VSZpSUQcVtQIC6RSWi7LliVJ5YknYM89s47GzCpdMW9F3kbSwJwD\nHZ3zvfrWHtDarmEQf9w4JxYzK235dIuNBH6dPkQJ8BYwMn398c+LFpl9yF13wapVyVQvZmalrDVP\n6O8MUK7vdin3brF165JB/KlTk5eBmZl1hKJ1i0n6qKTxwB+BP0q6tiHRtJWkn0laLGlR+q6Y3jnb\nxkpaLmmZpONzygdIWpJum5BT3lXStLR8vqR92xNbqbrssmSsxYnFzMpBPgP6dwNLgNtJ7hY7HfhE\nRHy5zQeVejTMUSbpHJK3XZ4lqS/JXWlHkcwG8ABQlU7zvwAYHRELJM0CboiI2ZJqgEMjokbSqcDJ\nETGsiWOWbcvlr3+FT30qGcTfa6+sozGzzqSYA/oHRsQlEfFcRDwbEbXAga2OMEejyS93BF5Ll4cA\nd0bEhoh4AXgGGChpT6BHRCxI600BTkqXB5MkPoC7gOPaE1upiYBzz4WxY51YzKx85JNc1kva/ESF\npM8A77T3wJIuk/Q34Ay23BiwF7Aip9oKkhZM4/KVbJnnbG/gJYCIqAfWSOrZ3vhKxfTp8NJLSYIx\nMysX+dwt9m1gSs44yxvAiJa+JGkOsEcTm8ZFxMyIuBi4WNJFwPVA0aeWqa2t3bxcXV1NdXV1sQ/Z\nLu+8AxdcAL/+NWy7bdbRmFlnUFdXR11dXbv306a7xSSdHxHXt/voyX73AWZFxKFpoiEirki3zQYu\nAV4E5kVEn7R8OHBsRIxK69RGxHxJXYBVEbFbE8cpuzGXH/0Ili9PXgZmZpaFjnify5qc25DHtPZA\nuSRV5awOARamyzOAYZK2k7Q/UAUsiIhXgLWSBkpquKngnpzvNLSkhgJz2xNbqXjmGZg4Ea65JutI\nzMxar7XvcymUn0s6GNgIPAuMAoiIpZJ+Cywlefq/Jqe5UUMyl1k3kpbO7LR8EjBV0nJgNcncZ2Wt\nYRD/wguhV6+sozEza728u8U+8CXppYjo3XLN0lFO3WL33AMXXQSLF8N222UdjZl1ZgWfFVnS2yRT\n7Tele2sPZPlZvx7OPx9uvdWJxczK11aTS0TsuLVtVjxXXAFHHglf+ELWkZiZtV2busXKUTl0iz37\nLAwcCAsXQu+y6nQ0s0pV9LvFrPjOPx++9z0nFjMrf1ndLWaN/Nd/JXOI3XVX1pGYmbWfk0sJePdd\nOO+85LkWD+KbWSVwt1gJuOoq6NcPjj++5bpmZuXAA/oZe/55OOooePxx2GefrKMxM/sgD+iXqQsu\nSD5OLGZWSTzmkqFZs+Cpp2DatKwjMTMrLCeXjLz7bjJ/2C9+AV27Zh2NmVlhuVssI9dcA4ceCiec\nkHUkZmaF5wH9DLz4IgwYAI8+Cvvtl3U0ZmZb5wH9MnLBBclzLU4sZlapPObSwe67D554Au64I+tI\nzMyKxy2XDvTee3DOOTBhAmy/fdbRmJkVj5NLB7ruOvjHf4QvfSnrSMzMissD+h3kb3+DI46ABQvg\ngAMyC8PMrFU8oF/ixoyB0aOdWMysc8gkuUj6maTFkhZJmiupd1q+n6T1khamn5tzvjNA0hJJyyVN\nyCnvKmlaWj5f0r5ZnFNzHngAHnsMfvCDrCMxM+sYWbVcroqIwyOiHzAduCRn2zMR0T/91OSUTwRG\nRkQVUCVpUFo+Elidlo8HruyIE8jX++8nLZbrr4du3bKOxsysY2SSXCLirZzVHYHXmqsvaU+gR0Qs\nSIumACely4OB29Plu4DjChhqu40fDwceCCeemHUkZmYdJ7PnXCRdBpwOvAMck7Npf0kLgTXADyPi\nz8DewIqcOivTMtKfLwFERL2kNZJ6RsTrxT6HlqxYAVdfDfPng1o9HGZmVr6KllwkzQH2aGLTuIiY\nGREXAxdLuoikO+sbwMtA74h4Q9IRwHRJhxQqptra2s3L1dXVVFdXF2rXTRozBmpq4KCDinoYM7OC\nqauro66urt37yfxWZEn7ALMi4tAmts0DxgCrgD9GRJ+0fDhwbESMkjQbqI2I+ZK6AKsiYrcm9tWh\ntyLPnQtnnZVMqd+9e4cd1sysoMrqVmRJVTmrQ4CFafmukrZJlw8AqoDnImIVsFbSQEki6U67J/3+\nDGBEujwUmNsBp9Cs999PnsQfP96Jxcw6p6zGXH4u6WBgI/AsMCotPxb4qaQNwCbg7Ih4M91WA0wG\nupG0dGan5ZOAqZKWA6uBYR1zClt3ww2w774wZEjWkZiZZSPzbrGO0lHdYi+/DJ/4BDz0EFRVtVzf\nzKyUtbVbzMmlwL72Ndh/f7jssqIfysys6JxcWtARyaWuDkaMgKVLYYcdinooM7MOUVYD+pVow4bk\nSfzrrnNiMTNzcimQG2+EvfaCL38560jMzLLnbrECWLUKDjsM/vIXOPjgohzCzCwTHnNpQTGTy9e/\nDr16wRVXFGX3ZmaZaWtyyWxusUrx4IPJZ+nSrCMxMysdHnNph/r6ZBD/2mthxx2zjsbMrHQ4ubTD\nTTfBxz4GQ4dmHYmZWWnxmEsbvfJKMoj/4IPQp0/BdmtmVlI8oN+CQieXESNg993hqqsKtkszs5Lj\nAf0O9Je/JFPqP/101pGYmZUmj7m0Un09fOc7cM010KNH1tGYmZUmJ5dWuuUW6NkTTj0160jMzEqX\nx1xa4f/+Dw45JJmg8pCCvXzZzKx0eUC/BYVILt/4RtJqufbaAgVlZlbiPKBfZA89BPff70F8M7N8\neMwlDxs3JoP4V18NO+2UdTRmZqUv0+QiaYykTZJ65pSNlbRc0jJJx+eUD5C0JN02Iae8q6Rpafl8\nSfsWOs5f/jJJKsOHF3rPZmaVKbPkIqk38EXgxZyyvsCpQF9gEHCzpIa+vonAyIioAqokDUrLRwKr\n0/LxwJWFjPPvf4fa2uR9LWp1r6OZWeeUZcvlOuDCRmVDgDsjYkNEvAA8AwyUtCfQIyIWpPWmACel\ny4OB29Plu4DjChnk2LFw2mlw6KGF3KuZWWXLZEBf0hBgRUQ8oQ82B/YC5uesrwD2Bjakyw1WpuWk\nP18CiIh6SWsk9YyI19sb58MPw6xZHsQ3M2utoiUXSXOAPZrYdDEwFjg+t3qx4mirhkH8K6+EnXfO\nOhozs/JStOQSEV9sqlzSocD+wOK01dILeEzSQJIWSe+c6r1IWiwr0+XG5aTb9gFeltQF2HlrrZba\n2trNy9XV1VRXV281/ttug27dkrdMmpl1FnV1ddTV1bV7P5k/RCnpeWBARLyeDujfARxN0t31AHBQ\nRISkh4FzgQXAvcANETFbUg1wWESMkjQMOCkihjVxnLwfoly9Gvr2TZ5rOfzwgpymmVlZKueHKDf/\nxo+IpZJ+CywF6oGanIxQA0wGugGzImJ2Wj4JmCppObAa+FBiaa1x45K5w5xYzMzaJvOWS0fJt+Xy\nyCMweHAyiP/Rj3ZAYGZmJaytLRc/oZ9j06ZkEP+KK5xYzMzaw8klx6RJsO22cPrpWUdiZlbe3C2W\nev116NMH7rsP+vXrwMDMzEqYp9xvQUvJZdQo2GabZJoXMzNLlPPdYpl77DH4/e/9JL6ZWaF0+jGX\nhkH8yy+HXXbJOhozs8rQ6ZPL5MnJzzPOyDIKM7PK0qnHXN54IxnEv/deGDAgo8DMzEqYB/Rb0FRy\nGT06maBy4sSMgjIzK3Ee0G+lhQvhd7/zIL6ZWTF0yjGXhkH8Sy+Fnj1brm9mZq3TKZPLlClJd9jI\nkVlHYmZWmTrdmMubbyaD+DNnwpFHZh2VmVlp84B+CxqSy7nnwnvvwS9/mXVEZmalzwP6eVi8GKZN\ng6VLs47EzKyydaoxl9Gj4ac/hX/4h6wjMTOrbJ0quaxfD2edlXUUZmaVr1Mll5tuSmY+NjOz4up0\nA/pmZpa/snzNsaQxkjZJ6pmu7ydpvaSF6efmnLoDJC2RtFzShJzyrpKmpeXzJe2bxbmYmdkWmSUX\nSb2BLwIvNtr0TET0Tz81OeUTgZERUQVUSRqUlo8EVqfl44Erix17Kaqrq8s6hKLy+ZWvSj43qPzz\na6ssWy7XARfmU1HSnkCPiFiQFk0BTkqXBwO3p8t3AccVMshyUen/g/v8ylclnxtU/vm1VSbJRdIQ\nYEVEPNHE5v3TLrE6SZ9Jy/YGVuTUWZmWNWx7CSAi6oE1Dd1sZmaWjaI9RClpDrBHE5suBsYCx+dW\nT3++DPSOiDckHQFMl3RIsWI0M7Pi6PC7xSQdCswF3kmLepG0RI6OiP9rVHceMAZYBfwxIvqk5cOB\nYyNilKTZQG1EzJfUBVgVEbs1cVzfKmZm1gZlMf1LRDwJ7N6wLul5YEBEvC5pV+CNiNgo6QCgCngu\nIt6UtFbSQGABcDpwQ7qLGcAIYD4wlCRxNXXcVv/HMTOztimFucVyWxTHAj+VtAHYBJwdEW+m22qA\nyUA3YFarw7tAAAAGzUlEQVREzE7LJwFTJS0HVgPDOiRqMzPbqk7zEKWZmXWcipv+RdIgScvShyp/\nsJU6N6TbF0vq39ExtkdL5yepWtKanAdRf5hFnG0h6d8kvSppSTN1yvnaNXt+ZX7tekuaJ+kpSU9K\nOncr9cry+uVzfmV+/baX9LCkRZKWSvr5Vurlf/0iomI+wDbAM8B+wLbAIqBPozr/QtKtBjAQmJ91\n3AU+v2pgRtaxtvH8Pgv0B5ZsZXvZXrs8z6+cr90eQL90eUfgfyvs314+51e21y+Nv3v6swvJGPZn\n2nP9Kq3lcjTJE/4vRMQG4D+AIY3qbH7oMiIeBj4qaXfKQz7nB1tu7S4rEfHfwBvNVCnna5fP+UH5\nXrtXImJRuvw28DSwV6NqZXv98jw/KNPrBxARDXfwbkfyh+zrjaq06vpVWnLZ/EBlagVbHrZsrk6v\nIsdVKPmcXwCfSputsyT17bDoiq+cr10+KuLaSdqPpIX2cKNNFXH9mjm/sr5+kj4iaRHwKjAvIhq/\nVrFV168U7hYrpHzvTmj810W53NWQT5yPkzyI+o6kE4DpwMeLG1aHKtdrl4+yv3aSdgT+Ezgv/Qv/\nQ1UarZfV9Wvh/Mr6+kXEJqCfpJ2B+yRVR0Rdo2p5X79Ka7msBHrnrPfmg9PGNFWn4SHOctDi+UXE\nWw3N24j4A7BtBU2HU87XrkXlfu0kbUsyv99vImJ6E1XK+vq1dH7lfv0aRMQa4F7gyEabWnX9Ki25\nPEoyY/J+krYDTiV5yDLXDOBfASQdA7wZEa92bJht1uL5SdpdktLlo0luN2/cd1quyvnataicr10a\n9yRgaURcv5VqZXv98jm/Mr9+u0r6aLrcjWTG+oWNqrXq+lVUt1hE1EsaDdxHMiA1KSKelnR2uv2X\nETFL0r9IegZYB3wjw5BbJZ/zI5mlYJSkepIpdsrmoVJJdwKfA3aV9BJwCcldcWV/7aDl86OMrx3w\naeDrwBOSGn4pjQP2gYq4fi2eH+V9/fYEbpf0EZJGx9SImNue351+iNLMzAqu0rrFzMysBDi5mJlZ\nwTm5mJlZwTm5mJlZwTm5mJlZwTm5mJlZwTm5mLWRpDpJAzrgOOem06BPLeIxXijHp8mtdFXUQ5Rm\nHazND4lJ6hIR9XlWHwUcFxEvt/V4efADb1ZQbrlYRUunynla0q/SlzzdJ2n7dNvmlkc6/cXz6fIZ\nkqZLul/S85JGS/qepMclPSRpl5xDnJ6+GGqJpKPS7++g5MVgD6ffGZyz3xmS5gJzmoj1u+l+lkg6\nLy27BTgAmC3p/Eb1z5B0j5KXWP1V0o9ztn09Pf5CSbekT14j6WZJj6T/LWqbiKGbpD9IGimpu6R7\nlbxAaomkr7bjUlgn4+RincFBwI0RcSjwJvCVtDzY+l/shwAnA0cBlwFrI+II4CHS+ZVIZojtFhH9\ngRrg39Lyi4G5ETEQ+Cfgaknd0239ga9ExOdzD5YmuTNI3tlzDPBNSYdHxLeBl4HqrcxpdRTwZeAT\nwCmSBkjqA3wV+FQa2ybgtIbYIuIo4HDgc5IOzdlXD5L5o/49IiYBJwArI6JfRBwGzN7KfyuzD3G3\nmHUGz0fEE+nyYyRv8mzJvIhYB6yT9CYwMy1fQvKLHJLEdCckLwKTtJOS6cqPB06U9L20XleSOagC\nmBMRbzZxvM8Ad0fEegBJdwPHAotbiPP+iHgj5zufATYCA4BH03kUuwGvpPVPlfRNkn/7ewJ9gSdJ\nEuU9wJURcWda9wngGklXAP8VEX9uIRazzZxcrDN4L2d5I7B9ulzPltb79nxQ7nc25axvovl/Nw0t\noS9HxPLcDZIGkkz4t7Xv5b4rQ7Q8DtJ4e+53bo+IcY2Ovz8wBjgyItZI+jVbzjuAP5O0VhoS5nIl\n70n/EnCppLkR8bMWYjID3C1mnVPDL/EX2PLOiqGt/G7D8qkAkj5DMgX5WpJZq8/dXCn5Bd34u439\nN3BSOuaxA3BSWtZSLF+UtEs6TfoQkgQxFxgqabf0+D0l7UPS7bUOWKvk9bQnNNrfj4E3JN2Ufm9P\n4N2I+HfgGuCIFuIx28wtF+sMGv+F37B+DfBbSd8ieTlS5GyPJuo33hbAu5IeJ/m3dGZa/jPgeklP\nkPwB9xzJ+8e3OsYTEQslTQYWpEW3RkRDl9jWWjCR1r+L5MVNUyPicQBJPwTuTwfyNwA1EbEgnS5+\nGcnraj/UzRUR56U3I1xJkqSulrQJeJ/krjWzvHjKfbMyJekMYEBEnJN1LGaNuVvMrHw1d7ebWabc\ncjEzs4Jzy8XMzArOycXMzArOycXMzArOycXMzArOycXMzArOycXMzAru/wNGWdhYWMwNkAAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105394438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(log_evidence[-4:])\n",
    "plt.ylabel('Log-evidence')\n",
    "plt.xlabel('number of peaks')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Bayes factor is related to the exponential of the difference between the log-evidence values.  Thus, 0 peaks is not very likely compared to 1 peak. But 1 peak is not as good as 2 peaks. 3 peaks is also not as good as 2 peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r01 = np.exp(log_evidence[-4] - log_evidence[-3])\n",
    "r12 = np.exp(log_evidence[-3] - log_evidence[-2])\n",
    "r23 = np.exp(log_evidence[-2] - log_evidence[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 1.10153937106e-45 12978005.6223\n"
     ]
    }
   ],
   "source": [
    "print(r01, r12, r23)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These numbers tell us that zero peaks is 0 times as likely as one peak. Two peaks is 9e44 times more likely than one peak. Two peaks is 1.3e7 times more likely than three peaks.  However, caution has to be taken with these values. The log-priors for this sampling are uniform but improper, i.e. they are not normalised properly.  Internally the lnprior probability is calculated as 0 if all parameters are within their bounds and `-np.inf` if any parameter is outside the bounds. The `lnprob` function defined above is the log-likelihood alone. Remember, that the log-posterior probability is equal to the sum of the log-prior and log-likelihood probabilities. Extra terms can be added to the lnprob function to calculate the normalised log-probability.  These terms would look something like:\n",
    "\n",
    "$\\log \\left(\\prod_i \\frac{1}{\\max_i - \\min_i}\\right)$\n",
    "\n",
    "where $\\max_i$ and $\\min_i$ are the upper and lower bounds for the parameter, and the prior is a uniform distribution. Other types of prior are possible. For example, you might expect the prior to be Gaussian."
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}