(This article was first published on

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.**R snippets**, and kindly contributed to R-bloggers)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

print**<-**replicate**(**1000000, sum**(**sample**(**deck, 26**)))****(**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

print**[**, 2**]****<-**result**[**, 2**]****/**sum**(**result**[**, 2**])****(**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

points**=**"HCP", ylab**=**"Probability"**)****(**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 his blog:**R snippets**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...