Bridge hand distribution: simulation vs exact calculation

[This article was first published on R snippets, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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.

To leave a comment for the author, please follow the link and comment on their blog: R snippets. offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)