###########
## CAS 6 ##
###########

# Ispricati o Indikatoru.
# Indikator ili Bernulijeva raspodela u oznaci Ber(p), gde paramerat p - verovatnoca da se ostvari dogadjaj.
# Neka je X~ Ber(p), tada vazi:
# P(X=0) = 1-p
# P(X=1) = p.
# E(X) = p
# D(X) = p*(1-p)

# Neka su X_1, X_2,..., X_n ~ Ber(p) i nezavisne. Tada vazi X_1+...+X_n = Y ~ Bin(n, p).
# Odakle mozemo lako da uocimo da zapravo i zbir dve nezavisne slucajne velicine koje
# obe imaju Binomnu raspodelu sa istim parametrom p, 
# je takodje Binomna sa istim parametrom p, dok se odgovarajuci parametri "n" sabiraju.

# 1) 
# (a) Simulirati binomnu slucajnu velicinu sa parametrima "n" i "p", ako je dat generator slucajnih brojeva na (0,1). 
# Naci matematicko ocekivanje i disperziju te slucajne velicine.

# (b) Simulirati geometrijsku slucajnu velicinu sa parametrom "p", ako je dat generator slucajnih brojeva na (0,1). 
# Naci matematicko ocekivanje i disperziju te slucajne velicine.


# Resenje:

# (a)
binomna <- function(n, p)
{
  # generisemo n brojeva izmedju 0 i 1 i gledamo koliko njih je upalo u [0, p]
  brojevi = runif(n)
  sum(brojevi <= p)
}
binomna(5, 1/2)

# Matematicko ocekivanje binomne raspodele sa parametrima n i p je E(X) = n*p

mean(replicate(10000, binomna(5, 1/2))) # ocekivano i vraca oko 2.5

# Dok je disperzija n*p*(1-p)

var(replicate(10000, binomna(5, 1/2))) # 5*1/4 = 1.25


# Sve je to zapravo vec ugradjeno u R-u i za to sluze funkcija:
rbinom(n = koliko_zelimo_slucajnih_brojeva, size =  parametar_n, prob = parametar_p)

# (b)
geometrijska <- function(p, maks_broj)
{
  # generisemo niz brojeva iz [0, 1]
  brojevi = runif(maks_broj)
  # i gledamo kada je prvi put verovatnoca p veca od niza generisanih brojeva,
  # jer which od tog izraza vraca prvu poziciju gde se desio FALSE.
  realizacije = which((p <= brojevi) == FALSE)[1]
  # Ako nikada nije, onda vratimo 1
  if(length(realizacije) == 0 ) return(1)
  return(realizacije)
}

geometrijska(1/8, 1000)
mean(replicate(100000, geometrijska(1/8, 1000) ) )

var(replicate(100000, geometrijska(1/8, 1000) ) )


# Matematicko ocekivanje geometrijske raspodele sa parametrom p je E(X) = 1/p, ako uzima vrednosti od 1,
# i (1-p)/p ako uzima od 0.

mean(replicate(10000, geometrijska(1/2))) 

# Dok je disperzija (1-p)/p^2 u oba slucaja. Zasto?

var(replicate(10000, geometrijska(1/2)))


# Sve je to zapravo vec ugradjeno u R-u i za to sluze funkcija:
rgeom(n = koliko_zelimo_slucajnih_brojeva, prob = parametar_p)
mean(rgeom(10000, 1/8) )+1


### Napomenuti funkcije: rbinom() i rgeom().


# Puasonova raspodela sa parametrom a>0 u oznaci P(a)
# P(X = k) = a^k*exp(-a)/k!
# E(X) = a
# D(X) = a

# generator pseudoslucajnih brojeva iz Puasonove raspodele:
rpois(n = 10 , lambda = 3)
# Funkcija raspodele: P(X<=x) = ppois(x, a)
a = 2; x = 2
ppois(q = x, lambda = a)
# Vrednost u tacki: P(X=x) = dpois(x, a)
a = 2; x = 2
dpois(x = x , lambda = a)


# 2) Slucajna velicina X ima Puasonovu P(a) raspodelu. 
# Izracunati ocekivanje i disperziju te slucajne velicine.

mean(rpois(n = 10000, lambda = 2))
var(rpois(n = 10000, lambda = 2))


# 3) Slucajna velicina X ima Puasonovu P(a) raspodelu, dok Y ima Puasonovu P(b).
# X i Y su nezavisne. Ako je Z = X + Y, naci raspodelu slucajne velicine Z i njeno ocekivanje.

simulacija <- function(a,b)
{
  pois_a = rpois(1, a)
  pois_b = rpois(1, b)
  
  raspodela = pois_a+pois_b
  return(raspodela)
}

# Generisanje broja iz Z
simulacija(2, 4)

# Ocena matematickog ocekivanja za Z.
mean(replicate(10000, simulacija(2, 4)))

# Teorema: zbir nezavisnih slucajnih velicina koje imaju Puasonovu raspodelu je takodje Puasonova raspodela, sa parametrom
# koji predstavlja zbir parametara pocetnih Puasonovih raspodela.

# 4) Baca se kockica za igru. Odrediti ocekivani broj bacanja do pojave svih brojeva.
kockica_svi_brojevi <- function()
{
  x = 6 + rgeom(1, 5/6) + rgeom(1, 4/6) + rgeom(1, 3/6) + rgeom(1, 2/6) + rgeom(1, 1/6) 
  return(x)
}

mean(replicate(100000, kockica_svi_brojevi())) 
# oko 14.7 treba da bude


# 5) Iz kutije u kojoj su cetiri cedulje numerisane brojevima {1, 2, 3, 4} izvlacimo,
# bez vracanja, dok ne izvucemo cedulju sa neparnim brojem.
# Ako je X - zbir izvucenih brojeva, a Y -  broj izvlacenja. 
# Izracunati ocekivanje i disperziju od X i Y.

# Resenje:
izvlacenja <- function()
{
  redosled = sample(1:4, 3, replace = FALSE)
  suma = cumsum(redosled)
  
  y = which(suma %% 2 == 1)[1] # to ce zapravo biti Y.
  
  x = suma[y]
  return(cbind(x, y))
}
izvlacenja()

rowMeans(replicate(10000, izvlacenja(), simplify = TRUE))

rezultat = replicate(10000, izvlacenja(), simplify = TRUE)

disperzija_X = var(rezultat[1, ])
disperzija_X

disperzija_y = var(rezultat[2, ])
disperzija_y

### Zakon velikih brojeva

# Cebisovljeva nejednakost:
# Teorema:
# Neka je X diskretna slucajna velicina sa ocekivanom vrednoscu 
# E(X)=m i epsilon>0. 
# Tada:
# P(|X-m|>=epsilon) <= D(X)/epsilon^2

# Zakon velikih brojeva:
# Teorema:
# Neka su X_1,..., X_n nezavisne i jednakoraspodeljene slucajne velicine 
# sa konacnim ocekivanjem E(X_i)=m i disperzijom D(X) = s^2.
# Neka je S_n = X_1+...+X_n. 
# Tada za svako eps>0 vazi
# P(|S_n/n -m|>=eps) -> 0, kada n -> infty.

# 6)
# Neka je X takva slucajna velicina da vazi E(X)=m, i D(X)=s^2. 
# Naci takav broj k da vazi:
# P{|X - m|>=k} <= 0.01.

# I pusticemo simulaciju za neku takvu slucajnu velicinu. 
# Neka je to X~Bin(n, p), td np=m, a np(1-p)=s^2
zadatak6 <- function(k,  N, n, p)
{
  br = 0 # broji slucajeve kada je |X - m| >= k
  vektor = rbinom(N, n, p) # vektor slucajnih brojeva iz trazene raspodele duzine N.
  br = sum(ifelse(abs(vektor - n*p) >= k, 1, 0) )
  return(br/N)
}
# za n = 100, p = 0.3
zadatak6(10, 10000, 100, 0.3) # mora manje k
zadatak6(30, 10000, 100, 0.3) # moze i vece k

zadatak6(13, 10000, 100, 0.3) # minimalno takvo k. 

# Dok teorema nam kaze
# 30*0.7/eps^2 = 0.01
# Odnosno:
eps = sqrt(30*0.7/0.01)
eps

zadatak6(eps, 10000, 100, 0.3)