Topics for Probability Presentation, Preferably Using Computational Solutions

Question

I'm looking for interesting applications of computing science for finding the probability of otherwise difficult problems. This is for an undergraduate level presentation in a probability course of around 30 minutes, from a comp sci major who has beginner to intermediate skills in both probability/statistics and comp sci.

I hope to have a topic that is interesting and can be explained without entirely relying on proofs and equations.

Answer 1

A couple of ideas:

Frustration Solitaire: Turn over cards in a well-shuffled deck. Call out 1 and turn over the first card, call out 2 and turn over the second card, and so on to the 13th card. Then start again calling out 1, 2, ..., 13 in sequence; repeat two more sequences of 13 to finish the deck. If, at any point, the denomination of your card matches the denomination you call out (a "hit"), you lose the game. Otherwise, you win. (The name of this game comes from the perhaps surprisingly low probability of winning.)

The exact probability of winning 0.01623 is found by an inclusion-exclusion procedure described at https://arxive.org/pdf/math/0703900.pdf. [Also, It has been suggested to approximate the distribution of the number $X$ of hits as $\mathsf{Pois}(\lambda = 52/13 = 4),$ so $P(Win) = $ $P(X = 0) \approx$ $e^{–4} = $ $0.018$ (which is a reasonable exercise using the Poisson distribution), but $\mathsf{Binom}(n = 52, p = 1/13)$ gives a better approximation $P(X = 0) \approx (12/13)^{52} = 0.0156,$ as suggested by the simulation of a million games below.]

Simulation using R statistical software:

set.seed(1234)
m = 10^6;  deckd = rep(1:13, times=4)
x = replicate( m, sum(deckd == sample(deckd, 52)) ) 
mean(x==0);  2*sd(x==0)/sqrt(m);  mean(x);  2*sd(x)/sqrt(m)
[1] 0.015997       # sim P(Win)
[1] 0.0002509272   # aprx margin of sim error for P(Win)
[1] 4.004366       # sim E(Hits)
[1] 0.003883274    # aprx margin of sim error for E(Hits)

In the figure below, the blue dots are binomial approximations and the red circles are Poisson approximations. For $P(X = 0),$ as well as for some other values, the binomial approximation is the better one.

mx = max(x);  cutp = (-1:mx)+.5
hist(x, prob=T, br=cutp, col="skyblue2", main="Simulated Frustration Hits")
k = 0:mx; pdfb = dbinom(k, 52, 1/13);  points(k, pdfb, pch=19, col="blue")
pdfp = dpois(k, 4);  points(k, pdfp, col="red")

Maybe you are not interested in the part about Poisson and binomial approximations, in which case it is easy enough to skip them.

If you want to demonstrate the game with a deck of cards. that may be feasible. Additional simulation shows that, if you don't win, then the expected number of cards turned over is less than 13.

Birthday Matches: This is a well-known problem in elementary probability with many links on this site. Here is a very brief summary. Suppose you have $n = 25$ randomly selected people in a room. Then the probability is about 0.57 that at least two of the 25 were born on the same day of the year.

Assuming that there are 365 equally likely birthdays in a year, the combinatorial solution is $$P(\text{No Match}) = \frac{P(365, 25)} {365^{25}} = \prod_{i=0}^{24}\left(1 - \frac{i}{365}\right) = 0.4313.$$

Most people are surprised that the probability of a match is so high. Maybe this is because they are thinking of the chances their own birthday will be matched, forgetting that there are ${25 \choose 2} = 300$ possible pairs of people. Or maybe it is because it takes 366 people to ensure a match, and 25 is much smaller than that. However, the probability of a match increases in a non-linear manner with the number of people in the room.

The following graph shows this non-linear relationship:

max.n = 60;  p = numeric(max.n)   # at the start, p is 60-vector of 0's
for (i in 1:max.n) {
   q = 1 - (0:(i-1))/365
   p[i] = 1 - prod(q)  }       # replace ith element of p with probability
n=1:max.n
plot(n, p, ylab="p = P(Match)", xlab="n = Nr of People in Room")
 lines(c(0, 23, 23), c(.5, .5, 0), col="blue", lty="dashed", lwd=2)
 abline(h=0:1, col="green2"); abline(v=0, col="green2")

The dotted lines point out that the probability of a match exceeds 1/2 for 23 people in the room.

The following simulation in R, for 100,000 rooms, verifies that the probability of no match is about 0.43 and also shows that the expected number of matches is about 0.8.

m = 10^5;  n = 25;  x = numeric(m)  
for (i in 1:m)  {
  b = sample(1:365, n, repl=T)
  x[i] = n - length(unique(b))  }

mean(x);  2*sd(x)/sqrt(m)
[1] 0.80434
[1] 0.005447805

mean(x==0);  2*sd(x==0)/sqrt(m)
[1] 0.43291
[1] 0.003133697

The histogram below shows the simulated distribution of the number of matches in a room with 25 people.

cp = (0:(max(x) + 1)) - 0.5
hist(x, breaks=cp, prob=T, col="skyblue2")

Note: In fact the distribution of births in the US is not exactly uniform. Actual 1997-99 data show that monthly averages range from a low of about 94.9% of uniform in January 1999 to a high of about 107.4% in September 1999. It is possible to run a simulation with the actual birthday distribution. The results for the probability of a match and the expected number of matches do not change by more than 1 in the second decimal place. So for practical purposes the assumption of equally likely birthdays leads to good results.

But you wouldn't want to try a demonstration of match probabilities at a convention of twins or at a meeting of the Sagittarian Society. There would be way too many matches.