#OK to post homework
#Yukun Yao, Feb. 1, Assignment 3



##Problem 3## The scaled k-th moments of a probability distribution is mk/ (m2)^(k/2), where mk is the k-th moment about the mean. Using MomsAM(f,t,K): Write a procedure ScaledMomsAM(f,t,K): that inputs a probability generating function f (of the variable t), and outputs the first K scaled moments about the mean. The first and second entries of ScaledMomsAM(f,t,K) should be the mean and variance, but starting at the 3rd it should be the scaled moments.

ScaledMomsAM:=proc(f,t,K) local L,ave,var,k:
L:=MomsAM(f,t,K):
ave:=L[1]:
var:=L[2]:
[ave, var, seq(L[k]/var^(k/2), k=3..K)]:
end:



##Problem 4## I let F:=n->ScaledMomsAM(((1+t)/2)^n, t, 14) and tried large n’s. I found that the 3rd through 14th scaled moments converge to a finite number when n goes to infinity, respectively. 



##Problem 5## By letting H:= n-> evalf(ScaledMomsAM(PQs(n, t), t, 4)), I got 
#H(50) = [258.9189445, 783.2000949, 1.007591274, 4.460496096]
#H(60) = [330.9441904, 1171.868087, .9940500268, 4.439700480]
#H(70) = [406.2628196, 1641.296046, .9830305181, 4.421410263]
#H(80) = [484.4076432, 2191.950937, .9738804886, 4.405517724]
#H(90) = [565.0278497, 2824.183446, .9661484404, 4.391689671]
#H(100) = [647.8502586, 3538.266678, .9595165948, 4.379588738]

#The estimated skewness is around 0.9 and the estimated kurtosis is around 4. 











###### Below is from C3.txt ######

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

#MFtoGF(M, t): inputs a prob. mass function of a finite and discrete prob. dist. given
#[  ... [outcome, Pr(outcome)], ...] and a variable name (of type symbol) and outputs
#the prob. gen. function
MFtoGF:=proc(M,t) local i:

add( t^M[i][1]*M[i][2], i=1..nops(M)):
end:

 #GFtoMF(f,t): inputs a  prob. gen. function in the variable t, and outputs the prob. mass function
 #for a finite discrete prob. dist.
 #as a list of pairs [outcome, Pr(outcome)]

GFtoMF:=proc(f,t) local i,M,c:
 if subs(t=1,f)<>1 then
    print(`Not prob. dist.`):
    RETURN(FAIL):
 fi:

  M:=[]:

  for i from ldegree(f,t) to degree(f,t) do
      c:=coeff(f,t,i):
       if c<0 then
          print(`OMG NOT a prob. `):
         RETURN(FAIL):
       elif c>0 then
         M:=[op(M),[i,c]]:
       fi:
od:
M:
end:
      
 #Moms(f,t,K): inputs a prob. gen. function and outputs the first K moments
 Moms:=proc(f,t,K) local L,i,F:
 if subs(t=1,f)<>1 then
    print(`Not prob. dist.`):
    RETURN(FAIL):
 fi:
 
 F:=f:

  L:=[]:

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

 end:

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