#OK to post homework
#Blair Seidler, 2021-02-28 Assignment 10

with(combinat):

Help:=proc():
print(`ExactGFinv(n,x), ExactGFmaj(n,x), ExactGFdes(n,x)`):
print(`EstGFinv(n,x,K) ,EstGFmaj(n,x,K) , EstGFdes(n,x,K)`):
print(`PnSeq(N,n)`):
print(`Gf(X,x), inv(L), desc(L), major(L)`):
end:

# 1. I used several utility functions to test things out, but I did not code
#    any of the functions in this exercise the direct way (i.e. by looking at
#    all of permute(n) and counting)

#ExactGFinv(n,x): inputs a positive integer n and a variable name n, and outputs the generating
#function whose coefficient of x^k is the exact number of permutations with inv=k
ExactGFinv:=proc(n,x) option remember:
if not (type(n,integer) and n>=0) then RETURN(FAIL): fi:
if n=1 then RETURN(1): fi:

expand(sum(x^i,i=0..n-1)*ExactGFinv(n-1,x)):

#I think that the reason this works is that I can think about inserting the nth element
#anywhere in the previous permutations. If I insert it at the beginning, it adds n-1 inversions.
#In position 2, it adds n-2 inversions ... at the end, it adds 0 inversions.

#So this formula essentially inserts the last element into all positions of all permutations
#of n-1 and counts the number of inversions in parallel.

#This may be the first moment in my entire life when I truly believe that a generating function
#made my life significantly easier!!!
end:

#Ndesc(n,k): outputs the number of permutations of n with exactly k descents
Ndesc:=proc(n,k) option remember:
if not (type(n,integer) and n>=0 and type(k,integer)) then RETURN(FAIL): fi:
if n=0 or k<0 or k>=n then RETURN(0): fi:
if k=0 or k=n-1 or n=1 then RETURN(1): fi:
(k+1)*Ndesc(n-1,k)+(n-k)*Ndesc(n-1,k-1):
end:

#ExactGFdes(n,x): inputs a positive integer n and a variable name n, and outputs the generating
#function whose coefficient of x^k is the exact number of permutations with des=k
ExactGFdes:=proc(n,x) option remember:
if not (type(n,integer) and n>=0) then RETURN(FAIL): fi:
if n=1 then RETURN(1): fi:

add(Ndesc(n,k)*x^k,k=0..n):

end:

#Note: I started playing around with major indices after midnight. I kept getting the
#same numbers as I did for inversions. I thought I was going crazy, so I looked on
#Wikipedia and found that MacMahon proved this result in 1913. So I'm not going to
#write another function for it.
#I am going to think about a combinatorial proof at some point before I look it up.
ExactGFmaj:=proc(n,x): ExactGFinv(n,x): end:

# 2. By using randperm(n), write procedures

#EstGFinv(n,x,K): samples K random permutations and finds the empircal pgf of inversions
EstGFinv:=proc(n,x,K):
Pgf([seq(inv(randperm(n)),i=1..K)],x):
end:

#EstGFmaj(n,x,K): samples K random permutations and finds the empircal pgf of major indices
EstGFmaj:=proc(n,x,K):
Pgf([seq(major(randperm(n)),i=1..K)],x):
end:

#EstGFdes(n,x,K): samples K random permutations and finds the empircal pgf of descents
EstGFdes:=proc(n,x,K):
Pgf([seq(desc(randperm(n)),i=1..K)],x):

end:

#Compare plotDist and SAc for up to the 8th moment, with n=100, and K=2000.
#I actually ran my trials with K=100000, because it was fast enough. Both the
#distributions and the moments are very close to the exact values.

# 3. 

#PnSeq(N,t): inputs a positive integer N and outputs the first N terms of Pn(t)
PnSeq:=proc(N,t) local f:
f := taylor(exp(t*(exp(x)-1)), x = 0, N+1):

[seq(i!* coeff(f,x,i),i=0..N)]:

end:

#(i) Using SAc, do they seem to be asymptotically normal?

# I looked at PlotDist, and they certainly seem to be asymptotically normal.

#(ii) Does the average number of parts and the variance divided by n tend to some constants?

#I tried to do some big runs, but they were taking too long to complete. I may try again later.

# 4. Using the procedure to generate a random set partition of {1,...,n} from a previous class,
#    write estimated versions and see how close they get to the exact values, for say, n=100.

EstPnSeq:=proc(n,t,K) local i,f,r:
#f:=0:
#for i from 1 to K do
#  r:=RandSPn(n):
#  f:=f+t^nops(r):
#od:
#f:
add(t^nops(RandSPn(n)),i=1..K):
end:

#I ran this for n=100 and K=20000, and the mean, variance, and moments 3 through 8 are pretty close.

#### Utility functions ####

#Gf(X,x): inputs a RV X (given as list of values where the value of
#Mr. or Ms. i is L[i]) outputs the ordinary generating function of X
#using the variable
Gf:=proc(X,x) local i:add(x^X[i],i=1..nops(X)):end:


inv:=proc(L) local i,j,s:
s:=0:
for i from 1 to nops(L)-1 do
for j from i+1 to nops(L) do
if L[i]>L[j] then s:=s+1: fi:
od: od:
s:
end:

desc:=proc(L) local i,s:
s:=0:
for i from 1 to nops(L)-1 do
if L[i]>L[i+1] then s:=s+1: fi:
od:
s:
end:

major:=proc(L) local i,s:
s:=0:
for i from 1 to nops(L)-1 do
if L[i]>L[i+1] then s:=s+i: fi:
od:
s:
end:

#### From C10.txt #####
#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,f1,L:
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:

L:=[seq([(i-mu)/sig,coeff(f/subs(x=1,f),x,i)*sig],i=trunc(mu-K*sig)..trunc(mu+K*sig))]:

plot(L):

end:

#Pgf(X,x): inputs a RV X (given as list of values where the value of
#Mr. or Ms. i is L[i]) outputs the prob. generating function of X
#using the variable
Pgf:=proc(X,x) local i:add(x^X[i],i=1..nops(X))/nops(X):end:


#### From Homework 9 ####
#f:=Sum(x^a[i], i=1..n);
#diff(f,x)= subs(x=1,Sum(a[i]*x^(i-1),i=1..n))/n;
#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:

#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:

#### From Homework 7 ####
#RandSPn(n): inputs a positive integer n and outputs a (uniformly) at random set
#partition of {1,...,n} using LD(L) and RandSPnk(n,k)
RandSPn:=proc(n) local c,k,L:
L:=[]:
for k from 1 to n do
  L:=[op(L),NuSPnk(n,k)]:
od:
c:=LD(L):
RandSPnk(n,c):
end:
#SPnk(n,k): Inputs pos. integers k and n and 1<=k<=n
#outputs the SET of set-partitions of {1,...,n} with exactly
#k parts
SPnk:=proc(n,k) local SP,SP1,sp,i:
option remember:
if n=1 then
  if k=1 then 
   RETURN({{{1}}}):
  else
    RETURN({}):
  fi:
fi:

#Case (i): n is a singleton
SP1:=SPnk(n-1,k-1):

SP:={seq(sp union {{n}},sp in SP1)}:

#Case (ii) n has friends
SP1:=SPnk(n-1,k):

for sp in SP1 do
 for i from 1 to nops(sp) do
  SP:=SP union {{op(1..i-1,sp), sp[i] union {n},op(i+1..nops(sp),sp)}}:
 od:
 od: 
SP:
end:  

#SPn(n): The set of all set-partitions of {1,...,n}
#(same as combinat[setpartition](n))
SPn:=proc(n) local k:  { seq(  op(SPnk(n,k)),k=1..n)}:end:

#NuSPnk(n,k): The number of set-partitions of {1,...,n}
#into exactly k parts (aka Strirling numbers of the second kind)
NuSPnk:=proc(n,k) option remember:

if n=1 then
  if k=1 then
   RETURN(1):
  else
   RETURN(0):
 fi:
fi:
NuSPnk(n-1,k-1)+ k*NuSPnk(n-1,k):
end:

NuSPn:=proc(n) local k:
option remember:
add(NuSPnk(n,k),k=1..n):
end:

#RandSPnk(n,k): a random set partition of {1,..,n} with k parts
RandSPnk:=proc(n,k) local c,SP1,k1:
if n=1 then
 if k=1 then
   RETURN({{1}}):
  else
   RETURN(FAIL):
 fi:
fi:
c:=LD([NuSPnk(n-1,k-1), k*NuSPnk(n-1,k)]):

if c=1 then
 SP1:=RandSPnk(n-1,k-1):
 RETURN(SP1 union {{n}}):
 else
 SP1:=RandSPnk(n-1,k):
  k1:=rand(1..k)():
  RETURN( { op(1..k1-1,SP1), SP1[k1] union {n}, op(k1+1..k,SP1)}):
  
fi:
end:


#### From C4.txt ####
#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:
#END FROM C4.txt