#OK to post homework
#George Spahn, 2/28/2021, Assignment 10

with(combinat):

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

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

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

#ExactGFinv:=proc(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
ExactGFinv:=proc(n,x) 
add(x^inv(p),p in PermuI(n)):
end:


#ExactGFdes:=proc(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 des=k
ExactGFdes:=proc(n,x) 
add(x^des(p),p in PermuI(n)):
end:

#ExactGFmaj:=proc(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 maj=k
ExactGFmaj:=proc(n,x) 
add(x^maj(p),p in PermuI(n)):
end:

# EstGFinv(n,x,K)  samples K random permutations and finds the empircal 
# prob. gen. function
EstGFinv:=proc(n,x,K) 
evalf(add(x^inv(randperm(n)),i=1..K)/K):
end:

# EstGFdes(n,x,K)  samples K random permutations and finds the empircal 
# prob. gen. function
EstGFdes:=proc(n,x,K) 
evalf(add(x^des(randperm(n)),i=1..K)/K):
end:

# EstGFmaj(n,x,K)  samples K random permutations and finds the empircal 
# prob. gen. function
EstGFmaj:=proc(n,x,K) 
evalf(add(x^maj(randperm(n)),i=1..K)/K):
end:


#PnSeq(n,t) inputs a positive integer n and outputs Pn(t)
PnSeq:=proc(n,t) 
coeff(taylor(exp(t*(exp(x) - 1)), x, n+1),x,n)*n!:
end:

# EstSPn(n,N,K) generates N random SPs of {1..n} and computes the moments of 
# the distribution of the number of parts.
EstSPn:=proc(n,N,K) 
SA([seq(nops(RandSPn(n)),i=1..N)],K):
end:

# It seems to be asymptotically normal

# I suspect the number of parts and the variance divided by n tend to 0.

# exact values for n=100
# [28.62528186, 5.745753494, 0.1250899875, 3.015649197, 1.251569630]

# approx values for n=100
# [28.55200000, 5.877296000, 0.1262598041, 2.946272800, 0.7406444972]


#PermuI(n): the list of all permutations of {1, ...,n} using inertion
PermuI:=proc(n)  local S,i,pi,pi1,L,j,i1:

option remember:
if n=0 then
 RETURN([[]]):
fi:


S:=PermuI(n-1):

L:=[]:

for i from 1 to nops(S) do
pi:=S[i]:

for j from 1 to n do
 pi1:=[op(1..j-1,pi),n,op(j..nops(pi),pi)]:
 L:=[op(L),pi1]:
od:
od:

L:
end:


#RandSPn(n)  Inputs a positive integer n and outputs a (uniformly) at random
# set partition of {1,...,n}
RandSPn:=proc(n) local Die,Die1:

Die1:=[seq(NuSPnk(n,k),k=1..n)]:
Die:=LD(Die1):
RandSPnk(n,Die):
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:

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


#Ave(X): inputs a RV X (given as list of values where the value of
#Mr. or Ms. i is L[i]) outputs the average (alias expectation, alias mean)
Ave:=proc(X) local i:add(X[i],i=1..nops(X))/nops(X):end:

#Mom(X,k): inputs a RV X (given as list of values where the value of
#Mr. or Ms. i is L[i]) outputs the k-th moment about the mean
Mom:=proc(X,k) local i,mu:
mu:=Ave(X):add((X[i]-mu)^k ,i=1..nops(X))/nops(X):
end:

#SA(X,K): inputs a discrete r.v. X (given as a list of numbers) and
#a pos. integer K outputs a list of length K whose
#(i) first entry is the average
#(ii) second entry variance
#The remaining entries are the scaled-third moment (about the mean) aka skeweness
#scaled fourth-moment (about the mean) aka kurtosis and so on
#up to the K-th moment
SA:=proc(X,K) local i,M,mu,sig2:
if not (type(X,list) and type(K,integer) and K>=2) then
 print(`Bad input`):
 RETURN(FAIL):
fi:
mu:=Ave(X,1):
sig2:=Mom(X,2):
M:=[mu,sig2]:

if sig2=0 then
 RETURN(M):
fi:

[op(M), seq(evalf(Mom(X,i)/sig2^(i/2)),i=3..K)]:

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: