【Puzzle】不公平的硬币3

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

352

2021-10-31

Puzzle28, 概率67

摘要: 不公平的硬币3

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

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

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

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

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

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

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

问题

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

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

思路参考

(1) 题目对比

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

为了解决这个新问题，我们在文章 Puzzle-不公平的硬币2 中先解决了一个简化版的问题，并总结了这类题的方法。为了方便对比，我们把简化版题目抄在下面。

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

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

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

可以看出，本题相比于我们已经解决的简化版问题，主要复杂在以下 3 点上

复杂点一：本题是 10 枚硬币，其中有 3 枚不公平硬币；简化版为 2 枚硬币，其中 1 枚不公平
复杂点二：本题不公平硬币抛出正面的概率是 0.8，简化版为 1.0
复杂点三：本题可以抛 3 次，简化版为抛 2 次

(2) 方法回顾

我们直接用解决简化版问题时总结的方法解决本题，下面我们把方法再梳理一下。

我们的目标是在各种可能的证据下，求各种可能的假设成立的条件概率，取其中条件概率最大的假设作为预测。进而计算整体的预测正确的概率。

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

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

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

$P(H_{j}|e_{i}) = \frac{P(e_{i}|H_{j})P(H_{j})}{P(e_{i})}$

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

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

$P(T|e_{i}) = \max\limits_{j}P(H_{j}|e_{i})$

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

$P(T) = \sum\limits_{i}P(e_{i})P(T|e_{i})$

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

$\begin{align} P(T) &= \sum\limits_{i}P(e_{i})P(T|e_{i}) \\ &= \sum\limits_{i}P(e_{i})\max\limits_{j}P(H_{j}|e_{i}) \\ &= \sum\limits_{i}P(e_{i})\max\limits_{j}\frac{P(e_{i}|H_{j})P(H_{j})}{P(e_{i})} \\ &= \sum\limits_{i}(\max\limits_{j}P(e_{i}|H_{j})P(H_{j})) \\ \end{align}$

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

$P(T) = \sum\limits_{i}(\max\limits_{j}P(e_{i}|H_{j})P(H_{j}))$

(3) 假设和证据

首先我们看都有哪些可能的证据。由于可以抛出 3 次(复杂点一)，因此观察到的证据可能有”正正正”、”正正反”、”正反正”、”正反反”, “反正正”、”反正反”、”反反正”、”反反反”这八种。分别记为 e1, e2, e3, e4, e5, e6, e7, e8。

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

下面我们基于这些可能的假设和证据，推导一下先验概率和似然值。

(4) 先验概率

由于本题是 10 枚硬币(复杂点二)，其中有 3 枚不公平硬币，因此

$P(H_{1}) = \frac{7}{10} \\ P(H_{2}) = \frac{3}{10} \\$

(5) 似然值

公平硬币抛出正面的概率为 p = 0.5，不公平硬币抛出正面的概率是 q = 0.8(复杂点三)。

如果硬币是公平硬币，则得到 8 种证据的概率均为 1/8，也就是 8 个似然值均为 1/8

$P(e_{i}|H_{1}) = p^{3} = \frac{1}{8}$

如果硬币是不公平硬币，则得到 8 种证据的概率分别为

$P(e_{1}|H_{2}) = q^{3} = \frac{64}{125}\\ P(e_{2}|H_{2}) = q^{2}(1 - q) = \frac{16}{125} \\ P(e_{3}|H_{2}) = q^{2}(1 - q) = \frac{16}{125} \\ P(e_{4}|H_{2}) = q(1 - q)^{2} = \frac{4}{125} \\ P(e_{5}|H_{2}) = q^{2}(1 - q) = \frac{16}{125} \\ P(e_{6}|H_{2}) = q(1 - q)^{2} = \frac{4}{125} \\ P(e_{7}|H_{2}) = q(1 - q)^{2} = \frac{4}{125} \\ P(e_{8}|H_{2}) = (1 - q)^{3} = \frac{1}{125} \\$

(6) 套公式

$\begin{align} P(T) &= \sum\limits_{i}(\max\limits_{j}P(e_{i}|H_{j})P(H_{j})) \\ &= \max(P(e_{1}|H_{1})P(H_{1}), P(e_{1}|H_{2})P(H_{2})) \\ & \quad\quad + \max(P(e_{2}|H_{1})P(H_{1}), P(e_{2}|H_{2})P(H_{2})) \\ & \quad\quad + \max(P(e_{3}|H_{1})P(H_{1}), P(e_{3}|H_{2})P(H_{2})) \\ & \quad\quad + \max(P(e_{4}|H_{1})P(H_{1}), P(e_{4}|H_{2})P(H_{2})) \\ & \quad\quad + \max(P(e_{5}|H_{1})P(H_{1}), P(e_{5}|H_{2})P(H_{2})) \\ & \quad\quad + \max(P(e_{6}|H_{1})P(H_{1}), P(e_{6}|H_{2})P(H_{2})) \\ & \quad\quad + \max(P(e_{7}|H_{1})P(H_{1}), P(e_{7}|H_{2})P(H_{2})) \\ & \quad\quad + \max(P(e_{8}|H_{1})P(H_{1}), P(e_{8}|H_{2})P(H_{2})) \\ &= \max(\frac{1}{8}\times\frac{7}{10}, \frac{64}{125}\times\frac{3}{10}) \\ & \quad\quad + \max(\frac{1}{8}\times\frac{7}{10}, \frac{16}{125}\times\frac{3}{10}) \\ & \quad\quad + \max(\frac{1}{8}\times\frac{7}{10}, \frac{16}{125}\times\frac{3}{10}) \\ & \quad\quad + \max(\frac{1}{8}\times\frac{7}{10}, \frac{4}{125}\times\frac{3}{10}) \\ & \quad\quad + \max(\frac{1}{8}\times\frac{7}{10}, \frac{16}{125}\times\frac{3}{10}) \\ & \quad\quad + \max(\frac{1}{8}\times\frac{7}{10}, \frac{4}{125}\times\frac{3}{10}) \\ & \quad\quad + \max(\frac{1}{8}\times\frac{7}{10}, \frac{4}{125}\times\frac{3}{10}) \\ & \quad\quad + \max(\frac{1}{8}\times\frac{7}{10}, \frac{1}{125}\times\frac{3}{10}) \\ &= \frac{64}{125}\times\frac{3}{10} + \frac{1}{8}\times\frac{7}{10} \times 7 \\ &= \frac{192}{1250} + \frac{49}{80} \\ &= \frac{7661}{10000} = 0.7661 \\ \end{align}$

可以看到，由于正常硬币的先验概率是 0.7，我们直接预测正常硬币可以得到 0.7 的正确概率。有了抛出三次的结果作为证据，我们可以得到 0.7661 的正确概率，有所提高，说明数据还是有用的。

蒙特卡洛模拟

模拟的核心是当得到证据 e 之后，如何给出预测结果，也就是代码中 model 实例中的 predict 方法。根据前面的推导，这取决于 16 个似然值。首先随机生成数据模拟出 8 个后验概率(代码中 model 实例的 train 方法)，再随机生成数据模拟预测正确的概率(代码中 model 实例的 test 方法)。

import numpy as np
from multiprocessing import Pool

class Coin:
    def __init__(self, p, t):
        self.p = p
        self.t = t # 硬币种类 id

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

class Model:
    def __init__(self, coins, n_flip):
        self.coins = coins
        self.n_flip = n_flip
        # mapping 记录每种证据下两种硬币种类的次数
        # 共有 16 个参数
        self.mapping = {}
        for i in range(1 << n_flip):
            e = [0] * n_flip
            for j in range(n_flip):
                e[n_flip - 1 - j] = i >> j & 1
            self.mapping[tuple(e)] = 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(1e5)
        for i in range(n_train):
            idx = np.random.randint(0, len(self.coins))
            coin = self.coins[idx]
            e = []
            for _ in range(self.n_flip):
                e.append(int(coin.flip()))
            e = tuple(e)
            self.update(e, coin.t)
        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, len(self.coins))
            coin = self.coins[idx]
            # 抛出三次，获取证据
            e = []
            for _ in range(self.n_flip):
                e.append(int(coin.flip()))
            e = tuple(e)
            # 根据证据猜硬币种类
            if self.predict(e) == coin.t:
                correct += 1
        return correct / T

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

模拟结果

证据 (0, 0, 0) 下的正常硬币有8775个, 特殊硬币有239个
证据 (0, 0, 1) 下的正常硬币有8857个, 特殊硬币有988个
证据 (0, 1, 0) 下的正常硬币有8803个, 特殊硬币有933个
证据 (0, 1, 1) 下的正常硬币有8708个, 特殊硬币有3729个
证据 (1, 0, 0) 下的正常硬币有8686个, 特殊硬币有953个
证据 (1, 0, 1) 下的正常硬币有8764个, 特殊硬币有3858个
证据 (1, 1, 0) 下的正常硬币有8685个, 特殊硬币有3774个
证据 (1, 1, 1) 下的正常硬币有8913个, 特殊硬币有15335个
===训练结束===
猜正确的概率是: 0.765926
猜正确的概率是: 0.765982
猜正确的概率是: 0.766050
猜正确的概率是: 0.766293
猜正确的概率是: 0.766128
猜正确的概率是: 0.765904
猜正确的概率是: 0.766045
猜正确的概率是: 0.766111

潮汐朝夕的生活实验室 \Doge 陪伴一个算法工程师的职业生涯

数概率