Investigating `scipy.optimize.curve_fit` covariance output

curve_fit.ipynb

@taldcroft

I had this gist pop up as a search result while I was looking for something else concerning curve_fit. However, I wonder, isn't the absolute_sigma flag exactly what you need (and the only thing you need) in the case you're describing above?

import numpy as np
import scipy.optimize

f = lambda x, a: a * x

a = 3
x = np.linspace(0, 10, 10)
s = np.random.normal(size=len(x))
y = f(x, a) + s

print("{:-^72s}".format(" sigma=s "))
popt, pcov = scipy.optimize.curve_fit(f, x, y, sigma=s)
perr = np.sqrt(np.diag(pcov))
print("popt =", popt)
print("pcov =", pcov)
print("perr =", perr)

print("{:-^72s}".format(" sigma=100*s "))
popt, pcov = scipy.optimize.curve_fit(f, x, y, sigma=100*s)
perr = np.sqrt(np.diag(pcov))
print("popt =", popt)
print("pcov =", pcov)
print("perr =", perr)

print("{:-^72s}".format(" sigma=s, absolute_sigma=True "))
popt, pcov = scipy.optimize.curve_fit(f, x, y, sigma=s, absolute_sigma=True)
perr = np.sqrt(np.diag(pcov))
print("popt =", popt)
print("pcov =", pcov)
print("perr =", perr)

print("{:-^72s}".format(" sigma=100*s, absolute_sigma=True "))
popt, pcov = scipy.optimize.curve_fit(f, x, y, sigma=100*s, absolute_sigma=True)
perr = np.sqrt(np.diag(pcov))
print("popt =", popt)
print("pcov =", pcov)
print("perr =", perr)

results for example in

------------------------------- sigma=s --------------------------------
popt = [ 2.98538192]
pcov = [[ 0.00015454]]
perr = [ 0.01243135]
----------------------------- sigma=100*s ------------------------------
popt = [ 2.98538192]
pcov = [[ 0.00015454]]
perr = [ 0.01243135]
--------------------- sigma=s, absolute_sigma=True ---------------------
popt = [ 2.98538192]
pcov = [[ 0.00016045]]
perr = [ 0.01266702]
------------------- sigma=100*s, absolute_sigma=True -------------------
popt = [ 2.98538192]
pcov = [[ 1.60453399]]
perr = [ 1.26670201]

As you can see, perr does get bigger, if I use sigma=100*s - but you have to set the absolute_sigma flag to True.

In the third call you can see that perr is (more or less) the same as in the first two calls to curve_fit. This is because the sigma argument's values are supposed to be weights in standard deviations of the y data and we're using np.random.normal to generate these. Therefore, absolute_sigma=True has no effect in this case because absolute and relative sigma weights are the same value. (The small remaining difference 0.01243135 vs 0.01266702 is due to the low number of points, which makes the standard deviation of the points deviate a little from the "real" standard deviation used to generate them.)

Conclusion: if one is using the sigma argument of curve_fit for absolute errors instead of weights (standard deviations) one should use absolute_sigma=True for perr to be appropriately scaled to these errors. More information in the docs for curve_fit as well https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html