#OK to post homework
#Wanying Rao, 02/21/2021, Assignment 9

with(combinat):

#1
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(f/(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), evalf(subs(x=1,f2)/sig2^(k/2))]:
od:
m:
end:


#2
#[n/2, n/4, 0., 16.*(0.1875000000*n^2 - 0.1250000000*n)/n^2, 0., 64.*(0.2500000000*n - #0.4687500000*n^2 + 0.2343750000*n^3)/n^3, 0., 256.*(-1.062500000*n + 2.296875000*n^2 - #1.640625000*n^3 + 0.4101562500*n^4)/n^4, 0., 1024.*(7.750000000*n - 6.152343750*n^4 + #0.9228515625*n^5 - 18.51562500*n^2 + 15.99609375*n^3)/n^5]
#lim:[inf,inf,0,3,0,15,0,105,0,945]
#for all k>=3
#if k is even, then k-th scaled moment is (k-1)!! = (k-1)(k-3)...1.
#if k is odd, then k-th scaled moment is 0.


#3
#n=100
#X := [seq(nops(RandSPn(100)), i = 1 .. 1000)]:
#SA(X,2);
#I got ave= 28.7, var=5.8
#mu := evalf(add(k*NuSnk(100, k), k = 0 .. 100)/bell(100));
#                                        -84
#                    mu := 1.331974547 10   

#sig2 := add((k - mu)^2*NuSPnk(100, k), k = 0 .. 200)/bell(100);
#                      sig2 := 825.1525145


#n=200
#X := [seq(nops(RandSPn(200)), i = 1 .. 1000)];
#SA(X, 2);
#ave=50,var=9.1
#mu := evalf(add(k*NuSnk(200, k), k = 0 .. 200)/bell(200));
#                                        -214
#                   mu := 2.572135990 10    

#sig2 := add((k - mu)^2*NuSPnk(200, k), k = 0 .. 200)/bell(200);
#                      sig2 := 2506.844431

#n=300
#X := [seq(nops(RandSPn(300)), i = 1 .. 1000)];
#SA(X, 2);
#ave=69.7, var=11.7
#mu := evalf(add(k*NuSnk(300, k), k = 0 .. 300)/bell(300));
#                                        -362
#                    mu := 3.184953145 10    

#sig2 := add((k - mu)^2*NuSPnk(300, k), k = 0 .. 300)/bell(300);
#                      sig2 := 4852.870541


#From hw8WanyingRao.txt
#By using LD(L) and RandSPnk(n,k) write a procedure
#RandSPn(n)
#that inputs a positive integer n and outputs a (uniformly) at random set partition of #{1,...,n}
RandSPn := proc(n) local i,L,r:
L := [seq(NuSPnk(n,i),i=1..n)]:

r := LD(L):
RandSPnk(n,r):
end: