In [1]: 
from pylab import *
%matplotlib inline
import statsmodels.api as sm

# nastavitve za izris grafov (http://matplotlib.org/1.3.1/users/customizing.html)
rc('text', usetex=True)
rc('font', size=12, family='serif', serif=['Computer Modern'])
rc('xtick', labelsize='small')
rc('ytick', labelsize='small')
rc('legend', fontsize='medium')
rc('figure', figsize=(5, 3))
rc('lines', linewidth=2.0)
rc('axes', color_cycle=['k'])

Linearna regresija

Linearna zveza \(y=kx\) je najpreprostejša in najpogostejša zveza med dvema fizikalnima količinama, zlasti še, ker lahko tudi druge funkcijske odvisnosti v ozkem intervalu aproksimiramo z linearno zvezo: \(\Delta y=k\Delta x\). Vemo, da sorazmernostni koeficient \(k\) za majhne \(\Delta x\) limitira k odvodu \(\mathrm{d}y/\mathrm{d}x\).

Kadar moramo določiti koeficient v linearni zvezi dveh količin \(x\) in \(y\), je s stališča merske tehnike koristno izmeriti več kot en par \((x,y)\). Če določimo vrednosti \(y_i\), ki ustrezajo celi vrsti izbranih \(x_i\), \(i = 1, 2, ..., N,\) lahko z njimi določimo \(k\) z večjo natančnostjo. Obenem preverimo, ali linearna zveza zares velja, pa tudi, ali naša merska naprava dobro deluje v širšem intervalu spremenljivk. Zlasti se lahko zgodi, da ima merilec te ali druge količine premaknjeno ničlo. Tedaj se naša linearna zveza ne pokaže kot premica skozi izhodišče koordinatnega sistema, pač pa je nekoliko premaknjena. Zato velja, da iz vrste meritev \(x_i\), \(y_i\) nikoli ne računamo koeficienta kot povprečje kvocientov \(y_i/x_i\), pač pa vedno kot naklon premice, ki jo potegnemo skozi izmerjene točke.

Problem najboljše premice skozi dane merske točke je mogoče definirati na mnogo načinov, ki pa se vsi prevedejo na eno od dveh matematičnih oblik. Če smo za meritve že določili korelacijski koeficient \(R_{x,y}\), zadošča vedeti, da gre najboljša premica vedno skozi težišče oblaka točk \((\bar x, \bar y)\) in da ima naklon \(R_{x,y}\sigma_y/\sigma_x\).

Sicer pa lahko koeficienta \(k\) in \(n\) najboljše premice \(y = kx + n\) določimo tudi z elementarnim računom. Zahtevamo, naj bo vsota kvadratov razdalj točk od premice najmanjša (načelo najmanjših kvadratov – least squares principle):

\[S=\sum\left(y_i-kx_i-n\right)^2=\mathrm{min}.\]

V zaresni merilni tehniki se spodobi, da so meritve opremljene s statistično napako, torej \(y_i= y_i\pm\epsilon_i.\) Kadar se \(\epsilon_i\) od izmerka do izmerka znatno spreminja, je teža posameznih točk različna. Statistično neoporečen rezultat dobimo, če sumande utežimo z \(\epsilon_i^{-2}\):

\[S=\sum\left(\frac{y_i-kx_i-n}{\epsilon_i}\right)^2=\mathrm{min}.\]

To je druga, popolnejša oblika za najboljšo premico. V njej vsoto \(S\), ki smo jo minimizirali, običajno označimo s \(\chi^2\) in jo lahko uporabljamo tudi za diagnozo kvalitete ujemanja. Za dobro ujemanje je pričakovana vrednost \(\chi^2/N=1\pm\sqrt{2/N}\). Govorimo o testu hi-kvadrat.

Zgornji postopek lahko uporabimo tudi pri nekaterih nelinearnih zvezah med dvema količinama. Če bi npr. skozi točke radi potegnili najboljšo eksponentno funkcijo \(y = Ae^{-\lambda x}\), le-to najprej z logaritmiranjem pretvorimo v linearno zvezo: \(\ln{y}=\ln{A}-\lambda x.\) Vidimo, da lahko \(A\) in \(\lambda\) določimo z običajno linearno regresijo za spremenljivki \(\ln{y}\) in \(x\).

Druga posplošitev je, da poskusimo podatkom prilagoditi linearno kombinacijo poljubnih funkcij \(f_i(x)\) s konstantnimi koeficienti \(c_i,\)

\[y=c_1 f_1(x)+c_2 f_2(x)+c_3 f_3(x)\ldots,\]

tako da minimiziramo vsoto

\[\chi^2=\sum\left(\frac{y_i-c_1 f_1(x_i)-c_2 f_2(x_i)-c_3 f_3(x_i)...}{\epsilon_i}\right)^2=\mathrm{min}.\]

Kot prvi primer bomo poiskali najboljšo premico za podatke, prikazane na naslednjih štirih grafih.

In [2]: 
podatki=loadtxt('anscombe.dat')
fig, ((axa, axc), (axb, axd)) = subplots(2, 2, figsize=(10, 6), sharex='col', sharey='row')
axa.set_xlim([3, 20])
axc.set_xlim([3, 20])
axa.set_ylim([2, 13])
axb.set_ylim([2, 13])
axa.set_ylabel(r'$y$')
axb.set_ylabel(r'$y$')
axb.set_xlabel(r'$x$')
axd.set_xlabel(r'$x$')

x1=podatki[:, 0]
y1=podatki[:, 1]
axa.plot(x1, y1, 'ko');
axa.grid(alpha=0.5)
axa.annotate('(a)', xy=(0.03, 0.92), xycoords='axes fraction')

x2=podatki[:, 2]
y2=podatki[:, 3]
axb.plot(x2, y2, 'ko');
axb.grid(alpha=0.5)
axb.annotate('(b)', xy=(0.03, 0.92), xycoords='axes fraction')

x3=podatki[:, 4]
y3=podatki[:, 5]
axc.plot(x3, y3, 'ko');
axc.grid(alpha=0.5)
axc.annotate('(c)', xy=(0.03, 0.92), xycoords='axes fraction')

x4=podatki[:, 6]
y4=podatki[:, 7]
axd.plot(x4, y4, 'ko');
axd.grid(alpha=0.5)
axd.annotate('(d)', xy=(0.03, 0.92), xycoords='axes fraction')

subplots_adjust(hspace=0.075, wspace=0.05)

Podatki so namenoma pripravljeni tako, da je najboljša premica za vse štiri sete enaka. Primerjava premice in podatkov na naslednjih štirih grafih pokaže, da je opis podatkov s premico primeren le za podatke (a). Podatki (c) in (d) vsebujejo točko, ki močno odstopa od obnašanja ostalih točk (outlier). Posledično je premica, ki jo določimo z linearno regresijo, zavajajoča. Podatkov (b) pa s premico ne moremo opisati. Linearna regresija s polinomom druge stopnje pa te podatke opiše pravilno.

In [4]: 
fig, ((axa, axc), (axb, axd)) = subplots(2, 2, figsize=(10, 6), sharex='col', sharey='row')
axa.set_xlim([3, 20])
axc.set_xlim([3, 20])
axa.set_ylim([2, 13])
axb.set_ylim([2, 13])
axa.set_ylabel(r'$y$')
axb.set_ylabel(r'$y$')
axb.set_xlabel(r'$x$')
axd.set_xlabel(r'$x$')

x=linspace(3, 20, 18)

fit1=sm.OLS(y1, sm.add_constant(x1, prepend=True)).fit()
f1=lambda x: fit1.params[0] + fit1.params[1] * x
axa.plot(x, f1(x), 'r');
axa.plot(x1, y1, 'ko');
axa.grid(alpha=0.5)
axa.annotate('(a)', xy=(0.03, 0.92), xycoords='axes fraction')

fit2=sm.OLS(y2, sm.add_constant(x2, prepend=True)).fit()
f2=lambda x: fit2.params[0] + fit2.params[1] * x
axb.plot(x, f2(x), 'r');
fit2x=sm.OLS(y2, sm.add_constant(column_stack((x2, x2 ** 2)), prepend=True)).fit()
f2x=lambda x: fit2x.params[0] + fit2x.params[1] * x + fit2x.params[2] * x ** 2
axb.plot(x, f2x(x), 'b');
axb.plot(x2, y2, 'ko');
axb.grid(alpha=0.5)
axb.annotate('(b)', xy=(0.03, 0.92), xycoords='axes fraction')
axb.legend([r'$y=kx+n$', r'$y=ax^2+kx+n$'], loc=8)

fit3=sm.OLS(y3, sm.add_constant(x3, prepend=True)).fit()
f3=lambda x: fit3.params[0] + fit3.params[1] * x
axc.plot(x, f3(x), 'r');
axc.plot(x3, y3, 'ko');
axc.grid(alpha=0.5)
axc.annotate('(c)', xy=(0.03, 0.92), xycoords='axes fraction')

fit4=sm.OLS(y4, sm.add_constant(x4, prepend=True)).fit()
f4=lambda x: fit4.params[0] + fit4.params[1] * x
axd.plot(x, f4(x), 'r');
axd.plot(x4, y4, 'ko');
axd.grid(alpha=0.5)
axd.annotate('(d)', xy=(0.03, 0.92), xycoords='axes fraction')

subplots_adjust(hspace=0.075, wspace=0.05)

Na naslednjih treh grafih smo prilagodili najboljšo premico podatkom z napakami in izračunali \(\chi^2/N\). Na grafu (a) so napake tako majhne, da podatkov s premico ne moremo pravilno opisati. To nam potrdi vrednost \(\chi^2/N\), ki močno presega \(1+\sqrt{2/N}\). Na grafu (c) je vrednost \(\chi^2/N\) premajhna, kar nakazuje, da smo mogoče precenili napake meritev. Na grafu (b) je \(\chi^2/N\) v mejah med \(1\pm\sqrt{2/N}\).

In [5]: 
def chi2(x, y, e, f):
   return sum(((y - f) / e) ** 2); 

fig, ((axa, axb, axc)) = subplots(3, 1, figsize=(5, 9), sharex='col', sharey='row')
axa.set_xlim([3, 20])
axa.set_ylim([2, 13])
axb.set_ylim([2, 13])
axc.set_ylim([2, 13])
axa.set_ylabel(r'$y$')
axb.set_ylabel(r'$y$')
axc.set_ylabel(r'$y$')
axc.set_xlabel(r'$x$')

x=linspace(3, 20, 18)

e1=[0.5 for i in x1]
fit1=sm.WLS(y1, sm.add_constant(x1, prepend=True), [i ** -2 for i in e1]).fit()
f1=lambda x: fit1.params[0] + fit1.params[1] * x
axa.plot(x, f1(x), 'r');
axa.errorbar(x1, y1, yerr=e1, fmt='ko');
axa.grid(alpha=0.5)
axa.annotate('(a)', xy=(0.03, 0.92), xycoords='axes fraction')
h1=chi2(x1, y1, e1, f1(x1)) / len(x1)
axa.annotate(r'$\chi^2/N=' + '%.2f' % h1 + r'$', xy=(0.4, 0.1), xycoords='axes fraction')

e2=[1.0 for i in x1]
fit2=sm.WLS(y1, sm.add_constant(x1, prepend=True), [i ** -2 for i in e2]).fit()
f2=lambda x: fit2.params[0] + fit2.params[1] * x
axb.plot(x, f2(x), 'r');
axb.errorbar(x1, y1, yerr=e2, fmt='ko');
axb.grid(alpha=0.5)
axb.annotate('(b)', xy=(0.03, 0.92), xycoords='axes fraction')
h2=chi2(x1, y1, e2, f2(x1)) / len(x1)
axb.annotate(r'$\chi^2/N=' + '%.2f' % h2 + r'$ ($\sqrt{2/N}=0.43$)', xy=(0.4, 0.1), xycoords='axes fraction')

e3=[2.0 for i in x1]
fit3=sm.WLS(y1, sm.add_constant(x1, prepend=True), [i ** -2 for i in e3]).fit()
f3=lambda x: fit3.params[0] + fit3.params[1] * x
axc.plot(x, f3(x), 'r');
axc.errorbar(x1, y1, yerr=e3, fmt='ko');
axc.grid(alpha=0.5)
axc.annotate('(c)', xy=(0.03, 0.92), xycoords='axes fraction')
h3=chi2(x1, y1, e3, f2(x1)) / len(x1)
axc.annotate(r'$\chi^2/N=' + '%.2f' % h3 + r'$', xy=(0.4, 0.1), xycoords='axes fraction')

subplots_adjust(hspace=0.075, wspace=0.05)

Poglejmo si še, kako poiščemo najboljšo eksponentno funkcijo \(y=Ae^{-\lambda x}\) za podatke na grafu (a) spodaj. Vrednosti \(y_i\) najprej logaritmiramo in poiščemo najboljšo premico \(\ln{y}=\ln(A)-\lambda x,\) graf (b). Iz koeficientov te premice preberemo \(A\) in \(\lambda\) in narišemo ustrezno eksponentno funkcijo (c).

In [6]: 
fig, ((axa, axb, axc)) = subplots(3, 1, figsize=(5, 9), sharex='col', sharey='row')
axa.set_xlim([0, 10])
axa.set_ylim([0, 200])
axb.set_ylim([1, 6])
axc.set_ylim([0, 200])
axa.set_ylabel(r'$y$')
axb.set_ylabel(r'$\ln{y}$')
axc.set_ylabel(r'$y$')
axc.set_xlabel(r'$x$')

N=50
x=linspace(0, 10, N)
y=exp(log(200.0 * exp(-0.3 * x)) + np.random.normal(0, 0.05, N))

axa.plot(x, y, 'ko')
axa.legend(['podatki'])
axa.annotate('(a)', xy=(0.03, 0.05), xycoords='axes fraction')

axb.plot(x, log(y), 'ko')
fit=sm.OLS(log(y), sm.add_constant(x, prepend=True)).fit()
f=lambda x: fit.params[0] + fit.params[1] * x
axb.plot(x, f(x), 'r')
axb.legend(['podatki', r'$\ln{y}=\ln{A}-\lambda x$'])
axb.annotate('(b)', xy=(0.03, 0.05), xycoords='axes fraction')

axc.plot(x, y, 'ko')
axc.plot(x, exp(f(x)), 'r')
axc.legend(['podatki', r'$y=Ae^{-\lambda x}$'])
axc.annotate('(c)', xy=(0.03, 0.05), xycoords='axes fraction')
axc.annotate(r'$A=' + '%.3g' % exp(fit.params[0]) + r'$', xy=(0.75, 0.5), xycoords='axes fraction')
axc.annotate(r'$\lambda=' + '%.3g' % -fit.params[1] + r'$', xy=(0.75, 0.4), xycoords='axes fraction')

subplots_adjust(hspace=0.075, wspace=0.05)

Naloge

  1. V datoteki Maribor.zip so zbrani vremenski podatki z merilne postaje Maribor-letališče za obdobje od 1. 1. 1977 do 31. 1. 2015. Poišči premico, ki najbolje opiše odvisnost med vlažnostjo in oblačnostjo v letu 2014.

  2. Teorija kemijske kinetike napove za sigmoidno krivuljo iz podatkov adrenalin.dat naslednjo odvisnost \(F/F_\textrm{max}=c/(a+c)\), kjer pomeni \(a\) koncentracijo s polovičnim maksimalnim učinkom. Določi koeficienta \(F_\textrm{max}\) in \(a\). Pretvori v linearno zvezo – ena pot je uvedba recipročnih spremenljivk \(1/F\) in \(1/c\), druga pa je uvedba spremenljivke \(c/F\).

  3. Iz povprečnih dnevnih temperatur v obdobju od 1. 1. 2000 do 31. 12. 2009 na merilni postaji Maribor-letališče (Maribor.zip) oceni povprečno temperaturo za vsak mesec v letu (povprečna temperatura v januarjih v tem obdobju, ...). Oceni še napake teh povprečij.

    Časovno odvisnost mesečnih povprečij modeliraj s funkcijama

    \[T(t)=T_0+T_1\cos\left(2\pi\frac{t}{t_0}\right)\]

    in

    \[T(t)=T_0+T_1\cos\left(2\pi\frac{t-t_1}{t_0}\right),\]

    kjer je \(t\) čas, merjen od začetka zime, in \(t_0\) eno leto. Določi parametre \(T_0\), \(T_1\) in \(t_1\) ter primerjaj vrednosti \(\chi^2\) za obe funkciji.