Histogram

Histogram je grafička reprezentacija raspodjele niza datih numerickih podataka. Koristi se za ocjenu funkcije gustine raspodjele.

Histogram se pravi na sledeći način: 1) Dobijeni podaci se sortiraju. 2) Odabere se dužina podeoka d. 3) Podijeli se cio interval (raspon podataka) na podintervale duzine d. 4) Na x-osi se oznache ti dobijeni intervali, a na odgovarajuća vrijednost na y-osi je broj elemenata iz uzorka koji su upali u taj interval.

U R-u postoji ugradjena funkcija hist(). Dva važna argumenta funkcije su: x- vektor čiji histogram želimo da prikažemo breaks - tu predajemo vektor sa krajevima svakog podintervala.

Kad pravimo histogram prvo treba da odaberemo koliko ćemo kategorija (odnosno tih podintervala) da imamo. Taj broj dobijamo iz formule: \[ k=[\log_2(N)]+1 \]

k-broj kategorija, N-obim uzorka (veličina tog vektora)

Odavde možemo naći širinu svakog intervala po formuli \[ d=\frac{R}{k} \]

gdje je R- raspon uzorka (razlika najvećeg i najmanjeg člana range(x))

Primjer. Tvrdi se da je prosječna minimalna cena bezolovnog benzina u Americi bila 1.35$. U reklamne svrhe kompanija želi da pokaze kako je njihova cijna niza. Da bi potkrijepili svoju tvrdnju, statističari iz firme su sakupili sledeće podatke na osnovu slučajnog uzorka:

cijene <- 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)

Izračunati u R-u uzoračku sredinu i medijanu, uzorački raspon, i nacrtati histogram nad zadatim podacima.

mean(cijene) # ili sum(cijene)/length(cijene)

## [1] 1.340167

median(cijene) 

## [1] 1.35

range(cijene) # ili c(min(cijene),max(cijene))

## [1] 1.20 1.49

# Pravimo histogram:
n <- length(cijene)
k <- floor(log(n, base = 2)) + 1
d <- diff(range(cijene)) / k
k

## [1] 6

d

## [1] 0.04833333

# Sortiramo vektor:

cijene <- sort(cijene)

# Pravimo podjelu na intervale:


podjela <- cijene[1] + 0:k * d
podjela

## [1] 1.200000 1.248333 1.296667 1.345000 1.393333 1.441667 1.490000

hist(cijene,
breaks = podjela,
col = 'coral',
main = "")

# Hoćemo da uporedimo histogram koji se dobija ako ne zadamo sami podjele:

par(mfrow = c(1, 2))
hist(cijene,
breaks = podjela,
col = 'coral' ,
main = "Histogram cijena")
hist(cijene, col = 'lightblue' , main = "Histogram cijena")

Primjer. Dati su sledeći rezultati ispita:

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)

Odrediti: histogram, uzoračku sredinu, uzoračku disperziju, uzoračko standardno odstupanje, medijanu, kvantile.

# Uzoracka sredina
mean(rezultati)

## [1] 13.77604

# Uzoracka disperzija 
var(rezultati)

## [1] 38.12563

# Standardno odstupanje
sd(rezultati)

## [1] 6.174596

# Medijana
median(rezultati)

## [1] 14

# Kvantili
quantile(rezultati)

##   0%  25%  50%  75% 100% 
##    0    9   14   18   28

# pokazuje koje vrijednosti statistika poretka imamo na 0%, 25%, 50% (medijana)
# 75% i 100%.

k <- floor(log(length(rezultati) , 2)) + 1
d <- diff(range(rezultati)) / k
podjela <- sort(rezultati)[1] + 0:k * d
hist(rezultati,
breaks = podjela,
col = 'coral',
main = "Histogram rezultata")

Poligon frekvencija

Grafičko predstavljanje podataka koje ima za cilja razumijevanje oblika raspodjele, na sličan način kao histogram.
Da bismo konstruisali poligon frekvencija, treba na početku da podijelimo podatke na neke klase (intervale), baš kao kod histograma. Potom se označi tačka na grafiku čija x koordinata ima vrijednost sredine intervala, a y koordinata je frekvencija odgovarajuće klase. Poligon dobijamo tako što dužima spojimo te označene tačke. Možemo uključiti jednu klasu prije najmanje vrijednosti među podacima i jednu posle najveće. Na taj način dobićemo da poligon dodiruje x osu sa obje strane.

Pokazaćemo poligon na primjeru podataka iz baze discoveries koja sadrži podatke o broju velikih izuma i naučnih otkrića između 1860. i 1959. godine.

discoveries

## Time Series:
## Start = 1860 
## End = 1959 
## Frequency = 1 
##   [1]  5  3  0  2  0  3  2  3  6  1  2  1  2  1  3  3  3  5  2  4  4  0  2
##  [24]  3  7 12  3 10  9  2  3  7  7  2  3  3  6  2  4  3  5  2  2  4  0  4
##  [47]  2  5  2  3  3  6  5  8  3  6  6  0  5  2  2  2  6  3  4  4  2  2  4
##  [70]  7  5  3  3  0  2  2  2  1  3  4  2  2  1  1  1  2  1  4  4  3  2  1
##  [93]  4  1  1  1  0  0  2  0

Po formuli dobijamo da treba da imamo otprilike 7 kategorija, pa za širinu intervala možemo uzeti 2 i dodajemo još dvije klase na krajevima.

range(discoveries)

## [1]  0 12

# Granice intervala

cut.points <- seq(-2 , 14 , by = 2)

# Tabele frekvencija

disc.cut <- cut(discoveries, cut.points, right = FALSE)
disc.freq <- table(disc.cut)
cbind(disc.freq)

##         disc.freq
## [-2,0)          0
## [0,2)          21
## [2,4)          46
## [4,6)          19
## [6,8)          10
## [8,10)          2
## [10,12)         1
## [12,14)         1

# Crtamo histogram i poligon frekvencija

hist(discoveries, 
     breaks=cut.points, 
     col="slategray3", 
     border = "dodgerblue4",
     right=FALSE,
     xlab = "x-osa", 
     main = "Histogram")

plot(disc.freq,
type = "b",
col = "orange",
main = "Poligon frekvencija")

Da li u primjeru sa rezultatima ispita možete prepotpostaviti koju raspodjelu imaju podaci?

Histogram kao ocjena gustine

Nekoliko primjera u kojima ćemo da generišemo uzorak iz poznate raspodjele i da poredimo oblik histograma sa teorijskom gustinom. (Koristimo histogram iz paketa lattice)

library ( lattice )
# Uzorak obima 500 iz standardne normalne raspodjele 
x <- rnorm (500)
hist (x, prob =TRUE , xlab ="",ylab ="",col ='coral1',border ='bisque', main="N(0,1)")
curve ( dnorm (x ) ,add =TRUE , lwd =3, col ='cornflowerblue ')

# Uzorak obima 1000 iz eksponencijalne raspodjele sa parametrom 4
y <- rexp (1000 ,4)
hist (y, prob =TRUE , xlab ="",ylab ="",col ='cornflowerblue',border ='bisque' , main="Exp(1)")
curve ( dexp (x ,4) ,add =TRUE ,lwd =3, col ='coral1')

# Uzorak iz uniformne U(0 ,1) raspodjele obima 500
z <- runif (500)
hist (z, prob =TRUE , xlab ="",ylab ="",col ='darksalmon',border ='bisque', main="U(0,1)")
curve ( dunif (x),add =TRUE , lwd =3, col ='cadetblue')

Primjer. Neka je \(S_n=X_1+...+X_n\), gdje su \(X_i, i=1,...,n\) nezavisne slučajne veličine sa \(\mathcal{E}(2)\) raspodjelom. Generisati slučajnu veličinu \(S_n\). Generisati \(N=1000\) brojeva iz iste raspodjele kao \(S_n\) i nacrtati njihov histogram. Da li je rezultat očekivan? (Podjsetite se CGT.)

s_n <- function(n) {
  x <- rexp(n, rate = 2)
  sum(x)
  
}

N <- 1000
n <- 100
s <- replicate(N, s_n(n))
hist(s,
probability = T,
col = 'lightblue',
main = "")

# Znamo da je EX=1/lambda, Dx=1/lambda^2


EX <- 1 / 2
DX <- 1 / 4

# Hocemo da standardizujemo podatke

s.z <- (s - n * EX) / sqrt(n * DX)
hist(s.z,
probability = T,
col = 'lightblue',
main = "")
curve(dnorm(x),
add = T,
lwd = 2,
col = 'coral')

Specijalan slučaj centralne granične teoreme je Muavr-Laplasova teorema:

Poznato nam je da se binomna slucajna veličina B(n,p) moze predstaviti kao zbir n nezavisnih slucajnih velicina sa Bernulijevom(p) raspodjelom (tj. zbir n nezavisnih indikatora) i takođe \(ESn=np, DSn=np(1-p)\). Ako je \(np>10\) možemo je aproksimirati normalnom.

# Generisac1emo B(n,p) kao zbir n nezavisnih Bernulijevih slucajnih velicina Ber(p)
binom <- function(n, p) {
  x <- sample(c(0, 1), n, replace = TRUE, prob = c(1 - p, p))
  b <- sum(x)
  return(b)
  
  
}
uzorak <- replicate(1000, binom(50, 0.7))
# standardizujemo podatke
st.uzorak <- (uzorak - 50 * 0.7) / sqrt(50 * 0.7 * 0.3)

plot(density(st.uzorak), lwd=2, col='coral', main = "")  # ocjena gustine iz uzorka
curve(dnorm(x), col='lightblue', lwd=2, add=T, from = -4, to =4) # teorijska kriva gustine N(0,1) raspodjele

Za neke poznate familije raspodjela u R-u postoji funkcija koja generiše uzorak iz te raspodjele , računa funkciju raspodjele u zadatoj tački, gustinu (kod apsolutno neprekidnih) ili vjerovatnoću iz zakona raspodjele (kod diskretnih) i kvantile raspodjele. Ove funkcije imaju prefikse r, p, d i q, redom, uz odgovorajuće ima za konkretnu raspodjelu. Više o tome:

 help("distribution")

## starting httpd help server ... done

Box-plot

Koristi se za detekciju autlajera.

boxplot(discoveries,col = "pink")

range(discoveries, na.rm = T)

## [1]  0 12

median(discoveries, na.rm = T)   # ! na.rm izostavlja NA vrijednosti iz uzorka

## [1] 3

q1 <- quantile(discoveries, na.rm = T)[2]
q3 <- quantile(discoveries, na.rm = T)[4]
qr <- IQR(discoveries, na.rm = T)
q1

## 25% 
##   2

q3

## 75% 
##   4

q1 - 1.5 * qr # f1

## 25% 
##  -1

q3 + 1.5 * qr # f3

## 75% 
##   7

q1 - 3 * qr # F1

## 25% 
##  -4

q3 + 3 * qr # F3

## 75% 
##  10

max(discoveries[discoveries < q3 + 1.5 * qr], na.rm = T) # a3

## [1] 6

Empirijska funkcija raspodjele

Neka je \((X_1, ..., X_n)\) uzorak iz raspodele F. \[F_n(x)=\frac{1}{n} \sum\limits_{k=1}^n I\{X_k \leq x\}\]

je empirijska funkcija raspodele. Primjer.

x <- rpois(20, lambda = 3)
table(x)

## x
## 1 2 3 4 5 6 
## 3 3 4 6 3 1

Fn <- ecdf(x)
plot(Fn, main = "Empirijska funkcija raspodjele")

Teorema [Glivenko-Kanteli]:

\[P\{\lim\limits_{N\to \infty} \sup\limits_{x \in \mathbb{R}}|F_N(x)-F(x)|=0\}=1\]

Primjer.

x<-runif(1000)
Fn<-ecdf(x)
plot(Fn, main="Empirijska funkcija raspodjele")
curve(punif(x), from = -1, to= 2, add = T, col="blue")

Zadatak: Nacrtati histogram datih podataka. Na osnovu njegovog izgleda procijeniti koja raspodjela najbolje opisuje podatke i dodati krivu gustine te raspodjele na histogram. Podaci:

 [1] 0.20461484 1.08674044 0.03856004 0.02560224 0.07928007 0.02340203 0.19352116 1.23953988
  [9] 1.19405939 1.06024342 0.17020231 0.09084827 0.28421320 0.04786419 0.43709257 0.02565908
 [17] 0.54787534 0.60735597 0.65346201 1.47148633 0.68876187 0.80113815 0.92430191 1.24893579
 [25] 0.38591012 0.43479177 1.10081436 0.70354737 0.23738834 1.79175686 0.13811964 0.71256143
 [33] 0.50764806 0.81095308 1.47837866 0.37798366 0.11811460 0.75100736 0.20178296 0.08794913
 [41] 0.06401065 0.28292047 0.23669375 0.57242330 0.73547956 0.61130502 0.79899670 1.30857987
 [49] 0.79412217 0.15405063 0.70493104 0.13028673 0.08695839 0.70228993 0.06374205 1.47505467
 [57] 0.32026473 0.55152157 0.31690261 0.75443969 0.51264233 0.07227614 0.44484945 1.18322093
 [65] 2.19167193 0.09054506 0.65836795 0.03997710 1.35456828 0.03298052 0.44563128 0.59318928
 [73] 0.21952265 0.26990897 0.33004873 0.60127290 1.18151449 0.02101827 0.38136168 0.37967533
 [81] 0.72016394 0.07404668 0.84811450 0.20616675 0.13154191 1.37544235 0.31510046 0.37547524
 [89] 0.93996441 0.54370520 0.20721877 0.75727264 0.20043624 0.24302641 0.02929712 0.22780797
 [97] 0.03895056 0.15394162 0.06309514 0.42370784

Simulacije

Zadatak 1: Simulirati bacanje regularne kockice za igru. Izracunati frekvenciju pojavljivanja šestice u 1000 izvodjenja eksperimenta. Nacrtati histogram frekvencija.

sample(6, 1, prob = rep(1/6,6))

## [1] 1

s<-sample(6, 1000, prob = rep(1/6,6), replace = TRUE)
mean(s==6)  # ekvivalentan poziv sum(s==6)/length(s)

## [1] 0.153

hist(s, main = "")

Zadatak 2: Bacaju se istovremeno kocka za igru i novčić.
a) Ispisati skup svih ishoda.

kocka<-c(1,2,3,4,5,6)
novcic<-c("pismo", "glava")

matrica.ishoda<-outer(kocka, novcic, FUN="paste", sep="")
matrica.ishoda

##      [,1]     [,2]    
## [1,] "1pismo" "1glava"
## [2,] "2pismo" "2glava"
## [3,] "3pismo" "3glava"
## [4,] "4pismo" "4glava"
## [5,] "5pismo" "5glava"
## [6,] "6pismo" "6glava"

2. Da li su svi ishodi jednako vjerovatni? Ocijeniti vjerovatnocu ishoda “5pismo”.

# Upisujemo sve ishode u vektor

vektor.ishoda<-as.vector(matrica.ishoda)

# Pozicija ishoda pismo5 u vektoru ishoda

pismo5<-which(vektor.ishoda=="5pismo")

# Generisemo vektor od hiljadu elemenata iz skupa 1:12 na slucajan nacin

simulacije<-sample(length(vektor.ishoda), 1000,  replace = TRUE)

# Ocjena vjerovatnoce na osnovu frekvencija

mean(simulacije==pismo5)

## [1] 0.09

Zadatak 3: IgraČ baca strelicu na tablu na kojoj se nalaze tri koncentrična kruga poluprečnika 5, 10 i 15 cm. Najveći krug je upisan u kvadrat stranice 30 cm. Smatra se da je svaki dio table ravnopravan.

1. Napisati funkciju strelica() koja simulira opisanu igru. (Funkcija treba da vrati 1-ako je pogođen prvi krug, 2 za drugi i 3 za treći, -1 ako je strelica van krugova)

strelica<-function(){

# Postavljamo koordinatni pocetak u centar kruga i na slucajan nacin biramo x i 
# y koor- dinatu za pogodjenu tacku, (x,y) su iz [-15,15]x[-15,15]
  
  poz<-runif(2, -15, 15)
  
# racunamo rastojanje tacke od centra
  
  r<-norm(poz, type="2")  # ekvivalentno: sqrt(sum(poz^2))

  if(r<=5)              return(1)
    else if(r>5 & r<=10)  return(2)
      else if(r>10 & r<=15) return(3)
        else                  return(-1)      

}

2. Ocijeniti vjerovatnocu da je da je pogodjen 1., odnosno 2. krug.

# Koristimo funkciju replicate cijim pozivom 
# 1000 puta pozivamo funkciju strelica() i rezultate tih poziva smjestamo 
# u jedan vektor

simulacije<-replicate(1000, strelica(), simplify = "array")


# Napomena: Ako fji replicate proslijedimo argument simplify= F, vratice listu

table(simulacije)

## simulacije
##  -1   1   2   3 
## 216  92 256 436

frekvencije<-table(simulacije)/1000
frekvencije

## simulacije
##    -1     1     2     3 
## 0.216 0.092 0.256 0.436

# Ako izracunamo vjerovatnoce teorijski, gledajuci kao gemetrijsku vjerovatnocu 
# po formuli: P(A)=m(A)/m(kvadrat) , dobicemo priblizno:

#  -1           1             2         3

#  0.21460    0.08727    0.26180    0.43633

3. Igrači A i B naizmjenično bacaju po tri strelice na tablu. Na početku igre imaju po 20 poena. Ako se pogodi prvi krug gubi se 1 poen, drugi 3, treci 5 i van kruga 10 poena. Igru započinje igrač A, a završava se kada jedan od igrača izgubi sve poene. Simulirati igru. Napisati funkciju pikado() koja vraća 1 ako je pobijedio igrač A, odnosno 2 ako je pobijedio igrač B.

 # Pravimo pomocnu fju koja racuna koliko treba oduzeti poena igracu nakon jedne
 # partije. Prosledjujemo joj vektor sa oznakama pogodjenih polja


poeni<-function(x){
  
  # racunamo koliko je puta pogodjeno svako polje
  
  krug.1<-sum(x==1)
  krug.2<-sum(x==2)
  krug.3<-sum(x==3)
  kvadrat<-sum(x==-1)
  
  return(krug.1 + 3*krug.2 + 5*krug.3 +10*kvadrat)
  
}



# Hocemo da napisemo funckiju pikado() koja vraca 1 ako je pobijedio igrac A, a
# 2 u suprotnom.


pikado<-function(){
  
  poeni.A<-20
  poeni.B<-20
  
# Napomena:  
# Indikator da je igrac A na redu: U aritmeticki izrazima T ce imati vr 1, a F 0
  
  A.na.potezu<-TRUE
  
  
  while(poeni.A*poeni.B>0){
    
  # Oduzimamo poeni igracu koji je na potezu
    
  poeni.A<-poeni.A - A.na.potezu*poeni(replicate(3,strelica()))
  poeni.B<-poeni.B - (!A.na.potezu)*poeni(replicate(3,strelica()))
    
  # Mijenjamo igraca
    
  A.na.potezu<-!A.na.potezu
    
    
  }
  
  # Ko prvi izgubi poene izgubio je igru.
       
  ifelse(poeni.A<=0 ,return (2), return(1))
      
  
}

# Ocijeniti vjerovatnocu da pobijedi igrac koji zapocinje igru, odnosno igrac A.
# Sta se ocekuje?


s<-replicate(1000,pikado())


mean(s==1)

## [1] 0.21

mean(s==2)

## [1] 0.79