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:

set.seed(1)

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:

library(caTools)

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] +

                                HCP.sums[i]

}

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.