概率DP

期望DP

「NOIP2016」换教室

又是一道当年没做过现在来补票的NOIP题目

定义DP状态为fi,j,0/1f_{i, j, 0 / 1}，其中i表示在第i个时间段，连同当前时间段总共用了j次换教室的机会，在这个时间段换（1）或者不换教室（0）的最小期望路程和。答案就是max⁡{fn,i,0,fn,i,1},i∈[0,m]\max \left\{f_{n, i, 0}, f_{n, i, 1}\right\}, i \in[0, m]

  • 如果当前阶段不换教室
    fi,j,0=min⁡(fi−1,j,0+wci−1,ci,fi−1,j,1+wdi−1,ci⋅pi−1+wci−1,ci⋅(1−pi−1))f_{i, j, 0}=\min \left(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\left(1-p_{i-1}\right)\right)
  • 如果当前阶段换教室
    fi,j,0=min⁡(fi−1,j,0+wci−1,ci⋅(1−pi)+wdi−1,ci⋅pi,(fi−1,j,1+wdi−1,ci⋅pi−1+wci−1,ci⋅(1−pi−1))⋅(1−pi)+(fi−1,j,1+wdi−1,di⋅pi−1+wci−1,di⋅(1−pi−1))⋅pi)f_{i, j, 0}=\min (f_{i-1, j, 0}+w_{c_{i-1}, c_{i}}\cdot (1-p_{i})+w_{d_{i-1}, c_{i}}\cdot p_{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}))\cdot (1-p_i) + (f_{i-1, j, 1}+w_{d_{i-1}, d_{i}} \cdot p_{i-1}+w_{c_{i-1}, d_{i}} \cdot(1-p_{i-1}))\cdot p_i)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#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;
}

HDU4035 Maze

一拿到题，心想这咋递推啊，不是有后效性吗，那看来得搞个高斯消元了。一看数据范围n<=1e4n<=1e4，高斯消元绝对炸。然后一看是树结构，于是开始考虑如何利用树的性质消元从而规避后效性。

首先推出式子：
dp[u]=kill[u]100⋅dp[1]+(1−kill[u]+esc[u]100)⋅(∑dp[v]+1dv+dp[fa]+1dv)dp[u]=\frac{kill[u]}{100}\cdot dp[1]+(1-\frac{kill[u]+esc[u]}{100})\cdot(\sum \frac{dp[v]+1}{dv} +\frac{dp[fa]+1}{dv})

对于一个节点，它的dp值由三个部分的贡献构成，即dp[1],dp[fa]+1,dp[v]+1dp[1],dp[fa]+1,dp[v]+1。我们可以把一个dp值扩展成一个结构体，结构体中有三个元素，分别存储它的常数项，父节点的贡献（父节点的dp值占的比例），根节点的贡献（根节点的dp值占的比例）。从而可以展开上面的式子进行化简，从而避免了后效性。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#include<bits/stdc++.h>
using namespace std;
typedef double db;
const int maxn=1e4+5;
const double eps=0.00000000001;
int n;
vector<int>e[maxn];
int kill[maxn],esc[maxn];
struct Node{
	double rt,fa,cons;//ratio : dp[1]  dp[father]  const
	Node(){
		rt=fa=cons=0.0;
	}
}dp[maxn];
int sz[maxn];
void dfs(int u,int p){
	sz[u]=1;
	double idx=0.0;
	double temp=1.0-(kill[u]+esc[u])/100.0;
	for(auto v:e[u]){
		if(v==p)continue;
		sz[u]++;
		dfs(v,u);
		idx+=dp[v].fa;
	}	
	if(u==1)return;
	idx=1.0-idx*temp/sz[u];
	for(auto v:e[u]){
		if(v==p)continue;
		dp[u].cons+=temp/sz[u]*(dp[v].cons+1.0)/idx;
		dp[u].rt+=temp/idx/sz[u]*dp[v].rt;
	}
	dp[u].fa=temp/sz[u]/idx;
	dp[u].cons+=temp/sz[u]/idx;
	dp[u].rt+=kill[u]/100.0/idx;
}
int main(){
	int T,t=0;scanf("%d",&T);
	while(++t<=T){
		scanf("%d",&n);
		for(int i=1;i<n;i++){
			int u,v;
			scanf("%d%d",&u,&v);
			e[u].push_back(v);
			e[v].push_back(u);
		}
		for(int i=1;i<=n;i++){
			scanf("%d%d",&kill[i],&esc[i]);
		}
		dfs(1,0);sz[1]--;
		double ans=0.0,idx=1.0;
		for(auto v:e[1]){
			idx-=(dp[v].rt+dp[v].fa)/sz[1];
			ans+=(dp[v].cons+1.0)/sz[1];
		}
		if(idx<=eps)printf("Case %d: impossible\n",t);
		else printf("Case %d: %.6lf\n",t,ans/idx);
		for(int i=1;i<=n;i++)e[i].clear();
		for(int i=1;i<=n;i++){
			dp[i].fa=dp[i].cons=dp[i].rt=0.0;
		}
	}


	return 0;
}

「SCOI2008」奖励关

这道题目并不难，但是具有一定的代表性，还是写一下总结吧。

该题目不能从前往后DP，我认为不仅仅是因为后效性（因为后效性其实可以通过给初始数组赋小于0的值再加上特判从而规避掉）。还有一个更重要的问题是，从前往后DP不能无重复地计算最终答案。

该问题的本质其实是对于总共nkn^k种不同的宝箱生成方式分别计算一个最优解，最后答案为各最优解的和，然后再除以nkn^k。

如果从前往后DP，不同的状态可能是由同一种宝箱生成方式得到的，显然不能直接加起来作为答案。

所以考虑从后往前DP，这样一来，任何被计入dp[1][0]dp[1][0]的贡献，都对应着一个独一无二的宝箱生成方式（可以用归纳法证明）。

引用在洛谷上看到的一句话：“概率是顺推，而期望需要逆推。”

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#include<bits/stdc++.h>
using namespace std;
int K,n,v[16],state[16];
double dp[105][1<<15];
inline double Max(register double a,register double b){
	return a>b?a:b;
}
int main(){
	scanf("%d%d",&K,&n);
	int tmp;
	for(int i=1;i<=n;++i){
		scanf("%d",&v[i]);
		while(~scanf("%d",&tmp)){
			if(!tmp)break;
			state[i]|=(1<<(tmp-1));
		}
	}
	for(int i=K;i>=1;--i){
		for(int j=0;j<(1<<n);++j){
			for(int k=1;k<=n;++k){
				if((j&state[k])==state[k])dp[i][j]+=Max(dp[i+1][j],dp[i+1][j|(1<<(k-1))]+v[k]);
				else dp[i][j]+=dp[i+1][j];
			}
			dp[i][j]/=n;
		}
	}
	printf("%.6lf",dp[1][0]);
	return 0;
}

有后效性DP

HDU Time Travel

WA了一下午终于AC了不得不说又是一道经典的“天坑”题目。

首先把点扩展成2∗n−22*n-2个，即0,1,2,...,n−2,n−1,n−2,...,2,10,1,2,...,n-2,n-1,n-2,...,2,1。这样一来，可以始终看成是从左往右移动，超出右边界则对2∗n−22*n-2取模即可。

然后重头戏来了。考虑题目中所说的impossible的情况。

通常的想法是，“这不很简单吗，无解的情况不就是高斯消元是系数为0的时候吗？”我们注意到，其实这种情况其实就对应着走进了死胡同，永远无法到达出口。只要有可能经过死胡同，那期望步数就直接成为无穷大，从而无解。但是只要从起点出发经过这些死胡同的概率为0 ，那么依然有解。

所以就算有时候高斯消元某个未知数的系数为0，仍然有可能算出起点对应的项的答案。

于是我们考虑BFS，标记一下从起点出发不可能经过的点，然后在高斯消元的时候进行排除。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#include<bits/stdc++.h>
using namespace std;
const int maxn=2e2+5;
const double eps=0.0000001;
int n,m,y,x,d,flag;
double e[maxn];
double a[maxn][maxn];
int vis[maxn];
void solve(){
	for(int i=0;i<2*n-2;i++){
		for(int j=i;j<2*n-2;j++){
			if(abs(a[j][i])>eps){
				for(int k=0;k<=2*n-2;k++)
					swap(a[j][k],a[i][k]);
				break;
			}
		}
		if(abs(a[i][i])<=eps){
			if(!vis[i])continue;
			flag=1;return;
		}
		double temp=a[i][i];
		for(int j=i;j<=2*n-2;j++)
			a[i][j]/=temp;
		for(int j=i+1;j<2*n-2;j++){
			double temp=a[j][i];
			for(int k=i;k<=2*n-2;k++){
				a[j][k]-=temp*a[i][k];
			}
		}
	}
	for(int i=n*2-3;i>=0;i--){
		if(!vis[i])continue;
		a[i][n*2-2]/=a[i][i];
		for(int j=i-1;j>=0;j--){
			a[j][n*2-2]-=a[i][n*2-2]*a[j][i];
		}
	}
}
void bfs(int u){
	queue<int>q;
	q.push(u);
	while(!q.empty()){
		int p=q.front();q.pop();
		for(int i=1;i<=m;i++){
			if(e[i]<eps)continue;
			int temp=(p+i)%(2*n-2);
			if(!vis[temp]){
				q.push(temp);
				vis[temp]=1;
			}
		}
	}
}
int main(){
	int T;scanf("%d",&T);
	while(T--){
		flag=0;
		scanf("%d%d%d%d%d",&n,&m,&y,&x,&d);
		for(int i=0;i<2*n-2;i++)vis[i]=0;
		if(d==1)x=2*(n-1)-x;
		for(int i=1;i<=m;i++){
			scanf("%lf",&e[i]);
			e[i]/=100;
		}		
		if(x==y){
			printf("0.00\n");continue;
		}
		vis[x]=1;bfs(x);
		for(int i=0;i<2*n-2;i++){
			for(int j=0;j<=2*n-2;j++)
				a[i][j]=0.0;
		}
		for(int i=0;i<2*n-2;i++){
			a[i][i]=1.0;
			if(i==y||i==(2*(n-1)-y))continue;
			for(int j=1;j<=m;j++){
				a[i][(j+i)%(2*n-2)]-=e[j];
				a[i][2*n-2]+=e[j]*j;
			}

		}
		solve();
		
		if(!flag)printf("%.2f\n",a[x][2*n-2]);
		else printf("Impossible !\n");
	}

	return 0;
}

「HNOI2013」游走

容易想到贪心，对期望访问次数最多的边赋尽量小的值。

考虑构建方程组。我刚开始做的时候就naive了，想的是对每条边都设一个未知数，然后用这些未知数去表达节点1的dp值。可是这么做只有60分，居然超时了。才反应过来这有可能是完全图，对每条边都设未知数然后高斯消元铁定会炸…

上面这种方法，dp[i]dp[i]定义为从i点出发还要多少代价到达n点。这次我们换一种方式，定义dp[i]dp[i]值为经过i点的期望次数，这样就能把系数矩阵控制在n2n^2范围内

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#include<bits/stdc++.h>
using namespace std;
const int maxn=5e2+5;
const double eps=1e-9;
int n,m;
int sz[maxn];
vector<int>e[maxn];
double a[maxn][maxn];
int e1[maxn*maxn],e2[maxn*maxn];
void solve(){
	for(int i=1;i<n;i++){
		if(abs(a[i][i])<eps){
			for(int j=i+1;j<n;j++){
				if(abs(a[j][i])>=eps){
					for(int k=1;k<=n;k++)
						swap(a[i][k],a[j][k]);
					break;
				}
			}
		}
		double temp=a[i][i];
		for(int j=i;j<=n;j++)a[i][j]/=temp;
		
		for(int j=i+1;j<n;j++){
			temp=a[j][i];
			for(int k=i;k<=n;k++){
				a[j][k]-=temp*a[i][k];
			}
		}
	}
	for(int i=n-1;i>=1;i--){
		for(int j=i-1;j>=1;j--){
			a[j][n]-=a[i][n]*a[j][i];
		}
	}
}
int main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=m;i++){
		int x,y;
		scanf("%d%d",&x,&y);
		e[x].push_back(y);
		e[y].push_back(x);
		e1[i]=x,e2[i]=y;
		sz[x]++;sz[y]++;
	}
	for(int i=1;i<n;i++){
		a[i][i]=1.0;
		for(auto j:e[i]){
			if(j==n)continue;
			a[i][j]-=1.0/sz[j];
		}
	}
	a[1][n]=1.0;
	solve();
	double ans=0.0;
	vector<double>vec;
	for(int i=1;i<=m;i++){
		double temp=0;
		if(e1[i]!=n)temp+=a[e1[i]][n]/sz[e1[i]];
		if(e2[i]!=n)temp+=a[e2[i]][n]/sz[e2[i]];
		vec.push_back(temp);
	}
	sort(vec.begin(),vec.end());
	for(int i=0,j=m;i<m;i++,j--){
		ans+=vec[i]*j;
	}
	printf("%.3lf\n",ans);
	return 0;
}