import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2 

modelová funkce: parabola $y = \theta_0 + \theta_1 x + \theta_2 x^2$

simulace dat

theta_true=([0.3,-5,1.5])
Ndata=100
x=np.linspace(0,20,Ndata)
y=np.polyval(theta_true,x)

rozšumění

ey=np.random.uniform(0,5,Ndata) #chyby y
y=y+np.random.normal(0,ey,Ndata) #rozšuměné hodnoty y

fit metodou nejmenších čtverců

theta_fit,cov=np.polyfit(x,y,2,w=1/ey,cov='unscaled')
print("nafitované parametry paraboly:")
print(f"theta_0 = {theta_fit[0]:5.3f} +/- {np.sqrt(cov[0,0]):5.3f}")
print(f"theta_1 = {theta_fit[1]:5.3f} +/- {np.sqrt(cov[1,1]):5.3f}")
print(f"theta_2 = {theta_fit[2]:5.3f} +/- {np.sqrt(cov[2,2]):5.3f}")

nafitované parametry paraboly:
theta_0 = 0.300 +/- 0.002
theta_1 = -4.995 +/- 0.024
theta_2 = 1.508 +/- 0.083

plt.errorbar(x,y,ey,lw=0,elinewidth=1,capsize=5,marker='o',label='data')
plt.plot(x,np.polyval(theta_fit,x),c='red',lw=3,label='fit')
plt.xlabel("x")
plt.ylabel("y")
plt.show()

$\chi^2$ test

chi2_exp=np.sum(((y-np.polyval(theta_fit,x))/ey)**2)
print(f"chi2 = {chi2_exp:8.2f}")
Ndf=Ndata-3 
print(f"počet stupňů volnosti = {Ndf}")
P=1-chi2.cdf(chi2_exp,Ndf)
print(f"P hodnota: {P}") 
print(f"redukované chi2 = {chi2_exp/Ndf:5.3f}")

chi2 =   106.07
počet stupňů volnosti = 97
P hodnota: 0.24850627209966392
redukované chi2 = 1.093