#April 15, 2013, starting gambling theory
Help:=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: