В своем последнем сообщении в блоге я показал, как получить апостериорные вероятности для проблемы костей Байеса 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 Марка Нидхэма в блоге Марка Нидхэма . |