【Puzzle】萨缪尔·佩皮斯的问题

  |  

摘要: 《概率50题》掷骰子的概率问题

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


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

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


问题描述

Isaac Newton Helps Samuel Pepys

Pepys wrote Newton to ask which of the following is most likely

a. A person rolls at least one 6 when six dice are rolled

b. A person rolls at least two 6s when twelve dice are rolled

c. A person rolls at least three 6s when eighteen dice are rolled

What is the answer?

萨缪尔佩皮斯写信问牛顿。以下三件事哪件的概率更大

  • 投掷 6 枚骰子至少一个是 6
  • 投掷 12 枚骰子至少两个是 6
  • 投掷 18 枚骰子至少三个是 6

思路参考

全概率公式 + 概率 DP

假设我们要求投掷 n 枚骰子,至少 m 枚是 6 的概率。

如果 m 是 0,则概率是 1,因为不论最终多少枚是 6,都可以算是至少 0 枚。

如果 m 大于 n,则概率是 0,因为不可能投掷出 6 的个数比骰子的个数还多。

当 0 < m < n 时,我们考虑其中一枚骰子的两种情况:

(1) 该骰子是 6,概率 1/6,则剩下的 n - 1 枚骰子需要至少 m - 1 枚是 6
(2) 该骰子不是 6,概率 5/6,则剩下的 n - 1 枚骰子需要至少 m 枚是 6

综合以上分析,再结合全概率公式

我们可以用概率DP的方式解决,算法如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
状态定义
dp[i][j] := 投掷 i 枚骰子,至少 j 个是 6 的概率

初始化
dp[i][0] = 1.0
dp[i][j] = 0 (j > i)

答案
dp[6][1]
dp[12][2]
dp[18][3]

状态转移
dp[i][j] = 1/6 * dp[i - 1][j - 1] + 5/6 * dp[i - 1][j]

下面编程计算 dp[n][m]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#include <vector>
#include <iostream>

using namespace std;

int main()
{
int n, m;
cin >> n >> m;
vector<vector<double>> dp(n + 1, vector<double>(m + 1, 0.0));
for(int i = 0; i <= n; ++i)
dp[i][0] = 1.0;
for(int i = 1; i <= n; ++i)
for(int j = 1; j <= min(i, m); ++j)
dp[i][j] = 1 / 6.0 * dp[i - 1][j - 1] + 5 / 6.0 * dp[i - 1][j];
cout << dp[n][m] << endl;
}

计算结果:

1
2
3
dp[6][1] = 0.665
dp[12][2] = 0.619
dp[18][3] = 0.597。

因此第一个事件的概率更大。


蒙特卡洛模拟

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np
import operator
from multiprocessing import Pool
from functools import reduce

def run(n, m):
l = np.random.randint(1, 7, n)
return np.where(l == 6)[0].shape[0] >= m

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

args = [(6, 1)] * 3 + [(12, 2)] * 3 + [(18, 3)] * 3
pool = Pool(8)
ts = pool.map(test, args)
for i in range(len(args)):
t = ts[i]
n, m = args[i]
print("P(at least {} 6s when {} dice are rolled): {:.6f}".format(n, m, t))

模拟结果

1
2
3
4
5
6
7
8
9
P(at least 6 6s when 1 dice are rolled): 0.665889
P(at least 6 6s when 1 dice are rolled): 0.664726
P(at least 6 6s when 1 dice are rolled): 0.664759
P(at least 12 6s when 2 dice are rolled): 0.618201
P(at least 12 6s when 2 dice are rolled): 0.618719
P(at least 12 6s when 2 dice are rolled): 0.618631
P(at least 18 6s when 3 dice are rolled): 0.597280
P(at least 18 6s when 3 dice are rolled): 0.597390
P(at least 18 6s when 3 dice are rolled): 0.596858

Share