Plotting an Unreasonable Function

Question

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!

Answer 1

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}]

Answer 2

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[88]:= With[{d = 52}, 
 Table[{n, Binomial[n - 1, d - 1]}, {n, d, d + 5}]]

Out[88]= {{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[0] -> 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[117]:= Table[p[52, n + 51]/(52!/52^52), {n, 1, 6}]

Out[117]= {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.