Игра Монти Холл в R с базовыми функциями

Question

Просто для удовольствия и для обучения R, я попытался доказать правило Игры Монти Холла (изменение вашего выбора после того, как открылись одни ворота, дает вам больше шансов на победу), я создал этот воспроизводимый код (объяснение каждого шага находится в пределах код):

## First I set the seed

set.seed(4)

## Then I modelize the presence of the prize as a random variable between gates 1,2,3


randomgates <- ceiling(runif(10000, min = 0, max = 3))

## so do I with the random choice.

randomchoice <- ceiling(runif(10000, min = 0, max = 3))

## As the opening of a gate is dependent from the gate you chose (the gate you chose cannot be opened)
## I modelize the opening of the gate as a variable which cannot be equal to the choice.

options <- c(1:3)

randomopen <- rep(1,10000)

for (i in 1:length(randomgates)) {
  realoptions <- options[options != randomchoice[i]]
  randomopen[i] <- realoptions[ceiling(runif(1,min = 0, max = 2))]
}

##Just to make data more easy to handle, I make a dataset

dataset <- cbind(randomgates, randomchoice, randomopen)

## Then I creat a dataset which only keeps the realization of the games in which we carry on (
## the opened gate wasn't the one with the price within)

steptwo <- dataset[randomopen != randomgates,]

## The next step is just to check if the probability of carry on is 2/3, which indeed is

carryon <- randomopen != randomgates

sum(carryon)/length(randomgates) 

## I format the dataset as a data frame

steptwo <- as.data.frame(steptwo)

## Now we check what happens if we hold our initial choice when game carries on

prizesholding <- steptwo$randomgates == steptwo$randomchoice

sum(prizesholding)

## creating a vector of changing option, dependant on the opened gate, in the dataset that
## keeps only the cases in which we carried on playing (the opened gate wasn't the one with the prize)

switchedchoice <- rep(1,length(steptwo$randomgates)) 

for (i in 1:length(steptwo$randomgates)) {
  choice <- options[options != steptwo$randomchoice[i]]
  switchedchoice[i] <- choice[ceiling(runif(1,min = 0, max = 2))]
}

## Now we check how many times you guess the prize gate when you switch your initial choice

prizesswitching <- steptwo$randomgates == switchedchoice

sum(prizesswitching)/length(steptwo$randomgates)

Когда я проверяю вероятность, не меняя первоначальный выбор в тех случаях, когда игра продолжалась (открытие ворот не совпадало с выигрышем), я получал то, что ожидал (закрыть 1/3 вероятности выигрыша приза), что относится к следующей инструкции:

carryon <- randomopen != randomgates

sum(carryon)/length(randomgates) 

Моя проблема возникает, когда я проверяю вероятность выигрыша приза после изменения своего выбора (условно, очевидно, чтобы не открыв дверь, в которой находится приз), вместо получения 1/2, как утверждает Монти Холл, я получаю 1/4, это относится к следующей инструкции:

prizesswitching <- steptwo$randomgates == switchedchoice

sum(prizesswitching)/length(steptwo$randomgates)

Я знаю, что делаю что-то плохое, потому что это уже более чем доказано, что Монти Холл держит, но я не могу обнаружить недостаток. Кто-нибудь знает, что это такое?

Если вы не знаете, в чем проблема Монти Холла, вы можете найти легкую для чтения информацию в википедии:

Игра Монти Холла

Редактировать: Как указывал @Dason, одной из проблем было то, что я вводил какую-то случайность при изменении первоначального выбора, что не имеет смысла, так как остался только один вариант.

Другая проблема заключалась в том, что я не подходил к проблеме, предполагая, что Монти Холл знает, где приз. Я изменил свой код с исходного на этот, и проблема решена:

# Prepare each variable for 10000 experiments

## First I set the seed

set.seed(4)

## Then I modelize the presence of the prize as a random variable between gates 1,2,3


randomgates <- ceiling(runif(10000, min = 0, max = 3))

## so do I with the random choice.

randomchoice <- ceiling(runif(10000, min = 0, max = 3))

## As the opening of a gate is dependent from the gate you chose (the gate you chose cannot be opened
##, neither the one with the prize does), I modelize the opening of the gate as a variable which cannot be equal to the choice.

options <- c(1:3)

randomopen <- rep(1,10000)

for (i in 1:length(randomgates)) {
  randomopen[i] <- options[options != randomchoice[i] & options != randomgates[i]]
}

##Just to make data more easy to handle, I make a dataset

dataset <- cbind(randomgates, randomchoice, randomopen)

## I format the dataset as a data frame

steptwo <- as.data.frame(dataset)

## Now we check what happens if we hold our initial choice when game carries on

steptwo$prizesholding <- steptwo$randomgates == steptwo$randomchoice

with(steptwo, sum(prizesholding))

## creating a vector of changing option, dependant on the opened gate, in the dataset that
## keeps only the cases in which we carried on playing (the opened gate wasn't the one with the prize)

steptwo$switchedchoice <- rep(1,length(steptwo$randomgates)) 

for (i in 1:length(steptwo$randomgates)) {
  steptwo$switchedchoice[i] <- options[options != steptwo$randomchoice[i] & options != steptwo$randomopen[i]]
}

## Now we check how many times you guess the prize gate when you switch your initial choice

steptwo$prizesswitching <- steptwo$randomgates == steptwo$switchedchoice

with(steptwo, sum(prizesswitching)/length(randomgates))

Answer 1

В каждом раунде есть приз за дверью и выбранная дверь. Монти Холл откроет дверь, которая не является pri__door или selected_door (setdiff между 1: 3 и вектором (выигрыш_door, selected_door), со случайным выбором между двумя, если setdiff состоит из двух элементов). Тогда дверь выключателя - это дверь, которая не выбрана или не открыта.

n <- 1e4
set.seed(2020)
df <- 
  data.frame(
    prize_door = sample(1:3, n, replace = TRUE),
    chosen_door = sample(1:3, n, replace = TRUE))

df$opened_door <- 
  mapply(function(x, y){
    available <- setdiff(1:3, c(x, y))
    available[sample(length(available), 1)]
  }, df$prize_door, df$chosen_door)

df$switch_door <- 
  mapply(function(x, y) setdiff(1:3, c(x, y)),
  df$chosen_door, df$opened_door)



with(df, mean(prize_door == chosen_door))
# [1] 0.3358
with(df, mean(prize_door == switch_door))
# [1] 0.6642

График вероятностей при увеличении n

probs <- 
  data.frame(
    chosen_p = with(df, cumsum(prize_door == chosen_door))/(1:n),
    switch_p = with(df, cumsum(prize_door == switch_door))/(1:n))

plot(probs$switch_p, type = 'l', ylim = c(0, 1))
lines(probs$chosen_p, col = 'red')
abline(h = 1/3)
abline(h = 2/3)

Answer 2

Похоже, это помогает:

n_iter <- 10000

set.seed(4)

doors <- 1:3
prizes <- sample.int(n = 3, size = n_iter, replace = TRUE)
your_pick <- sample.int(n = 3, size = n_iter, replace = TRUE)
open_door <- rep(0, n_iter)
switched_door <- rep(0, n_iter)

for (i in 1:n_iter) {
  remaining_choices <- setdiff(doors, c(your_pick[i], prizes[i]))

  if (length(remaining_choices) > 1) {
    open_door[i] <- sample(remaining_choices, size = 1)
  } else {
    open_door[i] <- remaining_choices
  }

  switched_door[i] <- setdiff(doors, c(your_pick[i], open_door[i]))
}

> mean(your_pick == prizes)
[1] 0.3305
> mean(switched_door == prizes)
[1] 0.6695

Базовые функции sample.int и sample помогают немного упростить ситуацию. Элемент remaining_choices содержит возможные двери, которые может открыть хозяин игрового шоу, длина которого составляет 1 или 2 в зависимости от вашего первоначального выбора. Если длина равна 2, мы выбираем из этого вектора, и если это 1, эта дверь автоматически открывается.