数概率

【Puzzle】选票盒

字数统计: 1.6k字   |   阅读时长: 7分
,

摘要: 《概率50题》选票盒

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

这是概率面试题连载第 24 期，我们继续看 《Fifty challenging problems in probability with solutions》 中的一道题。

往期的内容整理在这篇文章里；或者看这个 github 仓库。

本题是关于概率DP的，这个点之前有一个类似的题目我们解决过，可以看这篇文章 【Puzzle】萨缪尔·佩皮斯的问题，我们来看一下这个题。

问题

The-Ballot-Box

In an election, two candidates, Albert and Benjamin, have in a ballot box a and b votes respectively with a>b
Ballots are randomly drawn and tallied
What is the chance that at least once after the first tally, the candidates have the same number of tallies?

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

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

思路参考

假设抽过几轮之后，现在还剩 i 张 A 选票，j 张 B 选票，也就是之前已经抽出了 a - i 张 A 选票，以及 b - j 张 B 选票。这种情况下我们记【从下一次抽取选票开始，直到盒子中的剩余票抽取完，双方已宣布票数至少出现一次平局】的概率为 dp[i][j]。在这个定义下，我们要求的就是 dp[a][b]。

此时已经抽出的 A 的票数 a - i，B 的票数 b - j，A 比 B 多出的差值是 k = a - i - (b - j)

我们考虑下一次抽取选票的两种情况，第一种是抽到 a，记为事件 Xa，概率为 Pa = i / (i + j)，第二种是抽到 b，记为事件 Xb，概率为 Pa = j / (i + j)。

由全概率公式，我们可以写出 dp[i][j] 的转移方程

1
dp[i][j] = pa * dp[i - 1][j] + pb * dp[i][j - 1]

dp[i][j] 的边界值

我们首先看一下 dp[i][j] 的一些边界值

  • (1) i = 0 且 j = 0

当 i, j 均为 0 时，盒子中已经没有票可以抽出，当然也就无法在后续的抽票中达成平局，dp[i][j] = 0

  • (2) i = 0

当 i 为 0，后续的抽票只会抽出 B 选票。因此 dp[i][j] 取决于已经抽出的 A 比 B 多出的部分 k 与还剩的 B 的张数 j

如果 k > 0 且 j >= k，则 dp[i][j] = 1，否则 dp[i][j] = 0

  • (3) j = 0

当 j 为 0，后续的抽票只会抽出 A 选票。因此 dp[i][j] 取决于已经抽出的 A 比 B 多出的部分 k 与还剩的 A 的张数 i

如果 k < 0 且 i >= -k，则 dp[i][j] = 1，否则 dp[i][j] = 0

dp[i][j] 的转移方程

下面我们分别考虑转移方程中的 ans1 = pa * dp[i - 1][j] 和 ans2 = pb * dp[i][j - 1] 这两部分结果。

  • (1) 下一次抽出的是 A 的选票

如果 k = -1，那么这次抽到 a 后就出现了双方平局的结果，因此 ans1 = pa * 1；

否则就要继续算 ans1 = pa * dp[i - 1][j]。

  • (2) 下一次抽出的是 B 的选票

如果 k = 1，那么这次抽到 b 后就出现了双方平局的结果，因此 ans2 = pb * 1；

否则就要继续算 ans2 = pb * dp[i][j - 1]。

概率 DP

有了以上推导之后，我们就可以用概率 DP 算法计算答案了，算法如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
状态定义
dp[i][j] := 还剩 i 张 A 选票，j 张 B 选票时，后续抽取过程中会出现平局的概率

答案
dp[a][b]

初始化
dp[0][0] = 0
dp[i][0] = 1    if i >= -k
         = 0    other
dp[0][j] = 1    if j >= k
         = 0    other

状态转移
ans1 = pa                   if k == -1
       pa * dp[i - 1][j]    other
ans2 = pb                   if k == 1
       pb * dp[i][j - 1]    other
dp[i][j] = ans1 + ans2

编程计算答案

Python 代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
@lru_cache(maxsize=int(1e6))
def solve(i, j):
    if i == 0 and j == 0:
        return 0.0
    k = (a - i) - (b - 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 * solve(i - 1, j)
    if k == 1:
        ans2 = pb
    else:
        ans2 = pb * solve(i, j - 1)
    return ans1 + ans2

以 a = 5, b = 8 为例，计算结果为 0.769231

蒙特卡洛模拟

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
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 test(args):
    a, b = args
    box = Box(a, b)
    ans = box.solve()
    np.random.seed()
    T = int(1e7)
    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

args = [(5, 5), (5, 6), (5, 7), (5, 8), (5, 9), (5, 10), (5, 11), (5, 12)]
pool = Pool(8)
ts = pool.map(test, args)
for t in ts:
    a, b, ans, sim = t
    print("a={}, b={}".format(a, b))
    print("    计算结果:{:.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
a=5, b=5
    计算结果:1.000000
    模拟结果:1.000000
a=5, b=6
    计算结果:0.909091
    模拟结果:0.908938
a=5, b=7
    计算结果:0.833333
    模拟结果:0.833231
a=5, b=8
    计算结果:0.769231
    模拟结果:0.769221
a=5, b=9
    计算结果:0.714286
    模拟结果:0.714182
a=5, b=10
    计算结果:0.666667
    模拟结果:0.666567
a=5, b=11
    计算结果:0.625000
    模拟结果:0.624997
a=5, b=12
    计算结果:0.588235
    模拟结果:0.588035
Share