This story is about a sad man. The man is ordered by his doctor to take one pill each day in order to stay happy. The doctor gives him a pill box with $N$ pills. The sad man decides in his sadness that he wants only to be slightly happy. Therefore, he takes only half the dosage each day.

On day 1 he randomly grabs a whole pill from the box, he breaks it in half and consumes one half while putting the other half back in the box. On day 2 he might either grab one of the $N-1$ whole pills with probability $P(X_2=W) = \frac{N-1}{N}$, or the 1 half pill with probability $P(X_2=H) = \frac{1}{N}$. If he gets the half pill, he simply consumes it. Before the bottle is empty, $2\cdot N$ days will have passed.

The sad man is mathematically inclined and starts to wonder about how the contents of his pill box change over time. He wonders:

  1. If a lot of sad men around the world do this, what is the average number of days that will pass before only halves remain inside the bottle?
  2. The first day the pill will be a whole one. The last day the pill will be half. But how does the probability $P(X_t=H)$ behave for $1<t\leq 2N$?

Quick visualization through code

Before we analyze this problem mathematically, it is way too tempting to write some code and make some figures.

function experiment(w_arr::Array{Int, 2}, h_arr::Array{Int, 2}, day::Int, man::Int)
    prob_w = w_arr[day, man] / (h_arr[day, man] + w_arr[day, man])
    if rand() < prob_w   # Taken a whole: split it, and return one half
        w_arr[day + 1, man] = w_arr[day, man] - 1
        h_arr[day + 1, man] = h_arr[day, man] + 1
    else  # Taken a half: consume it
        w_arr[day + 1, man] = w_arr[day, man]
        h_arr[day + 1, man] = h_arr[day, man] - 1
    end
end

function main(N_pills::Int, N_men::Int)
    N_days = 2 * N_pills
    w_arr = zeros(Int, N_days, N_men)
    h_arr = similar(w_arr)

    w_arr[1, :] .= N_pills

    for i=1:N_men
        for j=1:N_days-1
            experiment(w_arr, h_arr, j, i)
        end
    end
    return w_arr, h_arr
end

We can track the development of the contents of the pill box with e.g. 100 pills and a million sad men by using

wholes, halves = main(100, Int(1e6));

Plotting these results gives the following interesting figures:

svg
The average (solid red) number of halves, standard deviation (ribbon), and 10 example trajectories.
svg
The same figure for wholes