FIVE Applications of Wilf-Zeilberger Theory to Enumeration and Probability ^{1}

Moa APAGODU and Doron ZEILBERGER

Dedicated to Sergei Abramov on becoming FIVE-dozen years old

Explicit Formulas vs. Algorithms

In the old days, when one had to find some sequence, a(n), there were two extremes. In the lucky case, one had an explicit formula. For example, the probability of tossing a fair coin 2n times and getting exactly n Heads, equals (2n)!/(2^{2n}n!^{2}). Sometimes, cheatingly, one considered as ‘explicit’ expressions in terms of sums (or multisums) or integrals (or multi-integrals). The other extreme was to just have a numerical algorithm, that for each (numeric!) input n, found the output. In that case the algorithm was rated by its efficiency.

Another compromise was an asymptotic formula, valid (approximately!) for large n.

But what’s a formula?, it is a kind of algorithm. Of course, it is more than that, theoretically, but from a practical point of view it should be judged by the efficiency of the implied algorithm.

The Holonomic Ansatz

Let’s look at the explicit formulas that are called ‘closed-form’, or more precisely hypergeometric sequences. A sequence a(n) is called hypergeometric if the ratio a(n+1)/a(n) is a rational function of n, i.e. a quotient P (n)/Q(n) where P (n) and Q(n) are polynomials. For example for the above-mentioned probability of getting exactly n Heads when tossing a fair coin 2n times, p(n) := (2n)!/(2^{2n}n!^{2}), we have p(n + 1)/p(n) = (2n + 1)/(2(n + 1)), or, by cross-multiplying

2(n + 1)p(n + 1) − (2n + 1)p(n) = 0 .

This is an example of a first-order linear recurrence equation with polynomial coefficients. Once you have the trivial value p(0) = 1 you can use it to compile a table of p(n) for n < N, for any desired N in O(N) operations.

The same is true for solutions of any linear recurrence equation with polynomial coefficients,

∑_{i=0}^{L} a_{i}(n)p(n + i) = 0 ,

of order L. The only difference is that we need L initial conditions, p(0), p(1), …, p(L − 1). We also assume that a_{L}(n) = 0 has no positive integer roots.

1 Oct. 18, 2006. Accompanied by Maple packages AppsWZ and AppsWZmulti downloadable from Zeilberger’s website.

Sample input and output can be viewed in:

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/ appswz.html . The research of the second author was supported in part by the NSF.

Such sequences were dubbed P-recursive by Richard Stanley in his seminal paper [St], but we prefer the name holonomic. Zeilberger[Z] famously showed that many sequences that arise in combinatorics, probability, and elsewhere are holonomic, and this was made into a full-fledged algorithmic theory by Wilf and Zeilberger [WilZ].

In our humble opinion, a holonomic representation of a sequence is to be considered explicit, since it is almost as good as a closed-form (i.e. hypergeometric). In the previous literature on WZ theory, there were few scattered examples of potential applications, but the focus was on the theory and the algorithms, not on specific applications.

Why this Paper?

The purpose of this paper is to fill this gap. We only list five such applications, but the reader can doubtless find many others. It is hoped that our implementation of these five applications will aid the reader to implement other ones that he or she might be interested in.

The Maple packages AppsWZ and AppsWZmulti

This article is accompanied by two Maple packages. AppsWZ that does applications of the singlevariable case ([AlZ]), and AppsWZmulti that does applications of the multi-variable case ([WilZ],[ApZ]).

Asymptotics

Another nice feature of being a solution of a linear recurrence equation with polynomial coefficients is that using the Birkhoff-Trijinski method (see [WimZ] for a lucid exposition), one can deduce the asymptotics to any desired order. This algorithm has been implemented by us in Maple, and is part of both packages.

First Application: Rolling a Die

If instead of tossing a coin n times, you roll a k-faced die, marked with positive or negative numbers, and you win the amount shown on the landed face (or lose, if it is a negative number). What is the probability that after n rolls, you break even? More generally, how likely are you to win exactly d dollars?

If the i^{th} face (i = 1 … k) shows the amount m_{i}, and lands with probability p_{i}, let the probability generating function of the die be defined by

P (x) := ∑_{i=1}^{k} p_{i}x^{m_{i}} .

It is very well known (and very very easy to see), that the probability, let’s call it a_{d}(n), of winding up with d dollars after n rolls is the constant term of

\frac{P (x)^{n}}{x^{d}} .

But that’s exactly grist for the Almkvist-Zeilberger mill! In the package AppsWZ this is accom-

panied by the commands RecProbVisit and RecProbVisitE.

For example, let a(n) be the probability of rolling a fair (standard) die 2n times and having the total score being exactly the expected value 7n. Then a(n) satisfies the following third-order linear recurrence equation with polynomial coefficients

−4(5 + 2n)(2n + 3)(2n + 1)(7n + 19)(5n + 11)(7n + 20)(7n + 13)(n + 2)(n + 1)a(n)

+4(7n+20)(5+2n)(2n+3)(n+2)(25480n^{5}+223496n^{4}+755066n^{3}+1223233n^{2}+946889n+279936)a(n+1)−

(5n + 6)(5 + 2n)(6 + 7n)(499359n^{6} + 6777015n^{5} + 38079431n^{4}
+113390385n^{3} + 188723986n^{2} + 166469280n + 60800544)a(n + 2)+

30(5n + 14)(5n + 13)(5n + 12)(7n + 12)(5n + 11)(5n + 6)(7n + 13)(6 + 7n)(n + 3)a(n + 3) = 0 ,

and the Birkhoff-Trijinksi method implies that the asymptotics is:

(.197833497170804)n^{−1/2}(1−(111/1400)/n−(12037/5488000)/n^{2}+(1367631/1097600000)/n^{3}+…) .

Readers can produce their own output for the scenario of their choice.

Second Application: How many ways to have r people chip in to pay a bill of n cents

In George Pólya’s classic ‘picture writing’ paper[Po], he considers the problem of figuring out how many ways can one person pay a bill of n cents using any number of coins. If the denominations are {d_{1}, …, d_{k}} ({1, 5, 10, 25, 50, 100} in the US), then the required number is the coefficient of x^{n} in the generating function

∏_{i=1}^{k} \frac{1}{1 − x^{d_{i}}} .

Calling this number a(n), this entails, trivially, a linear recurrence equation with constant coefficients of order ∑_{i=1}^{k} d_{i}. But by allowing polynomial coefficients, one can get, thanks to Almkvist-Zeilberger, a recurrence of order ≤ lcm(d_{1}, …, d_{k}).

More generally, of a_{r}(n) is the number of ways of breaking n cents with (up to) r people participating (or equivalently, one person with r pockets in his or her pants or dress), the generating function is

∏_{i=1}^{k} \frac{1}{(1 − x^{d_{i}})^{r}} ,

and applying Almkvist-Zeilberger to

\frac{1}{x^{n+1}} ∏_{i=1}^{k} \frac{1}{(1 − x^{d_{i}})^{r}} ,

produces a linear recurrence, still with polynomial coefficients (but now these are polynomials in both n and r) of the above order. On the other hand if you stick to constant coefficients then the order would be r(∑_{i=1}^{k} d_{i}), and would only make sense for (small!) numeric r, while within the holonomic ansatz, one can have symbolic r without any increase in the order. The procedure that takes care of this problem in the Maple package AppsWZ is TeamEffortMoneyChanging. For example, entering:

TeamEffortMoneyChanging({1,5,10,25},n,N,r);

would tell you that if a_{r}(n) is the number of ways of breaking n cents with (up to) r people chipping-in, using only pennies, nickels, dimes, and quarters, and denoting by N the shift operator in n: (Nx(n) := x(n + 1)), it turns out that a_{r}(n) is annihilated by the following 30^{th}-order linear recurrence operator with polynomial coefficients:

(−41r − n) + (−r)N + (−r)N^{2} + (−r)N^{3} + (−r)N^{4}+
(−47r − n − 5)N^{5} + (−2r)N^{6} + (−2r)N^{7} + (−2r)N^{8} + (−2r)N^{9} + (−22r)N^{10}
+(−2r)N^{11} + (−2r)N^{12} + (−2r)N^{13} + (−2r)N^{14}
+(−22r)N^{15} + (−2r)N^{16} + (−2r)N^{17} + (−2r)N^{18} + (−2r)N^{19} + (−22r)N^{20}+
(−2r)N^{21} + (−2r)N^{22} + (−2r)N^{23} + (−2r)N^{24} + (−6r + n + 25)N^{25}
+(−r)N^{26} + (−r)N^{27} + (−r)N^{28} + (−r)N^{29} + (n + 30)N^{30} .

Third Application: Hidden Markov Models

In contemporary bioinformatics (see for example the fascinating and lucid expository article by Lior Pachter and Bernd Sturmfels in the wonderful book that they edited([PaStu]), one has to estimate probabilities from output. The usual approach is maximum likelihood, that entails solving huge systems of polynomial equations that are handled via Buchberger’s Gröbner Basis Algorithm. But, following Laplace, we can also use a Bayesian approach and get quotients of integrals. These integrals can be handled via the multi-Almkvist-Zeilberger ([ApZ]) algorithm. For any specific instance one may want to use numerical integration. If one wants to compile a table anticipating all possible outputs (of course up to a certain length), then the recurrences supplied by Almkvist-Zeilberger are much more efficient.

Let’s illustrate it with the simple case of (only) two dice, that enables us to stay within the singlevariable Almkvist-Zeilberger ([AlZ]).

In the proverbial dishonest casino, suppose that there are two kinds of (identically looking) dice, with the same possible outcomes (faces), let’s call them {1, …, m}, each with a known probability distribution (p_{1}, …, p_{m}) and (q_{1}, …, q_{m}). What is unknown is the probability, x, of using the first die (and hence 1 − x for the other die). Judging from the sequence of outcomes (or rather by the relative frequencies of the landed faces), one has to estimate x. If the (hidden) probability of using the first die was x, (whatever it is), then the probability of the output distribution being a_{1}, …, a_{m} is (a_{1} + … + a_{m})!/(a_{1}! ⋯ a_{m}!) times

L(x) := ∏_{i=1}^{m}(xp_{i} + (1 − x)q_{i})^{a_{i}} .

The maximum-likelihood estimate is to maximize L(x) by solving L′(x) = 0 (with more dice one gets partial derivatives and systems of equations in several variables). But, following the more democratic approach of Laplace, that considers all scenarios (and that famously tells you that if so far you succeeded m times and failed n times, then your estimated probability of success in your next try is not m/(m+n) but rather (m+1)/(m+n+2)), we would have not the root of L′(x) = 0 but rather

\frac{∫_{0}^{1} x ∏_{i=1}^{m}(xp_{i} + (1 − x)q_{i})^{a_{i}} dx}{∫_{0}^{1} ∏_{i=1}^{m}(xp_{i} + (1 − x)q_{i})^{a_{i}} dx} .

Calling the top quantity T (a_{1}, …, a_{m}) and the bottom quantity B(a_{1}, …, a_{m}), we see that Almkvist-Zeilberger can manufacture pure linear recurrences (with polynomial coefficients) in each of the discrete variables a_{1}, …, a_{m} (in Nature m = 4, or m = 20). The multi-dice analog of this (which would employ the multi-variable Almkvist-Zeilberger , done in [ApZ]) may be of interest to bioinformaticians. So far we have only implemented the interface in the two-dice case.

The procedure that takes care of this problem in AppsWZ is ExpProbFirstDie(Lp1,Lp2,a0), where Lp1, and LP2 are the probability distributions and a0 is the list of respective outcome. The novelty here is that both top and bottom integrals are not computed directly but via the recurrences outputted by Almkvist-Zeilberger. With the option remember it should be much more efficient if one wants to pre-compute a table of estimated probability for each a0 of size less than some (lagre) pre-assigned value.

For example, if there are two coins, one fair and one loaded with probability of a Head being 1/3, and the outcome was 10 and 10 (which ML would say that the fair coin was used all the time),

evalf(ExpProbFirstDie([1/2,1/2],[1/3,2/3],[10,10])):

would estimate that the fair coin was used only %58.76 of the time.

Fourth Application: Lattice Paths Counting

We all know that the number of ways of walking from the origin (0, 0) to the point (m, n) in the square lattice, with unit northbound and unit eastbound steps {(0, 1), (1, 0)} is

F (m, n) = \frac{(m + n)!}{m!n!} .

This immediately implies (and is equivalent to) the fact that

\frac{F (m + 1, n)}{F (m, n)} = \frac{m + n + 1}{m + 1} , \frac{F (m, n + 1)}{F (m, n)} = \frac{m + n + 1}{n + 1} .

and cross-multiplying yields

(m + 1)F (m + 1, n) − (m + n + 1)F (m, n) = 0 , (n + 1)F (m, n + 1) − (m + n + 1)F (m, n) = 0 .

In other words the discrete function F (m, n) satisfies pure linear recurrences equations with polynomial coefficients, that happen, in this simple case, to be first-order. Recall that a recurrence is pure if only one of the variables changes at a time. For example, F (m, n) trivially satisfies the “mixed” (partial) recurrence F (m, n) = F (m − 1, n) + F (m, n − 1).

An amazing consequence of Wilf-Zeilberger theory ([WilZ]) is that this is still true for an arbitrary set of (positive) steps, and in arbitrary dimension. Of course the pure recurrences are no longer (usually) first-order, but as above, this is a minor computational disadvantage.

Indeed, if we are walking in the d-dimensional (hyper)cubic lattice, starting at the origin, and with a set of steps S (all with non-negative coordinates, excluding the step 0 [staying in place]), the generating function is trivially seen to be

∑_{m} F (m_{1}, …, m_{d})x_{1}^{m_{1}} ⋯ x_{d}^{m_{d}}

= \frac{1}{(1 − ∑_{(s_{1},...,s_{d})∈S} x_{1}^{s_{1}} ⋯ x_{d}^{s_{d}})} .

So it follows that our discrete function of interest, F (m_{1}, …, m_{d}) equals the formal residue of

\frac{x_{1}^{−m_{1}−1} ⋯ x_{d}^{−m_{d}−1}}{(1 − ∑_{(s_{1},...,s_{d})∈S} x_{1}^{s_{1}} ⋯ x_{d}^{s_{d}})} .

If you are only interested in getting to points on the diagonal, then f(n) := F (n, n, …, n) is given by the formal residue of

\frac{x_{1}^{−n−1} ⋯ x_{d}^{−n−1}}{(1 − ∑_{(s_{1},...,s_{d})∈S} x_{1}^{s_{1}} ⋯ x_{d}^{s_{d}})} ,

and once again thanks to Wilf-Zeilberger theory, satisfies a linear recurrence equation with polynomial coefficients.

This is implemented in the second Maple package accompanying this paper, AppsWZmulti that is “powered” by the Maple package MultiAlmkvistZeilberger that accompanied [ApZ]. The relevant

procedures (in AppsWZmulti) are LatticePaths and LatticePathsDiagonal and for the verbose
versions LatticePathsStory and LatticePathsDiagonalStory.

For example, if you type LatticePaths({[0,1],[1,0],[1,1]},m,M); you would get that the following two operators annihilate F (m_{1}, m_{2}), the number of ways of going from (0, 0) to (m_{1}, m_{2}) using the steps {(0, 1), (1, 0), (1, 1)} (where M_{1}, M_{2} are the shift operators in the m_{1}, m_{2} variables respectively)

[−(m_{1}+1)/(2+m_{1})−((1+2m_{2})/(2+m_{1}))M_{1}+M_{1}^{2}, −(m_{2}+1)/(2+m_{2})−((1+2m_{1})/(2+m_{2}))M_{2}+M_{2}^{2}] , which in everyday parlance means that F (m_{1}, m_{2}) satisfies

(m_{1} + 2)F (m_{1} + 2, m_{2}) − (2m_{2} + 1)F (m_{1} + 1, m_{2}) − (m_{1} + 1)F (m_{1}, m_{2}) = 0 ,
(m_{2} + 2)F (m_{1}, m_{2} + 2) − (2m_{1} + 1)F (m_{1}, m_{2} + 1) − (m_{2} + 1)F (m_{1}, m_{2}) = 0 .

If you type LatticePathsDiagonal({[0,1],[1,0],[1,1]},n,N); you would get that f(n), the

number of ways of getting from (0, 0) to (n, n) using the same set of steps is:

(n + 2)f(n + 2) − 3(2n + 3)f(n + 1) + (n + 1)f(n) = 0 ,

subject to the initial conditions f(0) = 1,f(1) = 3. Thanks to Birkhoff-Trijinski, its asymptotics is

C(3+2^{3/2})^{n}n^{−1/2}(1+((3/32)2^{1/2}−1/4)/n+(113/1024−(9/128)2^{1/2})/n^{2}+((1545/32768)2^{1/2}−245/4096)/n^{3}) ,

for some constant C.

A more interesting example are the (old-time) basketball numbers, which is the number of ways a basketball game that ended with the score n : n can proceed. Recall that in the old days (before 1961), an atom of basketball-scoring could be only of one or two points. Equivalently, this number is the number of ways of walking, in the square lattice, from (0, 0) to (n, n) using the atomic steps {(1, 0), (2, 0), (0, 1), (0, 2)}. Entering this into LatticePathsDiagonalStory yields that, calling this number F (n), it satisfies the third-order linear recurrence:

(16/5)(2n + 3)(11n + 26)(1 + n)/((n + 3)(2 + n)(11n + 15))F (n)
−(4/5)(121n^{3} + 649n^{2} + 1135n + 646)/((n + 3)(2 + n)(11n + 15))F (1 + n)
−(2/5)(176n^{2} + 680n + 605)/((11n + 15)(n + 3))F (2 + n) + F (n + 3) = 0 ,

subject to the initial conditions:

F (0) = 1, F (1) = 2, F (2) = 14 .

For the record, the first few terms are:

[1, 2, 14, 84, 556, 3736, 25612, 177688, 1244398, 8777612, 62271384, 443847648, 3175924636, 22799963576, 164142004184, 1184574592592, 8567000931404, 62073936511496, 450518481039956, 3274628801768744, 23833760489660324].

The asymptotics is:

(.37305616)(4 + 2(3^{1/2}))^{n}n^{−1/2}(1 + (67/1452)3^{1/2} − (119/484))/n

+((6253/117128)−(7163/234256)3^{1/2})/n^{2} +(−(32645/15460896)3^{1/2} +(129625/10307264))/n^{3}) ,

or in floating-point:

(.37305616)(7.464101616)^{n}n^{−1/2}·
(1. − .1659453140/n + .42398086 · 10^{−3}/n^{2} + .8918933381 · 10^{−2}/n^{3}) .

Fifth Application: Random Walk in Higher Dimensions This is the multivariate analog of the First Application. The relevant procedures in AppsWZmulti

are RandomWalkRecurrence and RandomWalkRecurrenceE. By now readers should be able to gen-

erate their own examples. A few sample input and output files are given in the webpage of this article. References

[AlZ] G. Almkvist and D. Zeilberger, The method of differentiating under the integral sign, J. Symbolic Computation 10, 571-591 (1990).

[ApZ] M. Apagodu and D. Zeilberger Multi-Variable Zeilberger and Almkvist-Zeilberger Algorithms and the Sharpening of Wilf-Zeilberger Theory , Adv. Appl. Math. 37 (2006), 139-152.

[PaStu] L. Pachter and B. Sturmfels, Statistics, in: “Algebraic Statistics for Computational Biology” [L. Pachter and B. Sturmfels, eds.], Cambridge University Press, 2005, 3-33.

[Po] G. Pólya , On picture writing, Amer. Math. Monthly 63 (1956), 689-697.

[St] R. Stanley, Differentiably finite power series, Europ. J. Combinatorics 1 (1980), 175-188.

[WilZ] H.S. Wilf and D. Zeilberger, An algorithmic proof theory for hypergeometric (ordinary and “q”) multisum/integral identities, Invent. Math. 108(1992), 575-633 .

[WimZ] J. Wimp and D. Zeilberger, Resurrecting the asymptotics of linear recurrences, J. Math. Anal. Appl. 111(1985), 162-177 .

[Z] D. Zeilberger, A Holonomic systems approach to special functions identities, J. of Computational and Applied Math. 32 (1990), 321-368 . —– Moa Apagodu, Department of Mathematics, Virginia Commonwealth University, Richmond, VA 23284-2014, USA. mapagodu at vcu dot edu. Doron Zeilberger, Department of Mathematics, Rutgers University (New Brunswick), Hill Center-Busch Campus, 110 Frelinghuysen Rd., Piscataway, NJ 08854-8019, USA. zeilberg at math dot rutgers dot edu .