\[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] 0
ber2 <- function(p, N) {
replicate(N, ber(p))
}
# Ocjenjujemo srednju vrijednost
mean(ber2(0.8, 1000))
## [1] 0.798
ber3 <- function(p, N) {
sample(c(0, 1),
size = N,
replace = TRUE,
prob = c(1 - p, p))
}
mean(ber3(0.3, 100))
## [1] 0.25
\[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.969
# Ocjena vjerovatnoce P{X=3}
mean(r==4)
## [1] 0.035
# Poredjenje sa teorijskim rezulatima:
10*0.7
## [1] 7
choose(10,4)*0.7^4*(1-0.7)^(10-4)
## [1] 0.03675691
# 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 5 8 5 9 6 9 8 7 5 7 8 8 7 7 6 7 10 7 5 5 9
## [24] 8 8 9 7 9 7 7 8 8 6 7 6 6 8 7 5 3 8 7 5 5 5 7
## [47] 7 8 7 7 7 9 9 7 6 6 8 7 8 2 8 9 8 8 7 8 7 6 7
## [70] 8 9 8 7 8 6 6 5 5 7 7 7 4 6 8 6 8 8 8 9 8 6 9
## [93] 5 7 5 7 7 6 6 8
# 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.
\[X \in \mathcal{G}(p), \quad P\{X=k\}=(1-p)^{k-1}p,\quad k=1,2,\dots \]
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.022
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] 1.96
1 - pgeom(4, 0.7)
## [1] 0.00243
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
sum(dgeom(seq(from = 1, to = 1000, by = 2),prob = 0.7))
## [1] 0.2307692
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
\[ 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
}
\[ 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<-function(n,lambda){
r<-replicate(n,eksponencijalna(lambda))
sum(r)
}
Koristeći samo generator uniformne raspodjele generisati uzorke iz sledećih raspodjela
\[ \left(\begin{matrix}
-1 & 0 & 1 & 2 \\
p/3 & p/6 & p & (1-3/2*p)
\end{matrix}\right)
\] za odgovarajuće vrijednosti parametra p.
b) [Vejbulova raspodjela] \[ F(x)=1-e^{-(x/\alpha)^\beta}, \ x \geq 0, \ \alpha>0, \ \beta >0.
\]
Parametri raspodjele treba da budu zadati kao argumenti funkcija. Za odabrane parametre (konkretne vrijednosti po izboru) generisati uzorak iz date raspodjele i nacrtati histogram. Naći ocjenu srednje vrijednosti i uporediti je sa teorijskom. Kod aposlutno neprekidnih raspodjela dodati na histogram teorijsku krivu gustine.
(Očekivanje Vejbulove raspodjele je \(\alpha \Gamma (1+1/\beta)\).)
.