算法动态规划

期望DP

字数统计: 1.8k字   |   阅读时长: 8分
,

摘要: 本文介绍期望 DP 的原理和例题

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

与概率DP类似(参考文章概率DP)，期望DP不是一种特殊的动态规划类型，而是一种常见的应用场景。类似的还有计数DP，博弈DP，等等。实际上在解决期望DP问题时候，最终用的DP类型可能是线性DP，树形DP，图上DP等等，与具体问题有关系。

期望 DP 有一些理论基础，主要是线性性质和全期望公式。

  • 期望的定义
  • 期望的线性性质

对于任意随机变量 X 和 Y

当 X, Y 独立时，E(XY) = E(X)(Y)

  • 全期望公式

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

其中 dp[i] 对应了全期望公式中的 E(Y)，p[i][j] 对应 P(xj)，dp[j] 对应 E(Y|xj)

本文我们详细拆解期望DP最经典的问题，涂砖块。

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

涂砖块1

n 个格子，每次随机涂一个，可能重复涂，求涂满 m 个格子的期望次数。

分析

dp[i] := 如果已经涂满了 i 个格子，再涂满 m - i 个格子所需次数的期望

现在执行一次涂色，有两种可能

一种是涂到了之前涂过的格子，概率 i / n，此时又回到了 dp[i] 的状态

另一种是涂到了之前没涂过的格子，概率 (n - i) / n，那么现在就涂满 i + 1 个格子了，到达 dp[i + 1] 的状态

状态转移方程如下：

1
dp[i] = 1 + i / n * dp[i] + (n - i) / n * dp[i + 1]

基于以上分析，写出期望DP如下

1
2
3
4
5
6
7
8
9
10
11
状态定义
dp[i] := 如果已经涂满了 i 个格子，再涂满 m - i 个格子所需次数的期望

初始化
dp[n] = 0

答案
dp[0]

状态转移
dp[i] = (n + (n - i) * dp[i + 1]) / (n - i);

答案验证

期望DP的实现代码如下

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
#include <vector>
#include <fstream>
#include <iostream>

using namespace std;

double solve(int i, const int n, const int m, vector<double>& dp)
{
    if(dp[i] > -0.5)
        return dp[i];
    if(i == m)
        return dp[i] = 0.0;

    return dp[i] = (n + (n - i) * solve(i + 1, n, m, dp)) / (n - i);
}

int main()
{
    ifstream fin("data.txt");
    vector<double> dp;
    int n, m;
    while((fin >> n) && (fin >> m))
    {
        cout << "n: " << n << ", m: " << m << endl;
        dp.assign(m + 1, -1.0);
        double ans = solve(0, n, m, dp);
        cout << ans << endl;
    }
}

模拟代码如下

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
import numpy as np
from multiprocessing import Pool

class Simulater:
    def __init__(self, n, m):
        self.n = n
        self.m = m

    def test(self, T):
        np.random.seed()
        N = 0
        for _ in range(T):
            i = 0
            j = 0
            visited = [False] * self.n
            while j < self.m:
                i += 1
                x = np.random.randint(0, self.n)
                if not visited[x]:
                    visited[x] = True
                    j += 1
            N += i
        return N / T

with open("data.txt") as f:
    for line in f.readlines():
        line = line.strip().split(" ")
        n = int(line[0])
        m = int(line[1])
        s = Simulater(n, m)
        T = int(5e5)
        pool = Pool(8)
        ts = pool.map(s.test, [T] * 8)
        print(sum(ts) / 8)

测试数据

1
2
10 5
34 19

测试结果

1
2
3
4
5
6
n: 10, m: 5
计算: 6.45635
模拟: 6.45649925
n: 34, m: 19
计算: 27.1994
模拟: 27.19897675

涂砖块2

n 个格子，每次随机涂一个，可能重复涂，求涂 m 次后，涂满格子的期望个数。

思路分析1

dp[i][j] := 已经涂了 i 次，其中涂满了 j 个格子，后续的 m - i 次可以新涂满格子的期望个数

现在执行一次新的涂色，有两种可能

一种是涂到了已经涂过的，概率 j / n，那么转移到 dp[i + 1][j]

另一种是涂到了之前没涂过的，概率 (n - j) / n，那么转移到 dp[i + 1][j + 1]，同时新涂满的格子数加 1

转移方程如下

1
dp[i][j] = j / n * dp[i + 1][j] + (n - j) / n * (1 + dp[i + 1][j + 1])

基于以上分析，整理出以下期望DP算法

1
2
3
4
5
6
7
8
9
10
11
状态定义
dp[i][j] := 已经涂了 i 次，其中涂满了 j 个格子，后续的 m - i 次可以新涂满格子的期望个数

初始化
dp[m][j] = 0

答案
dp[0][0]

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

思路分析2

dp[i] := 涂 i 次涂满格子的期望个数

在涂第 i 次时，有两种情况

第一种是涂到前面涂过的，概率为 dp[i - 1] / n

第二种是涂到前面没涂过的，概率为 (n - dp[i - 1]) / n

转移方程如下

1
2
3
dp[i] = dp[i - 1] / n * dp[i - 1] + (1 - dp[i - 1] / n) * (dp[i - 1] + 1)
      = dp[i - 1] + 1 - dp[i - 1] / n
      = 1 + (n - 1) / n * dp[i - 1]

基于以上分析，整理出以下期望DP算法

1
2
3
4
5
6
7
8
9
10
11
状态定义
dp[i] := 涂 i 次涂满格子的期望个数

初始化
dp[0] = 0

答案
dp[m]

状态转移
dp[i] = 1 + (n - 1) / n * dp[i - 1]

验证答案

思路 1 的代码

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
#include <vector>
#include <fstream>
#include <iostream>

using namespace std;

double solve(int i, int j, const int n, const int m, vector<vector<double>>& dp)
{
    if(dp[i][j] > -0.5)
        return dp[i][j];
    if(i == m)
        return dp[i][j] = 0.0;

    double ans = 0.0;
    ans += j / (double)n * solve(i + 1, j, n, m, dp);
    ans += (n - j) / (double)n * (1 + solve(i + 1, j + 1, n, m, dp));
    return dp[i][j] = ans;
}

int main()
{
    ifstream fin("data.txt");
    vector<vector<double>> dp;
    int n, m;
    while((fin >> n) && (fin >> m))
    {
        cout << "n: " << n << ", m: " << m << endl;
        dp.assign(m + 1, vector<double>(m + 1, -1.0));
        double ans = solve(0, 0, n, m, dp);
        cout << ans << endl;
    }
}

思路2代码：

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
#include <vector>
#include <fstream>
#include <iostream>

using namespace std;

double solve(int i, const int n, const int m, vector<double>& dp)
{
    if(dp[i] > -0.5)
        return dp[i];
    if(i == 0)
        return dp[i] = 0.0;
    return dp[i] = 1 + (n - 1) / (double)n * solve(i - 1, n, m, dp);
}

int main()
{
    ifstream fin("data.txt");
    vector<double> dp;
    int n, m;
    while((fin >> n) && (fin >> m))
    {
        cout << "n: " << n << ", m: " << m << endl;
        dp.assign(m + 1, -1.0);
        double ans = solve(m, n, m, dp);
        cout << ans << endl;
    }
}

模拟代码

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
import numpy as np
from multiprocessing import Pool

class Simulater:
    def __init__(self, n, m):
        self.n = n
        self.m = m

    def test(self, T):
        np.random.seed()
        N = 0
        for _ in range(T):
            N += len(set(np.random.randint(0, self.n, self.m)))
        return N / T

with open("data.txt") as f:
    for line in f.readlines():
        line = line.strip().split(" ")
        n = int(line[0])
        m = int(line[1])
        s = Simulater(n, m)
        T = int(5e5)
        pool = Pool(8)
        ts = pool.map(s.test, [T] * 8)
        print(sum(ts) / 8)

测试数据

1
2
10 5
34 19

测试结果

1
2
3
4
5
6
7
8
n: 10, m: 5
思路1计算结果: 4.0951
思路2计算结果: 4.0951
模拟结果: 4.09472625
n: 34, m: 19
思路1计算结果: 14.7183
思路2计算结果: 14.7183
模拟结果: 14.719252500000001
Share