# Counting and chance in Python: simulating the matching and birthday problems

The R section of Blitzstein and Hwang’s first chapter, rebuilt in Python in the style of CS61A: small named functions, doctests, and simulations passed around as values.

## 0. What this note is for

Chapter 1 of Blitzstein and Hwang’s *Introduction to Probability* ends with a
short section of R code: vectors, factorials, random sampling, and two
simulations that check results proved earlier in the chapter. This note covers
the same ground in Python.

The translation is not line by line. R is built around vectors; Python is
built around names, functions, and lists. So the code below follows the habits
of CS61A instead:

- every idea gets a **small named function** with a docstring;
- every docstring carries **doctests**, examples that Python can check;
- a simulation is a **function passed as a value**, not a special command.

Everything uses only the standard library (`math`, `random`, `collections`),
so it runs in any Python 3 and in the exercises on this page.

## 1. From vectors to lists

R’s `c(3,1,4,1,5,9)` becomes a Python list:

```python
>>> v = [3, 1, 4, 1, 5, 9]
>>> sum(v), max(v), min(v), len(v)
(23, 9, 1, 6)
```

The sequence $1,2,\dots,n$, written `1:n` in R, is `range(1, n + 1)` in
Python. A range stops *before* its second argument, which is why the `+ 1`
is there. Counting down from $m$ to $n$ needs a negative step:
`range(m, n - 1, -1)`.

### 1.1. Positions start at zero

The biggest change is indexing. R numbers positions from 1; Python numbers
them from 0. So `v[0]` is the first entry, and the $i$th entry in the book’s
sense is `v[i - 1]`.

### What would Python display? Indexing

*Fill in the blank · en*

What would Python display?

```python
>>> v = [3, 1, 4, 1, 5, 9]
>>> v[1]
```

Your answer: 1 (Accepted absolute error: 0)

#### Hints

1. Python numbers positions from 0, not from 1.

2. `v[0]` is `3`. Which entry comes next?

#### Reference solution

`1`. Python positions start at 0, so `v[1]` is the **second** entry. In R, `v[1]` would be the first entry, `3`.

Source: [Adapted from Blitzstein & Hwang, Introduction to Probability, §1.8 (R), rewritten in Python](http://stat110.net) · Original exercise; the book's R examples are paraphrased, not copied

Picking out the 1st, 3rd, and 5th entries, `v[c(1,3,5)]` in R, is a list
comprehension over positions. Removing the 2nd through 4th entries,
`v[-(2:4)]` in R, is two slices joined together:

```python
>>> [v[i - 1] for i in (1, 3, 5)]
[3, 4, 5]
>>> v[:1] + v[4:]
[3, 5, 9]
```

A negative index means something else in Python: `v[-1]` is the *last*
entry, not “everything except the first”.

### 1.2. Componentwise work is written out

R applies arithmetic to each entry of a vector automatically. In Python a
list comprehension says so explicitly. The vector
$\left(1,\frac{1}{2^2},\frac{1}{3^2},\dots,\frac{1}{100^2}\right)$ is

```python
>>> squares = [1 / k ** 2 for k in range(1, 101)]
>>> squares[:3]
[1.0, 0.25, 0.1111111111111111]
```

Powers use `**`. In Python `^` is bitwise exclusive or, so `2 ^ 3` is `1`,
not `8`. This is one of the easiest mistakes to make when porting R code.

R also *recycles* a shorter vector: `v + 3` adds 3 to every entry. Python
lists do not.

### What would Python display? No recycling

*Single choice · en*

What happens here?

```python
>>> v = [3, 1, 4]
>>> v + 3
```

- \[ ] added: `[6, 4, 7]`
- \[ ] appended: `[3, 1, 4, 3]`
- \[ ] error: A `TypeError`

Your answer: error

#### Hints

1. `+` between two lists joins them. Is `3` a list?

2. Python does not repeat a shorter value to fit a longer one.

#### Reference solution

Python raises `TypeError: can only concatenate list (not "int") to list`. For lists, `+` means concatenation, and it needs a list on both sides. Write `[x + 3 for x in v]` to add 3 to each entry, which gives `[6, 4, 7]`.

Source: [Adapted from Blitzstein & Hwang, Introduction to Probability, §1.8 (R), rewritten in Python](http://stat110.net) · Original exercise; the book's R examples are paraphrased, not copied

## 2. Factorials and binomial coefficients

The standard library has both: `math.factorial(n)` and `math.comb(n, k)`.

```python
>>> import math
>>> math.factorial(20)
2432902008176640000
>>> math.comb(10, 3)
120
```

Python integers have no fixed size, so `math.factorial(1000)` is an exact
2568-digit integer. The limit appears only when we convert to a
floating-point number: `float(math.factorial(171))` raises `OverflowError`,
the same point at which R returns `Inf`. For large $n$, work with logarithms.
Since $n!=\Gamma(n+1)$, the natural log $\log(n!)$ is `math.lgamma(n + 1)`,
and $\log\binom{n}{k}$ is `lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1)`.

Writing $\binom{n}{k}$ ourselves is a good first recursion. The base cases
come straight from the definition, and the recursive case is Pascal’s rule.

> **Theorem · 2.1 · Pascal’s rule**
>
> For integers $n\ge 1$ and $1\le k\le n-1$,
>
> $$
> \binom{n}{k}=\binom{n-1}{k-1}+\binom{n-1}{k}.
> $$
>
> Choose a particular person out of the $n$. A committee of size $k$ either
> includes that person, leaving $k-1$ more to choose from the other $n-1$, or
> it does not, leaving all $k$ to choose from the other $n-1$.

### Binomial coefficients by recursion

*Python · en*

Complete `choose(n, k)` **recursively**, using Pascal's rule $\binom{n}{k}=\binom{n-1}{k-1}+\binom{n-1}{k}$. Do not use `math.comb` or factorials. Return 0 when $k<0$ or $k>n$.

```python
def choose(n, k):
    """Return the number of ways to choose K items from N, using Pascal's rule.

    >>> choose(5, 2)
    10
    """
    pass
```

#### Public tests

Small case

```pycon
>>> choose(5, 2)
10
```

Middle of a row

```pycon
>>> choose(10, 5)
252
```

Choose none or all

```pycon
>>> choose(4, 0)
1
>>> choose(4, 4)
1
```

Too many

```pycon
>>> choose(3, 5)
0
```

#### Hints

1. There are two base cases: impossible choices, and choosing none or all.

2. If $k$ is 0 or equal to $n$, there is exactly one way.

3. Otherwise, return the sum of two smaller `choose` calls.

#### Reference solution

```python
def choose(n, k):
    if k < 0 or k > n:
        return 0
    if k == 0 or k == n:
        return 1
    return choose(n - 1, k - 1) + choose(n - 1, k)
```

The first base case handles impossible requests; the second is the edge of Pascal's triangle. Each recursive call lowers $n$, so the recursion ends. It makes many repeated calls, which is fine for small $n$ and a good reason to meet memoization later.

Source: [Adapted from Blitzstein & Hwang, Introduction to Probability, §1.8 (R), rewritten in Python](http://stat110.net) · Original exercise; the book's R examples are paraphrased, not copied

## 3. Random samples

The `random` module plays the role of R’s `sample`. Positions in the book
are the numbers $1,\dots,n$, so we sample from `range(1, n + 1)`.

| In R | In Python | Meaning |
| --- | --- | --- |
| `sample(n, k)` | `random.sample(range(1, n + 1), k)` | $k$ draws without replacement |
| `sample(n, k, replace=TRUE)` | `random.choices(range(1, n + 1), k=k)` | $k$ draws with replacement |
| `sample(n)` | `random.sample(range(1, n + 1), n)` | a random permutation |
| `sample(letters, 7)` | `random.sample(string.ascii_lowercase, 7)` | seven distinct letters |
| `sample(4, 3, replace=TRUE, prob=…)` | `random.choices(range(1, 5), weights=[0.1, 0.2, 0.3, 0.4], k=3)` | weighted, with replacement |

`random.choices` always samples with replacement; the standard library has no
weighted sampling without replacement.

These numbers are pseudo-random: a deterministic algorithm produces them.
That is useful. Calling `random.seed(110)` first makes a simulation
reproducible, so someone else running your code gets the same numbers.

## 4. A simulation is a function

R’s `replicate(times, expression)` evaluates an expression many times. In
Python we pass the experiment itself, a function of no arguments, and call it
as often as we like:

```python
def replicate(times, experiment):
    """Call EXPERIMENT, a function of no arguments, TIMES times.
    Return the results as a list.

    >>> replicate(3, lambda: 'hi')
    ['hi', 'hi', 'hi']
    """
    return [experiment() for _ in range(times)]
```

This is the CS61A idea of a **higher-order function**. `replicate` knows
nothing about cards or birthdays; it only knows how to repeat. The `lambda`
wraps an expression in a function so that it is evaluated later, and anew on
each call, rather than once.

## 5. The matching problem

Shuffle a deck of cards labeled $1,\dots,n$ and turn them over one at a time.
A **match** happens when card $i$ is the $i$th card turned over. The book
shows (Example 1.6.4) that the probability of at least one match is

$$
1-\frac{1}{2!}+\frac{1}{3!}-\cdots+(-1)^{n+1}\frac{1}{n!}\;\approx\;1-\frac{1}{e}\approx 0.632.
$$

To simulate it we need to count matches in a given deck. In Python, `True`
and `False` behave as `1` and `0` in arithmetic, just as `TRUE` and `FALSE`
do in R.

### What would Python display? Adding truth values

*Fill in the blank · en*

What would Python display?

```python
>>> deck = [2, 1, 3, 4]
>>> sum([card == position for position, card in enumerate(deck, start=1)])
```

Your answer: 2 (Accepted absolute error: 0)

#### Hints

1. The list inside is `[False, False, True, True]`.

2. In arithmetic, `True` counts as `1` and `False` as `0`.

#### Reference solution

`2`. Cards 3 and 4 are in their own positions, so the list is `[False, False, True, True]`, and its sum counts the `True` values.

Source: [Adapted from Blitzstein & Hwang, Introduction to Probability, §1.8 (R), rewritten in Python](http://stat110.net) · Original exercise; the book's R examples are paraphrased, not copied

### Count the matches

*Python · en*

Complete `count_fixed_points(deck)`. `deck` lists the cards in the order they are turned over; card $i$ matches when it is turned over $i$th. Return the number of matches.

```python
def count_fixed_points(deck):
    """Return how many cards in DECK sit in their own position.

    Positions are numbered from 1, like the cards.

    >>> count_fixed_points([1, 2, 3])
    3
    >>> count_fixed_points([2, 1, 3])
    1
    """
    pass
```

#### Public tests

Every card matches

```pycon
>>> count_fixed_points([1, 2, 3])
3
```

One match

```pycon
>>> count_fixed_points([2, 1, 3])
1
```

No matches

```pycon
>>> count_fixed_points([3, 1, 2])
0
>>> count_fixed_points([2, 1])
0
```

Empty deck

```pycon
>>> count_fixed_points([])
0
```

#### Hints

1. `enumerate(deck, start=1)` gives each position together with its card.

2. Compare each card with its position, then add up the results.

#### Reference solution

```python
def count_fixed_points(deck):
    return sum(card == position for position, card in enumerate(deck, start=1))
```

Each comparison is `True` or `False`, which `sum` counts as 1 or 0. An empty deck has no matches, so the sum is 0.

Source: [Adapted from Blitzstein & Hwang, Introduction to Probability, §1.8 (R), rewritten in Python](http://stat110.net) · Original exercise; the book's R examples are paraphrased, not copied

With `count_fixed_points` written, one shuffle and a full simulation are each
a line or two:

```python
import random

def shuffled_matches(n):
    """Shuffle cards 1..N and count the matches."""
    deck = random.sample(range(1, n + 1), n)
    return count_fixed_points(deck)

random.seed(110)
results = replicate(10 ** 4, lambda: shuffled_matches(100))
sum(r >= 1 for r in results) / 10 ** 4     # proportion with a match
```

> **Example · 5.1 · What one run gave**
>
> With seed 110 and $10^4$ shuffles of 100 cards, the proportion of decks with
> at least one match was $0.6267$, against $1-1/e\approx0.6321$. A different
> seed or Python version can give a slightly different proportion. The gap
> between the two shrinks roughly like $1/\sqrt{10^4}=0.01$.

## 6. The birthday problem

With $k$ people and 365 equally likely birthdays (ignoring February 29), the
probability that at least two share a birthday is

$$
1-\frac{365\cdot364\cdots(365-k+1)}{365^k}.
$$

`math.prod` computes the product directly, and exact integer arithmetic
keeps it accurate until the final division:

```python
>>> k = 23
>>> 1 - math.prod(range(365 - k + 1, 366)) / 365 ** k
0.5072972343239854
```

As a function with doctests, multiplying one factor at a time:

```python
def birthday_match(k, days=365):
    """Probability that at least two of K people share a birthday.

    >>> round(birthday_match(23), 3)
    0.507
    >>> birthday_match(1)
    0.0
    >>> birthday_match(366)
    1.0
    """
    no_match = 1.0
    for i in range(k):
        no_match *= (days - i) / days
    return 1 - no_match
```

R’s `qbirthday(p)` answers the reverse question: how many people are needed
for the probability to reach $p$? Python’s standard library has no
equivalent, so write it.

### How many people are needed?

*Python · en*

`birthday_match(k, days)` is provided. Complete `people_needed(p, days)`: return the smallest number of people $k\ge1$ for which `birthday_match(k, days)` is at least `p`. This is what R's `qbirthday` computes. Assume $0<p\le1$.

```python
def birthday_match(k, days=365):
    """Probability that at least two of K people share a birthday."""
    no_match = 1.0
    for i in range(k):
        no_match *= (days - i) / days
    return 1 - no_match


def people_needed(p, days=365):
    """Return the smallest number of people for which
    birthday_match reaches at least P.

    >>> people_needed(0.5)
    23
    """
    pass
```

#### Public tests

Even odds

```pycon
>>> people_needed(0.5)
23
```

Ninety percent

```pycon
>>> people_needed(0.9)
41
```

Ninety-nine percent

```pycon
>>> people_needed(0.99)
57
```

A shorter year

```pycon
>>> people_needed(0.5, days=2)
2
```

#### Hints

1. Start with one person and add people until the probability is large enough.

2. A `while` loop fits: keep going while `birthday_match(k, days) < p`.

#### Reference solution

```python
def people_needed(p, days=365):
    k = 1
    while birthday_match(k, days) < p:
        k += 1
    return k
```

The loop always stops: with `days + 1` people a shared birthday is certain, so the probability reaches 1.

Source: [Adapted from Blitzstein & Hwang, Introduction to Probability, §1.8 (R), rewritten in Python](http://stat110.net) · Original exercise; the book's R examples are paraphrased, not copied

### 6.1. Simulating instead

To simulate, draw 23 birthdays with replacement and ask how often the most
common birthday appears. `collections.Counter` does what R’s `tabulate` does:

```python
from collections import Counter

def most_shared(k, days=365):
    """Draw K birthdays; return the size of the largest group sharing one."""
    birthdays = random.choices(range(1, days + 1), k=k)
    return max(Counter(birthdays).values())

random.seed(110)
r = replicate(10 ** 4, lambda: most_shared(23))
sum(x >= 2 for x in r) / 10 ** 4           # at least one shared birthday
sum(x >= 3 for x in r) / 10 ** 4           # at least one triple
```

The simulation answers questions the formula does not. R’s
`pbirthday(23, coincident=3)` gives about $0.014$ for a *triple* birthday
among 23 people, and `qbirthday(0.5, coincident=3)` gives 88. Python has no
built-in for triples, but `most_shared` handles them with no extra work.

> **Example · 6.1 · What one run gave**
>
> With seed 110 and $10^4$ groups, 23 people shared a birthday in $0.5039$ of
> groups (exact: $0.5073$) and had a triple in $0.013$ of them. With 88 people,
> a triple appeared in $0.5073$ of groups, consistent with R’s answer of 88.

Unequal birthday probabilities make the formula hard but the simulation
easy: pass `weights=` to `random.choices`.

## 7. Comments and docstrings

A `#` starts a comment in both languages. Python adds the docstring: the
string at the top of a function, which says what the function does. With
doctests inside, the documentation is also a test. Run

```bash
python3 -m doctest simulation.py -v
```

and Python checks every `>>>` example in the file. Docstrings for what a
function promises, comments for why a line is written the way it is.

## References

- Joseph K. Blitzstein and Jessica Hwang, *Introduction to Probability*,
  Chapter 1, §1.6 (Example 1.6.4, the matching problem) and §1.8 (R).
  Course materials: [stat110.net](http://stat110.net).
- [CS61A](https://cs61a.org/) and John DeNero’s
  [*Composing Programs*](https://www.composingprograms.com/), for the style:
  higher-order functions, doctests, and “What would Python display?”
- Python documentation: [`random`](https://docs.python.org/3/library/random.html),
  [`math`](https://docs.python.org/3/library/math.html), and
  [`collections.Counter`](https://docs.python.org/3/library/collections.html#collections.Counter).