How to get the detailed results (probability tree) of executing a bernoulli experiment a large number of times

Question

Supposing the following experiment : Execute the same bernoulli trial (with the probability of success P) N number of times

I need the following information : All the possible sequences of success/failure with its probability to happen.

Example : A Bernouilli experiment with a probability of success P = 40% executed 3 times would yield the following results (S is a success, F is a failure) :

FFF 0.216

SFF 0.144

FSF 0.144

SSF 0.096

FFS 0.144

SFS 0.096

FSS 0.096

SSS 0.064

I tried to bruteforce it to obtain the results, but it chokes rapidly with only N = 25, I get an OutOfMemoryException...

using System;
using System.Linq;
using System.Collections.Generic;
using System.Text.RegularExpressions;

namespace ConsoleApplication
{
    class Program
    {
        static Dictionary<string, double> finalResultProbabilities = new Dictionary<string, double>();

        static void Main(string[] args)
        {
            // OutOfMemoryException if I set it to 25 :(
            //var nbGames = 25;
            var nbGames = 3;
            var probabilityToWin = 0.4d;

            CalculateAverageWinningStreak(string.Empty, 1d, nbGames, probabilityToWin);

            // Do something with the finalResultProbabilities data...
        }

        static void CalculateAverageWinningStreak(string currentResult, double currentProbability, int nbGamesRemaining, double probabilityToWin)
        {
            if (nbGamesRemaining == 0)
            {
                finalResultProbabilities.Add(currentResult, currentProbability);
                return;
            }

            CalculateAverageWinningStreak(currentResult + "S", currentProbability * probabilityToWin, nbGamesRemaining - 1, probabilityToWin);
            CalculateAverageWinningStreak(currentResult + "F", currentProbability * (1 - probabilityToWin), nbGamesRemaining - 1, probabilityToWin);
        }
    }
}

I need to be able to support up to N = 3000 in a timely manner (obtaining the result in less than 3 seconds for any P)

Is there a mathematical way to do this optimally?

Answer 1

Since you're interested in the length of the longest winning streak, at any point in the trials there are only two relevant facts about the history: (1) how long the longest winning streak is so far (m) and (2) how long the current winning streak is (k). The initial state is (0, 0). The transitions are (m, k) -> (max(m, k + 1), k + 1) on a win, and (m, k) -> (m, 0) on a loss. Once we know the probability of all of the final states, we can average.

EDIT: the revised version of this code has an optimization to cut down on the computation needed for very unlikely events. Specifically, we ignore all states with a streak of length greater than or equal to some cutoff. Given an absolute error budget abserr, we determine that we can exclude probability mass up to abserr / n from the final expected value computation, since the longest streak possible is n. There are less than n places to start a streak, and at each place, the probability of a streak of length cutoff is pwin**cutoff. Using a crude union bound, we need

n * pwin**cutoff <= abserr / n
pwin**cutoff <= abserr / n**2
log(pwin) * cutoff <= log(abserr / n**2)
cutoff >= log(abserr / n**2, pwin),

where the last inequality flips direction because log(pwin) < 0.

This revised code runs in less than three seconds despite a poor quality evaluator (i.e., the Python 3 interpreter on 2014-era hardware).

import math


def avglongwinstreak(n, pwin, abserr=0):
    if n > 0 and pwin < 1 and abserr > 0:
        cutoff = math.log(abserr / n**2, pwin)
    else:
        cutoff = n + 1
    dist = {(0, 0): 1}
    for i in range(n):
        nextdist = {}
        for (m, k), pmk in dist.items():
            winkey = (max(m, k + 1), k + 1)
            if winkey[0] < cutoff:
                nextdist[winkey] = nextdist.get(winkey, 0) + pmk * pwin
            losskey = (m, 0)
            nextdist[losskey] = nextdist.get(losskey, 0) + pmk * (1 - pwin)
        dist = nextdist
    return sum(m * pmk for ((m, k), pmk) in dist.items())


print(avglongwinstreak(3000, 0.4, 1e-6))

Answer 2

Here's a different approach, which is exact and fast enough being only quadratic. The expected value of the longest win streak is equal to

 n
sum Pr(there exists a win streak of length at least k).
k=1

We reason about the probability as follows. Either the record opens with a length-k win streak (probability pwin**k), or it opens with j wins for some j in 0..k-1 followed by a loss (probability pwin**j * (1 - pwin)), on which condition the probability is equal to the probability of a length-k win streak in n - (j + 1) tries. We use memoization to evaluate the recurrence that this logic implies in pwinstreak; the faster version in fastpwinstreak uses algebra to avoid repeated summations.

def avglongwinstreak(n, pwin):
    return sum(fastpwinstreak(n, pwin, k) for k in range(1, n + 1))


def pwinstreak(n, pwin, k):
    memo = [0] * (n + 1)
    for m in range(k, n + 1):
        memo[m] = pwin**k + sum(pwin**j * (1 - pwin) * memo[m - (j + 1)]
                                for j in range(k))
    return memo[n]


def fastpwinstreak(n, pwin, k):
    pwink = pwin**k
    memo = [0] * (n + 1)
    windowsum = 0
    for m in range(k, n + 1):
        memo[m] = pwink + windowsum
        windowsum = pwin * windowsum + (1 - pwin) * (memo[m] - pwink *
                                                     memo[m - k])
    return memo[n]


print(avglongwinstreak(3000, 0.4))

Version that allows error:

def avglongwinstreak(n, pwin, abserr=0):
    avg = 0
    for k in range(1, n + 1):
        p = fastpwinstreak(n, pwin, k)
        avg += p
        if (n - k) * p < abserr:
            break
    return avg


def fastpwinstreak(n, pwin, k):
    pwink = pwin**k
    memo = [0] * (n + 1)
    windowsum = 0
    for m in range(k, n + 1):
        memo[m] = pwink + windowsum
        windowsum = pwin * windowsum + (1 - pwin) * (memo[m] - pwink *
                                                     memo[m - k])
    return memo[n]


print(avglongwinstreak(3000, 0.4, 1e-6))