Лекція 7. Векторизовані обчислення в R з використанням apply-функцій

       Використання оператора циклів for(), що є характерним для мов низького рівня програмування типу Basic, вважається поганим стилем програмування на мові R, орієнтованою на векторизацію обчислень. Замість виконання послідовно скалярних операцій над кожним з елементів масиву, значно ефективніше виконувати паралельні обчислення, за яких програма обробляє одночасно весь масив (вектор) цілком або поелементно в кожен момент часу. Очевидно, що такий підхід потенційно може призвести до значного прискорення однотипних обчислень над великими масивами даних.

       Розглянемо найпростіший приклад векторизованих обчислень в R. Припустимо, у нас є вектор з 10 додатних чисел, і ми хочемо обрахувати квадратний корінь кожного із них. Замість написання циклу для почергового виконання цієї операції над кожним елементом, ми просто подаємо весь цей вектор на вхід функції sqrt(), яка повертає вектор із результатами обчислень:

x <- 1:10
sqrt(x)

 [1] 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751
 [8] 2.828427 3.000000 3.162278

       Принцип векторизованих обчислень застосовують до більш складних об’єктів R – матриць, списків і таблиць даних (для мови R немає різниці між останніми двома типами об’єктів: фактично таблиця даних є списком із декількох компонентів – векторів однакового розміру). У базовій комплектації R є ціле сімейство функцій, призначених для організації векторизованих обчислень над такими об’єктами. У назві всіх цих функцій є слово apply (з англ. застосувати), якому передує літера, що вказує на принцип роботи тієї чи іншої функції (див. докладніше в довідковому файлі ?apply). За такої умови:

• apply() на відміну від for() можна легко розпаралелити, перейменувавши функцію в її паралельну версію з пакету snow;
• алгоритми, записані без циклів, легше змінювати, містять менше помилок і легше набираються в командному рядку;
• результат роботи apply() може бути аргументом функції; крім того, будь-яка функція може бути вставлена в якості аргументу в apply(). Нижче надамо огляд apply-функцій із прикладами використання.

apply()

       Функція apply() використовується у випадках, коли необхідно застосувати будь-яку функцію до всіх рядків або стовпців матриці (або масивів більшої розмірності):

apply(x, MARGIN, FUN, ...)

де x – це перетворений об’єкт; MARGIN – індекс, що позначає напрямок процесу обчислень (стовпцями або рядками); FUN – функція, що використовуються для обчислень; … – це будь-які інші параметри функції FUN. Для матриці або таблиці даних MARGIN=1 позначає рядки, а MARGIN=2 – стовпці. Оскільки FUN означає будь-яку функцію R, то функція apply() – це потужний засіб модульної обробки даних. Створимо звичайну двовимірну матрицю:

M <- matrix(seq(1,16), 4, 4)
# Знайдемо мінімальні значення в кожному рядку матриці:
apply(M, 1, min)

[1] 1 2 3 4

# Знайдемо мінімальні значення в кожному стовпці матриці
apply(M, 2, max)

[1]  4  8 12 16

# Приклад із тривимірним масивом:
M <- array(seq(32), dim=c(4,4,2))
# Застосуємо функцію sum() до кожного елементу M[*,],
# тобто виконаємо додавання за вимірами 2 і 3:
apply(M, 1, sum)

[1] 120 128 136 144

# Результат – одновимірний вектор:
# Застосуємо функцію sum() до кожного елементу M[*, *,],
# – тобто виконаємо додавання за третім виміром:
apply(M, c(1,2), sum)

     [,1] [,2] [,3] [,4]
[1,]   18   26   34   42
[2,]   20   28   36   44
[3,]   22   30   38   46
[4,]   24   32   40   48

# Результат – матриця:

       У разі потреби обчислення суми і середніх значень за рядками або стовпчиках матриць використовувати дуже швидкі і спеціально оптимізовані для цього функції colSums(), rowSums(), colMeans() і rowMeans().

lapply()

       Функція lapply() використовується за потреби застосування певної функції до кожної компоненти списку. У такому випадку одержують результат у вигляді списку (буква “l” в назві lapply() означає list – “список”).

# Створимо список із трьома компонентами-векторами:
x <- list(a=1, b=1:3, c=10:100)
# З’ясуємо розмір кожної компоненти списку х (функція length()):
lapply(x, FUN=length)

$a
[1] 1

$b
[1] 3

$c
[1] 91

# Виконаємо додавання елементів у кожній компоненті списку х:
lapply(x, FUN=sum)

$a
[1] 1

$b
[1] 6

$c
[1] 5005

sapply()

       Функція sapply() використовується, якщо потрібно застосувати деяку функцію до кожної компоненти списку, і отримати резульат у вигляді вектору (буква “s” в назві sapply() означає simplify – “спростити”).

# Список із трьох компоненті:
x <- list(a=1, b=1:3, c=10:100)
# З’ясуємо розмір кожної компоненти списку х:
sapply(x, FUN=length)

 a  b  c 
 1  3 91 

# Додавання всіх елементів у кожній компоненті списку х:
sapply(x, FUN=sum)

   a    b    c 
   1    6 5005 

       У деяких більш “просунутих” випадках sapply() може видати результат у вигляді багатовимірного масиву. Наприклад, якщо застосовувана функція повертає вектори однакової довжини, sapply() об’єднає ці вектори в матрицю (по стовпцях):

sapply(1:5, function(x) rnorm(3, x))

         [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 1.858850 3.217806 3.908056 3.431791 6.384863
[2,] 1.453731 2.875744 3.221383 5.675227 6.743612
[3,] 1.485595 1.580141 4.426198 2.820517 5.828481

       Якщо застосовувана функція повертає матрицю, то sapply() перетворить кожну матрицю у вектор і об’єднає такі вектори в одну велику матрицю:

sapply(1:5, function(x) matrix(x, 2, 2))

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    1    2    3    4    5
[3,]    1    2    3    4    5
[4,]    1    2    3    4    5

       Роботу функції sapply(), що наведено в останньому прикладі, можна скасувати з допомогою аргументу simplify="array" – у цьому випадку матриці будуть об’єднані в один багатовимірний масив:

sapply(1:5, function(x) matrix(x,2,2), simplify="array")

, , 1

     [,1] [,2]
[1,]    1    1
[2,]    1    1

, , 2

     [,1] [,2]
[1,]    2    2
[2,]    2    2

, , 3

     [,1] [,2]
[1,]    3    3
[2,]    3    3

, , 4

     [,1] [,2]
[1,]    4    4
[2,]    4    4

, , 5

     [,1] [,2]
[1,]    5    5
[2,]    5    5

replicate()

       Функція replicate() є свого роду “обгорткою” для функції sapply() і дозволяє провести серію обчислень із метою генерації набору чисел відповідно до заданого алгоритму. Синтаксис функції має вигляд:

replicate(n, expr, simplify=TRUE)

де n – число повторів, expr – функція або група виразів, які треба повторити n раз, simplify=TRUE – необов’язковий параметр, який пробує спростити результат і представити його у вигляді вектору або матриці значень.

       Розглянемо приклад використання функції replicate() для перевірки статистичної гіпотези про рівність медіан двох вибірок бутстреп-методом (див. розділ 2.7). Створена нами функція boot_med() з кожної вихідної вибірки x і y вилучає по N псевдовибірок, використовуючи алгоритм “випадкового вибору з поверненням”, і знаходить різницю їхніх медіан. Тут функція sample.int() формує цілочисельні набори випадкових індексів indx і indy для кожного з порівнюваних векторів:

boot_med <- function(x, y, N=100) {
  replicate(N, {
    indx <- sample.int(length(x), length(x), replace=T)
    indy <- sample.int(length(y), length(y), replace=T)
    median(x[indx]) - median(y[indy])
    })
  }

       Для тестування функції перевіримо однорідність медіан для двох випадкових вибірок із нормального розподілу із середнім mean=10 і стандартним відхиленням sd=4:

Y <- rnorm(100, sd=4, mean=10)
X <- rnorm(100, sd=4, mean=10)
quantile(boot_med(Y, X, 10000), probs=c(0.025, 0.975))

      2.5%      97.5% 
-0.3757855  1.9518084 

       Оскільки довірчий інтервал різниці медіан включає 0, то нульову гіпотезу відкидати не слід.

vapply()

       Функція vapply() схожа з sapply(), але працює трохи швидше завдяки тому, що користувач однозначно вказує тип значень, що повертаються (буква “v” в назві vapply() означає velocity – “швидкість”). Такий підхід дозволяє також уникати повідомлень про помилки (і переривання обчислень), що виникають під час роботи з sapply() у деяких ситуаціях. Під час виклику vapply() користувач повинен навести приклад очікуваного типу значень, що повертаються. Для цього служить аргумент FUN.VALUE:

# Аргументу FUN.VALUE присвоєно логічне значення FALSE.
# Цим задається тип значень, що повертаються функцією
a <- vapply(NULL, is.factor, FUN.VALUE=FALSE)
# Функція sapply() застосована до того ж NULL-об’єкту:
b <- sapply(NULL, is.factor)
# Перевірка типу змінних:
is.logical(a)

[1] TRUE

is.logical(b)

[1] FALSE

mapply()

       Функція mapply() використовується за потреби поелементного застосування будь-якої функції одночасно до декількох об’єктів (наприклад, отримати суму перших елементів векторів, потім суму других елементів векторів, і т.ін.). Результат повертається у вигляді вектору або масиву іншої розмірності (див. приклади для sapply() вище). Буква “m” в назві mapply() означає multivariate – “багатовимірний” (мається на увазі одночасне виконання обчислень над елементами декількох об’єктів).

mapply(sum, 1:5, 1:5, 1:5)

[1]  3  6  9 12 15

       Функція rapply() використовується у випадках, коли необхідно застосувати якусь функцію до компонентів вкладеного списку (буква “r” в назві rapply() означає recursively – “рекурсивно”).

# Функція користувача, яка додає ! до елементу
# об’єкту x якщо цей елемент є текстовим виразом,
# або додає 1 якщо цей елемент є числом:
myFun <- function(x) {
  if (is.character(x)) {
    return(paste(x, "!", sep=""))
    }
  else {
    return(x + 1)
    }
  }

# Приклад вкладеного списку:
l <- list(a=list(a1="Boo", b1=2, c1="Eeek"),
          b=3, c="Yikes",
          d=list(a2=1, b2=list(a="Hey", b3=5)))

# Рекурсивне застосування функції myFun до списку l:
rapply(l, myFun)

    a.a1     a.b1     a.c1        b        c     d.a2   d.b2.a 
  "Boo!"      "3"  "Eeek!"      "4" "Yikes!"      "2"   "Hey!" 
 d.b2.b3 
     "6" 

# Якщо необхідно повернути результат у вигляді вкладеного списку,
# можна скористатися аргументом how=“replace”. У цьому випадку # вихідні значення в списку 1 будуть замінені на нові:
rapply(l, myFun, how="replace")

$a
$a$a1
[1] "Boo!"

$a$b1
[1] 3

$a$c1
[1] "Eeek!"


$b
[1] 4

$c
[1] "Yikes!"

$d
$d$a2
[1] 2

$d$b2
$d$b2$a
[1] "Hey!"

$d$b2$b3
[1] 6

tapply()

       Функція tapply() використовується у випадках, коли необхідно застосувати будь-яку функцію fun() до окремих груп елементів вектору x, що заданий відповідно до рівнів будь-якого фактору group:

tapply(x, group, fun, ...)

       Наприклад, у такому фрагменті коду функція sample() використовується для створення двох випадкових вибірок: з 50 значень цілих чисел від 1 до 4 і пов’язаних із ними міток чотирьох груп A-D. Функція tapply() підраховує суми х для кожного зі значень фактору:

x <- sample(1:4, size=50, replace=T)
gr <- as.factor(sample(c("A", "B", "C", "D"), size=50, replace=T))
tapply(x, gr, sum)

 A  B  C  D 
29 21 44 33 

by()

       Функція by() є свого роду аналогом функції tapply(), з тією різницею, що вона застосовується для таблиць. Таблиця data поділяється відповідно до заданого стовпця-фактору group на підмножину підтаблиць і для обробки кожної такої частини визначається функція fun():

by(data, group, fun, ...)

       Розглянемо приклад із таблицею Молюски, яка представлена в розділі 1.3, і розрахуємо середні значення довжини раковини ZMlength і чисельності інфузорій CAnumber для кожного з трьох обстежених озер Lake(Баторіно, Мястро й Нарочь):

outer()

       Функція outer() дозволяє виконати комбінаторну операцію fun() над елементами двох масивів або векторів x і у, не вдаючись до явного використання “подвійного” циклу:

outer(x, y, fun="*",...)

       За умовчанням здійснюється операція попарного перемноження

x <- 1:5
y <- 1:5
outer(x,y)

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    2    4    6    8   10
[3,]    3    6    9   12   15
[4,]    4    8   12   16   20
[5,]    5   10   15   20   25

       Використовуючи, наприклад, функцію outer() разом із функцією paste(), можна згенерувати всі можливі попарні комбінації “зв’язок” елементів символьного й цілочисельного векторів:

x <- c("A", "B", "C", "D")
y <- 1:10
outer(x, y, paste, sep="")

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] "A1" "A2" "A3" "A4" "A5" "A6" "A7" "A8" "A9" "A10"
[2,] "B1" "B2" "B3" "B4" "B5" "B6" "B7" "B8" "B9" "B10"
[3,] "C1" "C2" "C3" "C4" "C5" "C6" "C7" "C8" "C9" "C10"
[4,] "D1" "D2" "D3" "D4" "D5" "D6" "D7" "D8" "D9" "D10"