# Matthew Russell # Experimental Math # Homework 13 # I give permission to post this ##################################################################################### # Write a program # AveExp(L) # that 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 leaves. AveExp:=proc(L) local B: B:=BT(nops(L)): return add(ExpV(DressUp(T,L)), T in B)/nops(B): end: ##################################################################################### # Write a program # AveExpSeq(f,n,N) # that 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: return [seq(AveExp([seq(subs(n=i,f),i=1..m)]), m=1..N)]: end: ##################################################################################### # Can you conjecture explicit expressions in n for the following? (I don't know the answer) # AveExp([seq(i,i=1..n)]): # = 1/2*(n+1)/2 # AveExp([seq(i^2,i=1..n)]): # = # AveExp([seq(2^i,i=1..n)]): # = # AveExp([seq(3^i,i=1..n)]): # = ##################################################################################### # The n-th moment of a random variable X # is defined to be the expectation of X^n. # By elementary (discrete!) probability theory # if P(x) is the prob. generating function # (so we must have P(1)=1), # the average is P'(1) (let's call it a)). # The centralized p.g.f is F(x):=P(x)/x^a. # The r-th moment about the mean is m[r]:=(xd/dx)^r P(x) # (where d/dx is the differentiation operator with respect to x). # In particular m[2] is called the variance # and its square-root is called the standard deviation. # Finally the normalized moments are # N[r]:=m[r]/m[2]^(r/2) # Write a program # Stat(P,x,R) # that 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 me, P1, va, sd: if R<=2 then return FAIL: fi: if subs(x=1,P)<>1 then return FAIL: fi: me:=subs(x=1,diff(P,x)): P1:=P/x^me: va:=subs(x=1,m(P1,x,2)): sd:=va^(1/2): return [me, va, seq(subs(x=1,m(P1,x,r))/sd^r,r=3..R)]: end: m:=proc(P,x,r) option remember: if r=0 then return P: fi: return m1(m(P,x,r-1),x): end: m1:=proc(P,x) return x*diff(P,x): end: # evalf(Stat(((1+x)/2)^10000,x,10)); # = [5000., 2500., 0., 2.999800000, 0., 14.99700016, 0., 104.9580059, 0., 944.3701638] # evalf(Stat(((1+x)/2)^100000,x,10)); # = [50000., 25000., 0., 2.999980000, 0., 14.99970000, 0., 104.9958001, 0., 944.9370016] # evalf(Stat(((1+x)/2)^1000000,x,10)); # = [5.00000*10^5, 2.50000*10^5, 0., 2.999998000, 0., 14.99997000, 0., 104.9995800, 0., 944.9937000] # Can you conjecture a trend? # It appears that evalf(Stat(((1+x)/2)^(10^k),x,10)); gives # [1/2*10^k, 1/4*10^k, 3, 0, 3*5, 0, 3*5*7, 0, 3*5*7*11] ##################################################################################### # What is Stat(FGF(StPete(n),x),x,10) for n from 1 to 20? # Can you guess a formula for # (i) the average (easy) # (ii) variance (m[2]) # (iii) N[r]? (I don't know the answer) in terms of n?