如何获得大量执行伯努利实验的详细结果（概率树）

Question

假设如下实验：执行相同的伯努利试验（成功概率为P）N次

我需要以下信息：所有可能的成功/失败序列及其发生概率。

示例： 成功概率为 P = 40% 的 Bernouilli 实验执行 3 次将产生以下结果（S 为成功，F 为失败）：

FFF 0.216

SFF 0.144

FSF 0.144

SSF 0.096

FFS 0.144

SFS 0.096

FSS 0.096

SSS 0.064

我试图暴力破解它以获得结果，但它仅在 N = 25 时迅速窒息，我得到了 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);
        }
    }
}

我需要能够及时支持最多 N = 3000（任意 P 不到 3 秒即可获得结果）

有没有一种数学方法可以最佳地做到这一点？

Answer 1

由于您对最长连胜的长度感兴趣，因此在试验的任何时候，历史上只有两个相关事实：(1) 迄今为止最长的连胜有多长时间 (m) 和 ( 2）当前的连胜纪录是多久（k）。初始状态为 (0, 0)。转换是 (m, k) -> (max(m, k + 1), k + 1) 获胜， (m, k) -> (m, 0) 失败。一旦我们知道所有最终状态的概率，我们就可以进行平均。

编辑：此代码的修订版本进行了优化，以减少极不可能发生的事件所需的计算。具体来说，我们忽略长度大于或等于某个cutoff 的所有状态。给定绝对误差预算abserr，我们确定我们可以从最终的期望值计算中排除高达abserr / n 的概率质量，因为最长的连续可能是n。开始连胜的地方少于n，并且在每个地方，出现长度为cutoff 的连胜的概率为pwin**cutoff。使用粗略的联合界限，我们需要

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

最后一个不等式反转方向，因为log(pwin) < 0。

尽管评估器质量不佳（即 2014 年硬件上的 Python 3 解释器），但修改后的代码运行时间不到三秒。

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

这是一种不同的方法，它准确且足够快只是二次方。最长连胜的期望值等于

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

我们对概率的推理如下。要么以长度-k 连续获胜（概率pwin**k）打开记录，要么以j 获胜（概率为j in 0..k-1）然后输掉（概率pwin**j * (1 - pwin)），在这种情况下概率为等于在n - (j + 1) 尝试中连续获胜的概率-k。我们使用 memoization 来评估这个逻辑在pwinstreak 中暗示的重复； fastpwinstreak 中的更快版本使用代数来避免重复求和。

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))

允许出错的版本：

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))