The effect of each bet on your capital is multiplicative - your capital is multiplied by $1 + 2f$ on heads and by $1 - f$ on tails - so that your fortune at the end depends only on the total number of heads and tails, not on the order in which they occured.
In particular, for a given $f$, there will be a number $h(f)$ such that if the number of heads flipped is greater than or equal to $h(f)$, our final fortune will be greater than one billion.
import matplotlib.pyplot as plt
import numpy as np
f = 0.25
hs = range(1001)
final_amounts = [(1 + 2 * f)**h * (1 - f)**(1000 - h) for h in hs]
big_enough = [h for h, final_amount in zip(hs, final_amounts) if final_amount >= 10**9]
fig, ax = plt.subplots(figsize=(12, 8))
ax.semilogy(hs, final_amounts)
ax.axhline(1e9, color='red', linestyle='--')
ax.axvspan(min(big_enough), max(big_enough), color='green', alpha=0.2)
ax.set_xlim(0, 1000)
ax.set_ylabel('Final fortune')
ax.set_xlabel('Number of heads')
ax.set_title('{0} heads needed for f = {1}'.format(min(big_enough), f));
We can solve directly for the number of heads needed for a given $f$.
$$ (1 + 2f)^h (1 - f)^{1000 - h} = 10^9$$$$ h \log_{10} (1 + 2f) + (1000 - h) \log_{10} (1 - f) = 9 $$$$ h = \left \lceil{\frac{9 - 1000 \log_{10} (1 - f)}{\log_{10} (1 + 2f) - \log_{10} (1 - f)}}\right \rceil $$def min_h_greater_given_f(f, mult=2):
''' Returns the minimimum number of heads needed to end up with at least
10^9 pounds after 1000 flips if your capital is multiplied by
(1 + mult * f) for each head and by (1 - f) for each tail. If even
1000 heads are not enough, returns 1001.
'''
if f == 0:
# Impossible to end up above 1 if wagering 0.
h = 1001
elif f == 1:
# Need to never get a tail if wagering everything.
h = 1000
else:
exact = (9 - 1000 * np.log10(1 - f)) / (np.log10(1 + mult * f) - np.log10(1 - f))
h = min(1001, np.ceil(exact))
return h
An optimal $f$ minimizes the number of heads needed. Notice that because the number of heads is discrete, a range of $f$ values will be optimal.
fs = np.linspace(0, 1, 1000)
hs = [min_h_greater_given_f(f) for f in fs]
min_h = min(hs)
optimal_fs = [f for f, h in zip(fs, hs) if h == min_h]
fig, (full_ax, zoom_ax) = plt.subplots(1, 2, figsize=(12, 6))
full_ax.plot(fs, hs)
full_ax.set_ylim(0, 1001)
full_ax.set_xlabel('f')
full_ax.set_ylabel('Heads needed')
zoom_ax.plot(fs, hs)
zoom_ax.set_ylim(min_h - 10, min_h + 10)
zoom_ax.set_xlim(min(optimal_fs) - 0.1, max(optimal_fs) + 0.1);
zoom_ax.set_xlabel('f')
zoom_ax.set_ylabel('Heads needed');
The probability of achieving a given number of heads is given by a Binomial(1000, 0.5) distribution. We are interested in the probability of getting greater than or equal to a given number of heads. We can compute this ourselves from the formula...
from scipy.misc import comb
hs = range(1002)
binomial_coefficients = [comb(1000, h) for h in hs]
pdf = np.array(binomial_coefficients) / 2.**1000
p_geq = np.cumsum(pdf[::-1])[::-1]
fig, (pdf_ax, p_geq_ax) = plt.subplots(1, 2, figsize=(12, 6))
pdf_ax.plot(hs, pdf)
pdf_ax.axis('tight')
p_geq_ax.plot(hs, p_geq)
p_geq_ax.axis('tight');
... or use scipy
's built-in binomial distribution object.
from scipy.stats import binom
p_geq_alternate = [binom(1000, 0.5).sf(h - 1) for h in range(1002)]
print np.allclose(p_geq, p_geq_alternate)
True
The solution is the probability of getting at least the minimum number of heads.
print p_geq[min_h]
0.999992836186
The problem as stated multiplies your wager by 2 on a win. What is the smallest this factor can be while still leaving the probability of ending up above one billion greater than 0.5, assuming that you play with an optimal f given the value of the factor?
def best_p_given_mult(mult):
''' For a given value mult of the factor, returns the probability
of ending up a billionaire given optimal f.
'''
fs = np.linspace(0, 1, 1000)
min_h = min(min_h_greater_given_f(f, mult) for f in fs)
best_p = p_geq[min_h]
return best_p
full_mults = np.linspace(1, 2, 100)
full_ps = [best_p_given_mult(mult) for mult in full_mults]
zoom_mults = np.linspace(1.48, 1.52, 100)
zoom_ps = [best_p_given_mult(mult) for mult in zoom_mults]
fig, (full_ax, zoom_ax) = plt.subplots(1, 2, figsize=(16, 8))
full_ax.plot(full_mults, full_ps)
full_ax.axis('tight')
full_ax.set_ylabel('P(billionaire) given optimal f')
full_ax.set_xlabel('Factor on heads')
zoom_ax.plot(zoom_mults, zoom_ps)
zoom_ax.axis('tight')
zoom_ax.set_ylabel('P(billionaire) given optimal f')
zoom_ax.set_xlabel('Factor on heads');
We know that p < 0.5 if the factor is 1 and p > 0.5 if the factor is 2. We can repeatedly bisect this interval, maintaining the invariant that the left point of each new interval always produces a p < 0.5 and the right point always produces a p > 0.5.
low, high = 1., 2.
num_digits = 15
while (high - low) > 1 * 10**-num_digits:
mid = (low + high) / 2
mid_p = best_p_given_mult(mid)
if mid_p < 0.5:
low = mid
else:
high = mid
print 'A factor of {0:0.{num_digits}f} is high enough, but {1:0.{num_digits}f} is not.'.format(high, low, num_digits=num_digits)
A factor of 1.504676213769604 is high enough, but 1.504676213769603 is not.