Fitting signal with background with scipy curve_fit()

user_9314030 0
Viewed 1065

Below is a toy model of my current problem. I have a background with a shape of wide gaussian and a sharp signal peak that is slighly off-centered from the background mean. I want to extract the properties of the sharp peak (i.e. width, peak loc, and etc). Here is a link to the plot with failed fits to the peak:

graph of background + offset signal.

For initial guess of the parameters, weirdly, when a value much larger than the actual standard deviation of the distribution is used, the fit performs better...something is wrong, but can't figure out. I would appreciate for any tips on fitting peaks with backgrounds.

Below is what I tried.

#Fake Data
data  = np.random.normal(loc=3.25, scale=0.01, size=15000)
data2 = np.random.normal(loc=3.0, scale=0.3, size=25000)

#Bins
bins  = np.arange(0, 6.1, 0.1)

#Hitogram with its defined bins
data_entries_1, bins = np.histogram(data, bins=bins)
data_entries_2, bins = np.histogram(data2, bins=bins)

#Add two generated histograms - Final y data
data_entries = data_entries_1 + data_entries_2

#Cetner of each bins - Final x data
bin_centers = np.array([0.5*(bins[i] + bins[i+1]) for i in range(len(bins)-1)])

#fit func 1
def fit_func1(x, A, mu, sigma):
    #Define functions here
    first_func  = A*np.exp(-1.0*(x - mu)**2 / (2*sigma**2))
    return first_func

#fit func 2
def fit_func2(x, B, mu2, sigma2):
    #Define functions here
    second_func = B*np.exp(-1.0*(x - mu2)**2 / (2*sigma2**2))
    return second_func

#total fit function
def fit_func(x, A, mu, sigma, B, mu2, sigma2):
    #Define functions here
    first_func  = A*np.exp(-1.0*(x - mu)**2 / (2*sigma**2))
    second_func = B*np.exp(-1.0*(x - mu2)**2 / (2*sigma2**2))
    final_func  = first_func + second_func
    return final_func

#Fit it
popt1, pconv1 = curve_fit(fit_func1, xdata=bin_centers, ydata=data_entries_1, p0=[20000, 3.25, 1.])
popt2, pconv2 = curve_fit(fit_func2, xdata=bin_centers, ydata=data_entries_2, p0=[2000, 3.0, 0.3])
popt, pconv   = curve_fit(fit_func, xdata=bin_centers, ydata=data_entries, p0=[20000, 3.25, 1.,\
                                                                           2000, 3.0, 0.3])

# Generate enough x values to make the curves look smooth.
xspace = np.linspace(0, 6, 100)

# Plot the histogram and the fitted function.
plt.step(bin_centers, data_entries, label=r'Histogram entries')
plt.plot(xspace, fit_func1(xspace, *popt1), label='Fitted function1')
plt.plot(xspace, fit_func2(xspace, *popt2), label='Fitted function2')
plt.plot(xspace, fit_func(xspace, *popt), label='Fitted function', linestyle="--")
plt.xlim(1, 5)
plt.legend()
plt.show()

UPDATE Following the suggestions from the selected answer, bin width is reduced to be smaller than the sigma of the actual peak to be fitted. Also, to reduce the # of free parameters to fit, mu of the fitting gaussian function is fixed to a constant (actual mean of the distribution).

#Generate Fake Data
data  = np.random.normal(loc=3.25, scale=0.04, size=15000)
data2 = np.random.normal(loc=3.0, scale=0.3, size=25000)

#Bins
bins  = np.arange(0, 6.1, 0.02)

#Compute mean to use as contraints when fitting
mids = np.array([0.5*(bins[i] + bins[i+1]) for i in range(len(bins)-1)])
mean_1 = np.average(mids, weights=data_entries_1)
mean_2 = np.average(mids, weights=data_entries_2)

#fit func 1
def fit_func1(x, A, sigma):
    #Define functions here
    first_func  = A*np.exp(-1.0*(x - mean_1)**2 / (2*sigma**2))
    return first_func

#fit func 2
def fit_func2(x, B, sigma2):
    #Define functions here
    second_func = B*np.exp(-1.0*(x - mean_2)**2 / (2*sigma2**2))
    return second_func

#total fit function
def fit_func(x, A, sigma, B, sigma2):
    #Define functions here
    first_func  = A*np.exp(-1.0*(x - mean_1)**2 / (2*sigma**2))
    second_func = B*np.exp(-1.0*(x - mean_2)**2 / (2*sigma2**2))
    final_func  = first_func + second_func
    return final_func

enter image description here

curve-fitting scipy python
2 Answers

There are several problems.

  1. the plt.step function does not do what you think it does. It takes the edges of the bins, not the centers of the bin.
  2. The rarrow peak is a normal distribution with a sigma much smaller than the bin width. Essentially, you try to fit a three-parameter function (A, mu, sigma) on a single (x, y) value; this is bound to fail. The outcome of this fit is not reproducible on my system; rerunning the same code will sometimes even generate an error depending on the random generator output.

See below for an illustration of these points.

Close-up of multi-gauss fit

def normal(x, sigma):
    a = 0.5/sigma**2
    return np.sqrt(a/np.pi) * np.exp(-a*x**2)
    
def truefunc(x):
    return 1500*normal(x-3.25, 0.01) + 2500*normal(x-3.0, 0.3)

plt.close('all')
xspace = np.linspace(3, 3.5, 200)
plt.plot(bin_centers, data_entries, 'ko', label=r'Histogram entries')
plt.plot(xspace, truefunc(xspace), label='True distribution')
plt.plot(xspace, fit_func1(xspace, *popt1), label='Fitted function1')
plt.plot(xspace, fit_func2(xspace, *popt2), label='Fitted function2')
plt.plot(xspace, fit_func(xspace, *popt), label='Fitted function', linestyle="--")
plt.xlim(3, 3.5)
plt.legend()
plt.show()

So, you will need to decrease the bin size so that you can actually resolve the width of the narrow peak, or redefine fit_func1 to take two parameters: peak height and peak position (mu). Fix sigma to bin_width/6. You'll have to postprocess the fit result to get the area under the peak consistent with the histogram.

If you decrease the bin size, the histogram will be noisy. It's probably a good idea to use the sigma parameter of curve_fit; set it to

np.sqrt(data_entries + 0.25)

This is statistically not completely correct, but it's much better than assuming a fixed error on all histogram entries. (You can play with the 0.25 value; it should be > 0 and < 1).

user_62288910

If your goal is to analyze the peak properties you can go without fitting. To get the peak position just do:

peak_pos = bin_centers[data_entries.argmax()]

If you have multiple such peaks you could also use scipy.signal.find_peaks.

To get the peak width you can use scipy.signal.peak_widths. Optionally you can fit the background broad Gaussian and subtract before analyzing the peak.

user_37672390
Related