Iteratively replacing $3$ chocolates in a box of $10$
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<m$.
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.
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 [1990] (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. $$
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.
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*}
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) $$