{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Quickstart\n", "\n", "In this simple example we will generate some simulated data, and fit them with 3ML." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's start by generating our dataset:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Configuration read from /home/giacomov/.threeML/threeML_config.yml\n", "Plotter is MatPlotlib\n", "Using Gaussian statistic (equivalent to chi^2) with the provided errors.\n" ] } ], "source": [ "from threeML import *\n", "\n", "# Let's generate some data with y = Powerlaw(x)\n", "\n", "gen_function = Powerlaw()\n", "\n", "\n", "# Generate a dataset using the power law, and a\n", "# constant 30% error\n", "\n", "x = np.logspace(0, 2, 50)\n", "\n", "xyl_generator = XYLike.from_function(\"sim_data\", function = gen_function, \n", " x = x, \n", " yerr = 0.3 * gen_function(x))\n", "\n", "y = xyl_generator.y\n", "y_err = xyl_generator.yerr" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now fit it easily with 3ML:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using Gaussian statistic (equivalent to chi^2) with the provided errors.\n", "Best fit values:\n", "\n" ] }, { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
resultunit
parameter
source.spectrum.main.Powerlaw.K1.02 -0.08 +0.091 / (cm2 keV s)
source.spectrum.main.Powerlaw.index-2.002 +/- 0.031
\n", "
" ], "text/plain": [ " result unit\n", "parameter \n", "source.spectrum.main.Powerlaw.K 1.02 -0.08 +0.09 1 / (cm2 keV s)\n", "source.spectrum.main.Powerlaw.index -2.002 +/- 0.031 " ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Correlation matrix:\n", "\n" ] }, { "data": { "text/html": [ "\n", "\n", "\n", "
1.00-0.86
-0.861.00
" ], "text/plain": [ " 1.00 -0.86\n", "-0.86 1.00" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Values of -log(likelihood) at the minimum:\n", "\n" ] }, { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
-log(likelihood)
data19.354004
total19.354004
\n", "
" ], "text/plain": [ " -log(likelihood)\n", "data 19.354004\n", "total 19.354004" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Values of statistical measures:\n", "\n" ] }, { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
statistical measures
AIC42.963327
BIC46.532054
\n", "
" ], "text/plain": [ " statistical measures\n", "AIC 42.963327\n", "BIC 46.532054" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fit_function = Powerlaw()\n", "\n", "xyl = XYLike(\"data\", x, y, y_err)\n", "\n", "parameters, like_values = xyl.fit(fit_function)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot data and model:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEOCAYAAACetPCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl41dWdx/H394YbwhIWAVFZBAQXEBQSMLhVRREX1Coq\nmA5VEFEL2jrTjrV1m1qZaaeOG6go1KVxQWQU3LAWEXWIkqsWUUQpigShbBEiEshy5o+bQJab5N7k\nrr98Xs/DI/d3f8uJz+/Jl3O+53yPOecQEREJly/RDRARkdSiwCEiIhFR4BARkYgocIiISEQUOERE\nJCIKHCIiEhEFDhERiYgCh4iIRESBQ0REIqLAISIiEWmV6AbEQteuXV2fPn0S3QwRkZQSCAS2Oee6\nNXaeJwNHnz59KCgoSHQzRERSipmtD+c8DVWJiEhEFDhERCQiChwiIhIRT+U4zGwsMLZ///6JboqI\nJKHS0lIKCwspKSlJdFMSKiMjg549e+L3+5t0vXlxI6fs7Gyn5LiI1PbVV1+RmZlJly5dMLNENych\nnHNs376d4uJi+vbtW+M7Mws457Ibu4eGqkSkxSgpKWnRQQPAzOjSpUuzel0KHCLSorTkoFGluf8P\nFDjCcPkjy7n8keWJboaISA19+vRh27ZtzT4nUgocIiISEQWOMBSXlLLxuz0E1hcluikikuK+/vpr\njj76aK688kqOPPJIcnNzefPNNznppJMYMGAAH3zwATt27OCiiy5iyJAh5OTksHLlSgC2b9/O6NGj\nGTRoEFdffTXVJzf95S9/YcSIERx//PFMnTqV8vLymP0MnpqOGwuB9UV8vrmYCge5j+WTd3UOWYd3\nTnSzRCQa/nxe3WODLoIRU2DfD5B3ad3vj78ChubC7u0wb2LN7656JazHrl27lueff565c+cyfPhw\nnn76ad59910WLlzI3XffTa9evRg6dCgvvvgiS5YsYeLEiXz88cfceeednHzyydx222288sorzJkz\nB4DVq1fz3HPP8d577+H3+7n++uvJy8tj4sSJjbSkaRQ4qqnKYzw3deT+Y/nrtlNRGdRLyyrIX7dd\ngUNEmqVv374MHjwYgEGDBjFq1CjMjMGDB/P111+zfv16XnjhBQDOOOMMtm/fzq5du1i2bBkLFiwA\n4LzzzqNz5+Dvor/97W8EAgGGDx8OwJ49ezj44INj1v6kDxxm1g6YBewDljrn8uL5/Jx+XfAZVDjw\nt/KR069LPB8vIrHUUA8hvW3D37frEnYPo7bWrVvv/7vP59v/2efzUVZWFvHCPOccP/3pT5kxY0aT\n2hOphOQ4zGyumW0xs1W1jo8xszVmttbMbq48fDEw3zk3Bbgglu0KlcvIOrwzRx+SSc/ObTRMJSJx\nccopp5CXF/w38tKlS+natSsdOnTg1FNP5emnnwbgtddeo6go+Ltq1KhRzJ8/ny1btgCwY8cO1q8P\nq9BtkyQqOf44MKb6ATNLA2YC5wADgQlmNhDoCWyoPC1m2Z6qXEZh0R5yH8uvETwyM/z06NRGQUNE\n4uKOO+4gEAgwZMgQbr75Zp544gkAbr/9dpYtW8agQYNYsGABvXv3BmDgwIHcddddjB49miFDhnDW\nWWexadOmmLUvYSVHzKwP8LJz7tjKzyOBO5xzZ1d+/nXlqYVAkXPuZTN71jk3vrF7N6XkyMy31vLH\nxWsASKOCm07vzc/OPi6ie4hIclu9ejXHHHNMopuRFEL9v0jFkiM9ONCzgGDA6AEsAC4xs4eARfVd\nbGbXmFmBmRVs3bo14odX5TJ8VOCnlJwPpsN790HZ3ojvJSLiZckUOEJyzu12zl3lnLuuocS4c262\ncy7bOZfdrVujOx/WUZXLOKxzO/ImHEFWv0Pgr7fBg8Phu2+a9TOIiHhJMs2q2gj0qva5Z+WxsDW3\nrHpmhp/MDD9Zxx0Hx82DfyyBvz8LHXoGT/hhB7Q9qEn3FhHximTqcawABphZXzNLB8YDCyO5gXNu\nkXPumo4dO0anRUecARfPBp8vGDQeGAYLroGdB+KZ6liJpBYvbiURqeb+P0jUdNxngOXAUWZWaGaT\nnXNlwDRgMbAamOec+zTC+441s9k7d+6MfqPT0iHrKvj0RXggC5b8HvZ+H/3niEjMZGRksH379hYd\nPKr248jIyGjyPRIyVOWcm1DP8VeBV5tx30XAouzs7ClNvUe9WreHM2+H7KvgzTtg2R/gwyfJbHcf\nxb4o9XBEJKZ69uxJYWEhTZlA4yVVOwA2lad2AKyW45jy5ZdfxvZhGz6AL//KuZ+cwq6SMu47tztZ\ngwc1elmosiYiIskgFafjNlvUcxwN6TWCwBHXVy4a/IHcvDUE5twI2/8R+2eLiCSQpwJHvB0ogGiU\nkk7+N9/DzBPg9Vtgj0qwi4g3eSpwxDQ5HkLVokEAvz+NnPG3BEsuv/9QMIDs+6HONdrbQ0RSnacC\nR1yHqghRAPGYI+CC+2HqOzySNp7L//z34IkbVoBzDdbDEhFJFcm0ADAl7V80WL0A4iHHsqRtcfDv\nX78Hj58LfX9EfudfaW8PEUl5nupxJKVeI+CcP8DmleQU/Bw/ZUDdvT20kFBEUoWnAke8cxxhSfPD\nCVPhho/IOnE0eem/53r/y+RNGq7ehoikJE8FjnjnOBpSJwnepjOMuZu87r+ivOvRZPXtChXl8Pmr\n4KG1NCLifcpxNFOohXxVSfAKB7mP5dfYOXBzqx5sbtUjeOJnL8H8q6DncDoVT+TT8p4E1hepJyIi\nSc1TPY5kcWB9x4EkeJXnpo48EGwGXggXziKw1cfb29rzbdFuch9drtlWIpLUPBU4kiXHUWN9R60k\neA2+NBiaS372vezFTwU+SsvKyF+3LX6NFRGJkKcCR7LkOOqs72hk6CnnyEMx8wEOf6tW5PTrCnuL\n4cOngnmQCGmGlojEkqcCRzLJzPDTo1ObsPIVBwJNW/KmjAxes3IeLJzG+ruHcdf9M+PQYhGR8Chw\nJIk6gSZ7Elz2JBmuhN/uuAXyLoWta8K6l8qaiEgsKXAkKzMYeCE3dZvNXzInwzf5sPiWRi9TWRMR\niTVPTcdt7p7jyajM0lnU/lJ+8pNbYF/ljoPfbYDPXuKKlUMoN3+NKcGhZnRpeq+IRJOnehzJkhyP\niXZdoXOf4N8/XQBv/IZ7tl7D8JL3aiwgDHtGl4hIE3kqcKSyGus7qgmZrzjpRvjJC5RaOv9W9Dt4\n/Hz49mMg8hldIiKRUuBIYg3mK/qfya+6zuLRDtNh6+dQMHf/V5HM6BIRiZSnchxe01i+YufeCp4o\nPYNhF08lq0e74MHCAp476m04cXqd+2m/cxGJBvU4YqS+oadINJSvqNEbeXIVga2VJ659E5bOgAey\n4eNnoKKiWW0QEalNgSOJNZSvqLce1mk3w1WvQ2Z3ePFaePR0WN/wKnKtNBeRSChwJLn68hUNzp46\nfCRcvQR+PBt2b4VvgkFBCwNFJBo8FTiSpchhPDQ6e8rng+Muh2kFkHM9gfVFfLFpJxuLdpP7qBYG\nikjTeSpweHodRwhhzZ5Kbwv+DPLXbacccPgoLSslf9nrUB7cxlY9ERGJhKcCh9Qvp18XzHwYDr9V\nkPPln+ChEwksf0slSkQkIgocLUTV0FaPzm3Jm3oqWVfcCRWl5H+9s95Np0REQtE6jhYkM8NPZoaf\nrD4HAedB/7PIKfwe38rlOFeB3xw5h+rfEiLSMAWOJBfTxXqt0snqcxBHH5JJn53vM7l8PlkLNsKp\n/wonXAf+jNg9W0RSlgJHCotWUMnM8LM942SyLrkI3rgV3rwjWMLkggeg32lReYaIeIfGJeSArgPg\nimdh4kvQuiO/fWVdcGFgiNXnWjQo0nKpx9GChN1D6XcaXPsOX87OD35+5SbYtxtG3QadesWqeSKS\nIpI+cJhZP+A3QEfn3LhEt8eLQgYUq1yW7hy06wYfPw2rFwaLJ570c4pLStlVUkZgfVGNdSQqpCji\nfTEdqjKzuWa2xcxW1To+xszWmNlaM7u5oXs459Y55ybHsp3SADM44zcwvQCOGQvL/kjgnh+zZtNO\nrf0QaaFineN4HBhT/YCZpQEzgXOAgcAEMxtoZoPN7OVafw6OcfukATVWlHfqDZc8BpPfJD89h3KC\nPZLaaz+0Cl3E+2IaOJxzy4AdtQ6PANZW9iT2Ac8CFzrnPnHOnV/rz5ZYtk/qV+8mUr2GkzPuJnyV\nQ1l+ysj5x32wbW3DG0+JiGckYlZVD2BDtc+FlcdCMrMuZvYwMNTMft3AedeYWYGZFWzdujV6rW2h\n6i3bzoFV6L06ZZA3ciNZW/4XZp1A/ut5WoUu0gIk/XRc59x259y1zrkjnHMzGjhvtnMu2zmX3a1b\nt3g20ZMaLNtOcO3HYZ3bknXB9XDDhzD0J+RsfIIM9gXrYdW6RtN3RbwjEYFjI1B9TmfPymPN1pLK\nqsdao2Xbq2t/MIy9j6yfPc7N7V+lTyd/8JquFcFZWSLiKYkIHCuAAWbW18zSgfHAwmjcuKWVVY+1\nsMq2V9d9EK91m8TBnTuQ1bM9/HkMPHkBbP5ESXMRD4n1dNxngOXAUWZWaGaTnXNlwDRgMbAamOec\n+zRKz1OPI2kYjLgGNq8iMOsqvtxUpKS5iEfEdAGgc25CPcdfBV6NwfMWAYuys7OnRPveUlOjC/zS\nWsGIKTB4HPlPPUP5uprTd8PuxYhI0kn65Hgk1ONIvOemjqwZVNp0JufsCWA+fFQcSJpv/0fIGlgi\nkvw8FTiU40hOwUR7Bw7r3C6YND/Y4LEzYc6Z3PrAHM22EkkxngockrxqJNpbd4Cz74Zd3/K77Tdx\nY9HdULS+zjWawiuSnDwVODRUFV11hp2ixeeD4yfA9ADPt89l2N734cHhsOOr6D9LRKLOU4FDQ1Up\nJr0d8zP/hZ93mwNn3g4H9Q0e3/ABlJcltm0iUi9PBQ5JPcUlpawqbkfgsCuCB77bAH8+Fx45he7F\nq7T2QyQJKXBIXIQa9gpZFLFjTxg3l8DuLryxrQvfFu0m99HlCh4iScRTgUM5jtQSspCiGQy8gPxh\n/81e/FTgo7SsjPzPD9TFVNJcJLE8FTiU40gtDRVSzOnfHbPg6+lP85FzdGV5s9Uv08rti3dTRaSa\npN86VryrqpDirpIy7hs/tMZq8pDfbV4Fz+Xyp7RDycucDC7nwBa3IhI3nupxSOppqJBine8OORZ+\nsoCPKvqTsW0lgVmTYOOHcW6xiHiqx2FmY4Gx/fv3T3RTJEYC/mH8rOQ6nHO03lBK3uP/StYvX4H0\ntvvzHjFZeyIi+3mqx6Ech/dVJdQdRqmlkz/wt5DeFioqGL17Ia0rShLdRBHP81TgEO+rk1AfNiz4\nYf17TN41i3u3ToaPn1YBRZEY8tRQlXhLqCGnehPqfU/huvb3c+iuv3PegnvIev+RYD2sPidpCEsk\nytTjkJQTKqEeWF/E4u1dmbtvFLnldxD4ri0snKbSJSIx4Kkeh5LjqSdavYAaiwkrjPyse8g6zgdp\nrWhdUcIFu5+HPcdAm05ReZ5IS+apHoeS4y1XndzHgEOg21EAdPt+NZu/+57APRfD+7OhvDSBLRVJ\nfZ7qcUjLVV/uI7C+iP/dcTgVrjcZ35eR98pdZK14FEbfBQNGawGhSBMocEjKqW94KzPDT2aGv0bu\n48AQllGKn/zB/0HW5lvhvfu5fEkmmIW8X30JdSXaRTw2VCVSW50hrBEj4brlMG4OmNG5fBu88q/w\n/ZbENlQkhShwiKdVDWH17NwmuN/54Z2hVTpkHgLAMfs+gcDjcP8weOceKG14AWFxSan2CJEWT4FD\nPK++eljFJaU8vzeHwEVLoe8p8Lc7g1vYrloQMkCE3D+kESoBL16kwCEtUo0g8HwhgRNnwcSFkNGR\nwEeBkAEi5P4hjVAPRbzIU4FDGzm1bKF2GaxPyCDQ70cw9W3yD7niwHelZeR/tg5oeP+QUJrSQxFJ\nBZ4KHFrHIeGqNwj40sg58jB8BobDTyk5H0yHv/0HWYe0qpsvqRRqSKopPRSRVKDpuOJ5EdW8qv3d\n+f3J+nwQvPMn+PApzm01kXc6nRly/5DaqoJThQuvhyKSKjzV4xCJRFibSA06Gi55FK7+GxzUl4PL\n/xnyXqFyGSFndIl4gHocIuHomQ2TFnNZRRmXpfnhs4XB8u2jf0dgd1c+31xMhYPcx/JrBIlQixJF\nUp0Ch0i4zCDNH/x7yU74+l2YlUP+wXdQ4foBB3IZVYFCK8zFixQ4REJo9Bf+sH+BI8fAW78nZ8WT\nZHALe/Hjb5WmXIZ4nnIcIk3VvhuMvZesnz3BnW2f56T2m4LDVL06gHOJbp1IzChwiDRX94EsOHga\npV0HBYeo8mfBE2Nh098T3TKRmEiJoSozuwg4D+gAzHHOvZHgJokHRDX/UL08e5uD4J+fwiM/guNz\nYdSt+2tjiXhBzHscZjbXzLaY2apax8eY2RozW2tmNzd0D+fci865KcC1wOWxbK9Isw3NhRs+ghOn\nwcrnggUUP3460a0SiRpzMR6LNbNTge+BJ51zx1YeSwO+AM4CCoEVwAQgDZhR6xaTnHNbKq/7E5Dn\nnPuwoWdmZ2e7goKCqP4cIk2yYx389XbIuR4OHwmleyCtNfg0SizJx8wCzrnsxs6L+VCVc26ZmfWp\ndXgEsNY5tw7AzJ4FLnTOzQDOr30PMzPgP4HXGgsaIknloH5w+VMHPr9xK2wMwNl3BwOJSApK1D97\negAbqn0urDxWn+nAmcA4M7s21Almdo2ZFZhZwdatW6PXUpFo6jUCijfBn8fAvImw46tEt0gkYinR\nX3bO3e+cy3LOXeuce7iec2Y757Kdc9ndunWLdxNFwjPkMpgegNN+DV/+FWaOgJXzapyiPTwk2TUa\nOMxsuplFu17CRqBXtc89K481i8qqS0pIbwen3RwMIEMug57Dg8d3b4fysohvp0Aj8RZOj6M7sMLM\n5lXOhLJGr2jcCmCAmfU1s3RgPLCwuTdVWXVJKR0OgwtnwkF9g59fvBYePpnj9oae2KEAIcmi0cDh\nnPstMACYA1wJfGlmd5vZEeE8wMyeAZYDR5lZoZlNds6VAdOAxcBqYJ5z7tMm/gzVn6Ueh6Qm52DY\nRCgr4extTzDg25cIfPxRWJdql0GJt7ByHC44Z3dz5Z8yoDMw38z+EMa1E5xzhzrn/M65ns65OZXH\nX3XOHemcO8I59/tm/AzVn6Ueh6QmMzhmLIGxb3D5vtt4eu+J5D77FYG3FzV4mXYZlEQIJ8dxo5kF\ngD8A7wGDnXPXAVnAJTFuX0TU45BUl79+F2WkUUEapaSTv7dyGGvHV/ywZ0+dnoV2GZRECKfHcRBw\nsXPubOfc8865UgDnXAUh1lwkknockupqbGnrTyPnqB5QUU7giV/yxeadbCz6oUbPItJ90EWiIZwc\nx+3OufX1fLc6+k0S8b76Et0hdw30pZHfawqlpOEwSkvLyP/4k/rPF4mxlFjHES4NVYkXhNrSNmf4\nCLA0DIefMnICN8E/3qr3fJFY8lTg0FCVeFVVz6JH57bkTR5O1tn/An1OAaDfvi9oXVGS4BZKS5IS\nZdVFpNr+5QN6woAbgwdLS5hRcheYDz66DY6boAKKEnN6w0RSmT8DxudBp17w0vUw+0fw1TuAFgxK\n7Hiqx2FmY4Gx/fv3T3RTRJos4g2meo2AyX+FVS/Am3fAE+fDlLdi0jYR8FiPQzkOabHMYPA4mLYC\nLn4UegwDYGjJ+7Cn5qJA9USkuTzV4xDxsrB6Iv42wcKJQLuKYn7+3d1w/73BarzZkyDNH+NWSkvg\nqR6HSKqIR32pzfsyuNQ/i0CHM+G1X8GskbDm9WBdLJFm8FTg0DoOSQXxqC9V9YxPdrUh99tLCIx6\nBnDw7ARa/fBPFUWUZvFU4FCOQ1JBPOpL1XlG2ZFwfT6B0QtYvjU9GLRmv0dg9T+i/mzxPk8FDpFU\nEI/6UiGfkeYnv6TXgYBSXkH+s3fDsv+G0j0N3k8JdalOgUMkzuJRX6q+Z9QMKK3I6dUOlvwOHhwO\nn8xX/kPCollVIgmwfxV4DOtLhXpGVUDZVVLGfeOHknX4ufDVOFh8C7z8C65+tz3Fvo51ZnAVl5Sy\nq6SMwPoi1cQSb/U4lBwXaVydooh9T4Vr3oZJiyn2dQz2OpbcBd99A2izKKnLU4FDyXGRJvKlQfeB\nAPQqWw//9yA8kA1v3kn+F99qsyipwVOBQ0Sab4O/D0wvgEEXwbv3kLPiRlpRDmizKAlS4BCR/fYv\nTPyuHVw8G65eQlb3NB7OeJDDO6XXSeZrtlXLpOS4iEfVV6KkvuNVuYwKB7mP5VcGiSyY9DrzH3qZ\nQ9Lak3VYBrz4MzjpRuh2ZCybL0lMPQ4RARpYmGhGUVrX4N83fwKfvQSzcuDVX9K+YldiGisJpcAh\nIkCYCxN7jYAbPoKsn8KKx7hq011kbllB4Ott8W2sJJSGqkQECLXG40Auo8bwVvtucP7/EOg5kUnP\nrqdsXxrvzlkRs8WMknw81ePQOg6R5qmzxqMB+UWZ7MNPBb7g0Nbq9ZB3KXz7cYPXKaGe+jwVOLSO\nQyR+6gxtdSyCjQGYfRpL/2sc1856ObENlJjxVOAQkfipUw9r5OnB/MeJ0zlpz1Lu3ToJ3v6D6l95\nkAKHiDRZnaGtjI4w+nfc1G02H7ceDtvXBre1hf0BJB6bWElsKXCISNRtaXUo/9P5t3DhrOCBzZ/A\nY2cSeH+Z6l55gAKHiETd/l5FYXHwwO5tULyJ/EWPgqsAml/3Skn2xFHgEJH9nps6st6V5eEKWU33\niNNhWgE52SNIp5Q0yvFbRbPqXmnIK3EUOESkyUIFmnpXoKe3Jeui6RzXPZ0fZ3xE3vB1wdyIc1Be\nVu8zQvUsVOo9sRQ4RCSqGl2B3rYLGw4dTdZFNwQ/f/E6PHQifPnXsJ8Rj33bpX5JHzjM7Bgze9jM\n5pvZdYluj4g0LOytcatmW7XKgIpSyBsHT/0Y/vlZo8+Ix77tUr+YBg4zm2tmW8xsVa3jY8xsjZmt\nNbObG7qHc261c+5a4DLgpFi2V0SiI5IV6BxxOlz/Ppx9d3AB4cMnwZLf7/86VC4jHvu2S/1iXavq\nceBB4MmqA2aWBswEzgIKgRVmthBIA2bUun6Sc26LmV0AXAc8FeP2ikiMhUy+t0rn8pXDaN/xUeb0\nWQIH9QMg8NXWEKXeg0EiHvu2S2gxDRzOuWVm1qfW4RHAWufcOgAzexa40Dk3Azi/nvssBBaa2SvA\n06HOMbNrgGsAevfuHZX2i0h8fe/rAOf81/7P+e+8AS4TSNufy6gKFM2d/SVNl4jquD2ADdU+FwIn\n1HeymZ0GXAy0Bl6t7zzn3GxgNkB2drZqHEhS0y+98OQMOAz/Z99RBvipIKf9VqB/opvV4iV9WXXn\n3FJgaTjnmtlYYGz//nqxRLwga+Tp9H9/KV2//5wbfC+Q9coHsPvXcFqDqVGJsUQEjo1Ar2qfe1Ye\nazbn3CJgUXZ29pRo3E9EmiaaPar2bVpT0uY4sq7MhXfvhcMr58iU7ALzQev2UXuWhCcR03FXAAPM\nrK+ZpQPjgYUJaIeIJJFGV4K3zoRRt0LfU4Kfl86AB4bBh09BRXn8Gioxn477DLAcOMrMCs1ssnOu\nDJgGLAZWA/Occ59G6XnayEkkBTVpJfigi6FTb1g4DWb/CL5aFvuGCgDmPFgrPzs72xUUFCS6GSIS\npplvreWPi9cAkGZw0+ij+NnpYeQqnYNVL8Cbd8DODXDmHXDyL2LZVE8zs4BzLrux85I+OR4JJcdF\nUlPVSvAKF+FKcDMYPA6OPg+Wz4RjLgge3/VtcEV624Ni1+gWLOlLjkRCW8eKpKZmrwT3t4FT/w26\nHBH8/Oov4f6hkP8wlJdGv8EtnKcCh4ikrojKlDTg8keW88tt58Khx8Hr/w6zcmDNa+Cc9vCIEk8F\nDiXHRQTgG38/mPgSXDEPMHhmPKx4LNHN8gxPBQ4NVYnIfmZw5Nlw/XI497/h2EsA6FX6NRRvTmzb\nUpynAoeISB1pfhgxZX+ifOrOe+H+YbDsj1C6B9A2tJHyVODQUJWINKS4pJQbfL8mcMilsOQueCAb\nVj4fnNYrYfNU4NBQlYjUtwK9apHh6p1+cr8eQ+DsF4O9kAVXc3LJW1F5dkvpuXgqcIhIy9bQCvQ6\n282W9IJr3oZL5rA849TgF18tg6L1TX5+o2VTPEKBQ0Q8o6G9yENuN+vzweBxfLfX8W3RDwTm/xEe\nHB5ciV6yK6JnN6lsSopS4BARz2hoL/L6FhlW/cLf8F0JuTuvI9B7Erz7P8EFhAVzobysznNCDUk1\nFLS8RiVHRCQpRKMUe1Vw2FVSxn3jh9ZZTBhqu9kav/DLHfm9rybrzMtg8S3w8i+g0+HQf1Sjz25y\n2ZQU5Kkeh5LjIhLpCvSQvZQew+Cq1+Cni+CIM4Jfrnwetq6p9z7NLpuSQjzV4xARiVS9vRQz6FuZ\nNC8tgTd+A7u3wfDJlO85g817WxNYX1QjQITq0XiRAoeItHiN/sL3Z8B1/wdLZxDIf5tV+05kL2Xk\nPpZfo3fRUvaSV+AQkRYj0l/sVQnw56aOhHZd4bw/ke+Ws+/dbTh8+5PgXu9h1OapHIdWjotIrOUM\nPjq41znuQE7k3Xvh248S3bS48VTgUHJcRGLtQBK8bXCY6mCD5Q/C7NPgf68NbiIVgpdWlWuoSkRa\nvEiHsOrkRKYH4J0/Qf5D8OmLcNKNcNINkN4uBq1NPE/1OEREEiKjI5z1HzBtRbCU+/IHYe/3YV2a\nij0R9ThExFMSOrOpcx+47Ingfh+Z3YNVd1+5CY4dR3FJObtKyupM4U1F6nGIiNSjvqKFz00d2XCA\nyjwk+N9d38IXiwnMvZG1m3ZQWPSDJ+pYKXCIiIQQlaKFHXvAtALy+1xPGT7AKC0tI/+LTftPScWK\nup4KHJqOKyLRErWiheltyRn1YzAfPirwU0ZOv+BuhE0JTsmQE/FU4NB0XBGJloYq7UYqOIW3A4d1\nbkfe1SeQdcQhUFpC/oIHcJW7D9YOTskQIOrjqcAhIhIt0S5auL/4Yv/Dggd2FpJTHqA1+/BRjj/N\nwgpOyTB6nBJpAAAGpklEQVS0pcAhIlKPSCvtRqRrf7J+/jw/6fwZ01u9RJ7vdrJW3gH7fgBCB4hk\n2SxKgUNEJFFapfNJp1F80uNyskaOCpZtb5VRb4BIls2itI5DRCTBdvsy4Zz/hIpy8PnI//wbcOVA\nWo1CismyWZQCh4hIHIS1MNGXBkBO9wr8lFMG+Kkgp/0WoH+jOxzGiwKHiEgChQooWccPpf/SnXT9\n/nNu8L1A1isfwDeXwkUPJ8VmUQocIiJJqH2b1pS0OY6sK3ODZduLN0Faq2CgKS9NaNtSIjluZu3M\nrMDMzk90W0RE4qp1Joy6FS6cGfy8ZTXcOwQCTwRzIgkQ08BhZnPNbIuZrap1fIyZrTGztWZ2cxi3\n+ndgXmxaKSKSAqxyNaKrgE69YNEN8MiPYN3bcW9KrHscjwNjqh8wszRgJnAOMBCYYGYDzWywmb1c\n68/BZnYW8BmwJcZtFRFJft0HwaTFMO7PsHcnPHkBPH9lsBJvnMQ0x+GcW2ZmfWodHgGsdc6tAzCz\nZ4ELnXMzgDpDUWZ2GtCOYJDZY2avOucqYtluEZFEa3AWlhkcezEcdS68/xBUlB3okcRBIpLjPYAN\n1T4XAifUd7Jz7jcAZnYlsK2+oGFm1wDXAPTu3TtabRURSV7+DDj5F3F/bMrMqnLOPd7I97OB2QDZ\n2dnx67OJiGcldFOoCFUVRIxHmxMxq2oj0Kva556Vx5pNZdVFRGIvEYFjBTDAzPqaWTowHlgYjRur\nrLqISOzFejruM8By4CgzKzSzyc65MmAasBhYDcxzzn0ay3aIiEj0xHpW1YR6jr8KvBrt55nZWGBs\n//79o31rEZGkVlxSyq6SMgLri2JejiQlVo6HS0NVItISxXufDk8FDiXHRaQlivc+HZ4KHOpxiEhL\nFM390cORMus4REQktHjv0+GpHoeGqkSkpYrp/ui1eCpwaKhKRCT2PBU4REQk9jwVODRUJSISe54K\nHBqqEhGJPU8FDhERiT0FDhERiYinAodyHCIisWcujvvUxkt2drYrKChIdDNERFKKmQWcc9mNneep\nHoeIiMSeAoeIiEREgUNERCKiwCEiIhHxVODQrCoRkdjzVODQynERkdjzVOAQEZHYU+AQEZGIKHCI\niEhEPLly3Mx2Al82cEpHoL4MeldgW9QbFXsN/UzJ/Kym3qsp14V7TTjnNXSO3q/keVZz7hXptfF6\nvxr6vrnv1+HOuW6NnuWc89wfYHZTvwcKEt3+WPzMyfqspt6rKdeFe0045zXyDun9SpJnNedekV4b\nr/eroe/j9X55dahqUTO/T0Xx/Jmi+aym3qsp14V7TTjnNXSO3q/keVZz7hXptfF6vyJ5Vkx4cqiq\nOcyswIVR5EukKfR+SSzF6/3yao+jOWYnugHiaXq/JJbi8n6pxyEiIhFRj0NERCKiwCEiIhFR4BAR\nkYgocDTCzNqZ2RNm9qiZ5Sa6PeItZtbPzOaY2fxEt0W8x8wuqvzd9ZyZjY7WfVtk4DCzuWa2xcxW\n1To+xszWmNlaM7u58vDFwHzn3BTggrg3VlJOJO+Xc26dc25yYloqqSjC9+vFyt9d1wKXR6sNLTJw\nAI8DY6ofMLM0YCZwDjAQmGBmA4GewIbK08rj2EZJXY8T/vslEqnHifz9+m3l91HRIgOHc24ZsKPW\n4RHA2sp/Ae4DngUuBAoJBg9oof+/JDIRvl8iEYnk/bKg/wJec859GK026BfhAT040LOAYMDoASwA\nLjGzh/BmKQmJj5Dvl5l1MbOHgaFm9uvENE08oL7fX9OBM4FxZnZttB7WKlo38irn3G7gqkS3Q7zJ\nObed4PizSNQ55+4H7o/2fdXjOGAj0Kva556Vx0SiQe+XxFJc3y8FjgNWAAPMrK+ZpQPjgYUJbpN4\nh94viaW4vl8tMnCY2TPAcuAoMys0s8nOuTJgGrAYWA3Mc859msh2SmrS+yWxlAzvl4ociohIRFpk\nj0NERJpOgUNERCKiwCEiIhFR4BARkYgocIiISEQUOEREJCIKHCIiEhEFDhERiYgCh0gcmNlwM1tp\nZhmVu0p+ambHJrpdIk2hleMicWJmdwEZQBug0Dk3I8FNEmkSBQ6ROKksPrcCKAFOdM5pR0lJSRqq\nEomfLkB7IJNgz0MkJanHIRInZraQ4JaefYFDnXPTEtwkkSbRDoAicWBmE4FS59zTZpYG/J+ZneGc\nW5LotolESj0OERGJiHIcIiISEQUOERGJiAKHiIhERIFDREQiosAhIiIRUeAQEZGIKHCIiEhEFDhE\nRCQi/w/WmShJcXoq+wAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xyl.plot(x_scale='log', y_scale='log')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute the goodness of fit using Monte Carlo simulations (NOTE: if you repeat this exercise from the beginning many time, you should find that the quantity \"gof\" is a random number distributed uniformly between 0 and 1. That is the expected result if the model is a good representation of the data)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "cac624743e414b71adf0392598c79467" } }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "The null-hypothesis probability from simulations is 0.84\n" ] } ], "source": [ "gof, all_results, all_like_values = xyl.goodness_of_fit()\n", "\n", "print(\"The null-hypothesis probability from simulations is %.2f\" % gof['data'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The procedure outlined above works for any distribution for the data (Gaussian or Poisson). In this case we are using Gaussian data, thus the log(likelihood) is just half of a $\\chi^2$. We can then also use the $\\chi^2$ test, which give a close result without performing simulations:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The null-hypothesis probability from theory is 0.83\n" ] } ], "source": [ "import scipy.stats\n", "\n", "# Compute the number of degrees of freedom\n", "n_dof = len(xyl.x) - len(fit_function.free_parameters)\n", "\n", "# Get the observed value for chi2 \n", "# (the factor of 2 comes from the fact that the Gaussian log-likelihood is half of a chi2)\n", "obs_chi2 = 2 * like_values['-log(likelihood)']['data']\n", "\n", "theoretical_gof = scipy.stats.chi2(n_dof).sf(obs_chi2)\n", "\n", "print(\"The null-hypothesis probability from theory is %.2f\" % theoretical_gof)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are however many settings where a theoretical answer, such as the one provided by the $\\chi^2$ test, does not exist. A simple example is a fit where data follow the Poisson statistic. In that case, the MC computation can provide the answer." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "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.14" }, "toc": { "colors": { "hover_highlight": "#DAA520", "navigate_num": "#000000", "navigate_text": "#333333", "running_highlight": "#FF0000", "selected_highlight": "#FFD700", "sidebar_border": "#EEEEEE", "wrapper_background": "#FFFFFF" }, "moveMenuLeft": true, "nav_menu": { "height": "21px", "width": "254px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": false, "toc_section_display": "block", "toc_window_display": false, "widenNotebook": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false }, "widgets": { "state": { "b07a76dd4e07418b9722dc113ae5d245": { "views": [ { "cell_index": 8 } ] } }, "version": "1.2.0" } }, "nbformat": 4, "nbformat_minor": 1 }