Instalacija R/Rstudio

Prvo je potrebno instalirati programski jezik R, a onda Rstudio IDE. Uputstvo sa slicicama se moze videti na ovom linku.

Osnovni koraci

Deklaracija promjenljive

Operator dodjele u R-u je <- (ali moze da se koristi i =)

x <- 9

Za pokretanje linije koda (ako koristite RStudio): Ctrl+Enter
Postoji i operator -> za dodjelu sa druge strane

9 -> x

Numericke vrijednosti su podrazumijevano tipa double ( na primer 5/2 nije cjelobrojno dijeljenje već je jednako 2.5 )

Elementarne funkcije: exp(), cos(), sin(), log(), ..

Aritmetički operatori: +,-,*,/

Logički operatori: <=,>=, ==, !=, <,>, & , |, !

Vektori

U R-u postoji mnoštvo funkcija za rad sa vektorima i poželjno da se podaci čuvaju u formi vektora (ili matrica)

Neki od načina zadavanja vektora:

numeric.vec <- c(1, 10, 49)
character.vec <- c("a", "b", "c")
boolean.vec <- c(TRUE, FALSE)

vec <- rep(0, 10)
vec <- rep(c(1, 2, 3), 3)

vec <- seq(1, 5)
vec <- seq(0, 1, 0.1)
vec <- seq(9, 1, -1)
vec <- seq(1, 5, length.out = 4)

vec <- 1:9

for petlja

vec <- numeric(10)

for (i in 1:10) {
vec[i] <- i ^ 2

}

Operacije sa vektorima

a <- 1:5
b <- c(3, 4, 5, 6, 7)

Aritmetičke operacije

5 * a
## [1]  5 10 15 20 25
a + b
## [1]  4  6  8 10 12
a - b
## [1] -2 -2 -2 -2 -2
a * b
## [1]  3  8 15 24 35
a / b
## [1] 0.3333333 0.5000000 0.6000000 0.6666667 0.7142857

Logičke operacije

a == b
## [1] FALSE FALSE FALSE FALSE FALSE
a > b
## [1] FALSE FALSE FALSE FALSE FALSE
a <= 5
## [1] TRUE TRUE TRUE TRUE TRUE

Ako hoćemo da poredimo cijela dva vektora

all(a == b)
## [1] FALSE
all(a == a)
## [1] TRUE

Ako dva vektora nisu iste dužine, vrijednosti kraćeg vektora se ponavljaju ciklično dok se ne pogode iste dimenzije sa dužim vektorom

u<-c(10, 20, 30) 
v<-c(1, 2, 3, 4, 5, 6, 7, 8, 9) 
u + v 
## [1] 11 22 33 14 25 36 17 28 39

I funkcije u R-u su vektorske (primjenjuje se na svaki član vektora)

cos(a)
## [1]  0.5403023 -0.4161468 -0.9899925 -0.6536436  0.2836622
log(a)
## [1] 0.0000000 0.6931472 1.0986123 1.3862944 1.6094379

Pristupanje elementima vektora

Indeksiranje vektora počinje od jedinice!

v[3]  # element na 3. poziciji u vektoru v
## [1] 3

Možemo da izdvojimo i neke cijele vektore

v[c(1, 2)]
## [1] 1 2
v[c(3, 4, 5)]
## [1] 3 4 5
v[2:4]
## [1] 2 3 4
v[c(2, 1, 3)]
## [1] 2 1 3
v[v < 5]
## [1] 1 2 3 4
v[v != 8]
## [1] 1 2 3 4 5 6 7 9

Ako stavimo negativan predznak indeksu rezultat je originalni vektor umanjen za clan na poziciji koja odgovara apsolutnoj vrijednosti broja u zagradi

v[-4]
## [1] 1 2 3 5 6 7 8 9
v[-c(5, 7)]
## [1] 1 2 3 4 6 8 9
v[-(1:3)]  # Napomena: ":" ima manji prioritet od aritmetickih operacija 
## [1] 4 5 6 7 8 9

Novi vektor može se izdvojiti iz datog vektora uz pomoć vektora sa logičkim vrijednostima koji je iste dužine kao originalni. Njegovi elementi su TRUE (ili T) ako odgovarajući element u orignalnom vektoru treba da bude izdvojen, odnosno FALSE (ili F) u suprotnom.

s <- c("a", "b", "c", "d")
l <- c(T, F, T, F)
s[l]
## [1] "a" "c"

Možemo da dodijelimo imena članovima vektora

v <- c("Ivo", "Andric")
names(v) <- c("Ime", "Prezime")
v
##      Ime  Prezime 
##    "Ivo" "Andric"

Zatim mozemo da pristupimo elementima po imenu

v["Ime"]
##   Ime 
## "Ivo"

Neke funkcije nad vektorima

vec <- rev(vec) # obrce redosled clanova vektora
length(vec)
## [1] 10

Važne funkcije nad vektorima

sum(vec)   # zbir 
## [1] 385
min(vec)  
## [1] 1
max(vec)   
## [1] 100
mean(vec)  # aritmeticka sredina
## [1] 38.5

Sortiranje

sort(vec)
##  [1]   1   4   9  16  25  36  49  64  81 100
sort(vec, decreasing = TRUE)
##  [1] 100  81  64  49  36  25  16   9   4   1

Izmjene u vektoru

vec[2] <- 0
vec[vec < 5] <- 1
vec <- vec[1:4]

Matrice

M <- matrix(c(2, 4, 3,
              1, 5, 7),
            nrow = 2,
            byrow = TRUE)   # popunjavamo matrice po redovime (po default-u je po kolonama)
M[2, 3]
## [1] 7
M[2, ]       # drugi red
## [1] 1 5 7
M[, 3]       # treca kolona
## [1] 3 7
M[, c(1, 2)] 
##      [,1] [,2]
## [1,]    2    4
## [2,]    1    5

Funkcije

Sintaksa za pisanje funkcija

# func.name<- function (argument) {
  
#         statement
# }

function je rezervisana riječ za deklaraciju funkcije.
Izrazi u okviru vitičastih zagrada čine tijelo funkcije (zagrade su opcione ako tijelo sadrži samo jedan izraz.
Konačno, funkcija se poziva sa func.name(argumenti).

Primjer

pow <- function(x, y) {
  
  result <- x^y
  
  result  #  return(result)
  
}

Pozivanje funkcije

pow(2, 3)
## [1] 8

Ili..

pow(x = 2, y = 3)
## [1] 8
pow(y = 3, x = 2)  #nije bitan redosled argumenata u ovom slucaju
## [1] 8

if/else

x <- 5
if (x %% 2 == 0) {
print("x je paran")
} else
print("x je neparan")
## [1] "x je neparan"
# ili

ifelse( x %% 2 == 0, print("paran"), print("neparan"))
## [1] "neparan"
## [1] "neparan"

Slučajnost u R-u (pseudo-slučajnost)

Funkcija sample()
U svom najjednostavnijem obliku, funkcija sample permutuje elemente zadatog vektora na slučajan način

x <- 1:10
sample(x)
##  [1]  1  3  2  5  9  4  6  8  7 10
sample(10)  # je poziv ekvivalentan prethodnom 
##  [1]  3  4  1  7 10  2  6  8  9  5

Ova funkcija ima argument replace čija je vrijednost podrazumijevano FALSE, što znači da se simulira slučajno biranje elemenata nekog skupa bez ponavljanja.
Ako zadamo parametar replace=TRUE funkcija vraca elemente odabrane na slučajan način sa ponavljanjem (svaki put se biraju iz cijelog skupa koji je zadat)
Još jedan argument je prob (odnosno probability) i on je podrazumijevano NULL što odgovara zakonu diskretne ravnomjerne raspodjele, a možemo da zadamo neke težine (vjerovatnoće) elementima skupa.
Argument size određuje koliko elemenata biramo iz proslijeđenog vektora. Napomena: Zbog gore pomenuta dva načina zadavanja moze da dođe do zabune prilikom sledećih izraza

sample(x[x > 8])  #izlaz: vektor duzine 2
## [1]  9 10
sample(x[x > 9])  #izlaz: vektor duzine 10
##  [1]  8  7  3 10  4  2  6  5  9  1

Hocemo da simuliramo bacanje novcica

coin <- c("Heads", "Tails")
sample(coin, size = 1)
## [1] "Tails"

Sada želimo 100 simulacija ovog slučajnog eksperimenta

#sample(coin, size = 100)

Javlja grešku. Treba da promijenimo vrijednost argumenta replace

sample(coin, size = 1000, replace = TRUE)
##    [1] "Heads" "Tails" "Tails" "Heads" "Tails" "Tails" "Tails" "Tails" "Heads"
##   [10] "Heads" "Heads" "Heads" "Tails" "Heads" "Tails" "Tails" "Tails" "Heads"
##   [19] "Heads" "Heads" "Heads" "Heads" "Tails" "Tails" "Heads" "Tails" "Tails"
##   [28] "Heads" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails" "Heads" "Heads"
##   [37] "Heads" "Tails" "Tails" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails"
##   [46] "Tails" "Heads" "Tails" "Tails" "Heads" "Heads" "Heads" "Tails" "Tails"
##   [55] "Heads" "Heads" "Heads" "Tails" "Heads" "Heads" "Heads" "Tails" "Heads"
##   [64] "Tails" "Tails" "Tails" "Heads" "Heads" "Tails" "Tails" "Heads" "Tails"
##   [73] "Heads" "Heads" "Heads" "Tails" "Tails" "Heads" "Heads" "Heads" "Tails"
##   [82] "Heads" "Tails" "Tails" "Tails" "Heads" "Heads" "Tails" "Heads" "Heads"
##   [91] "Tails" "Tails" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails" "Tails"
##  [100] "Heads" "Heads" "Tails" "Heads" "Tails" "Tails" "Heads" "Tails" "Tails"
##  [109] "Heads" "Heads" "Heads" "Tails" "Heads" "Heads" "Tails" "Tails" "Tails"
##  [118] "Tails" "Heads" "Heads" "Tails" "Heads" "Tails" "Tails" "Heads" "Heads"
##  [127] "Tails" "Heads" "Heads" "Heads" "Heads" "Tails" "Tails" "Tails" "Heads"
##  [136] "Heads" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails" "Tails" "Heads"
##  [145] "Tails" "Tails" "Heads" "Heads" "Heads" "Heads" "Tails" "Heads" "Heads"
##  [154] "Heads" "Tails" "Heads" "Heads" "Heads" "Heads" "Heads" "Heads" "Heads"
##  [163] "Heads" "Heads" "Heads" "Heads" "Heads" "Tails" "Heads" "Heads" "Heads"
##  [172] "Tails" "Heads" "Tails" "Tails" "Heads" "Heads" "Tails" "Heads" "Heads"
##  [181] "Heads" "Tails" "Tails" "Tails" "Heads" "Tails" "Tails" "Tails" "Heads"
##  [190] "Tails" "Heads" "Tails" "Heads" "Tails" "Tails" "Tails" "Heads" "Heads"
##  [199] "Heads" "Heads" "Tails" "Heads" "Heads" "Heads" "Heads" "Heads" "Heads"
##  [208] "Tails" "Heads" "Tails" "Tails" "Heads" "Heads" "Heads" "Heads" "Heads"
##  [217] "Tails" "Heads" "Heads" "Heads" "Heads" "Heads" "Heads" "Heads" "Heads"
##  [226] "Tails" "Heads" "Tails" "Tails" "Heads" "Tails" "Heads" "Heads" "Heads"
##  [235] "Tails" "Heads" "Tails" "Tails" "Tails" "Heads" "Tails" "Tails" "Heads"
##  [244] "Heads" "Tails" "Heads" "Heads" "Tails" "Tails" "Tails" "Tails" "Heads"
##  [253] "Tails" "Tails" "Tails" "Tails" "Heads" "Tails" "Tails" "Tails" "Heads"
##  [262] "Heads" "Heads" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails" "Heads"
##  [271] "Tails" "Tails" "Heads" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails"
##  [280] "Heads" "Tails" "Tails" "Tails" "Tails" "Heads" "Tails" "Tails" "Tails"
##  [289] "Tails" "Tails" "Tails" "Tails" "Heads" "Heads" "Heads" "Heads" "Tails"
##  [298] "Tails" "Tails" "Tails" "Tails" "Heads" "Heads" "Tails" "Heads" "Heads"
##  [307] "Tails" "Heads" "Heads" "Tails" "Tails" "Heads" "Heads" "Heads" "Tails"
##  [316] "Heads" "Heads" "Heads" "Tails" "Heads" "Tails" "Heads" "Tails" "Tails"
##  [325] "Heads" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails" "Heads" "Tails"
##  [334] "Tails" "Tails" "Heads" "Tails" "Heads" "Tails" "Heads" "Tails" "Tails"
##  [343] "Tails" "Heads" "Tails" "Heads" "Tails" "Tails" "Tails" "Tails" "Heads"
##  [352] "Tails" "Heads" "Tails" "Tails" "Heads" "Tails" "Tails" "Heads" "Heads"
##  [361] "Heads" "Heads" "Tails" "Tails" "Heads" "Tails" "Heads" "Tails" "Heads"
##  [370] "Tails" "Tails" "Tails" "Tails" "Tails" "Heads" "Heads" "Tails" "Tails"
##  [379] "Tails" "Tails" "Heads" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails"
##  [388] "Tails" "Tails" "Tails" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails"
##  [397] "Tails" "Heads" "Heads" "Heads" "Tails" "Heads" "Heads" "Heads" "Tails"
##  [406] "Tails" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails"
##  [415] "Tails" "Tails" "Heads" "Heads" "Heads" "Heads" "Tails" "Tails" "Tails"
##  [424] "Heads" "Heads" "Tails" "Tails" "Heads" "Tails" "Tails" "Tails" "Tails"
##  [433] "Tails" "Tails" "Tails" "Tails" "Heads" "Heads" "Tails" "Tails" "Heads"
##  [442] "Tails" "Tails" "Heads" "Heads" "Heads" "Tails" "Tails" "Heads" "Heads"
##  [451] "Heads" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails"
##  [460] "Tails" "Heads" "Tails" "Tails" "Tails" "Tails" "Heads" "Heads" "Tails"
##  [469] "Tails" "Tails" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails" "Heads"
##  [478] "Heads" "Tails" "Tails" "Heads" "Heads" "Heads" "Tails" "Tails" "Heads"
##  [487] "Heads" "Tails" "Tails" "Heads" "Heads" "Heads" "Heads" "Heads" "Tails"
##  [496] "Heads" "Heads" "Heads" "Tails" "Heads" "Heads" "Heads" "Tails" "Heads"
##  [505] "Tails" "Heads" "Heads" "Heads" "Tails" "Tails" "Tails" "Heads" "Heads"
##  [514] "Tails" "Heads" "Tails" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails"
##  [523] "Heads" "Heads" "Heads" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails"
##  [532] "Heads" "Tails" "Heads" "Tails" "Tails" "Tails" "Tails" "Heads" "Tails"
##  [541] "Heads" "Tails" "Tails" "Tails" "Heads" "Heads" "Tails" "Tails" "Tails"
##  [550] "Tails" "Tails" "Heads" "Tails" "Tails" "Tails" "Heads" "Tails" "Heads"
##  [559] "Tails" "Heads" "Heads" "Heads" "Heads" "Tails" "Heads" "Tails" "Heads"
##  [568] "Heads" "Heads" "Tails" "Heads" "Heads" "Tails" "Heads" "Heads" "Tails"
##  [577] "Tails" "Tails" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails" "Tails"
##  [586] "Heads" "Heads" "Tails" "Heads" "Heads" "Tails" "Tails" "Tails" "Heads"
##  [595] "Tails" "Heads" "Heads" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails"
##  [604] "Heads" "Tails" "Heads" "Tails" "Tails" "Tails" "Heads" "Heads" "Heads"
##  [613] "Tails" "Tails" "Tails" "Heads" "Tails" "Tails" "Heads" "Tails" "Tails"
##  [622] "Heads" "Heads" "Heads" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails"
##  [631] "Tails" "Heads" "Tails" "Heads" "Tails" "Tails" "Tails" "Tails" "Tails"
##  [640] "Heads" "Heads" "Tails" "Heads" "Tails" "Tails" "Tails" "Tails" "Heads"
##  [649] "Heads" "Tails" "Tails" "Heads" "Heads" "Heads" "Heads" "Tails" "Tails"
##  [658] "Heads" "Heads" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails" "Heads"
##  [667] "Heads" "Tails" "Tails" "Tails" "Tails" "Heads" "Tails" "Tails" "Heads"
##  [676] "Heads" "Heads" "Heads" "Heads" "Heads" "Tails" "Tails" "Heads" "Heads"
##  [685] "Heads" "Heads" "Heads" "Heads" "Tails" "Tails" "Tails" "Heads" "Heads"
##  [694] "Heads" "Tails" "Heads" "Heads" "Heads" "Tails" "Heads" "Tails" "Heads"
##  [703] "Tails" "Heads" "Heads" "Heads" "Heads" "Tails" "Tails" "Heads" "Tails"
##  [712] "Tails" "Heads" "Tails" "Tails" "Heads" "Tails" "Tails" "Tails" "Tails"
##  [721] "Heads" "Tails" "Tails" "Heads" "Tails" "Heads" "Heads" "Heads" "Heads"
##  [730] "Tails" "Tails" "Tails" "Tails" "Tails" "Heads" "Heads" "Tails" "Tails"
##  [739] "Tails" "Heads" "Heads" "Tails" "Tails" "Heads" "Heads" "Tails" "Tails"
##  [748] "Heads" "Heads" "Heads" "Tails" "Tails" "Tails" "Tails" "Heads" "Heads"
##  [757] "Tails" "Tails" "Heads" "Heads" "Tails" "Heads" "Heads" "Tails" "Heads"
##  [766] "Heads" "Heads" "Heads" "Tails" "Heads" "Heads" "Tails" "Tails" "Tails"
##  [775] "Tails" "Heads" "Tails" "Tails" "Heads" "Tails" "Tails" "Tails" "Tails"
##  [784] "Heads" "Tails" "Tails" "Heads" "Heads" "Heads" "Tails" "Tails" "Tails"
##  [793] "Heads" "Heads" "Tails" "Heads" "Tails" "Tails" "Heads" "Heads" "Tails"
##  [802] "Heads" "Tails" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails" "Heads"
##  [811] "Heads" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails"
##  [820] "Tails" "Heads" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails" "Tails"
##  [829] "Tails" "Tails" "Tails" "Tails" "Heads" "Tails" "Tails" "Tails" "Heads"
##  [838] "Heads" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Heads" "Tails"
##  [847] "Heads" "Heads" "Tails" "Heads" "Tails" "Tails" "Tails" "Tails" "Tails"
##  [856] "Tails" "Tails" "Heads" "Tails" "Heads" "Tails" "Tails" "Heads" "Heads"
##  [865] "Heads" "Tails" "Heads" "Heads" "Heads" "Heads" "Heads" "Heads" "Tails"
##  [874] "Tails" "Heads" "Heads" "Heads" "Heads" "Heads" "Tails" "Tails" "Heads"
##  [883] "Tails" "Tails" "Heads" "Tails" "Tails" "Tails" "Heads" "Tails" "Heads"
##  [892] "Heads" "Tails" "Heads" "Heads" "Heads" "Tails" "Tails" "Heads" "Heads"
##  [901] "Heads" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails"
##  [910] "Heads" "Heads" "Heads" "Tails" "Heads" "Tails" "Tails" "Heads" "Heads"
##  [919] "Heads" "Heads" "Tails" "Heads" "Heads" "Tails" "Tails" "Heads" "Heads"
##  [928] "Heads" "Tails" "Heads" "Tails" "Tails" "Heads" "Heads" "Tails" "Heads"
##  [937] "Tails" "Heads" "Tails" "Heads" "Tails" "Heads" "Heads" "Heads" "Heads"
##  [946] "Tails" "Heads" "Tails" "Tails" "Tails" "Heads" "Tails" "Heads" "Heads"
##  [955] "Tails" "Heads" "Heads" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails"
##  [964] "Tails" "Tails" "Heads" "Heads" "Heads" "Heads" "Tails" "Tails" "Tails"
##  [973] "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Heads" "Heads" "Heads"
##  [982] "Tails" "Heads" "Heads" "Heads" "Tails" "Tails" "Tails" "Heads" "Heads"
##  [991] "Heads" "Tails" "Heads" "Tails" "Tails" "Heads" "Tails" "Heads" "Heads"
## [1000] "Heads"

Prikaz rezultata eksperimenta

table(sample(coin, size = 1000, replace = TRUE))
## 
## Heads Tails 
##   502   498

Simulirajmo bacanje faličnog novcica

sample(coin, size = 100, replace = TRUE, prob = c(1/3,2/3))
##   [1] "Heads" "Heads" "Heads" "Heads" "Heads" "Tails" "Tails" "Tails" "Tails"
##  [10] "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails"
##  [19] "Tails" "Tails" "Tails" "Heads" "Tails" "Tails" "Tails" "Heads" "Tails"
##  [28] "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Heads" "Tails" "Tails"
##  [37] "Tails" "Tails" "Heads" "Tails" "Tails" "Heads" "Tails" "Tails" "Heads"
##  [46] "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails"
##  [55] "Heads" "Heads" "Tails" "Heads" "Heads" "Tails" "Tails" "Tails" "Tails"
##  [64] "Tails" "Heads" "Heads" "Tails" "Tails" "Tails" "Heads" "Heads" "Tails"
##  [73] "Heads" "Tails" "Heads" "Heads" "Heads" "Tails" "Tails" "Heads" "Tails"
##  [82] "Tails" "Heads" "Tails" "Tails" "Tails" "Tails" "Tails" "Heads" "Tails"
##  [91] "Heads" "Tails" "Heads" "Tails" "Tails" "Tails" "Heads" "Tails" "Tails"
## [100] "Heads"
table(sample(coin, size = 100, replace = TRUE, prob = c(1/3,2/3)))
## 
## Heads Tails 
##    33    67

Jos jedna funkcija za generisanje slucajnosti je runif.
Ona generiše slučajan broj iz intervala (0,1) (po zakonu neprekidne uniformne raspodjele).
runif(size, min = a, max = b)

r <- runif(1000, 0, 1)

length(r[r < 0.5]) / 1000
## [1] 0.484
# ili

mean(r<0.5)
## [1] 0.484
# Kad primijenimo arimeticku funkciju na logicki vektor, vrijednosti T i F se tretiraju kao 1 i 0.

length(r[r > 2 / 3]) / 1000
## [1] 0.327
# ili

mean(r>2/3)
## [1] 0.327

Liste

Lista je generički vektor koji sadrži druge objekte.

n <- c(2, 3, 5)
s <- c("aa", "bb", "cc", "dd", "ee")
b <- c(TRUE, FALSE, TRUE, FALSE, FALSE)
lst <- list(n, s, b)   # lst sadrzi  kopije od n, s, b

lst[1] #niz n
## [[1]]
## [1] 2 3 5
lst[2] #niz s
## [[1]]
## [1] "aa" "bb" "cc" "dd" "ee"

Na ovaj način dobijamo kopije prvog, odnosno drugog člana liste
Međutim ako hoćemo da direktno pristupimo članu liste koristimo [[ ]]

lst[[1]]   
## [1] 2 3 5

Tada možemo da mijenjamo sadržaj liste

lst[[2]]<-c("a","b")
lst[[2]][1]<-"c"

Simulacije slučajnih događaja

U narednim primjerima, uz pomoc slucajnih (pseudo slucajnih) funkcija simuliracemo neke slucajne dogadjaje i racunati relativnu ucestalost (frekvenciju) tih dogadjaja (kolicnik broja koliko puta se realizovao taj dogadjaj i ukupnog broja simulacija istog). Intuitivno je da je takav broj aproksimacija vjerovatnoce tog dogadjaja i takodje je prirodno pretpostaviti da cemo za veci broj ponavljanja simulacija dobiti precizniju aproksimaciju. Ovaj pristup cemo teorijski moci da opravdamo pred kraj kursa nakon nekih vaznih teorema.

  • Primjer 1: Simulirati bacanje regularne kockice za igru. Izračunati frekvenciju pojavljivanja šestice u 1000 izvođenja eksperimenta.

sample(6, 1, prob = rep(1/6,6))
## [1] 3
s <- sample(6, 1000, prob = rep(1/6,6), replace = TRUE)
mean(s == 6)  # ekvivalentan poziv sum(s==6)/length(s)
## [1] 0.166

  • Primjer 2: Kockica za igru je tako napravljena da je vjerovatnoća padanja nekog broja proporcionalna broju tackica na toj strani. Odrediti vjerovatnoću da padne paran broj.

s <- sample(1:6, 10000, replace = TRUE, prob = c(1/21,2/21,3/21,4/21,5/21,6/21))

mean(s%%2 == 0)
## [1] 0.5762

  • Primjer 3: Neka je dat novčić kod kojeg sa vjerovatnoćom \(p > 0.5\) pada glava. Smisliti fer žrijeb izmedju dva igrača koristeći taj novčić.

# Ideja je da definisemo koji ishodi predstavljaju pobjedu za igraca 1 odnosno 2.
# Npr neka se baca novcic dva puta i kazemo da je igrac 1 pobijedio ako se dogo-
# dio ishod pismo-glava, a igrac 2 ako se dogodio ishod glava-pismo. Na ovaj na-
# cin vjerovatnoce pobjede oba igraca su jednake. Provjerite ovaj rezultat
# eksperimentalno. 


nefer.igra <- function(p) {
# neka 0 predstavlja pismo, a 1 glavu

  novcic <- sample(c(0, 1), 2, replace = TRUE, prob = c(p, 1 - p))
  
  if (novcic[1] == 0 & novcic[2] == 1)
    return(1)
  else if (novcic[1] == 1 & novcic[2] == 0)
    return(2)
  else
    return(nefer.igra(p)) # ako se nije desilo PG/GP, bacamo opet 2 puta
}
# Sada hocemo da provjerimo odnos pobjeda igraca 1 odnosno 2 u 10000 odigranih
# partija. Neka je p=0.1, 1-p=0.9

r <- replicate(10000, nefer.igra(0.1))
sum(r == 1) / sum(r == 2)
## [1] 1.041233
# Odnos je priblizno isti, sto znaci da smo napravili fer igru.

  • Primjer 4: Na turniru treba igrati 3 partije stonog tenisa protiv sampiona A i slabijeg igraca B, po jednoj od sema A-B-A ili B-A-B. Nagrada se dobija ako se pobedi u bar dve partije uzastopno. Koja sema je povoljnija?

# Vazne funkcije koji se koriste u resenju: sapply i replicate!

# Simuliramo 1 turnir po datoj semi
# Argumenti p_a i p_b odredjuju verovatnoce da u partiji pobedi A, odnosno B.
# Argument sema odredjuje po kojoj semi radimo "ABA" ili "BAB"
# Vracamo 1 ako smo dobili nagradu, a 0 ako nismo
turnir <- function(p_a, p_b, sema = "ABA") {
  sema_vec <- unlist(strsplit(sema, split = ""))
  
  vrv_seme <- ifelse(sema_vec == "A", p_a, p_b)
  
  nase_pobede <- sapply(vrv_seme, function(p) {
    sample(c(0, 1), 1, prob = c(p, 1-p))
    }
  )
  
  return(nase_pobede[1] + nase_pobede[2] == 2 ||
           nase_pobede[2] + nase_pobede[3] == 2)
}

# Simulirajmo semu ABA 10000 puta i vidimo verovatnocu da pobedimo
mean(replicate(1e4, turnir(0.7, 0.2, "ABA")))
## [1] 0.41
# Simulirajmo semu BAB 10000 puta i vidimo verovatnocu da pobedimo
mean(replicate(1e4, turnir(0.7, 0.2, "BAB")))#
## [1] 0.2908
# Veca je verovatnoca da dobijemo po prvoj semi! (Uporedite rezultate sa teorisjkim resenjem sa casa.)

  • Primjer 5: Svaka od 15 ispitnih cedulja sadrzi 2 razlicita pitanja. Student zna odgovor na 25 pitanja. Da bi polozio ispit, mora da zna oba sa cedulje koju prvu izvuce ili prvo sa druge cedulje. Odrediti verovatnocu da student polozi.

# Simuliramo 1 polaganje studenta.
ispit <- function() {
  pitanja <- sample(1:30) # permutovacemo pitanja, da ne bi postojao poznat redosled
  cedulje <- matrix(pitanja, ncol = 2) # svaka cedulja je jedan red u matrici
  
  poznata_pitanja <- 1:25
  
  # Moguci problemi: Ako su pitanja 1:30, a poznata 1:25 (bez sample), verovatnoca
  # ispada 1, jer postoji poznata struktura cedulja, pa je prvo pitanje uvek poznato
  # (jer matrix se puni po kolonama, pa je prvo pitanje uvek iz prvih 15)
  
  # izvlaci prvu cedulju i ako zna oba pitanja, polozio je
  izvucena_cedulja <- sample(1:15, 1)
  izvucena_pitanja <- cedulje[izvucena_cedulja, ]
  if (all(izvucena_pitanja %in% poznata_pitanja)) {
    return(TRUE)
  } else if (any(izvucena_pitanja %in% poznata_pitanja)){
    # Ako ne zna oba, ali zna bar jedno, izvlaci opet
    ostale_cedulje <- cedulje[-izvucena_cedulja, ]
    # PROBLEM: Kada ne bismo izbacili vec izvucenu cedulju, nego opet vukli iz
    # 1:15, dobili bismo vrv 0.928, sto je neodoljivo blizu tacnog resenja, ali
    # nije 0.935 koliko treba da bude!
    
    novo_pitanje <- ostale_cedulje[sample(1:14, 1), 1] # prvo sa nove cedulje
    # ako zna prvo, polozio je
    if (novo_pitanje %in% poznata_pitanja) {
      return(TRUE)
    }
  }
  # u svim ostalim sljucajevima pada
  return(FALSE)
}

# simuliramo ispit 10000 puta da ocenimo verovatnocu
mean(replicate(1e5, ispit()))
## [1] 0.93611

  • Primjer 6: U voz koji ima \(m\) vagona penje se \(n\) (\(n\geq m\)) putnika. Odrediti verovatnocu da u svaki vagon udje bar po 1 putnik.

# simuliramo jedan raspored putnika u vozu
voz <- function(m, n) {
  # broj ljudi po vagonu je vektor duzine m
  vagoni <- numeric(m)
  
  # svaki putnik koji naidje bira jedan vagon slucajno
  for (i in 1:n) {
    izabrani_vagon <- sample(1:m, 1)
    vagoni[izabrani_vagon] = vagoni[izabrani_vagon] + 1
  }
  
  return(vagoni)
}

mean(replicate(1e4, all(voz(5, 10) > 0)))
## [1] 0.5192
# teorijsko resenje sa casa:
resenje <- function(m, n) {
  k <- 0:m
  clanovi_sume <- choose(m, k) * (-1)^k * ((m-k)/m)^n
  sum(clanovi_sume)
}
resenje(5, 10)
## [1] 0.5225472
# priblizno isto

Korisne funkcije u vezi sa poznatim raspodjelama

Postoje ugrađene funkcije za funkciju raspodjele, gustinu, funkcije kvantila, i generisanje slučajnih brojeva iz zadate raspodjele.
Nazivi raspodjela: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 \(d(x)=P\{X=x\}\))
  • kvantili - q - (\(p\)-kvantil je najmanje x za koje je \(F(x) \geq p\))
  • 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

Slucajni brojevi iz binomne raspodele: rbinom(n,size,prob), gde je n obim uzorka

rbinom(20, 10, 0.5)
##  [1] 6 6 4 5 7 7 4 5 4 3 4 4 4 5 7 7 6 6 2 7
#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 vjerovatnoce 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 raspodjela N(m,sigma^2)
# vrijednost 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) raspodjele
## [1] 0.09666703
# vrednost funkcije raspodjele u tacki z
# pnorm(z, mean=0, sd=1)
pnorm(0)
## [1] 0.5
pnorm(3, mean = 4)
## [1] 0.1586553
pnorm(3,
mean = 4,
)
## [1] 0.1586553
# crtanje gustine:

curve(dnorm(x), from = -5, to = 5)

# 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.09856476  0.09765878  0.88680605 -0.94145258  0.80530352  0.76363860
##  [7] -1.04399641 -0.11248777 -1.33952603  0.01380600 -0.97600379 -0.44669711
## [13] -0.21570089  0.99214310  0.08077420 -0.08577377 -0.33053157 -0.48942222
## [19] -0.49136838 -2.03563853
#Uniformna raspodela: U[a,b]
# gustina raspodjele:    min=a,max=b
dunif(2)
## [1] 0
dunif(2, min = 3, max = 5)  # U[3,5]
## [1] 0
# funkcija raspodjele:
punif(1, min = 0, max = 3)
## [1] 0.3333333
# 0.5 kvantil raspodjele (Fi^-1(0.5))
qunif(0.5)
## [1] 0.5
# generisanje slucajnih brojeva: 
runif(10)  # 10 brojeva iz U[0,1]
##  [1] 0.12422099 0.55177639 0.40773330 0.39498476 0.94824515 0.15083128
##  [7] 0.06416432 0.98159625 0.14066038 0.89267725
#Eksponencijalna raspodjela: Exp(lambda)

dexp(2, rate = 2) # f(2), za Exp(2) raspodjelu
## [1] 0.03663128
pexp(qexp(0.4))  # F(F^-1(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 simulacija

# U ~ U(0,1)

n <- 10 ^ 5

u <- runif(n)

mean(u <= 0.3)  # P{U<=0.3}
## [1] 0.30151
mean(u > 0.4)   # P{U>0.4}
## [1] 0.59973
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.5018
pnorm(0, sd = sqrt(10))
## [1] 0.5
# Z ~ Pois(3)

z <- rpois(n, lambda = 3)

mean(z == 2) # P{Z=2}
## [1] 0.22664
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 moze izraziti analiticki
  • Eksponencijalna raspodjela
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)
}
y<-Exp_random_samp(0.2, 1000)

hist (y, prob =TRUE , xlab ="",ylab ="",col ='cornflowerblue',border ='bisque' , main="Exp(0.2)")
curve ( dexp (x ,0.2) ,add =TRUE ,lwd =3, col ='coral1', xlim = c(0.1,40))

  • Erlangova raspodela \(\gamma(n,\lambda), n \in \mathcal{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] 5.070928
#Jos neke apsolutno neprekidne raspodele koje se mogu modelirati ovom metodom su
#Vejbulova, Freseova i Gumbelova
#Domaci: Generisati slucajne brojeve iz ovih raspodjela prethodno opisanom metodom za neki konkretan izbor parametara.

Neki primjeri apsolutno neprekidnih raspodela

# Neka je X ~ Exp(a). Generisati slucajne brojeve iz sledecih raspodela i izracunati srednju vrijednost.
# 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.6362813
# 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.3630189
# 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.1548

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 <- 10^4
n <- 10^5
s <- replicate(N, s_n(n))
hist(s,
probability = T,
col = 'lightblue',
main = "", ylab="")

# 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 = "", ylab="")
curve(dnorm(x),
add = T,
lwd = 2,
ylim=c(0,0.8),
col = 'coral')

Aproksimacija binomne raspodele normalnom

Poznato nam je da se binomna slucajna veličina \(\mathcal{B}(n,p)\) moze predstaviti kao zbir n nezavisnih slucajnih velicina sa Bernulijevom(p) raspodelom (tj. zbir n nezavisnih indikatora) i takođe \(ES_n=np, \ DS_n=np(1−p)\). Ako je np>10 možemo je aproksimirati normalnom.

# Generisacemo 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)
  
  
}
n <- 150
p <- 0.4
uzorak <- replicate(10000, binom(n, p))
# standardizujemo podatke
st.uzorak <- (uzorak - n*p) / sqrt(n*p*(1-p))

hist(st.uzorak,
probability = T,
col = 'lightblue',
main = "", xlab = "", ylab = "")
curve(dnorm(x),
add = T,
lwd = 2,
col = 'coral', ylab = "")

Zadaci:

  1. Racunar u procesu sabiranja brojeva vrsi zaokruzivanje na najblizi ceo broj. Pretpostavlja se da su greske nastale zaokruzivanjem nezavisne i uniformno raspodeljene na segmentu [−0.5, 0.5].
  1. Ako se sabira 1500 brojeva , izracunati verovatnocu da apsolutna vrednost ukupne greske bude veca od 15.
  2. Koliko se najvise brojeva moze sabrati, pa da sa verovatnocom 0.9 vrednost ukupne greske bude manja od 10?
#1.a 

1-(pnorm(1.34)-pnorm(-1.34))
## [1] 0.1802453
2-2*pnorm(1.34)
## [1] 0.1802453
#1.b 
qnorm(0.9)
## [1] 1.281552