-
-
Save taldcroft/5014170 to your computer and use it in GitHub Desktop.
The code Eric gave on the mailing list also gives the expected result:
chi2 = np.sum(((yn-const(x, *popt))/sigma)**2)
print np.sqrt(np.diag(pcov)/(chi2/(x.shape[0]-1)))
FYI: You can see the factor that is multiplied to pcov in scipy.optimize.fit_curve
here.
@cdeil - thanks, that worked!
BTW, I did not get any notification. Looks like you need to give an explicit @taldcroft.
I prefer your version of the scaling correction as being a little more readable. :-)
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!!
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
If you write
instead of
you will get the expected result of 10.
(the notation in my and your example was different.)