#OK to post homework #Victoria Chayes, 2/21/21, Assignment 9 Help:=proc(): print(`SAc(f,x,K)`): end: with(combinat): #1: Convince yourself that if f(x) is the pgf of a random variable X, then μ=f'(1), Var(X)= (x d/dx)^2 (f(x)/xμ)|x=1, and more generally the kth coment about the mean is (x d/dx)^k (f(x)/xμ)|x=1. Using this, finish writing procedure SAc(f,x,K), that we barely started in class. Make sure that SAc(Pgf(X,x),x,8)=SA(X,8) for several X obtained by using RRV(100,1000) (say) SAc:=proc(f,x,K) local f1,L,mu,i: if subs(x=1,f)=0 then RETURN(0): fi: f1:=f/subs(x=1,f): mu:=subs(x=1,x*diff(f1,x)): f1:=f1/x^mu: L:=[]: for i from 1 to K do f1:=x*diff(f1,x): L:=[op(L), subs(x=1,f1)]: od: [mu,op(2..K,L)]: end: #Testing did not produce the same solution. I went back through class code from 2019 and 2020, and found that testing did produce the same solution to the code we wrote as a class then. We also had code for moments, not just moments about the mean, and that was producing different solutions than the Mom code, so that might be where the difference is coming from? But I could not for the life of me locate the mistake, in either my code/ our old code, or C9.txt. #2 We will soon see (and you probably already know that) that the pdf of tossing a fair coin n times is ((1+x)/2)^n. Using the above mentioned SAc that you wrote above, and leaving n symbolic, find closed-form expressions for the average, variance, and third through the 10th scaled moments. By using Maple's limit command, find the limits. Can you conjecture an explicit expression for the k-th scaled moment as an expression in k? # S:=SAc(((1+x)/2)^n,x,10) # [1 1 3 2 1 15 3 15 2 1 # S := [- n, - n, 0, -- n - - n, 0, -- n - -- n + - n, 0, # [2 4 16 8 64 32 4 # # 105 3 147 2 17 105 4 # - --- n + --- n - -- n + --- n , 0, # 64 64 16 256 # # 31 945 5 4095 3 1575 4 1185 2] # -- n + ---- n + ---- n - ---- n - ---- n ] # 4 1024 256 256 64 ] #However, this leads to #[seq(limit(S[i],n=infinity),i=1..10)] # [infinity, infinity, 0, infinity, 0, infinity, 0, infinity, 0, infinity] #Clearly, odd moments past 1 are 0, and even moments are polynomials of leading degree k/2, but I do not see a solution beyond that / my problem may be that my program SAc is indeed incorrect. #3 Using procedure RandSPn(n), with n=100,200,300 from C8.txt , 10000 times, and SA(X), estimate the average and variance of the random variable "number of parts" of a random set-partion of {1,...,n}. See how close it is to the exact values μ:=Sum(k*Snk(n,k),k=0..n)/Bell(n), sig2:=Sum((k-μ)2*NuSPnk(n,k),k=0..n)/Bell(n). #We first calculate # evalf(sum('k*NuSnk(100,k)',k=0..100)/bell(100)) # -84 # 1.331974547 10 # evalf(sum('k*NuSnk(200,k)',k=0..200)/bell(200)) # -214 # 2.572135990 10 # evalf(sum('k*NuSnk(300,k)',k=0..300)/bell(300)) # -362 # 3.184953145 10 # sum('(k-1.331974547*10^(-84))^2*NuSPnk(100,k)',k=0..100)/bell(100) # 825.1525147 # sum('(k-2.572135990*10^(-214))^2*NuSPnk(200,k)',k=0..200)/bell(200) # 2506.844431 # sum('(k-3.184953145*10^(-362))^2*NuSPnk(300,k)',k=0..300)/bell(300) # 4852.870540 # These looked a little bit wonky to me. If the original mu was supposed to be with NuSPnk too, then the values are # evalf(sum('k*NuSPnk(100,k)',k=0..100)/bell(100)) # 28.62528186 # evalf(sum('k*NuSPnk(200,k)',k=0..200)/bell(200)) # 49.97510229 # evalf(sum('k*NuSPnk(300,k)',k=0..300)/bell(300)) # 69.57332709 # sum('(k-28.62528186)^2*NuSPnk(100,k)',k=0..100)/bell(100) # 5.745753494 # sum('(k-49.97510229)^2*NuSPnk(200,k)',k=0..200)/bell(200) # 9.333583258 # sum('(k-69.57332709)^2*NuSPnk(300,k)',k=0..300)/bell(300) # 12.42269762 #Now to try to find these numbers with random variables: # X:=[seq(nops(RandSPn(100)),i=1..1000)]: # evalf(SA(X,2)) # [28.56100000, 5.806279000] # X:=[seq(nops(RandSPn(200)),i=1..1000)]: # evalf(SA(X,2)) # [50.11600000, 9.126544000] # X:=[seq(nops(RandSPn(300)),i=1..1000)]: # evalf(SA(X,2)) # [69.52300000, 12.26747100] #This makes me think that indeed there was a typo still in the question, as these are comparable to the second set of values, and are, in fact, very good matches to them. ########################### ### Programs From Class ### ########################### RRV:=proc(n,K) local i,ra: ra:=rand(0..K): [ seq(ra(),i=1..n)]: end: Ave:=proc(X) local i:add(X[i],i=1..nops(X))/nops(X):end: Mom:=proc(X,k) local i,mu: mu:=Ave(X):add((X[i]-mu)^k ,i=1..nops(X))/nops(X): end: SA:=proc(X,K) local i,M,mu,sig2: if not (type(X,list) and type(K,integer) and K>=2) then print(`Bad input`): RETURN(FAIL): fi: mu:=Ave(X,1): sig2:=Mom(X,2): M:=[mu,sig2]: if sig2=0 then RETURN(M): fi: [op(M), seq(evalf(Mom(X,i)/sig2^(i/2)),i=3..K)]: end: Pgf:=proc(X,x) local i:add(x^X[i],i=1..nops(X))/nops(X):end: Snk:=proc(n,k) local S1,S2,s: option remember: if k<0 or k>n then RETURN({}): fi: if n=0 then if k=0 then RETURN({{}}): else RETURN({}): fi: fi: S1:=Snk(n-1,k): S2:=Snk(n-1,k-1): S1 union {seq( s union {n}, s in S2)}: end: NuSnk:=proc(n,k) option remember: if k<0 or k>n then RETURN(0): fi: if n=0 then if k=0 then RETURN(1): else RETURN(0): fi: fi: NuSnk(n-1,k)+ NuSnk(n-1,k-1): end: RandSPn:=proc(n) local i,k,die1: die1:=[seq(NuSPnk(n,i),i=1..n)]: k:=LD(die1): RandSPnk(n,k): end: RandSPnk:=proc(n,k) local c,SP1,k1: if n=1 then if k=1 then RETURN({{1}}): else RETURN(FAIL): fi: fi: c:=LD([NuSPnk(n-1,k-1), k*NuSPnk(n-1,k)]): if c=1 then SP1:=RandSPnk(n-1,k-1): RETURN(SP1 union {{n}}): else SP1:=RandSPnk(n-1,k): k1:=rand(1..k)(): RETURN( { op(1..k1-1,SP1), SP1[k1] union {n}, op(k1+1..k,SP1)}): fi: end: LD:=proc(L) local n,i,su,r: n:=nops(L): r:=rand(1..convert(L,`+`))(): su:=0: for i from 1 to n do su:=su+L[i]: if r<=su then RETURN(i): fi: od: end: NuSPnF:=proc(n) local k: option remember: if n=0 then RETURN(1): fi: add( binomial(n-1,k)*NuSPnF(n-1-k) , k=0..n-1): end: NuSPnk:=proc(n,k) option remember: if n=1 then if k=1 then RETURN(1): else RETURN(0): fi: fi: NuSPnk(n-1,k-1)+ k*NuSPnk(n-1,k): end: