#OK to post homework
#Wanying Rao, 28/02/2021, Assignment 11


with(combinat):
#1
RandPnk := proc(n,k) local k1,L,c:

if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN(FAIL):
fi:
 
if k>n then
 RETURN(FAIL):
fi:

if k=n then
 RETURN([n]):
fi:

L:=[seq(0,i=1..min(n-k,k))]:

for k1 from min(n-k,k) by -1 to 1 do
 L[k1] := pnk(n-k,k1):
od:

c := LD(L):
[k,op(RandPnk(n-k,c))]:
 
end:

RandPn := proc(n) local  L,c,k:

if n=0 then
 RETURN([]):
fi:

L := [seq(pnk(n,k),k=1..n)]:
c := LD(L):

RandPnk(n,c):

end:


#2
EstStatAnalOfLargestPart := proc(n,K,N) local P,k:

P := Array([seq(0,i=1..K)]):

for k from 1 to K do
 P[k] := nops(RandPn(n)):
od:

SA(P,N):
end:

#4
pnSeqGFt := proc(N,t)local i,q,f:
f:=taylor(mul(1/(1-t*q^i),i=1..N),q=0, N+1 ):
[seq(coeff(f,q,i),i=0..N)]:
end:


#From C11.txt
#pnk(n,k):
pnk:=proc(n,k) local k1:
option remember:
if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN(0):
fi:
if n=k then
 RETURN(1):
else
 RETURN(add(pnk(n-k,k1),k1=1..min(k,n-k))):
fi:
end:

#pnSeqGF(N):
pnSeqGF:=proc(N) local i,q,f:
f:=taylor(mul(1/(1-q^i),i=1..N),q=0, N+1 ):
[seq(coeff(f,q,i),i=0..N)]:
end:

#From C10.txt
#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:

#From C9.txt
#RRV(n,K): A random random variable on {1, ... n}
#where the values are from 0 to K (uniformly)
RRV:=proc(n,K) local i,ra:
ra:=rand(0..K):
[ seq(ra(),i=1..n)]:
end:


#Ave(X): inputs a RV X (given as list of values where the value of
#Mr. or Ms. i is L[i]) outputs the average (alias expectation, alias mean)
Ave:=proc(X) local i:add(X[i],i=1..nops(X))/nops(X):end:

#Mom(X,k): inputs a RV X (given as list of values where the value of
#Mr. or Ms. i is L[i]) outputs the k-th moment about the mean
Mom:=proc(X,k) local i,mu:
mu:=Ave(X):add((X[i]-mu)^k ,i=1..nops(X))/nops(X):
end:

#SA(X,K): inputs a discrete r.v. X (given as a list of numbers) and
#a pos. integer K outputs a list of length K whose
#(i) first entry is the average
#(ii) second entry variance
#The remaining entries are the scaled-third moment (about the mean) aka skeweness
#scaled fourth-moment (about the mean) aka kurtosis and so on
#up to the K-th moment
#Changed type(X,list) to type(X,Array) or type(X,list)
SA:=proc(X,K) local i,M,mu,sig2:
if not ((type(X,Array) or type(X,list)) and type(K,integer) and K>=2) then
 print(`Bad input`):
 RETURN(FAIL):
fi:
mu:=Ave(X,1):
sig2:=Mom(X,2):
M:=[mu,sig2]:

if sig2=0 then
 RETURN(M):
fi:

[op(M), seq(evalf(Mom(X,i)/sig2^(i/2)),i=3..K)]:

end:

#From hw9
#Corrected: f2 should be x*diff(x*diff(f1/(x^mu),x),x).
SAc := proc(f,x,K) local f1,f2,mu,sig2,m,k:
if subs(x=1,f)=0 then
 RETURN(0):
fi:

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

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

m := [mu,sig2]:
for k from 3 to K do
f2:=x*diff(f2,x):
m := [op(m), subs(x=1,f2)/sig2^(k/2)]:
od:
m:
end: