Diskretne slučajne veličine

Metoda inverzne transformacije u diskretnom slučaju

Neka je \(X\) slučajna veličina sa vjerovatnoćama

\(P(X=x_i)=p_i, \quad i=0,1,... \quad \sum\limits_{i=0}^{\infty}p_i=1.\)

Da bi smo generisali jednu realizaciju slučajne veličine \(X\), prvo generišemo slučajan broj U iz intervala \((0,1)\) i potom kažemo da je \(X\) uzela vrijednost \(x_i\) ako

\[ \sum\limits_{j=0}^{i-1}p_j \leq U < \sum\limits_{j=0}^{i}p_j \]

Dakle, algoritam je sledeći:

-Generišemo slučajan broj U iz \(U(0,1)\) raspodjele

-Nađemo indeks \(I\) za koji važi

\[ \sum\limits_{j=0}^{I-1}p_j \leq U < \sum\limits_{j=0}^{I}p_j \]

Bernulijeva slučajna veličina

\[X \in Ber(p)\]

ber <- function(p) {
  # Generišemo slučajan broj iz (0,1)
  
  U <- runif(1)
  
  # Vratimo slučajan broj iz Ber(p)
  
  ifelse(U < p, 1, 0)
  
}

ber(0.5)
## [1] 1

Ako želimo da generešimo vektor sa N vrijednosti iz \(Ber(p)\) raspodjele, iskoristićemo već napravljenu funkciju ber() i ponoviti je N puta:

ber2 <- function(p, N) {
  replicate(N, ber(p))
  
  
}
# Ocjenjujemo srednju vrijednost
mean(ber2(0.8, 1000))
## [1] 0.782

Drugi način:

ber3 <- function(p, N) {
  sample(c(0, 1),
         size = N,
         replace = TRUE,
         prob = c(1 - p, p))
}
mean(ber3(0.3, 100))
## [1] 0.27

Binomna slučajna veličina

\[X \in \mathcal{B}(n,p)\]

binomna <- function(n, p) {
  k <- (-1):n
  vec <- choose(n, k) * p ^ k * (1 - p) ^ (n - k)
  
  x <- runif(1)
  
  which((cumsum(vec)[-length(vec)] < x) & (x <= cumsum(vec)[-1])) - 1
  
  
}
# Srednja vrijednost
r<-replicate(1000,binomna(10,0.7))
mean(r)
## [1] 6.932
# Ocjena vjerovatnoce P{X=3}
mean(r==4)
## [1] 0.042
# Poredjenje sa teorijskim rezulatima:
10*0.7
## [1] 7
choose(10,4)*0.7^4*(1-0.7)^(10-4)
## [1] 0.03675691

Ako iskoristimo poznatu činjenicu da zbir N nezavisnih slučajnih veličina sa Ber(p) raspodjelom (N nezavisnih indikatora) ima binomnu B(n,p) raspodjelu, možemo da generišemo vektor iz binomne raspodjele na jednostavniji način tako što iskoristimo neku od funkcija pomenutih gore.

binomna2 <- function(n, p) {
  vec <- sample(c(1, 0),
                size = n,
                replace = TRUE,
                prob = c(p, 1 - p))
  
  k <- sum(vec)
  
  return(k)
  
}
binomna2(10,0.7)
## [1] 7
mean(replicate(1000,binomna2(10,0.7)))
## [1] 7.059
Za vježbu: Napraviti funkciju koja ocjenjuje disperziju slučajne veličine X iz B(n,p) raspodjele.

U R-u postoji funkcija rbinom() koja generiše vektor iz B(n,p) raspodjele.

# Prvi argument je velicina vektora koju zelimo da dobijemo, dok su drugi i treci parametri binomne raspidjele n i p, redom.
rbinom(100, size = 10, prob = 0.7)
##   [1]  5  6  7  7  8  5  7  9  6  4  7  7  8  6  6  7  7  6  6  7  8  8  9
##  [24]  5  7  9  8  6  6  5  7  7  6  7  8  7  6  5  8  8  6  7  7  7  8  8
##  [47]  9  8  8  6  7  7  8  9  6  5  6  7  6  8  6  7  6  7 10  6  5  7  3
##  [70]  9  8  6  7  6  6  6  7  7  8  8  6  7  7  8  7  5 10  6  6  7  5  8
##  [93]  7  9 10  7  7  8  9  7
# funkcija raspodjele dobija se pozivom pbinom()
# npr F(2) za X iz B(10,0.7):

pbinom(2, 10, 0.7)
## [1] 0.001590386
# vjerovatnoce iz zakona raspodjele dobijaju se fjom dbinom()

# P{X=4}
dbinom(4, 10, 0.7)
## [1] 0.03675691
# Napomena: fje pbinom i dbinom daju prave (teorijske vjerovatnoce), a ne ocijenjene.

Geometrijska slučajna veličina

\[X \in \mathcal{G}(p), \quad P\{X=k\}=(1-p)^{k-1}p,\quad k=1,2,\dots \]

Simuliraćemo vrijednosti nezavisnih indikatora, kod kojih se uspjeh javlja sa vjerovatnoćom p. Potom gledamo kod kojeg indikatora se prvi put pojavila jedinica (odnosno uspijeh) i vratimo njegov redni broj. Pošto geometrijska slučajna veličina uzima beskonačno mnogo vrijednosti, mi taj broj moramo ograničiti pa ćemo gledati 1000 indikatora i ako su svi bili 0, proglasićemo da se nikad nije dogodio uspjeh, tj. smatramo da je broj pokušaja do prvog uspjeha beskonačan. Dakle, sve ishode koji su veći od 1000 proglasimo za ishod beskonačno koji se dešava sa vjerovatnoćom \(P\{X > 1000\}\).

geometrijska <- function(p) {
  vec <- sample(c(1, 0),
                size = 1000,
                replace = TRUE,
                prob = c(p, 1 - p))
  
  # Trazimo prvi element iz vektora koji je jednak jedinici
  
  
  if (all(vec == 0))
    return("infty")
  
  which.max(vec == 1)
  
  # Napomena: Ako je x logicki vektor (sto postizemo zadavanjem nekog uslova u 
  # zagradi)which.min(x) i which.max(x) vracaju indeks prvog FALSE ili TRUE, 
  # redom, posmatrajuci FALSE<TRUE. Kako mi trazimo najmanji indeks takav da nam
  # je uslov ispunjen, koristicemo which.max
}
# Primjer: 
geometrijska(0.5)
## [1] 2
mean(replicate(1000,geometrijska(0.5)))
## [1] 2.01

Drugi način (inverzna transformacija)

geometrijska2<-function(p){
  
  U<-runif(1)
  k<-1:1000
  # vektor vjerovatnoca iz zakona raspodjele
  prob<-c(0,(1-p)^(k-1)*p)
  which((cumsum(prob)[-length(prob)] < U) & (U <= cumsum(prob)[-1])) 
  
  
}
mean(replicate(1000,geometrijska2(0.5)))
## [1] 2.002
Funkcija rgeom() će generisati vektor iz geometrijske raspodjele, koja se malo razlikuje od prethodne definicije. Tačnije, to će biti geometrijska slučajna veličina koja broji neuspjehe do prvog uspjeha, tj. uzima vrijednosti 0,1,2… Dakle, ukoliko koristite ovu funkciju potrebno je korigovati vrijednosti za 1. Isto kao kod binomne, postoje funkcije pgeom() i dgeom().

Primjer: Vjerovatnoća da košarkaš pogodi koš je p=0.7. On gađa sve dok ne pogodi. Izračunati:

a) vjerovatnoću da košarkaš ima više od 5 pokušaja

b) vjerovatnoću da gađa između 3 i 7 puta

c) vjerovatnoću da je bio paran broj pokušaja

d) srednji broj pokušaja

a) \(P\{X>5\}=P\{X-1>5-1\}=P\{X-1>4\}=1-P\{X-1\leq 4\}\)

1 - pgeom(4, 0.7)
## [1] 0.00243

b) \(P\{3\leq X \leq 7\}=P\{X \leq 7\}-P\{X<3\}=P\{X-1\leq 6\}-P\{X-1\leq 1\}\)

pgeom(6,0.7)-pgeom(1,0.7)
## [1] 0.0897813
# ili P{3<=X<=7}=P{X-1=3-1}+P{X-1=4-1}+...+P{X-1=7-1}

sum(dgeom(2:6,0.7))
## [1] 0.0897813
c) Napomena: ovdje nam se javlja beskonačna suma pa cemo da nametnemo neku granicu. Npr. pgeom(1000,0.7) = 1 sto znači da je vjerovatnoća da X uzme vrijednost veću od 1000 0 (to se dešava zbog zaokruživanja, teorijski te vjerovatnoće su bliske nuli), pa za granicu uzimamo 1000.
Zbog toga sto gledamo fje za X-1, zapravo racunamo vjerovatnocu da bude neparan broj pokusaja
sum(dgeom(seq(from = 1, to = 1000, by = 2),prob = 0.7))
## [1] 0.2307692

d) Slicno kao u prethodnom primjeru, ogranicicemo broj vrijednosti na 1000

f <- function(k, p=0.7) {
  pk <- (1 - p) ^ (k - 1) * p
  return(k * pk)
  
}


ocekivanje <- sum(sapply(1:1000 , f))
ocekivanje 
## [1] 1.428571
Za vježbu: Baca se regularna kockica za igru. Ocijeniti očekivani broj do pojave svih brojeva.

Apsolutno neprekidne slučajne veličine

Metoda inverzne transformacije

Neka je X slučajna veličina sa neprekidnom i rastućom funkcijom raspodjele \(F\). Označimo sa \(F^{-1}\) inverz funkcije \(F\). Tada slučajnu veličinu X generišemo na sledeći način:

-Generišemo slučajan broj U iz \(U(0,1)\) raspodjele

-Vratimo vrijednost \(F^{-1}(U)\)

Na ovaj način dobijamo vrijednosti iz raspodjele X jer pri ovim uslovima važi \(F(X) \in U(0,1)\), odosno \(F^{-1}(U)\) ima istu raspodjelu kao X.

Ako F nije neprekidna ili rastuća, koristimo uopšteni inverz (ovo objašnjava zašto ovaj metod može da se primijeni i u diskretnom slučaju):

\[F^{-1}(u)=inf\{x:F(x)\geq u\}\]

Uniformna raspodjela na proizvoljnom intervalu (a,b)

\[ U \in \mathcal{U}(0,1) \rightarrow a+(b-a)U \in \mathcal{U}(a,b) \]

unif<-function(a,b){
U<-runif(1)
a+(b-a)*U
}

Eksponencijalna raspodjela

\[ U \in \mathcal{U}(0,1) \rightarrow -\frac{1}{\lambda}\log(1-U) \in \mathcal{E}(\lambda),\quad \lambda>0 \]

eksponencijalna<-function(lambda){

U <- runif(1)
-1/lambda * log(1-U)
# ili -1/lambda * log(U)
}

Gama raspodjela kao zbir nezavisnih eksponencijalnih

gama<-function(n,lambda){

r<-replicate(n,eksponencijalna(lambda))
sum(r)

}