> #ok to post
;
> #Yifan Zhang, 10/14/2020
;
> 
;
> 
;
> #Q1.
;
> #write a procedure
> #CPg(L)
> #that inputs a list of lists L=[L1,L2, ...,Lk] of "databases" (represented as lists L1, L2, ..., Lk) and outputs
> #CP(L1,L2, ..., Lk)
> #In particular, CPg([L1,L2]) should output the same as CP(L1,L2) for any two lists of integers L1 and L2
;
> #
;
> #
;
> #
;
> CP:=proc(L1,L2) local i,j,L:
> L:=[]:
> for i from 1 to nops(L1) do
>  for j from 1 to nops(L2) do
>   L:=[op(L),L1[i]+L2[j]]:
>  od:
> od:
> L:
> end:
> CPg:= proc(L) local num, result, i, j:
> option remember: 
> num := nops(L): result:=[]:
> for i from 1 to (num-1) do
> for j from i+1 to num do 
> result:= [op(result), CP(L[i], L[j])]:
> od:od:
> result:
> 
> end:
> #For example. 
;
> bigList:= [[1,2,3], [3,4,5],[213,24,2],[65,84,2],[24,6,8]];
Typesetting:-mprintslash([(bigList := [[1, 2, 3], [3, 4, 5], [213, 24, 2], [65,
84, 2], [24, 6, 8]])],[[[1, 2, 3], [3, 4, 5], [213, 24, 2], [65, 84, 2], [24, 6
, 8]]])

;
> CPg(bigList);
[[4, 5, 6, 5, 6, 7, 6, 7, 8], [214, 25, 3, 215, 26, 4, 216, 27, 5], [66, 85, 3,
67, 86, 4, 68, 87, 5], [25, 7, 9, 26, 8, 10, 27, 9, 11], [216, 27, 5, 217, 28,
6, 218, 29, 7], [68, 87, 5, 69, 88, 6, 70, 89, 7], [27, 9, 11, 28, 10, 12, 29,
11, 13], [278, 297, 215, 89, 108, 26, 67, 86, 4], [237, 219, 221, 48, 30, 32, 
26, 8, 10], [89, 71, 73, 108, 90, 92, 26, 8, 10]]

;
> nops(CPg(bigList))
10

;
> 
;
> 
;
> 
;
> #Q2. 
;
> #Procedures AveClever(L,x) and kthMomentClever(L,x) inputs a database L, and first finds its weight-enumerator. Suppose that you already know the weight-enumerator, call it f, in the variable x. Modify these procedures to write procedures
> #AveGF(f,x) and kthMomentGF(f,x)
> #that do the same things if you already have the weight-enumerator f (in terms of x).
> #In particular, for any "data-base" L, check that
> #AveGF(WtEn(L,x),x)= AveClever(L,x)
> #and
> #kthMomentGF(WtEn(L,x),x)= kthMomentClever(L,k)
;
> WtEn:=proc(L,x) local i:
> add(x^L[i],i=1..nops(L)):
> end:
> AveClever:=proc(L) local f,x:
> f:=WtEn(L,x):
> subs(x=1,diff(f,x))/subs(x=1,f):
> end:
> AveGF:= proc(f,x) option remember:
> subs(x=1,diff(f,x))/subs(x=1,f):
> end:
> evalb(AveGF(WtEn(L,x),x)=AveClever(L))
true

;
> 
;
> kthMomentClever:=proc(L,k) local x,mu,f,f1,i:
> f:=WtEn(L,x):
> mu:=AveClever(L):
> f1:=f/x^mu:
> for i from 1 to k do
> f1:=expand(x*diff(f1,x)):
> od:
> subs(x=1,f1)/nops(L):
> end:
> 
;
> kthMomentGF:=proc(f,k,x) local r,mu,f1,i:
> mu:= AveGF(f,x):
> f1:=f/x^mu:
> 
> for i from 1 to k do
> f1:= expand(x*diff(f1,x)):
> od:
> subs(x=1,f1):     #At this step, we should not divide by nops()
> end:
> 
> #Debug my procedure. 
;
> f:= WtEn([1,2,4,3],x); f;
Typesetting:-mprintslash([(f := x^4+x^3+x^2+x)],[x^4+x^3+x^2+x])
x^4+x^3+x^2+x

;
> mu:=AveGF(f,x)
Typesetting:-mprintslash([(mu := 5/2)],[5/2])

;
> f1:=f/x^mu; f1
Typesetting:-mprintslash([(f1 := (x^4+x^3+x^2+x)/x^(5/2))],[(x^4+x^3+x^2+x)/x^(
5/2)])
(x^4+x^3+x^2+x)/x^(5/2)

;
> f1:=expand(x*diff(f1,x)): f1
3/2*x^(3/2)+1/2*x^(1/2)-1/2/x^(1/2)-3/2/x^(3/2)

;
> f1:=expand(x*diff(f1,x)): f1
9/4*x^(3/2)+1/4*x^(1/2)+1/4/x^(1/2)+9/4/x^(3/2)

;
> subs(x=1,f1)/nops(f)
5/4

;
> evalb(kthMomentGF(WtEn(L,x),2,x)= kthMomentClever(L,2))
true

;
> 
;
> 
;
> 
;
> #Q3.
;
> #The scaled moment about the mean of a "random variable" X (in other words some numerical attribute on a combinatorial) set, defined on a "sample spaces" is defined by
> #m_k(X)/m_2(X)^(k/2)
> #, where m_k(X) is the k-th moment about the mean.
> #Write a procedure
> #ScaledMomentGF(f,x,k)
> #that compute it, given the weight-enumerator (alias generating-function)
;
> #
;
> #kthMomentGF(f,k,x)
;
> #
;
> ScaledMomentGF:= proc(f,x,k) : local a, b:
> a:= kthMomentGF(f,2,x):
> b:=kthMomentGF(f,k,x):
> b/a^(k/2):
> end:
> 
;
> 
;
> 
;
> 
;
> #Q4. 
;
> #We proved in class that the weight-enumerator of the set of words in the alphabet {0,1} of length n, according to their sum (equivalently, the number of 1-s) is (1+x)n
> #Leaving n as a symbol, use Maple to find explicit expressions for
> #ScaledMomentGF(((1+x)/2)^n,x,k), k=2,3,4,5,6,7,8,9,10
> #Then take limits as n goes to infinity (using the Maple command lim( ..., n=infinity)) Get a sequence of numbers. What is it?
> #Is it in the OEIS?
;
> #
;
> #
;
> 
;
> ScaledMomentGF(((1+x)/2)^n,x,2);  limit(%, n=infinity);
1
1

;
> ScaledMomentGF(((1+x)/2)^n,x,3); limit(%, n=infinity);
0
0

;
> ScaledMomentGF(((1+x)/2)^n,x,4); limit(%, n=infinity);
16*(-1/8*n+3/16*n^2)/n^2
3

;
> ScaledMomentGF(((1+x)/2)^n,x,5); limit(%, n=infinity);
0
0

;
> ScaledMomentGF(((1+x)/2)^n,x,6); limit(%, n=infinity);
64*(1/4*n-15/32*n^2+15/64*n^3)/n^3
15

;
> ScaledMomentGF(((1+x)/2)^n,x,7); limit(%, n=infinity);
0
0

;
> ScaledMomentGF(((1+x)/2)^n,x,8); limit(%, n=infinity);
256*(-17/16*n+147/64*n^2-105/64*n^3+105/256*n^4)/n^4
105

;
> ScaledMomentGF(((1+x)/2)^n,x,9); limit(%, n=infinity);
0
0

;
> ScaledMomentGF(((1+x)/2)^n,x,10); limit(%, n=infinity);
1024*(31/4*n-1185/64*n^2+945/1024*n^5+4095/256*n^3-1575/256*n^4)/n^5
945

;
> 
;
> ScaledMomentGF(((1+x)/2)^n,x,4); limit(%, n=infinity);
16*(-1/8*n+3/16*n^2)/n^2
3

;
> limit(%, n=infinity);
3

;
> #Bell shape 
;
> 
;
> #So I get a sequence, 1,0,3,0,15,0,105,0,945
;