Python Forum

I am trying to plot a simple curve in Python using matplotlib with a Gaussian fit which has both x and y errors. The best fit curve should take into account both errors. I have also built in a way of ignoring the baseline and to isolate the data to only a certain x range. My main issue is that I cant manage to get the Scipy ODR to work. It is reporting the wrong number of parameters (see in code comments). I know the lingress line is incorrect as this was taken from a straight line fir ODR but not sure what to put there. Any help would be appreciated. I will paste the full Python code below and see if I can upload the data file. I may have other issues once the parameters issue is solved! Thanks.

Here is the data file if it works - CSV DATA FILE

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from scipy import odr
from scipy.stats import linregress
# These two lines enable formatted printing of Pandas DataFrames
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

def GaussModelConstant(x, a, x0, sigma, const):
    y_gauss = a*np.exp(-(x-x0)**2/(2*sigma**2)) + const
    return y_gauss

# change the filename below as required
results = pd.read_csv('result-calibration.csv', skip_blank_lines=True, skipinitialspace = True)
results.head()

x = results['Channel'][260:280]
y = results['Counts'][260:280]
results['Counts_Error'] = np.sqrt((results['Counts']))
xerror = 0.5
yerror = results['Counts_Error'][260:280]

model = odr.Model(GaussModelConstant)
rdata = odr.RealData(x, y, sx = xerror, sy = yerror)
# This line doesnt seem right I have tried replacing with the one below which I need to add somewhere?
#p0 = [1400, 270, 10, 200]
init_guess = linregress(x, y)[0:2]
# Get error when I run - TypeError: GaussModelConstant() missing 3 required positional arguments: 'x0', 'sigma', and 'const'
odr = odr.ODR(rdata, model, beta0 = init_guess)
result_outputs = odr.run()
popt = result_outputs.beta
perr = np.sqrt(np.diag(result_outputs.cov_beta))

a_fit     = popt[0]
x0_fit    = popt[1]
sigma_fit = popt[2]
const_fit = popt[3]
perr = np.sqrt(np.diag(pcov))
a_err     = perr[0]
x0_err    = perr[1]
sigma_err = perr[2]
const_err = perr[3]

print('fit parameters with error estimates')
print('***************************************************')
print(f'A     = {a_fit: .3g} +/- {a_err: .3g}')
print(f'x0    = {x0_fit: .3g} +/- {x0_err: .3g}')
print(f'sigma = {sigma_fit: .3g} +/- {sigma_err: .3g}')
print(f'Constant \'base\' offset = {const_fit: .3g} +/- {const_err: .3g}')
print('***************************************************')

plt.figure(figsize=(15, 15))
plt.errorbar(x, y, yerr = yerror, xerr = xerror, fmt = 'o', ecolor = 'k', label = 'Errors')
plt.plot(x, GaussModelConstant(x, *popt), color = "red")

Looks like your GaussModelConstant function may not be properly defined. According to the scipy odr documentation (https://docs.scipy.org/doc/scipy/reference/odr.html) the fitting function has the following form:

def f(B, x):
    '''Linear function y = m*x + b'''
    # B is a vector of the parameters.
    # x is an array of the current x values.
    # x is in the same format as the x passed to Data or RealData.
    #
    # Return an array in the same format as y passed to Data or RealData.
    return B[0]*x + B[1]

If you correct that and use init_guess = [1400, 270, 10, 200] it works fine.

83dons

Tuxedo