【专题总结】概率&期望DP

下面是对一系列求概率&期望问题的总结：

问题的形式一般是求一个事件发生的概率或者期望，而概率和期望的求解形式是有些区别的，一般方法是dp&记忆化搜索，特殊情况可能用到人肉消元
& 高斯消元。

如果是概率类的问题，那么求的方法一般是利用这个事件之前的概率，dp[i]的含义是i事件发生的概率(当然这里不一定是一维数组，写一维是为了简单起见，思路是一样的，下同)

大概的流程是 dp[i] = sigma(p[j] * dp[j])，j是可能转移到i的事件

即当前事件发生的概率是通过“可能转移到当前事件的事件”推过来的

初始状态的dp[i] = 1

而对于期望类的问题，求的方法一般是利用这个事件之后的概率，dp[j]表示从j状态转移到终止状态的期望

那么要求的其实就是从起始状态转移到终止状态的期望，即dp[s]

转移和求概率相反，是从“当前事件可能转移到的后续事件”推到当前事件，也就是倒推的一个过程

也就是dp[i] = sigma(p[j] * (dp[j] + cost))，j是i的后续状态，由于是求期望，所以要加上转移所消耗的步数cost

初始状态的话就是终止状态的dp[t] = 0

这样一般情况的概率&期望dp都可以解决了，最后会放一些这类可以直接解决的问题。

但是很多时候直接dp是解决不了问题的，为什么呢？因为状态之间的关系可能会比较复杂

比如求概率的时候，你在求某一个状态的时候可能无法确定所有需要知道的状态的概率。比如A可以转移到B，B也可以转移到A，那么求A的概率的时候需要知道B的概率，求B的概率的时候需要知道A的概率，这下就两个都求不出来了。

解决的方法就是先设这些概率为变量，然后建立相互之间关系的方程，最后消元得出答案

大概思路就是：每个转移方程可以构成一个包含方程，变量为各个状态的概率，这样就是n个未知数n个方程，理论上可以用高斯消元求出每个未知数

那么这种转移方程嵌套的问题也就能解决了

之前说过还有一种方法叫做人肉消元：这种方法是在高斯消元复杂度过高时候或者变量间关系不那么复杂时候运用的方法，比高斯消元能更快的得到答案。

具体是怎么搞呢？一般dp方程的形式是这样的：

dp[i] = sigma(A * dp[j]) + B * dp[i] + C

这里A，B，C都是常数，然后j是除了i之外的其他转移状态

我们发现，如果知道了所有的dp[j]，那么dp[i] = (sigma(A * dp[j]) + C) / (1-B) 就解出来了

那问题就转换成怎么在之前所说转移嵌套的情况下得到dp[j]们

我们把某一个dp状态（比如dp[0]）作为变量，使得每个dp状态的表示形式变成这样：

dp[i] = A*dp[0] + B * dp[i] + C

这样就可以从dp[0]的后续状态们开始推出所有的状态的表示，然后用已知数值的状态（求概率时是初始状态，求期望时是末尾状态）就可以轻松得到dp[0]，那么根据这个表示法所有状态都可以求出来了

剩下要做的就是怎么维护A,B,C这三个常数了。这个在基本转移方程的基础上只要用这种表达形式稍微写一下很容易就能得出写法，就是在思维的过程可能会相对较多一些。

由于这个方法实质是通过人脑构造的方法进行消元，因此我称之为”人肉消元“（笑）

这种方法非常优美相比高斯消元也很快，但从做题速度的角度来说，思维的速度和敲模板的速度还是见仁见智吧。

下面给一些这类的题目：

高斯消元：

HDU 4418 Time Travel

#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<string.h>
#include<queue>
#include<math.h>
using namespace std;

#define CLR(a,b) memset(a,b,sizeof(a))
#define eps 1e-9
const int MAXN=220;
double a[MAXN][MAXN],x[MAXN];//方程的左边的矩阵和等式右边的值，求解之后x存的就是结果
int equ,var;//方程数和未知数个数

int Gauss()
{
int i,j,k,col,max_r;
for(k=0,col=0;k<equ&&col<var;k++,col++)
{
max_r=k;
for(i=k+1;i<equ;i++)
if(fabs(a[i][col])>fabs(a[max_r][col]))
max_r=i;
if(fabs(a[max_r][col])<eps)return 0;
if(k!=max_r)
{
for(j=col;j<var;j++)
swap(a[k][j],a[max_r][j]);
swap(x[k],x[max_r]);
}
x[k]/=a[k][col];
for(j=col+1;j<var;j++)a[k][j]/=a[k][col];
a[k][col]=1;
for(i=0;i<equ;i++)
if(i!=k)
{
x[i]-=x[k]*a[i][k];
for(j=col+1;j<var;j++)a[i][j]-=a[k][j]*a[i][col];
a[i][col]=0;
}
}
return 1;
}
const int N = 111;
const int M = 211;
int nn,m,s,t,d,n;
double sum;
double p[M];
int has[M],vis[M],k;

int bfs(int u)
{
CLR(has, -1);
CLR(a, 0);
CLR(vis ,0);
int v,flg = 0;
queue<int> q;
q.push(u);
k = 0 ;
has[u] = k ++;
while(!q.empty()){
u = q.front();
q.pop();
if(vis[u]) continue;
vis[u] = 1;
if(u == t || u == n-t)
{
a[has[u]][has[u]] = 1;
x[has[u]] = 0;
flg = 1;
continue;
}
a[has[u]][has[u]] = 1;
x[has[u]] = sum;
for(int i = 1; i <= m ;i ++){
if(fabs(p[i]) < eps) continue;
v = (u + i) % n;
if(has[v] == -1) has[v] = k ++;
a[has[u]][has[v]] -= p[i];
q.push(v);
}
}
return flg;
}
int main()
{
int T;
scanf("%d",&T);
while(T--){
scanf("%d%d%d%d%d",&nn,&m,&t,&s,&d);
sum = 0;
for(int i = 1 ; i <= m ; i++){
scanf("%lf",&p[i]);
p[i] /= 100;
sum += p[i] * i;
}
if(s == t){
printf("0.00\n");
continue;
}
n = 2 * (nn-1);
if(d > 0)s = (n-s) % n;
if(!bfs(s)){
puts("Impossible !");
continue;
}
equ = var = k;
Gauss();
printf("%.2f\n",x[has[s]]);
}
return 0;
}

人肉消元：
HDU 4035 Maze

非常棒的一道题，也很经典，写完这题就可以体会人肉消元的美妙之处了

#include <stdio.h>
#include <cstdlib>
#include <string>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <stack>
#include <cmath>
using namespace std;

typedef long long LL;
typedef pair<int,int> PII;
#define CLR(x,y) memset(x,y,sizeof(x));
#define PB push_back
#define MP make_pair
#define INF 0x3f3f3f3f

const double EPS = 1e-9;
const int N = 10000+20;
int n;
vector<int> edge
;
int fa
;
bool leaf
;
double k
,e
;
double dp
,xi
,r
;

bool ok;
void dfs(int v,int w)
{
fa[v] = w;
leaf[w] = 0;
for(int i = 0; i < edge[v].size() ; i ++){
int u = edge[v][i];
if(u == w)continue;
dfs(u, v);
}
}
void work(int s)
{
if(leaf[s]){
dp[s] = (1-k[s]-e[s]);
xi[s] = k[s];
r[s] = (1-k[s]-e[s]);
return;
}
double sumdp = (1-k[s]-e[s]), sumxi = k[s];
r[s] = (1-k[s]-e[s]) / edge[s].size();
double sumr = 0;
for(int i = 0 ;i < edge[s].size() ; i ++){
int v = edge[s][i];
if(v == fa[s])continue;
work(v);
sumdp += (1-k[s]-e[s]) / edge[s].size() * dp[v];
sumxi += (1-k[s]-e[s]) / edge[s].size() * xi[v];
sumr += (1-k[s]-e[s]) / edge[s].size() * r[v];
}
if(1-sumr < EPS){
ok = 0;
return;
}
dp[s] = sumdp / (1-sumr);
xi[s] = sumxi / (1-sumr);
r[s] = r[s] / (1-sumr);
}

void solve()
{
CLR(dp, 0);
CLR(xi, 0);
CLR(r, 0);
ok = 1;
work(1);
if(!ok || dp[1] < EPS || 1-xi[1] < EPS){
printf("impossible\n");
}else{
dp[1] = dp[1] / (1-xi[1]);
printf("%.6f\n",dp[1]);
}
}
int main()
{
int T,cas = 0;;
scanf("%d",&T);
while(T--){
cas ++;
printf("Case %d: ",cas);
scanf("%d",&n);
for(int i = 1; i <= n ;i ++)edge[i].clear();
CLR(fa, -1);
for(int i = 1; i < n ;i ++){
int x,y;
scanf("%d%d",&x,&y);
edge[x].push_back(y);
edge[y].push_back(x);
}
for(int i = 1 ;i <= n ;i ++){
scanf("%lf%lf",&k[i],&e[i]);
k[i]/=100;
e[i]/=100;
}
for(int i = 1; i <= n ; i++){
leaf[i] = 1;
}
dfs(1,-1);
solve();
}
return 0;
}

HDU 4870 Rating

2014年多校第一场的题， 需要消元。当时用人肉消元写的，不过似乎其他人都是直接高斯消元的。。。

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<iostream>
#include<iomanip>
#include<algorithm>
#include<vector>
#include<string>
#include<list>
#include<stack>
#include<queue>
#include<set>
#include<map>
using namespace std;

#define CLR(x,y) memset(x,y,sizeof(x));
#define PB push_back
#define MP make_pair
typedef long long LL;
typedef pair<int,int> PII;
const int INF=0x3f3f3f3f;

const double EPS = 1e-8;
const int N = 1000+20;
double A

,B

,C

;
bool chk

;
double p;

void solve()
{
for(int i = 999; i >= 0 ; i --){
A[i][0] = 1;
B[i][0] = 0;
C[i][0] = 1/p;
for(int j = 1 ; j <= min(i,100); j ++){
A[i][j] = p;
B[i][j] = 1 - p;
C[i][j] = 1;
}
for(int j = 101; j <= i ; j ++){
int b1 = max(0,j-50);
int b2 = max(0,j-100);
double pp = (1-p) * A[i][b1] * A[i][b2];
A[i][j] = p / (1-pp);
B[i][j] = (1-p) * (A[i][b2]*B[i][b1] + B[i][b2]) / (1-pp);
C[i][j] = ((1-p) * (A[i][b2] *C[i][b1] + C[i][b2]) + 1 ) / (1-pp);
}
for(int j = i ; j >= 0 ; j --){
if(j + 50 > i){
C[i][j] += A[i][j] * C[j+50][i];
}else{
B[i][j] += B[i][j+50] * A[i][j];
C[i][j] += C[i][j+50] * A[i][j];
}
}
C[i][0] = C[i][0] / (1-B[i][0]);
for(int j = 1; j <= i ; j ++){
C[i][j] += B[i][j] * C[i][0];
}
}
printf("%.7f\n",C[0][0]);
}
int main()
{
while(~scanf("%lf",&p)){
CLR(chk, 0);
CLR(A,0);CLR(B,0);CLR(C,0);
solve();
}
return 0;
}

下面是一些简单一些的题：

POJ 2096 Collecting Bugs

直接逆推期望

#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <stack>
#include <cmath>
using namespace std;

typedef long long LL;
typedef pair<int,int> PII;
#define CLR(x,y) memset(x,y,sizeof(x));
#define PB push_back
#define MP make_pair
#define INF 0x3f3f3f3f

const int N = 1000+20;
double n,s;
double dp

;

int main()
{
while(~scanf("%lf%lf",&n,&s)){
dp[(int)n][(int)s] = 0;
for(int i = n ; i >= 0 ;i --){
for(int j = s; j >= 0; j --){
if(i == n && j == s)continue;
dp[i][j] = (dp[i][j+1] * i * (s-j) + dp[i+1][j+1] * (n-i) * (s-j) + dp[i+1][j] * (n-i) * j + n*s) / (n*s - i*j);
}
}
printf("%.4f\n",dp[0][0]);
}
return 0;
}

HDU 4405 Aeroplane chess

建图直接dp，也是非常直接

#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <stack>
#include <cmath>
using namespace std;

typedef long long LL;
typedef pair<int,int> PII;
#define CLR(x,y) memset(x,y,sizeof(x));
#define PB push_back
#define MP make_pair
#define INF 0x3f3f3f3f

const int N = 1000+20;
double n,s;
double dp

;

int main()
{
while(~scanf("%lf%lf",&n,&s)){
dp[(int)n][(int)s] = 0;
for(int i = n ; i >= 0 ;i --){
for(int j = s; j >= 0; j --){
if(i == n && j == s)continue;
dp[i][j] = (dp[i][j+1] * i * (s-j) + dp[i+1][j+1] * (n-i) * (s-j) + dp[i+1][j] * (n-i) * j + n*s) / (n*s - i*j);
}
}
printf("%.4f\n",dp[0][0]);
}
return 0;
}

HDU 4098 Activation

这题需要维护一下状态自己在转移方程里面的系数，也就是在消元那部分表示法里面的B这个数

#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <stack>
#include <cmath>
using namespace std;

typedef long long LL;
typedef pair<int,int> PII;
#define CLR(x,y) memset(x,y,sizeof(x));
#define PB push_back
#define MP make_pair
#define INF 0x3f3f3f3f

const double EPS = 1e-8;
const int N = 2000+20;
int n,m,k;
double p1,p2,p3,p4;
double dp

;
double xi
;
void solve()
{
CLR(dp,0);
double p21 = p2 / (1-p1),p31 = p3 / (1-p1), p41 = p4 / (1-p1);
dp[1][1] = p41 / (1- p21);
for(int i = 2 ; i <= n ; i++){
CLR(xi,0);
dp[i][1] = p41;
xi[1] = p21;
for(int j = 2 ; j <= i ; j ++){
dp[i][j] = p21 * dp[i][j-1] + p31 * dp[i-1][j-1] ;
xi[j] = p21 * xi[j-1];
if(j <= k)dp[i][j] += p41;
}
dp[i][i] = dp[i][i] / (1-xi[i]);
xi[i] = 0;
for(int j = 1; j <= i ; j++){
dp[i][j] = dp[i][j] + xi[j] * dp[i][i];
xi[j] = 0;
}
}
printf("%.5f\n",round(dp
[m]*100000)*1.0/100000);
}

int main()
{
while(~scanf("%d%d%d",&n,&m,&k)){
scanf("%lf%lf%lf%lf",&p1,&p2,&p3,&p4);
if(p4 < EPS)printf("%.5f\n",0.0);
else solve();
}
return 0;
}

POJ 3744 Scout YYF I

需要稍微推一下公式

#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <stack>
#include <cmath>
using namespace std;

typedef long long LL;
typedef pair<int,int> PII;
#define CLR(x,y) memset(x,y,sizeof(x));
#define PB push_back
#define MP make_pair
#define INF 0x3f3f3f3f

double fpow(double n,int k)
{
double ans = 1,base = n;
while(k){
if(k & 1)ans *= base;
base *= base;
k >>= 1;
}
return ans;
}
const int N = 300;
int m
;
double p;
int n;
double work(int l)
{
if(l == 0)return 0;
// cout << l << "!" << endl;
double ret = (fpow(p-1,l)-1) / (p - 2);
return ret * (1-p);
}
int main()
{
while(~scanf("%d%lf",&n,&p)){
double ans = 1.0;
int last = 1;
for(int i = 0; i < n ;i ++){
scanf("%d",&m[i]);
}
sort(m, m + n);
for(int i = 0 ;i < n ;i ++){
int x = m[i];
ans *= work(x-last);
last = x+1;
}
printf("%.7f\n",ans);
}
return 0;
}

ZOJ 3329 One Person Game

维护一下自身系数，也是很直接

#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <stack>
#include <cmath>
using namespace std;

typedef long long LL;
typedef pair<int,int> PII;
#define CLR(x,y) memset(x,y,sizeof(x));
#define PB push_back
#define MP make_pair
#define INF 0x3f3f3f3f

const int N = 1000+20;

double dp
;
int n,k1,k2,k3,a,b,c;
double p[19];
double xi
;
void solve()
{
double sum = k1*k2*k3;
for(int i = 1 ; i<= k1 ; i++){
for(int j = 1; j <= k2 ;j ++){
for(int k = 1 ;k <= k3; k ++){
if(!(i == a && j == b && k == c))
p[i+j+k]+=1;
}
}
}
for(int i = n ; i >= 0 ; i--){
dp[i] = 1;
xi[i] = 1.0/sum;
for(int j = i+1 ; j <= min(n,i+k1+k2+k3) ; j ++){
dp[i] += p[j-i]/sum * dp[j];
xi[i] += p[j-i]/sum * xi[j];
}
}
dp[0] = dp[0] / (1-xi[0]);
printf("%.13lf\n",dp[0]);
}
int main()
{

int T;
scanf("%d",&T);
while(T--){
CLR(p, 0);
CLR(dp, 0);
CLR(xi, 0);
scanf("%d%d%d%d%d%d%d",&n,&k1,&k2,&k3,&a,&b,&c);
solve();
}
return 0;
}

Codeforces 148D Bag of mice
转移略有点不好想，注意哪些状态是可以获胜的

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<string>
#include<cmath>
#include<vector>
#include<map>
#include<queue>
#include<stack>
using namespace std;

#define CLR(a,b) memset(a,b,sizeof(a))
const int N = 2000+20;

int w,b;
double dp

;
inline int lef(int i){
return w + b - i/2 - i;
}
inline int lefb(int i, int j){
return lef(i) - (w - j);
}
inline int lefw(int i, int j){
return w - j;
}
int main()
{
while(~scanf("%d%d",&w,&b)){
CLR(dp,0);
double ans = 0;
dp[0][0] = 1;
for(int i = 1 ; lef(i-1) > 0 ; i ++){
for(int j = 0 ;j <= (i-1)/2+1 && j <= w ;j ++ ){
if(i & 1){
if(lefb(i-1,j) >= 0){
dp[i][j] = dp[i-1][j] * lefb(i-1,j) / lef(i-1);
ans += dp[i-1][j] * lefw(i-1 ,j) / lef(i-1);
}
}else{
if(lef(i-1) == 1)continue;
if(lefb(i-1,j) > 0)
dp[i][j] += dp[i-1][j] * lefb(i-1, j) / lef(i-1) * (lefb(i-1,j)-1) / (lef(i-1)-1);
if(lefb(i-1,j-1) >= 0)
dp[i][j] += dp[i-1][j-1] * lefb(i-1, j-1) / lef(i-1) * (lefw(i-1,j-1)) / (lef(i-1)-1);
}
}
}
printf("%.9f\n",ans);
}
return 0;
}

ZOJ 3380 Patchouli's Spell Cards

这题解决的问题是：m个位置，每个位置放1-n中的一个数，问同一个数不少于L个的方案有多少

比较机智的递推，dp[i][j]表示用前i个数放j个位置且每个数不超过L的方案数

由于要表示成分数，需要java的BigInteger

import java.util.Scanner;
import java.math.BigInteger;

/**
* Created by NExPlain on 14-7-25.
*/

public class Main {
public static int N = 111;
public static BigInteger[] fac = new BigInteger
;
public static BigInteger[][] dp = new BigInteger

;
public static BigInteger C(int n,int m){
if(m < 0 || m > n)return BigInteger.ZERO;
return fac
.divide(fac[m]).divide(fac[n - m]);
}
public static BigInteger GCD(BigInteger a, BigInteger b){
if(b.equals(BigInteger.ZERO))return a;
return GCD(b, a.mod(b));
}
public static void main(String[] args){
fac[0] = BigInteger.ONE;
for(int i = 1 ; i < N ; i ++){
fac[i] = fac[i-1].multiply(BigInteger.valueOf(i));
}
int m,n,l;
Scanner cin = new Scanner(System.in);
while(cin.hasNext()){
m = cin.nextInt();
n = cin.nextInt();
l = cin.nextInt();
if(m < l){
System.out.println("mukyu~");
}else{
for(int i = 0 ;i <= n ; i ++){
for(int j = 0 ;j <= m ;j ++){
dp[i][j] = BigInteger.ZERO;
}
}
dp[0][0] = BigInteger.ONE;
for(int i = 1; i <= n ; i++){
for(int j = 1 ;j <= m ; j ++){
for(int k = 0 ; k < l && k <= j ; k ++){
dp[i][j] = dp[i][j].add(dp[i-1][j-k].multiply(C(m-(j-k),k)));
}
}
}
BigInteger up = BigInteger.ZERO;
for(int i = 1; i <= n ; i++)
up = up.add(dp[i][m]);
BigInteger down = BigInteger.valueOf(n).pow(m);
up = down.subtract(up);
// System.out.println(up + "/" + down);
BigInteger gg = GCD(up,down);
up = up.divide(gg);
down = down.divide(gg);
System.out.println(up + "/" + down);
}
}
}
}

HDU 4336 Card Collector

状压dp转移

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<string>
#include<algorithm>
#include<cmath>
using namespace std;

#define CLR(a,b) memset(a,b,sizeof(a))
const double EPS = 1e-9;

const int N = 21;
double p
;
int n;
double dp[1<<N];
void solve()
{
CLR(dp, 0);
dp[(1<<n)-1] = 0;
for(int i = (1<<n)-2 ; i >= 0 ; i --){
double sum = 0;
double nul = 0;
for(int j = 0 ;j < n ; j++){
if((i & (1 << j)) == 0){
dp[i] += (dp[i + (1 << j)] + 1) * p[j];
sum += p[j];
}else{
sum += p[j];
nul += p[j];
}
}
nul += 1 - sum;
dp[i] += nul;
dp[i] = dp[i] / (1 - nul);
}
printf("%.9f\n",dp[0]);
}
int main()
{
while(~scanf("%d",&n)){
for(int i = 0 ;i < n ; i++)scanf("%lf",&p[i]);
solve();
}
return 0;
}