В своем последнем сообщении в блоге я показал, как получить апостериорные вероятности для проблемы костей Байеса 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 200.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 name1 0.00000000 42 0.03333333 63 0.02500000 84 0.01666667 125 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 weighted1 4 0.00000002 6 0.39215693 8 0.29411764 12 0.19607845 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 score1 6 4 0.25000000 0.000000002 6 6 0.16666667 0.166666673 6 8 0.12500000 0.125000004 6 12 0.08333333 0.083333335 6 20 0.05000000 0.05000000 l2 = likelihoods3(dice, c(6, 4, 8, 7, 7, 2))> l2 observation roll likelihood score1 6 4 0.25000000 0.000000002 4 4 0.25000000 0.250000003 8 4 0.25000000 0.000000004 7 4 0.25000000 0.000000005 7 4 0.25000000 0.000000006 2 4 0.25000000 0.250000007 6 6 0.16666667 0.166666678 4 6 0.16666667 0.166666679 8 6 0.16666667 0.0000000010 7 6 0.16666667 0.0000000011 7 6 0.16666667 0.0000000012 2 6 0.16666667 0.1666666713 6 8 0.12500000 0.1250000014 4 8 0.12500000 0.1250000015 8 8 0.12500000 0.1250000016 7 8 0.12500000 0.1250000017 7 8 0.12500000 0.1250000018 2 8 0.12500000 0.1250000019 6 12 0.08333333 0.0833333320 4 12 0.08333333 0.0833333321 8 12 0.08333333 0.0833333322 7 12 0.08333333 0.0833333323 7 12 0.08333333 0.0833333324 2 12 0.08333333 0.0833333325 6 20 0.05000000 0.0500000026 4 20 0.05000000 0.0500000027 8 20 0.05000000 0.0500000028 7 20 0.05000000 0.0500000029 7 20 0.05000000 0.0500000030 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))> l1Source: local data frame [5 x 2] roll weighted1 4 0.00000002 6 0.39215693 8 0.29411764 12 0.19607845 20 0.1176471 l2 = likelihoods3(dice, c(6, 4, 8, 7, 7, 2))> l2Source: local data frame [5 x 2] roll weighted1 4 0.0000000002 6 0.0000000003 8 0.9158452724 12 0.0804034265 20 0.003751302 |
Теперь мы получили тот же результат, что и с нашими вложенными циклами for, поэтому я думаю, что рефакторинг прошел успешно.
| Ссылка: | Р: Замена циклов фреймами данных от нашего партнера по JCG Марка Нидхэма в блоге Марка Нидхэма . |