Chi square to polynomial fit

Question

I have this polynomial fit to my data.

fig3 = plt.figure(3)

for dataset in [Bxfft]:
    dataset = np.asarray(dataset)
    freqs, psd = signal.welch(dataset, fs=266336/300, window='hamming', nperseg=8192)
    plt.semilogy(freqs, psd/dataset.size**0, color='r')

# Polynomial 5th grade
def line(freqs, a, b, c, d, e, f):
    return a*freqs**5 + b*freqs**4 + c*freqs**3 + d*freqs**2 + e*freqs + f

popt, pcov = curve_fit(line, freqs, np.log10(psd))
plt.semilogy(freqs, 10**line(freqs, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5]), 'black')

This is what I get:

I would like to count chi square, but to be honest I have no idea how to do it. I was able to do something like this, but I think this is wrong.

chisquare = chi(popt)
print chisquare
Power_divergenceResult(statistic=-0.4318298090941465, pvalue=1.0)

Answer 1

chi-square is generally defined as sum-of-squares of data-fit. For your example that would be:

 best_fit = 10**line(freqs, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5])
 chi_square = ((psd - best_fit)**2).sum()

Allow me to also suggest using lmfit (https://lmfit.github.io/lmfit-py/) for curve-fitting as an alternative to curve_fit that handles many such chores. Using lmfit, your example might look like this:

from lmfit import Model
def line(freqs, a, b, c, d, e, f):
    return a*freqs**5 + b*freqs**4 + c*freqs**3 + d*freqs**2 + e*freqs + f

# turn your model function into an lmfit Model
pmodel = Model(line)

# make parameters with initial guesses. note that parameters are
# named 'a', 'b', 'c', etc based on your `line` function, not numbered.
params = pmodel.make_params(a=1, b=-0.5, c=0, d=0, e=0, f=0)

# fit model to data with these parameters, specify independent variable 
result = pmodel.fit(np.log10(psd), params, freqs=freqs)

# this result has chi-square calculated:
print('Chi-square = ', result.chisqr)

# print report with fit statistics, parameter values and uncertainties
print(result.fit_report()) 

# plot, using best-fit in result
plt.semilogy(freqs, psd, color='r')
plt.semilogy(freqs, 10**(result.best_fit), color='k')
plt.show()