Краткое доказательство парадокса Монти Холла в р


Это попытка в краткой доказательство парадокса Монти Холла по теории вероятностей и статистики в р.

Для тех, кто незнаком, сценарий таков:

  • Существует кандидатом на игровое шоу. Этот участник выбирает между тремя дверями (А, B и C), за одним из которых является новый автомобиль.
  • Хозяин (Монти Холла) показывает, что за одной из двух других дверей.
  • Участнику предлагается либо с их же дверь или выключатель.
  • Вопрос: какой выбор предлагает самую высокую вероятность успеха (приз автомобиль)?

Для выполнения этой задачи, есть некоторые предпосылки:

  • Вероятность того, что автомобиль находится за дверью, прежде чем что-то случится, равномерно распределена.
  • Хозяин не откроет машину на первом открытии двери (это то, что компенсирует вероятность).

Предположим (произвольно), что участник выбирает дверь А. Если автомобиль находится позади двери, хозяин открыть дверь B или C с равным (50/50) вероятность - это интуитивное. Однако, если автомобиль находится позади двери Б, узел должен открыть дверь с целью предотвращения выявления автомобилей, таким образом, вероятность открытия дверей в нулю. Это предвидение местоположения автомобиля-это то, что искажает проблему.

Результатом является то, что участник заканчивает с 66% шансов на победу путем переключения выбранной ими двери, и 33%, если они придерживаться их нынешний выбор. В ниже код, он является обязательным для участника, чтобы переключиться - это было для повышения вычислительной эффективности:

doors <- c("A", "B", "C")
wins <- 0
n <- 10000

for (i in 1:n) {
  # Choose which door the car is behind
  car <- sample(doors)[1]

  # The contestant chooses a door
  contestant <- sample(doors)[1]

  # Monty opens a door that the car is not behind
  monty <- sample(doors[which((doors != car) & (doors != contestant))])[1]

  # Force the contestant to switch
  contestant <- sample(doors[which((doors != monty) & (doors != contestant))])[1]

  # Count the wins
  if(contestant == car){wins <- wins+1}
}

(wins/n)


193
4
задан 4 февраля 2018 в 02:02 Источник Поделиться
Комментарии
1 ответ

Фундаментальные идеи вашего кода являются твердыми, однако есть несколько областей, которые можно улучшить. Для начала, когда вы называете sample и вы только хотите один номер вернули, вы можете воспользоваться size аргумент. Для векторов только размер 3, как в вашем случае, особой разницы нет:

library(microbenchmark)
microbenchmark(sample(doors, size = 1),
sample(doors)[1], times = 10^4, unit = "relative")
Unit: relative
expr min lq mean median uq max neval
sample(doors, size = 1) 1.000000 1.00000 1.000000 1.000000 1.000000 1.000000 10000
sample(doors)[1] 1.104667 1.09216 1.111044 1.086809 1.082041 1.072216 10000

Впрочем, как и размер вашего входного вектора увеличивается, наблюдается заметное различие в методах.

doorsHuges <- rep(doors, 10^6)
length(doorsHuges)
[1] 3000000

microbenchmark(sample(doorsHuges, size = 1),
sample(doorsHuges)[1], unit = "relative")
Unit: relative
expr min lq mean median uq max neval
sample(doorsHuges, size = 1) 1.00 1.00 1.00 1.00 1.000 1.00 100
sample(doorsHuges)[1] 65693.19 31524.98 10609.52 9126.07 6347.235 4338.52 100

Далее, мы можем значительно улучшить нашу эффективность за счет того, что sample может выполнять много повторений за один раз при помощи replace аргумент. Например, вместо того, чтобы использовать цикл for и генерирует один сценарий, мы можем сделать следующее:

sample(doors, n, replace = TRUE)

Выигрыш в эффективности огромна:

microbenchmark(forLoop = for(i in 1:1000){sample(doors, 1)},
usingRep = sample(doors, 1000, replace = TRUE), unit = "relative")
Unit: relative
expr min lq mean median uq max neval
forLoop 211.0778 204.1896 127.7005 199.9415 191.6994 185.1182 100
usingRep 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 100

Теперь, давайте посмотрим на наши типы данных. Поскольку мы главным образом связаны с вероятностями, фактические значения наш главный вектор не имеет значения. Поскольку это так, мы должны смотреть в использовании integers когда мы можем. Мы можем добиться этого, используя seq_along(doors) или просто length(doors). Наблюдать:

intDoors <- seq_along(doors)
microbenchmark(intDoors == 2, doors == "B", times = 10^4, unit = "relative")
Unit: relative
expr min lq mean median uq max neval
intDoors == 2 1.000000 1.000000 1.00000 1.000000 1.000000 1.000000 10000
doors == "B" 1.449541 1.422819 1.35682 1.396226 1.390533 1.053855 10000

А мы на целые числа, which это удивительная функция, которая является очень эффективной и очень интуитивным. Однако, когда мы индексируем, это совершенно лишнее. Мы можем вместо этого использовать логические подмножества. Это чище и менее эффективным. Наблюдать:

 mySamp <- sample(10^6, 10^6)
boolTest <- rep(FALSE, 10^6)
boolTest[mySamp] <- TRUE
testIndex <- 1:10^6

microbenchmark(testIndex[which(boolTest)], testIndex[boolTest], unit = "relative")
Unit: relative
expr min lq mean median uq max neval
testIndex[which(boolTest)] 1.533577 1.350712 1.185725 1.36155 1.386107 0.1012212 100
testIndex[boolTest] 1.000000 1.000000 1.000000 1.00000 1.000000 1.0000000 100

Для того, чтобы решить monty вектор и вторая итерация contestantмы можем реализовать vapply вместе с length чтобы избежать вызовов sample где длина 1.

monty <- vapply(1:numReps, function(x) {
mySet <- intDoors[intDoors != car[x] & intDoors != contestant[x]]
if (length(mySet) > 1) sample(mySet, 1) else mySet}, 1L)

И, наконец, вместо повышения winмы можем сделать все это за один раз при помощи sum(contestant == car). Это имеет преимущество векторизации и внутреннего принуждения с логические для целых чисел. Положить все это вместе, мы получаем:

funImproved <- function(numReps) {
car <- sample(intDoors, numReps, replace = TRUE)
contestant <- sample(intDoors, numReps, replace = TRUE)

monty <- vapply(1:numReps, function(x) {
mySet <- intDoors[intDoors != car[x] & intDoors != contestant[x]]
if (length(mySet) > 1) sample(mySet, 1) else mySet}, 1L)

contestant <- vapply(1:numReps, function(x) {
mySet <- intDoors[intDoors != monty[x] & intDoors != contestant[x]]
if (length(mySet) > 1) sample(mySet, 1) else mySet}, 1L)

wins <- sum(contestant == car)
(wins/numReps)
}

Функция ОП является:

funOP <- function(numReps) {
wins <- 0

for (i in 1:numReps) {
# Choose which door the car is behind
car <- sample(doors)[1]

# The contestant chooses a door
contestant <- sample(doors)[1]

# Monty opens a door that the car is not behind
monty <- sample(doors[which((doors != car) & (doors != contestant))])[1]

# Force the contestant to switch
contestant <- sample(doors[which((doors != monty) & (doors != contestant))])[1]

# Count the wins
if(contestant == car){wins <- wins+1}
}

(wins/numReps)
}

И окончательное сравнение видит улучшение около 4,5 х.

microbenchmark(funOP(n), funImproved(n),
times = 50, unit = "relative")
Unit: relative
expr min lq mean median uq max neval
funOP(n) 5.060335 4.518107 4.54938 4.539024 4.183145 7.831481 50
funImproved(n) 1.000000 1.000000 1.00000 1.000000 1.000000 1.000000 50

Результаты очень похожи, а также:

## the average difference over 100 trials when n = 10^5
mean(replicate(100, funOP(10^5) - funImproved(10^5)))
[1] -0.0001978

В качестве последнего замечания, когда происходит выполнение какой-либо процедуры, которая требует проверки, это хорошая идея, чтобы использовать set.seed так что вы можете проверить свои результаты предсказуемо. Например:

sapply(1:100, function(x) {
set.seed(42)
funImproved(10^3)
})
[1] 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668
[16] 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668
[31] 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668
[46] 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668
[61] 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668
[76] 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668
[91] 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668 0.668

9
ответ дан 4 февраля 2018 в 09:02 Источник Поделиться