In [1]: 
from pylab import *
%matplotlib inline
import pandas
from pandas.tools.plotting import scatter_matrix
from scipy.stats import pearsonr
from statsmodels.tsa.stattools import ccf
from statsmodels.tsa.stattools import acf
import numpy

# 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'])

# število predalčkov v histogramu (pravilo Freedmana in Diakonisa)
def stevilo_predalckov(x):
    iqr = percentile(x, 75) - percentile(x, 25)
    sirina = 2.0 * iqr * len(x) ** (-1.0 / 3.0)
    return int((max(x) - min(x)) / sirina)

Korelacija

Korelacijski koeficient

Korelacijski koeficient med dvema količinama, \(y_1\) in \(y_2\), ki sta podani kot \(N\)-razsežna vektorja, je

\[R_{y_1,y_2}=\frac{1}{N-1}\sum_{i=1}^{N}\left(\frac{y_{1i}-\overline{y_1}}{\sigma_{y_1}}\right)\left(\frac{y_{2i}-\overline{y_2}}{\sigma_{y_2}}\right)\]

kjer sta \(\overline{y_1}\) in \(\overline{y_2}\) povprečji obeh količin, \(\sigma_{y_1}\) in \(\sigma_{y_2}\) pa njuna standardna odklona.

Korelacijski koeficient je brezdimenzijska količina, ki lahko zavzame vrednosti med -1 in 1. Pozitivne vrednosti pomenijo, da sta spremenljivki korelirani, se pravi, da ena količina v povprečju raste, ko raste druga. Negativne vrednosti \(R\) pomenijo antikorelacijo, ko vrednosti prve količine v povprečju padajo, ko druge naraščajo. Pri tem so vrednosti na obeh konceh intervala rezervirane za jasne linearne odvisnosti količin; tiste v sredi, recimo med -0.3 in 0.3, pa za blago, komaj razločno povezavo med gibanjem obeh spremenljivk.

Poglejmo si korelacijske koeficiente med količinami, prikazanimi na naslednjih štirih grafih:

In [2]: 
N=100
x=linspace(1, N, N)
y1=-0.8 + 1.6 * x / N + np.random.normal(0, 0.05, N)
y2=-0.2 + 0.3 * (x / N) ** 2 + np.random.normal(0, 0.1, N)
y3=0.0 + 0.0 * (x / N) + np.random.normal(0, 0.1, N)
y4=0.6 - 1.2 * (x / N) + np.random.normal(0, 0.2, N)
p = numpy.random.permutation(len(y1))
y1 = y1[p]
y2 = y2[p]
y3 = y3[p]
y4 = y4[p]

fig, (ax1, ax2, ax3, ax4) = subplots(4, 1, sharex='col', sharey='row', figsize=(3, 9))
ax1.plot(x, y1, "ko", markersize=4.0)
ax1.set_ylabel(r'$y_1$')
ax1.set_ylim(-1, 1)
ax2.plot(x, y2, "ko", markersize=4.0)
ax2.set_ylabel(r'$y_2$')
ax2.set_ylim(-1, 1)
ax3.plot(x, y3, "ko", markersize=4.0)
ax3.set_ylabel(r'$y_3$')
ax3.set_ylim(-1, 1)
ax4.plot(x, y4, "ko", markersize=4.0)
ax4.set_ylabel(r'$y_4$')
ax4.set_ylim(-1, 1)
ax4.set_xlabel(r'$i$')
subplots_adjust(hspace=0.1)

Poskusimo najprej grafično predstaviti korelacije med temi količinami. To naredimo tako, da narišemo grafe, kjer eno od količin nanesemo na vodoravno os, drugo pa na navpično os (izvendiagonalni grafi na sliki spodaj, na diagonali so prikazani histogrami posameznih količin).

In [3]: 
y={'$y_1$': y1, '$y_2$': y2, '$y_3$': y3,'$y_4$': y4}
df=pandas.DataFrame(y)
scatter_matrix(df, alpha=0.5, figsize=(8, 8), color='k', hist_kwds={'alpha': 0.5});
subplots_adjust(hspace=0.1, wspace=0.1);

Vidimo, da sta količini \(y_1\) in \(y_2\) korelirani, količini \(y_1\) in \(y_4\) ter \(y_2\) in \(y_4\) pa antikorelirani. Med količino \(y_3\) in ostalimi količinami pa korelacije ni opaziti. Ta opažanja potrdi izračun korelacijskih koeficientov:

In [4]: 
pandas.set_option('float_format', '{:,.2f}'.format)
df.corr()
Out[4]:
$y_1$ $y_2$ $y_3$ $y_4$
$y_1$ 1.00 0.68 0.11 -0.86
$y_2$ 0.68 1.00 0.04 -0.49
$y_3$ 0.11 0.04 1.00 -0.11
$y_4$ -0.86 -0.49 -0.11 1.00

Korelacijski koeficienti med količino \(y_3\) in ostalimi količinami so res majhni. Lahko na tej osnovi trdimo, ta korelacije med temi količinami ni? Da se pokazati, da mora biti korelacijski koeficient po absolutni vrednosti večji od \(\frac{2}{\sqrt{N}}\), da lahko s 95% verjetnostjo trdimo, da sta količini korelirani. V zgornjem primeru, kjer je \(N=100\), torej o korelacijah med \(y_3\) in ostalimi količinami ne moremo govoriti.

Spodnji graf prikazuje porazdelitev korelacijskih koeficientov neke količine z različnimi vektorji z naključnimi elementi. Če bi bilo število elementov v vsakem od vektorjev neskončno, bi moral korelacijski koeficient biti enak 0. Pri končnem številu elementov \(N\) pa so korelacijski koeficienti porazdeljeni okoli 0 tako, da je v 95% primerov korelacijski koeficient po absolutni vrednosti manjši od \(\frac{2}{\sqrt{N}}\) (območje med rdečima črtama).

In [6]: 
N=1000
def kor1():
    y1=linspace(0, 1, N)
    y2=np.random.normal(0, 1, N)
    r,z=pearsonr(y1, y2)
    return r

kor=[kor1() for i in xrange(N)]
hist(kor, stevilo_predalckov(kor), normed=True, alpha=0.5);
axvline(-2/sqrt(N), color='r', linewidth=2.0, ls='--')
axvline(2/sqrt(N), color='r', linewidth=2.0, ls='--')
xlabel(r'Korelacijski koeficient');
ylabel(r'Verjetnostna gostota');

Korelacijska funkcija

Pomembna karakteristika povezave signalov \(y_1(x)\) in \(y_2(x)\) je njuna korelacijska funkcija

\[R_{y_1,y_2}(\Delta x) = R_{y_1(x),\,y_2(x+\Delta x)},\]

torej korelacijski koeficient pri zakasnitvi \(\Delta x\). Iz nje zelo jasno razberemo morebitne fazne premike signalov: zgodi se, da je skalarni produkt signalov majhen, da pa zrase, ko enega od signalov premaknemo po osi \(x\). Seveda postaja z večanjem premika \(\Delta x\) interval, po katerem lahko računamo korelacijski koeficient, vse krajši, zato je definicijo smiselno popraviti tako, da računamo povprečno vrednost produkta – se pravi, da pri vsakem \(\Delta x\) delimo rezultat s številom točk, na katerih se signala sploh še prekrivata.

Poglejmo si korelacijsko funkcijo signalov, prikaznih na spodnji sliki.

In [19]: 
N=1000
x=linspace(0, 5, N)
y1=sin(2 * pi * x) + np.random.normal(0, 0.3, N)
y2=-0.4 + 0.5 * sin(2 * pi * x - 2.4) + np.random.normal(0, 0.5, N)
y3=0.2 * np.random.normal(0, 1, N)
plot(x, y1, 'ko', x, y2, 'ro', x, y3, 'yo', markersize=3.0, alpha=0.4);
xlabel(r'x');
ylabel(r'y');
legend([r'$y_1$', r'$y_2$', r'$y_3$']);

Korelacijska funkcija signalov \(y_1(x)\) in \(y_2(x)\) je prikazana s črno barvo, korelacijska funkcija signala \(y_1\) s signalom \(y_3\), ki vsebuje naključen šum, pa je prikazana z rumeno barvo. Z rdečima črtama je prikazana meja območja med \(\pm\frac{2}{\sqrt{N}}\). Če je vrednost korelacijske funkcije pri nekem \(\Delta x\) izven tega intervala, smo lahko z vsaj 95% verjetnostjo prepričani, da sta pri tem zamiku signala korelirana.

In [20]: 
c12=ccf(y1, y2)
c13=ccf(y1, y3)
plot(x, c12, 'k', x, c13, 'y');
xlabel(r'$\Delta x$')
ylabel('Korelacijska funkcija');
ylim(-1.1, 1.1)
axhline(-2/sqrt(N), color='r', linewidth=2.0, ls='--');
axhline(2/sqrt(N), color='r', linewidth=2.0, ls='--');
legend([r'$y_1$ in $y_2$', r'$y_1$ in $y_3$'], ncol=2, frameon=False, loc=2);

Avtokorelacijska funkcija

Tako funkcijo je zanimivo izračunati celo za en sam signal, takrat jo imenujemo avtokorelacijska funkcija:

\[R_{y,y}(\Delta x)=R_{y(x),\,y(x+\Delta x)}\]

Pri periodičnih signalih lahko iz ekstremov \(R_{y,y}(\Delta x)\) zelo natančno odčitamo periodo signala, pri neperiodičnih pa zvemo, kako daleč se signal "spominja" prejšnjega poteka.

Na spodnjem grafu sta prikazani avtokorelacijski funkciji signala \(y_1\) s črno bravo in naključnega šuma \(y_3\) z rumeno bravo. Z rdečima črtama je prikazana meja območja med \(\pm\frac{2}{\sqrt{N}}\).

In [23]: 
c11=acf(y1, nlags=N, unbiased=True)
c33=acf(y3, nlags=N, unbiased=True)
plot(x, c11, 'k', x, c33, 'y');
xlabel(r'$\Delta x$')
ylabel('Avtokorelacijska funkcija');
ylim(-1.3, 1.3)
axhline(-2/sqrt(N), color='r', linewidth=2.0, ls='--');
axhline(2/sqrt(N), color='r', linewidth=2.0, ls='--');
legend([r'$y_1$', r'$y_3$'], ncol=2, frameon=False, loc=2);

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. Izračunaj korelacijske koeficiente med minimalnimi in maksimalnimi dnevnimi temperaturami, oblačnostjo, vlažnostjo ter dnevno količino padavin v letu 2014. Katere od teh količin so korelirane?

  2. Pred leti so v okviru mednarodnega projekta v našem znanem gradbenem podjetju merili hitrost strjevanja betona. Ulili so nekaj metrov velik betonski blok, v katerega je bila vdelana cela vrsta termočlenov za sprotno merjenje temperature. Datoteka Beton.dat podaja izmerke v razdobju šestih dni v dveh merilnih točkah. Prva je blizu površine, druga globoko v notranjosti. (Prvi stolpec je zaporedna številka meritve – časovni interval med njimi lahko oceniš iz očitnih dnevnih nihanj temperature.) Določi efektivno zakasnitev med obema signaloma iz njune korelacijske funkcije.

  3. Iz podatkov v datoteki Maribor.zip izračunaj avtokorelacijsko funkcijo oblačnosti za obdobje od 1. 1. 2000 do 31. 12. 2009 in komentiraj njeno obnašanje pri majhnih in velikih zakasnitvah.