動態規劃之經典數學期望和概率DP

起因：在一場訓練賽上。有這麼一題沒做出來。

題目連結：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程式碼：

// Author : RioTian
// Time : 20/12/09
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e2 + 10;
double dp[N][N][N];
int main() {
    // freopen("in.txt", "r", stdin);
    ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    int a, b, c;
    cin >> a >> b >> c;
    for (int i = 99; i >= a; i--)
        for (int j = 99; j >= b; j--)
            for (int k = 99; k >= c; k--) {
                // 令 t = x + y + z，減少程式碼量
                double t = i + j + k;
                dp[i][j][k] = i / t * (dp[i + 1][j][k] + 1) +
                              j / t * (dp[i][j + 1][k] + 1) +
                              k / t * (dp[i][j][k + 1] + 1);
            }
    // 關於 C++ 的輸出控制可以在我的以前部落格找到
    cout << fixed << setprecision(9) << dp[a][b][c] << endl;
}

時間和空間複雜度：

\(O(n^3)\)

當然在訓練賽的題解我也提到也可以用蒙特卡洛方法模擬在\(O(n)\)解決，這裡就不再解釋了，有興趣的話可以去題解報告的那篇部落格查閱

做完上面那道概率DP題，算是入了門，但僅僅入門不夠的需要加強，所以開始學習各位大神的部落格。

下面先說一下9974dalao的總結：

很多概率題總逃不開用dp轉移。

期望題總是倒著推過來的，概率是正著推的，多做題就會理解其中的原因

有些期望題要用到有關 概率 或 期望的常見\(\color{red}公式或思想\)

​ 遇到dp轉移方程(組)中有環的，多半逃不出\(\color{red}高斯消元\)（手動 和 寫程式碼 兩種）

​ 這套題中還有道樹上的dp轉移，還用dfs對方程迭代解方程， 真是大開眼界了

​ 當然還有與各種演算法結合的題，關鍵還是要\(\color{red}學會分析\)

​ 當公式或計算時有除法時， 特別要注意\(\color{red}分母是否為零\)

下面幾個題型來自Oi wiki 和 kuangbin 的總結

DP 求概率

這類題目採用順推，也就是從初始狀態推向結果。同一般的 DP 類似的，難點依然是對狀態轉移方程的刻畫，只是這類題目經過了概率論知識的包裝。

"例題 Codeforces 148 D Bag of mice"

​ 題目大意：袋子裡有 w 只白鼠和 b 只黑鼠，公主和龍輪流從袋子裡抓老鼠。誰先抓到白色老鼠誰就贏，如果袋子裡沒有老鼠了並且沒有誰抓到白色老鼠，那麼算龍贏。公主每次抓一隻老鼠，龍每次抓完一隻老鼠之後會有一隻老鼠跑出來。每次抓的老鼠和跑出來的老鼠都是隨機的。公主先抓。問公主贏的概率。

設 \(f_{i,j}\) 為輪到公主時袋子裡有 \(i\) 只白鼠， \(j\) 只黑鼠，公主贏的概率。初始化邊界， \(f_{0,j}=0\) 因為沒有白鼠了算龍贏， \(f_{i,0}=1\) 因為抓一隻就是白鼠，公主贏。
考慮 \(f_{i,j}\) 的轉移：

• 公主抓到一隻白鼠，公主贏了。概率為 \(\frac{i}{i+j}\) ；
• 公主抓到一隻黑鼠，龍抓到一隻白鼠，龍贏了。概率為 \(\frac{j}{i+j}\cdot \frac{i}{i+j-1}\) ；
• 公主抓到一隻黑鼠，龍抓到一隻黑鼠，跑出來一隻黑鼠，轉移到 \(f_{i,j-3}\) 。概率為 \(\frac{j}{i+j}\cdot\frac{j-1}{i+j-1}\cdot\frac{j-2}{i+j-2}\) ；
• 公主抓到一隻黑鼠，龍抓到一隻黑鼠，跑出來一隻白鼠，轉移到 \(f_{i-1,j-2}\) 。概率為 \(\frac{j}{i+j}\cdot\frac{j-1}{i+j-1}\cdot\frac{i}{i+j-2}\) ；

考慮公主贏的概率，第二種情況不參與計算。並且要保證後兩種情況合法，所以還要判斷 \(i,j\) 的大小，滿足第三種情況至少要有 3 只黑鼠，滿足第四種情況要有 1 只白鼠和 2 只黑鼠。

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;
}

相關習題

• POJ3071 Football
• CodeForces 768 D Jon and Orbs

DP 求期望

" 例題 POJ2096 Collecting Bugs"

​ 題目大意：一個軟體有 s 個子系統，會產生 n 種 bug。某人一天發現一個 bug，這個 bug 屬於某種 bug 分類，也屬於某個子系統。每個 bug 屬於某個子系統的概率是 \(\frac{1}{s}\) ，屬於某種 bug 分類的概率是 \(\frac{1}{n}\) 。求發現 n 種 bug，且 s 個子系統都找到 bug 的期望天數。

令 \(f_{i,j}\) 為已經找到 \(i\) 種 bug 分類， \(j\) 個子系統的 bug，達到目標狀態的期望天數。這裡的目標狀態是找到 \(n\) 種 bug 分類， \(j\) 個子系統的 bug。那麼就有 \(f_{n,s}=0\) ，因為已經達到了目標狀態，不需要用更多的天數去發現 bug 了，於是就以目標狀態為起點開始遞推，答案是 \(f_{0,0}\) 。

考慮 \(f_{i,j}\) 的狀態轉移：

• \(f_{i,j}\) ，發現一個 bug 屬於已經發現的 \(i\) 種 bug 分類， \(j\) 個子系統，概率為 \(p_1=\frac{i}{n}\cdot\frac{j}{s}\)
• \(f_{i,j+1}\) ，發現一個 bug 屬於已經發現的 \(i\) 種 bug 分類，不屬於已經發現的子系統，概率為 \(p_2=\frac{i}{n}\cdot(1-\frac{j}{s})\)
• \(f_{i+1,j}\) ，發現一個 bug 不屬於已經發現 bug 分類，屬於 \(j\) 個子系統，概率為 \(p_3=(1-\frac{i}{n})\cdot\frac{j}{s}\)
• \(f_{i+1,j+1}\) ，發現一個 bug 不屬於已經發現 bug 分類，不屬於已經發現的子系統，概率為 \(p_4=(1-\frac{i}{n})\cdot(1-\frac{j}{s})\)

再根據期望的線性性質，就可以得到狀態轉移方程：

\[\begin{aligned} f_{i,j} &= p_1\cdot f_{i,j}+p_2\cdot f_{i,j+1}+p_3\cdot f_{i+1,j}+p_4\cdot f_{i+1,j+1} + 1\\ &=(p_2\cdot f_{i,j+1}+p_3\cdot f_{i+1,j}+p_4\cdot f_{i+1,j+1}+1)/(1-p_1) \end{aligned} \]

AC 程式碼

// Author : RioTian
// Time : 20/12/10
#include <cstdio>
#include <cstring>
#include <iostream>
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」換教室"

​ 題目大意：牛牛要上 \(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程式碼

#include <bits/stdc++.h>
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
• HDU4035 Maze
• 「SCOI2008」獎勵關

有後效性 DP

例題 "CodeForces 24 D Broken robot"

​ 題目大意：給出一個 \(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 列的矩陣，然後進行 高斯消元 即可解出答案。

AC程式碼：這個有點繞，程式碼博主沒有自己寫，copy下dalao的程式碼了（侵權刪）。

#include <bits/stdc++.h>
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
• 「HNOI2013」遊走

參考文獻

kuangbin 概率 DP 總結

ACM裡的期望和概率問題 從入門到精通

論文學習：

《資訊學競賽中概率問題求解初探》

《淺析競賽中一類數學期望問題的解決方法》

《有關概率和期望問題的研究 》

關於一些練習題：

概率dp總結