Vektorizácia je spôsob, ako zrýchliť vyčíslenie výrazov, v ktorých sa rovnaká operácia aplikuje veľakrát na rôzne prvky.
Uvažujme napríklad situáciu, kde chceme vypočítať súčet logaritmov 1000 čísel:
x <- abs(rnorm(1000))
Takýto problém vieme riešiť napríklad pomocou for cyklu:
log_sum <- 0
for(i in 1:length(x)) log_sum <- log_sum + log(x[i])
Lepšie je ale aplikovať funkciu sum
priamo na vektor x
:
log_sum <- sum(log(x))
Úloha: Otestujte rýchlosť oboch implementácií pre \(n=10^7\).
monte_carlo <- function(N) {
hits <- 0
for (i in seq_len(N)) {
u1 <- runif(1)
u2 <- runif(1)
if (u1 ^ 2 > u2)
hits <- hits + 1
}
return(hits / N)
}
Rýchlosť:
N <- 500000
system.time(monte_carlo(N))
## user system elapsed
## 1.18 0.00 1.19
monte_carlo_vec <- function(N) sum(runif(N)^2 > runif(N))/N
N <- 500000
system.time(monte_carlo_vec(N))
## user system elapsed
## 0.01 0.00 0.01
Funkcie typu apply
sú určené na repetitívnu manipuláciu s podmnožinami dát v maticiach, poliach a dátových rámcoch. Používame ich na to, aby sme sa vyhli, častokrát vnoreným, for cyklom. Funckia, ktorú týmto spôsobom aplikujeme, môže byť napríklad agregačná (priemer, súčet), transformačná alebo iná vektorizovaná funkcia.
K dispozícii máme funkcie apply()
, lapply()
, sapply()
, vapply()
, mapply()
, rapply()
a tapply()
. V nasledujúcich príkladoch sa pozrieme na niektoré z nich.
Pri každom príklade si pozrite help danej funkcie a všimnite si, ako sa medzi sebou líšia výstupy jednotlivých funkcií.
Ako by sme počítali riadkové alebo stĺpcové súčty matice?
X <- matrix(rnorm(30), nrow=5, ncol=6)
apply(X, 2 ,sum)
## [1] -1.4062491 0.1178578 -1.9623120 0.1005855 1.1808389 -1.7279902
Úloha: Vypočítajte riadkové priemery matice X.
Pre matice máme v R špeciálne funkcie, ktoré nám umožňujú robiť podobné často sa vyskytujúce operácie, napr. colSums
alebo rowMeans
. Pre polia väčších rozmerov a pre menej bežné funkcie je ale funkcia apply užitočná.
pole <- array(1:(2*3*4), dim = c(2,3,4))
apply(pole, 3, function(x) sum(x^2))
## [1] 91 559 1459 2791
apply(pole, 1:2, var)
## [,1] [,2] [,3]
## [1,] 60 60 60
## [2,] 60 60 60
apply(pole, c(3,1), max)
## [,1] [,2]
## [1,] 5 6
## [2,] 11 12
## [3,] 17 18
## [4,] 23 24
apply(pole, c(1,3), max)
## [,1] [,2] [,3] [,4]
## [1,] 5 11 17 23
## [2,] 6 12 18 24
Funkcia apply aplikovaná na objekty typu list.
l <- list(a = 1:10, b = 11:20)
lapply(l, mean)
## $a
## [1] 5.5
##
## $b
## [1] 15.5
lapply(l, sum)
## $a
## [1] 55
##
## $b
## [1] 155
Príklad použitia: ak chceme vykonať rovnakú operáciu na viacerých dátových rámcoch, spojíme ich do listu a použijeme lapply.
data(iris)
data(swiss)
lst <- list(iris=iris[,-5], swiss=swiss[,-1])
lapply(lst, colMeans)
## $iris
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.843333 3.057333 3.758000 1.199333
##
## $swiss
## Agriculture Examination Education Catholic
## 50.65957 16.48936 10.97872 41.14383
## Infant.Mortality
## 19.94255
l <- list(a = 1:10, b = 11:20)
l.mean <- sapply(l, mean)
class(l.mean)
## [1] "numeric"
l.mean[1]
## a
## 5.5
l.mean[['a']]
## [1] 5.5
Výstupom sapply
je v prípade, že oba listy sú rovnakej veľkosti, matica, nie list, čo môže často byť pohodlnejšie na následnú manipuláciu.
lst <- list(iris=iris[,-5], swiss=swiss[,-(1:2)])
sapply(lst, colMeans)
## iris swiss
## Sepal.Length 5.843333 16.48936
## Sepal.Width 3.057333 10.97872
## Petal.Length 3.758000 41.14383
## Petal.Width 1.199333 19.94255
Funkcia tapply
užitočná najmä pre dátové súbory. Vezmime si dáta obsahujúce mená detí v USA v rokoch 1880-2017. Kvôli prehľadnosti budeme uvažovať len roky 2015-2017.
library(babynames)
## Warning: package 'babynames' was built under R version 4.1.3
data(babynames)
bn <- babynames[babynames$year>=2015, ]
tapply(X=bn$n, INDEX=bn$sex, FUN=sum, na.rm=TRUE)
## F M
## 5254610 5633346
tapply(X=bn$n, INDEX=list(bn$sex, bn$year), FUN=sum, na.rm=TRUE)
## 2015 2016 2017
## F 1778883 1763916 1711811
## M 1909804 1889052 1834490
Uvažujme takúto dvojitú sumu: \[f(n)= \sum_{r=1}^n \sum_{s=1}^r \frac{s^2}{10+4r^3}\]
Ukážeme si, ako túto sumu vypočítať rôznymi (viac alebo menej efektívnymi) spôsobmi.
suma1 <- function(n){
sum <- 0
for (r in 1:n){
for (s in 1:r)
sum <- sum+s^2/(10+4*r^3)
}
sum
}
suma2 <- function(n){
mat <- matrix(0,ncol=n,nrow=n)
sum( (col(mat)^2) / (10+4 * row(mat) ^ 3) * (col(mat) <= row(mat)) )
}
sapply
suma3 <- function(n){
f <- function(r){ sum(((1:r)^2 / (10 + 4 * r^3)))}
sum(sapply(1:n, FUN=f))
}
suma4 <- function(n)
sum(outer(1:n, 1:n, FUN = function(r,s){ (s<=r)*(s^2)/(10+4*r^3) }))
Pomocou system.time
zistite, ktorá implementácia je najrýchlejšia a ktorá najpomalšia
system.time(suma1(1000))
## user system elapsed
## 0.09 0.00 0.09