###################################################################### ##TwoCrvs: Save this file as TwoCrvs # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read TwoCrvs # ##Then follow the instructions given there # ## # ##Written by Brian Nakamura and Doron Zeilberger, Rutgers University # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: June 4, 2013 print(`Created: June 4. 2013`): print(` This is TwoCrvs `): print(`It is a short package that accompanies the article `): print(` On the Asymptotic Statistics of The Number of Occurrences Of Multiple Permutation Pattterns `): print(`by Svante Janson, Brian Nakamura and Doron Zeilberger `): print(`and also available from Nakamura and 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: GP, GP1, Info1, Info1V, Moms1, Moms1One `): print(``): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: Alpha, AlphaSeq1 `): print(` BiNorm, BiNormD, BN, `): print(`EDist1, EDist2, MaxError1, Moms, MomsOne, MomSeq1,`): print(` MomsOneS , Sipur1, Sipur2, SymbMom `): print(` `): elif nops([args])=1 and op(1,[args])=Alpha then print(`Alpha(G,n): Given a list of lists such that`): print(`G[i+1][j+1] is the mixed (i,j) moment (about the mean)`): print(`of a certain pair of random variables, as polynomials in n,`): print(`outputs the asymptotic mixed-alpha coefficients`): elif nops([args])=1 and op(1,[args])=AlphaSeq1 then print(`AlphaSeq1(L,x,n,K): inputs a list of probability generating `): print(`functions L in the variable x, parametrized by`): print(`a discrete variable n, starting with n=1,`): print(`and returns`): print(` the list whose first two entries are the expectaction`): print(`and variance, followed by the asymptotic expansion`): print(`for as many alpha coefficients as possible. `): print(`But for the odd ones it gives the square `): print(`Try:`): print(`AlphaSeq1([seq((1+x)^i,i=1..10)],x,n,2);`): elif nops([args])=1 and op(1,[args])=BiNorm then print(`BiNorm(r,N1,N2): the list of of size N1+1 of lists of size N2+1`): print(`such that the (i+1,j+1) entry is the`): print(`the mixed (i,j) moment of the bivariate `): print(`normal distribution with correlation r,`): print(`Try: BiNorm(1/2,5,5);`): elif nops([args])=1 and op(1,[args])=BiNormD then print(`BiNormD(r,x,y): the Bi-Normal distribution with correllation r`): print(`Try: BiNormD(1/2,x,y);`): elif nops([args])=1 and op(1,[args])=BN then print(`BN(L,r,s,n,N1,N2): given a list L of combinatorial generating`): print(`functions L in the variables r and s corresponding to`): print(`two combinatorial random variables, returns the guessed`): print(`averages (expectations) followed by the `): print(`list of lists whose (i+1,j+1) entry is the mixed (i,j)`): print(`moment about the mean, as well as the asymptotic`): print(`list of corresponding alpha coefficients, and whether`): print(`it is asymptotically normal (with the found correlation,`): print(`the (2,2) entry of the alpha matrix)`): print(`try: BN([seq((1+r+s)^i,i=0..25)],r,s,n,2,2);`): elif nops([args])=1 and op(1,[args])=Sipur2 then print(`Sipur2(L,r,s,n,N1,N2,SHEM1,SHEM2): verbose version of`): print(` BN(L,r,s,n,N1,N2) `): print(`try: Sipur2([seq((1+r+s)^i,i=0..25)],r,s,n,2,2,NumUpSteps,NumRightSteps);`): elif nops([args])=1 and op(1,[args])=EDist1 then print(`EDist1(f,x,K): Given a prob. gen. function f(x) that is a `): print(`polynomial in x`): print(`for a discrete r.vs`): print(`gives the empirical dist. of`): print(` its normalized version (X-mu1)/sig1, between `): print(` mu1-K*sig and `): print(`mu1+K*sig.`): print(` For example, try: `): print(` EDist1((1+x)^100,x,2); `): elif nops([args])=1 and op(1,[args])=EDist2 then print(`EDist2(f,x,y,K): Given a prob. gen. function f(x,y) that is a `): print(`polynomial`): print(`for a pair of discrete r.vs`): print(`gives the empirical dist. of`): print(` its normalized version (X-mu1)/sig1, (Y-sig2)/sig2 between `): print(` [mu1-K*sig,mu1-K*sig] and `): print(`[mu1+K*sig,mu1+K*sig) and`): print(`it also returns the set of ratios to the binormal dist. `): print(`with the same correlation coeff.`): print(` For example, try: `): print(` EDist2((1+x+y)^40,x,y,5); `): elif nops([args])=1 and op(1,[args])=GP then print(`GP(L,x): The unique polynomial P(x)`): print(`such that P(i)=L[i+1] for i from 0 to nops(L)-1`): print(`Try: `): print(`GP([seq(i^3,i=0..5)],n);`): elif nops([args])=1 and op(1,[args])=GP1 then print(`GP1(L,x): The unique polynomial P(x)`): print(`such that P(i)=L[i] for i from 1 to nops(L)-1`): print(`Try: `): print(`GP1([seq(i^3,i=1..5)],n);`): elif nops([args])=1 and op(1,[args])=GPLL then print(`GPLL(G,n): given a list of lists of lists of numbers guesses`): print(`polynomials in n for them, calling the output H,`): print(`H[i][j](n) equals G[i][j][n+1] `): elif nops([args])=1 and op(1,[args])=Info1 then print(`Info1(L,x,N): inputs a sequence of combinatorial generating functions`): print(`(believed to have a limit distribution, after it is centralized and normalized)`): print(`returns (i) the sequence of averages (ii) the sequence of variances`): print(`(iii) the sequences of alpha coefficients`): print(` Try: `): print(` Info1([seq((1+x)^i,i=1..20)],x,8); `): elif nops([args])=1 and op(1,[args])=Info1V then print(`Info1V(L,x,N,SHEM1,SHEM2): verbose version for`): print(`Info1(L,x,N), where SHEM1 is the Sample space and SHEM is the random variable`): print(`(believed to have a limit distribution, after it is centralized and normalized)`): print(`returns (i) the sequence of averages (ii) the sequence of variances`): print(`Try:`): print(`Info1V([seq((1+x)^i,i=1..20)],x,8,"Tossing a coin n times", "Number of Heads");`): elif nops([args])=1 and op(1,[args])=MaxError1 then print(`MaxError1(f,x,K): Given a prob. gen. function f(x) that is a `): print(`polynomial`): print(`for for a discrete r.v.`): print(` gives the max difference between `): print(` its normalized version (X-mu1)/sig1 `): print(` and the normal distribution `): print(`For example, try: `): print(` MaxError1((1+x)^400,x,2); `): elif nops([args])=1 and op(1,[args])=Moms then print(`Moms(L,r,s,N1,N2): inputs a list of generating function `): print(`representing weight enumerators of a certain combinatorial`): print(`set according to variables r and s, returns the`): print(`list of lists of lists whose (i+1,j+1) entry is the mixed (i,j)`): print(`moment about the mean, try:`): print(`Moms([seq((1+x+y)^i,i=0..20)],x,y,2,2);`): elif nops([args])=1 and op(1,[args])=MomsOne then print(`MomsOne(L,r,N1): inputs a list of generating function `): print(`representing weight enumerators of a certain combinatorial`): print(`set according to the variable r, returns the`): print(`list of whose i+1 entry is the i-th`): print(`moment about the mean, try:`): print(`MomsOne([seq((1+x)^i,i=0..20)],x,10);`): elif nops([args])=1 and op(1,[args])=MomsOneS then print(`MomsOneS(Li,r,n): inputs a list of generating polynomials in the variable r, finds as many moments (about the mean)`): print(`as it can, try:`): print(`MomsOneS([seq((1+x)^i,i=0..20)],x,n);`): elif nops([args])=1 and op(1,[args])=Moms1 then print(`Moms1(P,r,s,N1,N2): inputs a generating function `): print(`representing weight enumerators of a certain combinatorial`): print(`set according to variables r and s, returns the`): print(`list of lists whose (i+1,j+1) entry is the mixed (i,j)`): print(`moment about the mean, try:`): print(`Moms1((1+x+y)^30,x,y,2,2);`): elif nops([args])=1 and op(1,[args])=Moms1One then print(`Moms1One(P,r,N1): inputs a generating function `): print(`representing weight enumerators of a certain combinatorial`): print(`set according to variable r, and returns the`): print(`a list (i+1) entry is the i-th`): print(`moment about the mean, try:`): print(`Moms1One((1+x)^10,x,4);`): elif nops([args])=1 and op(1,[args])=MomSeq1 then print(`MomSeq1(L,x,n): inputs a list of probability generating `): print(`functions L in the variable x, parametrized by`): print(`a discrete variable n, starting with n=1,`): print(`and returns`): print(` (i) the list whose first entry is the expectation `): print(` followed by as many moments as possible. Try: `): print(`MomSeq1([seq((1+x)^i,i=1..10)],x,n);`): elif nops([args])=1 and op(1,[args])=Sipur1 then print(`Sipur1(L,x,n,K,SHEM): inputs a list L of`): print(`combinatorial generating functions `): print(`starting at n=1, and phrased in`): print(`terms of a variable x, a discrete variable n`): print(`a positive integer L, and a name SHEM`): print(`returns an article with`): print(`explicitly rigorously derived expressions for`): print(`as many moments that it can find, phrased in terms of n`): print(`It also lists the alpha coefficients and investigates`): print(`whether it is asymptotically normal up to the`): print(`moments that it found. Try`): print(`Sipur1([seq((1+x)^i,i=1..10)],x,n,2,"Number of Heads on Tossing a Fair Coin n times");`): elif nops([args])=1 and op(1,[args])=SymbMom then print(`SymbMom(L,x,i,n): inputs a list,L, of generating polynomials in x and`): print(`outputs, the polynomial in n that is an exact expression for`): print(`the moment about the mean if i>1, and the average if i=1, try:`): print(`SymbMom([seq((1+x)^i,i=1..20)],x,2,n);`): else print(`There is no ezra for`,args): fi: end: ####guess polynomial Di:=proc(L) local i: [ seq( L[i+1]-L[i],i=1..nops(L)-1)]: end: #GP(L,x): The unique polynomial P(x) #such that P(i)=L[i+1] for i from 0 to nops(L)-1 GP:=proc(L,x) local i,L1,P: L1:=L: P:=0: for i from 1 to nops(L) do P:=P+L1[1]*binomial(x,i-1): L1:=Di(L1): if convert(L1,set)={0} then RETURN(factor(expand(P))): fi: od: FAIL: end: #GP1(L,x): The unique polynomial P(x) #such that P(i)=L[i] for i from 1 to nops(L) GP1:=proc(L,x) local gu: gu:=GP(L,x): if gu=FAIL then RETURN(FAIL): else factor(subs(x=x-1,gu)): fi: end: ###End from GP #Moms1(P,r,s,N1,N2): inputs a generating function #representing weight enumerators of a certain combinatorial #set according to variables r and s, returns the #list of lists whose (i+1,j+1) entry is the mixed (i,j) #moment about the mean, try: #Moms1([seq((1+x+y)^i,i=0..20)],x,y,2,2); Moms1:=proc(P,r,s,N1,N2) local P1,P1r,P1s,P2,P3,gu,i,j,gu1: P1:=P/subs({r=1,s=1},P): P1r:=subs({r=1,s=1},diff(P1,r)): P1s:=subs({r=1,s=1},diff(P1,s)): P1:=P1/r^P1r/s^P1s: gu:=[]: P2:=P1: for i from 0 to N1 do gu1:=[]: P3:=P2: for j from 0 to N2 do gu1:=[op(gu1),subs({r=1,s=1},P3)]: P3:=s*diff(P3,s): od: gu:=[op(gu),gu1]: P2:=r*diff(P2,r): od: gu: end: #Aves1(P,r,s,N1,N2): inputs a generating function #representing weight enumerators of a certain combinatorial #set according to variables r and s, returns the #two averages (for r and s) #Aves1([seq((1+x+y)^i,i=0..20)],x,y,2,2); Aves1:=proc(P,r,s,N1,N2) local P1,P1r,P1s,P2,P3,gu,i,j,gu1: P1:=P/subs({r=1,s=1},P): P1r:=subs({r=1,s=1},diff(P1,r)): P1s:=subs({r=1,s=1},diff(P1,s)): [P1r,P1s]: end: #Moms(Li,r,s,N1,N2), inputs a list of bivariate g.f. in (r,s) #and outputs the #list of lists of length N1+1,N2+1 such that L[i+1][j+1] #is the (i,j)-th mixed moment Moms:=proc(Li,r,s,N1,N2) local i: [seq(Moms1(Li[i],r,s,N1,N2),i=1..nops(Li))]: end: #Aves(Li,r,s), inputs a list of bivariate g.f. in (r,s) #and outputs the two lists of averages Aves:=proc(Li,r,s) local i: [seq(Aves1(Li[i],r,s),i=1..nops(Li))]: end: #GPLL(G,n): given a list of lists of lists of numbers guesses #polynomials in n for them GPLL:=proc(G,n) local i,a,b,gu,gu1,L: gu:=[]: for a from 1 to nops(G[1]) do gu1:=[]: for b from 1 to nops(G[1][1]) do L:=[seq(G[i][a][b],i=1..nops(G))]: gu1:=[op(gu1),GP(L,n)]: od: gu:=[op(gu),gu1]: od: gu: end: #Alpha(G,n): Given a list of lists such that #G[i+1][j+1] is the mixed (i,j) moment (about the mean) #of a certain pair of random variables, as polynomials in n, #outputs the asymptotic mixed-alpha coefficients Alpha:=proc(G,n) local gu,gu1,m02,m20,NoGood,i,j: m02:=G[1][3]: m20:=G[3][1]: gu:=[]: for i from 0 to nops(G)-1 do gu1:=[]: for j from 0 to nops(G[i+1])-1 do if degree(G[i+1][j+1],n)degree(m02,n)*i/2+degree(m20,n)*j/2 then gu1:=[op(gu1),NoGood]: else gu1:=[op(gu1), simplify(lcoeff(G[i+1][j+1],n)/lcoeff(m20,n)^(i/2)/lcoeff(m02,n)^(j/2))]: fi: od: gu:=[op(gu),gu1]: od: gu: end: #BiNorm(r,N1,N2): the list of lists of size (N1+1,N2+1) #such that the (i+1,j+1) entry is the #the mixed (i,j) moment of the bivariate #normal distribution with correlation r, #Try: BiNorm(1/2,5,3); BiNorm:=proc(r,N1,N2) local F,gu,x,y,gu1,i,j: F:=1/2/Pi/sqrt(1-r^2)*exp(-1/2/(1-r^2)*(x^2+y^2-2*r*x*y)): gu:=[]: for i from 0 to N1 do gu1:=[]: for j from 0 to N2 do gu1:=[op(gu1),int(int(F*x^i*y^j,x=-infinity..infinity),y=-infinity..infinity)]: od: gu:=[op(gu),simplify(gu1)]: od: gu: end: #BiNormD(r,x,y): the Bi-Normal distribution with correllation r #Try: BiNormD(1/2,x,y); BiNormD:=proc(r,x,y) : 1/2/Pi/sqrt(1-r^2)*exp(-1/2/(1-r^2)*(x^2+y^2-2*r*x*y)): end: #OtoDavar(L1,L2): are the lists L1 and L2 the same? OtoDavar:=proc(L1,L2) local i,j: if nops(L1)<>nops(L2) then RETURN(false): fi: if [seq(nops(L1[i]),i=1..nops(L1))]<>[seq(nops(L2[i]),i=1..nops(L2))] then RETURN(false): fi: for i from 1 to nops(L1) do for j from 1 to nops(L2) do if simplify(L1[i][j]-L2[i][j])<>0 then RETURN(false): fi: od: od: true: end: #BN(L,r,s,n,N1,N2): given a list L of combinatorial generating #functions L in the variables r and s corresponding to #two combinatorial random variables, returns the guessed #list of lists whose (i+1,j+1) entry is the mixed (i,j) #moment about the mean, as well as the aymptotic #list of corresponding alpha coefficients, and whether #it is asymptotically normal (with the found correlation, #the (2,2) entry of the alpha matrix) #try: BN(L15,r,s,n,2,2); BN:=proc(L,r,s,n,N1,N2) local gu,mu,lu,rho,lu1,ka,kan1,kan2: ka:=Aves(L,r,s): kan1:=GP([seq(ka[i][1],i=1..nops(ka))],n): if kan1=FAIL then RETURN(FAIL): fi: kan2:=GP([seq(ka[i][2],i=1..nops(ka))],n): if kan2=FAIL then RETURN(FAIL): fi: gu:=Moms(L,r,s,N1,N2): mu:=GPLL(gu,n): if member(FAIL,ListTools[Flatten](mu)) then print(`[N1,N2]=`, [N1,N2], `are too big, make them smaller .`): lprint(gu,mu): RETURN(FAIL): fi: lu:=Alpha(mu,n): rho:=lu[2][2]: lu1:=BiNorm(rho,N1,N2): if OtoDavar(lu,lu1) then RETURN([[kan1,kan2],mu,lu,true]): else RETURN([[kan1,kan2],mu,lu,lu1,false]): fi: end: #Sipur2(L,r,s,n,N1,N2,SHEM1,SHEM2): verbose version of BN(L,r,s,n,N1,N2) #where SHEM1 and SHEM2 are the names of the two combinatorial #random variables #try: Sipur2([seq((1+r+s)^i,i=0..10)],r,s,n,2,2); Sipur2:=proc(L,r,s,n,N1,N2,SHEM1,SHEM2) local gu,t0,i,j: t0:=time(): gu:=BN(L,r,s,n,N1,N2): if gu=FAIL then RETURN(gu): fi: print(`Evidence for the BiVariate Asymptotic Normality`): print(`for the Pair of Combinatorial Statistics `): print(SHEM1 , `and `, SHEM2): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`A certain family of combinatorial sets has two natural`): print(`statistics defined on its objects, and the first`, nops(L) ): print(`terms of the weight-enumerators, in the variabels, r,s,`): print(`is as inputted`): print(``): print(`The expectations of the two random variables for for the n-th`): print(`family are, respectively`): print(gu[1]): print(`and in Maple input format`): print(``): lprint(gu[1]): print(``): print(`Let's investigate all the mixed moments (about the mean)`): print(`(i,j) for i between 0 and`, N1, `and for j between 0 and `, N2): print(`The variance of the two random variables are`): print(gu[2][3][1], gu[2][1][3]): print(`and in Maple input format `): lprint(gu[2][3][1], gu[2][1][3]): print(``): print(`the covariance is`): print(``): print(gu[2][2][2]): print(``): print(`and in Maple input format `): print(``): lprint(gu[2][2][2]): print(``): for i from 0 to N1 do for j from 0 to N2 do print(``): print(`The mixed`, (i,j), `moment (about the mean) is`): print(``): print(gu[2][i+1][j+1]): print(``): print(`and in Maple input format `): print(``): lprint(gu[2][i+1][j+1]): od: od: print(`The asymptotic correlation is`): print(``): print(gu[3][2][2]): print(``): print(`and in Maple input format `): print(``): lprint(gu[3][2][2]): print(`Let's see whether it coincides with the moments of a bivariate`): print(`normal distribution with that correlation, i.e. `, gu[3][2][2]): if nops(gu)=4 and gu[4] then print(`Yes they agree!, there are the same, namely`): print(``): print(gu[3]): print(``): print(`and in Maple input format `): print(``): lprint(gu[3]): else print(`What a shame!, they don't. The asymptotic alpha matrix is`): print(``): lprint(gu[3]): print(``): print(`while it should be`): print(``): lprint(gu[4]): fi: print(`This ends this fascinating article, that took`, time()-t0, `seconds. `): end: #plotDistOld(f,x,K): Given a prob. gen. function f(x) that has a #Taylor series #for a discrete r.v. #plots its normalized version (X-mu)/sig between mu-K1*sig and #mu+K2*sig #For example, try: #plotDist((1+x)^40,x,5,5); plotDistOld:=proc(f,x,K1,K2) local mu,f1,lu,gu,sig,i,j1: f1:=f/subs(x=1,f): mu:=subs(x=1,diff(f1,x)): gu:=f1/x^mu: sig:=sqrt(subs(x=1,x*diff(x*diff(gu,x),x))): lu:=taylor(f1,x=0,trunc(mu)+K2*trunc(sig)+10): lu:=[seq([i,coeff(lu,x,i)],i=max(0,trunc(mu-K1*sig))..trunc(mu+K2*sig))]: lu:=evalf([seq([(lu[j1][1]-mu)/sig,lu[j1][2]*sig],j1=1..nops(lu))]): plot(lu,scaling=constrained): end: #plotDist(f,x,K): Given a polynomial prob. gen. function f(x) that has a #for a discrete r.v. #plots its normalized version (X-mu)/sig between mu-K1*sig and #mu+K2*sig #For example, try: #plotDist((1+x)^40,x,5,5); plotDist:=proc(f,x,K) local mu,f1,lu,gu,sig,i: f1:=f/subs(x=1,f): mu:=subs(x=1,diff(f1,x)): gu:=f1/x^mu: sig:=sqrt(subs(x=1,x*diff(x*diff(gu,x),x))): lu:=[]: for i from trunc(mu)-K*trunc(sig) to trunc(mu)+K*trunc(sig) do lu:=[op(lu), [(i-mu)/sig,coeff(f1,x,i)*sig]]: od: gu:=plot(lu): end: #EDist(f,x,K): Given a polynomial prob. gen. function f(x) that has a #Taylor series #for a discrete r.v. #plots its normalized version (X-mu)/sig between mu-K1*sig and #mu+K2*sig #For example, try: #plotDist((1+x)^40,x,5,5); EDist:=proc(f,x,K) local mu,f1,lu,gu,sig,i: f1:=f/subs(x=1,f): mu:=subs(x=1,diff(f1,x)): gu:=f1/x^mu: sig:=sqrt(subs(x=1,x*diff(x*diff(gu,x),x))): lu:=[]: for i from trunc(mu)-K*trunc(sig) to trunc(mu)+K*trunc(sig) do lu:=[op(lu), [(i-mu)/sig,coeff(f1,x,i)*sig]]: od: evalf(lu): end: #EDist2Ori(f,x,y,K): Given a prob. gen. function f(x,y) that is a #polynomial #for a pair of discrete r.vs #gives the empirical dist. of # its normalized version (X-mu1)/sig1, (Y-sig2)/sig2 between #[mu1-K*sig,mu1-K*sig) and #[mu1+K*sig,mu1+K*sig) and #For example, try: #EDist2((1+x+y)^40,x,y,5); EDist2Ori:=proc(f,x,y,K) local mu1,mu2,f1,lu,gu,sig1,sig2, i,j,f11: f1:=f/subs({x=1,y=1},f): mu1:=evalf(subs({x=1,y=1},diff(f1,x))): mu2:=evalf(subs({x=1,y=1},diff(f1,y))): gu:=f1/x^mu1/y^mu2: sig1:=evalf(sqrt(subs({x=1,y=1},x*diff(x*diff(gu,x),x)))): sig2:=evalf(sqrt(subs({x=1,y=1},y*diff(y*diff(gu,y),y)))): lu:=[]: for i from trunc(mu1)-K*trunc(sig1) to trunc(mu1)+K*trunc(sig1) do f11:=coeff(f1,x,i): for j from trunc(mu2)-K*trunc(sig2) to trunc(mu2)+K*trunc(sig2) do if coeff(f11,y,j)<>0 then lu:=[op(lu), [[(i-mu1)/sig1, (j-mu2)/sig2], coeff(f11,y,j)*sig1*sig2]]: fi: od: od: lu: end: #MaxError2(f,x,y,K): Given a prob. generating function, a polynomial f #in the variables x,y #finds how the maximum deviation from approximating it by a #Bi-Normal Distribution #with the same correlation #Try: #MaxError2((1+x+y)^400,x,y,5); MaxError2:=proc(f,x,y,K) local mu,r,F,gu,i: mu:=Moms([f],x,y,2,2)[1]: r:=evalf(mu[2][2]/sqrt(mu[1][3]*mu[3][1])): F:=BiNormD(r,x,y) : gu:=EDist2(f,x,y,K); max(seq(abs(evalf(subs({x=gu[i][1][1],y=gu[1][1][2]},F)-gu[i][2])),i=1..nops(gu))): end: #EDist1(f,x,K): Given a prob. gen. function f(x) that is a #polynomial #for for a discrete r.v. #gives the empirical dist. of # its normalized version (X-mu1)/sig1 #[mu1-K*sig,mu1+K*sig] #For example, try: #EDist1((1+x)^400,x,2); EDist1:=proc(f,x,K) local mu1,f1,lu,gu,sig1, i,f11: f1:=f/subs({x=1},f): mu1:=evalf(subs({x=1},diff(f1,x))): gu:=f1/x^mu1: sig1:=evalf(sqrt(subs({x=1},x*diff(x*diff(gu,x),x)))): lu:=[]: for i from trunc(mu1)-K*trunc(sig1) to trunc(mu1)+K*trunc(sig1) do f11:=coeff(f1,x,i): lu:=[op(lu), [(i-mu1)/sig1, f11*sig1]]: od: lu: end: #MaxError1(f,x,K): Given a prob. gen. function f(x) that is a #polynomial #for for a discrete r.v. #gives the max difference between # its normalized version (X-mu1)/sig1 #and the normal distrubution #For example, try: #MaxError1((1+x)^400,x,2); MaxError1:=proc(f,x,K) local gu,F,i: gu:=EDist1(f,x,K): F:=exp(-x^2/2)/sqrt(2*Pi): max(seq(abs(gu[i][2]-evalf(subs(x=gu[i][1],F))), i=1..nops(gu))): end: #MaxError1R(f,x,K): Given a prob. gen. function f(x) that is a #polynomial #for for a discrete r.v. #gives the max difference between 1 and the ratios of # its normalized version (X-mu1)/sig1 #and the normal distrubution #For example, try: #MaxError1R((1+x)^400,x,2); MaxError1R:=proc(f,x,K) local gu,F,i: gu:=EDist1(f,x,K): F:=exp(-x^2/2)/sqrt(2*Pi): max(seq(abs(1-gu[i][2]/evalf(subs(x=gu[i][1],F))), i=1..nops(gu))): end: #EDist2Plus(f,x,y,K): Given a prob. gen. function f(x,y) that is a #polynomial #for a pair of discrete r.vs #gives the empirical dist. of # its normalized version (X-mu1)/sig1, (Y-sig2)/sig2 between #[mu1-K*sig,mu1-K*sig) and #[mu1+K*sig,mu1+K*sig) and #For example, try: #EDist2((1+x+y)^40,x,y,5); EDist2:=proc(f,x,y,K) local mu1,mu2,f1,lu,gu,sig1,sig2, i,j,f11,r,F,ku: f1:=f/subs({x=1,y=1},f): mu1:=evalf(subs({x=1,y=1},diff(f1,x))): mu2:=evalf(subs({x=1,y=1},diff(f1,y))): gu:=f1/x^mu1/y^mu2: sig1:=evalf(sqrt(subs({x=1,y=1},x*diff(x*diff(gu,x),x)))): sig2:=evalf(sqrt(subs({x=1,y=1},y*diff(y*diff(gu,y),y)))): r:=evalf(subs({x=1,y=1},y*diff(x*diff(gu,x),y)))/sig1/sig2: F:=BiNormD(r,x,y) : lu:=[]: ku:={}: for i from trunc(mu1)-K*trunc(sig1) to trunc(mu1)+K*trunc(sig1) do f11:=coeff(f1,x,i): for j from trunc(mu2)-K*trunc(sig2) to trunc(mu2)+K*trunc(sig2) do if coeff(f11,y,j)<>0 then lu:=[op(lu), [[(i-mu1)/sig1, (j-mu2)/sig2], evalf(coeff(f11,y,j)*sig1*sig2), evalf(subs({x=(i-mu1)/sig1, y=(j-mu2)/sig2},F))]]: ku:=ku union {evalf(coeff(f11,y,j)*sig1*sig2)/ evalf(subs({x=(i-mu1)/sig1, y=(j-mu2)/sig2},F))}: fi: od: od: lu,ku: end: ######ONE VARIABLE #Moms1One(P,r,N1): inputs a generating function #representing weight enumerators of a certain combinatorial #set according to variable r, and returns the #a list (i+1) entry is the i-th #moment about the mean, try: #Moms1One([seq((1+x)^i,i=0..20)],x,4); Moms1One:=proc(P,r,N1) local P1,P1r,P2,gu,i: P1:=P/subs({r=1},P): P1r:=subs({r=1},diff(P1,r)): P1:=P1/r^P1r: gu:=[]: P2:=P1: for i from 0 to N1 do gu:=[op(gu),subs({r=1},P2)]: P2:=r*diff(P2,r): od: gu: end: #MomsOne(Li,r,N1), inputs a list of g.f.s in r #and outputs the #list of length N1+1 such that L[i+1] #is the it-th moment MomsOne:=proc(Li,r,N1) local i: [seq(Moms1One(Li[i],r,N1),i=1..nops(Li))]: end: #SymbMom(L,x,i,n): inputs a list,L, of generating polynomials in x and #outputs, the polynomial in n that is an exact expression for #the moment about the mean if i>1, and the average if i=1, try: #SymbMom([seq((1+x)^i,i=1..20)],x,2,n); SymbMom:=proc(L,x,i,n) local L1,L11,j,i1: L1:=[seq(L[j]/subs(x=1,L[j]),j=1..nops(L))]: if i=1 then L11:=subs(x=1,diff(L1,x)): RETURN(GP1(L11,n)): fi: L1:=[seq(L1[j]/x^subs(x=1,diff(L1[j],x)),j=1..nops(L1))]: for j from 2 to i do L1:=[seq(x*diff(L1[i1],x),i1=1..nops(L1))]: od: L11:=subs(x=1,diff(L1,x)): GP1(L11,n): end: #MomSeq1(L,x,n): inputs a list of probability generating #functions L in the variable x, parametrized by #a discrete variable n, starting with n=1, #and returns #(i) the list whose first entry is the expectation #followed by as many moments as possible MomSeq1:=proc(L,x,n) local gu,lu,i: option remember: gu:=[SymbMom(L,x,1,n)]: for i from 2 do lu:=SymbMom(L,x,i,n): if lu=FAIL then RETURN(gu): else gu:=[op(gu),lu]: fi: od: end: #Asy1(R,n,K): inputs a rational function R in n #and outputs the asymptotic to order K. Try: #Asy1((n+1)/(n+2),n,3); Asy1:=proc(R,n,K) local R1,x,i: R1:=normal(subs(n=1/x,R)): R1:=taylor(R1,x=0,K+3): add(coeff(R1,x,i)/n^i,i=0..K): end: #AlphaSeq1(L,x,n,K): inputs a list of probability generating #functions L in the variable x, parametrized by #a discrete variable n, starting with n=1, #and returns # the list whose first two entries are the expectaction #and variance, followed by the asymptotic expansion #for as many alpha coefficients as possible. #But for the odd ones it gives the square.Try: #AlphaSeq1([seq((1+x)^i,i=1..10)],x,n,2); AlphaSeq1:=proc(L,x,n,K) local gu,mu,i,lu: gu:=MomSeq1(L,x,n): if nops(gu)<2 or gu[2]=0 then RETURN(gu): fi: lu:=[gu[1],gu[2]]: for i from 3 to nops(gu) do mu:=gu[i]: if mu=0 then lu:=[op(lu),mu]: else if i mod 2=0 then mu:=gu[i]/gu[2]^(i/2): else mu:=gu[i]^2/gu[2]^i: fi: mu:=Asy1(mu,n,K): lu:=[op(lu),mu]: fi: od: lu: end: #Sipur1(L,x,n,K,SHEM): inputs a list L of #combinatorial generating functions #starting at n=1, and phrased in #terms of a variable x, a discrete variable n #a positive integer L, and a name SHEM #returns an article with #explicitly rigorously derived expressions for #as many moments that it can find, phrased in terms of n #It also lists the alpha coefficients and investigates #whether it is asymptotically normal up to the #moments that it found. Try #Sipur1([seq((1+x)^i,i=1..10)],x,n,2,"Number of Heads on Tossing a Fair Coin n times"); Sipur1:=proc(L,x,n,K, SHEM) local gu,mu,t0,i: t0:=time(): gu:=MomSeq1(L,x,n): if nops(gu)<2 then RETURN(FAIL): fi: mu:=AlphaSeq1(L,x,n,K): print(`Rigorous Expressions for the first`, nops(gu), ` moments `): print(`For the Combinatorial Statistic`, SHEM): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`Let X(n) be the combinatorial random variable `): print(SHEM): print(`The expected value is`): print(gu[1]): print(`and in Maple input notation`): lprint(gu[1]): print(``): print(`The variance is`): print(gu[2]): print(``): print(`and in Maple input notation`): lprint(gu[2]): print(``): for i from 3 to nops(gu) do print(`The `, i, `-th moment about the mean is`): print(gu[i]): print(``): print(`and in Maple input notation`): lprint(gu[i]): print(``): od: if nops(mu)>2 then print(`The square of the skewness to order `, K, `is `): print(mu[3]): print(``): print(`and in Maple input notation`): lprint(mu[3]): print(``): if coeff(mu[3],n,0)=0 then print(`it indeed tends to 0, as it should `): fi: fi: if nops(mu)>3 then print(`The kurtosis to order `, K, `is `): print(mu[4]): print(``): print(`and in Maple input notation`): print(``): lprint(mu[4]): print(``): if coeff(mu[4],n,0)=3 then print(`It indeed tends to`, 3, `as it should be! `): fi: fi: for i from 5 to nops(mu) do if i mod 2=1 then print(`The square of the `, i , ` -th alpha coeff. to order `, K, `is `): print(mu[i]): print(``): print(`and in Maple input notation`): print(``): lprint(mu[i]): print(``): if coeff(mu[i],n,0)=0 then print(`It indeed tends to`, 0, `as it should be! `): fi: else print(`The `, i , ` -th alpha coeff. to order `, K, `is `): print(``): print(mu[i]): print(``): print(`and in Maple input notation`): lprint(mu[i]): print(``): if coeff(mu[i],n,0)=i!/(i/2)!/2^(i/2) then print(`It indeed tends to`, i!/(i/2)!/2^(i/2) , `as it should be! `): fi: fi: od: print(`This ends this exciting article, that took`, time()-t0, `seconds to generate`): end: #Info1(L,x,N): inputs a sequence of combinatorial generating functions #(believed to have a limit distribution, after it is centralized and normalized) #returns (i) the sequence of averages (ii) the sequence of variances #(iii) the sequences of alpha coefficients #Try: #Info1([seq((1+x)^i,i=1..20)],x,8); Info1:=proc(L,x,N) local L1,i,Av1,Var1,M,A,j: if nops(L)<6 then print(L, `is too short`): RETURN(FAIL): fi: L1:=[seq(L[i]/subs(x=1,L[i]),i=1..nops(L))]: Av1:=subs(x=1,diff(L1,x)): L1:=[seq(L1[j]/x^Av1[j],j=1..nops(L1))]: L1:=x*diff(L1,x): L1:=x*diff(L1,x): Var1:=subs(x=1,L1): for i from 3 to N do L1:=x*diff(L1,x): M[i]:=subs(x=1,L1): od: for i from 3 to N do A[i]:=evalf([seq(M[i][j]/Var1[j]^(i/2),j=nops(M[i])-5..nops(M[i]))]): od: Av1,Var1,[seq(A[i],i=3..N)]: end: #Info1V(L,x,N,SHEM1,SHEM2): verbose version for #Info1(L,x,N), where SHEM1 is the Sample space and SHEM is the random variable #(believed to have a limit distribution, after it is centralized and normalized) #returns (i) the sequence of averages (ii) the sequence of variances #Try: #Info1V([seq((1+x)^i,i=1..20)],x,8),"Tossing a coin n times", "Number of Heads"); Info1V:=proc(L,x,N,SHEM1,SHEM2) local gu,j: if nops(L)<6 then print(L, `is too short`): RETURN(FAIL): fi: gu:=Info1(L,x,N): print(`Numerics for the combinatorial Random variable`, SHEM2, `defined in the samle space of`, SHEM1): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`The input is the sequence of generating functions for the random variable`, SHEM2, `defined over the sample space`, SHEM1): print(`for n from 1 to `, nops(L)): print(``): print(`The sequence of averages is`): print(``): lprint(gu[1]): print(``): print(`and in floating-point`): lprint(evalf(gu[1])): print(``): print(`The sequence of variances is`): print(``): lprint(gu[2]): print(``): print(``): print(`and in floating-point`): lprint(evalf(gu[2])): print(``): for j from 1 to nops(gu[3]) do print(`The `, j+2, ` -th alpha coefficients for n from`, nops(L)-5, `to `, nops(L) , `are : `): print(gu[3][j]): od: print(`This ends this numerical short article.`): end: