Automated Derivation of Limiting Distributions Of Combinatorial Random Variables Whose Generating Functions are Rational

By Doron ZEILBERGER

Abstract: The author’s methodology of automated derivation of explicit expressions for the average, variance, and higher moments of combinatorial random-variables that possess explicit generating functions is extended to those whose bi-variate generating function is an arbitrary rational function.

Maple Package and Sample Output

This article is accompanied by a Maple package, BiVariateMoms.txt, that is available, along with some sample input and output files, from

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/crv.html .

Ancient History

If you toss a fair coin n times, and are interested in the random variable ‘number of Heads’, let’s call it H_{n}, and plot the graph of

k → P rob(H_{n} = k) ,

and n is large enough, you very famously, get something close to a bell-shaped curve, and as n gets larger and larger, it would get closer and closer to it. The official names of the bell-shaped curve are ‘standard normal distribution’ and ‘Gaussian distribution’. More precisely, letting μ_{n} denote the expectation (aka mean, aka average), and σ_{n} the standard deviation (that happen to be equal to \frac{n}{2} and \frac{\sqrt{n}}{2} respectively), and defining the centralized and scaled version,

X_{n} := \frac{H_{n} − μ_{n}}{σ_{n}} ,

then the sequence of discrete random variables, X_{n}, ‘tend to’, in a certain precise sense, to the (continuous distribution) called standard normal distribution, whose probability density function is famously

\frac{1}{\sqrt{2π}} e^{− \frac{x^{2}}{2}} .

This means That

lim_{n→∞} P r(X_{n} ≤ x) = ∫_{−∞}^{x} \frac{1}{\sqrt{2π}} e^{−t^{2}/2} dt .

This was first proved, by de Moivre and Laplace, using fairly ad hoc methods, but has since been greatly generalized, and the collective name is Central Limit Theorems.

From Enumeration to Statistics Suppose that we have a finite set S, on which a certain numerical attribute, called random variable, X, (using the probability/statistics lingo), is defined. For any non-negative integer i, let’s define

N_{i} := ∑_{s∈S} X(s)^{i} .

In particular, N_{0}(X) is the number of elements of S. The expectation of X, E[X], denoted by μ, is, of course,

μ = \frac{N_{1}}{N_{0}} .

For i > 1, the i-th straight moment is

E[X^{i}] = \frac{N_{i}}{N_{0}} .

The i-th moment about the mean is

m_{i} := E[(X − μ)^{i}] = E[∑_{r=0}^{i} \binom{i}{r}(−1)^{r}μ^{r}X^{i−r}] = ∑_{r=0}^{i}(−1)^{r}\binom{i}{r}μ^{r}E[X^{i−r}]

= ∑_{r=0}^{i}(−1)^{r}\binom{i}{r} (\frac{N_{1}}{N_{0}} )^{r} \frac{N_{i−r}}{N_{0}}

= \frac{1}{N_{0}^{i}} ∑_{r=0}^{i}(−1)^{r}\binom{i}{r}N_{1}^{r}N_{0}^{i−r−1}N_{i−r} .

Finally, the most interesting quantity, statistically speaking, apart from the mean μ and variance m_{2} are the scaled-moments, also known as, alpha coefficients.

α_{i} := \frac{m_{i}}{m_{2}^{i/2}} .

Getting the Moments via Generating Functions

The generating function of our random variable, X, also called weight-enumerator, with respect to the variable t, is defined by

f(t) = f_{S,X}(t) := ∑_{s∈S} t^{X(s)} .

Note that when t = 1, we get the naive enumeration, |S|.

More generally we have, with the above notation,

N_{i} = (t\frac{t}{dt})^{i} f(t)|_{t=1} .

In many enumeration scenarios, we don’t have just one finite set, S, but an infinite sequence of sets S_{n}, and the generating function f_{n}(t) can be expressed explicitly in terms of both the ‘continuous’ variable, t, and discrete variable n.

For example, if the set S_{n} is the set {H, T }^{n} and X_{n}(s) is the “number of heads”, then by ‘independence’, we have a closed form formula for f_{n}(t):

f_{n}(t) = (1 + t)^{n} .

More generally, if S_{n} is the set of sequences of outcomes of rolling, n times, a k-faced die whose faces are marked with a_{1}, …, a_{k} dots, in other words the set {a_{1}, …, a_{k}}^{n}, and X_{n} is the random variable ‘total number of dots’ (i.e. the sum of the sequence of outcomes), then

f_{n}(t) = (∑_{j=1}^{k} t^{a_{j}} )^{n} .

So, for each non-negative integer i, we have a sequence of numbers, N_{i}(n), obtained by applying (t\frac{d}{dt})^{i} to f_{n}(t), and then plugging-in t = 1. We always get polynomials for the N_{i}(n). From them we can get polynomial expressions for the moments-about-the-mean, m_{i}(n). This is better delegated to computers. In many cases, one can also get the leading terms for symbolic i, and prove completely automatically, asymptotic normality up to any desired moment (i.e. that the scaled moments, α_{i}(n), tend, as n goes to infinity, to those of the standard normal distribution 1, 0, 3, 0, 5, 0, 15, 0, 105, . . .. In fact, it is even possible to prove it for an arbitrary moment (symbolic i), hence giving fully automated proofs of central limit theorems. This was described in [Z1] and [Z2].

The Grand Generating Function

Since we have an infinite sequence of polynomials, f_{n}(t), we can form the grand generating function

F (t, z) := ∑_{n=0}^{∞} f_{n}(t)z^{n} .

For example, for the coin-tossing example, we have

F (t, z) = \frac{1}{1 − (1 + t)z} ,

and for the more general die-rolling example, we have

F (t, z) = \frac{1}{1 − z (∑_{j=1}^{k} t^{a_{j}})} .

Notice that these are rational functions in both t and z.

But there are many examples in enumerative combinatorics (and statistical physics!), where there is no ‘closed-form’ expressions for the individual f_{n}(t), but the ‘grand-generating function’, F (t, z) is rational. These occur, for example, in tiling problems, see the beautiful article [Z3] with its accompanying Maple package TILINGS.

Note that in such cases, the enumerating sequence, {N_{0}(n)} is C-finite, i.e. satisfies a (homog.) linear recurrence equation with constant coefficients (see [Z4] for a brief overview of C-finite sequences, and the modern classic [KP] for an in-depth account.) Obviously,

∑_{n=0}^{∞} N_{0}(n)z^{n} = F (1, z) .

More generally,

∑_{n=0}^{∞} N_{i}(n)z^{n} = (t\frac{d}{dt})^{i} F (t, z)|_{t=1} ,

are always rational functions of z, and hence the sequences {N_{i}(n)}_{n=0}^{∞} are always C-finite.

Since F (t, z) is a rational function, we can write it as a quotient of polynomials

F (t, z) = \frac{P (t, z)}{Q(t, z)} .

It follows from the quotient rule from ‘calc1’, that, for each i,

∑_{n=0}^{∞} N_{i}(n)z^{n} = (t\frac{d}{dt})^{i} F (t, z)|_{t=1}= \frac{P_{i}(z)}{Q(1, z)^{i+1}} ,

for some polynomial P_{i}(z). Hence, by partial fraction decomposition, it follows that, denoting by L the order of the recurrence for N_{0}(n) (alias the degree, in z, of Q(1, z)), we are guaranteed that for each i, there exist polynomials A_{i,j}(n), 0 ≤ j ≤ L − 1, of degree i in n, such that

N_{i}(n) = ∑_{j=0}^{L−1} A_{i,j}(n) · N_{0}(n − j) .

Since we are guaranteed, a priori, for each i, that the L polynomials A_{i,0}(n), …, A_{i,L−1}(n) (of degree i in n), exist, we can ask our computers to find them empirically, by cranking-out sufficiently many terms of the sequence {N_{0}(n)} and {N_{i}(n)}, and using undetermined coefficients.

Once we have ‘explicit’ expressions (in the above sense) for N_{i}(n), we can get ‘explicit’ expressions (in the same sense), for the (straight) moments, and from them, the more informative moments about the mean, that, in turn, lead to the scaled moments α_{i}(n). Denoting by β the smallest positive root of Q(1, z), (a certain algebraic number) and noting that

lim_{n→∞} \frac{N_{0}(n − 1)}{N_{0}(n)} = β ,

and hence, more generally, for 0 ≤ j < L,

lim_{n→∞} \frac{N_{0}(n − j)}{N_{0}(n)} = β^{j} ,

we can get asymptotic expressions (with exponential ‘error’ in n) for all the quantities, in particular, for the expectation μ(n) and variance m_{2}(n) in terms of n and β. Also we can get asymptotic expressions (done automatically!) for the scaled moments α_{i}(n), and verify, each time, that as n goes to infinity, they converge to 0 if i is odd, and i!/(2^{i/2}(i/2)!) if i is even (i.e. the moments of the standard normal distribution). Even better, we can get more refined asymptotics, all done automatically, to any order in n.

It is possible to show, under mild hypothesis, that the sequence of random variable X_{n} whose grand-generating function is rational in t and z is always asymptotically normal, hence all we need is to find, in each case, the expectation, μ(n), and the variance m_{2}(n), in terms of N_{0}(n) and its shifts, and if all we care about is the limiting distribution, then it suffices to ask the computer to find the asymptotic expressions of these in terms of n and β.

The Maple Package BiVariateMoms.txt

Everything is implemented in the Maple package BiVariateMoms.txt. See the webpage of this article, where there are numerous input and output files.

The simplest non-trivial example

Let the set S_{n} be the set of sequences in {1, 2} that add-up to n, and let X_{n}(s) be the random variable ‘number of 2s’ in s. It is readily seen that

F (t, z) = \frac{1}{1 − z − t z^{2}} .

Here N_{0}(n) = F_{n+1}.Then the program says that the expectation, μ(n) equals exactly

\frac{2 n F_{n+1} − (n + 1) F_{n}}{5F_{n+1}} ,

while the variance, m_{2}(n), equals

m_{2}(n) = \frac{n(n + 3) F_{n+1}^{2} − (n^{2} − 1) F_{n+1}F_{n} − (n + 1)^{2} F_{n}^{2}}{25 F_{n+1}^{2}} ·

Let ϕ be the golden ratio, then N_{0}(n) = a(n) is asymptotic (with exponentially small error) to \frac{ϕ+1}{ϕ+2} · ϕ^{n}. The expectation is, asymptotically,

\frac{(2ϕ − 1)n − 1}{5ϕ} ,

and the variance, m_{2}(n), is, asymptotically,

\frac{(3ϕ + 1)n + (ϕ − 1)}{25(ϕ + 1)} .

For all the moments up to the sixth, see the output file

http://www.math.rutgers.edu/~zeilberg/tokhniot/oBiVariateMoms1.txt .

For numerous other examples, see the front of this article, mentioned above:

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/crv.html .

Acknowledgment: Many thanks to Tamar Zeilberger for insightful comments.

References

[KP] Manuel Kauers and Peter Paule, “The Concrete Tetrahedron”, Springer, 2011.

[Z1] Doron Zeilberger, The Automatic Central Limit Theorems Generator (and Much More!),

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/georgy.html ; Also appeared

in: “Advances in Combinatorial Mathematics: Proceedings of the Waterloo Workshop in Computer Algebra 2008 in honor of Georgy P. Egorychev”, chapter 8, pp. 165-174, (I.Kotsireas, E.Zima, eds. Springer Verlag, 2009.)

[Z2] Doron Zeilberger, HISTABRUT: A Maple Package for Symbol-Crunching in Probability theory, The Personal Journal of Shalosh B. Ekhad, Aug. 25, 2010,

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/histabrut.html . Also pub-
lished in https://arxiv.org/abs/1009.2984.

[Z3] Doron Zeilberger, Automatic CounTilings, The Personal Journal of Shalosh B. Ekhad and Doron Zeilberger, Jan. 20, 2006.

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/tilings.html .
[Z4] Doron Zeilberger, The C-finite Ansatz, http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/cfinite.html,

also published in the Ramanujan Journal 31 (2013), 23-32.

Doron Zeilberger, Department of Mathematics, Rutgers University (New Brunswick), Hill Center-Busch Campus, 110 Frelinghuysen Rd., Piscataway, NJ 08854-8019, USA;

DoronZeil at gmail dot com ; http://www.math.rutgers.edu/~zeilberg/ .

Exclusively published in The Personal Journal of Shalosh B. Ekhad and Doron Zeilberger

(http://www.math.rutgers.edu/~zeilberg/pj.html) .

Dec. 24, 2016