> ### 前言 > 你清茶園不是人待的地方！ > 裡面的個個都是人才，說話又好聽——就是我太菜了啥也聽不懂，這次期中還考的賊**爛，太讓人鬱悶了。 > 最近課上講這個馬爾科夫鏈蒙特卡羅方法，我也學得一塌糊塗。這時我猛然想起了自己的部落格園密碼（霧），來更個部落格吧。 **[Warning]** 本人數學水平差勁，下文用詞不嚴謹、缺少部分證明，請酌情閱讀。若出鍋，歡迎指正。 ## 啥是馬爾科夫鏈？ 馬爾科夫鏈(Markov Chain)，簡單來說就是一個用來隨機遊走的有向圖，每條邊(u, v)的邊權$p_{uv}$代表“當前在u，下一步走到v”的概率，顯然需要 $$p_{uv}\ge 0, \sum_{v}p_{uv}=1.$$ 下文中我們假設這個有向圖是**強連通的**，即任取兩個點u和v，都存在從u到v的、邊權都大於0的路徑（當然從v到u的路徑也要存在）。 馬爾科夫鏈（也就是這個隨機遊走過程）的美妙性質在於它**收斂**。怎麼個收斂法呢？ 設這個圖有n個點，令$n$維行向量$\mathbf p(t)$表示隨機遊走了t步之後的概率分佈（在時間t分別位於每個點的概率），$\mathbf p(t)_i$ 就是第t步到點i的概率。初始狀態$\mathbf p(0)$是隨便欽定的。 再定義一個量$\mathbf a(t)$，名叫“長期平均概率分佈(long-term average probability distribution)”， $$\mathbf a(t) = \frac1t \sum_{i=0}^{t-1} \mathbf p(i).$$ 顧名思義，就是把前$t$個時間點的概率分佈取個平均。 ### 定理（馬爾科夫鏈基本定理） - 存在唯一概率分佈向量$\mathbf \pi$使得$\mathbf \pi P = \mathbf \pi$。（這個P就是由$p_{xy}$構成的矩陣。） - $\lim_{t\rightarrow \infty} \mathbf a(t)$ 存在，而且就等於$\mathbf \pi$。 (證明待補充……放假再敲？） 這個定理就很令人開心了——不管欽定初始狀態$\mathbf p(0)$的你是歐皇還是非酋，只要遊走足夠多步，$\mathbf a(t)$肯定會收斂到唯一的答案。 誒？為啥要定義這個a(t)呢？p(t)自己它不收斂麼？ 還真不收斂。考慮一下$w_{12} = 1, w_{21}=1, p(0)=(1, 0)$。在這個圖上游走就是在兩個點之間反覆橫跳，p(t)是不會收斂的！但是a(t)就能收斂到$(\frac12, \frac12)$。 ## 啥是蒙特卡羅方法？ 蒙特卡羅(Monte Carlo)是誰？不是誰，這是個賭場的名字 (=_=|||)。取這個名字大概只是因為……它是隨機的，賭場也是隨機的？不得不說這個洋氣名字實在太勸退了。 其實蒙特卡羅方法就是……抽樣估計。小學學的撒豆子求面積啊，蒲豐投針計算圓周率啊，都可以視作蒙特卡羅演算法。 ## 啥是馬爾科夫鏈蒙特卡羅方法(MCMC)？ ~~I have a Markov Chain, I have a Monte Carlo, ah, Markov Chain Monte Carlo!~~ 我們已經知道，馬爾科夫鏈是個好東西，它保證能收斂到某個概率分佈。現在我們**已知一個概率分佈**，想要構造出相應的馬爾科夫鏈。這有什麼用呢？有些時候，概率分佈長得比較複雜，直接根據它生成隨機變數非常困難，例如想在一個高維空間中的凸多邊形中隨機取一個點——這時候我們就可以先構造出這個概率分佈對應的馬爾科夫鏈，然後在上面隨機遊走來取點。 不同的馬爾科夫鏈可能收斂到相同的概率分佈，但是收斂的速率有快有慢。在應用中，我們肯定希望馬爾科夫鏈收斂得快一點，少遊走幾步就能獲得和想要的概率分佈差不多的結果。下面是兩個比較好用的構造方法。 ### Metropolis-Hasting 演算法 有的時候又叫Metropolis演算法，反正是個一個名字挺長的演算法。以下簡稱MH演算法。 和它不好記的名字相反，MH演算法描述起來非常簡單、非常符合直覺！ 首先你要有一個**無向連通圖**。先不管這個圖是怎麼構建的，很有可能是出題人送給你的（假如你在做相關的練習題的話）。 設概率分佈是$\mathbf p$，點$i$的概率是$p_i$。假如你當前在點$i$： - 先從與點$i$相鄰的所有點中，等概率隨機取一個點$j$; - 如果$p_j \ge p_i$，則立刻走到$j$； - 如果$p_j < p_i$， 隨機一下，有$\frac{p_j}{p_i}$的概率走到$j$，否則待在原地不動。 這樣$p_j$更大的$j$更有可能被走到，很符合直覺吧！ 把上面的文字策略用公式表示一下，馬爾科夫鏈中點$i$到點$j$的邊權$p_{ij}$就是 $$p_{ij} = \frac1r \min(1, \frac{p_j}{p_i}),$$ 而剩下的概率都分給了$p_{ii}$， $$p_{ii} = 1 - \sum_{j\ne i} p_{ij}.$$ 那麼這個馬爾科夫鏈到底能不能收斂到$\mathbf p$這個概率分佈呢？ #### 引理 令$\mathbf \pi$為一個概率分佈，$P$為馬爾科夫鏈的邊權矩陣，如果對任意兩點$x,y$，均有 $$\pi_x p_{xy} = \pi_y p_{yx},$$ 則$\mathbf \pi$就是馬爾科夫鏈收斂到的概率分佈。 證明：若$\pi_x p_{xy} = \pi_y p_{yx}$，那麼列舉$y$，對等式兩邊分別求和，就有$\pi_x = \sum_y \pi_y p_{yx}$，即$\mathbf \pi = \mathbf \pi P$。根據上面的馬爾科夫鏈基本定理可知收斂到$\mathbf \pi$。 把MH演算法構造的馬爾科夫鏈代入這個引理，可知確實收斂到$\mathbf p$！ （課上講了一個應用MH演算法破譯密碼的有趣例子，我來不及敲完了，可以看看[這篇文章](https://math.uchicago.edu/~shmuel/Network-course-readings/MCMCRev.pdf)的開頭。） ### Gibbs 取樣 Gibbs 取樣(Gibbs Sampling)是另一種MCMC，它使用的是一種特殊的無向圖：每個點對應$d$維空間中的一個格點$\mathbf x = (x_1, x_2, \cdots, x_d), x_i\in\mathbb Z$，如果$\mathbf x$和$\mathbf y$只有一個個座標不同，則二者之間有連邊。這樣形成的就是一個類似$d$維網格（lattice）的結構，但每條（與座標軸平行的）直線上的點都形成一個**團**。  （一個用來Gibbs取樣的無向圖，圖片來自Blum, Hopcroft, Kannan, "Foundations of Data Science"） 對於相鄰的兩點$\mathbf x$和$\mathbf y$，不失一般性，假設它們第一維座標不同（$x_1 \ne y_1$）而剩下的座標都相同。那麼就令 $$p_{\mathbf{xy}} = \frac1d p(y_1|x_2, x_3, \cdots, x_d).$$ 而剩餘的概率留給$p_{\mathbf{xx}}$, $$p_{\mathbf{xx}} = 1-\sum_{\mathbf y \ne x}p_{\mathbf{xy}}.$$ 驗證： $$p(x_1|x_2, x_3, \cdots, x_d) p_{\mathbf{xy}} = p(y_1|x_2, x_3, \cdots, x_d) p_{\mathbf{yx}},$$ 故根據條件概率公式有 $$p(x_1) p_{\mathbf{xy}} = p(y_1) p_{\mathbf{yx}},$$ 所以概率分佈確實收斂到$\mathbf p.$ ## ε-混合時間(ε-mixing time) （我沒查到這個ε-mixing time怎麼翻譯成中文……有大佬知道的話評論區告訴我一下qaq） 儘管我們知道馬爾科夫鏈早晚會收斂的，但它到底早收斂還是晚收斂，對我們很重要。為了衡量收斂的快慢，定義ε-混合時間(ε-mixing time)為$t$的最小值，滿足：對任意初始狀態$\mathbf p(0)$,$\|\mathbf a(t) - \mathbf \pi(t)\|_1 < \epsilon$（這個範數就是曼哈頓距離那種，各維度距離直接相加的和）。為了保證演算法的效率，我們想知道這個ε-混合時間的一個上界。 一個上界是由導率(conductance)確定的。啥是導率呢？顧名思義，就和物理上的電導率差不多。對於節點集合$S$，令$\pi(S) = \sum_{x\in S} \pi_x$。令$\bar S$為$S$的補集，若$\pi(S) \le \pi(\bar S)$，那麼導率$\Phi(S)$就是 $$\Phi(S) = \frac{\sum_{(x, y)\in(S, \bar S)}\pi_xp_{xy}}{\pi(S)} = \sum_{x\in S}\frac{\pi_x}{\pi(S)}\sum_{y\in\bar S}p_{xy},$$ 即“已知當前在$S$中，下一步跳出$S$的概率”。 直覺上，如果$\Phi(S)$都很大，那麼我們可以在圖上來去自如，收斂就比較快；反之，如果某個$\Phi(S)$很小，一旦進入$S$就被困住、很難逃脫，收斂就可能很慢。 因此，定義整個馬爾科夫鏈的導率為 $$\Phi = \min_{S\subset V, S\ne \emptyset} \Phi(S).$$ ### 定理 任何無向圖隨機遊走的ε-混合時間有上界 $$O\left(\frac{\ln(1/\pi_{min})}{\Phi^2\epsilon^3}\right).$$ （證明……敲不完……下次也不一定） ### 例子 一條$n$個節點順序組成的鏈，$\mathbf \pi = (\frac1n, \frac1n, \cdots, \frac1n)$。 可以構造馬爾科夫鏈$p(1, 1) = p(n, n) = \frac12, p(i, i+1) = p(i + 1, i) = \frac12, \forall 1 \le i \le n-1$。 最難逃脫的$S$就由前半條鏈（後半條也行）構成的集合，$\pi(S) = \frac12$, $$\Phi = \Phi(S)= \frac{\frac1n \cdot \frac12}{\frac12} = \frac1n, $$ ε-混合時間的上界就是 $$O\left(\frac{\ln n}{\left(\frac1n\right)^2 \epsilon^3}\right) = O\left(\frac{n^2\ln n}{\epsilon^3}\right).$$ --- **寫作業去啦 未完待續_(:з」∠)_**