動態規劃之經典數學期望和概率DP
阿新 • • 發佈：2020-12-10
起因：在一場訓練賽上。有這麼一題沒做出來。
題目連結：http://acm.hdu.edu.cn/showproblem.php?pid=6829
題目大意：有三個人，他們分別有$X，Y，Z$塊錢（$1<=X,Y,Z<=1e6$)，錢數最多的（如果不止一個那麼隨機等概率的選一個）隨機等可能的選另一個人送他一塊錢。直到三個人錢數相同為止。輸出送錢輪數的期望，如果根本停不下來，輸出-1。
根據題目的意思，其實就是每次向包裡隨機加入一枚錢幣，直到包裡某種錢幣數量達到 100。本題的核心是如何計算期望。本題屬於標準的動態規劃求期望問題。直接套用模板即可。
> 一道”簡單“概率DP題，沒怎麼了解概率DP導致做不出
當理解基礎的知識以後發現的確比較簡單
**DP 陣列定義**
定義 $DP[i][j][k]$，表示有 i 枚金幣， j 枚銀幣， k 枚銅幣的期望。
**初值**
所有的期望都為零。
**遞推方法**
使用逆推。
**狀態轉移方程：**
$$
s = X + Y + Z\\
dp(i,j,k) = \frac{i}{s}*(dp(i + 1,j,k) + 1) + \frac{j}{s}*dp(i,j + 1,k) + 1) \\+ \frac{j}{s}*dp(i,j,k + 1) + 1))
$$
**AC程式碼：**
```cpp
// Author : RioTian
// Time : 20/12/09
#include 
#include
#include
using namespace std;
const int N = 1e3 + 10;
int n, s;
double dp[N][N];
int main() {
while (cin >> n >> s) {
dp[n][s] = 0;
for (int i = n; i >= 0; i--)
for (int j = s; j >= 0; j--) {
if (i == n && j == s) //跳過初始值
continue;
double p1, p2, p3, p4;
p1 = (double)i * j / (n * s); // dp[i][j];
p2 = (double)(n - i) * j / (n * s); // dp[i+1][j];
p3 = (double)i * (s - j) / (n * s); // dp[i][j+1];
p4 = (double)(n - i) * (s - j) / (n * s);
dp[i][j] = (1 + p2 * dp[i + 1][j] + p3 * dp[i][j + 1] +
p4 * dp[i + 1][j + 1]) /
(1 - p1);
}
printf("%.4f\n", dp[0][0]);
}
}
```
### "例題 [「NOIP2016」換教室](http://uoj.ac/problem/262)"
題目大意：牛牛要上 $n$ 個時間段的課，第 $i$ 個時間段在 $c_i$ 號教室，可以申請換到 $d_i$ 號教室，申請成功的概率為 $p_i$ ，至多可以申請 $m$ 節課進行交換。第 $i$ 個時間段的課上完後要走到第 $i+1$ 個時間段的教室，給出一張圖 $v$ 個教室 $e$ 條路，移動會消耗體力，申請哪幾門課程可以使他因在教室間移動耗費的體力值的總和的期望值最小，也就是求出最小的期望路程和。
對於這個無向連通圖，先用 Floyd 求出最短路，為後續的狀態轉移帶來便利。以移動一步為一個階段（從第 $i$ 個時間段到達第 $i+1$ 個時間段就是移動了一步），那麼每一步就有 $p_i$ 的概率到 $d_i$ ，不過在所有的 $d_i$ 中只能選 $m$ 個，有 $1-p_i$ 的概率到 $c_i$ ，求出在 $n$ 個階段走完後的最小期望路程和。
定義 $f_{i,j,0/1}$ 為在第 $i$ 個時間段，連同這一個時間段已經用了 $j$ 次換教室的機會，在這個時間段換（1）或者不換（0）教室的最小期望路程和，那麼答案就是 $max \{f_{n,i,0},f_{n,i,1}\} ,i\in[0,m]$ 。注意邊界 $f_{1,0,0}=f_{1,1,1}=0$ 。
考慮 $f_{i,j,0/1}$ 的狀態轉移：
- 如果這一階段不換，即 $f_{i,j,0}$ 。可能是由上一次不換的狀態轉移來的，那麼就是 $f_{i-1,j,0}+w_{c_{i-1},c_{i}}$ , 也有可能是由上一次交換的狀態轉移來的，這裡結合條件概率和全概率的知識分析可以得到 $f_{i-1,j,1}+w_{d_{i-1},c_{i}}\cdot p_{i-1}+w_{c_{i-1},c_{i}}\cdot (1-p_{i-1})$ ，狀態轉移方程就有
$$
\begin{aligned}
f_{i,j,0}=min(f_{i-1,j,0}+w_{c_{i-1},c_{i}},f_{i-1,j,1}+w_{d_{i-1},c_{i}}\cdot p_{i-1}+w_{c_{i-1},c_{i}}\cdot (1-p_{i-1}))
\end{aligned}
$$
- 如果這一階段交換，即 $f_{i,j,1}$ 。類似地，可能由上一次不換的狀態轉移來，也可能由上一次交換的狀態轉移來。那麼遇到不換的就乘上 $(1-p_i)$ ，遇到交換的就乘上 $p_i$ ，將所有會出現的情況都列舉一遍出進行計算就好了。這裡不再贅述各種轉移情況，相信通過上一種階段例子，這裡的狀態轉移應該能夠很容易寫出來。
AC程式碼
```cpp
#include
using namespace std;
const int maxn = 2010;
int n, m, v, e;
int f[maxn][maxn], c[maxn], d[maxn];
double dp[maxn][maxn][2], p[maxn];
int main() {
scanf("%d %d %d %d", &n, &m, &v, &e);
for (int i = 1; i <= n; i++)
scanf("%d", &c[i]);
for (int i = 1; i <= n; i++)
scanf("%d", &d[i]);
for (int i = 1; i <= n; i++)
scanf("%lf", &p[i]);
for (int i = 1; i <= v; i++)
for (int j = 1; j < i; j++)
f[i][j] = f[j][i] = 1e9;
int u, V, w;
for (int i = 1; i <= e; i++) {
scanf("%d %d %d", &u, &V, &w);
f[u][V] = f[V][u] = min(w, f[u][V]);
}
for (int k = 1; k <= v; k++)
for (int i = 1; i <= v; i++)
for (int j = 1; j < i; j++)
if (f[i][k] + f[k][j] < f[i][j])
f[i][j] = f[j][i] = f[i][k] + f[k][j];
for (int i = 1; i <= n; i++)
for (int j = 0; j <= m; j++)
dp[i][j][0] = dp[i][j][1] = 1e9;
dp[1][0][0] = dp[1][1][1] = 0;
for (int i = 2; i <= n; i++)
for (int j = 0; j <= min(i, m); j++) {
dp[i][j][0] =
min(dp[i - 1][j][0] + f[c[i - 1]][c[i]],
dp[i - 1][j][1] + f[c[i - 1]][c[i]] * (1 - p[i - 1]) +
f[d[i - 1]][c[i]] * p[i - 1]);
if (j != 0) {
dp[i][j][1] =
min(dp[i - 1][j - 1][0] + f[c[i - 1]][d[i]] * p[i] +
f[c[i - 1]][c[i]] * (1 - p[i]),
dp[i - 1][j - 1][1] +
f[c[i - 1]][c[i]] * (1 - p[i - 1]) * (1 - p[i]) +
f[c[i - 1]][d[i]] * (1 - p[i - 1]) * p[i] +
f[d[i - 1]][c[i]] * (1 - p[i]) * p[i - 1] +
f[d[i - 1]][d[i]] * p[i - 1] * p[i]);
}
}
double ans = 1e9;
for (int i = 0; i <= m; i++)
ans = min(dp[n][i][0], min(dp[n][i][1], ans));
printf("%.2lf", ans);
return 0;
}
```
比較這兩個問題可以發現，DP 求期望題目在對具體是求一個值或是最優化問題上會對方程得到轉移方式有一些影響，但無論是 DP 求概率還是 DP 求期望，總是離不開概率知識和列出、化簡計算公式的步驟，在寫狀態轉移方程時需要思考的細節也類似。
### 習題
- [HDU3853 LOOPS](http://acm.hdu.edu.cn/showproblem.php?pid=3853)
- [HDU4035 Maze](http://acm.hdu.edu.cn/showproblem.php?pid=4035)
- [「SCOI2008」獎勵關](https://www.luogu.com.cn/problem/P2473)
## 有後效性 DP
### 例題 "[CodeForces 24 D Broken robot](https://codeforces.com/problemset/problem/24/D)"
題目大意：給出一個 $n*m$ 的矩陣區域，一個機器人初始在第 $x$ 行第 $y$ 列，每一步機器人會等概率地選擇停在原地，左移一步，右移一步，下移一步，如果機器人在邊界則不會往區域外移動，問機器人到達最後一行的期望步數。
在 $m=1$ 時每次有 $\frac{1}{2}$ 的概率不動，有 $\frac{1}{2}$ 的概率向下移動一格，答案為 $2\cdot (n-x)$ 。
設 $f_{i,j}$ 為機器人機器人從第 i 行第 j 列出發到達第 $n$ 行的期望步數，最終狀態為 $f_{n,j}=0$ 。
由於機器人會等概率地選擇停在原地，左移一步，右移一步，下移一步，考慮 $f_{i,j}$ 的狀態轉移：
- $f_{i,1}=\frac{1}{3}\cdot(f_{i+1,1}+f_{i,2}+f_{i,1})+1$
- $f_{i,j}=\frac{1}{4}\cdot(f_{i,j}+f_{i,j-1}+f_{i,j+1}+f_{i+1,j})+1$
- $f_{i,m}=\frac{1}{3}\cdot(f_{i,m}+f_{i,m-1}+f_{i+1,m})+1$
在行之間由於只能向下移動，是滿足無後效性的。在列之間可以左右移動，在移動過程中可能產生環，不滿足無後效性。
將方程變換後可以得到：
- $2f_{i,1}-f_{i,2}=3+f_{i+1,1}$
- $3f_{i,j}-f_{i,j-1}-f_{i,j+1}=4+f_{i+1,j}$
- $2f_{i,m}-f_{i,m-1}=3+f_{i+1,m}$
由於是逆序的遞推，所以每一個 $f_{i+1,j}$ 是已知的。
由於有 $m$ 列，所以右邊相當於是一個 $m$ 行的列向量，那麼左邊就是 $m$ 行 $m$ 列的矩陣。使用增廣矩陣，就變成了 m 行 m+1 列的矩陣，然後進行 [高斯消元](../math/gauss.md) 即可解出答案。
AC程式碼：這個有點繞，程式碼博主沒有自己寫，copy下dalao的程式碼了（侵權刪）。
```cpp
#include
using namespace std;
const int maxn = 1e3 + 10;
double a[maxn][maxn], f[maxn];
int n, m;
void solve(int x) {
memset(a, 0, sizeof a);
for (int i = 1; i <= m; i++) {
if (i == 1) {
a[i][i] = 2;
a[i][i + 1] = -1;
a[i][m + 1] = 3 + f[i];
continue;
} else if (i == m) {
a[i][i] = 2;
a[i][i - 1] = -1;
a[i][m + 1] = 3 + f[i];
continue;
}
a[i][i] = 3;
a[i][i + 1] = -1;
a[i][i - 1] = -1;
a[i][m + 1] = 4 + f[i];
}
for (int i = 1; i < m; i++) {
double p = a[i + 1][i] / a[i][i];
a[i + 1][i] = 0;
a[i + 1][i + 1] -= a[i][i + 1] * p;
a[i + 1][m + 1] -= a[i][m + 1] * p;
}
f[m] = a[m][m + 1] / a[m][m];
for (int i = m - 1; i >= 1; i--)
f[i] = (a[i][m + 1] - f[i + 1] * a[i][i + 1]) / a[i][i];
}
int main() {
scanf("%d %d", &n, &m);
int st, ed;
scanf("%d %d", &st, &ed);
if (m == 1) {
printf("%.10f\n", 2.0 * (n - st));
return 0;
}
for (int i = n - 1; i >= st; i--) {
solve(i);
}
printf("%.10f\n", f[ed]);
return 0;
}
```
### 習題
- [HDU Time Travel](http://acm.hdu.edu.cn/showproblem.php?pid=4418)
- [「HNOI2013」遊走](https://loj.ac/problem/2383)
## 參考文獻
[kuangbin 概率 DP 總結](https://www.cnblogs.com/kuangbin/archive/2012/10/02/2710606.html)
[ACM裡的期望和概率問題 從入門到精通](https://www.cnblogs.com/idyllic/p/13460122.html)
**論文學習：**
**[《資訊學競賽中概率問題求解初探》](http://wenku.baidu.com/view/1c41152de2bd960590c677a8.html)**
[**《淺析競賽中一類數學期望問題的解決方法》**](http://wenku.baidu.com/view/90adb02acfc789eb172dc8a8.html)
[**《有關概率和期望問題的研究 》**](http://wenku.baidu.com/view/56147518a8114431b90dd81e.html)
關於一些練習題：
**[概率dp總結](https://www.zybuluo.com/01010101/note/13
AC 程式碼(建議獨立思考以後再閱讀)
// Author : RioTian
// Time : 20/12/10
#include
using namespace std;
typedef long long ll;
const int N = 1e3 + 10;
int w, b;
double dp[N][N];
int main() {
// freopen("in.txt", "r", stdin);
ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
cin >> w >> b;
memset(dp, 0.0, sizeof dp);
for (int i = 1; i <= w; ++i)
dp[i][0] = 1;
for (int j = 1; j <= b; ++j)
dp[0][j] = 0;
for (int i = 1; i <= w; ++i)
for (int j = 1; j <= b; ++j) {
dp[i][j] += (double)i / (i + j);
if (j >= 3)
dp[i][j] += (double)j / (i + j) * (j - 1) / (i + j - 1) *
(j - 2) / (i + j - 2) * dp[i][j - 3];
if (i >= 1 && j >= 2)
dp[i][j] += (double)j / (i + j) * (j - 1) / (i + j - 1) * i /
(i + j - 2) * dp[i - 1][j - 2];
}
cout << fixed << setprecision(9) << dp[w][b] << endl;
}