I've got a game that can be described by flipping coins over 7 days. Every day I get a number of coins described by a Poisson distribution. Each coin flips heads with probability $p$, and if I succeed, I get to keep the coin for the next day and flip again. Every Successful coin flip, I get to take a percentage of the reward pot which starts off with $1000. What is my expected reward for playing this game?

For example: On Monday, I receive 2 coins due to the poissson spawn rate. That day I flip the first coin. It comes up tails so I lose it. The second coin is heads. I get my 15% of the pot (so \$150) the first day and get to keep the coin for Tuesday. On Tuesday, I receive 0 coins, but I still have my coin from Monday. I get to flip that coin and it comes up heads again. My payoff this time is 15% of 1000-150=850 so \$127.5. If I receive no more coins, or all my coin tosses are tails for the other 5 days, my total reward is \$150 + \$127.5 = \$177.5

Im trying to solve this piece by piece and verify it with simulations, however, I cant get the two to align. Here's what I've got so far:

Let $\lambda$ be the parameter for the poissson number of coins per day. I would expect to receive a total of $7\lambda$ total coins over the course of a week. Each coin follows a geometric distribution, so it should last on average $1/p$ days. Therefore, I would expect to have a total of $7\lambda/p$ coin flips for the week. This is where my math doesn't align with my simulation. I know my math should be an over estimate since Im ignoring the truncation that happens when a coin spawns on the last day. If I receive a reward from it, it would only be for one day instead of $1/p$ days. However, my simulations give higher values than my math.

For reference, here is my code:

import numpy as np

trials = int(1e5)
spawn_rate = 0.08147
days = 7
persist = 0.75

toss_count = 0 
for _ in xrange(trials):
    # A vector of 7 days, with each element recording the number of coins spawned that day
    coins = np.random.poisson(spawn_rate, days)

    # Convert the spawned coins into the number of successful tosses per day
    tosses_per_day = [0] * days
    for idx,num_coins in enumerate(coins):
        for coin in range(num_coins):
            # For each coin today, add to the toss count as long as this coin persists
            offset = 0
            while np.random.uniform() < persist:
                    # The coin landed heads, so add it to the toss count for this day
                    tosses_per_day[idx + offset] += 1
                except IndexError:
                    # Do not count when a coin could continue past 7 days
                offset += 1
    toss_count += np.sum(tosses_per_day)

# return the average heads per week
print toss_count / float(trials)

Here is another program that produces the same results as above. This leads me to believe my original expression of $7\lambda/p$ is either wrong or the assumption of extending past a week does not hold well enough for the time frame.

import numpy as np

persist = 0.75
spawn_rate = 0.08147

def day(yesterday_coins, heads_to_date, total_days):
    if total_days <= 0:
        return heads_to_date
        today_coins = yesterday_coins + np.random.poisson(spawn_rate)
        tomorrow_coins = today_coins
        heads = 0;
        for _ in xrange(today_coins):
            if np.random.uniform() < persist:
                heads += 1
                tomorrow_coins -= 1
        return day(tomorrow_coins, heads_to_date + heads, total_days -1)

trials = int(1e5)
count = 0
for _ in xrange(trials):
    count += day(0,0,7)
print count / float(trials)

Keep in mind that Im looking for a mathematical expression for the expected value rather than a program to compute it for me.

  • $\begingroup$ +1 for validating math against code. However, It is unclear how the code computes what you described. Could you clean it up? (why xrange(2); why sum / 1e5; don't use sum as variable; replace xrange -> range; is p == 1 - persist?) $\endgroup$
    – psarka
    Jan 3, 2018 at 16:38
  • $\begingroup$ @psarka updated. The probability the coin lands heads is denoted with persist because that is also the probability I will keep the coin the next day. $\endgroup$ Jan 3, 2018 at 17:02
  • $\begingroup$ I feel like I'm missing something. If I run your program I get ~1.07. If I do the calculation days*spawn_rate/(1-persist) I get ~2.28 > 1.07 Also I'm not sure I trust your daily calculation. You might consider working on clarifying that! $\endgroup$ Jan 3, 2018 at 19:39
  • $\begingroup$ @KitterCatter You're right. I was dividing by persist as opposed to (1-persist). However, I thought the math would be a better approximation to the simulation which is why Im not convinced its correct. How can I approach this so these two methods match? $\endgroup$ Jan 3, 2018 at 19:55
  • $\begingroup$ I see two routes: 1) Check/refactor your code, comparing to a bad baseline won't really clarify things. 2) Make more realistic approximations: For example you know 7 days doesn't work because of carryover. Why don't you try truncating the 7th day? If that has a big effect maybe your approximations aren't that realistic afterall? $\endgroup$ Jan 3, 2018 at 20:02

1 Answer 1


I tried taking your description and turning it into some code that simulates a day, a week, and a number of trials at a time. This tracks number of coin flips per week to expect.

I should note that a correction to your approximation of $7\lambda/p$ to $6 \lambda/p +\lambda$ is more than 10% for the given parameters so I wouldn't expect that this simplification to hold all that well. I would expect your simplification would work better in the realm where persist is quite low, say 10%. And in fact that is the case:

0.62 vs 0.63

import numpy as np

def simulate_a_day(num_persist_coins, spawn_rate, persist_prob):
    num_coins = int(np.random.poisson(spawn_rate)+num_persist_coins)
    persist_coins = np.sum([i < persist_prob for i in np.random.uniform(0,1,num_coins)])
    return (num_coins,persist_coins)

def simulate_a_week(spawn_rate, persist_prob, num_days):
    running_total = 0
    num_persist_coins = 0
    for _ in range(num_days):
        num_coins, num_persist_coins = simulate_a_day(num_persist_coins,
        running_total += num_coins
    return running_total

def simulate_multiple_trials(num_trials, spawn_rate, persist_prob, num_days):
    weekly_total = 0
    for _ in range(num_trials):
        weekly_total += simulate_a_week(spawn_rate, persist_prob, num_days)
    return weekly_total/num_trials

if __name__ == '__main__':


    print("compared to:")

Since the problem description changed I will add a note: The maths that you do to get $n_{days} \lambda/p$ is for the number of flips. If you want the number of heads to compare with you would need to multiply by the probability of a heads: $n_{days}(1-p) \lambda/p$. Again this will work poorly when a coin is expected to persist at greater probability. A first correction can be obtained by fixing coins spawned on day 7 and only counting the heads, while letting other days' coins persist until they get a tails: \begin{equation} (n_{days}-1)\lambda (1-p)/p + (1-p)\lambda = (n_{days})\lambda (1-p)/p - \lambda(1-p)^2/p \end{equation} Thus as $p$ approaches 0 the relative error increases and you should not expect the calculation to match simulation. (For the case of p=1/4 the relative error is about 10% so you really shouldn't expect your calculation to hold up all that well.)

I think this series of approximations should lead you to a geometric series which is solvable: $$ \left(\sum_i i (1-p)^{1-i}\right) \lambda (1-p)^7 = \lambda(1-p)^{n_{days}+1}\frac{(n_{days}-(n_{days}+1)(1-p))(1-p)^{-n_{days}}+(1-p)}{p^2} $$

Using this sum I get results consistent with simulation.

  • $\begingroup$ Im looking for an expression for the expectation rather than simulation code. Im only using code to verify my math as I go along. Sorry if that wasn't clear in my original question. Having said that, I do like your code. The only issue is that every day, coins must pass their persistence check for the day they are spawned as well. So your simulate_day function should return persistent_coins for both values. $\endgroup$ Jan 3, 2018 at 20:30
  • $\begingroup$ Then your math doesn't work since you are calculating number of flips. $\endgroup$ Jan 3, 2018 at 20:35
  • $\begingroup$ If you want number of heads it should be $\lambda*persist*n_{days}/(1-persist)$ $\endgroup$ Jan 3, 2018 at 20:42
  • $\begingroup$ Can you explain how you arrived at the first correction? where did the $-\lambda(1-p)^2/p$ come from? $\endgroup$ Jan 3, 2018 at 20:57
  • $\begingroup$ Day 7 coins you only count heads, all other days are allowed to persist until you get a tails. $\endgroup$ Jan 3, 2018 at 21:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.