#OK to post homework
#Blair Seidler, 2021-01-24, Assignment 1

with(combinat):

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

# 3.

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

n:=nops(pi):
f:=0:

for i from 1 to n do
 if pi[i]=i then
  f:=f+1:
 fi:
od:
f:

end:


# 4.

#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. The method used here is
#similar to the approach used in class for RandDer(n), which is equivalent in functionality
#to RandAlmostDer(n,0)
RandAlmostDer:=proc(n,k) local pi:

pi:=randperm(n):

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

end:

# 5.

#EstimateAveFP(n,K): that inputs a positive integer n and a large positive integer K, 
#picks random permutations of {1,...,n} (using randperm(n) in the combinat package of Maple)
#and outputs the average number of fixed points.
EstimateAveFP:=proc(n,K) local pi,i,totfp:

totfp:=0:

for i from 1 to K do
 pi:=randperm(n):
 totfp:=totfp+FP(pi):
od:
totfp/K:

end: 

#By taking K fairly large, say K=10000, and various n, can you guess the exact value?

#Running several trials with n=40, 100, or 1000 and K=10000 or 100000 resulted in an average
#number of fixed points which was always very close to 1.

# 6. [Optional challenge]
# Can you find the exact value of the average? Can you prove it mathematically?

#The exact value of the average number of fixed points is 1.
#Proof: We can use linearity of expectation for this problem. Let X(i) be the event
#that a permutation of [n] chosen uniformly at random has fixed point i (with or without
#other fixed points). E(X(i)) can be calculated directly. There are n! total permutations.
#If we fix position i, there are (n-1)! permutations of the other elements. This means that
#E(X(i))=(n-1)!/n!=1/n. Let X be the random variable for the number of fixed points in a
#randomly chosen permutation of [n]. Then X=sum(from 1 to n) X(i), and by linearity of
#expectation E(X)=E(sum(from 1 to n) X(i))=sum(from 1 to n) E(X(i))=sum(from 1 to n) 1/n=1.
#Note that linearity of expectation is valid even if the events are dependent.

# 7.

#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). 
EstimateMomFP:=proc(n,r,K) local pi,i,tot:

tot:=0:
for i from 1 to K do
 pi:=randperm(n):
 tot:=tot+FP(pi)^r:
od:
tot/K:

end:

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

#Running several trials with n=10, 40, 100, or 1000 and K=100000 resulted in an average
#value of FF(pi)^2 of 2. Because I was curious, I ran similar trials for larger values of r which
#seem to indicate that the average value of FF(pi)^3 is 5 and the average value of FF(pi)^4 is 15.

# 8. [Bigger Optional challenge]

#Can you find the exact value of the average of random variable FP(pi)^2 over all permutations of
#length n? Can you prove it mathematically?

#The exact value of the random variable FP(pi)^2 is 2.
#Proof: We can again use linearity of expectation for this problem. Let X(i) be the event
#that a permutation of [n] chosen uniformly at random has fixed point i (with or without
#other fixed points). As in problem 6, we set X=sum(from 1 to n) X(i). Now E(X^2) can be
#calculated directly. E(X^2)=E((sum X(i))^2)=E(sum (over i) X(i)^2 + sum (over i=/=j) X(i)X(j))
# =n (1/n) + (n(n-1)) (1/(n(n-1)))=1+1=2.

#At this point, I became curious about the experimentally observed values E(FP(pi)^3)=5 and
#E(FP(pi)^4)=15. The technique for computing FP(pi)^2 generalizes to FP(pi)^r as follows: find all
#partitions of r. For each partition, multiply the multinomial coeeficient for the partition by the
#number of distinct permutations of the partition, then divide by the number of permutations of the
#partition. Summing those value for all partitions gives the exact value of E(FP(pi)^r). I was able
#to confirm the experimental values for r=3 and r=4. I also used this method to calculate that
#FP(pi)^5=52 and FP(pi)^6=203, and experimental results seem to confirm those values as well.