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

英勇黄铜

4.8仓促的决斗

在某一天的决斗场，有两个人会在 5 点到 6 点之间随机的某一分钟到来，然后 5 分钟后离开。

如果某一时刻这两个人同时在场，会发生决斗。

问：这两个人发生决斗的概率是多少？

思路参考

记两个人分别为 a 和 b，如图所示，横轴 x 是 a 到达的时间，纵轴 y 是 b 到达的时间，到达的时间范围均为 0 ~ 1。

(x, y) 会随机地在图中的正方形区域中产生，如果两个人的时间间隔小于 1/12 则会发生决斗，也就是 abs(x - y) < 1/12，这对应于图中的阴影部分。

(x, y) 落在阴影的概率是阴影面积与正方形面积的比值。

蒙特卡洛模拟

import numpy as np
from multiprocessing import Pool

def test(T):
    array_arrival_times = 60 * np.random.random(size=(T, 2))
    array_diff = np.abs(np.diff(array_arrival_times))
    return np.mean(array_diff <= 5)

T = int(1e7)
args = [T] * 8
pool = Pool(8)
ts = pool.map(test, args)
print("发生决斗的概率: {:.6f}".format(sum(ts) / len(ts)))

模拟结果

发生决斗的概率: 0.159805

4.9帮助赌徒戒赌

布朗沉迷于轮盘赌，并且强迫症地总是在 13 这个数字上押注 1 美元。

轮盘赌有 38 个等概率的数字，如果玩家押注的数字出现，则拿到 35 倍投注金额，加上原始投注金额，共 36 倍，如果押注的数字没出现，则损失原始投注金额。

好朋友为了帮助布朗戒赌，设计了一个赌注，与布朗等额投注 20 美元（好朋友和布朗均需要出 20 美元），赌以下事情：布朗在第 36 次游戏完成之后会亏钱（如果果真亏钱了，则朋友得到双方的筹码共 40 美元，否则布朗得到这 40 美元）。

问：这个戒赌方案是否有效（如果布朗能挣钱，说明无效，如果布朗会亏钱则说明有效）。

思路参考

布朗每次游戏的获胜概率为 p = 1/38，期望从赌场拿回的钱为 36 * p。

考虑 36 次游戏。布朗需要付出 36 美元与赌场的押注金额，以及 20 美元与朋友的押注金额。下面我们看布朗期望能够拿回多少钱。

布朗拿回的钱中，一部分来自赌场，另一部分来自朋友。

来自赌场的部分，布朗每次期望拿回 36 * p，36 次期望拿回 36 * 36 * p。

下面考虑来自朋友的部分:

每一次获胜，布朗会获得 36 美元。也就是说布朗只要赢一次就已经不亏了，也就意味着他只要不是 36 场全输了，就可以从朋友那里拿回 40 美元。

布朗在 36 次游戏中至少赢一次的概率为:1−(1−p) 36

因此布朗期望从好朋友那里拿回的金额为 40 * (1 - (1 - p) ^ 36)

把两部分加在一起，布朗在 36 次游戏中拿回金额的期望（平均可拿到金额）为E=40(1−(1−p) 36 )+36×36×p

这个数只要大于 36 + 20 = 56 美元，那么朋友设计的这个戒赌方案反而会助长布朗去赌博，戒赌方案就无效。下面我们编程计算这个数

def E(p):
    q = (1 - p) ** 36
    return 40 * (1 - q) + 36 * 36 * p

income = E(1 / 38)
print("Except Income: {:.6f}".format(income))

计算结果

Except Income: 58.790419

由于期望能拿回 58.790419 美元，比 56 美元多，布朗能够挣钱，朋友设计的戒赌方案无效。

蒙特卡洛模拟

from multiprocessing import Pool
import numpy as np
p = 1/38

def game1():
    # 1 次游戏，返回游戏结束拿回的钱
    if np.random.rand() < p:
        return 36
    return 0

def game36():
    # 36 次游戏，返回36次游戏结束拿回的钱
    w = 0
    for _ in range(36):
        w += game1()
    if w >= 36:
        return w + 40
    return w

def test(T):
    np.random.seed()
    total = 0
    for _ in range(T):
        total += game36()
    return total / T

T = int(1e6)
args = [T] * 8
pool = Pool(8)
incomes = pool.map(test, args)
for x in incomes:
    print("Average Income: {:.6f}".format(x))

模拟结果

Average Income: 58.894648
Average Income: 58.697580
Average Income: 58.741948
Average Income: 58.745660
Average Income: 58.757380
Average Income: 58.819856
Average Income: 58.800196
Average Income: 58.767288

4.10 较短的一节木棍

一根棍子，随机选个分割点，将棍子分为两根。

a. 较小的那一根的长度平均是多少？

b. 较小的那一根与较长的那一根的比例平均是多少？

思路参考

假设木棍长度为 1，分隔点 X 是 (0, 1) 上的均匀分布，也就是 X ~ U(0, 1)，概率密度函数 fX(x)=1, 0 < x < 1。

记较短的那一根的长度为 Y，y 与 x 的关系如下

下面我们来看 (a) 和 (b) 两问：

(a)

要求 Y 的期望，我们还需要两个知识点，一个是连续型随机变量的期望，一个是连续型随机变量函数的概率密度函数。

知识点复习

(1) 连续型随机变量的期望

(2) 随机变量函数的概率密度函数

X 的概率密度函数为 f X(x)，Y = g(X)，Y 的概率密度函数是多少？

解决这个问题可借助下面通用流程：

  • 根据 X 的范围（记作 Rx）求出 Y 的范围（记作 Ry）；
  • 对 Ry 中的每一个 y，求出分布函数 FY(y)。【一般用积分的方法】

FY(y)=P(Y≤y)=P(g(X)≤y)=P(X∈Ry)=∫ RyfX(x)dx

其中 x∈Ry 与 g(X)≤y 是相同的随机事件。Ry={x:g(x)≤y}

  • 对分布函数 FY(y) 求导得到密度函数 f Y(y)。

而如果 g(X) 是严格单调函数，x = h(y) 为相应的反函数，且可导，则 Y = g(X) 的密度函数可以直接给出

计算

有了以上两个知识点，我们就可以开始计算了，下面我们推导 Y 的概率密度函数。

g 是一个分段函数，在 (0, 1/2) 上，Y = X，单调递增，在 (1/2, 1) 上，Y = 1 - X，单调递减。因此我们可以对这两个范围分别来看。

(1) 对于 0 < x < 1/2

在 0 < x < 1/2 范围内，y = g(x) = x，因此 0 < y < 1/2，其反函数 x = h(y) = y，导数 h'(y) = 1。y 的密度函数如下

fY​(y)=fX​(h(y))∣h′(y)∣=1

y 的期望如下

(2) 对于 1/2 < x < 1

在 1/2 < x < 1 范围内，y = g(x) = 1 - x，因此 0 < y < 1/2，其反函数 x = h(y) = 1 - y，导数 h'(y) = -1。y 的密度函数如下

fY​(y)=fX​(h(y))∣h′(y)∣=1

y 的期望如下

将两部分加到一起即可得到所求的期望值，答案为 1/4。

(b)

要求较小的那一根与较长的那一根的比例的期望，也就是求 y / (1 - y) 的期望。

由 (a)，我们已经知道 Y 的范围是 (0, 1/2)，并且 f Y(y) 根据 X 的取值来源于两部分，第一部分是 1 < X < 1/2，第二部分是 1/2 < X < 1。将 X 取值在这两个范围时的 f Y(y) 相加，可以得到：

f Y (y)=2y∈(0,1/2)

这里 f Y(y)=2，但是 y 的取值是 (0, 1/2)，这是一个均匀分布。基于 Y，我们定义随机变量 Z = Y / (1 - Y)

g(y) 是严格单调的，其反函数为 h(z) = z / (z + 1)，其导数为 ℎ′(z)=(z+1) −2 , z 的范围为 (0, 1)，因此 z 的密度函数如下

fZ​(z)=fY​(h(z))∣h′(z)∣=2×(z+1)−2

Z 的期望如下

蒙特卡洛模拟

import numpy as np

def pyfunc(x):
    """
    0 < x < 1
    """
    if x > 0.5:
        x = 1 - x
    return x

def pyfunc2(x):
    """
    0 < x < 0.5
    """
    return x / (1 - x)

func = np.frompyfunc(pyfunc, 1, 1)
func2 = np.frompyfunc(pyfunc2, 1, 1)

def test(T):
    X = np.random.uniform(0, 1, T)
    X = func(X)
    Y = func2(X)
    return (np.mean(X), np.mean(Y))

T = int(1e7)
m1, m2 = test(T)
print("较短的木棍的期望长度: {:.6f}".format(m1))
print("较短的木棍与较长的木根长度比值的期望: {:.6f}".format(m2))

模拟结果

较短的木棍的期望长度: 0.249920
较短的木棍与较长的木根长度比值的期望: 0.386187

4.11完美手牌

从一副洗好的牌（52张）中抽 13 张，拿到完美手牌（拿到的 13 张是同花色）的概率是多少。

思路参考

52 张牌一共有 52! 中排列顺序。

完美手牌需要抽取的 13 张是同一花色，这 13 张同一花色的牌有 13! 种排列顺序，未被抽到的 39 张牌有 39! 中排列顺序。

因此对于指定花色，所有抽到的牌都是这一指定花色的概率为

一共有四种花色，因此拿到完美手牌的概率为

编程计算

from math import factorial as fac

p = 4.0 * fac(13) * fac(39) / fac(52)
print("{:.6e}".format(p))

计算结果

6.299078e-12

4.12双骰子赌博

双骰子赌博是美国的流行游戏，规则如下：

每次同时投掷两个骰子并且只看点数和。第一次投掷的点数和记为 x。若 x 是 7 或 11，玩家获胜，游戏结束；若 x 是 2、3 或 12，则输掉出局，其余情况则需要继续投掷。从第二次投掷开始，如果投掷出第一次的点数和 x，则获胜；如果投掷出 7，则输掉出局。重复投掷直至分出胜负。

问：获胜的概率是多少？

思路参考

首先考虑第一次投掷

第一次投掷一共有直接输，直接赢，需要继续投掷三种可能性。

(1) 直接输的情况是两枚骰子的和为 2、3、12，两个骰子各自的点数情况有 (1, 1), (2, 1), (1, 2), (6, 6) 这 4 中可能性。而两个骰子各自点数共有 6 * 6 = 36 种可能性，因此第一次投掷直接输的概率为 1/9。

(2) 直接赢的情况是两枚骰子的和为 7, 11，两个骰子各自的点数情况有 (1, 6), (2, 5), (3, 4), (4, 3), (5, 2), (6, 1), (5, 6), (6, 5) 共 8 种可能性，概率 2/9。

(3) 剩下 2/3 概率的情况，需要继续投掷，我们记此前的第 1 次投掷的点数和为 x。

考虑第二次以及后续可能的投掷

从第 2 次开始，以后每次投掷的输、赢、继续投掷的条件都一样。因此我们可以只考虑第 2 次投掷，如果投掷结果是需要继续投掷，递归处理即可。

记某次投掷直接赢的概率为 p(x)，与第一次的投掷结果 x 有关，直接输的概率为 q，与第一次的投掷结果 x 无关，则需要继续投掷的概率为 (1 - p(x) - q)，下面我们看输和赢的具体条件：

(1) 输的情况是投掷出 7，两个骰子各自的点数情况共有 (1, 6), (2, 5), (3, 4), (4, 3), (5, 2), (6, 1) 共 6 种情况，概率 q = 1 / 6。

(2) 赢的情况是投掷出 x。两个骰子各自的点数情况需要看 x。下面我们枚举 x 可能取值（4、5、6、8、9、10）：

  • x = 4: (1, 3), (2, 2), (3, 1)，概率为 1/12；
  • x = 5: (1, 4), (2, 3), (3, 2), (4, 1)，概率为 1/9；
  • x = 6: (1, 5), (2, 4), (3, 3), (2, 4), (5, 1)，概率为 5/36；
  • x = 8: (2, 6), (3, 5), (4, 4), (5, 3), (6, 2)，概率为 5/36；
  • x = 9: (3, 6), (4, 5), (5, 4), (6, 3)，概率为 1/9；
  • x = 10: (4, 6), (5, 5), (6, 4)，概率为 1/12。

总结一下，直接赢的概率 p(x) 与 x 的所有可能取值的对应关系如下

p(4) = 1/12
p(5) = 1/9
p(6) = 1/36
p(8) = 1/36
p(9) = 1/9
p(10) = 1/12

(3) 剩下的情况是需要继续投掷的，概率为 (1 - p(x) - q)

综合以上 3 条，从第二次投掷开始，最终能赢有两种情况，一种是第二次投掷直接赢，另一种是需要继续投掷，并在后续投掷中赢了。记 P(x) 为从第二次投掷开始，最终会赢的概率。

将第一次投掷直接赢和进入第二次投掷并最终赢这两种情况合起来，就是最终能赢的两种情况，概率记为 prob

下面编程计算 prob 这个数

mapping = {
    4: 1/12,
    5: 1/9,
    6: 5/36,
    8: 5/36,
    9: 1/9,
    10: 1/12,
}

q = 1/6

prob = 2/9
for x, p in mapping.items():
    P = p / (p + q)
    prob += p * P
print("Probability of win: {:.6f}".format(prob))

计算结果

Probability of win: 0.492929

蒙特卡洛模拟

from multiprocessing import Pool
import numpy as np

def throw2(x0):
    # 第二次以后的投掷
    x1 = np.random.randint(1, 7)
    x2 = np.random.randint(1, 7)
    x = x1 + x2
    if x == 7:
        return 0
    elif x == x0:
        return 1
    else:
        return throw2(x0)

def throw1():
    # 第一次投掷
    x1 = np.random.randint(1, 7)
    x2 = np.random.randint(1, 7)
    x = x1 + x2
    if x == 7 or x == 11:
        return 1
    elif x == 2 or x == 3 or x == 12:
        return 0
    else:
        return throw2(x)

def test(T):
    np.random.seed()
    n = 0
    for _ in range(T):
        n += throw1()
    return  n / T

T = int(1e7)
args = [T] * 8
pool = Pool(8)
probs = pool.map(test, args)
for prob in probs:
    print("Probability of win: {:.6f}".format(prob))

模拟结果

Probability of win: 0.492935
Probability of win: 0.492876
Probability of win: 0.492793
Probability of win: 0.492726
Probability of win: 0.493062
Probability of win: 0.493019
Probability of win: 0.493142
Probability of win: 0.492961

4.13 收集优惠券

每个盒子中有一张优惠券，优惠券共有 5 种，编号为 1 到 5。必须收集所有5张优惠券才能获得奖品。

问：收集到一套完整的优惠券平均需要打开多少盒子？

思路参考

单独考虑每个优惠券。

记 T1 为一个随机变量，表示拿到第一个编号所需盒子数。P(T1 = 1) = 1，E(T1) = 1

记 T2 为一个随机变量，表示拿到第二个编号所需要的额外盒子数。根据全期望公式

这里 X 表示新开盒子编号是否见过，共有两种情况：

  • 前面开盒子时见过，概率为 1/5，此时平均需要的额外盒子数（即期望）为 1 + E(T2)；【可参考$4-4 方法二、有向图总结】
  • 前面开盒子时没见过，概率为 4/5，此时平均需要的额外盒子数为 1。

记 T3 为一个随机变量，表示从第二个编号开始到拿到第三个新编号所需的额外盒子数。根据全期望公式

这里的 X 是新开的盒子的数字是否是见过的那 2 个。共有两种情况

  • 是前面见过的那2个，概率 2/5，此时额外的盒子数的期望是 1 + E(T3)
  • 不是前面见过的那2个，概率 3/5，此时额外的盒子数就是 1

记 T4 为一个随机变量，表示从第三个编号开始到拿到第四个新编号所需的额外盒子数。根据全期望公式

这里 X 表示新开盒子编号是否已出现过，分为两种情况：

  • 前面开盒子时见过，概率为 3/5，此时平均需要额外盒子数（期望）为 1 + E(T4)；
  • 前面开盒子时没见过，概率为 2/5，此时平均需要额外盒子数（期望）为 1。

记 T5 为一个随机变量，表示从第四个编号开始到拿到第五个新编号所需的额外盒子数。根据全期望公式

这里 X 表示新开盒子编号是否已出现过，分为两种情况：

  • 前面开盒子时见过，概率为 4/5，此时平均需要额外盒子数（期望）为 1 + E(T5)；
  • 前面开盒子时没见过，概率为 1/5，此时平均需要额外盒子数（期望）为 1。

收集到一套完整的优惠券需要打开的盒子个数的期望为 E1 + E2 + E3 + E4 + E5 = 137/12 = 11.416667

蒙特卡洛模拟

import numpy as np
import operator
from multiprocessing import Pool
from functools import reduce

def run():
    setting = set()
    n = 0
    while True:
        n += 1
        x = np.random.randint(1, 6)
        if x in setting:
            continue
        setting.add(x)
        if len(setting) == 5:
            break
    return n

def test(T):
    np.random.seed()
    y = (run() for _ in range(T))
    n = reduce(operator.add, y)
    return n / T

T = int(1e7)
args = [T] * 8
pool = Pool(8)
ts = pool.map(test, args)
for t in ts:
    print("{:.6f} coupons on average are required".format(t))
11.415423 coupons on average are required
11.417857 coupons on average are required
11.415529 coupons on average are required
11.417148 coupons on average are required
11.416793 coupons on average are required
11.413937 coupons on average are required
11.416149 coupons on average are required
11.416380 coupons on average are required

4.14 第 2 强的选手拿亚军

一场由 8 个人组成的淘汰赛，对阵图如上图。8 个人随机被分到图中的 1~8 号位置。

第二强的人总是会输给最强的人，同时也总是会赢剩下的人。决赛输的人拿到比赛的亚军。

问: 第二强的人拿到亚军的概率是多少？

思路参考

第二强的人拿到亚军，等价于最强的人和第二强的人在决赛遇到，等价于这两个人一个被分到 1~4，另一个被分到 5~8。

x1: 最强的人被分到 1~4

x2: 最强的人被分到 5~8

y1: 第二强的人被分到 1~4

y2: 第二强的人被分到 5~8

第二强的人获得亚军的概率为 2/7+2/7=0.571428

蒙特卡洛模拟

import numpy as np
import operator
from multiprocessing import Pool
from functools import reduce

ladder = range(8)

def run():
    l = np.random.permutation(ladder)
    x = int(np.where(l == 0)[0][0])
    y = int(np.where(l == 1)[0][0])
    return (x <= 3 and y >= 4) or (x >= 4 and y <= 3)

def test(T):
    np.random.seed()
    y = (run() for _ in range(T))
    n = reduce(operator.add, y)
    return n / T

T = int(1e7)
args = [T] * 8
pool = Pool(8)
ts = pool.map(test, args)
for t in ts:
    print("P(the second-best player wins the runner-up cup): {:.6f}".format(t))

模拟结果

P(the second-best player wins the runner-up cup): 0.571250
P(the second-best player wins the runner-up cup): 0.571030
P(the second-best player wins the runner-up cup): 0.571432
P(the second-best player wins the runner-up cup): 0.571021
P(the second-best player wins the runner-up cup): 0.571379
P(the second-best player wins the runner-up cup): 0.571495
P(the second-best player wins the runner-up cup): 0.571507
P(the second-best player wins the runner-up cup): 0.571226

4.15 孪生骑士

(a)

8 个骑士进行淘汰赛，对阵图与上图这样的网球对阵图类似。8 个骑士的水平一样，任意一组比赛双方的获胜概率均为 0.5。8 个骑士中有两个人是孪生兄弟。

问：这两个孪生兄弟在比赛过程中相遇的概率？

(b)

将 8 个骑士的淘汰赛改为 2^n 个人的淘汰赛，问这两个孪生兄弟相遇的概率？

思路参考

(a)

8 个骑士的比赛共有 3 轮，我们一轮一轮地考虑。

记事件 X1 为孪生骑士在第一轮相遇，若想在第一轮相遇，则它们必须同时被分到 (1, 2), (2, 1), (3, 4), (4, 3), (6, 5), (5, 6), (8, 7), (7, 8) 这八种情况之一。概率为

记事件 X2 为孪生骑士在第二轮相遇，若想在第二轮相遇，则需要它们两个被分到以下十六种情况的一种，同时它们两个还要都赢了

(1, 3), (1, 4), (2, 3), (2, 4), (5, 7), (5, 8), (6, 7), (6, 8) 
(3, 1), (4, 1), (3, 2), (4, 2), (7, 5), (8, 5), (7, 6), (8, 6)

概率为

记事件 X3 为孪生骑士在第三轮相遇，若想在第三轮相遇，需要它们两个其中一个被分到 [1..4] 中的一个，另一个被分到 [5..8] 中的一个。且他们在前两轮的共四场比赛都要赢，概率为

把 X1, X2, X3 这三种情况的概率加起来，得到孪生骑士相遇的概率

(b)

记事件 Xi 为孪生骑士在第 i 轮相遇，其中 i 取值为 1, 2, ..., n，共 T = 2^n 人。

对第 i 轮，将 1 ~ 2^n 连续地分为 B = (2^n)/(2^(i-1)) 个桶，每个桶 C = 2^(i-1) 个元素。孪生骑士需要在编号为 (2k+1, 2k+2) 的相邻桶中（k=0,1,...,(2^n)/(2^i) - 1），相邻桶的组数为 N = (2^n)/(2^i)，孪生骑士在满足要求的桶中的概率为

孪生骑士在第 1,2,...,i-1 轮中，共有 (i-1) * 2 场比赛。这些比赛还要赢，概率为

将两个概率相乘就是孪生骑士在第 i 轮相遇的概率

孪生骑士相遇的概率为

蒙特卡洛模拟

import numpy as np
import operator
from multiprocessing import Pool
from functools import reduce

def P(n):
    return 1 / (2 ** (n - 1))

def run(n):
    l = np.random.permutation(ladder)
    for _ in range(1, n + 1):
        for i in range(0, len(l), 2):
            if l[i] + l[i + 1] == 1:
                return 1
            if np.random.rand() < 0.5:
                l[int(i / 2)] = l[i]
            else:
                l[int(i / 2)] = l[i + 1]
        l = l[:int(len(l) / 2)]
    return 0

def test(n):
    T = int(1e5)
    np.random.seed()
    y = (run(n) for _ in range(T))
    return reduce(operator.add, y) / T

def main(n):
    print("n = {}: P(n) = {}".format(n, P(n)))
    args = [n] * 8
    pool = Pool(8)
    ts = pool.map(test, args)
    for t in ts:
        print("P(twin knight meet): {:.6f}".format(t))
    print()

for n in range(3, 7):
    global ladder
    ladder = range(2 ** n)
    main(n)

模拟结果

n = 3: P(n) = 0.25           
P(twin knight meet): 0.249970
P(twin knight meet): 0.249815
P(twin knight meet): 0.250011
P(twin knight meet): 0.249984
P(twin knight meet): 0.249990
P(twin knight meet): 0.250079
P(twin knight meet): 0.249889
P(twin knight meet): 0.250108 

n = 4: P(n) = 0.125
P(twin knight meet): 0.124822
P(twin knight meet): 0.124963
P(twin knight meet): 0.124807
P(twin knight meet): 0.125172
P(twin knight meet): 0.124867
P(twin knight meet): 0.124937
P(twin knight meet): 0.124887
P(twin knight meet): 0.124992

n = 5: P(n) = 0.0625
P(twin knight meet): 0.062401
P(twin knight meet): 0.062294
P(twin knight meet): 0.062492
P(twin knight meet): 0.062508
P(twin knight meet): 0.062588
P(twin knight meet): 0.062417
P(twin knight meet): 0.062393
P(twin knight meet): 0.062459

n = 6: P(n) = 0.03125
P(twin knight meet): 0.031231
P(twin knight meet): 0.031373
P(twin knight meet): 0.031206
P(twin knight meet): 0.031259
P(twin knight meet): 0.031189
P(twin knight meet): 0.031246
P(twin knight meet): 0.031311
P(twin knight meet): 0.031256

4.16有放回抽样还是无放回抽样

两个桶 A, B，桶内有红色和黑色的球，其中：

  • A 桶有 2 个红球和 1 个黑球。
  • B 桶有 101 个红球和 100 个黑球。

随机选出一个桶放到你面前，然后从中抽样两个球，基于这两个球的颜色猜该桶是 A 还是 B。

你可以选择有放回抽样或者无放回抽样，也就是你可以决定在抽取第二个球之前，是否将抽出的第一个球放回去再抽第二次。

问：为了使得猜中桶名的概率最大，应该有放回抽样还是无放回抽样，猜中的概率是多少？

思路参考

本题的目标是根据抽样两个球的结果，来猜当前是哪个桶。有两种获取证据的方式，一种是有放回抽样抽两个球，另一种是无放回抽样抽两个球。我们要选择使得猜中概率更大的那个抽样方式。

根据第 3 章总结，我们应先确定“证据是什么”、“假设有哪些”，然后计算似然值，再推公式计算。

(1) 假设和证据

首先我们看都有哪些可能的证据。由于抽样两个球，且球的种类有红和黑，因此无论有放回抽样还是无放回抽样，观察到的证据都可能有 ”两红”、”一红一黑”、”两黑” 这三种，分别记为 e1, e2, e3。

然后我们看都有哪些可能的假设。由于有 A 和 B 这两种桶，因此针对随机取出的桶有两种假设：桶名为 A 和 桶名为 B，分别记为 H1, H2。

不论抽样是否有放回，都有上述的三种证据和两种假设。

(2) 先验概率

由于一共两个桶，随机抽取一个，因此抽到 A 或 B 的先验概率为 1/2，与抽样是否有放回无关。

(3) 似然值

  • 有放回抽样

对于 A 桶，有放回抽样两个球，得到 e1, e2, e3 三种证据的概率分别为

对于 B 桶，有放回抽样两个球，得到 e1, e2, e3 三种证据的概率分别为

  • 无放回抽样

对于 A 桶，无放回抽样两个球，得到 e1, e2, e3 三种证据的概率分别为

对于 B 桶，无放回抽样两个球，得到 e1, e2, e3 三种证据的概率分别为

(4) 推导 P(T)

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

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

(5) 计算

我们分别计算在有放回抽样和无放回抽样 2 次时，预测正确的概率 P(T)

  • 有放回抽样

  • 无放回抽样

蒙特卡洛模拟

import numpy as np
from multiprocessing import Pool

class Box:
    def __init__(self, nR, nB, t):
        # 当前子在盒子中的
        self.nR = nR
        self.nB = nB
        # 抽出未放回的
        self.mR = 0
        self.mB = 0
        self.t = t

    def sample_with_replacement(self):
        if np.random.randint(0, self.nR + self.nB) < self.nR:
            return "R"
        else:
            return "B"

    def sample_without_replacement(self):
        if np.random.randint(0, self.nR + self.nB) < self.nR:
            self.nR -= 1
            self.mR += 1
            return "R"
        else:
            self.nB -= 1
            self.mB += 1
            return "B"

    def reset(self):
        self.nB += self.mB
        self.mB = 0
        self.nR += self.mR
        self.mR = 0

    def sample(self, method):
        if method == "with":
            return self.sample_with_replacement()
        else:
            return self.sample_without_replacement()

class Model:
    def __init__(self, boxs, method):
        self.boxs = boxs
        self.method = method
        # mapping 记录每种证据下两种盒子种类的次数
        # 共有 6 个参数
        self.mapping = {"RR": np.array([0, 0])
                       ,"RB": np.array([0, 0])
                       ,"BB": 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.boxs))
            box = self.boxs[idx]
            # 抽取两次，获取证据
            e = [box.sample(self.method), box.sample(self.method)]
            if e[0] == "B" and e[1] == "R":
                e[1], e[0] = e[0], e[1]
            e = "".join(e)
            self.update(e, box.t)
            box.reset()
        for k, v in self.mapping.items():
            print("证据 {} 下的A桶有{}个, B桶有{}个".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.boxs))
            box = self.boxs[idx]
            # 抽取 2 次，获取证据
            e = [box.sample(self.method), box.sample(self.method)]
            if e[0] == "B" and e[1] == "R":
                e[1], e[0] = e[0], e[1]
            e = "".join(e)
            # 根据证据猜硬币种类
            if self.predict(e) == box.t:
                correct += 1
            box.reset()
        return correct / T

class Solver:
    def __init__(self, method):
        self.method = method
        self.boxs = [Box(2, 1, 0), Box(101, 100, 1)]
        self.ans = -1

    def __call__(self):
        model = Model(self.boxs, self.method)
        model.train()
        T = int(1e7)
        args = [T] * 8
        pool = Pool(8)
        ts = pool.map(model.test, args)
        self.ans = sum(ts) / len(ts)

solver1 = Solver("with")
solver1()
print("有放回时，P(T): {:.6f}".format(solver1.ans))
solver2 = Solver("without")
solver2()
print("无放回时，P(T): {:.6f}".format(solver2.ans))

模拟结果

证据 RR 下的A桶有22280个, B桶有12755个
证据 RB 下的A桶有21997个, B桶有24808个
证据 BB 下的A桶有5561个, B桶有12599个
===训练结束===
有放回时，P(T): 0.596081
证据 RR 下的A桶有16452个, B桶有12755个
证据 RB 下的A桶有33433个, B桶有25112个
证据 BB 下的A桶有0个, B桶有12248个
===训练结束===
无放回时，P(T): 0.623070


阅读量：2035

点赞量：0

收藏量:0