# Iteratively replacing $3$ chocolates in a box of $10$

In the fridge there is a box containing 10 expensive high quality Belgian chocolates, which my mum keeps for visitors. Every day, when mum leaves home for work, I secretly pick 3 chocolates at random, I eat them and replace them with ordinary cheap ones, that have exactly the same wrapping. On the next day I do the same, obviously risking to eat also some of the cheap ones. How many days on average will it take for the full replacement of the expensive chocolates with cheap ones?

I would say $$10/3$$ but this is very simplistic. Also, the total number of ways to pick 3 chocolates out of 10 is $$\binom {10} 3=\frac {10!}{3!7!} = 120$$ which means that after 120 days I will have replaced all chocolates but I don't think it is correct.

Any help?

• It would take at least $4$ days to replace all $10$ chocolates, so the average cannot be $10/3$ (since that's less than $4$). – Minus One-Twelfth Jun 29 '19 at 22:03
• It looks like your problem is a "coupon collector's problem" where you select coupons in "batches" (i.e. select multiple coupons at a time rather than selecting just one coupon at a time). This page may be of interest to you, as well as searching online for things like "coupon collector's problem multiple coupons". – Minus One-Twelfth Jun 29 '19 at 22:07
• It's about $9.045127$, based on my simulation with python. – Botond Jun 29 '19 at 22:21
• I can't solve it analytically, but I was interested in the result so I simulated it :) – Botond Jun 29 '19 at 22:25
• +1 for the question. I can't think when a natural problem, amenable to backwards induction, caused me so much programming grief. That's the sign of a good problem, I figure. – lulu Jun 30 '19 at 0:30

Sketch (not a complete solution, but a road map towards one):

We proceed recursively. Let $$E_i$$ denote the expected number of day it takes given that you have exactly $$i$$ good ones left. The answer you want is $$E_{10}$$.

Let's compute $$E_1$$, for example. Each day the good one gets selected with probability $$\frac {3}{10}$$. Thus $$E_1=\frac {10}3$$

Now let's consider $$E_2$$. The first draw you gets $$0,1$$ or $$2$$ of the good ones. We quickly deduce that $$\binom {10}3 \times E_2=\binom 83\times (E_2+1)+\binom 82\times \binom 21 \times (E_1+1)+ \binom 81\times \binom 22\times (1)$$

This resolves to $$E_2=\frac {115}{24}\approx 4.7917$$

Similarly, for $$i>2$$ $$\binom {10}3\times E_i=\binom {10-i}3\times (E_i+1)+\binom {10-i}2\times \binom i1\times (E_{i-1}+1)$$ $$+\binom {10-i}1\times \binom i2\times (E_{i-2}+1)+\binom {10-i}0\times \binom {i}3\times (E_{i-3}+1)$$

And we can solve the system step by step. The computation requires a bit of attention since in the recursion we must define $$\binom nm=0$$ when $$n.

In the end I get about $$\boxed {9.046}$$ but as the comments will clearly indicate, a great many careless errors were made en route so I advise checking carefully.

• I don't understand your $\ldots$ dots. There are only four terms, aren't there? – TonyK Jun 29 '19 at 23:34
• Yes, the dots are very misleading. Also, the expression after "We quickly deduce" is not compatible with the expression after "Similarly": the first has $8, 7, 6$ where the second simply has $10-i$. – TonyK Jun 29 '19 at 23:38
• @PeterForeman I have made so many errors in the program that I gladly defer to your computation. My practice is "after $3$ idiotic errors, set the thing aside for a day". I passed that threshold some while ago. Thanks for doing a better job with it. – lulu Jun 29 '19 at 23:49
• @lulu You computed $E_2$ incorrectly as you wrote $\binom82$ at the end instead of $\binom81$! – Peter Foreman Jun 29 '19 at 23:54
• Both your and @PeterForeman's numerical evaluations have disappeared! I'm curious to see whether they seem to be converging on an answer :-) – TonyK Jun 30 '19 at 0:03

A self-contained, closed-form solution to the problem.

Suppose there are $$n$$ expensive chocolates and each day you secretly pick $$m$$ chocolates for replacement. Let $$T$$ be a random variable denoting the number of days it takes to replace all the expensive chocolates. Observe that given any set $$S$$ of chocolates, the probability that you have managed to avoid $$S$$ in each of the first $$i$$ days is $$\left[\binom{n-|S|}{m}/\binom{n}{m}\right]^i.$$ Indeed, each of the chosen sets of chocolates is independent of the others and has probability $$\binom{n-|S|}{m}/\binom{n}{m}$$ of avoiding the given set $$S$$, since there are $$\binom{n-|S|}{m}$$ $$m$$-element subsets avoiding $$S$$.

Thus, the probability that you have managed to avoid at least one chocolate during the first $$i$$ days is $$\mathbb P(\text{avoid chocolate number 1})+\cdots +\mathbb P(\text{avoid chocolate number n})$$ $$-\mathbb P(\text{avoid chocolates number 1 and 2})-\mathbb P(\text{avoid chocolates number 1 and 3})-\cdots -\mathbb P(\text{avoid chocolates number n-1 and n})$$ $$\text{+ similar sum over all 3 element subsets, etc}.$$ The previous calculation is using the principle of inclusion-exclusion: we start by overshooting the probability, then subtract off the double counting, then add back in the triple counting, and so on until we have nailed down the exact probability we are after.

Notice that each of the probabilities in the previous display depends only on the number of elements in the set. Thus, we can group them up by number of elements, noting that we will have $$\binom{n}{s}$$ sets contributing to the term for the $$s$$-element subsets. Therefore, $$\mathbb P(T>i)=\sum_{s=1}^n\binom{n}{s}\cdot (-1)^{s+1}\left[\binom{n-s}{m}/\binom{n}{m}\right]^i.$$

Using the tail sum for expectation formula which states that $$\mathbb ET=\sum_{i=0}^{\infty}\mathbb P(T>i),$$ we obtain that $$\mathbb ET=\sum_{i=0}^{\infty}\sum_{s=1}^n\binom{n}{s}(-1)^{s+1}\left[\binom{n-s}{m}/\binom{n}{m}\right]^i.$$ All that remains is to simplify this expression. Manipulating infinite alternating sums can be subtle, but in this we can rearrange the order of summation since the series converges absolutely. Therefore $$\mathbb ET=\sum_{s=1}^n\binom{n}{s}(-1)^{s+1}\sum_{i=0}^{\infty}\left[\binom{n-s}{m}/\binom{n}{m}\right]^i=\sum_{s=1}^n\binom{n}{s}(-1)^{s+1}\frac{1}{1-\binom{n-s}{m}/\binom{n}{m}},$$ where we have applied the geometric series formula $$\sum_{i=0}^{\infty}x^i=\frac{1}{1-x}$$ in the second equality.

Finally, we can multiply through in the fraction to obtain that $$\boxed{\mathbb ET=\binom{n}{m}\sum_{s=1}^n\frac{(-1)^{s+1}\binom{n}{s}}{\binom{n}{m}-\binom{n-s}{m}}}.$$

Substituting $$n=10$$ and $$m=3$$, we can evaluate the formula using this python snippet:

from math import factorial
from fractions import Fraction

def binom(n,k): # binomial coefficient n choose k
return 0 if n < k else factorial(n)/factorial(k)/factorial(n-k)

def e(n,m): # expected number of days to replace all n chocolates if replaced m at a time
return binom(n,m)*sum(Fraction(pow(-1,s+1)*binom(n,s),binom(n,m)-binom(n-s,m)) for s in range(1,n+1))

e(10,3) # Fraction(8241679, 911064), approximately 9.05 days


This is the general analytical solution to the problem. To see all the missteps and sources that inspired this solution, keep reading below for the previous iterations of this answer that led to here.

See the bottom of this post for a previous incarnation of this answer, which cited an incorrect formula and was therefore incorrect.

This problem has appeared many times in the literature. A closed-form solution is known and is valid with $$3$$ and $$10$$ replaced with arbitrary integers $$m$$ and $$n$$ satisfying $$0\leq m\leq n$$.

Theorem 2 in Stadje  (specifically equation (2.15) therein) with $$p=1$$ and $$l=s=n$$ states that the desired expectation equals $$\binom{n}{m}\sum_{j=0}^{n-1}\frac{(-1)^{n-j+1}\binom{n}{j}}{\binom{n}{m}-\binom{j}{m}}.$$ In this case $$n=10$$ and $$m=3$$ and therefore the expected number of days equals $$\frac{8241679}{911064}\approx 9.05.$$

A very nice and self-contained derivation of this formula appears at mathoverflow, although beware that the answer containing this derivation has a minor mistake, which I have corrected both in this post and over at the linked mathoverflow post as well.

Original answer, which cited an incorrect formula from mathoverflow. I have corrected that formula and posted an explanation of the mistake at mathoverflow.

This problem has been solved (with $$3$$ and $$10$$ replaced $$b$$ and $$n$$, respectively) over at math overflow (and in the classical literature). The general formula for the expected number of batches required is $$\sum_{s=1}^{n-b}\frac{(-1)^{s+1}\binom{n}{s}}{1-\binom{n-s}{b}/\binom{n}{b}}.$$ Plugging in $$n=10$$ and $$b=3$$ yields $$\frac{\binom{10}{1}}{1-\binom{9}{3}/\binom{10}{3}}-\frac{\binom{10}{2}}{1-\binom{8}{3}/\binom{10}{3}}+\frac{\binom{10}{3}}{1-\binom{7}{3}/\binom{10}{3}}-\cdots +\frac{\binom{10}{7}}{1-\binom{3}{3}/\binom{10}{3}},$$ which evaluates to $$\frac{41039983}{911064}\approx 45.$$ Compared to the other answer I see we agree in the denominator but not in the numerator...

[EDIT:] In fact, the incorrect answer is exactly $$36$$ greater than the correct answer. The incorrect answer equals $$\binom{10}{3}\sum_{j=3}^{9}\frac{(-1)^{11-j}\binom{10}{j}}{\binom{10}{3}-\binom{j}{3}}$$ whereas the correct answer equals $$\binom{10}{3}\sum_{j=0}^{9}\frac{(-1)^{11-j}\binom{10}{j}}{\binom{10}{3}-\binom{j}{3}}$$ and the difference is $$\binom{10}{3}\sum_{j=0}^{2}\frac{(-1)^{11-j}\binom{10}{j}}{\binom{10}{3}-\binom{j}{3}}=-\binom{10}{0}+\binom{10}{1}-\binom{10}{2}=-1+10-45=-36.$$

• I understand the answer on mathoverflow but this value cannot be correct as no simulation I have run takes anywhere near that number of days. – Peter Foreman Jun 30 '19 at 0:21
• Even if you do one at a time you only get to $\frac {100}3$ which is a lot lower. – lulu Jun 30 '19 at 0:23
• Yes, I am looking for the mistake... – pre-kidney Jun 30 '19 at 0:24
• Looking for the mistake is a recurrent theme for this question! – TonyK Jun 30 '19 at 0:27
• @PeterForeman I have performed a literature search and found the correct formula, which has a closed form and agrees with the other answers in this thread now. – pre-kidney Jun 30 '19 at 1:46

Following on from @lulu I wrote the following program to compute the expected number of days required to take $$n$$ 'expensive chocolates' given that you can only take $$k$$ in one day:

from fractions import Fraction

#n - number of 'expensive chocolates' at start
#g - number of 'expensive chocolates' in current state of recursion
#k - number of chocolates taken each turn
def expected(n,g,k):
if g<=0:
return Fraction()
if k>=n:
return Fraction(1,1)
ex=ncr(n-g,k)
for i in range(1,k+1):
ex+=ncr(n-g,k-i)*ncr(g,i)*(expected(n,g-i,k)+Fraction(1,1))
return ex/(ncr(n,k)-ncr(n-g,k))

def ncr(n,r):
if r>n:
return Fraction()
c=Fraction(1,1)
while r>0:
c*=Fraction(n,r)
n-=1
r-=1
return c


For example, running expected(10,10,3) gives Fraction(8241679, 911064) which is equivalent to $$\frac{8241679}{911064}\approx9.046212999\text{ days}$$ for the case in question.

• With $100,000$ I got $9.06933$ and $1.96$ standard deviations came out to about $.0206$ so we agree. – saulspatz Jun 30 '19 at 1:22

This is just a numerical supplement to @lulu's nice answer. Since I wanted to better see the development of $$E_k, 1\leq k\leq 10$$, I did the calculation somewhat detailed. At the end a recurrence relation is given.

We denote with $$E_k$$ the expected number of days it takes given that we have exactly $$k$$ good ones left and obtain following @lulu's approach:

\begin{align*} \binom{10}{3}E_1&=\binom{9}{3}\left(E_1+1\right)+\binom{9}{2}\binom{1}{1}\\ \binom{10}{3}E_2&=\binom{8}{3}\left(E_2+1\right)+\binom{8}{2}\binom{2}{1}\left(E_1+1\right) +\binom{8}{1}\binom{2}{2}\\ \binom{10}{3}E_3&=\binom{7}{3}\left(E_3+1\right)+\binom{7}{2}\binom{3}{1}\left(E_2+1\right) +\binom{7}{1}\binom{3}{2}\left(E_1+1\right)+\binom{7}{0}\binom{3}{3}\\ \binom{10}{3}E_4&=\binom{6}{3}\left(E_4+1\right)+\binom{6}{2}\binom{4}{1}\left(E_3+1\right) +\binom{6}{1}\binom{4}{2}\left(E_2+1\right)+\binom{6}{0}\binom{4}{3}\left(E_1+1\right)\tag{1}\\ \binom{10}{3}E_5&=\binom{5}{3}\left(E_5+1\right)+\binom{5}{2}\binom{5}{1}\left(E_4+1\right) +\binom{5}{1}\binom{5}{2}\left(E_3+1\right)+\binom{5}{0}\binom{5}{3}\left(E_2+1\right)\\ \binom{10}{3}E_6&=\binom{4}{3}\left(E_6+1\right)+\binom{4}{2}\binom{6}{1}\left(E_5+1\right) +\binom{4}{1}\binom{6}{2}\left(E_4+1\right)+\binom{4}{0}\binom{6}{3}\left(E_3+1\right)\\ \binom{10}{3}E_7&=\binom{3}{3}\left(E_7+1\right)+\binom{3}{2}\binom{7}{1}\left(E_6+1\right) +\binom{3}{1}\binom{7}{2}\left(E_5+1\right)+\binom{3}{0}\binom{7}{3}\left(E_4+1\right)\\ E_{10}&=E_7+1\\ \end{align*}

The terms in each of the lines above which are not multiplied by $$E_k$$ can be simplified using Vandermonde's identity. Taking for instance the line (1) we have \begin{align*} \binom{10}{3}E_4&=\binom{6}{3}\left(E_4+1\right)+\binom{6}{2}\binom{4}{1}\left(E_3+1\right) +\binom{6}{1}\binom{4}{2}\left(E_2+1\right)+\binom{6}{0}\binom{4}{3}\left(E_1+1\right)\\ &=\sum_{j=0}^3\binom{6}{3-j}\binom{4}{j}E_{4-j}+\color{blue}{\sum_{j=0}^3\binom{6}{3-j}\binom{4}{j}}\\ &=\sum_{j=0}^3\binom{6}{3-j}\binom{4}{j}E_{4-j}+\color{blue}{\binom{10}{3}}\\ \binom{10}{3}\left(E_4-1\right)&=\sum_{j=0}^3\binom{6}{3-j}\binom{4}{j}E_{4-j}\tag{2} \end{align*}

Setting $$E_0=1$$ we obtain the recurrence relation according to (2) \begin{align*} \color{blue}{\binom{10}{3}\left(E_k-1\right)}&\color{blue}{=\sum_{j=0}^{\min\{3,k\}}\binom{10-k}{3-j}\binom{k}{j}E_{k-j}\qquad\qquad 1\leq k\leq 7}\\ \color{blue}{E_0}&\color{blue}{=1}\\ \color{blue}{E_{10}}&\color{blue}{=E_7+1} \end{align*}

We obtain solving the recurrence relation: \begin{align*} E_1&=\frac{10}{3}\approx 3.333\\ E_2&=\frac{115}{24}\approx 4.792\\ E_3&=\frac{787}{136}\approx5.787\\ E_4&=\frac{6\,661}{1\,020}\approx 6.530\\ E_5&=\frac{15\,999}{2\,244}\approx 7.125\\ E_6&=\frac{330\,641}{43\,384}\approx 7.621\\ E_7&=\frac{7\,330\,615}{911\,064}\approx 8.046\\ \color{blue}{E_{10}}&\color{blue}{=\frac{8\,241\,679}{911\,064}\approx 9.046} \end{align*}

• All solutions are very nice - I had to pick one so I gave the green mark to @lulu. Thank you all! – Sal.Cognato Jul 1 '19 at 17:18
• @Sal.Cognato: Good decision and thanks for the nice comment. – Markus Scheuer Jul 1 '19 at 17:32

HINT for a Markov approach

If at a certain point you have $$A$$ ordinary and $$B$$ belgian chocolates, when you select the three chocolates you can select :
- $$[a,a,a]$$ (three ordinary ch.) with prob. $${A \over {A + B}}{{A - 1} \over {A - 1 + B}}{{A - 2} \over {A - 2 + B}} = {{A^{\,\underline {\,3\,} } } \over {\left( {A + B} \right)^{\,\underline {\,3\,} } }}$$
- $$[a,a,b]$$ with prob. $${A \over {A + B}}{{A - 1} \over {A - 1 + B}}{B \over {A - 1 + B}} = {{A^{\,\underline {\,2\,} } B^{\,\underline {\,1\,} } } \over {\left( {A + B} \right)^{\,\underline {\,2\,} } \left( {A - 1 + B} \right)}}$$
- $$[a,b,a]$$ with prob. $${A \over {A + B}}{B \over {A - 1 + B}}{{A - 1} \over {A - 1 + B - 1}} = {{A^{\,\underline {\,2\,} } B^{\,\underline {\,1\,} } } \over {\left( {A + B} \right)^{\,\underline {\,3\,} } }}$$
and so on, obtaining the development of $$\left( {A + B} \right)^{\,\underline {\,3\,} } = \sum\limits_j {\left( \matrix{ 3 \hfill \cr j \hfill \cr} \right)A^{\,\underline {\,3 - j\,} } B^{\,\underline {\,j\,} } }$$ The falling factorial at the numerator will ensure that we are not going to pick more pieces than available, and we shall ensure that the fraction be null, whenever the numerator is null.

At the end, we are going to restore three ordinary chocolates, i.e. bring $$A$$ to $$A+3$$ and restore the total to the original $$N=A+B$$

It is therefore possible to write the probabilities that after a cycle of picking and restoring the number of ordinary ch. passes from $$A$$ to $$A+0,\, A+1, \, A+2, \, A+3$$ (with $$B$$ being the complement to $$N$$).
That means that we can set up a 4-diagonal Markov Matrix for $$A$$, or else for $$B$$, and study the characteristics of that.

Performing the calculations as above the $$A$$-th row of the transition matrix will be $$\begin{array}{*{20}c} {} & | & \cdots & A & {A + 1} & {A + 2} & {A + 3} \\ \hline A & | & \ddots & {\frac{{A^{\,\underline {\,3\,} } }}{{N^{\,\underline {\,3\,} } }}} & {3\frac{{A^{\,\underline {\,2\,} } \left( {N - A} \right)^{\,\underline {\,1\,} } }}{{N^{\,\underline {\,3\,} } }}} & {3\frac{{A^{\,\underline {\,1\,} } \left( {N - A} \right)^{\,\underline {\,2\,} } }}{{N^{\,\underline {\,3\,} } }}} & {\frac{{\left( {N - A} \right)^{\,\underline {\,3\,} } }}{{N^{\,\underline {\,3\,} } }}} \\ \end{array}$$

Example

To contain the dimensions of the matrix, let's consider the case $$N=4$$, then the transition matrix is $$T\left( 4 \right) = \frac{1}{4} \, \left( {\matrix{ 0 & 0 & 0 & 1 & 0 \cr 0 & 0 & 0 & 3 & 1 \cr 0 & 0 & 0 & 2 & 2 \cr 0 & 0 & 0 & 1 & 3 \cr 0 & 0 & 0 & 0 & 1 \cr } } \right)$$