當(dāng)前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

Probability, Matringale,Markov Chain, MCMC

發(fā)布時間：2023/12/2 编程问答 45 豆豆

生活随笔收集整理的這篇文章主要介紹了 Probability, Matringale,Markov Chain, MCMC 小編覺得挺不錯的,現(xiàn)在分享給大家,幫大家做個參考.

一、基本知識

1. 條件概率

條件概率是指在某件事情已經(jīng)發(fā)生的前提下，另一件事情在此基礎(chǔ)上發(fā)生的概率，舉例來說P(A丨B)表示B發(fā)生的基礎(chǔ)上，A也發(fā)生的概率，基本公式為：

2. 條件期望

在上述概率下的期望我們稱之為條件期望。
計算：

3. 隨機(jī)過程

在概率論概念中，隨機(jī)過程是隨機(jī)變量的集合。若一隨機(jī)系統(tǒng)的樣本點(diǎn)是隨機(jī)函數(shù)，則稱此函數(shù)為樣本函數(shù)，這一隨機(jī)系統(tǒng)全部樣本函數(shù)的集合是一個隨機(jī)過程。實(shí)際應(yīng)用中，樣本函數(shù)的一般定義在時間域或者空間域。

顧名思義，它其實(shí)就是個過程，比如今天下雨，那么明天下不下雨呢？后天下不下雨呢？從今天下雨到明天不下雨再到后天下雨，這就是個過程。那么怎么預(yù)測N天后到底下不下雨呢？這其實(shí)是可以利用公式進(jìn)行計算的，隨機(jī)過程就是這樣一個工具，把整個過程進(jìn)行量化處理，用公式就可以推導(dǎo)出來N天后的天氣狀況，下雨的概率是多少，不下雨的概率是多少。說白了，隨機(jī)過程就是一些統(tǒng)計模型，利用這些統(tǒng)計模型可以對自然界的一些事物進(jìn)行預(yù)測和處理，比如天氣預(yù)報，比如股票，比如市場分析，比如人工智能。

隨機(jī)過程是一族時間函數(shù)的集合，隨機(jī)過程的每個樣本函數(shù)（一個隨機(jī)信號）是一個確定的時間函數(shù)x(t)（t=0~T），隨機(jī)過程在一個確定的時刻t1是一個隨機(jī)變量X(t1)。
從上述對隨機(jī)過程描述中可以看出，從信號分析的角度來看，一個隨機(jī)過程就是針對一個實(shí)驗(yàn)測量無數(shù)個的樣本，每個樣本就是一個隨機(jī)信號（要分清什么是隨機(jī)過程和隨機(jī)信號），而所有樣本測量的時間都是0~T（這看起來很麻煩和不可思議），然后，我們把這些樣本（每個隨機(jī)信號）放在一起擺放，見圖

其實(shí)，這就是一個隨機(jī)過程，然后，我們怎樣描述呢？我們從兩個角度來描述它（和上述定義對應(yīng)）。
（1）從圖中縱向的角度（0~T），一個隨機(jī)過程是你這次針對于某個任務(wù)進(jìn)行測量的所有樣本信號的集合。
（2）從橫向的角度（在0-T區(qū)間的任意時刻），你可以拿把“菜刀”在任意時刻橫著切一下，這時，你會看到所有你測量的樣本信號在這個時刻的幅度值都不一樣，它們是隨機(jī)的，為了表征這些幅度值的統(tǒng)計特性，我們用此時刻的隨機(jī)變量X(t)來表征，那么在0~T時間區(qū)間內(nèi)有多少隨機(jī)變量呢？當(dāng)然是無數(shù)個（相當(dāng)于你切了無數(shù)刀），這些無數(shù)個隨機(jī)變量的集合就是隨機(jī)過程。
如果所有工作都像這樣去做的話，那么，所有隨機(jī)信號分析就沒有辦法完成了，所以，人們將問題進(jìn)行簡化，盡量將隨機(jī)過程說成平穩(wěn)的并且具有各態(tài)歷經(jīng)的。

二、鞅（martingale）

在概率論中，鞅（martingale）是滿足下述條件的隨機(jī)過程：已知過去某一時刻 s 以及之前所有時刻的觀測值，若某一時刻 t 的觀測值的條件期望等于過去某一時刻 s 的觀測值，則稱這一隨機(jī)過程是鞅。而于博弈論中，鞅經(jīng)常用來作為公平博弈的數(shù)學(xué)模型。

鞅的原名 martingale 原指一類于18世紀(jì)流行于法國的投注策略，稱為加倍賭注法[1]。這類策略中最簡單的一種策略是為博弈設(shè)計的。在博弈中，賭徒會擲硬幣，若硬幣正面向上，賭徒會贏得賭本，若硬幣反面向上，賭徒會輸?shù)糍€本。這一策略使賭徒在輸錢后加倍賭金投注，為的是在初次贏錢時贏回之前輸?shù)舻乃绣X，同時又能另外贏得與最初賭本等值的收益。當(dāng)賭徒的財產(chǎn)和可用時間同時接近無窮時，他擲硬幣后硬幣正面向上的概率會接近1，由此看來，加倍賭注法似乎是一種必然能贏錢的策略。然而，賭金的指數(shù)增長最終會導(dǎo)致使用這一策略的賭徒破產(chǎn)。

鞅的概念首先是由保羅·皮埃爾·萊維于 1934 年提出的，但他只提出了離散時間的版本，而且沒有給予命名。直到 1939 年，約翰?維爾將此概念推廣到連續(xù)時間的情況，并且首次提出 martingale 這個名稱。約瑟夫·利奧·杜布（Joseph Leo Doob）等人在鞅的相關(guān)理論的初期發(fā)展做出重大貢獻(xiàn)，而完成這些工作的部分動機(jī)是為了表明成功的投注策略不可能存在。此外，伊藤清在分析應(yīng)用方面作出了重要的貢獻(xiàn)。從 1970 年代開始，鞅論就在純粹數(shù)學(xué)和應(yīng)用數(shù)學(xué)的很多領(lǐng)域中有廣泛的應(yīng)用，特別是在數(shù)學(xué)物理和金融數(shù)學(xué)中。

以下我們討論T個離散階段下，鞅的相關(guān)定義和性質(zhì)。

我們將一個隨機(jī)試驗(yàn)的所有基本結(jié)果稱為樣本點(diǎn)，記為：
稱其構(gòu)成的集合文件樣本空間：
其對應(yīng)的概率空間記為：
記每個階段的狀態(tài)為S0，S1, …, ST，我們定義以下集合：
Ft可以理解為t時刻前所有狀態(tài)的信息，即整個隨機(jī)過程到達(dá)當(dāng)前狀態(tài)的所有可能路徑。我們再定義:

這樣一來F就包含了所有時刻，所有狀態(tài)，所有路徑的信息。以下的性質(zhì)顯而易見：

因?yàn)楝F(xiàn)t2時刻的信息肯定是要多于t1時刻的。
我們繼續(xù)有以下定義：

接下來我們給出條件期望的基本性質(zhì)：

2.1 Martingale的定義

我們開始討論鞅的定義，鞅有以下三種等價的定義：

我們以下證明三個定義是等價的：
以上便是鞅的基本定義，為了更好的理解鞅，我們可以舉個例子。
假設(shè)賭徒在t時刻有 $M_t$ 的財富，如果賭局是公平的，那么他未來財富的期望應(yīng)當(dāng)是與當(dāng)前財富相同，同時不受之前賭博策略的影響（因?yàn)槲覀兦蟮氖菞l件期望）。事實(shí)上martingale理論最早便是用于賭博策略的研究，我們現(xiàn)在討論的鞅是公平賭局，也存在像上鞅和下鞅這樣不公平賭局。至于為什么中文會翻譯成鞅呢，因?yàn)樵诜ㄎ睦?#xff0c;martingale有兩層意思，一是“倍賭策略”（每輸一局就賭資加倍，只要贏一局就扳回所有損失），二是“馬韁繩”，將其翻譯成鞅（馬的韁繩）是用了它的第二層含義。

2.2 Martingale, supermartingale, submartingale

鞅： $E[X_t | F_s] = X_s （t>s）$ ；
上鞅： $E[X_t | F_s] < X_s （t>s）$ ；
下鞅： $E[X_t | F_s] > X_s （t>s）$ 。

通俗理解：Fs表示你在s時刻所掌握的所有信息，那么E[Xt | Fs]就表示你在s時刻所掌握的所有信息對X在未來t時刻預(yù)期的判斷，如何和X在s時刻是一致的說明游戲是公平的；上鞅說明對過程X是看衰的，下鞅則是看好。
鞅是公平賭博，下鞅是有利賭博，上鞅是不利賭博。

鞅直觀意義是一條（有固定趨勢的）線。
考慮一個函數(shù)f(t)，橫軸是時間軸，縱軸是函數(shù)值。這個函數(shù)可以是平的（就是 $f (t) = c$ ），也可以單調(diào)增，也可以單調(diào)減。因此從圖像上看是一條有固定趨勢的線。
如果僅僅是線，那么未免太乏味無趣。概率就是把固定取值的變量 $X = c$ 取代為可以同時取多個值、每個值都有一定權(quán)重的隨機(jī)變量 $X$ ，相當(dāng)于引入了多個可能狀態(tài)。隨機(jī)過程進(jìn)一步引入一個時間軸，考慮隨機(jī)變量的序列 $X_n$ （或者連續(xù)版本的 $X_t$ ），因此如果畫在剛剛那個圖里， $X_n(w)$ 可以看做某個大體上是線、局部有些抖動的點(diǎn)列，而 $X_n$ 就是這些點(diǎn)列的集合。
那么離散鞅是一種點(diǎn)列集（隨機(jī)過程），這個點(diǎn)列集大體趨勢是 $f (t) = c$ ，只是會有一些抖動，而且是考慮一個抖動圍繞這條線的點(diǎn)列集整體。類似的，上鞅、下鞅就是單調(diào)減、單調(diào)增函數(shù)的點(diǎn)列集。

2.3 Stopping time

關(guān)于隨機(jī)過程 X1，X2，X3，… 的停時是隨機(jī)變量 τ，這一隨機(jī)變量具有如下性質(zhì)：對于每一個時間，事件 τ = t 的發(fā)生與否僅取決于 X1，X2，X3，…，Xt 的取值。從定義中可以感受到的直覺是在任一特定時刻 t，我們都可以知道在這一時刻隨機(jī)過程是否到了停時。現(xiàn)實(shí)生活中停時的例子如賭徒離開賭桌的時刻，這一時刻可能是賭徒以前贏得錢財?shù)暮瘮?shù)（例如，僅當(dāng)他沒有錢時，他才可能離開賭桌），但是他不可能根據(jù)還未完成的博弈的結(jié)果來選擇離開還是留下。

上述停時定義滿足強(qiáng)條件，下面給出一個弱條件的停時定義：若事件 τ = t 的發(fā)生與否統(tǒng)計獨(dú)立于 Xt+1，Xt+2，… 但并不是完全決定于時刻 t 以及之前的過程歷史，則隨機(jī)變量 τ 是停時。雖然這是一個弱條件，但在需要用到停時的證明中的一些情況也算是足夠強(qiáng)的條件。

鞅的一個基本性質(zhì)是若 $X_{t})_{{t>0}}$ 是下\上鞅且 $τ\tau$ 是停時，由 $Xtτ:=Xmin?{τ,t}X_{t}^{\tau }:=X_{{\min\{\tau ,t\}}}$ 定義的對應(yīng)停止過程 $(Xtτ)t>0(X_{t}^{\tau })_{{t>0}}$ 也是下\上鞅。

停時鞅的概念引出了一系列定理，例如可選停止定理（又稱可選抽樣定理）：在特定條件下，停時的鞅的期望值等于其初始值。利用這一定理，我們可以證明對于一個壽命有限且房產(chǎn)有限的賭徒，成功的投注策略不可能存在。

2.3.1 Optional stopping theorem

In probability theory, the optional stopping theorem (or Doob’s optional sampling theorem) says that, under certain conditions, the expected value of a martingale at a stopping time is equal to its initial expected value. Since martingales can be used to model the wealth of a gambler participating in a fair game, the optional stopping theorem says that, on average, nothing can be gained by stopping play based on the information obtainable so far (i.e., without looking into the future). Certain conditions are necessary for this result to hold true. In particular, the theorem applies to doubling strategies.

The optional stopping theorem is an important tool of mathematical finance in the context of the fundamental theorem of asset pricing.

2.3.1.1 discrete-time version of the theorem

2.3.1.2 Applications

2.3.1.3 Proof

三、markov chain

3.1 Markov Property

在學(xué)習(xí)馬爾可夫鏈之前，我們首先了解一下什么是馬爾可夫性質(zhì)（或者叫馬爾可夫特性，Markov Property）。
假設(shè)你有一個系統(tǒng)，這個系統(tǒng)具有 $M$ 種可能的狀態(tài)，并且在狀態(tài)之間正在移動。真實(shí)世界中，我們常常見到這種例子，例如天氣從炎熱到溫和，或者是股票市場從熊市到牛市再到停滯的狀態(tài)（stagnant states）。

3.1.1 馬爾可夫性的定義

關(guān)于馬爾可夫性，我們給出了如下的Definition：

從上述的式子可以看出，t+1時刻的狀態(tài)包含了1,…,t時刻狀態(tài)的全部歷史信息，并且當(dāng)我們知道t時刻的狀態(tài)后，我們只關(guān)注于環(huán)境的信息，而不用管之前所有狀態(tài)的信息，這就是馬爾可夫性，當(dāng)論文中說某一狀態(tài)或其他信息符合馬爾可夫性時，我們也應(yīng)當(dāng)聯(lián)想到這個性質(zhì)。

3.1.2 State Transition Matrix狀態(tài)傳輸矩陣

如果一個過程具有馬爾可夫性質(zhì)，它也被成為馬爾可夫過程（Markov Process）。

3.2 Markov Processes馬爾可夫過程

在概率論及統(tǒng)計學(xué)中，馬爾可夫過程（英語：Markov process）是一個具備了馬爾可夫性質(zhì)的隨機(jī)過程，因?yàn)槎韲鴶?shù)學(xué)家安德雷·馬爾可夫得名。馬爾可夫過程是不具備記憶特質(zhì)的（memorylessness）。換言之，馬爾可夫過程的條件概率僅僅與系統(tǒng)的當(dāng)前狀態(tài)相關(guān)，而與它的過去歷史或未來狀態(tài)，都是獨(dú)立、不相關(guān)的。

3.3 markov chain馬爾可夫鏈

具備離散狀態(tài)的馬爾可夫過程，通常被稱為馬爾可夫鏈。

3.3.1 定義與解釋

3.3.1.1 定義舉例1

3.3.1.2 定義舉例2

Markov Chain 體現(xiàn)的是狀態(tài)空間的轉(zhuǎn)換關(guān)系，下一個狀態(tài)只決定于當(dāng)前的狀態(tài).通俗點(diǎn)說就是：你下一時刻在哪只跟你現(xiàn)在在哪有關(guān)，于你以前去過哪沒有關(guān)系. 如下圖：

這個狀態(tài)圖的轉(zhuǎn)換關(guān)系可以用一個轉(zhuǎn)換矩陣 T 來表示：

舉一個例子，如果當(dāng)前狀態(tài)為 u(x) = (0.5, 0.2, 0.3), 那么下一個矩陣的狀態(tài)就是 u(x)T = (0.18, 0.64, 0.18), 依照這個轉(zhuǎn)換矩陣一直轉(zhuǎn)換下去，最后的系統(tǒng)就趨近于一個穩(wěn)定狀態(tài) (0.22, 0.41, 0.37) 。而事實(shí)證明無論你從那個點(diǎn)出發(fā)，經(jīng)過很長的 Markov Chain 之后都會匯集到這一點(diǎn)。

滿足什么條件下經(jīng)過很長的 Markov Chain 后系統(tǒng)會趨近一個穩(wěn)定狀態(tài)呢，大概的條件如下：

Irreducibility. 即圖是聯(lián)通的，各個狀態(tài)之間都有過去的辦法，舉個不聯(lián)通的例子，比如爬蟲爬不到內(nèi)部局域網(wǎng)的網(wǎng)頁
Aperiodicity. 非周期性, 即圖中遍歷不會陷入到一個死圈里，進(jìn)去了再也出不來，有些網(wǎng)站為了防機(jī)器人，會專門設(shè)置這種陷阱
Detailed Balance，這是保證系統(tǒng)有穩(wěn)態(tài)的一個重要條件，詳細(xì)說明見下面。
假設(shè) p(x) 是最后的穩(wěn)態(tài)，那么 detailed balance 可以用公式表示為：

什么意思呢？假設(shè)上面狀態(tài)圖 x1 有 0.22 元， x2 有 0.41 元，x3 有 0.37 元，那么 0.22×1 表示 x1 需要給 x2 錢，以此類推，手動計算，可以發(fā)現(xiàn)下一個狀態(tài)每個人手中的錢都沒有變。值得說明的是，這里體現(xiàn)了一個很重要的特性，那就是從一個高概率狀態(tài) xi 向一個低概率狀態(tài) x(i-1) 轉(zhuǎn)移的概率等于從這個低概率狀態(tài)向高概率狀態(tài)轉(zhuǎn)移的概率（reversible，至于要不要轉(zhuǎn)移又是另外一回事)
ergodicity, 遍歷性，是指不管事物現(xiàn)在出于什么狀態(tài)，在較長時間內(nèi)，馬爾科夫過程逐漸趨于穩(wěn)定狀況，而且穩(wěn)定狀態(tài)與初始狀況無關(guān)。

3.3.1.3 通俗解釋

先說說我們村智商為0的王二狗，人傻不拉幾的，見人就傻笑，每天中午12點(diǎn)的標(biāo)配，仨狀態(tài)：吃，玩，睡。這就是傳說中的狀態(tài)分布。

你想知道他n天后中午12點(diǎn)的狀態(tài)么？是在吃，還是在玩，還是在睡？這些狀態(tài)發(fā)生的概率分別都是多少？
先看個假設(shè)，他每個狀態(tài)的轉(zhuǎn)移都是有概率的，比如今天玩，明天睡的概率是幾，今天玩，明天也玩的概率是幾幾，看圖更清楚一點(diǎn)。

這個矩陣就是轉(zhuǎn)移概率矩陣P，并且它是保持不變的，就是說第一天到第二天的轉(zhuǎn)移概率矩陣跟第二天到第三天的轉(zhuǎn)移概率矩陣是一樣的。（這個叫時齊，有興趣的同學(xué)自行百度）。
有了這個矩陣，再加上已知的第一天的狀態(tài)分布，就可以計算出第N天的狀態(tài)分布了。

S1 是4月1號中午12點(diǎn)的的狀態(tài)分布矩陣 [0.6, 0.2, 0.2]，里面的數(shù)字分別代表吃的概率，玩的概率，睡的概率。
那么
4月2號的狀態(tài)分布矩陣 S2 = S1 * P (倆矩陣相乘)。
4月3號的狀態(tài)分布矩陣 S3 = S2 * P (看見沒，跟S1無關(guān)，只跟S2有關(guān))。
4月4號的狀態(tài)分布矩陣 S4 = S3 * P (看見沒，跟S1，S2無關(guān)，只跟S3有關(guān))。
…
4月n號的狀態(tài)分布矩陣 Sn = Sn-1 * P (看見沒，只跟它前面一個狀態(tài)Sn-1有關(guān))。

總結(jié)：馬爾可夫鏈就是這樣一個任性的過程，它將來的狀態(tài)分布只取決于現(xiàn)在，跟過去無關(guān)！

就把下面這幅圖想象成是一個馬爾可夫鏈吧。實(shí)際上就是一個隨機(jī)變量隨時間按照Markov性進(jìn)行變化的過程。

附：S2 的計算過程 (沒興趣的同學(xué)自行略過)

3.3.2 n步轉(zhuǎn)移概率

除了一步轉(zhuǎn)移以外，馬爾科夫鏈還有n步轉(zhuǎn)移，即通過n次達(dá)到目標(biāo)狀態(tài)

其實(shí)是Markov Chain 的基本性質(zhì)“無后效性”，即事物將來的狀態(tài)及其出現(xiàn)的概率的大小，只取決于該事物現(xiàn)在所處的狀態(tài)，而與以前時間的狀態(tài)無關(guān)

3.3.2.1 C-K方程

查普曼-柯爾莫格洛夫方程(Chapman-Kolmogorov equation,C-K equation)給出了計算 [公式] 步轉(zhuǎn)移概率的一個方法：

3.3.3 狀態(tài)分類

①可達(dá)的：i到j(luò)的概率為正

②常返的：從i出發(fā)后還能再回到i

③非常返的：從i出發(fā)后不一定能回到j(luò)

常返和非常返的區(qū)別！：常返指出發(fā)后可以無限次回到自身；非常返并不是無法回到自身，而是只能有限次回到自身。

常返類：相互可達(dá)的狀態(tài)組。

馬爾科夫鏈的分解：

一個馬爾科夫鏈可以分解成一個或者多個常返類，加上可能的非常返類。
一個常返態(tài)從它所屬的類里任何一個狀態(tài)出發(fā)都是可達(dá)的，但從其他類里的常返狀態(tài)出法是不可達(dá)的。
從任何一個常返狀態(tài)出發(fā)都不可到達(dá)非常返狀態(tài)
從一個非常返狀態(tài)出發(fā)，至少有一個常返態(tài)是可達(dá)的

常返類的周期：

若一個常返類能被分成幾個組，這幾個組之間的下一個狀態(tài)轉(zhuǎn)移必定不在自己組。則這個常返類被稱為有周期的。

命題得證。

很顯然，暫態(tài)也是一個類性質(zhì)。而利用上述性質(zhì)可以得到：有限馬爾可夫鏈的所有狀態(tài)不可能都是暫態(tài)，有限不可約馬爾可夫鏈的所有狀態(tài)都是常返態(tài)。

3.3.4 穩(wěn)態(tài)

3.3.4.1 穩(wěn)態(tài)概率舉例1

3.3.4.2 穩(wěn)態(tài)概率舉例2

舉個具體的例子。社會學(xué)家把人按其經(jīng)濟(jì)狀況分為3類：下層，中層，上層，我們用1,2,3表示這三個階層。社會學(xué)家發(fā)現(xiàn)決定一個人的收入階層最重要的因素就是其父母的收入階層。如果一個人的收入屬于下層類別，則它的孩子屬于下層收入的概率為0.65，屬于中層收入的概率為0.28，屬于上層收入的概率為0.07。從父代到子代，收入階層轉(zhuǎn)移概率如下

我們用P表示這個轉(zhuǎn)移矩陣，則

假設(shè)第1代人的階層比例為
則前10代人的階層分布如下

我們可以看到，在相同的轉(zhuǎn)移矩陣作用下，狀態(tài)變化最終會趨于平穩(wěn)。對于第n代人的階層分布，我們有
從表達(dá)式上我們可以看到，π是一維向量，P是兩維矩陣，P進(jìn)行足夠多次自乘后，值趨于穩(wěn)定。

3.3.5 馬氏鏈平穩(wěn)分布

在轉(zhuǎn)移矩陣P作用下達(dá)到的平穩(wěn)狀態(tài)，我們稱之為馬氏鏈平穩(wěn)分布。對于這個特性，有如下精彩定理

我在這里直觀的解釋一下上面定理

條件

（1）非周期馬氏鏈：馬氏鏈轉(zhuǎn)移要收斂，就一定不能是周期性的。不做特別處理，我們處理的問題基本上都是非周期性的，在此不做多余解釋。
（2）存在概率轉(zhuǎn)移矩陣P，任意兩個狀態(tài)是連通的：這里的連通可以不是直接相連，只要能夠通過有限次轉(zhuǎn)移到達(dá)即可。比如對于a, b, c狀態(tài)，存在a->b, b->c，則我們認(rèn)為a到c是可達(dá)的。

結(jié)論

（1）不論初始狀態(tài)是什么，經(jīng)過足夠多次概率轉(zhuǎn)移后，會存在一個穩(wěn)定的狀態(tài)π。
（2）概率轉(zhuǎn)移矩陣自乘足夠多次后，每行值相等。即

3.3.5.1 馬爾科夫鏈平穩(wěn)狀態(tài)定理的物理解釋

我們再用一個更加簡單的例子來闡明這個定理的物理含義。假設(shè)城市化進(jìn)程中，農(nóng)村人轉(zhuǎn)移為城市人的概率為0.5，城市人轉(zhuǎn)移為農(nóng)村人的概率為0.1。

假設(shè)一開始有100個農(nóng)村人，0個城市人，每代轉(zhuǎn)移人數(shù)如下

可以看到，城市化進(jìn)程中馬爾科夫平穩(wěn)狀態(tài)就是農(nóng)村人轉(zhuǎn)移為城市人的速度等于城市人轉(zhuǎn)移為農(nóng)村人的速度。對于上述轉(zhuǎn)移矩陣P，平穩(wěn)分布為農(nóng)村人17%，城市人83%。如果我們可以得到當(dāng)前中國城市化轉(zhuǎn)移矩陣P，我們就可以算出中國最終城市化率大概為多少（這里不考慮P的變化）。同時如果我們知道了中國城市化人口比例，我們就能知道城市化進(jìn)程還可以持續(xù)多少代人。

3.3.5.2 長期頻率解釋

3.3.6 Birth-Death process 生滅過程

3.3.6.1 平穩(wěn)分布

3.3.6.2 性質(zhì)

3.3.6.3 例子

3.3.6.4 Probability balance equations 概率平衡方程

3.3.7 吸收概率和吸收期望時間

3.3.7.1 平均時間

平均吸收時間：i的平均吸收時間為從狀態(tài)i開始直到到達(dá)吸收態(tài)所需要的期望步數(shù)

平均吸收時間計算：解平均吸收時間方程組

3.3.8 為什么馬爾可夫鏈這么重要

因?yàn)樗撵o態(tài)分布（或叫平穩(wěn)性分布，stationary distribution）性質(zhì)。
我們用一下Wiki上給出的股票市場的例子（見下圖），這是一個具有靜態(tài)空間的馬爾可夫鏈（a continuous-time Markov chain with state-space），包含了牛市、熊市、停滯市場三種狀態(tài)（Bull market, Bear market, Stagnant market）。

它各狀態(tài)之間由轉(zhuǎn)移概率組成的轉(zhuǎn)移矩陣（ transition rate matrix）如下，記為 Q：

你可以嘗試從其他狀態(tài)出發(fā)，但最終的靜態(tài)分布一定是相同的。那么得出這個靜態(tài)分布有什么用處呢？

靜止?fàn)顟B(tài)分布的重要性在于它可以讓你在隨機(jī)的時間里，對于一個系統(tǒng)的任意一個狀態(tài)確定其概率。（The stationary state distribution is important because it lets you define the probability for every state of a system at a random time.）

對于上述例子，你可以說整個系統(tǒng)，62.5%的概率系統(tǒng)將處于牛市，31.25%的概率將處于熊市，而6.25%的概率將處于停滯態(tài)。

直觀上說，你可以想想成在一個狀態(tài)鏈上進(jìn)行隨機(jī)游走。你在一個狀態(tài)出，然后你決定你前往下一個狀態(tài)的方法是在給定當(dāng)前狀態(tài)下觀察下一個狀態(tài)的概率分布。

3.3.9 小結(jié)

馬爾科夫鏈?zhǔn)且粋€數(shù)學(xué)對象，包含一系列狀態(tài)以及狀態(tài)之間的轉(zhuǎn)移概率，如果每個狀態(tài)轉(zhuǎn)移到其他狀態(tài)的概率只與當(dāng)前狀態(tài)相關(guān)，那么這個狀態(tài)鏈就稱為馬爾科夫鏈。有這樣一個馬爾科夫鏈之后，我們可以任取一個初始點(diǎn)，然后根據(jù)狀態(tài)轉(zhuǎn)移概率進(jìn)行隨機(jī)游走。假設(shè)能夠找到這樣一個馬爾科夫鏈，其狀態(tài)轉(zhuǎn)移概率正比于我們想要采樣的分布（如貝葉斯分析中的后驗(yàn)分布），采樣過程就變成了簡單地在該狀態(tài)鏈上移動的過程。

3.4 連續(xù)時間馬爾科夫鏈CTMC

3.4.1 定義

3.4.1.1 定義1:一般定義

和起始時間t無關(guān)的話，我們稱這是時間齊次的馬爾科夫鏈。這個轉(zhuǎn)移矩陣和離散時間不同的是，離散時間給出的是一步轉(zhuǎn)移概率，但是連續(xù)馬爾科夫鏈的轉(zhuǎn)移概率給出的是和時間相關(guān)的。

3.4.1.2 定義2

我們考慮連續(xù)時間馬爾科夫鏈從一個狀態(tài) i開始，到狀態(tài)發(fā)生變化，比如變成j所經(jīng)過的時間，由于馬爾科夫鏈的馬爾科夫性，這個時間是具有無記憶性的，所以這個時間是服從指數(shù)分布的。這和離散時間馬爾科夫鏈?zhǔn)敲芮邢嚓P(guān)的，離散時間馬爾科夫鏈中的時間是離散時間，因?yàn)橛蔁o記憶性，所以是服從幾何分布的。

這樣我們就可以這樣定義連續(xù)時間馬爾科夫鏈。馬爾科夫鏈?zhǔn)沁@樣的一個過程。

(i)在轉(zhuǎn)移到不同的狀態(tài) $i$ 前，它處于這個狀態(tài)的時間是速率為 $v_i$ 的指數(shù)分布。
(ii)當(dāng)離開狀態(tài) $i$ 時，以某種概率 $P_{ij}$ 進(jìn)入下一個狀態(tài) $j$ ,當(dāng)然 $P_{ij}$ 滿足

對比于半馬爾科夫鏈，我們可以發(fā)現(xiàn)，連續(xù)時間馬爾科夫鏈?zhǔn)且环N特殊的半馬爾科夫鏈，在一個狀態(tài)所待的時間是只不過是一個具體的分布–指數(shù)分布，而半馬爾科夫鏈只是說所待的時間是任意的一個隨機(jī)時間。

3.4.1.3 定義3

3.4.1.4 舉例

3.4.2 兩個定義(1和2)之間的關(guān)系

3.4.3 CTMC的極限概率

3.4.4 CTMC極限分布的應(yīng)用

3.4.5 時間可逆性

3.4.6 截止的馬爾科夫鏈

四、martingale與markov chain的關(guān)系

二者都是隨機(jī)過程中的一種過程。

Martingale的詞本意是指馬車上馬夫控制馬前進(jìn)的韁繩（如果我記得沒錯的話），所以從詞源來看刻畫了一種過程前進(jìn)（未來）與現(xiàn)在出發(fā)點(diǎn)關(guān)系的含義。具體來說韁繩的作用是使得馬車的前進(jìn)方向與現(xiàn)在所沖的方向一致，所以在概率上來解釋就是未來發(fā)生的所有路徑可能性的期望值與現(xiàn)在的出發(fā)點(diǎn)一致。

從這個意義上來說Matingale的核心是說明了一個過程在向前演化的過程中的穩(wěn)定性性質(zhì)。但它并沒有說明這個過程是如何到達(dá)這一時間點(diǎn)的（是否由上一個時間點(diǎn)所在的位置決定，matingale并沒有說明）。再用馬車的例子來說，Martingale告訴了你馬車在未來是怎么向前走的，中間會有左右的波動（比如馬、車夫走神了，路上有坑馬要繞開，etc.），但整體來說馬是沿著一條直線向前走的。

而馬爾科夫過程的核心在于點(diǎn)明了過程的演化是無記憶性的。還拿馬車來舉栗子，假設(shè)車夫喝醉了，他沒有意識并在一個很大的廣場，馬車下一刻前進(jìn)的方向并不需要是一條直線（經(jīng)過車夫與馬的直線，這種情況下韁繩是繃直的，是martingale），或者說韁繩由于車夫沒繃緊是松垮的，這種情況下馬車在下一刻可以去任何一個方向，整體上來說前進(jìn)方向也不必須有什么穩(wěn)定性規(guī)律可循，但整個過程唯一的共性是馬邁出前腿的時候，能夠到達(dá)的所有可能范圍，是由它的后腿（你現(xiàn)在所在的位置）決定的（但馬可能扭一扭屁股，身子彎曲一下，所以不必須走直線，不必需走直線，不必需走直線！），而并沒有由上一步馬所在的位置決定，這也就是所謂的無記憶性。

所以從這兩個角度來理解，兩個名詞

有共性：

都從一個過程的全生命角度描畫了一個過程的演進(jìn)性質(zhì)，

有重疊：
當(dāng)還是馬車?yán)拥臅r候,martingale也是一個markov（因?yàn)殡m然走直線，但下一刻的位置還是只由現(xiàn)在決定，只是馬身子不能扭曲，不能改變方向），這個例子在概率上最熟悉的模型就是brownian motion了；而反過來，馬車未來位置由現(xiàn)在決定，又走直線，所以此時markov process 也是一個martingale （例子還是brownian motion）；

但更重要的是兩個過程本質(zhì)上不是在講一回事：比如還是馬車車夫，喝醉了但走在一個三維空間，這時候它是一個markov process,但是由于方向不確定，此時已經(jīng)不是martingale而變成了一個local martingale；而反過來，假設(shè)有一個錯幀宇宙，空間共享但時間差一天，這時候同一個馬車走在不同的宇宙里（但行走軌跡獨(dú)立），韁繩拉直，此時兩架馬車都走之前，兩架馬車組成的系統(tǒng)是一個martingale，但是由于下一時刻前進(jìn)的方向與宇宙1中的此時有關(guān)，也與宇宙2中的昨天有關(guān)，所以兩架馬車組成的系統(tǒng)就不再是一個markov了。

總結(jié)一下，brownian motion (wiener process)既是markov process 又是 martingale; 而markov process 與martingale是相交而非包含與反包含關(guān)系。只能說你中有我我中有你，但你不屬于我我也不屬于你

五、Monte Carlo

5.1 蒙特卡羅方法引入

蒙特卡洛（Monte Carlo）方法來自于摩納哥的蒙特卡洛賭場，許多紙牌類游戲需要計算其勝利的概率。用它作為名字大概是因?yàn)槊商乜_方法是一種隨機(jī)模擬的方法，這很像賭博場里面的扔骰子的過程。最早的蒙特卡羅方法都是為了求解一些不太好求解的求和或者積分問題。比如積分：
如果我們很難求解出f(x)的原函數(shù)，那么這個積分比較難求解。當(dāng)然我們可以通過蒙特卡羅方法來模擬求解近似值。(我們可以將蒙特卡洛理解為簡單的模擬，通過模擬的情景來計算其發(fā)生的概率。)
如何模擬呢？假設(shè)我們函數(shù)圖像如下圖:

5.2 概率分布采樣

上一節(jié)我們講到蒙特卡羅方法的關(guān)鍵是得到x的概率分布。如果求出了x的概率分布，我們可以基于概率分布去采樣基于這個概率分布的n個x的樣本集，帶入蒙特卡羅求和的式子即可求解。但是還有一個關(guān)鍵的問題需要解決，即如何基于概率分布去采樣基于這個概率分布的n個x的樣本集。

5.3 接受-拒絕采樣

5.4 蒙特卡羅方法小結(jié)

使用接受-拒絕采樣，我們可以解決一些概率分布不是常見的分布的時候，得到其采樣集并用蒙特卡羅方法求和的目的。但是接受-拒絕采樣也只能部分滿足我們的需求，在很多時候我們還是很難得到我們的概率分布的樣本集。比如：

從上面可以看出，要想將蒙特卡羅方法作為一個通用的采樣模擬求和的方法，必須解決如何方便得到各種復(fù)雜概率分布的對應(yīng)的采樣樣本集的問題。而馬爾科夫鏈就是幫助找到這些復(fù)雜概率分布的對應(yīng)的采樣樣本集的白衣騎士。

5.5 應(yīng)用場景

通常情況下，許多概率的計算非常復(fù)雜甚至是無法計算的，但是我們可以通過計算機(jī)對期待發(fā)生的場景進(jìn)行大量的模擬，從而計算出其發(fā)生的概率。簡單的公式表達(dá)為：

學(xué)習(xí)過概率的同學(xué)應(yīng)當(dāng)知道這就是典型的頻率學(xué)派的觀點(diǎn)。

對于蒙特卡洛方法，經(jīng)典的例子就是計算 $π\(zhòng)pi$ 值。在（-1，1）之間隨機(jī)取兩個數(shù)，如果在單位圓內(nèi)，則記一次，在圓外則不計入次數(shù)。

import random import numpy as np from math import sqrt, pinum1 = 1000 # simulation times freq1 = 0 x1 = [] y1 = [] for i in range(1,num1+1):x_i, y_i = random.uniform(-1, 1), random.uniform(-1, 1)square_distance = x_i**2+y_i**2x1.append(x_i)y1.append(y_i)if square_distance <= 1: # 圓的方程x^2+y^2=1freq1 += 1simulation_pi1 = 4*freq1/num1;num2 = 10000 # simulation times freq2 = 0 x2 = [] y2 = [] for i in range(num2):x_i, y_i = random.uniform(-1, 1), random.uniform(-1, 1)square_distance = x_i**2+y_i**2x2.append(x_i)y2.append(y_i)if square_distance <= 1: # 圓的方程x^2+y^2=1freq2 += 1simulation_pi2 = 4*freq2/num2;% matplotlib inline import matplotlib.pyplot as pltfig, ax = plt.subplots(figsize=(8,8))theta = np.linspace(0,2*pi,500) x,y = np.cos(theta), np.sin(theta) ax.plot(x, y, color='red', linewidth=1.0)ax.scatter(x1,y1,color='blue',alpha=0.75) ax.set_xlim(-1,1); ax.set_ylim(-1,1);

結(jié)果為：

真實(shí)的pi值： 3.141592653589793
蒙特卡洛方法下得出的pi值（模擬1000次）： 3.188
相對誤差（模擬1000次）： 1.4771917153924754 %
蒙特卡洛方法下得出的pi值（模擬10000次）： 3.1292
相對誤差（模擬10000次）： 0.39447041536821975 %

5.6 小結(jié)

蒙特卡洛就是一種模擬事情發(fā)生的方法。

六、MCMC

從名字我們可以看出，MCMC由兩個MC組成，即蒙特卡羅方法（Monte Carlo Simulation，簡稱MC）和馬爾科夫鏈（Markov Chain ，也簡稱MC）。
Monte Carlo （蒙特卡羅）的核心是尋找一個隨機(jī)的序列。

蒙特卡洛馬爾可夫鏈 Markov Chain Monte Carlo簡稱MCMC，是一個抽樣方法，用于解決難以直接抽樣的分布的隨機(jī)抽樣模擬問題。用來在概率空間，通過隨機(jī)采樣估算興趣參數(shù)的后驗(yàn)分布。

那MCMC到底是什么呢？《告別數(shù)學(xué)公式，圖文解讀什么是馬爾可夫鏈蒙特卡羅方法》里面這樣解釋：MCMC方法是用來在概率空間，通過隨機(jī)采樣估算興趣參數(shù)的后驗(yàn)分布。
說的很玄，蒙特卡羅本來就可以采樣，馬爾科夫鏈可以采樣，為啥要將他們合在一起？下面給出兩個動機(jī)，后面將從蒙特卡羅開始一直推到gibbs采樣，來深入了解為什么需要MCMC。
再次感謝劉建平MCMC，他是網(wǎng)上寫的最詳細(xì)的——將整個脈絡(luò)梳理出來了，看完收獲很多。本文幾乎涵蓋了它所有內(nèi)容，因此只能算一篇讀書筆記。

簡而言之，就是用馬爾科夫鏈方法來抽樣。

抽樣算法的主要任務(wù)是找到符合給定分布p(x)的一系列樣本。對于簡單的分布，可以通過基本的抽樣算法進(jìn)行抽樣。大多數(shù)分布都是不容易直接抽樣的，馬爾可夫鏈蒙特卡羅算法解決了不能通過簡單抽樣算法進(jìn)行抽樣的問題，是一種重要的實(shí)用性很強(qiáng)的抽樣算法。

馬爾可夫鏈蒙特卡羅算法（簡寫為MCMC）的核心思想是找到某個狀態(tài)空間的馬爾可夫鏈，使得該馬爾可夫鏈的穩(wěn)定分布就是我們的目標(biāo)分布p(x)。這樣我們在該狀態(tài)空間進(jìn)行隨機(jī)游走的時候，每個狀態(tài)x的停留時間正比于目標(biāo)概率p(x)。在用MCMC進(jìn)行抽樣的時候，我們首先引進(jìn)一個容易抽樣的參考分布q(x)，在每步抽樣的過程中從q(x)里面得到一個候選樣本y, 然后按照一定的原則決定是否接受該樣本，該原則的確定就是要保證我們得到的原本恰好服從p(x)分布。

在了解什么是MCMC的時候，我們需要考慮一個情況。舉例，我們已經(jīng)知道一個分布（例如beta分布）的概率密度函數(shù)PDF，那么我們怎么樣從這個分布中提取樣本呢？

MCMC給我們提供了一種方式來從概率分布中進(jìn)行采樣，而這種方法在貝葉斯統(tǒng)計中的后驗(yàn)概率進(jìn)行采樣十分方便。

6.1 動因

MCMC通常用于解決高維度的積分和最優(yōu)化問題，這兩種問題也是機(jī)器學(xué)習(xí)、物理、統(tǒng)計、計量等領(lǐng)域的基礎(chǔ)。比如：

6.2 MCMC的定義

Wikipedia的解釋如下:

Markov Chain Monte Carlo (MCMC) methods are a class of algorithms for sampling from a probability distribution based on constructing a Markov chain that has the desired distribution as its stationary distribution. The state of the chain after a number of steps is then used as a sample of the desired distribution. The quality of the sample improves as a function of the number of steps.
馬爾可夫鏈蒙特卡羅方法是一類以期望分布為平穩(wěn)分布的馬爾可夫鏈為基礎(chǔ)，對概率分布進(jìn)行抽樣的算法。經(jīng)過一系列步驟之后，鏈的狀態(tài)將用作期望分布的樣本。樣品的質(zhì)量隨著步驟數(shù)量的增加而提高。

MCMC方法提供了可以創(chuàng)建一個以Beta分布為其平穩(wěn)分布（stationary distribution）的馬爾科夫鏈的算法，從而使取樣變得簡單，因?yàn)槲覀兛梢詮囊粋€均勻分布中取樣(這相對容易)。

如果我們基于一些算法，重復(fù)地從一個隨機(jī)的狀態(tài)開始，遍歷到下一個狀態(tài)，我們將會創(chuàng)建一個馬爾可夫鏈，一個以beta分布作為其平穩(wěn)分布的馬爾可夫鏈。

這類算法，一個典型的代表就是Metropolis-Hastings Algorithm算法。

MCMC由梅特羅波利斯(Metropolis)于1949年基于馬爾可夫鏈的基本性質(zhì)提出，下面介紹一下與馬爾可夫鏈相關(guān)的性質(zhì)。

6.3 相關(guān)性質(zhì)

一、穩(wěn)定分布

穩(wěn)定分布是指當(dāng)我們給定狀態(tài)空間的一個初始分布 $π0\pi_0$ 以后，按照轉(zhuǎn)移矩陣進(jìn)行跳轉(zhuǎn)最終達(dá)到的穩(wěn)定狀態(tài)。

即每個狀態(tài)的流出概率等于該狀態(tài)注入的概率，該方程稱為全局平衡方程。滿足該方程的分布稱為轉(zhuǎn)移矩陣A的穩(wěn)定分布。

二、細(xì)致平衡

對于一個穩(wěn)定分布，我們可以給其一個更強(qiáng)的條件限制，使得任意一個狀態(tài) $i$ 滿足如下條件，

下面我們介紹基本的梅特羅波利斯算法(Metropolis algorithm)。假設(shè)p(x)是我們的目標(biāo)分布函數(shù)，我們想得到一系列服從該分布的樣本。我們考慮樣本的產(chǎn)生過程構(gòu)成一個馬爾可夫鏈，并且讓p(x)是該馬爾可夫鏈的穩(wěn)定分布，那么該樣本序列就服從p(x)?，F(xiàn)在p(x)是已知的，問題的關(guān)鍵在于如何構(gòu)造這個產(chǎn)生過程，也即如何構(gòu)造馬爾可夫鏈。

即目標(biāo)概率p(x)滿足細(xì)致平衡方程，因此p(x)是轉(zhuǎn)移概率的穩(wěn)定分布。

最后，回顧一下梅特羅波利斯抽樣的主要思想。我們首先構(gòu)造了一個馬爾可夫鏈，接著證明了p(x)滿足其細(xì)致平衡方程，進(jìn)而說明p(x)是該鏈的穩(wěn)定分布。然后將抽樣的過程看成是在馬爾可夫鏈狀態(tài)空間進(jìn)行跳轉(zhuǎn)的過程，跳轉(zhuǎn)的候選狀態(tài)由參考分布q(x)產(chǎn)生。最后得到一個的跳轉(zhuǎn)序列，該序列在每個狀態(tài)x的停留時間與p(x)成比，即服從p(x)分布。

6.4 Metropolis-Hastings (MH) Algorithm算法

梅特羅波利斯－黑斯廷斯算法（英語：Metropolis–Hastings algorithm）是統(tǒng)計學(xué)與統(tǒng)計物理中的一種馬爾科夫蒙特卡洛（MCMC）方法，用于在難以直接采樣時從某一概率分布中抽取隨機(jī)樣本序列。得到的序列可用于估計該概率分布或計算積分（如期望值）等。梅特羅波利斯－黑斯廷斯或其他MCMC算法一般用于從多變量（尤其是高維）分布中采樣。對于單變量分布而言，常會使用自適應(yīng)判別采樣（adaptive rejection sampling）等其他能抽取獨(dú)立樣本的方法，而不會出現(xiàn)MCMC中樣本自相關(guān)的問題。

該算法的名稱源于美國物理學(xué)家尼古拉斯·梅特羅波利斯與加拿大統(tǒng)計學(xué)家W·K·黑斯廷斯

為了更形象地理解這個算法，我們用下面這個例子來類比。假設(shè)我們想知道某個湖的水容量以及這個湖中最深的點(diǎn)，湖水很渾濁以至于沒法通過肉眼來估計深度，而且這個湖相當(dāng)大，網(wǎng)格近似法顯然不是個好辦法。為了找到一個采樣策略，我們請來了兩個好朋友小馬和小萌。經(jīng)過討論之后想出了如下辦法，我們需要一個船（當(dāng)然，也可以是竹筏）和一個很長的棍子，這比聲吶可便宜多了，而且我們已經(jīng)把有限的錢都花在了船上。

（1）隨機(jī)選一個點(diǎn)，然后將船開過去。

（2）用棍子測量湖的深度。

（3）將船移到另一個地點(diǎn)并重新測量。

（4）按如下方式比較兩點(diǎn)的測量結(jié)果。

如果新的地點(diǎn)比舊的地點(diǎn)水位深，那么在筆記本上記錄下新的測量值并重復(fù)過程（2）。
如果新的地點(diǎn)比舊的地點(diǎn)水位淺，那么我們有兩個選擇：接受或者拒絕。接受意味著記錄下新的測量值并重復(fù)過程（2）；拒絕意味著重新回到上一個點(diǎn)，再次記錄下上一個點(diǎn)的測量值。

如何決定接受還是拒絕新的測量值呢？這里的一個技巧便是使用Metropolis-Hastings準(zhǔn)則，即接受新的測量值的概率正比于新舊兩點(diǎn)的測量值之比。

事實(shí)上，理論保證了在這種情形下，如果我們能采樣無數(shù)次，最終能得到完整的后驗(yàn)。幸運(yùn)地是，實(shí)際上對于很多問題而言，我們只需要相對較少地采樣就可以得到一個相當(dāng)準(zhǔn)確的近似。

6.5 數(shù)學(xué)表達(dá)

6.6 算法介紹

6.7 MCMC 采樣

6.7.1 背景

給定一個的概率分布 P(x), 我們希望產(chǎn)生服從該分布的樣本。

前面介紹過一些隨機(jī)采樣算法（如拒絕采樣、重要性采樣）可以產(chǎn)生服從特定分布的樣本，但是這些采樣算法存在一些缺陷（如難以選取合適的建議分布，只適合一元隨機(jī)變量等）。

下面將介紹一種更有效的隨機(jī)變量采樣方法：MCMC 和 Gibbs采樣，這兩種采樣方法不僅效率更高，而且適用于多元隨機(jī)變量的采樣。

6.7.2 隨機(jī)矩陣

在MCMC采樣中先隨機(jī)一個狀態(tài)轉(zhuǎn)移矩陣Q，然而該矩陣不一定能滿足細(xì)致平穩(wěn)定理，一次會做一些改進(jìn)，具體過程如下

6.7.3 算法具體流程

MCMC采樣算法的具體流程如下

6.7.4 MCMC: Metropolis-Hastings algorithm

然而關(guān)于MCMC采樣有收斂太慢的問題，所以在MCMC的基礎(chǔ)上進(jìn)行改進(jìn)，引出M-H采樣算法

M-H 算法的具體流程如下

M-H算法在高維時同樣適用

6.7.4.1 基本概念

從一個概率分布（目標(biāo)分布 $P (x)$ ）中得到隨機(jī)樣本序列
這個序列可以用于：a) 近似估計分布 $P (x)$ ; b) 計算積分(期望)
用于高維分布取樣
缺陷: MCMC固有缺點(diǎn), 樣本自相關(guān)性

6.7.4.2 優(yōu)勢

可以從任意的概率分布 $P (x)$ 中取樣，只要滿足條件：函數(shù)f(x)成比例于 $P (x)$ 的密度。
更寬松的要求： $f (x)$ 僅需要與 $P (x)$ 的密度成比例。

6.7.4.3 要點(diǎn)

生成樣本值序列；樣本值產(chǎn)生得越多，這些值的分布就越近似于 $P (x)$
迭代產(chǎn)生樣本值：下一個樣本的分布僅僅取決于當(dāng)前樣本值（馬爾科夫鏈特性）
接受/拒絕概率：接受計算出的值為下一個樣本值/拒絕并重復(fù)使用當(dāng)前樣本值；基于 $P (x)$ ，接受概率通過比較 $f(x_t)$ （當(dāng)前值）和 $f (x^{'})$ （備選值）得出

6.7.5 Metropolis Algorithm (對稱分布)

Input: $f (x)$ ，與目標(biāo)分布 $P (x)$ 成比例的函數(shù)

6.7.5.1 初始化：

6.7.5.2 迭代 $t$

6.7.6 缺陷

樣本自相關(guān)：相鄰的樣本會相互相關(guān)，雖然可以通過每 $n$ 步取樣的方式來減少相關(guān)性，但這樣的后果就是很難讓樣本近似于目標(biāo)分布 $P (x)$ 。

自相關(guān)性可通過增加步調(diào)長度（jumping width，與jumping function

g (x ∣ y)

的方差有關(guān)）來控制，但同時也增加了拒絕備選樣本的幾率

α\alpha

。

過大或過小的jumping size會導(dǎo)致slow-mixing Markov Chain, 即高度自相關(guān)的一組樣本，以至于我們需要得到非常大的樣本量

n

才能得到目標(biāo)分布

P (x)

。

初始值的選擇：盡管Markov Chain最后都會收斂到目標(biāo)分布，然而初始值的選擇直接影響到運(yùn)算時間，尤其是把初始值選在在了“低密度”區(qū)域。因此，選擇初值時，最好加入一個“burn-in period”（預(yù)燒期，預(yù)選期）。

6.7.7 優(yōu)勢

抗“高維魔咒”（curse of dimensionality）：維度增加，對于rejection sampling方法來說，拒絕的概率就是呈指數(shù)增長。而MCMC則成為了解決這種問題的唯一方法
多元分布中，為了避免多元初始值以及 $g (x ∣ y)$ 選擇不當(dāng)而導(dǎo)致的問題，Gibbs sampling是另外一個更適合解決多元分布問題的MCMC 方法。Gibbs sampling從多元分布的各個維度中分別選擇初始值，然后這些變量分別同時取樣。

6.7.8 衍生算法

6.7.9 小結(jié)

一般來說M-H采樣算法較MCMC算法應(yīng)用更廣泛，然而在大數(shù)據(jù)時代，M-H算法面臨著兩個問題：

1）在高維時的計算量很大，算法效率很低，同時存在拒絕轉(zhuǎn)移的問題，也會加大計算量

2）由于特征維度大，很多時候我們甚至很難求出目標(biāo)的各特征維度聯(lián)合分布，但是可以方便求出各個特征之間的條件概率分布（因此就思考是否能只知道條件概率分布的情況下進(jìn)行采樣）。

6.8 Gibbs 采樣

6.8.1 二維的流程

因此可以得出在二維的情況下Gibbs采樣算法的流程如下

6.8.2 多維

而在多維的情況下，比如一個n維的概率分布π(x1, x2, …xn)，我們可以通過在n個坐標(biāo)軸上輪換采樣，來得到新的樣本。

對于輪換到的任意一個坐標(biāo)軸xi上的轉(zhuǎn)移，馬爾科夫鏈的狀態(tài)轉(zhuǎn)移概率為 P(xi|x1, x2, …, xi?1, xi+1, …, xn)，即固定n?1個坐標(biāo)軸，在某一個坐標(biāo)軸上移動。

而在多維的情況下Gibbs采樣算法的流程如下

6.8.3 小結(jié)

由于Gibbs采樣在高維特征時的優(yōu)勢，目前我們通常意義上的MCMC采樣都是用的Gibbs采樣。

當(dāng)然Gibbs采樣是從M-H采樣的基礎(chǔ)上的進(jìn)化而來的，同時Gibbs采樣要求數(shù)據(jù)至少有兩個維度，一維概率分布的采樣是沒法用Gibbs采樣的，這時M-H采樣仍然成立。

6.9 其他算法

除了最常見的MH那幾個算法，后來還有很多新的比較驚艷的算法出現(xiàn)，比如說slice sampling，elliptical slice sampling，generalized elliptical slice sampling，上面說的BPS， forward event chain MC，還有和神經(jīng)網(wǎng)絡(luò)結(jié)合的NNGHMC，A-Nice-MC，以及利用了batch optimization思想的stochastic gradient HMC以及stochastic gradient Langevin dynamic等。

6.10 參考資料

以下為列表，鏈接見原文

統(tǒng)計之都-MCMC

HANS-MCMC 算法及其應(yīng)用

知乎-MCMC 專欄

機(jī)器學(xué)習(xí)之MCMC算法

知乎-MCMC 算法

知乎-MCMC 算法中接受概率是什么意思

MCMC 和 Metropolis–Hastings 算法

馬爾可夫鏈蒙特卡洛(MCMC)算法

CSDN-MCMC

MCMC相關(guān)算法介紹及代碼實(shí)現(xiàn)

算法資料
http://civs.stat.ucla.edu/MCMC/MCMC_tutorial.htm

http://www.soe.ucsc.edu/classes/cmps290c/Winter06/paps/mcmc.pdf

http://public.lanl.gov/kmh/talks/maxent00b.pdf

http://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo

google keywords： MCMC tutorial

MCMC preprint service:

http://www.statslab.cam.ac.uk/~mcmc/

David MacKay’s book (electronic version availiable):

http://www.inference.phy.cam.ac.uk/mackay/itila/

Radford M. Neal’s review: Probabilistic Inference using Markov Chain Monte Carlo Methods

http://www.cs.toronto.edu/~radford/review.abstract.html

6.11 代碼實(shí)現(xiàn)MCMC采樣beta分布

import random # Lets define our Beta Function to generate s for any particular state. # We don't care for the normalizing constant here.def beta_s(w,a,b): # beta distribution pdfreturn w**(a-1)*(1-w)**(b-1)# This Function returns True if the coin with probability P of heads comes heads when flipped. def random_coin(p):unif = random.uniform(0,1)if unif>=p:return Falseelse:return True# This Function runs the MCMC chain for Beta Distribution. def beta_mcmc(N_hops,a,b):states = []cur = random.uniform(0,1)for i in range(0,N_hops):states.append(cur)next = random.uniform(0,1)ap = min(beta_s(next,a,b)/beta_s(cur,a,b),1) # Calculate the acceptance probabilityif random_coin(ap):cur = nextreturn states[-1000:] # Returns the last 1000 states of the chain

在定義了相關(guān)函數(shù)后，我們來檢查一下我們采樣的結(jié)果和實(shí)際beta概率分布之間的差距

import numpy as np from matplotlib import pylab as pl import scipy.special as ss % matplotlib inline plt.rcParams['figure.figsize'] = (17.0, 4.0)# Actual Beta PDF. def beta(a, b, i):e1 = ss.gamma(a + b)e2 = ss.gamma(a)e3 = ss.gamma(b)e4 = i ** (a - 1)e5 = (1 - i) ** (b - 1)return (e1/(e2*e3)) * e4 * e5# Create a function to plot Actual Beta PDF with the Beta Sampled from MCMC Chain. def plot_beta(a, b):Ly = []Lx = []i_list = np.mgrid[0:1:100j]for i in i_list:Lx.append(i)Ly.append(beta(a, b, i))pl.plot(Lx, Ly, label="Real Distribution: a="+str(a)+", b="+str(b))plt.hist(beta_mcmc(100000,a,b),density=True,bins =25, histtype='step',label="Simulated_MCMC: a="+str(a)+", b="+str(b))pl.legend()pl.show()plot_beta(0.1, 0.1) plot_beta(1, 1) plot_beta(2, 3)

相關(guān)圖形如下所示：

從上面的采樣值可以看出，我們的采樣值還是很接近beta分布本身。雖然我們用了beta分布最為例子，但實(shí)際上其他的分布也是可以類似來進(jìn)行采樣。

針對那些我們無法直接利用共軛分布得到的后驗(yàn)分布，特別是高維的分布，MCMC方法更顯的尤為重要。

6.12 MH Algorithm的進(jìn)一步說明

從貝葉斯的角度看，Metropolis-Hastings算法使得我們能夠從任意分布中以概率p(x)得到采樣值，只要我們能算出某個與p(x)成比例的值。這一點(diǎn)很有用，因?yàn)樵陬愃曝惾~斯統(tǒng)計的許多問題中，最難的部分是計算歸一化因子，也就是貝葉斯定理中的分母。Metropolis-Hastings算法的步驟如下：

最后，我們會得到一連串記錄值，有時候也稱采樣鏈或者跡。如果一切都正常進(jìn)行，那么這些采樣值就是后驗(yàn)的近似。在采樣鏈中出現(xiàn)次數(shù)最多的值就是對應(yīng)后驗(yàn)中最可能的值。該過程的一個優(yōu)點(diǎn)是：后驗(yàn)分析很簡單，我們把對后驗(yàn)求積分的過程轉(zhuǎn)化成了對采樣鏈所構(gòu)成的向量求和的過程。

強(qiáng)烈建議閱讀以下博文，作者用一個簡單的例子實(shí)現(xiàn)了metropolis方法，并將其用于求解后驗(yàn)分布，文中用非常好看的圖展示了采樣的過程，同時簡單討論了最初選取的步長是如何影響結(jié)果的。
https://github.com/findmyway/Bayesian-Analysis-with-Python/blob/master/MCMC-sampling-for-dummies.ipynb

6.13 一句話總結(jié)MCMC

MCMC generates samples from the posterior distribution by constructing a reversible Markov-chain that has as its equilibrium distribution the target posterior distribution.
https://twiecki.io/blog/2015/11/10/mcmc-sampling/

6.14 MCMC如何解決具體問題

尋找p(x)最大值（simulated annealing）
轉(zhuǎn)換核的混搭和循環(huán)
Gibbs 抽樣
Monte Carlo EM（期望-最大化）
輔助變量抽樣（Auxiliary variable samplers)
可逆跳躍MCMC（Reversible jump MCMC）

-----------------------------------**
鏈接：
https://zhuanlan.zhihu.com/p/83314877
https://www.zhihu.com/question/26887292/answer/310904395
https://zh.wikipedia.org/wiki/%E9%9E%85_(%E6%A6%82%E7%8E%87%E8%AE%BA)
https://www.zhihu.com/question/52778636/answer/133133860
https://www.zhihu.com/question/31173033/answer/113202955
https://zhuanlan.zhihu.com/p/116725922
https://zhuanlan.zhihu.com/p/86995916
https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm
https://zhuanlan.zhihu.com/p/37121528
https://houbb.github.io/2020/01/28/math-05-mcmc
https://zhuanlan.zhihu.com/p/21112618

總結(jié)

以上是生活随笔為你收集整理的Probability, Matringale,Markov Chain, MCMC的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。

上一篇： MATLAB 将不同长度的一维数组存入二
下一篇：马尔可夫蒙特卡罗 MCMC 原理及经典实