#OK to post homework
#Blair Seidler, 2021-02-21 Assignment 9

with(combinat):

Help:=proc(): print(`SAc(f,x,K), FindCoinMoments()`): end:

#f:=Sum(x^a[i], i=1..n);
#diff(f,x)= subs(x=1,Sum(a[i]*x^(i-1),i=1..n))/n;
#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:

#### Yes, SAc(Pgf(X,x),x,8)=SA(X,8) for a variety of X:=RRV(something,something)
#### I removed the evalf before trying to use SAc for problem 2

# 2.

#FindCoinMoments()
FindCoinMoments:=proc() local f,M,K:
K:=10:
f:=((1+x)/2)^n:

M:=simplify(SAc(f,x,K)):

print(M):

[seq(limit(M[i],n=infinity),i=3..K)]:

end:

#### I ran this procedure with K=20 (so looking at the limits of the first
#### 20 moments), and found the limits in the OEIS:
#A123023: a(n) = (n-1)*a(n-2), a(0)=1, a(1)=0.
#One of the descriptions is the n-th moment of the standard normal distribution,
#so I think that I'm in the right place.

# 3.

#ChkParts()
ChkParts:=proc() local N,n,Xsize,X,i,musig,mu,sig2,k:
N:=[100,200,300]:
Xsize:=10000:

for n in N do
  X:=[seq(nops(RandSPn(n)),i=1..Xsize)]:
  musig:=SA(X,2):
  printf(`Experimental: n=%d, mu=%f, sig2=%f\n`,n,musig[1],musig[2]):
  mu:=add(seq(k*NuSPnk(n,k),k=0..n))/NuSPn(n):
  sig2:=add(seq((k-mu)^2*NuSPnk(n,k),k=0..n))/NuSPn(n):
  printf(`Exact:        n=%d, mu=%f, sig2=%f\n`,n,mu,sig2):
od:
RETURN():
end:

#### Results: (pretty close to the exact values:
#Experimental: n=100, mu=28.596800, sig2=5.662630
#Exact:        n=100, mu=28.625282, sig2=5.745753

#Experimental: n=200, mu=49.969200, sig2=9.164651
#Exact:        n=200, mu=49.975102, sig2=9.333583

#Experimental: n=300, mu=69.667200, sig2=12.430044
#Exact:        n=300, mu=69.573327, sig2=12.422698


#### From C8.txt ####
#RandSPn(n): (from Yuxuan Yang's hw7): A randon set-partition of {1, ..,n}
RandSPn:=proc(n) local i,k,die1:
die1:=[seq(NuSPnk(n,i),i=1..n)]:
k:=LD(die1):
RandSPnk(n,k):
end:

#NuSPnk(n,k): The number of set-partitions of {1,...,n}
#into exactly k parts (aka Strirling numbers of the second kind)
NuSPnk:=proc(n,k) option remember:

if n=1 then
  if k=1 then
   RETURN(1):
  else
   RETURN(0):
 fi:
fi:
NuSPnk(n-1,k-1)+ k*NuSPnk(n-1,k):
end:

NuSPn:=proc(n) local k:
option remember:
add(NuSPnk(n,k),k=1..n):
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:

#RandSPnk(n,k): a random set partition of {1,..,n} with k parts
RandSPnk:=proc(n,k) local c,SP1,k1:
if n=1 then
 if k=1 then
   RETURN({{1}}):
  else
   RETURN(FAIL):
 fi:
fi:
c:=LD([NuSPnk(n-1,k-1), k*NuSPnk(n-1,k)]):

if c=1 then
 SP1:=RandSPnk(n-1,k-1):
 RETURN(SP1 union {{n}}):
 else
 SP1:=RandSPnk(n-1,k):
  k1:=rand(1..k)():
  RETURN( { op(1..k1-1,SP1), SP1[k1] union {n}, op(k1+1..k,SP1)}):
  
fi:
end:

#### From C9.txt ####

Help9:=proc(): print(`RRV(n,K), Ave(X), Mom(X,k) , SA(X,K), Pgf(X,x) `): end:

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

#Pgf(X,x): inputs a RV X (given as list of values where the value of
#Mr. or Ms. i is L[i]) outputs the prob. generating function of X
#using the variable
Pgf:=proc(X,x) local i:add(x^X[i],i=1..nops(X))/nops(X):end: