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.