# Generisanje slucajnih velicina


# Binomna raspodjela

# X~B(n,p)  broj uspjeha u n pokusaja


# Na proslom casu smo vidjeli da ako zelimo da generisemo neki uzorak iz binomne
# raspodje sa parametrima n i p, odnosno da generisemo slucajnu velicinu X iz 
# binomne B(n,p) raspodjele, za to imamo gotovu funkciju rbinom(velicina
# vektora, n, p).


# Hocemo da generisemo binomnu slucajnu velicinu bez koriscenja fje rbinom.

# 1. nacin:


# Objasnjenje: na casu



binomna1 <- 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
  
  
}

# 2. nacin

# Neka  X1, X2,...,Xn  imaju Bernulijevu raspodjelu Ber(p), tj. to su indikatori
# nekog dogadjaja cija je vjerovatnoca p. Tada ce zbir X=X1+X2+...+Xn imati
# binomnu B(n,p) raspodjelu.


# Kako nam za binomnu raspodjelu treba n indikatora, generisacemo vektor od n 0
# ili 1 koje ce da predstavljaju uspjehe odnosno neuspjehe. Zbir takvog vektora
# ce da da bas broj uspjeha u n nezavisnih izvodjenja eksperimenta, sto bas
# predstavlja binomna slucajna velicina



binomna2 <- function(n, p) {
  vec <- sample(c(1, 0),
                size = n,
                replace = TRUE,
                prob = c(p, 1 - p))
  
  k <- sum(vec)
  
  return(k)
  
}

# Sa proslog casa znamo da je ocekivanje X iz B(n,p) EX=np.

# Ocjenjivanje matematickog ocekivanja EX:

# Ocjena za ocekivanje neke slucajne velicine je srednja vrijednost zbira 
# nezavisnih slucajnih velicina iz iste raspodjele.

# Dakle, da bismo ga ocijenili, generismo vektor iz nase raspodjele duzine N i
# nadjemo njegovu srednju vrijednost.

# Mozemo da iskoristimo bilo koju od funkcija rbinom(), binomna1() ili
# binomna2() za generisanje takvog vektora


ocekivanje<-function(n,p){
  
  vec<-replicate(10000, binomna1(n,p))
  mean(vec)
  
}
ocekivanje(10,0.7) # 10*0.7=7


# Disperzija DX=np(1-p)
# Rekli smo da je ocjena za disperziju var(x)

disperzija<-function(n,p){
  
  vec<-replicate(10000, binomna1(n,p))
  var(vec)
  
}

disperzija(10,0.7) 
10*0.7*0.3


# Geometrijska raspodjela

# X~G(p) broj pokusaja do prvog uspjeha

# Na proslom casu smo vidjeli da postoji gotova funkcija rgeom() koja ce da 
# vrati uzorak iz geometrijske raspodjele, medjutim ne bas one koju smo mi 
# definisali, nego pomjerene, koja uzima vrijednosti od 0, tj. broji neuspjehe 
# do prvog uspjeha
# Mi cemo generisati bas slucajnu velicinu koja broji sve pokusaje do prvog 
# uspjeha.

# Algoritam za generisanje broja iz geometrijske raspodjele sa parametrom p.

# 1. Generisemo nezavisne Bernulijeve slucajne velicine X1,X2,...sa parametrom p
# (ogranicemo broj na 1000 zbog veoma male terojske vjerovatnoce da se prvi 
# uspjeh ne desi prije 1000. pokusaja)

# 2. Trazimo najmanji indeks i takav da je Xi uzela vrijednost 1, tj trazimo prvi
# elemenat u vektoru koji je jednak jedinici

# 3. Trazeni broj pokusaja je bas i.



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
  # Ako su svi nula, smatramo da nikad nije doslo do uspjeha
  
  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
}

# Teorijski izracunato, ocekivanje EX=1/p za X~G(p).

# Provjeravamo simulacijama, analogno kao za binomnu raspodjelu.


ocekivanje<-function(p){
  
  vec<-replicate(10000, geometrijska(p))
  mean(vec)
  
}
ocekivanje(0.4) 
1/0.4


# Disperzija DX=(1-p)/p^2


disperzija<-function(p){
  
  vec<-replicate(10000, geometrijska(p))
  var(vec)
  
}
disperzija(0.4)
(1-0.4)/0.4^2


# Napomena: Neka je X ~ G(p), a neka je Y pomjerena geometrijska sa istim 
# parametrom, odnosna ona kojoj odgovata fja rgeom(p). Tada Y mozemo da 
# predstavimo preko X kao: Y=X-1. Odatle, EY=E(X-1)=EX-1 -> EY= 1/p - 1
# Disperzija ce biti ista.

ocekivanje<-function(p){
  
  vec<-rgeom(1000,p)
  mean(vec)
  
}
ocekivanje(0.4) 
1/0.4-1


# Puasonova raspodjela 

# broj realizacija nekog dogadjaja na zadatom intervalu, npr. u jedincici
# vremena


# X ~ P(lambda), lambda >0

# X={0,1,2,3,...}

# Prefiksi r, p, d i q oznavaju iste funkcije kao do sad:

dpois(x, lambda) # P{X=x}

ppois(x, lambda) # F(x)=P{X<=x}-funckija raspodjele

rpois(n, lambda) # Vektor velicine n sa slucajnim uzorkom iz Puasonove raspodjele


# Ocekivanje teorijsko EX=lambda, kao i disperzija DX=lambda

# Mozemo da ga ocijenimo kao sva prethodna ocekivanja.

ocekivanje <- function(lambda) {
  vec <- rpois(1000, lambda)
  mean(vec)
  
}
ocekivanje(3)


# Takodje i disperziju:


disperzija <- function(lambda) {
  vec<- rpois(1000, lambda)
  var(vec)
  
}
disperzija(3)


# Primjer:

# Broj poziva za hotelske rezervacije je u prosjeku 3 poziva po satu. Izracunati
# vjerovatnocu u toku jednog sata stignu bar 2 poziva.

# X~P(3) Trazi se: P{X>=2}=1-P{X<2}=1-P{x<=1}=1-F(1)

#  resenje:


1 - ppois(1, lambda = 3)

# Primjer:

# Slucajna velicina X ~ P(3), dok Y~P(5) i nezavisne su. Ocijeniti ocekivanje
# slucajne velicine Z=X+Y.

# Hocemo da generisemo vektor velicine N iz raspodjele Z:

Z<-function(N){
  
  X<-rpois(N, lambda=3)
  Y<-rpois(N, lambda=5)
  
  X+Y
  
}

mean(Z(1000))