I recently stumbled onto “The Riddler” series from Five Thirty Eight. Each week, they present a puzzle to solve that involves math and probability, and since I don’t already have enough to do :grimmace:, I thought I’d give this week’s a shot:

Where Will The Seven Dwarfs Sleep Tonight?Each of the seven dwarfs sleeps in his own bed in a shared dormitory. Every night, they retire to bed one at a time, always in the same sequential order, with the youngest dwarf retiring first and the oldest retiring last. On a particular evening, the youngest dwarf is in a jolly mood. He decides not to go to his own bed but rather to choose one at random from among the other six beds. As each of the other dwarfs retires, he chooses his own bed if it is not occupied, and otherwise chooses another unoccupied bed at random.

- What is the probability that the oldest dwarf sleeps in his own bed?
- What is the expected number of dwarfs who do not sleep in their own beds?

I had some ideas of how to begin, but as it’s been a while since I worked on this kind of logic puzzle, I thought the safest thing to do at the outset would be to determine the answer(s) emprically. I can definitely write an algorithm that simulates the conditions above, and then it’s just a matter of running it a few thousand times and checking how many times we get the outcomes mentioned.

## Bootstrapping probabilities

As is my want, I did the work in julia. First, I made a `Beds`

type that
holds the information for which beds are occupied and by whom. It’s really just
a wrapper around an array, where the index for the array = the dwarf who’s bed
it is (`1`

for the youngest, `7`

for the oldest), and the value at each index is
either `0`

(the bed is empty) or a number `1:7`

indicating which dwarf is
sleeping there.

```
struct Beds
beds::Vector{Int}
# a constructor that builds a `Beds` object with `n` zeros
function Beds(n::Int)
new(zeros(1:n))
end
end
# these functions return a list of indicies for beds that are vacant or occupied
vacant(b::Beds) = findin(b.beds, 0)
occupied(b::Beds) = findin(x -> x!=0, b.beds)
# this function returns `true` if theres no one in the bed
isvacant(b::Beds, i::Int) = b.beds[i] == 0
```

I also made a convenience function that adds a particular dwarf to a particular bed, but spits an error if the bed already has someone in it, or if that dwarf is already sleeping somewhere else.

```
function sleep!(b::Beds, bed::Int, dwarf::Int)
isoccupied(b, bed) && error("That bed's already taken!")
in(dwarf, b.beds) && error("That dwarf already has a bed!")
b.beds[bed] = dwarf
end
```

Now we have our constructs, we need to implement the logic of the puzzle. So I made a function that shows what happens when a particular dwarf gets sleepy:

- If it’s
`dwarf1`

(the youngest), he picks a random bed other than his own - For the other dwarves, if their bed is vacant, they just grab their own bed
- If their bed is
*not*vacant, they grab a random bed from among the vacant ones

```
function sleepy(dwarf::Int, b::Beds)
if dwarf == 1
bed = rand(vacant(b)[2:end])
elseif isvacant(b, dwarf)
bed = dwarf
else
bed = rand(vacant(b))
end
sleep!(b, bed, dwarf)
end
```

At bedtime, we iterate through the dwarves from `1:7`

, and they go to bed
according to the logic above. The I also made a function that will iterate the
go-to-bed process some number of times, and returns the answers to the two
questions in the riddle.

```
function bedtime(numbeds::Int)
b = Beds(numbeds)
for d in 1:numbeds
sleepy(d, b)
end
# 1. Does the eldest dwarf sleep in his own bed?
b.beds[end] == numbeds ? eldest=true : eldest=false
# 2. What proportion of dwarves are not in their own beds?
wrongbed = sum([i != b.beds[i] for i in 1:numbeds]) / numbeds
return eldest, wrongbed
end
function iterbed(numbeds::Int, itr::Int)
# store the outcomes of each run
eldest = []
wrongbed = []
for _ in 1:itr
(x, y) = bedtime(numbeds)
push!(eldest, x)
push!(wrongbed, y)
end
return mean(eldest), mean(wrongbed)
end
```

Put these functions into a file called `dwarves.jl`

(or download the one I made
here), then head into a julia terminal:

```
julia> include("static/files/dwarves.jl");
julia> iterbed(7, 1000)
(0.41, 0.41199999999999853)
```

In other words:

- What is the probability that the oldest dwarf sleeps in his own bed? $0.41$
- What is the expected number of dwarfs who do not sleep in their own beds? $0.41 * 7 ≈ 3$

The nice thing about this code is that it’s generic enough that we can do the same caluclations for any number of dwarves.

```
julia> iterbed(2, 1000)
(0.0, 1.0)
julia> iterbed(3, 1000)
(0.244, 0.7503333333333289)
julia> iterbed(4, 1000)
(0.337, 0.60625)
julia> iterbed(5, 1000)
(0.361, 0.5231999999999994)
julia> iterbed(100, 1000)
(0.495, 0.05149999999999994)
```

Or, we can plot it:

The answers clearly converge on ~$0.5$ for `Q1`

and on $0$ for `Q2`

for
large numbers of dwarves. But how do we get the exact numbers?

## Exact Values

This one took me a while to puzzle out. For question 1, I suspected it would be
easier to find the solution to `eldest not in own bed`

, and then subtract that
from 1 to find the solution to `eldest in own bed`

. And it was easier to start
with smaller numbers of dwarves.

The smallest number we can have is 2. Under that circumstance, the youngest
always takes the oldest’s bed, so the answers to `Q1`

and `Q2`

are $0.0$ and $1$
respectively. With 3 dwarves, there’s a $1\/2$ chance that `dwarf1`

takes the
eldest’s bed, but we also have to account for what happens if he takes bed 2
($1\/2$ chance) **and** dwarf 2 takes the eldest’s bed ($^{1}⁄_{2}$ chance, since there
are 2 beds remaining). So the overall probability is:

$$^{1}⁄_{2} + (^{1}⁄_{2})*(^{1}⁄_{2}) = 0.75$$

For reasons that should become obvious, it was easier to think of this as:

$$\frac{1 + ^{1}⁄_{2}}{2}$$

What about for 4 dwarves?

You might notice there seems to be some repeating pattern… if not, it should be really clear from the n = 5 example:

The pattern seems to be the following recursive function:

$$f(n) = 1 + \sum_{i=2}^{n-1} \frac{f(i)}{i}$$

In otherwords,

```
function f(n::Int)
n < 2 && error("cannot calculate for less than 2 dwarves")
if n == 2
return 1
else
return 1 + sum([f(x) / x for x in 2:n-1])
end
end
```

Then the probability $P$ of the eldest being in his own bed is

$$1 - \frac{f(n)}{n-1}$$

```
p_eldest(n::Int) = 1 - f(n) / (n-1)
```

There’s probably a more efficient way to do this (the recursion basically makes the code unrunnable on my system for more than ~25 dwarves), but it goes pretty nicely with our simulations:

For `Q2`

, it’s a similar idea, but we can look at it a different way. Say there
are 3 dwarves - `dwarf1`

can either choose `bed2`

or `bed3`

. There’s a $^{1}⁄_{2}$
chance he’ll pick `bed3`

, in which case ^{2}⁄_{3} dwarves will be in the wrong bed.
Or, there’s a $^{1}⁄_{2}$ chance he’ll pick `bed2`

, in which case there’s a 50-50
split between all 3 dwarves being in the wrong bed and $^{2}⁄_{3}$ dwarves being in the
wrong bed.

As the number of dwarves increases, the decision tree gets bigger and bigger,
but it follows a similar pattern. There’s a $\frac{1}{n-1}$ chance that `dwarf1`

picks the last bed, and two dwarves will be in the wrong. There’s a
$\frac{1}{n-1}$ chance that he’ll pick the second to last bed, which leaves a
50-50 chance of 2 dwarves in the wrong bed and 3 dwarves in the wrong bed. Then,
there’s a $\frac{1}{n-1}$ chance of picking the 3rd to last bed, which leaves a
$^{1}⁄_{3}$ chance of 2 in the wrong, a 1/3rd chance of 3 in the wrong (if the 3rd
from the last picks the 5th bed), and 1/3rd chance of a 50-50 split between 3 in
the wrong and 4 in the wrong. Again, it keeps iterating in a predictable
pattern.

I’m a bit ashamed to say that I spent the better part of the weekend figuring out what that pattern was and how to code it in an alorithm. The issue is that each term in the equation has to be modified by how many dwarves are already in the wrong bed (imagine 1 picks 2, 2 picks 3, 3 picks 4 and so on).

I’m not even sure I can explain the logic, but basically, if you’re the eldest dwarf, and your bed is taken, you’ve only got 1 choice, and the number of wrong beds will be 2 + however many other dwarves are in the wrong beds. If you’re the second oldest, you have 2 choices, and the number wrong will be 2 or 3 + however many and so on. Each level down, you get one more possible choice.

I thought about this in the number of terms (if there are 3 beds open, there are three terms) and the number of levels of fractions below you). For 4 dwarves, the equation is:

$$P_(wrongbed,4) = \frac{2 + \frac{2+3}{2} + \frac{2+3+\frac{3+4}{2}}{3}}{3}$$

The general formula I figured out (again, there may be one way more simple, this is what I got to) is:

$$g_(t,l) = \frac{l + \sum_{i=1}^{t-1} g_(i,t+1)}{t}$$

Where $t$ is the term number and $l$ is the level. For $n$ dwarves then, the number that should be in the wrong bed is

$$N_(wrongbed, n) = \frac{\sum_{i=1}^{n-1} g(i,2)}{n-1}$$

and the proportion is just that divided by $n$.

```
function g(t, l)
if t == 1
return Float64(l)
else
return (l + sum([g(i,l+1) for i in 1:t-1])) / t
end
end
p_wrongbeds(n) = sum([g(i,2) for i in 1:n-1]) / (n - 1) / n
```

Finally, let’s plot it to see if it lines up with the simulation:

Nice!