SIC ITUR AD ASTRA

馬可夫過程 (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\) 上述的算法在第2步就是去篩出隨機生成的\(N\)個點中，哪些是符合半徑等於\(1\)的狀況 (即這些點在紅色的圓上)，如果不等於\(1\)則將之刪去。我不做計算，但我可以打包票，上面的算法接受率 (acceptance rate) 是超級低，生成數以萬計的點可能都不見得有一個是符合半徑等於\(1\)的條件。 　但換個角度來看生成的點，點座標\((x,y)\)其實也代表著向量，既然單位圓有半徑為\(1\)，我們只要將每個點的向量都做歸一化的動作，我們就可以得到單位向量\((x^\prime,y^\prime)\) (即向量距原點的長度為\(1\))，所以可以保證歸一化後的單位向量一定會在圖1的紅色圓上，而未歸一化前的點座標則用來決定在極座標上的角度\(\phi\)，即 \[ \phi=\tan^{-1}\left(\frac{y}{x}\right). \] 示意圖如圖2。這樣不僅保證所有的點都在單位圓上，又由於點是隨機生成的，在圓周上的位置\(\phi\)也是在\(0\)和\(2\pi\)之間均勻分佈。 圖2. 歸

Mathematica®生成400根針的隨機投放示意圖 布豐投針問題 (Buffon's needle problem) 是這樣敘述的： 有一房間的地板由數片平行並列的木板所構成，板與板之間溝紋的距離皆相同，如上圖棕色隔線的間距。試問隨機扔一針，該針觸及溝紋的機率為何？ 這項問題主要由18世紀的布豐伯爵所提出的，並由義大利數學家Mario Lazzarini在1901年實地進行拋針實驗，共拋出了3,408根針(🥴)並計算出\(\pi\)的近似值約為\(355/113\)，參見 維基百科 。在我們繼續追本溯源這個歷史問題以前，讓我們先看看當代最常見估計圓周率的方法。 一、灑點估計法 這個方法概念非常的簡單，我們將產生\(N\)個座標\((x,y)\)隨機散佈在\(0\)和\(1\)之間的點，所這些點將被包在以原點為中心，邊長為2的正方形的第一象限中 (first quarter)。接著我們計算出每個點距原點的長度\(d=\sqrt{x^2+y^2}\)，若\(d<1\)則紀錄，反之則丟棄，最後我們紀錄了\(n\)個點，這些點都被包在以原點為中心，半徑為\(1\)的\(1/4\)圓內，見圖1。 圖1. 在第一象限中藍色正方形的面積為\(1\)而紅色區域的面積為圓的\(1/4\) 　現在我們相信這些處於單位圓內點的數量和全部點的數量的比值應等於那\(1/4\)圓面積和正方形面積之比，即 \[ \frac{n}{N}=\frac{A_圓}{A_方}.\tag{1} \] 而\(A_正 = 1^2\)和\(A_圓=\pi 1^2/4\)，所以我們便得到了 \[ \pi = 4\frac{n}{N},\tag{2} \] 正是\(\pi\)的估計。我們預測若\(N\)愈大，則其估計值愈接近真正的圓周率。另外\(\pi\)可以解析表示成 \[ \pi = \int_{-1}^{1}\frac{dx}{\sqrt{1-x^2}}. \tag{3} \] 上式是由現代分析學之父卡爾·魏爾斯特拉斯 (Karl T. W. Weierstraß, 1815-1897) 在1841年所定義的。 　我們接著可以透過 rand 來實現這個簡單的蒙地卡羅模擬，腳本如下： import numpy as np import matp

藝術化的元胞自動機 (cellular automaton)，取自 Softology's blog 在今年初我在自己的這個部落格開了一個關於機率與統計 (P&S) 的系列，是包含了我自己學習古典機率的心得，後來再加上由於計算機效能長足進步後日漸興起的貝氏推論 (Bayesian inference) 的應用。在研習的過程當中，我感到許多的內容與統計力學和計算機理論是緊密相連的 (必然如此，這些學問都是嘗試去理解大自然的方法，也許從不同的立足點出發，但最終它們想攻剋的都是同一做山頭)，所以不僅僅重新讀了以前懵懵懂懂的統計力學，還對資訊理論以及這些算法做了不少功課，我看見所有的學問都是一脈相承、自有而擁有的，包含最近這幾年不斷引起大眾耳目的人工智慧，以及在推論領域中已穩紮穩打許久的資料科學。 　所以雖然P&S我還沒認為應該結束這個系列，但我想再開一個資訊與算法 ( I nformation & A lgorithm, IA ) 的專題。不像PS系列，在這裡我應該會著重在我研讀這類資料的時候遇到的問題如何解決，以及已經解決的問題，我能夠用簡單的程式腳本來實現，所以這個系列每篇都如同一個小專題一樣，而不是我循序漸進的學習過程。 　這類的問題像是圓周率\(\pi\)的估計，除了耳熟能詳的用computer隨機產生一堆點然後計算有多少落在半徑為1的圓內來估計外，我們也要提及幾乎被遺忘的18世紀法國布豐伯爵 (Georges-Louis Leclerc, Comte de Buffon, 1707-1788)，他發明了布豐投針 (Buffon's needles) 估計\(\pi\)的蒙地卡羅方法 (Monte-Carlo method)。不像是傳統邏輯符號的推論，蒙地卡羅模擬 (Monte-Carlo simulation) 在於能用大自然處理事情的方式來估計未知的參數，或是推論可能發生的情況，而當代硬體計算效能之強大，即使在一台簡便的筆記型電腦上也足以處理這些運算，甚至得出精細的結果，所以電腦是實現蒙地卡羅模擬非常好的工具。除此之外，這個系列也包含了我遇到一些有趣的演算法，除了實踐它們外，我希望能逆向找出這類模擬背後的數理邏輯。 　至於原本的P&S系列，我最想談的還是Bayesian analysis，

當我們說處在海平面為0的時候，亦即以海的平面為基準當作高度為0 (基平面)，例如玉山高度3,925公尺便是高於海面3,925公尺。不過當我們望向一望無際的海面，所有的水分子都處在這個為0的基平面上嗎？原則上這些水分子會以這個基平面為準做上下的震盪，而若我們把所有的水分子在某一瞬間偏離這個平面的高度 (含低於基平面的高度，取負值) 總和的平均，我們會得到非常接近0的值。所以所謂的海平面是一個平均的概念，而真實的水分子會因為熱運動不斷上下偏離這個平均，但其總和仍為0，而這個偏離的行為我們便稱作 漲落 (fluctuation)。 　我們回到杯水的例子，給定宏觀的狀態\(N\)、\(V\)和\(T\)後雖然水分子的平均能量是\(\bar{E}=3k_B T/2\)，但每個水分子真正的能量\(E\)卻是是落在\(0\)到\(\infty\)之間，所以其能量的分佈可用正則分佈來描述。而在眾多的水分子當中，即便都是能量\(E\)也可能有不同的內在狀態 (intrinsic properties)，像是自旋等等，對於這種多個分子共享相同能量確有不同內在狀態的情況我們稱做 簡併 (degenerate)。考慮了簡併後的正則分佈可以改寫成 \[ p(E)=\frac{g(E)e^{-\beta E}}{\mathcal{Z}},\tag{7.1} \] 其中\(\mathcal{Z}=\int g(E) e^{-\beta E}dE\)為配分函數而\(g(E)\)表示在能量為\(E\)時有多少水分子處於簡併的狀態，它是能量\(E\)的函數，也叫做 態密度 (density of state)。這裡用積分表示由於分子數量眾多，所以能量間隔可以考慮成連續譜的狀態。 　原則上\(e^{\beta E}\)是個隨能量遞減的單調函數，而態密度\(g(E)\)則隨能量遞增，一個遞增一個遞減相乘在一起會產生一個極值，而極值的位置則由 \[ \left. \frac{\partial \{g(E)e^{-\beta E}\}}{\partial E}\right|_{E=E^*}=0 \to \left. \frac{\partial \ln g(E)}{\partial E}\right|_{E=E^*}=\beta \] 決定。另外當\(E\)給定後，這個系統的熵正比於簡併態

高爾頓板，或稱做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

在前面的幾篇文章裡面我討論到一些機率的問題時提到了Shannon entropy來決定least informative prior的技巧，在這篇雜記裡我想要對Shannon entropy做一些探討。Shannon entropy源自於C. E. Shannon在資訊理論上的工作，特別是引入了統計力學對一個物理系統微觀狀態上的觀點。所以在這裡，我想從統計力學的角度推導出Shannon entropy或者稱作資訊量\(I\)的關係式\[I = -\sum_i p_i \ln p_i\]並且探討一些關於資訊理論的簡單問題。 一、系統的狀態 統計力學主要在是從微觀的角度來推論出巨觀的所有現象學上的熱力學特性，這些巨觀的特性諸如壓力\(P\)、內能\(U\)、熵\(S\)和定容或定壓比熱\(C_P\)、\(C_V\)...等等。我們這裡並不去討論如何推導這些特性，我們只要知道為了得出這些特性，我們必須了解系統現在到底處在什麼樣的狀態上。 　舉例來說，桌上裝滿水的杯子，杯子體積是\(V\)且有\(N\)個水分子 (已知水的質量便可以除上水分子的莫爾質量乘上Avogadro常數得到總水分子數\(N\))，假設它已經和室溫\(T\)達成熱平衡，那麼我從中取出任意個水分子，試問該水分子的動能為何？下方我們皆以能量\(E\)來代稱動能。學過普通物理學也許會認為每個水分子的能量是\[\bar{E}=\frac{3}{2}k_B T \tag{7.1}\]其中\(k_B\)為Boltzmann constant。答案是的沒錯，但注意到我們取的是\(\bar{E}\)即平均能量，表示若我測量了無數個水分子的能量，那麼平均來講每個水分子會帶有\(\bar{E}\)，不過對任意個水分子而言，他所帶的能量是在0到\(\infty\)中隨機賦值，這表示所有測量的結果\(E\)是一組統計上的隨機變量 (random variable)。所以這個問題便回到了在一給定的杯水的巨觀狀態\(N\)、\(V\)和\(T\)後，任意取出一個水分子的能量\(E\)的值最有可能為何？而這可以由\(E\)的機率分佈函數來決定。熟悉統計力學的人便可以馬上領會到，這個杯水對應到的是 正則系綜 1 (canonical ensemble)，而每個水分子則是構成該系綜的成員 (members) 之一

Stirling 近似公式 (Stirling's approximation)，或者表示成\[\ \ln n! \simeq n\ln n - n \]當\(n\)極大的時候，Stirling近似和精確值如圖1所示。 圖1.  \(\ln \lambda !\)精確值 (藍)、Stirling 近似公式 (橘) 與利用 Poisson 分佈所推導含NLO項的值 (綠) 　關於這條近似公式我們有一個簡單的推導，可以從 Poisson 分佈著手，其表示式如下\[ P(n|\lambda)=e^{-\lambda}\frac{\lambda^n}{n!}. \]而當\(\lambda\)很大 (\(>10\)) 的時候，Poisson分佈在\(n\sim \lambda\)附近可以被一個平均值和變異數皆為\(\lambda\)的高斯分佈來近似，即\[ e^{-\lambda}\frac{\lambda^n}{n!}\simeq \frac{1}{\sqrt{2\pi \lambda}}e^{-\frac{(n-\lambda)^2}{2\lambda}} \]取\(n=\lambda\)我們有\[\begin{eqnarray*} e^{-\lambda}\frac{\lambda^{\lambda}}{\lambda!}&\simeq&\frac{1}{\sqrt{2\pi\lambda}}\\\Rightarrow\lambda!&\simeq&\lambda^\lambda e^{-\lambda}\sqrt{2\pi \lambda} \end{eqnarray*} \]等式兩邊取\(\ln\)後我們得到\[ \ln \lambda! \simeq \lambda \ln \lambda -\lambda + \frac{1}{2}\ln 2\pi \lambda \]即 Stirling 近似公式！其中最後一項\(\frac{1}{2}\ln 2\pi \lambda\)為 next-to-leading-order 修正 (NLO correction)，我們從 Poisson 分佈出發的解在\(\lambda \gtrsim 0.2\)時似乎比原來的近似公式還要來的更接近\(\l