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 hocemo da izracunamo:

• funkcija raspodele - p
• gustina - d (kod diskretnih raspodela ovo 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 trazimo P{X<=q}, da smo to postavili na FALSE, racunalo bi se P{X>q}
#log.p=T znacilo bi da vraca log(p) umesto pog.p=T znacilo bi da vraca 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 6 5 6 6 5 8 1 6 4 4 4 8 4 4 4 3 4 4 3

#crtanje verovatnoca:
#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 verovatnoce predstavimo stubicima:
plot(x,
dbinom(x, size = 50, prob = 0.4),
type = "h",
main = "") #h od histogram

#mozemo 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 tacki 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 nacin:
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)) #vraca bas 2, ovo su inverzne f-je

## [1] 2

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

##  [1]  0.916856930  0.161056830 -1.043932920  0.029444169 -1.108827620
##  [6]  0.495227985 -0.001098055  0.214943793 -1.275785227 -0.224672256
## [11] -0.098605205  1.020913184 -0.422764621 -0.624737345  0.555244089
## [16]  0.789973195  1.016657271 -0.672296333 -0.430541216 -0.253345658

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

## [1] 0

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.58695384 0.15599255 0.50377994 0.31202325 0.02344991 0.69467740
##  [7] 0.68462019 0.76073674 0.14427789 0.88482455

#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.30086

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

## [1] 0.59967

1 - punif(0.4)  # P{U>0.4}=1-P{U<=0.4} Tacna vrijednost!

## [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.50063

pnorm(0, sd = sqrt(10))

## [1] 0.5

# Z ~ Pois(3)

z <- rpois(n, lambda = 3)

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

## [1] 0.22298

dpois(2, lambda = 3) # P{Z=2} Tacna vrijednost!

## [1] 0.2240418

Generisanje apsolutno neprekidnih slucajnih velicina (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 generisemo slucajan broj U iz U[0,1]

#ovo tvrdjenje nam daje mogucnost da generisemo sl brojeve iz onih
#apsolutno neprekidnih raspodela za koje se inverz funkcije raspodele se moze izraziti analiticki

• Eksponencijalna raspodela

Exp_random <- function(lambda)
{
u <- runif(1)
x <- (-1 / lambda) * log(1 - u)
return(x)
}
# Mozemo da napravimo funkciju koja generise uzorak obima N iz eksponencijalne raspodjele.

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

##  [1] 6.1906379 2.3634118 1.0231841 1.3464902 0.3614683 0.1735427 2.4090706
##  [8] 0.5612781 9.1145657 4.5270305

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

Erlang_random <- function(n, lambda)
{
x <- rexp(n, lambda)
sum(x)
}
Erlang_random(10, 2)

## [1] 4.444395

#Jos neke apsolutno neprekidne raspodele koje se mogu modelirati ovom metodom su
#Vejbulova, Freseova i Gumbelova
#Pogledajte njihove funkcije raspodele i napisite funkcije za
#generisanje slucajnih brojeva iz tih raspodela.
#Ove tri raspodele su vazne jer su jedine tri nedegenerisane granicne raspodele
#normiranog (nekim konstantama) niza maksimuma nezavisnih slucajnih velicina.

Neki primjeri apsolutno neprekidnih raspodela

# Neka je X ~ Exp(a). Generisati uzorke iz sledecih raspodela i izracunati 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.6382001

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

vrednost = simulacija_b(10000, 2)

# srednja vrednost
mean(vrednost)

## [1] 0.3595116

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

# srednja vrednost
mean(vrednost)

## [1] 0.1584

Aproksimacije binomne raspodele

# Posmatrajmo eksperiment bacanja lopte u kos. Neka imamo ukupno 300 bacanja 
# i verovatnoca pogotka je 1/5.
# Izracunati verovatnocu da je bilo najvise 70 pogotka.

# Neka je X~Bin(300, 1/5). Trazi se P{X<=100}=
# 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 da u nekom gradu stanovnik ima bicikl sa verovatnocom 0.02,
# a motor sa verovatnocom 0.01, s tim sto niko nema i bicikl i motor. 
# Izracunati verovatnocu da od 100 slucajno odabranih stanovnika broj
# onih koji poseduju bar jedno od ova dva prevozna sredstva bude izmedju
# 2 i 6 ukljucujuci i te brojeve.

# Zapravo prava vrednost je da aproksimiramo 
# bas tu binomnu B~Bin(100, 0.03) sa Puasovom(3), koju cemu oznaciti 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 novcica n puta, pri cemu u svakom bacanju pismo pada sa 
# verovatnocom p, nezavisno od ostalih bacanja. Neka je Xi slucajna velicina 
# koja uzima vrednost 1 ako je u i-tom bacanju palo pismo, 0 inace, za sve 
# i=1,..n. Hocemo da proverimo (simulacijom) da li za ovakav niz slucajnih 
# velicina vazi slabi zakon velikih brojeva.

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

# ESn=nEX1=n*p, sto 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.0356

# Verovatnoca dosta varira u zavisnosti do izbora eps zbog nepreciznosti samih
# ocena, a kako znamo da su ove slucajne velicine nezavisne i imaju konacno ocekivanje, na osnovu teoreme sa casa znamo da za takav niz vazi zakon velikih brojeva.

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

CENTRALNA GRANIČNA TEOREMA

Neka je \(S_n=X_1+...+X_n\), gdje 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 raspodjele 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 domaci: Uraditi ovaj zadatak za razlicite 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.

# 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