4.概率面试题（连载）

4.17圆周上取 3 个点构成锐角三角形

在一个圆环上随机取 3 点，求这 3 个点组成一个锐角三角形的概率。

思路参考

如图，在圆周上取不同的两个点 X, Y，当 X, Y 确定后，再在圆周上取点 Z，要 ΔXYZ 构成锐角三角形，需要 X 在圆弧 AB 上，如果超出圆弧 AB 的范围，有以下几种情况，哪一种都会使得 ΔXYZ 不成为锐角三角形。

Z 如果在圆弧 AX 上，则角 X 就会变成钝角
Z 如果在圆弧 YB 上，则角 Y 就会变成钝角
Z 如果在圆弧 XY 上，则角 Z 就会变成钝角

当弧 XY 确定后，对应的圆心角为 θ，该圆心角的范围是 [0, pi)，其概率密度函数为

p(θ)=π/1

此时再取点 Z 落在弧 AB 上的概率为 AB 对应的圆心角，也就是 θ 与整个圆的角度范围的比值

P(A∣θ)= θ/2π

由连续型随机变量的全概率公式

蒙特卡洛模拟

import numpy as np
import math
from multiprocessing import Pool

def dist(x, y):
    return (x[0] - y[0]) ** 2 + (x[1] - y[1]) ** 2

def test(T):
    n = 0
    np.random.seed()
    for _ in range(T):
        theta_x = np.random.uniform(0.0, 2 * math.pi)
        theta_y = np.random.uniform(0.0, 2 * math.pi)
        theta_z = np.random.uniform(0.0, 2 * math.pi)
        x = (math.cos(theta_x), math.sin(theta_x))
        y = (math.cos(theta_y), math.sin(theta_y))
        z = (math.cos(theta_z), math.sin(theta_z))
        d_xy = dist(x, y)
        d_yz = dist(y, z)
        d_zx = dist(z, x)
        if not d_xy + d_yz > d_zx:
            continue
        if not d_yz + d_zx > d_xy:
            continue
        if not d_zx + d_xy > d_yz:
            continue
        n += 1
    return n / T

T = int(1e6)
args = [T] * 8
pool = Pool(8)
ts = pool.map(test, args)
for t in ts:
    print("锐角三角形的概率: {:.6f}".format(t))

模拟结果

锐角三角形的概率: 0.249651
锐角三角形的概率: 0.250418
锐角三角形的概率: 0.249653
锐角三角形的概率: 0.249947
锐角三角形的概率: 0.250386
锐角三角形的概率: 0.250427
锐角三角形的概率: 0.250429
锐角三角形的概率: 0.250860

4.18选票盒

一场选举两个候选人 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] 的转移方程

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

dp[i][j] 的边界值

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

i = 0 且 j = 0

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

i = 0, j != 0

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

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

j = 0, i != 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] 这两部分结果。

下一次抽出的是 A 的选票

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

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

下一次抽出的是 B 的选票

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

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

概率 DP

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

状态定义
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 代码如下

@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

蒙特卡洛模拟

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

模拟结果

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

4.19 向正方形区域扔硬币

在一个游戏中，一个玩家从 5 英尺远的地方把一分钱硬币扔到一张 1 12英寸正方形的桌子上，硬币直径为 3/4 英寸。

如果硬币完全落在正方形内，则玩家赢。

如果硬币落在桌上，他赢的机会有多大？

思路参考

如果硬币可以落在桌子上，那么硬币的重心在正方形范围内。

如果硬币完全落在正方形内，硬币重心到正方形四条边的距离都要大于硬币的半径。

硬币半径为 3/8 英寸。设硬币重心为 (x, y)，那么要满足以上条件，需要

求出 x，y 的范围均是 [3/8, 5/8]。根据几何概型概率公式，当硬币落到桌子上时，赢的机会有 (5/8 - 3/8) * (5/8 - 3/8) / 1 = 1/16 = 0.0625。

蒙特卡洛模拟

import numpy as np

def check(x, y):
    return x >= 3/8 and (1 - x) >= 3/8 and y >= 3/8 and (1 - y) >= 3/8

def test():
    T = int(1e7)
    s = 0
    for t in range(T):
        x = np.random.rand()
        y = np.random.rand()
        if check(x, y):
            s += 1
    print("his chance to win: {:.6f}".format(s / T))

for i in range(10):
    test()

模拟结果

his chance to win: 0.062507
his chance to win: 0.062518
his chance to win: 0.062433
his chance to win: 0.062596
his chance to win: 0.062261
his chance to win: 0.062555
his chance to win: 0.062552
his chance to win: 0.062652
his chance to win: 0.062406
his chance to win: 0.062344

4.20 不公平的硬币 2

在第三章贝叶斯公式例题中，我们解决了初级的「不公平的硬币」问题。

不公平的硬币 2，不公平的硬币 3，不公平的硬币 4 均为该题的变种，考察点也是贝叶斯公式，并且解题过程也类似。

下面是不公平的硬币 2 的描述。

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

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

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

思路参考

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

原问题是 10 枚硬币，其中有 1 枚不公平硬币；本题简化为 2 枚，其中 1 枚不公平；
原问题不公平硬币抛出正面的概率是 0.8，本题简化为 1.0；
原问题可以抛 3 次，本题简化为抛 2 次。

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

具体地分为 3 个步骤:

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

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

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 个步骤:

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

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

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

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

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

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

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

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

蒙特卡洛模拟

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

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

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

模拟结果

训练进度: 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

4.21不公平的硬币 3

在第三章的贝叶斯公式贝叶斯公式的例题中，我们解决了基础版的不公平的硬币问题。

不公平的硬币 2，不公平的硬币 3，不公平的硬币 4 均为该题的变种，考察点也是贝叶斯公式，并且解题过程也类似。

下面是不公平的硬币 3 的描述。

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

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

思路参考

(1) 题目对比

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

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

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

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

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

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

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

(2) 方法回顾

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

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

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

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

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

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

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

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

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

(3) 假设和证据

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

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

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

(4) 先验概率

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

(5) 似然值

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

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

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

(6) 套公式

记 P(T) 为猜测正确的概率，P(T|e) 为当获取到证据为 e 时，猜测正确的概率。

可以看到，由于正常硬币的先验概率是 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

4.22不公平的硬币 4

在第三章的贝叶斯公式的例题中，我们解决了基础版的不公平的硬币问题。

不公平的硬币 2，不公平的硬币 3，不公平的硬币 4 均为该题的变种，考察点也是贝叶斯公式，并且解题过程也类似。

下面是不公平的硬币 4 的描述。

我们在以上不公平的硬币 3 的变种的基础上，再换一种问法，本题由于比较复杂，需要用有向图来对整体流程进行建模。下面我们来看这个题。

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

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

目前已经抛出了第 1 次。

(1) 如果第一次抛出的结果是正面，则三次抛出完成后预测正确的概率是多少？

(2) 如果第一次抛出的结果是反面，则三次抛出完成后预测正确的概率是多少？

(1) 问题对比

本题与文章不公平的硬币 3 解决的问题相比，仅仅是增加了一个已知第一次抛出的结果，求的是已知第一次抛出结果后，预测正确的条件概率。为了方便对比，我们首先把上一期的题目抄一遍如下

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

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

该题求的是未知第一次抛出结果时的预测正确的概率，答案是 0.7661，计算过程可以参考上一期文章。

对于本题，我们还是先把假设和证据分别是什么列清楚，然后把各个假设的先验概率以及各个假设下得到各个证据的似然值梳理清楚，再分析已知第一次抛出的结果后，预测正确的条件概率的计算。

(2) 假设和证据

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

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

(3) 先验概率

由于本题是 10 枚硬币，其中有 3 枚不公平硬币，因此

(4) 似然值

公平硬币抛出正面的概率为 p = 0.5，不公平硬币抛出正面的概率是 q = 0.8。

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

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

(5) 第一次抛出后的情况

第一次抛出可能有两种结果：正面和反面，分别记为 X1 和 X2。如果是正面，最终预测对的概率就是 P(T|X1)，如果是反面，最终预测对的概率就是 P(T|X2)。

P(T|X1) 和 P(T|X2) 正是我们要求的值。这两个值在整个获取证据流程中是中间的位置，上游有整体的正确概率 P(T)；下游有三次抛出完成后，各个证据下的正确概率 P(T|e)，还是有点复杂的，我们借助有向图来理解会比较直观。

(6) 有向图建模

X1 对应 e1, e2, e3, e4 这四种证据，因此问当第一次抛出结果为 X1 时预测正确的概率，相当于问 3 次抛出结果后取得的证据为 e1, e2, e3, e4 这四种之一时预测正确的概率。

X2 对应 e5, e6, e7, e8 这四种证据，因此问当第一次抛出结果为 X2 时预测正确的概率，相当于问 3 次抛出结果后取得的证据为 e5, e6, e7, e8 这四种之一时预测正确的概率。

为了理解以上的公式，我们把整个流程抽象成一个有向图，如下：

三层节点的含义如下

根节点（第一层节点）表示一次还没有抛出，还可以抛出 3 次的状态，该状态的值表示整体的预测正确的概率。
第二层节点表示还剩下两次没有抛出的状态，根据第一次的结果可以分为两种状态，状态的值表示在该状态下预测正确的概率。
第三层节点表示三次抛出均完成的状态，根据三次抛出的结果可以分为八种状态（8 种可能证据），状态的值表示在该证据下下预测正确的概率。

有向边表示从一个状态跳向另一个状态的概率。要计算某节点表示的状态对应的预测正确概率，只需找到该节点的所有出边对应的下一个节点持有的概率值，与边权表示的状态转移概率进行加权求和后就可以得到

例如在上图中，如果我们要求 P(T)，只需要找到两个出边对应的下游节点 P(T|X1) 和 P(T|X2)，以及转移概率 P(X1) 和 P(X2)，即可得到

如果要求 P(T|X1)，只需要找到四个出边对应的下游节点 P(T|e1), P(T|e2), P(T|e3), P(T|e4)，以及转移概率 P(e1|X1), P(e2|x1), P(e3|X1), P(e4|x1)，即可得到

P(T|X_{1}) = P(e_{1}|X_{1})P(T|e_{1}) + P(e_{2}|X_{1})P(T|e_{2}) \\ \nonumber & \quad\quad + P(e_{3}|X_{1})P(T|e_{3}) + P(e_{4}|X_{1})P(T|e_{4}) \\

(7) 推导 P(T|X)

其中 P(ei|Hj) 是我们之前算好的似然值，P(Hj) 是之前算好的先验概率。

由于 X1 包含 ei，即 ei 成立则 X1 一定也成立，所以

因此 P(T|X1) 和Ｐ(T|X2) 的最终计算公式如下

这里 P(X1) 和 P(X2) 需要额外算一下，用到全概率公式

(8) 计算 P(T|X)

将先验概率，似然值，和刚刚算好的 P(X1), P(X2) 代入公式计算

蒙特卡洛模拟

模拟的逻辑与 $4-21 中的基本一样：模拟的核心是当得到证据 e 之后，如何给出预测结果，也就是代码中 model 实例中的 predict 方法。根据前面的推导，这取决于 16 个似然值。首先随机生成数据模拟出 8 个后验概率(代码中 model 实例的 train 方法)，再随机生成数据模拟预测正确的概率（代码中 model 实例的 test 方法）。只是在预测的时候，增加了对条件概率 P(T|X1) 和 P(T|X2) 的模拟，验证前面的推导过程。

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()
        n_X1 = 0
        n_X2 = 0
        n_X1_corrent = 0
        n_X2_corrent = 0
        for _ in range(T):
            # 随机选硬币
            idx = np.random.randint(0, len(self.coins))
            coin = self.coins[idx]
            # 抛出 n_flip 次，获取证据
            e = []
            # 抛出 n_flip 次
            for _ in range(self.n_flip):
                e.append(int(coin.flip()))
            e = tuple(e)
            if e[0] == 1:
                n_X1 += 1
            else:
                n_X2 += 1
            # 根据证据猜硬币种类
            if self.predict(e) == coin.t:
                correct += 1
                if e[0] == 1:
                    n_X1_corrent += 1
                else:
                    n_X2_corrent += 1
        return correct / T, n_X1_corrent / n_X1, n_X2_corrent / n_X2

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)
ans = [0, 0, 0]
for t in ts:
    ans[0], ans[1], ans[2] = t
print("P(T): {:.6f}".format(ans[0]))
print("P(T|X1): {:.6f}".format(ans[1]))
print("P(T|X2): {:.6f}".format(ans[2]))

模拟结果

证据 (0, 0, 0) 下的正常硬币有8849个, 特殊硬币有216个
证据 (0, 0, 1) 下的正常硬币有8689个, 特殊硬币有908个
证据 (0, 1, 0) 下的正常硬币有8718个, 特殊硬币有988个
证据 (0, 1, 1) 下的正常硬币有8765个, 特殊硬币有3892个
证据 (1, 0, 0) 下的正常硬币有8747个, 特殊硬币有926个
证据 (1, 0, 1) 下的正常硬币有8846个, 特殊硬币有3902个
证据 (1, 1, 0) 下的正常硬币有8709个, 特殊硬币有3810个
证据 (1, 1, 1) 下的正常硬币有8669个, 特殊硬币有15366个
===训练结束===
无额外决策的结果
P(T): 0.766069
P(T|X1): 0.705227
P(T|X2): 0.853590

4.概率面试题（连载）—下-灵析社区

4.17圆周上取 3 个点构成锐角三角形

4.18选票盒

4.19 向正方形区域扔硬币

4.20 不公平的硬币 2

4.21不公平的硬币 3

4.22不公平的硬币 4