# Thomas Sznigir # Experimental Math Spring 2013 # Homework 24 # I hereby give permission to post this # Problem 1: Write a procedure OptimalGamblingTerse(N,i,p,T) # that does not tell you the story but returns TWO outputs: # true (or fals) if you came out alive (or dead) # A list of the intermediate probabilities (of length less # than T in the lucky case, and of length T in the sad case) OptimalGamblingTerse:=proc(N,i,p,T) local time,strategy,coin,money,L,result: time:=T: money:=i: L:=[]: while money>0 and money0 do strategy:=BestSDD(N,money,p,time): result:=strategy[1]: if result = Pray then L:=[op(L),0$time]: break: fi: L:=[op(L),1.0*strategy[2]]: coin:=LC(p): if coin then money:=money+result: else money:=money-result: fi: time:=time-1: if money=0 then L:=[op(L),0$time]: fi: od: evalb(money=N),L: end: # Problem 2: Write a procedure ManyOptimalGambling(N,i,p,T,K) # that runs OptimalGamblingTerse(N,i,p,T) K times and returns # TWO NUMBERS (in floating point): # the ratio of times the gambler came out alive # the "closest call", the lowest ever (intermediate) # probability of survival that ever showed up in the games # that ended happily. ManyOptimalGambling:=proc(N,i,p,T,K) local k,result,rec,challenge,success: rec:=1: success:=0: for k from 1 to K do: result:=OptimalGamblingTerse(N,i,p,T): if result[1] then challenge:=min(op(result[2])): if challenge < rec then rec:=challenge: fi: success:=success+1: fi: od: if success > 0 then 1.0*success/K,rec: else 0,0: fi: end: # Problem 3: Try out ManyOptimalGambling(20,9,11/20,10,5000); # How close is the first output to BestSDD(20,9,11/20,10)[2];? # What is the second output? # We got 0.5700000000 for the first output, compared to # 0.5651360230 for BestSDD. # The second output was 0.7297623438e-1 # Problem 4: Write a procedure # BestKellyFactor(N,i,p,TimeIsMoney,resolution) where # N is the ultimate desired capital # i is your current capital # p is the probability of a successful one round (usually larger than 1/2) # TimeIsMoney is the cost that you have to pay by playing a single round # resolution is a small positive number (say 1/200) # that tries out all the possible Kelly strategies, starting with # f=resolution, and then f=2*resolution,..., and so on, all the way to the # bold scenario f=1, and finds the f that maximizes # (VGD(p,N,TiS(N))[i]-VGD(p,N,KeS(N,f)[i]))*TimeIsMoney # -(VGW(p,N,TiS(N))[i]-VGW(p,N,KeS(N,f)[i]))*N # In other words you maximize the expected savings from exiting the casino # sooner than you would with timid play, take away the expected extra loss # due to your impatience and unwillingness to play the boring timid # strategy. BestKellyFactor:=proc(N,i,p,TimeIsMoney,resolution) local f,champ,rec,challenge: f:=resolution: champ:=f: rec:=0: while f <= 1 do: challenge:= (VGD(p,N,TiS(N))[i]-VGD(p,N,KeS(N,f))[i])*TimeIsMoney -(VGW(p,N,TiS(N))[i]-VGW(p,N,KeS(N,f))[i])*N: if challenge > rec then: rec:=challenge: champ:=f: fi: f:=f + resolution: od: # Depending on choice of resolution, the case f=1 might not be tested, # so we include it here: challenge:= (VGD(p,N,TiS(N))[i]-VGD(p,N,KeS(N,1))[i])*TimeIsMoney -(VGW(p,N,TiS(N))[i]-VGW(p,N,KeS(N,1))[i])*N: if challenge > rec then: rec:=challenge: champ:=1: fi: champ, 1.0*rec: end: # Write a procedure WhatKellyCharges(p,N,i) # that finds out how much is TimeIsMoney when f is # Kelly's recommended value f=2*p-1 # WhatKellyCharges(p,N,i,tol:=0.001): p, N, and i are # given as above. tol is the tolerance. If the user # does not specify a tolerance, it gets the default # value 0.001. Estimates TimeIsMoney for which # f is 2*p-1 within tol of the true value. # This is done by finding an approximate zero of # BestKellyFactor(N,i,p,TimeIsMoney,resolution)[1]-(2*p-1) # using Bisection method. We take the resolution:=tol/10. # Note that we can assume TimeIsMoney<=1 since that is where # the timid strategy breaks down and so # BestKellyFactor ceases to work properly. WhatKellyCharges:=proc(p,N,i,tol:=0.01) local x0,x1,x,y0,y1,y,resolution,kelly,j: kelly:=2*p-1: x0:=0: x1:=1: j:=0: resolution:=1/200: y0:=BestKellyFactor(N,i,p,x0,resolution)[1]-kelly: y1:=BestKellyFactor(N,i,p,x1,resolution)[1]-kelly: if y0*y1 > 0 then FAIL: fi: while abs(y0)>tol and j<20 do x:=(x1-x0)/2: y:=BestKellyFactor(N,i,p,x,resolution)[1]-kelly: if y<=0 then x0:=x: y0:=y: elif y>0 then x1:=x: y1:=y: fi: j:=j+1 od: x0,y0,j: end: # Since we are using a discrete method to find f, # we cannot always find the desired zero to within # tolerance. Using a finer resolution could work, # but then the procedure takes a very long time to run. # We limited the number of iterations to 20 and # took the best result so obtained. # With this in mind, we have the following results: # WhatKellyCharges(11/20,20,10) = 0.0208 where f is # within 1/25 of the Kelly recommendation # WhatKellyCharges(11/20,40,32) = 0.0417 where f is # within 7/200 of the Kelly recommendation # Old Stuff #C24.txt: April 22, 2013, concluding Gambling theory #(How to gamble if you're in a hurry) Help:=proc(): print(`TiS(N), BoS(N), KeS(N,f) `): print(`BestSDD(N,i,p,T) , OptimalGamblingStory(N,i,p,T)`) : end: #TiS(N): the timid strategy TiS:=proc(N) [1$(N-1)]:end: #BoS(N): the bold strategy BoS:=proc(N) local i: [seq( min(i,N-i), i=1..N-1)]:end: #KeS(N,f): The Kelly Strategy with fraction f of your capital KeS:=proc(N,f) local i: [seq( min(ceil(i*f),N-i) ,i=1..N-1)]: end: #BestSDD(N,i,p,T): How much should you bet if you #currently have i dollars, but need to come up #with N dollars in <=T rounds (or else you would die) #and your chance of winning one round is p #Outputs the stake that would maximize your chances #of staying alive, followed by your chance of doing so BestSDD:=proc(N,i,p,T) local x,champ, rec,cand: option remember: if i=N then RETURN(GiveTheSharkTheMoney, 1): fi: if i=0 then RETURN(Pray, 0): fi: if T=0 then RETURN(Pray,0): fi: champ:=Pray: rec:=0: for x from 1 to min(i,N-i) do cand:=(1-p)*BestSDD(N,i-x,p,T-1)[2]+p*BestSDD(N,i+x,p,T-1)[2]: if cand>rec then champ:=x: rec:=cand: fi: od: champ,rec: end: #slightly edited and embellished after class OptimalGamblingStory:=proc(N,i,p,T) local t,jake,coin,cc: print(``): print(`You are at a casino where a single round is winning with`): print(`probability `, p ): print(`You currently have`, i, `dollars, but `): print(`you must come up with `, N,` dollars, that you owe to a loan shark`): print(`before `, T, `rounds are up `): print(`or else the loan shark would shoot you.`): cc:=i: t:=T: while cc>0 and cc0 do jake:=BestSDD(N,cc,p,t): print(`You currently have`, cc, `dollars `): if jake[2]=0 then print(`Pray for your life `): break: fi: print(`Your best chance of staying alive is`, evalf(jake[2]) ): print(`provided that you bet `, jake[1], `dollars `): coin:=LC(p): print(`Let's toss the coin `): if coin then print(`Yea! You won the stake of`, jake[1], `dollars! `): cc:=cc+jake[1]: else print(`Oh no! You lost the stake of`, jake[1], `dollars! `): cc:=cc-jake[1]: fi: t:=t-1: od: if cc=1 and L[i]<=min(i,N-i)),i=1..N-1)}<>{true} then RETURN(FAIL): fi: ca:=s: du:=0: while ca>0 and ca=N), du: end: #ManyGGgames(p,s,N,HowMany,L): repeats GGgame(p,s,N,L) #HowMany times, and returns the percentage of times #the gambler won, and the average duration if it won, #followed by the average duration if it lost followed #by the average duration ManyGGgames:=proc(p,s,N,HowMany,L) local DUw, DUl, DU, W,i,jake: W:=0: DUw:=0: DUl:=0: DU:=0: for i from 1 to HowMany do jake:=GGgame(p,s,N,L): if jake[1] then W:=W+1: DUw:=DUw+jake[2]: DU:=DU+jake[2]: else DUl:=DUl+jake[2]: DU:=DU+jake[2]: fi: od: if W=0 or W=HowMany then print(`This is never going to happen`): RETURN(FAIL): fi: evalf(W/HowMany), evalf(DUw/W), evalf(DUl/(HowMany-W)), evalf(DU/HowMany): end: #VGW(p,N,L): inputs prob. of winning one round, p #and the exit capital N, outputs the List of #size N-1 whose i-th entry is the EXACT prob. #of winning if you currently posses i dollars #following the strategy L VGW:=proc(p,N,L) local w,eq,unk,i,sol: if {seq(evalb(L[i]>=1 and L[i]<=min(i,N-i)),i=1..N-1)}<>{true} then RETURN(FAIL): fi: #w[i]: prob. of ultimately winning with current capital i #in Gambler's ruin name unk:={ seq(w[i],i=0..N)}: eq:={w[0]=0, w[N]=1, seq(w[i]=p*w[i+L[i]]+(1-p)*w[i-L[i]], i=1..N-1)}: sol:=solve(eq,unk): subs(sol, [seq(w[i],i=1..N-1)]): end: #VGD(p,N,L): inputs prob. of winning one round, p #and the exit capital N, outputs the List of #size N-1 whose i-th entry is the EXACT EXpected Duration #of winning if you currently posses i dollars #and follow strategy L VGD:=proc(p,N,L) local Du,eq,unk,i,sol: if {seq(evalb(L[i]>=1 and L[i]<=min(i,N-i)),i=1..N-1)}<>{true} then RETURN(FAIL): fi: #Du[i]: expected duration of #a gambler's ruin with max. cap N and prob. p of winning #a dolllar unk:={ seq(Du[i],i=0..N)}: eq:={Du[0]=0, Du[N]=0, seq(Du[i]=p*Du[i+L[i]]+(1-p)*Du[i-L[i]]+1, i=1..N-1)}: sol:=solve(eq,unk): subs(sol, [seq(Du[i],i=1..N-1)]): end: #IsAB(L1,L2) is the list L1 always >= list L2 IsAB:=proc(L1,L2) local i: evalb({seq(evalb(L1[i]>=L2[i]),i=1..nops(L1))}={true}): end: #AllVecsLess(L): The SET of all the vectors of length nops(L) #of pos. integers such that 1<=M[i]<=L[i] AllVecsLess:=proc(L) local n, L1, Phil,i, tom: option remember: n:=nops(L): if n=0 then RETURN({[]}): fi: L1:=[op(1..n-1,L)] : Phil:=AllVecsLess(L1): {seq(seq([op(tom),i], tom in Phil), i=1..L[n])}: end: #AllStra(N): The set of all possible strategies with max capital #N AllStra:=proc(N) local i: AllVecsLess([seq(min(i,N-i),i=1..N-1)]): end: #BestStra(p,N) BestStra:=proc(p,N) local champs,cu,S,cu1,i: S:=AllStra(N): champs:={S[1]}: cu:=VGW(p,N,S[1]): for i from 2 to nops(S) do cu1:=VGW(p,N,S[i]): if cu1=cu then champs:=champs union {S[i]}: elif IsAB(cu1,cu) then champs:={S[i]}: cu:=cu1: fi: od: champs: end: ###################old stuff from C22.txt #April 15, 2013, starting gambling theory Help22:=proc(): print(` LC(p) , GRgame(p,s,N) `): print(`ManyGames(p,s,N,HowMany), VW(p,N) , VD(p,N) `): end: #LC(p): simulating a loaded coin whose #that returns true with prob. p, p rational LC:=proc(p) local a,b,ra: a:=numer(p): b:=denom(p): ra:=rand(1..b)(): evalb(ra<=a): end: #GRgame(p,s,N):inputs prob. p of success (a rational number) #starting capital, and exit capital N (if you are a winner) #once you have 0 dollars you must leave the casino #Returns, true (false) if you won (or lost) followed #by the number of rounds it took GRgame:=proc(p,s,N) local ca,du,toss: ca:=s: du:=0: while ca>0 and ca=N), du: end: #ManyGames(p,s,N,HowMany): repeats GRgame(p,s,N) #HowMany times, and returns the percentage of times #the gambler won, and the average duration if it won, #followed by the average duration if it lost followed #by the average duration ManyGames:=proc(p,s,N,HowMany) local DUw, DUl, DU, W,i,jake: W:=0: DUw:=0: DUl:=0: DU:=0: for i from 1 to HowMany do jake:=GRgame(p,s,N): if jake[1] then W:=W+1: DUw:=DUw+jake[2]: DU:=DU+jake[2]: else DUl:=DUl+jake[2]: DU:=DU+jake[2]: fi: od: if W=0 or W=HowMany then print(`This is never going to happen`): RETURN(FAIL): fi: evalf(W/HowMany), evalf(DUw/W), evalf(DUl/(HowMany-W)), evalf(DU/HowMany): end: #VW(p,N): inputs prob. of winning one round, p #and the exit capital N, outputs the List of #size N-1 whose i-th entry is the EXACT prob. #of winning if you currently posses i dollars VW:=proc(p,N) local w,eq,unk,i,sol: #w[i]: prob. of ultimately winning with current capital i #in Gambler's ruin name unk:={ seq(w[i],i=0..N)}: eq:={w[0]=0, w[N]=1, seq(w[i]=p*w[i+1]+(1-p)*w[i-1], i=1..N-1)}: sol:=solve(eq,unk): subs(sol, [seq(w[i],i=1..N-1)]): end: #VD(p,N): inputs prob. of winning one round, p #and the exit capital N, outputs the List of #size N-1 whose i-th entry is the EXACT EXpected Duration #of winning if you currently posses i dollars VD:=proc(p,N) local Du,eq,unk,i,sol: #Du[i]: expected duration of #a gambler's ruin with max. cap N and porb. p of winning #a dolllar unk:={ seq(Du[i],i=0..N)}: eq:={Du[0]=0, Du[N]=0, seq(Du[i]=p*Du[i+1]+(1-p)*Du[i-1]+1, i=1..N-1)}: sol:=solve(eq,unk): subs(sol, [seq(Du[i],i=1..N-1)]): end: