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()
No description has been provided for this image
In [ ]: