数概率

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

字数统计: 902字   |   阅读时长: 4分
,

摘要: 《概率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