How To Return The Fit Error In Python Curve_fit

I'm trying to fit function to a data set of an experiment using python. I can get a really good approximation and the fit looks pretty good, but the error given for the parameters

Solution 1:

If you had printed out the full fit report from lmfit (or properly untangled to components of the covariance matrix from curve_fit) you would see that the parameters a and b are 100% correlated.

Basically, this is the fitting algorithm telling you that your data is not described well by you model and that you don't need that many parameters (or perhaps these parameters and this model) to describe your data.

Indeed, if you plot the data, there is a gentle slope up. Your function is the sum of three different terms:

   (4.9+a*x+b)*x
 + (a*x+b)*((mo/q)-x)*log(1-x/(mo/q)) 
 - 0.5*9.8*x*x

There are a couple of things to note:

1. mo and q only appear together, and as mo/q. They will not be independent.
2. a and b only appear together, and in the same form in multiple places
3. there are 2 purely quadratic terms in x, one of them hardwired.
4. the logarithmic term also has a quadratic prefactor. Importantly, the x data does not vary by more than 1 order of magnitude, so that log term won't really vary by much, giving mostly a third quadratic term. (as an aside: if you're taking logs, you should ensure that the argument is actually positive -- log(1-x*a) is asking for trouble)

To summarize: your model is too complicated for your data.

I very much doubt you need to log term at all. As it turns out, you can get a pretty good fit with simple parabolic model:

import numpy as np
import matplotlib.pyplot as plt
from lmfit.models import ParabolicModel

x = np.array([0.0333, 0.0667, 0.1, 0.133, 0.167, 0.2, 0.233 , 0.267, 0.3 ,
          0.333, 0.367, 0.4, 0.433, 0.467, 0.5, 0.533 , 0.567 , 0.6 ])

y = np.array([0.104, 0.249 , 0.422, 0.6, 0.791, 1.0, 1.23, 1.47, 1.74,
          2.02, 2.33, 2.64, 2.99, 3.34, 3.71, 4.08, 4.47, 4.85 ])

qmodel = ParabolicModel()
result = qmodel.fit(y, x=x, a=1, b=2, c=0)
print(result.fit_report())
fitlabel = "fit: a=%5.3f, b=%5.3f, c=%5.3f" % (result.params['a'].value,
                                               result.params['b'].value,
                                               result.params['c'].value)
plt.plot(x, y, label='Daten')
plt.plot(x, result.best_fit, label=fitlabel)
plt.xlabel("t [s]")
plt.ylabel("z [m]")
plt.legend(loc='upper left')
plt.title("Anpassung vor Zeitpunkt T (Model: a*x^2+b*x+c)")
plt.show()

which will give a report of

[[Model]]
    Model(parabolic)
[[Fit Statistics]]
    # fitting method   = leastsq
    # functionevals   = 9
    # datapoints      = 18
    # variables        = 3
    chi-square         = 0.00298906
    reducedchi-square = 1.9927e-04
    Akaikeinfocrit   = -150.657052
    Bayesianinfocrit = -147.985936
[[Variables]]
    c: -0.02973853 +/- 0.01120090 (37.66%) (init = 0)
    b:  3.67707491 +/- 0.08142567 (2.21%) (init = 2)
    a:  7.51540814 +/- 0.12492370 (1.66%) (init = 1)
[[Correlations]] (unreported correlations are < 0.100)
    C(b, a) = -0.972
    C(c, b) = -0.891
    C(c, a) =  0.785

and a plot of

Solution 2:

The data and the formula came from a physics experiment so I had to use the function. I could determine the actual value of mo and I also had to think about the meaning of the different variables. First in the formula "vo" and "ve" have to be equal, because the data "starts" after a period of acceleration and also because in the original fit "ve" was assumed to be a linear function but physically it should be a constant velocity. That's why in the following code vo+ve=ve. That left the fit function with two parameters which could be determined easily and with a small error.

def func(x, ve, q):
    return (2*ve*x) + ve*((0.65/q) - x)*log(1 - (q*x)/0.65) - 0.5*9.8*x*x

popt, pcov = curve_fit(func,x,y,bounds=((-100, 0), (1500, 1.5)))
plt.plot(x, func(x, *popt), color='orange', label='Anpassung: $v_e$=%5.3f, q=%5.3f' % 
tuple(popt))
plt.plot(x,y,':', label='Daten')

print(pcov)

plt.grid(True)
plt.legend(loc='upper left')
plt.xlabel("t [s]")
plt.ylabel("z [m]")
plt.title('Anpassung vor Zeitpunkt T', )
plt.savefig('anpass.pdf')
plt.savefig('anpass.png')
plt.show()

Which gives the covariance matrix:

[[ 0.01807039 -0.00305967]
 [-0.00305967  0.00065509]]

And the Fit:

Fitted Function

Which are realistic values for the function.

Thank you for your help and I hope this post is useful for someone someday!