Python gaussian fit on simulated gaussian noisy data
Question: I need to interpolate data coming from an instrument using a gaussian fit. To
this end I thought about using the `curve_fit` function from `scipy`. Since
I'd like to test this functionality on fake data before trying it on the
instrument I wrote the following code to generate noisy gaussian data and to
fit it:

    
    
    from scipy.optimize import curve_fit
    import numpy
    import pylab
    
    # Create a gaussian function
    
    def gaussian(x, a, b, c):
        val = a * numpy.exp(-(x - b)**2 / (2*c**2))
        return val
    
    # Generate fake data.
    zMinEntry = 80.0*1E-06
    zMaxEntry = 180.0*1E-06
    zStepEntry = 0.2*1E-06
    
    x = numpy.arange(zMinEntry, 
                     zMaxEntry, 
                     zStepEntry, 
                     dtype = numpy.float64)    
    
    n = len(x)                 
    
    meanY = zMinEntry + (zMaxEntry - zMinEntry)/2
    sigmaY = 10.0E-06
    
    a = 1.0/(sigmaY*numpy.sqrt(2*numpy.pi))
    y = gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))
    
    
    # Fit
    popt, pcov = curve_fit(gaussian, x, y)
    
    # Print results
    print("Scale =  %.3f +/- %.3f" % (popt[0], numpy.sqrt(pcov[0, 0])))
    print("Offset = %.3f +/- %.3f" % (popt[1], numpy.sqrt(pcov[1, 1])))
    print("Sigma =  %.3f +/- %.3f" % (popt[2], numpy.sqrt(pcov[2, 2])))
    
    
    pylab.plot(x, y, 'ro')
    pylab.plot(x, gaussian(x, popt[0], popt[1], popt[2]))
    pylab.grid(True)
    pylab.show()
    

Unfortunately this does not work properly, the output of the code is the
following:

    
    
    Scale =  6174.816 +/- 7114424813.672
    Offset = 429.319 +/- 3919751917.830
    Sigma =  1602.869 +/- 17923909301.176
    

And the plotted result is (blue is the fit function, red dots is the noisy
input data):

[![enter image description
here](http://i.stack.imgur.com/MCqXY.png)](http://i.stack.imgur.com/MCqXY.png)

I also tried to look at
[this](http://stackoverflow.com/questions/19206332/gaussian-fit-for-python)
answer, but couldn't figure out where my problem is. Am I missing something
here? Or am I using the `curve_fit` function in the wrong way? Thanks in
advance!


Answer: I agree with Olaf in so far as it is a question of scale. The optimal
parameters differ by many orders of magnitude. However, scaling the parameters
with which you generated your toy data does not seem to solve the problem for
your actual application.
[`curve_fit`](http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.optimize.curve_fit.html)
uses
[`lestsq`](http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.optimize.leastsq.html#scipy.optimize.leastsq),
which numerically approximates the Jacobian, where numerical problems arise
because of the differences in scale (try to use the `full_output` keyword in
`curve_fit`).

In my experience it is often best to use
[`fmin`](http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.optimize.fmin.html)
which does not rely on approximated derivatives but uses only function values.
You now have to write your own least-squares function that is to be optimized.

Starting values are still important. In your case you can make sufficiently
good guesses by taking the maximum amplitude for `a` and the corresponding
x-values for `b`and `c`.

In code, it looks like this:

    
    
    from scipy.optimize import curve_fit,fmin
    import numpy
    import pylab
    
    # Create a gaussian function
    
    def gaussian(x, a, b, c):
        val = a * numpy.exp(-(x - b)**2 / (2*c**2))
        return val
    
    # Generate fake data.
    zMinEntry = 80.0*1E-06
    zMaxEntry = 180.0*1E-06
    zStepEntry = 0.2*1E-06
    
    x = numpy.arange(zMinEntry, 
                     zMaxEntry, 
                     zStepEntry, 
                     dtype = numpy.float64)    
    
    n = len(x)                 
    
    meanY = zMinEntry + (zMaxEntry - zMinEntry)/2
    sigmaY = 10.0E-06
    
    a = 1.0/(sigmaY*numpy.sqrt(2*numpy.pi))
    y = gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))
    
    print a, meanY, sigmaY
    # estimate starting values from the data
    a = y.max()
    b = x[numpy.argmax(a)]
    c = b
    
    # define a least squares function to optimize
    def minfunc(params):
        return sum((y-gaussian(x,params[0],params[1],params[2]))**2)
    
    # fit
    popt = fmin(minfunc,[a,b,c])
    
    # Print results
    print("Scale =  %.3f" % (popt[0]))
    print("Offset = %.3f" % (popt[1]))
    print("Sigma =  %.3f" % (popt[2]))
    
    
    pylab.plot(x, y, 'ro')
    pylab.plot(x, gaussian(x, popt[0], popt[1], popt[2]),lw = 2)
    pylab.xlim(x.min(),x.max())
    pylab.grid(True)
    pylab.show()
    

[![enter image description
here](http://i.stack.imgur.com/tk92z.png)](http://i.stack.imgur.com/tk92z.png)

