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

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

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

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

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

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-12
using 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-12
using 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-12
using 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]]);

    }
}