In [1]: 
from pylab import *
from scipy import stats
%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'])

# š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)   

# Gaussova porazdelitev
gauss=lambda x, m, s: (1 / (s * sqrt(2 * pi)) * exp(-0.5 * ((x - m) / s) ** 2))

Povprečja

Porazdelitev vrednosti spremenljivke je lahko zelo pestra in bogata, kot smo videli na naših zgledih. Kadar si hočemo to sliko poenostaviti na nekaj številskih parametrov, zelo pogosto uporabljamo povprečja.

Kot zgled najprej naračunamo 1000 dogodkov (dogodek je spet vsota 100 metov kocke) in narišemo histogram verjetnostne gostote dogodkov:

In [2]: 
# dogodek1 = vsota 100 metov kock
def dogodek1():
    # izračunaj vsoto 100 naključnih celih števil med 1 in 6
    return sum(randint(1, 7, 100))

# N dogodkov1
def N_dogodkov1(N):
    # izračunaj N naključnih dogodkov
    return [dogodek1() for i in xrange(N)]

dogodki1=N_dogodkov1(1000)

figure(figsize=(10,3));
subplot(121);
plot(dogodki1, '.');
xlabel(r'Zaporedna \v st. dogodka');
ylabel(r'Vsota 100 metov kocke');
subplot(122);
hist(dogodki1, stevilo_predalckov(dogodki1), normed=True, alpha=0.5);
xlabel(r'Vsota 100 metov kocke');
ylabel(r'Verjetnostna gostota');
tight_layout();
Porazdelitev bomo opisali s tremi parametri (\(N\) je število izmerjenih vrednosti, \(x_i\) so posamezne izmerjene vrednosti):

  • oceno povprečja spremenljivke (ang. sample mean), to je aritmetično sredino izmerjenih vrednosti spremenljivke: \(\bar x=\frac{1}{N}\sum_{i=1}^N x_i\),

  • oceno standardnega odklona ali standardne deviacije (ang. sample standard deviation), ki opiše razsutost ali širino porazdelitve: \(\sigma=\sqrt{C_\sigma\frac{1}{N}\sum_{i=1}^N(x_i-\bar x)^2}\); \(C_\sigma=\frac{N}{N-1}\).

  • in oceno poševnosti (ang. sample skewness), ki meri asimetrijo porazdelitve: \(\gamma=C_\gamma\frac{1}{N}\sum_{i=1}^N\left(\frac{x_i-\bar x}{\sigma}\right)^3\); \(C_\gamma=\frac{N^2}{(N-1)(N-2)}\). Pri negativni poševnosti je levi rep porazdelitve daljši. Obratno je pri pozitivni poševnosti daljši desni rep porazdelitve.

Koeficienta \(C_\sigma\) in \(C_\gamma\) poskrbita za popravek, ker obravnavamo samo naključni vzorec (ang. sample) \(N\)-tih vrednosti iz porazdelitve. Če imamo na voljo celotno porazdelitev, sta \(C_\sigma=C_\gamma=1\). V tem primeru dobimo prave vrednosti povprečja, standardnega odklona in poševnosti porazdelitve (ang. population mean, ...)

Za spremenljivko, ki smo ji določili \(\bar x\) in \(\sigma\), si v najbolj grobem približku predstavljamo, da je porazdeljena po Gaussovi (normalni) porazdelitvi \(G(x)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{\left(x-\bar x\right)^2}{2\sigma^2}}\). Ta je simetrična glede na vrh pri \(x=\bar x\), njena poševnost je torej enaka 0. Kljub temu lahko pri končnem številu meritev ocena poševnosti odstopa od te teoretične vrednosti tudi za spremenljivko, katere porazdelitev je Gaussova. Teorija pove, da lahko z veliko gotovostjo (približno 95%) trdimo, da porazdelitev ni Gaussova, če je \(\left|\gamma\right|>2\sigma_\gamma\), kjer je \(\sigma_\gamma=\sqrt\frac{6N(N-1)}{(N-2)(N+1)(N+3)}\).

Podobno je s povprečjem porazdelitve. Ocene povprečja vzorca iz Gaussove porazdelitve so okoli pravega povprečja porazdeljene, standardni odklon te porazdelitve je \(\sigma_\bar x=\frac{\sigma}{\sqrt{N}}\). Ta izraz lahko uporabimo kot oceno napake izračunanega povprečja (ang. standard error of the mean). Vidimo, da kvaliteta ocene povprečja narašča s povečevanjem števila meritev.

Izračunajmo torej parametre porazdelitve naših dogodkov in narišimo pripadajočo Gaussovo porazdelitev (rdeča črta):

In [3]: 
N=len(dogodki1)
print "Ocena povprečja:", stats.tmean(dogodki1)
print "Ocena napake povprečja:", stats.sem(dogodki1)
print "Ocena standardnega odklona:", stats.tstd(dogodki1)
print "Ocena poševnosti:", stats.skew(dogodki1, bias=False)
print "Standardni odklon porazdelitve poševnosti za Gaussovo porazdelitev:", sqrt(6.0 * N * (N - 1) / (N - 2) / (N + 1) / (N + 3))

ax=subplot()
hist(dogodki1, stevilo_predalckov(dogodki1), normed=True, alpha=0.5);
gauss=lambda x, m, s: (1 / (s * sqrt(2 * pi)) * exp(-0.5 * ((x - m) / s) ** 2))
x=linspace(100, 600, 1000);
ax.set_autoscalex_on(False)
plot(x, gauss(x, stats.tmean(dogodki1), stats.tstd(dogodki1)), 'r')
xlabel(r'Vsota 100 metov kocke');
ylabel(r'Verjetnostna gostota');

Ocena povprečja: 349.62
Ocena napake povprečja: 0.533948219119
Ocena standardnega odklona: 16.8849252501
Ocena poševnosti: -0.0386095275642
Standardni odklon porazdelitve poševnosti za Gaussovo porazdelitev: 0.0773438156196

Na kratko: povprečje porazdelitve je \(349.6\pm0.5\), njen standardni odklon pa \(16.9\). Porazdelitev ni poševna.

Poglejmo si še, kako so porazdeljena povprečja in poševnosti 1000 vzorcev (vsak vzorec vsebuje po 1000 dogodkov). Vidimo, da je širino obeh porazdelitev dobro opišeta zgoraj izračunani vrednosti.

In [4]: 
povprecja=[stats.tmean(N_dogodkov1(1000)) for i in xrange(1000)];
posevnosti=[stats.skew(N_dogodkov1(1000)) for i in xrange(1000)];
figure(figsize=(10, 3));
subplot(121);
hist(povprecja, stevilo_predalckov(povprecja), normed=True, alpha=0.5);
xlabel(r'Povpre\v cje vzorca 1000 dogodkov');
ylabel(r'Verjetnostna gostota');
subplot(122);
hist(posevnosti, stevilo_predalckov(posevnosti), normed=True, alpha=0.5);
xlabel(r'Po\v sevnost vzorca 1000 dogodkov');
ylabel(r'Verjetnostna gostota');
tight_layout();

Ponovimo zdaj račun za porazdelitev drugačnih dogodkov (harmonično povprečje 10 metov kocke):

In [5]: 
# dogodek2 = harmonično povprečje 10 metov kocke
def dogodek2():
    # izračunaj hamonično povprečje 10 naključnih celih števil med 1 in 6
    return stats.hmean(randint(1, 7, 10))

# N dogodkov2
def N_dogodkov2(N):
    # izračunaj N naključnih dogodkov
    return [dogodek2() for i in xrange(N)]

dogodki2=N_dogodkov2(1000)

figure(figsize=(10,3));
subplot(121);
plot(dogodki2, '.');
xlabel(r'Zaporedna \v st. dogodka');
ylabel(r'Harm. povpr. 10 metov kocke');
subplot(122);
hist(dogodki2, stevilo_predalckov(dogodki2), normed=True, alpha=0.5);
xlabel(r'Harm. povpr. 10 metov kocke');
ylabel(r'Verjetnostna gostota');
tight_layout();
In [6]: 
N=len(dogodki2)
print "Ocena povprečje:", stats.tmean(dogodki2)
print "Ocena napake povprečja:", stats.sem(dogodki2)
print "Ocena standardnega odklona:", stats.tstd(dogodki2)
print "Ocena poševnosti:", stats.skew(dogodki2, bias=False)
print "Standardni odklon porazdelitve poševnosti za Gaussovo porazdelitev:", sqrt(6.0 * N * (N - 1) / (N - 2) / (N + 1) / (N + 3))

ax=subplot()
hist(dogodki2, stevilo_predalckov(dogodki2), normed=True, alpha=0.5);
gauss=lambda x, m, s: (1 / (s * sqrt(2 * pi)) * exp(-0.5 * ((x - m) / s) ** 2))
x=linspace(1, 6, 1000);
ax.set_autoscalex_on(False)
plot(x, gauss(x, stats.tmean(dogodki2), stats.tstd(dogodki2)), 'r')
xlabel(r'Harm. povpr. 10 metov kocke');
ylabel(r'Verjetnostna gostota');

Ocena povprečje: 2.59144556051
Ocena napake povprečja: 0.0190889471981
Ocena standardnega odklona: 0.603645512808
Ocena poševnosti: 0.985425840089
Standardni odklon porazdelitve poševnosti za Gaussovo porazdelitev: 0.0773438156196

Vidimo, da je ocenjena poševnost porazdelitve prevelika, da bi porazdelitev lahko bila Gaussova. Pozitivna poševnost se kaže v tem, da je desni rep porazdelitve daljši. Na kratko: povprečje porazdelitve je \(2.59\pm0.02\), njen standardni odklon je \(0.60\), poševnost pa \(1.0\).

Naloge

  1. Oceni povprečje, napako povprečja, standardni odklon in poševnost za spremenljivko v podatkih Ozadje.dat, ki prikazujejo rezultate meritve absorpcije rentgenskih žarkov (logaritem razmerja vpadnega in prepuščenega toka, drugi stolpec) brez merjenca, tako da pričakujemo konstantne ali skoraj konstantne vrednosti. Primerjaj histogram verjetnostne gostote z grafom Gaussove porazdelitve z enakim povprečjem in standardnim odklonom. Lahko trdiš, da porazdelitev ni Gaussova?

  2. V datoteki Maribor.zip so zbrani vremenski podatki z merilne postaje Maribor-letališče za obdobje od 1. 1. 1977 do 31. 1. 2015. Ponovi postopek iz 1. naloge za povprečno dnevno temperaturo (datoteka TG_STAID003331.txt) v tem obdobju.

  3. Oceni povprečje in standardni odklon porazdelitve vsote vrednosti pri metu dveh kock tako, da z generatorjem naključnih števil generiraš \(10^6\) metov dveh kock. Parametra poskusi oceniti brez shranjevanja vrednosti v tabelo. Namig: \(\frac{1}{N}\sum_{i=1}^N(x_i-\bar x)^2=\frac{1}{N}\sum_{i=1}^N x_i^2-{\bar x}^2\). Koliko oceni odstopata od teoretičnih vrednosti?