3

I currently want to fit data with errors in x and y and im using the scipy.odr package to get my results. Im just wondering about the correct use of errors sx and sy. So here is an example.

Lets assume im measuring a voltage V for different currents I. So I have measured 1V with an error of +- 0.1V and so on. So if I'm assuming no error in current measurement I could use spipy.curve_fit as follows. Absolute_sigma is set True because my absolute error of +- 0.1V.

So I'm getting:

from scipy.optimize import curve_fit
y=[1,2.5,3,4]
x=[1,2,3,4]
yerr=[0.1,0.1,0.1,0.1]

def func(x, a, b):
    return a*x+b

popt, pcov = curve_fit(func, x, y,sigma=yerr,absolute_sigma=True)

print(popt)
print(np.sqrt(np.diag(pcov)))

[0.95 0.25]
[0.04472136 0.12247449]

In a second step I want to use the odr-package with errors in both, current and voltage. According to the documentation it should be used as follows: sx and sy are the errors for my measurement data. So I should assume to get similar results to curve_fit if I'm using a very small error for sx.

from scipy.odr import *
x_err = [0.0000000001]*x.__len__()
y_err = yerr

def linear_func(p, x):
   m, c = p
   return m*x+c

linear_model = Model(linear_func)

data = RealData(x, y, sx=x_err, sy=y_err)

odr = ODR(data, linear_model, beta0=[0.4, 0.4])

out = odr.run()

out.pprint()

Beta: [0.94999996 0.24999994]
Beta Std Error: [0.13228754 0.36228459]

But as you can see, the result is different from the curve_fit above with absolute_sigma=True.

Using the same data with curve_fit and absolute_sigma=False leads to the same results as the ODR-Fit . 

popt, pcov = curve_fit(func, x, y,sigma=yerr,absolute_sigma=False)

print(popt)
print(np.sqrt(np.diag(pcov)))

[0.95 0.25]
[0.13228756 0.36228442]

So I guess the ODR-Fit is not really taking care of my absolute errors as curve_fit and absolute_sigma=True does. Is there any way to do that or am I missing something?

eax
  • 71
  • 1
  • 4

1 Answers1

4

The option absolute_sigma=True in curve_fit() gives out the real covariance matrix in that sense that np.sqrt(np.diag(pcov)) yields the 1-sigma standard deviation as error as defined in, e.g., Numerical Recipes, i.e., it could have a unit, like meter or so. A very helpful summary comes with kmpfit

ODR, however, gives out the standard error derived from a scaled covariance matrix. See here, How to compute standard error from ODR results? or the example below.

This scaling by ODR is performed such that a reduced chi2 - calculated with the input weights being scaled in the same manner - yields approx 1. Or from the curve_fit docs:

This constant is set by demanding that the reduced chisq for the optimal parameters popt when using the scaled sigma equals unity.

The remaining question is now: What does sd_beta really mean? It's the standard error, see e.g., Standard error of mean versus standard deviation (there may exist conditions where the magnitude of both are the same. See comments below) see also the preceding discussion here: https://mail.python.org/pipermail/scipy-user/2013-February/034196.html

Now, via scaling pcov with the reduced chi2 one obtains for the parameters errors the same output of a) curve_fit(..., absolute_sigma=False) and b) ODR which are a priory relative errors:

# calculate chi2
residuals = (y - func(x, *popt))
chi_arr =  residuals / yerr
chi2_red = (chi_arr**2).sum() / (len(x)-len(popt))
print('red. chi2:\t\t\t\t', chi2_red)
print('np.sqrt(np.diag(pcov) * chi2_red):\t', np.sqrt(np.diag(pcov) * chi2_red))

yielding:

red. chi2:                           8.75
np.sqrt(np.diag(pcov) * chi2_red):   [0.13228756 0.36228442]

Note, that, however, for curve_fit(..., absolute_sigma=True) and ODR the covariance matrices pcov from curve_fitand cov_beta from the ODR output are still the same. Here, the disadvantage of back- and forth rescaling becomes perceivable!

The relative error now of course implies, that if the input errors are scaled, the magnitude of the relative errors remains:

# scaled uncertainty
yerr = np.asfarray([0.1, 0.1, 0.1, 0.1]) * 2

popt, pcov = curve_fit(func, x, y, sigma=yerr, absolute_sigma=False)

print('popt:\t\t\t', popt)
print('np.sqrt(np.diag(pcov)):\t', np.sqrt(np.diag(pcov)))

with the same output of:

popt:                    [0.95 0.25]
np.sqrt(np.diag(pcov)):  [0.13228756 0.36228442]
red-isso
  • 305
  • 1
  • 7
  • I am not sure your conclusion is right. If you use ODR in leastsq mode to fit a constant with some normal noise you will see that sd_beta matches the standard deviation of the data (the same for curve_fit with absolute_sigma=False). Have a look here: https://stackoverflow.com/a/44686138/2588210 - sd_beta in fact is the standard deviation of the parameter. – Christian K. Jul 13 '20 at 21:21
  • Yes, ODR in leastsq mode will give the same `sd_beta` as `np.sqrt(diag(pcov))` (s.a.) for `curve_fit(.., absolute_sigma=False)`. For both the cov. matrix is scaled by the red. chi2: For `curve_fit(...=False/ True)` the output of `pcov` changes, ODR gives the *bare* cov. matrix and does the scaling before `sd_beta` is calculated (see entry on VCVI p. 85 in the odrpack_guide.pdf. However, this is the standard *error* (not dev.) (see outp. of `out.pprint()`. The ODR docs sometimes also talk about the sd (see vector SDI). I think it's not really a bug but more an unlucky mingling of conventions. – red-isso Jul 14 '20 at 20:32
  • As a quick addon to @ChristianK. 's comment, yes they could actually reflect the real standard deviation but only if the sample size is much larger than the degree's of freedom. I admit, I would have to do the full math to work this out completely... Maybe someone here, has already done this ;) – red-isso Jul 15 '20 at 06:02
  • It is basically all explained in here: https://en.wikipedia.org/wiki/Ordinary_least_squares For small degrees of freedom one has to consider an additional increase of the uncertainty according to Student's statistics. – red-isso Jul 17 '20 at 07:42
  • So what should I do if I want ODR to produce the same results as curve_fit(..., absolute_sigma=True)? – Astral Cai Aug 03 '23 at 14:25