# Experimental Math rocks! # Homework 24 # Phil Benjamin # Free to post ... HelpHW24 := proc() print(`Problem 1: OptimalGamblingTerse(N,i,p,T)`); print(`Problem 2: ManyOptimalGambling(N,i,p,T,K)`); print(`Problem 3: [experiment: problematic result??]`); print(`Problem 4: BestKellyFactor(N,i,p,TimeIsMoney,resolution)`); print(`Problem 5: WhatKellyCharges(p,N,i)`); end proc: # Problem 1 # OptimalGamblingTerse(N,i,p,T) # Return: # true (alive) or false (dead) # list of the intermediate probabilities OptimalGamblingTerse := proc(N,i,p,T) local t,BSDD,BSDDbet,BSDDprob,coin,cc; local probList; # Initialize cc := i; t := T; probList := []; while cc>0 and cc0 do BSDD := BestSDD(N,cc,p,t); BSDDbet := BSDD[1]; BSDDprob := BSDD[2]; probList := [op(probList),BSDDprob]; if BSDDprob=0 then break; end if; coin := LC(p); if coin then cc := cc + BSDDbet; else cc := cc - BSDDbet; end if; t := t-1; end do; return [evalb(cc 1/2) # TimeIsMoney - cost of playing a single round # resolution - small positive number (say 1/200) # Try out all possible Kelly strategies, that is: # f = resolution,2*resolution,..,1 # to maximize: # (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 BestKellyFactor := proc(N,i,p,TimeIsMoney,resolution) local f, bestFact, gameValue, bestValue; # Initialize bestFact := -1; bestValue := -1; for f from resolution to 1 by resolution do gameValue := TimeIsMoney * (VGD(p,N,TiS(N))[i] - VGD(p,N,KeS(N,f))[i]) - N * (VGW(p,N,TiS(N))[i] - VGW(p,N,KeS(N,f))[i]); if gameValue > bestValue then bestValue := gameValue; bestFact := f; end if; end do; return bestFact; end proc: # Sample run: BestKellyFactor(20,5,3/5,1,1/20) --> 19/20 # Problem 5 # WhatKellyCharges(p,N,i) # Find TimeIsMoney when f=2*p-1 # Experiments: # WhatKellyCharges(11/20,20,10) --> 0.04 # WhatKellyCharges(11/20,40,21) --> 0.03 # WhatKellyCharges(11/20,100,40) --> 0.01 # Note: original values # WhatKellyCharges(11/20,400,150) did not terminate # Strategy: Use formula from Problem 4 with f=2*p-1 but # TimeIsMoney unknown (solve for it) # Assume Kelly is very fair and wants to make no profit! WhatKellyCharges := proc(p,N,i) local f,TimeIsMoney; f := 2*p - 1; TimeIsMoney := N * (VGW(p,N,TiS(N))[i] - VGW(p,N,KeS(N,f))[i]) / (VGD(p,N,TiS(N))[i] - VGD(p,N,KeS(N,f))[i]); return evalf(TimeIsMoney); end proc: ##################################################### # C24.txt #C24.txt: April 22, 2013, concluding Gambling theory #(How to gamble if you're in a hurry) HelpC24:=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: