数概率

【Puzzle】不公平的硬币2

字数统计: 3.7k字   |   阅读时长: 15分
,

摘要: 不公平的硬币2

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

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

目前的题目来源主要是《Fifty challenging problems in probability with solutions》这本书以及朋友在面试中遇到的题。对于每道题，首先给出思路参考，然后用蒙特卡洛模拟的方式验证结果。

在上一期，我们看了一个贝叶斯公式里面比较简单但是很有代表性的一个题。今天我们在那道题的基础上做一点改变形成一个新的题目来看看。为了方便对比，我们首先把上一期的题目抄一遍如下

有 10 枚硬币，其中有 1 枚是不公平的，随机抛出正面的概率为 0.8，另外 9 枚都是公平的，也就是随机抛出正面的概率是 0.5。

现在从 10 枚硬币中随机取出 1 枚，随机地抛了 3 次，结果为”正正反”。求该硬币为不公平硬币的概率。

如果我们把问题的问法改一下，题目的其它条件不变，就可以变成一个新的题目，如下

有 10 枚硬币，其中有 1 枚是不公平的，随机抛出正面的概率为 0.8，另外 9 枚都是公平的，也就是随机抛出正面的概率是 0.5。

现在从 10 枚硬币中随机取出 1 枚，可以随机地抛 3 次，根据结果来预测该硬币是否为不公平硬币。求预测正确的概率。

原问题是给定一个特定的证据下，求特定的假设成立的条件概率。修改后的新问题是在各种可能的证据下，求各种可能的假设成立的条件概率，进而取其中条件概率最大的假设作为预测。

这个问题还是挺复杂的，本期我们暂时不解决这个问题，而是先看一个简化后的问题找找感觉。下一期我们再基于简化问题的思路解决上面这个我们最终要解决的问题。

问题

有 2 枚硬币，一枚是普通硬币，也就是随机抛出正面的概率是 0.5，另一枚是特殊硬币，两面都是正面，也就是随机抛出正面的概率是 1。

现在随机取出一枚硬币，允许抛出 2 次并根据两次抛出结果来预测该硬币是普通硬币还是特殊硬币。

问预测正确的概率是多少。

思路参考

本题与最终要解决的问题整体上是一样的，只是做了以下三点简化:

  • 原问题是 10 枚硬币，其中有 1 枚不公平硬币；本题简化为 2 枚，其中 1 枚不公平
  • 原问题不公平硬币抛出正面的概率是 0.8，本题简化为 1.0
  • 原问题可以抛 3 次，本题简化为抛 2 次
我们的目标是在各种可能的证据下，求各种可能的假设成立的条件概率，取其中条件概率最大的假设作为预测。进而计算整体的预测正确的概率。

具体地分为 3 个步骤:

  • step1: 用贝叶斯公式求出在各种可能的证据下，各种可能的假设成立的条件概率
  • step2: 对每一种可能的证据，取各种可能的假设中条件概率最大的假设作为预测，该条件概率就是拿到该证据时预测正确的概率
  • step3: 用全概率公式计算整体的预测正确的概率

下面我们分别看这三个步骤

step1:

首先我们看都有哪些可能的证据。由于可以抛出 2 次，因此观察到的证据可能有”正正”、”正反”、”反正”、”反反”这四种。分别记为 e1, e2, e3, e4。

然后我们看都有哪些可能的假设。由于有公平和不公平这两种硬币，因此针对随机取出的硬币有两种假设，”公平” 和 “不公平”，分别记为 H1, H2。

我们目标所求的在各种可能的证据下，求各种可能的假设成立的条件概率，就是 P(H1|e1), P(H1|e2), P(H1|e3), P(H1|e4), P(H2|e1), P(H2|e2), P(H2|e3), P(H2|e4) 这八个条件概率。

由于 H1 和 H2 是对立事件，因此 P(H2|e) = 1 - P(H1|e)，我们实际上只需要求 4 个条件概率，P(H1|e1), P(H1|e2), P(H1|e3), P(H1|e4)。

我们在上一期解决的问题是给定一个特定的证据下，求特定的假设成立的条件概率。相当于求这里的 4 个条件概率中的 1 个。我们用同样的方法分别求出 4 个条件概率即可。

下面我们以 P(H1|e1) 为例，首先写出用贝叶斯公式计算 P(H1|e1) 的公式

由于硬币是 2 枚中有 1 枚不公平，因此随机取出 1 枚的两种假设的先验概率分别为 P(H1) = P(H2) = 0.5。

下面我们看两种假设下，抛两次产生 e1 这个证据的似然值

  • 如果取出的是正常硬币(H1)，抛两次产生 e1 这个证据的概率是 P(e1|H1) = 1/4
  • 如果取出的是特殊硬币(H2)，抛两次产生 e1 这个证据的概率是 P(e1|H2) = 1

有了以上先验概率和似然值，我们就可以用贝叶斯公式计算条件概率了。

P(H1|e2), P(H1|e3), P(H1|e4) 也可以类似地算出，首先我们把所需的似然值写出来。

  • 如果取出的是正常硬币(H1)，抛两次产生 e2 这个证据的概率是 P(e2|H1) = 1/4
  • 如果取出的是特殊硬币(H2)，抛两次产生 e2 这个证据的概率是 P(e2|H2) = 0
  • 如果取出的是正常硬币(H1)，抛两次产生 e3 这个证据的概率是 P(e3|H1) = 1/4
  • 如果取出的是特殊硬币(H2)，抛两次产生 e3 这个证据的概率是 P(e3|H2) = 0
  • 如果取出的是正常硬币(H1)，抛两次产生 e4 这个证据的概率是 P(e4|H1) = 1/4
  • 如果取出的是特殊硬币(H2)，抛两次产生 e4 这个证据的概率是 P(e4|H2) = 0

step2:

对证据 e1, e2, e3, e4，分别取各种可能的假设中条件概率最大的假设作为预测，该条件概率就是拿到该证据时预测正确的概率

step3:

用全概率公式计算整体的预测正确的概率

总结

这里我们对这类问题做个总结，首先我们的目标是在各种可能的证据下，求各种可能的假设成立的条件概率，取其中条件概率最大的假设作为预测。进而计算整体的预测正确的概率。

基于这个大方向，具体做的时候分为 3 个步骤:

  • step1: 用贝叶斯公式求出在各种可能的证据下，各种可能的假设成立的条件概率

具体地，在证据 ei 下，假设 Hj 成立的条件概率为

  • step2: 对每一种可能的证据，取各种可能的假设中条件概率最大的假设作为预测，该条件概率就是拿到该证据时预测正确的概率

具体地，在证据 ei 下，预测正确的概率为

  • step3: 用全概率公式计算整体的预测正确的概率

逐层代入之后，会发现 P(T) 的最终公式中 P(ei) 是被削掉了，如下

其中 P(Hj) 是先验概率，也就是说我们其实只需要求出所有似然值 P(ei|Hj) 即可套下面的最终公式了，而无需做一步用贝叶斯公式求后验概率了，因此证据因子 P(ei) 不用算了。

前面为了思路的连贯性，还是按步骤推的答案。如果套最终公式，计算过程就是下面这样

蒙特卡洛模拟

模拟的核心是当得到证据 e 之后，如何给出预测结果，也就是代码中 model 实例中的 predict 方法。根据前面的推导，这取决于 8 个后验概率，这 8 个后验概率可以直接根据计算结果手动写进去，也可以用数据模拟出这 8 个后验概率，类似于用训练数据学出模型的 8 个参数，然后再用模型和另一批数据模拟预测正确的概率。

代码中实现的是后者，也就是先随机生成数据模拟出 8 个后验概率(代码中 model 实例的 train 方法)，再随机生成数据模拟预测正确的概率(代码中 model 实例的 test 方法)。

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
import numpy as np
from multiprocessing import Pool

class Coin:
    def __init__(self, p):
        self.p = p

    def flip(self):
        """
        如果抛出正面，返回 True
        """
        return np.random.rand() < self.p

class Model:
    def __init__(self, coins):
        self.coins = coins
        # mapping 记录每种证据下两种硬币种类的次数
        # 共有 8 个参数
        self.mapping = {(0, 0): np.array([0, 0])
                       ,(0, 1): np.array([0, 0])
                       ,(1, 0): np.array([0, 0])
                       ,(1, 1): np.array([0, 0])
                       }

    def update(self, e, idx):
        """
        根据数据 (e, idx) 更新 mapping 的参数
        数据的含义是证据 e 对应硬币种类 idx
        """
        self.mapping[e][idx] += 1

    def train(self):
        """
        随机生成数据训练
        """
        n_train = int(1e7)
        for i in range(n_train):
            idx = np.random.randint(0, 2)
            coin = self.coins[idx]
            r1 = coin.flip()
            r2 = coin.flip()
            e = (int(r1), int(r2))
            self.update(e, idx)
            if (i + 1) % int(2e6) == 0:
                print("训练进度: {}/{}".format((i + 1), n_train))
                for k, v in self.mapping.items():
                    print("    证据 {} 下的正常硬币有{}个, 特殊硬币有{}个".format(k, v[0], v[1]))
        print("===训练结束===")

    def predict(self, e):
        """
        根据证据 e 猜硬币种类，返回 0 或 1
        """
        return np.argmax(self.mapping[e])

    def test(self, T):
        correct = 0
        np.random.seed()
        for _ in range(T):
            # 随机选硬币
            idx = np.random.randint(0, 2)
            coin = coins[idx]
            # 抛出两次，获取证据
            r1 = coin.flip()
            r2 = coin.flip()
            e = (int(r1), int(r2))
            # 根据证据猜硬币种类
            if self.predict(e) == idx:
                correct += 1
        return correct / T

coins = [Coin(0.5), Coin(1.0)]
model = Model(coins)
model.train()
T = int(1e6)
args = [T] * 8
pool = Pool(8)
ts = pool.map(model.test, args)
for t in ts:
    print("猜正确的概率是: {:.6f}".format(t))

模拟结果

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
训练进度: 2000000/10000000
    证据 (0, 0) 下的正常硬币有250160个, 特殊硬币有0个
    证据 (0, 1) 下的正常硬币有250105个, 特殊硬币有0个
    证据 (1, 0) 下的正常硬币有249805个, 特殊硬币有0个
    证据 (1, 1) 下的正常硬币有251115个, 特殊硬币有998815个
训练进度: 4000000/10000000
    证据 (0, 0) 下的正常硬币有499945个, 特殊硬币有0个
    证据 (0, 1) 下的正常硬币有499713个, 特殊硬币有0个
    证据 (1, 0) 下的正常硬币有499823个, 特殊硬币有0个
    证据 (1, 1) 下的正常硬币有501512个, 特殊硬币有1999007个
训练进度: 6000000/10000000
    证据 (0, 0) 下的正常硬币有749591个, 特殊硬币有0个
    证据 (0, 1) 下的正常硬币有750072个, 特殊硬币有0个
    证据 (1, 0) 下的正常硬币有750249个, 特殊硬币有0个
    证据 (1, 1) 下的正常硬币有750605个, 特殊硬币有2999483个
训练进度: 8000000/10000000
    证据 (0, 0) 下的正常硬币有999616个, 特殊硬币有0个
    证据 (0, 1) 下的正常硬币有1000364个, 特殊硬币有0个
    证据 (1, 0) 下的正常硬币有1000514个, 特殊硬币有0个
    证据 (1, 1) 下的正常硬币有1000868个, 特殊硬币有3998638个
训练进度: 10000000/10000000
    证据 (0, 0) 下的正常硬币有1249086个, 特殊硬币有0个
    证据 (0, 1) 下的正常硬币有1250621个, 特殊硬币有0个
    证据 (1, 0) 下的正常硬币有1251428个, 特殊硬币有0个
    证据 (1, 1) 下的正常硬币有1249869个, 特殊硬币有4998996个
===训练结束===
猜正确的概率是: 0.874737
猜正确的概率是: 0.875113
猜正确的概率是: 0.874732
猜正确的概率是: 0.874939
猜正确的概率是: 0.874946
猜正确的概率是: 0.874690
猜正确的概率是: 0.875392
猜正确的概率是: 0.874712
Share