In [37]:
p=peri[2]
x,y=array(dset_min[p]).T
plot(x,y,'+')
amat=array([ones_like(x),x,exp(-p3fix*x)]).T
acov=linalg.inv(amat.T.dot(amat))
pars=acov.dot(amat.T.dot(y))
plot(x,amat.dot(pars))
xlabel("cas [s]")
ylabel("amplit. [mV]")
pars
Out[37]:
Problém vychází z úlohy demonstrující elektromagnetickou indukci (č. 1 z fyzikálního praktika 2). V přiložených souborech jsou hodnoty maxima a minim kmitů tlumeného kyvadla (v abs. hodnotě - amplituda $A$) s příslušnými časy (levý sloupec, v sekundách) pro různé tlumící odpory. S jejich pomocí nalezněte parametry odpovídající kombinovanému modelu $$A(t)=p_0 + p_1 t + p_2 \exp (-p_3 t)$$ pro všechny sady v zadání; nelineární parametr volte fixní $p_3=6.4\ 10^{-3} s^{-1}$. Rozhodněte, zda body s největšími (absolutními) hodnotami reziduí jsou hrubé chyby, ty případně vyřaďte a regresi opakujte.
Stanovte nejistoty a korelační koeficienty mezi parametry. Můžete se pokusit zlepšit fit modelu uvolněním parametru $p_3$ (pokud máte nástroje na nelineárním optimalizaci); problémem je zde značná korelace mezi lineárním a exponenciálním modelem, který vede k nárůstu nejistot a fluktuacím parametru $p_1$, který je potřeba dále.
Parametr $p_1$ dle teorie má být nepřímo úměrný součtu zátěžového odporu $R_p$ a vnitřního odporu cívky $R_c$. Z lineární regrese mezi $R_p$ a $1/p_1$ (s příslušnými nejistotami) určete 95% konfidenční interval pro odpor $R_c$
případně stanovte jeho horní mez, za předpokladu, že $R_c$ musí být nezáporné (pro určení mezí ve fyzikálně ohraničených oblastech lze těžit z postupů v druhé části kap.9 knihy F. Jamese, použít můžete ale i nejjednodušší Bayesovský přístup demonstrovaný v zápisníku zde).
Jaká je reziduální suma čtverců (minimalizovaná hodnota) u této závěrečné regrese a na co z ní lze usuzovat?
Data najdete na adrese http://physics.muni.cz/praktika/static/docs/TPX/data/2018/zap_xxx.dat kde xxx je vaše id. číslo (v opačném případě se ozvěte opět na munz@physics.muni.cz).
%matplotlib inline
from matplotlib.pyplot import *
from numpy import *
from scipy import stats
indir="C:/Users/Admin/Documents/Vyuka/2018/"
#indir="/home/limu/Dokumenty/Vyuka/MMZM/zapocet/2018/"
adata=open(indir+"zap_162880.dat").readlines()
def read_exper(adata):
dset_min={}
dset_max={}
cmin=True
for a in adata:
if a.find('#')>=0:
if a.find('ohm')>0:
cat=int(a[2:a.find('ohm')])
elif a.find('max')>0:
cmin=False
dset_max[cat]=[]
else:
cmin=True
dset_min[cat]=[]
elif a.find('.')>0:
val=[float(b) for b in a.split()]
if len(val)>1:
if cmin: dset_min[cat].append(val[:2])
else: dset_max[cat].append(val[:2])
return dset_min,dset_max
dset_min,dset_max=read_exper(adata)
r_zatez=list(dset_min.keys())
p3fix=6.4e-3
r_zatez
p=peri[2]
x,y=array(dset_min[p]).T
plot(x,y,'+')
amat=array([ones_like(x),x,exp(-p3fix*x)]).T
acov=linalg.inv(amat.T.dot(amat))
pars=acov.dot(amat.T.dot(y))
plot(x,amat.dot(pars))
xlabel("cas [s]")
ylabel("amplit. [mV]")
pars