关于概率dp的个人理解与总结

原文来自：/article/3708792.html

首先，概率dp主要解决的是关于概率问题和期望问题的求解。

难点和普通dp一样在于dp[i][j][k]的维数控制和含义，其实就是转移方程的构建。

然后一般地，求概率是正推、求期望是逆推。（开始的很多状态不可能发生概率为0，最后的状态出口期望为0）

对于求概率当前点的概率是由前面能到达当前点的点乘上到达当前点的概率得到的。

也就是dp[i]=Σ(dp[j]*p[j][i]) i是当前点、j是前面的点。

对于求期望当前点的期望是由当前点所能到达的点得到的。（注意下，对应的概率是当前点到达之后点的概率，因为你要到达后面的点后面的点才能传递给你期望！）

也就是E[i]=Σ((E[j]+k)*p[i][j]) i是当前点、j是当前点所能到达的点，k是所需的期望值。

然后就是对于期望的话，如果成环的话数据范围小的话可以用高斯消元解决（之前的博文提到），如果范围大就要推导公式了（后面的题目有提到）。

poj 2096 Collecint Bugs

dp求期望入门题。

题意，思路，公式引用别人的：/content/3596333.html

dp求期望
逆着递推求解
题意：（题意看题目确实比较难道，n和s都要找半天才能找到）
一个软件有s个子系统，会产生n种bug
某人一天发现一个bug,这个bug属于一个子系统，属于一个分类
每个bug属于某个子系统的概率是1/s,属于某种分类的概率是1/n
问发现n种bug,每个子系统都发现bug的天数的期望。

求解：
dp[i][j]表示已经找到i种bug,j个系统的bug，达到目标状态的天数的期望
dp
[s]=0;要求的答案是dp[0][0];
dp[i][j]可以转化成以下四种状态:
dp[i][j],发现一个bug属于已经有的i个分类和j个系统。概率为(i/n)*(j/s);
dp[i][j+1],发现一个bug属于已有的分类，不属于已有的系统.概率为 (i/n)*(1-j/s);
dp[i+1][j],发现一个bug属于已有的系统，不属于已有的分类,概率为 (1-i/n)*(j/s);
dp[i+1][j+1],发现一个bug不属于已有的系统，不属于已有的分类,概率为 (1-i/n)*(1-j/s);
整理便得到转移方程

代码：

[cpp] view
plain copy

#include"cstdlib"

#include"cstdio"

#include"cstring"

#include"cmath"

#include"queue"

#include"algorithm"

#include"iostream"

#define eps 1e-5

using namespace std;

double dp[1020][1020];

int main()

{

int n,s;

while(cin>>n>>s)

{

memset(dp,0,sizeof(dp));

for(int i=n; i>=0; i--)

{

for(int j=s; j>=0; j--)

{

if(i==n&&j==s) continue;

double tep=1.0;

if(i<n) tep+=(n*1.0-i)*j/(n*s)*dp[i+1][j];

if(j<s) tep+=(s*1.0-j)*i/(n*s)*dp[i][j+1];

if(i<n&&j<s) tep+=(n*1.0-i)*(s*1.0-j)/(n*s)*dp[i+1][j+1];

dp[i][j]=tep/(1-(i*j*1.0)/(s*n));

}

}

printf("%.4f\n",dp[0][0]);

}

return 0;

}

类似的简单题：

hdu 3853 注意停留原地的情况，应该跳过不处理。还有一点走一步的期望是2。

hdu 4405 处理一下连续移动的情况

SGU 495 用概率求期望

CF 148D 考虑下边界就OK了，方程很好推

zoj 3640 Help Me Escape

题意：

一只吸血鬼，有n条路给他走，每次他随机走一条路，

每条路有个限制，如果当时这个吸血鬼的攻击力大于等于某个值，那么就会花费t天逃出去。

否则，花费1天的时间，并且攻击力增加，问他逃出去的期望天数。

思路：

方程很好像 dp[i] 有i点战斗力逃出去的期望

那么如果 x>c[i]时 dp[i]=(1.0/n)*Σt[i]

否则 dp[i]=(1.0/n)*Σ(dp[i+c[i]]+1)

因为很多点会被重复用，所以写成记忆化搜索。

还有一个要注意的，t[i]是向下取整的。

代码：

[cpp] view
plain copy

#include"cstdlib"

#include"cstdio"

#include"cstring"

#include"cmath"

#include"queue"

#include"algorithm"

#include"iostream"

#define eps 1e-12

using namespace std;

#define N 10020

double dp
;

int c[123];

int t[123];

int MAX;

int n,f;

double dfs(int x)

{

if(x>MAX) x=MAX; //小优化 超过最大值 就为最大值

if(fabs(dp[x]+1)>eps) return dp[x];

double tep=0.0;

for(int i=1;i<=n;i++)

{

if(x>c[i]) tep+=1.0/n*t[i];

else tep+=1.0/n*(dfs(x+c[i])+1);

}

return dp[x]=tep;

}

int main()

{

while(cin>>n>>f)

{

int i;

MAX=-1;

memset(dp,-1,sizeof(dp));

for(i=1;i<=n;i++)

{

scanf("%d",&c[i]);

MAX=max(MAX,c[i]+1);

t[i]=(int)((1+sqrt(5.0))/2*c[i]*c[i]); //这里要注意了 t[i]是一个向下取整的！

}

dp[f]=dfs(f);

printf("%.3f\n",dp[f]);

}

return 0;

}

poj 3744 Scout YYF I

其实概率方程很简单，就是有地雷的时候那个点是不能得到其他点的概率要是0.

然后就是分段的思想了，按区间分步求，起点是1，终点MAX+1。

然后需要注意的是：

1、有可能1是雷，你直接爆炸~ bomb :) ~

2、范围很大，需要用矩阵快速幂优化

3、雷不是有序的，要记得排序

代码：

[cpp] view
plain copy

#include"cstdlib"

#include"cstdio"

#include"cstring"

#include"cmath"

#include"queue"

#include"algorithm"

#include"iostream"

#define eps 1e-12

using namespace std;

#define N 10020

struct matrix

{

double mat[3][3];

};

matrix matmul(matrix a,matrix b,int n)

{

int i,j,k;

matrix c;

memset(c.mat,0,sizeof(c.mat));

for(i=0; i<n; i++) for(j=0; j<n; j++) for(k=0; k<n; k++) c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];

return c;

}

matrix matpow(matrix a,int k,int n)

{

matrix b;

int i;

memset(b.mat,0,sizeof(b.mat));

for(i=0; i<n; i++) b.mat[i][i]=1;

while(k)

{

if(k&1) b=matmul(a,b,n);

a=matmul(a,a,n);

k>>=1;

}

return b;

}

int main()

{

int n;

double p;

while(scanf("%d%lf",&n,&p)!=-1)

{

double q=1-p;

int i;

int mine[23];

for(i=0; i<n; i++) scanf("%d",&mine[i]);

sort(mine,mine+n);

int k=0;

matrix ans,a,b,c;

memset(a.mat,0,sizeof(a.mat));

memset(b.mat,0,sizeof(b.mat));

memset(c.mat,0,sizeof(c.mat));

if(mine[k]==1)

{

puts("0.0000000");

continue;

}

else if(mine[k]==2)

{

a.mat[0][0]=0;

a.mat[0][1]=1;

k++;

}

else

{

a.mat[0][0]=p;

a.mat[0][1]=1;

}

int yuan=2;

b.mat[0][0]=p;

b.mat[1][0]=q;

b.mat[0][1]=1;

c=b;

for(i=k;i<n;i++)

{

c=b;

ans=matmul(a,matpow(c,mine[i]-yuan-1,2),2);

a.mat[0][0]=0;

a.mat[0][1]=ans.mat[0][0];

yuan=mine[i];

}

ans=matmul(a,matpow(c,mine[n-1]+1-yuan,2),2);

printf("%.7f\n",ans.mat[0][0]);

}

return 0;

}

还有几道中等的题目：

poj 2151 dp[i][j][k]第i个人在前j题中做出了k题的概率 最后所有人至少作出一道题的概率减去所有人都做1~N-1题的概率 即为所求

poj 3071 全概率公式 巧妙确定范围的方法

[cpp] view
plain copy

for(i=0; i<(1<<n); i++) dp[0][i]=1;

for(i=1; i<=n; i++)

{

for(j=0; j<(1<<n); j++)

{

int tep=j/(1<<(i-1));

tep^=1;

int s=tep*(1<<(i-1));

int e=tep*(1<<(i-1))+(1<<(i-1));

for(k=s; k<e; k++)

dp[i][j]+=dp[i-1][j]*dp[i-1][k]*p[j][k];

}

}

以下几题都比较难了：

hdu 4089 Activation

题意：有n个人排队等着在官网上激活游戏。Tomato排在第m个。
对于队列中的第一个人。有一下情况：
1、激活失败，留在队列中等待下一次激活（概率为p1)
2、失去连接，出队列，然后排在队列的最后（概率为p2）
3、激活成功，离开队列（概率为p3）
4、服务器瘫痪，服务器停止激活，所有人都无法激活了。
求服务器瘫痪时Tomato在队列中的位置<=k的概率

解析：
概率DP；
设dp[i][j]表示i个人排队,Tomato排在第j个位置，达到目标状态的概率(j<=i)
dp
[m]就是所求
j==1:    dp[i][1]=p1*dp[i][1]+p2*dp[i][i]+p4;
2<=j<=k: dp[i][j]=p1*dp[i][j]+p2*dp[i][j-1]+p3*dp[i-1][j-1]+p4;
k<j<=i:  dp[i][j]=p1*dp[i][j]+p2*dp[i][j-1]+p3*dp[i-1][j-1];
化简：
j==1:    dp[i][1]=p*dp[i][i]+p41;
2<=j<=k: dp[i][j]=p*dp[i][j-1]+p31*dp[i-1][j-1]+p41;
k<j<=i:  dp[i][j]=p*dp[i][j-1]+p31*dp[i-1][j-1];

其中:
p=p2/(1-p1);
p31=p3/(1-p1)
p41=p4/(1-p1)

可以循环i=1->n 递推求解dp[i].在求解dp[i]的时候dp[i-1]就相当于常数了。
在求解dp[i][1~i]时等到下列i个方程
j==1:   dp[i][1]=p*dp[i][i]+c[1];
2<=j<=k:dp[i][j]=p*dp[i][j-1]+c[j];
k<j=i:  dp[i][j]=p*dp[i][j]+c[j];
其中c[j]都是常数了。上述方程可以解出dp[i]了。
首先是迭代得到 dp[i][i].然后再代入就可以得到所有的dp[i]了。

注意特判一种情况。就是p4<eps时候，就不会崩溃了，应该直接输出0

代码：

[cpp] view
plain copy

#include"cstdlib"

#include"cstdio"

#include"cstring"

#include"cmath"

#include"queue"

#include"algorithm"

#include"iostream"

#define eps 1e-12

using namespace std;

double dp[2022][2022];

int main()

{

int n,m,k;

double p1,p2,p3,p4;

while(cin>>n>>m>>k>>p1>>p2>>p3>>p4)

{

memset(dp,0,sizeof(dp));

int i,j;

double p22,p33,p44;

p22=p2/(1-p1);

p33=p3/(1-p1);

p44=p4/(1-p1);

if(fabs(p4)<=eps)

{

puts("0.00000");

continue;

}

dp[1][1]=p44/(1-p22);

for(i=2;i<=n;i++)

{

double p=1.0,tep=0.0;

for(j=i;j>=1;j--,p*=p22) //中间就是迭代求出其中一个，因为五个方程五个未知数，形式相等。

{

if(j==1) tep+=p*p44;

else if(j>=2&&j<=k) tep+=p*(p33*dp[i-1][j-1]+p44);

else tep+=p*p33*dp[i-1][j-1];

}

dp[i][i]=tep/(1-p);

for(j=1;j<i;j++)

{

if(j==1) dp[i][j]=p22*dp[i][i]+p44;

else if(j>=2&&j<=k) dp[i][j]=p22*dp[i][j-1]+p33*dp[i-1][j-1]+p44;

else dp[i][j]=p22*dp[i][j-1]+p33*dp[i-1][j-1];

}

}

printf("%.5f\n",dp
[m]);

}

return 0;

}

zoj 3329 One Person Game

题意：有三个骰子，分别有k1,k2,k3个面，人物开始分数为零。

每次掷骰子，如果三个面分别是a,b,c的话分数归零，否则分数加上三个骰子的点数和。

分数超过n则游戏结束。求游戏结束的期望掷骰子次数。

思路： 设投出a,b,c 的概率为p0，投出k分的概率为p[k]

则对于每个分数 dp[i]=p0*dp[0]+Σ(dp[i+k]*p[k])+1 ①

又设 dp[i]=A[i]dp[0]+B[i] ②，则dp[i+k]=A[i+k]dp[0]+B[i+k]

代入① 得： dp[i]=p0*dp[0]+Σ((A[i+k]dp[0]+B[i+k])*p[k])+1

化简得： dp[i]=(p0+Σ(A[i+k]*p[k]))*dp[0] + Σ(B[i+k]*p[k]) +1

对比②得： A[i]=Σ(A[i+k]*p[k])+p0

B[i]=Σ(B[i+k]*p[k]) +1

当i+k>n时A[i+k]=B[i+k]=0

所以dp[0]=B[0] / (1-A[0])可求出

代码：

[cpp] view
plain copy

#include"cstdlib"

#include"cstdio"

#include"cstring"

#include"cmath"

#include"queue"

#include"algorithm"

#include"iostream"

#define eps 1e-12

using namespace std;

int main()

{

int t;

cin>>t;

while(t--)

{

int n,k[4],a,b,c;

int i,j,l;

scanf("%d%d%d%d%d%d%d",&n,&k[1],&k[2],&k[3],&a,&b,&c);

double p[22];

memset(p,0,sizeof(p));

for(i=1;i<=k[1];i++) for(j=1;j<=k[2];j++) for(l=1;l<=k[3];l++) p[i+j+l]+=1.0/(k[1]*k[2]*k[3]); //求出掷出k点的概率

p[a+b+c]-=1.0/(k[1]*k[2]*k[3]); //单独考虑下分别a,b,c点的情况

p[0]=1.0/(k[1]*k[2]*k[3]);

double A[600],B[600];

memset(A,0,sizeof(A));

memset(B,0,sizeof(B));

for(i=n;i>=0;i--) //递推求解

{

A[i]=p[0];

B[i]=1;

if(i==n) continue;

for(j=1;j<=k[1]+k[2]+k[3];j++)

{

A[i]+=p[j]*A[i+j];

B[i]+=p[j]*B[i+j];

}

}

printf("%.15f\n",B[0]/(1-A[0]));

}

return 0;

}

总结下这类概率DP:

既DP[i]可能由DP[i+k]和DP[i+j]需要求的比如DP[0]决定

相当于概率一直递推下去会回到原点

比如

(1):DP[i]=a*DP[i+k]+b*DP[0]+d*DP[i+j]+c;

但是DP[i+k]和DP[0]都是未知

这时候根据DP[i]的方程式假设一个方程式：

比如:

(2):DP[i]=A[i]*DP[i+k]+B[i]*DP[0]+C[i];

因为要求DP[0],所以当i=0的时候但是A[0],B[0],C[0]未知

对比(1)和(2)的差别

这时候对比(1)和(2)发现两者之间的差别在于DP[i+j]

所以根据(2)求DP[i+j]然后代入(1)消除然后对比(2)就可以得到A[i],B[i],C[i]

然后视具体情况根据A[i],B[i],C[i]求得A[0],B[0],C[0]继而求DP[0]

然后下面那题可以好好感受一下这个！

hdu 4035 Maze

经典的概率dp，关于成环的概率dp的求解，公式的推导。

题意，思路，公式引用别人的：/content/3596333.html

dp求期望的题。
题意：
有n个房间，由n-1条隧道连通起来，实际上就形成了一棵树，
从结点1出发，开始走，在每个结点i都有3种可能：
1.被杀死，回到结点1处（概率为ki）
2.找到出口，走出迷宫 （概率为ei）
3.和该点相连有m条边，随机走一条

求：走出迷宫所要走的边数的期望值。
设 E[i]表示在结点i处，要走出迷宫所要走的边数的期望。E[1]即为所求。
叶子结点：
E[i] = ki*E[1] + ei*0 + (1-ki-ei)*(E[father[i]] + 1);
= ki*E[1] + (1-ki-ei)*E[father[i]] + (1-ki-ei);
非叶子结点：（m为与结点相连的边数）
E[i] = ki*E[1] + ei*0 + (1-ki-ei)/m*( E[father[i]]+1 + ∑( E[child[i]]+1 ) );
= ki*E[1] + (1-ki-ei)/m*E[father[i]] + (1-ki-ei)/m*∑(E[child[i]]) + (1-ki-ei);
设对每个结点：E[i] = Ai*E[1] + Bi*E[father[i]] + Ci;
对于非叶子结点i，设j为i的孩子结点，则
∑(E[child[i]]) = ∑E[j]
= ∑(Aj*E[1] + Bj*E[father[j]] + Cj)
= ∑(Aj*E[1] + Bj*E[i] + Cj)
带入上面的式子得
(1 - (1-ki-ei)/m*∑Bj)*E[i] = (ki+(1-ki-ei)/m*∑Aj)*E[1] + (1-ki-ei)/m*E[father[i]] + (1-ki-ei) + (1-ki-ei)/m*∑Cj;
由此可得
Ai =        (ki+(1-ki-ei)/m*∑Aj)   / (1 - (1-ki-ei)/m*∑Bj);
Bi =        (1-ki-ei)/m            / (1 - (1-ki-ei)/m*∑Bj);
Ci = ( (1-ki-ei)+(1-ki-ei)/m*∑Cj ) / (1 - (1-ki-ei)/m*∑Bj);
对于叶子结点
Ai = ki;
Bi = 1 - ki - ei;
Ci = 1 - ki - ei;
从叶子结点开始，直到算出 A1,B1,C1;
E[1] = A1*E[1] + B1*0 + C1;
所以
E[1] = C1 / (1 - A1);
若 A1趋近于1则无解...

代码：

[cpp] view
plain copy

#include"cstdlib"

#include"cstdio"

#include"cstring"

#include"cmath"

#include"queue"

#include"algorithm"

#include"iostream"

#define eps 1e-12 //注意开1e-8会WA

#include"vector"

#define N 10010

using namespace std;

struct node

{

int v,w;

node(int vv)

{

v=vv;

}

};

vector<node>edge
;

double A
,B
,C
,Asum
,Bsum
,Csum
;

double p
,e
;

void dfs(int u,int f) //最简单的对树的搜索

{

int n=edge[u].size();

if(n==1&&u!=1) //如果是叶子节点

{

A[u]=Asum[u]=p[u];

B[u]=Bsum[u]=1-p[u]-e[u];

C[u]=Csum[u]=1-p[u]-e[u];

return ;

}

Asum[u]=Bsum[u]=Csum[u]=0;

for(int i=0;i<n;i++) //不是叶子节点的 按公式求值

{

if(edge[u][i].v==f) continue;

dfs(edge[u][i].v,u);

Asum[u]+=A[edge[u][i].v];

Bsum[u]+=B[edge[u][i].v];

Csum[u]+=C[edge[u][i].v];

A[u]=(p[u]+(1-p[u]-e[u])/n*Asum[u]) /(1-(1-p[u]-e[u])/n*Bsum[u]);

B[u]=((1-p[u]-e[u])/n) /(1-(1-p[u]-e[u])/n*Bsum[u]);

C[u]=((1-p[u]-e[u])+(1-p[u]-e[u])/n*Csum[u]) /(1-(1-p[u]-e[u])/n*Bsum[u]);

}

return ;

}

int main()

{

int t,cas=1;

cin>>t;

while(t--)

{

int n,i;

cin>>n;

for(i=0;i<=n;i++)

edge[i].clear();

for(i=1;i<=n-1;i++)

{

int a,b;

scanf("%d%d",&a,&b);

edge[a].push_back(node(b));

edge[b].push_back(node(a));

}

for(i=1;i<=n;i++) //注意读进来的数是百分比

{

scanf("%lf%lf",&p[i],&e[i]);

p[i]*=0.01;

e[i]*=0.01;

}

dfs(1,1);

printf("Case %d: ",cas++);

if(fabs(1-A[1])<=eps) //不断被杀死 逃不出去 T_T

{

puts("impossible");

continue;

}

printf("%.6f\n",C[1]/(1-A[1]));

}

return 0;

}