#It's OK to post.  Tong Cheng, Feb 3rd, 2019, Assignment 3

Help:= proc(): ( ` ScaledMomsAM ` ) : end:

##From C2.txt
#PQs(n,t):the prob. generating function of the running time of Quicksort
 PQs:=proc(n,t) local k:
   option remember:
  if n=0 or n=1 then
   RETURN(1):
  else:

  RETURN(expand(t^(n-1)*add(1/n*PQs(k-1,t)*PQs(n-k,t),k=1..n))):
 fi:

 end:

###end from C2.txt

#MomsAM(f,t,K): inputs a prob. gen. function and outputs the first K moments about the mean
 MomsAM:=proc(f,t,K) local L,i,F,ave:
 if subs(t=1,f)<>1 then
    print(`Not prob. dist.`):
    RETURN(FAIL):
 fi:
 
 F:=f:

  L:=[]:

   F:=t*diff(F,t):
    ave:=subs(t=1,F):

    F:=f/t^ave:

   for i from 1 to K do
    F:=t*diff(F,t):
    L:=[op(L),subs(t=1,F)]:
   od:
   [ave,op(2..nops(L),L)]:

 end:


# inputs a probability generating function f(t)
# outputs the first K scaled moments about the mean
ScaledMomsAM:= proc(f, t, K)  local a:

 a:=MomsAM(f,t,K)/( MomsAM(f,t,2)[2]^(K/2) ):
 a:

end:

# problem5 procedure

Skew_Kurtosis:=proc(t, K) local  i, g, seq:
global Sk, Ku:
Sk:=[]:
Ku:=[]:
g:='g':
for i from 50 by 10 to 100 do
g:=PQs(i,t):
seq:=ScaledMomsAM(g, t, K):
Sk:=[ op(Sk),  seq[3] ]:
Ku:=[ op(Ku), seq[4] ]:



od:

evalf(Sk);
Ku;
end:
#Then do "evalf(Sk), evalf(Ku)"