Функциональное программирование в R

Элементы функционального программирования

R сочетает в себе несколько парадигм программирования. Если попытаться выделить какую-то главную, то пожалуй - это элементы функционального проаграммирования. Эта тема относится к продвинутому уровню. Это нормально чувствовать себя на этом уровне некомфортно и неуверенно

Объектно-ориентированные системы

Когда мы говорим об объектах языка R таких как вектор, функция и т.д., то это не то же самое когда речь идет об объектно-ориентированном программировании.

Что входит в понятие объектино-ориентированного программирования?:
1. Классы

2. Поля классов

3. Методы, определенные для классов

Классы

В R можно создавать классы. Стандартно их сразу 3 в R:
* S3
- нет формальной деклараци класса
- функция может иметь разное поведение в зависимости от класса (method dispatch)
- такие функции называются generic
Вторые два класса встречаются значительно реже, прочитать про них подробнее можно в книге Advanced R * S4
- строгое определеие класса и его полей
- больше возможностей для поведения методов
* Reference classes

Generic функции

Generic функции следует разобрать подробнее, т.к. всречаются повсеместно. Можно проверить является ли функция print generic функцией

length(methods(print))

[1] 217

Если мы получили много вариантов, значит функция может себя по-разному вести в зависимости от того, чтобы приходит ей в качестве аргумента. Какие могут быть примеры?

Если берем print(x), то если:
- x - дата-фрейм, то вызывается print.data.frame(x);
- если x - функция, то print.function(x) и так далее
- если ни один из методов не подходит, то print.default(x)

Функции без сторонних эффектов

В R нет указателей на объекты, все объекты передаются “по значению” (хотя есть нюансы!)

f <- function(k) {
  k <- k+1
  a <- a + k*2
  a
}
k <- 5
f(k)

[1]  13  14  54 111

Примеры функций, отражающие функциональное программирование в языке R

Функция **replicate*
В задачах моделирования случайных процессов приходится встречаться с запросом повторить один и тот же выхов несколько раз, при этом каждый новый вызов может иметь новый результат (зависящий, например, датчика случайных чисел)

Фунции семейтсва apply apply, sapply. lapply, tapply

mapply - многомерная версия sapply. Она приниает на вход функцию и произвольное количество списков аргументов, которые необходиом перебрать. Например:

mapply(seq, from = 1:4, to = 2:5, by = 1/(1 + 1:4))

[[1]]
[1] 1.0 1.5 2.0

[[2]]
[1] 2.000000 2.333333 2.666667
[4] 3.000000

[[3]]
[1] 3.00 3.25 3.50 3.75 4.00

[[4]]
[1] 4.0 4.2 4.4 4.6 4.8 5.0

Получаем тоже самое, что получили бы при последовательном вызове функции seq для каждого первого элемента аргументов и далее каждого n-го элемента аргументов

outer

перебор всех возможных комбинаций аргументов.
Допустим, хочу применить paste0 ко всем возможным аргументам массива LETTERS в нижнем и в ВЕРХНЕМ регистре. Тогда можно воспользоваться функцией outer

m <- outer(letters, LETTERS, paste0)

Теперь можно увидеть, что m - это матрица размером 26*26. При этом диагональные элементы будут содержать в себе одинаковые буквы

diag(m)

 [1] "aA" "bB" "cC" "dD" "eE"
 [6] "fF" "gG" "hH" "iI" "jJ"
[11] "kK" "lL" "mM" "nN" "oO"
[16] "pP" "qQ" "rR" "sS" "tT"
[21] "uU" "vV" "wW" "xX" "yY"
[26] "zZ"

m[1:5,1:4]

     [,1] [,2] [,3] [,4]
[1,] "aA" "aB" "aC" "aD"
[2,] "bA" "bB" "bC" "bD"
[3,] "cA" "cB" "cC" "cD"
[4,] "dA" "dB" "dC" "dD"
[5,] "eA" "eB" "eC" "eD"

Vectorize

Не все функции векторизированы по умолчанию, то их можно сделать такими с помощью функции Vectorize

do.call

Последняя функция из списка функционального прогрммирования. Она осуществляет вызов некоторой функции на списке аргументов. Например, пришел дата-фрейм, но он сегментирован на 3 отдельных файла, которые нужно объединить между собой. Можно их соединить вместе с помощью rbind, а можно с помощью do.call

df1 <- data.frame(ld = 1:2, value = rnorm(2))
df2 <- data.frame(ld = 3:4, value = runif(2))
df3 <- data.frame(ld = 222, value = 7)  
do.call(rbind, list(df1,df2, df3))

Но зачем тогда вообще do.call, если в этом примере можно было обойтись без него просто точно так же склеив дата фреймы через rbind.
Функция нужна, когда мы не знаем сколько именно объектов нужно будет склеить

do.call(rbind,lapply(list.files(), function(file) read.csv(file)))

cannot open file 'ASUS': Permission deniedError in file(file, "rt") : cannot open the connection

Глосаррий

• S3, S4, reference class
  
• generic function
  
• copy-on-modify semantics
  
• ?replicate
  
• ?mapply
  
• ?outer
  
• Vectorize
  
• do.call

Моделирование траектории случайного процесса

Упражнение может иметь реальное применение в физике элементарных частиц.

# Random walk with absorption
simulate_walk <- function(lower = -10, upper = 10, n_max = 200, p = 1e-3) {
   # на каждом шаге есть позиция, начинаем из центра интервала
  current_position <- (lower + upper) / 2
  for (i in 1:n_max) {
    #сначала определяем вероятность поглощения частицы
    is_absorbed <- rbinom(1, 1, p)
    #если произошло - симуляция заканчивается
    if (is_absorbed) return(list(status = "Absorbed", 
                                 position = current_position, 
                                 steps = i))
    #если поглощение не произошло - совершаем 1 случайный шаг
    current_position <- current_position + rnorm(1)
    if (current_position < lower) return(list(status = "Left breach", 
                                              position = current_position, 
                                              steps = i))
    if (current_position > upper) return(list(status = "Right breach", 
                                              position = current_position, 
                                              steps = i))
  }
  #возвращаем такое сообщение при достижении максимального количества шагов
  return(list(status = "Max steps reached", 
              position = current_position,#позиция, на которой завершилась симуляция
              steps = n_max)) #количество шагов, на которой завершилась симуляция
}

Есть одномерная область от -10 до 10, в которой может находится элементарная частица. Частица может прыгать либо вправо, либо влево и на каждом шаге есть вероятность того, что она поглотится (p = 1e-3). Общее количество переходов ограничиваем 200 шагами (n_max = 200).
Почему использовали цикл for, хотя для R - это один из медленных способов, к которому нужно относится с осторожностью. Дело в том, что на каждом из шагов симуляция может быть прервана и закончится - каждое определение перехода это какая-то сложная вычислительная заадача. Если генерировать такую ситуацию через векторизацию, то мы генерируем сразу 200 переходов (потенциальный максимум), тогда как цикл может прерваться уже на первой итерации в случае поглощения частицы. Т.е. в такой ситуации цикл for позволяет сэкономить те переходы, которые выпадают за границы симуляции

Теперь можно с помощью функции replicate повторить n-ое число таких симуляций

# Simulate results
result <- replicate(1000, simulate_walk(), simplify = FALSE)
result <- data.frame(
  status = sapply(result, function(x) x$status),
  position = sapply(result, function(x) x$position),
  steps = sapply(result, function(x) x$steps)
)

Приводим результат в том виде, который нам удобен. В данной конкретной ситуации превращаем его в дата-фрейм

Можно посмотреть, что примерно получилось

str(result)

'data.frame':   1000 obs. of  3 variables:
 $ status  : chr  "Left breach" "Absorbed" "Left breach" "Left breach" ...
 $ position: num  -10.17 -6.95 -10.65 -11.12 10.04 ...
 $ steps   : num  56 27 84 167 20 87 100 79 26 121 ...

Теперь, когда у нас есть какая-то разумная статистика, то можно ее определеить, например, с помощью функции tapply

tapply(result$position, result$status, length)

         Absorbed       Left breach Max steps reached      Right breach 
               99               389               120               392 

tapply(result$steps, result$status, mean)

         Absorbed       Left breach Max steps reached      Right breach 
         77.54545          79.48329         200.00000          82.80612 

Больший интерес представляют задачи с блужданием по плоскости, то есть в размерности 2.
Возьмите написанную мной функцию и измените её так, чтобы блуждание начиналось в центре координат (0, 0), а все переходы по координатам x и y были бы независимы и имели стандартное нормальное распределение. Если вас пугают эти слова, то это то же самое, что делал я, только отдельно по x и по y.

• Процесс обрывается в момент выхода за границу круга с центром в (0, 0) и радиусом 6
  
• Вероятность поглощения на каждом шаге равна 0.01
  
• Максимальное количество шагов — 100
  
• Расстояние, разумеется, евклидово: расстояние от точки (x, y) до (0, 0) есть
  
• Один шаг процесса подразумевает изменение обеих координат одновременно!

В ответе укажите вероятность вылета частицы в процентах, с точностью до целых процентов, в виде XX (например, 14, без указания значка процентов). Вероятность должна получиться больше 50%.

simulate_walk <- function(n_max = 100, p = 1e-2) {
  x<-0
  y<-0
  current_position <- 0 
  for (i in 1:n_max) {
    is_absorbed <- rbinom(1, 1, p)
    if (is_absorbed) return(list(status = "Absorbed", 
                                 position = current_position, 
                                 steps = i))
    x<-x+rnorm(1)
    y<-y+rnorm(1)
    current_position <- sqrt(x**2 + y**2)
    if (current_position > 6) return(list(status = "breach", 
                                              position = current_position, 
                                              steps = i))
  }
  return(list(status = "Max steps reached", 
              position = current_position,
              steps = n_max))
}

# Simulate results
result <- replicate(100000, simulate_walk(), simplify = FALSE)
result <- data.frame(
  status = sapply(result, function(x) x$status),
  position = sapply(result, function(x) x$position),
  steps = sapply(result, function(x) x$steps)
)

# Inspect results
tapply(result$position, result$status, length)

         Absorbed            breach 
            18934             80997 
Max steps reached 
               69 

a

[1] 81.05293

Общий глоссарий для этого урока:

S3, S4, reference classes
Generic function
Copy-on-modify semantics
?replicate
?mapply
?outer
?Vectorize
?do.call