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.
