Saturday, December 8, 2012

Bridge hand distribution: simulation vs exact calculation

Recently I played bridge with my friends. Being frustrated with several consecutive poor hand distributions we asked ourselves a question what is the probability of having a hand good enough for a small slam. A well known rule of thumb is that you need 33+ HCP for 6NT. But we could not find information about the probability of such an event. So we decided to calculate it.
And the question was whether to calculate it exactly or simulate the result. With GNU R you can do both. However - not surprisingly - simulation is much easier and faster to perform and exact calculation is coding error prone. Have a look at the codes.

First start with simulated result. The code is clean and simple:

deck <- rep(c(1:4, rep(0, 9)), 4)
points <- replicate(1000000, sum(sample(deck, 26)))
print(2 * mean(points > 32))

and we get the following result:

[1] 0.00687

On the other hand exact result requires careful counting of all possible card combinations:

result <- data.frame(HCP=0:40, exact = 0)
result[1, 2] <- choose(16, 0) * choose(36, 26)
for (i in 1:16) {
    HCP.sums <- table(rowSums(combs(rep(1:4,4), i))) *
                choose(36, 26 - i)
    levels <- as.integer(names(HCP.sums)) + 1
    for (i in 1:length(HCP.sums))
        result[levels[i], 2] <- result[levels[i], 2] +
result[, 2] <- result[, 2] / sum(result[, 2])
print(2 * sum(result$exact[result$HCP > 32]))

It took: (a) some thinking before coding, (b) 10x more time to code and (c) one mistake in the process to get the desired results. The output of the procedure is the following:

[1] 0.0069593

So the conclusion is that some pair gets a hand good enough for 6NT on the average once in 150 games.

I wanted check sure whether the simulation and exact results are similar so I decided to plot the resulting HCP distributions using the code:

par(mar=c(4, 4, 1, 1))
result$sim <- sapply(0:40, function(x) { mean(points == x) })
plot(result$HCP, result$exact, pch = 4,
     xlab = "HCP", ylab = "Probability")
points(result$HCP, result$sim, pch = 3, col = "red")

Tthe result given below shows that simulation approximates exact result pretty well.


  1. Not bad -- now write a version for suited slams, where you must (at the very least) throw in point counts for short (<2 cards) and long( >6 cards) suits. :-)

  2. Nice. Another extension could be to check for control of all four suits. For example, either all four Aces or three Aces plus a protected King for a small slam in No Trump. In suited slams, you'd need to include singletons and voids. Much higher level of programming difficulty, I guess.

  3. Another thought: you've calculated the chances of getting *at least* 6NT. Properly speaking, you should calculate the chance of getting *exactly* 6NT by subtracting the chance of getting 7NT.

  4. print(2 * mean(points > 32))

    Why multiply by 2? Is it because you want the probability of either partnerships having as slam hand as opposed to the probability of just you and your partner have a slam hand?

    1. Yes - I count the probability that either partnership has at least 32 points.

      You could alternatively write:
      print(mean(pmax(points, 40 - points) > 32))

      This is equivalent because the distribution is symmetrical and the events are mutually exclusive.