进阶数学问题

这里称作进阶数学问题，其实着实是因为将线性方程组和期望的问题混合到一起总结，有些不太好起名字，所以采用了进阶数学问题。说是进阶，我觉得期望可能用的稍微多一点，还算是进阶，线性方程组直接用的应该挺少的，不过这里会给出一道利用线性方程组计算解的问题，个人感觉还是挺厉害的一道题。
一、线性方程组
关于线性方程组的解法，线性代数里面应该学了挺多了，有常见的高斯消元法，还有一些更快的算法，这里就不再多说了，如果有兴趣的话可以详细地去查一查。下面给出一般情况下都够用的高斯消元法：

const double EPS = 1E-8;
typedef vector<double> vec;
typedef vector<vec> mat;

// A是系数矩阵，b是结果向量
vec gauss_jordan(const mat& A,const vec& b){
    int n = A.size();
    mat B(n,vec(n+1));     //用矩阵B来存储系数矩阵和结果向量，方便之后统一消元
    for(int i = 0;i < n;i++){
        for(int j = 0;j < n;j++){
            B[i][j] = A[i][j];
        }
    }
    for(int i = 0;i < n;i++)B[i][n] = b[i];
    for(int i = 0;i < n;i++){
        /*这里pivot是想找到当前要消掉的这一列中，系数绝对值最大的那个，
        一方面保证当前的主元(就是用来消的系数)不为0，另一方面有人证明这
        样可以更稳定。*/
        int pivot = i;
        for(int j = i;j < n;j++){
            if(abs(B[j][i]) > abs(B[pivot][i]))pivot = j;
        }
        swap(B[i],B[pivot])  //如果不放心swap函数，可以自己写一个，不复杂
        
        //<EPS相当于系数是0了，经过上面的pivot，主元找最大的都是0，那么要不
        //就是无解，要么就是无穷个解了
        if(abs(B[i][i]) < EPS)return vec();
        for(int j = i+1;j <= n;j++)B[i][j] /= B[i][i];
        for(int j = 0;j < n;j++){
            if(i != j){
                for(int k = i+1;k <= n;k++)B[j][k] -= B[j][i] * B[i][k];
            }
        }
    }
    vec x(n);
    for(int i = 0;i < n;i++)x[i] = B[i][n];
    return x;
}

一道例题，算是线性方程组向期望过渡的题：
有一个N*M大小的格子，从(0,0)出发，每一步朝着上下左右4个各自中可以移动的格子等概率移动。另外有一些格子中有石头，因此无法移至这些格子。求第一次到达(N-1,M-1)格子的期望步数。题目假定至少存在一条从(0,0)出发到(N-1,M-1)的路径。

书中给出的期望公式： $E(x,y) = \tfrac{1}{4}E(x-1,y) + \tfrac{1}{4}E(x+1,y) + \tfrac{1}{4}E(x,y-1) + \tfrac{1}{4}E(x,y+1) + 1$
这里的公式其实来的很奇怪，这里的(x,y)与(x-1,y)等周围四个格子并没有明显的等式关系，甚至这里是两个变量的期望，为什么可以这样求？一种想法可能是，E(x,y)用p1(x1,y1)+p2(x2,y2)...这种，概率乘取值这种方法来求，但是这里显然不是，因为E(x-1,y)必然不是取值，它也是一个期望。上面也说了，它们之间也没有明显的加和关系，所以也不能直接拆开。综上，这个公式究竟是哪儿来的？
我翻了很多博客，没有人对这里有过讲解，仿佛二维图表中这样算期望已经被默认了。我后来也想了很多，只能有一个不敢保证对的方法，帮助自己和大家理解，说明这种方法确实是一个可以推广的方法。
对于这里，上下左右是等可能的，所以各是 $\tfrac{1}{4}$ ，就算到不了，它也占着 $\tfrac{1}{4}$ 的概率。接下来，设X₁₁，X₁₂...都是代表可以到达X₁的路径，X₁代表实际的(x-1,y)，然后 X₂等也有相应的含义。那么，E(X)= [p₁₁ * (X₁₁+1) + p₁₂ * (X₁₂+1)+...+p_1n * (X₁+1)+p₂₁ * (X₂₁+1)+...+p_2n * (X_2n+1)+...p₃₁ * (X₃₁+1)+...+p_4n * (X_4n+1)。

上面的公式确实挺乱。。所以我也不打算再敲公式了，不如这里口糊。。首先，这里肯定可以整体提个 $\tfrac{1}{4}$ 出来，于情，它们最后都要从各自的位置，转移到(x,y)，这一步转移的概率为 $\tfrac{1}{4}$ ；于理，任何数都可以提个 $\tfrac{1}{4}$ 出来啊。。因为提个 $\tfrac{1}{4}$ 就相当于乘了个4。接下来，提完以后的部分，先把+1都拿出来，可以凑出4(p₁₁+...+p_1n)这种式子，把4补进去，其实就是从起点出发，通过各条路径到达(x-1,y)的概率和，而概率和这个东西一定是等于1的，不管怎么计算，概率之和为1这一条铁定不会变，所以这一块凑起来是1；然后剩下的部分，其实很明显就是到达(x-1,y)的期望了，（概率*取值的形式）所以就变成了公式的样子。

上面的口糊也比较严重，不过我相信是可以理解的，并且最终是可以得到在格子上运动时，期望的公式的。所以这时我们把这个公式当作已知条件，接着往下探索。此时没有任何一个点的期望是已知的，所以直接求期望也不太现实，在这里，如果将每个E(x,y)都当成一个变量，也就是有N * M个变量的话，其实每一个变量都能列出一个方程。

对于是石头的格子，E(x,y) = 0。可以到达终点的格子，可以列出E(x,y) = E(x-1,y) + E(x+1,y) + E(x,y-1) + E(x,y+1) + 1的等式。如果想要优化的话，其实还可以提前把无法到达终点的格子找出来，然后它们的E(x,y)其实也是0，保险起见最好提前求一下。列出线性方程组之后，交给前面的gauss_jordan来处理就好了，返回值就是每个点的期望。完结撒花，接下来上代码：

char grid[max_n][max_m+1];    //+1是为了存字符串末尾的'\0'
int n,m;

bool can_goal[max_n][max_m];      //记录每个点是否可以到达终点
int dx[4] = {-1,1,0,0}, dy[4] = {0,0,-1,1};

// 搜索可以到达终点的点，记录
void dfs(int x,int y){
    can_goal[x][y] = true; //这里敢直接标成true,其实是因为搜索是从终点开始的逆向搜索
    for(int i = 0;i < 4;i++){
        int nx = x+dx[i],ny = y+dy[i];
        if(0 <= nx && nx < N && 0 <= ny && ny < M && !can_goal[nx][nyh] && grid[nx][ny] != '#'){
            dfs(nx,ny);
        }
    }
}

void solve(){
    mat A(N*M, vec(N*M, 0));     //mat是gauss_jordan定义的类型
    vec b(N*M, 0);
    for(int x = 0;x < N;x++){
        for(int y = 0;y < M;y++){
            can_goal[x][y] = false;      //初始化数组
        }
    }
    dfs(N-1,M-1);
    
    /* A是一个(N*M,N*M)的矩阵，也就是系数矩阵。第x*M+y行的方程代表的是(x,y)
    位置建立的方程。如果E(x,y) = 0，那么该行第x*M+y个变量t,t=0即可。如果E(x,y)
    不是0，就用期望的那个等式，将该行相邻位置的变量系数都设置一下，并且更新B向量.*/
    for(int x = 0;x < N;x++){
        for(int y = 0;y < M;y++){
            if((x == N-1 && y == M-1) || !can_goal[x][y]){
                A[x*M+y][x*M+y] = 0;
                continue;
            }
            int move_grid = 0;
            for(int k = 0;k < 4;k++){
                int nx = x+dx[k],ny = y+dy[k];
                if(0 <= nx && nx < N && 0 <= ny && ny < M && !can_goal[nx][nyh] && grid[nx][ny] == '.'){
                    A[x*M+y][nx*M+ny] = -1; //4E(x,y) - E(x-1,y) - E(x+1,y) - E(x,y-1) - E(x,y+1) = 4
                    move_grid++;
                }
            }
            b[x*M+y] = A[x*M+y][x*M+Y] = move_grid;
        }
    }
    vec res = gauss_jordan(A,b);
    printf("%.8f\n",res[0]);
}

这道题最难的，其实我个人觉得就是推导二维格子上的期望公式。如果有了这个，想到了列线性方程组求解，那么这个题就可以做了。

二、期望总结
在一些考数学的算法题目中，还比较喜欢考一些期望问题。这里我将期望问题，按照自己的想法分为三类。第一类是很偏数学的考点，就是考察对于期望的理解，以及期望公式的运用，比如E(x) = p₁ * x₁ + p₂ * x₂ + ... + p_n * x_n，这是期望最基本的定义；再比如说E(x+y) = E(x) + E(y)，两个随机变量独立时，E(xy) = E(x) * E(y)。第二类就是上面那道题目的类型，在二维格子中的期望公式，这个真的是需要有见过这类型的题目，知道这里的期望公式到底是怎么用的，才能做出来，一般情况下，都是列出线性方程组来求解的。第三类，是这里想要重点说一下的类型，是算法导论上的指示器随机变量，用它来求解很多期望问题，真的有种神兵利器的感觉。

首先先来一个问题，来感受一下指示器随机变量是什么：
抛n次硬币，求正面向上的次数？
ok这道题并不难，甚至可以说非常简单。求法也比较通俗易懂，用定义的话，正面向上i次的概率p_i = $\tfrac{1}{2^n}$ * Cⁱ_n，所以E(x) = $\tfrac{1}{2^n}$ * $sum_{i=1}^n$ Cⁱ_n * i，右边的累加，打开看一下的话，很容易发现其实就是 $sum_{i=1}^n$ Cⁱ_n-1，也就是2^n-1，与左边相乘，最后的结果是 $\tfrac{n}{2}$ 。

如果再换一种思路呢？这里其实很容易想到，每次抛掷硬币，其实是个独立重复的过程，相互之间并不影响，那么把每一次的期望加起来，不就是n次的期望了嘛。接下来，每一个抛掷就是一个两点分布，所以E(x) = p * 1 + (1-p) * 0 = p。所以最终E(X) = n*p = $\tfrac{n}{2}$ 。

其实分析到这里，指示器随机变量的两个核心问题已经都分析出来了，1. 将一个分n次进行的随机变量，分解为n个进行一次的随机变量。2. 每一个进行一次的随机变量服从两点分布，所以E(x) = p。合起来，E(X) = n*p。

到此，指示器随机变量的核心已经说完了。不过还不能完结撒花，我们需要了解为什么可以这样用，以及究竟怎么用，才能真正的在需要的时候用出来。如果只是用它来解决抛硬币的问题，这个方法有点过于鸡肋。。

关于指示器随机变量的第一个问题，为什么E(x) = n * E(x)。这里其实是需要看X这个随机变量，是否可以拆成n个子过程的。如果X本身可以拆成n个子过程，那么X = n*x。（这样表达不准确，其实应该是X = $sum_{i=1}^n$ x_i。然后代入期望的公式中，E(X) = E( $sum_{i=1}^n$ x_i) = $sum_{i=1}^n$ E(x_i)。其实我认为这才是指示器随机变量真正的核心，后面的两点分布，并不是每一次都能转化的出来的，但是一旦有了它，很多问题的难度就可以直接降一个level。

下一步的两点分布，通常定义为，发生的话，取值为1，未发生，取值为0，这里其实也没有太多要说的。好了，接下来看一些例题，这是比较有意思的地方。（例题参考了人家的博客，如有冒犯，请联系我，我会马上整改）

例一.
抛掷一枚骰子n次，求骰子和的期望。
因为这n次，其实可以拆成n个一次，与抛硬币一样，所以E(X) = $sum_{i=1}^n$ E(x_i)，每一次投掷骰子，很容易计算期望是3.5，所以结果是3.5n。这里就是一个，E(x_i)不是两点分布的例子。

例二.
一个随机的n排列，求所有排列中，满足i = a[i]的元素数目的期望。

这道题就没那么简单了，首先第一个问题，能不能把X拆分？这里的拆分方案，可能是Aⁿ_n种排列，计算每一种排列的期望。但是其实拿出每一次排列的话，这个时候就不是在算期望了，而是能确切地得出有几个i = a[i]了，而且这样计算也挺难的，我们研究指示器随机变量的做法。

i = a[i]，无非有n种可能，第一个位置上的值是1，第二个位置上的值是2...第n个位置上的值是n。这样可以顺利地将X拆成n个小的随机变量。根据公式，求出每个子随即变量的期望，就可以累加出X的期望。接下来的问题是，求每个位置上a[i] = i的期望。很明显这个问题要比前面的问题都难很多，这个问题的n个位置之间，是有相关性的，而不像前面的抛硬币，掷骰子这些，是重复独立实验。

计算各个位置a[i] = i的相关性，如果开始考虑位置之间的相关性，那这个题基本就做出来没啥希望了。。这里要用期望的眼光来分析这件事，我们这里算的，是当前位置E(a[i] = i)，并没有说求E(x|y)这种考虑条件的东西（事实上这种东西我也没见过），所以我们完全可以单独拿出i这个位置来分析，a[i] = i时，有n-1的全排列种情况，一共有n的全排列种情况，所以期望 = 1*(n-1)! / n!，也就是1/n。然后把n个位置累加起来，其实最后的结果就是1.（这里真的挺绕的）

例三.
有n种礼券，每次发一张，拿到每张的概率都是一样的，设X次可以领到n种不同礼券的话，求X的期望。

如果按照上面那道题的想法的话，有一个直接的猜测，计算每一张礼券的期望，然后再合起来，就是所有礼券的期望（还是相当于n张礼券是没有任何关联的）。但是这样行不通！这是这道题最大的特点，它事件的总数是不确定的。(这里拿一次礼券就算一次事件）第一次拿礼券，肯定是一张没拿到过的礼券；然后第二次拿，拿到之前拿到过的，这次相当于是没用的；理论上来说，是存在一直拿不完n张不同的礼券，所以要一直拿的情况的。所以这里说，事件的总数是不确定的。

上面那道题，可以不考虑任何顺序的问题，分析到i时，就只考虑i的情况，本质上是因为总共只有n!种情况，然后可以直接拿出来处理；但是这里是不可能直接拿出当前的情况与总的情况来处理的，所以一定不能单独分析第i张礼券拿的时候究竟是什么情况的，这里一定需要考虑之前的影响，来一步步地计算到n的。

这里说一下正确的思路，要拿到第一张不同的礼券，必然只需要1次；然后拿第二张，E(x₂) = 1 * $\tfrac{n-1}{n}$ + 2 * $\tfrac{1}{n}$ * $\tfrac{n-1}{n}$ + 3 * $(\tfrac{1}{n})^2$ * $\tfrac{n-1}{n}$ + m * $(\tfrac{1}{n})$ ^m-1 * $\tfrac{n-1}{n}$ 。要求解这个式子，其实只要知道一下几何分布的结论：1 * (1-p)⁰ * p + 2 * (1-p)¹ * p + ... + m * (1-p)^m-1 * p = $\tfrac{1}{p}$ 。上面那个式子其实就是一个几何分布，所以E(x₂) = $\tfrac{n}{n-1}$ ，...E(x_n) = $\tfrac{n}{1}$ ，综上E(X) = E(x₁) + E(x₂) + ... + E(x_n) = n( $\tfrac{1}{1} + \tfrac{1}{2} + ... + \tfrac{1}{n}$ )，这个式子，如果是 $10^8$ 以内的话，可以直接遍历一遍得到结果。如果n是10¹⁸这个数量级的话，那肯定不能遍历了。这个问题一般是用 $\tfrac{1}{1} + \tfrac{1}{2} + ... + \tfrac{1}{n}$ = ln(n+1) + r来实现。r是欧拉常数，一般取0.5772156649够了。

其他还有指示器随机变量的例题，比如说雇佣问题，上楼梯问题，坐公交车问题等，花样很多，这里就不一一说明了。个人认为，如果考期望问题，逃不过上面这三类；如果考指示器随机变量，大体就上面两类，即总方案数确定和总方案数不确定两种。