In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
data
In [3]:
x=np.array([0.1,2.1,3.3,4.9,5.6,6.5,7.7,8,10.5,13]) #data x
ex=np.array([0.9,0.3,0.1,0.1,0.5,0.5,0.1,0.2,0.9,2]) #chyby x
y=np.array([19,23,25,31,35,43,46,44,52,56]) #data y
ey=np.array([2,2,2,2,1,5,3,4,2,1]) #chyby y
In [4]:
Ndata=np.size(x)
print(f"počet naměřených hodnot: {Ndata}")
počet naměřených hodnot: 10
modelová funkce
In [6]:
def primka(x,a,b):
return a*x+b
$\chi^2$, které budeme minimalizovat (minimilizujeme kvadrát funkce chisq)
In [8]:
def chisq(theta,x,y,ex,ey):
return (y-theta[0]-theta[1]*x)/np.sqrt(ey**2+(ex*theta[1])**2)
In [9]:
#počáteční odhady parametrů a,b
theta_initial=([1,1])
#minimalizace kvadrátu funkce chisq nelineární metodou nejmenších čtverců
res=leastsq(chisq,theta_initial,args=(x,y,ex,ey),full_output=True)
#vektor res[0] obsahuje odhady parametrů a,b
a=res[0][1]
b=res[0][0]
#matice res[1] je kovarianční matice
e_a=np.sqrt(res[1][1][1]) #chyba parametru a
e_b=np.sqrt(res[1][0][0]) #chyba parametru a
cov_ab=res[1][0][1] #cov(a,b)
print("modelová funkce y=ax+b")
print(f"a = {a:3.1f} +/- {e_a:3.1f}")
print(f"b = {b:3.0f} +/- {e_b:3.0f}")
print(f"cov(a,b) = {cov_ab:3.1f}")
modelová funkce y=ax+b a = 3.6 +/- 0.4 b = 15 +/- 2 cov(a,b) = -0.6
In [10]:
y_fit=a*x+b
extrapolace
In [12]:
x0=-b/a
#obecná metoda přenosu chyb
e_x0=np.sqrt((b/a**2)**2*e_a**2+(1/a)**2*e_b**2-2*(b/a**2)*(1/a)*cov_ab)
print(f"průsečík s osou x: {x0:3.1f} +/- {e_x0:3.1f}")
průsečík s osou x: -4.2 +/- 0.9
In [13]:
plt.errorbar(x,y,ey,ex,lw=0,elinewidth=2,capsize=5,marker='o')
plt.plot(x,y_fit,c='red',lw=3)
plt.errorbar(x0,0,0,e_x0,lw=0,elinewidth=2,capsize=5,marker='*',c='green')
plt.plot([x0,x[0]],[0,y_fit[0]],c='green',ls='dashed',lw=2)
plt.xlabel("x")
plt.ylabel("y")
plt.show()
In [ ]: