【算法學習筆記】機率與指望DP

本文學習自 Sengxian 學長的博客

以前也在CF上寫了一些機率DP的題並作過總結

建議閱讀完本文再去接着閱讀這篇文章：Here

前言

單純只用到機率的題並非不少，從現有的 OI/ACM 比賽中來看，大多數題目須要機率與指望結合起來（指望就是用機率定義的），因此本文主要講述指望 DP。

指望 DP 有一些固定的方法，這裏分多種方法來說述。

---

講解

例一

#3036. 綠豆蛙的歸宿

題意：

給定一個起點爲 \(1\)，終點爲 \(n\) 的有向無環圖。到達每個頂點時，若是有 \(K\) 條離開該點的道路，能夠選擇任意一條道路離開該點，而且走向每條路的機率爲 \(\frac 1 K\)。問你從 \(1\) 出發走到 \(n\) 的路徑指望總長度是多少。

方法一：直接定義指望狀態

這道題的終點很明確，那就是走到 \(n\) 即中止。對於指望 DP，咱們通常採用逆序的方式來定義狀態，即考慮從當前狀態到達終點的指望代價。由於在大多數狀況下，終點不惟一，而起點是惟一的。
咱們定義 \(dp(i)\)爲從 \(i\) 出發走到終點 \(n\) 的路徑指望總長度，根據全指望公式，獲得（設​ \(G_i\)爲從 \(i\) 的邊的集合）：

\[dp(i) = \sum\limits_{e\in G_i}\frac{dp(e_{to}) + e_{const}}{|G_i|} \]

由於這是一個有向無環圖，每一個點須要其能到達的點的狀態，故咱們採用拓撲序的逆序進行計算便可。

【AC Code】

const int N = 100000, M = 2 * N;
int n, m;
struct node { int v, w;};
vector<node>e[N];
int d[N];    // 出度
double f[N]; // dp

double dfs(int u) {
    if (f[u] >= 0)return f[u];
    f[u] = 0;
    for (auto [v, w] : e[u]) { // tuple 需開啓 C++17
        f[u] += (w + dfs(v)) / d[u];
    }
    return f[u];
}

int main() {
    cin.tie(nullptr)->sync_with_stdio(false);
    memset(f, -1, sizeof(f));
    cin >> n >> m;
    for (int i = 0; i < m; ++i) {
        int u, v, w;
        cin >> u >> v >> w;
        e[u].push_back(node{v, w});
        d[u]++;//出度++
    }
    cout << fixed << setprecision(2) << dfs(1);
}

方法二：利用指望的線性性質

根據指望的線性性質，\(E[aX +bY]=aE[X] + bE[Y]\)。因此另一種求指望的方式是分別求出每一種代價產生的指望貢獻，而後相加獲得答案。在本題中，路徑指望總長度等於每條邊產生的指望貢獻之和。而每條邊產生又等於通過這條邊的指望次數乘這條邊的代價。因此，咱們只須要算每條邊的指望通過次數便可。

邊 \((u,v,w)\) 的指望通過次數是很差直接算的，但若是咱們能算得點 \(u\) 的指望通過次數爲 \(dp(u,v)\)，那麼邊 \((u,v,w)\) 的指望通過次數是 \(dp(u)*\frac1{|G_u|}\) ，對答案的貢獻就是 \(w*dp(u)*\frac1{|G_u|}\)

如何計算點 \(u\) 的指望通過次數 \(dp(u)\)呢？咱們依然考慮 DP 的方式，首先有 \(dp(u) = 1\)，轉移採起刷表的方式：

\[dp(e_{to})\leftarrow dp(u)*\frac1{|G_u|},e\in G_u \]

在用邊 \(e\) 刷表的同時，邊 \(e\) 的貢獻就能夠計算了，十分簡潔。由於這種方法計算答案十分的便捷，並且適用範圍廣，因此這種『利用指望的線性性質，單獨計算貢獻的方法』是咱們計算指望首選的方法。

【AC Code】這裏貼 Sengxian 學長的代碼

typedef long long ll;
inline int readInt() {
    static int n, ch;
    n = 0, ch = getchar();
    while (!isdigit(ch)) ch = getchar();
    while (isdigit(ch)) n = n * 10 + ch - '0', ch = getchar();
    return n;
}

const int MAX_N = 100000 + 3, MAX_M = MAX_N * 2;
struct edge {
    edge *next;
    int to, cost;
    edge(edge *next = NULL, int to = 0, int cost = 0): next(next), to(to), cost(cost) {}
} pool[MAX_M], *pit = pool, *first[MAX_N];
int n, m, deg[MAX_N], outDeg[MAX_N];
double f[MAX_N];

void solve() {
    static int q[MAX_N];
    int l = 0, r = 0;
    q[r++] = 0;

    double ans = 0;
    f[0] = 1.0;
    while (r - l >= 1) {
        int u = q[l++];
        for (edge *e = first[u]; e; e = e->next) {
            f[e->to] += f[u] / outDeg[u];
            ans += f[u] * e->cost / outDeg[u];
            if (--deg[e->to] == 0) q[r++] = e->to;
        }
    }

    printf("%.2f\n", ans);
}

int main() {
    n = readInt(), m = readInt();
    for (int i = 0, u, v, w; i < m; ++i) {
        u = readInt() - 1, v = readInt() - 1, w = readInt();
        first[u] = new (pit++) edge(first[u], v, w);
        deg[v]++, outDeg[u]++;
    }
    solve();
}

例二

接着咱們考慮 Codeforces 518D 這道題，以便體會方法二的好處。

題意：有 \(n\) 我的排成一列，每秒中隊伍最前面的人有 \(p\) 的機率走上電梯（一旦走上就不會下電梯），或者有 \(1-p\) 的機率不動。問你 \(T\) 秒事後，在電梯上的人的指望。

方法一

在本題這樣一個狀況中，方法一是用不了的，由於咱們的結束狀態不明確。

方法三：使用指望的定義計算

若是 \(X\) 是離散的隨機變量，輸出值爲 \(x_1,x_2,...\)，輸出值相應的機率爲 \(p_1,p_2,...\)，那麼指望值是一個無限數列的和（若是不收斂，那麼指望不存在）：

\[E[x] =\sum\limits_ip_ix_i \]

在本題中，若是設 \(dp(i,j)\) 爲 \(i\) 秒事後，電梯上有 \(j\) 我的的機率，那麼答案是：

\[\sum\limits_{0\le k\le n}dp(T,K)*K \]

因此咱們只須要求 \(dp(i, j)\) 就能夠了，初始值 \(dp(0, 0) = 1\) 就能夠了，仍然是採用刷表法：

\[dp(i + 1,j + 1) \leftarrow dp(i,j)*p\\ dp(i + 1,j)\leftarrow dp(i,j) * (1 - p) \]

【AC Code】

const int N = 2e3 + 10;
double p, dp[N][N];
int main() {
    cin.tie(nullptr)->sync_with_stdio(false);
    int n, t;
    cin >> n >> p >> t;
    dp[0][0] = 1;
    for (int i = 0; i < t; ++i) {
        dp[i + 1][n] += dp[i][n];
        for (int j = 0; j < n; ++j) if (dp[i][j] > 1e-10) {
                dp[i + 1][j + 1] += dp[i][j] * p;
                dp[i + 1][j] += dp[i][j] * (1 - p);
            }
    }
    double ans = 0;
    for (int i = 0; i <= n; ++i) ans += i * dp[t][i];
    cout << fixed << setprecision(6) << ans;
}

方法二

那麼以前提到的適用範圍廣的方法二，是否能在這裏用呢？答案是確定的。

延續方法三的 DP，咱們不妨將狀態之間的轉移抽象成邊，只不過只有 \(dp(i, j)\) 到 \(dp(i + 1, j + 1)\) 的邊纔有爲 \(1\) 的邊權，其他都爲 \(0\)。由於這個 DP 涵蓋了全部可能出現的狀況，因此咱們仍然能夠利用指望的線性性質，在刷表的過程當中進行計算答案。

本題中，沒有直觀的邊的概念，可是咱們能夠將狀態之間的轉移抽象成邊，因爲 \(dp(i, j)\)到 \(dp(i + 1, j + 1)\) 這一個轉移是對答案有 \(1\) 的貢獻的，因此咱們將它們之間的邊權賦爲 \(1\)。

這一題將方法二抽象化了，實際上大多數題並不是是直觀的，而是這種抽象的形式。

const int N = 2e3 + 10;
double p, dp[N][N];
int main() {
    cin.tie(nullptr)->sync_with_stdio(false);
    int n, t;
    cin >> n >> p >> t;
    dp[0][0] = 1;
    double ans = 0;
    for (int i = 0; i < t; ++i) {
        dp[i + 1][n] += dp[i][n];
        for (int j = 0; j < n; ++j) if (dp[i][j] > 1e-10) {
                dp[i + 1][j + 1] += dp[i][j] * p;
                dp[i + 1][j] += dp[i][j] * (1 - p);
                ans += dp[i][j] * p;
            }
    }
    cout << fixed << setprecision(6) << ans;
}

例三

BZOJ 3450

題意：給定一個序列，一些位置未肯定（是 \(\texttt{o}\) 與 \(\texttt{x}\) 的概率各佔 \(\%50%\)）。對於一個 \(\texttt{ox}\) 序列，連續 \(x\) 長度的 \(\texttt{o}\) 會獲得 \(x^2\) 的收益，請問最終獲得的序列的指望收益是多少？

分析

這個題若是一段一段的處理，實際上並非很好作。咱們觀察到 \((x + 1) ^ 2 - x ^ 2 = 2x + 1\)，那麼根據指望的線性性質，咱們能夠單獨算每個字符的貢獻。咱們設 \(dp_i\) 爲考慮前 ii 個字符的指望得分，\(l_i\) 爲以 \(i\) 爲結尾的 comb 的指望長度，\(Comb_i\) 爲第 \(i\)個字符，那麼有 3 種狀況：

1. \(s_i = o\) ，則 \(dp_i = dp_{i - 1} + l_{i - 1} * 2 + 1,l_i = l_{i - 1} + 1\)
2. \(s_i = x\) ，則 \(dp_i = dp_{i - 1}\)
3. \(s_i =\ ?\)， 則 \(dP_i = dp_{i - 1} + \frac{l_i*2 + 1}{2},l_i = \frac{l_{i - 1} + 1}{2}\)

對於前兩種狀況，實際上是很是直觀的，對於第三種狀況，其實是求了一個平均長度。例如 ?oo，兩種狀況的長度 \(l_i\) 分別爲 \([0,1,2]\) 和 \([1,2,3]\) ，可是求了平均以後，長度 \(l_i\) 變成了 \([0.5,1.5,2.5]\) ，這樣因爲咱們的貢獻是一個關於長度的一次多項式 \((2x + 1)\) ，因此長度平均以後，貢獻也至關於求了一個平均，天然可以求得正確的得分指望。

【AC Code】

const int N = 3e5 + 10;
double dp[N], Comb[N];
int main() {
    cin.tie(nullptr)->sync_with_stdio(false);
    int n; string s;
    cin >> n >> s;
    for (int i = 0; i < n; ++i) {
        if (s[i] == 'o') {
            dp[i] = dp[i - 1] + Comb[i - 1] * 2 + 1;
            Comb[i] = Comb[i - 1] + 1;
        } else if (s[i] == 'x') {
            dp[i] = dp[i - 1];
            Comb[i] = 0;
        } else {
            dp[i] = dp[i - 1] + (Comb[i - 1] * 2 + 1) / 2;
            Comb[i] = (Comb[i - 1] + 1) / 2;
        }
    }
    cout << setprecision(4) << fixed << dp[n - 1];
}

思考：若是長度爲 \(a\) 的 comb 的貢獻爲 \(a^3\) 時該如何解決？

Tips：因爲 \((a + 1)^3 - a^3 = 3a^3 + 3a + 1\) ，因此咱們要維護 \(a^2\) 和 \(a\) 的指望，注意 \(E_{a^2} \not= E^2_a\)，因此維護 \(a^2\) 的指望是必要的。

---

例四

BZOJ 4318

題意：給定一個序列，每一個位置 \(o\) 的概率爲 \(p_i\) ，爲 \(x\) 的概率爲 \(1-p_i\) 。對於一個 \(\texttt{ox}\) 序列，連續 \(x\) 長度的 \(\texttt{o}\) 會獲得 \(x^3\) 的收益，請問最終獲得的 \(ox\) 序列的指望收益是多少？

分析

延續例三的思路，咱們仍是分別求每個位置的貢獻。根據 \((a + 1)^3 - a^3 = 3a^3 + 3a + 1\)，咱們只須要維護 \(l(i)\)爲以 \(i\) 爲結尾的 comb 的指望長度，\(l_2(i)\)爲以 \(i\) 爲結尾的 comb 的指望長度的平方。注意 \(E[a^2] \not =E^2[a]\)，因此維護 \(a^2\) 的指望是必要的。

int main() {
    cin.tie(nullptr)->sync_with_stdio(false);
    int n;
    double p, l1 = 0, l2 = 0, ans = 0;
    cin >> n;
    for (int i = 0; i < n; ++i) {
        cin >> p;
        ans += (3 * l2 + 3 * l1 + 1) * p;
        l2 = (l2 + 2 * l1 + 1) * p;
        l1 = (l1 + 1) * p;
    }
    cout << fixed << setprecision(1) << ans;
}

總結

指望 DP 通常來講有它固定的模式，一種模式是直接 DP，定義狀態爲到終點指望，採用逆序計算獲得答案。一種模式是利用指望的線性性質，對貢獻分別計算，這種模式通常要求咱們求出每種代價的指望使用次數，而每種代價每每體如今 DP 的轉移之中。最後的兩個例題是典型的分離變量，用指望的線性性質計算答案的例子，若是狀態過於巨大，那麼就得考慮分離隨機變量了。

本總結只是解釋了機率與指望 DP 的冰山一角，它能夠變化無窮，但那些實際上並不僅屬於機率與指望 DP，真正核心的內容，仍是逃不出咱們幾種方法。

想要深刻了解一些機率的DP的請閱讀這篇文章：Here