# The coupon collector's earworm

[EDITED mostly to report on the answer by Kevin Costello (and to improve the gp code at the end)]

I thank Nicolas Dupont for the following question (and for permission to disseminate it further):

I have a playlist with, say, $N$ pieces of music. While using the shuffle option (each such piece is played randomly at each step), I realized that, generally speaking, I have to hear quite a lot of times the same piece before the last one appears. It makes me think of the following question:

At the moment the last non already heard piece is played, what is the max, in average, of number of times the same piece has already been played?

I have not previously enountered this variant of the coupon collector's problem. If it is new, then (thinking of the original real-world context origin of the problem) I propose to call it the "coupon collector's earworm".

Is this in fact a new question? If not, what is known about it already?

Let us call this expected value $M_N$. Nicolas observes that $M_1=1$ and $M_2=2$ (see below), and asks:

I doubt there is such an easy formula for general $N$ - would you be able to find some information on it, e.g. its asymptotic behavior?

[Nicely answered by Kevin Costello: $M_N$ is asymptotic to $e \log N$ as $N \rightarrow \infty$. Moreover, the maximal multiplicity is within $o(\log N)$ of $e \log N$ with probability $\rightarrow 1$. I don't recall any other instance of a naturally-arising asymptotic growth rate of $e \log N$...]

Recall that the standard coupon collector's problem asks for the expected value of the total number of shuffles until each piece has appeared at least once. It is well known that the answer is $N H_N$ where $H_N = \sum_{i=1}^N 1/i$ is the $N$-th harmonic number. Hence the expected value of the average number of plays per track is $H_N$, which grows as $\log N + O(1)$. The expected maximum value $M_N$ must be at least as large $-$ in fact larger once $N>1$, because one track is heard only once, so the average over the others is $(N H_N-1) / (N-1)$. One might guess that $M_N$ is not that much larger, because typically it's only the last few tracks that take most of shuffles to appear. But it doesn't seem that easy to get at the difference between the expected maximum and the expected average, even asymptotically. Unlike the standard coupon collector's problem, here an exact answer seems hopeless once $N$ is at all large (see below), so I ask:

How does the difference $M_N - H_N$ behave asymptotically as $N \rightarrow \infty$?

[Looks like this was a red herring: "One might guess" plausibly that $M_N - H_N$ is dwarfed by $H_N$ for large $N$, but by Kevin Costello's answer $M_N - H_N$ asymptotically exceeds $H_N$ by a factor $e - 1$, and that factor is more complicated than $\lim_{N\rightarrow\infty} M_N/H_N = e$, so analyzing the difference $M_N-H_N$ is likely not a fruitful approach.]

@ For each $N>1$, the expected value of the maximum count is given by the convergent $(N-1)$-fold sum $$M_N = \sum_{a_1,\ldots,a_{N-1} \geq 1} N^{-\!\sum_i \! a_i} \Bigl(\sum_i a_i\Bigr)! \frac{\max_i a_i}{\prod_i a_i!}.$$ Indeed, we may assume without loss of generality that the $N$-th track is heard last; conditional on this assumption, the probability that the $i$-th track will be heard $a_i$ times for each $i<N$ is $N^{-\!\sum_i \! a_i}$ times the multinomial coefficient $(\sum_i a_i)! \big/ \prod_i a_i!$. Numerically, these values are $$2.00000, \quad 2.84610+, \quad 3.49914-, \quad 4.02595\!-$$ for $N=2,3,4,5$.

@ A closed form for $M_N$ is available for $N \leq 3$ and probably not beyond. Trivially $M_1 = 1$; and N.Dupont already obtained the value $M_2 = 2$ by evaluating $M_2 = \sum_{a \geq 1} a/2^a$. But for $N=1$ and $N=2$ the problem reduces to the classical coupon collector's problem. Already for $N=3$ we have a surprise: $M_3 = 3/2 + (3/\sqrt{5})$, which has an elementary but somewhat tricky proof. For $N=4$, I get $$M_4 = \frac73 - \sqrt{3} + \frac4\pi \int_{x_0}^\infty \frac{(2x+1) \, dx}{(x-1) \sqrt{4x^3-(4x-1)^2}}$$ where $x_0 = 3.43968\!+$ is the largest of the roots (all real) of the cubic $4x^3-(4x-1)^2$. I don't expect this to simplify further: the integral is the period over an elliptic curve of a differential with two simple poles that do not differ by a torsion point. In general one can reduce the $(N-1)$-fold sum to an $(N-2)$-fold one (which is one route to the value of $M_3$ and $M_4$), or to an $(N-3)$-fold integral, but probably not beyond.

@ It's not too hard to simulate this process even for considerably larger $N$. In GP one can get a single sample of the distribution from the code

try(N) = v=vector(N); while(!vecmin(v),v[random(N)+1]++); vecmax(v)


[turns out that one doesn't need to call vecmin each turn:

try(N, m,i)= v=vector(N); m=N; while(m, i=random(N)+1; v[i]++; if(v[i]==1,m--)); vecmax(v)


does the same thing in $\rho+O(1)$ operations per shuffle rather then $\rho+O(N)$, where $\rho$ is the cost of one random(N) call.]

So for example

sum(n=1,10^4,try(100)) / 10000.


averages 1000 samples for $M_{100}$; this takes a few seconds, and seems to give about $11.7$.

• I feel like I've seen this tackled in discussions of the Coupon Collector's problem before, but I can't remember specific details - have you checked The Art Of Computer Programming? That's the first place I'd look (and in particular, in the relevant exercises) for a result... – Steven Stadnicki Aug 26 '15 at 5:02
• Related is Schteuzow's "Asymptotics for the maximum in the coupon collector's problem", Math. Sci. 27 (2002), 85--90. MR1950287. But there the computation is giving $N$ dollars to $N$ people, and enquiring concerning the expected maximum wealth. The technique of "Poissonization" should be relevant though. – kantelope Aug 26 '15 at 16:07
• In sampling $M_{1000}$, it took 8 minutes for 10000 trials and obtained about $17.8$, which is approximately $\log(1000)/\log(100)$ bigger than the $M_{100}$ sample. – kantelope Aug 26 '15 at 17:29
• @Steven Stadnicki - I hadn't checked The Art Of Computer Programming; that's a good idea, but I cannot locate the problem in either the main text or the exercises. Were you able to find it there? – Noam D. Elkies Aug 31 '15 at 17:51
• @kantelope try the new faster code. – Noam D. Elkies Aug 31 '15 at 17:51

For the asymptotic case: Let $t_1=n \log n - Cn$ and $t_2 = n \log n + Cn$, where $C$ is slowly tending to infinity. It is a classic result that as $C$ tends to infinity the probability all coupons are collected at time $t_1$ tends to $0$, and the probability all coupons are collected at time $t_2$ tends to $1$.

So the number being asked for in the original post is with high probability sandwiched between the most collected coupon at time $t_1$ and the most collected coupon at time $t_2$. This has been studied a fair amount, though often in language of load-balancing and balls-in-bins ("if you toss $m$ balls in $n$ bins, what is the typical number of balls in the bin with the most balls"). In particular, there's a nice analysis of Raab and Steger that uses the second moment method to give a tight concentration on this. It follows from Theorem $1$ in their paper that at time $c n \log n$, the most common coupon has almost surely been collected $(d_c+o(1)) \log n$ times, where $d_c$ is the larger real root of $$1+x(\log c - \log x +1)-c=0$$

In our case, we have $d_1=e$, so the most common song will have been heard $(e+o(1))\log n = (e+o(1))H_n$ times.

In response to the comment: The first moment part of the calculation also comes with come concentration estimates for free. At time $t_1$, the probability that we've seen some coupon at least $C \log n$ times can be bounded by \begin{eqnarray*} & & \sum_{j=C \log n}^{t_1} n \binom{t_1}{j} n^{-j} (1-\frac{1}{n})^{t_1-j} \\ &\leq& \sum_{j=C \log n}^{t_1} n \left(\frac{e t_1}{ nj}\right)^j (1-\frac{1}{n})^{n \log n + o(\log n)} (1-\frac{1}{n})^{-j} \\ &=& n^{o(1)} \sum_{j=C \log n}^{t_1} \left(\frac{e \log n(1+o(1))}{j}\right)^j \\ &\leq& n^{o(1)} \sum_{j=C \log n}^{t_1} \left(\frac{e}{C}+o(1)\right)^j \end{eqnarray*}

Which should decay rapidly enough in $C$ to guarantee the mean lies close to the concentration.

The linked bounds for the probability of the completion time lying outside $[t_1, t_2]$ also decay quite rapidly (e.g. cardinal's answer, and David's subsequent comments on it).

• Thanks for this analysis; the appearance of $e$ as the asymptotic coefficient of $\log N$ is unexpected, and the results you quote almost give the asymptotic behavior of $M_N$ en route to describing the distribution near $e \log N$. I write "almost" because the distribution could be asymptotically concentrated at $e \log N$ but still have a heavy enough tail to shift the mean $M_N$ above $(e+o(1)) \log N$. Usually such "concentration" results come with error estimates that force the mean and mode to be quite close; is this the case here too? – Noam D. Elkies Aug 29 '15 at 14:30
• Seems so. I've added a bit on the concentration. – Kevin P. Costello Aug 29 '15 at 18:27
• Thanks for this as well. I expect that I'll "accept" your answer soon (I already "up-voted" it), and will edit the question to reflect your answer. – Noam D. Elkies Aug 29 '15 at 18:36

Let $X$ be the total number of plays and let $Y$ be the number of plays for the track with the most plays.

I think a good strategy is to estimate

$$E [ Y- X | X = X_0 ]$$

It seems like this involves a simpler problem, and that this will be very flat as a function of $X$. So if we get a good asymptotic for this, we probably don't need to know much about the probability distribution of $X$ to calculate your expectation.

What is the distribution of $Y$ conditional on a certain value of $X$? Well it's just the random model where $X_0-1$ balls fall into $N-1$ bins and you take the max of the number of balls in each bin, conditional on each bin having at least one ball.

Now we want to prove some bounds relating $E [ Y- X | X = X_0 ]$ to $\sqrt{X_0/N-1} \sqrt{\log (N-1)}$ and $\sqrt{X_0/N } \sqrt{\log (N-1)}$ and stuff. I'm just thinking of the number of balls in a fixed bin as approximately Gaussian, and these Gaussians as approximately independent, to get these kinds of figures. I'lll try to actually justify them later or at least get an upper bound.

If we get some asymptotic for $E [ Y- X | X = X_0 ]$ that is roughly linear in $\sqrt{X_0}$, then we only have to estimate $E [ \sqrt{X} ]$ which shouldn't be incredibly hard because that function is just barely concave so it should be quite close to $\sqrt{ E[X] }$ .