概率与期望计算

字数统计: 6.9k字 | 阅读时长: 32分

2022-01-04

摘要: 概率与期望计算若干例题

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

本文整理自《算法竞赛进阶指南》0x38，首先介绍概率与期望的基本概念和公式，然后有若干例题。

除了第一题是直接求数学期望之外，其余的题目的解法都是根据概率或期望的公式构造动态规划的转移方程解决。

基本概念

1. 概率

样本空间为 $\Omega$，对于 $\Omega$ 中的每一个随机事件 A，都存在实值函数 P(A)，满足

$\begin{align} &P(A) >= 0 \\ &P(\Omega) = 1 \\ &\forall A_{1}, A_{2}, ..., A_{k}, 有 \sum\limits_{i}P(A_{i}) = P(\bigcup\limits_{i}A_{i}) \\ \end{align}$

则称 P(A) 是随机事件 A 的概率。

2. 期望

若随机变量 X 的取值有 x1, x2, …, 一个随机事件可以表示为 $X = X_{i}$，其概率为 $P(X = X_{i}) = p_{i}$

则称 $E(X) = \sum\limits_{i}p_{i}x_{i}$ 为随机变量 X 的数学期望。

期望的性质

数学期望是新型函数

$E(aX + bY) = aE(X) + bE(y)$

例题

217. 绿豆蛙的归宿 (全期望公式, 期望DP)

给出一个有向无环的连通图，起点为 1，终点为 N，每条边都有一个长度。

数据保证从起点出发能够到达图中所有的点，图中所有的点也都能够到达终点。

绿豆蛙从起点出发，走向终点。

到达每一个顶点时，如果有 K 条离开该点的道路，绿豆蛙可以选择任意一条道路离开该点，并且走向每条路的概率为 1/K。

现在绿豆蛙想知道，从起点走到终点所经过的路径总长度的期望是多少？

输入格式：
第一行: 两个整数 N，M，代表图中有 N 个点、M 条边。
第二行到第 1+M 行: 每行 3 个整数 a,b,c，代表从 a 到 b 有一条长度为 c 的有向边。

输出格式：
输出从起点到终点路径总长度的期望值，结果四舍五入保留两位小数。

数据范围：
1<=N<=1e5,
1<=M<=2N

输入样例：
4 4
1 2 1
1 3 2
2 3 3
3 4 4
输出样例：
7.00

算法(期望DP)

用邻接表 vector<vector<Edge>> 建图，其中 Edge 中存有某条边指向的点，以及这条边的长度。具体如下

struct Edge
{
    int v;
    int l;
    Edge(int v, int l):v(v),l(l){}
};

我们要求的是从 1 到 N 的总长度期望值，首先如果已经在 N 点，那么期望长度 E(N) 当然是 0。

然后考虑在 i 点的情况，假设 i 有 K 条出边，这 K 条出边分别连接到 g[i][1].v, g[i][2].v, …, g[i][K].v 这 K 个点，分别记为 i1, i2, …, iK，那么从 i 走到这 K 个点的概率均为 1/K，如果我们知道从 i1, i2, …, iK 走到 N 的期望长度，那么由全期望公式

$\begin{align} E_{x}(x) &= E_{y}(E_{x|y}(x|y)) \\ &= \sum\limits_{y}p(y=y) \times E_{x|y=y}(x|y=y) \\ \end{align}$

我们可以写出从 E(N1), E(N2), …, E(NK) 推出 E(N) 的公式

$E(N) = \sum\limits_{i=1}\limits^{K}E(Ni)\frac{1}{K}$

下面问题变成了求 E(N1), E(N2), …, E(NK)，这些期望值的求解方法均与 E(N) 一样：找到它们的入边及其对应的前一个点，然后通过全期望公式计算，继续向前推导，直至推导到起始点。

这个过程会有很多重复计算，因此需要动态规划保存中间结果，动态规划推导的顺序是按照反图的顺序从 N 到 1。算法如下

状态定义
dp[i] := 从 i 走到 N 点的期望长度

初始值
dp[N] = 0

答案
dp[1]

状态转移
dp[i] = (1 / g[i]) * sum(g[i][j].l + dp[g[i][j].v])
其中 g 是邻接表

代码 (C++)

#include <iostream>
#include <queue>
#include <vector>
#include <iomanip>

using namespace std;

struct Edge
{
    int v;
    int l;
    Edge(int v, int l):v(v),l(l){}
};

int main()
{
    int N, M;
    cin >> N >> M;
    vector<vector<Edge>> g(N + 1);
    vector<vector<Edge>> rg(N + 1);
    for(int i = 0; i < M; ++i)
    {
        int a, b, c;
        cin >> a >> b >> c;
        rg[b].push_back(Edge(a, c));
        g[a].push_back(Edge(b, c));
    }

    vector<int> outdegree(N + 1);
    for(int i = 1; i <= N; ++i)
        outdegree[i] = g[i].size();

    queue<int> q;
    q.push(N);

    vector<double> dp(N + 1);
    while(!q.empty())
    {
        int i = q.front();
        q.pop();
        int n = g[i].size();
        for(int j = 0; j < n; ++j)
            dp[i] += g[i][j].l + dp[g[i][j].v];
        if(n > 0)
            dp[i] /= n;
        for(const Edge &son: rg[i])
        {
            --outdegree[son.v];
            if(outdegree[son.v] == 0)
                q.push(son.v);
        }
    }
    cout << std::fixed << std::setprecision(2);
    cout << dp[1] << endl;
}

218. 扑克牌 (期望DP)

Admin 生日那天，Rainbow 来找 Admin 玩扑克牌。

玩着玩着 Rainbow 觉得太没意思了，于是决定给 Admin 一个考验。

Rainbow 把一副扑克牌(54 张)随机洗开，倒扣着放成一摞。

然后 Admin 从上往下依次翻开每张牌，每翻开一张黑桃、红桃、梅花或者方块，就把它放到对应花色的堆里去。

Rainbow 想问问 Admin，得到 A 张黑桃、B 张红桃、C 张梅花、D 张方块需要翻开的牌的张数的期望值 E 是多少？

特殊地，如果翻开的牌是大王或者小王，Admin 将会把它作为某种花色的牌放入对应堆中，使得放入之后 E 的值尽可能小。

由于 Admin 和 Rainbow 还在玩扑克，所以这个程序就交给你来写了。

输入格式
输入仅由一行，包含四个用空格隔开的整数，A,B,C,D。

输出格式
输出需要翻开的牌数的期望值 E，四舍五入保留 3 位小数。

如果不可能达到输入的状态，输出 -1.000。

数据范围
0<=A,B,C,D<=15

输入样例：
1 2 3 4
输出样例：
16.393

算法(期望DP)

初始的时候，我们手上没有牌，经过从牌堆取牌，我们的目标是得到 A 张黑桃，B 张红桃，C 张梅花，D 张方块。

假设当前的状态如下

(1) 四种花色状态: 已经取到了 a 张黑桃，b 张红桃，c 张梅花，d 张方块
(2) 大王的状态为 x，有视为黑桃，红桃，梅花，方块，未翻到这 5 种，分别编码为 0,1,2,3,4
(3) 小王的状态为 y，同样有视为黑桃，红桃，梅花，方块，未翻到这 5 种，分别编码为 0,1,2,3,4
(4) 牌堆中剩余牌张数为 S = 54 - (a + b + c + d + (x != 4) + (y != 4))

在当前状态下，考虑我们取到的下一张牌，取到的花色和概率如下

黑桃: (13 - a) / S
红桃: (13 - b) / S
梅花: (13 - c) / S
方块: (13 - d) / S
大王: 1 / S
小王: 1 / S

如果我们取到的是四种花色牌，则状态的转移比较简单，以取到黑桃牌为例，新状态如下

(1) 已经取到了 a + 1 张黑桃，b 张红桃，c 张梅花，d 张方块
(2) 大小王的状态不变

红桃，梅花，方块的情况类似。

如果我们取到的是大小王，则状态转移稍微复杂一点，由于大小王可以视为任意花色，所以四种可能的花色都要考虑，这是四种不同的状态。以抽到大王，视为黑桃为例，新状态如下

(1) 四种花色状态不变
(2) 大王状态改为 0
(3) 小王状态不变

抽到小王，以及视为红桃，视为梅花，视为方块的情况类似。

四种状态中，我们要取对应的期望最小的那一种状态。

基于以上分析，我们可以抽象出期望 DP 算法，具体如下

状态定义
dp[a][b][c][d][x][y] := 已经取到了 a 张黑桃，b 张红桃，c 张梅花，d 张方块，大王状态为 x，小王状态为 y 时，达到目标的 A 黑桃，B 红桃，C 梅花，D 方块所需的步数期望

初始值
dp[a][b][c][d][x][y] = 0.0, 其中:
    a + (x == 0) + (y == 0) >= A
    b + (x == 1) + (y == 1) >= B
    c + (x == 2) + (y == 2) >= C
    d + (x == 3) + (y == 3) >= D
dp[a][b][c][d][x][y] = -1.0, 表示当前状态无法取到目标状态，具体判断条件如下
    diff = max(A - 13 - (x == 0) - (y == 0), 0)
           + max(B - 13 - (x == 1) - (y == 1), 0)
           + max(C - 13 - (x == 2) - (y == 2), 0)
           + max(D - 13 - (x == 3) - (y == 3), 0)
    diff > (x == 4) + (y == 4)

答案
dp[0][0][0][0][4][4]

状态转移
dp[a][b][c][d][x][y] = 
    1 + (13 - a) / S * dp[a+1][b][c][d][x][y]
    + (13 - b) / S * dp[a][b+1][c][d][x][y]
    + (13 - c) / S * dp[a][b][c+1][d][x][y]
    + (13 - d) / S * dp[a][b][c][d+1][x][y]
    + 1 / S * min(
        dp[a][b][c][d][0][y]
        ,dp[a][b][c][d][1][y]
        ,dp[a][b][c][d][2][y]
        ,dp[a][b][c][d][3][y]
        )
    + 1 / S * min(
        dp[a][b][c][d][x][0]
        ,dp[a][b][c][d][x][1]
        ,dp[a][b][c][d][x][2]
        ,dp[a][b][c][d][x][3]
        )

代码(Python)

由于状态的维度过高，需要五维数组，所以用 Python 写记忆化会舒服不少

import functools

class Solution:
    def __init__(self, A, B, C, D):
        self.A = A
        self.B = B
        self.C = C
        self.D = D

    @functools.lru_cache(maxsize=int(2e6))
    def solve(self, a, b, c, d, x, y):
        reached = True
        reached = reached and (a + (x == 0) + (y == 0) >= self.A)
        reached = reached and (b + (x == 1) + (y == 1) >= self.B)
        reached = reached and (c + (x == 2) + (y == 2) >= self.C)
        reached = reached and (d + (x == 3) + (y == 3) >= self.D)
        if reached:
            return 0.0
        diff = max(self.A - 13 - (x == 0) - (y == 0), 0)
        diff += max(self.B - 13 - (x == 1) - (y == 1), 0)
        diff += max(self.C - 13 - (x == 2) - (y == 2), 0)
        diff += max(self.D - 13 - (x == 3) - (y == 3), 0)
        if diff > (x == 4) + (y == 4):
            return 1e9
        S = 54 - (a + b + c + d + (x != 4) + (y != 4))
        ans = 1.0
        if 13 - a > 0:
            ans += (13 - a) / S * self.solve(a + 1, b, c, d, x, y)
        if 13 - b > 0:
            ans += (13 - b) / S * self.solve(a, b + 1, c, d, x, y)
        if 13 - c > 0:
            ans += (13 - c) / S * self.solve(a, b, c + 1, d, x, y)
        if 13 - d > 0:
            ans += (13 - d) / S * self.solve(a, b, c, d + 1, x, y)
        if x == 4:
            ans += 1 / S * min(self.solve(a, b, c, d, 0, y)
                              ,self.solve(a, b, c, d, 1, y)
                              ,self.solve(a, b, c, d, 2, y)
                              ,self.solve(a, b, c, d, 3, y)
                              )
        if y == 4:
            ans += 1 / S * min(self.solve(a, b, c, d, x, 0)
                              ,self.solve(a, b, c, d, x, 1)
                              ,self.solve(a, b, c, d, x, 2)
                              ,self.solve(a, b, c, d, x, 3)
                              )
        return ans

if __name__ == "__main__":
    tmp = input()
    A, B, C, D = tmp.split(" ")
    A = int(A)
    B = int(B)
    C = int(C)
    D = int(D)
    s = Solution(A, B, C, D)
    ans = s.solve(0, 0, 0, 0, 4, 4)
    if ans > 55:
        print("{:.3f}".format(-1.0))
    else:
        print("{:.3f}".format(round(ans, 3)))

代码(C++)

可惜 Python 代码有点慢，通过不了。还是补一个 C++ 的代码吧。

不过数组就不能用 vector 写了，因为嵌套的太多了，但是像下面这样直接声明一个 int dp[16][16][16][16][5][5] 然后 memset(dp, -1, sizeof(dp)) 初始化是可以的。

#include <iostream>
#include <cstring>
#include <iomanip>

using namespace std;

double dp[16][16][16][16][5][5];
int A, B, C, D;

double solve(int a, int b, int c, int d, int x, int y)
{
    if(dp[a][b][c][d][x][y] >= 0.0)
        return dp[a][b][c][d][x][y];
    bool reached = true;
    reached = reached && (a + (x == 0) + (y == 0)) >= A;
    reached = reached && (b + (x == 1) + (y == 1)) >= B;
    reached = reached && (c + (x == 2) + (y == 2)) >= C;
    reached = reached && (d + (x == 3) + (y == 3)) >= D;
    if(reached)
        return dp[a][b][c][d][x][y] = 0.0;
    int diff = max(A - 13 - (x == 0) - (y == 0), 0);
    diff += max(B - 13 - (x == 1) - (y == 1), 0);
    diff += max(C - 13 - (x == 2) - (y == 2), 0);
    diff += max(D - 13 - (x == 3) - (y == 3), 0);
    if(diff > (x == 4) + (y == 4))
        return dp[a][b][c][d][x][y] = 1e9;
    double ans = 1.0;
    double S = 54 - (a + b + c + d + (x != 4) + (y != 4));
    if(13 - a > 0)
        ans += (13 - a) / S * solve(a + 1, b, c, d, x, y);
    if(13 - b > 0)
        ans += (13 - b) / S * solve(a, b + 1, c, d, x, y);
    if(13 - c > 0)
        ans += (13 - c) / S * solve(a, b, c + 1, d, x, y);
    if(13 - d > 0)
        ans += (13 - d) / S * solve(a, b, c, d + 1, x, y);
    if(x == 4)
    {
        double tmp = 1e9;
        tmp = min(tmp, solve(a, b, c, d, 0, y));
        tmp = min(tmp, solve(a, b, c, d, 1, y));
        tmp = min(tmp, solve(a, b, c, d, 2, y));
        tmp = min(tmp, solve(a, b, c, d, 3, y));
        ans += 1 / S * tmp;
    }
    if(y == 4)
    {
        double tmp = 1e9;
        tmp = min(tmp, solve(a, b, c, d, x, 0));
        tmp = min(tmp, solve(a, b, c, d, x, 1));
        tmp = min(tmp, solve(a, b, c, d, x, 2));
        tmp = min(tmp, solve(a, b, c, d, x, 3));
        ans += 1 / S * tmp;
    }
    return dp[a][b][c][d][x][y] = ans;
}

int main()
{
    memset(dp, -1, sizeof(dp));
    cin >> A >> B >> C >> D;
    double ans = solve(0, 0, 0, 0, 4, 4);
    cout << std::fixed << std::setprecision(3);
    if(ans > 55)
        cout << -1.0 << endl;
    else
        cout << ans << endl;
}

232. 守卫者的挑战 (概率DP)

打开了黑魔法师 Vani 的大门，队员们在迷宫般的路上漫无目的地搜寻着关押 applepi 的监狱的所在地。

突然，眼前一道亮光闪过，“我，Nizem，是黑魔法圣殿的守卫者。如果你能通过我的挑战，那么你可以带走黑魔法圣殿的地图……”。

瞬间，队员们被传送到了一个擂台上，最初身边有一个容量为 K 的包包。

擂台赛一共有 N 项挑战，各项挑战依次进行。

第 i 项挑战有一个属性 ai，如果 ai≥0，表示这次挑战成功后可以再获得一个容量为 ai 的包包；如果 ai=−1，则表示这次挑战成功后可以得到一个大小为 1 的地图残片。

地图残片必须装在包包里才能带出擂台，包包没有必要全部装满，但是队员们必须把获得的所有的地图残片都带走（没有得到的不用考虑，只需要完成所有 N 项挑战后背包容量足够容纳地图残片即可），才能拼出完整的地图。

并且他们至少要挑战成功 L 次才能离开擂台。

队员们一筹莫展之时，善良的守卫者 Nizem 帮忙预估出了每项挑战成功的概率，其中第 i 项挑战成功的概率为 pi%。

现在，请你帮忙预测一下，队员们能够带上他们获得的地图残片离开擂台的概率。

输入格式
第一行三个整数 N,L,K。
第二行 N 个实数，第 i 个实数 pi 表示第 i 项挑战成功的概率的百分比。
第三行 N 个整数，第 i 个整数 ai 表示第 i 项挑战的属性值。

输出格式
一个实数，表示所求概率，四舍五入保留 6 位小数。

数据范围
0<=K<=2000,
0<=N<=200,
−1<=ai<=1000,
0<=L<=N,
0<=pi<=100

输入样例：
3 1 0
10 20 30
-1 -1 2

输出样例：
0.300000

算法(概率DP)

状态定义
dp[i][j][k] := 至少成功 j 次，净容量 k，从 i 到 N 的成功概率

初始化
dp[N + 1][j][k] := 1  k >= 0 and j <= 0
                   0  other

答案
dp[1][L][K]

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

代码(Python)

由于 dp[i][j][k] 中的 j, k 可能出现负数，用 Python 写记忆化方便一点

import functools

class Solution:
    def __init__(self, N, p, a):
        self.N = N
        self.p = p
        self.a = a

    @functools.lru_cache(maxsize=int(1e7))
    def solve(self, i, j, k):
        if i == self.N + 1:
            if j <= 0 and k >= 0:
                return 1.0
            else:
                return 0.0
        return self.p[i] * self.solve(i + 1, j - 1, k + self.a[i]) \
                    + (1 - self.p[i]) * self.solve(i + 1, j, k)

if __name__ == "__main__":
    tmp = input()
    N, L, K = [int(i) for i in tmp.split(" ")]
    tmp = input()
    tmp = tmp.split(" ")
    p = [0.0] * (N + 1)
    for i in range(N):
        p[i + 1] = float(tmp[i]) / 100
    tmp = input()
    tmp = tmp.split(" ")
    a = [0] * (N + 1)
    for i in range(N):
        a[i + 1] = int(tmp[i])
    s = Solution(N, p, a)
    ans = s.solve(1, L, K)
    ans = round(ans, 6)
    print("{:.6f}".format(ans))

代码(C++)

可惜 Python 还是被卡时间了，还是补一下 C++ 吧。但是要注意以下两点，这也是不愿意写 C++ 的原因。

(1) dp[i][j][k] 中，j, k 加个 n 再访问 dp 中的对应位置，避免负数
(2) 当 k > N 的时候，将 k 视为 N 再去访问 dp，但是转移的时候还是用 k 去转移

#include <iostream>
#include <iomanip>
#include <vector>

using namespace std;

vector<int> a;
vector<double> p;
int N;
vector<vector<vector<double>>> dp;

double solve(int i, int j, int k)
{
    int k_tmp = min(k, N);
    if(dp[i][j + N][k_tmp + N] >= -0.5)
        return dp[i][j + N][k_tmp + N];
    if(i == N + 1)
    {
        if(j <= 0 && k >= 0)
            return dp[i][j + N][k_tmp + N] = 1.0;
        else
            return dp[i][j + N][k_tmp + N] = 0.0;
    }
    return dp[i][j + N][k_tmp + N] = p[i] * solve(i + 1, j - 1, k + a[i])
            + (1 - p[i]) * solve(i + 1, j, k);
}

int main()
{
    int L, K;
    cin >> N >> L >> K;
    dp.assign(N + 2, vector<vector<double>>(N * 2 + 1, vector<double>(N * 2 + 1, -1.0)));
    a.assign(N + 1, 0);
    p.assign(N + 1, 0.0);
    for(int i = 0; i < N; ++i)
    {
        double tmp;
        cin >> tmp;
        p[i + 1] = tmp / 100.0;
    }
    for(int i = 0; i < N; ++i)
    {
        int tmp;
        cin >> tmp;
        a[i + 1] = tmp;
    }
    double ans = solve(1, L, K);
    cout << std::fixed << std::setprecision(6);
    cout << ans << endl;
}

233. 换教室

对于刚上大学的牛牛来说，他面临的第一个问题是如何根据实际情况申请合适的课程。

在可以选择的课程中，有 2n 节课程安排在 n 个时间段上。

在第 i (1≤i≤n) 个时间段上，两节内容相同的课程同时在不同的地点进行，其中，牛牛预先被安排在教室 ci 上课，而另一节课程在教室 di 进行。

在不提交任何申请的情况下，学生们需要按时间段的顺序依次完成所有的 n 节安排好的课程。

如果学生想更换第 i 节课程的教室，则需要提出申请。

若申请通过，学生就可以在第 i 个时间段去教室 di 上课，否则仍然在教室 ci 上课。

由于更换教室的需求太多，申请不一定能获得通过。

通过计算，牛牛发现申请更换第 i 节课程的教室时，申请被通过的概率是一个已知的实数 ki，并且对于不同课程的申请，被通过的概率是互相独立的。

学校规定，所有的申请只能在学期开始前一次性提交，并且每个人只能选择至多 m 节课程进行申请。

这意味着牛牛必须一次性决定是否申请更换每节课的教室，而不能根据某些课程的申请结果来决定其他课程是否申请。

牛牛可以申请自己最希望更换教室的 m 门课程，也可以不用完这 m 个申请的机会，甚至可以一门课程都不申请。

因为不同的课程可能会被安排在不同的教室进行，所以牛牛需要利用课间时间从一间教室赶到另一间教室。

牛牛所在的大学有 v 个教室，有 e 条道路。

每条道路连接两间教室，并且是可以双向通行的。

由于道路的长度和路况不同，通过不同的道路耗费的体力可能会有所不同。

当第 i (1≤i≤n−1) 节课结束后，牛牛就会从这节课的教室出发，选择一条耗费体力最少的路径前往下一节课的教室。

现在牛牛想知道，申请哪几门课程可以使他因在教室间移动耗费的体力值的总和的期望值最小，请你帮他求出这个最小值。

输入格式
第一行四个整数 n,m,v,e，n 表示这个学期内的时间段的数量; m 表示牛牛最多可以申请更换多少节课程的教室; v 表示牛牛学校里教室的数量; e 表示牛牛的学校里道路的数量。
第二行 n 个正整数，第 i 个正整数表示 ci，即第 i 个时间段牛牛被安排上课的教室；保证 1<=ci<=v。
第三行 n 个正整数，第 i 个正整数表示 di，即第 i 个时间段另一间上同样课程的教室；保证 1<=di<=v。
第四行 n 个实数，第 i 个实数表示 ki，即牛牛申请在第 i 个时间段更换教室获得通过的概率；保证 0<=ki<=1。
接下来 e 行，每行三个正整数aj,bj,wj,表示有一条双向道路连接教室 aj,bj ,通过这条道路需要耗费的体力值是 wj ;保证1<=aj,bj<=v,1<=wj<=100。
保证 1<=n<=2000,0<=m<=2000,1<=v<=300,0<=e<=90000。
保证通过学校里的道路，从任何一间教室出发，都能到达其他所有的教室。
保证输入的实数最多包含 3 位小数。

输出格式
输出一行，包含一个实数，四舎五入精确到小数点后恰好 2 位，表示答案。
你的输出必须和标准输出完全一样才算正确。
测试数据保证四舎五入后的答案和准确答案的差的绝对值不大于 4*10−3 。 (如果你不知道什么是浮点误差，这段话可以理解为: 对于大多数的算法， 你可以正常地使用浮点数类型而不用对它进行特殊的处理)

输入样例：
3 2 3 3
2 1 2
1 2 1
0.8 0.2 0.5 
1 2 5
1 3 3
2 3 1

输出样例：
2.80

算法(期望DP)

本题主要分为两步，第一步针对 v 个节点 e 条边，用 Floyd 算法预处理出所有点之间的最短路径，dist[i][j] 表示 i 到 j 的最短路径。

然后依次考虑 N 个时间段，假设当前在第 i 个时间段，在教室 c[i] 或 d[i]，可以申请的额度还剩 j 个。

考虑第 i + 1 个时间段是否换，如果 j = 0 没得考虑，只能不换，如果 j > 0，就可以换或者不换。

如果不换，那么总期望距离加上从当前点到 c[i + 1] 的距离，然后考虑第 i + 1 个时间段。

如果换，并且以概率 p[i + 1] 成功了，那么总期望距离加上从当前点到 d[i + 1] 的距离，申请额度减 1，然后考虑第 i + 1 个时间段。

同时也有 (1 - p[i + 1]) 概率失败，那么总期望加上从当前点到 c[i + 1] 的距离，申请额度减 1，然后考虑第 i + 1 个时间段。

基于以上分析，总结一下，可以写出以下期望 DP 算法

注意 DP 状态的推导方向与分析过程是相反的，主要是为了实现方便。

状态定义
dp1[i][j] := 从第 1 个时间段到第 i 个时间段，申请了 j 次，第 i 个时间段没申请的期望距离
dp2[i][j] := 从第 1 个时间段到第 i 个时间段，申请了 j 次，第 i 个时间段申请了的期望距离

初始化
dp1[1][0] = 0
dp2[1][1] = 0

答案
min(dp1[N][j], dp2[N][j])

状态转移
dp1[i][j] = min(dist[c[i - 1]][c[i]] + dp1[i - 1][j]
               ,p[i - 1] * (dist[d[i - 1]][c[i]] + dp2[i - 1][j])
                   + (1 - p[i - 1]) * (dist[c[i - 1]][c[i]] + dp2[i - 1][j])
               )
          = min(dist[c[i - 1]][c[i]] + dp1[i - 1][j]
               ,dp2[i - 1][j] + dist[d[i - 1]][c[i]] * p[i - 1]
                              + dist[c[i - 1]][c[i]] * (1 - p[i - 1])
               )
dp2[i][j] = min(p[i] * (dist[c[i - 1]][d[i]] + dp1[i - 1][j - 1])
                   + (1 - p[i]) * (dist[c[i - 1]][c[i]] + dp1[i - 1][j - 1])
               ,p[i] * (p[i - 1] * (dist[d[i - 1]][d[i]] + dp2[i - 1][j - 1])
                        + (1 - p[i - 1]) * (dist[c[i - 1]][d[i]] + dp2[i - 1][j - 1])
                ) + (1 - p[i]) * (p[i - 1] * (dist[d[i - 1]][c[i]] + dp2[i - 1][j - 1])
                        + (1 - p[i - 1]) * (dist[c[i - 1]][c[i]] + dp2[i - 1][j - 1])
                )
               )
          = min(dp1[i - 1][j - 1] + dist[c[i - 1]][d[i]] * p[i]
                                  + dist[c[i - 1]][c[i]] * (1 - p[i])
               ,dp2[i - 1][j - 1] + dist[d[i - 1]][d[i]] * p[i] * (p[i - 1]
                                  + dist[c[i - 1]][d[i]] * p[i] * (1 - p[i - 1])
                                  + dist[d[i - 1]][c[i]] * (1 - p[i]) * p[i - 1]
                                  + dist[c[i - 1]][c[i]] * (1 - p[i]) * (1 - p[i - 1])
               )

代码(C++)

#include <iostream>
#include <iomanip>
#include <vector>
#include <cstring>

using namespace std;

const double INF = 1e9;
const int M = 310;

int dist[M][M];
vector<vector<double>> dp1, dp2;
vector<int> c, d;
vector<double> p;
int n;

int main()
{
    int m, v, e;
    cin >> n >> m >> v >> e;
    c.assign(n + 1, -1);
    for(int i = 0; i < n; ++i)
        cin >> c[i + 1];
    d.assign(n + 1, -1);
    for(int i = 0; i < n; ++i)
        cin >> d[i + 1];
    p.assign(n + 1, -1);
    for(int i = 0; i < n; ++i)
        cin >> p[i + 1];
    memset(dist, -1, sizeof(dist));
    for(int i = 1; i <= v; ++i)
        dist[i][i] = 0;
    for(int i = 0; i < e; ++i)
    {
        int a, b, w;
        cin >> a >> b >> w;
        if(dist[a][b] == -1)
            dist[a][b] = dist[b][a] = w;
        else
            dist[a][b] = dist[b][a] = min(dist[a][b], w);
    }
    // 遍历所有节点，其中 k 是用于松弛的点
    for(int k = 1; k <= v; ++k)
        for(int i = 1; i <= v; ++i)
            for(int j = 1; j <= v; ++j)
                // 使用 k 松弛 i > j 的最短路径
                if(dist[i][k] != -1 && dist[k][j] != -1)
                {
                    if(dist[i][j] != -1)
                        dist[i][j] = min(dist[i][j], dist[i][k] + dist[k][j]);
                    else
                        dist[i][j] = dist[i][k] + dist[k][j];
                }
    dp1.assign(n + 1, vector<double>(n + 1, INF));
    dp2.assign(n + 1, vector<double>(n + 1, INF));
    dp1[1][0] = 0;
    dp2[1][1] = 0;
    for(int i = 2; i <= n; ++i)
    {
        dp1[i][0] = dist[c[i - 1]][c[i]] + dp1[i - 1][0];
        for(int j = 1; j <= min(i, m); ++j)
        {
            dp1[i][j] = min(dist[c[i - 1]][c[i]] + dp1[i - 1][j]
                           ,p[i - 1] * (dist[d[i - 1]][c[i]] + dp2[i - 1][j])
                               + (1 - p[i - 1]) * (dist[c[i - 1]][c[i]] + dp2[i - 1][j])
                           );
            dp2[i][j] = min(dp1[i - 1][j - 1] + dist[c[i - 1]][d[i]] * p[i]
                                              + dist[c[i - 1]][c[i]] * (1 - p[i])
                           ,dp2[i - 1][j - 1] + dist[d[i - 1]][d[i]] * p[i] * p[i - 1]
                                              + dist[c[i - 1]][d[i]] * p[i] * (1 - p[i - 1])
                                              + dist[d[i - 1]][c[i]] * (1 - p[i]) * p[i - 1]
                                              + dist[c[i - 1]][c[i]] * (1 - p[i]) * (1 - p[i - 1])
                           );
        }
    }
    double ans = INF;
    for(int j = 0; j <= m; ++j)
        ans = min(ans, min(dp1[n][j], dp2[n][j]));

    cout << std::fixed << std::setprecision(2);
    cout << ans << endl;
}

潮汐朝夕的生活实验室 \Doge 陪伴一个算法工程师的职业生涯

数概率

概率与期望计算

基本概念

1. 概率

2. 期望

期望的性质

例题

216. Rainbow的信号 (求数学期望)

217. 绿豆蛙的归宿 (全期望公式, 期望DP)

算法(期望DP)

代码 (C++)

218. 扑克牌 (期望DP)

算法(期望DP)

代码(Python)

代码(C++)

232. 守卫者的挑战 (概率DP)

算法(概率DP)

代码(Python)

代码(C++)

233. 换教室

算法(期望DP)

代码(C++)