# Osmi dvocas

# STATISTIKA

# Generalna uvodna prica:

# Zadatke koje smo razmatrali u teoriji verovatnoca su oblika
# (a) Dat je prostor verovatnoca i verovatnoca P
# (b) Odrediti verovatnocu nekog dogadjaja A.

# U praksi je uglavnom(!) lako uociti skup mogucih ishoda i izdvojiti klasu podskupova koje i zovemo dogadjajima.
# Medjutim verovatnoca P nije unapred data. Cesto se (pomocu predznanja i intuicije) moze pretpostaviti
# da P pripada nekoj klasi P* - dopustivih verovatnosnih mera.

# Def: Uredjenu trojku (Omega, sigma-algebra A, P*), gde je Omega= skup mogucih ishoda, A - sigma-algebra nad Omega,
# P*= familija dopustivih verovatnoca, zovemo Statisticka struktura (model).

# Taj objekat i jeste od osnovni u matematickoj statistici.

# Primer: 
# Nesimetrican novcic se baca N puta, ako sa 1 oznacimo pojavu pisma, a sa 0 pojavu glave, onda je prostor ishoda:
# Omega={omega=(a_1, a_2, ..., a_N) | a_k=0 ili a_k=1, za k od 1 do N}.
# P*={P_u | 0<=u<=1}m gde je u - nepoznata verovatnoca pojave pisma.
# P_u=u^(a_1+...+a_N) * (1-u)^(N-a_1-...-a_N)

# Generalni zadatak Matematicke statistike je:
# da se osnovu registrovanih podataka, dakle na osnovu realizovanog dogadjaja, odabere ona verovatnoca
# iz klase dopustivih, koja ce najbolje opisivati slucajnu pojavu koja se razmatra ili da se ta klasa sto vise suzi.

# Statistike poretka
# Neka je (X_1, X_2,..., X_n) uzorak iz raspodele F i neka je (x_1, ..., x_n) jedna realizacija tog uzorka.
# Poredjamo brojeve x_1, ..., x_n po velicini. Dobijamo niz x_(1)<=x_(2)<=...<=x_(n). Definisemo X_(k)=x_(k) za svako k
# i niz X_(i) zovemo statistikama poretka.

# Primetimo da je:
# X_(1)=min(X_1,..., X_n)
# X_(n)=max(X_1,..., X_n).

# Empirijska funkcija raspodele:
# Def: Neka je (X_1, ..., X_n) uzorak iz raspodele F. Funkcija F_n:R-> Hi, gde je Hi- skup slucajnih velicina, data sa
# F_n(x)=1/n sum_{k=1}^n I{X_k<=x}
# zove se Empirijska funkcija raspodele.


# Histogrami:
# Histogram - graficka reprezentacija raspodele niza datih numerickih podataka.
# Stoga predstavlja ocenu-prikaz funkcije gustine.
# Prvi je to koristio Karl Pirson (Karl Pearson).

# Histogram se pravi na sledeci nacin:
# 1) Dobijene podatke sortirati.
# 2) Odabrati duzinu podela, nazovemo tu duzinu d.
# 3) Podelimo ceo interval (raspon (range) podataka) i podintervale duzine d.
# 4) Plotujemo dobijene intervale - na X-osi - vrednosti, na Y-osi broj X_k-ova koji su upali u dati interval.

# U R-u postoji ugradjena funkcija "hist()" i vec smo je koristili, ako se secate.
# Neki bitni argumenti funkcije su
# x- vektor koji plotujemo: hist(x)
# breaks="Sturges" - tu predajemo vektor sa krajevima svakog podintervala.
# I mnogi drugi argumenti postoje, ali nama su od kljucnog znacaja ova dva.
# Ko zeli da pogleda ostale argumente : help(hist).

# Jedno od pitanja koje prirodno nastaje: Kako odabrati duzine podintervala, sto smo oznacili sa d?
# Uglavnom se prvo odabere broj kategorija, tj koliko podintervala cemo imati.
# A to se racuna po sledecoj formuli
# k=[log_2(N)]+1
# gde je [] - ceo deo broja, N - broj opservacija, tj koliko imamo podataka.
# A dalje d se dobija prema sledecoj formuli:
# d=R/k, gde je R-range, odnosno raspon vektora.

# Napomena: u opstem slucaju rastojanje "d" ne mora da bude jednako na svakom podintervalu!

# 1)
# Tvrdi se da je prosecna minimalna cena bezolovnog benzina u Americi bila 1.35$. U reklamne svrhe
# kompanija zeli da pokaze kako je njihova cena niza. Da bi potkrepili svoju tvrdnju, statisticari iz
# firme su sakupili sledece podatke na osnovu slucajnog uzorka:
c(1.22, 1.37, 1.27, 1.20, 1.42, 1.41, 1.22, 1.24,
  1.28, 1.42, 1.48, 1.32, 1.40, 1.26, 1.39, 1.45,
  1.44, 1.49, 1.47, 1.47, 1.24, 1.34, 1.27, 1.35,
  1.34, 1.45, 1.49, 1.45, 1.23, 1.20, 1.42, 1.34,
  1.43, 1.21, 1.49, 1.36, 1.24, 1.20, 1.45,
  1.23, 1.25, 1.24, 1.35, 1.23, 1.39, 1.38,
  1.46, 1.48, 1.26, 1.36, 1.22, 1.46, 1.39,
  1.22, 1.29, 1.47, 1.24, 1.35, 1.21, 1.21)
# Napisati program u R-u koji racuna uzoracku sredinu i medijanu,
# uzoracki raspon, i iscrtava histogram nad zadatim podacima.

# Resenje
# uzoracka sredina iliti srednja vrednost, za to sluzi funkcija mean() - vec smo i koristili
cene=c(1.22, 1.37, 1.27, 1.20, 1.42, 1.41, 1.22, 1.24,
       1.28, 1.42, 1.48, 1.32, 1.40, 1.26, 1.39, 1.45,
       1.44, 1.49, 1.47, 1.47, 1.24, 1.34, 1.27, 1.35,
       1.34, 1.45, 1.49, 1.45, 1.23, 1.20, 1.42, 1.34,
       1.43, 1.21, 1.49, 1.36, 1.24, 1.20, 1.45,
       1.23, 1.25, 1.24, 1.35, 1.23, 1.39, 1.38,
       1.46, 1.48, 1.26, 1.36, 1.22, 1.46, 1.39,
       1.22, 1.29, 1.47, 1.24, 1.35, 1.21, 1.21)
mean(cene)

# Medijana - centralna statistika poretka u slucaju neparnog N
# u slucaju parnog - poluzbir centralnih statistika poretka 
length(cene)
# Vidimo da je paran broj, zato poluzbir centralnih:
cene_sort=sort(cene)
medijana=(cene_sort[length(cene)/2]+cene_sort[length(cene)/2+1])/2
medijana
# Ili smo mogli pozivom funkcije median() da dobijemo to isto.
median(cene)

# Uzoracki raspon pomocu funkcije range()
r=range(cene)
r
r[2]-r[1]
# ili pesacki
cene_sort[length(cene_sort)]
cene_sort[1]
raspon=cene_sort[length(cene_sort)]- cene_sort[1]
raspon

# Histogram:
k=floor(log(length(cene), base=2) ) +1
k
# to je broj kategorija
# sada duzine podela
d=raspon/k
d
podele=cene_sort[1]+0:k * d
podele
# Sada crtamo histogram
hist(cene_sort, breaks=podele )
# A mogli smo i nesortiran vektor da mu prosledimo
hist(cene, breaks = podele)

# A kada bismo pozvali hist() bez prosledjenih podela, mogao bi skroz drugacije da odradi.
hist(cene)
# Sto i jeste slucaj ovde
# Da bismo sve videli na jednom plotu:

# pomocu sledeceg reda zadajemo da u plot-prozoru smesti odmah dva plota.
par(mfrow=c(1,2))
hist(cene, breaks=podele)
hist(cene)


# 2)
# Prilikom proucavanja rasta dece, posmatra se obim glave deteta pri rodenju izrazeno u cm. Dobijeni
# su sledeci podaci:
c(33.1, 33.7, 33.7, 33.8, 33.4,
  33.9, 33.6, 33.4, 34.1, 34.2,
  34.5, 34.2, 34.6, 34.9, 34.8,
  34.0, 34.5, 34.2, 34.2, 34.7,
  34.7, 34.6, 34.3, 34.3, 34.2,
  35.1, 36.0, 35.8, 35.2, 35.6,
  36.1, 35.1, 35.3, 35.2)
# Napisati program u R-u koji iscrtava histogram i stablo-lisce dijagram nad zadatim podacima.

# resenje:
obim=c(33.1, 33.7, 33.7, 33.8, 33.4,
       33.9, 33.6, 33.4, 34.1, 34.2,
       34.5, 34.2, 34.6, 34.9, 34.8,
       34.0, 34.5, 34.2, 34.2, 34.7,
       34.7, 34.6, 34.3, 34.3, 34.2,
       35.1, 36.0, 35.8, 35.2, 35.6,
       36.1, 35.1, 35.3, 35.2)

length(obim)
k=floor(log(length(obim) , 2)) +1
# funkcije diff() vraca razlike susednih clanova vektora, u nasem slucaju vektor je duzine dva, pa stoga bas raspon vraca
d=diff(range(obim))/k

podele=sort(obim)[1]+0:k*d
hist(obim, breaks = podele)

# Stablo- lisce dijagram
stem(obim)

# Za domaci pronadjite uzoracku srednju vrednost i medijanu. Kao i kvantile.


# 3)
# Dati su sledeci rezultati ispita:
c(28, 27, 26, 25, 24, 23, 21, 21, 20, 19, 19, 18, 18, 18, 17, 17, 17, 17,
  16, 16, 16, 15.5, 15, 15, 15, 15, 14, 13, 13, 13, 13, 12, 12, 11, 11, 
  11, 11, 11, 10, 10, 10, 9, 9, 8, 7, 6, 5, 4, 0, 0, 25, 23, 21, 21, 21, 
  21, 20, 19.5, 19, 19, 18, 18, 17, 17, 17, 17, 16, 15, 15, 15, 14, 14, 
  14, 13.5, 13, 13, 12, 12, 10, 10, 9, 9, 9, 9, 8,  7, 7, 7, 7, 5, 5, 5, 
  5, 4, 3, 2)
# Histogram, uzoracku sredinu, uzoracku disperziju, uzoracko standardno odstupanje, medijanu, stablo-lisce, kvantile?
# Da li prepoznajete rezultate? :)

# resenje:
rezultati=c(28, 27, 26, 25, 24, 23, 21, 21, 20, 19, 19, 18, 18, 18, 17, 17, 17, 17,
            16, 16, 16, 15.5, 15, 15, 15, 15, 14, 13, 13, 13, 13, 12, 12, 11, 11, 
            11, 11, 11, 10, 10, 10, 9, 9, 8, 7, 6, 5, 4, 0, 0, 25, 23, 21, 21, 21, 
            21, 20, 19.5, 19, 19, 18, 18, 17, 17, 17, 17, 16, 15, 15, 15, 14, 14, 
            14, 13.5, 13, 13, 12, 12, 10, 10, 9, 9, 9, 9, 8,  7, 7, 7, 7, 5, 5, 5, 
            5, 4, 3, 2)

mean(rezultati)
median(rezultati)
range(rezultati)



k=floor(log(length(rezultati) , 2)) +1
d=diff(range(rezultati))/k
podele=sort(rezultati)[1]+0:k*d
podele
# valja proveriti da li sve vrednosti upadaju u te intervale.
hist(rezultati, breaks=podele)

stem(rezultati)

# Uzoracka disperzija (popravljena):
var(rezultati)

# Uzoracko standardno odstupanje:
sqrt(var(rezultati))


# Kvantili: 
quantile(rezultati)
# pokazuje koje vrednosti statistika poretka imamo na 0%, 25%, 50% (medijana)
# 75% i 100%.


# Specijalno mozemo da vidimo da rezultati prolaze i testove normalnosti 
# (kasnije cemo detaljnije o njima pricati)
ks.test(rezultati, "pnorm", mean(rezultati), sqrt(var(rezultati)))
shapiro.test(rezultati)
# Velike su p-vrednosti, zato mozemo da prihvatimo hipotezu o normalnosti podataka.

# S tim sto ne bi smo bas smeli direktno ks.test da primenimo vec neku drugaciju vrstu njega, jer smo ocenili parametre.
# Zato instaliramo sledeci paket:
install.packages("nortest")
library(nortest)
# i u tom paketu trazimo funkciju lillie.test
help(lillie.test)
lillie.test(rezultati)
# Takodje i tu dobijamo velike p-vrednosti, zato mozemo da prihvatimo hipotezu o normalnosti nasih podataka.