WnioskowanieStatystyczne/ Regresja liniowa i test chi2: Różnice pomiędzy wersjami

Z Brain-wiki
 
(Nie pokazano 11 wersji utworzonych przez 3 użytkowników)
Linia 7: Linia 7:
 
{{Solution | title = Przykład: rozkłady ''Y''  dla każdego punktu ''X'' |text =  
 
{{Solution | title = Przykład: rozkłady ''Y''  dla każdego punktu ''X'' |text =  
 
<source lang= python>
 
<source lang= python>
 +
# -*- coding: utf-8 -*-
 
import scipy.stats as st
 
import scipy.stats as st
 
import pylab as py
 
import pylab as py
Linia 44: Linia 45:
 
</math>
 
</math>
 
:<math>
 
:<math>
  \hat b_0=\overline{y}-b\overline{x}
+
  \hat b_0=\overline{y}-\hat b_1\overline{x}
 
</math>
 
</math>
 
Z tymi współczynnikami otrzymujemy równanie opisujące prostą regresji:
 
Z tymi współczynnikami otrzymujemy równanie opisujące prostą regresji:
Linia 61: Linia 62:
 
Funkcję estymującą parametry i ich standardowe odchylenia można zaimplementować w pythonie następująco:
 
Funkcję estymującą parametry i ich standardowe odchylenia można zaimplementować w pythonie następująco:
 
<source lang= python>
 
<source lang= python>
 +
# -*- coding: utf-8 -*-
 
import scipy.stats as st
 
import scipy.stats as st
 
import pylab as py
 
import pylab as py
Linia 121: Linia 123:
 
(b0, b1, s_b0, s_b1, residua ) = regresja_liniowa(X,Y)
 
(b0, b1, s_b0, s_b1, residua ) = regresja_liniowa(X,Y)
  
print 'Równanie prostej: y = b0 + b1*x'
+
print('Równanie prostej: y = b0 + b1*x')
print 'dopasowane współczynniki: b0 = %.3f, b1 = %.3f' %(b0, b1)
+
print('dopasowane współczynniki: b0 = %.3f, b1 = %.3f' %(b0, b1))
print 's_b0 = %.4f, s_b1= %.4f '%(s_b0, s_b1)
+
print('s_b0 = %.4f, s_b1= %.4f '%(s_b0, s_b1))
 
 
py.errorbar(X,Y,sigma, fmt = None)
+
py.errorbar(X,Y,sigma)
 
Y_reg = b0 + b1*X  
 
Y_reg = b0 + b1*X  
 
py.plot(X,Y_reg)
 
py.plot(X,Y_reg)
Linia 154: Linia 156:
 
<source lang= python>
 
<source lang= python>
 
R2 = 1 - ss_err/ss_tot  
 
R2 = 1 - ss_err/ss_tot  
print 'R2 = %.2f' %(R2)
+
print('R2 = %.2f' %(R2))
 
</source>
 
</source>
  
 
W przypadku regresji liniowej <math>ss_{reg} + ss_{err} = ss_{tot}</math>. Możemy to sprawdzić [http://brain.fuw.edu.pl/edu/STAT:Regresja_liniowa#Interpretacja_wsp.C3.B3.C5.82czynnika_korelacji analitycznie] i numerycznie:
 
W przypadku regresji liniowej <math>ss_{reg} + ss_{err} = ss_{tot}</math>. Możemy to sprawdzić [http://brain.fuw.edu.pl/edu/STAT:Regresja_liniowa#Interpretacja_wsp.C3.B3.C5.82czynnika_korelacji analitycznie] i numerycznie:
 
<source lang= python>
 
<source lang= python>
print 'ss_tot = %.3f' %(ss_tot)
+
print('ss_tot = %.3f' %(ss_tot))
print 'ss_reg + ss_err =%.3f'%(ss_reg+ss_err)
+
print('ss_reg + ss_err =%.3f'%(ss_reg+ss_err))
 
</source>
 
</source>
 
czyli  
 
czyli  
Linia 194: Linia 196:
 
F = (ss_reg *(N-2))/ss_err
 
F = (ss_reg *(N-2))/ss_err
 
p_F = 1-st.f.cdf(F,1,N-2)
 
p_F = 1-st.f.cdf(F,1,N-2)
print 'F = %.2f, p_F = %.2f'%(F, p_F)
+
print('F = %.2f, p_F = %.2f'%(F, p_F))
 
</source>
 
</source>
 
''Wnioskowanie:'' Jeśli <tt>p_F</tt> jest duże to nie mamy powodu aby odrzucić hipotezę zerową. Jeśli zaś jest ono mniejsze niż ustalony poziom istotności to odrzucamy hipotezę zerową i przyjmujemy alternatywną.
 
''Wnioskowanie:'' Jeśli <tt>p_F</tt> jest duże to nie mamy powodu aby odrzucić hipotezę zerową. Jeśli zaś jest ono mniejsze niż ustalony poziom istotności to odrzucamy hipotezę zerową i przyjmujemy alternatywną.
Linia 221: Linia 223:
 
b1_l = b1 + s_b1*t_kryt
 
b1_l = b1 + s_b1*t_kryt
 
b1_h = b1 - s_b1*t_kryt
 
b1_h = b1 - s_b1*t_kryt
print '%.1f procentowe przedziały ufności parametrów:'%((1-alpha)*100)
+
print('%.1f procentowe przedziały ufności parametrów:'%((1-alpha)*100))
print 'b0: [%.2f %.2f ] '%(b0_l, b0_h)
+
print('b0: [%.2f %.2f ] '%(b0_l, b0_h))
print 'b1: [%.2f %.2f ] '%(b1_l, b1_h)
+
print('b1: [%.2f %.2f ] '%(b1_l, b1_h))
 
</source>
 
</source>
  
Linia 264: Linia 266:
 
<source lang= python>
 
<source lang= python>
 
W, p =st.shapiro(residua)
 
W, p =st.shapiro(residua)
print 'Test normalności residuów: p = %.3f'%(p)
+
print('Test normalności residuów: p = %.3f'%(p))
 
</source>
 
</source>
 
Jeśli tak to zmienna:
 
Jeśli tak to zmienna:
Linia 277: Linia 279:
 
else:
 
else:
 
p_chi2 = 1 - st.chi2.cdf(chi2, N-2)
 
p_chi2 = 1 - st.chi2.cdf(chi2, N-2)
print 'chi2 = %.2f, p_chi2 = %.3f' %(chi2, p_chi2)
+
print('chi2 = %.2f, p_chi2 = %.3f' %(chi2, p_chi2))
 
</source>
 
</source>
Czasem używamy zredukowanego, czyli podzielonego przez ilość stopni swobody <math>\chi^2</math>:
+
Czasem używamy zredukowanego, czyli podzielonego przez liczbę stopni swobody <math>\chi^2</math>:
 
*  Jeśli jest on znacząco większy niż 1 to model nie pasuje do danych, lub nie doszacowaliśmy standardowego odchylenia <math>\sigma</math>.  
 
*  Jeśli jest on znacząco większy niż 1 to model nie pasuje do danych, lub nie doszacowaliśmy standardowego odchylenia <math>\sigma</math>.  
 
*  Jeśli jest sporo mniejszy niż 1 to prawdopodobnie oszacowane przez nas <math>\sigma</math> jest większe niż rzeczywiste.
 
*  Jeśli jest sporo mniejszy niż 1 to prawdopodobnie oszacowane przez nas <math>\sigma</math> jest większe niż rzeczywiste.
 
To jakościowe porównanie można uściślić szacując prawdopodobieństwo zaobserwowania wartości <math>\chi^2_{zred}</math> bardziej ekstremalnych niż otrzymane w dopasowaniu. Zmienna <math>\chi^2_{zred}</math> podlega innemu rozkładowi prawdopodobieństwa niż <math>\chi^2</math>, możemy go jednak łatwo wyznaczyć w drodze symulacji:
 
To jakościowe porównanie można uściślić szacując prawdopodobieństwo zaobserwowania wartości <math>\chi^2_{zred}</math> bardziej ekstremalnych niż otrzymane w dopasowaniu. Zmienna <math>\chi^2_{zred}</math> podlega innemu rozkładowi prawdopodobieństwa niż <math>\chi^2</math>, możemy go jednak łatwo wyznaczyć w drodze symulacji:
<source lang= python><source lang= python>
+
<source lang= python>
 
chi2_zred = chi2/(N-2)
 
chi2_zred = chi2/(N-2)
 
# potrzebny jest nam rozkład chi2_zred:
 
# potrzebny jest nam rozkład chi2_zred:
Linia 292: Linia 294:
 
else:
 
else:
 
p_chi2_zred = np.sum(dist_chi2_zred<=chi2_zred)/float(N_dist)
 
p_chi2_zred = np.sum(dist_chi2_zred<=chi2_zred)/float(N_dist)
print 'chi2_zred = %.2f, p_chi2_zred = %.3f' %(chi2_zred, p_chi2_zred)
+
print('chi2_zred = %.2f, p_chi2_zred = %.3f' %(chi2_zred, p_chi2_zred))
 
</source>
 
</source>
  
 
=== Dopasowanie krzywej do danych gdy wariancje dla poszczególnych punktów pomiarowych są różne ===
 
=== Dopasowanie krzywej do danych gdy wariancje dla poszczególnych punktów pomiarowych są różne ===
Często w fizyce potrzebujemy dopasować jakąś bardziej skomplikowaną zależność niż prosta.  Często też potrafimy oszacować błędy pomiarowe dla różnych wartości zmiennej niezależnej, przy czym może się zdarzyć, że błędy te nie są jednakowe dla różnych wartości zmiennej niezależnej. Do dopasowania współczynników używamy zasady [[WnioskowanieStatystyczne/MLF|największej wiarygodności]], która prowadzi do procedur minimalizacji ważonego średniego błędu kwadratowego. Możemy wówczas użyć standardowych procedur minimalizacji gradientowej. Należy jednak pamiętać, że metody gradientowe znajdują najbliższe ''minimum lokalne'' analizowanej funkcji. W przypadku funkcji nieliniowych skutkiem tego  jest zależność wyniku od punktu startu minimalizacji. Poniżej rozważymy przykład dopasowania zależności wykładniczej.  
+
Często w fizyce potrzebujemy dopasować jakąś bardziej skomplikowaną zależność niż prosta.  Często też potrafimy oszacować błędy pomiarowe dla różnych wartości zmiennej niezależnej, przy czym może się zdarzyć, że błędy te nie są jednakowe dla różnych wartości zmiennej niezależnej. Do dopasowania współczynników używamy zasady [[WnioskowanieStatystyczne/MLF|największej wiarygodności]], która prowadzi do procedur minimalizacji ważonego średniego błędu kwadratowego. Możemy wówczas użyć standardowych procedur minimalizacji gradientowej. Należy jednak pamiętać, że metody gradientowe znajdują najbliższe ''minimum lokalne'' analizowanej funkcji. W przypadku funkcji nieliniowych skutkiem tego  jest zależność wyniku od punktu startu minimalizacji.
 +
 
 +
====Dopasowanie dowolnej funkcji====
 +
Poniżej rozważymy przykład dopasowania zależności wykładniczej.  
  
 
<source lang= python>
 
<source lang= python>
Linia 304: Linia 309:
 
import pylab as py
 
import pylab as py
 
import numpy as np
 
import numpy as np
 
+
 
# funkcja używana do symulowania danych
 
# funkcja używana do symulowania danych
 
def zanik(x, amp, wykladnik, blad_wzgledny):
 
def zanik(x, amp, wykladnik, blad_wzgledny):
Linia 314: Linia 319:
 
y += st.norm.rvs(size=num_points) * sigma       
 
y += st.norm.rvs(size=num_points) * sigma       
 
return (y, sigma)
 
return (y, sigma)
+
 
 
# Funkcja, którą chcemy dopasować do danych:
 
# Funkcja, którą chcemy dopasować do danych:
 
+
def funkcja_do_fitowania(x,a,b):
def funkcja_do_fitowania(b,x):
+
y = a*x**b
y = b[0]*x**b[1]
 
 
return y
 
return y
 
+
def funkcja_bledu( b, x, y, err):
+
def funkcja_bledu(x, y, funkcja, params, err):
'''Suma kwadratów tej funkcji jest minimalizowana w procesie optymalizacji parametrów'''
+
'''Suma kwadratów tej funkcji jest minimalizowana w procesie optymalizacji parametrów.
y_fit = funkcja_do_fitowania(b, x) # aktualne wartości y z dopasowania
+
Nam przyda się do obliczenia residuów.'''
residuum = y-y_fit  
+
y_fit = funkcja(x, *params) # aktualne wartości y z dopasowania
        # residua wchodzą do sumy kwadratów z wagą odwrotnie proporcjonalną do standardowego odchylenia
+
residuum = y-y_fit # residua wchodzą do sumy kwadratów z wagą odwrotnie proporcjonalną do standardowego odchylenia
 
residuum_wazone = residuum/ err  
 
residuum_wazone = residuum/ err  
 
return residuum_wazone
 
return residuum_wazone
 
  
 
# Generujemy punkty z szumem
 
# Generujemy punkty z szumem
Linia 335: Linia 337:
 
X = np.linspace(1.1, 10.1, num_points)
 
X = np.linspace(1.1, 10.1, num_points)
 
Y, sigma = zanik(X, 10.0, -2.0, 0.1)    # symulowane dane  
 
Y, sigma = zanik(X, 10.0, -2.0, 0.1)    # symulowane dane  
 
+
 
 
 
# Dopasowujemy parametry  
 
# Dopasowujemy parametry  
# musimy podać wartości startowe dla procedury minimalizacji
+
# nie musimy podawać wartości startowych (params_init) dla procedury minimalizacji (wtedy funkcja zakłada wartości startowe równe 1)
b_init = [2.0, -1.0]
+
# jednak zazwyczaj dobrze jest podpowiedzieć algorytmowi, gdzie powinien zacząć
 +
# nie musimy również podawać wartości sigma, ale jeśli są one różne dla różnych punktów, to podanie ich sprawi, że algorytm będzie się bardziej troszczył
 +
# o dopasowanie do punktów pomiarowych zmierzonych z dobrą dokładnością, a bardziej swobodnie podejdzie do tych o dużych niepewnościach
 +
params_init = [2.0, -1.0]
 +
params_final, covar = opt.curve_fit(funkcja_do_fitowania,X,Y,params_init,sigma)
  
# funkcja opt.leastsq tak zmienia parametry b
+
print("Dopasowane parametry",params_final)
# aby suma kwadratów wartości zwracanych przez funkcja_bledu była jak najmniejsza
+
print("Macierz kowariancji\n",covar)
  
out = opt.leastsq(funkcja_bledu, b_init, args=(X, Y, sigma), full_output=1)
+
# dopasowane parametry
 +
amp=params_final[0]
 +
wykladnik=params_final[1]
  
# odczytujemy wyniki minimalizacji
 
b_final = out[0]
 
# Funkcja lestsq zwraca nie macierz kowariancji tylko według dokumentacji:
 
# cov_x is a Jacobian approximation to the Hessian of the least squares
 
# objective function. This matrix must be multiplied by the residual
 
# variance to get the covariance of the parameter estimates - see curve_fit
 
# zatem obliczamy wariancję resztkową:
 
s_sq = (funkcja_do_fitowania(b_final,X)**2).sum()/(len(Y)-len(b_init))
 
# i obliczamy kowariancję:
 
covar = out[1]*s_sq
 
print  b_final
 
print covar
 
# dopasowane parametry
 
amp=b_final[0]
 
wykladnik=b_final[1]
 
 
# standardowe błędy dopasowania
 
# standardowe błędy dopasowania
amp_err = np.sqrt( covar[0][0] )  
+
amp_err = np.sqrt(covar[0][0])
wykladnik_err = np.sqrt(covar[1][1] )
+
wykladnik_err = np.sqrt(covar[1][1])
 
+
 
# test chi2 dobroci dopasowania.  
 
# test chi2 dobroci dopasowania.  
 
# Jeśli znamy wariancję błędu pomiarowego można zastosować test chi2 do oceny jakości dopasowania.
 
# Jeśli znamy wariancję błędu pomiarowego można zastosować test chi2 do oceny jakości dopasowania.
 
# Po pierwsze powinniśmy przetestować czy residua mają rozkład normalny
 
# Po pierwsze powinniśmy przetestować czy residua mają rozkład normalny
 
+
residua = funkcja_bledu( b_final, X, Y, sigma)# tym razem residua już są podzielone przez standardowe odchylenie, każde przez swoje
+
residua = funkcja_bledu(X, Y, funkcja_do_fitowania, params_final, sigma)# tym razem residua już są podzielone przez standardowe odchylenie, każde przez swoje
 
W, p =st.shapiro(residua)
 
W, p =st.shapiro(residua)
print 'Test normalności residuów: p = %.3f'%(p)
+
print('Test normalności residuów: p = %.3f'%(p))
 
# jeśli tak to zmienna:
 
# jeśli tak to zmienna:
 
chi2 = np.sum(residua**2)
 
chi2 = np.sum(residua**2)
 
# podlega rozkładowi chi-kwadrat o N - n  ilości stopni swobody (n - ilość fitowanych parametrów), czyli u nas N-2
 
# podlega rozkładowi chi-kwadrat o N - n  ilości stopni swobody (n - ilość fitowanych parametrów), czyli u nas N-2
 
# możemy zbadać jakie jest prawdopodobieństwo zaobserwowania takiej, bądź bardziej ekstremalnej wartości chi2:
 
# możemy zbadać jakie jest prawdopodobieństwo zaobserwowania takiej, bądź bardziej ekstremalnej wartości chi2:
#
+
 
 
N = len (X)
 
N = len (X)
if chi2 < N-2:
+
liczba_stopni_swobody = N-len(params_final) # liczba punktów - liczba parametrów
p_chi2 = st.chi2.cdf(chi2, N-2)
+
if chi2 < liczba_stopni_swobody:
 +
p_chi2 = st.chi2.cdf(chi2, liczba_stopni_swobody)
 
else:
 
else:
p_chi2 = 1 - st.chi2.cdf(chi2, N-2)
+
p_chi2 = st.chi2.sf(chi2, liczba_stopni_swobody) # równoważne 1-st.chi2.cdf(chi2, N-2), ale sf ma lepszą dokładność dla małych wartości
print 'chi2 = %.2f, p_chi2 = %.3f' %(chi2, p_chi2)
+
print('chi2 = %.2f, p_chi2 = %.3f' %(chi2, p_chi2))
 
# czasem używamy zredukowanego chi2, czyli podzielonego przez ilość stopni swobody
 
# czasem używamy zredukowanego chi2, czyli podzielonego przez ilość stopni swobody
chi2_zred = chi2/(N-2)
+
chi2_zred = chi2/liczba_stopni_swobody
 
# jeśli jest on znacząco większy niż 1 to model nie pasuje do danych, lub nie doszacowaliśmy sigmy,  
 
# jeśli jest on znacząco większy niż 1 to model nie pasuje do danych, lub nie doszacowaliśmy sigmy,  
 
# jeśli jest sporo mniejszy niż 1 to prawdopodobnie oszacowane przez nas sigma jest większe niż rzeczywiste
 
# jeśli jest sporo mniejszy niż 1 to prawdopodobnie oszacowane przez nas sigma jest większe niż rzeczywiste
 
# potrzebny jest nam rozkład chi2_zred:
 
# potrzebny jest nam rozkład chi2_zred:
 
N_dist = 100000
 
N_dist = 100000
dist_chi2_zred = np.sum(st.norm.rvs(size=(N-2,N_dist))**2  ,0)/(N-2)
+
dist_chi2_zred = np.sum(st.norm.rvs(size=(liczba_stopni_swobody,N_dist))**2  ,0)/liczba_stopni_swobody
 
p_chi2_zred = np.sum(dist_chi2_zred>=chi2_zred)/float(N_dist)
 
p_chi2_zred = np.sum(dist_chi2_zred>=chi2_zred)/float(N_dist)
print 'chi2_zred = %.2f, p_chi2_zred = %.3f' %(chi2_zred, p_chi2_zred)
+
print('chi2_zred = %.2f, p_chi2_zred = %.3f' %(chi2_zred, p_chi2_zred))
+
 
+
 
 
 
##########
 
##########
# rysunek
+
# wykres
 
##########
 
##########
 
+
py.subplot(2,1,1)
py.clf()
+
py.plot(X, funkcja_do_fitowania(X,amp,wykladnik))    # Fit
py.plot(X, funkcja_do_fitowania(b_final,X))    # Fit
 
 
py.errorbar(X, Y, yerr=sigma, fmt='k.')  # Dane i błędy
 
py.errorbar(X, Y, yerr=sigma, fmt='k.')  # Dane i błędy
 
py.text(5, 6.5, 'amplituda = %5.2f +/- %5.2f' % (amp, amp_err))
 
py.text(5, 6.5, 'amplituda = %5.2f +/- %5.2f' % (amp, amp_err))
Linia 408: Linia 399:
 
py.ylabel('Y')
 
py.ylabel('Y')
 
py.xlim(1, 11)
 
py.xlim(1, 11)
 +
py.subplot(2,1,2)
 +
py.plot(X, residua) # residua
 +
py.xlabel('X')
 +
py.ylabel('dY')
 +
py.title(u'Wykres residuów')
 
py.show()
 
py.show()
 
</source>
 
</source>
  
==== Przykład z pracowni ====
+
====Dopasowanie wielomianu====
[[TI:Programowanie_z_Pythonem/Zadania_dodatkowe#Wczytywanie_danych_i_wykres| Przykład w odnośniku ]] przedstawia dopasowanie parametrów krzywej opisującej dochodzenie ciała do równowagi termicznej z otoczeniem.
+
Poniżej rozważymy przykład dopasowania zależności wielomianowej.
 +
 
 +
<source lang=python>
 +
# -*- coding: utf-8 -*-
 +
import scipy.stats as st
 +
import pylab as py
 +
import numpy as np
 +
 +
# funkcja używana do symulowania danych
 +
def wielomian_z_szumem(x, wspolczynniki,blad_wzgledny):
 +
'''Definicja funkcji wielomianowej. Użyjemy jej do wytworzenia danych'''
 +
W = np.poly1d(wspolczynniki) # funkcja zwracająca obiekt wielomianu o zadanych wspolczynnikach
 +
#można go używać tak, jak zwykłej funkcji, ale obsługuje też działania na wielomianach
 +
y = W(X)# idealne dane
 +
sigma = blad_wzgledny * y # zakładamy, że stały jest błąd względny pomiaru
 +
# przeliczamy go na standardowe odchylenie symulowanego  błędu
 +
# symulujemy szum z obliczonym odchyleniem standardowym i dodajemy go do danych idealnych
 +
y += st.norm.rvs(size=num_points) * sigma     
 +
return (y, sigma)
 +
 
 +
def funkcja_bledu_dla_wielomianow(x, y, wspolczynniki, err):
 +
'''Suma kwadratów tej funkcji jest minimalizowana w procesie optymalizacji parametrów.
 +
Nam przyda się do obliczenia residuów.'''
 +
W = np.poly1d(wspolczynniki)
 +
y_fit = W(x) # aktualne wartości y z dopasowania
 +
residuum = y-y_fit # residua wchodzą do sumy kwadratów z wagą odwrotnie proporcjonalną do standardowego odchylenia
 +
residuum_wazone = residuum/ err
 +
return residuum_wazone
 +
 
 +
# Generujemy punkty z szumem
 +
num_points = 20
 +
X = np.linspace(-4, 6, num_points)
 +
wspolczynniki_wielomianu= (0.3,1,-2,4)
 +
stopien_wielomianu=len(wspolczynniki_wielomianu)-1
 +
blad_wzgledny_pomiaru=0.1
 +
Y, sigma = wielomian_z_szumem(X, wspolczynniki_wielomianu, blad_wzgledny_pomiaru)    # symulowane dane
 +
 +
# Dopasowujemy parametry
 +
# tym razem skorzystamy z funkcji np.polyfit, która nie potrzebuje parametrów początkowych, ani zdefiniowanej funkcji, którą ma dopasować
 +
# podajemy jej tylko nasze dane oraz stopień wielomianu, który ma dopasować oraz opcjonalne wagi
 +
# UWAGA! Tym razem wagi muszą być odwrotnością odchyleń standardowych (1/sigma, a nie sigma, jak w curve_fit)
 +
# funkcja ta domyślnie zwraca tylko dopasowane parametry (wspolczynniki wielomianu), a nie zwraca macierzy kowariancji,
 +
# jeśli jest nam ona potrzebna, to musimy jej zarządać poprzez dodanie opcji cov=True (full=False, ale to jest domyślnie)
 +
params_final, covar=np.polyfit(X, Y, deg=stopien_wielomianu, w=1/sigma, cov=True)
 +
 
 +
print("Dopasowane wspolczynniki wielomianu",params_final)
 +
print("Macierz kowariancji\n",covar)
 +
 
 +
# standardowe błędy dopasowania
 +
niepewnosci=[]
 +
for i in range(len(params_final)):
 +
niepewnosci.append(np.sqrt(covar[i][i]))
 +
print(niepewnosci)
 +
 
 +
 
 +
# test chi2 dobroci dopasowania.
 +
# Jeśli znamy wariancję błędu pomiarowego można zastosować test chi2 do oceny jakości dopasowania.
 +
# Po pierwsze powinniśmy przetestować czy residua mają rozkład normalny
 +
 +
residua = funkcja_bledu_dla_wielomianow(X, Y, params_final, sigma)# tym razem residua już są podzielone przez standardowe odchylenie, każde przez swoje
 +
W, p =st.shapiro(residua)
 +
print('Test normalności residuów: p = %.3f'%(p))
 +
# jeśli tak to zmienna:
 +
chi2 = np.sum(residua**2)
 +
# podlega rozkładowi chi-kwadrat o N - n  ilości stopni swobody (n - ilość fitowanych parametrów), czyli u nas N-2
 +
# możemy zbadać jakie jest prawdopodobieństwo zaobserwowania takiej, bądź bardziej ekstremalnej wartości chi2:
 +
 
 +
 
 +
N = len (X)
 +
liczba_stopni_swobody = N-len(params_final) # liczba punktów - liczba parametrów
 +
if chi2 < liczba_stopni_swobody:
 +
p_chi2 = st.chi2.cdf(chi2, liczba_stopni_swobody)
 +
else:
 +
p_chi2 = st.chi2.sf(chi2, liczba_stopni_swobody) # równoważne 1-st.chi2.cdf(chi2, N-2), ale sf ma lepszą dokładność dla małych wartości
 +
print('chi2 = %.2f, p_chi2 = %.3f' %(chi2, p_chi2))
 +
 
 +
 
 +
 
 +
# czasem używamy zredukowanego chi2, czyli podzielonego przez ilość stopni swobody
 +
chi2_zred = chi2/liczba_stopni_swobody
 +
# jeśli jest on znacząco większy niż 1 to model nie pasuje do danych, lub nie doszacowaliśmy sigmy,
 +
# jeśli jest sporo mniejszy niż 1 to prawdopodobnie oszacowane przez nas sigma jest większe niż rzeczywiste
 +
# potrzebny jest nam rozkład chi2_zred:
 +
N_dist = 100000
 +
dist_chi2_zred = np.sum(st.norm.rvs(size=(liczba_stopni_swobody,N_dist))**2  ,0)/liczba_stopni_swobody
 +
p_chi2_zred = np.sum(dist_chi2_zred>=chi2_zred)/float(N_dist)
 +
print('chi2_zred = %.2f, p_chi2_zred = %.3f' %(chi2_zred, p_chi2_zred))
 +
 +
 +
##########
 +
# wykres
 +
##########
 +
py.subplot(2,1,1)
 +
W=np.poly1d(params_final)
 +
py.plot(X, W(X))    # Fit
 +
py.errorbar(X, Y, yerr=sigma, fmt='k.')  # Dane i błędy
 +
py.title(u'Dopasowanie metodą najmniejszych kwadratów')
 +
py.text(-4.6, 92, u'dopasowane współczynniki = '+str(np.round(params_final,3)))
 +
py.text(-4.6, 86, u'niepewności współczynników = '+str(np.round(niepewnosci,3)))
 +
py.text(-4.6, 80, u'prawdziwe współczynniki = '+str(np.round(wspolczynniki_wielomianu,3)))
 +
py.xlabel('X')
 +
py.ylabel('Y')
 +
py.xlim(X.min()-1, X.max()+1)
 +
py.subplot(2,1,2)
 +
py.plot(X, residua) # residua
 +
py.xlabel('X')
 +
py.ylabel('dY')
 +
py.title(u'Wykres residuów')
 +
py.show()
 +
py.show()</source>

Aktualna wersja na dzień 14:57, 5 cze 2019

Wstęp

Załóżmy, że mamy dwie zmienne losowe ciągłe [math]X[/math] i [math]Y[/math]. Chcielibyśmy wykorzystać wiedzę o wartościach zmiennej [math]X[/math] do przewidywania wartości zmiennej [math]Y[/math]. Mówimy, że zmienna [math]X [/math] jest niezależna, a zmienna [math]Y[/math] zależna. W fizyce taką wiedzę opisujemy przy pomocy równań. Równania fizyczne często wyrażają związki przyczynowo-skutkowe. W takim wypadku, która zmienna jest zależna, a która niezależna ma głębszy sens. Jednak nie zawsze tak musi być. Wartości dwóch zmiennych mogą zależeć od trzeciej nieobserwowanej zmiennej. W tej sytuacji wiedza o wartości jednej z tych zmiennych może być wykorzystana do przewidywania wartości drugiej, ale nie ma między nimi związku przyczynowo-skutkowego.

Regresja

W ogólności, dla każdej wartości zmiennej [math]X[/math] mamy rozkład wartości zmiennej [math]Y[/math].


Regresja liniowa

Dalej będziemy rozważać regresję liniową, tzn. założymy, że punkty [math](X,Y)[/math] są generowane przez model liniowy o następującym równaniu:

[math] y = b_0 + b_1x + \epsilon [/math]

współczynniki [math]b_0[/math] i [math]b_1[/math] można wyestymować stosując metodę największej wiarygodności:

[math]\hat b_1=\frac{\underset{i=1}{\overset{N}{\sum }}(x_{i}-\overline{x})(y_{i}- \overline{y})}{\underset{i=1}{\overset{N}{\sum }}(x_{i}-\overline{x})^{2}}, [/math]
[math] \hat b_0=\overline{y}-\hat b_1\overline{x} [/math]

Z tymi współczynnikami otrzymujemy równanie opisujące prostą regresji: [math]\hat y = \hat b_0 + \hat b_1 x[/math]

Zakłądając, że [math]\epsilon [/math] pochodzi z rozkładu normalnego o wariancji [math]\sigma^2[/math] estymowane współczynniki są zmiennymi losowymi pochodzącymi z rozkładów normalnego o średniej takiej jak wyestymowany współczynnik i wariancji odpowiednio:

[math] v_{b_1} = \frac{\sigma^2}{\sum_i (x_i - \bar x)^2} [/math]
[math] v_{b_0} = \frac{\sigma^2}{n} + \frac{\bar x ^2}{\sum_i (x_i - \bar x)^2} \sigma^2 [/math]

Wariancję [math]\sigma^2 = E[S^2][/math] można estymować przez:

[math]S^2 = \frac{1}{n-2}\sum_{i=1}^n(y_i - \hat y_i)^2[/math]

Warto tu zwrócić uwagę na prosty fakt, że niepewność oszacowania współczynników można zmniejszyć zwiększając zakres zmiennej [math]x[/math].

Funkcję estymującą parametry i ich standardowe odchylenia można zaimplementować w pythonie następująco:

# -*- coding: utf-8 -*-
import scipy.stats as st
import pylab as py
import numpy as np

def regresja_liniowa(X,Y):
	'''równanie dopasowywanej prostej to y = b0 + b1*x
	argumenty:
	  X - zmienna niezależna
	  Y - zmienna zależna
	funkcja zwraca:
	  b0, b1, - estymaty parametrów
	  s_b0, s_b1, - estymaty standardowego odchylenia parametrów 
	  residua - różnice między punktami pomiarowymi a punktami na dopasowanej prostej
	'''
	N = len(X)
	x_sr = np.mean(X)
	y_sr = np.mean(Y)
	# estymatory parametrów 
	# korzystamy z tego że numpy wykonuje odejmowania i potęgowania dla każdego elementu tablicy X i Y
	b1 = np.sum((X-x_sr)*(Y-y_sr))/np.sum((X-x_sr)**2)
	b0 = y_sr - b1*x_sr
	
	# teraz liczymy kilka rzeczy przydatnych do oceny jakości modelu
	Y_reg = b0 + b1*X # wartości Y przewidywane przez model
	residua = Y - Y_reg # residua, czyli zmienność Y nie wynikająca z modelu

	
	sse = np.sum(residua**2)
	# estymator wariancji residuów, bywa nazywany średnim błędem kwadratowym regresji :
	v_e = sse/(N-2)
	# estymatory standardowych błędów parametrów
	s_b0 = np.sqrt(v_e) * np.sqrt(1.0/N + x_sr**2/np.sum( (X-x_sr)**2))
	s_b1 = np.sqrt(v_e) * np.sqrt( 1.0/np.sum( (X -x_sr)**2 ))
	return (b0, b1, s_b0, s_b1, residua )

Przykład: Dopasowanie prostej do punktów (zakładamy jednakową wariancję Y dla każdego X)

Wytwórzmy dane zgodnie z modelem:

[math]y = -13 + 3 x + \epsilon[/math]
i [math]\sigma_\epsilon =19[/math]:
# symulowana zależność ma następującą postać y = b0 + b1*x 
# wartości parametrów 
b0 = -13.0
b1 = 3.0


X = np.arange(30, 70, 0.5)
sigma = 19.0
n = 1
Y = np.zeros(len(X))
for i in range(len(X)):
	Y[i] =  b0 + b1*X[i] + st.norm.rvs(size = n, loc=0, scale = sigma)

Korzystając ze zdefiniowanej powyżej funkcji regresja_liniowa estymujemy parametry i ich odchylenia standardowe:

(b0, b1, s_b0, s_b1, residua ) = regresja_liniowa(X,Y)

print('Równanie prostej: y = b0 + b1*x')
print('dopasowane współczynniki: b0 = %.3f, b1 = %.3f' %(b0, b1))
print('s_b0 = %.4f, s_b1= %.4f '%(s_b0, s_b1))
	
py.errorbar(X,Y,sigma)
Y_reg = b0 + b1*X 
py.plot(X,Y_reg)

py.show()

Ocena jakości dopasownia

Współczynnik [math]R^2[/math]

Aby wyrazić współczynnik [math]R^2[/math] potrzebujemy następujących wyrażeń - sum kwadratów (ss). Są one miarą zmienności.

  • [math]ss_{tot}[/math] - całkowita suma kwadratów - proporcjonalna do wariancji próby,
  • [math]ss_{reg}[/math] - suma kwadratów regresji - zwana też wyjaśnioną sumą kwadratów,
  • [math]ss_{err}[/math] - suma kwadratów residuów - niewyjaśniona suma kwadratów.

Poszczególne składniki wymienionych powyżej sum kwadratów są zilustrowane na poniższym rysunku.

Plik:Regresja1.svg
Dla wybranego punktu [math](x_i,y_i)[/math] zaznaczono różnice będące składnikami poszczególnych sum kwadratów

Implementacja:

y_sr = np.mean(Y)
ss_tot = np.sum( (Y -  y_sr)**2 ) 
ss_reg = np.sum( (Y_reg - y_sr)**2 )
ss_err = np.sum( (residua)**2 )

mając te sumy [math]R^2[/math] definiujemy jako:

[math]R^2 = 1 - \frac{ss_{err}}{ss_{tot}}[/math]
R2 = 1 - ss_err/ss_tot 
print('R2 = %.2f' %(R2))

W przypadku regresji liniowej [math]ss_{reg} + ss_{err} = ss_{tot}[/math]. Możemy to sprawdzić analitycznie i numerycznie:

print('ss_tot = %.3f' %(ss_tot))
print('ss_reg + ss_err =%.3f'%(ss_reg+ss_err))

czyli

[math]R^2 = \frac{ss_{reg}}{ss_{tot}} [/math],

można więc interpretować [math]R^2[/math] jako frakcję zmienności Y tłumaczoną przez model. W przypadku regresji liniowej współczynnik [math]R^2[/math] równy jest kwadratowi współczynnika korelacji [math]\rho[/math]

[math]R^2 = \rho^2[/math]

(dowód)

Test F dla hipotezy o braku korelacji

Często interesujące jest zweryfikowanie hipotezy o istotności zależności między Y a X (proszę nie mylić tego z istnieniem związku przyczynowo-skutkowego). Matematycznie równoważne jest to postawieniu hipotezy:

[math]H_0: b_1 = 0 [/math]
[math]H_1: b_1 \ne 0 [/math]

albo:

[math]H_0: \rho = 0 [/math]
[math]H_1: \rho \ne 0 [/math]

Wykorzystamy do tego test równości wariancji oparty o rozkład F. Jeśli zgodnie z [math]H_0[/math] [math]b_1 = 0[/math] to prosta regresji jest pozioma i wariancja wyjaśniona przez regresję (proporcjonalna do [math]ss_{reg}[/math]) jest równa wariancji niewyjaśnionej (proporcjonalna do [math]ss_{err}[/math]). Wariancje te można estymować dzieląc odpowiednie sumy kwadratów zdefiniowane w poprzednim paragrafie przez odpowiadającą im liczbę stopni swobody. Jeśli mamy N punktów danych, to:

  • liczba stopni swobody dla [math]ss_{tot}[/math] jest [math]N-1[/math], poniważ jeden stopień swobody jest tracony na obliczenie średniej,
  • liczba stopni swobody dla [math]ss_{err}[/math] jest [math]N-2[/math], ponieważ do policzenia tej sumy kwadratów musimy wyznaczyć dwa parametry prostej,
  • liczba stopni swobody odpowiadająca [math]ss_{reg}[/math] jest 1, bo jest [math]ss_{reg}[/math] związana jest z poprzednimi sumami kwadratów równaniem, czyli swobody jest tyle ile wynosi różnica w stopniach swobody tamtych sum.

Zatem:

  • estymator wariancji wyjaśnionej:
[math]s_{reg} = \frac{ss_{reg}}{1}[/math]
  • estymator wariancji niewyjaśnionej:
[math]s_{err} = \frac{ss_{err}}{N-2}[/math]

Wielkość

[math] F = \frac{ss_{reg}(N-2)}{ss_{err}} [/math] podlega rozkładowi F o [math](1,N-2)[/math] stopniach swobody.

W naszym przykładzie:

# test F
N = len(X)
F = (ss_reg *(N-2))/ss_err
p_F = 1-st.f.cdf(F,1,N-2)
print('F = %.2f, p_F = %.2f'%(F, p_F))

Wnioskowanie: Jeśli p_F jest duże to nie mamy powodu aby odrzucić hipotezę zerową. Jeśli zaś jest ono mniejsze niż ustalony poziom istotności to odrzucamy hipotezę zerową i przyjmujemy alternatywną.

Przedziały ufności dla parametrów

Przedziały ufności dla parametrów [math]b_0[/math] i [math]b_1[/math] pokazują zakres, w jakim z zadanym prawdopodobieństwem znajdują się ich "prawdziwe" wartości.

Jeśli residua mają rozkład normalny, to estymatory parametrów [math]b_0[/math] i [math]b_1[/math] również będą miały rozkład normalny. Zmienne:

[math]t = \frac{\hat b_0 - b_0}{s_{\hat b_0}}\ \sim\ t_{N-2},[/math]
[math]t = \frac{\hat b_1 - b_1}{s_{\hat b_1}}\ \sim\ t_{N-2},[/math]

podlegają rozkładowi t z (N−2) stopniami swobody.

Używając powyższych statystyk t można skonstruować przedziały ufności w standardowy sposób (porównaj z przykładem). Jeśli przedział ma mieć poziom ufności [math]1 - \alpha[/math] to potrzebna nam będzie wartość krytyczna z rozkładu [math]t^*_{N-2}[/math] taka, że prawdopodobieństwo zaobserwowania wartości t nie większej od niej jest [math]\alpha/2[/math]. Wówczas:

[math] b_1\in \Big[\ \hat b_1 - s_{\hat b_1} t^*_{N-2},\ \hat b_1 + s_{\hat b_1} t^*_{N-2}\ \Big] [/math]

oraz

[math] b_0 \in \Big[\ \hat b_0 - s_{\hat b_0} t^*_{N-2},\ \hat b_0 + s_{\hat b_0} t^*_{N-2}\ \Big] [/math]

Implementacja:

# przedziały ufności:
alpha = 0.05  # zakładam 95% przedział ufności
# wartość krytyczna w rozkładzie t
t_kryt = st.t.ppf(alpha/2, N-2)
b0_l = b0 + s_b0*t_kryt
b0_h = b0 - s_b0*t_kryt
b1_l = b1 + s_b1*t_kryt
b1_h = b1 - s_b1*t_kryt
print('%.1f procentowe przedziały ufności parametrów:'%((1-alpha)*100))
print('b0: [%.2f %.2f ] '%(b0_l, b0_h))
print('b1: [%.2f %.2f ] '%(b1_l, b1_h))
Przedziały ufności dla modelu

Widzieliśmy, że parametry dopasowanej prostej nie są wyznaczone dokładnie. Tzn. jeśli dostalibyśmy inne realizacje danych (X,Y) to ta sama procedura regresji zwraca nieco inne parametry modelu. Jak widzieliśmy powyżej można wyznaczyć przedziały ufności wewnątrz których parametry te znajdują się z określonym prawdopodobieństwem. Różnym parametrom odpowiadają różne proste. Proste te wyznaczają na płaszczyźnie (x,y) pewien obszar. Obszar ten to przedział ufności dla modelu. Jego granice można wyznaczyć obliczając dla każdej wartości x błąd standardowy regresji ze wzoru:

[math]s_{reg}(x_i) = \sqrt{\frac{ss_{err}}{N-2}} \cdot \sqrt{\frac{1}{N} + \frac{(x_i - \bar X)^2}{\sum_{j=1}^N(x_j - \bar X)^2}} [/math]

odległość krzywej wyznaczającej obszar ufności od prostej regresji znajdujemy mnożąc ten błąd standardowy przez odpowiednią wartość krytyczną z rozkładu [math]t_{N-2}[/math]:

[math]d_i = t^*_{N-2}s_{reg}(x_i) [/math]

Implementacja:

# Przedział ufności modelu:
alpha = 0.05  # zakładam 95% przedział ufności
# wartość krytyczna w rozkładzie t
t_kryt = st.t.ppf(alpha/2, N-2)
sse = np.sum(residua**2)
# estymator wariancji residuów, bywa nazywany średnim błędem kwadratowym regresji :
v_e = sse/(N-2)
x_sr = np.mean(X)
# Odległość brzegów przedziału ufności od prostej regresji
d = t_kryt*np.sqrt(v_e)*np.sqrt(1.0/N + (X- x_sr)**2/np.sum((X-x_sr)**2))
# Ilustracja: dla każdego X cieniujemy obszar pomiędzy Y_reg-d,Y_reg+d i nadajemy mu przezroczystość 0.5
py.fill_between(X,Y_reg-d,Y_reg+d,alpha=0.5)
Przedziały ufności dla obserwacji

Przedział zmienności dla modelu nie mówi nam wiele o tym jak daleko od wyznaczonej prostej mogą pojawiać się nowe obserwacje (x,y). Aby zobrazować obszar, w którym z określonym prawdopodobieństwem mogą wystąpić nowe obserwacje potrzebujemy przedziału ufności dla obserwacji. Jego granice można wyznaczyć obliczając dla każdej wartości x błąd standardowy ze wzoru:

[math]s_{reg}(x_i) = \sqrt{\frac{ss_{err}}{N-2}} \cdot \sqrt{1+\frac{1}{N} + \frac{(x_i - \bar X)^2}{\sum_{j=1}^N(x_j - \bar X)^2}} [/math]

odległość krzywej wyznaczającej obszar ufności od prostej regresji znajdujemy mnożąc ten błąd standardowy przez odpowiednią wartość krytyczną z rozkładu [math]t_{N-2}[/math]:

[math]d_i = t^*_{N-2}s_{reg}(x_i) [/math]
# przedział ufności na obserwacje
d = t_kryt*np.sqrt(v_e)*np.sqrt(1+1.0/N + (X- x_sr)**2/np.sum((X-x_sr)**2))
py.fill_between(X,Y_reg-d,Y_reg+d, facecolor='gray',alpha=0.5)
Test [math]\chi^2[/math]

Jeśli znamy wariancję błędu pomiarowego można zastosować test [math]\chi^2[/math] do oceny jakości dopasowania. Po pierwsze powinniśmy przetestować czy residua mają rozkład normalny

W, p =st.shapiro(residua)
print('Test normalności residuów: p = %.3f'%(p))

Jeśli tak to zmienna:

[math]\chi_{fit}^2 = \sum_{i=1}^N {\left( \frac{y_i-y_{reg}}{\sigma} \right)^2 }[/math]

podlega rozkładowi [math]\chi^2[/math] o [math]N - n[/math] ilości stopni swobody (n - ilość estymowanych parametrów), czyli u nas N-2. Możemy zbadać jakie jest prawdopodobieństwo zaobserwowania takiej ([math]\chi_{fit}^2[/math]), bądź bardziej ekstremalnej wartości [math]\chi^2[/math]:

chi2 = np.sum(residua**2)/sigma**2
N = len(X)
if chi2 < N-2:
	p_chi2 = st.chi2.cdf(chi2, N-2)
else:
	p_chi2 = 1 - st.chi2.cdf(chi2, N-2)
print('chi2 = %.2f, p_chi2 = %.3f' %(chi2, p_chi2))

Czasem używamy zredukowanego, czyli podzielonego przez liczbę stopni swobody [math]\chi^2[/math]:

  • Jeśli jest on znacząco większy niż 1 to model nie pasuje do danych, lub nie doszacowaliśmy standardowego odchylenia [math]\sigma[/math].
  • Jeśli jest sporo mniejszy niż 1 to prawdopodobnie oszacowane przez nas [math]\sigma[/math] jest większe niż rzeczywiste.

To jakościowe porównanie można uściślić szacując prawdopodobieństwo zaobserwowania wartości [math]\chi^2_{zred}[/math] bardziej ekstremalnych niż otrzymane w dopasowaniu. Zmienna [math]\chi^2_{zred}[/math] podlega innemu rozkładowi prawdopodobieństwa niż [math]\chi^2[/math], możemy go jednak łatwo wyznaczyć w drodze symulacji:

chi2_zred = chi2/(N-2)
# potrzebny jest nam rozkład chi2_zred:
N_dist = 100000
dist_chi2_zred = np.sum(st.norm.rvs(size=(N-2,N_dist))**2  ,0)/(N-2)
if chi2_zred>1:
	p_chi2_zred = np.sum(dist_chi2_zred>=chi2_zred)/float(N_dist)
else:
	p_chi2_zred = np.sum(dist_chi2_zred<=chi2_zred)/float(N_dist)
print('chi2_zred = %.2f, p_chi2_zred = %.3f' %(chi2_zred, p_chi2_zred))

Dopasowanie krzywej do danych gdy wariancje dla poszczególnych punktów pomiarowych są różne

Często w fizyce potrzebujemy dopasować jakąś bardziej skomplikowaną zależność niż prosta. Często też potrafimy oszacować błędy pomiarowe dla różnych wartości zmiennej niezależnej, przy czym może się zdarzyć, że błędy te nie są jednakowe dla różnych wartości zmiennej niezależnej. Do dopasowania współczynników używamy zasady największej wiarygodności, która prowadzi do procedur minimalizacji ważonego średniego błędu kwadratowego. Możemy wówczas użyć standardowych procedur minimalizacji gradientowej. Należy jednak pamiętać, że metody gradientowe znajdują najbliższe minimum lokalne analizowanej funkcji. W przypadku funkcji nieliniowych skutkiem tego jest zależność wyniku od punktu startu minimalizacji.

Dopasowanie dowolnej funkcji

Poniżej rozważymy przykład dopasowania zależności wykładniczej.

# -*- coding: utf-8 -*-
import scipy.stats as st
import scipy.optimize as opt
import pylab as py
import numpy as np
 
# funkcja używana do symulowania danych
def zanik(x, amp, wykladnik, blad_wzgledny):
	'''Definicja funkcji zaniku wykładniczego. Użyjemy jej do wytworzenia danych'''
	y = amp * (x**wykladnik)       # idealne dane
	sigma = blad_wzgledny * y      # zakładamy, że stały jest błąd względny pomiaru
				       # przeliczamy go na standardowe odchylenie symulowanego  błędu 
	# symulujemy szum z obliczonym odchyleniem standardowym i dodajemy go do danych idealnych										
	y += st.norm.rvs(size=num_points) * sigma       
	return (y, sigma)
 
# Funkcja, którą chcemy dopasować do danych:
def funkcja_do_fitowania(x,a,b):
	y = a*x**b
	return y
 
def funkcja_bledu(x, y, funkcja, params, err):
	'''Suma kwadratów tej funkcji jest minimalizowana w procesie optymalizacji parametrów.
	Nam przyda się do obliczenia residuów.'''
	y_fit = funkcja(x, *params) # aktualne wartości y z dopasowania
	residuum = y-y_fit # residua wchodzą do sumy kwadratów z wagą odwrotnie proporcjonalną do standardowego odchylenia
	residuum_wazone = residuum/ err 
	return residuum_wazone

# Generujemy punkty z szumem
num_points = 20
X = np.linspace(1.1, 10.1, num_points)
Y, sigma = zanik(X, 10.0, -2.0, 0.1)     # symulowane dane 
 
# Dopasowujemy parametry 
# nie musimy podawać wartości startowych (params_init) dla procedury minimalizacji (wtedy funkcja zakłada wartości startowe równe 1)
# jednak zazwyczaj dobrze jest podpowiedzieć algorytmowi, gdzie powinien zacząć
# nie musimy również podawać wartości sigma, ale jeśli są one różne dla różnych punktów, to podanie ich sprawi, że algorytm będzie się bardziej troszczył
# o dopasowanie do punktów pomiarowych zmierzonych z dobrą dokładnością, a bardziej swobodnie podejdzie do tych o dużych niepewnościach
params_init = [2.0, -1.0]
params_final, covar = opt.curve_fit(funkcja_do_fitowania,X,Y,params_init,sigma)

print("Dopasowane parametry",params_final)
print("Macierz kowariancji\n",covar)

# dopasowane parametry
amp=params_final[0]
wykladnik=params_final[1]

# standardowe błędy dopasowania
amp_err = np.sqrt(covar[0][0])
wykladnik_err = np.sqrt(covar[1][1])
 
# test chi2 dobroci dopasowania. 
# Jeśli znamy wariancję błędu pomiarowego można zastosować test chi2 do oceny jakości dopasowania.
# 	Po pierwsze powinniśmy przetestować czy residua mają rozkład normalny
 
residua = funkcja_bledu(X, Y, funkcja_do_fitowania, params_final, sigma)# tym razem residua już są podzielone przez standardowe odchylenie, każde przez swoje
W, p =st.shapiro(residua)
print('Test normalności residuów: p = %.3f'%(p))
# jeśli tak to zmienna:
chi2 = np.sum(residua**2)
# podlega rozkładowi chi-kwadrat o N - n  ilości stopni swobody (n - ilość fitowanych parametrów), czyli u nas N-2	
# możemy zbadać jakie jest prawdopodobieństwo zaobserwowania takiej, bądź bardziej ekstremalnej wartości chi2:

N = len (X)
liczba_stopni_swobody = N-len(params_final) # liczba punktów - liczba parametrów
if chi2 < liczba_stopni_swobody:
	p_chi2 = st.chi2.cdf(chi2, liczba_stopni_swobody)
else:
	p_chi2 = st.chi2.sf(chi2, liczba_stopni_swobody) # równoważne 1-st.chi2.cdf(chi2, N-2), ale sf ma lepszą dokładność dla małych wartości
print('chi2 = %.2f, p_chi2 = %.3f' %(chi2, p_chi2))
# czasem używamy zredukowanego chi2, czyli podzielonego przez ilość stopni swobody
chi2_zred = chi2/liczba_stopni_swobody
# jeśli jest on znacząco większy niż 1 to model nie pasuje do danych, lub nie doszacowaliśmy sigmy, 
# jeśli jest sporo mniejszy niż 1 to prawdopodobnie oszacowane przez nas sigma jest większe niż rzeczywiste
# potrzebny jest nam rozkład chi2_zred:
N_dist = 100000
dist_chi2_zred = np.sum(st.norm.rvs(size=(liczba_stopni_swobody,N_dist))**2  ,0)/liczba_stopni_swobody
p_chi2_zred = np.sum(dist_chi2_zred>=chi2_zred)/float(N_dist)
print('chi2_zred = %.2f, p_chi2_zred = %.3f' %(chi2_zred, p_chi2_zred))
 
 
##########
# wykres
##########
py.subplot(2,1,1)
py.plot(X, funkcja_do_fitowania(X,amp,wykladnik))     # Fit
py.errorbar(X, Y, yerr=sigma, fmt='k.')  # Dane i błędy
py.text(5, 6.5, 'amplituda = %5.2f +/- %5.2f' % (amp, amp_err))
py.text(5, 5.5, u'wykładnik = %5.2f +/- %5.2f' % (wykladnik, wykladnik_err))
py.title(u'Dopasowanie metodą najmniejszych kwadratów')
py.xlabel('X')
py.ylabel('Y')
py.xlim(1, 11)
py.subplot(2,1,2)
py.plot(X, residua) # residua
py.xlabel('X')
py.ylabel('dY')
py.title(u'Wykres residuów')
py.show()

Dopasowanie wielomianu

Poniżej rozważymy przykład dopasowania zależności wielomianowej.

# -*- coding: utf-8 -*-
import scipy.stats as st
import pylab as py
import numpy as np
 
# funkcja używana do symulowania danych
def wielomian_z_szumem(x, wspolczynniki,blad_wzgledny):
	'''Definicja funkcji wielomianowej. Użyjemy jej do wytworzenia danych'''
	W = np.poly1d(wspolczynniki) # funkcja zwracająca obiekt wielomianu o zadanych wspolczynnikach
								 #można go używać tak, jak zwykłej funkcji, ale obsługuje też działania na wielomianach
	y = W(X)# idealne dane
	sigma = blad_wzgledny * y	# zakładamy, że stały jest błąd względny pomiaru
								# przeliczamy go na standardowe odchylenie symulowanego  błędu 
	# symulujemy szum z obliczonym odchyleniem standardowym i dodajemy go do danych idealnych										
	y += st.norm.rvs(size=num_points) * sigma       
	return (y, sigma)

def funkcja_bledu_dla_wielomianow(x, y, wspolczynniki, err):
	'''Suma kwadratów tej funkcji jest minimalizowana w procesie optymalizacji parametrów.
	Nam przyda się do obliczenia residuów.'''
	W = np.poly1d(wspolczynniki)
	y_fit = W(x) # aktualne wartości y z dopasowania
	residuum = y-y_fit # residua wchodzą do sumy kwadratów z wagą odwrotnie proporcjonalną do standardowego odchylenia
	residuum_wazone = residuum/ err 
	return residuum_wazone

# Generujemy punkty z szumem
num_points = 20
X = np.linspace(-4, 6, num_points)
wspolczynniki_wielomianu= (0.3,1,-2,4)
stopien_wielomianu=len(wspolczynniki_wielomianu)-1
blad_wzgledny_pomiaru=0.1
Y, sigma = wielomian_z_szumem(X, wspolczynniki_wielomianu, blad_wzgledny_pomiaru)     # symulowane dane 
 
# Dopasowujemy parametry
# tym razem skorzystamy z funkcji np.polyfit, która nie potrzebuje parametrów początkowych, ani zdefiniowanej funkcji, którą ma dopasować
# podajemy jej tylko nasze dane oraz stopień wielomianu, który ma dopasować oraz opcjonalne wagi
# UWAGA! Tym razem wagi muszą być odwrotnością odchyleń standardowych (1/sigma, a nie sigma, jak w curve_fit)
# funkcja ta domyślnie zwraca tylko dopasowane parametry (wspolczynniki wielomianu), a nie zwraca macierzy kowariancji,
# jeśli jest nam ona potrzebna, to musimy jej zarządać poprzez dodanie opcji cov=True (full=False, ale to jest domyślnie)
params_final, covar=np.polyfit(X, Y, deg=stopien_wielomianu, w=1/sigma, cov=True)

print("Dopasowane wspolczynniki wielomianu",params_final)
print("Macierz kowariancji\n",covar)

# standardowe błędy dopasowania
niepewnosci=[]
for i in range(len(params_final)):
	niepewnosci.append(np.sqrt(covar[i][i]))
print(niepewnosci)


# test chi2 dobroci dopasowania. 
# Jeśli znamy wariancję błędu pomiarowego można zastosować test chi2 do oceny jakości dopasowania.
# 	Po pierwsze powinniśmy przetestować czy residua mają rozkład normalny
 
residua = funkcja_bledu_dla_wielomianow(X, Y, params_final, sigma)# tym razem residua już są podzielone przez standardowe odchylenie, każde przez swoje
W, p =st.shapiro(residua)
print('Test normalności residuów: p = %.3f'%(p))
# jeśli tak to zmienna:
chi2 = np.sum(residua**2)
# podlega rozkładowi chi-kwadrat o N - n  ilości stopni swobody (n - ilość fitowanych parametrów), czyli u nas N-2	
# możemy zbadać jakie jest prawdopodobieństwo zaobserwowania takiej, bądź bardziej ekstremalnej wartości chi2:


N = len (X)
liczba_stopni_swobody = N-len(params_final) # liczba punktów - liczba parametrów
if chi2 < liczba_stopni_swobody:
	p_chi2 = st.chi2.cdf(chi2, liczba_stopni_swobody)
else:
	p_chi2 = st.chi2.sf(chi2, liczba_stopni_swobody) # równoważne 1-st.chi2.cdf(chi2, N-2), ale sf ma lepszą dokładność dla małych wartości
print('chi2 = %.2f, p_chi2 = %.3f' %(chi2, p_chi2))



# czasem używamy zredukowanego chi2, czyli podzielonego przez ilość stopni swobody
chi2_zred = chi2/liczba_stopni_swobody
# jeśli jest on znacząco większy niż 1 to model nie pasuje do danych, lub nie doszacowaliśmy sigmy, 
# jeśli jest sporo mniejszy niż 1 to prawdopodobnie oszacowane przez nas sigma jest większe niż rzeczywiste
# potrzebny jest nam rozkład chi2_zred:
N_dist = 100000
dist_chi2_zred = np.sum(st.norm.rvs(size=(liczba_stopni_swobody,N_dist))**2  ,0)/liczba_stopni_swobody
p_chi2_zred = np.sum(dist_chi2_zred>=chi2_zred)/float(N_dist)
print('chi2_zred = %.2f, p_chi2_zred = %.3f' %(chi2_zred, p_chi2_zred))
 
 
##########
# wykres
##########
py.subplot(2,1,1)
W=np.poly1d(params_final)
py.plot(X, W(X))     # Fit
py.errorbar(X, Y, yerr=sigma, fmt='k.')  # Dane i błędy
py.title(u'Dopasowanie metodą najmniejszych kwadratów')
py.text(-4.6, 92, u'dopasowane współczynniki = '+str(np.round(params_final,3)))
py.text(-4.6, 86, u'niepewności współczynników = '+str(np.round(niepewnosci,3)))
py.text(-4.6, 80, u'prawdziwe współczynniki = '+str(np.round(wspolczynniki_wielomianu,3)))
py.xlabel('X')
py.ylabel('Y')
py.xlim(X.min()-1, X.max()+1)
py.subplot(2,1,2)
py.plot(X, residua) # residua
py.xlabel('X')
py.ylabel('dY')
py.title(u'Wykres residuów')
py.show()
py.show()