Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Short script that loads a three-column data file and fits the data to a Gaussian model. Times the execution & plots the result.
import numpy as np
from scipy.optimize import leastsq
import time
# Load data from file
rawdata = np.genfromtxt('./gaussianData.xye')
x = rawdata[:, 0]
y = rawdata[:, 1]
sy = rawdata[:, 2]
# Define model & distance function
def gaussian(p, x):
return p[0]*np.exp(-0.5*((x-p[1])/p[2])**2) + p[3]
def residuum(p, x, y, sy):
return (y - gaussian(p, x))/sy
# Initialize start parameter
p = np.array([20, 3, 0.2, 0])
# Also measure the execution time: starting now
t0 = time.time()
# Perform LM fit
p_fit, cov_x, info_dict, mesg, ier = leastsq(residuum, p, args=(x,y,sy), full_output=True)
# Finish time measurement
t1 = time.time()
# Extract standard errors from covariance matrix
std_errs = np.sqrt(np.diag(cov_x))
# Print execution time & fit results
print(f"Execution time: {(t1 - t0)*1e3} ms")
for i in range(len(p_fit)):
print(f"{p_fit[i]} +/- {std_errs[i]}")
# Print additional information of fit performance
print(f"#Function eval: {info_dict['nfev']}")
print(f"#Chi2: {np.sum(residuum(p_fit, x, y, sy)**2)}")
print(f"#Red. Chi2: {np.sum(residuum(p_fit, x, y, sy)**2) / (len(x) - len(p))}")
# Plot the result
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(3, 2))
ax = fig.add_axes([0.15, 0.2, 0.8, 0.7])
ax.errorbar(x,y,sy, ls='None', color='black')
ax.plot(x, gaussian(p_fit, x), marker='None', color='#0EA8DF')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
fig.savefig('gaussianFit_scipy.png', dpi=200)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.