Automatic Proofs of Asymptotic ABNORMALITY (and much more!) of Natural Statistics Defined on Catalan-Counted Combinatorial Families

Shalosh B. EKHAD and Doron ZEILBERGER^{1}

Preliminary Sermon: Humans will be Humans; The Medium is the Message

The famous Catalan numbers (see [Sl1]), count zillions of combinatorial families (see [St]) and many humans have fun trying to find ‘nice’ bijections between family A and family B. While this may be fun for a while, sooner or later this game gets old, especially since the real reason Catalan numbers are so ubiquitous is their simplicity, and that humans can only grasp simple things.

Indeed, (see [Z]), the reason for the ubiquity of the sequence of Catalan numbers, {c_{n}}, is that their generating function

C(z) := ∑_{n=0}^{∞} c_{n}z^{n} ,

satisfies the simplest possible (genuinely!) algebraic equation, namely

C(z) = 1 + zC(z)^{2} ,

that is equivalent to the quadratic recurrence satisfied by the Catalan numbers themselves, namely:

c_{n} = ∑_{k=1}^{n} c_{k−1}c_{n−k} , c_{0} = 1 .

Often, the members of the combinatorial family in question posses natural statistics, for example for Dyck paths, the number of ‘inversions’ (D (not necessarily immediately) ahead of U), or for 132-avoiding permutations, the number of occurrences of some given pattern, then it may happen that two different statistics ‘amazingly’ have the same average! Wow!, Let’s find a bijection! See, e.g., the humanly-generated article [B] (by human Miklós Bóna), that appeared in the very prestigious (and very selective!) Electronic journal of Combinatorics, that does its best [alas, not always successfully] to only accept the best papers. [It often errs on both sides, rejecting truly seminal papers, and accepting quite a few trivial ones.]

Humans can, with some effort, find closed-form expressions for the average (aka expectation, aka first moment), of a given combinatorial statistic (aka random variable), and if they try really hard, may be able to find the variance (aka second moment [about the mean]), but beyond that they should enlist their much superior silicon brethrern, and develop algorithms for discovering (and proving!) closed-form expressions for as many as possible moments. In addition to its intrinsic interest, this activity would also indicate whether the combinatorial statistic in question seems to be asymptotically normal (if the standardized moments, starting with the third, converge, as n goes to ∞, to 0, 3, 0, 15, 0, 105, …, the famous moments of the normal distribution), or whether it is (rigorously-) provably not normal (if the expression for the skewness (alias standardized 3rd moment) does not tend to zero, we are done!) .

1 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 , http://www.math.rutgers.edu/~zeilberg/ . First version March 21, 2014. This version: March 29, 2014. Accompanied by the Maple packages http://www.math.rutgers.edu/~zeilberg/tokhniot/AlgFunEq , and http://www.math.rutgers.edu/~zeilberg/tokhniot/Cheyne .

Supported in part by the NSF. Exclusively published in the Personal Journal of Shalosh B. Ekhad and Doron Zeilberger, and arxiv.org .

In the present article, a collaboration between a human (DZ) and computer (SBE) we do just that! The human wrote a Maple package available, free of charge from:

http://www.math.rutgers.edu/~zeilberg/tokhniot/AlgFunEq ,

that once written, can handle zillions of possible statistics defined on Catalan-families, and surprisesurprise, find zillions of Bóna-style ‘surprises’, and of course, prove them all fully rigorously. More importantly, it can prove asymptotic abnormality, by deriving (and proving!) closed-form expressions for the skewness, as an expression in n, and having done that, (automatically!) take the limit as n → ∞, and realize that that limit is not zero. Since for any asymptotically normal sequence of random variables, that limit should be zero, this constitutes a fully rigorous proof of asymptotic abnormality. But why stop with the skewness? Our program also finds closed-form expressions for the kurtosis, aka standardized fourth moment (and proves that its limit, as n → ∞, is not 3), as well as expressions for higher moments.

It is true that, with great effort, very smart humans, like Svante Janson ([J]), can do it by entirely human means (and can even handle all moments, at least recursively), but they can only do the leading asymptotics! Not even Svante Janson can find, just by hand, e.g., an exact closed-form expression, in n, for the sixth moment of the random variable ‘number of occurrences of the pattern 213’ in the set of 132-avoiding permutations of length n.

But do we really care about the 6th moment of some stupid statistic defined on some stupid family of sets? Of course not! The Medium is the Message! This article is but a case-study in human-computer collaboration. The human teaching the computer how to solve every conceivable problem in a wide class of combinatorial problems, by the human designing algorithms, then implementing them (in our case in Maple), and then letting the computer execute them. Once we get better and better at this kind of collaboration, we would be ready for the big time! Stand by (in 100 years or less) for a computer-generated proof of RH and NP ≠ P.

Maple packages and Sample Input and Output Files

As usual in the ongoing collaboration between the authors of the present article, the most important part is not the article, but the Maple packages that come with it, that can be dowloaded, free of charge, from the front of this article

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

and the numerous input and output files also available there. In particular, the file

http://www.math.rutgers.edu/~zeilberg/tokhniot/oAlgFunEq1BonaRedux ,

reproduces, in 0.117 seconds, all the results (all rigorously proved!) of the first section of [B], that handled averages of the number of occurrences of patterns of length ≤ 3 (see below for details). Before describing our new algorithms, we challenge our readers (both humans and machines), with a neat conjecture, based on ample evidence outputted by our Maple package AlgFunEq. Lots and Lots of Bóna-Style Surprises and a Conjecture Let AV_{n}(132) be the set of 132-avoiding permutations of {1, …, n}, and for any permutation π and pattern p, let a_{p}(π) be the number of occurrences of the pattern p in the permutation π (in other words, if π has length n and p has length k, the number of k-tuples

1 ≤ i_{1} < i_{2} < … < i_{k} ≤ n ,

such that π_{i_{1}}π_{i_{2}} … π_{i_{k}} is ‘order-isomorphic’ to p). For each pattern p, define the sequence

A_{p}(n) := ∑_{π∈AV_{n}(132)} a_{p}(π) .

In [B], Bóna observed that trivially A_{231}(n) = A_{312}(n) (since 231 and 312 are inverses of each other, and the class of 132 permutations is closed under taking inverse), but that surprisingly (at least to him), both are also equal to A_{213}(n). Hence we have the following facts.

• For k = 1 there is 1 Bóna class (of course!)

• For k = 2 there are 2 Bóna classes (of course!)

• For k = 3 there are 3 Bóna classes Can you spot a pattern?(pun intended!). Hint: the term for k = 4 is not 4. The first-named author

(See http://www.math.rutgers.edu/~zeilberg/tokhniot/oAlgFunEq1 ,
http://www.math.rutgers.edu/~zeilberg/tokhniot/oAlgFunEq1a9 and
http://www.math.rutgers.edu/~zeilberg/tokhniot/oAlgFunEq1a10)

rigorously proved that the numbers of distinct Bóna classes (i.e. distinct sequences {A_{p}(n)} as p ranges over all 132-avoiding patterns of length k) for 1 ≤ k ≤ 10, starting with k = 1, are as follows:

1, 2, 3, 5, 7, 11, 15, 22, 30, 42, ….

Going to the indispensable OEIS, we immediately realized that these are the first ten values of [Sl2], and this naturally leads to the following intriguing conjecture.

Conjecture ($100 donation to the OEIS in honor of the prover or disprover): For every k ≥ 1, the number of distinct sequences A_{p}(n), as p ranges over all the (2k)!/(k!(k +1)!) 132-avoiding patterns of length k, is exactly p(k), the number of integer partitions of k.

[Of course the letter p in ‘p(k)’ has no relation whatsoever to the letter p in “A_{p}(n)”, except that both words, ‘partition’ and ‘pattern’, happen to start with it.]

Added March 24, 2014: It was just brought to our attention that Kate Rudolph ([R]) has already made this conjecture, and proved it, using elegant bijections, for k ≤ 7. Note that we proved it (of course, not as elegantly, but then again, we gave explicit closed-forms for all averages!) for k ≤ 10. We could have easily gone (with our modest computational resources) to k = 13.

Added March 31, 2014: Cheyne Homberger has just pointed out to us that our “conjecture” was already a theorem before it was made. Kate Rudolph’s conjecture was brilliantly solved by Lynn Chua and Krishanu Roy Sankar, in [CS]. We made the promised donation to the OEIS.

Ideally one would like to have not just explicit expressions for the average (expectation) of the random variable ‘number of occurrences of the pattern p’, i.e. A_{p}(n)/c_{n}, but explicit expressions (or failing this, efficient algorithms for generating many terms) for computing as many as possible higher moments.

Higher Moments

Every infinite sequence of sets, let’s call it {C_{n}}_{n=0}^{∞}, counted by the Catalan numbers, i.e. such that |C_{n}| = c_{n} is (most probably) so because of a natural structure-bijection [sometimes obvious (e.g. binary trees, Dyck paths), sometimes less so (e.g. 123-avoiding permutations)]

C_{n} ↔ ⋃_{k=1}^{n} C_{k−1} × C_{n−k} ,

leading immediately to the famous recurrence

c_{n} = ∑_{k=1}^{n} c_{k−1}c_{n−k} , c_{0} = 1 .

Many times, the members of our Catalan family have personalities, and posses numerical attributes (usually positive integers, but not necessarily), interchangeably called statistics and random variables. Let’s call such a statistic s → i(s). Then a natural question is

• What is the average, getting a brand-new numerical sequence

a_{n} := \frac{1}{c_{n}} ∑_{s∈C_{n}} i(s) .

But why stop here? For each power r, we may be interested in the r-th moment, getting yet another numerical sequence, one for each r,

m_{n}^{(r)} := \frac{1}{c_{n}} ∑_{s∈C_{n}}(i(s))^{r} .

For statisticians (and probabilists), more insightful sequences are moments about the mean,

M_{n}^{(r)} := \frac{1}{c_{n}} ∑_{s∈C_{n}}(i(s) − a_{n})^{r} ,

that of course, are easily computed from m_{n}^{(r)}, using the binomial theorem. For r > 2, in fact, the most interesting quantities are the standardized moments, aka alpha coefficients,

α_{n}^{(r)} = \frac{M_{n}^{(r)}}{(M_{n}^{(2)})^{r/2}} ,

that almost always converge, as n → ∞, to a sequence of real numbers, let’s call them β^{(r)}. When that happens, there is a limiting distribution, that in many cases (but not for Catalan families!) is the good-old Gaussian (aka normal) distribution, and that happens when β^{(r)} = 0 for r odd and β^{(r)} = 1 · 3 ⋯ (r − 1) = r!/(2^{r/2}(r/2)!), for r even.

Weighted-Counting of Catalan Objects According to a Statistic

Knowing all the moments is equivalent to knowing explicitly the generating function (aka weightenumerator) according to the statistic i(s), using the indeterminate t, getting a family of polynomials

P_{n}(t) := ∑_{s∈C_{n}} t^{i(s)} .

Once we know P_{n}(t) we can expand it in terms of t − 1

P_{n}(t) = ∑_{r=0}^{∞} \frac{f_{n}^{(r)}}{r!} (t − 1)^{r} ,

immediately getting

f_{n}^{(r)} := ∑_{s∈C_{n}} r!\binom{i(s)}{r} ,

from which the factorial moments can be gotten upon dividing by P_{n}(1) = c_{n}.

From the factorial moments, the usual moments can be easily computed (using Stirling numbers of the second kind).

Unlike the (numerical) enumerating sequence, c_{n}, that has a lovely closed-form, namely the famous

c_{n} = \frac{(2n)!}{n!(n + 1)!} ,

it is (usually) too much to hope for a closed-form expression for P_{n}(t), and in fact, their generating function is usually not even algebraic, i.e. the ‘grand-generating function’, w.r.t to z, say

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

usually does not satisfy an analogous algebraic equation to C(z) = 1 + zC(z)^{2}.

But it so happens (in many cases!), that the generating functions for the average (times c_{n}) and the r-th factorial (and hence actual) moments (again times c_{n}), for each specific (i.e. numeric) r are algebraic! In fact the polynomial equations satisfied by those generating functions often happen to be of degree 2, just like the one for the Catalan numbers, but, of course, with much more complicated coefficients. How can me find them?

Functional Recurrences to the Rescue

For purely pedagogical reasons, let’s first consider a very simple example. As in [B] and [J], our population is the set of 132-avoiding permutations, i.e. the set of permutations, π, of {1, …, n}, such that you never have 1 ≤ i_{1} < i_{2} < i_{3} ≤ n with π_{i_{1}} < π_{i_{3}} < π_{i_{2}}. Let’s first convince ourselves that this is indeed a Catalan family.

Take a typical such permutation, π, and look for the location of the largest entry, n. Suppose n stands at the k-th place, i.e. π_{k} = n. Then it is easy to see that all the entries standing to the left of n, i.e. {π_{1}, …, π_{k−1}} are all larger than all the entries standing to the right of n, namely {π_{k+1}, …, π_{n}}, or else a forbidden 132 pattern would emerge with the n playing the role of the ‘3’ in 132.

Hence every such permutation can be written as

π = π_{1}nπ_{2} ,

where π_{1} is a permutation of the set {n−k+1, …, n−1} and π_{2} is a permutation of {1, 2, …, n−k}. Of course, both π_{1} and π_{2} are 132-avoiding on their own right, and the map is a bijection. Hence the number of 132-avoiding permutations of {1, …, n} with π_{k} = n equals c_{k−1}c_{n−k}, and summing over 1 ≤ k ≤ n yields the Catalan recurrence, c_{n} = ∑_{k=1}^{n} c_{k−1}c_{n−k}, for the cardinality of the set of 132-avoiding permutations.

But now let’s consider the simple statistic ‘number of 21 patterns’.

Let a_{21}(π) be the number of 21 patterns of π. Using the above decomposition π = π_{1}nπ_{2}, we clearly have

a_{21}(π_{1}nπ_{2}) = a_{21}(π_{1}) + a_{21}(π_{2}) + k(n − k) ,

since a 21 pattern may either be entirely contained in π_{1}, entirely contained in π_{2}, or the ‘2’ may belong to π_{1} (k − 1 possibilities) and the ‘1’ may belong to π_{2} (n − k possibilities), so altogether (k −1)(n−k) possibilities, and of course, the 2 may be the ‘n’, and that gives n−k extra scenarios, so altogether we have (k − 1)(n − k) + (n − k) = k(n − k) additional occurrences of the patter 21.

Let’s define the weight-enumerator,

P_{n}(t) := ∑_{π∈AV_{132}(n)} t^{a_{21}(π)} .

For 1 ≤ k ≤ n, let AV_{132}^{(k)}(n) be the subset of AV_{132}(n) for which π_{k} = n, then of course

P_{n}(t) := ∑_{π∈AV_{132}(n)} t^{a_{21}(π)} = ∑_{k=1}^{n} ∑_{π∈AV_{132}^{(k)}(n)} t^{a_{21}(π)}

= ∑_{k=1}^{n} ∑_{{π_{1}∈AV_{132}(k−1) \atop π_{2}∈AV_{132}(n−k)}} t^{a_{21}(π_{1})+a_{21}(π_{2})+k(n−k)} = ∑_{k=1}^{n} t^{k(n−k)} ∑_{{π_{1}∈AV_{132}(k−1) \atop π_{2}∈AV_{132}(n−k)}} t^{a_{21}(π_{1})}t^{a_{21}(π_{2})}

= ∑_{k=1}^{n} t^{k(n−k)} (  ∑_{π_{1}∈AV_{132}(k−1)} t^{a_{21}(π_{1})}) (  ∑_{π_{2}∈AV_{132}(n−k)} t^{a_{21}(π_{2})})

= ∑_{k=1}^{n} t^{k(n−k)}P_{k−1}(t)P_{n−k}(t) ,

and hooray!, we found the (non-linear) recurrence equation

P_{n}(t) = ∑_{k=1}^{n} t^{k(n−k)}P_{k−1}(t)P_{n−k}(t) , P_{0}(t) = 1 ,  (NLR)

from which one can immediately get the first one hundred (or whatever) terms, but alas, no closed form.

Nevertheless, one can easily get explicit expressions for both the generating functions, and the sequences themselves, for the average and higher moments, m_{n}^{(r)} (times c_{n}), for numeric r, up to any desired r. Of course as r gets larger, the ‘explicit’ expressions would get more and more complicated, and there is (probably) no hope to get a symbolic expression in r, but we do what we can.

One way to derive expressions for higher moments is empirical. After you crank out the first 100 terms of the sequence P_{n}(t), guess (using, e.g. the built-in Maple package gfun developed by Salvy and Zimmerman[SaZ]) explicit expressions for the numerical sequences {P′_{n}(1)}, {P′′_{n}(1)}, etc. But one can proceed purely ‘rigorously’ as follows.

Suppose that we are only interested in the moments up to r ≤ R, then write,

P_{n}(1 + z) = ∑_{r=0}^{R} \frac{1}{r!}f_{n}^{(r)}z^{r} + O(z^{R+1}) .

Now plug this in into the above non-linear recurrence (NLR), with t replaced by 1 + z

∑_{r=0}^{R} \frac{1}{r!}f_{n}^{(r)}z^{r}+O(z^{r+1}) . = ∑_{k=1}^{n}(1+z)^{k(n−k)} (∑_{r=0}^{R} \frac{1}{r!}f_{k−1}^{(r)}z^{r} + O(z^{R+1}))) (∑_{s=0}^{R} \frac{1}{s!}f_{n−k}^{(s)}z^{s} + O(z^{R+1}))) .

Now use the binomial theorem to expand (1 + z)^{k(n−k)} to order R:

(1 + z)^{k(n−k)} = 1 + k(n − k)z + … + \binom{k(n − k)}{R}z^{R} + O(z^{R+1}) ,

and compare the coefficients of z^{0}, z, …, z^{R} on both sides, getting non-linear numerical recurrences.

Comparing the coefficient of z^{0} we get the good-old recurrence for c_{n} = f_{n}^{(0)}. Comparing the coefficients of z leads to the non-linear (numerical) recurrence for the sequence f_{n}^{(1)} that assumes that you already know f_{n}^{(0)} (as indeed you do, it being equal to the Catalan number c_{n}).

f_{n}^{(1)} = ∑_{k=1}^{n} k(n − k)f_{k−1}^{(0)}f_{n−k}^{(0)} + ∑_{k=1}^{n} f_{k−1}^{(1)}f_{n−k}^{(0)} + ∑_{k=1}^{n} f_{k−1}^{(0)}f_{n−k}^{(1)} ,

and so on and so forth.

From here we have a rigorous proof that a priori, the sequences f_{n}^{(1)}, f_{n}^{(2)}, etc. all have algebraic generating functions. One way is to teach the computer how to translate the recurrences into a system of algebraic equations for the generating functions, and then solve it, but a much more reasonable way is to use the non-linear recurrences (that now only involve numbers) to crank out sufficiently many terms to guess the algebraic generating functions, that can be justified a posteriori by plugging-in.

Analogously, for the pattern 12 one has the non-linear recurrence

P_{n}(t) = ∑_{k=1}^{n} t^{k−1}P_{k−1}(t)P_{n−k}(t) , P_{0}(t) = 1 ,

(why?).

Patterns of length 3

For the two patterns of length 2, namely 12 and 21 we got away with simple (non-linear) recurrences, but for patterns of length 3, we have two new concepts. The first is (non-linear) functional recurrence equation, and the second one is catalytic variable.

Let’s try and find an analogous recurrence for the weight-enumerators

P_{n}(t) := ∑_{π∈AV_{132}(n)} t^{a_{231}(π)} ,

(Note that this new P_{n}(t) is not the same as in the previous section, it is local notation).

So let’s try to express

a_{231}(π) = a_{231}(π_{1}nπ_{2}) ,

in terms of a_{231}(π_{1}) and a_{231}(π_{2}), where

π_{1} ∈ AV_{132}(k − 1)[{n − k + 1, …, n − 1}] ,

and

π_{2} ∈ AV_{132}(n − k)[{1, …, n − k}] .

Here are the possible scenarios for an occurrence of the pattern 231 in π = π_{1}nπ_{2}.

• it is completely immersed in π_{1} : a_{231}(π_{1}) ways.

• it is completely immersed in π_{2} : a_{231}(π_{2}) ways.

• the ‘23’ part of the pattern 231 belongs to π_{1} and the ‘1’ part belongs to π_{2}: a_{12}(π_{1}) · (n − k) ways.

• the ‘2’ part of the pattern 231 belongs to π_{1}, the ‘3’ is n, and the ‘1’ part belongs to π_{2}: (k − 1)(n − k) ways. Hence

a_{231}(π_{1}nπ_{2}) = a_{231}(π_{1}) + a_{231}(π_{2}) + a_{12}(π_{1}) · (n − k) + (k − 1)(n − k) .

Alas, we have an uninvited guest, a_{12}(π), so we need to figure out how to express it in terms of a_{12}(π_{1}) and a_{12}(π_{2}), but that’s easy

a_{12}(π_{1}nπ_{2}) = a_{12}(π_{1}) + a_{12}(π_{2}) + k − 1 .

In order to figure out a recurrence for P_{n}(t), we need to introduce a catalytic variable, q, that takes care of the r.v. a_{12}(π) and define:

Q_{n}(t, q) := ∑_{π∈AV_{132}(n)} t^{a_{231}(π)}q^{a_{12}(π)} .

At the end of the day, once we have Q_{n}(t, q), we would plug-in q = 1 and get our desired P_{n}(t) = Q_{n}(t, 1), but until then we would have to put-up with q. It would be convenient to define a weight

W t(π)(t, q) := t^{a_{231}(π)}q^{a_{12}(π)} .

So we have W t(π_{1}nπ_{2})(t, q) := t^{a_{231}(π)}q^{a_{12}(π)} = t^{a_{231}(π_{1})+a_{231}(π_{2})+a_{12}(π_{1})·(n−k)+(k−1)(n−k)}q^{a_{12}(π_{1})+a_{12}(π_{2})+k−1}

= t^{(k−1)(n−k)}q^{k−1} · ( t^{a_{231}(π_{1})+a_{12}(π_{1})·(n−k)}q^{a_{12}(π_{1})} ) · ( t^{a_{231}(π_{2})}q^{a_{12}(π_{2})} )
= t^{(k−1)(n−k)}q^{k−1} · ( t^{a_{231}(π_{1})}(qt^{n−k})^{a_{12}(π_{1})} ) · ( t^{a_{231}(π_{2})}q^{a_{12}(π_{2})} )
= t^{(k−1)(n−k)}q^{k−1} · W t(π_{1})(t, t^{n−k}q) · W t(π_{2})(t, q) ,

leading to the functional (non-linear) recurrence

Q_{n}(t, q) = ∑_{k=1}^{n} t^{(k−1)(n−k)}q^{k−1} Q_{k−1}(t, t^{n−k}q) Q_{n−k}(t, q) , Q_{0}(t, q) = 1 .

Similarly, for the other four 132-avoiding patterns, we have:

• a_{123}(π) (with catalytic variable q corresponding to a_{12}(π)):

Q_{n}(t, q) = ∑_{k=1}^{n} q^{k−1} Q_{k−1}(t, tq) Q_{n−k}(t, q) , Q_{0}(t, q) = 1 .

• a_{321}(π) (with catalytic variable q corresponding to a_{21}(π)):

Q_{n}(t, q) = ∑_{k=1}^{n} q^{k(n−k)} Q_{k−1}(t, t^{n−k}q) Q_{n−k}(t, t^{k}q) , Q_{0}(t, q) = 1 .

• a_{213}(π) (with catalytic variable q corresponding to a_{21}(π)):

Q_{n}(t, q) = ∑_{k=1}^{n} q^{k(n−k)} Q_{k−1}(t, tq) Q_{n−k}(t, q) , Q_{0}(t, q) = 1 .

• a_{312}(π) (with catalytic variable q corresponding to a_{12}(π)):

Q_{n}(t, q) = ∑_{k=1}^{n} q^{k−1}Q_{k−1}(t, q)Q_{n−k}(t, t^{k}q) , Q_{0}(t, q) = 1 .

All the above can be (and has been!) ‘taught’ to the computer, and the computer can automatically derive such functional equations. The above verbose derivation was only for the benefit of explaining to humans the algorithms that would eventually be executed by computers.

Unfortunately, one can’t get closed-form expressions for the Q_{n}(t, q), and not even a closed-form expression for their generating function, but the above functional recurrences are fairly efficient for generating quite a few terms, and by plugging-in q = 1 and taking successive derivatives with respect to t and then plugging-in t = 1 one can generate quite large beginnings of the moment-sequences, and have a guessing program (e.g. gfun [SaZ]) guess either closed-forms or recurrences.

But for those who abhor guessing, one can get, completely automatically, algebraic equations satisfied by the generating functions for the factorial moments, like we did above for the P_{n}(t) of a_{21}(π). Now we need multi-variable Taylor expansions, and need to put up with more multi-indexed sequences, but so what? The computer does not mind!

How to represent Functional Recurrences in Maple?

The beauty of mathematics, (and computers!) is that we can generalize and consider a very general class of functional recurrences that include all the above as very spacial cases.

Here goes:

P_{n}(t_{1}, …, t_{m}) =

∑_{k=1}^{n} c(t_{1}, …, t_{m}, k, n)P_{k−1}(a_{1}(t_{1}, …, t_{m}), …, a_{m}(t_{1}, …, t_{m}))·

P_{n−k}(b_{1}(t_{1}, …, t_{m}), …, b_{m}(t_{1}, …, t_{m})) , P_{0}(t_{1}, …, t_{m}) = 1 .

Here the a_{i}’s and b_{i}’s are arbitrary polynomials in their variables t_{1}, …, t_{m}, and c(t_{1}, …, t_{m}, k, n) is a polynomial in t_{1}, …, t_{m}, k, m but where the powers of the t_{i}’s (but not n, k) are allowed to be polynomials in n, k. For example c(t_{1}, t_{2}, k, n) = n^{2}t_{1}^{n+k^{2}}t_{2}^{n^{3}} + k^{3}t_{2}^{k^{3}+n} is quite acceptable, but c(t_{1}, t_{2}, k, n) = k^{k} is not.

In the Maple package AlgFunEq such a Functional Equation is represented by the following data structure

[c, [a_{1}, …, a_{m}], [b_{1}, …, b_{m}]], [t_{1}, …, t_{m}]] .

From this one can get, automatically, an efficient scheme for computing the mixed factorial moments (and hence the pure ones). This is accomplished by procedure FAscheme.

More General Functional Recurrences

As general as the above form of the functional Catalan-type recurrences is, it is not general enough to consider many patterns of length larger than 3. Let’s take for example a_{4321}.

We have

a_{4321}(π_{1}nπ_{2}) = a_{4321}(π_{1}) + a_{4321}(π_{2}) + ka_{321}(π_{2}) + a_{21}(π_{1})a_{21}(π_{2}) + (n − k)a_{321}(π_{1}) .

It is easy to see, because of the product term a_{21}(π_{1})a_{21}(π_{2}), that an analogous derivation to the one for a_{231}(π) carried above does not lead to such a functional recurrence.

But what we can do is introduce “generalized products”. First defining it on pairs of monomials, that the computer defines automatically, according to its needs, and then extends it by bi-linearity to apply to any pair of polynomials. Calling this product F, we have to handle functional equations of the more general form:

P_{n}(t_{1}, …, t_{m}) =

∑_{k=1}^{n} c(t_{1}, …, t_{m}, k, n) · F(P_{k−1}(a_{1}(t_{1}, …, t_{m}), …, a_{m}(t_{1}, …, t_{m})) , P_{n−k}(b_{1}(t_{1}, …, t_{m}), …, b_{m}(t_{1}, …, t_{m}))) ,

P_{0}(t_{1}, …, t_{m}) = 1 .

We confess that we were too lazy to implement these more general types in order to compute higher moments, but for the special cases of averages (and that’s all that Bóna did for patterns of length

3) it is fully implemented for any pattern, see procedure MilonK.

Note that even in this more general setting, we are guaranteed that the generating function of each specific moment is algebraic, and hence it justifies the empirical guessing.

Using these, we were able to find closed-form expressions for the averages for all patterns of length ≤ 10, mentioned at the beginning of this article.

Rigorous Proofs of Asymptotic Abnormality

Using the output of http://www.math.rutgers.edu/~zeilberg/tokhniot/oAlgFunEq3 we see

(rigorously!) that the random variables π → a_{p}(π) are never asymptotically normal, for patterns of length ≤ 3, in accordance with Janson’s humanly-generated paper ([J]) (who proved it for patterns of all lengths). That file also contains asymptotic expressions, to order 4 of all standardized moments up to the sixth.

What About The Number of Pattern-Occurrences in 123-Avoiding Permutations?

The Catalan Structure of 123-avoiding permutations is a bit more subtle. For any 123-avoiding permutation π, let U(π) be the permutation obtained by finding the right-to-left maxima, circling them, introducing an empty slot right before the first entry, sliding all the non-left-to-right maxima one unit to the left, then increasing all the entries by 1, and finally sticking a 1 at the remaining open slot.

It is (almost) readily seen that any permutation π of {1, …, n} can be written, for some k, 1 ≤ k ≤ n

π = π_{1}U(π_{2}) ,

where π_{1} is a permutation of {k + 1, …, n} and π_{2} is a permutation of {1, …, k − 1}.

Cheyne Homberger[H] found generating functions, and explicit expressions, for the averages of the random variables “number of occurrences of pattern p” defined on the set of 123-avoiding permutations, for all patterns p of length ≤ 3.

Let’s be more general and try to find a functional recurrence equation for the weight-enumerator

P_{n}(t) := ∑_{π∈AV_{123}(n)} t^{a_{213}(π)} .

It turns out that we need two auxiliary statistics, and hence two catalytic variables: σ_{1}(π) := a_{213}(U(π)) − a_{213}(π) ,

σ_{2}(π) := σ_{1}(U(π)) − σ_{1}(π) ,

and let’s define the sequence of polynomials

Q_{n}(t, s_{1}, s_{2}) := ∑_{π∈AV_{123}(n)} t^{a_{213}(π)}s_{1}^{σ_{1}(π)}s_{2}^{σ_{2}(π)} ,

then the reader is welcome to prove the (slightly) ‘weird’ functional recurrence

Q_{n}(t, s_{1}, s_{2}) = s_{2} Q_{n−1}(t, s_{1}, s_{2}) + ∑_{k=2}^{n} Q_{n−k}(t, s_{1}, s_{2}) Q_{k−1}(t, ts_{1}, s_{1}s_{2}) ,

subject to the initial condition Q_{0}(t, s_{1}, s_{2}) = 1. From this functional recurrence, one can compute quite a few terms. Then setting s_{1} = 1, s_{2} = 1, we get P_{n}(t) = Q_{n}(t, 1, 1). See:

http://www.math.rutgers.edu/~zeilberg/tokhniot/oCheyne2 ,

where one can find the list of P_{n}(t) for 0 ≤ n ≤ 35. Using this, we easily reproduced (a more compact version of) Homberger’s explicit expression for what he called a_{n} in [H], that is:

∑_{π∈AV_{123}(n)} a_{213}(π) .

The expression is:

−\frac{3}{8} 4^{n} + \frac{1}{2} (n + 2) (2 n − 1) c_{n−1} ,

(recall that c_{n} are the Catalan numbers). Going beyond, we found an explicit expression for the second moment (times c_{n}), i.e. for:

∑_{π∈AV_{123}(n)}(a_{213}(π))^{2} ,

that turned out to be:

−\frac{9}{128} (3 n^{2} + 7 n + 6) 4^{n} + (\frac{19}{60} n^{4} + \frac{57}{20} n^{3} + \frac{67}{30} n^{2} + \frac{1}{10} n − 1) c_{n−1} ,

but to our disappointment, the third moment turned out to be not nice at all, and all we could find was a linear recurrence equation of order 4 with coefficients that are polynomials of degree 4, see the output file

http://www.math.rutgers.edu/~zeilberg/tokhniot/oCheyne3  .

All this was generated with the Maple package Cheyne available from

http://www.math.rutgers.edu/~zeilberg/tokhniot/Cheyne  .

It should be be possible to expand this Maple package to handle larger patterns, but enough is enough.

Empirical Coda

Recall that at the beginning of this article we conjectured that the number of ‘Bóna classes’ for 132-avoiding permutations, i.e. the number of different sequences that show up as averages of the random variables ‘number of occurrences of the pattern p’ for all 132 patterns of length k is the number of integer partitions of k, and we proved it rigorously for k ≤ 10.

The output file

http://www.math.rutgers.edu/~zeilberg/tokhniot/oCheyne1

contains all the beginnings (through n = 9) of the analogous sequences for 123-avoiding permutations, and it is very possible that the sequence ‘total number of Bóna classes’ for patterns of length k defined on the set of 123-avoiding permutations, starting at k = 1 is:

1, 2, 3, 6, 12, 32, … ,

but we can’t see a pattern. Can you?

References

[B] Miklós Bóna, Surprising symmetries in objects counted by Catalan numbers, The Electronic Journal of Combinatorics 19(1) (2012), #P62.

http://www.combinatorics.org/ojs/index.php/eljc/article/download/v19i1p62/pdf

[CS] Lynn Chua and Krishanu Roy Sankar, Equipopularity Classes of 132-Avoiding Permutations, The Electronic Journal of Combinatorics 21(1)(2014), #P59

http://www.combinatorics.org/ojs/index.php/eljc/article/download/v21i1p59/pdf

[H] Cheyne Homberger, Expected Patterns in Permutation Classes, The Electronic Journal of Combinatorics 19(3)(2012), #P43

http://www.combinatorics.org/ojs/index.php/eljc/article/download/v19i3p43/pdf

[J] Svante Janson, Patterns in random permutations avoiding the pattern 132, arxiv.org, 22 Jan.

2014, http://arxiv.org/abs/1401.5679

[R] Kate Rudolph, Pattern Popularity in 132-Avoiding Permutations, The Electronic Journal of Combinatorics 20(1)(2013), #P8

http://www.combinatorics.org/ojs/index.php/eljc/article/download/v20i1p8/pdf

[SaZ] Bruno Salvy and Paul Zimmerman, GFUN: a Maple package for the manipulation of generating and holonomic functions in one variable, ACM Transactions on Mathematical Software 20(1994), 163-177 .

[Sl1] Neil Sloane, Sequence #108, http://oeis.org/A000108 .

[Sl2] Neil Sloane, Sequence #41, http://oeis.org/A000041 .

[St] Richard Stanley, Catalan Addendum, http://www-math.mit.edu/~rstan/ec/catadd.pdf .
[Z] Doron Zeilberger, Opinion #49, http://www.math.rutgers.edu/~zeilberg/Opinion49.html