Comparison of curve_fit and scipy.odr - absolute sigma

Question

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?

Answer 1

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]