VEROVATNOĆA U R-U

Postoje ugrađene funkcije za funkciju raspodele, gustinu, funkcije kvantila, i generisanje slučajnih brojeva iz zadate raspodele.
Nazivi raspodela:unif,norm,pois,beta,gamma,binom,cauchy,chisq,exp...

Određeni prefiksi se dodaju na ime raspodele, u zavisnosti od toga šta hoćemo da izračunamo:

• funkcija raspodele - p
• gustina - d (kod diskretnih raspodela ovo je verovatnoća iz zakona raspodele da slučajna veličina uzme zadatu vrednost)
• kvantili - q - inverzna od funkcije raspodele (za verovatnoću p vraca q tako da je \(p=P\{X\geq q\}\))
• generisanje slučajnih brojeva iz neke raspodele - r (od random)

# tako za npr. normalnu raspodelu imamo pnorm,dnorm,qnorm,rnorm`

Binomna raspodela:B(n,p)

Za računanje pojedinačnih verovatnoća \(P\{X=x\}\): dbinom(x,size,prob,log=FALSE), gde je size=n, prob=p

dbinom(0, size = 10, prob = 0.5)

## [1] 0.0009765625

• Funkcija raspodele: pbinom(q,size,prob,lower.tail=T,log.p=F)

# lower.tail podrazumeva da tražimo P{X<=q}, da smo to postavili na FALSE, računalo bi se P{X>q}
#log.p=T značilo bi da vraća log(p) umesto p
pbinom(8, size = 10, prob = 0.3)

## [1] 0.9998563

• Funkcija kvanitla: q, takvo da je \(P\{X<=q\}=p\): qbinom(p,size,prob,lower.tail=T,log.p=F)

qbinom(0.99, size = 10, prob = 0.5)

## [1] 9

Uzorak elemenata iz binomne raspodele: rbinom(n,size,prob), gde je n obim uzorka

rbinom(20, 10, 0.5)

##  [1] 5 7 4 5 5 8 3 6 5 5 6 2 4 4 5 5 7 4 5 4

#crtanje verovatnoća:
#npr B(50, 0.4):
x <- 0:50
plot(
x,
dbinom(x, size = 50, prob = 0.4),
main = "Binomna raspodela B(50, 0.4)",
xlab = "k",
ylab = "binomne verovatnoce"
)

#bolje se vidi ako verovatnoće predstavimo stubićima:
plot(x,
dbinom(x, size = 50, prob = 0.4),
type = "h",
main = "") #h od histogram

#možemo dodati naslov, podnaslov, font i sve ostalo na plot:
plot(
x,
dbinom(x, size = 50, prob = 0.4),
type = "h",
main = "Binomna raspodela B(50, 0.4)",
xlab = "k",
ylab = "binomne verovatnoce"
)

#binomni koeficijent: choose(n,k)
choose(4, 2)

## [1] 6

#Normalna raspodela N(m,sigma^2)
# vrednost gustine f(x)
#dnorm(x,mean=0,sd=1,log=F), mean=m, sd=sigma)
dnorm(3, mean = 2, sd = 4) # f(3), gde je f gustina N(2,4^2) raspodele

## [1] 0.09666703

# vrednost funkcije raspodele u tački z
# pnorm(z,mean=0,sd=1,lower.tail=T,log.p=F)
pnorm(0)

## [1] 0.5

pnorm(3, mean = 4, lower.tail = F)

## [1] 0.8413447

pnorm(3,
mean = 4,
lower.tail = F,
log.p = T)

## [1] -0.1727538

# crtanje gustine:
x <- seq(-5, 5, 0.01)
plot(x, dnorm(x), type = "l", main = "", xlab = "", ylab = ""
)

# drugi način:
curve(dnorm(x), from = -5, to = 5)

# crtanje f-je raspodele:
plot(x, pnorm(x), type = "l")

# kvantili:
# qnorm(p,mean=0,sd=1)
qnorm(0.9)

## [1] 1.281552

qnorm(pnorm(2)) #vraća baš 2, ovo su inverzne f-je

## [1] 2

#uzorak pseudoslučajnih brojeva iz N(m,sigma^2) raspodele:
#rnorm(size,mean=0,sd=1)
rnorm(20) #iz N(0,1)

##  [1]  0.197373361  1.496193840 -0.584803680 -0.645306118 -0.835235202
##  [6] -1.465491249  0.798979558 -1.098964950 -1.418150439  0.157138598
## [11]  1.298926643 -0.566736084  0.676489192  0.439335900  0.576869277
## [16] -0.009352715  0.319715809 -0.540779175 -0.055474444  0.310518705

#Uniformna raspodela: U[a,b]
#dunif(x,min=0,max=1,log=F), min=a,max=b
dunif(2)

## [1] 0

dunif(0.5)

## [1] 1

dunif(2, min = 3, max = 5)

## [1] 0

#punif(q,min=0,max=1,lower.tail=TRUE,log.p=FALSE)
punif(1, min = 0, max = 3)

## [1] 0.3333333

#qunif(p,min=0,max=1,lower.tail=TRUE,log.p=FALSE)
qunif(0.5)

## [1] 0.5

runif(10)#iz U[0,1]

##  [1] 0.65919378 0.41849826 0.40960594 0.19110584 0.07928669 0.33748221
##  [7] 0.18046183 0.33267440 0.76666183 0.85790404

#Eksponencijalna raspodela: Exp(lambda)
#dexp(x,rate=1,log=F) rate=lambda
#pexp(q,rate=1,lower.tail=T,log.p=F)
#qexp(p,rate=1,lower.tail=T,log.p=F)
dexp(2, rate = 2) # f(2), za Exp(2) raspodelu

## [1] 0.03663128

pexp(qexp(0.4))

## [1] 0.4

x <- seq(0, 25, 0.01)
y <- dexp(x, rate = 3)
plot(x,
y,
main = "Gustina Exp(3) raspodele",
ylab = "f(x)",
type = "l")

Aproksimacija verovatnoća na osnovu uzorka

# U ~ U(0,1)

n <- 10 ^ 5

u <- runif(n)

mean(u <= 0.3)  # P{U<=0.3}

## [1] 0.29902

mean(u > 0.4)   # P{U>0.4}

## [1] 0.60176

1 - punif(0.4)  # P{U>0.4}=1-P{U<=0.4} Tačna vrednost!

## [1] 0.6

mean(u == 1)    # P{U=1}

## [1] 0

# X ~ N(0,10)

x <- rnorm(n, sd = sqrt(10))
mean(x <= 0)   # P{X<=0}

## [1] 0.50052

pnorm(0, sd = sqrt(10))

## [1] 0.5

# Z ~ Pois(3)

z <- rpois(n, lambda = 3)

mean(z == 2) # P{Z=2}

## [1] 0.22622

dpois(2, lambda = 3) # P{Z=2} Tačna vrednost!

## [1] 0.2240418

Generisanje apsolutno neprekidnih slučajnih veličina (Metod inverzne transformacije)

Zadatak: Neka je slučajna veličina \(X\) apsolutno neprekidnog tipa sa funkcijom raspodele \(F\). Tada slučajna veličina \(F(X)\) ima \(\mathcal{U}[0,1]\) raspodelu. Takođe, ako \(U\) ima \(\mathcal{U}[0,1]\) raspodelu i postoji inverz \(F^{-1}\), tada \(F^{-1}(U)\) ima funkciju raspodele \(F\).

#pretpostavljamo da znamo da generišemo slučajan broj U iz U[0,1]

#ovo tvrđenje nam daje mogućnost da generišemo slučajne brojeve iz onih
#apsolutno neprekidnih raspodela za koje se inverz funkcije raspodele ne može izraziti analitički

• Eksponencijalna raspodela

Exp_random <- function(lambda)
{
u <- runif(1)
x <- (-1 / lambda) * log(1 - u)
return(x)
}
# Možemo da napravimo funkciju koja generiše uzorak obima N iz eksponencijalne raspodele.

Exp_random_samp <- function(lambda, N)
{
u <- runif(N)
x <- (-1 / lambda) * log(1 - u)
return(x)
}
Exp_random_samp(0.2, 10)

##  [1] 5.8868444 1.6683383 3.3203804 2.4558121 3.2079470 4.5254784 4.3837781
##  [8] 6.1809803 0.1688166 3.5565010

• Erlangova raspodela \(\gamma(\lambda,n)\) (zbir \(n\) nezavisnih eksponencijalnih slučajnih veličina sa parametrom \(\lambda>0\))

Erlang_random<-function(n,lambda)
{
  x<-0
  for(i in 1:n)
  {
    gama<-runif(1)
    x<-x+(-1/lambda)*log(1-gama)
  }
  return(x)
}
Erlang_random(10,2)

## [1] 4.289066

#Još neke apsolutno neprekidne raspodele koje se mogu modelirati ovom metodom su
#Vejbulova, Frešeova i Gumbelova
#Pogledajte njihove funkcije raspodele i napišite funkcije za
#generisanje slučajnih brojeva iz tih raspodela.
#Ove tri raspodele su važne jer su jedine tri nedegenerisane granične raspodele
#normiranog (nekim konstantama) niza maksimuma nezavisnih slučajnih veličina.

Neki primeri apsolutno neprekidnih raspodela

# Neka je X ~ Exp(a). Generisati uzorke iz sledećih raspodela i izračunati srednju vrednost.
# a) Y = |1 - X|
# b) Z = min{X, X^2}
# v) T = [X]


# a)
simulacija_a <- function(n, a)
{
  x = rexp(n, a)
  y = abs(1-x)
  return(y)
}

vrednost = simulacija_a(10000, 2)

# srednja vrednost
mean(vrednost)

## [1] 0.6342722

# b)
simulacija_b <- function(n, a)
{
  x = rexp(n, a)
  # pmin() - uzima član po član minimume.
  z = pmin(x, x^2)
  
  return(z)
}

vrednost = simulacija_b(10000, 2)

# srednja vrednost
mean(vrednost)

## [1] 0.3612592

# v)
simulacija_v <- function(n, a)
{
  x = rexp(n, a)
  # floor - isto što i ceo deo.
  t = floor(x)
  
  
  return(t)
}
vrednost = simulacija_v(10000, 2)

# srednja vrednost
mean(vrednost)

## [1] 0.157

Aproksimacije binomne raspodele

# Posmatrajmo eksperiment bacanja lopte u koš. Neka imamo ukupno 300 bacanja 
# i verovatnoća pogotka je 1/5.
# Izračunati verovatnoću da je bilo najviše 70 pogotka.

# Neka je X~Bin(300, 1/5). Traži se P{X<=70}=
# P{(X- E(X))/sqrt(D(X)) <= (70- E(X))/ D(X)}=
# P{X*<= (70- 60)/sqrt(48)} , gde je sada  X* ~ N(0, 1)
vrednost <- (70 - 300 * 1 / 5) / (sqrt(300 * 1 / 5 * 4 / 5))
pnorm(vrednost, 0, 1) # ovde smo dobili vrednost nakon aproksimacije

## [1] 0.9255427

pbinom(70, 300, 1 / 5) # ovde prava vrednost, bez aproksimacije.

## [1] 0.9330246

# Relativno su bliske.

# Dok u slucaju n>30, a n*p<=10 aproksimiramo sa Puasonovom 
# raspodelom sa parametrom Poiss(n*p)

# Poznato je da u nekom gradu stanovnik ima bicikl sa verovatnoćom 0.02,
# a motor sa verovatnoćom 0.01, s tim što niko nema i bicikl i motor. 
# Izračunati verovatnoću da od 100 slučajno odabranih stanovnika broj
# onih koji poseduju bar jedno od ova dva prevozna sredstva bude između
# 2 i 6 uključujući i te brojeve.

# Zapravo prava vrednost je da aproksimiramo 
# bas tu binomnu B~Bin(100, 0.03) sa Puasovom(3), koju ćemu označiti sa
# B*, tj:
# P{2<=B<=6} =~ P{2<=B*<=6} = P{B*<=6} - P{B*<2} 
# = P{B*<=6} - P{B*<=1} = ppois(6, 3) - ppois(1, 3)

# Apsoksimirana vrednost
ppois(6, 3) - ppois(1, 3)

## [1] 0.7673432

# Prava vrednost.
pbinom(6, 100, 0.03) - pbinom(1, 100, 0.03)

## [1] 0.7741504

Zakon velikih brojeva

# Primer: Bacanje novčića n puta, pri čemu u svakom bacanju pismo pada sa 
# verovatnoćom p, nezavisno od ostalih bacanja. Neka je Xi slučajna veličina 
# koja uzima vrednost 1 ako je u i-tom bacanju palo pismo, 0 inače, za sve 
# i=1,..n. Hoćemo da proverimo (simulacijom) da li za ovakav niz slučajnih 
# veličina važi slabi zakon velikih brojeva.

# Treba da proverimo da li P{|Sn/n-p|>=eps}~0 ako je n veliko.

# ESn=nEX1=n*p, što znaci da je ESn/n=p

zvb <- function(n, p, eps = 0.1) {
  s <- sample(c(1, 0), n, replace = TRUE, prob = c(p, 1 - p))
  
  ifelse((abs(mean(s) - p)) >= eps, 1, 0)
  
}
mean(replicate(10000, zvb(100, 0.5)))

## [1] 0.0353

# Verovatnoća dosta varira u zavisnosti od izbora eps zbog nepreciznosti samih
# ocena, a kako znamo da su ove slučajne veličine nezavisne i imaju konačno očekivanje, na osnovu teoreme sa časa znamo da za takav niz važi zakon velikih brojeva.

Histogram kao ocena gustine

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

library ( lattice )
# Uzorak obima 500 iz standardne normalne raspodele 
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 raspodele 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) raspodele 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')

CENTRALNA GRANIČNA TEOREMA

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

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

# Za domaći: Uraditi ovaj zadatak za različite primere raspodela koje zadovoljavaju
# uslove CGT.

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) raspodelom (tj. zbir n nezavisnih indikatora) i takođe \(ESn=np, DSn=np(1-p)\). Ako je \(np>10\) možemo je aproksimirati normalnom.

# Generisaćemo B(n,p) kao zbir n nezavisnih Bernulijevih slučajnih veličina 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 = "")  # ocena gustine iz uzorka
curve(dnorm(x), col='lightblue', lwd=2, add=T, from = -4, to =4) # teorijska kriva gustine N(0,1) raspodele