В моем продолжении чтения Think Bayes следующая проблема, которую необходимо решить, — это проблема локомотива, которая определяется следующим образом:
Железная дорога нумерует свои локомотивы в порядке 1..N.
Однажды вы видите локомотив с номером 60. Оцените, сколько локомотивов у железной дороги.
Интересная вещь в этом вопросе заключается в том, что изначально кажется, что у нас недостаточно информации, чтобы найти какой-либо ответ. Тем не менее, мы можем получить оценку, если мы придумали до работы с.
Простейший принцип состоит в том, чтобы предположить, что существует один железнодорожный оператор с, например, от 1 до 1000 железных дорог с равной вероятностью каждого размера.
Затем мы можем написать такой же код, как и в случае с игрой в кости, чтобы обновить предыдущее на основе поездов, которые мы видели.
Сначала мы создадим фрейм данных, который фиксирует произведение «количества локомотивов» и наблюдений локомотивов, которые мы видели (в данном случае мы видели только один локомотив с номером «60» 🙂
library(dplyr)
possibleValues = 1:1000
observations = c(60)
l = list(value = possibleValues, observation = observations)
df = expand.grid(l)
> df %>% head()
value observation
1 1 60
2 2 60
3 3 60
4 4 60
5 5 60
6 6 60Затем мы хотим добавить столбец, который представляет вероятность того, что наблюдаемый локомотив мог прибыть из определенного парка. Если количество железных дорог меньше 60, то вероятность 0, в противном случае имеем 1 / numberOfRailroadsInFleet:
prior = 1 / length(possibleValues)
df = df %>% mutate(score = ifelse(value < observation, 0, 1/value))
> df %>% sample_n(10)
value observation score
179 179 60 0.005586592
1001 1001 60 0.000999001
400 400 60 0.002500000
438 438 60 0.002283105
667 667 60 0.001499250
661 661 60 0.001512859
284 284 60 0.003521127
233 233 60 0.004291845
917 917 60 0.001090513
173 173 60 0.005780347Чтобы найти вероятность каждого размера флота, мы пишем следующий код:
weightedDf = df %>%
group_by(value) %>%
summarise(aggScore = prior * prod(score)) %>%
ungroup() %>%
mutate(weighted = aggScore / sum(aggScore))
> weightedDf %>% sample_n(10)
Source: local data frame [10 x 3]
value aggScore weighted
1 906 1.102650e-06 0.0003909489
2 262 3.812981e-06 0.0013519072
3 994 1.005031e-06 0.0003563377
4 669 1.493275e-06 0.0005294465
5 806 1.239455e-06 0.0004394537
6 673 1.484400e-06 0.0005262997
7 416 2.401445e-06 0.0008514416
8 624 1.600963e-06 0.0005676277
9 40 0.000000e+00 0.0000000000
10 248 4.028230e-06 0.0014282246Давайте построим фрейм данных, чтобы увидеть, как изменяется вероятность для каждого размера парка:
library(ggplot2)
ggplot(aes(x = value, y = weighted), data = weightedDf) +
geom_line(color="dark blue")
Наиболее вероятный выбор — размер флота 60, основанный на этой диаграмме, но альтернативой было бы найти среднее значение апостериорного значения, которое мы можем сделать так:
> weightedDf %>% mutate(mean = value * weighted) %>% select(mean) %>% sum()
[1] 333.6561Теперь давайте создадим функцию со всем этим кодом, чтобы мы могли поиграть с некоторыми другими априорами и наблюдениями:
meanOfPosterior = function(values, observations) {
l = list(value = values, observation = observations)
df = expand.grid(l) %>% mutate(score = ifelse(value < observation, 0, 1/value))
prior = 1 / length(possibleValues)
weightedDf = df %>%
group_by(value) %>%
summarise(aggScore = prior * prod(score)) %>%
ungroup() %>%
mutate(weighted = aggScore / sum(aggScore))
return (weightedDf %>% mutate(mean = value * weighted) %>% select(mean) %>% sum())
}Если мы обновим наши наблюдаемые железные дороги на номера 60, 30 и 90, мы получим следующие средства для постеров, предполагающих разные приоры:
> meanOfPosterior(1:500, c(60, 30, 90))
[1] 151.8496
> meanOfPosterior(1:1000, c(60, 30, 90))
[1] 164.3056
> meanOfPosterior(1:2000, c(60, 30, 90))
[1] 171.3382На данный момент функция предполагает, что мы всегда хотим иметь одинаковый априор, т. Е. Каждый вариант имеет равную возможность быть выбранным, но мы могли бы изменить априор, чтобы увидеть, как различные допущения влияют на апостериор.
Мы можем реорганизовать функцию для получения значений и априоров вместо вычисления априоров в функции:
meanOfPosterior = function(values, priors, observations) {
priorDf = data.frame(value = values, prior = priors)
l = list(value = priorDf$value, observation = observations)
df = merge(expand.grid(l), priorDf, by.x = "value", by.y = "value") %>%
mutate(score = ifelse(value < observation, 0, 1 / value))
df %>%
group_by(value) %>%
summarise(aggScore = max(prior) * prod(score)) %>%
ungroup() %>%
mutate(weighted = aggScore / sum(aggScore)) %>%
mutate(mean = value * weighted) %>%
select(mean) %>%
sum()
}Теперь давайте проверим, получаем ли мы одинаковые апостериорные средства для одинаковых приоров:
> meanOfPosterior(1:500, 1/length(1:500), c(60, 30, 90))
[1] 151.8496
> meanOfPosterior(1:1000, 1/length(1:1000), c(60, 30, 90))
[1] 164.3056
> meanOfPosterior(1:2000, 1/length(1:2000), c(60, 30, 90))
[1] 171.3382Теперь, если вместо единообразного априора мы будем использовать степенной закон, в котором предполагается, что меньшие флоты более вероятны:
> meanOfPosterior(1:500, sapply(1:500, function(x) x ** -1), c(60, 30, 90))
[1] 130.7085
> meanOfPosterior(1:1000, sapply(1:1000, function(x) x ** -1), c(60, 30, 90))
[1] 133.2752
> meanOfPosterior(1:2000, sapply(1:2000, function(x) x ** -1), c(60, 30, 90))
[1] 133.9975
> meanOfPosterior(1:5000, sapply(1:5000, function(x) x ** -1), c(60, 30, 90))
[1] 134.212
> meanOfPosterior(1:10000, sapply(1:10000, function(x) x ** -1), c(60, 30, 90))
[1] 134.2435Теперь мы получаем очень похожие апостериорные средства, которые сходятся на 134, и это наше лучшее предсказание.