# Joey Reichert # 640:640 Spring 2013 # March 11, 2013 # hw14.txt # I give permission to post Help:=proc() : print("PrWinBreakEvenLose(m,p);"): end: # Problem #1 # According to Larry Shepp, (beta(p)-p)/sqrt(2*p) # tends to a certain fixed constant, that we may # call the Shepp constant. Estimate it. # Using p = 1200, # (beta(p)-p)/sqrt(2*p) is approx. 0.8369089956 # Problem #2 # Write a program PrWinBreakEvenLose(m,p) that inputs # two positive integers m and p and outputs a list of # length 3 whose first entry is the probability of # exiting with a positive amount, zero amount, and # negative amouts, respectively. By examinining # PrWinBreakEvenLose(beta(p),p) for many p, decide # whether it is a good idea to play the Shepp "maximize # expectation strategy", or whether it is too risky, if # you hate to lose. PrWinBreakEvenLose:=proc(m,p) local i,prW,prL,prBE: option remember: if m = 0 and p = 0 then return [0,0,0]: elif m = 0 then return [1,0,0]: elif p = 0 then return [0,0,1]: fi: prW:=p/(2*(m+p))( PrWinBreakEvenLose(m,p-1)[1] + PrWinBreakEvenLose(m-1,p)[1] ); prL:=m/(2*(m+p))( PrWinBreakEvenLose(m,p-1)[3] + PrWinBreakEvenLose(m-1,p)[3] ); prBE:=1-prW-prL: return [prW,prBE,prL]: end: # evalf(seq(PrWinBreakEvenLose(beta(p),p),p=1..1000)); # returns results ranging from [0.2,0.5,0.3] to # [0.2455, 0.5, 0.2545], which suggests in the best case, using # the Shepp strategy will break even 50% of the time, but # it is more likely you will lose money than win, in the long run. # Hence it is not too risky to play if you hate to lose, # as long as you don't mind breaking even. But if you want to WIN, # it is too risky. # Problem #3 # Verify empirically the inequalities proved in # William Boyce's beautiful article for m,p up to 200. # {seq(seq(evalb(Larry(m,p+1)-Larry(m,p) >= 0),m=0..200),p=0..200)}; # {seq(seq(evalb(Larry(m+1,p)-Larry(m,p) <= 0),m=0..200),p=0..200)}; # {seq(seq(evalb(Larry(m+1,p+1)-Larry(m,p) >= 0),m=0..200),p=0..200)}; # All output {true}. # Problem #4 # In William Boyce's beautiful article, it is proved # (Lemma 3.6) that V(1,p)=p^2/(p+1). Conjecture # expressions for V(2,p) and V(3,p), and if possible, # prove them by induction. # To verify: # {seq(evalb(Larry(1,p) = p^2/(p+1)),p=1..200)}; # Returns {true}. # V(2,p) = (p-1)*p/(p+1) # {seq(evalb((p-1)*p/(p+1) = Larry(2,p)),p=1..200)}; # Return {true}. # I am sure that V(3,p) can be easily found, but # I am omitting it here for time. ### C14.txt stuff ### #March 11, 2013, Larry Shepp Urn Help14:=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 #beta(p)-p is asymtotically sqrt(p) times a (explicit) Constant beta:=proc(p) local m: option remember: for m from p while LarryE(m,p)>0 do od: m-1: end: #LarryPGF(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: