IA3 - 再看$n$維球上的點，淺談均勻性與不對稱性

哈伯深空圖像 (NASA/ESA/S. Beckwith(STScI) and The HUDF Team)

根據宇宙學原理 (cosmological principle)，一般相信宇宙在大尺度底下 (>Mpc) 是均勻 (homogeneous) 且各向同性 (isotropic) 的，亦即宇宙不存在中心，就像在球的表面上行走一樣，我們不會找到哪個地方在球面上有特別突出的重要性一樣，另外在宇宙各處的物質密度也應該是均等的。這樣均勻、公正且不偏頗的特性，同樣的也在統計取樣中，特別是沒有任何先驗資訊底下扮演了重要的角色。這裡我們回顧上個專題I&A2中生成任意維度單位球面上點的例子，我們以為生成的點是均勻地散布在球面上，但實際上這是一個錯覺，為什麼我們生成$(x,y)$的座標明明就是均勻的分佈，但最後會導致不均勻，或不對稱性的發生呢？

一、均勻性

讓我們回想一下上一個專題談到的均勻性，我們用二維單位圓上的點的例子來舉例比較容易想像。見下圖1，為了產生在單位圓上的點，我們可以隨機生成$(x,y)$分別在$(-1,1)$之間的點並正歸化它們，以形成單位圓上的點，而向量$(x,y)$則用於決定在極座標上的角度$\phi$。可是從圖1會很明顯的發現紫色的點它們都是對應到相同的$\phi$角，而綠色的則是角度為$\pi/2$，但由於$\phi$的位置的點可以超出單位圓，所以產生在$\phi$這個角度的單位圓上的點，實際上是比角$\pi/2$上的綠點來得多。亦即I&A2的算法2其實生成的並不是均勻分布的點，比如說在$\phi=(\frac{\pi}{4},\frac{3\pi}{4},\frac{5\pi}{4},\frac{7\pi}{4})$處的點 (二維的例子剛好在這裡的點是最多的) 會多於$\phi=(0,\frac{\pi}{2},\pi,\frac{3\pi}{2})$的點 (這裡的點是最少的)。

圖1. 藉由平均分佈生成單位原上的點，以及各個角度$\phi$所對應點出現的機率

　我們可以重新run一次I&A2的算法2，取${\tt dim}=2$，然後利用 $$\phi = \tan^{-1}\left(\frac{x}{y}\right)$$ 得到在單位圓上的點所對應的角度。圖2即是生成$10000$個點後對每個點的$\phi$作直方圖 (histogram)，不難發現在$\phi=(\frac{\pi}{4},\frac{3\pi}{4},\frac{5\pi}{4},\frac{7\pi}{4})$處的點的數量最多，而$\phi=(0,\frac{\pi}{2},\pi,\frac{3\pi}{2})$最少，完全不是均勻分佈的狀態！

圖2. 生成$10000$個二維單位圓上的點，將每個點所對應的$\phi$作直方圖

　那麼要如何避免這種不均勻呢，從圖1的概念上，也許讀者會想說，那我可以像布豐投針裡的那個例子，我丟掉任何$x^2+y^2>1$的點，我只保留那些在粉紅圓內所生成的點。是的，沒錯，這樣的話的確會保證在各個$\phi$上的點的數量是均勻的。
　但是這樣會衍生出一個問題，以圖1的狀況，假設我一開始生成$10000$個隨機分佈在正方形內的點，不過為了保持均勻性，我必須reject掉藍色區域的點，所以不是每個生成的點都能有做有效率的利用。那麼如果以三維的例子，則是立方體和單位球，由於生成的點是佈滿立方體的，但僅有那些分佈在單位球內的點會被拿來做計算，不在球內的點會被reject掉。推到任意維度的情況，我們便可以推論生成$N$個佈滿$n$維超立方體 (hyper-cube) 的點，僅有那些屬於超球 (hyper-sphere) 內的點是有用的，故接受率 (acceptance rate) $\eta$ 便是超球的體積和超立方體體積之比，即 $$\eta=\frac{V_{超球}}{V_{超立方體}}\leq 1.\tag{1}$$ 而數學家已經幫我們計算出在任意維度的球體積 $$V_{超球}(R)=\frac{\pi^{n/2}R^n}{\Gamma(\frac{n}{2}+1)}\tag{2}$$ 其中$R$是球半徑而$n$是維度，另外立方體的體積 $$V_{超立方體}(r)=r^n.\tag{3}$$ 不過這裡立方體的邊長$r$是兩倍於球的半徑$R$，可以從圖1看得出來，對單位圓而言是內切一個邊長為$2$的正方形。所以我們可以寫下 $$\eta = \frac{\pi^{n/2}}{2^n \Gamma(\frac{n}{2}+1)},$$ 而隨著維度$n$的增加，接受率$\eta$會愈來愈小，當$n\to\infty$時$\eta\to 0$，見圖3。

圖3. 接受率$\eta$隨著維度$n$的增加而迅速減少，注意$n$是整數

　這張圖告訴我們，為了均勻性的產生在單位圓或單位球上的點，如果採用拋棄$x^2+y^2>1$或$x^2+y^2+z^3>1$這種方式，那麼最終有很大一部分生成的點都會被reject掉。以二維和三維為例，$\eta$分別是$0.7$和$0.5$，即在這兩個維度上有3和5成的點會被拋棄。如果我們想產生一個$10$維球面上的點，其對應的$\eta=0.0024$，即每生成$1000$個隨機的點，僅僅有$2$顆是成功分佈在$10$維球的球面上，所以為了完整的覆蓋住$10$維球面，我們可能得要生成非常非常多的點，這不僅僅沒有效率，而且非常浪費硬體的資源，所以我們要來看一個更好的做法，不但要保持均勻性，更要減少reject的數量，最好每個點都不會被reject！

二、從高斯分佈生成各向均勻數量的點

所以從上面的例子 (特別是圖1的示意圖) 可以看出，我們之所以會得到單位圓上的點在不同$\phi$角上的分佈不均勻是由於一開始生成隨機向量$(x,y)$在某些角方向的點多於其它角的方向，為了避免這樣的不均勻性，又不想要上面提到$\eta$在高維空間急速減少的情況，如果我們可以在各個方向灑點的數量都是均勻的那就太好了！比如說，生成的點像圖4這樣。

圖4. 在$xy$平面上以原點為中心在各個方向$\phi$的點數量都是均勻的

　如果我們把在圖4幾個特定$\phi$方向的點作PDF，我們會得到圖5，我們會發現不管在哪個$\phi$的方向上，它們的PDF分佈圖都是對稱的，亦即在各個$\phi$上點的數量都是均勻的，不會有像之前一樣不對稱的情況而導致某些特定的$\phi$上產生比較多的點。

圖5. 同圖4，但將不同$\phi$方向上的點作PDF

　所以我們現在得到了，如果要各個方向產生點的數量是均勻的，則在任意維度上點的座標是以原點為中心的高斯分佈，對於多維度的高斯分 (multivariate Gaussian distribution) 佈我們有 $$p(x_1,x_2,\cdots,x_n)=\frac{1}{\sqrt{(2\pi)^n |\boldsymbol{\Sigma}|}}\exp\left(-\frac{1}{2}\mathbf{x}^T\boldsymbol{\Sigma}^{-1}\mathbf{x}\right),\tag{4}$$ 而$\boldsymbol{\Sigma}$稱作協方差矩陣 (covariance matrix)，它是用來描述不同軸之間的關聯性。由於不同軸 (例如$x$和$y$之間) 是無關聯的，所以$\boldsymbol{\Sigma}$是對角化的且其行列式 (determinant) 滿足 $$\det (\boldsymbol{\Sigma})=n^{-1}.\tag{5}$$ 至於為什麼是這個值是由於將${\tt rand}(-1,1)$轉換到高斯分佈時所作的測度變換，不過這裡不在這裡深究背後的數學，記得協方差矩陣的行列式即維度$n$的倒數。為了得到均勻分佈在$n$維球面上的點，我們必須將原來任意維度上的座標${\tt rand}(-1,1)$改成均值為$0$的高斯分佈${\tt gauss}(\sigma)$，而標準差$\sigma=\sqrt{\det (\boldsymbol{\Sigma})}=1/\sqrt{n}$，I&A2的算法2可以改寫如下：

1. 用${\tt gauss}$生成任意維度上的座標
   $\sigma=1/\sqrt{n}$
   $x_i ={\tt gauss}(\sigma),~i=1,2,\cdots,n$
2. 計算歸一化常數$r=\sqrt{\sum_{i=1}^n x_i^2}$
3. ${\tt return}:~(x_1/r,x_2/r,\cdots,x_n/r)$

上述patched過的I&A2算法2可以用python實現如下

def nsphere_patch(N,dim=2):
    """
    Generating points on the unit sphere in arbitrary dimensions.
    N: How many points we want to generate?
    dim: Dimension of the sphere.
    """
    
    # Step 1
    X = np.random.normal(0,1/np.sqrt(dim),size=(dim,N))
    # Step 2: Normalization constant
    r = np.linalg.norm(X,axis=0)
    
    return X/r

一樣以二維的例子為例，我們產生$N=10000$個隨機$(x,y)$，然後看它們在各個$\phi$上的分佈如圖6，圖6與圖2相似只是現在不對稱的情況不再，在各個角度點的數量都是接近一致的。而在這個改進的算法2裡，我們也不存在上一節提到rejection的問題，因為每個點都被利用到了！

圖6. 利用改進的算法得到在單位圓上分佈的點，在各個$\phi$方向上點的數量是接近均勻的

三、小結

原則上，要生成任意$n$維球面上的點，再不考慮各向點的數量是均勻分佈的條件下，可以用任何方式去生成隨機的點，例如I&A2的${\tt rand}$或任何自訂分佈，因為只要最後我們做了歸一化的動作，這些隨機的$(x,y,\cdots)$都會回到單位球的球面上，但是由於任意自訂的分佈，也許在不同$\phi$上會產生較多的點，而導致最後在球面上某個方向會聚集較多的點，這就是所謂的不對稱性。
　不對稱性或均勻性都有其重要的價值在。例如在Bayesian analysis裡，我們要對後驗機率分佈 (posterior probability distribution) 做取樣，對於高密度的區域，我們應該提高取樣的數量，對於其它低密度的地方，則應該降低取樣的數量。這時候不對稱性就扮演了重要的角色，因為我們不能平等的對待後驗機率分佈的每個地方，對那些不重要的地方應該降低或忽視，以避免硬體花費太多時間在徒勞無功的地方著墨。
　而當我們想要測定一種新藥對人體的影響時，在什麼先驗資訊都沒有的狀態，我們不應該偏好特定的種族、性別或年齡，而應該以均勻隨機採樣的方式做投藥，我們才能得出整體性的概念，之後若發現藥效特別對某個項目反應結果高於平均，我們才應該進一步做不對稱的分析以避免在研究初始就產生了偏誤 (bias)。


Source on github: I&A3_homogeneity.ipynb