关于高斯消元求期望的个人理解

看了一天的高斯消元求期望。

但是貌似因为自己的概率论学的很差，所以总结了一套自己独特的推断。

没有证明，纯粹是个人想法。

个人感觉，做高斯期望求概率的题目分为以下几步：

1.  遍历所有出现的状态，找到所有能达到的状态并重新将其进行编号，可以存在一个num[]。

2. 判断出口（终点）是否被编号，意思就是判断能不能出去（到达），不能到达直接输出不能到达。

3.  遍历所有出现的状态，找到所有标过号的状态 E[num[i]]，那么每个E[num[i]]都满足：

     E[num[i]]=(E[num[0]]+step[0])*p[0]+…(E[num[j]]+step[j])*p[j]    其中E[num[0~j]] 为所有E[num[i]] 能达到的状态，step[0~j]为所需要的步数(比如在迷宫图中，step[0~j]=1代表走一步能到)，p[0~j] 代表到达这种状态的概率（这个概率不能为零，要在之前给排除，不然会让方程无解！！）

    假设可达状态有cont个，则所有的状态就构成了cont*(cont+1)的增广矩阵。

    这时候再令所有的出口（终点）的期望为0（因为如果当前就在终点就不用走路了，所以期望是0）

4. 运用高斯消元法求解实数增广矩阵（注意精度每次用绝对值最大的行去消元），求解出起点编号所对应的解 x[num[start]] 就是我们要的答案了。

注意点：

1.不能达到的状态一定要舍去，不然会让方程无解。

2.数据记得初始化！

3.为了减少精度丢少，每次消元的时候用绝对值最大的进行消元。

4.注意控制状态数，高斯消元是O(n^3)的算法，若状态过多，则尝试转换需找状态的思路。（例如：要在E(i,j)中i==j的状态，其实可以当做是找状态E(i-j=0)=E(0)的状态）

几道题目：

hdu 2262 Where is the canteen 

题意：n*m的迷宫图，”.“是空地、”@“是起点（仅有一个）、”X“是墙、”$“是终点（不止一个），问从起点到终点的步数期望。

其实很简单就是bfs跑一遍图，然后把点重新编号，不能到的点去掉。

然后在外面遍历一遍n*m建立增广矩阵（其实在bfs里面直接建也可以，但是个人不喜欢这样），然后判断能否走到终点，能就进行高斯消元。

然后求出期望。这题同前面一个题解的公式就不列了。

#include"cstdlib"#include"cstdio"#include"cstring"#include"cmath"#include"queue"#include"algorithm"#include"iostream"#define eps 1e-12using namespace std;struct node{    int x,y;};int f;int n,m;int equ,var;char map[20][20];double a[300][300];double x[300];int num[20][20];int move[4][2]= {{1,0},{0,1},{-1,0},{0,-1}};void debug(){    for(int i=0;i<equ;i++)    {        for(int j=0;j<=var;j++) printf("%3.2f ",a[i][j]);        puts("");    }    puts("");}int ok(int x,int y){    if(x<0||y<0||x>=n||y>=m) return 0;    if(map[x][y]=='#') return 0;    return 1;}int bfs(int xx,int yy){    int i,id=0;    node cur,next;    queue<node>q;    cur.x=xx;    cur.y=yy;    q.push(cur);    while(!q.empty())    {        cur=q.front();        q.pop();        if(num[cur.x][cur.y]!=-1) continue;        if(map[cur.x][cur.y]=='$') f=1;        num[cur.x][cur.y]=id++;        for(i=0; i<4; i++)        {            next.x=cur.x+move[i][0];            next.y=cur.y+move[i][1];            if(ok(next.x,next.y)) q.push(next);        }    }    return id;}void gauss(){    int i,j;    int row,col;    for(row=0,col=0;row<equ&&col<var;row++,col++)    {        int maxr=row;        for(i=row+1;i<equ;i++) if(fabs(a[i][col])>fabs(a[maxr][col])) maxr=i; //为找绝对值最大的行        for(i=0;i<=var;i++) swap(a[row][i],a[maxr][i]);        for(i=row+1;i<equ;i++)        {            if(fabs(a[i][col])>eps)            {                double ta=a[i][col]/a[row][col];                for(j=col;j<=var;j++) a[i][j]-=ta*a[row][j];             }        }    }    //debug();    for(i=row-1;i>=0;i--)    {        double tep=a[i][var];        for(j=i+1;j<var;j++)        {            if(fabs(a[i][j])<=eps) continue;            tep-=a[i][j]*x[j];        }        x[i]=tep/a[i][i];    }}int main(){    while(cin>>n>>m)    {        int i,j,k;        for(i=0; i<n; i++) cin>>map[i];        for(i=0; i<n; i++)        {            for(j=0; j<m; j++)                if(map[i][j]=='@') break;            if(j!=m) break;        }        memset(num,-1,sizeof(num));        memset(a,0,sizeof(a));        f=0; //用f判断能不能到达终点        var=equ=bfs(i,j);        if(!f)        {            puts("-1");            continue;        }        for(i=0;i<n;i++)        {            for(j=0;j<m;j++)            {                if(num[i][j]==-1) continue;                int tep=num[i][j];                int cnt=0;                for(k=0;k<4;k++)                {                    int xx=i+move[k][0];                    int yy=j+move[k][1];                    if(xx<0||yy<0||xx>=n||yy>=m) continue;                    if(map[xx][yy]=='#'||num[xx][yy]==-1) continue;                    cnt++;                    a[tep][num[xx][yy]]-=1;                }                a[tep][tep]=cnt;                a[tep][var]=cnt;                if(map[i][j]=='$')                {                    memset(a[num[i][j]],0,sizeof(a[num[i][j]]));                    a[tep][tep]=1;                }            }        }        //debug();        gauss();        printf("%.6f\n",x[0]);    }}

zjut  1317 掷飞盘

http://cpp.zjut.edu.cn/ShowProblem.aspx?ShowID=1317

这题我们以两个飞盘的距离为状态进行转移。
那么E[n]=E[n+2]/4+E[n-2]/4+E[n]/2+1，化简成:2E[n]-E[n+2]-E[n-2]=4。
首先对于两个飞盘给定的起始距离n，我们可以先搜索一下可否到达状态0，如果不行，则直接输出INF。在搜索的过程中，顺便把重新标号也进行了。为什么这题也要重新标号呢？
因为该题中，假设m是偶数，那么对于任意的n，n+1和n-1都是不可达的状态。请一定记得，如果方程组中有不可达的点的话，就会导致无解的情况！

#include"cstdlib"#include"cstdio"#include"cstring"#include"cmath"#include"queue"#include"algorithm"#include"iostream"#define exp 1e-12using namespace std;int equ,var;int num[505];double x[500];double a[500][500];int move[4][2]={{1,1},{1,-1},{-1,1},{-1,-1}};int m,n;void debug(){    for(int i=0; i<equ; i++)    {        for(int j=0; j<=var; j++) printf("%3.2f ",a[i][j]);        puts("");    }    puts("");}struct node{    int x,y;};int sign(int x){    return x<0?-x:x;}void gauss(){    int i,j;    int row,col;    for(row=0,col=0; row<equ&&col<var; row++,col++)    {        int maxr=row;        for(i=row+1; i<equ; i++) if(fabs(a[i][col])>fabs(a[maxr][col])) maxr=i;        for(i=0; i<=var; i++) swap(a[row][i],a[maxr][i]);        for(i=row+1; i<equ; i++)        {            if(fabs(a[i][col])>exp)            {                double ta=a[i][col]/a[row][col];                for(j=col; j<=var; j++) a[i][j]-=a[row][j]*ta;            }        }    }    for(i=row-1; i>=0; i--)    {        double tep=a[i][var];        for(j=i+1; j<var; j++)        {            if(fabs(a[i][j])<=exp) continue;            tep-=a[i][j]*x[j];        }        x[i]=tep/a[i][i];    }}int bfs(int x,int y)   //重标号  标记所有可到达的点{    int id=0;    node cur,next;    queue<node>q;    cur.x=x;    cur.y=y;    q.push(cur);    while(!q.empty())    {        cur=q.front();        q.pop();        if(num[min(sign(cur.x-cur.y),sign(min(cur.x,cur.y)+m-max(cur.x,cur.y)))]!=-1) continue;        num[min(sign(cur.x-cur.y),sign(min(cur.x,cur.y)+m-max(cur.x,cur.y)))]=id++;        for(int i=0;i<4;i++)        {            next.x=cur.x+move[i][0];            next.y=cur.y+move[i][1];            if(next.x==0) next.x=m;            if(next.y==0) next.y=m;            if(next.x==m+1) next.x=1;            if(next.y==m+1) next.y=1;            q.push(next);        }    }    return id;}int main(){    while(cin>>m>>n)    {        if(n>m/2)             // 注意一下 很关键！            n=m-n;        memset(num,-1,sizeof(num));        memset(a,0,sizeof(a));        equ=var=bfs(1,1+n);        if(num[0]==-1)       //到达不了直接输出INF        {            puts("INF");            continue;        }        int i;        for(i=1;i<=m;i++)    //注意一下 i-2和i+2之后的距离变化就好了        {            if(num[i]==-1) continue;            int tep=num[i];            a[tep][num[sign(i-2)]]-=0.25;            a[tep][num[min(i+2,m-i-2)]]-=0.25;            a[tep][tep]-=0.5;            a[tep][tep]+=1;            a[tep][var]=1;        }        a[num[0]][num[0]]=1;        gauss();        printf("%.2f\n",x[num[n]]);    }}

hdu 4870 Rating

题意：一个人注册两个账号，初始rating都是0，他每次拿低分的那个号去打比赛，赢了加50分，输了扣100分，胜率为p，他会打到直到一个号有1000分为止，问比赛场次的期望

思路：f(i, j)表示i >= j，第一个号i分，第二个号j分时候，达到目标的期望，那么可以列出转移为f(i, j) = p f(i', j') + (1 - p)f(i'' + j'') + 1
f(i', j')对应的是赢了加分的状态，f(i'', j'')对应输的扣分的状态，跑一遍找出可达状态，利用高斯消元去求解方程组，解出f(0, 0)就是答案。

很明显此题的终点只有一个 f(950,1000)

#include"cstdlib"#include"cstdio"#include"cstring"#include"cmath"#include"queue"#include"algorithm"#include"iostream"#define exp 1e-12using namespace std;int equ,var;double x[320];double a[300][300];int num[1200][1200];struct node{    int x,y;};void debug(){    for(int i=0; i<equ; i++)    {        for(int j=0; j<=var; j++) printf("%3.2f ",a[i][j]);        puts("");    }    puts("");}void gauss(){    int i,j;    int row,col;    for(row=0,col=0; row<equ&&col<var; row++,col++)    {        int maxr=row;        for(i=row+1; i<equ; i++) if( fabs(a[i][col])>fabs(a[maxr][col])) maxr=i;        for(i=0; i<=var; i++) swap(a[row][i],a[maxr][i]);        for(i=row+1; i<equ; i++)        {            if(fabs(a[i][col])>exp)            {                double ta=a[i][col]/a[row][col];                for(j=col; j<=var; j++) a[i][j]-=a[row][j]*ta;            }        }    }    for(i=row-1; i>=0; i--)    {        double tep=a[i][var];        for(j=i+1; j<var; j++)        {            if(fabs(a[i][j])<=exp) continue;            tep-=a[i][j]*x[j];        }        x[i]=tep/a[i][i];    }}int bfs(int x,int y){    queue<node>q;    node cur,next;    cur.x=x;    cur.y=y;    q.push(cur);    int id=0;    while(!q.empty())    {        cur=q.front();        q.pop();        if(num[cur.x][cur.y]!=-1) continue;        num[cur.x][cur.y]=id++;        next.x=cur.x+50;        next.y=cur.y;        if(next.x<=950)        {            if(next.x>next.y) swap(next.x,next.y);            q.push(next);        }        next.x=max(0,cur.x-100);        next.y=cur.y;        if(next.x<=950)        {            if(next.x>next.y) swap(next.x,next.y);            q.push(next);        }    }    return id;}int main(){    double p;    memset(num,-1,sizeof(num));    int cont=bfs(0,0);    num[950][1000]=cont++;    equ=var=cont;    while(cin>>p)    {        memset(a,0,sizeof(a));        double q=1-p;        int i,j,k;        for(i=0; i<=950; i+=50)        {            for(j=i; j<=950; j+=50)            {                if(num[i][j]==-1) continue;                int tep=num[i][j];                a[tep][tep]+=1;                a[tep][var]=1;                int xx,yy;                xx=max(0,i-100);                yy=j;                if(num[xx][yy]!=-1) a[tep][num[xx][yy]]-=q;                xx=i+50;                yy=j;                if(xx>yy) swap(xx,yy);                if(num[xx][yy]!=-1) a[tep][num[xx][yy]]-=p;            }        }        memset(a[num[950][1000]],0,sizeof(a[num[950][1000]]));        a[num[950][1000]][num[950][1000]]=1;        gauss();        printf("%.6f\n",x[num[0][0]]);    }}