<jupyter_start><jupyter_text>Zápočtová úloha do F7270 Veronika Dočkalová Naměřené hodnoty signálu detektoru částic v čase budeme fitovat exponenciální funkcí $y=A\exp(−t\lambda)$. Pro účely lineární regrese provedeme logaritmickou transformaci naměřených dat i fitované funkce $\ln(y)=\ln(A)-t\lambda$ a budeme hledat hodnoty parametrů $\ln(A)$ a $\lambda$. Nejdříve naimportujeme všechny potřebné balíčky. Do proměnných x,y uložíme naměřené hodnoty, z nichž ovšem odstraníme všechny nulové body, a spočítáme $\ln(y)$<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
mereni=np.loadtxt("sada_424108_300ks.txt")
pole = []
for x in mereni:
    if x[1] != 0:
        pole.append(x)
pole = np.array(pole)
x,y=data=np.transpose(pole)
ln_y=np.log(y)<jupyter_output><empty_output><jupyter_text>Dále vytvoříme sloupcovou modelovou matici M (2 x $n$), jejíž první řádek tvoří jedničky a druhý hodnoty -$x$. Z ní zkonstruujeme Hessián H=M$^{T}$M.<jupyter_code>M=np.array([np.ones(len(x)), -x]).T
H=M.T.dot(M)
H<jupyter_output><empty_output><jupyter_text>Nyní můžeme dosadit do vzorce pro výpočet matice parametrů $\hat{\theta}$=H$^{-1}$M$^T$Y, kde Y je matice (vektor) hodnot $\ln(y)$<jupyter_code>H_inv=np.linalg.inv(H)
theta=H_inv.dot(M.T.dot(ln_y))
theta<jupyter_output><empty_output><jupyter_text>a spočítat hodnoty parametrů $A$ a $\lambda$<jupyter_code>A=np.exp(theta[0])
L=theta[1]
A,L<jupyter_output><empty_output><jupyter_text>Nakonec vykreslíme data a výslednou fitovanou funkci do grafu<jupyter_code>plt.plot(x,y, "x", label="mereni")
plt.plot(x,A*np.exp(-x*L), label="fit")
plt.xlabel("$t$")
plt.ylabel("signal detektoru")
plt.legend(loc=0)<jupyter_output><empty_output><jupyter_text>Nejistoty parametrů určíme z kovarianční matice cov($\hat{\theta}$)=$\sigma^2$H$^{-1}$. Odhad hodnoty $\sigma^2$ získáme ze vztahu $\sigma^2=Q^2_{min}/(n-r)$, kde $Q^2_{min}=\Sigma_i(y_{m_i}-y_{f_i})^2$ je suma reziduí v minimu ($y_m$, $y_f$ jsou naměřené hodnoty, resp. hodnoty fitu) a $n-r$ je počet stupňů volnosti, tj. počet hodnot $-$ 2 (parametry).<jupyter_code>Q_min=np.sum((ln_y-(theta[0]-x*theta[1]))**2)
Q_min
sigma=Q_min/(len(x)-2)
cov_theta=sigma*H_inv
cov_theta<jupyter_output><empty_output><jupyter_text>Na diagonále kovarianční matice jsou rozptyly parametrů $\sigma^2(\ln A)$, $\sigma^2(\lambda)$, z nichž po odmocnění dostaneme nejistoty. Mimo diagonálu se nachází člen $\rho\sigma^2(\ln A)\sigma^2(\lambda)$, ze kterého můžeme určit korelační koeficient.<jupyter_code>sigma_lnA=np.sqrt(cov_theta[0,0])
sigma_L=np.sqrt(cov_theta[1,1])
corel=cov_theta[0,1]/(sigma_lnA*sigma_L)
sigma_lnA, sigma_L, corel<jupyter_output><empty_output><jupyter_text>Nejistotu $\sigma(\ln A)$ převedeme pomocí zákona o šíření nejistot na nejistotu, ze kterého je možné odvodit $\sigma(A)=A\sigma(\ln A)$. Finální hodnoty fitovaných parametrů jsou<jupyter_code>sigma_A=A*sigma_lnA
print("A = {} +- {}".format(A, sigma_A))
print("L = {} +- {}".format(L, sigma_L))<jupyter_output>A = 89.78324291461347 +- 3.875849500200264
L = 0.4778960556560432 +- 0.008439898343748392<jupyter_text>Výše uvedený postup je správný pro konstantní hodnoty nejistot měřených bodů. V případě, že každý bod má jinou nejistotu, je nutné do výpočtu zahrnout navíc váhovou matici W, která tuto skutečnost reflektuje. Váhová matice je diagonální matice, jejíž prvky jsou tvořeny převrácenou hodnotou rozptylu jednotlivých naměřených bodů.  a) Konstantní nejistoty Jelikož v průběhu zpracování měření provádíme logaritmicku transformaci, budeme muset pracovat s nekonstantními nejistotami i tehdy, když původní data mají nejistotu konstantní. V případě a) platí mezi rozptylem transformaovaných dat a původními hodnotami úměra  $\sigma^2 \approx y^{-2}$. Zkonstruujeme tedy váhovou matici W_a, na jejíž diagonále budou hodnoty $y^2$.<jupyter_code>W_a=np.eye(len(x))*(y**2)<jupyter_output><empty_output><jupyter_text>Váhovou matici je nutné zahrnout jak do výpočtu Hessiánu H=M$^{T}$WM, tak do výpočtu matice parametrů $\hat{\theta}$=H$^{-1}$M$^T$WY. Abychom nemuseli opisovat celý výše uvedený postup, definujeme si funkci, která provede celou analýzu dat automaticky.<jupyter_code>def analyza(data, W, title):
    x=data[0]
    y=data[1]
    ln_y=np.log(y)
    M=np.array([np.ones(len(x)), -x]).T
    H=M.T.dot(W.dot(M))
    H_inv=np.linalg.inv(H)
    theta=H_inv.dot(M.T.dot(W.dot(ln_y)))
    A=np.exp(theta[0])
    L=theta[1]
    Q_min=np.sum((ln_y-(theta[0]-x*theta[1]))**2)
    sigma=Q_min/(len(x)-2)
    cov_theta=sigma*H_inv
    sigma_lnA=np.sqrt(cov_theta[0,0])
    sigma_L=np.sqrt(cov_theta[1,1])
    corel=cov_theta[0,1]/(sigma_lnA*sigma_L)
    sigma_A=A*sigma_lnA
    plt.plot(x,y, "x")
    plt.plot(x,A*np.exp(-x*L), label=title)
    plt.xlabel("$t$")
    plt.ylabel("signal detektoru")
    plt.legend(loc=0)
    print("ln(A) = {} +- {}".format(theta[0], sigma_lnA))
    print("korel.koef. = {}".format(corel))
    print("Q_min= {}".format(Q_min))
    print("A = {} +- {}".format(A, sigma_A))
    print("L = {} +- {}".format(L, sigma_L))
    return A,L
title_a="pripad a)"
a=analyza(data, W_a, title_a) # do proměnné "a" si uložíme hodnoty vektoru theta<jupyter_output>ln(A) = 4.531650623564943 +- 0.0010375505321566893
korel.koef. = 0.6937799893165646
Q_min= 50.12905949567006
A = 92.91179698148564 +- 0.0964006844017747
L = 0.44895602717997285 +- 0.0006814487251427288<jupyter_text>Dále získáme rozdělení reziduí<jupyter_code>reziduum_a=y-(A*np.exp(-x*L))
plt.hist(reziduum_a, bins=20)<jupyter_output><empty_output><jupyter_text>a porovnáním šikmosti (skewness) a špičatosti (kurtosis) otestujeme hypotézu, zda rezidua odpovídají normálně rozděleným datům. Koeficient šikmosti spočítáme pomocí podílu třetího centrálního momentu a třetí mocniny směrodatné odchylky$$\mathrm{skew}=\frac{\mu_3}{\sigma^3}=\frac{E\left[ (x-\mu)^3 \right]}{\left( E\left[ (x-\mu)^2 \right]\right)^{3/2}}$$a její nejistotu jako $\sqrt{6/n}$. Koeficient špičatosti je dán vzorcem$$\mathrm{kurt}=\frac{\mu_4}{\sigma^4}-3=\frac{E\left[ (x-\mu)^4 \right]}{\left( E\left[ (x-\mu)^2 \right]\right)^2}-3$$s odchylkou $\sqrt{24/n}$<jupyter_code>expvalue_a=np.mean(reziduum_a)
CM4_a=np.mean((reziduum_a-expvalue_a)**4)
CM3_a=np.mean((reziduum_a-expvalue_a)**3)
CM2_a=np.mean((reziduum_a-expvalue_a)**2)
skew_a=np.sqrt(CM3_a**2/CM2_a**3)
sigma_skew_a=np.sqrt(6/len(x))
kurt_a=np.sqrt(CM4_a**2/CM2_a**4)-3
sigma_kurt_a=np.sqrt(24/len(x))
print("skew = {} +- {}".format(skew_a, sigma_skew_a))
print("kurt = {} +- {}".format(kurt_a, sigma_kurt_a))
print("2sigma = {} ".format(2*np.sqrt(CM2_a)))
print("3sigma = {} ".format(3*np.sqrt(CM2_a)))<jupyter_output>skew = 2.2812664058225205 +- 0.1426148065624585
kurt = 15.183455555907255 +- 0.285229613124917
2sigma = 10.640374915767218 
3sigma = 15.960562373650827<jupyter_text>Vidíme, že šikmost i špičatost se od nuly liší o více než $2\sigma$, tudíž zamítáme nulovou hypotézu, že rezidua mají normální rozdělení. Z měření tedy odděláme odlehlé hodnoty zatížené hrubou chybou, tj. takové, jejichž reziduum se od střední hodnoty reziduí liší o více než $3\sigma$, a postup zopakujeme.<jupyter_code>sigma3=3*np.sqrt(CM2_a)
pole_a1 = []
for z in pole:
   if (expvalue_a-sigma3) < (z[1] -(a[0]*np.exp(-z[0]*a[1]))) < expvalue_a+sigma3:
        pole_a1.append(z)
pole_a1 = np.array(pole_a1)
x1,y1=data1=np.transpose(pole_a1)
W_a1=np.eye(len(x1))*(y1**2)
title_a1="pripad a) bez odlehlych hodnot"
a1=analyza(data1, W_a1, title_a1)<jupyter_output>ln(A) = 4.509360443028548 +- 0.0011298384142490972
korel.koef. = 0.7156923771919719
Q_min= 49.14578134568504
A = 90.86368741869804 +- 0.10266128450596744
L = 0.4448481406870539 +- 0.0007175048418888641<jupyter_text>Dále definujeme funkci, která provede výpočet a analýzu reziduí.<jupyter_code>def rezidua(data, W):
    title="title"
    x=data[0]
    y=data[1]
    ln_y=np.log(y)
    M=np.array([np.ones(len(x)), -x]).T
    H=M.T.dot(W.dot(M))
    H_inv=np.linalg.inv(H)
    theta=H_inv.dot(M.T.dot(W.dot(ln_y)))
    reziduum=y-(np.exp(theta[0])*np.exp(-x*theta[1]))
    plt.hist(reziduum, bins="auto")
    expvalue=np.mean(reziduum)
    CM4=np.mean((reziduum-expvalue)**4)
    CM3=np.mean((reziduum-expvalue)**3)
    CM2=np.mean((reziduum-expvalue)**2)
    skew=np.sqrt(CM3**2/CM2**3)
    sigma_skew=np.sqrt(6/len(x))
    kurt=np.sqrt(CM4**2/CM2**4)-3
    sigma_kurt=np.sqrt(24/len(x))
    sigma2=2*np.sqrt(CM2)
    sigma3=3*np.sqrt(CM2)
    print("skew = {} +- {}".format(skew, sigma_skew))
    print("kurt = {} +- {}".format(kurt, sigma_kurt))
    print("2sigma = {} ".format(sigma2))
    print("3sigma = {} ".format(sigma3))
    return expvalue, sigma3
r_a1=rezidua(data1,W_a1) # do proměnné "r_a" uložíme hodnoty expvalue a sigma3<jupyter_output>skew = 0.4128300486622038 +- 0.14408763192842225
kurt = 2.259824194493061 +- 0.2881752638568445
2sigma = 8.399400566576297 
3sigma = 12.599100849864445<jupyter_text>Vidíme, že šikmost i špičatost se od nuly liší o méně než $2\sigma$, tudíž nulovou hypotézu, že rezidua mají normální rozdělení, nemůžeme zamítnout. Rezidua tedy odpovídají normálně rozděleným datům.  b) Rozptyl úměrný naměřeným hodnotám V případě b), kdy je rozptyl úměrný naměřeným hodnotám, bude po transformaci platit úměra $\sigma^2 \approx y^{-1}$. Na diagonále váhové matice se tedy budou nacházet hodnoty $y$.<jupyter_code>W_b=np.eye(len(x))*(y)
title_b="pripad b)"
b=analyza(data, W_b, title_b)
r_b=rezidua(data,W_b)<jupyter_output>skew = 2.1313234669615806 +- 0.1426148065624585
kurt = 15.326658180133009 +- 0.285229613124917
2sigma = 10.675018291482376 
3sigma = 16.012527437223564<jupyter_text>Opět vidíme, že rozdělení rezudií není normální, proto odděláme odlehlé hodnoty.<jupyter_code>pole_b1 = []
for z in pole:
   if (r_b[0]-r_b[1]) < (z[1] -(b[0]*np.exp(-z[0]*b[1]))) < r_b[0]+r_b[1]:
        pole_b1.append(z)
pole_b1 = np.array(pole_b1)
x2,y2=data2=np.transpose(pole_b1)
W_b1=np.eye(len(x2))*(y2)
title_b1="pripad b) bez odlehlých hodnot"
b1=analyza(data2, W_b1, title_b1)
r_b1=rezidua(data2,W_b1)<jupyter_output>skew = 0.304628577570488 +- 0.14334554477024897
kurt = 2.399426489469051 +- 0.28669108954049793
2sigma = 8.676805472612967 
3sigma = 13.01520820891945<jupyter_text>Nyní již jsou šikmost i špičatost měnší než $2\sigma$, tudíž jsme ověřili, že i v tomto případě mají rezidua normální rozdělení. Pro zajímavost můžeme nakonec vykreslit všechny tři případy, tj. výpočet bez váhové matice, případ a) a b) (oba bez odlehlých hodnot).<jupyter_code>a1=analyza(data1, W_a1, title_a1)
b1=analyza(data2, W_b1, title_b1)
plt.plot(x,A*np.exp(-x*L), label="bez W")
plt.legend(loc=0)
plt.rcParams['figure.figsize'] = 12, 8<jupyter_output>ln(A) = 4.509360443028548 +- 0.0011298384142490972
korel.koef. = 0.7156923771919719
Q_min= 49.14578134568504
A = 90.86368741869804 +- 0.10266128450596744
L = 0.4448481406870539 +- 0.0007175048418888641
ln(A) = 4.47948977477985 +- 0.007448733794281786
korel.koef. = 0.7426812855700488
Q_min= 44.31615345217824
A = 88.18966460346387 +- 0.6569013350381976
L = 0.45069001937841313 +- 0.002728624510593743
{'fraction_non_alphanumeric': 0.15443244714635695, 'fraction_numerical': 0.11332909899283186, 'mean_word_length': 5.7953144266337855, 'pattern_counts': {'":': 0, '<': 65, '<?xml version=': 0, '>': 61, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '460669', 'license': 'MIT', 'n_tokens_mistral': 5984, 'n_tokens_neox': 5328, 'n_words': 874, 'path': 'analysis/dockalova_lis1.ipynb', 'repo_name': 'limu007/galaxytea', 'size': '11021'}