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)
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()