圆内均匀随机地取三个点，三点确定的新圆在旧圆内的概率是多少？ 第1页

先给结论：

在n维的单位球内独立、（按体积）均匀地取n+1个点，它们确定的新球被原单位球所包含的概率为

前几项为

数值验证了一下，应该是正确的

至于过程嘛，你看这答案就知道过程不会太简单了，等我有空再来写...

1.总体思路

在 维空间中取 个点，总的样本空间是 维空间中的某个区域 ，其上的概率密度为常数 ，其中 是 维单位球的体积。满足题中条件的样本构成样本空间的一个子集 ，我们只要求出其 维体积即可得到对应的概率。现在的问题是要直接在原来的坐标系下描述这个集合比较困难，想求体积就必须换一个方式来刻画（换元），从而得到较为规则的积分区域和较为容易计算的积分。我的方法与 @Ao Sun 的方法类似，但又稍有不同：

注意到任意两个 维球之间都构成位似关系，且一般有两个位似中心，对应的位似比一正一负。考虑正的那个位似中心，我们发现两球为包含关系当且仅当该位似中心在大球内部。因此我们可以用球内的位似中心 ，位似比 以及球面上的 个点 来表示满足条件的点：

其中 . 除去测度为零的退化情况外，该表示方法恰好一一对应了所有满足条件的情形，因此我们可以将 的体积写成

2.初遇难关

积分区域现在已经非常规则，但问题转移到了被积函数——体积元之比上。二维的情形较为简单，可以以弧长/极角作为坐标，有 ，直接计算雅可比行列式得

而三维以上的表达式就颇为恐怖了，因此如何计算这个体积元之比成为了我们遇到的第一个难题。由于时间有限，就先更新到这里，有兴趣的读者可以自行尝试去找到一般的规律并得到其表达式。实际上，我一开始也是猜出来的。

3.大胆猜想

为了得到一般的规律，我们首先硬着头皮算一下三维情形的结果。用球坐标 表示球面上的点，结果为：

这当然不是手算出来的，否则我已经吐血翻白眼不省人事了。观察二维和三维的两个式子，很容易发现带 的因子有很简单的规律：前 行都带 ，后面有 行带 ，于是贡献的因子为 。然而剩下的部分才是重头戏所在，我们看这个因子

这个式子看起来如此规律，其是否有什么简洁明了的意义呢？答案是有！其恰好就是圆上三个点构成三角形的面积的2倍！于是我们不禁思考：其它维度的表达式是否也有类似的意义？代入几个特殊值验证一下，我们很容易发现，三维情形的表达式恰好就是球面四个点构成四面体的体积的6倍！

一个大胆的猜想浮上心头： 维情形下的这个因子，是否恰好是球面上 个点构成的 维单形（三角形、四面体等在 维空间中的对应物）的 维体积的 倍？写成数学语言就是

用数学软件验证了几个 ，发现都是正确的！我们几乎可以确定，这就是一般情况下的表达式！接下来的任务就是给出一个证明。由于我最近比较忙，并且后面还有很多步骤，所以这篇回答可能会更很久很久。有兴趣的读者可以自行尝试往下推导，说不定你就能得到一个比我更简洁的解法呢？

4.小心求证

要证明上面的猜想，说简单不简单，说难其实也不难。由于我们要求的是两个体积元之比，它们各自都有良好的定义，因此这个比值与我们采用的坐标系统无关。只要我们使用 个参数来描述，并且在某点附近不退化，那么尽管对应的雅可比行列式会随着坐标系选取而改变，该点处两个体积元的比值却是一个不变量。因此，我们可以选择一个用起来最舒服的坐标系来计算。

如上图所示，在每个 处，我们选取包含 本身的一组标准正交基。其中 是径向单位矢量，其余的 则是切向的单位矢量。对每个 ，我们采取 处的基底。如此我们对球面上的几个点采用了不同的坐标系，但上面已经论证过这不会影响最终的结果。

现在事情变得简单了不少：显然有

由于我们的基底本身就是按照球面面元的方向选取的，才能把面元写成如此简单的形式。现在

这是一个雅可比行列式，里面的每个元素都可以简单计算得到：

这玩意看起来令人晕乎乎的，实际上只是不好排版而已，算起来十分简单——下面那些行每行都仅有一个非零元素1，因此可以直接划掉它们所在行和列，最后只剩下这么个东西：

再注意到 ，代入即得：

又由于 ，最后一行里的点积部分可以写成前面几行的线性组合，因此最后有

我们终于证明了体积元之比的表达式：

将其代入最初的积分式，把所有变量积掉，我们就得到了一个看起来比较简单的式子：

其中 为 维单位球的体积和表面积， 为单形 的平均体积。看起来是不是有点结果那味了？有一部分因子已经出来了。

于是乎，我们将原问题转化为了一个全新的问题：在 维单位球表面独立、（按面积）均匀地取 个随机点，它们构成的 维单形的体积期望是多少？这个问题看起来简化了一些，毕竟"体积"比之"外接球球心"来说要好处理得多，但依然不容易，解决之路依然漫长。

5.准备工作

在介绍新问题的解决方法之前，我们需要先做一些准备工作，以便接下来的过程能够被顺畅地阅读。

首先定义Beta函数，其中 为正实数：

再定义Beta分布的概率密度函数，这是一个 上的分布：

之所以要特别提出Beta分布，是因为其与球内的分布有着十分密切的联系，以至于接下来我们会不止一次遇到这个分布。

一个Beta分布密度函数在另一个Beta分布下的期望值为：

在 维单位球内随机取一个点，其半径 是取值在 上的随机变量，累积分布函数为 ，概率密度函数为

6.叠 叠 乐

现在我们已经做好了准备工作，可以开始向着答案进军了。作为第一步，考虑一个更简化的问题：将其中一个点固定到原点处，则形成的单形的体积期望是多少？这个问题比原问题的简单之处在于，有 条棱是经过球心的，因此它们张成的平面等也经过球心，使得计算有章可循。

考虑依次在球面上选取 个随机点，现在已经选取了 个的情况。记这 个点为 ，则 张成一个 维的子空间 。现在考虑第 个点，我们将其位矢分解为垂直和平行分量 ，其中 。

那么 维单形 的 维体积为

其中 为 维单形 的 维体积。这其实就是底乘高的高维推广。注意到分子中的两项是独立的：无论前 个点如何分布，由于球对称性，都不会影响 的分布。对于独立的随机变量，积的期望等于期望的积，因此我们有

下面我们来求 的值。为此我们首先考虑 在球内而不是球面上均匀分布的情况。

对于任意 的情形，其 满足 ，且体积元 ，因此若取总体为 上的均匀分布，则

从而有 的概率密度

于是 的概率密度为

至于常数系数我们不关心，因为从上面的式子我们已经足以看出 服从Beta分布，即

于是有

接下来我们想办法将其转化为 在球面上均匀分布的情况。对于函数 ，如果

即 的径向和角向是完全独立的，则有

对于 来说， ，

现在，我们可以用递推的方法得到所求的体积期望

是不是又看到了一些熟悉的因子？我们已经解决了简化版的问题，现在只要将简化版问题的答案和原问题的答案联系起来，就能得到最终的结果了！

7.建立联系

那么我们要如何在两个问题之间建立联系呢？不妨观察一下，原图形的体积和简化图形的体积之间有什么关系？我们很容易发现，单形 的体积可以写成 个（简化的）单形体积的代数和：

其中每一项的正负号是我们主要关心的点。在什么情况下这个符号会变成负的？稍微想象一下就能得到答案： 的符号是负的，当且仅当 在 构成的超平面异侧。我们可以以此作为出发点。

将 构成的平面作为 的平面，则我们关心的是 在垂直方向上的坐标，即 。若 ，则 贡献为负号，反之为正号。现在我们考虑 的对踵点，也就是以球心为中心的对称点。显然 和 处的概率密度是相同的，而对于两种情况的平均体积来说，若 ，则 的贡献正负抵消，为零，否则为一。若能得到 和 的分布，以及"底面积" 的概率分布，就能计算出 对体积期望的总贡献。而 个单形对总体积的贡献必然是相等的，因此我们可以写出体积期望的表达式：

同样的，我们可以由此写出 的体积期望：

两个式子看起来非常相似，不过我们仍需进行进一步的计算才能发现两者之间的联系。

球面上 的面元为

从这个形式中我们能发现

于是

这玩意看起来有点复杂，而剩下 和 也不像是什么简单东西，而我们还需要对这些复杂的东西做积分。这看上去非常得难，不过做到最后我们会发现，其实我们是在自己吓自己。

8.故技重施

那么我们要怎么去求剩下的那些分布函数呢？回顾一下本回答的开头，我们采用了一套新的坐标体系来描述这些点，并且通过体积元之比将两者联系起来。下面我们用同样的方法，选取一套包含了 的坐标系来描述球面上的 个点：

首先，选定随机变量 ，其分布函数未知；

接下来选定球面上的随机单位向量 ，其分布为球面上均匀分布；

在赤道面 上选定 个随机点 ，其分布未知，是否独立未知；

通过某个正交旋转，将赤道面变换为与 垂直的大圆，得到上面 个点的像

将大圆沿 的方向平行移动距离 ，且同时随着球面而缩小，得到

如此，我们就得到了球面上的 个随机点，且它们构成的超平面与原点的距离为 。

而我们的任务就是找到这样一个分布，使得上面的采样方法与直接在球面上独立地取 个随机点等价。容易证明的是，除去退化情形外，上述表示方法恰好构成一一对应，因此我们所要求的概率密度实际上就是两个坐标空间中的体积元之比。

对于正交旋转的选择，一个很自然的想法是只进行 面内的旋转，而保持所有其它分量不变，其中 为垂直于赤道面的单位矢量。这个旋转的表达式不难求得，为

其中 分别为 的位矢。最终得到的 点位矢为

这个变换的性质非常不错，可以帮我们大大减少计算难度。首先，其是一个仿射变换，因此位矢的变化保持线性关系；其次，其在超平面内是“正交”的：

因此我们可以直接由赤道面上 处的一组正交基变换得到超平面内 处的正交基，并由此确定一个坐标系。每个 点处由此得到 个基矢量，剩余的最后一个基矢量为 面内的球面切向量，我们规定其指向靠近 的方向。

将该方向坐标记作 ，超平面内基矢确定的坐标记为 ，则有

即 。也就是说 点处某一个坐标的变化，仅仅会改变 处对应的一个坐标，而对其它坐标全无影响。回顾一下第4节的内容，我们发现这其实就是雅可比行列式中一行只有一个非零值的情形，我们可以直接“划掉”对应的行和列，只保留剩余部分：

变量一下子就减少到了 个呢！接下来我们就仔细地研究一下：

注意到含 的部分和不含 的部分已经完全分离开来，我们可以直接写出

其中 是与 无关的函数。从中我们可以很容易地看出 的概率密度为

而由正交旋转的性质我们有 ，因此

我们还没有求解 的具体形式，所以后面的常数因子我们没办法直接算出来。但是没关系！还记得我们前面特意算的简化问题吗？我们可以从那里面得到这个常数因子，然后再用这个因子去计算最终的结果：

两者相除就可以消掉未知常数：

现在，只要把后面的两个积分算出来，我们就成功打通了通往原问题答案的一整条道路。剩下的事情就是做一些乘法的化简而已。万里长征就快到达终点了，有没有感到很兴奋呢？

9.大结局！

分母的积分是简单的：

而分子的积分则稍显复杂一些，因为 本身就是一个变上限积分，并且无法简单地表示为 的函数。

此时我们需要利用一个很经典的公式：令 ，则有

代入分子的积分式即得

将两个积分以及第6节的结果再代入上一节末尾的式子得

把这个结果再代入第4节末尾的式子，得到

最后，利用伽马函数的倍元公式 化简得

Bingo！大功告成！我们得到了 维情形下的通式！

由于 函数在整数时取值为整数，半整数时取值为 的有理数倍，因此当 为偶数时 为有理数， 为奇数时则含有 的因子。 随着 的增大而下降，当 时有

这个问题有点意思，我以为结果会是一个很不常规的数字，但我用NumPy拿条件概率和极坐标的方法分别做了两次蒙特卡洛，得到的结果都是非常接近0.4。我最后算的解析解也是 ，而不是 。这个回答里首先我们证明这两种蒙特卡洛的方法是符合要求的，然后再用换元积分法求解析解。

我们先用数学的语言确定一些符号：让 表示单位圆盘， 表示中心在原点 的四边与坐标轴平行且长为2的一个正方形。设 ，然后 则表示 三点的外接圆。以后积分的时候这三个点对应的面积元素我也会写成 ，为方便我们简化符号。

条件概率的方法

1. 生成 ，并且独立的生成 ;
2. 如果 ，那么保留点 ；反之，则返回第一步。

假设现在我们有三个二维随机变量 由上述过程独立生成，下面我们要证明下式中的第一个等号：

上面(1)式中的等号成立，因为如果三个点构成的外接圆在单位圆内，那么这三个点必须也在单位圆内。考虑等式的左边，因为相同的理由：

(2)中的等号是对 的联合密度函数积分获得的。化简(2)中的式子后，可以把它和(1)中的表达式对比后发现是相等的。

极坐标的方法

1. 生成一个随机变量 , 独立的生成另一个随机变量 ;
2. 让 ，那么 服从单位圆盘上的均匀分布。

因为 ，所以换元之后求雅可比矩阵的行列式，我们可以得到 。

用换元法求解析解

由条件概率的部分，我们可以得到我们只需要计算：

令 外接圆圆心的极坐标为 ，外接圆半径为 ，现在 , , 三点的坐标可以表示成：

，其中 ，并且 。示性函数里的事件可以简单的表示为 ，我也是在用这个条件做蒙特卡洛模拟的时候，想到这种换元的方法。所以我们有

我用Mathematica算了一下，如下：

                f                   =                   {         r1         *         Cos         [         theta         ]                   +                   r2         *         Cos         [         a1         ],                   r1         *         Sin         [         theta         ]                   +                   r2         *         Sin         [         a1         ],                                r1         *         Cos         [         theta         ]                   +                   r2         *         Cos         [         a2         ],                   r1         *         Sin         [         theta         ]                   +                   r2         *         Sin         [         a2         ],                                r1         *         Cos         [         theta         ]                   +                   r2         *         Cos         [         a3         ],                   r1         *         Sin         [         theta         ]                   +                   r2         *         Sin         [         a3         ]};                            x                   =                   {         r1         ,                   r2         ,                   theta         ,                   a1         ,                   a2         ,                   a3         };                            b                   =                   D         [         f         ,                   {         x         }];                            Det         [         b         ]            

得到的结果是：

       r1 r2^3 (Cos[a3] (Sin[a1] - Sin[a2]) + Cos[a1] (Sin[a2] - Sin[a3]) +     Cos[a2] (-Sin[a1] + Sin[a3]))     

人工再化简一下后带入(4)式，我们有

再把(5)式代入(3)式，我们有：

做到这一步，我感觉希望很大了，因为变量分离的很好。现在先把(6)中的 积掉剩下 ，这里细节我就不写了，被积函数很漂亮很容易就积出来了，最后得到：

现在考虑(7)中的定积分（不包括最前面那个常数），令 , , ，再换元求雅可比式，(7)中的积分可以写成：

接下来我们对 的符号分类讨论。

1. 如果 非负，那么有两种情况：

• 即 ，此时有 或者 。
• 即 ，此时有 或者 。

2. 如果 为负，我们可以由 ，可得：

• ，此时有 或者 。
• ，此时有 或者 。

显然以上四种情况的两种子情况是对称的因为 。在下面的等式中反复使用被积函数的对称性： 结合(7), (8)和(9)式，我们可以得到想要求的概率是 。最后我附上蒙特卡洛模拟的代码，验证答案。

蒙特卡洛的NumPy代码

       import numpy as np  def circum_inside(coord):  ax = coord[0, 0]; ay = coord[1, 0] # Coordinates of the point A  bx = coord[0, 1]; by = coord[1, 1] # Coordinates of the point B  cx = coord[0, 2]; cy = coord[1, 2] # Coordinates of the point C   d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))   # Coordinates of the circumcenter  ux = ((ax * ax + ay * ay) * (by - cy) + (bx * bx + by * by) * (cy - ay) + (cx * cx + cy * cy) * (ay - by)) / d  uy = ((ax * ax + ay * ay) * (cx - bx) + (bx * bx + by * by) * (ax - cx) + (cx * cx + cy * cy) * (bx - ax)) / d    # Radius of the circle  r = np.sqrt((ax - ux)**2 + (ay - uy)**2)   # Determine whether or not the circumcircle is covered by the unit circle  if np.sqrt(ux**2+uy**2) + r <= 1:   return True  def is_inside(a):  return a[0]**2 + a[1]**2 <= 1   N = 10**5 # Number of Monte Carlo samples in both methods  # Generate samples by rejecting/accepting  count = 0 # Initialize the number of desired groups k = 0; data_inside = np.zeros((2, 3))  while k <= 3*N:  pt = np.random.uniform(-1, 1, 2) # draw one point  if is_inside(pt): # Determine if the point lies inside the unit circle   s = k % 3; data_inside[0, s] = pt[0]; data_inside[1, s] = pt[1]    if k % 3 == 2 and circum_inside(data_inside): count+=1   k+=1  per_rej = count / N  # Generate samples by polar coordinates ct = 0 # Initialize the number of desired groups for i in range(N):  data_polar = np.array([np.sqrt(np.random.random(3)), np.random.uniform(-np.pi, np.pi, 3)])  # Tranform polar coornates to cartesian coordinates  sp = data_polar.copy()  sc = np.array([sp[0, :]*np.cos(sp[1, ]), sp[0, :]*np.sin(sp[1, ])])  if circum_inside(sc): ct+=1  per_polar = ct / N print(per_rej, per_polar)      

两种方法分别输出的结果如下：

       0.40116 0.40151