from pylab import *
%matplotlib inline
# 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', frameon=False, fontsize='medium')
rc('figure', figsize=(5, 3))
rc('lines', linewidth=2.0)
rc('axes', color_cycle=['k'])
Histogrami¶
Histogram je vrsta grafa, ki se uporablja v statistiki za prikaz porazdelitve neke količine. Histogram narišemo tako, da na vodoravno os nanašamo predalčke, tj. različne vrednosti količine ali pa intervale teh vrednosti, na navpično os pa s pokončnimi stolpci predstavimo ustrezne frekvence - tj. kako pogosto vrednosti padejo v vsakega od predalčkov.
Najprej definirajmo nekaj funkcij, s katerimi bomo naračunali podatke, ki jih bomo prikazali na histogramu.
# met ene kocke
def met():
# izračunaj naključno celo število med 1 in 6
return randint(1, 7)
# met N kock
def N_metov(N):
# izračunaj N naključnih celih števil med 1 in 6
return randint(1, 7, N)
# dogodek = vsota metov N kock
N=100
def dogodek():
# izračunaj vsoto N naključnih celih števil med 1 in 6
return sum(randint(1, 7, N))
# N dogodkov
def N_dogodkov(N):
# izračunaj N naključnih dogodkov
return [dogodek() for i in xrange(N)]
print 'met:', met()
print '10 metov:', N_metov(10)
print 'dogodek:', dogodek()
print '10 dogodkov:', N_dogodkov(10)
Naračunajmo 1000 dogodkov,
dogodki=N_dogodkov(1000)
plot(dogodki, '.');
xlabel(r'Zaporedna \v st. dogodka')
ylabel(r'Vsota 100 metov kocke');
in narišimo histograma z različnima številoma predalčkov, 4 in 100. Opazimo, da je v prvem primeru predalčkov premalo, da bi uspešno prikazali porazdelitev, v drugem pa jih je preveč in se posledično porazdelitev izgubi v šumu.
figure(figsize=(5, 6))
subplot(211)
hist(dogodki, 4, alpha=0.5);
title(r'Premalo predal\v ckov')
xlabel(r'Vsota 100 metov kocke');
ylabel(r'\v Stevilo dogodkov');
subplot(212)
hist(dogodki, 100, alpha=0.5);
title(r'Preve\v c predal\v ckov')
xlabel(r'Vsota 100 metov kocke');
ylabel(r'\v Stevilo dogodkov');
tight_layout()
Število predalčkov¶
- Sturgesova formula: \(N_\mathrm{p}=[1+\log_2(N)]\), kjer je \(N_\mathrm{p}\) število predalčkov, \(N\) število podatkov in \([x]\) prvo celo število, ki je večje od \(x\).
- Pravilo Freedmana in Diaconisa: \(w=\frac{2\mathrm{IQR(x)}}{\sqrt[3]{N}}\), kjer je \(w\) širina predalčka, \(N\) število podatkov in \(\mathrm{IQR}(x)\) razlika med zgornjim in spodnjim kvartilom podatkov \(x\). Spodnji kvartil je vrednost, pri kateri je 25% podatkov manjših od te vrednosti. Zgornji kvartil pa je vrednost, pri kateri je 75% podatkov manjših od te vrednosti. \(\mathrm{IQR}(x)\) torej podaja značilno širino porazdelitve vrednosti podatkov.
Poglejmo si, kakšno število predalčkov predlagata ti dve pravili:
n=len(dogodki)
print 'Število podatkov:', n
# Sturgesova formula
np=int(ceil(1 + log2(n)))
print "Število predalčkov po Sturgesovi formuli:", np
# pravilo Freedmana in Diaconisa
q1=percentile(dogodki, 25)
q3=percentile(dogodki, 75)
iqr = q3 - q1
print 'IQR(x):', iqr
sirina = 2 * iqr * n ** (-1. / 3.)
np=int((max(dogodki) - min(dogodki)) / sirina)
print "Število predalčkov po pravilu Freedmana in Diaconisa:", np
Narišimo histogram s številom predalčkov po pravilu Freedmana in Diakonisa:
hist(dogodki, np, alpha=0.5);
xlabel(r'Vsota 100 metov kocke');
ylabel(r'\v Stevilo dogodkov');
Namesto števila dogodkov v vsakem od predalčkov lahko s histogramom prikažemo verjetnostno gostoto. V tem primeru je ploščina pod celotnim histogramom enaka 1, ploščina histograma na nekem intervalu pa predstavlja verjetnost, da je vsota metov znotraj tega intervala. Pri tem načinu prestavitve je bolj smiselno, da histogram narišemo s črto namesto s stolpci.
hist(dogodki, np, normed=True, histtype='step');
xlabel(r'Vsota 100 metov kocke');
ylabel(r'Verjetnostna gostota');
Kumulativni histogram¶
V kumulativnem histogramu je vrednost v predalčku enaka vsoti števila dogodkov v tem predalčku in v vseh predalčkih levo od njega. Vrednost v zadnjem predalčku je torej enaka številu vseh dogodkov.
n, bins, patches=hist(dogodki, np, cumulative=True, histtype='step');
patches[0].set_xy(patches[0].get_xy()[:-1]);
xlabel(r'Vsota 100 metov kocke');
ylabel(r'Kumulativno \v stevilo dogodkov');
Tudi s kumulativnim histogramom lahko predstavimo verjetnostno gostoto. Vrednost v predalčku je v tem primeru enaka verjetnosti, da bo dogodek v tem predalčku ali pa v enem od predalčkov levo od njega. Vrednost v zadnjem predalčku je enaka 1.
n, bins, patches=hist(dogodki, np, normed=True, cumulative=True, histtype='step');
patches[0].set_xy(patches[0].get_xy()[:-1]);
ylim([0,1])
xlabel('Vsota 100 metov kocke');
ylabel('Kumulativna verjetnost');
Kumulativni histogram je še posebej priročen za ocenjevanje verjetnosti, da je vrednost prikazane količine v nekem intervalu. Na spodnjem primeru iz kumulativnega histograma (a) brez težav ocenimo verjetnost, da je vsota 100 metov manjša od 330, na približno 20%. Ustrezna količina na običajnem histogramu (b) je ploščina histograma na tem intervalu in je vizualno ne moremo tako natančno oceniti.
figure(figsize=(5, 6))
subplot(211)
n, bins, patches=hist(dogodki, np, normed=True, cumulative=True, histtype='step');
patches[0].set_xy(patches[0].get_xy()[:-1]);
xlabel(r'Vsota 100 metov kocke');
ylabel(r'Kumulativna verjetnostna gostota');
grid();
axvline(330, color='r', linewidth=2.0, ls='--')
annotate(r'(a)', xy=(0.9,0.9), xycoords='axes fraction')
subplot(212)
hist(dogodki, np, normed=True, histtype='step');
xlabel(r'Vsota 100 metov kocke');
ylabel(r'Verjetnostna gostota');
grid();
axvline(330, color='r', linewidth=2.0, ls='--')
annotate(r'(b)', xy=(0.9,0.9), xycoords='axes fraction')
tight_layout()
Naloge¶
-
V datoteki Maribor.zip so zbrani vremenski podatki z merilne postaje Maribor-letališče za obdobje od 1. 1. 1977 do 31. 1. 2015. Nariši histogram povprečne dnevne temperature (datoteka TG_STAID003331.txt) v tem obdobju.
-
V enakomernih korakih tabeliraj funkcijo \(y=(1/\tau)\exp(-t/\tau)\), kjer je \(\tau=2{,}2\) sekunde. Funkcija podaja verjetnostno gostoto radioaktivnih razpadov atoma z življenskim časom \(\tau\). Nariši kumulativno verjetnost kot funkcijo časa. S pomočjo histograma oceni, v kakšnem času razpade polovica vseh atomov.
- Nariši histogram porazdelitve podatkov (logaritem razmerja vpadnega in prepuščenega toka, drugi stolpec) iz datoteke Fe_Co.dat, ki vsebuje nekoliko predelan absorpcijski spekter EXAFS (Extended X-ray Absorption Fine Structure = drobna struktura rentgenskih absorpcijskih robov) mešanega železovo-kobaltovega oksida. Čeprav je neodvisna spremenljivka (energija fotona v eV, prvi stolpec) netrivialna in pomembna fizikalna količina, je zanimivo pogledati podatke tudi neodvisno od energije. Izberi primerno število predalčkov. Porazdelitev ima vrhove pri skoraj konstantnih vrednostih absorpcije med robovi. — Zanimiva je še ena modifikacija: koraki v energiji, pri katerih smo merili absorpcijo, niso enako razmaknjeni (ekvidistantni). Zato v predalčenju vse točke niso enako pravično obravnavane. Predstavljajmo si, da opravimo meritev še enkrat z zelo drobnim ekvidistantnim energijskim korakom. Vidimo, da je pravična teža vsake točke velikost koraka, točneje interval, ki obsega pol desnega in pol levega energijskega koraka. S tem smo določili matematično mero vsake točke. Primerjaj porazdelitev po prejšnjem in novem načelu. Če je prva normirana s številom vseh točk, je druga normirana z dolžino energijskega intervala celega spektra.