本文主要摘抄整理自Rickjin的《LDA數(shù)學(xué)八卦》
統(tǒng)計模擬中一個很重要的問題就是給定一個概率分布$p(x)$,我們?nèi)绾卧谟嬎銠C中生成它的樣本比庄。一般而言，均勻分布$$Uniform(0,1)$$的樣本是相對容易生成的：我們可以通過隨機數(shù)生成器生成偽隨機數(shù)序列帝簇，作為平均分布的樣本。一些常見的概率分布，可以基于$$Uniform(0,1)$$的樣本生成派撕，比如正態(tài)分布可以通過著名的Box-Muller變化得到Box-Muller變換如果隨機變量$$U_1,U_2$$獨立且$$U_1,U_2~Uniform[0,1]$$,令$$Z_0=\sqrt{-2ln U_1}cos(2\pi U_2)$$
$$Z_1=\sqrt{-2ln U_1}sin(2\pi U_2)$$則$$Z_0,Z_1$$獨立且服從標(biāo)準正態(tài)分布隧熙。其它幾個著名的連續(xù)分布包括指數(shù)分布片挂、Gamma分布、t分布贞盯、F分布音念、Beta分布、Dirichlet分布等等躏敢，也都可以通過類似的數(shù)學(xué)變換得到闷愤。但是，我們并不總是那么幸運件余，有的時候p(x)的形式十分復(fù)雜讥脐，或者是個高維分布的時候，樣本的生成就十分困難了啼器，此時就需要用到我們下面介紹的MCMC(Markov Chain Monte Carlo)和Gibbs采樣算法了旬渠。而要了解這兩個算法，我們需要對馬氏鏈的平穩(wěn)分布有所認識镀首。##馬氏鏈及其平穩(wěn)分布馬爾科夫鏈是隨機變量$$X_1,X_2,X_3,...$$的一個數(shù)列坟漱，在這個數(shù)列中$$X_{n+ 1}$$的分布僅跟$$X_n$$有關(guān),即$$P(X_{n+1}=x|X_1=x_1,X_2=x_2,...,X_n=x_n) = P(X_{n+1}=x|X_n=x_n)$$ 我們來看一個馬氏鏈的具體例子。社會學(xué)家經(jīng)常把人按照經(jīng)濟狀況分成下更哄、中芋齿、上三層（我們用1、2成翩、3表示）觅捆，他們發(fā)現(xiàn)，決定一個人的收入階層的最重要的因素就是其父母的收入階層：如果一個人的收入屬于下層類別麻敌，那么它的孩子屬于下層收入的概率是0.65,屬于中層收入的概率是0.28栅炒，屬于上層收入的概率是0.07.其轉(zhuǎn)移概率矩陣如下：
如上圖，其狀態(tài)轉(zhuǎn)移矩陣為$$\begin{equation}P=\left[\begin{array}{ccc} 0.65 & 0.28 & 0.07 \ 0.15 & 0.67 & 0.18 \ 0.12 & 0.36 & 0.52 \end{array} \right]\end{equation}$$假設(shè)第一代人處在下术羔、中赢赊、上的概率是$$\pi_0=[\pi_0(1),\pi_0(2),\pi_0(3)]$$,那么他們的子女分布比例將是$$\pi 2 = \pi 1 P = \pi 0 P^2,......,$$第n代子孫的收入分布比例將是$$\pi n=\pi {n-1}P=\pi 0P^n$$。假設(shè)初始概率分布為$$\pi_0 = [0.21,0.68,0.11]$$,則我們計算前n帶人的分布狀況如下
population income distritution transfer
我們發(fā)現(xiàn)從第7代人開始级历，這個分布就穩(wěn)定不變了释移，這是偶然的嗎？如果你還一個初始的概率分布寥殖，就會發(fā)現(xiàn)玩讳，經(jīng)過幾代人的迭代之后涩蜘，分布又收斂了。更為奇特的是熏纯，雖然兩次給定的是不同的初始概率分布同诫，但最終都會收斂到相同的概率分布$$\pi = [0.286,0.489,0.225]$$,也就是說收斂的結(jié)果和初始概率分布$$\pi_0$$無關(guān)。這說明收斂是由概率轉(zhuǎn)移矩陣P決定的樟澜。 $$\begin{equation}P^{20} = P{21}=...=P{100}=...=\left[ \begin{array}{ccc}0.286 & 0.489 & 0.225 \0.286 & 0.489 & 0.225\0.286 & 0.489 & 0.225 \end{array}\right]\end{equation}$$ 我們發(fā)現(xiàn)误窖，當(dāng)n足夠大的時候，這個$$P^n$$矩陣的每一行都是穩(wěn)定的收斂到$$\pi=[0.286,0.489,0.225]$$這個概率分布往扔。事實上贩猎，關(guān)于馬氏鏈的收斂我們有如下定理：馬氏鏈收斂定理: 如果一個非周期的馬氏鏈具有轉(zhuǎn)移矩陣P，且它的任何兩個狀態(tài)是連通的萍膛，那么$$\lim{n \to \infty}P{ij}^n$$存在且與i無關(guān)吭服，記$$\lim{n \to \infty}P{ij}^n = \pi(j)$$,我們有1. $$\begin{equation}\lim{n \to \infty}P^n= \left[ \begin{array}{ccccc}\pi(1)&\pi(2)&...&\pi(j)&...\\pi(1)&\pi(2)&...&\pi(j)&...\...&...&...&...&...\pi(1)&\pi(2)&...&\pi(j)&...\...&...&...&...&...\end{array}\right]\end{equation}$$;2. $$\pi(j)=\sum{i=0}^\infty \pi(i)P_{ij}$$;3. $$\pi$$是方程$$\pi P=\pi$$的唯一非負解其中，$$\pi = [\pi(1),\pi(2),...,\pi(j),...], \quad \sum_{i=0}^\infty = 1$$$$\pi$$稱為馬氏鏈的平穩(wěn)分布蝗罗。馬氏鏈收斂定理非常重要艇棕，所有的MCMC(Markov Chain Monte Carlo)方法都是以這個定理作為理論基礎(chǔ)的。幾個補充說明：1. 雖然我們在示例中用到的馬氏鏈的狀態(tài)空間是有限的串塑，但是沼琉，收斂定理其實并不要求狀態(tài)空間有限，這意味著桩匪，對于連續(xù)型變量打瘪，收斂定理依然有效。2. 兩個狀態(tài)i,j是連通并非指i可以直接一步轉(zhuǎn)移到j(luò)($$P_{ij}>0$$),而是指i可以通過有限的n步轉(zhuǎn)移到達j傻昙。也就是說闺骚，馬氏鏈的任何兩個狀態(tài)連通是指存在一個n，使得矩陣$$P^n$$中的任何一個元素的數(shù)值都大于0.3. 我們不打算證明整個定理妆档，但可以小證一下第二個結(jié)論僻爽，柿子專檢軟的捏`(∩_∩)′ $$\begin{align}P(X_{n+1}=j)&=\sum_{i=0}^\infty P(X_n = i)\cdot P(X_{n+1}=j|X_n = i)\&=\sum_{i=0}^\infty P(X_n =i)P_{ij}\end{align}$$ 上式兩邊取極限即得 $$\pi(j)=\sum_{i=0}^\infty \pi(i)P_{ij}$$由馬氏鏈的收斂定理，概率分布$$\pi_i(x)$$將收斂到平穩(wěn)分布$$\pi(x)$$,假設(shè)到第n步的時候馬氏鏈收斂$$\begin{align}X_0 &\sim \pi_0(x)\X_1 &\sim \pi_1(x)\ &\cdots \X_n &\sim \pi_n(x) = \pi(x) \X_{n+1} &\sim \pi(x)\X_{n+2} &\sim \pi(x)\ & \cdots \end{align}$$所以$$X_n,X_{n+1},X_{n+2},...\sim \pi(x)$$都是同分布(但是并不獨立)的隨機變量贾惦，如果我們從一個具體的初始狀態(tài)$$x_0$$開始胸梆，沿著馬氏鏈按照概率轉(zhuǎn)移矩陣做跳轉(zhuǎn)，那么我們得到一個轉(zhuǎn)移序列$$x_0,x_1,x_2,...,x_n,x_{n+1},...$$,其中须板，$$x_n,x_{n+1},...$$都將是平穩(wěn)分布$$\pi(x)$$的樣本碰镜。##Markov Chain Monte Carlo馬氏鏈能收斂到平穩(wěn)分布，而我們的目的要生成概率分布p(x)的樣本习瑰，一個很直觀的想法就是: 如果我們能夠構(gòu)造一個平穩(wěn)分布$$\pi(x)=p(x)$$的馬氏鏈绪颖，那么我們從任何一個初始狀態(tài)$$x_0$$出發(fā)，得到一個序列$$x_0,x_1,x_2,...,x_n,x_{n+1}...杰刽，$$如果馬氏鏈在第n步已經(jīng)收斂了菠发，于是我們就得到了$$\pi(x)$$的樣本$$x_n,x_{n+1},......$$馬氏鏈的收斂性質(zhì)主要由轉(zhuǎn)移矩陣P決定，所以基于馬氏鏈做采樣的關(guān)鍵問題是如何構(gòu)造轉(zhuǎn)移矩陣P贺嫂，使得平穩(wěn)分布$$\pi (x)=p(x)$$滓鸠。我們使用馬氏鏈的細致平穩(wěn)條件來構(gòu)造轉(zhuǎn)移矩陣P。細致平穩(wěn)條件：如果非周期馬氏鏈的轉(zhuǎn)移矩陣P和分布$$\pi (x)$$滿足 $$\pi(i)P_{ij}=\pi(j)P_{ji}, \quad for \ all \ i,j$$ 則$$\pi(x)$$是馬氏鏈的平穩(wěn)分布第喳，上式被稱為細致平穩(wěn)條件(detailed balance condition)糜俗。這個定理的物理意義是顯而易見的: 對于任何兩個狀態(tài)i,j, 從i轉(zhuǎn)移出去到j(luò)丟失的概率，恰好會被從j轉(zhuǎn)移回i的概率補充回來曲饱，所以狀態(tài)i上的概率$$\pi(i)$$是穩(wěn)定的悠抹。數(shù)學(xué)上的證明也十分簡單: $$\sum_{i=1}{\infty}\pi(i)P_{ij}=\sum_{i=1}{\infty}\pi(j)P_{ji}=\pi(j)\sum_{i=1}{\infty}=\pi(j)$$ 由于$$\pi$$是方程$$\pi P=\pi$$的解，所以$$\pi$$是平穩(wěn)分布扩淀。假設(shè)我們已經(jīng)有一個轉(zhuǎn)移矩陣為Q的馬氏鏈楔敌，$$q(i,j)$$(或$$q(j|i)$$或者$$q(i \to j)$$表示從狀態(tài)i轉(zhuǎn)移到狀態(tài)j的概率)，顯然驻谆，通常情況下 $$p(i)q(i,j) \neq p(j)q(j,i)$$ 也就是細致平穩(wěn)條件不成立卵凑，所以p(x)不太可能是這個馬氏鏈的平穩(wěn)分布。我們可否對馬氏鏈做一個改造胜臊，使得細致平穩(wěn)條件成立呢勺卢？譬如，我們引入一個$$\alpha(i,j)$$,我們希望 $$p(i)q(i,j)\alpha(i,j) = p(j)q(j,i)\alpha(j,i)$$ 取什么樣的$$\alpha(i,j)$$能滿足上述條件呢象对？最簡單的黑忱，按照對稱性，我們可以取 $$\alpha(i,j) = p(j)q(j,i), \ \ \alpha(j,i) = p(i)q(i,j)$$ 于是我們有 $$p(i)\underbrace{q(i,j)\alpha(i,j)}{Q'(i,j)} = p(j)\underbrace{q(j,i)\alpha(j,i)}{Q'(j,i)}$$ 于是我們把原來具有轉(zhuǎn)移矩陣Q的一個很普通的馬氏鏈勒魔，改造為了具有轉(zhuǎn)移矩陣Q'的馬氏鏈甫煞，而Q'滿足細致平穩(wěn)條件，因此馬氏鏈Q(jìng)'的平穩(wěn)分布就是p(x)(即馬氏鏈的初始狀態(tài)分布沥邻！)危虱。在改造Q的過程中引入的$$\alpha(i,j)$$稱為接受率，它的物理意義是從狀態(tài)i以q(i,j)的概率跳轉(zhuǎn)到狀態(tài)j的時候唐全，我們以$$\alpha(i,j)$$的概率接受這個轉(zhuǎn)移埃跷，以$$1-\alpha(i,j)$$的概率留在原狀態(tài)。利用馬爾科夫鏈采樣概率分布p(x)的MCMC算法如下所示：* 初始化馬爾科夫鏈初始狀態(tài)$$X_0=x_0$$;* 對$$t=0,1,2,...,$$循環(huán)以下過程就行采樣1. 第t個時刻馬氏鏈狀態(tài)為$$X_t=x_t$$, 采樣$$y \sim q(x|x_t)$$2. 從均勻分布采樣$$u \sim Uniform[0,1]$$ 3. 如果$$u < \alpha(x_t,y) = p(y)q(x_t|y)$$,則接受轉(zhuǎn)移$$x_t \to y$$,即$$X_{t+1} = y$$;否則不接受轉(zhuǎn)移邮利，$$X_{t+1}=x_t$$##Gibbs Sampling以上的MCMC算法已經(jīng)能夠work了弥雹，但是對于高維的情形，由于接受率$$\alpha$$的存在(通常\alpha <1),以上算法的效率不夠高延届，能夠找到一個轉(zhuǎn)移矩陣Q使得接受率$$\alpha =1$$呢剪勿？我們先看二維的情形，假設(shè)有一個概率分布p(x,y),考察x坐標(biāo)相同的的兩個點$$A(x_1,y_1),B(x_1,y_2)$$,我們發(fā)現(xiàn)$$p(x_1,y_1)p(y_2 | x_1)=p(x_1)p(y_1 |x_1)p(y_2 | x_1) \p(x_1,y_2)p(y_1 | x_1)=p(x_1)p(y_2 |x_1)p(y_1 | x_1)$$也就是說方庭，細致平穩(wěn)條件成立厕吉！事實上酱固，在$$x=x_1$$這條平行于y軸的直線上，任何兩個點之間的轉(zhuǎn)移滿足細致平穩(wěn)條件头朱。同樣的运悲，如果我們在$$y=y_1$$這條直線上任意取兩點$$A(x_1,y_1),C(x_2,y_1)$$,這兩個點之間的轉(zhuǎn)移也滿足細致平穩(wěn)條件。
平面上馬氏鏈轉(zhuǎn)移矩陣的構(gòu)造
于是我們可以按照如下規(guī)則構(gòu)造平面上任意兩點之間的轉(zhuǎn)移概率矩陣Q$$\begin{align} Q(A\rightarrow B) & = p(y_B|x_1) & \text{如果} \quad x_A=x_B=x_1 & \ Q(A\rightarrow C) & = p(x_C|y_1) & \text{如果} \quad y_A=y_C=y_1 & \ Q(A\rightarrow D) & = 0 & \text{其它} & \end{align}$$有了如上的轉(zhuǎn)移矩陣Q项钮，我們?nèi)菀昨炞C平面上任意兩點$$X,Y$$,都滿足細致平穩(wěn)條件$$p(X)Q(X \to Y)=p(Y)q(Y \to X)$$于是這個二維空間上的馬氏鏈將收斂到平穩(wěn)分布$$p(x,y)$$班眯。這個算法就稱為Gibbs算法，是由物理學(xué)家Gibbs首先提出來的烁巫。二維Gibbs Sampling算法1. 隨機初始化$$X_0=x_0,Y_0=y_0$$;2. 對$$t=0,1,2,...$$循環(huán)采樣* $$y_{t+1} \sim p(y|x_t)$$;* $$x_{t+1} \sim p(x|y_{t+1})$$署隘。如上圖所示，馬氏鏈的轉(zhuǎn)移只是輪換的沿著坐標(biāo)軸x軸和y軸做轉(zhuǎn)移亚隙，于是得到樣本$$(x_0,y_0),(x_0,y_1),(x_1,y_2),(x_2,y_2),...$$馬氏鏈收斂后磁餐，最終得到的樣本就是$$p(x,y)$$的樣本。我們上邊的算法雖然是坐標(biāo)軸輪換采樣的恃鞋，但是這其實是不強制要求的;一般的情形是崖媚，在t時刻，在x軸和y軸之間隨機的選一個坐標(biāo)軸恤浪，然后按條件概率做轉(zhuǎn)移畅哑，馬氏鏈也一樣收斂。由以上過程我們可以看出水由，如果$$x_1$$變?yōu)槎酁榍樾?$\vec x_1$$荠呐，細致平穩(wěn)條件同樣成立$$p(x_1,y_1)p(y_2 |x_1)=p(x_1,y_2)p(y_1|x_1)$$此時轉(zhuǎn)移矩陣Q由條件分布$$p(y|x_1)$$定義，所以n為空間中對于概率分布$$p(x_1,x_2,...,x_n)$$可以如下定義轉(zhuǎn)移矩陣1. 如果當(dāng)前狀態(tài)為$$(x_1,x_2,...,x_n)$$,馬氏鏈轉(zhuǎn)移的過程中砂客，只能沿著坐標(biāo)軸轉(zhuǎn)移泥张。沿著$$x_i$$這跟坐標(biāo)軸做轉(zhuǎn)移的時候，轉(zhuǎn)移概率由條件概率$$p(x_i|x_1,...,x_{i-1},x_{i+1},...,x_n)$$定義鞠值；2. 其它無法沿著單根坐標(biāo)軸進行的跳轉(zhuǎn)媚创，轉(zhuǎn)移概率都設(shè)置為0.于是我們可以把Gibbs Sampling算法從二維推廣到n維，當(dāng)以上算法收斂后彤恶，得到的就是概率分布$$p(x_1,x_2,...,x_n)$$的樣本钞钙，當(dāng)然這些樣本并不獨立，但是我們要求的是采樣得到的樣本符合給定的概率分布声离，并不要求獨立芒炼。##后記##1. MCMC算法和Gibbs Sampling都能得到指定的概率分布p(x)的分布樣本，但是樣本之間并不獨立;2. 由于MCMC算法存在接受率$$\alpha$$,使得在高維的時候效率低下术徊，Gibbs Sampling采用坐標(biāo)軸輪換的采樣方法,巧妙的避開了接受率$$\alpha$$(其實相當(dāng)于接受率為1)本刽。
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
