Problems with curve_fit fitting highly correlated data

Question

For my bachelor thesis, I am working on a project where I want to perform a fit to some data. The problem is a bit more complex, but I tried to minimize the problem here:

We have three data points (very little theory data is available), but these points are highly correlated.

Using curve_fit to fit these points, we get a horrible fit result, as you can see in this picture. (The fit could be easily improved by altering the fit parameters by hand).

Our fit results with correlations (blue) and with neglected correlations (orange):

The results get better when we use more parameters (as the fits essentially behave like solves by then).

My question: Why does this behaviour happen? (We use our own least-squares algorithm for our specific problem, but it suffers from the same problem). Is it a numerical problem, or is there any good reason for curve_fit to show this solution?

I would be very happy to have a good explanation to say why we can't use "only 2" parameters to fit these highly correlated 3 datapoints.

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
plt.rcParams['lines.linewidth'] = 1

y = np.array([1.1994, 1.0941, 1.0047])
w = np.array([1, 1.08, 1.16])
cor = np.array([[1, 0.9674, 0.8812],[0.9674, 1, 0.9523], [0.8812, 0.9523, 1]])
s = np.array([0.0095, 0.0104, 0.0072])

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

cov = np.zeros((3,3))
for i in range(3):
    for j in range(3):
        cov[i,j] = cor[i,j] * s[i] * s[j]

A1, B1 = curve_fit(f, w, y, sigma=cov)
A2, B2 = curve_fit(f, w, y)

plt.plot(w, f(w, *A1))
plt.plot(w, f(w, *A2))

plt.scatter(w, y)
plt.show()

Answer 1

It is not a numerical problem. The "problem" is that the off-diagonal terms of your covariance matrix are all positive and relatively large. These determine the correlations among the errors in the fit, so if all the terms are positive, you are saying that all the errors are positively correlated. If one is large, the others will tend to also be large with the same sign.

Here's an example similar to yours, with the covariance matrix

        [2.0  1.3  0.0]
sigma = [1.3  2.0  1.3]
        [0.0  1.3  2.0]

(The condition number of this matrix is 23.76, so we shouldn't expect any numerical problems.)

While the covariance between the first and third points is 0, it is 1.3 between the first and second, and between the second and third, and 1.3 is a relatively large fraction of the variances, which are all 2. So it will not be surprising if all the errors in the fitted model have the same sign.

This script does a fit of three points and plots the data and the fitted line.

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt


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


x = np.array([1, 2, 3])
y = np.array([2, 0.75, 0])
sig = np.array([[2.0, 1.3, 0.0],
                [1.3, 2.0, 1.3],
                [0.0, 1.3, 2.0]])

params, pcov = curve_fit(f, x, y, sigma=sig)

y_errors = f(x, *params) - y

plt.plot(x, y, 'ko', label="data")
plt.plot(x, f(x, *params), linewidth=2.5, label="fitted curve")
plt.vlines(x, y, f(x, *params), 'r')

for k in range(3):
    plt.annotate(s=r"$e_{%d}$" % (k+1), xy=(x[k]-0.05, y[k]+0.5*y_errors[k]), ha='right')

plt.xlabel('x')
plt.ylabel('y')
plt.axis('equal')
plt.grid()
plt.legend(framealpha=1, shadow=True)
plt.show()

As you can see in the plot, all the errors have the same sign.

We can confirm this reasoning by considering another covariance matrix,

        [ 2.0   1.3  -1.0]
sigma = [ 1.3   2.0  -1.3]
        [-1.0  -1.3   2.0]

In this case, all the off-diagonal terms are relatively large in magnitude. The covariance between the first and second errors is positive, and it is negative between the second and third and between the first and third. If these off-diagonal terms are large enough relative to the variances, we should expect the signs of the errors of the first two points to be the same, while the third error will have the opposite sign of the first two.

Here's the plot generated by the script when sig is changed to the above matrix:

The errors show the expected pattern.