#OK to post homework
#Wanying Rao, 02/28/2021, Assignment 10


with(combinat):
#1
#ExactGFinv(n,x),ExactGFmaj(n,x),ExactGFdes(n,x):inputs a positive integer n and 
#a variable name x, and outputs the generating function whose coefficient of x^k is 
#the exact number of permutations with inv=k, maj=k, des=k.

Inv := proc(L) local i,n:
n := nops(L):
if n = 0 then
 RETURN(0):
else
 i := 1:
 while L[i]<>n do
  i := i+1:
 od:
 RETURN(n-i+Inv([op(1..i-1,L),op(i+1..n,L)])):
fi: 
end:

Maj := proc(L) local i,n,m:
n := nops(L):
m := 0:
if n = 0 then
 RETURN(0):
else
 for i from 1 to n-1 do
  if L[i]>L[i+1] then m := m+i: fi:
 od:
fi:
m:
end:

Des := proc(L) local i,n,m:
n := nops(L):
m := 0:
if n = 0 then
 RETURN(0):
else
 for i from 1 to n-1 do
  if L[i]>L[i+1] then m := m+1: fi:
 od:
fi:
m:
end:

ExactGFinv := proc(n,x) local p,k,c,C:
C := Array([seq(0,i=0..n*(n-1)/2)]):
for p in permute(n) do
 c := Inv(p):
 C[c+1] := C[c+1]+1:
od:
add(C[j+1]*x^j,j=0..n*(n-1)/2):
end:

ExactGFmaj := proc(n,x) local p,k,c,C:
C := Array([seq(0,i=0..n*(n+1)/2)]):
for p in permute(n) do
 c := Maj(p):
 C[c+1] := C[c+1]+1:
od:
add(C[j+1]*x^j,j=0..n*(n+1)/2):
end:


ExactGFdes := proc(n,x) local p,k,c,C:
C := Array([seq(0,i=0..n)]):
for p in permute(n) do
 c := Des(p):
 C[c+1] := C[c+1]+1:
od:
add(C[j+1]*x^j,j=0..n):
end:


#2
EstGFinv := proc(n,x,K) local p,k,c,C:
C := Array([seq(0,i=0..n*(n-1)/2)]):
for k from 1 to K do
 p := randperm(n):
 c := Inv(p):
 C[c+1] := C[c+1]+1:
od:
add(C[j+1]*x^j,j=0..n*(n-1)/2):
end:

EstGFmaj := proc(n,x,K) local p,k,c,C:
C := Array([seq(0,i=0..n*(n+1)/2)]):
for k from 1 to K do
 p := randperm(n):
 c := Maj(p):
 C[c+1] := C[c+1]+1:
od:
add(C[j+1]*x^j,j=0..n*(n+1)/2):
end:

EstGFdes := proc(n,x,K) local p,k,c,C:
C := [seq(0,i=0..n)]:
for k from 1 to K do
 p := randperm(n):
 c := Des(p):
 C[c+1] := C[c+1]+1:
od:
add(C[j+1]*x^j,j=0..n):
end:

#3
PnSeq := proc(N,t) 
taylor(exp(t*exp(x)-1),x=0,N):
end:

#4
EstSet := proc(n,K) local L,i,s:
L := [seq(0,i=1..n)]:

for i from 1 to K do
 s := nops(RandSPn(n)):
 L(s):= L(s)+1:
od:
L:
end:

#From C10.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:

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

#From hw8
#By using LD(L) and RandSPnk(n,k) write a procedure
#RandSPn(n)
#that inputs a positive integer n and outputs a (uniformly) at random set partition of #{1,...,n}
RandSPn := proc(n) local i,L,r:
L := [seq(NuSPnk(n,i),i=1..n)]:

r := LD(L):
RandSPnk(n,r):
end:

#From hw9
#Corrected: f2 should be x*diff(x*diff(f1/(x^mu),x),x).
SAc := proc(f,x,K) local f1,f2,mu,sig2,m,k:
if subs(x=1,f)=0 then
 RETURN(0):
fi:

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

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

m := [mu,sig2]:
for k from 3 to K do
f2:=x*diff(f2,x):
m := [op(m), subs(x=1,f2)/sig2^(k/2)]:
od:
m:
end: