###################################################################### ##MVPoisson: Save this file as MVPoisson # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read MVPoisson # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: Feb.-April, 2009 #Version of March 15, 2010, thanks to useful comments of Arthur Hipke! print(`Created: Feb.-Apr. 2009`): print(`This version: March 15, 2010, using helpful comments of Arthur Hipke`): print(` This is MVPoisson `): print(`It the Maple package that accompanies the article `): print(`A symbolic Computation Approach to a Problem`): print(`involving Multivariate Poisson Distrubutions`): print(`by Eduardo Sontag and Doron Zeilberger`): print(`and also available from Sontag's and Zeilberger's websites`): 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: Ave,Ave2, Ave2f, Avef,`): print(` CheckRecsG, Cor, Cor2f, Corf, Cov, Cov2f, Covf `): print(` gtoc2, Hatkhala, Opes , SeqFromRec, ValueFromRec1 `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: `): print(`AllAve,AllAve2, AllAve2f, AllAvef, C,Cf,`): print(` C2,C2f, CorM, CorMf, CorM2f`): print(` FM1, GF1,Recs, RecsV, RecsG, SipurD, SipurDf, SipurD2f,Tay `): print(`ValueFromRecs, Var, Var2, Varf, Var2f, Zinn, ZinnC `): elif nops([args])=1 and op(1,[args])=AllAve then print(`AllAve(A,g,b): All the averages of x[j] in the `): print(`for j from 1 to n:=nops(A[1]), in `): print(`Eduardo Sontag prob. dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) g[i]'s`): print(`or numericl g-s (that has to be a list of length n`): print(`Warning: it is much faster, for some reason, to have g as rationals`): print(`or integers, not as floats`): print(` For example, try`): print(`AllAve(A1,g,[2,2,2,2]);`): print(`AllAve(A1,[1$8],[2,2,2,2]);`): elif nops([args])=1 and op(1,[args])=AllAvef then print(`AllAvef(A,g,b): Like AllAve(A,g,b), but using Almkvist-Zeilberger`): print(`Warning: for large matrices it is infeasible `): print(`Warning: for large matrices it is infeasible `): print(` For example, try`): print(`AllAvef(A2,g,[4,4]);`): print(`AllAvef(A2,[1,2,3,4,5],[4,4]);`): elif nops([args])=1 and op(1,[args])=AllAve2f then print(`AllAve2f(A,g,b): Like AllAve2(A,g,b), but using the Zeilberger`): print(`recurrence implied by the single sum`): print(` For example, try`): print(`AllAve2f(A2,g,[4,4]);`): print(`AllAve2f(A2,[1,2,3,4,5],[4,4]);`): elif nops([args])=1 and op(1,[args])=AllAve2 then print(`AllAve2(A,g,b): Like AllAve(A,g,b), but only applicable`): print(`to 0-1 two-rowed matrices, using direct summation `): print(` For example, try`): print(`AllAve2(A2,g,[4,4]);`): elif nops([args])=1 and op(1,[args])=Ave then print(`Ave(A,g,b,j): The average of x[j] in the Eduardo `): print(`Sontag prob. dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) g[i]'s`): print(`or numericl g-s`): print(`(1<=i<=n).`): print(` For example, try`): print(`Ave([[1,1]],g,[2],1);`): elif nops([args])=1 and op(1,[args])=Ave2 then print(`Ave2(A,g,b,j): like Ave(A,g,b,j) (using generating functions)`): print(`but a faster version (using single-summation) only applicable`): print(` for two rows `): print(` For example, try`): print(`Ave2(A2,g,[2,2],1);`): elif nops([args])=1 and op(1,[args])=Ave2f then print(`Ave2f(A,g,b,j): like Ave2(A,g,b,j) `): print(`but a faster version (using the recurrence) only applicable`): print(` for two rows `): print(` For example, try`): print(`Ave2f(A2,g,[2,2],1);`): elif nops([args])=1 and op(1,[args])=Avef then print(`Avef(A,g,b,j): Fast version of Ave(A,g,b,j)`): print(`The average of x[j] in the Eduardo `): print(`Sontag prob. dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) g[i]'s`): print(`or numericl g-s`): print(`(1<=i<=n).`): print(` For example, try`): print(`Ave([[1,1]],g,[2],1);`): elif nops([args])=1 and op(1,[args])=C then print(`C(A,b,g): The normalization constant, computed directly`): print(`via the generating function`): print(` For example, try`): print(`C([[1,2,1],[2,1,3]],[2,4],g);`): elif nops([args])=1 and op(1,[args])=Cf then print(`Cf(A,b,g): The normalization constant, computed `): print(`using a recurrence found by the Apagodu-Zeilberger`): print(`multivariable extension of Almkvist-Zeilberger.`) : print(` For example, try`): print(`Cf([[1,2,1],[2,1,3]],[2,4],g);`): elif nops([args])=1 and op(1,[args])=C2 then print(`C2(A,b,g): The normalization constant for`): print(`the case where A has two rows AND is a 0-1 matrix`): print(` For example, try`): print(`C2([[1,1,0],[1,0,1]],[2,4],g);`): elif nops([args])=1 and op(1,[args])=C2f then print(`C2f(A,b,g): fast version of C2(A,B,g), the`): print(`The normalization constant for`): print(`the case where A has two rows AND is a 0-1 matrix`): print(` For example, try`): print(`C2f([[1,1,0],[1,0,1]],[2,4],g);`): elif nops([args])=1 and op(1,[args])=CheckRecsG then print(`CheckRecsG(F,z,L,kama): Checks RecsG(F,a,b,B) for`): print(`up to L-degree terms`): print(`For example, try:`): print(`CheckRecsG(exp(x+y),[x,y],10,20);`): elif nops([args])=1 and op(1,[args])=Cor then print(`Cor(A,g,b,i,j): The correlation `): print(`of x[i] and x[j] in the Eduardo Sontag problem of a dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) `): print(`or numeric g[i]'s`): print(`(1<=i<=n).`): print(` For example, try: `): print(`Cor([[1,1]],g,[3],1,2);`): elif nops([args])=1 and op(1,[args])=CorM then print(`CorM(A,g,b): The correlation matrix`): print(`in the Eduardo Sontag problem of a dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) `): print(`or numeric g[i]'s`): print(`(1<=i<=n).`): print(` For example, try: `): print(`CorM([[1,1]],g,[3]);`): elif nops([args])=1 and op(1,[args])=CorMf then print(`CorMf(A,g,b): The correlation matrix, using Cf`): print(`in the Eduardo Sontag problem of a dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) `): print(`or numeric g[i]'s`): print(`(1<=i<=n).`): print(` For example, try: `): print(`CorMf(A2,g,[3,4]);`): elif nops([args])=1 and op(1,[args])=CorM2f then print(`CorM2f(A,g,b): The correlation matrix`): print(`in the Eduardo Sontag problem of a dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) `): print(`or numeric g[i]'s`): print(`(1<=i<=n).`): print(` For example, try: `): print(`CorM2f(A2,g,[3,4]);`): elif nops([args])=1 and op(1,[args])=Corf then print(`Corf(A,g,b,i,j): The correlation, using Cf `): print(`of x[i] and x[j] in the Eduardo Sontag problem of a dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) `): print(`or numeric g[i]'s`): print(`(1<=i<=n).`): print(` For example, try: `): print(`Corf([[1,1]],g,[3],1,2);`): elif nops([args])=1 and op(1,[args])=Cor2f then print(`Cor2f(A,g,b,i,j): The correlation, using C2f, `): print(`of x[i] and x[j] in the Eduardo Sontag problem of a dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) `): print(`or numeric g[i]'s`): print(`(1<=i<=n).`): print(` For example, try: `): print(`Cov2f([[1,1]],g,[3],1,2);`): elif nops([args])=1 and op(1,[args])=Cov2f then print(`Cov2f(A,g,b,i,j): The covarinace, using C2f. Only applicable`): print(`to two-rowed matrices ,`): print(`of x[i] and x[j] in the Eduardo Sontag problem of a dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) `): print(`or numeric g[i]'s`): print(`(1<=i<=n).`): print(` For example, try: `): print(`Covf([[1,1]],g,[3],1,2);`): elif nops([args])=1 and op(1,[args])=Covf then print(`Covf(A,g,b,i,j): The covarinace, using the recurrences `): print(`of x[i] and x[j] in the Eduardo Sontag problem of a dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) `): print(`or numeric g[i]'s`): print(`(1<=i<=n).`): print(` For example, try: `): print(`Covf([[1,1]],g,[3],1,2);`): elif nops([args])=1 and op(1,[args])=FM1 then print(`FM1(A,g,b,r,j): The j-th factorial`): print(`moment of x[r] in the Eduardo Sontag prob. dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) `): print(`or numeric g[i]'s (1<=i<=n). `): print(` For example, try`): print(`FM1([[1,1]],g,[2],2,1);`): elif nops([args])=1 and op(1,[args])=GF1 then print(`GF1(A,g,z): Given a (numeric or symbolic) matrix`): print(`A (say m by n) and a symbol g and a symbol z`): print(`finds the multi-variable generating function`): print(`in z[1], ..., z[m] of the normalizing constant`): print(`For example, try:`): print(`GF1([[1,1]],g,z);`): elif nops([args])=1 and op(1,[args])=gtoc2 then print(`gtoc2(A,g): given a 0-1 matrix A with 2 rows and a vector`): print(`g of numbers (or symbols) of size nops(A[1]) (=nops(A[2]))`): print(`outputs the vector [c11,c10,c01].`): print(`gtoc2([[1,1,0],[1,0,1]],[2,3,5]);`): elif nops([args])=1 and op(1,[args])=Hatkhala then print(`Hatkhala(F,z,K): The list of lists of lists etc. such that`): print(`L[i1][i2][i3]... is the coeff. of z1^i1*z2^i2*z3^i3...`): print(`for 1<=i1<=K[1], 1<=i2<=K[2] `): print(`in the multivarite expanson of F(z[1],z[2],z[3], ...)`): print(`For example try:`): print(`Hatkhala(exp(z1+z2+z1*z2),[z1,z2],[4,5]);`): elif nops([args])=1 and op(1,[args])=Opes then print(`Opes(F,z,b,B): the recurrence operators annihilating `): print(`F/z[1]^(b[1]+1)/.../z[n]^(b[n]+1) `): print(`using B[i] as shifts (n=nops(b)=nops(z)=nops(B)) For example`): print(`Opes(exp(x+y),[x,y],[a,b],[A,B]);`): elif nops([args])=1 and op(1,[args])=Recs then print(`Recs(A,g,b,B): Given a (numeric or symbolic) matrix`): print(`A (say m by n) and a symbol g and a symbol B`): print(`g may also be a list of numbers of length nops(A[1]) `): print(`[Note: it is much faster to have a numeric g rather than a`): print(`symbolic g], `): print(`finds the recurrences for the appropriate F(b[1],b[2], ...)`): print(`of the paper`): print(`in terms of b[1], ..., b[m] and B[1], ..., B[m].`): print(`For example, try:`): print(`Recs(A2,g,b,B);`): print(`Recs(A2,[1,2,3,4,5],b,B);`): elif nops([args])=1 and op(1,[args])=RecsV then print(`RecsV(A,g,b): A verbose version of Recs(A,g,b,B)`): print(`For example, try:`): print(`RecsV(A2,g,b);`): print(`RecsV(A2,[1,2,3,4,5],b);`): elif nops([args])=1 and op(1,[args])=RecsG then print(`RecsG(F,z,b,B): Given a (numeric or symbolic) matrix`): print(`A (say m by n) and a symbol g and a symbol z`): print(`finds the recurrences for the appropriate C`): print(`in terms of b[1], ..., b[m] and B[1], ..., B[m]`): print(`For example, try:`): print(`RecsG(exp(x+y),[x,y],b,B);`): elif nops([args])=1 and op(1,[args])=SeqFromRec then print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`): print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`): print(`satisfied by the recurrence ope(n,N)f(n)=0`): print(`extends it to the first K values`): elif nops([args])=1 and op(1,[args])=SipurD then print(`SipurD(A,g,j,K,n): The diagonal story for the matrix A`): print(`Poisson parameters g, the variable x[j], up to`): print(`[K, ..., K] using the variable n for asymptotics`): print(`using the slow procedure C. For example, try:`): print(`SipurD(A1,[1$nops(A1[1])],1,20,n);`): elif nops([args])=1 and op(1,[args])=SipurDf then print(`SipurDf(A,g,j,K,n): `): print(`like SipurD(A,g,j,K,n) `): print(`but using the Almkvist-Zeilberger speed-up. For example, try:`): print(`SipurDf(A2,[1$nops(A2[1])],1,40,n);`): elif nops([args])=1 and op(1,[args])=SipurD2f then print(`SipurD2f(A,g,j,K,n): `): print(`like SipurD(A,g,j,K,n) `): print(`but using the Zeilberger recurrence. Only valid for `): print(`two rowed-matrices . For example, try:`): print(`SipurD2f(A2,[1$nops(A2[1])],1,40,n);`): elif nops([args])=1 and op(1,[args])=Tay then print(`Tay(A,g,z): Given a (numeric or symbolic) matrix`): print(`A (say m by n) and a symbol g and a symbol z`): print(`finds the multi-variable Taylor exapansion up to z[1]^K...z[m]^K`): print(`of `): print(`exp [`): print(`g[1]*z[1]^A[1][1]*z[2]^A[2][1]*z[3]^A[3][1]*...*z[m]^A[m][1]+`): print(`(g[2]*z[1]^A[1][2]*z[2]^A[2][2]*z[3]^A[3][2]*...*z[m]^A[m][2]+`): print(`+....`): print(`(g[n]*z[1]^A[1][m]*z[2]^A[2][m]*z[3]^A[3][2]*...*z[m]^A[m][n])]`): print(`For example, try:`): print(`Tay([[1,1]],g,z,4);`): elif nops([args])=1 and op(1,[args])=ValueFromRec1 then print(`ValueFromRec1(ope,n,N,Ini,K): Given the first L`): print(`terms of the sequence Ini=[f(1), ..., f(L)]`): print(`satisfied by the recurrence ope(n,N)f(n)=0`): print(`finds f(K)`): print(`For example, try:`): print(`ValueFromRec1(N^2-N-1,n,N,[1,1],10);`): elif nops([args])=1 and op(1,[args])=ValueFromRecs then print(`ValueFromRecs(opes,n,N,Inis,K): given a list of operators`): print(`opes, a list of discrete variables n, a list of the corresponding`): print(`Initial Values, and a vector the same length as opes (n and N)`): print(`computes the value of the multisequence at K.`): print(`For example, try:`): print(`ValueFromRecs([N1-2,N2-3],[n1,n2],[N1,N2],[[1]],[4,5]);`): elif nops([args])=1 and op(1,[args])=Var then print(`Var(A,g,b,j): The variance of`): print(`of x[j] in the Eduardo Sontag prob. dist. given`): print(`by the matrix A and the vector b, with generic (symbolic) `): print(`or numeric g[i]'s`): print(`(1<=i<=n).`): print(` For example, try`): print(`Var(A2,g,[3,3],1);`): elif nops([args])=1 and op(1,[args])=Var2 then print(`Var2(A,g,b,j): like Var(A,g,b,j), but a faster`): print(`version for 0-1 two-rowed matrices`): print(`For example, try:`): print(`Var2(A2,g,[3,3],1);`): elif nops([args])=1 and op(1,[args])=Var2f then print(`Var2f(A,g,b,j): like Var2(A,g,b,j), but a faster`): print(`version using the recurrence`): print(`For example, try:`): print(`Var2f(A2,g,[2,2],1);`): elif nops([args])=1 and op(1,[args])=Varf then print(`Varf(A,g,b,j): like Var(A,g,b,j), but a faster`): print(`version using the recurrence given by Almkvist-Zeilberger`): print(`For example, try:`): print(`Varf(A2,g,[2,2],1);`): elif nops([args])=1 and op(1,[args])=Zinn then print(`Zinn(L), given a list of numbers L, returns theta followed`): print(`by mu such that conjecturally.`): print(`L[n] is appx. mu^n*n^theta `): elif nops([args])=1 and op(1,[args])=ZinnC then print(`ZinnC(resh,eps,V): inputs a sequence of numbers V`): print(`and tries to conjecture an approximate asymptotic expression of`): print(`form resh[V]=C*mu^V*V^theta, `): print(`with confidence eps. `): print(`For example, try:`): print(`ZinnC([seq(2*i/i,i=1..100)],.1,n);`): else print(`There is no ezra for`,args): fi: end: ###sample input A1:= [[1, 0, 0, 0, 1, 1, 0, 0], [0 , 1, 0, 0, 1, 0, 0, 1], [0, 0, 1, 0, 0, 1, 1, 0], [0, 0, 0, 1, 0, 0, 1, 1]]; #WARNING: IT COULD WELL BE QUITE TRIVIAL - all I have done is extract the matrix #from a biological network that I found in the literature and which satisfies #the needed assumptions for the theory to apply: #(for my own memory: #given in: #Gilles Gnacadja, Alex Shoshitaishvili, Michael J. Gresser, Brian Varnum, David # Balaban, Mark Durst, Chris Vezina, Yu Li #Monotonicity of interleukin-1 receptor$,1rs(Bligand binding with respect to #antagonist in the presence of decoy receptor #Journal of Theoretical Biology #Volume 244, Issue 3, 7 February 2007, Pages 478-488 #) #(2) from a paper of mine, there is also this simpler example; I prefer to put #it here so I don't forget, even if easy: A2:= [[1, 0, 0, 1, 1], [0, 1, 1, 1, 1]]: #(Madalena Chaves, Eduardo D. Sontag, Robert J. Dinerstein #Steady-states of receptor$,1rs(Bligand dynamics: a theoretical framework #Journal of Theoretical Biology #Volume 227, Issue 3, 7 April 2004, Pages 413-428) #(3) Here is another example: A3:= [[1, 0, 1, 0, 1, 1], [0, 1, 1, 1, 0, 1]]: #(computed, perhaps with rrors --to be checked!-- from: #Eric Batchelor, Mark Goulian #Robustness and the cycle of phosphorylation and dephosphorylation in a #two-component regulatory system #Proceedings of the National Academy of Sciences #Vol. 100, No. 2. (21 January 2003), pp. 691-696. #) ###end sample input #Tay(A,g,z): Given a (numeric or symbolic) matrix #A (say m by n) and a symbol g and a symbol z #finds the multi-variable Taylor exapansion up to z[1]^K...z[m]^K #of exp #(g[1]*z[1]^A[1][1]*z[2]^A[2][1]*z[3]^A[3][1]*...*z[m]^A[m][1]+ #(g[2]*z[1]^A[1][2]*z[2]^A[2][2]*z[3]^A[3][2]*...*z[m]^A[m][2]+ #+.... #(g[n]*z[1]^A[1][m]*z[2]^A[2][m]*z[3]^A[3][2]*...*z[m]^A[m][n]) #For example, try: #Tay([[1,1]],g,z,4); Tay:=proc(A,g,z,K) local m,n,gu,i,j: option remember: m:=nops(A): n:=nops(A[1]): gu:=add(g[j]*mul(z[i]^A[i][j],i=1..m),j=1..n): gu:=exp(gu): for i from 1 to m do gu:=taylor(gu,z[i]=0,K+1): gu:=add(expand(coeff(gu,z[i],j))*z[i]^j,j=0..K): od: gu: end: #C(A,b,g): The normalization constant. For example, try #C([[1,2,1],[2,1,3]],[2,4],g); C:=proc(A,b,g) local gu,i,z,K,m: option remember: m:=nops(A): K:=max(op(b))+1: gu:=Tay(A,g,z,K): for i from 1 to m do gu:=coeff(gu,z[i],b[i]): od: gu: end: #GF1(A,g,z): Given a (numeric or symbolic) matrix #A (say m by n) and a symbol g and a symbol z #finds the multi-variable generating function #in z[1], ..., z[m] of the normalizing constant #For example, try: #GF1([[1,1]],g,z); GF1:=proc(A,g,z) local m,n,gu,i,j: option remember: m:=nops(A): n:=nops(A[1]): gu:=add(g[j]*mul(z[i]^A[i][j],i=1..m),j=1..n): exp(gu): end: #Ave(A,g,b,j): The average of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) g[i]'s #(1<=i<=n). # For example, try #Ave([[1,1]],g,[2],1); Ave:=proc(A,g,b,j) local b1,i1,mu: b1:=[seq(A[i1][j],i1=1..nops(A))]: b1:=[seq(b[i1]-b1[i1],i1=1..nops(b))]: mu:=C(A,b,g): if mu=0 then FAIL: else normal(g[j]*C(A,b1,g)/mu): fi: end: #Ave2(A,g,b,j): Like Ave but a faster version applicable #to two-rowed matrices A #The average of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) g[i]'s #(1<=i<=n). # For example, try #Ave2([[1,1]],g,[2],1); Ave2:=proc(A,g,b,j) local b1,i1,mu: if nops(A)<>2 then RETURN(FAIL): fi: b1:=[seq(A[i1][j],i1=1..nops(A))]: b1:=[seq(b[i1]-b1[i1],i1=1..nops(b))]: mu:=C(A,b,g): if mu=0 then FAIL: else normal(g[j]*C(A,b1,g)/mu): fi: end: #gtoc2(A,g): given a 0-1 matrix A with 2 rows and a vector #g of numbers (or symbols) of size nops(A[1]) (=nops(A[2])) #outputs the vector [c11,c10,c01]. #gtoc2([[1,1,0],[1,0,1]],[2,3,5]); gtoc2:=proc(A,g) local c11,c10,c01,j: option remember: if nops(A)<>2 or nops(A[1])<>nops(A[2]) then RETURN(FAIL): fi: if {op(A[1]),op(A[2])} minus {0,1}<>{} then print(`The matrix A should be a 0-1 matrix`): RETURN(FAIL): fi: c11:=0: c10:=0: c01:=0: for j from 1 to nops(A[1]) do if A[1][j]=1 and A[2][j]=1 then c11:=c11+g[j]: elif A[1][j]=1 and A[2][j]=0 then c10:=c10+g[j]: elif A[1][j]=0 and A[2][j]=1 then c01:=c01+g[j]: fi: od: [c11,c10,c01]: end: C2old:=proc(A,b,g) local c11,c10,c01,j,k: if nops(A)<>2 or nops(b)<>2 then RETURN(FAIL): fi: c11:=0: c10:=0: c01:=0: for j from 1 to nops(A[1]) do if A[1][j]=1 and A[2][j]=1 then c11:=c11+g[j]: elif A[1][j]=1 and A[2][j]=0 then c10:=c10+g[j]: elif A[1][j]=0 and A[2][j]=1 then c01:=c01+g[j]: fi: od: expand( add(c11^k*c10^(b[1]-k)*c01^(b[2]-k)/k!/(b[1]-k)!/(b[2]-k)!,k=0..min(op(b)))): end: #C2(A,b,g): fast version of C(A,B,g) for the special case #of A being a 0-1 matrix with two rows C2:=proc(A,b,g) local c11,c10,c01,k,gu: if nops(b)<>2 then RETURN(FAIL): fi: gu:=gtoc2(A,g): if gu=FAIL then RETURN(FAIL): fi: c11:=gu[1]: c10:=gu[2]: c01:=gu[3]: expand( add(c11^k*c10^(b[1]-k)*c01^(b[2]-k)/k!/(b[1]-k)!/(b[2]-k)!,k=0..min(op(b)))): end: #C2f(A,b,g): fast version of C(A,B,g) for the special case #of A being a 0-1 matrix with two rows C2f:=proc(A,b,g) local c11,c10,c01,gu,b1,b2: option remember: gu:=gtoc2(A,g): b1:=b[1]: b2:=b[2]: if gu=FAIL then RETURN(FAIL): fi: c11:=gu[1]: c10:=gu[2]: c01:=gu[3]: if b1<=1 and b2<=1 then RETURN(C2(A,b,g)): fi: if b2>=2 then RETURN( expand( (1/c10/b2*c11*b1+1/c10/b2*c11-1/c10*c11+1/b2*c01)* C2f(A,[b1,b2-1],g)+1/c10/b2*c11*c01*C2f(A,[b1,b2-2],g))): else RETURN(C2(A,b,g)): fi: end: with(SumTools[Hypergeometric]): FindOpe2:=proc(b1,b2,B1,B2,c01,c10,c11) local F,k,ope1,ope2: F:=c11^k*c10^(b1-k)*c01^(b2-k)/k!/(b1-k)!/(b2-k)!: ope1:=Zeilberger(F,b1,k,B1)[1]: ope2:=Zeilberger(F,b2,k,B2)[1]: ope1:=expand(subs(b1=b1-2,ope1)/B1^2): ope2:=expand(subs(b2=b2-2,ope2)/B2^2): ope1:= - expand((coeff(ope1,B1,-1)*B1^(-1)+coeff(ope1,B1,-2)*B1^(-2))/coeff(ope1,B1,0)): ope2:= - expand((coeff(ope2,B2,-1)*B2^(-1)+coeff(ope2,B2,-2)*B2^(-2))/coeff(ope2,B2,0)): ope1,ope2: end: #ope1:=-1/c01/B1*c11+1/c01/b1/B1*c11+1/c01/b1/B1*c11*b2+1/b1/B1*c10+ #1/c01/b1/B1^2*c11*c10: #ope2:= #1/c10/b2/B2*c11*b1+1/c10/b2/B2*c11-1/c10/B2*c11+1/b2/B2*c01+ #1/c10/b2/B2^2*c11*c01: #C2f(A,b,g): fast version of C(A,B,g) for the special case #of A being a 0-1 matrix with two rows C2fOld:=proc(A,b,g) local c11,c10,c01,gu,b1,b2: option remember: gu:=gtoc2(A,g): b1:=b[1]: b2:=b[2]: if gu=FAIL then RETURN(FAIL): fi: c11:=gu[1]: c10:=gu[2]: c01:=gu[3]: if b1<=1 and b2<=1 then RETURN(C2(A,b,g)): fi: if b2>=2 then RETURN( expand( (1/c10/b2*c11*b1+1/c10/b2*c11-1/c10*c11+1/b2*c01)* C2f(A,[b1,b2-1],g)+1/c10/b2*c11*c01*C2f(A,[b1,b2-2],g))): else RETURN( expand( (-1/c01*c11+1/c01/b1*c11+1/c01/b1*c11*b2+1/b1*c10)*C2f(A,[b1-1,b2],g) +1/c01/b1*c11*c10*C2f(A,[b1-2,b2],g))): fi: end: #Ave2(A,g,b,j): The average of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) g[i]'s #(1<=i<=n). # For example, try #Ave2([[1,1]],g,[2],1); Ave2:=proc(A,g,b,j) local b1,i1,mu: if nops(A)<>2 then print(`The matrix should have two rows`): RETURN(FAIL): fi: b1:=[seq(A[i1][j],i1=1..nops(A))]: b1:=[seq(b[i1]-b1[i1],i1=1..nops(b))]: mu:=C2(A,b,g): if mu=0 then FAIL: else normal(g[j]*C2(A,b1,g)/mu): fi: end: #Ave2f(A,g,b,j): The average of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) g[i]'s #(1<=i<=n). # For example, try #Ave2f([[1,1]],g,[2],1); Ave2f:=proc(A,g,b,j) local b1,i1,mu: if nops(A)<>2 then print(`The matrix should have two rows`): RETURN(FAIL): fi: b1:=[seq(A[i1][j],i1=1..nops(A))]: b1:=[seq(b[i1]-b1[i1],i1=1..nops(b))]: mu:=C2f(A,b,g): if mu=0 then FAIL: else normal(g[j]*C2f(A,b1,g)/mu): fi: end: #####stuff from MultiAlmkvistZeilberger ####################################################################### ## MultiAlmkvistZeilberger Save this file as MultiAlmkvistZeilberger # ## same directory, get into Maple (by typing: maple ) # ## and then type: read MultiAlmkvistZeilberger # ## Then follow the instructions given there # ## # ## Written by Mohamud Mohammed and Doron Zeilberger, # #Rutgers University , # ## [mohamudm, zeilberg] at math dot rutgers dot edu # ####################################################################### ezraMAZ:=proc() if nargs=0 then print(`MultiAlmkvistZeilberger`): print(`A Maple package for finding recurrences for multi-integrals of `): print(` Hypergeometric/Hyperexponential Integrands.`): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(`Contains procedures: MAZ `): elif nargs=1 and args[1]=MAZ then print(`MAZ(POL,H,rat,x,n,N,pars): Given a polynomial, POL, an hyper-exponential function, H, and`): print(`a rational function, rat, of the continuous variables in the list x, and a discrete variable`): print(`n, and the shift-operator in n, N, and a set of auxiliary parameters, par`): print(`ouputs a recurrence operator, let's call it ope(n,n), and a multi-certificate, let's call it R`): print(`such that the function F(n;x):=POL*H*rat^n satisfies`): print(` ope(n,N)F=sum(D_{x[i]} (R[i]H*rat^n),i=1..nops(R)`): print(`In particular, try:`): print(`MAZ(1,((1-1/x)*(1-y))^b*((1-1/x/y)*(1-1/y))^c/x/y,(1-x)*(1-x*y),[x,y],a,A,{b,c});`): print(`MAZ(1,exp(-x^2/2-y^2/2),(x-y)^2,[x,y],n,N,{});`): else print(`There is no help for`, args): fi: end: ezOld:=proc():print(`MAZ1(POL,H,rat,x,n,N,L,pars), MAZ(POL,H,rat,x,n,N,pars)`): print(`FindDeg(POL,H,rat,x,n,L)`): print(`CheckMAZ1(POL,H,rat,x,n,N,L,pars)`): end: _EnvAllSolutions:=true: ####From MultiZeilberger #IV1(d,n): all the vectors of non-negative integres of length d whose #sum is n IV1:=proc(d,n) local gu,i,gu1,i1: if d=0 then if n=0 then RETURN({[]}): else RETURN({}): fi: fi: gu:={}: for i from 0 to n do gu1:=IV1(d-1,n-i): gu:=gu union {seq([op(gu1[i1]),i],i1=1..nops(gu1))}: od: gu: end: yafe:=proc(ope,N) local i: add(factor(coeff(ope,N,i))*N^i,i=ldegree(ope,N)..degree(ope,N)):end: #IV(d,n): all the vectors of non-negative integres of length d whose #sum is <=n IV:=proc(d,n) local i: {seq(op(IV1(d,i)),i=0..n)}: end: #GenPol(kList,a,deg): The generic polynomial in kList of #degree deg, using the indexed variable a, followed by the set #of coeffs. GenPol:=proc(kList,a,deg) local gu,i,i1: gu:=IV(nops(kList),deg): convert([seq(a[i]*convert([seq(kList[i1]^gu[i][i1],i1=1..nops(kList))],`*`), i=1..nops(gu))],`+`),{seq(a[i],i=1..nops(gu))}: end: #Extract1(POL,kList): extracts all the coeffs. of a POL in kList Extract1:=proc(POL,kList) local k1,kList1,i: k1:=kList[nops(kList)]: kList1:=[op(1..nops(kList)-1,kList)]: if nops(kList)=1 then RETURN({seq(coeff(POL,k1,i),i=0..degree(POL,k1))}): else RETURN({seq( op(Extract1(coeff(POL,k1,i),kList1)), i=0..degree(POL,k1))}): fi: end: ###########End from MultiZeilberger FindDeg:=proc(POL,H,rat,x,xSet,n,L) local s,t,h,e,i,Hbar,q,r,gu: s:=numer(rat): t:=denom(rat): h:=convert([seq(e[i]*subs(n=n+i,POL)*s^i*t^(L-i),i=0..L)],`+`): Hbar:=H*s^n/t^(n+L): gu:=normal(simplify(diff(Hbar,x)/Hbar)): q:=numer(gu): r:=denom(gu): max(degree(h,xSet)-degree( diff(r,x)+q,xSet)+degree(r,xSet) ,degree(h,xSet)-degree(r,xSet)+1+degree(r,xSet)): end: #MAZ1(POL,H,rat,x,n,N,L,pars): Given a polynomial POL of the discrete variable n #and the continuous variables in x #and a pure-hyperexponential function H of the variables in x, and a rational function #Rat of x, and a shift operator N, and a non-neg. integer L #returns FAIL or, if there is one, an operator ope(N,n) of order L and #a list of rational functions R, such that F:=POL*H*rat^n satisfies #ope(N,n)F=sum(diff(F*R[i],x[i]),i=1..nops(R)) . #For example, try MAZ1(1,exp(-x),x,[x],n,N,1,{}) MAZ1:=proc(POL,H,rat,x,n,N,L,pars) local deg,deg1,m,gu,i,i1: deg:=[]: for i from 1 to nops(x) do deg:=[op(deg),FindDeg(POL,H,rat,x[i],{seq(x[i1],i1=1..nops(x))},n,L)]: od: m:=min(op(deg)): for i from 0 to m do deg1:=[seq(deg[i1]-m+i,i1=1..nops(deg))]: gu:=MAZ1deg(POL,H,rat,x,n,N,L,pars,deg1): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: MAZ:=proc(POL,H,rat,x,n,N,pars) local gu,L: for L from 0 do gu:=MAZ1(POL,H,rat,x,n,N,L,pars): if gu<>FAIL then RETURN(gu): fi: od: end: CheckMAZ:=proc(POL,H,rat,x,n,N,pars) local F,gu,luL,luR,R,ope,i: F:=H*rat^n: gu:=MAZ(POL,H,rat,x,n,N,pars): if gu=FAIL then RETURN(FAIL): fi: ope:=gu[1]: R:=gu[2]: luL:=normal(add(subs(n=n+i,POL)*coeff(ope,N,i)*normal(simplify(subs(n=n+i,F)/F)),i=0..degree(ope,N))): luR:=normal(add( normal(diff(R[i]*F,x[i])/F ), i=1..nops(R))): evalb(normal(luL/luR)=1): end: #MAZ1deg(POL,H,rat,x,n,N,L,pars,deg): Given a polynomial POL of the discrete variable n #and the continuous variables in x #and a pure-hyperexponential function H of the variables in x, and a rational function #Rat of x, and a shift operator N, and a non-neg. integer L #returns FAIL or, if there is one, an operator ope(N,n) of order L and #a list of rational functions R, such that F:=POL*H*rat^n satisfies #ope(N,n)F=sum(diff(F*R[i],x[i]),i=1..nops(R)) . #For example, try MAZ1(1,exp(-x),x,[x],n,N,1,{}) MAZ1deg:=proc(POL,H,rat,x,n,N,L,pars,deg) local s,t,h,e,i,Hbar,gu,X,a,var,q,r,eq,ope,var1,ku,i1,eqN,var1N,opeN,bilti,meka: s:=numer(rat): t:=denom(rat): ope:=add(e[i]*N^i,i=0..L): h:=convert([seq(e[i]*subs(n=n+i,POL)*s^i*t^(L-i),i=0..L)],`+`): Hbar:=H*s^n/t^(n+L): gu:=[]: for i from 1 to nops(x) do gu:=[op(gu),normal(simplify(diff(Hbar,x[i])/Hbar))]: od: q:=[]: r:=[]: for i from 1 to nops(gu) do q:=[op(q),numer(gu[i])]: r:=[op(r),denom(gu[i])]: od: var:={}: X:=[]: for i from 1 to nops(x) do ku:=GenPol(x,a[i],deg[i]): X:=[op(X),ku[1]]: var:=var union ku[2]: od: var:=var union {seq(e[i],i=0..L)}: gu:=h: for i from 1 to nops(x) do gu:=expand(normal(gu-(diff(r[i],x[i])+q[i])*X[i]/r[i]-r[i]*diff(X[i]/r[i],x[i]))): od: eq:=Extract1(numer(gu),x): eqN:=subs(n=9/17,eq): eqN:=subs({seq(pars[i]=1/(i^2+3),i=1..nops(pars))},eqN): var1N:=solve(eqN,var): opeN:=subs(var1N,ope): if opeN=0 then RETURN(FAIL): fi: var1:=solve(eq,var): ope:=subs(var1,ope): X:=subs(var1,X): if ope=0 then RETURN(FAIL): fi: bilti:={}: for i from 1 to nops(var1) do if op(1,op(i,var1))=op(2,op(i,var1)) then bilti:=bilti union {op(1,op(i,var1))}: fi: od: for i from 1 to nops(bilti) do if coeff(ope,bilti[i],1)<>0 then ope:=coeff(ope,bilti[i],1): X:=[seq(coeff(X[i1],bilti[i],1),i1=1..nops(X))]: meka:=coeff(ope,N,degree(ope,N)): ope:=expand(normal(ope/meka)): ope:=yafe(ope,N): X:=[seq(normal(X[i1]/meka/t^L),i1=1..nops(X))]: RETURN(ope,X): fi: od: end: ############end stuff from MultiAlmkvistZeilberger #Opes(F,z,b,B): the recurrence operators annihilating #F/z[1]^(b[1]+1)/.../z[n]^(b[n]+1) #using B[i] as shifts (n=nops(b)=nops(z)=nops(B)) For example #Opes(exp(x+y),[x,y],[a,b],[A,B]); Opes:=proc(F,z,b,B) local i,j: option remember: if nops(z)<>nops(b) then RETURN(`the second and third arguments must have the same length`): fi: [seq( MAZ(1,F/mul(z[j]^(b[j]+1),j=1..i-1)/mul(z[j]^(b[j]+1),j=i+1..nops(z))/z[i], 1/z[i],z,b[i],B[i], {seq(b[j],j=1..i-1), seq(b[j],j=i+1..nops(b))})[1], i=1..nops(b))]: end: #Hatkhala(F,z,K): The list of lists of lists etc. such that #L[i1][i2][i3]... is the coeff. of z1^i1*z2^i2*z3^i3... #for 1<=i1<=K1, 1<=i2<=K2, ... #in the multivarite expanson of F(z[1],z[2],z[3], ...) #For example try: #Hatkhala(exp(z1+z2+z1*z2),[z1,z2],[3,4]); Hatkhala:=proc(F,z,K) local n,gu,i,gu1,lu: n:=nops(z): if nops(K)<>n then print(`Bad input`): RETURN(FAIL): fi: gu:=taylor(F,z[1]=0,K[1]+2): if n=1 then RETURN([seq(coeff(gu,z[1],i),i=1..K[1])]): fi: lu:=[]: for i from 1 to K[1] do gu1:=coeff(gu,z[1],i): lu:=[op(lu),Hatkhala(gu1,[op(2..n,z)],[op(2..nops(K),K)])]: od: lu: end: #RecsOld(A,g,b,B): Given a (numeric or symbolic) matrix #A (say m by n) and a symbol g and a symbol z #finds the recurrences for the appropriate C #in terms of b[1], ..., b[m] and B[1], ..., B[m] #For example, try: #RecsOld([[1,1]],g,b,B); RecsOld:=proc(A,g,b,B) local m,n,gu,i,j,z,opes: option remember: m:=nops(A): n:=nops(A[1]): gu:=add(g[j]*mul(z[i]^A[i][j],i=1..m),j=1..n): gu:=exp(gu): opes:=Opes(gu,[seq(z[i],i=1..m)],[seq(b[i],i=1..m)],[seq(B[i],i=1..m)]): opes,Hatkhala(gu,[seq(z[i],i=1..m)],[seq(degree(opes[i],B[i]),i=1..m)]): end: Yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: #SeqFromRec(ope,n,N,Ini,K): Given the first L-1 #terms of the sequence Ini=[f(1), ..., f(L-1)] #satisfied by the recurrence ope(n,N)f(n)=0 #extends it to the first K values SeqFromRec:=proc(ope,n,N,Ini,K) local ope1,gu,L,n1,j1: ope1:=Yafe(ope,N)[2]: L:=degree(ope1,N): if nops(Ini)<>L then ERROR(`Ini should be of length`, L): fi: ope1:=expand(subs(n=n-L,ope1)/N^L): gu:=Ini: for n1 from nops(Ini)+1 to K do gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)), j1=1..L)]: od: gu: end: #ValueFromRec1(ope,n,N,Ini,K): Given the first L #terms of the sequence Ini=[f(1), ..., f(L)] #satisfied by the recurrence ope(n,N)f(n)=0 #finds f(K) #For example, try: #ValueFromRec1(N^2-N-1,n,N,[1,1],10); ValueFromRec1:=proc(ope,n,N,Ini,K) local ope1,gu,L,n1,j1: ope1:=Yafe(ope,N)[2]: L:=degree(ope1,N): if nops(Ini)<>L then ERROR(`Ini should be of length`, L): fi: if K<=nops(Ini) then RETURN(Ini[K]): fi: ope1:=expand(subs(n=n-L,ope1)/N^L): gu:=Ini: for n1 from nops(Ini)+1 to K do gu:=[op(2..nops(gu),gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)), j1=1..L)]: od: expand(gu[nops(gu)]): end: #RecsG(F,z,b,B): Given a (numeric or symbolic) matrix #A (say m by n) and a symbol g and a symbol z #finds the recurrences for the appropriate C #in terms of b[1], ..., b[m] and B[1], ..., B[m] #For example, try: #RecsG(exp(x+y),[x,y],b,B); RecsG:=proc(F,z,b,B) local m,i,opes: option remember: m:=nops(z): opes:=Opes(F,z,[seq(b[i],i=1..m)],[seq(B[i],i=1..m)]): opes,Hatkhala(F,[seq(z[i],i=1..m)],[seq(degree(opes[i],B[i]),i=1..m)]): end: #Recs(A,g,b,B): Given a (numeric or symbolic) matrix #A (say m by n) and a symbol g and a symbol z #finds the recurrences for the appropriate C #in terms of b[1], ..., b[m] and B[1], ..., B[m] #For example, try: #Recs([[1,1]],g,b,B); Recs:=proc(A,g,b,B) local m,n,gu,i,j,z: option remember: m:=nops(A): n:=nops(A[1]): gu:=add(g[j]*mul(z[i]^A[i][j],i=1..m),j=1..n): gu:=exp(gu): RecsG(gu,[seq(z[i],i=1..m)],[seq(b[i],i=1..m)],[seq(B[i],i=1..m)]): end: #ValueFromRecs(opes,n,N,Inis,K): given a list of operators #opes, a list of discrete variables n, a list of the corresponding #Initial Values, and a vector the same length as opes (n and N) #computes the value of the multisequence at K #For example, try: #ValueFromRecs([N1-2,N2-3],[n1,n2],[N1,N2],[[1]],[4,5]); ValueFromRecs:=proc(opes,n,N,Inis,K) local k,opes1,n1, ope1,lu,j: k:=nops(opes): if nops(n)<>k or nops(N)<>k or nops(K)<>k then print(`Bad input`): RETURN(FAIL): fi: if k=1 then RETURN(ValueFromRec1(opes[1],n[1],N[1],Inis,K[1])): fi: lu:=[]: for n1 from 1 to nops(Inis) do opes1:=subs(n[1]=n1,[op(2..k,opes)]): lu:= [op(lu), ValueFromRecs(opes1,[op(2..k,n)],[op(2..k,N)],Inis[n1],[op(2..k,K)]) ]: od: ope1:=subs({seq(n[j]=K[j],j=2..k)},opes[1]): ValueFromRec1(ope1,n[1],N[1],lu,K[1]): end: #CheckRecsG(F,z,L,kama): Checks RecsG(F,a,b,B) for #up to L-degree terms #For example, try: #CheckRecsG(exp(x+y),[x,y],10); CheckRecsG:=proc(F,z,L,kama) local lu,i,ra,mu,n,v,c,ka1,ka2,b,B,kv: n:=nops(z): mu:=Hatkhala(F,z,[L$n]): ra:=rand(1..L): lu:=RecsG(F,z,b,B): kv:={}: for c from 1 to kama do v:=[seq(ra(),i=1..n)]: ka1:=ValueFromRecs(lu[1], [seq(b[i],i=1..nops(z))], [seq(B[i],i=1..nops(z))], lu[2],v): ka2:=mu: for i from 1 to n do ka2:=ka2[v[i]]: od: kv:=kv union {ka2-ka1}: od: evalb(kv={0}): end: #Cf(A,b,g): The normalization constant. For example, try #Cf([[1,2,1],[2,1,3]],[2,4],g); Cf:=proc(A,v,g) local b,B,lu,i,m: option remember: m:=nops(v): lu:=Recs(A,g,b,B): ValueFromRecs(lu[1],[seq(b[i],i=1..m)], [seq(B[i],i=1..m)],lu[2],v): end: #Avef(A,g,b,j): The average of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) g[i]'s #(1<=i<=n). # For example, try #Avef([[1,1]],g,[2],1); Avef:=proc(A,g,b,j) local b1,i1,mu: b1:=[seq(A[i1][j],i1=1..nops(A))]: b1:=[seq(b[i1]-b1[i1],i1=1..nops(b))]: mu:=Cf(A,b,g): if mu=0 then FAIL: else normal(g[j]*Cf(A,b1,g)/mu): fi: end: sn:=proc(resh,n): -1/log(op(n+1,resh)*op(n-1,resh)/op(n,resh)^2): #evalf(%): end: Zinn:=proc(resh) local s1,s2,n: n:=nops(resh)-2: s1:=sn(resh,n): s2:=sn(resh,n-1): evalf(2*(s1+s2)/(s1-s2)^2), evalf(sqrt(op(n+1,resh)/op(n-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))): end: #ZinnC(resh,eps,V): Like Zinn(resh), but with the constant #and with confidence eps. and in terms of the variable V #returns the conjectured asymptotics phrased in terms of V #For example, try: #ZinnC([seq(2*i/i,i=1..100)],.1); ZinnC:=proc(resh,eps,V) local s1,s2,n,mu,theta,C1,C2: n:=nops(resh)-2: s1:=sn(resh,n): s2:=sn(resh,n-1): theta:=evalf(2*(s1+s2)/(s1-s2)^2): mu:=evalf(sqrt(op(n+1,resh)/op(n-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))): C1:=resh[n+2]/(n+2)^theta/mu^(n+2): C2:=resh[n+1]/(n+1)^theta/mu^(n+1): if abs(C1/C2-1)>eps then RETURN(FAIL): fi: C1*mu^V*V^theta: end: #FM1(A,g,b,r,j): The j-th factorial #moment of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #FM1([[1,1]],g,[2],2,1); FM1:=proc(A,g,b,j,r) local b1,i1,mu: if r>=min(op(b)) then print(`The minimum of the vector`, b, `should be larger than`, r): RETURN(FAIL): fi: b1:=[seq(A[i1][j],i1=1..nops(A))]: b1:=[seq(b[i1]-r*b1[i1],i1=1..nops(b))]: mu:=C(A,b,g): if mu=0 then FAIL: else normal(g[j]^r*C(A,b1,g)/mu): fi: end: #FM1f(A,g,b,r,j): The r-th factorial #moment of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #FM1f([[1,1]],g,[2],2,1); FM1f:=proc(A,g,b,j,r) local b1,i1,mu: if r>=min(op(b)) then RETURN(FAIL): fi: b1:=[seq(A[i1][j],i1=1..nops(A))]: b1:=[seq(b[i1]-r*b1[i1],i1=1..nops(b))]: mu:=Cf(A,b,g): if mu=0 then FAIL: else normal(g[j]^r*Cf(A,b1,g)/mu): fi: end: #FM12(A,g,b,r,j): The r-th factorial #moment of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #FM12([[1,1]],g,[2],2,1); FM12:=proc(A,g,b,j,r) local b1,i1,mu: if r>=min(op(b)) then RETURN(FAIL): fi: b1:=[seq(A[i1][j],i1=1..nops(A))]: b1:=[seq(b[i1]-r*b1[i1],i1=1..nops(b))]: mu:=C2(A,b,g): if mu=0 then FAIL: else normal(g[j]^r*C2(A,b1,g)/mu): fi: end: #FM12f(A,g,b,r,j): The r-th factorial #moment of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #FM12f([[1,1]],g,[2],2,1); FM12f:=proc(A,g,b,j,r) local b1,i1,mu: if r>=min(op(b)) then RETURN(FAIL): fi: b1:=[seq(A[i1][j],i1=1..nops(A))]: b1:=[seq(b[i1]-r*b1[i1],i1=1..nops(b))]: mu:=C2f(A,b,g): if mu=0 then FAIL: else normal(g[j]^r*C2f(A,b1,g)/mu): fi: end: #Var(A,g,b,j): The variance of #of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #Var([[1,1]],g,[2],1); Var:=proc(A,g,b,j) local mu: mu:=Ave(A,g,b,j): if 2>=min(op(b)) then print(`The minimum of the vector`, b, `should be larger than`, 2): RETURN(FAIL): fi: FM1(A,g,b,j,2)+mu-mu^2: end: #Varf(A,g,b,j): The variance of #of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #Varf([[1,1]],g,[2],1); Varf:=proc(A,g,b,j) local mu: mu:=Avef(A,g,b,j): FM1f(A,g,b,j,2)+mu-mu^2: end: #Var2f(A,g,b,j): The variance of #of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #Var2f([[1,1]],g,[2],1); Var2f:=proc(A,g,b,j) local mu: if nops(A)<>2 then print(`The matrix should have two rows`): RETURN(FAIL): fi: mu:=Ave2f(A,g,b,j): FM12f(A,g,b,j,2)+mu-mu^2: end: #Var2(A,g,b,j): The variance of #of x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #Var2([[1,1]],g,[2],1); Var2:=proc(A,g,b,j) local mu,i1,j1: if nops(A)<>2 then print(`The matrix should have two rows`): RETURN(FAIL): fi: if {seq(seq(A[i1][j1],j1=1..nops(A[i1])),i1=1..nops(A))} minus {0,1}<>{} then ERROR(`Not a 0-1 matrix`): fi: mu:=Ave2(A,g,b,j): if 2>=min(op(b)) then print(`The minimum of the vector`, b, `should be larger than`, 2): RETURN(FAIL): fi: FM12(A,g,b,j,2)+mu-mu^2: end: #SipurD(A,g,j,K,n): The diagonal story for the matrix A #Poisson parameters g, the variable x[j], up to #[K, ..., K] using the variable n for asymptotics #using the slow procedure C. For example, try: #SipurD(A1,[1$nops(A1[1])],1,20,n); SipurD:=proc(A,g,j,K,n) local gu,i,x: print(`This is the story of the averages, variances `): print(`of the random variable`, x[j]): print(`under the constraint Ax=b`): print(`for the matrix`): print(convert(A,matrix)): print(`with the Poisson parameters given by`, g): print(`for b from `, [1$nops(A)], `to `, [K$nops(A)] ): gu:=evalf([seq(Ave(A,g,[i$nops(A)],j),i=1..K)]): print(`The sequence of averages is`): print(gu): print(`This seems to be asymptotically`): print(ZinnC(gu,.2,n)): gu:=evalf([seq(sqrt(Var(A,g,[i$nops(A)],j)),i=3..K)]): print(`The sequence of standard deviations is`): print(evalf(gu)): print(`This seems to be asymptotically`): ZinnC(gu,.2,n): end: #SipurDf(A,g,j,K,n): Like SipurD(A,g,K,n) but using #the Almkvist-Zeilberger recurrences #For example, try: #SipurDf(A2,[1$nops(A2[1])],1,20,n); SipurDf:=proc(A,g,j,K,n) local gu,i,x: print(`This is the story of the averages and variance`): print(`of the random variable`, x[j]): print(`under the constraint Ax=b`): print(`for the matrix`): print(convert(A,matrix)): print(`with the Poisson parameters given by`, g): print(`for b from `, [3$nops(A)], `to `, [K$nops(A)] ): gu:=evalf([seq(Avef(A,g,[i$nops(A)],j),i=3..K)]): print(`The sequence of averages is`): print(gu): print(`This seems to be asymptotically`): print(ZinnC(gu,.1,n)): gu:=evalf([seq(sqrt(Varf(A,g,[i$nops(A)],j)),i=3..K)]): print(`The sequence of standard deviations is`): print(evalf(gu)): print(`This seems to be asymptotically`): print(ZinnC(gu,.1,n)): end: #SipurD2f(A,g,j,K,n): Like SipurD(A,g,K,n) but using #the Zeilberger recurrence, and only applicable to two-rowed matrices #For example, try: #SipurD2f(A2,[1$nops(A2[1])],1,20,n); SipurD2f:=proc(A,g,j,K,n) local gu,i,x,C: if nops(A)<>2 then print(`The matrix should have two rows`): RETURN(FAIL): fi: print(`This is the story of the averages and variance`): print(`of the random variable`, x[j]): print(`under the constraint Ax=b`): print(`for the matrix`): print(convert(A,matrix)): print(`with the Poisson parameters given by`, g): print(`for b from `, [3$nops(A)], `to `, [K$nops(A)] ): gu:=evalf([seq(Ave2f(A,g,[i$nops(A)],j),i=3..K)]): print(`The sequence of averages is`): print(gu): print(`This seems to be asymptotically`): print(ZinnC(gu,.1,n)): if abs(Zinn(gu)[1]-1/2)<.1 then print(`If we assume that this grows like a constant times`, sqrt(n)): print(`then it is reasobale that the asymptotics of the averages is`): C:=evalf(gu[nops(gu)]/sqrt(K)): print(C*sqrt(n)): fi: gu:=evalf([seq(sqrt(Var2f(A,g,[i$nops(A)],j)),i=3..K)]): print(`The sequence of standard deviations is`): print(evalf(gu)): print(`This seems to be asymptotically`): print(ZinnC(gu,.1,n)): if abs(Zinn(gu)[1]-1/4)<.1 then print(`If we assume that this grows like a constant times`, n^(1/4)): print(`then it is reasobale that the asymptotics of the averages is`): C:=evalf(gu[nops(gu)]/K^(1/4)): print(C*n^(1/4)): fi: end: #AllAve(A,g,b): all the nops(A[1]) averages AllAve:=proc(A,g,b) local j,gu,i,b1: gu:=[seq(Ave(A,g,b,j),j=1..nops(A[1]))]: b1:=normal([seq(add(A[i][j]*gu[j],j=1..nops(gu)),i=1..nops(A))]): if {seq(evalb(evalf(abs(b[i]-b1[i]))<=10^(-Digits/2)),i=1..nops(b))} <>{true} then print(`Something is wrong, the image of the average point`): print(`in not`, b, `but rather `, b1): fi: gu: end: #AllAvef(A,g,b): all the nops(A[1]) averages AllAvef:=proc(A,g,b) local j,b1,gu,i: gu:=[seq(Avef(A,g,b,j),j=1..nops(A[1]))]: b1:=normal([seq(add(A[i][j]*gu[j],j=1..nops(gu)),i=1..nops(A))]): if b<>b1 then print(`Something is wrong, the image of the average point`): print(`in not`, b, `but rather `, b1): fi: gu: end: #AllAve2(A,g,b): all the nops(A[1]) averages AllAve2:=proc(A,g,b) local j,b1,gu,i: gu:=[seq(Ave2(A,g,b,j),j=1..nops(A[1]))]: b1:=normal([seq(add(A[i][j]*gu[j],j=1..nops(gu)),i=1..nops(A))]): if b<>b1 then print(`Something is wrong, the image of the average point`): print(`in not`, b, `but rather `, b1): fi: gu: end: #AllAve2f(A,g,b): all the nops(A[1]) averages AllAve2f:=proc(A,g,b) local i,j,b1,gu: gu:=[seq(Ave2f(A,g,b,j),j=1..nops(A[1]))]: b1:=normal([seq(add(A[i][j]*gu[j],j=1..nops(gu)),i=1..nops(A))]): if b<>b1 then print(`Something is wrong, the image of the average point`): print(`in not`, b, `but rather `, b1): fi: gu: end: #RecsV(A,g,b): a Verbose version of Recs(A,g,b,B) #RecsV(A2,g,b); RecsV:=proc(A,g,b) local gu, B,m,n,F,i,j,j1: m:=nops(A): n:=nops(A[1]): if type(g,list) and nops(g)<>n then print(`Bad input`): RETURN(FAIL): fi: gu:=Recs(A,g,b,B): print(`Let `, F(seq(b[i],i=1..m)) , `be the (unnormalized) weight`): print(` of the range of`, A, `being `, [seq(b[i],i=1..m)] ): print(`It satisfies the following pure linear recurrece equations`): print(`with polynomial coefficients `): for i from 1 to nops(gu[1]) do print( add(coeff(gu[1][i],B[i],j)* F(seq(b[j1],j1=1..i-1),b[i]+j,seq(b[j1],j1=i+1..m)), j=0..degree(gu[1][i],B[i]))=0) : print(): od: print(`These recurrences enable to compute, in linear time any`): print(`desired value of `, F(seq(b[i],i=1..m))): print( `subject to the inital conditions`): print(gu[2]): 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: ###Covariance moments #Cov(A,g,b,i,j): The covarinace #of x[i] and x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #Cov([[1,1]],g,[2],1,2); Cov:=proc(A,g,b,i,j) local b1,b2,i1,mu: if i=j then RETURN(Var(A,g,b,i)): fi: b1:=[seq(A[i1][j],i1=1..nops(A))]: b2:=[seq(A[i1][i],i1=1..nops(A))]: b1:=[seq(b[i1]-b1[i1]-b2[i1],i1=1..nops(b))]: mu:=C(A,b,g): if mu=0 then FAIL: else normal(g[i]*g[j]*C(A,b1,g)/mu -Ave(A,g,b,i)*Ave(A,g,b,j)): fi: end: #Covf(A,g,b,i,j): The covarinace using Cf #of x[i] and x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #Covf([[1,1]],g,[2],1,2); Covf:=proc(A,g,b,i,j) local b1,b2,i1,mu: if i=j then RETURN(Varf(A,g,b,i)): fi: b1:=[seq(A[i1][j],i1=1..nops(A))]: b2:=[seq(A[i1][i],i1=1..nops(A))]: b1:=[seq(b[i1]-b1[i1]-b2[i1],i1=1..nops(b))]: mu:=Cf(A,b,g): if mu=0 then FAIL: else normal(g[i]*g[j]*Cf(A,b1,g)/mu -Avef(A,g,b,i)*Avef(A,g,b,j)): fi: end: #Cov2f(A,g,b,i,j): The covarinace using C2f #of x[i] and x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #Cov2f([[1,1]],g,[2],1,2); Cov2f:=proc(A,g,b,i,j) local b1,b2,i1,mu: if nops(A)<>2 then print(`A must have two rows`): RETURN(FAIL): fi: if i=j then RETURN(Var2f(A,g,b,i)): fi: b1:=[seq(A[i1][j],i1=1..nops(A))]: b2:=[seq(A[i1][i],i1=1..nops(A))]: b1:=[seq(b[i1]-b1[i1]-b2[i1],i1=1..nops(b))]: mu:=C2f(A,b,g): if mu=0 then FAIL: else normal(g[i]*g[j]*C(A,b1,g)/mu -Ave2f(A,g,b,i)*Ave2f(A,g,b,j)): fi: end: ####end covariances ####start correlation #Cor(A,g,b,i,j): The correlation #of x[i] and x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #Cor([[1,1]],g,[2],1,2); Cor:=proc(A,g,b,i,j): Cov(A,g,b,i,j)/sqrt(Var(A,g,b,i)*Var(A,g,b,j)): end: #Corf(A,g,b,i,j): The correlation using the recurrences #of x[i] and x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #Corf([[1,1]],g,[2],1,2); Corf:=proc(A,g,b,i,j): Covf(A,g,b,i,j)/sqrt(Varf(A,g,b,i)*Varf(A,g,b,j)): end: #Cor2f(A,g,b,i,j): The correlation using the recurrences #of x[i] and x[j] in the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #Cor2f(A2,g,[4,5],1,2); Cor2f:=proc(A,g,b,i,j): if nops(A)<>2 then print(`A must have two rows`): RETURN(FAIL): fi: Cov2f(A,g,b,i,j)/sqrt(Var2f(A,g,b,i)*Var2f(A,g,b,j)): end: #CorM(A,g,b): The correlation matrix #of the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #CorM([[1,1]],g,e[2]); CorM:=proc(A,g,b) local i,j: [seq([seq(Cor(A,g,b,i,j),j=1..nops(A[1]))],i=1..nops(A[1]))]: end: #CorMf(A,g,b): The correlation matrix using Cf #of the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #CorMf([[1,1]],g,e[2]); CorMf:=proc(A,g,b) local i,j: [seq([seq(Corf(A,g,b,i,j),j=1..nops(A[1]))],i=1..nops(A[1]))]: end: #CorM2f(A,g,b): The correlation matrix using C2f #of the Eduardo Sontag prob. dist. given #by the matrix A and the vector b, with generic (symbolic) #or numeric g[i]'s #(1<=i<=n). # For example, try #CorM2f([[1,1]],g,e[2]); CorM2f:=proc(A,g,b) local i,j: [seq([seq(Cor2f(A,g,b,i,j),j=1..nops(A[1]))],i=1..nops(A[1]))]: end: