概率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(fi1,j,0+wci1,ci,fi1,j,1+wdi1,cipi1+wci1,ci(1pi1))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(fi1,j,0+wci1,ci(1pi)+wdi1,cipi,(fi1,j,1+wdi1,cipi1+wci1,ci(1pi1))(1pi)+(fi1,j,1+wdi1,dipi1+wci1,di(1pi1))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]100dp[1]+(1kill[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了不得不说又是一道经典的“天坑”题目。

首先把点扩展成2n22*n-2个,即0,1,2,...,n2,n1,n2,...,2,10,1,2,...,n-2,n-1,n-2,...,2,1。这样一来,可以始终看成是从左往右移动,超出右边界则对2n22*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了,想的是对每条边都设一个未知数,然后用这些未知数去表达节点1dp值。可是这么做只有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;
}