在前面的一系列文章中，我们用数学方法解决了这样一个问题：在 维空间中的 维超球面上随机、均匀地取 个点，它们位于同一个半超球面的概率是多少？答案是 。

　　不知道大家注意到没有，在知乎上比较复杂的概率问题下面，除了用数学方法解答的答案以外，总有些答案会借助编程模拟来提出猜想或者验证结果。现在咱们也来试着编个程序，验证一下上面的结果吧！

　　验证结果的程序，需要做两件事情：

　　1. 在高维超球面上随机、均匀地取点；
　　2. 对于给定的一组点，判断它们是否位于同一个超球面上。

这两件事情似乎都不简单。

　　件事其实是有「标准答案」的：从高维标准正态分布中随机采样，然后把取到的点归一化到单位超球面上去。从高维标准正态分布中随机采样，其底层使用的是 Box-Muller 变换，这也是一个非常巧妙的算法。不过我们现在倒不需要深入研究其原理，因为许多编程语言已经提供了「从高维标准正态分布中随机采样」的功能，比如 Matlab 的 randn 函数。利用这个函数，我们可以写出如下的程序框架，来验证我们算出的概率：

function prob = simulate(n, d)  % 估算d维空间中d-1维超球面上n个点位于同一个半超球面的概率
    total = 10000;              % 总共模拟10000次
    success = 0;                % n个点位于同一个半超球面的次数

    for i = 1:total
        % 从d维标准正态分布中随机抽取n个点（每行为一个点的坐标）
        points = randn(n, d);

        % 把所有点都归一化到单位超球面上
        % 在较老版本的Matlab上，要写成：
        % points = bsxfun(@rdivide, points, sqrt(sum(points .^ 2, 2)));
        points = points ./ sqrt(sum(points .^ 2, 2));

        % same_hemisphere判断多个点是否位于同一超球面上
        success = success + same_hemisphere(points);
    end

    prob = success / total;     % 返回估算出的概率
end

　　框架中的 same_hemisphere 函数，就负责来做第二件事 —— 判断多个点是否位于同一超球面上。这件事又该怎么做呢？

　　观察上图。以 A、B、C、D 为极点的四个半球面有交集（黑色区域，称为「可行域」），所以 A、B、C、D 点可以被同一个半球面所包含，这个半球面的极点可以在黑色区域内任取。要判断一组点是否位于同一个半（超）球面上，其实就是判断以它们为极点的半（超）球面是否有交集。在二维中，这样的交集是一段圆弧；在三维中，这样的交集是上图中那样的球面多边形。

　　如果你之前学过计算几何，也许现在就已经跃跃欲试了。计算几何里有一个经典问题是「半平面交」，我们在这里要做的是「半球面交」。「半平面交」交出来的结果是平面凸多边形，「半球面交」交出来的，就是球面凸多边形了。这样的图形可以用顶点的序列来表示……先不要着急，休息一会儿。「半球面交」的程序固然能写，但写起来会很麻烦；并且，「半球面交」只能解决三维情况，更高维的情况还解决不了。我们还需要更巧妙的办法！

　　再观察图 5.1。如果黑色区域存在，那么它一定会有顶点（忽略两个大圆重合等概率为零的特殊情况），这些顶点都是两个大圆的交点。我们可以检验大圆两两之间的交点，只要发现有一个交点使得以它为极点的半圆能够把输入的所有点都包含在内，那么就可以说，输入的所有点都位于同一个半球面了。

　　首先我们来枚举大圆的交点。以红圈和绿圈为例，设它们的两个交点为 X、X'。容易发现，向量 、 与向量 、 都垂直。在三维空间中，可以计算向量 、 的叉积并归一化，归一化后的坐标及其关于球心的对称点，就是 X 和 X'。下面就要检验 A、B、C、D 四点是否都在以 X（或 X'）为极点的半球面内了，这其实就是要看 、、、 这几个点积是否都同号（包括等于 0）。由于向量 本身就是由 、 叉乘而得的，所以点积 、 肯定等于 0，只需看 与 、 的点积是否同号。在图 5.1 中，这两个点积均为正，所以 A、B、C、D 位于同一个半球面上。

　　把上面的方法推广到高维。在 维空间中，要确定可行域的一个顶点 X，需要指定 个输入点，超球心 O 到这些输入点的向量都与 垂直。Matlab 提供了一个 null 函数，输入一个矩阵，可以求出它的零空间的一组基。把 个输入点排成一个 的矩阵，在非特殊情况下，这个矩阵的秩就应该是 ，其零空间只有 1 维，基向量就是 。枚举从 个输入点中选取 个的所有可能。只要在其中一种情况下，发现 与其余输入点的点积都同号（包括等于 0），就可以断定所有的输入点位于同一个半超球面上了。

　　在此基础上加上一些特殊情况的判断，就可以写出 same_hemisphere 函数的代码了：

function result = same_hemisphere(points)   % 若所有点位于同一个半超球面上，返回true
    [n d] = size(points);                   % n：点数，d：空间维数
    
    r = rank(points);                       % 所有输入点的秩
    if r < d                                % 若秩小于d，则所有输入点位于一个d-1维子空间里
        result = true;                      % 过超球心作一个d-1维超平面与这个子空间平行，
                                            % 则输入点均在其同侧，故在同一个半超球面上
        if r < n, disp('你中奖了！'); end   % 通常情况下秩应该等于点数，小于点数的概率为0
        return;                             
    end
    
    permutations = nchoosek(1:n, d-1);      % 枚举从n个点中选d-1个点的所有可能
    for i = 1:length(permutations)
        mask = false(n, 1);
        mask(permutations(i, :)) = true;
        X = null(points(mask, :));          % 求选中的d-1个点排成的矩阵的零空间的基
        if size(X, 2) > 1                   % 通常情况下零空间应该只有一维，高于一维的概率为0
            disp('你中奖了！');
            continue;                       % 如果真的高于一维了，则确定不了顶点X，跳过
        end
        prod = points(~mask, :) * X;        % 求未选中的点与X的点积
        if all(prod >= 0) || all(prod <= 0) % 若所有点积都同号（包括等于0）
            result = true;                  % 则可断定所有输入点位于同一个半超球面上
            return;
        end
    end
    
    result = false;
end

　　运行 simulate(6, 4)，就可以估算 4 维空间中的 6 个点位于同一个半超球面的概率。程序需要枚举在 6 个点中选 3 个点的所有可能，复杂度较高，检验 10000 组点需要花费 3 秒左右。运行 10 次得到的结果如下：

0.8119  0.8121  0.8114  0.8098  0.8175  0.8177  0.8065  0.8088  0.8158  0.8132

其平均值为 0.81247，十分接近我们算出的结果 。

　　代码中有两处「你中奖了」，是用来检测一些概率为 0 的特殊情况的。我参加信息学竞赛时省队的一个队友，在后来上「数值分析」课上机实验时，就遇到过概率为 0 的特殊情况。然而我把 simulate(6, 4) 运行了 1000 次，也就是检验了 1 千万组点，都没有中过奖。

参考文献

• 第 2 篇中 3blue1brown 的视频链接：[YouTube] [Bilibili]
• 第 3、4 篇中的方法主要来自这个网页，其中求解递推式使用的是多项式拟合法。
• 第 5 篇中判断多个点是否位于同一个半超球面上的算法来自这个网页，它只考虑了三维情况。

　　本系列共有 5 篇文章，以下是传送门：(1) (2) (3) (4) (5)