Let’s kick it off with an important thought experiment: what % crime rate among blond-haired males would make it ok to put them all into Trolley A, heading the opposite way toward Trolley B full of innocent baby goats with wave function Ψ = … just kidding!
I’ll share a math problem that has been kicking my ass lately.
A gambler with bankroll b will repeatedly bet on coinflips against an infinitely bankrolled casino. The gambler will wager $1 on each flip. Their probability of winning a flip is p and a winning flip pays $w, such that:
p ≠ 1/2
w ≠ 1
pw > 1-p, which is to say the gambler’s EV is positive.
If the gambler will play forever unless they bust, what is the probability of busting?
There are two ways to define ruin. If the gambler can’t decrease their stake, then busting is when b falls below $1. If the gambler can bet smaller, then busting is reaching $0.
What has kicked my ass is trying to solve it combinatorially with the help of lattice path visualizations; for now I’ve given up. There is a paper whose authors who took that approach and came up with a solution that only works for integer-valued w. There is another paper whose authors present a technique that sounds like it would apply to this problem and work for all rational w, but I haven’t digested it enough to know.
Most other authors only try to find estimates, most of which aren’t particularly good. One author, Katriel, seems to have found an exact solution for rational w where ruin is b<1; I haven’t tried it out yet.
When ruin is defined as reaching zero, there is an accurate upper-bound approximation which, to my knowledge, is original to Chen & Ankenman in The Mathematics of Poker. It’s also easy to program and fast to execute. I’ve coded it in Julia:
function rorexp(p::Vector{Float64}, w::Vector{Int64}, b::Number, tol=2^-10)
minimum(w)≥0 && return 0
sum(p.*w)≤0 && return 1
len = length(p)
nearone = 1+tol # Technically (1-tol, 1+tol), but >1 gives a more accurate answer.
goal = 1+tol/2
step = -1/16
x = step
prevx = step
prevdiff = Inf
zerohi = true
foundlow = false
while true
sum = 0.0
for k=1:len begin sum += p[k]*exp(w[k]*x) end end
diff = abs(sum-goal)
if diff >= prevdiff
if zerohi && step>0
zerohi = false
elseif !foundlow && step<0
foundlow = true
!zerohi && begin step/=2 end
end
x = prevx + (step /= -2)
elseif 1 < sum < nearone
return exp(b*x)
else
zerohi = (prevdiff==Inf)
prevx = x
prevdiff = diff
foundlow ? x+=(step/=2) : x+=step
end
end
end
function rorexp(p::Vector{Float64}, w::Vector{Float64}, b::Number, tol=2^-10)
minimum(w)≥0 && return 0
sum(p.*w)≤0 && return 1
len = length(p)
nearone = 1+tol
step = -1/16
x = step
while true
sum = 0.0
for k in 1:len begin sum += p[k]*exp(w[k]*x) end end
if (sum>nearone && step<0) || (sum<1 && step>0)
step /= -2
elseif 1 < sum < nearone
return exp(b*x)
end
x += step
if x>=0
x += (step /= -2)
while x>=0 begin x+=(step/=2) end end
end
end
end
(You need both function methods. One only accepts integer values into the w vector and the other only accepts decimals. The first method’s solver algorithm is a little more efficient, but doesn’t work with decimal payouts. Writing my own solver algo pretty much amounts to reinventing the wheel, but it was more fun than using whatever package I could have used.)
Example: 70% chance of winning each flip, laying 3:2 odds, starting with $2. This is the same as having $6 and betting $3 per flip, so it’s easier to type:
rorexp([.7, .3], [2, -3], 6)
Result: 34.83%
It works for games with more than two outcomes as well, for instance the SNG example in MOP:
rorexp([.12, .13, .13, .62], [785, 285, 185, -215], 1000)
Result: 54.32%
If you can’t move down in stakes when below one unit, there is a fast brute-force algorithm called Gauss-Seidel. It iterates a recurrence relation until the desired precision is reached:
function seidelgauss(p::Vector{Float64}, w::Vector{Int64}, b::Int64, T::Int64, tol=10^-9)
tol = max(tol, 4*eps())
dig = ceil(Int64, abs(log(10,tol))-2) # This many digits will be accurate.
pwlen = length(p)
maxloss = abs(minimum(w))
r = zeros(Float64, maximum(w)+T) # T is the dollar target at which we assume RoR ≈ 0
r[1 : maxloss] .= 1.0
while true
prev = r[b+1]
for j=(maxloss+1):T
r[j] = 0.0
for k=1:pwlen begin r[j] += p[k]*r[j+w[k]] end end
end
diff = abs(r[b+1] - prev)
diff<tol && break
end
return round(r[b+1]; digits=dig)
end
Same two examples:
seidelgauss([.7, .3], [2, -3], 6, 1000)
seidelgauss([.12,.13,.13,.62],[785,385,185,-215],1000,20000)
Results: 42.75% and 58.48%
Finally, here is a Monte Carlo sim applying to more variants of the problem:
# For better accuracy, use integer W & L when practical.
# For no grossTarget (ie to play against an infinite bankroll)), put Inf.
# Movedown == whether or not you have the ability to move down in stakes if your bankroll falls below the max possible loss.
function rorsim(p::Vector{Float64}, w::Vector{T}, initCapital::Number, grossTarget::Number,
sims::Int64, movedown=true, gameLimit=Inf, showprogress=false) where T<:Number
eternal = (gameLimit==Inf)
grossTarget==Inf && eternal && error("need a grossTarget or a gameLimit")
len = length(p)
maxloss = abs(minimum(w))
busts = 0
showprogress && begin qtr1, qtr2, qtr3 = ceil(Int64,sims/4), ceil(Int64,sims/2), ceil(Int64, 3*sims/4) end
for n=1:sims
!eternal && begin g=0 end
b = initCapital
while b<grossTarget && (b ≥ maxloss || (movedown && b>0)) && (eternal || g<gameLimit)
r = rand()
sum = p[1]
for k=1:len
if r<=sum || k==len
(b≥maxloss || w[k]≤0) ? b+=w[k] : b += w[k]*b/maxloss
break
else sum += p[k+1] end
end
!eternal && begin g+=1 end
end
(b≤0 || (b<maxloss && !movedown)) && begin busts+=1 end
showprogress && (n==qtr1 || n==qtr2 || n==qtr3) && begin println("n = ",n) end
end
return busts/sims
end
Example:
rorsim([.7, .3], [2,-3], 6, 100, 1000000)
Result: 34.28%
Side note: apparently there’s a plugin discourse-math which, if enabled, allows the use of LaTeX (via MathJax) and KaTeX. If not too much trouble, that would probably come in handy for this thread and potentially the Poker Hands thread.