#March 11, 2013, Larry Shepp Urn Help:=proc(): print(` Larry(m,p), LarryE(m,p) `): print(`LarryTale(N),beta(p) , LarryPGF(m,p,x)`): print(`PrDist(m,p) `): print(` Stat(P,x,K) `): end: #LarryE(m,p): input non.neg. integers m and p #representing a box (urn) with m (-1)dollar bills #and p 1-one dollar bills, outputs the #expected gain of the game (where at any given time) #you pick at random one or the other dollars bills #Pr. of picking a dollar bill is p/(m+p) and #Pr. of picking a -1 dollar bill m/(m+p) #And you quit whenever you want, LarryE:=proc(m,p) option remember: max(Larry(m,p),0): end: #Larry(m,p): input non.neg. integers m and p #representing a box (urn) with m (-1)dollar bills #and p 1-one dollar bills, outputs the #expected gain of the game (where at any given time) #you pick at random one or the other dollars bills #Pr. of picking a dollar bill is p/(m+p) and #Pr. of picking a -1 dollar bill m/(m+p) #And you quit whenever you want Larry:=proc(m,p) option remember: if m=0 then RETURN(p): fi: if p=0 then RETURN(-1): fi: p/(m+p)*(1+LarryE(m,p-1))+ m/(m+p)*(-1+LarryE(m-1,p)): end: #LarryTable(N): LarryTable:=proc(N) local m,p: for p from N by -1 to 0 do lprint(seq(evalf(LarryE(m,p),3),m=0..N)): od: end: #beta(p): The largest m for which LarryE(m,p) is >0 #beat(p)-p is asymtotically sqrt(p) times a (explicit) Constant beta:=proc(p) local m: for m from p while LarryE(m,p)>0 do od: m-1: end: #LarrtPGF(m,p,x): The generalized "polynomial" #whose exponets are numbers (rational) and #where the coeff. of x^a is the prob. of exactly getting #a dollars when you end the game (but following the #stupid max. exp. strategy) LarryPGF:=proc(m,p,x) option remember: if LarryE(m,p)=0 then RETURN(1): fi: if m=0 then RETURN(x^p): fi: expand(p/(m+p)*x*LarryPGF(m,p-1,x)+ m/(m+p)/x*LarryPGF(m-1,p,x)): end: #PrDist(m,p): the list of pairs [i,p(i)] where #the prob. of getting i is p(i) PrDist:=proc(m,p) local olivia,x,i: olivia:=LarryPGF(m,p,x): [seq([i, coeff(olivia,x,i)], i=ldegree(olivia,x)..degree(olivia,x))]: end: #Stat(P,x,K): inputs a prob. g.f. P in x and a pos. integer K #outputs a list of length K where the first entry if the #average, the second the variance, and the r-th is the s-called #alpha coefficiient m[r]/m[2]^(r/2) where m[r] is the #r-th moment about the mean #try: #evalf(Stat(((1+x)/2)^1000,x,10)); #Taken from homework problem, not done in class Stat:=proc(P,x,K) local av,va,L,P1,r: P1:=P/subs(x=1,P): av:=subs(x=1,diff(P1,x)): P1:=expand(P1/x^av): P1:=expand(x*diff(P1,x)): P1:=expand(x*diff(P1,x)): va:=subs(x=1,P1): L:=[av,va]: for r from 3 to K do P1:=expand(x*diff(P1,x)): L:=[op(L),subs(x=1,P1)/va^(r/2)]: od: L: end: