#Nathan Fox
#Homework 13
#I give permission for this file to be posted online

##Read old files
read(`C13.txt`):

#Help procedure
Help:=proc():
print(` AveExp(L) , AveExpSeq(f,n,N) , Stat(P,x,R) `):
end:

##Problem 1
#AveExp(L): inputs a list L, of different numbers and outputs
#the average of the expected value of DressUp(T,L) over all
#trees T with n=nops(L) leaves.
AveExp:=proc(L) local T,ret:

ret:=0:

for T in BT(nops(L)) do
    ret:=ret + ExpV(DressUp(T,L)):
od:

ret / C(nops(L)):

end:

##Problem 2
#AveExpSeq(f,n,N): inputs an expression f in n, and outputs
#the list of length N, whose m-th entry is AveExp([f(1), ..., f(m)])
AveExpSeq:=proc(f,n,N) local i,m:

[seq(AveExp([seq(subs(n=i,f),i=1..m)]),m=1..N)]:

end:

##Problem 3
#Conjecture: AveExp([seq(i,i=1..n)])=(n+1)/2

##Problem 4
#Stat(P,x,R): inputs an expression P in x, and a pos. integer
#R larger than 2, and checks that indeed P(1)=1, and returns
#a list of length R, whose first entry is the average, its
#second entry is the variance, and for r=3..R, the r-th entry
#is N[r].
Stat:=proc(P,x,R) local ret,Pp,F,i:

if subs(x=1,P)<>1 or R<=2 then
    return FAIL:
fi:

Pp:=diff(P,x):
ret:=[subs(x=1,Pp)]:

F:=P/x^(ret[1]):
Pp:=diff(diff(F,x),x):
ret:=[op(ret),subs(x=1,Pp)]:

for i from 3 to R do
    Pp:=diff(Pp,x):
    ret:=[op(ret),subs(x=1,Pp)/(ret[2]^(i/2))]:
od:

ret:

end:

##Problem 5
#evalf(Stat(((1+x)/2)^10000,x,10))=[5000., 2500., -0.06000000000, 3.004200000, -0.6003600000, 15.09903720, -6.316384535, 106.8928692, -76.11462140, 981.0475136]
#evalf(Stat(((1+x)/2)^100000,x,10))=[50000., 25000., -0.01897366596, 3.000420000, -0.1897480438, 15.00990037, -1.992752921, 105.1890287, -23.92307745, 948.5923741]
#evalf(Stat(((1+x)/2)^1000000,x,10))=[500000., 250000., -0.006000000000, 3.000042000, -0.06000036000, 15.00099000, -0.6300163800, 105.0189003, -7.560514085, 945.3591137]
#It appears that the odd moments (after the first) are going to
#zero and that the even moments (after the second) are stabilizing
#to the product of the first n/2 odd numbers (e.g. 4th moment going
#to 3).  The first and second moments are increasing tenfold each
#time, as they should be.

##Problem 6
#evalf(seq(Stat(PGG(StPete(n),x),x,10),n=2..7)) =
#[3., 1., -3., 12., -60., 360., -2520., 20160., -181440., 0.1814400e7 ],
#[4., 6., -0.4082482906, 1.833333333, -4.082482906, 11.66666667, -38.10317379, 140., -571.5476066, 2566.666667],
#[5., 21., 0.9663935137, 2.666666667, 2.790813834, 5.597667638, 2.731434817, 8.992549401, -8.499682253, 30.76495904],
#[6., 58., 2.159764430, 6.893579073, 19.19244145, 53.30938128, 139.2612693, 348.2052678, 821.7802262, 1836.243155],
#[7., 141., 3.481491879, 15.43503848, 67.32488035, 292.7185691, 1253.232774, 5271.068351, 21741.86819, 87874.63634],
#(the rest of the high moments were too large to output nicely)

#The average is always n+1
#The variance appears to be 6*2^n-n^2-4n-6