# What is a good library in Python for correlated fits in both the $x$ and $y$ data?

I have $$x$$ and $$y$$ data, both of which have their own covariance matrices. scipy.optimize.curve_fit will accept a covariance matrix for the $$y$$ data, called sigma, but not for $$x$$ data. I would like to pass a covariance matrix for my $$x$$ data also, since it is also correlated.

I found scipy.odr, however the 2d-array option for we and wd, as well as for covx and covy do not work (I have submitted an issue on github).

I have written my own ODR code in python, however it's quite slow. The ODRPACK that scipy.odr uses I imagine is quite fast since it is FORTRAN.

I thought I would probe to see if people knew of a nice option, in place of me using the slow code I wrote.

Is there an alterntative to scipy.odr for fitting correlated $$x$$ and $$y$$ data using python?

Thanks.

• Do you have some references to do that mathematically? – David Aug 9 '19 at 4:25
• @David I will write a answer and post my code. Maybe someone can do better/find something better. – kηives Aug 9 '19 at 20:05
• @David I posted it, and the reference for ODRPACK. – kηives Aug 16 '19 at 20:23

Here is the current code I am using to do correlated fits in the $$x$$ and $$y$$ directions. I wrote it from this reference. I adapted it from section 1.A.i and 4.B

#!/anaconda/bin/python

import numpy as np
from scipy.optimize import least_squares

class Fitter(object):
"""
Orthogonal distance regression (ODR) fitter object.
"""

def __init__(self,):
pass

def Model(self, function):
"""
Set the model to fit to.
"""
self.model = function

def Beta0(self, init_beta):
"""
Set the initial parameters.
"""
self.beta0 = np.array(init_beta)
self.bs = len(init_beta)

def Data(self, xdata, ydata, sigmax=None, sigmay=None, covx=None, covy=None):
"""
Input the x and y data, as well as the errors if they exist.
"""
if (len(xdata) != len(ydata)):
raise ValueError("x and y data must have same shape")
self.x = xdata
self.y = ydata
self.xs = len(xdata)
self.ys = len(ydata)
self.delta = np.zeros(shape=(self.xs,))
if (sigmax is not None) and (sigmay is not None):
Omega = np.block([[np.diag(1./sigmay**2), np.zeros((self.ys, self.ys))],
[np.zeros((self.xs, self.xs)), np.diag(1./sigmax**2)]])
elif (sigmax is not None) and (sigmay is None):
Omega = np.block([[np.eye(self.xs), np.zeros((self.ys, self.ys))],
[np.zeros((self.xs, self.xs)), np.diag(1./sigmax**2)]])
elif (sigmax is None) and (sigmay is not None):
Omega = np.block([[np.diag(1./sigmay**2), np.zeros((self.ys, self.ys))],
[np.zeros((self.xs, self.xs)), np.eye(self.xs)]])
elif (covx is not None) and (covy is not None):
Omega = np.block([[np.linalg.inv(covy), np.zeros((self.ys, self.ys))],
[np.zeros((self.xs, self.xs)), np.linalg.inv(covx)]])
elif (covx is None) and (covy is not None):
Omega = np.block([[np.linalg.inv(covy), np.zeros((self.ys, self.ys))],
[np.zeros((self.xs, self.xs)), np.eye(self.xs)]])
elif (covx is not None) and (covy is None):
Omega = np.block([[np.eye(self.xs), np.zeros((self.ys, self.ys))],
[np.zeros((self.xs, self.xs)), np.linalg.inv(covx)]])
else:
Omega = np.eye(2*self.xs)

self.L = np.linalg.cholesky(Omega)
assert np.allclose(Omega, self.L.dot(self.L.transpose()))

def Residuals(self, z):
"""
Calculates the sum of residuals for ODR.
"""
para = z[:self.bs]
arguments = z[self.bs:]

rx = arguments
ry = (self.model(self.x + rx, *para) - self.y)
r = np.hstack((ry, rx))

return np.dot(r, self.L)

def Run(self,):
"""
Runs the fitter.
"""
self.whole = np.append(self.beta0, self.delta)
self.out = least_squares(self.Residuals, self.whole, method='lm')
self.chi2 = np.sum(self.out.fun**2) / float(len(self.x)-len(self.beta0))


One passes in a function as the fitting form to Model, sets the initial guess with Beta0, passes in the data and errors (or covariances) with Data, then runs to find the best fit parameters with Run.

Mathematically what I think ODR does, and what I wrote here is trying to minimize weighted least squares. If $$C^{(x)}$$ and $$C^{(y)}$$ are covariance matrices for $$x$$ and $$y$$ data (or their diagonals for the squared errors) we try to minimize the expression, $$\sum_{i,j} (f(\bar{x}_i, \beta_{k}) - y_i) {C^{(y)}}^{-1}_{ij} (f(\bar{x}_j, \beta_{k}) - y_j) + (\bar{x}_{i} - x_{i}) {C^{(x)}}^{-1}_{ij} (\bar{x}_{j} - x_{j})$$ where $$\beta_{k}$$ are the fit parameters, $$x_{i}$$ and $$y_{i}$$ are the the data, $$f$$ is some fit function, and $$\bar{x}_{i}$$ is some ideal $$x$$ value that minimizes this expression. So in some sense one is solving for $$\beta$$ and $$\bar{x}$$. Alternatively we can replace $$\delta_{i} \equiv (\bar{x}_{i} - x_{i})$$, which gives $$\bar{x}_{i} = \delta_{i} + x_{i}$$, giving a new expression, $$\sum_{i,j} (f(\delta_{i} + x_{i}, \beta_{k}) - y_i) {C^{(y)}}^{-1}_{ij} (f(\delta_j + x_j, \beta_{k}) - y_j) + \delta_i {C^{(x)}}^{-1}_{ij} \delta_j$$ Putting the $$x$$ and $$y$$ inverse covariances in a block-diagonal matrix and taking Cholesky decomposition (the covariance matrices have to be Hermitian at least), This is dotted with the $$x$$ and $$y$$ residuals. This final vector when squared and summed (or dotted with itself) gives the sum of the squares which is what the scipy.least_squares wants.

• Could you please add fake data to test your code? – David Aug 20 '19 at 5:21