Fitting Gaussian curve to data in python

Question

I'm trying to fit and plot a Gaussian curve to some given data. This is what I have so far:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Generate data
mu, sigma = 0, 0.1
y, xe  = np.histogram(np.random.normal(mu, sigma, 1000))
x = .5 * (xe[:-1] + xe[1:])

def gauss (x, y):
    p = [x0, y0, sigma]
    return p[0] * np.exp(-(x-p[1])**2 / (2 * p[2]**2))

p0 = [1., 1., 1.]

fit = curve_fit(gauss, x, y, p0=p0)
plt.plot(gauss(x, y))
plt.show()

When I run the code I get this error:

TypeError: gauss() takes exactly 2 arguments (4 given)

I don't understand where I have given my function 4 arguments. I'm also not convinced I'm using the curve function correctly, but I'm not sure exactly what I'm doing wrong. Any help would be appreciated.

Edit

Here's the Traceback:

Traceback (most recent call last):
  File "F:\Numerical methods\rw893  final assignment.py", line 21, in <module>
    fitE, fitI = curve_fit(gauss, x, y, p0=p0)
  File "F:\Portable Python 2.7.5.1\App\lib\site-packages\scipy\optimize\minpack.py", line 515, in curve_fit
    res = leastsq(func, p0, args=args, full_output=1, **kw)
  File "F:\Portable Python 2.7.5.1\App\lib\site-packages\scipy\optimize\minpack.py", line 354, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
  File "F:\Portable Python 2.7.5.1\App\lib\site-packages\scipy\optimize\minpack.py", line 17, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
  File "F:\Portable Python 2.7.5.1\App\lib\site-packages\scipy\optimize\minpack.py", line 427, in _general_function
    return function(xdata, *params) - ydata
TypeError: gauss() takes exactly 2 arguments (4 given)

Answer 1

Probably your callback is called in curve_fit with a different number of parameters. Have a look at the documentation where it says:

The model function, f(x, ...). It must take the independent variable as the first argument and the parameters to fit as separate remaining arguments.

To make sure this works out you might want to take *args after the first argument and have a look at what you get.

Answer 2

Check the first scipy documentation docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.optimize.curve_fit.html:

scipy.optimize.curve_fit

scipy.optimize.curve_fit(f, xdata, ydata, p0=None, sigma=None, **kw)

Use non-linear least squares to fit a function, f, to data.

Assumes ydata = f(xdata, *params) + eps

Explaining the idea

The function to be fitted should take only scalars (not: *p0). I want to remind you that you hand over the initialization parameters x0, y0, sigma to the function gauss during the call of curve_fit. You call the initialization p0 = [x0, y0, sigma].

The function gauss returns the value y = y0 * np.exp(-((x - x0) / sigma)**2). Therefore the input values need to be x, x0, y0, sigma. The first parameter x is the data you know together with the result of the function y. The later three parameters will be fitted - you hand over them as initialization parameters.

Working example

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Create data:
x0, sigma = 0, 0.1
y, xe  = np.histogram(np.random.normal(x0, sigma, 1000))
x = .5 * (xe[:-1] + xe[1:])

# Function to be fitted
def gauss(x, x0, y0, sigma):
    p = [x0, y0, sigma]
    return p[1]* np.exp(-((x-p[0])/p[2])**2)

# Initialization parameters
p0 = [1., 1., 1.]
# Fit the data with the function
fit, tmp = curve_fit(gauss, x, y, p0=p0)

# Plot the results
plt.title('Fit parameters:\n x0=%.2e y0=%.2e sigma=%.2e' % (fit[0], fit[1], fit[2]))
# Data
plt.plot(x, y, 'r--')
# Fitted function
x_fine = np.linspace(xe[0], xe[-1], 100)
plt.plot(x_fine, gauss(x_fine, fit[0], fit[1], fit[2]), 'b-')
plt.savefig('Gaussian_fit.png')
plt.show()

Answer 3

from numpy import loadtxt
import numpy as np
from scipy import *
from matplotlib import *
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a, b, c, d, x0):
    return a*np.exp(-(x-x0)**2/(2*d**2)) + c

x = np.linspace(0,4,50)
y = func(x, 2.5, 1.3, 0.5, 1.0, 2.0)
yn = y + 0.2*np.random.normal(size=len(x))

p = [1,1,1,1,1]

popt, pcov = curve_fit(func, x, yn, p0=p)


plt.plot(x,func(x,popt[0],popt[1],popt[2],popt[3],popt[4]))
plt.plot(x,yn,'r+')
plt.show()

This should help. This can also be extended to a 3d Gaussian, then the input array 'x' should be a k-dimensional array for the (x,y) values and 'yn' should be the z-values.