概率期望 DP 总结

关于概率期望 DP，首先要了解概率及期望的一些基本定义概率 & 期望。

我们需要用到的比较重要的公式就是期望的定义和期望的线性性质。

定义

E(x)=\sum _{\omega\in \Omega}P(\omega)W(\omega)

即期望等于概率乘权值。需要注意的是， $x$ 是一个随机变量而并非一个事件，你可以把随机变量理解成一个函数，函数的自变量是样本点，即事件的某种可能的结果，定义域是一个样本空间，而函数的值域则是 $\mathbb{R}$ ，即随机变量相当于事件的某个结果往权值的一个映射。因此， $E(x)$ 是对 $x$ 的定义域内的所有自变量，即所有结果的概率和映射过来的权值的积求和，其中 $\Omega$ 就是定义域，一个样本空间了， $\omega$ 则是自变量，事件的某个结果， $P(\omega)$ 则是导致这个结果发生的概率， $W(\omega)$ 则是随机变量 $x$ 把 $\omega$ 当作自变量映射过来的权值，一个实数。

例如，现在我们掷两枚骰子，我们求掷出的点数之和的期望。那么随机变量 $x$ 的值域则是 $2\sim 12$ ，定义域则是不同的事件结果，即掷出的点数之和是多少，假设 $\omega$ 指掷出点数之和为 $2$ ，那么只有两次都掷出 $1$ 点才行， $P(\omega)=\dfrac{1}{6}\times \dfrac{1}{6}=\dfrac{1}{36}$ 。因此，可得：

\begin{aligned} E(x)& = \dfrac{1}{36}\times(2+12)+\dfrac{1}{18}\times(3+11)+\dfrac{1}{12}\times(4+10)\\ & +\dfrac{1}{9}\times(5+9)+\dfrac{5}{36}\times(6+8)+\dfrac{1}{6}\times7=7 \end{aligned}

性质

期望是线性函数，满足：

E(aX+bY)=a\times E(X)+b\times E(Y)

首先要了解，随机变量和函数不同的地方，函数相加定义域取交集，而两个随机变量相加它们的定义域样本空间应该是取笛卡尔积。

如此，我们便能证明这条性质：

\begin{aligned} E(X+Y) & =\sum _{x}\sum _{y}(x+y)P(X=x,Y=y)\\ & =\sum _{x}\sum _{y}xP(X=x,Y=y)+\sum _{x}\sum _{y}yP(X=x,Y=y)\\ & =\sum _{x} x\sum _{y}P(X=x,Y=y)+\sum _{y}y\sum _{x}P(X=x,Y=y)\\ & =\sum _{x} xP(X=x)+\sum _{y}yP(Y=y)\\ & =E(x)+E(y) \end{aligned}

其中第三步到第四步将 $P(X=x,Y=y)$ 拆开用到了一个概率公式，即对于两个不相关的事件 $X,Y$ ，满足 $P(XY)=P(X)P(Y)$ ，然后由于 $\sum _{x}P(X=x)=\sum _{y} P(Y=y)=1$ ，因此可以转成第四步。

好了，现在开始讲概率 DP，有很多存在后效性的 DP 也是概率 DP。

例 1

一个 $n\times m$ 的地图，从 $(1,1)$ 出发，每次随机向上下左右 $4$ 个方向移动到一个无障碍物的格子，求移动到 $(n,m)$ 的期望步数。

HDU2262 Where is the canteen

解析：这是一个图中求到某一点的期望步数的套路题，我们尝试从定义严格来推式子。

假设 $E _{u}$ 为点 $u$ 到终点 $t$ 的期望步数。

那么可得出：

\begin{aligned} E _{u} & =\sum _{i} ^{u\to t}p _{i}w _{i} \end{aligned}

其中 $i$ 枚举的是 $u\to t$ 的一条路径， $p _{i}$ 为走这条路径的概率， $w _{i}$ 为这条路径的长度。

接着推：

\begin{aligned} & =\sum _{j} ^{u\to v}\sum _{i} ^{v\to t}p _{i}p _{j}(w _{i}+w(u,v)) \end{aligned}

这里把路径 $u\to t$ 拆成了 $u\to v\to t$ ，然后概率由乘法原理变为 $p _{i}p _{j}$ ，权值也同样拆开。

\begin{aligned} & =\sum _{j} ^{u\to v}p _{j}\sum _{i} ^{v\to t}p _{i}w _{i}+p _{i}w(u,v) \end{aligned}

发现 $p _{j}$ 只与 $j$ 有关，把它跟 $j$ 的枚举放在一起，然后把括号用分配律打开。

\begin{aligned} & = \sum _{j} ^{u\to v}p _{j}(E _{v}+\sum _{i} ^{v\to t}p _{i}w(u,v)) \end{aligned}

发现 $\sum _{i} ^{v\to t}p _{i}w _{i}$ 就是 $E _{v}$ ，把它替换出来。

\begin{aligned} & = \sum _{j} ^{u\to v}p _{j}E _{v}+\sum _{i} ^{v\to t}p _{i}\sum _{j} ^{u\to v}p _{j}w(u,v) \end{aligned}

分配律，然后发现 $p _{i}$ 只与 $i$ 有关， $p _{j},w(u,v)$ 与 $j$ 有关，把它们各自放到一边。

\begin{aligned} & =\sum _{j} ^{u\to v}p _{j}(E _{v}+w(u,v)) \end{aligned}

可以发现 $\sum _{i} ^{v\to t}p _{i}=1$ ，然后再合并同类项即可。

于是得到了期望步数的套路式子：

\begin{aligned} E _{u} & =\sum _{j} ^{u\to v}p _{j}(E _{v}+w(u,v))\\ & =\sum ^{u\to v} \dfrac{(E _{v}+w(u,v))}{\deg _{u}} \end{aligned}

对于这道题 $p _{j}=\deg _{u}$ ，然后由于这道题能到达 $4$ 个不同的方向，即后面的状态也会影响到之前的状态，即 DP 具有后效性，于是采用高斯消元即可。需要注意的是先要 bfs 预处理出能到达终点的所有位置，对这些能到达的点才列出方程消元。同时注意高斯消元第一步消成的矩阵并不是上三角矩阵，而是阶梯形矩阵，所以回代时要小心。时间复杂度为 $O(n ^{6})$ 。

/*
 * @Author: clorf 
 * @Date: 2020-09-07 21:45:46 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-07 22:47:27
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<queue>
#include<cassert>
#define INF 1e9
using namespace std;
const int maxn=310;
const double Pi=acos(-1.0);
// 数组一定要算好开！！！！边和点的关系
// 记得看开不开long long！！！
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
// n和m不要搞混！！！任何循环时候都要注意！！
int n,m,deg[maxn];
int sx,sy,dir[4][2]={{0,1},{0,-1},{1,0},{-1,0}};
bool mark[maxn][maxn];
char s[maxn][maxn];
double a[maxn][maxn],ans[maxn],eps=1e-8;
inline int point(int x,int y)
{
    return (x-1)*m+y;
}
queue<pair<int,int> > q;
inline void bfs()
{
    while(!q.empty())
    {
        int nowx=q.front().first;
        int nowy=q.front().second;
        q.pop();
        for(int i=0;i<=3;i++)
        {
            int tx=nowx+dir[i][0];
            int ty=nowy+dir[i][1];
            if(tx<1||tx>n||ty<1||ty>m||mark[tx][ty]||s[tx][ty]=='#')
                continue;
            mark[tx][ty]=1;
            q.push(make_pair(tx,ty));
        }
    }
}
inline bool gauss()
{
    int i=1,j=1;
    while(i<=point(n,m)&&j<=point(n,m))
    {
        int maxx=i;
        for(int k=i;k<=point(n,m);k++)
            if(fabs(a[k][j])>fabs(a[maxx][j]))
                maxx=k;
        if(fabs(a[maxx][j])>=eps)
        {
            swap(a[i],a[maxx]);
            for(int k=i+1;k<=point(n,m);k++)
            {
                double num=a[k][j]/a[i][j];
                for(int l=j;l<=point(n,m)+1;l++)
                    if(fabs(a[i][l])>=eps)
                        a[k][l]-=num*a[i][l];
            }
            i++;
        }
        j++;
    }
    for(int k=i;k<=point(n,m);k++)
        if(fabs(a[k][point(n,m)+1])>=eps)
            return 0;
    i--;
    for(int k=i;k>=1;k--)
    {
        for(j=1;j<=point(n,m);j++)
            if(fabs(a[k][j])>=eps)
                break;
        ans[j]=a[k][point(n,m)+1]/a[k][j];
        for(int l=1;l<k;l++)
            a[l][point(n,m)+1]-=ans[j]*a[l][j];
    }
    return 1;
}
int main()
{
    while(scanf("%d%d",&n,&m)!=EOF)
    {
        memset(mark,0,sizeof(mark));
        memset(a,0,sizeof(a));
        memset(ans,0,sizeof(ans));
        memset(deg,0,sizeof(deg));
        memset(s,0,sizeof(s));
        for(int i=1;i<=n;i++)
        {
            scanf("\n");
            for(int j=1;j<=m;j++)
            {
                scanf("%c",&s[i][j]);
                if(s[i][j]=='@')
                {
                    sx=i;
                    sy=j;
                }
                if(s[i][j]=='$')
                {
                    mark[i][j]=1;
                    q.push(make_pair(i,j));
                }
            }
        }
        bfs();
        for(int i=1;i<=n;i++)
            for(int j=1;j<=m;j++)
            {
                if(!mark[i][j]||s[i][j]=='#')
                    continue;
                if(s[i][j]=='$')
                {
                    a[point(i,j)][point(i,j)]=1;
                    a[point(i,j)][point(n,m)+1]=0;
                    continue;
                }
                for(int k=0;k<4;k++)
                {
                    int tx=i+dir[k][0];
                    int ty=j+dir[k][1];
                    if(!mark[tx][ty]||s[tx][ty]=='#'||tx<1||tx>n||ty<1||ty>m)
                        continue;
                    deg[point(i,j)]++;
                    a[point(i,j)][point(tx,ty)]=1.0;
                }
                a[point(i,j)][point(i,j)]=-deg[point(i,j)];
                a[point(i,j)][point(n,m)+1]=-deg[point(i,j)];
            }
        bool flag=gauss();
        if(mark[sx][sy]&&flag)
            printf("%.6lf\n",ans[point(sx,sy)]);
        else
            printf("-1\n");
    }
    return 0;
}
/*
2 2
@#
.$
*/

例 2

二维方格，给出每个点向四方走的概率，保证和为 $1$ ，并且向边界外走的概率为 $0$ ，求从 $(1,1)$ 走到 $(n,m)$ 的期望步数。

UVA12177 First Knight

解析：就是上道题那个套路的模板题，直接推式子高斯消元。

/*
 * @Author: clorf 
 * @Date: 2020-09-11 09:22:02 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-11 11:50:49
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#define INF 1e9
using namespace std;
const int maxn=45;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int n,m,dir[4][2]={{1,0},{0,1},{-1,0},{0,-1}},r;
double p[maxn][maxn][4],a[maxn*maxn][maxn*maxn];
inline int point(int x,int y)
{
    if(x<1||x>n||y<1||y>m)
        return n*m+2;
    return (x-1)*m+y;
}
inline void gauss(int m,int n)
{
    int i=m,j=n;
    while(i>=1&&j>=1)
    {
        if(fabs(a[i][j]))
        {
            int downi=max(i-r,1);
            for(int k=i-1;k>=downi;k--)
            {
                double num=a[k][j]/a[i][j];
                int downj=max(j-r,1);
                for(int l=j;l>=downj;l--)
                    a[k][l]-=a[i][l]*num;
                a[k][n+1]-=a[i][n+1]*num;
            }
            i--;
        }
        j--;
    }
}
int main()
{
    while(scanf("%d%d",&n,&m)!=EOF&&n+m)
    {
        r=m;
        memset(a,0,sizeof(a));
        for(int k=0;k<4;k++)
            for(int i=1;i<=n;i++)
                for(int j=1;j<=m;j++)
                    scanf("%lf",&p[i][j][k]);
        for(int i=1;i<=n;i++)
            for(int j=1;j<=m;j++)
            {
                int s=point(i,j);
                a[s][s]=1.0;
                if(i==n&&j==m)
                    continue;
                a[s][point(n,m)+1]=1.0;
                for(int k=0;k<4;k++)
                    a[s][point(i+dir[k][0],j+dir[k][1])]=-p[i][j][k];
            }
        gauss(point(n,m),point(n,m));
        printf("%lf\n",a[1][n*m+1]/a[1][1]);
    }
    return 0;
}

例 3

跟上题相同，把步长改成 $2$ ，且一个格子只能向右向下走，还有概率停在本格子。

HDU3853 LOOPS

解析：跟上题套路一样。。。只不过把边权变成 $2$ ，然后相应的列式子：

E _{i,j}=p _{i,j,0}E _{i,j}+p _{i,j,1}E _{i,j+1}+p _{i,j,2}E _{i+1,j}

由于只有向右和向下的方向，所以这时前面的(左上的)状态都可以由后面的(右下的)状态得到，直接递推即可，不用高斯消元。

例 4

无向带权图，求从 $1$ 号节点到 $n$ 号节点的期望 XOR 和路径。

Luogu3211 XOR和路径

解析：自然而然地想到用刚才地定义式子一样来推。但是有一个问题，我们会在这一步地时候卡住：

\begin{aligned} & \sum _{j} ^{u\to v}\sum _{i} ^{v\to t}p _{i}p _{j}(w _{i}~\mathrm{xor}~w(u,v)) \end{aligned}

我们原套路中是加法，因此可以用分配律分开，但在这里是异或运算，异或并不满足分配律。简单来说，即异或的期望不等同于期望的异或。

这时候我们就需要换状态了。由于异或各位之间不互相影响，我们设 $f _{u,i}$ 表示 $u$ 到 $n$ 的路径期望异或和第 $i$ 位为 $1$ 的概率，这样样本点就转化到了只有 $2$ 个，为 $1$ 或为 $0$ ，而不是之前多条路径多个样本点。这样就可以列出式子：

\begin{aligned} f _{u,i} & =\dfrac{\sum _{w(u,v,i)=0}f _{v,i}+\sum _{w(u,v,i)=1}(1-f _{v,i})}{\deg _{u}}\\ & \Downarrow\\ \deg _{u}f _{u,i} & =\sum _{w(u,v,i)=0}f _{v,i}+\sum _{w(u,v,i)=1}(1-f _{v,i})\\ & \Downarrow\\ \deg _{u}f _{u,i}-\sum _{w(u,v,i)=0}f _{v,i}+\sum _{w(u,v,i)=1}f _{v,i} & =\sum _{w(u,v,i)=1} 1 \end{aligned}

然后对每一位高斯消元即可，答案即为 $\sum _{i}2 ^{i}f _{1,i}$ 。注意自环的情况，需要特判。

/*
 * @Author: clorf 
 * @Date: 2020-09-08 11:36:59 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-08 20:38:53
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#define INF 1e9
using namespace std;
const int maxn=210;
const double Pi=acos(-1.0);
// 数组一定要算好开！！！！边和点的关系
// 记得看开不开long long！！！
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
// n和m不要搞混！！！任何循环时候都要注意！！
int n,m;
int head[maxn],cnt,maxx,deg[maxn];
double a[maxn][maxn],ans[maxn],sum,eps=1e-8;
struct node
{
    int next;
    int to;
    int num;
}e[20010];
void add(int from,int to,int num)
{
    e[++cnt].next=head[from];
    e[cnt].to=to;
    e[cnt].num=num;
    head[from]=cnt;
    deg[from]++;
}
inline void connect(int x)
{
    a[n][n]=1;
    for(int u=1;u<n;u++)
    {
        a[u][u]=deg[u];
        for(int j=head[u];j;j=e[j].next)
        {
            int v=e[j].to;
            if((e[j].num&x)==0)
                a[u][v]--;
            else
            {
                a[u][v]++;
                a[u][n+1]++;
            }
        }
    }
}
inline void gauss()
{
    int i=1,j=1;
    while(i<=n&&j<=n)
    {
        int maxx=i;
        for(int k=i+1;k<=n;k++)
            if(fabs(a[k][j])>fabs(a[maxx][j]))
                maxx=k;
        if(fabs(a[maxx][j]))
        {
            swap(a[maxx],a[i]);
            for(int k=i+1;k<=n;k++)
            {
                double num=a[k][j]/a[i][j];
                for(int l=j;l<=n+1;l++)
                    a[k][l]-=num*a[i][l];
            }
            i++;
        }
        j++;
    }
    for(i=n;i>=1;i--)
    {
        for(j=i+1;j<=n;j++)
            a[i][n+1]-=ans[j]*a[i][j];
        ans[i]=a[i][n+1]/a[i][i];
    }
}
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=1;i<=m;i++)
    {
        int a,b,c;
        scanf("%d%d%d",&a,&b,&c);
        add(a,b,c);
        if(a!=b)
            add(b,a,c);
        maxx=max(maxx,c);
    }
    for(int i=1;i<=maxx;i<<=1)
    {
        memset(a,0,sizeof(a));
        memset(ans,0,sizeof(ans));
        connect(i);
        gauss();
        sum+=ans[1]*i;
    }
    printf("%.3lf",sum);
    return 0;
}

例 5

现在有一个 $n$ 个人的队列，你在第 $m$ 个位置，接下来会发生 $4$ 种不同的事件，每种事件都有一个概率：

队列保持不变，有 $p _{1}$ 的概率发生；
队列第一个人移到队尾，其他人前进一个位置，有 $p _{2}$ 的概率发生；
第一个人离开队列，其他人前进一个位置，有 $p _{3}$ 的概率发生；
队列冻结，永远保持现在状态不变，有 $p _{4}$ 的概率发生。

满足 $p _{1}+p _{2}+p _{3}+p _{4}=1$ 。

现在要求在你离开队列之前，队列已冻结，并且你在第前 $k$ 个位置的概率。

HDU4089 Activation

解析：设 $f _{i,j}$ 代表队列你前面有 $i$ 个人总共有 $j$ 个人转移到目标状态的概率。

稍微分类讨论一下，很容易推出转移式：

f _{i,j}=\begin{cases} p _{1}f _{i,j}+p _{2}f _{j-1,j}+p _{3}\times 0+p _{4}\times 1 && i=0\\ p _{1}f _{i,j}+p _{2}f _{i-1,j}+p _{3}f _{i-1,j-1}+p _{4}\times 1 && 1\le i\le k-1\\ p _{1}f _{i,j}+p _{2}f _{i-1,j}+p _{3}f _{i-1,j-1}+p _{4}\times 0 && i>k-1 \end{cases}

我们发现，如果确定了 $j$ ，那么对于每一个 $i\in [0,j-1]$ ，都会有一个 $f _{i,j}$ 的方程，即有一个 $j$ 个方程 $j$ 个未知数的方程组，但是如果用高斯消元的话很明显会超时。我们于是采取代入消元法，把 $f _{i,j}$ 的式子代入 $f _{i+1,j}$ 进去，这样整理后会改变的只有 $p _{1}$ 这个未知数的系数和 $p _{2},p _{3}$ 两个已知数的系数( $p _{4}$ 对应那项只有可能是 $p _{4}$ 和 $0$ )，于是这样代入 $n$ 次后，就是一个关于 $f _{j-1,j}$ 的一个一元一次方程，且最多只有 $4$ 个项，直接暴力求解，然后再将解出的 $f _{j-1,j}$ 回代回每一个方程，就能在 $O(j)$ 的时间内解出所有未知数，高斯消元则要 $O(j ^{3})$ 。最后答案即为 $f _{m-1,n}$ 。注意要开滚动数组优化空间。

/*
 * @Author: clorf 
 * @Date: 2020-09-09 09:30:25 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-09 11:05:44
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#define INF 1e9
using namespace std;
const int maxn=2010;
const double Pi=acos(-1.0);
// 数组一定要算好开！！！！边和点的关系
// 记得看开不开long long！！！
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
// n和m不要搞混！！！任何循环时候都要注意！！
int n,m,k;
double p1,p2,p3,p4,f[maxn][2];
double a,b[maxn],c[maxn],A,B,C,eps=1e-8;
int main()
{
    while(scanf("%d%d%d",&n,&m,&k)!=EOF)
    {
        memset(f,0,sizeof(f));
        scanf("%lf%lf%lf%lf",&p1,&p2,&p3,&p4);
        if(p4<=eps)
        {
            printf("0.00000\n");
            continue;
        }
        f[0][1]=p4/(p3+p4);
        for(int j=1;j<=n;j++)
        {
            memset(b,0,sizeof(b));
            memset(c,0,sizeof(c));//系数
            a=A=B=C=0.0;
            a=p2/(1-p1);
            b[0]=0.0;
            c[0]=p4/(1-p1);
            for(int i=1;i<=j;i++)
            {
                b[i]=p3/(1-p1)*f[i-1][(j%2)^1];
                if(i<=k-1)
                    c[i]=p4/(1-p1);
                else
                    c[i]=0.0;
            }
            A=a;
            B=b[j-1];
            C=c[j-1];
            for(int i=j-2;i>=0;i--)
            {
                B+=A*b[i];
                C+=A*c[i];
                A*=a;
            }
            f[j-1][j%2]=(B+C)/(1-A);
            f[0][j%2]=a*f[j-1][j%2]+b[0]+c[0];
            for(int i=1;i<=j-2;i++)
                f[i][j%2]=a*f[i-1][j%2]+b[i]+c[i];
        }
        printf("%.5lf\n",f[m-1][n%2]);
    }
    return 0;
}

例 6

一个无向图，节点 $1$ 有一个炸弹，在每个单位时间内，有可能在这个节点炸掉，也有 $p$ 的概率随机选择一条出去的路到其他的节点上。问最终炸弹在每个节点上爆炸的概率。

Luogu2973 Driving Out the Piggies G

解析：容易发现这道题和上道题都有一个共同点，它们都是起始状态唯一，但终止状态有多个的问题，但这道题最大的不同点是到达终止状态后立马结束，并且能结束程序的只有终止状态，而上道题在终止状态内还可以继续转移，并且在任何一个状态都可以结束程序。对于这类题，求在每个终止节点终止的概率，也是一类套路问题，我们一般设 $E _{u}$ 代表从起点出发直到到达终点结束，途中 $u$ 节点的期望经过次数作为状态，由于终点只会经过 $1$ 次便停止，所以终点的期望经过次数就是在这个点终止的概率(权值，即次数位 $1$ )。现在我们尝试从定义推期望经过次数的式子。

\begin{aligned} E _{u} & =\sum _{i} ^{s\to t}p _{i}w _{i} \end{aligned}

其中 $i$ 枚举的是从起点 $s$ 到终点 $t$ 的一条路径， $p _{i}$ 是走这条路径的概率， $w _{i}$ 是这条路径经过 $u$ 的次数。

我们发现到这里我们不能像之前那样把路径拆开处理了，我们需要借助其他方法。

我们设 $f _{k,u}$ 代表从起点 $s$ 出发，走 $k$ 条边正好到达 $u$ 的概率，即从 $s$ 出发走 $k$ 条边恰好到达 $u$ 的路径数除以从 $s$ 出发走 $k$ 条边的总路径数，注意这里概率的分子分母和 $p _{i}$ 的定义不一样。那么可以推出：

f _{k,u}=\sum \limits_{i,s\mathop{\to}\limits ^{k} u} ^{s\to t} p _{i}

左式和右式都是求的是 $s$ 走 $k$ 条边到达 $u$ 的概率，如果一条路径 $i$ 多次经过 $u$ ，那它也会多次被统计进不同的 $f _{k,u}$ 里面，因此相等。

于是我们可以继续处理式子：

\begin{aligned} & =\sum _{k=0} ^{\infty}\sum \limits_{i,s\mathop{\to}\limits ^{k} u} ^{s\to t} p _{i}\\ & =\sum _{k=0} ^{\infty}f _{k,u} \end{aligned}

发现我们上面那个式子转化的困难点就是不知道 $w _{i}$ ，即每条路径经过 $u$ 的次数是多少。于是我们转而枚举从起点 $s$ 出发到 $u$ 走的步数，这样依然能枚举完所有路径。这里的权值被消掉了，其实就是相当于把一条 $s\to t$ 的路径拆成了多条路径。假如一条 $s\to t$ 的路径多次经过 $u$ ，那么在这个式子里相当于把每一次经过 $u$ 都新建一条 $s\to t$ 的路径拆开，这样新建出来的路径对次数的贡献是 $1$ ，即只算这一次经过 $u$ ，之前经过 $u$ 不算在其中，如此就可以把 $w _{i}$ 分到多条路径中去统计，就可以把它消掉了。

接下来的难点就是 $f _{k,u}$ 的化简。它的状态跟之前做过的某题有点类似(HDU6331 Walking Plan)，但又不完全一样。我们可以通过举例或推导得出下面的式子：

f _{k+1,v}=\sum _{u=1} ^{n}f _{k,u}\times A _{v,u}

其中 $A _{i,j}$ 代表 $j\to i$ 的概率。

于是接着化简期望式子：

\begin{aligned} & = f _{0,u}+\sum _{(v,u)}A _{u,v}\sum _{k=0} ^{\infty}f _{k,v}\\ & =f _{0,u}+\sum _{(v,u)}A _{u,v}E _{v} \end{aligned}

枚举所有转移到 $u$ 的边 $(v,u)$ ，然后用 $f _{k,v}$ 来更新 $f _{k,v}$ ，注意 $f _{0,u}$ 要单独拆出来。然后发现 $\sum _{k=0} ^{\infty}f _{k,v}=E _{v}$ ，替换掉它即可。

在这道题中 $A _{u,v}$ 代表的是 $v\to u$ 的概率，其实就是 $\deg _{v}$ 。 $f _{0,u}$ 代表从起点 $s$ 出发，走 $0$ 步到达 $u$ 的概率，只有 $u=s$ 时概率为 $1$ ，其他情况为 $0$ ，因为走 $0$ 步只能到达它自己 $s$ 。

于是得出结论式：

\begin{aligned} E _{u} & =[u=s]+\sum _{(v,u)}\dfrac{E _{v}}{\deg _{v}} \end{aligned}

于是这道题就直接套式子，高斯消元即可。对每个终点，注意不仅要在那一点停止，还要爆炸，最后在乘以爆炸的概率即可。

/*
 * @Author: clorf 
 * @Date: 2020-09-11 15:10:06 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-11 15:39:11
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<vector>
#define INF 1e9
using namespace std;
const int maxn=310;
const double Pi=acos(-1.0);
const double eps=1e-8;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int n,m,P,Q;
double p,a[maxn][maxn],ans[maxn];
vector<int> e[maxn];
inline void gauss(int m,int n)
{
    int i=1,j=1;
    while(i<=m&&j<=n)
    {
        int maxx=i;
        for(int k=i+1;k<=m;k++)
            if(fabs(a[k][j])>fabs(a[maxx][j]))
                maxx=k;
        if(fabs(a[maxx][j])>eps)
        {
            swap(a[maxx],a[i]);
            for(int k=i+1;k<=m;k++)
            {
                double num=a[k][j]/a[i][j];
                for(int l=j;l<=n+1;l++)
                    a[k][l]-=num*a[i][l];
            }
            i++;
        }
        j++;
    }
    for(i=n;i>=1;i--)
    {
        for(j=i+1;j<=n;j++)
            a[i][n+1]-=a[i][j]*ans[j];
        ans[i]=a[i][n+1]/a[i][i];
    }
}
int main()
{
    scanf("%d%d%d%d",&n,&m,&P,&Q);
    p=1.0*P/Q;
    for(int i=1;i<=m;i++)
    {
        int a,b;
        scanf("%d%d",&a,&b);
        e[a].push_back(b);
        e[b].push_back(a);
    }
    for(int u=1;u<=n;u++)
    {
        a[u][u]=1.0;
        for(int i=0;i<e[u].size();i++)
        {
            int v=e[u][i];
            a[u][v]=(p-1.0)/e[v].size();
        }
        if(u==1)
            a[u][n+1]=1.0;
    }
    gauss(n,n);
    for(int i=1;i<=n;i++)
    {
        ans[i]*=p;
        printf("%.9lf\n",ans[i]);
    }
    return 0;
}

例 7

无向连通图，有两个人分别在 $s,t$ ，每一分钟若一个人在点 $u$ 时，则他有 $p _{u}$ 的概率不动，否则随机前往一个相邻的结点，求在每个点相遇的概率。

BZOJ3270 博物馆

解析：跟上题一样的套路题。分 $4$ 种情况照样推式子：

\begin{aligned} E _{u _{1},u _{2}} & =\dfrac{1-p _{v _{1}}}{\deg _{v _{1}}}p _{u _{2}}E _{v _{1},u _{2}}+\dfrac{1-p _{v _{2}}}{\deg _{v _{2}}}p _{u _{1}}E _{u _{1},v _{2}}\\ & +\dfrac{1-p _{v _{1}}}{\deg _{v _{1}}}\dfrac{1-p _{v _{2}}}{\deg _{v _{2}}}E _{v _{1},v _{2}}+p _{u _{1}}p _{u _{2}}E _{u _{1},u _{2}} \end{aligned}

同样的，起点 $E _{s,t}$ 要加 $1$ 。然后高斯消元。

/*
 * @Author: clorf 
 * @Date: 2020-09-09 21:27:03 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-10 21:31:53
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#define INF 1e9
using namespace std;
const int maxn=22;
const double Pi=acos(-1.0);
// 数组一定要算好开！！！！边和点的关系
// 记得看开不开long long！！！
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
// n和m不要搞混！！！任何循环时候都要注意！！
int n,m,a,b,cnt;
int deg[maxn];
bool dis[maxn][maxn];
double p[maxn],map[maxn*maxn][maxn*maxn],ans[maxn*maxn];
inline int point(int x,int y)
{
    return (x-1)*n+y;
}
inline void gauss(int m,int n)
{
    int i=1,j=1;
    while(i<=m&&j<=n)
    {
        int maxx=i;
        for(int k=i+1;k<=m;k++)
            if(fabs(map[k][j])>fabs(map[maxx][j]))
                maxx=k;
        if(fabs(map[maxx][j]))
        {
            swap(map[maxx],map[i]);
            for(int k=i+1;k<=m;k++)
            {
                double num=map[k][j]/map[i][j];
                for(int l=j;l<=n+1;l++)
                    map[k][l]-=num*map[i][l];
            }
            i++;
        }
        j++;
    }
    for(i=n;i>=1;i--)
    {
        for(int j=i+1;j<=n;j++)
            map[i][n+1]-=map[i][j]*ans[j];
        ans[i]=map[i][n+1]/map[i][i];
    }
}
int main()
{
    scanf("%d%d%d%d",&n,&m,&a,&b);
    for(int i=1;i<=m;i++)
    {
        int x,y;
        scanf("%d%d",&x,&y);
        dis[x][y]=dis[y][x]=1;
        dis[x][x]=dis[y][y]=1;
        deg[x]++;
        deg[y]++;
    }
    for(int i=1;i<=n;i++)
        scanf("%lf",&p[i]);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
        {
            int s=point(i,j);
            map[s][s]=-1;
            for(int t=1;t<=n*n;t++)
            {
                int y=(t-1)%n+1;
                int x=(t-y)/n+1;
                if(x==y)
                    continue;
                if(x==i&&y==j)
                    map[s][t]+=p[i]*p[j];
                else if(x==i&&dis[j][y])
                    map[s][t]+=p[i]*(1-p[y])/deg[y];
                else if(y==j&&dis[i][x])
                    map[s][t]+=p[j]*(1-p[x])/deg[x];
                else if(dis[i][x]&&dis[j][y])
                    map[s][t]+=(1-p[y])*(1-p[x])/(deg[x]*deg[y]);
            }
        }
    map[point(a,b)][point(n,n)+1]=-1;
    gauss(n*n,n*n);
    for(int i=1;i<=n;i++)
        printf("%.6lf ",ans[point(i,i)]);
    return 0;
}

例 8

$n$ 个点 $m$ 条边的无向图，初始在 $1$ 号节点，到达 $n$ 号节点结束，每一步随机走向下一个点，每走一条边获得这个边编号(编号从 $1\sim m$ )的分数，总分为获得的所有分数之和。现在对 $m$ 条边编号使总分期望值最小。

Luogu3232 游走

解析：很明显，根据倒序和小于乱序和小于正序和的定理，边的期望经过次数多的编号应该越小。我们将边的期望经过次数从大到小排序，然后依次按从 $1\sim m$ 编号即可。

边的期望经过次数由点的期望经过次数得到：

E _{(u,v)}=\dfrac{E _{u}}{\deg _{u}}+\dfrac{E _{v}}{\deg _{v}}

然后用套路高斯消元求出点的期望经过次数即可。

/*
 * @Author: clorf 
 * @Date: 2020-09-10 20:29:56 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-10 22:24:38
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<vector>
#define INF 1e9
using namespace std;
const int maxn=510,maxm=125010;
const double eps=1e-7;
const double Pi=acos(-1.0);
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
vector<int> p[maxn];
int n,m;
double a[maxn][maxn],ans[maxn],sum;
struct edge
{
    int x;
    int y;
    double E;
}e[maxm];
bool cmp(edge a,edge b)
{
    return a.E>b.E;
}
inline void gauss(int m,int n)
{
    int i=1,j=1;
    while(i<=m&&j<=n)
    {
        int maxx=i;
        for(int k=i+1;k<=m;k++)
            if(fabs(a[k][j])>fabs(a[maxx][j]))
                maxx=k;
        if(fabs(a[maxx][j])>eps)
        {
            swap(a[maxx],a[i]); 
            for(int k=i+1;k<=m;k++)
            {
                double num=a[k][j]/a[i][j];
                for(int l=j;l<=n+1;l++)
                    a[k][l]-=num*a[i][l];
            }
            i++;
        }
        j++;
    }
    for(i=n;i>=1;i--)
    {
        for(j=i+1;j<=n;j++)
            a[i][n+1]-=a[i][j]*ans[j];
        ans[i]=a[i][n+1]/a[i][i];
    }
}
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=1;i<=m;i++)
    {
        scanf("%d%d",&e[i].x,&e[i].y);
        p[e[i].x].push_back(e[i].y);
        p[e[i].y].push_back(e[i].x);
    }
    for(int i=1;i<n;i++)
    {
        a[i][i]=1.0;
        for(int j=0;j<(int)p[i].size();j++)
        {
            if(p[i][j]==n)
                continue;
            a[i][p[i][j]]=-(1.0/((double)p[p[i][j]].size()*1.0));
        }
        if(i==1)
            a[i][n]=1.0;
    }
    gauss(n-1,n-1);
    for(int i=1;i<=m;i++)
        e[i].E=ans[e[i].x]/(double)(p[e[i].x].size()*1.0)+ans[e[i].y]/(double)(p[e[i].y].size()*1.0);
    sort(e+1,e+m+1,cmp);
    for(int i=1;i<=m;i++)
        sum+=e[i].E*i;
    printf("%.3lf",sum);
    return 0;
}