read `EKHAD.txt`:
read `AsyRec.txt`:

with(combinat):

Help:=proc()
if args=NULL then
print(`RecToLiElt(rec, li, n, N)`, `RecToLi(rec, li,M,n,N )`):
print(`MixedMomentRecurrence(S,n,N,power)`):
print(`MixedMomentCount(S,M,power)`, `StraightMoment(S,M,i,p)`):
print(`ScaledMoment(S,M,i,p)`, `MixedMomentAboutMean(S,M,a1,a2,p1,p2)`):
print( `NormalMoment(S,M,a1,a2,p1,p2)`, `Correlation(S,M,a1,a2)`):
print(`ScaledMixedMomentAboutMean(S,M,a1,a2,p1,p2)`):
print(`Story2(S,M,s1,s2,po), Story(S,M,p1,p2)`):

elif args[1] = 'MixedMomentRecurrence' then
print(`MixedMomentRecurrence(S,n,N,power) finds a recurrence for the mixed moment of X_i, a randomly`):
print(`chosen ordered tree whose nodes have a number of children in S.`):
print(`Precisely, the recurrence gives E[mul(X_i^power[i+1], i in S)]`):
print(`This function also finds enough initial terms of the sequence to define it uniquely`):
print(`For example, MixedMomentRecurrence({0,1,3},n,N,[1,0,0,2]) gives a list of two elements.`): 
print(`The first is a recurrence for the expectation of the product of the number of leaves of a (uniformly)`):
print(`randomly chosen tree with the square of the number of nodes with 3 children.`):
print(`The second is the first 5 terms of the sequence`):
print(`For Example: MixedMomentRecurrence({0,1,3},n,N,[1,0,0,2]);`):

elif args[1] = 'MixedMomentCount' then
print(`MixedMomentCount(S,M,power) gives the first M terms of the sequence defined by `):
print(`MixedMomentRecurrence(S,n,N,power)`):
print(`For Example: MixedMomentCount({0,3,4,5},20,[1,0,0,2,3,1]);`):


elif args[1] = 'StraightMoment' then
print(`StraightMoment(S,M,i,p) computes the first M terms of the sequence of E[(X_{i,n}-mu)^p] `):
print(`where X_{i,n} is the number of nodes with i children in a random ordered tree chosen (uniformly)`):
print(`from all trees on n nodes whose nodes all have a number of children in S and mu is `):
print(`the average value of X_{i,n}`):
print(`For Example: StraightMoment({0,1,4},30,1,3);`):

elif args[1] = 'ScaledMoment' then
print(`ScaledMoment(S,M,i,p) computes StraightMoment(S,M,i,p)/StraightMoment(S,M,i,2)^(p/2)`):
print(`For Example: ScaledMoment({0,1,4},30,1,3);`):

elif args[1] = 'MixedMomentAboutMean' then
print(`MixedMomentAboutMean(S,M,a1,a2,p1,p2) computes the first M terms of the sequence `):
print(`E[(X_{a1,n}-mu1)^p1(X_{a2,n}-mu2)^p2] where X_{i,n} is the number of nodes with i `):
print(`children in a random ordered tree chosen (uniformly) from all trees on n nodes whose nodes all have`):
print(`a number of children in S, mu1 is the average value of X_{a1,n} and mu2 is the average value of X_{a2,n}`):
print(`For Example: MixedMomentAboutMean({0,1,4},10,0,4,3,5);`):

elif args[1] = 'ScaledMixedMomentAboutMean' then
print(`ScaledMixedMomentAboutMean(S,M,a1,a2,p1,p2) computes `):
print(`MixedMomentAboutMean(S,M,a1,a2,p1,p2)/(MixedMomentAboutMean(S,M,a1,a2,2,0)^(p1/2)*MixedMomentAboutMean(S,M,a1,a2,0,2)^(p2/2))`):
print(`For Example: ScaledMixedMomentAboutMean({0,1,4},10,0,4,3,5);`):

elif args[1] = 'NormalMoment' then
print(`NormalMoment(S,M,a1,a2,p1,p2) computes the first M terms of the sequence E[Y_{a1,n}^p1Y_{a2,n}^p2] where the Ys jointly follow a normal distribution with mean 0, variance 1, and correlation`):
print(`equal to the correlation of X_{a1,n} and X_{a2,n}`):
print(`For Example: NormalMoment({0,1,4},10,0,4,3,5);`):


elif args[1] = 'RecToLiElt' then
print(`RecToLiElt(rec,li,n,N) gives the next term of the sequence defined by the recurrence rec with initial terms li`):
print(`For Example: RecToLiElt(-1-N+N^2,[1,1],n,N);`):

elif args[1] = 'RecToLi' then
print(`RecToLi(rec,li,M n,N) gives the Mth term of the sequence defined by the recurrence rec with initial terms li`):
print(`For Example: RecToLi(-1-N+N^2,[1,1],10,n,N);`):

elif args[1] = 'Correlation' then
print(`Correlation(S,M,a1,a2) gives the correlation between the number nodes with a1`):
print(`and the number with a2 childrenin a random tree on M vertices with child counts in S`):
print(`For Example: Correlation({0,1,2,4},30,2,4);`):

elif args[1] = 'Story2' then
print(`Story1(S,M,s1,s2,po) first finds the expectation and variance of X_s1 and X_s2. Then, it computes all scaled mixed moments about the mean of X_s1 and X_s2 of order (p1, p2) for p1,p2 < po. It also computes what these moments should be if X_s1 and X_s2 were jointly normal with means 0 and variances 1 and their empirically observed correlation. It prints all of this information in a verbose manner.`):
print(`For Example: Story2({0,1,2,4},30,0,2,3);`):

elif args[1] = 'Story' then
print(`Story(S,M,p) finds every possible combination of s1 and s2 from S. Then, this procedure runs Story2(S,M,p1,p2) for every such choice.`):
print(`For Example: Story({0,1,2,4},30,3);`):
fi:
end:

RecToLiElt:=proc(rec,li,n,N) local f,terms,i,co,finalco,trec,nval:
trec := rec:
#terms := nops(rec):
i:=0:
co:=[]:
while trec <> coeff(trec,N,0) do
 i:=i+1:
 co:=[op(co),coeff(trec,N^i)]:
 trec := trec - coeff(trec,N^i)*N^i:
od:
co:=[rec-(add(co[i]*N^i,i=1..nops(co))),op(co)]:
nval:=nops(li)-nops(co)+2:
simplify(subs(n=nval,add(li[-j]*co[-(j+1)],j=1..nops(co)-1)/(-co[-1]))):
end:

RecToLi:=proc(rec,li,M,n,N) local elt,temp_li:
if nops(li)>=M then
 return(li[1..M]):
end:
temp_li:=li:
while nops(temp_li) < M do
 elt := RecToLiElt(rec,temp_li,n,N):
 temp_li:=[op(temp_li), elt]:
od:
temp_li:
end:

#MixedMoment(power,y,poly) applies the operator prod(y_i d/dy_i)^power[i+1] to the polynomial
#poly. Note: power MUST have a term corresponding to each y_i (terms may be zero)
MixedMoment:=proc(power,y,poly) local n,temp_poly,i,j:
n:=nops(power):
temp_poly := poly:
for i from 0 to n-1 do
 for j from 1 to power[i+1] do
  temp_poly:=diff(temp_poly,y[i])*y[i]:
 od:
od:
subs({seq(y[i]=1,i=0..n-1)}, temp_poly):
end:

#MixedMomentRecurrence(S,n,N,power) finds a recurrence for the mixed moment of X_i, a 
#randomly chosen ordered tree whose nodes have a number of children in S
MixedMomentRecurrence:=proc(S,n,N,power) local y,A,i,rec,deg,init,m,Roots,max_root,r:
A:=1/n*(add(y[i]*x^i, i in S union {0}))^n/x^n:
rec:=AZd(MixedMoment(power,y,A),x,n,N)[1]:
deg := degree(rec,N):
Roots:=roots(coeff(rec, N, degree(rec, N))):
max_root:=max(seq(r[1],r in Roots)):
init:=[seq(MixedMoment(power,y,GeneralInitTreeCount(S,m,y)),m=1..max(deg,ceil(max_root)+deg+1))]:
[rec,init]:
end:

#MixedMomentCount(S,M,power) returns the first M terms of the sequence of mixed
#moments of X_{i,n} where X_i is defined as above for trees on n vertices
MixedMomentCount:=proc(S,M,power) local rec,n,N:
rec:=MixedMomentRecurrence(S,n,N,power):
RecToLi(rec[1],rec[2],M,n,N):
end:

#StraightMoment(S,M,i,p) X_i is the number of nodes with i children in a random tree on M nodes
#and you want to know the pth moment
StraightMoment:=proc(S,M,i,p) local mo,power,n,N:
power:=[0$(max(S)+1)]:
power[i+1]:=1: 
mo:=[seq(MixedMomentCount(S,M,power*j)[-1],j=0..p)]:
evalf((1/mo[1]^p)*add((-1)^r*binomial(p,r)*mo[2]^r*mo[1]^(p-r-1)*mo[p-r+1],r=0..p)):
end:

#ScaledMoment(S,M,i,p)
ScaledMoment:=proc(S,M,i,p) local sec, mo,n,N,y:
sec:=StraightMoment(S,M,i,2):
mo:=StraightMoment(S,M,i,p):
evalf(mo/sec^(p/2)):
end:


MixedMomentAboutMean:=proc(S,M,a1,a2,p1,p2) local n,N,r,l,moments,power, numtreesrec,numtrees,mean1rec,mean1,mean2rec,mean2,mmrec,mmlist,newmm:
power:=[0$(max(S)+1)]:
numtreesrec:=MixedMomentRecurrence(S,n,N,power):
numtrees:=RecToLi(numtreesrec[1],numtreesrec[2],M,n,N)[-1]:
power[a1+1]:=1:
mean1rec:=MixedMomentRecurrence(S,n,N,power):
mean1:=evalf(RecToLi(mean1rec[1],mean1rec[2],M,n,N)[-1]/numtrees):
power[a1+1]:=0:
power[a2+1]:=1:
mean2rec:=MixedMomentRecurrence(S,n,N,power):
mean2:=evalf(RecToLi(mean2rec[1],mean2rec[2],M,n,N)[-1]/numtrees):
moments:=[]:
power[a1+1]:=0:
power[a2+1]:=0:
for r from 0 to p1 do
 moments:=[op(moments),[]]:
 for l from 0 to p2 do
  mmrec:=MixedMomentRecurrence(S,n,N,power):
  moments[r+1]:=[op(moments[r+1]),RecToLi(mmrec[1],mmrec[2],M,n,N)[-1]]:
  power[a2+1]:=power[a2+1]+1:
 od:
 power[a1+1]:=power[a1+1]+1:
 power[a2+1]:=power[a2+1]-(p2+1):
od:
newmm:=add(add(binomial(p1,r)*binomial(p2,l)*(-1)^(r+l)*mean1^r*mean2^l*moments[p1-r+1][p2-l+1]*numtrees^(-1),l=0..p2),r=0..p1):
end:

ScaledMixedMomentAboutMean:=proc(S,M,a1,a2,p1,p2) local m,m2:
m:=MixedMomentAboutMean(S,M,a1,a2,p1,p2):
m2:=[MixedMomentAboutMean(S,M,a1,a2,2,0),MixedMomentAboutMean(S,M,a1,a2,0,2)]:
m/(m2[1]^(p1/2)*m2[2]^(p2/2)):
end:

Correlation:=proc(S,M,a1,a2) local cov, v1, v2:
cov:=MixedMomentAboutMean(S,M,a1,a2,1,1):
v1:=MixedMomentAboutMean(S,M,a1,a2,2,0):
v2:=MixedMomentAboutMean(S,M,a1,a2,0,2):
cov/sqrt(v1*v2):
end:


f:=proc(x,y,c):
exp(-x^2/2-y^2/2+c*x*y)/(2*Pi/sqrt(1-c^2)):
end:

Mom:=proc(i,j,c) local x,y:
int(int(x^i*y^j*f(x,y,c),x=-infinity..infinity),y=-infinity..infinity):
end:

MM:=proc(i,j,c)
Mom(i,j,c)/( Mom(2,0,c)^((i+j)/2)):
end:

NormalMoment:=proc(S,M,a1,a2,p1,p2) local cor,c,expr,c1:
cor:=Correlation(S,M,a1,a2):
if abs(cor)>1.01 then
 return("Correlation too high"):
fi:
if cor >1 then
 cor:=1:
elif cor < -1 then
 cor:=-1:
fi:
expr:=simplify(MM(p1,p2,c) assuming csgn(-1/2+1/2*c^2)= -1):
c1:=1/2*(-1+sqrt(4*cor^2+1))/cor:
subs(c=cor,expr):
end:

AsyMeanVar:=proc(S,i) local n,N, power,mo0rec,mo1rec,asymo0,asymo1,toge,num,den,mean,mo2rec,asymo2,var:
power:=[0$(max(S)+1)]:
mo0rec:=MixedMomentRecurrence(S,n,N,power):

power[i+1]:=1:
mo1rec:=MixedMomentRecurrence(S,n,N,power):
asymo0:=AsyC(mo0rec[1],n,N,10,mo0rec[2],1000):
asymo1:=AsyC(mo1rec[1],n,N,10,mo1rec[2],1000):
toge:=expand(asymo1[2]/asymo0[2]):
num:=numer(toge):
den:=denom(toge):
mean:=lcoeff(num)/lcoeff(den)*asymo1[1]/asymo0[1]:

power[i+1]:=2:
mo2rec:=MixedMomentRecurrence(S,n,N,power):
asymo2:=AsyC(mo2rec[1],n,N,5,mo2rec[2],1000):
toge:=expand(asymo2[2]/asymo0[2]):
num:=numer(toge):
den:=denom(toge):
var:=lcoeff(num)/lcoeff(den)*asymo2[1]/asymo0[1]-(mean)^2:
return([evalf(mean),evalf(var)]):
end:

#GeneralInitTreeCount(S,n) returns b_n where b_n is a polynomial in y[S] with each
#monomial corresponding to a tree on n vertices where the number of vertices of degree i
#is the power of y[i]. It does not use the EKHAD package, and so it is suitable
#for finding initial conditions.
GeneralInitTreeCount:=proc(S,n,y) local poly,parts,s,m,part,S_temp:
option remember:
S_temp:=S union {0}:
if n = 0 then
 RETURN(1):
fi:
if n = 1 then
 return(y[0]):
fi:
poly:=0:
for s in S_temp do
 parts:=Part(s, n-1):
 for part in parts do
  poly:= poly + y[s]*mul(GeneralInitTreeCount(S_temp,m,y), m in part):
 od:
od:
poly:
end:

#Part(k,n) returns all possible partitions of n into k pieces (not in any particular order)
Part:=proc(k,n) local alone,together,s,tog,al,s_alone,s_together:
option remember:
if k=1 then
 return([n]):
fi:
if k>n or k =0 then
 return([]):
fi:
alone:=Part(k-1,n-1):
together:=Part(k,n-1):
s_alone:={seq(seq([op(1..i, al), 1, op(i+1..-1, al)], i =0..nops(al)),al in alone)}:
s_together:={seq(seq([op(1..i-1,tog), tog[i]+1, op(i+1..nops(tog),tog)], i=1..nops(tog)), tog in together)}:
[op(s_alone union s_together)]:
end:

#Story1 considers trees on M vertices with child counts in S. It writes a section of a paper
#analyzing these trees; it finds the expectation and variance of the number of nodes with 
#each possible number of children. For each pair of possible child numbers, it outputs the
#correlation as well as the scaled mixed moment where each number is raised to the power p1 
#and p2 in turn
Story1:=proc(S,M,p1,p2) local power, num_trees, averages, variances,s, mixedmoments, normalmoments, corrs, newmm, newnm, newc,S1, s1,i:
power:=[0$(max(S)+1)]:
num_trees:=MixedMomentCount(S,M,power)[M]:
if num_trees = 0 then
 return():
fi:
averages:=[]:
variances:=[]:
for s in S do
 power[s+1]:=1:
 averages:=[op(averages), 1.0*MixedMomentCount(S,M,power)[M]/num_trees]:
 power[s+1]:=0:
 variances:=[op(variances), StraightMoment(S,M,s,2)]:
od:
mixedmoments:=[]:
normalmoments:=[]:
corrs:=[]:
for s in S do
 newmm:=[]:
 newnm:=[]:
 newc:=[]:
 for s1 in S minus {s} do
  newmm:=[op(newmm), ScaledMixedMomentAboutMean(S,M,s,s1,p1,p2)]:
  newnm:=[op(newnm), NormalMoment(S,M,s,s1,p1,p2)]:
  newc:=[op(newc), Correlation(S,M,s,s1)]:
 od:
 mixedmoments:=[op(mixedmoments),newmm]:
 normalmoments:=[op(normalmoments), newnm]:
 corrs:=[op(corrs), newc]:
od:
S1:=[]:
for i from 1 to nops(S) do
 S1 := [op(S1), S minus {S[i]}]:
od:
i:='i': 
 

print(`This section will consider trees on`, M, `vertices where the number of children of each vertex lies in`, S, `. First, we compute the expectation E[X_i] where X_i is the number of vertices with i children in a tree chosen uniformly at random. We find the following values`, seq(E(X[S[i]]) = averages[i], i=1..nops(S))):
print(`Next, we compute the variance Var(X_i) for each i in S, and obtain`, seq(E(X[S[i]]) = variances[i], i=1..nops(S))):
print(`With these calculations finished, we begin demonstrating the fact that all the X_i are pairwise jointly normally distributed. In order to do so, we calculate`, m[i,j]=E((X[i]-mu[i])^p[1]*(X[j]-mu[j])^p[2]/(Var(X[i])^(p[1]/2)*Var(X[j])^(p[2]/2)))):
print(`After performing this calculation, we compute n_{i,j}, what the value would be if X_i and X_j were jointly distributed according to a bivariate normal distribution with means 0, variances 1, and correlation equal to whatever the correlation between X_i happens to be. If the values are similar that provides evidence that X_i and X_j are, in fact, jointly normal. Doing the calculations, we obtain:`):
print(seq(seq(op([`The empirical moment`, m[S[i],S1[i][j]] = mixedmoments[i][j], `while the corresponding normal moment is`, n[S[i],S1[i][j]] = normalmoments[i][j], `and the correlation`, Cor(X[S[i]],X[S1[i][j]]) = corrs[i][j]]), j=1..nops(S)-1), i=1..nops(S))):

return():
end:

#Story2 considers trees on M vertices with child counts in S. It writes a section of a paper analyzing these trees. For s1, s2 in S, it outputs the average and variance of the number of nodes with s1 and s2 children along with the correlation as well as all the scaled (p1,p2) mixed moments where p1,p2 <=p

Story2:=proc(S,M,s1,s2,po) local power, num_trees, averages, variances,s, mixedmoments, normalmoments, corrs, newmm, newnm, newc,i,p1,p2,f,enum_rec,n,N, num_trees_all, first_mom_rec, first_mom_rec2, first_moments:
if (nops(S) <= 1 or (not member(0, S)) or S = {0,1}) then
 return(FAIL):
fi:
power:=[0$(max(S)+1)]:
enum_rec:=MixedMomentRecurrence(S,n,N,power)[1]:
power[s1+1]:=1:
first_mom_rec:=MixedMomentRecurrence(S,n,N,power)[1]:
power[s2+1]:=1:
power[s1+1]:=0:
first_mom_rec2:=MixedMomentRecurrence(S,n,N,power)[1]:
power[s2+1]:=0:
num_trees_all:=MixedMomentCount(S,max(M,30),power):
num_trees:=num_trees_all[M]:
if num_trees = 0 then
 return():
fi:
averages:=[]:
variances:=[]:
first_moments:=[]:
for i from 1 to 2 do 
 s := [s1,s2][i]:
 power[s+1]:=1:
 first_moments:=[op(first_moments), MixedMomentCount(S,max(M,30),power)]:
 averages:=[op(averages), 1.0*first_moments[i][M]/num_trees]:
 power[s+1]:=0:
 variances:=[op(variances), StraightMoment(S,M,s,2)]:
od:
mixedmoments:=[]:
normalmoments:=[]:
for p1 from 1 to po do
 newmm:=[]:
 newnm:=[]:
 for p2 from 1 to po do
  newmm:=[op(newmm), ScaledMixedMomentAboutMean(S,M,s1,s2,p1,p2)]:
  newnm:=[op(newnm), NormalMoment(S,M,s1,s2,p1,p2)]:
 od:
 mixedmoments:=[op(mixedmoments),newmm]:
 normalmoments:=[op(normalmoments), newnm]:
od:
s:='s':
i:='i':
p1:='p1':
p2:='p2':

printf(`This section will consider trees on %d vertices where the number of children of each vertex lies in %a. `, M, S):

printf(`These trees have a generating function given by the functional equation %a; the coefficient on the %a term is the number of trees with %a vertices. `, f(x) = add(x*f(x)^s, s in S), x^n, n):

printf(`To extract the relevant coefficient, we can use the Lagrange inversion theorem, which states the the coefficient of %a in %a is the same as the coefficient of %a in the expression %a. `, x^n, f, z^(n-1), 1/n*(add(z^s, s in S))^n):

printf(`Finally, we can use the Almkvist-Zeilberger algorithm to find that these coefficients follow the recurrence: %a. `, enum_rec):

printf(`Using this recurrence we can generate the first 30 terms of the enumeration sequence: \n\n %a \n\n`, num_trees_all[1..30]):

printf(`We are interested in more than just the number of such trees, though. We are investigating the distribution of the statistics X_%d and X_%d, the numbers of vertices with %d and %d children respectively in a tree chosen uniformly at random. Therefore, we must generalize %a to a function of %a which satisfies the similar functional equation %a. \n\n`, s1,s2,s1,s2,f, [x,seq(y[s], s in S)], f(x,seq(y[s],s in S)) = x*add(y[s]*f(x,seq(y[s],s in S))^s, s in S)):

printf(`Our next step is to compute E[X_%d] and E[X_%d]. We do this by finding the total number of vertices with %d and %d children over all trees with %d vertices. Earlier, we used the Lagrange inversion theorem to extract the coefficients of %a, but this time we will first apply the operator %a to multiply every monomial by the number of nodes with %d children in the corresponding tree, and similarly for %d. `, s1, s2, s1, s2, M, %f, y[s1]*d/dy[s1], s1, s2):

printf(`As before, we combine Lagrange inversion with the Almkvist-Zeilberger algorithm to find the following recurrence for the total number of vertices with %d children over all trees with %a vertices: %a. Then, this recurrence gives the following sequence of vertex counts: \n\n %a \n\n`, s1, n, first_mom_rec, first_moments[1][1..30]):

printf(`Doing the same for vertices with %d children, we obtain the recurrence %a and the sequence \n\n %a \n\n`, s2, first_mom_rec2, first_moments[2][1..30]):

printf(`Armed with this information, we can compute both the expectation E[X_%d] = %f and E[X_%d] = %f. \n\n`, s1, averages[1], s2,  averages[2]):

printf(`In order to calculate the variance of X_%d and X_%d, we will need to calculate the second moments of these two variables, rather than just the first moment. Fortunately, this is simply a matter of applying the operator %a twice rather than once. Doing so gives Var(X_%d)= %f, and making a similar calculation for %d gives Var(X_%d)= %f. \n\n`, s1,s2, y[s1]*d/dy[s1], s1,variances[1],s2, s2, variances[2]):

printf(`With these calculations finished, we begin demonstrating the fact that X_%d and X_%d are jointly normally distributed. In order to do so, we calculate %a for all %a,%a. `, s1,s2, E((X[s1]-averages[1])^p1*(X[s2]-averages[2])^p2/(variances[1]^(p1/2)*variances[2]^(p2/2))), 1<=p1,p2 <= p):

printf(`This calculation proceeds much as the previous ones; we begin by applying the operators %a %a times and %a %a times. We use the Lagrange inversion theorem and Almkvist-Zeilberger algorithm to find a recurrence for the coefficient of %a, find the %d th term, and plug in 1 for all y. `, y[s1]*d/dy[s1], p1,y[s2]*d/dy[s2], p2, x^n, M):

printf(`Since our goal is to determine whether X_%d and X_%d are asymptotically normally distributed, we need to compare these moments to what the moments would be if they were normal. If they were normal, then X_%d and X_%d would have the pdf %a where c is the correlation between X_%d and X_%d. We can compute this correlation, and then compute moments based on this pdf. \n\n`, s1, s2, s1, s2, f(x,y)=1/2*e^(-1/2*x^2-1/2*y^2+c*x*y)*sqrt(1-c^2)/Pi, s1,s2):

printf(`After computing both the actual moments and what the moments would be if X_%d and X_%d were jointly normal, we find the following results: \n\n`, s1,s2):

for p1 from 1 to po do
 for p2 from 1 to po do
  printf(`The (%d,%d) scaled moment about the mean for X_%d and X_%d is %f, while the moment would be %f if the random variables were normally distributed. \n\n`,p1,p2,s1,s2,mixedmoments[p1][p2], normalmoments[p1][p2]):
 od:
od:

#print(`With these calculations finished, we begin demonstrating the fact that`, X[s1], ` and`, X[s2], `are jointly normally distributed. In order to do so, we calculate`, m[p[1],p[2]]=E((X[s1]-averages[1])^p[1]*(X[s2]-averages[2])^p[2]/(Var(X[i])^(p[1]/2)*Var(X[j])^(p[2]/2))),`for all`, 1<=p[1],p[2] <= p):

#print(`After performing this calculation, we compute n_{i,j}, what the value would be if`, X[s1],` and`, X[s2],` were jointly distributed according to a bivariate normal distribution with means 0, variances 1, and correlation equal to whatever the correlation between the two X_i happens to be. If the values are similar, that provides evidence that X_i and X_j are, in fact, jointly normal. Doing the calculations, we obtain:`):
#print(seq(seq(op([`The empirical moment`, m[i,j] = mixedmoments[i][j], `while the corresponding normal moment is`, n[i,j] = normalmoments[i][j]]), j=1..po), i=1..po)):

return():
end:


Story:=proc(S,M,po) local s1, s, S1:
printf(`On Statistics of Ordered Trees Where the Number of Children a Node Can Have is Restricted to %a. \n\n`, S):
printf(`By Lumpy (AKA Yonah's Computer) \n\n`):

for s in S do
 S1:= S minus {s}:
 for s1 in S1 do
  if s < s1 then
   Story2(S,M,s,s1,po):
  fi:
 od:
od:
return():
end: