關於概率dp的個人理解與總結
首先,概率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求期望入門題。
dp求期望 逆著遞推求解 題意:(題意看題目確實比較難道,n和s都要找半天才能找到) 一個軟體有s個子系統,會產生n種bug 某人一天發現一個bug,這個bug屬於一個子系統,屬於一個分類 每個bug屬於某個子系統的概率是1/s,屬於某種分類的概率是1/n 問發現n種bug,每個子系統都發現bug的天數的期望。 求解: dp[i][j]表示已經找到i種bug,j個系統的bug,達到目標狀態的天數的期望 dp[n][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); 整理便得到轉移方程程式碼:
#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]是向下取整的。
程式碼:
#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[N];
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、雷不是有序的,要記得排序
程式碼:
#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 全概率公式 巧妙確定範圍的方法
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[n][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
程式碼:
#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[n][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])可求出
程式碼:
#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的求解,公式的推導。
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則無解...
程式碼:
#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[N];
double A[N],B[N],C[N],Asum[N],Bsum[N],Csum[N];
double p[N],e[N];
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;
}