11588 一道「小黃鴨」機率題及其有趣擴充套件 (5)
在前面的一系列文章中，我們用數學方法解決了這樣一個問題：在
維空間中的
維超球面上隨機、均勻地取
個點，它們位於同一個半超球面的機率是多少？答案是
。
不知道大家注意到沒有，在知乎上比較複雜的機率問題下面，除了用數學方法解答的答案以外，總有些答案會藉助程式設計模擬來提出猜想或者驗證結果。現在咱們也來試著編個程式，驗證一下上面的結果吧！
驗證結果的程式，需要做兩件事情：
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 函式，就負責來做第二件事 —— 判斷多個點是否位於同一個半超球面上。這件事又該怎麼做呢？
圖 5。1：以 A、B、C、D 為極點的四個半球面有交集（黑色區域）
觀察上圖。以 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)
