4
\$\begingroup\$

This is an attempt at a concise proof of the Monty Hall problem in probability and statistics in R.

For those unfamiliar, the scenario is this:

  • There exists a contestant on a game show. This contestant chooses between three doors (A, B, and C), behind one of which is a new car.
  • The host (Monty Hall) reveals what is behind one of the other two doors.
  • The contestant is offered to either stick with their same door, or switch.
  • The question is: which choice offers the highest probability of success (winning the car)?

To complete this problem, there is some necessary background:

  • The probability that the car is behind any door before anything happens is uniformly distributed.
  • The host will never reveal the car on the first opening of a door (this is what offsets the probability).

Assume (arbitrarily) that the contestant chooses door A. If the car is behind door A, the host will open door B or C with equal (50/50) probability - this is intuitive. If the car is behind door B however, the host must open door C to prevent revealing the car, thus the probability of opening door B is zero. This foreknowledge of the location of the car is what skews the problem.

The result is that the contestant ends up with a 66% chance of winning by switching their chosen door, and a 33% if they stick with their current choice. In the code below, it is compulsory for the contestant to switch - this was to improve computational efficiency:

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)
\$\endgroup\$
9
\$\begingroup\$

The fundamental ideas of your code are solid, however there are a couple of areas that could be improved. For starters, when you call sample and you only want one number returned, you can make use of the size argument. For vectors of only size 3, as in your case, there isn't much difference:

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

However, as the size of your input vector increases, there is a noticeable difference in methods.

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

Next, we can greatly improve our efficiency by taking advantage of the fact that sample can perform many replications in one go by making use of the replace argument. For example, instead of using a for loop and generating a single scenario, we can do the following:

sample(doors, n, replace = TRUE)

The gains in efficiency are tremendous:

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

Now, let's look at our data types. Since we are mainly concerned with probabilities, the actual values of our main vector doesn't matter. Since this is the case, we should look to use integers when we can. We can achieve this by using seq_along(doors) or simply length(doors). Observe:

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

While we are on integers, which is an awesome function which is very efficient and very intuitive. However, when we are indexing, it is completely unnecessary. We can use logical subsetting instead. It is cleaner and just as efficient. Observe:

 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

In order to resolve the monty vector and the second iteration of contestant, we can implement vapply along with length to avoid calls to sample where the length is 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)

And finally, instead of incrementing win, we can get it all at once using sum(contestant == car). This takes advantage of vectorization and the internal coercion from logicals to integers. Putting it all together, we get:

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)
}

The OP's function is:

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)
}

And the final comparison sees an improvement of about 4.5x.

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 results are very similar as well:

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

As a final note, when one is performing any sort of procedure that calls for verification, it is a good idea to use set.seed so you can predictably verify your results. For example:

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
\$\endgroup\$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.