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

留言

張貼留言

這個網誌中的熱門文章

物理雜記3 - Ising模型與模擬退火

圖片
Ising模型主要是用來描述物質的鐵磁性 (ferromagnetism)，是由德國物理學家Ernst Ising在1920左右提出了處理一維的解析解。該模型包含了每個晶格上原子的自旋 (spin) 參數，兩兩相鄰的原子彼此間可以是偏好處在共同的態上 (都是自旋向上)  或是處在相反的態上 (一上一下)，而外部的磁場會對該晶格上的自旋狀態施加能量，導致其偏向某一個方向。在這篇文章裡，我們不重複Ising當年複雜的數學解析解，而且該解也只能適用在1D的狀態，對於2D以上的系統目前並沒有解析解，所以我們將要另闢蹊徑，採用MCMC模擬的方式來決定一個系統終態每個晶格上自旋的排列情形。 一、Ising模型簡介 考慮一個由圖1所構成的晶格系統，晶格上的點都代表一顆原子，而黑白則表原子的自旋狀態，我們用黑表示向下而白表示向上，所以這個binary系統可以用更簡潔的\(\pm 1\)來表示每個原子可能的自旋 (即這個系統的樣本空間\(\Omega=\{+1,-1\}\))，另外為了方便，我們採用物理學的習慣將自旋表示為希臘字母\(\sigma_{n}\)，下標\(n\)可以想像成在 n 處原子的自旋態。原則上討論的都是一個由費米子 (fermions) 構成的系統，所以自旋應該用\(\pm \frac{1}{2}\)來表述，但在計算機上我們用\(\pm 1\)比較乾脆且簡潔。另外補充一下冷知識自然界已知的物質幾乎由費米子構成，僅有媒介作用力的粒子是玻色子 (bosons)，它們的差別在於費米子有半整數 (即不能約化成整數的分數) 自旋而玻色子則是整數自旋。 圖1. 晶格系統，黑白點分別代表該處原子的自旋狀態，黑表向下，白表向上。紅色的箭頭指出中心原子的4位鄰居\(m\) 　Ising模型假定了每個晶格上的原子的自旋都會受其鄰居 (neighbors \(\mathcal{N}\)) 的狀態影響，例如圖1中間白色的原子雖然他在這一刻是處在\(+1\)的狀態，但根據其週邊四個鄰居 (分別是3白1黑) 的況，它仍有一定的機會會在下一刻翻轉自旋的狀態。 　對一個給定位置 n 的原子，它和它的周圍鄰居構成的局部區域 (local area) 形成的一個系統的能量由Ising模型給出 \[ E_n = -\sum_{m\in \mathcal{N}} J
閱讀完整內容

PS7 - 小專題: 為何高斯分佈擁有最大熵?

圖片
高爾頓板，或稱做quincunx，是由Galton爵士 (1822-1911) 所發明用來驗證中央極限定理 (CLT) 的裝置。它也是第一位討論均值回歸 (regression to the mean) 現象的科學家，但同時亦是優生學 (eugenics) 這樣錯誤觀念的推廣者 在 上篇雜記 的結尾我們提到了對於一個趨於穩態的系統擁有最大的熵，而這也是統計力學的直接應用，用物理的角度來講便是當系統平衡時 Gibbs 自由能不再變化 (\(\Delta G =0\))，而在此時會對應到系統的熵有全域極大值 (global maximum)。故大自然很多現象都已經處於穩態的狀況下，我們對這個現象的母體做抽樣，所得到的觀測樣本應也有最大熵，則描述這些觀測樣本最好的分佈便是符合該條件下擁有最大熵的機率分佈。而指數族的分佈符合最大熵的特點，我們在這個小專題內便嘗試證明高斯與二項式分佈在連續和離散的隨機變量下的條件擁有最大熵的性質。 一、機率熵 在 上篇雜記 裡裡我們探討了夏儂熵，或者稱作資訊量\(I\)的意含，所以一個系統內 (我們這裡就不再用系綜這種彆扭的歷史名詞了) 每個成員出現的機率是\(p_i\)，則對該系統的熵可以用夏儂熵公式來表示 \[ I(p) = -\sum_i p_i \ln p_i. \tag{7.1} \] 所以對某一組 連續 機率分佈函數\(p(x)\)，其中\(x\)是隨機變量，則將(7.1)的求和符號改成積分後有 \[ I(p) = -\int p(x)\ln p(x) dx, \tag{7.2} \] 稱為機率分佈\(p(x)\)的 機率熵 ，而積分的上下限則看隨機變量的範圍，例如高斯分佈就是\(\pm \infty\)。如果今天我們有兩組機率分佈\(p(x)\)和\(q(x)\)，我們便可以定義它們之間的 相對熵 ，或稱作 KL散度 (Kullback-Leibler divergence)，定義為 \[ D_{KL}(q,p)\equiv\int q(x)\ln \left(\frac{q(x)}{p(x)}\right) dx = I(q,p)-I(q)\geq 0, \tag{7.3} \] 是 正定的 (positively defined)，其中 \[ I(q,p)=-\int q(x)\ln p(x) d
閱讀完整內容

IA9a - K-means分群與EM算法: 理論

圖片
市場內分門別類的產品 原則上 K-means分群 是 EM算法 (expectation-maximization algorithm) 的一個特例，相較於EM從似然函數 (likelihood) 的方式出發並導出不同資料點對不同群 (或稱聚類、簇，cluster) 的貢獻度，K-means單純用歐氏距離來分類，雖然這加快了K-means算法的計算效率，但顯然的這麼naïve的方法會產生明顯的缺陷。但仍讓我們會先從簡單的K-means開始，並探討K-means算法的結果在某些情況下和人類肉眼辨識的結論會產生明顯的歧義，接著再推進到通用的EM算法，除了 python 的demo外，我也會證明在數學上如何從EM得到K-means的極限 (僅限高斯混合的情況)，而在這個part內，讓我們先從數學理論的角度來面。 一、K-means分群 今天如果我有一大坨資料共120個點如圖1所示，從資料的散佈趨勢來看，很明顯的分成兩群，那我們是不是可以有一個演算法來自動化將資料點如同我們臆測的一樣分成兩袋？一個最簡單的發法便是使用K-means分群，而這個大寫K表示群的數量，以圖1而言，我們猜測K=2，所以是個2-means分群問題。 圖1. 資料散佈圖 (普魯士藍) 與初始化的紅綠兩點 　所以為了達成分群，我們需要先初始化兩個點，這兩個點代表兩群的中心。事實上，初始化的位置不必要一定剛好在兩群的中心，可以是任意挑選的地方，但為了計算方便，我們會將起始化的兩點開始在靠近資料點的平均值附近，如同圖1中紅綠兩點。我們將這兩個初始化的點叫作\(p_k\)其中\(k=1,2\)表示第\(k\)群。另外將資料點的集合用\(\mathbf{D}=\{d_1,d_2\cdots d_n \cdots d_{120}\}\)表示，所以我們計算每個點\(d_n=(x_{d_n},y_{d_n})\)和\(p_k=(x_{p_k},y_{p_k})\)的歐氏距離 \[ \ell_{n,k}=\sqrt{(x_{d_n}-x_{p_k})^2+(y_{d_n}-y_{p_k})^2}. \tag{1} \] 若該點和\(p_1\)的距離小於和\(p_2\)的距離 (\(\ell_{n,1}<\ell_{n,2}\))，則我們說\(d_n\)比較靠近\(p_1\)且應該
閱讀完整內容