I think your n_subjects is too low. You need that to be high enough or you'll miss those low-probability winners that bring up the average.
Looking at the Julia code, I think what he is doing wrong is making all wins worth $.50 and all losses worth $.40, but the bet computes a win or loss based on your current wealth, not your starting wealth. His formula would work if you were always betting $1 no matter what your bankroll was, but that isn't what the actual post stipulates.
If you don't trust the Julia code, try running with the same parameters in Python.
My code is right up there and you can run it. You can even just run the OP's notebook that he provided but increase the number of trials. Change the "num_flips_per_sim" parameter he provides in cell 6 to anything over 500 and you will always get sum(count_lose_capital) == everyone.
The appropriate response to that is introspection, not repetition.
More precisely there is a finite time after which no-one ever passes above $0.0000000000000001.
That is a mathematical theorem.
This doesn’t depend on the number of test subjects, and you can add as many zeroes as you want.
Therefore in the long run the mean outcome is 0.
Forgive me if I have misinterpreted what you are are trying to say.
Edit: I’ve just realized that I have indeed missed your point.
It’s a mathematical theorem. (I would be curious to see a proof of your theorem, by the way.)
t=0 mean(w) = 1
t=1 mean(w) = 1/2*1.5 + 1/2*0.6 = 1.05
t=2 mean(w) = 1/4*1.5*1.5 + 1/2*1.5*0.6 + 1/4*0.6*0.6 = 1.1025
....
t mean(w) = 1.05^t
As for the proof of my theorem, By taking logarithms, the process becomes an additive random walk with negative drift (log 1.6 + log 0.5 < 0). This is well known to converge to negative infinity almost surely. After exponentiating to undo the logarithm, this is exactly the statement I made.
It does not matter how many test subjects there are ( as long as there’s finitely many) because, informally speaking , you can just wait for each of them to become irrevocably bankrupt in turn.
- for a fixed sample size we can find a time large enough that the probability of the sample mean being above $1 is as low as we want
- for a fixed time we can find a sample size large enough that the probability of the sample mean being below $1 is as low as we want
- when both the sample size and the horizon grow without limit which effect dominates will depend on how we make it happen
Adding "almost surely" to "everyone goes bankrupt and will never recover" or "there is a finite time after which no-one ever passes above $0.0000000000000001" is a subtle change but it's enough to allow for someone to go to infinity with infinitesimal probability.
This is why the distribution mean can grow exponentially, it wouldn't be possible if the everyone and no-one in those quotes were strictly true.
Just to confirm, I am using 'almost surely' in the technical sense, which means 'with probability 1.'
Consider the following statement:
If you keep flipping a fair coin every day, it is almost sure that after some day you will have gotten a tails.
This is the same 'almost surely' that I am referring to.
The point was that you didn't specify "almost surely" previously, that's why I asked for a proof to understand what did you mean exactly when you said that "everyone goes bankrupt and will never recover" and "or "there is a finite time after which no-one ever passes above $0.0000000000000001".
The mean of a random variable that is close to zero is close to zero, the mean of a random variable that is almost surely close to zero can be anything.
After 500 trials, you need 279 heads to stay above $1 net wealth. 1.5^278 + 0.6^222 = 0.50 and 1.5^279 + 0.6^221 = 1.26, so that's your breakeven point. The probability of getting at least 279 heads in 500 coin flips is 0.005364, so with 1,000 participants, you expect to see about 5 still above water.
At 1000 trials, the breakeven point becomes 558 and the probability of getting at least that many heads in 1000 flips is 0.00013614. So the expected number of people who stay above water in a pool of 1000 participants is 0. Out of 1,000,000, it is 13, so you're right, there are some, but at that point it's not nearly enough and we're not sampling the ones whose wealth is enough to actually bring the mean back up, so it keeps trending to 0 in any sample of a practical trial size.
This is a pretty interesting property of this problem, really. It's not related to ergodicity, but just the relative proportion of probability mass represented by above 1 and below one itself trending asymptotically toward 0 even though the analytical expectation trends toward infinity. I don't know that there is even a word for that, but seemingly which of those moves faster toward its limit would determine what sample ensemble average you really see when the number of realized states is far less than the number of possible states.
This probably has some implications for Pascal's Mugger type problems in decision theory. If some course of action has potentially infinite future payoff and destroys expected utility calculations because of that, but the expected number of possible universes in which a positive outcome happens at all trends toward 0 faster than the expectation trends toward infinity, that gives a decision rule. In this specific case, don't take this bet, at least not in an indefinitely repeating form.
A shorter run (say 100 steps) would be more likely to capture enough realisations to produce a reasonable estimate. You could assess this behaviour yourself, for very low step numbers, by calculating the variability in a sampled ensemble average, relative to the exhaustive (i.e. true) ensemble average.
This particular problem is another consequence of the properties dynamical system being examined, but not quite the same as the issues caused by its non-ergodicity.
import numpy as np
import itertools
from matplotlib import pyplot as plt
def ensemble_mean(outcomes):
# Assume we are given a (K, T) array of outcomes, and compute the ensemble average
# for T+1 time steps, starting with 1 wealth.
K, T = outcomes.shape
X = np.ones((K, T+1), dtype=np.float64)
X[:, 1:] = np.where(outcomes, 1.5, 0.6)
Z = np.cumprod(X, axis=1)
return Z.mean(axis=0)
time_steps = 20
all_outcomes = np.array(list(itertools.product([0, 1], repeat=time_steps-1)))
exhaustive_mean = ensemble_mean(all_outcomes)
ensemble_size = 100
ensemble_samples = 10000
ensemble_means = np.zeros((time_steps, ensemble_samples))
for i in range(ensemble_samples):
print(i)
# generate ensembles as though we were sampling (i.e. with replacement)
J = np.random.choice(all_outcomes.shape[0], size=ensemble_size, replace=True)
ensemble_means[:, i] = ensemble_mean(all_outcomes[J, :])
plt.hist(ensemble_means[-1], bins=1000, histtype='step')
plt.axvline(exhaustive_mean[-1])
plt.title("Modal sampled ensemble mean is below true ensemble mean")
plt.show()Sounds like a good strategy whether in Vegas or Wall Street.