# OK to post homework
# Lucy Martinez, 04-09-2026, Assignment 20

with(combinat):


# Question 1:
# Lucy Martinez asked a very good question:
#  In one bootstrap sample with n data entries,
#   why do we sample (with replacement) exactly n times. 
# Write a procedure GenBSsamp(L,k) that generates a bootstrap sample
#  of length k*nops(L), and generalize BSse(L,B) to compute GenBSse(L,B,k)
# With L:=CH(),
#  What is GenBSse(L,1000,1/2)?
#   ANSWER: It is approximately 3.713292516
#  What is GenBSse(L,1000,2)?
#   ANSWER: It is approximately 1.839265103
#  What is GenBSse(L,1000,1)? (alias BSse(L,1000)
#   ANSWER: It is approximately 2.656733622
#  Which of them is closest to SE(L)? 
#   ANSWER: SE(L)=2.635231383 so GenBSse(L,1000,1) is the closest one.
#    Wow! We definitely want to sample exactly n times!
GenBSsamp:=proc(L,k) local n,ra,L1,i:
n:=nops(L):
ra:=rand(1..n):
L1:=[seq(L[ra()],i=1..k*n)]:
L1:
end:


GenBSse:=proc(L,B,k) local M,i,mu:
M:=[seq(Av(GenBSsamp(L,k)),i=1..B)]: 

#here we estimate the standard-error
mu:=Av(M):
evalf(sqrt(add((M[i]-mu)^2,i=1..B)/(B-1))):
end:



# Question 2:
# Recall (from a previous class) that the number of inversions
#  of a permutation pi, denoted by inv(pi), is the number of pairs
#  1 <= i < j <= n such that pi[i] > pi[j].
# Another permutation statistic is the major index,
#  defined as the sum of the places 1 <= i <= n-1 such that pi[i] > pi[i+1].
# By picking 1000 random permutations of length 100,
#  estimate the average of inv(pi)*maj(pi) over all permutations of length 100
#  (BTW, I have no clue what it is supposed to be, there is no explicit formula
#   (that I know of). All I know is that it is between 0 and binomial(100,2)^2 )
# What is the estimated standard error?
# ANSWER:
# Some averages followed by the standard error are the following:
# 6107450.335, 19335.35735
# 6121392.480, 18904.93446
# 6108728.513, 19123.24171
# 6142230.834, 18492.90604

AvSEInvMaj:=proc(n,K) local L,i,pi:
L:=[]:
for i from 1 to K do
 pi:=randperm(n):
 L:=[op(L),invmaj(pi)]:
od:
Av(L), SE(L):
end:


invmaj:=proc(pi):
inv(pi)*maj(pi):
end:

#maj(pi): The major index : The sum of the places where
#pi[i]>pi[i+1]
maj:=proc(pi) local n,i,co:
co:=0:
n:=nops(pi):
for i from 1 to n-1 do
 if pi[i]>pi[i+1] then  
   co:=co+i:
 fi:
od:
co:
end:


#inv(pi): The number of inversions of the permutation pi
#For example inv([1,2,3])=0, inv([3,2,1])=3
inv:=proc(pi) local n,i,j,co:
n:=nops(pi):
co:=0:
for i from 1 to n do 
 for j from i+1 to n do   
  if pi[i]>pi[j] then 
   co:=co+1:
  fi:
 od:
od:
co:
end:

#########################From previous classes:
# C19.txt, April 6, 2026

Help19:=proc(): 
 print(`Av(L), CH(), SE(L)`):
 print(`BSsamp(L), BSse(L,B)`):
 print(`FP(pi), GFfpBF(n,x), GFfp(n,x), PerS(n,K)`):
end:


# (1+x)^100 models tossing a coin 100 times
# If L:=SMoms((1+x)^n,x,6):
# then 
#  limit(L[4], n=infinity)=6 (central limit theory)
#  limit(L[6], n=infinity)=15 (central limit theory)
#  limit(L[8], n=infinity)=105 (central limit theory)
# This is the normal distribution


#Today: We will learn about the bootstrapping method


#CH(): The list of the heights in Math640 (Sp. 2026) (w/o Pablo and Jike)
CH:=proc():
[61,63,71,68,54,73]:
end:


#Av(L): Given a list L of numbers outputs the average (aka mean or expectation)
Av:=proc(L) local n,i:
n:=nops(L):
evalf(add(L[i],i=1..n)/n):
end:


#SE(L): The classical textbook formula for the standard-error
SE:=proc(L) local n,mu,i:
n:=nops(L):
mu:=Av(L):
#add the squares of the deviations from the average: add((L[i]-mu)^2,i=1..n)
evalf(sqrt(add((L[i]-mu)^2,i=1..n))/n):
end:


#BSsamp(L): ONE bootstrap sample from the data list L
# Then, we can run this many times
BSsamp:=proc(L) local n,ra,L1,i:
n:=nops(L):
ra:=rand(1..n):
L1:=[seq(L[ra()],i=1..n)]:
end:


#BSse(L,B): The Bootstrap estimate of the standard error of the data set L
BSse:=proc(L,B) local M,i,mu:
#Here we are taking B bootstrap samples
# and then evaluating a desired statistic for each bootstrap list (B samples)
M:=[seq(Av(BSsamp(L)),i=1..B)]: 

#here we estimate the standard-error
mu:=Av(M):
evalf(sqrt(add((M[i]-mu)^2,i=1..B)/(B-1))):
end:



#FP(pi): The number of fixed points of the permtutation pi
FP:=proc(pi) local n,i,co:
n:=nops(pi):
co:=0:
for i from 1 to n do
 if pi[i]=i then
  co:=co+1:
 fi:
od:
co:
end:



#GFfpBF(n,x): The weight-enumerator of S_n according to x^FP(pi)
GFfpBF:=proc(n,x) local S,pi:
S:=permute(n):
add(x^FP(pi),pi in S):
end:

#Av(F)=Sum(Av({pi[i]=i},i=1..n)}:
# Probability that pi[i]=i is 1/n and probability that pi[i]<>i is (n-1)/n
#E[{pi[i]=i}]=1*1/n + 0*(n-1)/n= 1/n
# Now, add for each i from 1 to n: 1/n+...+1/n = 1 (n times)


#GFfp: (the actual formula) The weight-enumerator of S_n according to x^FP(pi)
GFfp:=proc(n,x) local k:
expand(add(binomial(n,k)*(x-1)^k*(n-k)!,k=0..n)):
end:



#PerS(n,K): A random sample of K permutations of length n 
# returns the list of the number of fixed points of each sample
PerS:=proc(n,K) local i:
[seq(FP(randperm(n)),i=1..K)]:
end:




######################################################################
# C18.txt, March 30, 2026

Help18:=proc(): 
 print(`OT(L), RandSYT(L), Moms(f,x,r) `):
 print(`WtE(S,f,x), Jike(n,x), SMoms(f,x,r), PlotDist(f,x)`):
end:

#Notes:
#Want to see how many rows does the Tableaux from the RSK give young
# For example if pi:=[1,2,4,3] then RS(pi)[2]=[[1,2,5],[3,4]]
#
#Turns out the number of rows of RS(pi)[1] is equal to the number
# of the largest decreasing subsequence

#Statistics Notes:
# If f=Sum(Pr(X=i), all i) is the probability generating function then
#  the average is defined as
#  Sum(Pr(X=i)*i, all i) 
# The r-th moment:
#  m_r(X)=add((X(s)-av)^r, s in S)/nops(S)
# The scaled moments called alpha is defined as
#  alpha_r=m_r/((s.d)^r) where s.d=standard deviation
#
# alpha_4=m_4/sig^4 is called the Kurtosis
# THE STANDARD DISTRIBUTION:
# m_r=0 if r is odd AND
# m_(2r)= 1*3*5*7*...*(2*r-1)=(2*r)!/2^r*r! (the even ones)


#WtE(S,f,x): the weight-enumerator of the set S according to the statistic f(s)
# WtE(permute(5), pi->nops(RS(pi)[1]),x):
WtE:=proc(S,f,x) local s:
add(x^f(s),s in S):
end:


#Moms(f,x,r): 
Moms:=proc(f,x,r) local mu,L,f1,sig,i:
f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
L:=[mu]:

#g.f of X-mu
# Centralized pgf:
f1:=f1/x^mu:

#f1 derived twice:
f1:=x*diff(f1,x):
f1:=x*diff(f1,x):
sig:=evalf(sqrt(subs(x=1,f1))):
L:=[mu,sig]:

for i from 3 to r do
 f1:=x*diff(f1,x):
 L:=[op(L),subs(x=1,f1)]:
od:
L:
end:



#SMoms: The scaled moments
SMoms:=proc(f,x,r) local L,sig,i,mu:
L:=Moms(f,x,r):
mu:=L[1]:
sig:=L[2]:
#L is the list of moments
# L[i]/sig^i are the scaled moments
[mu,sig,seq(L[i]/sig^i, i=3..r)]:
end:


#PlotDist(f,x): plots the prob. distribution inspired by the weight enumerator f in the variable x after it is
#turned to a prob. distribution (after dividing by subs(x=1,f)
PlotDist:=proc(f,x) local f1,i:
f1:=f/subs(x=1,f):
plot([seq([i,coeff(f1,x,i)],i=ldegree(f1,x)..degree(f1,x))]):
end:


Jike:=proc(n,x) local S,pi:
S:=permute(n) :
add(x^nops(RS(pi)[1]), pi in S):
end:



OT:=proc(L) local n,c:
 n:=convert(L,`+`):
 c:=Cells(L)[rand(1..n)()]:
while HLc(L,c)>1 do
 c:=Hook(L,c)[rand(1..HLc(L,c))()]:
od:
c:
end:


#RandSYT(L): a random SYT of shape L
RandSYT:=proc(L) local k,c,i,L1,n,Y1:
if L=[1] then
  RETURN([[1]]):
fi:

n:=convert(L,`+`):
k:=nops(L):
c:=OT(L):
i:=c[1]:

L1:=[op(1..i-1,L),L[i]-1,op(i+1..k,L)]:
if L1[-1]=0 then
 L1:=[op(1..nops(L1)-1,L1)]:
fi:

Y1:=RandSYT(L1):
if nops(L1)=k then
 [op(1..i-1,Y1),[op(Y1[i]),n],op(i+1..k,Y1)]:
else
 [op(1..k-1,Y1),[n]]:
fi:
end:



##################################from before
# C17.txt, March 23, 2026

Help17:=proc():
 print(`Cells(L), Conj(L), Hook(L,c), HL(L,c), NuSYTc(L)`):
 print(`HLc(L,c), NuSYTcc(L)`):
end:

#S:=SYT([3,3,3]):
# nops(S)=42
#One way that we could get a random Young Tableaux is to call in Maple: 
# S[rand(1..nops(S))()];
#However, SYT(L) is constructing all Young Tableaux - not efficient to do this
#We will set up for next week so we can construct a way to generate random Young Tableaux
#
#Robinson-Frame-Thrall
#Suppose we have placed the following boxes:
# ***
# ***
# ***
#NuSYT(L)=n!/mul(nops(Hook(c)), c a cell in the shape L)
# Let L be of shape:
# ***
# ***
# ***
# Here, n=9 (9 cells)
# To count the hook of any cell: count the the number of 
#  stars in the same column and to the right of it (row)
# So now, the hook length for cell=(1,1) is 5 
#  the hook length for cell (1,2) is 4
#  the hook length for cell (2,2) is 3
# If we do it all, we get
#  9!/(5*4*3*4*3*2*3*2*1)=42


#Cells(L): the set of n (=sum(L)) cells [i,j] in the shape L
Cells:=proc(L) local k,i,j:
k:=nops(L):
{seq(seq([i,j], j=1..L[i]),i=1..k)}:
end:

#Def: Conjugate of a partition: each column (left to right) becomes a row
# from top to bottom
# Example: 
# Say we have the following diagram
# 331
# ***
# ***
# *
# Then the conjugate is 322
# ***
# **
# **

#Conj(L): Given a partition L outputs the conjugate of L
Conj:=proc(L) local k,L1,i1,C1,i:
option remember:

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

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


#Hook(L,c): The set of the cells in the hook corresponding
# to the cell c=[i,j]
# i.e. the set of cells to the right and to the bottom of call
Hook:=proc(L,c) local k,i,j,C,i1,j1,i2:
k:=nops(L):
i:=c[1]:
j:=c[2]:
if not (i>=1 and i<=k) then
 return(FAIL):
fi:

if not (j>=1 and j<=L[i]) then
 return(FAIL):
fi:

C:={seq([i,j1], j1=j..L[i])}:

for i2 from i to k while j<=L[i2] do od:
i2:=i2-1:
C:=C union {seq([i1,j],i1=i..i2)}:
C:
end:


#HL(L,c): the hook-length of the cell c in the shape L
HL:=proc(L,c) local n,C,i:
nops(Hook(L,c)):
end:

#HLc(L,c): the hook-length of the cell c in the shape L
# done cleverly
HLc:=proc(L,c) local i,j,L1:
L1:=Conj(L):
i:=c[1]:
j:=c[2]:
L[i]-j+L1[j]-i+1:
end:


#NuSYTc(L): implementing the Frame-Robinson-Thrall hook lenght formula
NuSYTc:=proc(L) local n,C,i,c:
n:=add(L[i],i=1..nops(L)):
C:=Cells(L):
n!/mul(HL(L,c),c in C):
end:


#NuSYTcc(L): implementing the Frame-Robinson-Thrall hook lenght formula
# done cleverly
NuSYTcc:=proc(L) local n,C,i,c:
n:=add(L[i],i=1..nops(L)):
C:=Cells(L):
n!/mul(HLc(L,c),c in C):
end:


#######################################
Help16:=proc():
 print(`RS1(Y,i), RS1left(pi), RS(pi)`):
 print(` `):
end:



#RS1(Y,i): Inputs a partial Young tableaux and another integer i NOT yet in Y
# and places it in the right place, by a bumping process
# it returns a tableaux with one more box followed by the row 
RS1:=proc(Y,i) local k, NewY,lucy,bumpee,i1,j:
if Y=[] then
 return([[i]],1):
fi:
k:=nops(Y):
lucy:=RS11(Y[1],i):
NewY[1]:=lucy[1]:
bumpee:=lucy[2]:

if bumpee=0 then 
 return([NewY[1], op(2..k,Y)],1):
fi:

for i1 from 2 to k while bumpee<>0 do
 lucy:=RS11(Y[i1],bumpee):
 NewY[i1]:=lucy[1]:
 bumpee:=lucy[2]:
 if bumpee=0 then
  return( [seq(NewY[j],j=1..i1), op(i1+1..k,Y)],i1):
 fi:
od:
[seq(NewY[j],j=1..k),[bumpee]],k+1:
end:



#RS1left(pi): inputs a permutation pi (of size nops(pi)) and outputs ONE
# STY of whatever shape (with n boxes)
RS1left:=proc(pi) local Y,i:
Y:=[]:
for i from 1 to nops(pi) do
 Y:=RS1(Y,pi[i])[1]:
od:
Y:
end:



#RS(pi): inputs a permutation pi of ({1,2,...,n:=nops(pi)) and outputs
# a pair of SYT of the SAME shape (with n boxes)
#The Robinson–Schensted algorithm
RS:=proc(pi) local Y1, i,Yr,eaea,p:
if pi=[] then
 return([[],[]]):
fi:
Y1:=[[pi[1]]]:
Yr:=[[1]]:

for i from 2 to nops(pi) do
 eaea:=RS1(Y1,pi[i]):
 Y1:=eaea[1]:
 p:=eaea[2]:
 if p<=nops(Yr) then
  Yr:=[op(1..p-1,Yr),[op(Yr[p]),i],op(p+1..nops(Yr),Yr)]:
 else
  Yr:=[op(Yr),[i]]:
 fi:
od:
[Y1,Yr]:
#Yr is usually what they call Q
end:



############################
# C14.txt, March 09, 2026

Help14:=proc():
 print(`NuSYT(L), SYTpairs(n), NuSYTpairs(n)`):
 print(`RS11(a,i)`):
end:


#NuSYT(L): The NUMBER of Standard Young tableaux of shape L
NuSYT:=proc(L) local i, k,S,L1:
option remember:
k:=nops(L):

if k=0 then
  RETURN(1):
fi:
S:=0:

#now we look for all the legal rows where the element n can be placed
for i from 1 to k-1 do
if L[i]>L[i+1] then
 L1:=[op(1..i-1,L),L[i]-1,op(i+1..k,L)]:
 S:=S+NuSYT(L1): 
fi:
od:

if L[k]>1 then
L1:=[op(1..k-1,L),L[k]-1]:

S:=S+NuSYT(L1):
else
L1:=[op(1..k-1,L)]:
S:=S+NuSYT(L1):
fi:
S:
end:


#Pairs of SYT of the same shape
# We are interested in pairs (Y1,Y2) such that
# shape of Y1 equals shape of Y2
# For n=2:
#  SYT([1,1])={[ [1], [2] ] } so there is only 1 of these STY of shape [1,1]
#
#  SYTpairs(2)= { [[[1, 2]], [[1, 2]]], [ [[1], [2]], [[1], [2]] ]}

#SYTpairs(n):
SYTpairs:=proc(n) local S,P,p,S1,s1,s2:
option remember:
P:=Par(n):
S:={}:
for p in P do
 S1:=SYT(p): #the set of Young Tableaux of the partition p (shape p)
 S:=S union {seq(seq([s1,s2], s1 in S1),s2 in S1)}:
od:
S:
end:

#NuSYTpairs(n):
NuSYTpairs:=proc(n) local S,P, p:
option remember:
P:=Par(n):
S:=0:
for p in P do
S:=S+NuSYT(p)^2:
od:
S:
end:



#RSK: Find a mapping from the set of permutations to the set of STYpairs(n)
#Input: permutation of {1,2,...,n}
#Output: a pair of SYTs of the SAME shape (with n boxes)
#Example: Given [1,4,5,9] then where can we place 10 so that
# it is still increasing? --> put it at the end
# What if I want to place the number 3 using RSK such that if I place it
#  somewhere other than at the end, we must bump the number that was in it
#  So, [1,4,5,9] becomes [1,3,5,9]
# 

#RS11(a,i): inputs an INCREASING list of positive integers a and
# another positive integer i, outputs a pair a1,j where 
# (usually a1 is of the same length as a) and i is placed where it belongs
# and j is the entry that it bumped, unless i is larger than all the members of a
# (i.e. larger than a[-1]) then a1 is [op(a), i], and j is 0
# For example: RS1([11,12,13],4)= [4,12,13], 11
#  where the second output is the element that got bumped out
RS11:=proc(a,i) local k,j:
k:=nops(a):
for j from 1 to k while a[j]<i do
od:
if j=k+1 then
 return([op(a),i],0):
fi:
[op(1..j-1,a),i,op(j+1..k,a)], a[j]:
end:


#######################################
# C13.txt, March 05, 2026

Help13:=proc():
 print(`ParN(n,k), PFG(L), SYT(L), PSYT(Y)`):
end:


#ParN(n,k): outputs the partitions of n for which the first part equal to k
ParN:=proc(n,k) local s,S,T:
S:=Par(n):
T:={}:
for s in S do
 if s[1]=k then
  T:=T union {s}:
 fi:
od:
T:
end:


#PFG(L): The Ferrers graph of the partition L
PFG:=proc(L) local i:
for i from 1 to nops(L) do
 lprint(1$L[i]):
od:
end:



#SYT(L): the set of standard young tableaux of shape L
SYT:=proc(L) local i,s1,k,n,S,L1,S1:
option remember:
k:=nops(L):
n:=add(L[i],i=1..k):

if k=0 then
 return({[]}):
fi:

S:={}:

#now we look for all the legal rows where the element n can be placed
for i from 1 to k-1 do
 if L[i]>L[i+1] then
  L1:=[op(1..i-1,L),L[i]-1,op(i+1..k,L)]:
  S1:=SYT(L1):
  S:=S union {seq([op(1..i-1,s1) , [op(s1[i]),n], op(i+1..k,s1)] ,s1 in S1)}:
 fi:
od:

if L[k]>1 then 
 L1:=[op(1..k-1,L),L[k]-1]:
 S1:=SYT(L1):
 S:=S union {seq([op(1..k-1,s1) , [op(s1[k]),n]],s1 in S1)}:
else
 L1:=[op(1..k-1,L)]:
 S1:=SYT(L1):
 S:=S union {seq([op(1..k-1,s1) , [n] ],s1 in S1)}:
fi:
S:
end:


#PSYT(Y): prints the SYT Y
PSYT:=proc(Y) local i:
for i from 1 to nops(Y) do
 lprint(op(Y[i])):
od:
end:

###################################
# C12.txt, March 02, 2026

Help12:=proc():
 print(`Park(n,k), Par(n)`):
end:

#Park(n,k): The set of partitions of n into exactly k parts
Park:=proc(n,k) local S,k1,S1,s1:
option remember:

if n<k then
   RETURN({}):
fi:

if k=0 then
 if n=0 then
   RETURN({[]}):
   else
    RETURN({}):
fi:
fi:

S:={}:
for k1 from 0 to k do
 S1:=Park(n-k,k1):
 S:= S union {seq([op(s1+[1$k1]),1$(k-k1)], s1 in S1)}:
od:
S:
end:


#Par(n): The set of partitions of n, the clever way.
Par:=proc(n) local k:
{seq(op(Park(n,k)), k=1..n)}:
end: