I’m trying to make a data-fit with odr, bit it’s not showing the fit. Any ideas why?
f,ax = plt.subplots(1)
ax.errorbar(starttijden,afstanden,xerr=sig_tijd,yerr=sig_afstand,fmt='ko',label='Datapunten')
#ax.show()
A_start=0
B_start=40
## (1) Definieer een Python-functie die het model bevat, in dit geval een
## rechte lijn
## B is een vector met parameters, in dit geval twee (A = B[0], B = B[1])
## x is de array met x-waarden
def f(B, x):
return B[0] + B[1]*x
## (2) Definieer het model-object om te gebruiken in odr
odr_model = odr.Model(f)
## (3) Definieer een RealData object
## Een RealData-object vraagt om de onzekerheden in beide richtingen.
## !! De onzekerheid in de x-richting mag ook nul zijn (dan mag je
## sx=0 weglaten), maar dan moet bij onderdeel (4)/(4a) wel gekozen
## worden voor een kleinste-kwadratenfit !!
odr_data = odr.RealData(starttijden,afstanden,sx=sig_tijd,sy=sig_afstand)
## (4) Maak een ODR object met data, model en startwaarden
## Je geeft startwaarden voor parameters mee bij keyword beta0
odr_obj = odr.ODR(odr_data,odr_model,beta0=[A_start,B_start])
## (4a) Stel in op kleinste kwadraten (optioneel)
## Als de onzekerheden in de x-richting gelijk zijn aan nul, dan faalt de
## default-aanpak van ODR. Je moet dan als methode een kleinste-kwadraten-
## fit instellen. Dit gaat met het volgende commando (nu uit):
#odr_obj.set_job(fit_type=2)
## (5) Voer de fit uit.
## Dit gebeurt expliciet door functie .run() aan te roepen.
odr_res = odr_obj.run()
## (6) Haal resultaten uit het resultaten-object:
# (6a) De beste schatters voor de parameters
par_best = odr_res.beta
# (6b) De (EXTERNE!) onzekerheden voor deze parameters.
par_sig_ext = odr_res.sd_beta
# (6c) De (INTERNE!) covariantiematrix.
par_cov_ext = odr_res.cov_beta
# (6d) De chi-kwadraat en de gereduceerde chi-kwadraat van deze fit.
chi2 = odr_res.sum_square
chi2red = odr_res.res_var
# (6e) Een compacte weergave van de belangrijkste resultaten als output.
odr_res.pprint()
# Hier plotten we ter controle de fit met de dataset.
xplot=np.arange(-7.5,-5.5,100)
ax.plot(xplot,par_best[0] + par_best[1]*xplot,'r-',label='Beste fit')
ax.axes.set(xlabel='Starttijd', ylabel='Afstand tot epicentrum (km)',title="Starttijd van de aardbeving tegen de afstand met de beste fit")
ax.grid(alpha=0.7)
ax.legend()
Don’t pay attention to the comments