###################################################################### ##AsymptoticMoments: Save this file as AsymptoticMoments # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read AsymptoticMoments # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: Nov. 12, 2008 print(`Created: Nov. 12, 2008`): print(` This is AsymptoticMoments. `): print(`It is one of the two packages accompanying the article `): print(`The Automatic Central Limit Theorems Generator (and Much More!) `): print(`by Doron Zeilberger`): print(`It finds explicit expressions for the General Moments of`): print(`Large families of Discrete Probabiliy Distributions`): print(`leading immediately to Automatic Proofs of Generalized`): print(`Limit Theorems, and refinements thereof`) : print(`by Doron Zeilberger.`): print(`and also available from Zeilberger's website`): print(``): print(`Please report bugs to zeilberg at math dot rutgers dot edu`): print(``): print(`The most current version of this package and paper`): print(` are available from`): print(`http://www.math.rutgers.edu/~zeilberg/ .`): print(`For a list of the procedures type ezra();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): with(combinat): ezra1:=proc() if args=NULL then print(` The supporting procedures are: `): print(` Centralize, CheckMs, CheckRec, FM1, FMs, FMsAnsatz,`): print(`GuessPol, M1, MekT, Ms,MsAnsatz, NorFMsE, NorFMsO,`): print(`NorMsE, NorMsO,PGF, `): print(`CheckAsyFMsEempir,CheckAsyFMsOempir,`): print(`CheckAsyMsEempir,CheckAsyMsOempir,`): print(` Stir2 `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: `): print(` AsyFMsE, AsyFMsO,AsyMsE, AsyMsO, AsyNorFMsE, AsyNorFMsO`): print(`AsyNorMsE,AsyNorMsO, CheckAsyFMsEmpir, CheckAsyMsEmpir, PPk `): print(` Sipur, SipurFairDie, SipurPP`): elif nops([args])=1 and op(1,[args])=AsyMsE then print(`AsyMsE(P,q,n,r,s): The first s+1 terms `): print(`of the conjectured asymptotic expansion, in n, for the `): print(`2r-th moment`): print(`of the (centralized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`AsyMsE(p*q+(1-p),q,n,r,3);`): print(`AsyMsE((1-q^(n+1))/(1-q),q,n,r,3);`): print(`AsyMsE((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3);`): elif nops([args])=1 and op(1,[args])=AsyMsO then print(`AsyMsO(P,q,n,r,s): The first s+1 terms `): print(`of the conjectured asymptotic expansion, in n, for the `): print(`(2r+1)-th factorial moment`): print(`of the (centralized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`AsyMsO(p*q+(1-p),q,n,r,3);`): print(`AsyMsO((1-q^(n+1))/(1-q),q,n,r,3);`): print(`AsyMsO((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3);`): elif nops([args])=1 and op(1,[args])=AsyFMsE then print(`AsyFMsE(P,q,n,r,s): The first s+1 terms `): print(`of the conjectured asymptotic expansion, in n, for the `): print(`2r-th factorial moment`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`AsyFMsE(p*q+(1-p),q,n,r,3);`): print(`AsyFMsE((1-q^(n+1))/(1-q),q,n,r,3);`): print(`AsyFMsE((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3);`): elif nops([args])=1 and op(1,[args])=AsyFMsO then print(`AsyFMsO(P,q,n,r,s): The first s+1 terms `): print(`of the conjectured asymptotic expansion, in n, for the `): print(`(2r+1)-th factorial moment`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`AsyFMsO(p*q+(1-p),q,n,r,3);`): print(`AsyFMsO((1-q^(n+1))/(1-q),q,n,r,3);`): print(`AsyFMsO((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3);`): elif nops([args])=1 and op(1,[args])=AsyNorFMsE then print(`AsyNorFMsE(P,q,n,r,s): The first s+1 terms `): print(`of the conjectured asymptotic expansion, in n, for the `): print(`Normalized 2r-th factorial moment`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`The normalization is obtained by dividing by `): print(`the r-th power of the variance divided by (2r)!/r!/2^r .`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`AsyNorFMsE(p*q+(1-p),q,n,r,3);`): print(`AsyNorFMsE((1-q^(n+1))/(1-q),q,n,r,3);`): print(`AsyNorFMsE((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3);`): elif nops([args])=1 and op(1,[args])=AsyNorFMsO then print(`AsyNorFMsO(P,q,n,r,s): The first s+1 terms `): print(`of the conjectured asymptotic expansion, in n, for the `): print(`Normalized (2r+1)-th factorial moment`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`The normalization is obtained by dividing by `): print(`the r-th power of the variance divided by (2r)!/r!/2^r .`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`AsyNorFMsO(p*q+(1-p),q,n,r,3);`): print(`AsyNorFMsO((1-q^(n+1))/(1-q),q,n,r,3);`): print(`AsyNorFMsO((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3);`): elif nops([args])=1 and op(1,[args])=AsyNorMsO then print(`AsyNorMsO(P,q,n,r,s): The first s+1 terms in the`): print(`asymptotic expansion of the normalized (2r+1)-th moment`): print(`of the r.v. induced by P(q,q^n).`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`AsyNorMsO(p*q+(1-p),q,n,r,2);`): print(`AsyNorMsO((1-q^(n+1))/(1-q),q,n,r,3);`): print(`AsyNorMsO((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3);`): elif nops([args])=1 and op(1,[args])=AsyNorMsE then print(`AsyNorMsE(P,q,n,r,s): The first s+1 terms `): print(`of the conjectured asymptotic expansion, in n, for the `): print(`2r-th moment`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`AsyNorMsE(p*q+(1-p),q,n,r,3);`): print(`AsyNorMsE((1-q^(n+1))/(1-q),q,n,r,3);`): print(`AsyNorMsE((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3);`): elif nops([args])=1 and op(1,[args])=Centralize then print(`Centralize(P,q,n): Given a rational function`): print(`P, depending on q and q^n, centralizes and normalizes it`): print(`(see the paper for details). For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`Centralize(p*q+(1-p),q,n);`): print(`Centralize((1-q^(n+1))/(1-q),q,n);`): print(`Centralize((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n);`): print(``): elif nops([args])=1 and op(1,[args])=CheckMs then print(`CheckMs(P,q,n,R,S): checks MS(P,q,n,R) againt`): print(`the first S values of M1`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`CheckMs(p*q+(1-p),q,n,10,15);`): print(`CheckMs((1-q^(n+1))/(1-q),q,n,4,10);`): print(`CheckMs((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6,10);`): elif nops([args])=1 and op(1,[args])=CheckRec then print(`CheckRec(P,q,n,R): checks the recurrence for the factorial`): print(`moments. For example, try:`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`CheckRec(p*q+(1-p),q,n,10);`): print(`CheckRec((1-q^(n+1))/(1-q),q,n,4);`): print(`CheckRec((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6);`): elif nops([args])=1 and op(1,[args])=FM1 then print(`FM1(P,q,n,n1,R): The list of the first R factorial moments`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`for the case n1, i.e. where the pgf is`): print(` (a centralized-normalized version`): print(` of) P(q,q^1)* ...* P(q,q^n1) `): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`FM1(p*q+(1-p),q,n,5,10);`): print(`FM1((1-q^(n+1))/(1-q),q,n,3,10);`): print(`FM1((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,2, 12);`): elif nops([args])=1 and op(1,[args])=FMsAnsatz then print(`FMsAnsatz(P,q,n,R): The list of the first R factorial moments`): print(`as polynomials expressions in the variable n`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`Using the polynomial ansatz way (see also FMs that is much faster)`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`FMsAnsatz(p*q+(1-p),q,n,10);`): print(`FMsAnsatz((1-q^(n+1))/(1-q),q,n,4);`): print(`FMsAnsatz((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6);`): elif nops([args])=1 and op(1,[args])=FMs then print(`FMs (P,q,n,R): The list of the first R factorial moments`): print(`as polynomials expressions in the variable n`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`Using the recurrence way (much faster than FMsAnsatz`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`FMs(p*q+(1-p),q,n,10);`): print(`FMs((1-q^(n+1))/(1-q),q,n,4);`): print(`FMs((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6);`): elif nops([args])=1 and op(1,[args])=GuessPol then print(`GuessPol(L,n): guesses a polynomial of degree d in n for`): print(` the list L, such that P(i)=L[i] for i=1..nops(L) for example, try: `): print(`GuessPol([seq(i,i=1..10)],n);`): elif nops([args])=1 and op(1,[args])=MekT then print(`MekT(P,q,n,R): The first R Taylor coefficients in`): print(`P1(1+z) where P1 is the centralized-normalized version`): print(`of P, a rational function in q and q^n describing`): print(`a discrete family of discrete prob. distribution`): print(`(see the article for details).`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`MekT(p*q+(1-p),q,n,10);`): print(`MekT((1-q^(n+1))/(1-q),q,n,4);`); print(`MekT((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6);`): elif nops([args])=1 and op(1,[args])=Ms then print(`Ms(P,q,n,R): The list of the first R moments (about the mean)`): print(`as polynomials expressions in the variable n`): print(`of the (normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`using the expressions for the Factorial Moments`): print(`produced by procedure FMs, that uses the recurrence`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`Ms(p*q+(1-p),q,n,10);`): print(`Ms((1-q^(n+1))/(1-q),q,n,4);`): print(`Ms((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6);`): elif nops([args])=1 and op(1,[args])=MsAnsatz then print(`MsAnsatz(P,q,n,R): The list of the first R moments (about the mean)`): print(`as polynomials expressions in the variable n`): print(`of the (normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`using the polynomial ansatz`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`MsAnsatz(p*q+(1-p),q,n,10);`): print(`MsAnsatz((1-q^(n+1))/(1-q),q,n,4);`): print(`MsAnsatz((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6);`): elif nops([args])=1 and op(1,[args])=NorFMsE then print(`NorFMsE(P,q,n,R): The list of the first Normalized R even `): print(`factorial moments`): print(`as polynomials expressions in the variable n`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`The r-th element in the output is the 2r-th moment divided by`): print(`[the r-th power of the variance times (2r)!/r!/2^r ].`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`NorFMsE(p*q+(1-p),q,n,10);`): print(`NorFMsE((1-q^(n+1))/(1-q),q,n,4);`): print(`NorFMsE((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6);`): elif nops([args])=1 and op(1,[args])=NorFMsO then print(`NorFMsO(P,q,n,R): The list of the first Normalized R odd`): print(`factorial moments`): print(`as polynomials expressions in the variable n`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`The r-th element in the output is the (2r+1)-th moment divided by`): print(`[the r-th power of the variance times (2r)!/r!/2^r ].`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`NorFMsO(p*q+(1-p),q,n,10);`): print(`NorFMsO((1-q^(n+1))/(1-q),q,n,4);`): print(`NorFMsO((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6);`): elif nops([args])=1 and op(1,[args])=NorMsE then print(`NorMsE(P,q,n,R): The list of the first Normalized R even `): print(`moments`): print(`as polynomials expressions in the variable n`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`The r-th element in the output is the 2r-th moment divided by`): print(`[the r-th power of the variance times (2r)!/r!/2^r ].`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`NorMsE(p*q+(1-p),q,n,10);`): print(`NorMsE((1-q^(n+1))/(1-q),q,n,4);`): print(`NorMsE((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6);`): elif nops([args])=1 and op(1,[args])=NorMsO then print(`NorMsO(P,q,n,R): The list of the first Normalized R odd`): print(` moments`): print(`as polynomials expressions in the variable n`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`The r-th element in the output is the (2r+1)-th moment divided by`): print(`[the r-th power of the variance times (2r)!/r!/2^r ].`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`NorMsO(p*q+(1-p),q,n,10);`): print(`NorMsO((1-q^(n+1))/(1-q),q,n,4);`): print(`NorMsO((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6);`): elif nops([args])=1 and op(1,[args])=PGF then print(`PGF(P,q,n,n1): the probability generating function`): print(`induced by the rational function P(q,q^n), namely:`): print(`P(q,q^1)*...*P(q,q^n1)`): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`PGF(p*q+(1-p),q,n,5);`): print(`PGF((1-q^(n+1))/(1-q),q,n,5);`): print(`PGF((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,5);`): elif nops([args])=1 and op(1,[args])=M1 then print(`M1(P,q,n,n1,R): The list of the first R moments`): print(`of the (centralized and normalized) probability distribution`): print(`induced by the rational function P(q,q^n)`): print(`for the case n1, i.e. where the pgf is`): print(` (a centralized-normalized version`): print(` of) P(q,q^1)* ...* P(q,q^n1) `): print(`For example, try:`): print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`M1(p*q+(1-p),q,n,5,10);`): print(`M1((1-q^(n+1))/(1-q),q,n,3,10);`): print(`M1((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,2, 12);`): elif nops([args])=1 and op(1,[args])=PPk then print(`PPk(k,q,n): The P(q^n,q) corresponding to the generating function`): print(`of plane partitions whose 3D Ferrers diagram is inside`): print(` an n by n by k box`): print(`Here k is a pos. integer, and q and n are symbolic.`): print(`For example, try:`): print(`PPk(2,q,n)`): elif nops([args])=1 and op(1,[args])=CheckAsyFMsEmpir then print(`CheckAsyFMsEmpir(P,q,n,s,K): Emprically yet Rigorously `): print(`if K is chosen large enough, `): print(`proves the`): print(`conjectured asymptotic expression given by both AsyFMsO(P,q,n,r,s)`): print(`and AsyFMsE(P,q,n,r,s), for the odd and even factorial moments`): print(`respectively. For example, try:`); print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`CheckAsyFMsEmpir(p*q+(1-p),q,n,2,7);`): print(`CheckAsyFMsEmpir((1-q^(n+1))/(1-q),q,n,3,7);`): print(`CheckAsyFMsEmpir((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,3,7);`): elif nops([args])=1 and op(1,[args])=CheckAsyFMsEempir then print(`CheckAsyFMsEempir(P,q,n,s,K): Emprically yet Rigorously `): print(`if K is chosen large enough, `): print(`proves the`): print(`conjectured asymptotic expression given by AsyFMsO(P,q,n,r,s)`): print(`For example, try:`); print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`CheckAsyFMsOempir(p*q+(1-p),q,n,2,7);`): print(`CheckAsyFMsOempir((1-q^(n+1))/(1-q),q,n,3,7);`): print(`CheckAsyFMsOempir((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,3,7);`): elif nops([args])=1 and op(1,[args])=CheckAsyFMsOempir then print(`CheckAsyFMsOempir(P,q,n,s,K): Emprically yet Rigorously `): print(`if K is chosen large enough, `): print(`proves the`): print(`conjectured asymptotic expression given by AsyFMsO(P,q,n,r,s)`): print(`For example, try:`); print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`CheckAsyFMsOempir(p*q+(1-p),q,n,2,7);`): print(`CheckAsyFMsOempir((1-q^(n+1))/(1-q),q,n,3,7);`): print(`CheckAsyFMsOempir((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,3,7);`): elif nops([args])=1 and op(1,[args])=CheckAsyMsEmpir then print(`CheckAsyMsEmpir(P,q,n,s,K): Emprically yet Rigorously `): print(`if K is chosen large enough, `): print(`proves the`): print(`conjectured asymptotic expression given by both AsyMsO(P,q,n,r,s)`): print(`and AsyMsE(P,q,n,r,s), for the odd and even factorial moments`): print(`respectively. For example, try:`); print(`(for tossing a coin, the Mahonian statistics, and q-Catalan resp.):`): print(`CheckAsyMsEmpir(p*q+(1-p),q,n,2,7);`): print(`CheckAsyMsEmpir((1-q^(n+1))/(1-q),q,n,3,7);`): print(`CheckAsyMsEempir((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,3,7);`): elif nops([args])=1 and op(1,[args])=Sipur then print(`Sipur(n,r,s): A verbose description of `): print(`aymptotic expansions for the general (2r and 2r+1)`): print(`moment (about the mean) of`): print(`some important combinatorial distributions`): print(`proving, in particular, that they asymptotically normal`): print(`Do:`): print(`Sipur(n,r,2);`): elif nops([args])=1 and op(1,[args])=SipurFairDie then print(`SipurFairDie(n,r,K,s): A verbose description of `): print(`aymptotic expansions, to s+1 terms, of the general `): print(`moment (about the mean) of the`): print(`discrete probability distributions `): print(`obtained by rolling a fair k-sided die for k from 2 to K`): print(`or example, try`): print(`SipurFairDie(n,r,6,2);`): elif nops([args])=1 and op(1,[args])=SipurPP then print(`SipurPP(n,r,K,s): A verbose description of `): print(`aymptotic expansions, to s+1 terms, of the general `): print(`moment (about the mean) of the`): print(`discrete probability distributions `): print(`for plane partitions bounded in an n by n by k`): print(`box for fixed k (k=1..K),`): print(`Do:`): print(`SipurPP(n,r,3,2);`): elif nops([args])=1 and op(1,[args])=Stir2 then print(`Stir2(r,i): The polynomial in r that is Stirling2(r,r-i)`): print(`For example try: Stir2(r,2):`): else print(`There is no ezra for`,args): fi: end: #Centralize(P,q,n): Given a rational function #P, depending on q and q^n, centralizes and normalizes it #(see the paper for details). #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #Centralize(p*q+(1-p),q,n); #Centralize((1-q^(n+1))/(1-q),q,n); #Centralize((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n); Centralize:=proc(P,q,n) local gu0,gu1,P1: gu0:=normal(limit(P,q=1)): P1:=P/gu0: gu1:=normal(limit(normal(diff(P1,q)),q=1)): P1/q^gu1: end: #PGF(P,q,n,n1): the probability generating function #induced by the rational function P(q,q^n), namely: #P(q,q^1)*...*P(q,q^n1) #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.`): #PGF(p*q+(1-p),q,n,5); #PGF((1-q^(n+1))/(1-q),q,n,5); #PGF((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,5); PGF:=proc(P,q,n,n1) local P1: option remember: P1:=Centralize(P,q,n): if n1=0 then RETURN(1): else normal(PGF(P,q,n,n1-1)*subs(n=n1,P1)): fi: end: ##begin guessing polynomials #GuessPol1(L,d,n): guesses a polynomial of degree d in n for # the list L, such that P(i)=L[i+1] for i=1..nops(L) #For example, try: #GuessPol1([seq(i,i=1..10)],1,n); GuessPol1:=proc(L,d,n) local P,i,a,eq,var: if d>nops(L)-2 then ERROR(`the list is too small`): fi: P:=add(a[i]*n^i,i=0..d): var:={seq(a[i],i=0..d)}: eq:={seq(subs(n=i,P)-L[i],i=1..nops(L))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: subs(var,P): end: #GuessPol(L,n): guesses a polynomial of degree d in n for # the list L, such that P(i)=L[i] for i=1..nops(L) for example, try: #GuessPol([seq(i,i=1..10)],n); GuessPol:=proc(L,n) local d,gu: for d from 0 to nops(L)-2 do gu:=GuessPol1(L,d,n): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: ##end guessing polynomials #FM1(P,q,n,n1,R): The list of the first R factorial moments #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #for the case n1, i.e. where the pgf is (a centralized-normalized version #of) P(q,q^1)* ...* P(q,q^n1) #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #FM1(p*q+(1-p),q,n,5,10); #FM1((1-q^(n+1))/(1-q),q,n,3,10); #FM1((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,2, 12); FM1:=proc(P,q,n,n1,R) local gu,lu,r: option remember: gu:=PGF(P,q,n,n1): lu:=[]: for r from 1 to R do gu:=normal(diff(gu,q)): lu:=[op(lu),subs(q=1,gu)]: od: lu: end: #M1(P,q,n,n1,R): The list of the first R moments (about the mean) #of the probability distribution #induced by the rational function P(q,q^n) #for the case n1, i.e. where the pgf is (a centralized-normalized version #of) P(q,q^1)* ...* P(q,q^n1) #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #M1(p*q+(1-p),q,n,5,10); #M1((1-q^(n+1))/(1-q),q,n,3,10); #M1((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,2, 12); M1:=proc(P,q,n,n1,R) local gu,lu,r: option remember: gu:=PGF(P,q,n,n1): lu:=[]: for r from 1 to R do gu:=normal(q*diff(gu,q)): lu:=[op(lu),subs(q=1,gu)]: od: lu: end: #FMsAnsatz(P,q,n,R): The list of the first R factorial moments #as polynomials expressions in the variable n #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #Here we use guessing by the polynomial ansatz. #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #FMs(p*q+(1-p),q,n,10); #FMs((1-q^(n+1))/(1-q),q,n,4); #FMsAnsatz((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6); FMsAnsatz:=proc(P,q,n,R) local gu,lu,r,mu,n1: option remember: lu:=[seq(FM1(P,q,n,n1,R),n1=1..2*R+2)]: gu:=[]: for r from 1 to R do mu:=GuessPol([seq(lu[n1][r],n1=1..2*R+2)],n): if mu=FAIL then RETURN(FAIL,gu): else gu:=[op(gu),mu]: fi: od: gu: end: #MsAnsatz(P,q,n,R): The list of the first R moments #as polynomials expressions in the variable n #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #MsAnsatz(p*q+(1-p),q,n,10); #MsAnsatz((1-q^(n+1))/(1-q),q,n,4); #MsAnsatz((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6); MsAnsatz:=proc(P,q,n,R) local gu,lu,r,mu,n1: option remember: lu:=[seq(M1(P,q,n,n1,R),n1=1..2*R+2)]: gu:=[]: for r from 1 to R do mu:=GuessPol([seq(lu[n1][r],n1=1..2*R+2)],n): if mu=FAIL then RETURN(FAIL,gu): else gu:=[op(gu),mu]: fi: od: gu: end: sed1:=proc(f,x,a) local i,f1: f1:=f: if normal(subs(x=a,f1))<>0 then RETURN(0): fi: for i from 1 do f1:=diff(f1,x): if normal(subs(x=a,f1))<>0 then RETURN(i): fi: od: end: #MekT(P,q,n,R): The first R Taylor coefficients in #P1(1+z) where P1 is the centralized-normalized version #of P, a rational function in q and q^n describing #a discrete family of discrete prob. distribution #(see the article for details). #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #MekT(p*q+(1-p),q,n,10); #MekT((1-q^(n+1))/(1-q),q,n,4); #MekT((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6); MekT:=proc(P,q,n,R) local P1,P1a,z,i,shulaim,mone,mekh: P1:=normal(Centralize(P,q,n)): P1a:=subs(n=13,P1): mone:=numer(P1a): mekh:=denom(P1a): shulaim:=max( sed1(mone,q,1),sed1(mekh,q,1)): P1:=subs(q=1+z,P1): P1:=taylor(P1,z=0,R+shulaim+1): [seq(normal(coeff(P1,z,i)),i=1..R)]: end: #CheckRec(P,q,n,R): checks the recurrence for the factorial #moments. For example, try: #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #CheckRec(p*q+(1-p),q,n,10); #Checkec((1-q^(n+1))/(1-q),q,n,4); #CheckRec((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6); CheckRec:=proc(P,q,n,R) local gu,lu,r,s,ka,fu: gu:=FMs(P,q,n,R): lu:=MekT(P,q,n,R): fu:=[]: for r from 1 to R do ka:=gu[r]-subs(n=n-1,gu[r])- add(r!/(r-s)!*lu[s]*subs(n=n-1,gu[r-s]),s=2..r-2)-r!*lu[r]: ka:=expand(ka): fu:=[op(fu),ka]: od: evalb(fu=[0$nops(fu)]): end: #FMs(P,q,n,R): The list of the first R factorial moments #as polynomials expressions in the variable n #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #Here we use the recurrence #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #FMs(p*q+(1-p),q,n,10); #FMs((1-q^(n+1))/(1-q),q,n,4); #FMs((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6); FMs:=proc(P,q,n,R) local gu,lu,gu1,r,s,i: option remember: gu:=[0]: lu:=MekT(P,q,n,R): for r from 2 to R do gu1:=expand(add(r!/(r-s)!*lu[s]*subs(n=n-1,gu[r-s]) ,s=2..r-1)+r!*lu[r]): gu1:=expand(sum(subs(n=i,gu1),i=1..n)): gu:=[op(gu),gu1]: od: gu: end: #Stir2Ansatz(r,i): The polynomial in r that is Stirling2(r,r-i) #For example try: Stir2(r,2): Stir2Ansatz:=proc(r,i) local gu,j: gu:=[seq(stirling2(j,j-i),j=1..3*i+2)]: factor(GuessPol(gu,r)): end: #Ms(P,q,n,R): The list of the first R moments #as polynomials expressions in the variable n #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #Ms(p*q+(1-p),q,n,10); #Ms((1-q^(n+1))/(1-q),q,n,4); #Ms((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6); Ms:=proc(P,q,n,R) local gu,r1,r,i1: gu:=FMs(P,q,n,R): [seq(expand(add(subs(r=r1,Stir2(r,i1))*gu[r1-i1],i1=0..r1-1)+ subs(r=r1,Stir2(r,r1))),r1=1..R)]: end: #Stir2(r,i): The polynomial in r that is Stirling2(r,r-i) #Here r is stynbolic and i is a (numeric) non-neg. integer #For example try: Stir2(r,2): Stir2:=proc(r,i) local gu,r1: option remember: if i=0 then RETURN(1): fi: gu:=(r-i)*subs(r=r-1,Stir2(r,i-1)): expand(sum(subs(r=r1,gu),r1=1..r)): end: #NorFMsE(P,q,n,R): The list of the first Normalized R Even factorial moments #as polynomials expressions in the variable n #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #The normalization is obtained by dividing by #the r-th power of the variance divided by (2r)!/r!/2^r . #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #NorFMsE(p*q+(1-p),q,n,10); #NorFMsE((1-q^(n+1))/(1-q),q,n,4); #NorFMsE((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6); NorFMsE:=proc(P,q,n,R) local gu,i: gu:=FMs(P,q,n,2*R): gu:=normal([seq(gu[2*i]/gu[2]^i/(2*i)!*i!*2^i,i=1..R)]): end: #NorFMsO(P,q,n,R): The list of the first Normalized R odd factorial moments #starting with the first moment (=0) #as polynomials expressions in the variable n #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #the r-th power of the variance divided by (2r)!/r!/2^r . #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #NorFMsO(p*q+(1-p),q,n,10); #NorFMsO((1-q^(n+1))/(1-q),q,n,4); #NorFMsO((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6); NorFMsO:=proc(P,q,n,R) local gu,i: gu:=FMs(P,q,n,2*R+1): gu:=normal([seq(gu[2*i+1]/gu[2]^i/(2*i)!*i!*2^i,i=1..R)]): end: #NorMsE(P,q,n,R): The list of the first Normalized R Even moments #as polynomials expressions in the variable n #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #The normalization is obtained by dividing by #the r-th power of the variance divided by (2r)!/r!/2^r . #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #NorMsE(p*q+(1-p),q,n,10); #NorMsE((1-q^(n+1))/(1-q),q,n,4); #NorMsE((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6); NorMsE:=proc(P,q,n,R) local gu,i: gu:=Ms(P,q,n,2*R): gu:=normal([seq(gu[2*i]/gu[2]^i/(2*i)!*i!*2^i,i=1..R)]): end: #NorMsO(P,q,n,R): The list of the first Normalized R odd moments #starting with the first moment (=0) #as polynomials expressions in the variable n #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #the r-th power of the variance divided by (2r)!/r!/2^r . #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #NorMsO(p*q+(1-p),q,n,10); #NorMsO((1-q^(n+1))/(1-q),q,n,4); #NorMsO((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,6); NorMsO:=proc(P,q,n,R) local gu,i: gu:=Ms(P,q,n,2*R+1): gu:=normal([seq(gu[2*i+1]/gu[2]^i/(2*i)!*i!*2^i,i=1..R)]): end: #AsyNorFMsE(P,q,n,r,s): The first s+1 terms #of the conjectured asymptotic expansion, in n, for the #Normalized 2r-th factorial moment #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #The normalization is obtained by dividing by #the r-th power of the variance divided by (2r)!/r!/2^r . #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #AsyNorFMsE(p*q+(1-p),q,n,r,2); #AsyNorFMsE((1-q^(n+1))/(1-q),q,n,r,3); #AsyNorFMsE((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3); AsyNorFMsE:=proc(P,q,n,r,s) local gu,i,z,mu,lu,p1,i1: gu:=NorFMsE(P,q,n,3*s+4): gu:=normal(subs(n=1/z,gu)): gu:=[seq(taylor(gu[i],z=0,s+2),i=1..nops(gu))]: mu:=1: for i from 1 to s do lu:=[seq(coeff(gu[i1],z,i),i1=1..nops(gu))]: p1:=GuessPol(lu,r): if p1=FAIL then RETURN(FAIL, mu): else mu:=mu+factor(p1)/n^i: fi: od: mu: end: #AsyNorFMsO(P,q,n,r,s): The first s+1 terms #of the conjectured asymptotic expansion, in n, for the #Normalized (2r+1)-th factorial moment #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #The normalization is obtained by dividing by #the r-th power of the variance divided by (2r)!/r!/2^r . #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #AsyNorFMsO(p*q+(1-p),q,n,r,2); #AsyNorFMsO((1-q^(n+1))/(1-q),q,n,r,3); #AsyNorFMsO((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3); AsyNorFMsO:=proc(P,q,n,r,s) local gu,i,z,mu,lu,p1,i1: gu:=NorFMsO(P,q,n,3*s+4): gu:=normal(subs(n=1/z,gu)): gu:=[seq(taylor(gu[i],z=0,s+2),i=1..nops(gu))]: mu:=0: for i from 0 to s do lu:=[seq(coeff(gu[i1],z,i),i1=1..nops(gu))]: p1:=GuessPol(lu,r): if p1=FAIL then RETURN(FAIL, mu): else mu:=mu+factor(p1)/n^i: fi: od: mu: end: #AsyFMsE(P,q,n,r,s): The first s+1 terms #of the conjectured asymptotic expansion, in n, for the #2r-th factorial moment #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #The normalization is obtained by dividing by #the r-th power of the variance divided by (2r)!/r!/2^r . #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #AsyFMsE(p*q+(1-p),q,n,r,2); #AsyFMsE((1-q^(n+1))/(1-q),q,n,r,3); #AsyFMsE((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3); AsyFMsE:=proc(P,q,n,r,s) local gu,lu: gu:=AsyNorFMsE(P,q,n,r,s): lu:=FMs(P,q,n,2)[2]: lu^r*(2*r)!/r!/2^r*gu: end: #AsyFMsO(P,q,n,r,s): The first s+1 terms #of the conjectured asymptotic expansion, in n, for the #2r+1-th factorial moment #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #The normalization is obtained by dividing by #the r-th power of the variance divided by (2r)!/r!/2^r . #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #AsyFMsO(p*q+(1-p),q,n,r,2); #AsyFMsO((1-q^(n+1))/(1-q),q,n,r,3); #AsyFMsO((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3); AsyFMsO:=proc(P,q,n,r,s) local gu,lu: gu:=AsyNorFMsO(P,q,n,r,s): lu:=FMs(P,q,n,2)[2]: lu^r*(2*r)!/r!/2^r*gu: end: #AsyMsE(P,q,n,r,s): The first s+1 terms #of the conjectured asymptotic expansion, in n, for the #2r-th moment #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #AsyMsE(p*q+(1-p),q,n,r,2); #AsyMsE((1-q^(n+1))/(1-q),q,n,r,3); #AsyMsE((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3); AsyMsE:=proc(P,q,n,r,s) local gu,lu: gu:=AsyNorMsE(P,q,n,r,s): lu:=Ms(P,q,n,2)[2]: lu^r*(2*r)!/r!/2^r*gu: end: #AsyMsE1(P,q,n,r,s): The first s+1 terms #of the conjectured asymptotic expansion, in n, for the #2r-th moment #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #AsyMsE1(p*q+(1-p),q,n,r,2); #AsyMsE1((1-q^(n+1))/(1-q),q,n,r,3); #AsyMsE1((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3); AsyMsE1:=proc(P,q,n,r,s) local gu,lu: gu:=AsyNorMsE(P,q,n,r,s): lu:=Ms(P,q,n,2)[2]: [lu,(2*r)!/r!/2^r,gu]: end: #AsyMsO(P,q,n,r,s): The first s+1 terms #of the conjectured asymptotic expansion, in n, for the #2r+1-th moment #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #AsyMsO(p*q+(1-p),q,n,r,2); #AsyMsO((1-q^(n+1))/(1-q),q,n,r,3); #AsyMsO((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3); AsyMsO:=proc(P,q,n,r,s) local gu,lu: gu:=AsyNorMsO(P,q,n,r,s): lu:=Ms(P,q,n,2)[2]: lu^r*(2*r)!/r!/2^r*gu: end: #AsyMsO1(P,q,n,r,s): The first s+1 terms #of the conjectured asymptotic expansion, in n, for the #2r+1-th moment #of the (centralized and normalized) probability distribution #induced by the rational function P(q,q^n) #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #AsyMsO1(p*q+(1-p),q,n,r,2); #AsyMsO1((1-q^(n+1))/(1-q),q,n,r,3); #AsyMsO1((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,3); AsyMsO1:=proc(P,q,n,r,s) local gu,lu: gu:=AsyNorMsO(P,q,n,r,s): lu:=Ms(P,q,n,2)[2]: [lu,(2*r)!/r!/2^r,gu]: end: #########################Proving part ###for factorial moments #CheckAsyFMsEempir(P,q,n,s,K): #Emprically yet rigorously proves the #conjectured asymptotic expression given by AsyFMsE(P,q,n,r,s) #by checking the first K terms #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #CheckAsyFMsEempir(p*q+(1-p),q,n,2,20); #CheckAsyFMsEempir((1-q^(n+1))/(1-q),q,n,3,20); #CheckAsyFMsEempir((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,3,20); CheckAsyFMsEempir:=proc(P,q,n,s,K) local gu,lu,r1,i1,r,i,y: gu:=AsyNorFMsE(P,q,n,r,s,K): lu:=NorFMsE(P,q,n,K): lu:=normal(subs(n=1/y,lu)): lu:=[seq(taylor(lu[i],y=0,s+2),i=1..nops(lu))]: evalb({0}= normal( {seq(expand(add(coeff(lu[r1],y,i1)/n^i1,i1=0..s)-subs(r=r1,gu)),r1=1..K)} )): end: #CheckAsyFMsOempir(P,q,n,s,K): #Emprically yet rigorously proves the #conjectured asymptotic expression given by AsyFMsO(P,q,n,r,s) #by checking the first K terms #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #CheckAsyFMsOempir(p*q+(1-p),q,n,2,20); #CheckAsyFMsOempir((1-q^(n+1))/(1-q),q,n,3,20); #CheckAsyFMsOempir((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,3,20); CheckAsyFMsOempir:=proc(P,q,n,s,K) local gu,lu,r1,i1,r,i,y: gu:=AsyNorFMsO(P,q,n,r,s,K): lu:=NorFMsO(P,q,n,K): lu:=normal(subs(n=1/y,lu)): lu:=[seq(taylor(lu[i],y=0,s+2),i=1..nops(lu))]: evalb({0}= normal( {seq(expand(add(coeff(lu[r1],y,i1)/n^i1,i1=0..s)-subs(r=r1,gu)),r1=1..K)} )): end: CheckAsyFMsEmpir:=proc(P,q,n,s,K): CheckAsyFMsEempir(P,q,n,s,K) and CheckAsyFMsOempir(P,q,n,s,K): end: ###end checking part for Factorial moments ###begin checking part for actual moments #CheckAsyMsEempir(P,q,n,s,K): #Emprically yet rigorously proves the #conjectured asymptotic expression given by AsyMsE(P,q,n,r,s) #by checking the first K terms #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #CheckAsyMsEempir(p*q+(1-p),q,n,2,20); #CheckAsyMsEempir((1-q^(n+1))/(1-q),q,n,3,20); #CheckAsyMsEempir((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,3,20); CheckAsyMsEempir:=proc(P,q,n,s,K) local gu,lu,r1,i1,r,i,y: gu:=AsyNorMsE(P,q,n,r,s,K): lu:=NorMsE(P,q,n,K): lu:=normal(subs(n=1/y,lu)): lu:=[seq(taylor(lu[i],y=0,s+2),i=1..nops(lu))]: evalb({0}= normal( {seq(expand(add(coeff(lu[r1],y,i1)/n^i1,i1=0..s)-subs(r=r1,gu)),r1=1..K)} )): end: #CheckAsyMsOempir(P,q,n,s,K): #Emprically yet rigorously proves the #conjectured asymptotic expression given by AsyMsO(P,q,n,r,s) #by checking the first K terms #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #CheckAsyMsOempir(p*q+(1-p),q,n,2,20); #CheckAsyMsOempir((1-q^(n+1))/(1-q),q,n,3,20); #CheckAsyMsOempir((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,3,20); CheckAsyMsOempir:=proc(P,q,n,s,K) local gu,lu,r1,i1,r,i,y: gu:=AsyNorMsO(P,q,n,r,s): lu:=NorMsO(P,q,n,K): lu:=normal(subs(n=1/y,lu)): lu:=[seq(taylor(lu[i],y=0,s+2),i=1..nops(lu))]: evalb({0}= normal( {seq(expand(add(coeff(lu[r1],y,i1)/n^i1,i1=0..s)-subs(r=r1,gu)),r1=1..K)})): end: CheckAsyMsEmpir:=proc(P,q,n,s,K): CheckAsyMsEempir(P,q,n,s,K) and CheckAsyMsOempir(P,q,n,s,K): end: ###end checking part for genuine moments ####end checking part #AsyNorMsE(P,q,n,r,s): The first s+1 terms in the #asymptotic expansion of the normalized (2r)-th moment #of the r.v. induced by P(q,q^n). #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #AsyNorMsE(p*q+(1-p),q,n,r,2); #AsyNorMsE((1-q^(n+1))/(1-q),q,n,r,2); #AsyNorMsE((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,2); AsyNorMsE:=proc(P,q,n,r,s) local gu,i,lu,i1,lu1,lu2,y,deg,vu: lu:=Ms(P,q,n,2)[2]: deg:=degree(lu,n): lu1:=normal(lu/n^deg): lu2:=normal(subs(n=1/y,1/lu1)): lu2:=taylor(lu2,y=0,s+2): lu:=add(coeff(lu2,y,i1)/n^i1,i1=0..s): gu:=0: vu:=AsyNorFMsE(P,q,n,r,s): for i from 0 to s do gu:=gu+ subs(r=2*r,Stir2(r,2*i))*subs(r=r-i,vu)* simplify((2*r-2*i)!/(2*r)!)*simplify(r!/(r-i)!)*2^i*lu^i/n^(i*deg): od: vu:=AsyNorFMsO(P,q,n,r,s): for i from 0 to s-1 do gu:=gu+subs(r=2*r,Stir2(r,2*i+1))*subs(r=r-i-1,vu) *simplify((2*r-2*i-2)!/(2*r)!)*simplify(r!/(r-i-1)!)*2^(i+1)*lu^(i+1)/ n^((i+1)*deg): : od: add(factor(coeff(gu,n,-i))/n^i,i=0..s): end: #AsyNorMsO(P,q,n,r,s): The first s+1 terms in the #asymptotic expansion of the normalized (2r+1)-th moment #of the r.v. induced by P(q,q^n). #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #AsyNorMsO(p*q+(1-p),q,n,r,2); #AsyNorMsO((1-q^(n+1))/(1-q),q,n,r,2); #AsyNorMsO((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,r,2); AsyNorMsO:=proc(P,q,n,r,s) local gu,i,lu,vu,deg,lu1,lu2,y,i1: lu:=Ms(P,q,n,2)[2]: deg:=degree(lu,n): lu1:=normal(lu/n^deg): lu2:=normal(subs(n=1/y,1/lu1)): lu2:=taylor(lu2,y=0,s+2): lu:=add(coeff(lu2,y,i1)/n^i1,i1=0..s): gu:=0: vu:=AsyNorFMsO(P,q,n,r,s): for i from 0 to s do gu:=gu+subs(r=2*r+1,Stir2(r,2*i))* subs(r=r-i,vu) *simplify((2*r-2*i)!/(2*r)!)*simplify(r!/(r-i)!)*2^i*lu^i/n^(i*deg): : gu:=expand(gu): od: vu:=AsyNorFMsE(P,q,n,r,s): for i from 0 to s do gu:=gu+subs(r=2*r+1,Stir2(r,2*i+1))*subs(r=r-i,vu) *simplify((2*r-2*i)!/(2*r)!)*simplify(r!/(r-i)!)*2^i*lu^i/n^(i*deg): gu:=expand(gu): od: add(factor(coeff(gu,n,-i))/n^i,i=0..s): end: #PPk(k,q,n): The P(q^n,q) corresponding to the generating function #of plane partitions whose 3D Ferrers diagram is inside an n by n by k box #Here k is a pos. integer, and q and n are symbolic. #For example, try: #PPk(2,q,n) PPk:=proc(k,q,n) local i,x,gu: gu:=mul((1-x^2*q^(i+1))*(1-x^2*q^(i+2)),i=0..k-1)/ mul(1-x*q^(i),i=1..k)^2: gu:=normal(gu): subs(x=q^n,gu): end: #Sipur(n,r,s): A verbose description of #aymptotic expansions for the general (2r and 2r+1) #moment (about the mean) of #some important combinatorial distributions #proving, in particular, that they asymptotically normal #Do: #Sipur(n,r,2); Sipur:=proc(n,r,s) local P,p,q,lu1,lu2,i,KA: print(`Some Interesting Probability Limit Theorems`): print(`by Shalosh B. Ekhad`): print(`Theorem 1: Consider the Binomial distribution B(n,p)`): print(`that is the probability distribution of the r.v.`): print(`Number of Heads upon Tossing a loaded coin (Pr(H)=p)`): print(` n times`): print(`The even, 2r-th, moment, about the mean is (asymptotically, ignoring`): print(1/n^(s+1), `and beyond): `): lu1:=AsyMsE1(p*q+(1-p),q,n,r,s); print(lu1[1]^r*lu1[2]*lu1[3]): print(`The variance (i.e. r=1 case) is`, lu1[1]): print(`After dividing by it to the power r (i.e. normalizing it)`): print(`we get that the normalized 2r-th moment is`): print(lu1[2]*lu1[3]): print(`The odd, (2r+1)-th, moment, about the mean is`): lu2:=AsyMsO1(p*q+(1-p),q,n,r,s); print(lu2[1]^r*lu2[2]*lu2[3]): print(`The variance (again) is`, lu1[1]): print(`After dividing by that variance to the power (r+1/2)`): print(` (i.e. normalizing it)`): print(`we get that the normalized 2r+1-th moment is`): print(lu2[2]*lu2[3]/lu2[1]^(1/2)): print(`that obviously`): print(`goes to 0`): print(`In particular, it follows that the normalized`): print(`moments of the Binomial Distribution tend`): print(`to the famous moments`, (2*r)!/r!/2^r, `(for the even (2r)) and `): print(`0, for the odd ones) of the Gaussian Distribution`): print(exp(-x^2/2)): print(`and hence is asymptotically normal.`): print(`But we proved much more: we found not just the leading term`): print(`of the asymptotics, but its expansion up to term`, 1/n^s ,` . `): print(``): print(`Theorem 2: Consider the Mahonian Distribution, whose p.g.f. is`): print(Product(1-q^i,i=1..n)/(1-q)^n/n!): print(`This p.g.f. describes both the prob. that a random permutation`): print(`of size n will have a certain number of inversions, or a certain`): print(`major index (in addition to some other more obscure permutation `): print(` statistics). The average is, of course`, n(n-1)/4, `.Centralizing`): print(`we get that the 2r-th moment is`): P:=(1-q^(n+1))/(1-q); lu1:=AsyMsE1(P,q,n,r,s); print(lu1[1]^r*lu1[2]*lu1[3]): print(`After dividing by the variance to the power r`): print(lu1[1]^r ,`(i.e. normalizing it)`): print(`we get that the normalized 2r-th moment is`): print(lu1[2]*lu1[3]): print(`In particular, it follows that the Mahonian Distribution `): print(`on permutations is asymptotically normal. A well-knwon`): print(`fact, that can be found, e.g., in : Feller v. I, 3rd ed., p. 257`): print(`But he only talked abiut the leading term, while we have established`): print(`the asymptotic to`, s+1, `terms .`): print(`Of course, by symmetry, all the odd moments are zero`): print(`but we can derive this directly, getting that indeed`): print(`the odd moments eqaul`): print(AsyMsO(P,q,n,r,s)): print(``): print(`Theorem 3: Consider the q-Catalan Distribution whose p.g.f. is`): KA:=(1-q)/(1-q^(n+1))/((2*n)!/n!/(n+1)!): print(KA*Product((1-q^(n+i))/(1-q^i),i=1..n)): print(`This p.g.f. describes both the distribution of`): print(`so called Catalan words according to the major index.`): print(` statistics). The average is, of course`, n(n-1)/2, `.Centralizing`): print(`we get that the 2r-th moment is`): P:=(1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)); lu1:=AsyMsE1(P,q,n,r,s); print(lu1[1]^r*lu1[2]*lu1[3]): print(`The variance is`, lu1[1]): print(`After dividing by the variance to the power r (i.e. normalizing it)`): print(`we get that the normalized 2r-th moment is`): print(lu1[2]*lu1[3]): print(`In particular, it follows that the q-Catalan Distribution `): print(`is asymptotically normal. As recently proved by`): print(`W.Y.C. Chen, C.J. Wang, and L. X.W. Wang`): print(`in their article: "The Limiting Distribution of the Coefficients`): print(`of the q-Catalan Numbers .`): print(`But their result only implies the leading term, of the asymptotics`): print(`and we have established`): print(`the asymptotic to`, s+1, `terms .`): print(`Of course, by symmetry, all the odd moments are zero`): print(`but we can derive this directly, getting that indeed`): print(`all the odd moments eqaul`): print(AsyMsO(P,q,n,r,s)): end: #SipurPP(n,r,K,s): A verbose description of #aymptotic expansions, to s+1 terms, of the general #moment (about the mean) of the #discrete probability distributions #for plane partitions bounded in an n by n by k #box for fixed k (k=1..K), #Do: #SipurPP(n,r,3,2); SipurPP:=proc(n,r,K,s) local P,q,lu1,k: print(`Some Interesting Probability Limit Theorems About Plane Partitions`): print(`of Bounded Height`): print(`by Shalosh B. Ekhad`): for k from 1 to K do print(`Theorem Number `, k,` : Consider the set of plane`): print(`partitions whose 3D Ferrers diagram is inside`): print(`an n by n by `, k, ` box `): print(`Let the r.v. be the: the number of cells`): print(`The even, 2r-th, moment, about the mean is (asymptotically, ignoring`): print(1/n^(s+1), `and beyond): `): P:=PPk(k,q,n): lu1:=AsyMsE1(P,q,n,r,s); print(lu1[1]^r*lu1[2]*lu1[3]): print(`The variance is`, lu1[1]): print(`After dividing by the variance to the power r (i.e. normalizing it)`): print(`we get that the normalized 2r-th moment is`): print(lu1[2]*lu1[3]): print(`In particular, it follows that this discrete prob. dist.`): print(`is asymptotically normal, but it does much more!`): print(`We have established`): print(`the asymptotic to`, s+1, `terms .`): print(`Of course, by symmetry, all the odd moments are zero`): print(`but we can derive this directly, getting that the odd moments eqaul`): print(AsyMsO(P,q,n,r,s)): od: end: #CheckMs(P,q,n,R,S): Checks Ms(P,q,n,R,S) against M1 #For example, try: #(for tossing a coin, the Mahonian statistics, and q-Catalan resp.): #CheckMs(p*q+(1-p),q,n,10,20); #CheckMs((1-q^(n+1))/(1-q),q,n,10,20); #CheckMs((1-q^(2*n+1))*(1+q^(n+1)) /(1-q^(n+1)),q,n,10,20); CheckMs:=proc(P,q,n,R,S) local gu,n1: gu:=Ms(P,q,n,R): evalb({seq(evalb(expand(subs(n=n1,gu))=expand(M1(P,q,n,n1,R))), n1=1..S)}={true}): end: #SipurFairDie(n,r,K,s): A verbose description of #aymptotic expansions, to s+1 terms, of the general #moment (about the mean) of the #discrete probability distributions #obtained by rolling a fair k-sided die for k from 2 to K #For example, try #SipurFairDie(n,r,6,2); SipurFairDie:=proc(n,r,K,s) local P,q,lu1,k,i: print(`Generalization of the Central Limit Theorem for`): print(`Repeated Throws of a fair k-faced die, for k between 2(coin)`): print(` and `, K): print(``): print(`by Shalosh B. Ekhad`): print(``): for k from 2 to K do print(`Theorem Number `, k,` : Consider the r.v.`): print(`total number of dots`): print(`obtained on rolling n times a fair die, with`, k, ` faces `): print(`such that the faces have 1, 2, ..., `, k, `dots `): print(`The even, 2r-th, moment, about the mean is (asymptotically, ignoring`): print(1/n^(s+1), `and beyond): `): P:=add(q^i,i=1..k)/k: lu1:=AsyMsE1(P,q,n,r,s); print(lu1[1]^r*lu1[2]*lu1[3]): print(`The variance is`, lu1[1]): print(`After dividing by the variance to the power r (i.e. normalizing it)`): print(`we get that the normalized 2r-th moment is`): print(lu1[2]*lu1[3]): print(`In particular, it follows that this discrete prob. dist.`): print(`is asymptotically normal, but it does much more!`): print(`We have established`): print(`the asymptotic to`, s+1, `terms .`): print(`Of course, by symmetry, all the odd moments are zero`): print(`but we can derive this directly, getting that the odd moments eqaul`): print(AsyMsO(P,q,n,r,s)): od: end: