Is there a closed form expression for $\sum_{1\le i_1 \le i_2 \le \dots \le i_k \le n} i_1 i_2 \dots i_k$? These sums showed up in a probability problem I was working on. They're not quite the Stirling numbers of the first kind since it's possible to have e.g. $i_1 = i_2$. Denoting the sum by $(k\mid n)$ we have the recurrence relation
$(k\mid n) = n(k-1\mid n) + (k\mid n-1)$
I thought maybe the average term of the sum would be similar to the average product over a random list of $k$ numbers from $1$ to $n$. But it seems to always be a few times greater.
I'm beginning to despair of there being a closed form expression for $(k\mid n)$ in terms of factorials and powers. Does anyone know how to analyze the sum further and perhaps find a good approximation?
 A: Here I will show by induction that the Stirling numbers of the second kind give the general answer for the multisubset sum: with $Q_{n,k}$ being the set of all $k$-multisubsets of $[n]=\{1,\dots,n\}$,
$$\sum_{s\in Q_{n,k}}\prod_{i=1}^ks_i=S(n+k,n)$$
For $n=1$ the equation holds since the only multisubset in $Q_{1,k}$ is $\{1,1,\dots,1\}$, leading to a sum of $1$, while $S(1+k,1)=1$.
For $k=1$ the equation holds because the multisubset sum reduces to a sum of the numbers in $[n]$:
$$\sum_{i=1}^ni=\binom{n+1}2=S(n+1,n)$$
Now for $n,k>1$ the multisubsets in $Q_{n,k}$ can be separated into those with the largest possible element $n$ and those without it. The first subset of multisubsets is $Q_{n,k-1}$ with $n$ added to each multisubset; the second subset is $Q_{n-1,k}$. Thus
$$\sum_{s\in Q_{n,k}}\prod_{i=1}^ks_i=n\sum_{s\in Q_{n,k-1}}\prod_{i=1}^ks_i+\sum_{s\in Q_{n-1,k}}\prod_{i=1}^ks_i$$
Assuming the relation is true for $(n,k-1)$ and $(n-1,k)$ – inducting on $n+k$:
$$\sum_{s\in Q_{n,k}}\prod_{i=1}^ks_i=nS(n+k-1,n)+S(n+k-1,n-1)$$
But $S(n,k)=kS(n-1,k)+S(n-1,k-1)$, so the left-hand side is $S(n+k,n)$ by substituting $(n,k)\leftarrow(n+k,n)$. This proves the induction, so the multinomial sum is $S(n+k,n)$ for all $n,k\ge1$.
A: Parcly has given a proof that
$$
\sum_{1\le i_1\le \dots \le i_k\le n}i_1\cdots i_k={n+k \brace n}\tag{$\star$}
$$
Surely, such a beautiful equation must have a combinatorial proof? Indeed it does!
Given a partition of $\{1,2,\dots,n+k\}$ into $n$ unlabeled parts, define the "leader" of each part to be its smallest element. Every other element of $\{1,\dots,n+k\}$ is a "follower." Let $F_1,\dots,F_k$ be the list of followers, sorted in increasing order, so $F_1<\dots<F_k$.
Then, define the "signature" of a partition to be the sequence $(i_1,\dots,i_k)$, where for each $r\in \{1,\dots,k\}$, $i_r$ is the number of leaders which are less than $F_r$. Here is an example:
$$
\color{red}{1},2,5,6 \mid \color{red}{3}, 9 \mid \color{red}{4}\mid  \color{red}7,8
$$
The leaders are in red. The signature of this partition is $(1,3,3,4,4)$, because the followers are $2,5,6,8,9$, and there is $1$ leader which is numerically less than $2$ (only the leader $1$ is less than $2$), there are $3$ leaders less than $5$ (namely, $1,3$ and $4$), there are $3$ leaders less than $6$, and there are four leaders less than both $8$ and $9$.
Note that a signature must be a weakling increasing sequence. I claim that for each potential signature $(i_1,i_2,\dots,i_k)$, there are exactly $i_1\times \dots \times i_k$  partitions with that signature. Indeed, imagine choosing a partition with signature $(i_1,\dots,i_k)$ one element at a time, starting at $1$ and proceeding in numerical order.

*

*The signature $(i_1,\dots,i_k)$ completely determines the set of leaders. There is only one way to choose the part for a leader, since they must be put in a new part, and all the currently empty parts are identical.


*When choosing the $r^\text{th}$ follower, it must be placed in a currently nonempty part, with some other leader. There are currently $i_r$ nonempty parts, so there are $i_r$ ways to place this follower.
Summing over all possible signatures, we have proved $(\star)$.
A: We get from first principles the generating function where we mark
multisets containing some number of instances of $q$ giving the factor
being contributed and the exponent of $z$ how often  it occurs:
$$[z^k] \prod_{q=1}^n (1+qz+q^2z^2+q^3z^3+\cdots)
= [z^k] \prod_{q=1}^n \frac{1}{1-qz}
\\ = \frac{1}{n!} [z^k] \prod_{q=1}^n \frac{1}{1/q-z}
= (-1)^n \frac{1}{n!} [z^k] \prod_{q=1}^n \frac{1}{z-1/q}.$$
Now using partial fractions by residues we have
$$(-1)^n \frac{1}{n!} [z^k]
\sum_{p=1}^n \frac{1}{z-1/p}
\prod_{q=1}^{p-1} \frac{1}{1/p-1/q}
\prod_{q=p+1}^n \frac{1}{1/p-1/q}
\\ = (-1)^n \frac{1}{n!} [z^k]
\sum_{p=1}^n \frac{1}{z-1/p}
\prod_{q=1}^{p-1} \frac{pq}{q-p}
\prod_{q=p+1}^n \frac{pq}{q-p}
\\ = (-1)^n \frac{1}{n!} [z^k]
\sum_{p=1}^n \frac{p^{n-1}}{z-1/p}
\prod_{q=1}^{p-1} \frac{q}{q-p}
\prod_{q=p+1}^n \frac{q}{q-p}
\\ = (-1)^n [z^k]
\sum_{p=1}^n \frac{p^{n-2}}{z-1/p}
\prod_{q=1}^{p-1} \frac{1}{q-p}
\prod_{q=p+1}^n \frac{1}{q-p}
\\ = (-1)^n [z^k]
\sum_{p=1}^n \frac{p^{n-2}}{z-1/p}
\frac{(-1)^{p-1}}{(p-1)!} \frac{1}{(n-p)!}
\\ = (-1)^n \frac{1}{n!} [z^k]
\sum_{p=1}^n \frac{p^{n-1}}{z-1/p} {n\choose p} (-1)^{p-1}
\\ = (-1)^n \frac{1}{n!} [z^k]
\sum_{p=1}^n \frac{p^{n}}{pz-1} {n\choose p} (-1)^{p-1}
\\ = (-1)^n \frac{1}{n!} [z^k]
\sum_{p=1}^n \frac{p^{n}}{1-pz} {n\choose p} (-1)^p
\\ = (-1)^n \frac{1}{n!}
\sum_{p=1}^n  {n\choose p} (-1)^p p^{n+k}.$$
We may lower $p$ to zero and obtain
$$\frac{1}{n!}
\sum_{p=0}^n  {n\choose p} (-1)^{n-p} p^{n+k}
= {n+k \brace n}.$$
A: Consider the expansion
$$X^N = \sum_{k=0}^N {N\brace k} (x)_k$$
where $(x)_k = x(x-1)\cdots (x-k+1)$, and ${N\brace k}$ is the Sterling number of the second kind. It follows that for every $n$ between $0$ and $N$
$$\sum_{k=0}^n {N\brace k} (x)_k$$
is the Lagrange interpolation  polynomial for the function $f(x) = x^N$ and nodes $0$, $1$, $\ldots$, $n$.
The leading term ${N \brace n}$ of this polynomial can also be calculated with Cramer's rule as a quotient of determinants
$$\frac{\Delta'(x_0, x_1, \ldots, x_n)}{\Delta(x_0, x_1, \ldots, x_n)}$$
evaluated at $(x_0, x_1, \ldots, x_n) = (0, 1, \ldots, n)$.  Here, $\Delta$ is the usual Vandermonde determinant, while $\Delta'$ is obtained from $\Delta$ by having in last column the exponent $N$ instead of $n$.
Now, it is a fact that the above expression in $x_0$, $\ldots$, $x_n$ is the complete symmetric polynomial of degree $N-n$ in $x_0$, $\ldots$, $x_n$.
We are done.
