BZOJ 3270: 博物馆

最新推荐文章于 2020-01-13 20:23:43 发布

转载最新推荐文章于 2020-01-13 20:23:43 发布 · 68 阅读

0 ·

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/LLTYYC/p/10837689.html

本文探讨了在概率图模型中，如何通过状态转移方程和矩阵优化来求解两人在随机图上相遇的概率问题。介绍了状态转移的概念，详细解析了状态转移方程的建立与求解过程，包括矩阵优化的暴力转移方法和使用高斯消元解方程组的精确解法。

传送门

显然可以状态转移：

设 $f[k][x][y]$ 表示第 $k$ 时刻，第一个人在 $x$ ，第二个人在 $y$ 时的概率

那么转移显然：

$f[k][x][y]+=\sum_{u}\sum_{v}f[k-1][u][v]*(1-P_u)(1-P_v)/du[u]/du[v]$

其中 $u$ 和 $x$ 有边相连，$v$ 和 $y$ 有边向连，$du[i]$ 表示节点 $i$ 的度数，并且 $u!=v$

当然这只是一种情况的转移，还有三种情况：

$f[k][x][y]+=\sum_{u}f[k-1][u][y]*(1-P_u)P_y/du[u]$

$f[k][x][y]+=\sum_{v}f[k-1][x][v]*(1-P_v)P_x/du[v]$

$f[k][x][y]+=f[k-1][x][y]*P_x*P_y$

然后显然可以矩阵优化暴力转移，复杂度 $O(n^6log_k)$，$k$ 是转移步数，精度玄学

正解是考虑设 $f[x][y]$ 表示两人从起点到终点，经过状态 $(x,y)$ 即第一个人在 $x$，第二个人在 $y$ 的期望次数

状态转移的方程好像也差不多

$f[x][y]+=\sum_{u}\sum_{v}f[u][v]*(1-P_u)(1-P_v)/du[u]/du[v]$

$f[x][y]+=\sum_{u}f[u][y]*(1-P_u)P_y/du[u]$

$f[x][y]+=\sum_{v}f[x][v]*(1-P_v)P_x/du[v]$

$f[x][y]+=f[x][y]*P_x*P_y$

因为两人在同一个点时就不会继续走了，所以从起点出发两人在 $i$ 点相遇的概率就是两人从起点到终点，经过状态 $(i,i)$ 的期望次数（$f[i][i]$）

这个方程转移有环，可以列出所有转移方程然后用高斯消元解方程组

复杂度 $O(n^6)$

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
typedef long long ll;
typedef double db;
inline int read()
{
    int x=0,f=1; char ch=getchar();
    while(ch<'0'||ch>'9') { if(ch=='-') f=-1; ch=getchar(); }
    while(ch>='0'&&ch<='9') { x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
    return x*f;
}
const int N=507;
int fir[N],from[N<<1],to[N<<1],cntt;
inline void add(int a,int b) { from[++cntt]=fir[a],fir[a]=cntt,to[cntt]=b; }
int n,m,nn,pa,pb,id[N][N],du[N];
db P[N],A[N][N],ans[N];
void Gauss()//高斯消元解方程组
{
    int pos;
    for(int i=1;i<=nn;i++)
    {
        pos=0;
        for(int j=i;j<=nn;j++) if(fabs(A[j][i])>fabs(A[pos][i])||!pos) pos=j;
        swap(A[i],A[pos]);
        for(int j=i+1;j<=nn;j++)
        {
            db w=A[j][i]/A[i][i];
            for(int k=i;k<=nn+1;k++) A[j][k]-=w*A[i][k];
        }
    }
    for(int i=nn;i;i--)
    {
        for(int j=i+1;j<=nn;j++) A[i][nn+1]-=ans[j]*A[i][j];
        ans[i]=A[i][nn+1]/A[i][i];
    }
}
int main()
{
    n=read(),m=read(),pa=read(),pb=read();
    nn=n*n; int now=0,a,b;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++) id[i][j]=++now;
    for(int i=1;i<=m;i++)
    {
        a=read(),b=read();
        add(a,b); add(b,a);
        du[a]++; du[b]++;
    }
    for(int i=1;i<=n;i++) scanf("%lf",&P[i]);
    for(int i=1;i<=n;i++)//枚举u
        for(int j=1;j<=n;j++)//枚举v
        {
            if(i==j) continue;//注意如果i=j就不能转移
            for(int k=fir[i];k;k=from[k])//枚举x
                for(int l=fir[j];l;l=from[l])//枚举y
                    A[ id[to[k]][to[l]] ][id[i][j]]-=(1.0-P[i])*(1.0-P[j])/du[i]/du[j];
            for(int k=fir[i];k;k=from[k]) A[ id[to[k]][j] ][id[i][j]]-=(1.0-P[i])*P[j]/du[i];
            for(int k=fir[j];k;k=from[k]) A[ id[i][to[k]] ][id[i][j]]-=(1.0-P[j])*P[i]/du[j];
            A[id[i][j]][id[i][j]]-=P[i]*P[j];
            //构造矩阵
        }
    for(int i=1;i<=nn;i++) A[i][i]++;
    A[id[pa][pb]][nn+1]=1;//起点初始为1
    Gauss();
    for(int i=1;i<=n;i++) printf("%.6lf ",ans[id[i][i]]);
    return 0;
}