IA2 - \(n\)維球上的點
馬可夫過程 (Markov process) 生成三維球面上的點，此時點的位置不是均勻分佈，而是猶如醉漢走路 (drunkard's walk) 的情況。馬可夫過程在針對自然界運作或是社會結構演變的模擬中扮演著核心的角色！
一、二維圓盤的case
這次的專題非常簡單，我們先從二維的情況開始，二維球面上的點即為是分散在圓周上的點。我們知道如何產生散布在\((x,y)\)上的點 (圖1)，但如何將之收縮到單位圓的圓周上呢？ |
| 圖1. 均勻散布在邊長為1的正方形內的藍點以及單位圓 (紅) 的範圍 |
在圖1我們產生了1000個點，以及對應單位圓圈起來的面積。這次我們並不是要估計\(\pi\)，而是要找出哪些點分佈在紅色的圓周上，亦即產生出散布在單位圓上的點。你可能想說可以用以下的方式來實現 (以二維的case為例)：
- 產生\(N\)的點的座標分別在\((-1,1)\)之間
\(x={\tt rand}(-1,1)\)
\(y={\tt rand}(-1,1)\) - \({\tt if}~(x^2+y^2=1):~{\tt return}~(x,y)\)
\({\tt else:~return}~0\)
但換個角度來看生成的點，點座標\((x,y)\)其實也代表著向量，既然單位圓有半徑為\(1\)，我們只要將每個點的向量都做歸一化的動作，我們就可以得到單位向量\((x^\prime,y^\prime)\) (即向量距原點的長度為\(1\))，所以可以保證歸一化後的單位向量一定會在圖1的紅色圓上，而未歸一化前的點座標則用來決定在極座標上的角度\(\phi\)，即 \[ \phi=\tan^{-1}\left(\frac{y}{x}\right). \] 示意圖如圖2。這樣不僅保證所有的點都在單位圓上，又由於點是隨機生成的，在圓周上的位置\(\phi\)也是在\(0\)和\(2\pi\)之間均勻分佈。
 |
| 圖2. 歸一化前與歸一化後的點，以及在極座標上的角度\(\phi\) |
所以我們改寫上述的方式得到了算法1：
- 產生\(N\)個點的座標分別在\((-1,1)\)之間
\(x={\tt rand}(-1,1)\)
\(y={\tt rand}(-1,1)\) - 計算歸一化常數\(r=\sqrt{x^2+y^2}\)
- \({\tt return}~(x/r,y/r)\)
def unitcircle(N):
"""
Generating points on the unit circle.
N: How many points we want to generate?
"""
# Step 1
X,Y = np.random.uniform(-1,1,size=(2,N))
# Step 2: Normalization constant
r = np.sqrt(X**2+Y**2)
return X/r,Y/r
 |
| 圖3. 透過算法1生成100在單位圓上的點 |
二、\(n\)維球面的case
原則上由算法1的方式，我們可以很簡變得拓展到三維球面的狀況，只需要在算法1的第1步多加上一個\(z\)軸的位置，並且讓第2步歸一化常數變成計算\(r=\sqrt{x^2+y^2+z^2}\)即可，結果如圖4。 |
| 圖4. 單位球球面上均勻分佈的紫色點點 |
所以與其再寫一遍三維的算法，我們直接推廣到任意維的case，我們會發現一維 (點)、二維 (線) 和三維 (面) 都可以直接從這個算法導出。我們用\(x_i\)其中\(i=1,2,3,\cdots,d\)表示\(d\)維空間的第一個到第\(d\)個維的卡氏座標，所以算法2如下：
- 生成\(N\)個\(d\)維的點座標
\(x_1={\tt rand}(-1,1)\)
\(x_2={\tt rand}(-1,1)\)
\(\qquad\quad\vdots\)
\(x_d={\tt rand}(-1,1)\) - 計算歸一化常數\(r=\sqrt{\sum_{i=1}^d x_i^2}\)
- \({\tt return}~(x_1/r,x_2/r,\cdots,x_d/r)\)
def nsphere(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.uniform(-1,1,size=(dim,N))
# Step 2: Normalization constant
r = np.linalg.norm(X,axis=0)
# The following returns:
# 1st row: all N points' positions on x-axis
# 2nd row: all N points' positions on y-axis
# ...
# d-th row: all N points' positions on d-axis
# It can be projected to traditional (x,y,z,...,d) coordinate by transpose T
return X/r
三、小結
在這裡我們提供一個簡便生成任意維度單位球面上點的算法，而且也藉由python便利的性質，雖然我們用到許多駭人的術語，像是高維球之類的😩，但實際上操作起來並沒有很困難，此份章節我們也幾乎不用到數學，完全體現了第一原理模擬，從最簡單的概念出發，我們可以逼近自然界運作的方式，當然生成球面上的點和自然運作沒有什麼關聯就是😅。
除了均勻分佈外，如本專題一開始透過馬可夫過程生成球面上的點，馬可夫過程表示每一個步驟所得到的結果會與它的上一步有關，但僅僅如此而已，與更早的過程完全無關，是無記憶性的 (memory-less)。但由於必須給定初始條件 (可以透過隨機生成初始條件)，我們看到點並不是均勻隨機分佈在球面上，而是會呈現一幅醉漢走路的狀態。馬可夫過程在許多模擬自然或社會運作的計算中扮演著核心的角色，有機會我們再回到這個議題上。
Source on github: I&A2_nSphere.ipynb
Source on github: I&A2_nSphere.ipynb
留言
張貼留言