Concise Proof of the Monty Hall Problem in R
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 &lt;- 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 &lt;- 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 &lt;- sample(10^6, 10^6)
 boolTest &lt;- rep(FALSE, 10^6)
 boolTest[mySamp] &lt;- TRUE
 testIndex &lt;- 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 &lt;- vapply(1:numReps, function(x) {
                mySet &lt;- intDoors[intDoors != car[x] &amp; intDoors != contestant[x]]
                if (length(mySet) &gt; 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 &lt;- function(numReps) {
    car &lt;- sample(intDoors, numReps, replace = TRUE)
    contestant &lt;- sample(intDoors, numReps, replace = TRUE)

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

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

    wins &lt;- sum(contestant == car)
    (wins/numReps)
}


The OP's function is:

funOP &lt;- function(numReps) {
    wins &lt;- 0

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

        # The contestant chooses a door
        contestant &lt;- sample(doors)[1]

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

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

        # Count the wins
        if(contestant == car){wins &lt;- 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