#Please do not post homework #Guy Adami, 2026-04-01, Assignment 18 #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]: f1:=f1/x^mu: #g.f. of X-mu #Centralized pgf f1:=x*diff(f1,x): f1:=x*diff(f1,x): sig:=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:=proc(f,x,r) local L,sig,i: L:=Moms(f,x,r): sig:=L[2]: [op(1..2,L),seq(L[i]/sig^i , i=3..r)]: end: #------------------------------------------------------------ Conj:=proc(L) local k,C1,i1,L1,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: #_____________________________________ 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: #_____________________________________ #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: #_____________________________________ 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)}: end: #_____________________________________ RandSYT:= proc(L) local n, Y, L1, corners,Tcell,c, C, row, col, i,j: n:= add(i, i in L): L1:= L: Y:= []: #Makes a Young tableux of all zeros to be filled in later for i from 1 to nops(L) do Y := [op(Y), [seq(0, j = 1 .. L[i])]]: od: for i from n by -1 to 1 do corners:= []: #Get all current corners C:= Cells(L1): for c in C do if HLc(L1,c) = 1 then corners:= [op(corners),c]: fi: od: #Pick a random corner Tcell := corners[rand(1..nops(corners))()]; row:= Tcell[1]: col:= Tcell[2]: #Fill in that element into Y Y[row][col]:=i: #Update the shape L1[row]:=L1[row]-1: #remove row if = 0 if L1[row] = 0 then L1:= subsop(row=NULL, L1): fi: od: RETURN(Y); end: #------------------------------------------------------------ #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: #------------------------------------------------------------ #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: #------------------------------------------------------------ #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 put 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 RS11:=proc(a,i) local k,j: k:=nops(a): for j from 1 to k while a[j]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: #_____________________________________ #RS(pi): inputs a permutation pi of ({1, ..., n:=nops(pi)) and outputs a pair of SYT of the SAME shape (with n boxes) #The Robinson-Schenstead algorithm RS:=proc(pi) local Yl, i, Yr,eaea, p: if pi=[] then RETURN([[],[]]): fi: Yl:=[[pi[1]]]: Yr:=[[1]]: for i from 2 to nops(pi) do eaea:=RS1(Yl,pi[i]): Yl:=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: [Yl, Yr]: end: #_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_- #Question 1 nRolls:= (x + x^2 + x^3 + x^4 + x^5 + x^6)^n: #Performing Moms using this generating function: #Moms(nRolls,x,3): --> #mu = (7n)/2: #sig:= (sqrt(105) * sqrt(n) ) / 6 #---- fnorm:= x^(7*n/2) * exp(35*n/24 * (ln(x))^2); evalb(SMoms(nRolls,x,20)=Smoms(fnorm,x,20); #Returned false #Conjecture for average [2,2] # evalf(4 + sum(n*(1/2)^(n - 3), n = 4 .. N + 3)) #This conjecture comes from building a SYT from cell [1,1] similar to RandSYT, but the other direction #After the first 3 are filled because the min value of [2,2] is 4, call it a 50% chance to add to row 1 or add to row 2 #Adding to row 3 can only be done one time and then another added to row 3 before 2 will no longer make a SYT #This value converges to 9 for sufficient N which shouldn't happen (I am pretty sure) #_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_- #Question 2 with(Statistics): AppxAve22:=proc(n,K) local L, i, Cell_22, c : Cell_22:=[]: for i from 1 to K do L:=RandSYT([n,n,n]): c:= L[2][2]: Cell_22:=[op(Cell_22), c]: od: Mean(Cell_22); end: #AppxAve22(30,1000) = 7.55300000000000 #AppxAve22(30,10000) = 7.69790000000000 #These two values are pretty close #These values are close-ish to the conjectured value of 9 #_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_- #Question 3 with(combinat): AppxJike:=proc(n,x,K) local N,k, pi: N:=[]: for k from 1 to K do: N:=[op(N),randperm(n)]: od: WtE(N, pi->nops(RS(pi)[1]), x); end: PlotDist(AppxJike(100,x,1000),x); PlotDist(AppxJike(200,x,1000),x); #These plots have the same shape (normal looking distribution) with n=100 centered ~16-17 whereas n=200 centered ~24