In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
definice modelové funkce
In [3]:
def primka(theta,x):
return theta[0]+theta[1]*x
definice funkce, která se bude minimizovat.
Procedura leastsq minimizuje kvadrát této funkce
In [5]:
def chi(theta,x,y,ex,ey):
return (y-primka(theta,x))/np.sqrt(ey**2+(ex*theta[1])**2)
Načtení dat
In [7]:
data=np.loadtxt('data_ex_ey.txt')
x=data[:,0]
ex=data[:,1]
y=data[:,2]
ey=data[:,3]
fit
In [9]:
theta_initial=np.array([0,1]) #počáteční hodnoty parametrů přímky
res=leastsq(chi,theta_initial,args=(x,y,ex,ey),full_output=True)
print("primka y = ax + b")
print(f"a = {res[0][0]} +/- {np.sqrt(res[1][0][0])}")
print(f"b = {res[0][1]} +/- {np.sqrt(res[1][1][1])}")
primka y = ax + b a = -18.157039602263165 +/- 6.545554913976975 b = 3.9076737666627337 +/- 0.09500114466485858
Vykreslení výsledků
In [11]:
plt.errorbar(x,y,xerr=ex,yerr=ey,lw=0,elinewidth=1,capsize=5,marker='o',label='data')
plt.plot(x,primka(res[0],x),c='red',label='fit')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()
In [ ]: