-
-
Save taldcroft/5014170 to your computer and use it in GitHub Desktop.
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
I have a question regarding the sigma; How do you define the value of sigma that you will use?
You write above,
"chi = (yn - const(x, *popt)) / sigma"
but what is the correct value for sigma and how can you caculate it from your data?
Thank you very much!!