The binomial distribution $\mathrm{Binom}(n,p)$, which counts the total number of successes within n independent trials each succeeding with probability p, is of historical and fundamental importance for probability theory and applied statistical inference. In particular, it appears in quantitative variants of the central limit theorem [16, 8, 25, 29], and is broadly used in statistical modelling [17, 1, 9, 7, 26] also as a building block of more advanced models [27]; this includes performing A/B tests on conversion rates popular in e-commerce [11].

Despite this large body of work on approximate inference, little is known about the exact higher moments of the binomial distribution. Besides being of natural interest, such formulas are needed to develop various concentration bounds; for example, binomial moment bounds were proven useful for analyzing random projections [23].

While the textbooks usually cover only the variance, sometimes also the skewness and kurtosis), there have been only few research papers discussing formulas for binomial moments of order $d\gt 4$. The first recursion formula for binomial moments appeared in [3] for the special case of $p=\frac{1}{2}$. The case of general p was handled in [13] by means of recursions utilizing Stirling numbers of the first kind. This was subsequently simplified by means of moment generating functions in [15], and resulted in a more compact formula involving Stirling numbers of the second kind. Very recently, a recursion-free derivation of raw moments has been presented in [19]. Overall, the common idea is to see the moments as (more or less explicit) polynomials in n and p and group terms, in order to simplify the formula or establish other desired properties such as nonnegativity of terms or boundedness of recursion depth (as in the recent result on Beta moments [24]).

The discussed approaches still do not offer a satisfactory answer, as the formulas are not handy enough to be directly applicable. The author of the most general formula in [15] didn’t manage to obtain nonnaive bounds on the binomial moments: the bound $O{(nq)^{d}}$ with $q=1-p$ on the d-th central moment [15] valid for $p\lt \frac{1}{2}$ is trivial as the centered binomial random variable is bounded between $-pn$ and $qn$ – no extra formulas are needed; this bound is far from the true behaviour $O{(nq)^{d/2}}$ when $nq\to \infty $ (obtained by the Central Limit Theorem). The main formula in [15] is actually a mixture of positive and negative contributions, which makes its numerical convergence problematic and theoretical analysis very difficult (as seen above). Moreover, all the prior works do not fully exploit the symmetry and produce rather complicated formulas in terms of p; one could expect simpler expressions on the central moments with the variance $n{\sigma ^{2}}=np(1-p)$ as a variable. Lastly, the discussed prior works are rather scarce in their presentation of related works and techniques, in particular they seem to have overlooked that the formulas with the appearance of Stirling numbers follow easier by the established approach of factorial moments [14, 2].

Addressing the aforementioned issues with approaches in prior works, this paper offers the following novel contributions to computing binomial moments: In summary, when compared to prior works, these results bring a broader scope of the techniques, as well as lead to arguably more handy formula; another added value is the contributed algorithm and its Python implementation.1

The remainder of the paper is organized as follows: the necessary background is given in Section 2, proofs appear in Section 4, the results are presented in Section 3, the sharp asymptotics are discussed in Section 5, and Section 6 concludes the work. The Python implementation is presented in Appendix A.

A random variable S follows the binomial distribution with parameters n and p, denoted as $S\sim \mathrm{Binom}(n,p)$, when the probability density function is

Let d be a positive integer. The raw moment of order d of a random variable S is defined as $\mathbb{E}[{S^{d}}]$, while the central moment of order d of S equals $\mathbb{E}[{(S-\mathbb{E}[S])^{d}}]$. We also use the factorial moment, defined as $\mathbb{E}[{S^{\underline{d}}}]$ where ${x^{\underline{d}}}=x(x-1)\cdots (x-(d-1))$ is called the d-th falling factorial [12].

To state the results, we need Stirling numbers of the second kind. The symbol $\left\{\genfrac{}{}{0.0pt}{}{n}{k}\right\}$ stands for the number of ways of partitioning an n element set into k nonempty subsets. We also need multinomial coefficients defined as $\left(\genfrac{}{}{0.0pt}{}{d}{{d_{1}}\dots {d_{k}}}\right)=\frac{d!}{{d_{1}}!\cdots {d_{k}}!}$ when ${\textstyle\sum _{i=1}^{k}}{d_{i}}=d$ and ${\min _{i}}{d_{i}}\geqslant 0$ and 0 otherwise, which extend the binomial coefficients. By the multinomial theorem, we have that ${({x_{1}}+\cdots +{x_{n}})^{d}}={\textstyle\sum _{{d_{1}},\dots ,{d_{n}}}}\left(\genfrac{}{}{0.0pt}{}{d}{{d_{1}}\dots {d_{n}}}\right){x^{{d_{1}}}}\cdots {x_{n}^{{d_{n}}}}$, extending the binomial formula.

To work out the desired polynomial formulas, we need some standard algebraic notation. By $\mathbb{Z}[{x_{1}},\dots ,{x_{k}}]$ we denote polynomials with integer coefficients in variables ${x_{1}},\dots ,{x_{k}}$. A polynomial is symmetric if it is invariant under interchanging variables, and antisymmetric when it gets negated under exchange of any two variables. The fundamental theorem of symmetric polynomials states that any symmetric polynomial from $\mathbb{Z}[{x_{1}},\dots ,{x_{k}}]$ can be written as a polynomial in the elementary symmetric functions ${s_{j}}({x_{1}},\dots ,{x_{n}})={\textstyle\sum _{1\leqslant {i_{1}}\dots \leqslant {i_{j}}}}{x_{{i_{1}}}}\cdots {x_{{i_{j}}}}$ for $j=1,\dots ,k$, with integer coefficients. Moreover, antisymmetric polynomials can be written as a product of a symmetric polynomial and Vandermonde’s determinant ${\textstyle\prod _{1\leqslant i\lt j\leqslant k}}({x_{i}}-{x_{j}})$ (see, for example, [22, 28]).

The first result is the derivation of the closed-form formula for raw binomial moments. While this formula appears in prior works [3, 15, 13], the arguments developed here are simpler, as they avoid the machinery of generating functions. The two alternative proof techniques are, respectively: a) linking central and factorial moments b) a direct counting argument.

The two proofs appear respectively in Section 4.1 and Section 4.2. Below, in Table 1, we list the explicit expressions for the first 10 moments.

Symmetric structure. While in prior works the formulas are derived in terms of p, we go beyond that exploiting the symmetry and showing that the formulas can be written in terms of the variance ${\sigma ^{2}}=p(1-p)$, which makes them much simpler. The following theorem proves what can be conjectured by inspection of known formulas for small-order moments.

Table 2 illustrates this result, providing explicit moments of order $d=2,\dots ,10$. The usefulness of the formula in Theorem 2 is its simplicity when compared to the representation in terms of p alone. The result is intuitive, but not straightforward to prove; two arguments are given, based on a) the theory of symmetric functions, see Section 4.3 and b) the novel combinatorial formula, see Section 4.4. The algorithm deriving the exact formulas is discussed later.

Positive polynomial representation. As mentioned in the introduction, the only closed-form formula due to [15] is an alternating sum with no readable leading term, which makes it hard to use; in particular the discussion in [15] fails to give nontrivial bounds on binomial moments. The novelty of this work is a formula consisting of positive terms. This makes it more stable for numerical computations and more handy in theoretical analysis.

To illustrate how useful is our positive representation established in Theorem 3, we derive the sharp bounds on (normalized) central binomial moments. This problem has remained open so far; some recent works used ad hoc upper bounds to estimate the binomial moments and tails (works on random projections, particularly [23]).

The estimate is uniformly sharp in all parameters; for the special case when $n{\sigma ^{2}}\to \infty $ and d is fixed, the d-th central moment grows, up to a constant, as ${(n{\sigma ^{2}})^{d/2}}$ which matches the central limit theorem combined with the explicit formulas for moments of the normal distribution [20]. In contrast to Theorem 4, the formula in [15] gives in this setup only much worse $O{(nq)^{d}}$, which anyway follows trivially since S is bounded by n. The proof uses Theorem 3.

We have seen that the variance-based representation in Theorem 2 is useful, yet it is not immediate how to efficiently compute this polynomial. To this end, we develop two different algorithms, both implemented in the popular Python symbolic algebra package Sympy [18].

Algorithm 1 takes advantage of the fundamental theorem on symmetric polynomials (see, for example, [10]). Specifically, there is an explicit procedure for converting any symmetric polynomial in $p,q$ into a polynomial in variables $p+q,pq$ (the basic symmetric polynomials in two variables); substituting $p+q=1$ we are left with a polynomial in ${\sigma ^{2}}=pq$. For even d we start with a symmetric polynomial and obtain a polynomial in ${\sigma ^{2}}$. In turn, for odd d we start with an antisymmetric polynomial and apply this procedure to its symmetric part, which results in a polynomial in ${\sigma ^{2}}$ plus the factor $q-p=1-2p$.

In turn, Algorithm 2 uses the power of elimination theory, to recover the representation in ${\sigma ^{2}}$ from any formula written in terms of p. Essentially, it simplifies the input polynomial in p with respect to the polynomial ${\sigma ^{2}}-p(1-p)$ leveraging the elimination properties of Groebner bases (see, for example, [5]). The output is a polynomial in ${\sigma ^{2}}$ (plus the factor $1-2p$ for odd d).

The proof is based on the fact that the factorial moments of the binomial distribution are easy to compute. Namely (see [21]), we have the following proposition.

Then it remains to connect factorial moments to standard moments, or in other terms: factorial powers to powers. It is well known (see, for example, the discussion in [4]) that this base change is given in terms of the Stirling numbers of the second kind. We state this fact formally below.

Now Theorem 1 follows by applying Proposition 2 to $x:=S$, and then using Proposition 1 to compute the expectation of ${S^{\underline{k}}}$ for $k=0,\dots ,d$.

Here we take a direct approach, writing $S={\textstyle\sum _{i=1}^{n}}{X_{i}}$ where ${X_{i}}{\sim ^{iid}}\mathrm{Bern}(p)$. Using the multinomial expansion and the independence of ${X_{i}}$ we obtain We now group the expressions in the above sum, depending on the number of nonzero elements in ${({d_{i}})_{i}}$. Denoting $\| ({d_{i}}){\| _{0}}=\mathrm{\# }\{i:{d_{i}}\gt 0\}$ and using the fact that $\mathbb{E}{X_{i}^{{d_{i}}}}=p$ when ${d_{i}}\gt 0$, we obtain By the symmetry of multinomial coefficients this equals Finally, we observe that the expression ${\textstyle\sum _{{d_{1}}\dots {d_{k}}\gt 0}}\left(\genfrac{}{}{0.0pt}{}{d}{{d_{1}}\dots {d_{k}}}\right)$ counts the number of ways of partitioning $\{1,\dots ,d\}$ into k nonempty labeled subsets; thus, this numbers equals $k!\cdot \left\{\genfrac{}{}{0.0pt}{}{d}{k}\right\}$ which finishes the proof.

From Equation (1) we obtain Replacing k by $n-k$ and using the symmetry of binomial coefficients we obtain When d is even, comparing these two equivalent expressions we see that they are symmetric as polynomials in p and q. By the fundamental theorem of symmetric polynomials, this can be written as a polynomial in $pq$ and $p+q$; in our case $p+q=1$ and the claim follows. If d is odd then ${(nq-k)^{d}}=-{(k-nq)^{d}}$ and we get antisymmetric polynomials in $p,q$ which can be written as a product of $p-q$ and a symmetric polynomial. The latter, by the fundamental theorem, is a polynomial in $p+q$ and $pq$; since $p+q=1$ the result follows.

By inspecting the products ${\textstyle\prod _{i=1}^{k}}({q^{{d_{i}}-1}}-{(-p)^{{d_{i}}-1}})$ that appear in Theorem 3, it can be seen that each of them is symmetric in $p,q$ when ${\textstyle\sum _{i}}{d_{i}}=d$ is even, and antisymmetric when ${\textstyle\sum _{i}}{d_{i}}=d$ is odd. This is because $p,q\to {q^{{d_{i}}-1}}-{(-p)^{{d_{i}}-1}}$ is symmetric when ${d_{i}}$ is even, and antisymmetric otherwise. The claim now follows.

As in the proof of Theorem 1 we arrive at Denote $x=1-\frac{1}{1-p}$, then $\mathbb{E}{[({X_{i}}-\mathbb{E}[{X_{i}}])]^{{d_{i}}}}=p{(1-p)^{{d_{i}}}}(1-{x^{{d_{i}}-1}})$ and thus With some further simplifications and grouping we can write or equivalently This finishes the proof. In addition to that, in what follows, we discuss how to further group terms and speed up computations. We can write and since $\left(\genfrac{}{}{0.0pt}{}{d}{{d_{1}}\dots {d_{k}}}\right)=\left(\genfrac{}{}{0.0pt}{}{{d_{1}}+\cdots +{d_{j}}}{{d_{1}}\dots {d_{j}}}\right)\cdot \left(\genfrac{}{}{0.0pt}{}{d-({d_{1}}+\cdots +{d_{j}})}{{d_{j+1}}\dots {d_{k}}}\right)\cdot \left(\genfrac{}{}{0.0pt}{}{d}{{d_{1}}+\cdots +{d_{j}}}\right)$ we obtain and thus where ${S_{2}}(n,k)$ denotes the number of ways of partitioning an n-element set into k subsets of cardinality at least 2 (a variation on Stirling numbers of the second kind). This can be used to develop an equivalent, but faster to compute, formula.

Throughout the proof, we will use the elementary estimates ${(m/\mathrm{e})^{m}}\leqslant m!\leqslant {m^{m}}$ (the second inequality is obvious, the first one follows by rearranging to ${\mathrm{e}^{m}}\geqslant \frac{{m^{m}}}{m!}$ and Taylor’s expansion) and ${(n/k)^{k}}\leqslant \left(\genfrac{}{}{0.0pt}{}{n}{k}\right)\leqslant {(n\mathrm{e}/k)^{k}}$ (see, for example, [6]).

By applying Remark 1 we can assume $p\leqslant \frac{1}{2}$. Then $q\geqslant p$, and thus $0\leqslant {q^{{d_{i}}-1}}-{(-p)^{{d_{i}}-1}}\leqslant p+q=1$ for any ${d_{i}}\geqslant 2$. In view of Theorem 3, we obtain the bound Since we have we further obtain Denoting ${\sigma ^{2}}=pq$ and using the bound $\left(\genfrac{}{}{0.0pt}{}{n}{k}\right)\leqslant {(n\mathrm{e}/k)^{k}}$, we finally obtain so that (noticing that ${(d/2)^{1/d}}\leqslant {(5/2)^{1/5}}$ for $d\geqslant 2$)

We now move on to the lower bound. The idea is to fix k and look at terms such that ${d_{1}},\dots ,{d_{k}}$ are nearly equal, to make $\left(\genfrac{}{}{0.0pt}{}{d}{{d_{1}}\dots {d_{k}}}\right)$ possibly large.

Let us write $d=k\cdot r+\ell $ with an integer $r\geqslant 2$ and nonnegative integer ℓ such that $0\leqslant \ell \lt k$; this is possible as $k\leqslant d/2$, by writing $d=k\cdot r+\ell $ where $0\leqslant \ell \lt k$. Define ${d_{i}}=r+1$ when $1\leqslant i\leqslant \ell $ and ${d_{i}}=r$ when $\ell \lt i\leqslant k$, so that ${\textstyle\sum _{i}}{d_{i}}=d$. For this tuple $({d_{i}})$ we have Since $d!\geqslant {(d/\mathrm{e})^{d}}$, $r+1\leqslant \frac{3}{2}r$, and $d/r\geqslant k$, it follows that

When ${d_{i}}$ are even, we have ${q^{{d_{i}}-1}}-{(-p)^{{d_{i}}-1}}={q^{{d_{i}}-1}}+{p^{{d_{i}}-1}}\geqslant 2\cdot \frac{1}{{2^{{d_{i}}-1}}}$ by Jensen’s inequality applied to the function $u\to {u^{{d_{i}}-1}}$ and $p+q=1$. Since in the summation we consider ${({d_{i}})_{i}}$ such that ${\textstyle\sum _{i}}{d_{i}}=d$ and ${d_{i}}\geqslant 2$, we obtain

The above two bounds, in view of Theorem 3, imply that Denoting ${\sigma ^{2}}=pq$ and using the bound $\left(\genfrac{}{}{0.0pt}{}{n}{k}\right)\geqslant {(n/k)^{k}}$, we finally obtain which finishes the proof.