#Yuxuan Yang, Jan 23rd, Assignment 1 

Help:=proc(): print(`FP(pi),RandAlmostDer(n,k),EstimateAveFP(n,K),EstimateMomFP(n,r,K)`): end:
with(combinat):

#FP(pi): inputs a permutation pi of {1,...n}, (where n=nops(pi)) and outputs an integer between 0 and n indicating the number of fixed points, i.e. the number of places i such that pi[i]=i. Note that if pi is a derangement, then FP(pi)=0.
FP:=proc(pi) local n,i,j:
n:=nops(pi):
i:=0:
for j from 1 to n do 
 if pi[j]=j then 
  i:=i+1
 fi:
od:
i:
end:

#RandAlmostDer(n,k): inputs a positive integer n, another positive integer k, (between 0 and n) and outputs a random permutation whose number of fixed points is <=k. Note that RandAlmostder(n,0) does the same thing as RandDer(n)
RandAlmostDer:=proc(n,k) local pi:

if not (k>=0 and k<=n) then 
 RETURN(FAIL):
fi:

pi:=randperm(n):

while not FP(pi)<=k do 
pi:=randperm(pi):
od:
pi:
end:

#EstimateAveFP(n,K): inputs a positive integer n, and a large positive integer K, picks randompermutations of {1,...,n} (using randperm(n) in the combinat package of Maple) and outputs the average number of fixed points. By taking K fairly large, say K=10000, and various n, can you guess the exact value?
EstimateAveFP:=proc(n,K) local m,i,pi:
m:=0:
for i from 1 to K do 
pi:=randperm(n):
m:=m+FP(pi):
od:
m/K:
end:
#This average should be about 1.
#With linearity of expectation, each digit gives 1/n. The expection of number of fixed points is 1.

#EstimateMomFP(n,r,K): inputs a positive integer n, another positive integer, r, and a large positive integer K, picks randompermutations of {1,...,n} (using randperm(n) in the combinat package of Maple) and outputs the average number of FP(pi)^r. Note that EstimateMomFP(n,1,K) does the same thing as EstimateAve(n,K). By taking K fairly large, say K=10000, r=2, and various n, can you guess the exact value of the average of FF(pi)^2?
EstimateMomFP:=proc(n,r,K) local m,i,pi:
m:=0:
for i from 1 to K do 
pi:=randperm(n):
m:=m+FP(pi)^r
od:
m/K:
end:
#This 2nd moment should be about 2. Here is the proof.
#Let e_i be the indicator of i-th digit to be fixed, then the 2nd moment is the square of the sum of e_i.
#E(e_i ^2)=E(e_i)=1/n and E(e_i * e_j)=(n-2)!/n!=1/(n*(n-1)), so E(FP(pi)^2)=n*1/n+n(n-1)/(n(n-1))=2.