算法动态规划

概率DP

字数统计: 2.2k字   |   阅读时长: 9分
,

摘要: 本文介绍概率 DP 的原理和例题

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

概率DP并不是一种特殊的动态规划类型，而是一种常见的应用场景。类似的还有计数DP，博弈DP，等等。实际上在解决概率DP问题时候，最终用的DP类型，根据对问题建模后状态之间的关系，可能是线性DP，树形DP，图上DP等等，与具体问题有关系。

概率 DP 的理论基础主要是全概率公式

  • 全概率公式

概率 DP 的转移方程一般是下面这样，当然会根据问题的不同产生很多变化，比如会需要附加状态等等。

其中 dp[i] 对应了全概率公式中的 P(A)，p[i][j] 对应 P(Bj)，dp[j] 对应 P(A|Bj)

本文我们详细拆解两个概率DP的经典问题，锦标赛和招标。

另外 【连载】概率面试题 和 概率与期望计算 也有一些期望 DP 的问题。

问题1: 锦标赛

有 2^n 个队伍，排成一列，相邻的之间两两打淘汰赛。此后每轮比赛前各个队伍均按照升序排列，然后相邻的之间两两比赛。

有一个概率矩阵 P，Pij 为 i, j 比赛时，i 取胜的概率。

问哪个队最终获胜的概率最大。

测试数据，第一行为 n，此后 2^n 行，每行 2^n 个数为 P 矩阵。1 <= n <= 7

1
2
3
4
5
6
7
8
输入:
2
0.0 0.1 0.2 0.3
0.9 0.0 0.4 0.5
0.8 0.6 0.0 0.6
0.7 0.5 0.4 0.0
输出
2

思路(概率DP)

如果有 2^n 个队伍，需要比 n 轮。

以 n = 3 为例，共 8 个队伍，编号 0 ~ 7，三轮比赛如下

j = 0,1,2 为比赛轮次，i = 0,1,…,7 为队伍编号。

假设最终 i 队伍获得冠军，那么 i 队伍需要先到第 j = 2 轮比赛，然后在该轮比赛中获胜。

其中涉及到两件事，第一件是 i 队伍需要打到第 j = 2 轮；第二件事是 i 队伍赢下该轮他的对手。

第一件事比较好解决，它只需要递归地变成 i 队伍打到第 j = 1 轮并赢了该轮的对手。

对于第二件事，又可以拆成两个子问题，第一个子问题是 j = 2 轮的对手是哪个，第二个子问题是 i 队伍赢下该对手的概率是多少。

其中第二个子问题比较好解决，直接查矩阵 P 即可。

现在就剩下一个问题：队伍 i 在第 j 轮的对手可能是谁。

要直接说出队伍 i 在第 j 轮的所有可能对手会有点复杂，但是给定一个队伍 k，判断 i 和 k 是否有可能在第 j 轮称为对手就比较简单了。

具体做法是将 i 和 k 向右移 j 位，然后看移位后的 i 和 k 异或是否为 1。

以 n = 3, j = 0, 1, 2 为例，8 个队伍的二进制分别为

1
2
3
4
5
6
7
8
000
001
010
011
100
101
110
111

这步判断的代码如下

1
2
3
4
5
6
7
bool check(int i, int j, int k)
{
    // 判断 i, k 是否可能在第 j 轮成为对手
    i >>= j;
    k >>= j;
    return i ^ k == 1;
}

基于以上分析，我们写出概率 DP 算法如下

1
2
3
4
5
6
7
8
9
10
11
状态定义
dp[i][j] := 第 j 轮 i 获胜

答案
max(dp[n - 1][i])

初始化
dp[i][0] = p[i][k], if check(i, 0, k)

状态转移
dp[i][j] = sum(dp[i][j - 1] * dp[k]dp[j - 1] * p[i][k])

代码(C++)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#include <cmath>
#include <iostream>
#include <vector>

using namespace std;

bool check(int i, int j, int k)
{
    // 判断 i, k 是否可能在第 j 轮成为对手
    i >>= j;
    k >>= j;
    return i ^ k == 1;
}

int main()
{
    int n;
    cin >> n;
    int N = pow(2, n);
    vector<vector<double>> P(N, vector<double>(N, -1.0));
    for(int i = 0; i < N; ++i)
        for(int j = 0; j < N; ++j)
            cin >> P[i][j];
    vector<vector<double>> dp(N, vector<double>(n, -1.0));
    for(int i = 0; i < N; ++i)
    {
        if(i & 1)
            dp[i][0] = P[i][i - 1];
        else
            dp[i][0] = P[i][i + 1];
    }
    for(int j = 1; j < n; ++j)
    {
        for(int i = 0; i < N; ++i)
            for(int k = 0; k < N; ++k)
                if(check(i, j, k))
                    dp[i][j] = dp[i][j - 1] * dp[k][j - 1] * P[i][k];
    }
    double max_p = -1.0;
    int ans = -1;
    for(int i = 0; i < n; ++i)
        if(dp[i][n - 1] > max_p)
        {
            max_p = dp[i][n - 1];
            ans = i;
        }
    cout << ans + 1 << endl;
}

问题2: 招标

客户把 M 个需求打包到一个项目，进行招标，能解决最多问题的中标。

共 T 个供应商竞标。

第 i 个供应商能解决第 j 个问题的概率为 P[i][j]

问：每个竞标供应商都能解决至少一个问题，且中标者最少解决 N 个问题的概率。

测试数据，第一行为 M, T, N，接下来是 T 行 M 列，为 P 矩阵

1
2
3
4
5
6
输入
2 2 2
0.9 0.9
1 0.9
输出
0.972

思路(概率DP)

对于特定的一个供应商，他最终能解决的问题可能是 0, 1, …, M，共 M + 1 种可能。

我们首先需要用概率 DP 的方式求出这些概率。对于特定供应商 0 <= t <= T-1，解决每个问题的概率为 P[t][0..M-1]

1
2
3
4
5
6
7
8
9
10
11
12
状态定义
dp[i][j] := 第 [0..i] 这些问题，总共解决 j 个问题的概率

初始化
dp[0][0] = (1 - P[t][0])
dp[0][1] = P[t][0]

答案
sum(dp[M - 1][0..N-1])

状态转移
dp[i][j] = P[t][i] * dp[i - 1][j - 1] + (1 - P[t][i]) * dp[i - 1][j]

对于特定 t, 以上 DP 流程走完后，dp[M - 1][0..M] 就是供应商 t 最终能解决 0, 1, …, M 个问题的概率。

我们要求的是所有供应商都至少解决一个问题，且至少有一个供应商解决至少 N 个问题。

记以下两个事件符号

A = 供应商都至少解决一个问题
B = 至少有一个供应商解决至少 N 个问题

所求为 P(AB)，但由于 A, B 并不独立，不能直接将 P(A)P(B) 作为答案。因此需要用条件概率 P(AB) = P(A)P(B|A)。

P(A) 比较简单。对每个供应商，将 dp[M - 1][1..M] 加起来，就是该供应商至少解决一个问题的概率，然后将所有供应商的概率相乘即可。

P(B|A) 比较复杂，我们首先看没有 A 这个条件会怎样。

如果是 P(B)，对于每个供应商，将 dp[M - 1][0..N-1] 相加就是该供应商未能解决至少 N 个问题的概率将每个供应商的概率相乘即可得到所有供应商都未能解决至少 N 个问题的概率，然后用 1 减一下得到 P(B)

增加 A 条件后，对每个供应商，都不会有 dp[M - 1][0] 这个选项，因此算概率时候的分母就要将这一项减掉。

具体地，对于特定供应商，sum(dp[M - 1][1..N-1]) / (1 - dp[M - 1][0]) 即为该供应商在至少一个问题的条件下未能解决至少 N 个问题的概率。

将所有供应商的 sum(dp[M - 1][1..N-1]) / (1 - dp[M - 1][0]) 相乘即为所有供应商都在至少解决一个问题的条件下未能解决至少 N 个问题的概率。再用 1 减一下即得 P(B|A)

代码(C++)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include <vector>
#include <iostream>
#include <iomanip>

using namespace std;

int main()
{
    int M, T, N;
    cin >> M >> T >> N;
    vector<vector<double>> P(T, vector<double>(M, -1.0));
    for(int i = 0; i < T; ++i)
        for(int j = 0; j < M; ++j)
            cin >> P[i][j];
    vector<vector<double>> dp(M, vector<double>(M + 1, -1.0));
    double pa = 1.0;
    double pb_a = 1.0;
    for(int t = 0; t < T; ++t)
    {
        dp.assign(N, vector<double>(M + 1, -1.0));
        dp[0][0] = (1 - P[t][0]);
        dp[0][1] = P[t][0];
        for(int i = 1; i < M; ++i)
        {
            dp[i][0] = (1 - P[t][0]) * dp[i - 1][0];
            for(int j = 1; j <= i; ++j)
                dp[i][j] = P[t][i] * dp[i - 1][j - 1] + (1 - P[t][i]) * dp[i - 1][j];
        }
        pa *= 1 - dp[M - 1][0];
        double tmp = 0.0;
        for(int j = 1; j < N; ++j)
            tmp += dp[M - 1][j];
        pb_a *= tmp / (1 - dp[M - 1][0]);
    }
    pb_a = 1 - pb_a;
    double ans = pa * pb_a;
    cout << std::fixed << std::setprecision(4);
    cout << ans << endl;
}
Share