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.340167median(cijene) ## [1] 1.35range(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] 6d## [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.490000hist(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 0Po 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 ... doneBox-plot
Koristi se za detekciju autlajera.
boxplot(discoveries,col = "pink")
range(discoveries, na.rm = T)## [1] 0 12median(discoveries, na.rm = T) # ! na.rm izostavlja NA vrijednosti iz uzorka## [1] 3q1 <- quantile(discoveries, na.rm = T)[2]
q3 <- quantile(discoveries, na.rm = T)[4]
qr <- IQR(discoveries, na.rm = T)
q1## 25%
## 2q3## 75%
## 4q1 - 1.5 * qr # f1## 25%
## -1q3 + 1.5 * qr # f3## 75%
## 7q1 - 3 * qr # F1## 25%
## -4q3 + 3 * qr # F3## 75%
## 10max(discoveries[discoveries < q3 + 1.5 * qr], na.rm = T) # a3## [1] 6Ocjenjivanje parametara
Podsjetite se zadatka 18. Imali smo \(X_1,...,X_n\) iz uniformne \(\mathcal{U}[0,\theta]\) raspodjele. Za nepoznati parametar su bile predložene dvije ocjene: \(\hat{\theta}_1=2\overline{X}_n\) i \(\hat{\theta}_2=\frac{n+1}{n}X_{(n)}\). Pokazali smo da su obje ocjene nepristrane pa treba da uporedimo njihove disperzije. Poznato je da je popravljena uzoračka disperzija nepristrasna ocjena za disperziju, pa ćemo simulacijama da ocjenimo disperzije za \(\hat{\theta}_1\) i \(\hat{\theta}_2\).
# Uzmimo primjer uzorka iz U(0,5) raspodjele. Generisemo n brojeva:
n <- 1000
x <- runif(n, 0, 5)
# Ocjene parametra theta na osnovu uzorka x
theta1 <- 2 * mean(x)
theta2 <- (n + 1) / n * max(x)
# Hocemo da ocijenimo disperzije ovih ocjena
ocjene <- function(n) {
x <- runif(n, 0, 5)
theta1 <- 2 * mean(x)
theta2 <- (n + 1) / n * max(x)
return(c(theta1, theta2))
}
r <- replicate(1000, ocjene(100))
# Ocjene ocekivanje ocjena theta1 i theta2
rowMeans(r)## [1] 5.007768 4.997899# Obje su blizu 5 sto je ocekivano jer su ocjene nepristrasne
# Ocjene disperzija
var(r[1, ])## [1] 0.08026403var(r[2, ])## [1] 0.0028284# Druga je mnogo manja, pa je druga ocjena bolja
# Uporedite sa teorijskim disperzijama dobijenim na casu.Metod maksimalne vjerodostojnosti
Primjer. Naći ocjenu za nepoznati parametar \(\lambda\) kod Puasonove \(\mathcal{P}(\lambda)\) raspodjele metodom maksimalne vjerodostojnosti.
# Na casu smo dobila da je ocjena uzoracka sredina
x <- rpois(n, lambda = 3)
lambda.hat <- mean(x)
lambda.hat## [1] 2.975# ocjenu MMV mozemo dobiti i pozivom funkcije fitdistr() iz paketa MASS
#install.packages("MASS")
library(MASS)## Warning: package 'MASS' was built under R version 3.2.5fitdistr(x, "Poisson")## lambda
## 2.97500000
## (0.05454356)# Osim za Puasonovu , pozivom fitdistr() mozemo dobiti ocjenu i za druge klase raspodjela:
# "beta", "cauchy", "chi-squared", "exponential", "f", "gamma", "geometric", "log-normal", "lognormal", "logistic", "negative binomial", "normal", "Poisson", "t" , "weibull"Intervali povjerenja
- Naći 95% interval povjerenja za nepoznati parametar \(m\) na osnovu uzorka obima \(n=100\) iz \(\mathcal{N}(m,64)\) raspodjele.
# Prvo cemo generisati neki uzorak pa cemo posle smatrati m nepoznatim
x <- rnorm(100, 5, sd = 8)
sigma <- 8
beta <- 0.95
n <- 100
c <- qnorm((1 + beta) / 2)
int.pov <- c(mean(x) - c * sigma / sqrt(n), mean(x) + c * sigma / sqrt(n))
int.pov## [1] 3.031091 6.167034- U primjeru 1 pretpostaviti da je \(\sigma\) nepoznato.
# Prvo cemo generisati neki uzorak pa cemo posle smatrati m i sigma nepoznatim
x <- rnorm(100, 5, sd = 8)
# Sigma treba da ocijenimo
beta <- 0.95
n <- 100
# t-raspodjela ili Studentova raspodjela
c <- qt((1 + beta) / 2, n - 1)
int.pov <-
c(mean(x) - c * sqrt(var(x) / (n - 1)), mean(x) + c * sqrt(var(x)) / sqrt(n -1))
int.pov## [1] 3.323741 6.207091- Za isti primjer iz prethodna dva zadatka odrediti 95% interval za \(\sigma^2\) (gornji i donji).
x <- rnorm(100, 5, sd = 8)
beta <- 0.95
n <- length(x)
c <- qchisq(beta, n - 1)
gornji <- c((n - 1) * var(x) / c, "infty")
gornji## [1] "48.5341087895297" "infty"c <- qchisq(1 - beta, n - 1)
donji <- c(0, (n - 1) * var(x) / c)
donji## [1] 0.00000 77.62376Empirijska 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
## 0 1 2 3 4 5 6
## 2 2 1 3 5 5 2Fn <- 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")