Last active
June 7, 2018 14:24
-
-
Save taldcroft/5014170 to your computer and use it in GitHub Desktop.
Investigating `scipy.optimize.curve_fit` covariance output
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "metadata": { | |
| "name": "curve_fit" | |
| }, | |
| "nbformat": 3, | |
| "nbformat_minor": 0, | |
| "worksheets": [ | |
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Investigating `scipy.optimize.curve_fit` covariance output\n", | |
| "-----------------------------------------------------------\n", | |
| "\n", | |
| "In Aug 2011 there was a thread [Unexpected covariance matrix from scipy.optimize.curve_fit](http://mail.scipy.org/pipermail/scipy-user/2011-August/030412.html) where Christoph Deil reported that \"scipy.optimize.curve_fit returns parameter errors that don't scale with sigma, the standard deviation of ydata, as I expected.\" Today I independently came to the same conclusion.\n", | |
| "\n", | |
| "This thread generated some discussion but seemingly no agreement that the covariance output of `curve_fit` is not what would be expected. I think the discussion wasn't as focused as possible because the example was too complicated. With that I provide here about the simplest possible example, which is fitting a constant to a constant dataset, aka computing the mean and error on the mean. Since we know the answers we can compare the output of `curve_fit`." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "import numpy as np\n", | |
| "from scipy.optimize import curve_fit" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 1 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Generate a normally distributed dataset of 100 points with $\\mu = 0.0$ and $\\sigma_y = 1.0$" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "x = np.arange(100, dtype=np.float)\n", | |
| "y = np.zeros_like(x)\n", | |
| "yn = y + np.random.normal(size=len(y))\n" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 2 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Define a constant function" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "def const(x, a):\n", | |
| " return a * np.ones_like(x)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 3 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Fit a constant to a set of normally distributed points\n", | |
| "\n", | |
| "This is equivalent to determining the mean and uncertainty on the mean.\n", | |
| "We know the answer analytically:\n", | |
| "\n", | |
| "- $\\mu = 0.0$\n", | |
| "- $\\sigma_\\mu = \\sigma_y / N^{1/2}$\n", | |
| "\n", | |
| "If no `sigma` parameter is provided to `curve_fit` then we get the expected results:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "popt, pcov = curve_fit(const, x, yn)\n", | |
| "print 'Best fit constant =', popt[0]\n", | |
| "print 'Uncertainty on constant =', np.sqrt(pcov[0, 0])\n", | |
| "plot(x, yn, '.')\n", | |
| "plot(x, const(x, popt[0]))\n", | |
| "grid()" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "Best fit constant = -0.114460255805\n", | |
| "Uncertainty on constant = 0.0990467843199\n" | |
| ] | |
| }, | |
| { | |
| "output_type": "display_data", | |
| "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD9CAYAAACoXlzKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGOJJREFUeJzt3X1snWX9x/HPHlAyJBRE51xrWsZmR9fRjtRgZPZM1iWt\nTnFYxhIlHUlLApE0kBD/0KyAGRB/iQSID3PJFkzqAj6EiLNKsp4tFAdq1nXOkBFYlTZrzNymzE3A\n7f79cei5266P5366Ht6vZBmndKdXv+c+33Pfn+s615kXBEEgAIAz5mc9AABAvGjsAOAYGjsAOIbG\nDgCOobEDgGNo7ADgmEiNfWRkRLlcTrW1tVqxYoW6urpiGhYAoFTzoqxjP3nypE6cOKHa2lqdO3dO\na9as0c6dO3XLLbfEOUYAwBwsjPKPr732Wl177bWSpEWLFmnVqlUaGRmJZWAAgNJEOmMfa3BwULfc\ncouOHDmiq6++unDn8+bFcdcA4J0orTmWydNz587pjjvu0FNPPVVs6qOCIOBPEGjbtm2Zj8GUP9TC\n7Vq0twdqbAzU3Bzo9Gm/a1Hqn6giN/YLFy5o8+bNuvPOO7Vp06bIA3LV4ODgnL6/o0PK5aSWFunM\nmUSGlJm51sJlLtbi2DFp/37pt78tHMez5WItshK5sbe3t6uqqkoPPPBAHOPBB0p9cgBZW7So8HdD\ng7RjR7Zj8VWkxt7X16fdu3frwIEDqq+vV319vV588cW4xuaUtra2OX2/y0+OudbCZS7Wortbam2V\nfv97qaxs9v/OxVpkJbbJ00nvfN68WPIiH505UzhT37Fjbk8OAPaL2jt552lK8vn8nL6/rEx67jk3\nm/pca+Ey02qR5dyOabWwGY0dQBFzO24gigFQ1NJSaOoNDXPPyBGfqL2Txg6giLkdM5CxW4L8MEQt\nQqbVIsu5HdNqYTMaOwA4higGAAxDFAMAGIfGnhLywxC1CFGLELWID40dABxDxg4AhiFjBwCMQ2NP\nCflhiFqEqEWIWsSHxm45lz+QA0BpyNgtl8sVNm2SCntgP/dcpsMBEAMyds+5/IEcAEpDY09JUvlh\nqZ9WkyWy1BC1CFGL+CzMegCIZnTTJgAYRcYOGKijo/ChF4sWFa7KbLkaQzzI2AEH8UlGiILGnhLy\nwxC1CE1VCx8nxTku4mNUY2dNNlBg46Q4zGFUxs6abABwLGP38fITAOJmVGN3+fKT/DBELULUIkQt\n4mPUOnbWZANAdEZl7AAABzJ2VsIAPA8Qr8wbuy9vxCA/DFGL0GgtfHkeTIfjIj6ZN3ZWwsBkaZ1J\n8zwIcfUSXeYZ+5kzhQdyxw73VsLAfmm9t4LnQYj3s0TP2DNfFcNKGJgsrTNpngchrl6iyzyK8QX5\nYcimWiT93gqbapG00Vq4/H6WtGR+xg6kodRtcDmTTt9sa87WxlPLPGMH0kBu6x6XH9NM17Hffffd\nWrx4saqqqqLcDZA4clv38JhOLVJj37p1q3p6euIai9PIUkNZ1MLU3JbjIjTXWpj6mJogUsa+du1a\nDQ4OxjQUIDlk5e4hi59a4pOnbW1tqqyslCSVlZWprq5OuVxOUvgK7cPtXC5n1Hi4bc7tUaaMZ7Lb\nHR3Sa6/l9eEPS7/7XU5lZfH/vNGvxT3+Y8dyH2Txed12m5TPx1+fqLfz+bx2794tScV+GUXkydPB\nwUGtW7dOx48fv/TOmTwFnGDzRGVLS2GrhoYGe2Ib6zcB88XEs7PJ+PJW6tnUwhe21CKNicqkauFj\nFk9jN4hvG0H58kLmApub42gWb9u4o4gUxbS2tqqvr08nT57U4sWL9eCDD6qzszO8c6KYObHxkjEK\nmy/vgSRF7Z28Qckgrm4ENdWqBN9eyIDZImO3xGzyQ1cvGSdGTKO1MPHyPu14yJaMPQ3UIj40diRu\nqok3E1/IfJvnSJMPcyqm/I5EMUicTRET8VByfJhTiet3tDKKMeVVDekw8cx8KibGQ67wYW8XU37H\nTBq7j5e75Ichk9f0p/0i5NNxMdOLpi21mO7YNOXEIJP92E15VYO5Rl/8pcITycXLdt+4sl/PdMem\nKb9jJhm7TZkrskHWDVOlcWyyjh1O4sXfbD7umDgqjWPTyslTH9mSH6bB9TX9c5kfsPW4SGKezOS5\nl7FsODZp7EDMfFgcEPc8WUeH1Nk5c8P2obZxIIoBYubD/EDcccRs139PV1uX4iEy9hK5dBDALMwP\nzN1sXwynq61Lb4AiYy9R2pd0tmapSXC9FnPJYF2vxWx1d0uNjfkZr3Cmq22Sy6iny/ZNyP0n8rax\ns5YeMEdZmdTVFe0Kp5Q3B822KU93Imhi7u9tFMPlMoA4sv0k5lTI2AGgRHFk+0mcJJKxW2I0SzUx\nj0sbuXKIWoSyqMVs45vpsn0T17VnsleMz9gDBTCHKXu7xI0oJmU+rHEGEA0Zu2VsmLRljT/iwrFU\nGjJ2S4zmhybmcRMlvXyLXDnkei3mciy5Xos00dhxCdb4Iy4cS5NLehEFUQwuYUNcBDtwLE1upvXz\nZOwAnOZiTj/TIgoy9gz5sO92EqhFiFqEpqqFiW/Zjyrpz0ZlHXsErEmHjWw7A3Yxp096/TxRTASs\nSYeNbNve1secnow9Qz4ecLAfJyTmI2PPEPtul4ZahEzeHyVtHBfxMTpjty0LtBV19our+6MgZHQU\nY1sWaCvqDJjF6SjGxdlwE1Fns7HVM+bK6MZuahZYCpPzw7TrbHIt0jabWri4jnsyHBfxMTpjJwtM\nB3U2WxxXVMyj+MXojB1Igm1NLo5ltcyj2CVq7zT6jB1Igm3vGI7jiop5FLNNPNmIKnLG3tPTo1Wr\nVmnlypV64oknoo/IUeSHoaxrYVKTS6sWNsxXZX1cZCnueZRIjf3dd99Ve3u79u7dq4GBAXV3d+vQ\noUPRRwUkyIYmFzcbPuDFZ3GfbETK2A8cOKCuri7t27dPkvTwww9rwYIF+va3v124czJ2pMi27BwY\nNXEeJdOMfWhoSEuXLi3eLi8v1x//+Mdx39PW1qbKykpJUllZmerq6pTL5SSFl14+3/6//5POns1p\n0SLp3nvz+shHzBqfTbdfey2vw4clKaeOjkI9TRoft7k91e3+/rwWLdqtzk4V+2UkQQTd3d3B17/+\n9eLtnTt3Bvfcc0/xdsS7d0pvb++kX29sDAKp8Ke1NdUhZWaqWkTV3FyoY0NDEJw+nciPiF1StbAR\ntQhF7Z2RMvby8nINDw8Xbw8NDamioiLiS41fTJrIs52P2TkwmUgZ+3//+19df/316uvr05IlS9TQ\n0KBdu3ZpzZo1hTs3LGM3MYNl618AE2WasV9++eX6yU9+oi9+8Yu6cOGC7rrrrmJTN5GJ65d51yeA\nuCX+zlN1JXXvlhmUVJnxGEwxKGoxalDUYtSgvK9FsK3QjvkEpTnIMvbI5/PF2XDfUYvQXGvh8tYA\nHBchGjvgET7Wzg80dhghzYlpEyfB08Jkux+c/qCN2bLhgwhG34zgqrnsdRG1FlntT57EcTbXWri8\nNYDrz5E0OdHYffkggrlK8wUvzfX4Wa395ziDLZyIYsgdJ5fmRFuaEUFWcYSLx5nPsZbJjM/Ym5uD\nxA8YcsfJudiIkjDb5ubicebyKhubGZ+xp3HZakPumEV+aOpb7E3LUmcbsSRxnGVdC5O2tMi6Fi5J\nvLGbcMD4ZGyuLpn/gmcCk5pb2kx98Uc0iUcxp08HXh4wWWWXXFrPnYsRC+xm/Gee+vpEyWpfGp/P\nPkvFfj1wjRPLHcfKak37xJ87scGmlR/acGlNlhqiFiFqEZ/Ez9jTltWZ8sSf292dzeU9Z58AnFjH\nPlZWS/xYWgggLsavY0+7sWc1EZb0z+WNJIA/jF/Hnras1rTP9HNt3R8lCXOthQ17AZUq6nHhUm3I\n2OPjXGN3lc+rXVx6UYsbtcFknItiXOXzWmvmL6ZGbdxExg7n+fyiNhNq4yYy9pSVmmmSH4Zc3IM8\nq+PChtrMFs+R+KTS2F2a4CHTxGQ4LmCSVKKYLPcviXuZIJkmJsNxgThZEcVkuaIj6pnUxKsNG96y\nPxWXrpxMM9vjgscge148BkGCRu/+9OkgaG0t/J225uYgkIKgoWH6n9/eHgSNjYXvH/t9jY2Ffy8V\nfodS9fb2lv6PYxLX7xKVCbXIysTHwOdaTJRWLUx5HkwnamtO5Yw9ywme2Z5JTXVm79L6cZd+l1G2\nnX25+BjYxovHIKYXmEklfPexmurMPsurjbi59LuMsuHsaywXHwPb2PAYRO2drGP/AOuB7cSkpb3Y\n/2hqVkye2iDpuIg1uqE4a2HzZLbk93ExMf70uRZxo7EjVR0dUmdnfJm4S2/QiZMNcw9eZN0ZIYpJ\nAZecIT6TNR021Nm3+HMufYAoxgK8KzHEWVo6bKizb1dbafYBGnsKCk+yvNFPsrR0d0uNjXlrM/G4\nJZUr2zj3YFrGHnecleaLLY09BYVmZteTLCllZVJXF3VImm9nw0mI+ww7zRdbMnYAmESWS2nZjx0A\nEpDl5C6Tp5YwLT/MErUIUYuQabWwOc4qubE///zzqqmp0YIFC7R/dF0VACBzJUcxr7/+uubPn697\n7rlHDz/8sD7/+c9feudEMQAwZ1F758JS/2F1dXXJPxQAkJySG/tstbW1qbKyUpJUVlamuro65XI5\nSWGm5sPtsfmhCePJ8vbo10wZT5a3+/v71dnZacx4srz95JNPet0fdu/eLUnFfhnFtFFMU1OTRkZG\nLvn69u3btXHjRknSunXriGJmIZ/PFx9Q31GLELUIUYtQ5ssdaezA1NgnyB4mPVZGLHekeWMsG3YW\nTAv7BNnDpccq0nLHiooKHTx4UK2trfrc5z4X57icMzZfjsr0xjnTEyTOWmRlto/BTPuDuFCLuGRd\nCxs2Tputkht7a2ur3n77bZ0/f17/+Mc/1NfXF+e4MA3TzyxceoJMZbaPgY2bcfnKpceKLQUsZPrH\nwfmwz7bpjwHslvnk6bR3TmNPhKuN06TJq5m4+hjADEZMnromiQw7zvzQ5j0spKlrYXrENFZcj0HW\nubJJqEV8aOyTsKnBuMSHbB5IA1HMJHzIT02MPYg3gAIy9gT40GBs+LBjwFdk7AlIIsM2LT/MMvYw\nrRZZohYhahEfGrunXFqzC2A8ohgAMAxRDABgHBp7SsgPQ9QiRC1C1CI+NHYAcAwZOwAYhowdADAO\njT0l5IchahGiFiFqER8aOwA4xpqM3cS9TQAgCd5k7Oy4CACzY01jt31LV/LDUFq1MP2zYSWOi7F8\n+lzgpFnT2NnbJHu2PVm4yvOX74+9NRk7smfbVr8+7KuPydn+2HuTsSN7tsVhXOX5y/fHnsaeEhey\n1LieLGnVwobPhnXhuIgLnwscn4VZDwD2GH2yAGOxFHlyWdaFjB1AJLbNvaQlSl3I2AFkyra5l7Rk\nWRdrG7ttS+/IUkPUIuRCLWybe0lLlhO41mbso+tUpUKT5/IvfWSrkJh7mUqWdbE2Y7d9naotpmve\nZKtAMrzN2H1fp5qW6d7BR7YKmMnaxm7bOlVb88PpmnepL6621iIJ1CJELeJjbWNHOqZr3ra9uAK+\nsDZjBwBXeZuxA0iWbUuKEaKxp4T8MEQtQibXIu2tb02uhW1Kbuzf+ta3VFNTo5qaGm3YsEEjIyNx\njgtAxlj1ZK+SM/b9+/dr7dq1mj9/vr7zne/oxIkT2rlz5/g7J2MHrHXmTOFMfccOJsjTllnG3tjY\nqPnzC//8M5/5jE6cOFHyIACYh1VP9oplS4Ef/ehHuv322yf9f21tbaqsrJQklZWVqa6uTrlcTlKY\nqflwe2x+aMJ4srw9+jVTxpPl7f7+fnV2dhoznixvP/nkk173h927d0tSsV9GMW0U09TUNGl2vn37\ndm3cuFGS9Pjjj+vVV1/Vr371q0vvnCimKJ/PFx9Q31GLELUIUYtQ1N4ZaR37z372Mz399NPat2+f\nLr/88tgHBwA+ito7S45iXnrpJX33u99VPp+ftKkDALJR8uTpfffdp3feeUcbNmxQfX29OtJY6Gqx\nsfmy76hFiFqEqEV8Sj5jP3bsWJzjAADEhL1iAMAw7BUDABiHxp4S8sMQtSjo6JDq6vJssvUBjov4\n0NiBjBw7Jh0+nN4mW/AHGTuQET63F1PJ9A1KM945jR2YEptsYSpMnlqC/DBELQrKyqR7783T1D/A\ncREfGjsAOIYoBgAMQxQDABiHxp4S8sMQtQhRixC1iA+NHQAcQ8YOAIYhYwcAjENjTwn5YYhahKhF\niFrEh8YOAI4hYwcAw5CxA5BU2HcmlxPbAIPGnhbywxC1CMVZi2PHpP377d0GmOMiPjR2wBGLFhX+\nbmgo7BgJf5GxA45gG2B3sB87ADiGyVNLkB+GqEWIWoSoRXxo7ADgGKIYADAMUQwAYBwae0rID0PU\nIkQtQtQiPjR2AHAMGTsAGIaMHQAwDo09JeSHIWoRohYhahEfGjsAOIaMHQAMQ8YOABiHxp4S8sMQ\ntQhRixC1iA+NPSX9/f1ZD8EY1CJELULUIj4lN/Yf/vCHuvHGG7V69WrV19frlVdeiXNczjnDZ5UV\nUYsQtQhRi/iU3NjvuusuHT58WAMDA3rkkUf0zW9+M85xAQBKVHJjv+KKK4r/ffbsWS1dujSWAblq\ncHAw6yEYg1qEqEWIWsQn0nLHH//4x/re976ns2fPqq+vT8uWLRt/5/PmRR4gAPgosY/Ga2pq0sjI\nyCVf3759uzZu3Fi8vWvXLj377LPq7e0teSAAgHjE8galc+fO6aMf/ajOnz8fx5gAABGUnLG/+eab\nxf9+4YUXVFtbG8uAAADRLCz1Hz722GN69dVXdfHiRX384x/XT3/60zjHBQAoUcln7Dt37tSRI0d0\n9OhR9fb26tOf/vS4/9/T06NVq1Zp5cqVeuKJJyIP1CYjIyPK5XKqra3VihUr1NXVJUk6deqUmpqa\ntHr1am3YsMGrdbsXL15UQ0OD1q1bJ8nfWpw+fVqbNm3S6tWrdcMNN2hgYMDbWjz00ENavny5qqur\n9dWvflX//ve/dfz4cX32s5/VqlWrdOedd+r999/PepiJuPvuu7V48WJVVVUVvzbdcXD//ffrhhtu\n0Jo1a3To0KEZ7z+Rd56+++67am9v1969ezUwMKDu7u5ZDcYVCxcu1NNPP60jR46ov79fe/bs0csv\nv6xt27Zp/fr1GhgY0Be+8AVt27Yt66Gm5plnntHy5cuLK6V8rUVHR4daWlo0MDCgI0eO6LrrrvOy\nFocOHdLPf/5zHT16VK+//ro+9KEPadeuXbr//vvV2dmpv/zlLyorK9MzzzyT9VATsXXrVvX09Iz7\n2lTHwS9+8Qu98cYb+utf/6of/OAHamtrm/kHBAnYv39/sG7duuLtrq6u4NFHH03iR1nh9ttvD557\n7rnguuuuC/7+978HQRAEg4ODwbJlyzIeWTqGh4eD9evXB/v27QtyuVwQBIGXtTh58mSwZMmSS77u\nYy1OnDgRLF++PDh16lTw/vvvB1/60peCF198MbjyyiuD//3vf0EQBEE+nw9uvfXWjEeanOPHjweV\nlZXF21MdB1u3bg2effbZ4vdVVlYGQ0ND0953ImfsQ0ND496wVF5erqGhoSR+lPEGBwd18OBB3Xrr\nrePqsnTpUm9q8uCDD+rxxx/X/Pnh4eZjLd544w0tXrxYW7ZsUU1Njb7xjW/onXfe8bIWn/jEJ/TQ\nQw/pU5/6lD75yU/q6quv1po1a3TVVVdpwYIFkvypxaipjoPh4eE599NEGjtvTCo4d+6cWltb9dRT\nT+maa67JejiZ6Onp0VVXXaWbbrrJ+735L168qMOHD+vee+/V0aNHdcUVV+jRRx/NeliZePPNN/XY\nY4/prbfe0vDwsE6ePKmXXnop62E5I5HGXl5eruHh4eLtoaEhVVRUJPGjjHXhwgVt3rxZW7Zs0aZN\nmySNf6UdHh5WeXl5lkNMxSuvvKLf/OY3qqqq0pYtW3Tw4EF95Stf8bIWFRUVuuaaa7R27VpJ0m23\n3abDhw+roqLCu1q89tpramho0Mc+9jFddtll+vKXv6wDBw7oX//6ly5cuCDJn1qMmuo5MfEMfTZ1\nSaSxNzQ06NixY/rb3/6m9957T7/85S/V3NycxI8yVnt7u6qqqvTAAw8Uv9bS0qLu7m5JUnd3t1pa\nWrIaXmoeeeQRvf322zp+/Lj27Nmjm2++WS+88IKXtaioqFBFRYUGBgYkSb29vaqurlZzc7N3tbj+\n+uv1pz/9Sf/5z38UBIH27dun6upqNTY26vnnn5fkTy1GTfWcaGlp0Z49eyRJf/jDH3TllVfOvDdX\n7DMCH9i7d29QU1MTVFdXB9u3b0/qxxjp5ZdfDubNmxfceOONQV1dXVBXVxf8+te/Dv75z38G69ev\nD2pra4Ompqbg9OnTWQ81Vb29vcVJdV9r0d/fH9x0003BypUrg+bm5uDUqVPe1mLbtm3BsmXLghUr\nVgSbN28Ozp8/H7z11lvBzTffHNTU1AR33HFH8N5772U9zER87WtfC5YsWRJcdtllQXl5efD9739/\n2uPgvvvuC1auXBnU19cHf/7zn2e8/0Q/8xQAkD4+QQkAHENjBwDH0NgBwDE0dgBwDI0dABxDYwcA\nx/w/TXO1eMNTrhAAAAAASUVORK5CYII=\n" | |
| } | |
| ], | |
| "prompt_number": 15 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "But now setting `sigma=100` doesn't change the output covariance at all:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "sigma = 100\n", | |
| "popt, pcov = curve_fit(const, x, yn, sigma=sigma)\n", | |
| "print 'Best fit constant =', popt[0]\n", | |
| "print 'Uncertainty on constant =', np.sqrt(pcov[0, 0])" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "Best fit constant = -0.114460248457\n", | |
| "Uncertainty on constant = 0.099046785584\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 22 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Try to fix the covariance matrix using the method suggested by @cdeil.\n", | |
| "\n", | |
| "This gives the expected result." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "chi = (yn - const(x, *popt)) / sigma\n", | |
| "chi2 = (chi ** 2).sum()\n", | |
| "dof = len(x) - len(popt)\n", | |
| "factor = (chi2 / dof)\n", | |
| "pcov_sigma = pcov / factor\n", | |
| "print chi2\n", | |
| "print dof\n", | |
| "print factor\n", | |
| "print np.sqrt(pcov_sigma[0, 0])" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "0.00971216310043\n", | |
| "99\n", | |
| "9.81026575801e-05\n", | |
| "9.99999998803\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 25 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Use Sherpa to do the fit and compute covariance\n", | |
| "\n", | |
| "[Sherpa](http://cxc.harvard.edu/sherpa/) is a modeling and fitting package used in astronomy. There is a relatively steep learning curve and installation isn't so easy, but the optimization methods are robust." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "import sherpa.ui as ui # sherpa user interface" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 6 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "**Load dataset 1 with the same `x` and `yn` dataset, with `sigma=100`**" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "ui.load_arrays(1, x, yn, 100.0 * np.ones_like(yn))" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 7 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "**Set the model for dataset 1 to a constant model named `c1`**" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "ui.set_model(1, ui.const1d.c1)\n", | |
| "c1.c0.min = -100 # set the parameter bounds to -100 to +100\n", | |
| "c1.c0.max = 100" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 8 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "**Fit the dataset using Levenberg-Marquart minimization and $\\chi^2$ statistics**\n", | |
| "\n", | |
| "The answer is exactly the same as `curve_fit`." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "ui.fit()" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "Dataset = 1\n", | |
| "Method = levmar\n", | |
| "Statistic = chi2gehrels\n", | |
| "Initial fit statistic = 0.0221324\n", | |
| "Final fit statistic = 0.00971216 at function evaluation 4\n", | |
| "Data points = 100\n", | |
| "Degrees of freedom = 99\n", | |
| "Probability [Q-value] = 1\n", | |
| "Reduced statistic = 9.81027e-05\n", | |
| "Change in statistic = 0.0124202\n", | |
| " c1.c0 -0.11446 \n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 9 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "ui.plot_fit() # this shows the large error bars" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "display_data", | |
| "png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEICAYAAABxiqLiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHX6B/DPCAhekJESJS6CigJiCoJaecFrK5aal7ys\nlpritq3+1Lasn78C0rxtW1Zuaym5Zc2apJXlJa1Ek7xleKXCFEhISZQyFQXh+/vj6XBmGEA8AoP5\neb9e85qZc/me5zzne77POTODmpRSCkRERNepnqMDICKimxMLCBERGcICQkREhrCAEBGRISwgRERk\nCAsIEREZUusFZNKkSWjevDkCAwNLp02YMAFBQUEIDw9HeHg4jhw5AgBQSmH69OkIDQ1FREQEUlNT\naztcIiKqQK0XkIkTJ2Lz5s0200wmExITE5GamorU1FSEhYUBANatW4djx44hLS0Nr732GiZMmFDb\n4RIRUQVqvYD06NEDTZs2tZte3t8zbty4EWPHjgUAdOvWDefPn0dOTk6Nx0hERNfm7OgANI899hhK\nSkrQv39/LF68GK6ursjOzoaPj0/pMr6+vnbTALmDISKi63cj/xhJnfgSfdGiRThy5AhSU1ORm5uL\nefPmXXcbSik+lEJcXJzDY6grD+aCuWAuKn/cqDpRQLy8vAAArq6uGD9+PPbu3QtAv+PQ5OTkwNfX\n1yExEhGRrTpRQM6ePQsAKCkpwfr16xEaGgoAiImJwerVqwEAu3btgru7u93HV0RE5Bi1/h3IyJEj\nkZKSgry8PPj5+WHWrFnYsWMHsrKycP78eYSHh2Px4sUAgOHDh2Pbtm0IDQ2Fm5sbVq5cWdvh3nSi\no6MdHUKdwVzomAsdc1F9TKo6PghzMJPJVC2f5xER3UpudOysEx9hERHRzYcFhIiIDGEBISIiQ1hA\niIjIEBYQIiIyhAWEiIgMYQEhIiJDWECIiMgQFhAiIjKEBYSIiAxhASEiIkNYQIiIyBAWECIiMoQF\nhIiIDGEBISIiQ1hAiIjIEBYQIiIyhAWEiIgMYQEhIiJDWECIiMgQFhAiIjKEBYSIiAxhASEiIkNY\nQIiIyBAWECIiMoQFhIiIDGEBISIiQ1hAiIjIkFovIJMmTULz5s0RGBhYOu3cuXPo378/7rzzTgwY\nMAC//PJL6bzp06cjNDQUERERSE1Nre1wiYioAiallKrNDX755Zdo3Lgxhg0bhoyMDADAtGnT4Ovr\ni9mzZ2PhwoU4deoUXn75ZaxduxYrVqzApk2bsHv3bkydOhUHDx603wmTCbW8Gze95GR5lBUdLQ+i\n6lYdfY79tnrd6NhZ6wUEADIzM9G7d+/SAtK6dWskJyfDz88PWVlZ6Nu3L3744QdMmjQJvXv3xvjx\n4wEAgYGB2LlzJ3x8fGx3wioJZTtYZqY8AgKAy5eBn34C6tUDWrSQh9lcM53PUR29su1q88tq1Ajw\n9gbGjau8vXPngEuXAF9f4/tR3Xkxsr/XmnejcSgFmEzlt2e9nNY3AcDfH/DxAVxdrx2DkRxar5Oe\nDpw8CTg5yTYbNpTn2hyEL1wAnn8eWLDAPj7rvAQEyAOwj2/3bsBiATw9r72eto2yy/n4AIGBgLNz\n5cu5ugJXrly7bWvX088q2n83Nxm3tO1q861jKK+9qvpDFBBXV1cUFBSgXr16uHr1Kho3bozLly/j\n3nvvxezZs9GnTx8AQI8ePfDCCy+ga9eutjthMqFlyzgAUhA6dYpGQEA0IiKAf/8b2LRJluvTB/i/\n/5PnJ54APv8cGDwYKCkBvv0WaN++8gMLAMePy8lnfTAB+4Pr4SGxmEzSdoMGtgccqLzDlm3v2DEp\nfE5OVe9UubnA+fNAUBBQXAzk5wO33y7b/cc/gOeeAzp3BqZNA/btA/70J/v4rHPx5JPAZ59Jzqzj\nDQiQHP74ox5TRSfYyZOAu7vkJjBQ8vPQQ7gm6+NQWAi4uMi61jEeOiT70ru3/fq5uUBoqMwva8AA\nIDhYBiLrY9CyJXD1KpCTY3+8W7QA7rgDaNLEPk+RkbKtVq2AtDRZp2FD+7ysWwcMGQL07Qu89x6w\nbBnQq5fMv3pVBjVAcmX1qW6p6GjgX/8C3npL2q9qYZk1Sy4CZs2S4z5+PDB6tMz76iugoED6mbOz\nDJpOTlUfkP38gDNnpH8GBABZWdJG/fq2bYSHA488AuTl2cdrsUgcS5fKsR42DPjkE/vl1qwB3n9f\nngEgLk4e9eoBf/6zTAsKss/FG28AK1ZIzu6/X6Z37iznysmT+jiwciXw4ovAbbcBjz8ux/vxx2UM\nWL4cWLjQPqbJkyUX3bvbb3f7dmDiRLlYqOxY7dghx6NtW2DtWjmvZs6Uojt7thxzAPjLX+TY3XNP\n5RfOWVnAxYuyv/XqJePCheTS/vjWWwk3VECcDa9Zx6xZE48uXYCkJOD0aRkocnOBKVPKX75dOxlQ\n4+MlubfdJh0nOVk6R3a2DFBlT5yEBLnCBIBHHwVmzACGDpUOv2kT8J//yLxGjYCff5bnyZOBbt3k\n+eBB4IUXgNatZVvHj0v7Tk4Sw5AhwAMPyMCzZIneKVq1ksG7VSu5amvQQLadliYn4q5dstxzz8l7\nHx9g/nzpdPHxQEaGFM7fazaeegooKpLXbdrISRcfD5w9K51/1CjZ9ptvAl5eQOPGsmx4uCx36ZKc\nCNoAMGSIdPIhQ+REeecdiSE5WbapfeWVmSnxjR0LfPyx7MuJE/q8qlzxPfCA5LJRI+Ddd2UA8fIC\nmjeXwS8+XpZ9+mmJyctLTrQPPpD9u3JFjm/r1tLe6dNyTO68U459x46yXmamzNdimj1b8hYcLP0r\nJQUYOFD28bPP9KvYpk1le/37S/527pTnggKgWTM5JoAc399+0/ukl5ceu8kEPPusPB86JIXKbJYL\ngB9/BBYvluViYoB582SAzsqSQaJZMxnIP/tM+nnZweXkSSkgGnd3fbtt2gCbN8vzM89Iu888I/10\n1izgo4/09bRCsnGj9IcRIySGnj3lGQDuvltye/fdwKuvyt3Pq69Kvykq0rdrraBAzklALky2bCl/\nubLj3rx5kjMtf3/6k9xVp6ZKG7Nny7zBg6VAA1IUIiKAqVMlR/Hx+rbGjtX7oLX8fOlLbm7yvrhY\nzh+TSQpez54Sx48/yl3W66/LcpMnA6dOycVHbq7E4OMjxWv/fn27c+bIRUXbtrL88eMy/coVYPVq\nvYDs3y952rrVNr7u3aWIa8d82DDpr8OGAUlJ0VizJhpJSTLvrbcS7HfwOtSJAuLr64vs7Gz4+/sj\nJycHvr/3bm26xnpeWVpnysnRT/jr4eSkH8DBg+VKY/Bg4MMPpfO99pr9Or/8Ip0dkBOinK9n7OTl\nycdoq1bJey8v4MgReZ4xQz/xfvut4vYuXJATC5Bn60EoKUkG2DKf8lVZcbEMslou+vYFHn5YnhMT\npVBqLl0qv41TpyQm68GwpESex461Xfa22/Tl7rxTCmGjRlL4u3SR56++kkHnv/+V5erXB/7+dykY\nM2dKIZs5E9i7V64aNWvWyEnr5SV3E+3ayR3o0aPAgw/qx8BaXp5ciZZnxw65WADkytZkktjPnJE7\njjNnZN7zz5e/PiD5rQqTyfZqOipKnj/5RAZ4jZOTFEp3dyA2VopfbCywbZsURi1nc+bIANqkiRSC\n63XlCnDggH6sMjJk8HV1lf7s6SkF5Hq4uOjtTZsmg33XrnIMf/5ZX067uAEkB6NHy53DmjVyd38t\nP/0kx04rINXBuugGBUkRDQrSCxggg/uOHfr7xo3lrqFTJ2DRIvlIOD5eCs1jj5VfJCvj6Sk5u/de\nuWD78kspVleuAIMG6YWmJtWJn/HGxMTAYrEAACwWC2JiYkqnr169GgCwa9cuuLu7233/UdMKCuSK\ng2reiRP2V5WAfBySm1v78fyRJCZWXPCrqlEj/Qr9wAEZyOPj5a74Rn3/PfDrr9de7tIl/e6BdFev\nln+3VNNq/Q5k5MiRSElJQV5eHvz8/PD4448jISEBo0aNgsViQYsWLbDm9w81hw8fjm3btiE0NBRu\nbm5YaX15SUREDlXrBSRJ+/CtjK1lP8j73dKlS2syHCIiMqhOfIRFREQ3HxYQIiIyhAWEiIgMYQEh\nIiJDWECIiMgQFhAiIjKEBYSIiAxhASEiIkNYQIiIyBAWECIiMoQFhIiIDGEBISIiQ1hAiIjIEBYQ\nIiIyhAWEiIgMYQEhIiJDWECIiMgQFhAiIjKEBYSIiAxhASEiIkNYQIiIyBAWECIiMoQFhIiIDGEB\nISIiQ1hAiIjIEBYQIiIyhAWEiIgMYQEhIiJD6lQBqV+/PsLDwxEeHo7hw4cDADIyMnDXXXchLCwM\no0ePRlFRkYOjJCIioI4VEB8fH6SmpiI1NRVr164FAEyfPh0zZszAkSNHYDabsXTpUgdHSUREQB0r\nIGVdvXoV27dvx4gRIwAAY8aMwYYNGxwcFRERAYCzowOwlpubi86dO8NkMuHpp5/GXXfdBQ8PDzg5\nOQGQO5Ts7Oxy112xIh6bNgG7dwONG0cDiK61uImIbgZHjyYjLS0Z8fHV016dKiBZWVlo1qwZjh07\nhp49e+K9996r8rqTJ8eja1dgyRIgM7PGQiQiumm1bx+No0ejSwtIQkLCDbVXpz7CatasGQAgKCgI\nPXr0QFZWFn799VcUFxcDAHJycuDr6+vIEImI6Hd1poBcuHABhYWFAOSjrN27d6NDhw7o1asXkpKS\nAAAWiwUxMTGODJOIiH5XZwpIRkYGunbtio4dO6JHjx544okn0KlTJ7zyyit4+eWXERYWhl9++QXT\npk1zdKhERIQ69B1Ihw4dkJqaajc9MDAQu3btckBERERUmTpzB0JERDcXFhAiIjKEBYSIiAxhASEi\nIkNYQIiIyBAWECIiMoQFhIiIDGEBISIiQ1hAiIjIEBYQIiIyhAWEiIgMYQEhIiJDWECIiMgQFhAi\nIjKEBYSIiAxhASEiIkNYQIiIyBAWECIiMoQFhIiIDGEBISIiQ1hAiIjIEBYQIiIyhAWEiIgMYQEh\nIiJDWECIiMgQFhAiIjKEBYSIiAxhASEiIkNuigKyefNmhIWFISQkBIsWLXJ0OEREBMDZ0QFcy5Ur\nVzBlyhSkpKTA29sbkZGRGDBgAMLDw6+rndhYIDUVePppICQE2LULuHQJ+OUXwMWlhoKvo2JjgbQ0\nYNYsoG1bICUFKCoCEhIqX2fHDuDiRclZ/fq1F+/NJDYW2L8fmDMHiIqq+jopKcC5c5Jbs7lm49u0\nCWjcGDh0CDh8GPjhh5rfbtkYjh4Fzp+3367Wzy5ckHlubhW3caM5i40FPvlE2hk1qvLltJw98kjl\ny23ZIuPJ9OnXH8/NqMIC8uabb2L06NFo2LBhbcZjZ8+ePQgKCoK/vz8AYNiwYdiwYYNdAZk5E9i4\nUX+vddL8fOlg6enyvHcvkJEBnDkjy0VEAL6+wOXLlXfE2FhpA5DlnnwS+OIL4NtvgYED7Ze9fBkY\nNgzw9gY+/RTYswcYMcJYDmJjgVOngEmTgIAAYNs2wNUVOHZM9jEzs/wT8bPPAJMJ+Okn4LvvgNxc\nPRe//SZF9Icf9FzExgJLl5YfQ3o68P33tjm7ckXPRUoKkJ0N9OpV+X58/rkUry1bZKDVYqrqABAb\nK8d07Fhg7Vrb6fv3A8ePl5+LTz+VnP3tb/btnTgB/PWvMphYTz90qPz4YmPlGBQVAc89Z5+n/Hxg\n3z7JU/36MhhOmCDbKSzUc7ZzJ3DyJNCggeREa9tsBpQCBg0C/vvfynNRUAAMHw6sWVO1/KWny3EC\ngLNnbY+92Qzk5MggGRgIJCdL/NOmVa3t2FgpSD//rO/j0aPA448D7dpJHyksBPz8gK++0texjr28\nfqblzPoYpKfb56ykBLjvPqB5c2DrVonlvvvsY7x0Sc7FS5fkvDp1SrbVpAnw44967GfPAuPGyfa1\nnEVEAE2bll9009OBrCx9O9YXZLGxss7f/mbbz8rL4ZEjwN//bpuzCxdk/d9+0+Pbtw945hmga9fK\n2/vySxkj1q+Xdc6erb4LBpNSSpU3Y8aMGVi/fj169+6NiRMnonv37je+NQMsFgs2bdqEVatWAQAS\nExOxb98+LFu2rHQZk8kERFutFPD7g4iIdJkA/hOH0FBg5EggISEBFZSAqlGVKCwsVGvXrlWDBw9W\noaGhav78+SonJ6eyVaqdxWJR48aNK32/YsUKNXXqVJtlAKiQEKXuuUcpuXZTqkULeXZ2Vio/Xx7N\nmin10UfyOjJSqXHjlBo4UJarV0+m9+qltxEYqFRIiFLe3kr166dPz8/X12vVyna7I0fq8yIi9PUC\nAuyXs95Ws2blv7Zur0MHvT0fH/21m5t97Nr+W8dev76eC09PpbZskdedOik1dart+mXjGzJEcvbn\nP+vxmEy2uWjTpvJ91GJq1UqPycND2pgyRY5B//62bQQGKtWxo8Rrva1Onco/3g0bVpwLH5+Kj1X7\n9uW3p+WsvPa8vOzby8+X6evW6W3ffru+v+XlLDNTqW7dlBo6VPqzNi8y0rb9ivpFRETly02ZopSL\ni1J9+si2goKUmjdP4ujbV3Jpvd2y/azsPk6ZolSjRjK9spxp7QUH6+15eUkMQ4Yo1bSpbHfKFHkd\nFSXzoqKUGjvWvp+V7Y/duslz2Zxp2woMrPzczMyUvvjii/r0xo1tY+/YUZZr21apuXP16eX1syFD\nZLx46in7c6miflbRMS2bM+21k5NtfG3bVq291q31Npo2lTa0sfNGXHPt77//Xs2ePVu1atVKjRs3\nToWFhak5c+bc0Eavx44dO1Tv3r1L38fHx6t58+bZLANAbd2qJ0vrpIMHy8mr6d1bqc8/l9fLlyv1\nyCOSyGHDlGrQQKZXNBgOGSIHRMt3fr5Sfn7SjraONhjm58uBzs6W14GBSr38sr6cu7v9YKgdXD8/\n/XWjRnp7DRsqlZoqr4ODlZozR14PGCDvrWNv3172PyxMqZkzZbmYGKX8/fVcdO2q1K5d8nrJEqWm\nT9fX14qu9r5dO3m/YoVSkybpOXNz03Nxxx1Kvfuuvo7WSbX32uDg56fUG2/IvHvuUeree6WN8gbo\n8opufr5Srq5Kff+93rafn7Tdt69S4eG2uejQQeaFhCj1v/+rT9eKbn6+Uk2aKPXll/o8b29ZZ9Ag\nadu6vZAQmdepk1KxsfY5U0qOn1acIyNlMNRy5upqm7N33pH3q1cr9eCD+jxAqbNnbS9Uyhbd/HzZ\n7o8/6su1bKkvp/WzsoPZzJlK/fOfsq29eyVGbbuNGim1f7+8DglR6umn7XNW3rFq0ULPma+v3p6n\np1KbN+sXKlOmyLwzZ5S67Tb7Yz9ypFJvvqnUhAmyzvDhUpDKnpv5+Uq9954sb52zvDx57e+v1LJl\n5Z+bzs5KZWXJerGxslx+vlLR0dLftPbc3JRKS5P3s2Yp9cILMr1fPyksZftZfr5SzzyjVEJC+f3M\nw0OpHTv0eS1a6MfKxUVf7rbblNq0SV6Hhys1ebK8HjJEKbNZj695c6Xef19v77bb7C9UtH729tvy\n+q67ZFzU3GgBqfBXWImJiejevTsmT56M4OBgHD58GKtWrcKBAwewpqofuFaDqKgopKenIysrC4WF\nhVi3bh0Glv3SAYC7O2CxAJ06yXcPLVsCb7wB1LvG78zMZuDtt+W7AkDa8PaWzy+bNJFpnp7Af/5j\n+1mt2Qzccw/QqJGsc889wF13yXSzWT5v11736QM0bCjL9eolMZrN8r5+fWk3KQkICgIefVRe9+0r\nn4FqbTRvLvGYzbJ/DRrI63/+E3By0mNv0gR4/XXZ/zFj5EtIs1m+27hWLiwWYPBgwMNDj8/LC1i0\nyPbzUrMZWLVKz5nZDERG6rno1g2IjtbbAOQ7iJYtge7dZTmzGZg9W/8Bg/ZVW0QEsHu3fAY/Y4Z+\nDNzd5XiazRKfh4e0HRQE/OUv0vb8+YCzs74vjRoBb70l87ScWSzAvffKNC23LVvKdiwWICxMvlBt\n2dI2ZxaL9IMlS2TexIlyjLWcacembN969FG9L7zzjm3OoqIkxvL6pMmk58/fX3KVlCT569pVj107\nvhaLnrOkJOlnHTvKPC23nTpJDitiNgMtWuj7Mny4bc78/W3ba99ejlWHDnrO/vUvPWdms/Rhd3d5\nPWmS5KIsrb22bW3jK5sziwW44w4gLq78XNerp+elRw+9P3bvDnTpos/Tzp2y68fF6f3HbJbvOjw8\n7JdbsMC2nzVuDKxcadumljM/P327AQH6ONWhA/Dgg3KsBg2S80xbrrycmc1AYqJtbjt2lG1bLNKX\nBgyQ9oYPl/i09rp21c+5mTOr9wcwFQ4p6enpSExMxI4dOzBhwoTSL9OdnJxg0UaFWuDm5obly5dj\n0KBB6NixI0aPHo2IiIhylzWbgYcfLr+TVpU2GGoHpmtXfTCsbJ0nn7z2r7nMZvnSS1vOejA0m4GY\nGL1zP/+83kmvJ3Z/f33QvV5mM7B8uW0nDQuTXFxPG//zP3on1fJ2rS/sLBYphB99JANR3762RVcb\nDMtua9AgyVl5cWhFt+x066Jbdt7YseX/8sdslgHO3d1++ooV1y7ORpnNFRfdssv16yc5KzsYWixy\nTqxebeyLU7MZePFF2wsVDw9g2TI5VhXlrKq0C5UFC659nnXpUn7RrWydp56quV9aWhfdstNfeqni\nfvbnP+vF/7XXbqz/mM3yI5D69eX1u+/qRbemVThEVfb3FpGRkTUSTEUGDhxY7l1HTdMGw/Xra33T\ntxzrq+my0599Fpg3zzFx/RFoV6Jlr6ZvpL2WLe2L6Y20p11N083lpvhDQiIiqntYQIiIyBAWECIi\nMoQFhIiIDGEBISIiQ1hAiIjIEBYQIiIyhAWEiIgMYQEhIiJDWECIiMgQFhAiIjKEBYSIiAxhASEi\nIkNYQIiIyBAWECIiMoQFhIiIDGEBISIiQ1hAiIjIEBYQIiIyhAWEiIgMYQEhIiJDWECIiMgQFhAi\nIjKEBYSIiAxhASEiIkNYQIiIyBAWECIiMoQFhIiIDKkTBWTChAkICgpCeHg4wsPDceTIEQCAUgrT\np09HaGgoIiIikJqa6uBIiYhI4+zoAADAZDIhMTERPXv2tJm+bt06HDt2DGlpadi9ezcmTJiAgwcP\nOihKIiKyVifuQAC52yhr48aNGDt2LACgW7duOH/+PHJycmo7NCIiKkeduAMBgMceewwlJSXo378/\nFi9eDFdXV2RnZ8PHx6d0GV9fX7tpmhUr4rFpE7B7N9C4cTSA6FqLnYjoZnD0aDLS0pIRH1897dVa\nAenfvz9Onz5tN/3555/HokWL4OXlhStXruDhhx/GvHnzMHfu3Otqf/LkeHTtCixZAmRmVk/MRER/\nJO3bR+Po0ejSApKQkHBD7dVaAdm6des1l3F1dcX48ePxyiuvANDvODQ5OTnw9fWtsRiJiKjq6sR3\nIGfPngUAlJSUYP369QgNDQUAxMTEYPXq1QCAXbt2wd3dvdyPr4iIqPbVie9AJk+ejKysLJw/fx7h\n4eFYvHgxAGD48OHYtm0bQkND4ebmhpUrVzo4UiIi0tSJAvLBBx9UOG/p0qW1GAkREVVVnfgIi4iI\nbj4sIEREZAgLCBERGcICQkREhrCAEBGRISwgRERkCAsIEREZwgJCRESGsIAQEZEhLCBERGQICwgR\nERnCAkJERIawgBARkSEsIEREZAgLCBERGcICQkREhrCAEBGRISwgRERkCAsIEREZwgJCRESGsIAQ\nEZEhLCBERGQICwgRERnCAkJERIawgBARkSEsIEREZAgLCBERGcICQkREhrCAEBGRIbVaQJKSktC+\nfXs4OTlh+/btNvMWLFiAkJAQhIWFYcuWLaXTN2/ejLCwMISEhGDRokW1GS4REVXCuTY31qFDB3zw\nwQeYOnUqTCZT6fT9+/dj9erVOHz4MHJyctC9e3ecOHECJSUlmDJlClJSUuDt7Y3IyEgMGDAA4eHh\ntRn2H9L588CKFcDmzcCePUB+PhAfD3TsaLtcfj7w9tvAl18CqanAyZOyXLdujoiabmbHjgF790r/\nuXRJHvHxQHR01dvIywPeew84eBA4ehRIS7v+Nqj61GoBCQ4OLnf6hg0bMGLECDg7O6Nly5Zo06YN\n9uzZg5KSEgQFBcHf3x8AMGzYMGzYsKHSAvLDD8DXX0ununhRBsroaCAgADh7Fpg9G2jUCCgo0Of5\n+gJFRdXbEYuLgfnzARcX4JtvgJ9/BrKzgSZNqt7G+fMSEyAD+csvA02bAllZQEhI+etcvAj8+9+A\nl5cM+oWF0ka7dsDly3p7V68CTk7y+tFHgagowNtb4rTWtCnw0ENA377A4cPAqVPAgAFy8l+PhAR5\n3r0bOHAAeOMN2X529o3lPTsb+PZb4NdfgZ9+kvi19goKgFdekX3YtQtwdpZ5AQGSJy0XubmSs+bN\ngZwc4Pbbrz8Orf8AEsuaNUBKin2eSkr05VJTgePHgSVLgAYN9GN1rVycPau3UVgILFgA1K8v+1+R\nwkLghReAxo3lgmHrVmD9eqBZM6BePf0cOXdOcubpKf0sKOg6EwHbfnbyJJCYCGgfKnTpIvOKioDh\nw4GuXWX63/8OrFoFfPUVkJwseYmOBvz8pK9qebn9dmDUKOCBB+T82LIFOHJE1ikpkX6qbbegQMYD\nNzfbnBUUAIMGSf/LzZW+8d//ynLffafn4pdfgBEjJLe5uYDJBHz8MeDuLrm27mevvio527tX8hkf\nDwQG6kUSAE6fBpYtA1q0kPPIbC4/Z8eOAdOmyVjx88/Sn44ckfPzwgU9vhMngHfekfgvX7Y9Btb9\n7NtvpegeOnT9x/JaarWAVCQnJwddunQpfe/r64vs7GwopeDj42Mzfd++feW2MWxYPFxcpGN6eEQD\niEZMjJxcztfYy6IiIDISmDFD3r/4oj4vLU0OnnYwmjXTD2Bqqgyqr78unTE9XZ/Xrp10OAAYOlQf\nSA4eBDIy9M5XWAj84x9S1LZskfcffign+m+/yTrR0dJ+QYG8/+47fRDw9JTiom334kU99u7dpfMB\nEmvr1noY7CsPAAAMO0lEQVR72v4AcvK9/rq8vnRJtqO1l5Gh34GYzXJSffWVnNRubvpyx49Lmy+9\nBCiln7DaYBgXJ/mIigJatZICqBXViIjKjw+g3yEBst3hw+W4njkjgwwAjB0rj+RkeVifVPffL7ED\nwL59+vGIjgbuuQfYuVPPU16enJTNm0vM1vv42GNyYru6ykASHw906iQXClp80dGyfUD6T0aGPhhq\nxansMdCOg/a4/XagTx+Z/v33Etfy5bLPZ8/qbdx5p/QpQPrfP/8JWCySo8xM/RjUry+DdIsW9tu1\nNmOG7BsgcX/9NfD553Lsz53T27twQQbD5s1lUD1zRi4O7rhDv4vVtpufL+9/+kn6o9aGdaH09ATG\nj5eLE2tKSZ/S8vbyy/q8Jk1kgB8xQt5HRQExMXLMtDwCco6eO6fH9MgjgI+PDPTWcnNlkO3fX94v\nXCjnu3ahpbl8WfqIr6+8b9BABvygINvtlu1nrVrJsQSA/ful/+/ZI9u4ckVf7uGHJa8NG8rYpJTk\n8epV4Nlny794fOopuTCLj5eLV6X09po1k2IDAN98k4xjx5JL+/SNMimlbap69O/fH6dPn7abPn/+\nfNx///0AgN69eyMhIQE9e/YEAEydOhVdunTBI488AgAYP348Bg0aBKUUNm7ciFWrVgEAEhMTsW/f\nPixbtsx2J0wmaLtRVCTVVzsJqsr6wH/wgXSaJk2kk50/D9x3n32nt1ZYKI/GjSvfzuXLwKefyoAA\nyFVqVJR0kC5dpP2GDWUfCgtlILgeR45IR23Y8PrWq0hSklwl/n4TWCWXL8ugo9X6rVv1k7JsDq3z\nvngx0Lmz5D4nR47h7bfLCZqRAYwebb/+zp2S806d7OM4ehRo08a+L5w5I8c4Nrbq+wTIAGM2S3vW\ncV++DKxbJ8Wrsj5i1KZNcjUbHCyDw9Wr5ffvHTtksAgJkb7z6acyUAFy9/GXv0iujMZovc/Llsld\nQPPmUjC6dpXiYbS9VauAu++Wi5zK+siGDZKH8paryIULUrisrkVvmHVM778vha9JE9uYrl6Vi5by\n7mhTU6Xf3HWXPB86JOe/0Ri+/VbOmX79gJ49AQ8POZeuxXrsNKLaC0hVlC0gc+fOhVIKzz77LACg\nT58+mDt3LkpKShAXF4cvvvgCAJCQkABnZ2fMmTPHpr0bTQLRH92+fTLQ169fPe298QYwZIgUkOpQ\nXCx3DmXvCqhm3ejY6bDDZR10TEwM1q5di6KiImRmZiI9PR1dunRBVFQU0tPTkZWVhcLCQqxbtw4D\nBw50VMhENy3tLre6xMZWX/EA5K6TxePmU6vfgSQlJWHWrFnIy8vDyJEjERQUhJSUFHTu3BmjRo3C\nnXfeCScnJyQmJsLFxQUuLi5Yvnw5Bg0ahOLiYjz00EOIqMoH5kREVOMc8hFWdeNHWERE1++m/QiL\niIhubiwgRERkCAsIEREZwgJCRESGsIAQEZEhLCBERGQICwgRERnCAkJERIawgBARkSEsIEREZAgL\nCBERGcICQkREhrCAEBGRISwgRERkCAsIEREZwgJCRESGsIAQEZEhLCBERGQICwgRERnCAkJERIaw\ngBARkSEsIEREZAgLCBERGcICQkREhrCAEBGRISwgRERkCAsIEREZwgLyB5OcnOzoEOoM5kLHXOiY\ni+pTqwUkKSkJ7du3h5OTE7Zv3146PTk5GZ6enggPD0d4eDgWLFhQOm/z5s0ICwtDSEgIFi1aVJvh\n3pR4cuiYCx1zoWMuqo9zbW6sQ4cO+OCDDzB16lSYTKbS6SaTCUOHDsWbb75ps/yVK1cwZcoUpKSk\nwNvbG5GRkRgwYADCw8NrM2wiIipHrd6BBAcHo23btnbTlVJQStlN37NnD4KCguDv7w8XFxcMGzYM\nGzZsqI1QiYjoWpQDREdHq+3bt5e+T05OVi1atFBhYWGqX79+6tChQ0oppd599101bty40uVWrFih\npk6datceAD744IMPPgw8bkS1f4TVv39/nD592m76/Pnzcf/995e7TufOnZGRkQE3Nzd8+OGHuP/+\n+5GRkVHlbapy7l6IiKhmVXsB2bp163Wv07hx49LXQ4cORWxsLHJzc+Hn54ecnJzSednZ2fDz86uW\nOImI6MY47Ge81ncN586dK329c+dOmEwmeHl5ISoqCunp6cjKykJhYSHWrVuHgQMHOiJcIiIqo1Z/\nhZWUlIRZs2YhLy8PI0eORFBQEFJSUrB+/XosWbIEV69ehYuLC9577z3Uq1cPbm5uWL58OQYNGoTi\n4mI89NBDiIiIqM2QiYioIjf0DUodsGnTJtW+fXsVHBysFi5c6OhwatWpU6dUr169VFhYmAoKClJx\ncXFKKaXOnj2r+vXrpzp06KD69++v8vPzHRtoLSkuLlaRkZEqOjpaKXXr5uHcuXPqgQceUB06dFAh\nISHq4MGDt2wunnjiCdWmTRvVrl07NXToUPXrr7+qEydOqG7duqn27durUaNGqcLCQkeHWSMmTpyo\nvLy8VEBAQOm0yvrBtGnTVEhIiAoPD1fffPNNlbZxU/8luvZ3Ihs3bsShQ4dgsViQmprq6LBqjbOz\nM1599VUcPnwYBw4cwOrVq7Fz507ExcWhX79+OHToEPr06YO4uDhHh1orli5diqCgoNK/MbpV8xAb\nG4uYmBgcOnQIhw8fRqtWrW7JXKSmpuL999/H0aNH8d1336F+/fpYuXIlpk+fjhkzZuDIkSMwm81Y\nunSpo0OtERMnTsTmzZttplXUD9auXYtjx44hLS0Nr732GiZMmFC1jVR72atF27dvV7179y59Hx8f\nr+bOnevAiBxr+PDhas2aNapVq1bqxx9/VEoplZmZqVq3bu3gyGpeTk6O6tevn/riiy9K70BuxTzk\n5eUpb29vu+m3Yi5OnTqlgoKC1Llz51RRUZG677771CeffKLc3d3V1atXlVLyJwR9+/Z1cKQ1JyMj\nw+YOpKJ+MHHiRPX222+XLhcQEKCys7Ov2f5NfQeSnZ0NHx+f0ve+vr7Izs52YESOk5mZid27d6Nv\n3742efHx8bklcvL4449j4cKFqFdP79K3Yh6OHTuG5s2bY8yYMWjfvj3Gjx+P33777ZbMRYsWLfDk\nk0/C398fd9xxB5o2bYqIiAh4eHjAyckJwK2TC01F/SAnJ8fQWHpTFxDrfw7lVnbp0iWMHDkSr7zy\nCjw9PR0dTq3bvHkzPDw80Llz51v+b4JKSkpw8OBB/PWvf8XRo0fRqFEjzJ0719FhOcTx48exYMEC\nnDhxAjk5OcjLyzP0ZwZUsZu6gPj6+t7yfydSXFyMUaNGYcyYMRg2bBgA26uHnJwc+Pr6OjLEGvfV\nV19hw4YNCAwMxJgxY7B7924MGTLklssDAPj5+cHT0xM9evQAIH9XdfDgQfj5+d1yudi7dy+ioqLQ\nrFkzuLi4YPDgwdixYwd+/fVXFBcXA7h1cqGp6Jwoe8dR1bzc1AWEfycCTJkyBYGBgZg1a1bptJiY\nGFgsFgCAxWJBTEyMo8KrFc899xxOnjyJjIwMrF69Gt26dcNHH310y+UBkALi5+eHQ4cOAQC2bduG\n4OBgDBw48JbLRZs2bfD111/j4sWLUErhiy++QHBwMHr16oWkpCQAt04uNBWdEzExMVi9ejUAYNeu\nXXB3d7f5SKtC1fqNjQNs3Lix9Ge88+fPd3Q4tWrnzp3KZDKpjh07qk6dOqlOnTqpjz/++Jb9yaZS\nSm3btq30hxW3ah4OHDigOnfurEJCQtTAgQPVuXPnbtlcxMXFqdatW6u2bduqUaNGqYKCApuf8T74\n4IN/2J/xjhgxQnl7eysXFxfl6+urXnrppUr7wWOPPVb6M979+/dXaRsmpW7xD42JiMiQm/ojLCIi\nchwWECIiMoQFhIiIDGEBISIiQ1hAiGrAvn370LFjR1y5cgUXL15EWFgY0tLSHB0WUbXir7CIasgz\nzzyDy5cvo6CgAH5+fpg9e7ajQyKqViwgRDWkqKgIkZGRaNCgAXbt2sV/eof+cPgRFlENycvLw8WL\nF3HhwgUUFBQ4Ohyiasc7EKIaMnjwYIwdOxYnTpzAqVOn8Oqrrzo6JKJqVav/pS3RreLtt9+Gq6sr\nRo8ejZKSEtx9991ITk5GdHS0o0Mjqja8AyEiIkP4HQgRERnCAkJERIawgBARkSEsIEREZAgLCBER\nGcICQkREhvw/ShnkwwJ1mscAAAAASUVORK5CYII=\n" | |
| } | |
| ], | |
| "prompt_number": 10 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "**Compute the covariance**\n", | |
| "\n", | |
| "The diagonal element of the covariance is automatically converted to an effective $1-\\sigma$ confidence interval on the fitted parameter. In this case, we used $\\sigma_y=100$, so the answer of $\\sigma_\\mu = \\sigma_y / N^{1/2} = 100 / 100^{1/2} = 10$ is correct." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "ui.covar()" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "Dataset = 1\n", | |
| "Confidence Method = covariance\n", | |
| "Iterative Fit Method = None\n", | |
| "Fitting Method = levmar\n", | |
| "Statistic = chi2gehrels\n", | |
| "covariance 1-sigma (68.2689%) bounds:\n", | |
| " Param Best-Fit Lower Bound Upper Bound\n", | |
| " ----- -------- ----------- -----------\n", | |
| " c1.c0 -0.11446 -10 10\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 11 | |
| } | |
| ], | |
| "metadata": {} | |
| } | |
| ] | |
| } |
The code Eric gave on the mailing list also gives the expected result:
chi2 = np.sum(((yn-const(x, *popt))/sigma)**2)
print np.sqrt(np.diag(pcov)/(chi2/(x.shape[0]-1)))
FYI: You can see the factor that is multiplied to pcov in scipy.optimize.fit_curve
here.
@cdeil - thanks, that worked!
BTW, I did not get any notification. Looks like you need to give an explicit @taldcroft.
I prefer your version of the scaling correction as being a little more readable. :-)
I have a question regarding the sigma; How do you define the value of sigma that you will use?
You write above,
"chi = (yn - const(x, *popt)) / sigma"
but what is the correct value for sigma and how can you caculate it from your data?
Thank you very much!!
I had this gist pop up as a search result while I was looking for something else concerning curve_fit. However, I wonder, isn't the absolute_sigma flag exactly what you need (and the only thing you need) in the case you're describing above?
import numpy as np
import scipy.optimize
f = lambda x, a: a * x
a = 3
x = np.linspace(0, 10, 10)
s = np.random.normal(size=len(x))
y = f(x, a) + s
print("{:-^72s}".format(" sigma=s "))
popt, pcov = scipy.optimize.curve_fit(f, x, y, sigma=s)
perr = np.sqrt(np.diag(pcov))
print("popt =", popt)
print("pcov =", pcov)
print("perr =", perr)
print("{:-^72s}".format(" sigma=100*s "))
popt, pcov = scipy.optimize.curve_fit(f, x, y, sigma=100*s)
perr = np.sqrt(np.diag(pcov))
print("popt =", popt)
print("pcov =", pcov)
print("perr =", perr)
print("{:-^72s}".format(" sigma=s, absolute_sigma=True "))
popt, pcov = scipy.optimize.curve_fit(f, x, y, sigma=s, absolute_sigma=True)
perr = np.sqrt(np.diag(pcov))
print("popt =", popt)
print("pcov =", pcov)
print("perr =", perr)
print("{:-^72s}".format(" sigma=100*s, absolute_sigma=True "))
popt, pcov = scipy.optimize.curve_fit(f, x, y, sigma=100*s, absolute_sigma=True)
perr = np.sqrt(np.diag(pcov))
print("popt =", popt)
print("pcov =", pcov)
print("perr =", perr)results for example in
------------------------------- sigma=s --------------------------------
popt = [ 2.98538192]
pcov = [[ 0.00015454]]
perr = [ 0.01243135]
----------------------------- sigma=100*s ------------------------------
popt = [ 2.98538192]
pcov = [[ 0.00015454]]
perr = [ 0.01243135]
--------------------- sigma=s, absolute_sigma=True ---------------------
popt = [ 2.98538192]
pcov = [[ 0.00016045]]
perr = [ 0.01266702]
------------------- sigma=100*s, absolute_sigma=True -------------------
popt = [ 2.98538192]
pcov = [[ 1.60453399]]
perr = [ 1.26670201]
As you can see, perr does get bigger, if I use sigma=100*s - but you have to set the absolute_sigma flag to True.
In the third call you can see that perr is (more or less) the same as in the first two calls to curve_fit. This is because the sigma argument's values are supposed to be weights in standard deviations of the y data and we're using np.random.normal to generate these. Therefore, absolute_sigma=True has no effect in this case because absolute and relative sigma weights are the same value. (The small remaining difference 0.01243135 vs 0.01266702 is due to the low number of points, which makes the standard deviation of the points deviate a little from the "real" standard deviation used to generate them.)
Conclusion: if one is using the sigma argument of curve_fit for absolute errors instead of weights (standard deviations) one should use absolute_sigma=True for perr to be appropriately scaled to these errors. More information in the docs for curve_fit as well https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
If you write
instead of
you will get the expected result of 10.
(the notation in my and your example was different.)