#OK to post homework
#George Spahn, 2/28/2021, Assignment 11

#RandPnk(n,k): random integer partition of n with largest part k
RandPnk:=proc(n,k) local Die,d,P:

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:

Die:=[seq(pnk(n-k,k1),k1=1..min(n-k,k))]:
d:=LD(Die):
P:=RandPnk(n-k,d):
[k, op(P)]:
end:


#RandPn(n): random integer partition of n
RandPn:=proc(n) local Die,d:
Die:=[seq(pnk(n,k),k=1..n)]:
d:=LD(Die):
RandPnk(n,d):
end:



# EstStatAnalOfLargestPart(n,K,N) Inputs a positive integer n
# and (using previous functions) with N trials, estimates the expectation, variance,
# and the 3rd- through K-th scaled moments, of the r.v. "number of parts" when 
# running over the sample space of (integer) partitions of n.
EstStatAnalOfLargestPart:=proc(n,K,N) local L:
L:=[seq(nops(RandPn(n)),i=1..N)]:
SA(L,K):
end:

# EstStatAnalOfLargestPart(100, 6, 10000) 
# gives ~ [21, 67, 1.1, 5.2, 17, 75]

# EstStatAnalOfLargestPart(200, 6, 10000) 
# gives ~ [34, 137, 1.1, 5.2, 16, 69]

#The third and fourth scaled moments are indeed similar



# pnSeqGFt(N,t) Inputs N, and a variable t, and outputs the first N+1 
# (starting at n=0) of the sequence Pn(t) qn

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

#The above doesn't work for getting the taylor approximation,
# I didn't figure out how to fix it.


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



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

#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
SA:=proc(X,K) local i,M,mu,sig2:
if not (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: