Статьи

R: замена циклов на фреймы данных

В своем последнем сообщении в блоге я показал, как получить апостериорные вероятности для проблемы костей Байеса Think :

Предположим, у меня есть коробка с кубиками, которая содержит 4-сторонний кубик, 6-сторонний кубик, 8-сторонний кубик, 12-сторонний кубик и 20-сторонний кубик. Если вы когда-нибудь играли в Dungeons & Dragons, вы знаете, о чем я говорю.

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

Напомним, это было мое окончательное решение:

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
likelihoods = function(names, observations) {
  scores = rep(1.0 / length(names), length(names)) 
  names(scores) = names
  
  for(name in names) {
      for(observation in observations) {
        if(name < observation) {
          scores[paste(name)]  = 0
        } else {
          scores[paste(name)] = scores[paste(name)] *  (1.0 / name)
        }       
      }
    
  return(scores)
}
  
dice = c(4,6,8,12,20)
l1 = likelihoods(dice, c(6))
  
> l1 / sum(l1)
        4         6         8        12        20
0.0000000 0.3921569 0.2941176 0.1960784 0.1176471

Хотя это работает, мы вложили в циклы, которые не очень идиоматичны, поэтому давайте попробуем избавиться от них.

Первое, что мы хотим сделать, это вернуть фрейм данных, а не вектор, поэтому мы настраиваем первые две строки следующим образом:

1
2
scores = rep(1.0 / length(names), length(names)) 
df = data.frame(score = scores, name = names)

Затем мы можем избавиться от внутреннего цикла for и заменить его вызовом ifelse, заключенным в вызов dplyr mutate:

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
library(dplyr)
likelihoods2 = function(names, observations) {
  scores = rep(1.0 / length(names), length(names)) 
  df = data.frame(score = scores, name = names)
  
  for(observation in observations) {
    df = df %>% mutate(score = ifelse(name < observation, 0, score * (1.0 / name)) )
  }
  
  return(df)
}
  
dice = c(4,6,8,12,20)
l1 = likelihoods2(dice, c(6))
  
> l1
       score name
1 0.00000000    4
2 0.03333333    6
3 0.02500000    8
4 0.01666667   12
5 0.01000000   20

Наконец, мы подведем итоги, чтобы они были относительно взвешенными:

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
likelihoods2 = function(names, observations) {
  scores = rep(1.0 / length(names), length(names)) 
  df = data.frame(score = scores, name = names)
  
  for(observation in observations) {
    df = df %>% mutate(score = ifelse(name < observation, 0, score * (1.0 / name)) )
  }
  
  return(df %>% mutate(weighted = score / sum(score)) %>% select(name, weighted))
}
  
dice = c(4,6,8,12,20)
l1 = likelihoods2(dice, c(6))
  
> l1
  name  weighted
1    4 0.0000000
2    6 0.3921569
3    8 0.2941176
4   12 0.1960784
5   20 0.1176471

Теперь мы подошли к циклу. Избавиться от этого немного сложнее. Сначала мы создадим фрейм данных, содержащий строку для каждой пары (наблюдения, кости), имитирующую вложенные циклы:

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
likelihoods3 = function(names, observations) {
  l = list(observation = observations, roll = names)
  obsDf = do.call(expand.grid,l) %>%
    mutate(likelihood = 1.0 / roll,
           score = ifelse(roll < observation, 0, likelihood))  
  
  return(obsDf)
}
  
dice = c(4,6,8,12,20)
l1 = likelihoods3(dice, c(6))
  
> l1
  observation roll likelihood      score
1           6    4 0.25000000 0.00000000
2           6    6 0.16666667 0.16666667
3           6    8 0.12500000 0.12500000
4           6   12 0.08333333 0.08333333
5           6   20 0.05000000 0.05000000
  
l2 = likelihoods3(dice, c(6, 4, 8, 7, 7, 2))
> l2
   observation roll likelihood      score
1            6    4 0.25000000 0.00000000
2            4    4 0.25000000 0.25000000
3            8    4 0.25000000 0.00000000
4            7    4 0.25000000 0.00000000
5            7    4 0.25000000 0.00000000
6            2    4 0.25000000 0.25000000
7            6    6 0.16666667 0.16666667
8            4    6 0.16666667 0.16666667
9            8    6 0.16666667 0.00000000
10           7    6 0.16666667 0.00000000
11           7    6 0.16666667 0.00000000
12           2    6 0.16666667 0.16666667
13           6    8 0.12500000 0.12500000
14           4    8 0.12500000 0.12500000
15           8    8 0.12500000 0.12500000
16           7    8 0.12500000 0.12500000
17           7    8 0.12500000 0.12500000
18           2    8 0.12500000 0.12500000
19           6   12 0.08333333 0.08333333
20           4   12 0.08333333 0.08333333
21           8   12 0.08333333 0.08333333
22           7   12 0.08333333 0.08333333
23           7   12 0.08333333 0.08333333
24           2   12 0.08333333 0.08333333
25           6   20 0.05000000 0.05000000
26           4   20 0.05000000 0.05000000
27           8   20 0.05000000 0.05000000
28           7   20 0.05000000 0.05000000
29           7   20 0.05000000 0.05000000
30           2   20 0.05000000 0.05000000

Теперь нам нужно перебрать фрейм данных, сгруппировав его по «рулонам», чтобы мы получили по одной строке для каждого.

Мы добавим новый столбец, в котором хранится апостериорная вероятность для каждой кости. Это будет рассчитано путем умножения предыдущей вероятности на произведение записей «оценка».

Вот как выглядит наша новая функция правдоподобия:

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
likelihoods3 = function(names, observations) {
  l = list(observation = observations, roll = names)
  obsDf = do.call(expand.grid,l) %>%
    mutate(likelihood = 1.0 / roll,
           score = ifelse(roll < observation, 0, likelihood))  
  
  return (obsDf %>%
    group_by(roll) %>%
    summarise(s = (1.0/length(names)) * prod(score) ) %>%
    ungroup() %>%
    mutate(weighted = s / sum(s)) %>%
    select(roll, weighted))
}
  
l1 = likelihoods3(dice, c(6))
> l1
Source: local data frame [5 x 2]
  
  roll  weighted
1    4 0.0000000
2    6 0.3921569
3    8 0.2941176
4   12 0.1960784
5   20 0.1176471
  
l2 = likelihoods3(dice, c(6, 4, 8, 7, 7, 2))
> l2
Source: local data frame [5 x 2]
  
  roll    weighted
1    4 0.000000000
2    6 0.000000000
3    8 0.915845272
4   12 0.080403426
5   20 0.003751302

Теперь мы получили тот же результат, что и с нашими вложенными циклами for, поэтому я думаю, что рефакторинг прошел успешно.

Ссылка: Р: Замена циклов фреймами данных от нашего партнера по JCG Марка Нидхэма в блоге Марка Нидхэма .