(This article was originally published at R snippets, and syndicated at StatsBlogs.)

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.

**Please comment on the article here:** **R snippets**