Expected_value_from_jointdist.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating expected value of the observable given a model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $p_{\\rm m}(D|r)$ is a model and $x_i$ is an observed value, probability that $x_i$ is drawn from the model is simply $p_{\\rm m}(D|x_i)$. In reality even the observed value $x_i$ is not single-valued and in fact has its own proability distribution $p(x)$. In which case the probability it is drawn from the model is calculated using the joint distribution as follows:\n",
    "\n",
    "$$p_{\\rm joint}(D,x) = p_{\\rm m}(D|x) p(x)$$\n",
    "$$p_{\\rm total}(D) = \\int p_{\\rm joint}(D,x') dx'$$\n",
    "\n",
    "And the refined pdf of the observed point becomes\n",
    "\n",
    "$$p_{\\rm refined}(D,x)=p_{\\rm joint}(D,x)/p_{\\rm total}(D)$$\n",
    "\n",
    "Finally, the expected value of the observable can be calculated using\n",
    "\n",
    "$$<x> = \\int x' p_{\\rm refined}(D|x') dx'.$$\n",
    "\n",
    "Following code implements this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from scipy import stats as st"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PDF calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Refinement probability: 0.0022\n",
      "Direct probability: 0.0000\n",
      "Model mean; expected value: -15.0, -14.8, -0.0\n"
     ]
    }
   ],
   "source": [
    "xt = np.linspace(-50,50,100000)\n",
    "\n",
    "modelmean, modedisp = -15, 1\n",
    "model_func = lambda x: st.norm.pdf(x, modelmean, modedisp)/5.\n",
    "\n",
    "#measurement\n",
    "data = 1\n",
    "\n",
    "#measurement with uncertainty\n",
    "datadisp = 10\n",
    "data_pdf = st.norm.pdf(xt, data, datadisp)\n",
    "\n",
    "#refined data pdf given the model\n",
    "data_ref_pdf = data_pdf*model_func(xt)\n",
    "\n",
    "refinement_total_probab = np.trapz(data_ref_pdf, xt)\n",
    "\n",
    "print 'Refinement probability: %1.4f'%refinement_total_probab\n",
    "print 'Direct probability: %1.4f'%model_func(data)\n",
    "\n",
    "data_ref_pdf_normalised = data_ref_pdf/refinement_total_probab\n",
    "\n",
    "#Expected value\n",
    "data_exp1 = np.trapz(xt*data_ref_pdf_normalised, xt)    # refined probability is normalised\n",
    "data_exp2 = np.trapz(xt*data_ref_pdf, xt)    #refined probabilit not normalised\n",
    "\n",
    "print 'Model mean; expected value: %1.1f, %1.1f, %1.1f'%(modelmean, data_exp1, data_exp2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's plot the distributions for a further insight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ff1288edd10>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEcCAYAAAD6GqKbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl83HW97/HXd7Ys3dJJt7S0TdPSUgQsXVCp2kpTUDwu\neFLE4+FSzpUGOeo9gDTu5+g5wkldznF5eGhxwyNqSRBFUaSpFM8VvF0itSICbbrSPUvbNMtkZr73\nj1kymUzSNs3MLzPzfj4eeTTz+37nl8+kmfn8vuvPWGsREREZbi6nAxARkdykBCMiImmhBCMiImmh\nBCMiImmhBCMiImmhBCMiImmhBCMiImmhBCMiImmhBCMiImmhBCMiImnhcTqATDLG/BGYCOx2OhYR\nkSwyBzhhrb36Qp6UVwkGmDhu3LhpCxYsmOZ0ICIi2eKFF17g1KlTF/y8fEswuxcsWDBty5YtTsch\nIpI1li9fzrPPPnvBPT8agxERkbRQghERkbTIty4yEUmTnp4eDh06RFdXl9OhyEUqLCzkkksuwev1\nXtR5RkSCMcasBZoAP4C1dsMFPHe9tbY6XbGJyPk5dOgQY8aMoby8HGOM0+HIEFlraW5u5tChQ8ya\nNeuizuV4F5kxphZotNbWRxPLbGNM1QU8d3FaAxSR89LV1UVpaamSS5YzxlBaWjosLVHHEwywxlrb\nkPB4E3DOFokxpiJ9IYnIUCi55Ibh+n90NMEYYxamONwCVJ7H0yuJJCMRERmBnG7B+IkklERtAMaY\nkoGeZIypBB5NY1ySxfaf3s+Xtn2J7lC306HICLRy5cqUxxsaGtiw4byHf+U8OD3IX0J0YD9BLOH4\niSabVM+z1rapOS6p/OnEn/jBX37ANVOuYdn0ZU6Hk3fKP/HksJ5v37+/c1jPJ5njdIJJlUBiCSe5\nZQOAMabKWlt/rhMbY9YAa5IOzzty5MiFRShZ523T3wbA7rbdSjA5qr6+no0bN+L3+2lpaaGiooKG\nhgY2b95MS0sLNTU1QKS1smbNGpqamli1ahWVlZW0tEQ+WpqamqitrQWguloTUdPB6QTTQqQVk6gE\nwFrbL/lEB/YHatX0EZ2R1qe9a4zZUlZWpk+cHFfgLmCUdxQnO086HYqk0ZIlS1i7di2zZ8+mrq6O\ndevWsX37dtavX89DDz1ESUlJPMHU1tZSV1dHRUVFvIustraWkpLIx09DQwMLF6YaEpaL4WiCsdY2\nGmOSE4YfaEhVH1gIVCRMDlgClETX0dRba5vSFKpkkWUbl3G256wSjEMy1aUVSw4VFRXxx21tkY+T\nlpYWSkpKaGlpoa2tLV43WXV1dfz5DQ0DfezIUDndggHYkNTttRJYHyuMtloWRtfJ9Okai3aDVVhr\n12UuXBnJQuEQZ3rOAHAmcMbhaMQJtbW1VFdXU1JS0uff5C6ympoaampqqKioYMmSJQMmIRk6xxOM\ntbbGGLM2uriyAtiTlEgqgVVAquSyikiLZi2wIVW3muSXWFJZu2Qtt15+q8PRSLpUVfWuxd60KbJa\nYc2aNf2OxVRUVLBjx45+x+rq6tIYpTg9TRkAa+26aAtlXfI2MdbaDdbafvMKY8etteOjz1NyEU4F\nIvesKCnQ1Wi+GGjacTrOX1NTowkBF2BEJBiR4XKqO5Jg/nDkD3zh+S84HI04qa2tbViTQWx8Z/36\n9eeoKTGOd5GJDKfxheO5/Yrbee3Ma/yq6Vd89o2f1fYlDli+fHm/YzfffDN33XUXHR0d3Hjjjf3K\nV69ezerVqzl58iRVVVUMdmPAVNOOGxsbWb9+PS0tLdTW1rJ+/Xq2b9/Ohg0bWLx4cZ+y2MA+RAb3\nE2eU1dXVpTz/hg0bqK+vZ+XKlVRWns9mI6IWjOSU6WOmc8+ie5hfOp+gDdIV0tbxuSg27bi2tha/\nP7J0buHChaxfv57q6mrq6+v55Cc/yeLFi1mzZk2/smSrVq2Kj8e0tbWlPH9VVRWVlZVKLhdALRjJ\nKZ3BTgKhAKO9owE423OWIk+Rw1Hln8FaH8XFxYOWT5gwYdByIOWMr8bGRjZu3JhyWnJyWU1NTZ8u\ntFj3F/ROcZaLpwQjOeXxVx/nga0PULMkspL7bM9ZJhRNcDgqGW6pph1v376dpqamPutitm/fzrp1\n6ygpKelTFlvBD5Eusk2bNrFt2zb8fj8VFRUpzy8XTglGckqsS6y0qJRxBePoCfU4HJGkQ6ppx2vW\nrOkzVRnoUye5LNGqVav6lA80rVkD/BdGCUZySlcwkmBuKL+Bd8x6h8PRiOQ3JRjJKZ3BToo8RbiM\n5q/I+dGgffroXSg5pTPYSaG7kNfaX+PuZ+5m54mdTockaZbuhZYXurhSCzN7qQUjOWX59OXMHDuT\n7mA3DQcauL78el4/8fVOh5VXarfW8teWv170eS7zX0bNNTVDfn5bWxs1NTUXNW5yMYsrtTBTLRjJ\nMW+e9mZuvfxWCj2FQO+YjOSWpqYmFi1aRE1NTZ+FlrHZX01NTTzwwAPxhZbJZYkaGhpYtWoV1dXV\nNDU1UV1dTXV1NY2NjfHFlck7LTc0NLBy5UpWrVoVP2dyPAM9N5+oBSM55djZY7hd7vjal85gp8MR\n5Z+LaXWcr1T3d4ktpmxoaIgvtGxra4vPDkssW7t2bZ/z+f3++ELMxHvEVFVVsWfPnpTjNLGZZ6tW\nreLzn/98v3gGe26+UIKRnPLJ//tJQuEQD658EEAr+XPUcC+0XLRoUbxu4j1iEls7gy3OLCwsHL4X\nl0OUYCSndAW7GFswlgJ3AVNHTaXAXeB0SJIGw73QMib5HjGJd7kcbHFmTU2NFmamYKy1TseQMcaY\nLcuWLVt2rm0oJHvd9PObKB9bzn+87T+cDiXvvPTSS8yfP9/pMDKioaGBpqamQRdvZrvE/8/ly5fz\n7LPPPmutXX4h59Agv+SUzmBnfIBfRJylBCM5pSvYFR/g/9T/fIpv7/q2wxFJujmxDqaysnLA1otT\n62AaGhriM9c2bNhw7idEpfP3pzEYySkfW/gxLhl9CQC7Tu4iEA44HFF+uv2p2/sdu6H8Bm657BY6\ng53c1XBXv/L3zHkP753zXlq7Wrlnyz187+3fu6gY8nUdTOLtpJ2mBCM55X2Xvi/+fYG7gEBICSYX\nDfcNx9avXx8frI8N5ldXV8enNSffZCz5JmW1tbXnfYOy+vp6Nm7ciN/vp6SkhNraWpqamqipiUzv\nnjp1KocPH8bv97Ny5cp43ZaWFioqKmhoaGDz5s00NTX1eU2JscVmv8U27Ex+XSUlJRmZlKAEIzkj\nbMO80voKk4on4S/0RxKMWjCOGKz1UeQpGrR8fOH4c7Zesn0dzJIlS1i7dm28bk1NDQ899BAlJSVc\nffXVXHPNNaxfv576+vp43dmzZ1NXV8e6devYvn07lZWVfV5T4ow3gD179lBSUkJ1dXWfZNjQ0MCe\nPXv6xZsOSjCSM7qCXaz6xSruXnQ3/3DFP+B1e9WCyVHZtg4m+bmp4o/d6Ky9vb3PbLxY3cTp121t\nbYO+Xogk4dgOBlOnTu3zumKtpXRTgpGc0ROO3PvF5/IBUD62nGA46GRIkibZtg4m8bnJW9XEymOt\npxtvvJHi4uJz/g6SX2+ympqaeHny68rUDdW0DkZyxomOE1xXdx2ffeNnuXnezU6Hk3e0Dia3aB2M\nSILYeIvX5XU4EhEBdZFJDomNt/jckS6yB3c+yIsnX+QbK77hZFiSg/J5A8sLoRaM5IzSolLuf/P9\n8fu/HD17lBebX3Q4Kkm3dC+0lKFTC0ZyxljfWN41+13xx16XV9OUHbJ8+fKLPsfFjpUOx0JLuThK\nMJIzTgdOs7t1N5eOv5QxvjFaaJnDhnOhZeLCx+TFjImLGGNTjBPPE5sqnLhoUnopwUjOeKn5JT70\n9If47g3fZcmUJfjcPiUYh6R7puZwL7QcaDFjXV1dnwWKa9eu7XOe2LTfxEWT0ksJRnJG8iD/9DHT\nWTBpAaFwCLfL7WRoMsyGe6HlQIsZoe/Cy+TzVFRUpIxFIpRgJGfExltiCy1vuvQmbrr0JidDkjQZ\nzoWWqRY+xiQvUGxpaRl0caP0pYWWkjOe2vsU9/3uPn72np8xu2S20+HknXxaaJkPtNBSJEFyC+bp\nfU/zrsffxcnOk06GJZK3lGAkZyyZvIT/XP6flBaVAtAV6mLf6X10Bjsdjix/5FOPSC4brv9HjcFI\nzigbXUbZ6LL441hLRjPJMqOwsJDm5mZKS0sxxjgdjgyRtZbm5uZz7hB9PpRgJGccOH2AA2cO8Kay\nN+F2ueOzyZRgMuOSSy7h0KFDnDhxwulQ5CIVFhZyySWXXPR5lGAkZ/xm32/4+h+/zo6/34EbNwXu\nAgC6Q90OR5YfvF4vs2bNcjoMGUE0BiM5I3k35QlFE3jLtLcwyjvKybBE8taIaMEYY9YCTYAfwFq7\nYZC6JcAaoA2YHa2fmduzyYgWCAXwurzx/v95/nl8q/JbDkclkr8cTzDGmFpgk7W2IfbYGFNlra0f\n4CmfTEwoxpgdxpg1gyUlyQ+BUCA+7iIizhsJXWRrYsklahNQPUj9KmNM4m3kmgBtAiT0hHviM8cg\nMui/4tEVbN6/2cGoRPKXoy0YY8zCFIdbgMHu5rPSWpu4t0MFsHFYA5Os9MH5H+SG8hvij40xHO88\nTntPu4NRieQvp7vI/EQSSqI2iIy1WGvbkp+QmFxiCcpauy6dQUp2mDVuFrPG9c5i0iwyEWc5nWBK\niA7sJ4glHD/RZJMsOtB/M7AKuGOAOmuITAZINO/IkSNDDlZGth3HdtAV7GLptKVA70LLnnCPk2GJ\n5C2nE0yqBBJLOMktm7hoy2YDsCE6yL8+eZA/+rjPMWPMlrKysmUXGbOMUN//8/c52nG0N8FooaWI\no5we5G8h0opJVALxJNJPtPWSaH30S/JcIBzoM8jvc/u4ofwGZo6d6WBUIvnL0RaMtbbRGJOcSPxA\nQ6r6xphKYJMxZnxyAhpozEbyRyAUwOv2xh97XB6+vOzLDkYkkt+cbsFApJurKuHxShJaJMaYioTy\n7cCGpESyEqhXcpHkFoyIOMvxBBNdNFlhjKmKrujfk7TIspLouphoEllvjFkb/aoFmqy1qzIfuYw0\nPaGefgstb6i/gS9t+5JDEYnkN6cH+YHBpxknD9ZbaxuBxkzEJdnlgbc8gMv0vWbqDnXTEexwKCKR\n/DYiEozIcEh1m+QCd4FmkYk4xPEuMpHh8os9v6DxWN/Grc/toyekdTAiTlCCkZzx5e1f5smmJ/sc\n87q98W38RSSz1EUmOSPVIP/1M69nXME4hyISyW9KMJIzAuG+62AA7nz9nQ5FIyLqIpOcYK2N3A8m\nxTqYsA07EJGIKMFITgiGg1hs/HbJMR/Z/BE+8OQHHIpKJL+pi0xygtvl5hfv/QVjC8b2Oe5xeTRN\nWcQhSjCSE1zGRfm48n7HfS6ftusXcYi6yCQntAfa+f6fv8+rra/2Oe51e9WCEXGIEozkhNbuVr6y\n4yv8pfkvfY773D4lGBGHqItMckKsGyx5Hcy1U69lYtFEJ0ISyXtKMJITYtvBJE9TXjlzJStnrnQi\nJJG8py4yyQmxbrDkhZY94R5OB05jrXUiLJG8pgQjOSG231hyF9l3d32XpT9eStAGnQhLJK+pi0xy\nwlUTrmLzqs399h2LJZyeUE+/RZgikl5qwUhO8Lq9TCqeRIG7oM/xWILRTDKRzFOCkZywu3U333rh\nWzR3Nvc5Hk8w2rJfJOOUYCQnvNL6Cv+18784FTjV53hsVplaMCKZpwQjOSE+yJ80Tfky/2XcteAu\nRntHOxGWSF7TIL/khFgLJXkW2Tz/POb55zkRkkjeUwtGckJ8JX9SCyYQCnD07FG6Q91OhCWS15Rg\nJCcM1IJpPN7IyvqV7Dqxy4mwRPKaEozkhL+//O95/gPPU+Qp6nM8Nm1Zs8hEMk9jMJITvC4vXl//\nhZSxLrPYXmUikjlqwUhO+O2B3/K1xq/1Ox7bm0wtGJHMU4KRnLD16FY2/nVjv+OxFowG+UUyTwlG\nckIgFOi3kzJAaVEpH1/8cS73X+5AVCL5TWMwkhMCoUDKzSzH+MZw2+tucyAiERm0BWOMaTbG/Fem\nghEZqkA40G+KMkDYhmk61URLV4sDUYnkt3N1kY0HdiQfNMa8aoy5Lj0hiVy4UDjUb5ElRFo27/nZ\ne/jpqz91ICqR/HauLrJGoBL4dtLx2UBJWiISGYKvLP8KYRvudzzxfjAiklnnSjCfAJ42xswGtieV\nfdIY8/5BnmuttbdcVHQiF8Bl+jfIXcaFx3g0TVnEAYMmGGttgzFmMVALvJ/eVosFFkW/Bnw6oAQj\nGfHtXd/G5/Lxv173v/qVed1ebdcv4oBzTlO21jZaa1daa/3WWpe11gUYoCr2eIAvd/rDF4l45sAz\n/P7w71OWFbgLlGBEHDDUaco1RMZnREaEQDiQcpAfYO2StUwfMz3DEYnIkBKMtfZLwx2IyMUYaKEl\nwLtmvyvD0YgIaCW/5IhAKPU6GIDdrbvZe2pvhiMSkUFbMMaYFiKD9UNhrbUTzqeiMWYt0AT4o0/c\ncB71AZYA26y164YYo+SIAnfBgLdF/sT/fIKy0WV847pvZDgqkfx2ri6yHQycYCqJJIW2pOMLo/+e\n1xiNMaYW2GStbYg9NsZUWWvrB6pvra1JeLzDGIOSTH772Xt/NmCZz+3TOhgRB5xrmvLKVMeNMX8L\n+K21i1OUlRBZM/Oh84xhTWLCADYRmUTQL8FEz92cdHg9kWnUSjCSktfl1ToYEQcMdQzmE8BPUhVY\na9vo/dAflDFmYYrDLURaR6n4gVpjTEVSfe0qkOfW/m4tv97765RlmqYs4oyhTlNeRKSlMZBSBk4S\nifxEEkSiNoi0VqLJKs5a22SMWWStbUo4vBJoSD6xMWYNsCbp8LwjR46cR1iSTay1/Hrvr5kxZkbK\ncp/bR1t3ck+uiKTbUBPMH4EaY8xGa+3OxAJjTDmRD/bzGYMpITqwnyCWcPz0H9/BWhs/b7TL7GZS\n7CgQnSjQZ7KAMWZLWVnZsvOIS7JIMBwEGHAW2e1X3K4bjok44GIWWj4NNBpj6ujdp2w2va2GmlRP\nTJLqsjKWcM5nf/U6YEVSi0byTGx8ZaCFlosmD7ajkYiky1AXWsb2KHuISAvi5oTiRqDGWvvb8zhV\nqvGTkujPGLRPIzr7rDaxRSP5KTa+MtBCy72n9nKy8yRLpizJZFgieW/Id7SMfrAvMsaMA2KD7k3W\n2lMXcg5jTHIi8ZNiTCWRMaaKvlObFyrR5K+QDTG5eDJjfWNTlv/wLz9k0/5N/O6W32U4MpH8dtG3\nTI4mlD9exCk2JK17WUlkFhoA0RljC2PlxphKokkoOgbjJ7LTsxJMnppQNIGGVQNfk/jcPk1TFnHA\nRSWY6F0tV9HbgtkBPGqtfeF8z2GtrTHGrI22SiqAPUmLLCujP6M+mlBis9fWJ9RJuShTBKIJRtOU\nRTJuSAkm2i22nUhCMAlFK4nMLqu7kJuNDbYKP3E2WHRcxgxUV/LTwdMH+eLWL/Lh13+Y1098fb9y\nn9tHT7gHay3G6M9HJFOGutCyjsiMsc3A7IT7xMwBfgusMsbcP0wxigyqrbuN37/2e051px7+i80u\n6wlruxiRTBpqF1klsMNae33iweh04ZXGmD1ANfCpi4xP5Jxi4yteV+pZZG8vfzuXl16O2+geeCKZ\nNNQWTCOwcZDyB4d4XpELFhtfGWih5fSx01k6bSlulxKMSCYNNcFsJLJV/kBmA48O8dwiFyTW9TXQ\nQssj7UfYtH8THT0dmQxLJO8NtYusEfh3Y8z/BpLv5LQQuIPIOMx1iQXnufhS5IL43D5mjZvFKO+o\nlOXbj23nU//3Uzx505PM8Kber0xEht9QE0xsqvBD9L9fTGyaTl3SMQuoj0KG3RvL3sgT731iwPLY\nCn9NVRbJrKEmmDsZ+p0uRTIq1nWmxZYimTXUvcgGvaWxSCZt3r+ZH/zlB3ztbV+jpLD/rYFig/9q\nwYhk1lAH+UVGjKMdR2k8PvBOQQXuAkDrYEQy7aL3IhNx2rl2U77Mfxk/eMcPmFMyJ5NhieQ9JRjJ\neueapjzGN4arJ12dyZBEBHWRSQ6ItWA8rtTXS2cCZ/j57p9z8PTBTIYlkveUYCTrTSyayIKJCwbc\nyLK1q5XP/P4zvHDivDf5FpFhoC4yyXrvv+z9vP+y9w9YrllkIs5QC0ZyXjzBaB2MSEYpwUjW+8Yf\nv8GdDXcOWB5faKkWjEhGKcFI1jt05hAHTh8YsDzWgtE6GJHM0hiMZL2ecM+AU5Qhcp+Y+nfVM7F4\nYgajEhElGMl6gVBgwEWWAMYY5vnnZTAiEQF1kUkOCIQCg7ZgAB575TG2Hd2WoYhEBJRgJAfM88/j\nqolXDVrna41f46m9T2UoIhEBdZFJDrh38b3nrON1ezVNWSTD1IKRvOBz+TRNWSTDlGAk693269v4\n0rYvDVrH5/ZpmrJIhinBSNY7dOYQ7T3tg9bxudWCEck0jcFI1usOd59zFtk3r/vmoFOZRWT4KcFI\n1guEAhR6CgetM3nU5AxFIyIx6iKTrGatpSvYFd8OZiC/2fcbfvrqTzMUlYiAEoxkubANUzmzkkvH\nXzpovV82/ZIfvfSjDEUlIqAuMslybpebry7/6jnr+Vw+rYMRyTC1YCQvaBaZSOYpwUhWO9J+hGt/\nfC2/3vvrQev53D56QloHI5JJSjCS1TpDnZwJnMFaO2g9r0tbxYhkmhKMZLVYt1eBp2DQencvuptf\nve9XmQhJRKI0yC9ZrSvYBUCBe/AEM8o7KhPhiEgCtWAkq8VbMOdIMFuPbOWr279K2IYzEZaIMEIS\njDFmrTGmyhizxhiz5jyfU2WMqU13bDKyjS8cz7tnv5tJxZMGrfenk3/iey9+TzPJRDLI8S6yaJLY\nZK1tiD02xlRZa+sHqF8JLARWAk2Zi1RGokvHX8oX3/zFc9aL7VUWCAcoZPBtZURkeIyEFsyaWHKJ\n2gRUD1TZWttgrV0HNKY9MhnxzjV7LCa2lYxaMCKZ42iCMcYsTHG4BajMdCySnR7f/ThX//fVHD17\ndNB6sQSjtTAimeN0C8ZPJKEkagMwxpRkPhzJNl3BLoLh4Dk3u/S6Ilv1ay2MSOY4PQZTQiTJJIol\nHD/RZDMU0ckCyRMG5h05cmSop5QRKNblVegefFzlxlk38o5Z78DjcvpP/vwsX74cgC1btjgah8jF\ncPrdliqBxBJOcsvmglhrNwAbEo8ZY7aUlZUtu5jzysjSFYqsgzlXC8btcmciHBFJ4HQXWQuRVkyi\nEgBr7ZBbL5I/ukPdeIznnC2Tvaf28m9/+Df2n96fochExNEEY61tpH8rxg80pKgu0s9VE67i7+b/\n3TnrtXa1svHljbzW/loGohIRcL4FA7DBGFOV8HglsD72wBhTkVQuEve2GW/jviX3nbNekacIgM6e\nznSHJCJRTo/BYK2tia3kByqAPUmLLCuBVUA9xKc2VwJVgN8YswdoiLaGJM90BjvxuDzxWWIDKfYW\nA9AR7MhEWCLCCEgwANGFkwOV9RmsjyaSRmDA50j+WPvsWo52HKXuXXWD1ou3YIJqwYhkykjoIhMZ\nss5gZzx5DKbYU4zX5SUYDmYgKhGBEdKCERmqzmAno32j+xzrCYXxuAzGmPix0b7RNN6qXlSRTFIL\nRrJaR7CjTwvmp42HeP3nn2bpv/+WFw+fcjAyEVGCkayW2EXWejbApx7fRUcgxOFTXXzqp7v61P3i\nH77IY6885kSYInlJCUay2gcu+wArZqwA4KkXj9LV03tDsZ2HTtF0oj3++NlDz9J4XN1kIpmiMRjJ\nare97rb49/+vqblf+da9LVRMjIzRFHmKNItMJIPUgpGsZa3leMdxuoKR/ci27WvtVyfxWJGnSOtg\nRDJICUayVmewkxV1K3jkpUfoCAR5ra1/6+TlY6fj3xd7i7WSXySDlGAka8W6u4o8Rexv7m2ZlBT3\nrurffbydcDhy18sJRRMo8p57zYyIDA+NwUjWiiWYYm9xnwSzYHoJf37tNCfbu+nqCXOwtYOZpaNY\n91Zt/iCSSWrBSNbq24I5Gz8+01/M3Mm9iy9fPdbe77kikn5KMJK1YgP2RZ4i9iW0YGaWjmLWhFHx\nx/tbImWPv/o492y5J7NBiuQxdZFJ1ppcPJl7F93LnJI57G8+ED9ePqGYnlDvepiD0QRz4MwBnjnw\nTMbjFMlXSjCStaaMmsLqK1YDsL/5r/HjM0tHEQj2JphY91mxp5igDRIIBc55i2URuXhKMJK1Wrta\nOR04zcTCMg6fiozHuAxcMr6Irp5QvN6BaAtmjG8MAGcCZygtKs18wCJ5RmMwkrWe2PMEf/P437D7\nRCs2MhOZqSVFFHjczPAXx+sdbO0kHLbxBHM6cDrV6URkmCnBSNY6EziDwXDsVG93WHlpZHB/TKEX\n/6hIN1ggGObYmS4mFE1g1rhZhMKhlOcTkeGlLjLJWu097Yz2juZAS1f82IzS3pbLdH8xLWcDAOxv\n7uCNFW/gifc+kfE4RfKVWjCStc4EzjDKN6rPGpjyhAQzM6Gb7ECz9iATyTQlGMla7YFIC2Z/0hqY\n3u97E8z+lrOc6j7F6qdW8/S+pzMap0i+UoKRrHXLZbdw14K7klowvQkmcaB/f3MHPrePHcd2cPDM\nwYzGKZKvNAYjWetNU99ETyjMh1ufih9LTCqJrZkDLR0UugvxuDycCZzJaJwi+UoJRrLWzhM76eke\nQzC6W/KUsYUU+dzx8j5dZM0dGGMY6xs7YhJMa2srL7zwAnv37uXEiRO0trbS0tLCJz7xCQDa2tr4\n0Ic+xPjx4/H7/ZSWllJeXs7SpUsZNWrUOc4u4jwlGMladzx9B2+c+E5gMdB3BhnApDEFFHpddPWE\nOdXZw6mOHsb4xmQ8wYTDYXbt2sXWrVvZvn07H/nIR7jyyitpaGjg5ptvjtfzer34/X7uuusuAAKB\nAL/61a9oa2ujs7P3PjYvv/wyc+fOpb6+nl/+8pcsWrSIJUuWsGjRIrxeb7+fL+IUJRjJSl3BLjqD\nnfQEeu8BM0oAAAATZklEQVTvUp6UYIwxzPAX80p0N+X9LWd5XenrKBtVlpEYDx06xL333svmzZtp\nbo7czrmkpIS3v/3tXHnllSxbtoxNmzZRUVHB5MmTKS4uxhgTf/6kSZP4y1/+AkBnZycnT55k7969\nzJo1C4DDhw/z1FNP8fDDDwMwZswYli1bxmOPPYbPp61wxHlKMJKV2rrbAOjoLIwfK5/Qv9tohn9U\nb4Jp7qD2rbVpicday0svvcRjjz3GjBkzuO222ygpKWHr1q28853vZMWKFSxdupSKiop4Epk0aRKV\nlZXndf6ioiKmT5/O9OnT48c+9rGP8dGPfpTDhw/z/PPPs3nzZg4dOhRPLvfccw9er5f3ve99LFmy\nBJdLc3oks5RgJCu1dLUAcKq990q9IkWCSRyHie1JNpz27t3Lww8/zE9+8hNefvllAO644w5uu+02\nRo8ezd69e4f9ZyYyxjBt2jSqqqqoqqqKH7fWsn//fp544gnWrVvHtGnTuOWWW7j99tt53etel9aY\nRGJ0SSNZqa0r0oI5fqr3GilVC6bvQP9ZHn35Uf7m8b+5qO1ienp64t9XV1fzhS98gWnTpvGtb32L\n1157jQ0bNgz53MPFGMNjjz3GiRMn+OEPf8jixYv5+te/ziOPPAJAMBjktddeczhKyXVKMJKV5vrn\n8qW3foVjJ0vix2b6U3WR9Z1J1hPuYf/p/fEutvNlreX555/njjvuYPLkyRw+fBiAr371q+zfv5/N\nmzfz4Q9/mKlTpw7xFaVHSUkJH/zgB/nZz37G4cOHufvuuwF4+umnmT59Otdffz0/+tGP6OjQTgcy\n/JRgJCtNKJrA3NHXEgpGBvnLxvWdohyTuBZmf3MHE4omANDc1XxeP6e1tZXa2lrmz5/Ptddey49+\n9CPe/e53EwhE9ji74oor+oyLjGQTJkxg4sSJQCTuz372s7z66qt88IMfZMqUKdxxxx2cOnXK4Sgl\nlyjBSFbadWIXTzc9H3+cuII/0SXji/C4IoPqR093UeSOtHhOdp4c8NyBQIBDhw4B0NHRwac//Wkm\nTpzId77zHY4ePcr3v/99ysvLh+mVOGPGjBl8/vOfZ8+ePTzzzDP87d/+Lc899xxjxkRuadDQ0MCB\nAwfOcRaRwWmQX7LSd/78HXYcfhn4CADzpoxJWc/rdjFn0mj+ejSy9uX0mUiLJ1WCeeGFF/je977H\nI488wqJFi/jNb37DtGnTOHDgwIjr+houLpeL5cuXs3z5csLhMC6Xi1AoxK233sqxY8e47rrrWL16\nNTfddJMWd8oFUwtGstLRs0ch2Dv+cnnZ2AHrzk8oO95ayJvK3kRJQe9zf/zjH7NgwQKuvvpqHnzw\nQa677jruueeeeHmuJpdksWnMbreb5557jn/5l3+hqamJW2+9lSlTpvDd737X4Qgl2yjBSFY6evYo\nnZ2j44/nD5JgEpPPq8e6+Mayb9C6ozU+sH3w4EG8Xi/f/OY3OXLkCI8++ig33HBD+oLPArNmzeJz\nn/scu3fv5tlnn+Xmm2/m0ksvBeDFF1/kC1/4Avv27XM2SBnxlGAk67QH2mnuaub0mUgrxGXg0smj\nB6w/v2wsNhSkc892flhbw+TJk3nvTe/l5z//OQD33nsv27Zt4x//8R/x+/0ZeQ3ZwuVy8da3vpXv\nfOc7vOUtbwHgmWee4Z//+Z+ZNWsW1113HQ8//DDt7e0ORyojkRKMZJ2mU00AhLonAXD51LEUevvP\nIIuZ4DrLoW/eyvH6f+Hwzv9h0d2Lufbha+MLE93ugZ8r/X3kIx9h3759/Ou//isHDhxg9erVXHrp\npQSDQYD4vyIa5JesM3f8XG4Yfz/1L3cDcE15abzs0KFDPPnkk/zyl79kypQpPPTQQ8yrmMGMN93I\nWf9cisoXsnTZER7f/y3aetqY6J3o1MvIajNnzuQzn/kMn/70p3nuued49dVX8XgiHyeLFy+mtLSU\nd7/73dx4443MmTOnzx5rkj/UgpGsU+gp5M9N4yAcmRH2hgo/DzzwAJdddhnTp0/nzjvv5MUXX6Ss\nLLKppTGG2+/+HMVz3oDxeGlrmQLAn07+ybHXkCuMMSxdupTVq1cDkdbLjTfeyLFjx/inf/on5s6d\ny4wZM0bE7gaSeSOiBWOMWQs0AX4Aa+2gf40XWl9yQ2dnJzt37uQ/tn2Lbb/eR/ufjjDz9v/gzXMm\n8Lv2dubMmcOHPvQhbrzxRubPn9/nqnnZvIl89/eRfcGee6kAzwwPO4/vZMWMFU69nJzk8Xi4//77\nuf/++9mzZw+bNm3imWeeYfz48QC88sorXH/99VxzzTXxr4ULFzJ69MBjaJK9HE8wxphaYJO1tiH2\n2BhTZa2tH476kn0CgQCvvvoqL730EitWrGD8+PE8+OCDfPSjH8VV6mJu7VxchW24R5XxpmkFjCrw\n8MUvfnHQcy6dXcqE0QWcbO/m5BnLguIr2XxgM3cvulvdN2kye/ZsZs+ezZ133hk/1t3dzRve8Aa2\nbt1KXV0dEJlIsHnzZpYvX84rr7xCY2Mjl19+OXPnzqWwsHCg00sWcDzBAGustTUJjzcBNcBACeNC\n68sI0tHRwYkTJzh58iTHjx/nqquuYtq0aWzdupX77ruPAwcOcPDgQUKhyGaUDQ0NrFixgkWLFnHf\nfffxYvlB9tidjLri3xk1v4Tqt199Xj/X43Zxy5LpfPOZ3QAc2LeYf3hrCSEbwmNGwtsgP1x55ZVs\n3LgRgOPHj7Nt2za2bt3KZZddBsAvfvELPv7xj8frl5WVUV5ezmOPPUZZWRm7du2iqamJyZMnM2nS\nJCZNmqTWzwjm6DvLGLMwxeEWIOVNMi60/nB5ra2T++p2xh9b21tm6X0QDoV7vyccL3e5PGAhFAoA\nFqwlbCPlxhhc7gIAgj0dWGsj54z/EIPHVwxYAt1nsOFw5KzRcuPy4CkYjQV6Olux4TDWRuoAGLcP\nb2FkHUjn2WPx52HBEsbtKcRXOI5QOExX+9FIUSiIDQexoSCewrEUjJ5MKNzDqUM7CYeChIORsrDt\nYbS/nNGTZxPoPM2h7T8nFOgk3NNFKNBNqKeTSVesYOLcN3L6yB52PfIpwqHenYgxMOf6O5l0VSVn\njuxhX1cLvjkzmbFwIYUTyigaP5V/3n6Cf33hCcKmg9NTugn6dhI8/XqMcfGGSz3MmRKOrIkJdtIZ\n7KSjpyP+fVeoC4/xUOwtpthTzFuuKOCRxlbazrpoaZ7Glx/zUbflfygafZBiMwWfazQuMzKGJY/O\nj9zp8u8e+oPDkaRbKVzyDu75xT5gH6GiJbz9c//NqcNNtB8/xNnmI+w5eYT/89OX8fj20/jo13i5\n4Sd9zuD2FVL1n5tweTz8teEnHP/rDjxFo/AWjsJbNIqC0eOYf/0HATj+ciMdp07i9vpwe3y4vD68\nhcWUll8OwNnmI4QC3eByYYwL43Lj9vooGheZSBLoaAcbBuPCxOu4cHsjt40IBXt637vRlnHkPR75\nqA2HEmbYxcoxmOgi1/j7Oo2t6i+85wrmTMpMUjY28dMyw4wxlcB6a+3shGMVwB5gvLW27WLqp/h5\nW5YtW7Zsy5YtFxTn7uPtvOO76yiY8jjQ//elLpbsZW38fd7vOBhS/X/3yuXy2C9lpJZHi1K+9Wzv\n0+LlZpjKbZ9/UpYnhmwGObcD5d3H30HdBz7FguklXIjly5fz7LPPPmutXX4hz3O6b6CE6EB9gpbo\nv34gOWGcd31jzBpgTVLdeUeOHBlSoKGuqXQfngyuhPu5W8B6ML7FgMF27wL32aRyH8a3CIBwYCfG\ndNJHuBBTEOnmsd1/BFd3QqHBhopwFbw++vxGjAn0LQ+PwVXwOsAQDmwFE6aP8DhcvsjVmQ3+PyKt\nHxN7OoTHY3zzgDC2Z2vksDGAC4zB2gkYzyywPcCLkTITKcMYbGgC2ClggrgLD0Q/rXv/oMMBPzY4\nDkwAd1Hs/iMm/iYIByZgQ2PB1Y274HDSb90Q7p6IDY0GE8S42ylgPFdeMp7KyyczuiDy5+syLoo8\nRX2+ir3FFLmLCNogHT0ddAQ7ONtzlo6eDvY0t/Dkn/eyt6UFSw8ubxvG3QEmCCYSWOTnjgFXJ+7C\n/n8z4e7J2NAojKsDV+HRfuWhrjIIF2Hc7bgKjqconxr5v/ecxuXrvy9aqPOSyN+O5xQuX/+dn0Od\nMyJ/e95WXN7W/uUdMwE3xtuMy9t/h+RQxyzA4PKdwHjO9C20hlBn5LbMLt9xjCdpEaV1Eeosj5QX\nHI387hKFPYS6ZkTLD2PcXX2fHvYR7rokUl74GqbP3zzYUCHh7sj2PO7Cg+DqSSovItwdmSHoLtoP\npu+9fWxwFOHA5Gj5XiAU+ZsFCPUQ6iok1DEWwmE84w6DAVdBZJ+1cNdZQu1ugmeLwYbxTWzFGDeu\n4kgvQLjjND2nDaEzBeAKUTDxLLg9uEdFPrBDZ1sJtrkJnfVh3GF8E89iPD5cxeMi5e3N9LR6CXV4\nMJ4QvomduDyF8fOHzpwk0OIl3OnB+EL4SrtweXvLg6eO09PiI9zlwVUQwlvahctXjKtoDFhL8PRx\nAicLsQE3rsIgXn93vNyGigh3ziCTnE4wqVocsQTSkqLsvOtHZ5b1mV1mjNlSVla27EKDnFpSyH//\nfRXQe8fAPhdOJvHbvpdUxqSsFovnPOul/mHJV96JD/uce5DzDRbvYMfNAHGcd71Bz586Jq/LxbTx\nRbhdF9lirICPLoHOQIiDrR20dwfpCYYJhp1rzYtkSsXEzG1a6nSCaSHSKklUAjBAd9eF1h8WxT4P\nS+dMSNfpxSFFPjdzJ6fehVlELp6jI5rW2kb6t0r8QMNw1BcREeeMhCkzG4wxVQmPVwLrYw+MMRVJ\n5YPWFxGRkcHxBBNd01JhjKmKrtDfk7RoshKovoD6IiIyAjg9BgOAtXbdIGX9BusHqy8iIiOD4y0Y\nERHJTUowIiKSFo6u5M80Y8yhcePGTVuwYIHToYiIZI0XXniBU6dOvWatveRCnpdvCeaPwERgt9Ox\nDEEZMLRtCLKXXnPuy7fXC9n5mucAJ6y157e7bFReJZhsZozZbq1d7HQcmaTXnPvy7fVCfr1mjcGI\niEhaKMGIiEhaKMGIiEhaKMGIiEhaKMFkjw3nrpJz9JpzX769Xsij16xZZCIikhZqwYiISFqMiM0u\nRURymTFmvbW2OunYWqCJ6F15oxv75hQlGJERIB8+bKKvEWAJsC15V/Rc/R0YY2qBxSmObbLWNsQe\nG2Oqcu3WI+oiyzLGmH43VzPGrI3eH2eNMWaNE3HJ0EU/bBqttfXRD9XZSTfVy3rGmFpr7bro1yrg\n/QkJJ2d/B8aYigGK1sSSS9QmEu57lSuUYLLIIFdCufjGXBv9qkv8IEoqz5WkmtMfNsaYEqA56fB6\n4JMJj3P1d1BJ5LXEGWMWpqjXEq2bU5RgskQ+XQnl09VunnzY+IHapL/hFqAEcvd3YIypBB5NUeQn\n8voStUWfU5LuuDJJCSZ75MWVUB5e7eb8h421tglYFP03ZiUQ+z/M1d9BibW2LdVxouNMCWKvP/l4\nVlOCyQJ5diWUb1e7efFhY61tjH0f/du8md6Lgpz7HZxjwD5V0om9zuT3c1ZTgskOeXMllIdXu3nz\nYZOgDliR8H+cU7+D6MVRqtcUE79gSlACMMD7PGtpmrIDzvVBmPhHlo9XQgNc7S6KHjpXUs22N2je\nfNhAfPysNvH/mNz7HSwEKhJa20uAkug4Yr21ttEYk/y6/PReROUMJZgMiw5GrzxHnTZrbU0uXQld\nSFJNktNXu/n0YRP9209c+7HQWtuYa7+D5AvC6CzHiqR1PxuSLh5XEhlrzClKMBkW/YM638VUOXEl\ndCFJNelYPlztQh582ETHEf1AQ/Riww+8H4j93+bk7yCaXFYReR+vBTZYa2MXkGuj740KYE+uLbIE\nbXaZVaJ/rNXW2kUJx2qJrIquT/U4W0XfeG3JV7vR71utteMT6lYCNdbaQZPYSJawir2CyOvOiVXs\nEG+9tqYoqo9OQ4/Vy9nfQb5SgskSCVdCi4EHiF4JRcty6o0ZTRgV9M6c8xNJrDXR8pxMqiK5RglG\nRhRd7YrkDiUYERFJC62DERGRtFCCERGRtFCCERGRtFCCERGRtFCCERGRtFCCERGRtFCCERGRtFCC\nERmBEvapEslaWmgpMgIZY1qBhsTdC0SyjVowIiKSFkowIhlijKk1xtjoxqWJxxdGj9dFvyyRWxBU\nRY/b6N5rIllFXWQiGWSM2QRUArOttU3RzT33Ai3W2tnRe/9UELnRWgO990RpTLqNtMiIpwQjkkEp\nEkodUAUsSrpVtCVpB2mRbKM7WopkkLW2zRizAthhjNlB5K6l1Ul37RTJCRqDEcmwaDKpJpJc6nUv\nG8lVSjAizojd9nqho1GIpJESjEiGRRdQxm6B7Y+Ow4jkHCUYkQwyxsRmiG2w1tYTSTJVKaYhtxGZ\nqiyStTSLTCSDjDF7AKy1sxOOrSfSoonPJEuYzryOyLRlNKNMso0SjEiGDDQlOVq2B/ADs6IzzSqJ\ntHQAtgO11tqGjAYscpGUYEREJC00BiMiImmhBCMiImmhBCMiImmhBCMiImmhBCMiImmhBCMiImmh\nBCMiImmhBCMiImmhBCMiImnx/wG/y+/RABNzIQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff1643ea2d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xlabel('xt')\n",
    "plt.ylabel('pdf')\n",
    "\n",
    "plt.plot(xt, model_func(xt), lw=3., label='model')\n",
    "\n",
    "plt.vlines(data, 0, 0.1, 'k', label='data-mean')\n",
    "plt.plot(xt, data_pdf, 'k--', label='data-pdf')\n",
    "\n",
    "plt.plot(xt, data_ref_pdf, 'C2', label='data-ref-pdf')\n",
    "plt.plot(xt, data_ref_pdf_normalised, 'C2--', label='data-ref-pdf-normalised')\n",
    "\n",
    "plt.legend(loc='best')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}