DP训练——概率期望DP

2020-02-25 本文已影响0人云中翻月

概率期望DP

POJ2096
题意
一个软件有 $s$ 个子系统，会产生 $n$ 种bug。
某人一天发现一个bug，这个bug属于某种bug，发生在某个子系统中。
求找到所有的 $n$ 种bug，且每个子系统都找到bug，这样所要的天数的期望。
需要注意的是：bug的数量是无穷大的，所以发现一个bug，出现在某个子系统的概率是 $1/s$ ，属于某种类型的概率是 $1/n$ 。
题解
状态定义： $d[i,j]$ 表示目前已经收集了 $i$ 种bug， $j$ 种子系统的期望的剩余时间
边界： $d[0,0]$
目标： $d[n,s]=0$
转移方程：
分四种情况讨论。
$d[i,j]$ ：发现一个bug属于已经找到的 $i$ 种bug和 $j$ 个子系统中
$d[i+1,j]$ ：发现一个bug属于新的一种bug，但属于已经找到的 $j$ 种子系统
$d[i,j+1]$ ：发现一个bug属于已经找到的 $i$ 种bug，但属于新的子系统
$d[i+1,j+1]$ ：发现一个bug属于新的一种bug和新的一个子系统
以上四种的概率分别为：
$p0 = i*j / (n*s)$
$p1 = (n-i)*j / (n*s)$
$p2 = i*(s-j) / (n*s)$
$p3 = (n-i)*(s-j) / (n*s)$
$d[i,j] = p1*d[i,j] + p2*d[i+1,j] + p3*d[i,j+1] + p4*d[i+1,j+1] + 1$
整理后得到 $d[i,j] = ( n*s + (n-i)*j*d[i+1,j] + i*(s-j)*d[i,j+1] + (n-i)*(s-j)*d[i+1,j+1] )/( n*s - i*j )$
代码如下

/*

*/
#define method_1
#ifdef method_1
/*

*/
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<set>
#include<map>
#include<queue>
#include<stack>
#include<vector>
#include<cstring>
#include<cstdlib>
#include<iomanip>
#include<ctime>
#include<string>
#include<bitset>
#define D(x) cout<<#x<<" = "<<x<<"  "
#define E cout<<endl
using namespace std;
typedef long long ll;
typedef pair<int,int>pii;
const int maxn=1000+5;
const int maxs=1000+5;
const int INF=0x3f3f3f3f;
int n,s;
double d[maxn][maxs];
void dp(){
    d[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=(n-i)*j;
            double p2=i*(s-j);
            double p3=(n-i)*(s-j);
            d[i][j]=p1*d[i+1][j]+p2*d[i][j+1]+p3*d[i+1][j+1]+n*s;
            d[i][j]/=(n*s-i*j);
        }
    }
}
int main() {
//  ios::sync_with_stdio(false);
    //freopen("POJ2096.in","r",stdin);
    cin>>n>>s;
    dp();
    printf("%.4lf",d[0][0]);
    return 0;
}
#endif
#ifdef method_2
/*

*/

#endif
#ifdef method_3
/*

*/

#endif

CF 113D英文
 CF 113D中文
题意
$n(n<=22)$ 个点， $m(m<=n(n-1)/2)$ 条边的无向图。两个人开始分别位于 $a$ 点和 $b$ 点。
每个时刻，对于处在第 $i$ 个点的人来说，有 $p[i]$ 的概率不动， $1-p[i]$ 的概率去相邻的房间。
求两人在每个房间相遇的概率。
题解
状态定义： $f[i,j]$ 表示一个人在点 $i$ ，另一个人在点 $j$ 的期望相遇时间。
目标： $f[i,i]$
边界： $f[a,b]=0$
转移方程： $f[i][j] = p_i * p_j * f[i][j] + p_i * \sum \frac{1 - p_k}{du[k]} *f[i][k] + p_j * \sum \frac{1 - p_k}{du[k]} * f[k][j] + \sum \sum \frac{1 - p_k}{du[k]} \frac{1 - p_h}{du[h]} * f[k][h]$
关于DP方程的高斯消元：
因为状态定义是两维的，所以总共有 $n^2$ 个未知状态，需要 $n^2$ 个方程来描述。
具体地说，对于增广矩阵的第 $i$ 行，描述了所有与 $hash$ 后值为 $i$ 的状态有关的变量。
题解链接 https://www.cnblogs.com/mrclr/p/10885912.html
代码如下

/*

*/
#define method_1
#ifdef method_1
/*

*/
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<set>
#include<map>
#include<queue>
#include<stack>
#include<vector>
#include<cstring>
#include<cstdlib>
#include<iomanip>
#include<ctime>
#include<string>
#include<bitset>
#define D(x) cout<<#x<<" = "<<x<<"  "
#define E cout<<endl
using namespace std;
typedef long long ll;
typedef pair<int,int>pii;
const int maxn=22+2;
const int maxm=maxn*maxn;
const int INF=0x3f3f3f3f;
const double eps=1e-8;
int n,n2,m,a,b; //n2=n*n 为方程数
double p[maxn];
double f[maxn*maxn][maxn*maxn+1],ans[maxn*maxn];
//f为高斯矩阵，总共n^2行，每行对应一个未知数。共n^2+1列，第n^2+1列为方程右边的常数。
int head[maxn],tot=1;
struct node{
    int from,to;
}edge[maxm];
int deg[maxn];//出度
void add(int from,int to){
    deg[from]++,edge[++tot].from=head[from],head[from]=tot,edge[tot].to=to;
}
int _hash(int i,int j){
    return (i-1)*n+j;
}
void build(){
    f[_hash(a,b)][n2+1]=1.0;
    for(int x=1;x<=n;x++){
        for(int y=1;y<=n;y++){
            int u=_hash(x,y);
            f[u][u]=1.0; //确定最初的系数
            if(x!=y) f[u][u]-=p[x]*p[y];
            for(int i=head[x];i;i=edge[i].from){
                int l=edge[i].to;
                if(l==y) continue;
                int v=_hash(l,y);
                f[u][v]-=p[y]*(1-p[l])/deg[l];
            }
            for(int i=head[y];i;i=edge[i].from){
                int l=edge[i].to;
                if(l==x) continue;
                int v=_hash(x,l);
                f[u][v]-=p[x]*(1-p[l])/deg[l];
            }
            for(int i=head[x];i;i=edge[i].from){
                for(int j=head[y];j;j=edge[j].from){
                    int k=edge[i].to,l=edge[j].to;
                    if(k==l) continue;
                    int v=_hash(k,l);
                    f[u][v]-=(1-p[k])*(1-p[l])/deg[k]/deg[l];
                }
            }
        }
    }
}
void solve(){
    for(int i=1;i<=n2;i++){
        int pos=i;
        for(int j=i+1;j<=n2;j++) if(fabs(f[j][i])>fabs(f[pos][i])) pos=j;
        if(pos!=i){
            for(int j=1;j<=n2+1;j++) swap(f[pos][j],f[i][j]);
        }
        if(f[i][i]<eps) continue;
        double tp=f[i][i];
        for(int j=i;j<=n2+1;j++) f[i][j]/=tp;
        for(int j=i+1;j<=n2;j++){
            double tp=f[j][i];
            for(int k=i;k<=n2+1;k++) f[j][k]-=tp*f[i][k];
        }
    }
    for(int i=n2;i>=0;i--){
        ans[i]=f[i][n2+1];
        for(int j=i-1;j>=0;j--) f[j][n2+1]-=f[j][i]*ans[i];
    }
}
int main() {
//  ios::sync_with_stdio(false);
    //freopen("CF113D.in","r",stdin);
    cin>>n>>m>>a>>b;
    n2=n*n;
    int from,to;
    for(int i=1;i<=m;i++){
        cin>>from>>to;
        add(from,to);add(to,from);
    }
    for(int i=1;i<=n;i++) cin>>p[i];
    build();
    solve();
    for(int i=1;i<=n;i++) printf("%.6lf ",ans[_hash(i,i)]);
    return 0;
}
#endif
#ifdef method_2
/*

*/

#endif
#ifdef method_3
/*

*/

#endif

BZOJ1415
题意
$n(n<=1000)$ 个点， $e(e<=1000)$ 条边的无向图。聪聪开始在 $c$ 点，可可开始在 $m$ 点。
可可每次等概率的向周围点移动一格。聪聪每次朝着到可可距离缩小的点中编号最小的点移动，如果移动一次后没有遇到可可，就会再移动一次。
求聪聪遇到可可的期望时间。
题解
首先预处理出 $to[i,j]$ 表示聪聪从i向j移动时，两步后会抵达的位置。
然后开始dp。（因为dp的起点不确定，所以用记忆化搜索实现）
状态定义： $d[i,j]$ 表示聪聪在 $i$ 点，可可在 $j$ 点的期望移动时间
目标： $d[c,m]$
边界： $d[i,i]=0,d[i,j]=1$ (如果 $i$ 到 $j$ 的移动次数小于等于两次)
转移方程： $d[i,j]=\sum( \frac{d[to[i,j],k]}{out[j]+1})+\frac{d[to[i,j],j]}{out[j]+1}+1$ (其中 $j,k$ 之间存在边连接且 $k!=j$ )
代码如下

/*

*/
#define method_1
#ifdef method_1
/*

*/
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<set>
#include<map>
#include<queue>
#include<stack>
#include<vector>
#include<cstring>
#include<cstdlib>
#include<iomanip>
#include<ctime>
#include<string>
#include<bitset>
#define D(x) cout<<#x<<" = "<<x<<"  "
#define E cout<<endl
using namespace std;
typedef long long ll;
typedef pair<int,int>pii;
const int maxn=1000+5;
const int maxe=1000+5;
const int INF=0x3f3f3f3f;
int n,e;
int c,m;
int G[maxn][maxn];
int tot=1,head[maxn];
struct node{
    int from,to;
}edge[maxe<<1];
void add(int from,int to){
    edge[++tot].from=head[from],head[from]=tot,edge[tot].to=to;
}
int vis[maxn][maxn];
double d[maxn][maxn];
int to[maxn][maxn];
int dis[maxn][maxn];
int out[maxn];
queue<int>q;
void pre(){ //预处理i向j移动时下一步会移动到的位置
    for(int i=1;i<=n;i++){
        while(q.size()) q.pop();
        q.push(i);
        for(int j=1;j<=n;j++) dis[j][i]=INF;
        dis[i][i]=0;
        while(q.size()){
            int x=q.front();q.pop();
            for(int k=head[x];k;k=edge[k].from){
                int y=edge[k].to;
                if(dis[y][i]!=INF) continue;
                dis[y][i]=dis[x][i]+1;
                q.push(y);
            }
        }
    }
    /*
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            D(i);D(j);
            cout<<dis[i][j]<<" ";
        }
        cout<<endl;
    }
    */
    for(int x=1;x<=n;x++){
        for(int y=1;y<=n;y++){
            if(x==y) continue;
            if(G[x][y]) continue;
            int fs=INF,ss=INF;
            for(int i=head[x];i;i=edge[i].from){
                int s1=edge[i].to;
                if(dis[s1][y]+1!=dis[x][y]) continue;
                fs=min(fs,s1);
            }
            for(int i=head[fs];i;i=edge[i].from){
                int s2=edge[i].to;
                if(dis[s2][y]+1!=dis[fs][y]) continue;
                ss=min(ss,s2);
            }
            to[x][y]=ss;
        }
    }
    /*
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            D(i);D(j);
            cout<<to[i][j]<<endl;
        }
        cout<<endl;
    }
    */
}
double dfs(int x,int y){
    if(vis[x][y]) return d[x][y];
    if(x==y) return 0.0; //在同一个位置
    if(G[x][y]) return 1.0; //走一步到达
    if(to[x][y]==y) return 1.0; //走两步到达（注意，需要考虑这个边界，因为如果不考虑，则会多费一次可可的移动
    double ret=0.0;
    for(int i=head[y];i;i=edge[i].from){
        int z=edge[i].to;
        ret+=dfs(to[x][y],z)/(out[y]+1.0); //聪聪(x)先向着可可(y)移动到(to[x][y])，然后可可(y)再移动到(z)
    }
    ret+=dfs(to[x][y],y)/(out[y]+1.0);
    ret+=1.0;
    vis[x][y]=1;
    return d[x][y]=ret;
}
int main() {
//  ios::sync_with_stdio(false);
    //freopen("BZOJ1415.in","r",stdin);
    scanf("%d%d%d%d",&n,&e,&c,&m);
    int from,to;
    for(int i=1;i<=e;i++){
        scanf("%d%d",&from,&to);
        G[from][to]=G[to][from]=1; //记录两点是否相邻
        add(from,to);add(to,from);
        out[from]++;out[to]++;
    }
    pre();
    double ans=dfs(c,m);
    printf("%.3lf",ans);
    return 0;
}
#endif
#ifdef method_2
/*

*/

#endif
#ifdef method_3
/*

*/

#endif

BZOJ2707
题意
$n(n<=10000)$ 个点， $m(m<=1000000)$ 条边的有向图，保证所有 $SCC$ 大小不超过 $100$ 。
一个人从 $S$ 点出发，每次等概率的随即走向一个相邻点，求走到终点的期望步数。
题解
注意到题目保证所有 $SCC$ 大小不超过 $100$ ，因此可以考虑缩点。
缩点后的图为 $DAG$ ，且根据 $tarjan$ 算法的退栈顺序， $SCC$ 编号满足从后往前的拓扑序，可以直接进行递推计算。
而在每个 $SCC$ 内的所有点可以使用高斯消元。
这样理论最大时间复杂度为 $O(100*100^3)$ 。
状态定义： $f[i]$ 表示从 $i$ 走到终点的期望时间
目标： $f[S]$
边界： $f[T]=0$
转移方程： $f[i]=\sum_{存在i到j的有向边}(\frac{f[j]}{out[i]})+1$
PS：
为了在每个 $SCC$ 内进行高斯消元，点的编号需要离散化。
在计算高斯消元系数矩阵时，如果 $f[j]$ 已知（ $i,j$ 属于不同的 $SCC$ ），则作为常数加入 $B$ 矩阵。
否则（ $j$ 所在的 $SCC$ 中的所有点已计算完毕）作为变量系数加入 $A$ 矩阵。
题解链接 https://www.luogu.com.cn/problemnew/solution/P6030
代码如下

/*

*/
#define method_1
#ifdef method_1
/*

*/
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<set>
#include<map>
#include<queue>
#include<stack>
#include<vector>
#include<cstring>
#include<cstdlib>
#include<iomanip>
#include<ctime>
#include<string>
#include<bitset>
#define D(x) cout<<#x<<" = "<<x<<"  "
#define E cout<<endl
using namespace std;
typedef long long ll;
typedef pair<int,int>pii;
const int maxn=1e4+5;
const int maxs=100+5;
const int maxm=1e6+5;
const double INF=1e18;
const double eps=1e-6;
int head[maxn],tot=1;
struct node{
    int from,to;
}edge[maxm];
void add(int from,int to){
    edge[++tot].from=head[from],head[from]=tot,edge[tot].to=to;
}
int num=0,cnt=0,dfn[maxn],low[maxn],top=0,st[maxn],inx[maxn],c[maxn];
vector<int>scc[maxn];
int id[maxn]; //为每个点重新编号 id[i]表示i点在所处的scc中的编号
void tarjan(int x){
    dfn[x]=low[x]=++num;
    inx[x]=1;
    st[++top]=x;
    for(int i=head[x];i;i=edge[i].from){
        int y=edge[i].to;
        if(!dfn[y]){
            tarjan(y);
            low[x]=min(low[x],low[y]);
        }
        else if(inx[y]){
            low[x]=min(low[x],dfn[y]);
        }
    }
    if(low[x]==dfn[x]){
        int y;
        cnt++;
        do{
            y=st[top];
            top--;
            inx[y]=0;
            c[y]=cnt;
            scc[cnt].push_back(y);
            id[y]=scc[cnt].size();
        }while(x!=y);
    }
}
double f[maxn];
double A[maxs][maxs],B[maxs],X[maxs];
int out[maxn];
void GuassInit(){
    memset(A,0,sizeof(A));
    memset(B,0,sizeof(B));
    memset(X,0,sizeof(X));
}
void Guass(int n){
    for(int i=1;i<=n;i++){
        int p=i;
        for(int j=i+1;j<=n;j++){
            if(fabs(A[j][i])>fabs(A[p][i])) p=j;
        }
        if(p!=i){
            //swap(A[p],A[i]); //在本地不允许这种交换
            for(int j=1;j<=n;j++) swap(A[p][j],A[i][j]);
            swap(B[p],B[i]);
        }
        if(fabs(A[i][i])<eps) return; //存在自由元
        for(int j=i+1;j<=n;j++){
            double d=A[j][i]/A[i][i];
            B[j]-=d*B[i];
            for(int k=i;k<=n;k++) A[j][k]-=d*A[i][k];
        }
    }
    for(int i=n;i>=1;i--){
        for(int j=i+1;j<=n;j++) B[i]-=X[j]*A[i][j];
        X[i]=B[i]/A[i][i];
    }
}
int n,m,s,t;
void solve(){
    tarjan(s);
    if(!dfn[t]){ //t不可达
        printf("INF"); 
        return;
    }
    //按照tarjan退栈顺序，1~cnt就是scc拓扑序从后往前的顺序
    for(int i=1;i<=cnt;i++){
        GuassInit();
        for(int j=0;j<scc[i].size();j++){
            int x=scc[i][j];
            if(x==t){
                A[id[x]][id[x]]=1.0; //f[t]=0
                continue;
            }
            A[id[x]][id[x]]=1.0;
            B[id[x]]=1.0;
            for(int k=head[x];k;k=edge[k].from){
                int y=edge[k].to;
                if(c[x]==c[y]) A[id[x]][id[y]]-=1.0/out[x]; //处在同一个强连通分量里
                else B[id[x]]+=1.0/out[x]*f[y]; //y已经计算出结果
            }
            if(head[x]==0) B[id[x]]=INF;
        }
        Guass(scc[i].size());
        for(int j=0;j<scc[i].size();j++){
            f[scc[i][j]]=X[id[scc[i][j]]];
            if(f[scc[i][j]]>1e10) f[scc[i][j]]=INF;
        }
    }
    if(f[s]>1e10){
        printf("INF"); //存在岔路
        return;
    }
    printf("%.3lf",f[s]);
}
int main() {
//  ios::sync_with_stdio(false);
    //freopen("BZOJ2707.in","r",stdin);
    scanf("%d%d%d%d",&n,&m,&s,&t);
    int from,to;
    for(int i=1;i<=m;i++){
        scanf("%d%d",&from,&to);
        add(from,to);
        out[from]++;
    }
    solve();
    return 0;
}
#endif
#ifdef method_2
/*

*/

#endif
#ifdef method_3
/*

*/

#endif

BZOJ2337
题意
$n(n<=100)$ 个点的无向连通图，边有边权，每次移动等概率的移动向相邻点。
求从 $1$ 到 $n$ 路径的期望异或和。
题解
首先，考虑到异或不进位的特性，可以将每一位分开考虑dp，再在最后进行合并。
状态定义： $f[i]$ 表示从点 $i$ 到 $n$ 的路径期望异或和中，当前位为 $1$ 的概率（则 $1-f[i]$ 为当前位为 $0$ 的概率）
目标： $f[1]$
边界： $f[n]=0$
转移方程： $f[i]=\frac{1}{deg[i]}(\sum_{w(i,j)=0}f[j]+\sum_{w(i,j)=1}(1-f[j]))$
转移方程解释：若 $j$ 到 $n$ 的异或和为 $0$ ，那么 $i$ 到 $j$ 到 $n$ 的异或和即取决于 $<i,j>$ 在这一位上的边权。
PS：根据期望方程设置高斯消元矩阵的方式
考虑关于 $i$ 的方程，将所有变量移到方程左侧，方程右侧的常量即为 $B[i]$ 。
左侧各个变量 $j1,j2,j...k$ 的系数即为 $A[i,j]$ 的值。
于是上述方程可变形为 $-deg[i]*f[i]+\sum_{w(i,j)=0}f[j]-\sum_{w(i,j)=1}f[j]=\sum_{w(i,j)=1}1$
那么 $A[i,i]$ 就是 $-deg[i]$ ， $A[i][j]$ 就 $+1$ 或者 $-1$ （在本题中，这取决于 $w(i,j)$ 的值）
题解链接 https://www.luogu.com.cn/problemnew/solution/P3211
代码如下

/*

*/
#define method_1
#ifdef method_1
/*

*/
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<set>
#include<map>
#include<queue>
#include<stack>
#include<vector>
#include<cstring>
#include<cstdlib>
#include<iomanip>
#include<ctime>
#include<string>
#include<bitset>
#define D(x) cout<<#x<<" = "<<x<<"  "
#define E cout<<endl
using namespace std;
typedef long long ll;
typedef pair<int,int>pii;
const int maxn=100+5;
const int maxm=1e4+5;
const int maxlog=30;
const int INF=0x3f3f3f3f;
const double eps=1e-6;
int n,m;
struct node{
    int from,to,v;
}edge[maxm<<1];
int head[maxn],tot=1;
void add(int from,int to,int v){
    edge[++tot].from=head[from],head[from]=tot,edge[tot].to=to,edge[tot].v=v;
}
int deg[maxn];
double A[maxn][maxn],B[maxn],X[maxn];
void GuassInit(){
    memset(A,0,sizeof(A));
    memset(B,0,sizeof(B));
    memset(X,0,sizeof(X));
}
void Guass(){
    for(int i=1;i<=n;i++){
        int p=i;
        for(int j=i+1;j<=n;j++){
            if(fabs(A[j][i])>fabs(A[p][i])) p=j;
        }
        if(p!=i){
            //swap(A[p],A[i]); //在本地不允许这种交换
            for(int j=1;j<=n;j++) swap(A[p][j],A[i][j]);
            swap(B[p],B[i]);
        }
        if(fabs(A[i][i])<eps) return; //存在自由元
        for(int j=i+1;j<=n;j++){
            double d=A[j][i]/A[i][i];
            B[j]-=d*B[i];
            for(int k=i;k<=n;k++) A[j][k]-=d*A[i][k];
        }
    }
    for(int i=n;i>=1;i--){
        for(int j=i+1;j<=n;j++) B[i]-=X[j]*A[i][j];
        X[i]=B[i]/A[i][i];
    }
}
void solve(){
    double ans=0.0;
    for(int k=maxlog;k>=0;k--){
        GuassInit();
        for(int x=1;x<=n;x++){
            if(x==n){
                A[n][n]=1.0; 
                B[n]=0.0; //1*f[n]=0
                continue;
            }
            A[x][x]=-deg[x];
            B[x]=0.0;
            for(int i=head[x];i;i=edge[i].from){
                int y=edge[i].to,v=edge[i].v;
                int valV=(v>>k)&1;
                if(valV){
                    A[x][y]-=1.0;
                    B[x]-=1.0;
                }
                else{
                    A[x][y]+=1.0;
                }
            }
        }
        Guass();
        ans+=(1<<k)*X[1];
    }
    printf("%.3lf",ans);
}
int main() {
//  ios::sync_with_stdio(false);
    //freopen("BZOJ2337.in","r",stdin);
    scanf("%d%d",&n,&m);
    int from,to,v;
    for(int i=1;i<=m;i++){
        scanf("%d%d%d",&from,&to,&v);
        if(from==to){ //自环算一次
            deg[from]++;
            add(from,to,v);
        }
        else{
            deg[from]++,deg[to]++;
            add(from,to,v),add(to,from,v);
        }
    }
    solve();
    return 0;
}
#endif
#ifdef method_2
/*

*/

#endif
#ifdef method_3
/*

*/

#endif

DP训练——概率期望DP

概率期望DP

猜你喜欢

热点阅读