#OK to post
#Yuxuan Yang, Feb 27th, Assignment 10

with(combinat):

ExactGFinv:=proc(n,x) local i:
if n=0 then 
 RETURN(1):
fi:
add(x^i,i=0..n-1)*ExactGFinv(n-1,x):
end:

#permu with j at last digit
ExactGFdesj:=proc(n,x,j) local i:
option remember:
if n=0 or n=1 then
 RETURN(1):
fi:
expand(add(ExactGFdesj(n-1,x,i) ,i=1..j-1)+add(ExactGFdesj(n-1,x,i) ,i=j..n-1)*x):
end:


ExactGFdes:=proc(n,x) local i:
expand(add(ExactGFdesj(n,x,j),j=1..n)):
end:


ExactGFmajj:=proc(n,x,j) local i:
option remember:
if n=0 or n=1 then
 RETURN(1):
fi:
expand(add(ExactGFmajj(n-1,x,i) ,i=1..j-1)+add(ExactGFmajj(n-1,x,i) ,i=j..n-1)*x^(n-1)):

end:


ExactGFmaj:=proc(n,x) local i:
expand(add(ExactGFmajj(n,x,j),j=1..n)):
end:


EstGFinv:=proc(n,x,K) local i,gf, perm,j,k,inv:
gf:=0:
inv:=0:
for i from 1 to K do
 perm:=randperm(n):
 for j from 1 to n do
  for k from 1 to n do    
    if k<j and perm[k]>perm[j] then
      inv:=inv+1:
    fi:
  od:
 od:
 gf:=gf+x^inv:
 inv:=0:
od:
gf:
end:

EstGFdes:=proc(n,x,K) local i,gf, perm,j,num:
gf:=0:
num:=0:
for i from 1 to K do
 perm:=randperm(n):
 for j from 1 to n-1 do    
    if perm[j]>perm[j+1] then
      num:=num+1:
    fi:
 od:
 gf:=gf+x^num:
 num:=0:
od:
gf:
end:

EstGFmaj:=proc(n,x,K) local i,gf, perm,j,num:
gf:=0:
num:=0:
for i from 1 to K do
 perm:=randperm(n):
 for j from 1 to n-1 do    
    if perm[j]>perm[j+1] then
      num:=num+j:
    fi:
 od:
 gf:=gf+x^num:
 num:=0:
od:
gf:
end:

#3
PnSeq:=proc(N,t,n) local k,f,x,pnt:
f:=exp(t*(exp(x)-1) ):
pnt:=taylor(coeff(taylor(f,x,n+1),x,n)*(n!),t,N):
[seq(coeff(pnt,t,k),k=0..N-1)]:
end:

#the procedure above has some problems, no idea how to fix it and do problem 4.



















#C10.txt: Maple code for Lecture 10
Help10:=proc(): print(` SAc(f,x,K), LD(L) ,PlotDist(f,x,K) , OneTrial(p,n)`): 
print(`ManyTrials(p,n,K,x)`):
end:

#Taken from Blair Seidler's homework (who adapted a previous
#code from Dr. Z.)
#SAc(f,x,K): The SA(X,x,K) of the rv X such whose Pgf is f of x
SAc:=proc(f,x,K) local i,mu,sig2,M,f1:
if subs(x=1,f)=0 then
 RETURN(0):
fi:

f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
f1:=f1/x^mu:
f1:=x*diff(f1,x):
f1:=x*diff(f1,x):
sig2:=subs(x=1,f1):

M:=[mu,sig2]:

if sig2=0 then RETURN(M): fi:
for i from 3 to K do
  f1:=x*diff(f1,x):
#  M:=[op(M), evalf(subs(x=1,f1)/sig2^(i/2))]:
  M:=[op(M), subs(x=1,f1)/sig2^(i/2)]:
od:
M:
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:

#PlotDist(f,x,K): Inputs a prob. generating function of
#some r.v. X outputs a graph of its SCALED distribution
#A graph of the pdf of (X-mu)/sig. It draws it
#from -K*sig to K*sig

PlotDist:=proc(f,x,K) local i,mu,sig,M,f1:
if subs(x=1,f)=0 then
 RETURN(0):
fi:

f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
f1:=f1/x^mu:
f1:=x*diff(f1,x):
f1:=x*diff(f1,x):
sig:=sqrt(subs(x=1,f1)):
if sig=0 then RETURN(FAIL): fi:

#corrected after class, following Blair's suggestion
plot([seq([(i-mu)/sig,coeff(f,x,i)/subs(x=1,f)*sig],i=trunc(mu-K*sig)..trunc(mu+K*sig))]):

end:

#OneTrial(p,n): The outcome of ONE trial of tossing a coin n times
#where Pr(Head)=p. Returns the number of Heads
#p=1/3
OneTrial:=proc(p,n) local L,p1:
p1:=convert(p,rational):
L:=[denom(p1)-numer(p1), numer(p1)]:
add(LD(L)-1,i=1..n):
end:

#ManyTrials(p,n,K,x): the EMPIRICAL PGF  in the variable x, of repeating K times
#OneTrial(p,n)
ManyTrials:=proc(p,n,K,x) local i:
add( x^OneTrial(p,n) , i=1..K)/K:
end: