数概率

概率DP的解析解

字数统计: 2.7k字   |   阅读时长: 12分

摘要: 选票盒问题和系列赛问题的解析解

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

这是概率面试题连载第 28 期，往期的内容整理在这篇文章里；或者看这个 github 仓库。

今天我们回炉两道之前解决过的问题，这两道题都是概率DP的问题，当时我们通过分析问题找到概率DP的算法并编程求解的，本文我们推一下解析解，这两个题是递进的关系，第二题用到了第一题的结论，下面是这两个题的文章链接。

下面我们分别看这两个题

问题1: 选票盒

一场选举两个候选人 A 和 B，选票盒里有 A, B 两个候选人的选票，其中 A 有 a 张票，B 有 b 张票，a > b。选票被随机抽出并宣布。

问：从宣布第一张票开始直至盒子中的票宣布完，双方已宣布票数至少出现一次平局的概率

思路参考

首先我们从简单的情况开始看，假设 a = 3, b = 2，宣布选票的等可能序列如下，共十种

  • AAABB
  • AABAB
  • ABAAB
  • BAAAB
  • AABBA
  • ABABA
  • BAABA
  • ABBAA
  • BABAA
  • BBAAA

标红的是形成平局的序列。因此 a = 3, b = 2 时，至少出现一次平局的概率为 8/10

下面我们考虑一般情况，设当第2n张票宣布的时候出现第一次平局，由于 a > b，因此 n <= b。

对于每一个满足要求的序列，有两种可能情况

  • 第1类序列: A 的票数总是领先，直至第 2n 张票
  • 第2类序列: B 的票数总是领先，直至第 2n 张票

这两类序列的个数是一样多的，因为对于第1类的任意一个序列，将 A 改成 B，将 B 改成 A，就得到第 2 类的对应的序列。

由于这两类序列的个数是一样多，这两类序列出现的概率也就是一样的，因此我们不妨研究 B 一直领先直至第 2n 张票的情况。

B 想要一致领先，第一张票就一定是 B，而当 B 是第一张票的时候，无论后续宣布的票是什么顺序，在耗尽之前一定会形成平局，因为 a > b。

因此【B 一直领先直至第 2n 张票】的概率，就等于是第一张票宣布到 B 的概率。而这个概率由古典概型就可以得到 b / (a + b)

那么两类序列加起来总的概率就是 2b / (a + b)

验证结果

下面的代码对比了 b / (a + b) 的结果，以及概率DP的计算结果和蒙特卡洛模拟的结果。

其中概率DP的推导过程可以看 【Puzzle】选票盒

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
import numpy as np
from functools import lru_cache
from multiprocessing import Pool
import sys
sys.setrecursionlimit(int(1e5))

class Box:
    def __init__(self, a, b):
        self.na = a
        self.nb = b
        self.ma = self.mb = 0

    def solve(self):
        self.reset()
        return self._solve(self.na, self.nb)

    @lru_cache(maxsize=int(1e6))
    def _solve(self, i, j):
        if i == 0 and j == 0:
            return 0.0
        k = (self.na - i) - (self.nb - j)
        if i == 0:
            if k > 0 and j >= k:
                return 1.0
            return 0.0
        if j == 0:
            if k < 0 and i >= -k:
                return 1.0
            return 0.0
        pa = i / (i + j)
        pb = 1.0 - pa
        if k == -1:
            ans1 = pa
        else:
            ans1 = pa * self._solve(i - 1, j)
        if k == 1:
            ans2 = pb
        else:
            ans2 = pb * self._solve(i, j - 1)
        return ans1 + ans2

    def reset(self):
        self.na += self.ma
        self.nb += self.mb
        self.ma = self.mb = 0

    def draw(self):
        if np.random.randint(0, self.na + self.nb) < self.na:
            self.na -= 1
            self.ma += 1
        else:
            self.nb -= 1
            self.mb += 1
        return self.ma == self.mb

    def empty(self):
        return self.na + self.nb == 0

    def calc(self):
        return 2 * self.nb / (self.na + self.nb)

def test(args):
    a, b = args
    box = Box(a, b)
    ans = box.solve()
    x = box.calc()
    np.random.seed()
    T = int(1e5)
    n_tie = 0
    for _ in range(T):
        while not box.empty():
            if box.draw():
                n_tie += 1
                break
        box.reset()
    return a, b, ans, n_tie / T, x

args = [(5, 5), (6, 5), (7, 5), (8, 5), (9, 5), (10, 5), (11, 5), (12, 5)]
pool = Pool(8)
ts = pool.map(test, args)
for t in ts:
    a, b, ans, sim, x = t
    print("a={}, b={}".format(a, b))
    print("    概率DP计算结果:{:.6f}".format(ans))
    print("    模拟结果:{:.6f}".format(sim))
    print("    解析解:{:.6f}".format(x))

模拟结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
a=5, b=5
    概率DP计算结果:1.000000
    模拟结果:1.000000
    解析解:1.000000
a=6, b=5
    概率DP计算结果:0.909091
    模拟结果:0.910470
    解析解:0.909091
a=7, b=5
    概率DP计算结果:0.833333
    模拟结果:0.835590
    解析解:0.833333
a=8, b=5
    概率DP计算结果:0.769231
    模拟结果:0.768290
    解析解:0.769231
a=9, b=5
    概率DP计算结果:0.714286
    模拟结果:0.712870
    解析解:0.714286
a=10, b=5
    概率DP计算结果:0.666667
    模拟结果:0.667780
    解析解:0.666667
a=11, b=5
    概率DP计算结果:0.625000
    模拟结果:0.626520
    解析解:0.625000
a=12, b=5
    概率DP计算结果:0.588235
    模拟结果:0.589470
    解析解:0.588235

问题2: 系列赛中不出现平局

A 和 B 进行一场 N 轮的系列赛，并记录每一轮胜负，双方在一轮中的获胜概率均为 0.5

问：第1场比赛之后，直到第 N 场比赛比完之后，双方始终没有出现平局的概率

思路参考

一共比赛 N 轮，胜负序列共有 $2^{N}$ 种。

首先我们讨论一下 N 的奇偶问题。

如果在胜负序列中出现平局，则平局一定出现在第 2, 4, …, 2n 局中。

因此如果 N 为奇数的时候，包含平局的胜负序列取决于前 N-1 局的胜负序列，如果前 N-1 局的胜负序列没有出现平局，则第 N 局结束时也肯定不会出现平局。所以当 N 为奇数的时候，出现平局的概率就是 N-1 局比赛中出现平局的概率。

下面我们只需要考虑 N 为偶数时，胜负序列含平局的概率，记 N = 2n。

我们假设其中 A 赢了 a 场，B 赢了 b = N - a 场。这种情况出现的概率记为 P(a)，可以由古典概型给出

现在我们要问：A 赢了 a 场，B 赢了 b 场的情况下，胜负序列会出现平局的概率有多少，这个概率记为 P(T|a)。

我们可以直接用刚刚的选票盒那个题的结论：

  • 如果 a <= b，即 a <= n，概率是 2a / (a + b) = 2a/N
  • 如果 a > b，即 a > n，概率是 2b / (a + b) = 2(N-a)/N

由全概率公式，整体会出现平局的概率如下

将组合数用阶乘展开，然后化简

由上面的式子的形式，我们可以联想到二项式定理(注意中间少了一项)

我们可以得到当 N 为偶数 (N=2n) 时，胜负序列含有平局的概率 P(T) 的最终结果。

最终答案：不含平局的概率如下(下面的化简基于 N=2n)

前面我们论证过当 N 为奇数时，结果就是上面公式中 N 取 N - 1 即可。

验证结果

下面的代码对比了 的结果，以及概率DP的计算结果和蒙特卡洛模拟的结果。

其中概率DP的推导过程可以看 【Puzzle】系列赛中不出现平局

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
import numpy as np
from functools import lru_cache
from multiprocessing import Pool
import sys
sys.setrecursionlimit(int(1e5))

class Series:
    def __init__(self, N):
        self.N = N
        self.ma = self.mb = 0

    def solve(self):
        self.reset()
        return self._solve(0, 0)

    @lru_cache(maxsize=int(1e6))
    def _solve(self, i, j):
        if self.N - i == 0:
            return 1.0
        k = j * 2 - i
        if abs(k) > self.N - i:
            return 1.0
        pa = pb = 0.5
        if k == -1:
            ans1 = 0.0
        else:
            ans1 = pa * self._solve(i + 1, j + 1)
        if k == 1:
            ans2 = 0.0
        else:
            ans2 = pb * self._solve(i + 1, j)
        return ans1 + ans2

    def reset(self):
        self.N += self.ma + self.mb
        self.ma = self.mb = 0

    def match(self):
        if np.random.rand() < 0.5:
            self.ma += 1
        else:
            self.mb += 1
        self.N -= 1
        return self.ma == self.mb

    def finish(self):
        return self.N == 0

    def calc(self):
        from scipy.special import comb
        if self.N % 2 == 0:
            return comb(self.N, self.N / 2) / (2 ** self.N)
        else:
            return comb(self.N - 1, (self.N - 1) / 2) / (2 ** (self.N - 1))

def test(N):
    series = Series(N)
    ans = series.solve()
    res = series.calc()
    np.random.seed()
    T = int(1e6)
    n_tie = 0
    for _ in range(T):
        while not series.finish():
            if series.match():
                n_tie += 1
                break
        series.reset()
    return N, ans, res, (T - n_tie) / T

args = [6, 7, 8, 9, 12, 13, 18, 19]
pool = Pool(8)
ts = pool.map(test, args)
for t in ts:
    N, ans, res, sim = t
    print("N={}".format(N))
    print("    解析解:{:.6f}".format(res))
    print("    概率DP计算结果:{:.6f}".format(ans))
    print("    模拟结果:{:.6f}".format(sim))

模拟结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
N=6
    解析解:0.312500
    概率DP计算结果:0.312500
    模拟结果:0.313608
N=7
    解析解:0.312500
    概率DP计算结果:0.312500
    模拟结果:0.312354
N=8
    解析解:0.273438
    概率DP计算结果:0.273438
    模拟结果:0.272503
N=9
    解析解:0.273438
    概率DP计算结果:0.273438
    模拟结果:0.273682
N=12
    解析解:0.225586
    概率DP计算结果:0.225586
    模拟结果:0.225618
N=13
    解析解:0.225586
    概率DP计算结果:0.225586
    模拟结果:0.225854
N=18
    解析解:0.185471
    概率DP计算结果:0.185471
    模拟结果:0.185324
N=19
    解析解:0.185471
    概率DP计算结果:0.185471
    模拟结果:0.185529
Share