WnioskowanieStatystyczne/ Regresja liniowa i test chi2

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].

Przykład: rozkłady Y dla każdego punktu X

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

# symuowana zależność ma następującą postać y = b0 + b1*x 
# wartości parametrów 
b0 = 1 
b1 = 3

X = np.arange(0, 10,0.5)
# będę symulował zbieranie n wartości Y dlakażdego X[i], zakładam to samo odchylenie standardowe 
odch_std = 1
n = 30
Y = np.zeros((n,len(X)))
for i in range(len(X)):
	Y[:,i] =  b0 + b1*X[i] + st.norm.rvs(size = n, loc=0, scale = odch_std)

# narysujmy ten zbiór punktów 
for j in range(len(X)):
	py.plot(X, Y[j,:],'b,')
# wyróżnimy średnie
py.plot(X,np.mean(Y,0),'ro')
#  i odchykenia standardowe:
py.errorbar(X,np.mean(Y,0),odch_std,ecolor = 'k',elinewidth = 8)
py.show()

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()