# Plotting an Unreasonable Function

Without getting into too much detail, the following (very complicated) function recently appeared as a solution to a combinatorics problem I've been thinking about:

$$P(n) = \frac{52!}{52^{52}} \cdot \sum_{1=i_{1} < i_{2} < \ldots < i_{51} < i_{52} \le n} \left[ \prod_{k=1}^{51} \left( \frac{k}{52} \right)^{i_{k+1}-i_{k}-1}\right]$$

I'd like to plot this sucker in Mathematica but as it stands I don't see how that's going to happen. If there's a way to simplify this expression or some sort of efficient way to evaluate it I'd really appreciate any help.

I know that $P(n)=0$ for $1 \le n < 52$, $P(52)=\frac{52!}{52^{52}}$, and $P(n)$ should grow quickly towards 1 and then plateau, so it should look something like the CDF of the normal distribution centered around 52 when all's said and done.

Thanks!

• There's no way this can look like a CDF centered around $52$ when it only starts at $52$! – whuber Nov 6 '12 at 5:55
• Whoops, of course! I meant to say that I expected it to be centered around 100 (for no good reason). It looks like it's really centered around 208. I'm happy to see that my intuition about the general shape of the curve was surprisingly accurate. Thanks so much Sasha and whuber, still reading your answers but they look awesome so far. If you guys are interested this is basically the solution to the following problem: drawing from an infinite deck of cards, what is the probability that one will obtain a complete deck after n draws. – Jackson Nov 6 '12 at 6:15
• This is known as the Coupon Collector's Problem. The generating function is derived at en.wikipedia.org/wiki/…. For some asymptotic estimates see stats.stackexchange.com/questions/7774/…. The distribution is decidedly non-Normal, by the way: plot the $Q(n)$ to see how skewed it is. – whuber Nov 6 '12 at 6:20

Letting $j_k = i_{k+1}-i_k-1$ and writing $$Q(n) = P(n) - P(n-1) = C\sum_{0 \le j_1, j_2, \cdots, j_{51}\vert j_1+\cdots+j_{51}=n-52} \prod_{k=1}^{51}\left(\frac{k}{52}\right)^{j_k}\,,$$ with $C$ a constant, exhibits the $P(n)$ as cumulative sums of the $Q(n)$ and shows that $Q(n)$ is the coefficient of $x^{n-52}$ in the formal power series

$$q(x) = \frac{52!}{52^{52}} \prod_{k=1}^{51} \frac{1}{1 - \frac{k x}{52}}.$$

A Mathematica implementation (using a generic variable m for $52$) is

q[x_, m_] := Product[1/(1 - k x / m), {k, 1, m - 1}] m! / m^m;
p = Accumulate[CoefficientList[Normal[Series[q[x, 52.], {x, 0, 640}]], x]]


($0.0156$ seconds).

Here is a plot:

ListPlot[p, DataRange -> {52, Length[p] + 51}] • Congrats. I was getting on to same ideas, but you figured it out first. Notice that q[x,d]==d!/FactorialPower[d,d,x]. This shows that $q[1,d]$ equals 1 exactly. Of course it follows from the explicit product representation readily. – Sasha Nov 6 '12 at 6:02

Let's replace 52 with $d$. Then we are seeking to compute $$p_d(n) = \frac{(d-1)!}{d^{d-1}} \sum_{1 = i_1 < i_2 < \cdots < i_{d-1} < i_{d} \leqslant n} \left[ \prod_{k=1}^{d-1} \left( \frac{k}{d} \right)^{i_{k+1}-i_k-1} \right] \tag{1}$$ The set $\{i_2, \ldots,i_d\}$ is a length $d-1$ subset of natural consecutive numbers from 2 to $n$. The number such possible subsets is given by binomial expression: $$S_d(n) = \binom{n-1}{d-1}$$ The number of such subsets, and thus the number of combinations over which we need to sum grows rather fast:

In:= With[{d = 52},
Table[{n, Binomial[n - 1, d - 1]}, {n, d, d + 5}]]

Out= {{52, 1}, {53, 52}, {54, 1378}, {55, 24804}, {56,
341055}, {57, 3819816}}


Such a counting allows to build reasonable expectation about what is computable.

Here is the code implementing the $p_d(n)$:

Clear[p];
p[d_Integer, n_Integer] := Module[{set, i},
If[n < d, 0,
set = Subsets[Range[2, n], {d - 1}];
Total[(Function @@ {Product[(k/d)^(i[k] - i[k - 1] - 1), {k, 1,
d - 1}] /. {i -> 1} /. i -> Slot}) @@@ set] d!/d^d
]]


Here is how the sequence $p_d(n)$ looks for $d=8$: Now we can turn to the $d=52$:

In:= Table[p[52, n + 51]/(52!/52^52), {n, 1, 6}]

Out= {1, 53/2, 74889/208, 1390455/416, 13409133507/562432, \
156710118411/1124864}


As you can see $p_{52}(57) \approx 6 \cdot 10^{-17}$, that is quite small, and making progress on the problem requires more theoretical insight. Introducing $j_k = i_{k+1}-i_{k}-1$ we can rewrite ${\rm eq. } (1)$ as follows: $$p_d(n) = \frac{(d-1)!}{d^{d-1}} \sum_{j_1=0}^{n-d} \sum_{j_{2}=0}^{n-d} \cdots \sum_{j_{d-1}=0}^{n-d} \left( [j_1+\cdots+j_{d-1} \leqslant n-d ] \prod_{k=1}^{d-1} \left[\left(\frac{k}{d}\right)^{j_k} \right] \right)$$where $[j_1+\cdots+j_{d-1} \leqslant n-d ]$ stands for Iverson bracket. Introducing $j_d = n - d - \sum\limits_{k=1}^{d-1} j_k$ we rewrite: p_d(n) = \frac{(d-1)!}{d^{n-1}} \sum_{\begin{align}{c} j_1 \geqslant 0, \ldots, j_d \geqslant 0 \cr j_1 + \cdots + j_d = n-d \end{align}} \prod_{k=1}^d k^{j_k} leading to @whuber's answer, who beat me to it.

• I think you figured it out first, but worked harder to provide a detailed answer (+1). I was so glad to see our answers almost agree!--it suggests the method is good. Some details still need to be worked out, though. E.g., the constant term in your equation (1) doesn't seem to agree with the OP's constant term and that may slightly have affected the constant term in your final expression. – whuber Nov 6 '12 at 6:08
• It's more basic than that: I think the numerator and denominator got reversed in the first equation... – whuber Nov 6 '12 at 6:24
• Thanks! I have edited to correct the factor. – Sasha Nov 6 '12 at 6:28