#C12.txt: Maple Code for Lecture 12 of Dr. Z.'s (Doron Zeilberger) Experimental Mathematics class, Spring 2021. Delivered March 1, 2021

Help12:=proc(): print(` RandPnk(n,k), RandPn(n), pnSeqGFt(N,t)`): 
print(`TSprintf() , Ferr(L), Conj(L)  `):
end:

#Ferr(L): The Ferrers graph of the partition L
Ferr:=proc(L) local i,k:

k:=max(L):
for i from 1 to nops(L) do
print(TSprintf("*"$L[i], "  "$(k-L[i]))):

od:

end:

#Conj(L):Conjugate partition
Conj:=proc(L) local L1,k,i:

if L=[] then
  RETURN([]):
fi:

k:=nops(L):

L1:=L-[1$k]:

for i from 1 to nops(L1) while L1[i]<>0 do  od:
L1:=[op(1..i-1,L1)]:
[k, op(Conj(L1))]:

end:


#Nice print (found by Victoria Chayes)
TSprintf:= proc()
uses T= Typesetting;
     T:-mrow(seq(`if`(e::string, T:-mn(e), T:-Typeset(T:-EV(e))), e= [args]))
end proc:

#Start from Yuxuan Yang's hw11
#RandPnk(n,k): a uniformly at random INTEGER partition of
#n with largest part k
RandPnk:=proc(n,k) local die1, k1,k2:
if n=1 then 
RETURN([1]):
fi:
if n=k then 
RETURN([k]):
fi:

die1:=[seq(pnk(n-k,k1),k1=1..min(n-k,k))]:
k2:=LD(die1):
[k, op(RandPnk(n-k,k2))]:

end:

#RandPn(n): A random partition of n
RandPn:=proc(n) local k,i,die1,k1:

die1:=[seq(pnk(n,k),k=1..n)]:
k1:=LD(die1):
RandPnk(n,k1):
end:


EstStatAnalOfLargestPart:=proc(n,K,N) local i, nu:
nu:=[]:
for i from 1 to K do 
 nu:=[op(nu),nops(RandPn(n))]:
od:
SA(nu,N):
end:


#EstStatAnalOfLargestPart(100,10000,6)
#[21.77990000, 67.17725599, 1.144286217, 5.111565874, 15.86000751, 66.25212843]
#[21.63340000, 64.72900444, 1.099627897, 4.968728406, 15.48344264, 67.08132648]

#EstStatAnalOfLargestPart(200,10000,6)
#[34.37140000, 147.3746620, 1.174089079, 5.396169297, 18.52936230, 87.70007364]


#pnSeqGFt(N,t): the first N+1 weight-enumetors according to the r.v.
#"number of parts" 
pnSeqGFt:=proc(N,t) local q,f,i,ta:
f:=mul(1/(1-t*q^i),i=1..N+1):
ta:=taylor(f,q,N+1):
[seq(expand(coeff(ta,q,i)),i=0..N)]:
end:

#END from Yuxuan Yang's hw11

#C11.txt: Maple Code for Lecture 11 of Dr.Z's Experimental Mathematics class, Spring 2021, taught ar Rutgers University

Help11:=proc(): print(`PnClass(n), Pnk(n,k), Pn(n),pnk(n,k), pn(n) `): 
print(`pnSeq(N), pnSeqGF(N) `):
end:

#PnClass(n): THE SET of integer partitions of n
PnClass:=proc(n) local n1,L,L1,j:
option remember:
if n=0 then 
 RETURN({[]}):
elif n<0 then
 RETURN({}):
else
L:={}:

  for n1 from 1 to n do
   L1:=PnClass(n-n1):
   for j from 1 to nops(L1) do
     L:= L union {sort([op(L1[j]),n1],`>`)}:
   od:
 od:
   RETURN(L):
fi:
end:


#Pnk(n,k): The LIST of integer partitions of n with largest part k
Pnk:=proc(n,k) local k1,L,L1:
option remember:

if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN([]):
fi:
 
if k>n then
 RETURN([]):
fi:

if k=n then
 RETURN([[n] ]):
fi:

L:=[]:

for k1 from min(n-k,k) by -1 to 1 do
 L1:=Pnk(n-k,k1):
 L:=[op(L), seq([k, op(L1[j])],j=1..nops(L1))]:
od:

L:

end:

#Pn(n): The list of integer partitions of n in rev. lex. order
Pn:=proc(n) local k:option remember:[seq(op(Pnk(n,n-k+1)),k=1..n)]:end:

#pnk(n,k):
pnk:=proc(n,k) local k1:
option remember:
if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN(0):
fi:
if n=k then
 RETURN(1):
else
 RETURN(add(pnk(n-k,k1),k1=1..min(k,n-k))):
fi:
end:

pn:=proc(n) local k: option remember: 
if n=0 then
 RETURN(1):
fi:

add(pnk(n,k),k=1..n):end:

pnSeq:=proc(N) local n: [seq(pn(n),n=0..N)]:end:

#pnSeqGF(N):
pnSeqGF:=proc(N) local i,q,f:
f:=taylor(mul(1/(1-q^i),i=1..N),q=0, N+1 ):
[seq(coeff(f,q,i),i=0..N)]:
end:

#From C10.txt
#C10.txt: Maple code for Lecture 10
Help10:=proc(): print(` SAc(f,x,K), LD(L) ,PlotDist(f,x,K) , OneTrial(p,n)`): 
print(`ManyTrials(p,n,K,x)`):
end:

#Taken from Blair Seidler's homework (who adapted a previous
#code from Dr. Z.)
#SAc(f,x,K): The SA(X,x,K) of the rv X such whose Pgf is f of x
SAc:=proc(f,x,K) local i,mu,sig2,M,f1:
if subs(x=1,f)=0 then
 RETURN(0):
fi:

f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
f1:=f1/x^mu:
f1:=x*diff(f1,x):
f1:=x*diff(f1,x):
sig2:=subs(x=1,f1):

M:=[mu,sig2]:

if sig2=0 then RETURN(M): fi:
for i from 3 to K do
  f1:=x*diff(f1,x):
#  M:=[op(M), evalf(subs(x=1,f1)/sig2^(i/2))]:
  M:=[op(M), subs(x=1,f1)/sig2^(i/2)]:
od:
M:
end:

#LD(L): Inputs a list of positive integers L (of n:=nops(L) members)
#outputs an integer i from 1 to n with the prob. of i being
#proportional to L[i]
#For example LD([1,2,3]) should output 1 with prob. 1/6
#output 2 with prob. 1/3
#output 3 with prob. 3/6=1/2
LD:=proc(L) local n,i,su,r:
n:=nops(L):
r:=rand(1..convert(L,`+`))():
su:=0:

for i from 1 to n do
 su:=su+L[i]:
  if r<=su then
    RETURN(i):
  fi:
od:
end:

#PlotDist(f,x,K): Inputs a prob. generating function of
#some r.v. X outputs a graph of its SCALED distribution
#A graph of the pdf of (X-mu)/sig. It draws it
#from -K*sig to K*sig

PlotDist:=proc(f,x,K) local i,mu,sig,M,f1:
if subs(x=1,f)=0 then
 RETURN(0):
fi:

f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
f1:=f1/x^mu:
f1:=x*diff(f1,x):
f1:=x*diff(f1,x):
sig:=sqrt(subs(x=1,f1)):
if sig=0 then RETURN(FAIL): fi:

#corrected after class, following Blair's suggestion
plot([seq([(i-mu)/sig,coeff(f,x,i)/subs(x=1,f)*sig],i=trunc(mu-K*sig)..trunc(mu+K*sig))]):

end:

#OneTrial(p,n): The outcome of ONE trial of tossing a coin n times
#where Pr(Head)=p. Returns the number of Heads
#p=1/3
OneTrial:=proc(p,n) local L,p1:
p1:=convert(p,rational):
L:=[denom(p1)-numer(p1), numer(p1)]:
add(LD(L)-1,i=1..n):
end:

#ManyTrials(p,n,K,x): the EMPIRICAL PGF  in the variable x, of repeating K times
#OneTrial(p,n)
ManyTrials:=proc(p,n,K,x) local i:
add( x^OneTrial(p,n) , i=1..K)/K:
end: