Bridge hand distribution: simulation vs exact calculation

December 8, 2012

(This article was first published on R snippets, and kindly contributed to R-bloggers)

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 on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.


Mango solutions

RStudio homepage

Zero Inflated Models and Generalized Linear Mixed Models with R

Quantide: statistical consulting and training


CRC R books series

Contact us if you wish to help support R-bloggers, and place your banner here.

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)