Најједноставнији график

x <- seq(2, 16, 1)
y <- c(5, 6, 5, 9, 3, 11, 6, 15, 19, 5, 7, 14, 8, 12, 1)
plot(x, y)

Додатна подешавања при цртању графика

главни и споредни наслов, називи оса: main, sub, xlab, ylab
дужину оса можемо и експлицитно навести параметрима xlim и ylim
додавање тачака на график points
додавање текста на график text
линије и текст се могу додати у оквиру позива функцијеplot или их можемо накнадно додати
За додавање линије на график користи се функција abline која црта график функције \(y=a+bx\).
Ако је агрумент \(h=...\) или \(v=...\), црта хоризонталну, односно вертикалну линију са датом координатом.

plot(x, y, main = "Main title", sub = "subtitle", xlab = "x-osa", ylab = "y-osa")

plot(x, y, xlim = c(0,30), ylim = c(0,30))
points(25, 25) # Zadajemo koordinate tacke koju dodajemo
text(3, 3, "ubaceni tekst") # Zadajemo koordinate i tekst koji se ispisuje 
abline(0, 1) # Dobili smo grafik f-je y=x
abline(h = 20)
abline(v = 20) # Isto kao da smo napisali i abline(h = 20,v = 20)

Ако хоћемо да график буде веома сложен, често је добра идеја да на почетк изоставимо све елементе и после их појединачно додајемо.

plot(x, y, type = "n", xlab = "", ylab = "", axes = F) # Ne crta ose, ne obelezava ih, type="n" znaci da se tacke izostavljaju
points(x, y) # Dodaje sve tacke na grafik
axis(1) # Dodaje prvu osu
axis(2, at = seq(0, 20, 5)) # Dodaje drugu osu, ali je obelezava na mestima naznacenim vektorom seq(..)
box() # Dodaje okvir
title(main = "Glavni naslov", sub = "podnaslov", xlab = "x-osa", ylab = "y-osa")

График функције \(f(x)=x^3\)

x <- seq(-2,2,0.01)
plot(x, x^3, type = "l", main = "y=x^3", xlab ="x-osa", ylab ="y-osa", col = "blue")

curve(x^3, from = -2, to = 2, col = "blue", main = "y=x^3", xlab = "x-osa", ylab = "y-osa")

Функција par се користи за напреднија подешавања свих графичких параметара. Нека од ових подешавања могу и преко plot функције. Разлика је што их функција par трајно подешава, док се не врате на старо. Неке од опција су боја, дељина линије, тип линије, подела оса, величина графика, подела на више површина, итд.

Дељење платна за цртање

Параметар mfrow означава колико се графика одједном може налазити на платну. То је вектор димензијa

par(mfrow = c(1, 2)) # platno se deli na jednu vrstu i dve kolone (dva grafika jedan pored drugog)
plot(x, y)
plot(y, x)

# Vracamo sve na default
par(mfrow = c(1, 1), col = "black", bg = "white")
plot(x, y)

Комбиновање више графика један преко другог

x1 <- rnorm(20) #  Grafici uzorka normalne raspodele
par(mfrow = c(2,2)) 
plot(x1, type = "p", main = "points", ylab = "y-osa", xlab = "x-osa", col = "red")
plot(x1, type = "l", main = "lines", ylab = "y-osa", xlab = "x-osa", col = "orange")
plot(x1, type = "b", main = "both", ylab = "y-osa", xlab = "x-osa", col = "blue")
plot(x1, type = "o", main = "both overplot", ylab = "y-osa", xlab = "x-osa", col = "green")

Хистограм и функција густине

Да ли је расподела симетрична?

data(airquality)
hist(airquality$Temp)

# granice histograma
broj_kat <- ceiling(log(length(airquality$Temp), 2))+1
d <- (max(airquality$Temp) - min(airquality$Temp)) / broj_kat
granice <- min(airquality$Temp) - 1/2 + (0:broj_kat) * (d+1)
hist(airquality$Temp, breaks = granice, prob = TRUE)

attach(airquality)
hist(Temp, col = "darkmagenta", freq = FALSE)

Колико је процентуално времена температура била мања од 70? Већа од 90?

H <- hist(Temp, prob = TRUE, plot = FALSE)

## Warning in hist.default(Temp, prob = TRUE, plot = FALSE): argument 'probability'
## is not made use of

cumsum(H$density*diff(H$breaks))

## [1] 0.05228758 0.11764706 0.21568627 0.33986928 0.55555556 0.77777778 0.90849673
## [8] 0.98692810 1.00000000

# 21.57% ; 9.15%
hist(Temp, col = "darkmagenta", freq = FALSE, breaks = seq(50,110,5))
lines(density(Temp))

hist(Ozone, col = "magenta", freq = FALSE)

hist(Ozone, col = "magenta", freq = FALSE, breaks = seq(0,180,5))
lines(density(na.omit(Ozone)))

hist(Solar.R, col = "pink3")

hist(Solar.R, col = "pink3", freq = FALSE, breaks = seq(0,360,50))
lines(density(na.omit(Solar.R)))

hist(Wind, col = "lightblue")

hist(Wind, col = "lightblue", freq = FALSE, breaks = seq(0,22,2))
lines(density(Wind))

detach(airquality)
attach(faithful)
hist(eruptions, breaks = seq(1.4, 5.2, 0.2), prob = T, col = "orange")
lines(density(eruptions))

detach(faithful)
attach(mtcars)
hist(mpg, main = "Potrosnja goriva", col = "aquamarine")

hist(mpg, prob = TRUE, main = "Potrosnja goriva", col = "aquamarine")
lines(density(mpg), col = "red", lwd = 2)

mean(mpg)

## [1] 20.09062

median(mpg)

## [1] 19.2

diff(range(mpg))

## [1] 23.5

var(mpg)

## [1] 36.3241

sd(mpg)

## [1] 6.026948

quantile(mpg)

##     0%    25%    50%    75%   100% 
## 10.400 15.425 19.200 22.800 33.900

IQR(mpg)

## [1] 7.375

plot(mpg)
abline(mean(mpg), 0, col = "pink", lwd = 3)

hist(hp, main = "Konjska snaga", col = "cyan")

hist(hp, prob = TRUE, main = "Konjska snaga", col = "cyan", breaks = seq(45,350,10))

hist(hp, prob = TRUE, main = "Konjska snaga", col = "cyan", breaks = seq(45,350,20))

hist(hp, prob = TRUE, main = "Konjska snaga", col = "cyan")
lines(density(hp), col = "red", lwd = 2)

mean(hp)

## [1] 146.6875

median(hp)

## [1] 123

diff(range(hp))

## [1] 283

var(hp)

## [1] 4700.867

sd(hp)

## [1] 68.56287

quantile(hp)

##    0%   25%   50%   75%  100% 
##  52.0  96.5 123.0 180.0 335.0

plot(hp)
abline(mean(hp), 0, col = "pink", lwd = 3)

detach(mtcars)

Кутијасти дијаграм - BOXPLOT

Обједињује мере положаја и расејања и открива аутлајере

boxplot(mtcars$mpg, main = "Potrosnja goriva", col = "aquamarine")
points(mean(mtcars$mpg), col = 'red', pch = 16)

boxplot(mtcars$hp, main = "Konjska snaga", col = "cyan")
points(mean(mtcars$hp), col = 'red', pch = 16)

boxplot(mtcars$mpg ~ mtcars$vs, main = "Potrosnja goriva u zavisnosti od vrste motora" , col ="lightblue2")

boxplot(mtcars$mpg ~ mtcars$am, main = "Potrosnja goriva u zavisnosti od vrste menjaca" , col = "lightblue")

data(InsectSprays)
attach(InsectSprays)
boxplot(count ~ spray, data = InsectSprays, col = "lightgray")
means <- tapply(count, spray, mean)
points(means, col = "red", pch = 18)

boxplot(count ~ spray, data = InsectSprays, horizontal = T, col = "lightblue")

Који аутлајери су благи, а који екстремни?

Табеларни приказ и тракасти дијаграм - BARPLOT

summary(mtcars)

##       mpg             cyl             disp             hp       
##  Min.   :10.40   Min.   :4.000   Min.   : 71.1   Min.   : 52.0  
##  1st Qu.:15.43   1st Qu.:4.000   1st Qu.:120.8   1st Qu.: 96.5  
##  Median :19.20   Median :6.000   Median :196.3   Median :123.0  
##  Mean   :20.09   Mean   :6.188   Mean   :230.7   Mean   :146.7  
##  3rd Qu.:22.80   3rd Qu.:8.000   3rd Qu.:326.0   3rd Qu.:180.0  
##  Max.   :33.90   Max.   :8.000   Max.   :472.0   Max.   :335.0  
##       drat             wt             qsec             vs        
##  Min.   :2.760   Min.   :1.513   Min.   :14.50   Min.   :0.0000  
##  1st Qu.:3.080   1st Qu.:2.581   1st Qu.:16.89   1st Qu.:0.0000  
##  Median :3.695   Median :3.325   Median :17.71   Median :0.0000  
##  Mean   :3.597   Mean   :3.217   Mean   :17.85   Mean   :0.4375  
##  3rd Qu.:3.920   3rd Qu.:3.610   3rd Qu.:18.90   3rd Qu.:1.0000  
##  Max.   :4.930   Max.   :5.424   Max.   :22.90   Max.   :1.0000  
##        am              gear            carb      
##  Min.   :0.0000   Min.   :3.000   Min.   :1.000  
##  1st Qu.:0.0000   1st Qu.:3.000   1st Qu.:2.000  
##  Median :0.0000   Median :4.000   Median :2.000  
##  Mean   :0.4062   Mean   :3.688   Mean   :2.812  
##  3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.:4.000  
##  Max.   :1.0000   Max.   :5.000   Max.   :8.000

mtcars$cyl <- as.factor(mtcars$cyl)
mtcars$vs <- as.factor(mtcars$vs)
mtcars$am <- as.factor(mtcars$am)
mtcars$gear <- as.factor(mtcars$gear)
mtcars$carb <- as.factor(mtcars$carb)
summary(mtcars)

##       mpg        cyl         disp             hp             drat      
##  Min.   :10.40   4:11   Min.   : 71.1   Min.   : 52.0   Min.   :2.760  
##  1st Qu.:15.43   6: 7   1st Qu.:120.8   1st Qu.: 96.5   1st Qu.:3.080  
##  Median :19.20   8:14   Median :196.3   Median :123.0   Median :3.695  
##  Mean   :20.09          Mean   :230.7   Mean   :146.7   Mean   :3.597  
##  3rd Qu.:22.80          3rd Qu.:326.0   3rd Qu.:180.0   3rd Qu.:3.920  
##  Max.   :33.90          Max.   :472.0   Max.   :335.0   Max.   :4.930  
##        wt             qsec       vs     am     gear   carb  
##  Min.   :1.513   Min.   :14.50   0:18   0:19   3:15   1: 7  
##  1st Qu.:2.581   1st Qu.:16.89   1:14   1:13   4:12   2:10  
##  Median :3.325   Median :17.71                 5: 5   3: 3  
##  Mean   :3.217   Mean   :17.85                        4:10  
##  3rd Qu.:3.610   3rd Qu.:18.90                        6: 1  
##  Max.   :5.424   Max.   :22.90                        8: 1

table(mtcars$carb)

## 
##  1  2  3  4  6  8 
##  7 10  3 10  1  1

table(mtcars$am)

## 
##  0  1 
## 19 13

table(mtcars$am, mtcars$cyl)

##    
##      4  6  8
##   0  3  4 12
##   1  8  3  2

table(mtcars$carb, mtcars$cyl)

##    
##     4 6 8
##   1 5 2 0
##   2 6 0 4
##   3 0 0 3
##   4 0 4 6
##   6 0 1 0
##   8 0 0 1

barplot(table(mtcars$carb), col = "deepskyblue", main = "Broj karburatora")

barplot(table(mtcars$carb, mtcars$cyl), col = rainbow(6), main = "Broj karburatora u zavisnosti od broja cilindara", xlab = "broj cilindara")

barplot(t(table(mtcars$carb, mtcars$cyl)), col = c("deeppink", "darkviolet","limegreen"), main = "Broj cilindara u zavisnosti od broja karburatora", xlab = "broj karburatora")

barplot(t(table(mtcars$carb, mtcars$cyl)), beside = TRUE, main = "Broj cilindara u zavisnosti od broja karburatora", xlab = "broj karburatora")

Кружни дијаграм - PIECHART

par(mfrow = c(2,2))
pie.sales <- c(0.12, 0.3, 0.26, 0.16, 0.04, 0.12)
names(pie.sales) <- c("Blueberry", "Cherry",
                      "Apple", "Boston Cream", "Other", "Vanilla")
pie(pie.sales, main = "Obicna pitica")
pie(pie.sales, col = gray(seq(0.4, 0.9, length = 6)), clockwise = TRUE, main = "Nijanse sive")
pie(pie.sales, col = rainbow(6), clockwise = TRUE, main = "Boje duge")
library(plotrix)
pie3D(pie.sales, main = "3D pitica")

Погледати пакете lattice, ggplot2, plotrix,…

Функција sample()

У најједноставнијем облику функција sample пермутује задати низ на случајан начин.

x <- 1:10
sample(x)

##  [1]  4  1  2  5  3  9  6  8  7 10

sample(10)

##  [1]  7  5  6  9  3  2  8  4 10  1

Ова функција има аргумент replace који може узимати вредности FALSE или TRUE. Подразумевано, ако није другачије наведено, узима вредност FALSE, што значи да се симулира случајно бирање бројева без понављања. Ако проследимо вредност TRUE функција враћа случајан узорак са понављањем.

Још један аргумент је prob (односно probability) и он је подразумевано NULL, али можемо и да задамо неке тежине (вероватноће) елементима скупа.

Аргумент size одређује колико елемената бирамо из прослеђеног вектора.

НАПОМЕНА: Може доћи до забуне приликом следећих израза. Зашто?

sample(x[x > 8])

## [1] 10  9

sample(x[x > 9])

##  [1]  5 10  2  1  9  4  7  6  3  8

Симулирајмо бацање новчића

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

## [1] "Heads"

Желимо 100 симулација овог случајног експеримента.

# sample(coin, size = 100)

Зашто јавља грешку?

sample(coin, size = 100, replace = TRUE)

##   [1] "Tails" "Tails" "Heads" "Heads" "Heads" "Heads" "Heads" "Tails" "Heads"
##  [10] "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Heads" "Heads" "Tails"
##  [19] "Tails" "Heads" "Heads" "Tails" "Heads" "Tails" "Heads" "Heads" "Tails"
##  [28] "Heads" "Heads" "Heads" "Heads" "Tails" "Tails" "Tails" "Tails" "Heads"
##  [37] "Heads" "Heads" "Heads" "Tails" "Tails" "Heads" "Tails" "Tails" "Tails"
##  [46] "Heads" "Tails" "Tails" "Tails" "Heads" "Heads" "Heads" "Tails" "Tails"
##  [55] "Heads" "Heads" "Tails" "Heads" "Heads" "Tails" "Heads" "Tails" "Tails"
##  [64] "Heads" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Heads" "Heads"
##  [73] "Heads" "Tails" "Tails" "Heads" "Tails" "Heads" "Tails" "Heads" "Tails"
##  [82] "Heads" "Heads" "Heads" "Tails" "Heads" "Heads" "Heads" "Tails" "Heads"
##  [91] "Heads" "Heads" "Tails" "Tails" "Tails" "Heads" "Tails" "Tails" "Heads"
## [100] "Tails"

Приказ резултата експеримента

table(sample(coin, size = 1000, replace = TRUE))

## 
## Heads Tails 
##   523   477

Симулација бацања фаличног новчића

sample(coin, size = 100, replace = TRUE, prob = c(1/3,2/3))

##   [1] "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Heads" "Heads" "Heads"
##  [10] "Heads" "Tails" "Heads" "Heads" "Heads" "Tails" "Tails" "Tails" "Heads"
##  [19] "Tails" "Heads" "Heads" "Tails" "Heads" "Tails" "Tails" "Tails" "Tails"
##  [28] "Tails" "Heads" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails"
##  [37] "Heads" "Tails" "Tails" "Heads" "Tails" "Tails" "Heads" "Heads" "Tails"
##  [46] "Heads" "Tails" "Tails" "Heads" "Tails" "Tails" "Tails" "Tails" "Heads"
##  [55] "Tails" "Heads" "Heads" "Tails" "Tails" "Tails" "Tails" "Tails" "Heads"
##  [64] "Heads" "Tails" "Tails" "Tails" "Tails" "Tails" "Heads" "Tails" "Heads"
##  [73] "Heads" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails" "Tails"
##  [82] "Tails" "Heads" "Heads" "Heads" "Tails" "Tails" "Tails" "Tails" "Tails"
##  [91] "Tails" "Tails" "Tails" "Heads" "Tails" "Heads" "Tails" "Tails" "Tails"
## [100] "Tails"

table(sample(coin, size = 100, replace = TRUE, prob = c(1/3,2/3)))

## 
## Heads Tails 
##    38    62

Уочава се разлика у фреквенцијама.

Функција runif()

функција runif генерише случајан број из интервала \((0,1)\) по закону униформне расподеле. Ова функција може да се користи за генерисање случајности.

r <- runif(1000, 0, 1)
length(r[r < 0.5]) / 1000

## [1] 0.503

length(r[r > 2 / 3]) / 1000

## [1] 0.321

Задатак 1

Симулирати бацање регуларне коцкице за игру. Израчунати фреквенцију појављивања 6 у 1000 бацања.

sample(6, 1, prob = rep(1/6, 6))

## [1] 4

s <- sample(6, 1000, prob = rep(1/6, 6), replace = TRUE)
mean(s == 6)

## [1] 0.185

Задатак 2

Коцкица за игру је таква да је вероватноћа падања неког броја пропорционална броју тачкица на тој страни. Одредити вероватноћу да падне паран број.

s <- sample(1:6, 10000, replace = TRUE, prob = 1:6/sum(1:6))
mean(s%%2 == 0)

## [1] 0.5801

Задатак 3

Дат је фаличан новчић, где је вероватноћа да падне глава једнака \(p>0.5\). Осмислити фер игру са два играча, тј. игру у којој је једнако вероватно да победи сваки од играча.

Како дефинисати исходе тако да игра буде фер?

Новчић се баца два пута и кажемо да је победио играч 1 ако се догодио исход писмо-глава, а победио је играч 2 ако се догодио исход глава-писмо.

На овај начин су вероватноће победе оба играча једнаке. Проверимо овај резултат експериментално.

igra <- function(p) {
novcic <- sample(c("G","P"), 2, replace = TRUE, prob = c(p, 1 - p))
if (novcic[1] == "G" & novcic[2] == "P")
return (1)
else if (novcic[1] == "P" & novcic[2] == "G")
return (2)
else
return (0)
}

Проверимо однос победа играча 1, односо играча 2 у 10000 одиграних партија.

Нека је, на пример, \(p=0.1, 1-p=0.9\).

r <- replicate(10000, igra(0.1))
sum(r == 1) / sum(r == 2)

## [1] 0.9644779

Однос партија је приближно исти, близак 1, што значи да смо направили фер игру.

Задатак 4

Играч има почетни капитал \(X\) хиљада динара и он баца коцкицу. Ако добије 1, 3 или 5, његов капитал се увећа за 1000 динара, а у супротном он губи 1000 динара. Нацртати трајекторију његовог капитала након \(N\) бацања новчића (дозвољено је да играч иде у минус).

trajektorija <- function(X, N) {
  a <- numeric(N)
  s <- sample(6, N, replace = TRUE)
  a[1] <- X
  for (i in 1:N) # ifelse(test, yes, no)
  ifelse(s[i] %% 2 == 1, a[i + 1] <- a[i] + 1, a[i + 1] <- a[i] - 1)
    return(a)
}

Цратамо трајекторију за 20 бацања коцкице и почетни капитал од \(X=5\) хиљада динара.

plot(trajektorija(5, 20), type = "o", col = "blue")

Други начин: Ако желимо да избегнемо петљу.

Кумулативна сума вектора \(x=(x_1,x_2,...,x_n)\) је вектор \(x=(x_1,x_1+x_2, ... , x_1+x_2+...+x_n)\).

trajektorija2 <- function(X, N) {
s <- sample(6, N, replace = T)
i <- ifelse(s %% 2 == 1, 1,-1)
trajektorija <- cumsum(i) + X
return(trajektorija)
}
plot(trajektorija2(5, 20), type = "o", col = "red")

Приметимо да је други начин ефикаснији.

system.time(trajektorija(1, 10000))

##    user  system elapsed 
##    0.02    0.00    0.01

system.time(trajektorija2(1, 10000))

##    user  system elapsed 
##       0       0       0

Оценити вероватнићу да играч након \(N\) партија има бар једну хиљаду.

plus <- function(X, N) {
  a <- trajektorija2(X, N)
  ifelse(a[N] > 0, return(1), return(0))
}

frekv <- function(X, N, broj.sim = 1000) {
simulacije <- replicate(broj.sim, plus(X, N))
f <- mean(simulacije)
return(f)
}

frekv(5, 20)

## [1] 0.869

frekv(5, 100)

## [1] 0.712

frekv(100, 1000)

## [1] 1

Графичко приказивање података и симулације у R-у

Бојана Тодић