While answering this question I encountered the following probability distribution: $$T_n\sim p(t~|~\lambda,n)=n\lambda e^{-\lambda t}(1-e^{-\lambda t})^{n-1}$$ This is a PDF for the random variable $T_n$ , which has to do with waiting times at a bank, or restaurant, or whatever. See the linked question for further details. Anyway, it was of some importance to compute the first and second non-central moments of this random variable. This is detailed in my answer to the linked question. But, then I started wondering about the $m$th non-central moment of this distribution. $$\mathrm{E}({T_n}^m)=n\lambda \int_0^\infty t^m\cdot e^{-\lambda t}(1-e^{-\lambda t})^{n-1}\mathrm{d}t$$ Using a binomial expansion, $$\mathrm{E}({T_n}^m)=n\lambda \int_0^\infty t^m\cdot e^{-\lambda t}\sum_{k=0}^{n-1}{}_{(n-1)}\mathrm{C}_k(-e^{-\lambda t})^{n-1-k}\mathrm{d}t$$ Simplifying and switching the integral and summation signs (fingers crossed!) $$\mathrm{E}({T_n}^m)=n\lambda\sum_{k=0}^{n-1}(-1)^{n-1-k}{}_{(n-1)}\mathrm{C}_k\int_0^\infty t^me^{-(n-k)\lambda t}\mathrm{d}t$$ Now using $\tau=\lambda(n-k)t$: $$\mathrm{E}({T_n}^m)=n\lambda\sum_{k=0}^{n-1}(-1)^{n-1-k}{}_{(n-1)}\mathrm{C}_k\int_0^\infty \left(\frac{\tau}{\lambda(n-k)}\right)^m e^{-\tau}\frac{1}{\lambda(n-k)}\mathrm{d}\tau$$ $$\mathrm{E}({T_n}^m)=\frac{n}{\lambda^m}\sum_{k=0}^{n-1}\frac{(-1)^{n-1-k}{}_{(n-1)}\mathrm{C}_k}{(n-k)^{m+1}}\int_0^\infty \tau^m e^{-\tau}\mathrm{d}\tau$$ Great! This is a familiar integral and we can use the properties of the Gamma function. $$\mathrm{E}({T_n}^m)=\frac{n\cdot m!}{\lambda^m}\sum_{k=0}^{n-1}\frac{(-1)^{n-1-k}{}_{(n-1)}\mathrm{C}_k}{(n-k)^{m+1}}$$ But what can we say about the sum $$S(n,m)=\sum_{k=0}^{n-1}\frac{(-1)^{n-1-k}{}_{(n-1)}\mathrm{C}_k}{(n-k)^{m+1}}$$ For $m=0$, Wolfram finds $$S(n,0)=\frac{\Gamma(n)}{\Gamma(n+1)}=\frac{1}{n}$$ Which is consistent, because since $$\mathrm{E}({T_n}^m)=\frac{n\cdot m!}{\lambda^m}S(n,m)$$ Then $$\mathrm{E}({T_n}^0)=\frac{n\cdot 0!}{\lambda^0}\frac{1}{n}=1.$$ For $m=1$ it finds the striking $$S(n,1)=\frac{H_n}{n}$$ With $H_n$ being the $n$th harmonic number. For $m=2$ it finds the even more mysterious, $$S(n,2)=\frac{6{H_n}^2-6\digamma'(n+1)+\pi^2}{12n}$$ With $\digamma'$ being the first derivative of the digamma function. For higher $m$ it only continues to get more complicated. For $m=3$, the even crazier $$S(n,3)=\frac{2{H_n}^3+(\pi^2-6\digamma'(n+1))H_n+2(\digamma''(n+1)+2\zeta(3))}{12n}$$ With $\zeta$ being the Riemann zeta function. And for $m=4$ and onwards, Wolfram fails to find any representations.
Can anyone make any sense of this sum for higher $m$? Or perhaps even provide some justification for the lower order identities $$\sum_{k=0}^{n-1}\frac{(-1)^{n-1-k}{}_{(n-1)}\mathrm{C}_k}{n-k}=\frac{1}{n}$$ And $$\sum_{k=0}^{n-1}\frac{(-1)^{n-1-k}{}_{(n-1)}\mathrm{C}_k}{(n-k)^2}=\frac{H_n}{n}$$
