12

Is it possible to obtain the value of the chi squared as a direct output of scipy.optimize.curve_fit()?

Usually, it is easy to compute it after the fit by squaring the difference between the model and the data, weighting by the uncertainties and summing all up. However, it is not as direct when the parameter sigma is passed a 2D matrix (the covariance matrix of the data) instead of a simple 1D array.

Are really the best-fit parameters and its covariance matrix the only two outputs that can be extracted from curve_fit()?

Stefano
  • 359
  • 1
  • 5
  • 16

1 Answers1

11

It is not possible to obtain the value of chi^2 from scipy.optimize.curve_fit directly without manual calculations. It is possible to get additional output from curve_fit besides popt and pcov by providing the argument full_output=True, but the additional output does not contain the value of chi^2. (The additional output is documented e.g. at leastsq here).

In the case where sigma is a MxM array, the definition of the chi^2 function minimized by curve_fit is slightly different. In this case, curve_fit minimizes the function r.T @ inv(sigma) @ r, where r = ydata - f(xdata, *popt), instead of chisq = sum((r / sigma) ** 2) in the case of one dimensional sigma, see the documentation of the parameter sigma. So you should also be able to calculate chi^2 in your case by using r.T @ inv(sigma) @ r with your optimized parameters.

An alternative would be to use another package, for example lmfit, where the value of chi square can be directly obtained from the fit result:

from lmfit.models import GaussianModel

model = GaussianModel()

# create parameters with initial guesses:
params = model.make_params(center=9, amplitude=40, sigma=1)  

result = model.fit(n, params, x=centers)
print(result.chisqr)
jdamp
  • 1,379
  • 1
  • 15
  • 22
  • `lmfit` is indeed a great package. Unfortunately, it does not offer the option to pass an MxM matrix to `sigma` but only a 1D array. Hope this can be implemented in the future. – Stefano Oct 01 '18 at 14:26