日韩性视频-久久久蜜桃-www中文字幕-在线中文字幕av-亚洲欧美一区二区三区四区-撸久久-香蕉视频一区-久久无码精品丰满人妻-国产高潮av-激情福利社-日韩av网址大全-国产精品久久999-日本五十路在线-性欧美在线-久久99精品波多结衣一区-男女午夜免费视频-黑人极品ⅴideos精品欧美棵-人人妻人人澡人人爽精品欧美一区-日韩一区在线看-欧美a级在线免费观看

歡迎訪問 生活随笔!

生活随笔

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

编程问答

Probability, Matringale,Markov Chain, MCMC

發(fā)布時(shí)間:2023/12/2 编程问答 33 豆豆
生活随笔 收集整理的這篇文章主要介紹了 Probability, Matringale,Markov Chain, MCMC 小編覺得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.

一、基本知識

1. 條件概率

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

2. 條件期望

在上述概率下的期望我們稱之為條件期望。
計(jì)算:

3. 隨機(jī)過程

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

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

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

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

二、鞅(martingale)

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

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

鞅的概念首先是由保羅·皮埃爾·萊維于 1934 年提出的,但他只提出了離散時(shí)間的版本,而且沒有給予命名。直到 1939 年,約翰?維爾將此概念推廣到連續(xù)時(shí)間的情況,并且首次提出 martingale 這個(gè)名稱。約瑟夫·利奧·杜布(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個(gè)離散階段下,鞅的相關(guān)定義和性質(zhì)。

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

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

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

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

2.1 Martingale的定義

我們開始討論鞅的定義,鞅有以下三種等價(jià)的定義:

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

2.2 Martingale, supermartingale, submartingale

鞅:E[Xt∣Fs]=Xs(t>s)E[X_t | F_s] = X_s (t>s)E[Xt?Fs?]=Xs?t>s
上鞅:E[Xt∣Fs]<Xs(t>s)E[X_t | F_s] < X_s (t>s)E[Xt?Fs?]<Xs?t>s
下鞅:E[Xt∣Fs]>Xs(t>s)E[X_t | F_s] > X_s (t>s)E[Xt?Fs?]>Xs?t>s。

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

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

2.3 Stopping time

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

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

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

停時(shí)鞅的概念引出了一系列定理,例如可選停止定理(又稱可選抽樣定理):在特定條件下,停時(shí)的鞅的期望值等于其初始值。利用這一定理,我們可以證明對于一個(gè)壽命有限且房產(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è)你有一個(gè)系統(tǒng),這個(gè)系統(tǒng)具有MMM 種可能的狀態(tài),并且在狀態(tài)之間正在移動。真實(shí)世界中,我們常常見到這種例子,例如天氣從炎熱到溫和,或者是股票市場從熊市到牛市再到停滯的狀態(tài)(stagnant states)。

3.1.1 馬爾可夫性的定義

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

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

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


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

3.2 Markov Processes馬爾可夫過程

在概率論及統(tǒng)計(jì)學(xué)中,馬爾可夫過程(英語:Markov process)是一個(gè)具備了馬爾可夫性質(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)系,下一個(gè)狀態(tài)只決定于當(dāng)前的狀態(tài).通俗點(diǎn)說就是:你下一時(shí)刻在哪只跟你現(xiàn)在在哪有關(guān),于你以前去過哪沒有關(guān)系. 如下圖:

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

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

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

  • Irreducibility. 即圖是聯(lián)通的,各個(gè)狀態(tài)之間都有過去的辦法,舉個(gè)不聯(lián)通的例子,比如爬蟲爬不到內(nèi)部局域網(wǎng)的網(wǎng)頁
  • Aperiodicity. 非周期性, 即圖中遍歷不會陷入到一個(gè)死圈里,進(jìn)去了再也出不來,有些網(wǎng)站為了防機(jī)器人,會專門設(shè)置這種陷阱
  • Detailed Balance,這是保證系統(tǒng)有穩(wěn)態(tài)的一個(gè)重要條件,詳細(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 錢,以此類推,手動計(jì)算,可以發(fā)現(xiàn)下一個(gè)狀態(tài)每個(gè)人手中的錢都沒有變。值得說明的是,這里體現(xiàn)了一個(gè)很重要的特性,那就是從一個(gè)高概率狀態(tài) xi 向一個(gè)低概率狀態(tài) x(i-1) 轉(zhuǎn)移的概率等于從這個(gè)低概率狀態(tài)向高概率狀態(tài)轉(zhuǎn)移的概率(reversible,至于要不要轉(zhuǎn)移又是另外一回事)
  • ergodicity, 遍歷性,是指不管事物現(xiàn)在出于什么狀態(tài),在較長時(shí)間內(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ā)生的概率分別都是多少?
先看個(gè)假設(shè),他每個(gè)狀態(tài)的轉(zhuǎn)移都是有概率的,比如今天玩,明天睡的概率是幾,今天玩,明天也玩的概率是幾幾,看圖更清楚一點(diǎn)。

這個(gè)矩陣就是轉(zhuǎn)移概率矩陣P,并且它是保持不變的,就是說第一天到第二天的轉(zhuǎn)移概率矩陣跟第二天到第三天的轉(zhuǎn)移概率矩陣是一樣的。(這個(gè)叫時(shí)齊,有興趣的同學(xué)自行百度)。
有了這個(gè)矩陣,再加上已知的第一天的狀態(tài)分布,就可以計(jì)算出第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 (看見沒,只跟它前面一個(gè)狀態(tài)Sn-1有關(guān))。

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

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

附:S2 的計(jì)算過程 (沒興趣的同學(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),而與以前時(shí)間的狀態(tài)無關(guān)

3.3.2.1 C-K方程

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

3.3.3 狀態(tài)分類

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

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

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

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

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

馬爾科夫鏈的分解:

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

常返類的周期:

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



命題得證。

很顯然,暫態(tài)也是一個(gè)類性質(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

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

我們用P表示這個(gè)轉(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)分布。對于這個(gè)特性,有如下精彩定理

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

條件

(1)非周期馬氏鏈:馬氏鏈轉(zhuǎn)移要收斂,就一定不能是周期性的。不做特別處理,我們處理的問題基本上都是非周期性的,在此不做多余解釋。
(2)存在概率轉(zhuǎn)移矩陣P,任意兩個(gè)狀態(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)移后,會存在一個(gè)穩(wěn)定的狀態(tài)π。
(2)概率轉(zhuǎn)移矩陣自乘足夠多次后,每行值相等。即

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

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

假設(shè)一開始有100個(gè)農(nóng)村人,0個(gè)城市人,每代轉(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的變化)。同時(shí)如果我們知道了中國城市化人口比例,我們就能知道城市化進(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 吸收概率和吸收期望時(shí)間


3.3.7.1 平均時(shí)間

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

平均吸收時(shí)間計(jì)算:解平均吸收時(shí)間方程組

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

因?yàn)樗撵o態(tài)分布(或叫平穩(wěn)性分布,stationary distribution)性質(zhì)。
我們用一下Wiki上給出的股票市場的例子(見下圖),這是一個(gè)具有靜態(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)分布一定是相同的。那么得出這個(gè)靜態(tài)分布有什么用處呢?

靜止?fàn)顟B(tài)分布的重要性在于它可以讓你在隨機(jī)的時(shí)間里,對于一個(gè)系統(tǒng)的任意一個(gè)狀態(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.)

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

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

3.3.9 小結(jié)

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

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

3.4.1 定義

3.4.1.1 定義1:一般定義


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

3.4.1.2 定義2

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

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

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

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

3.4.1.3 定義3

3.4.1.4 舉例


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

3.4.3 CTMC的極限概率

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

3.4.5 時(shí)間可逆性

3.4.6 截止的馬爾科夫鏈


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

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

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

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

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

所以從這兩個(gè)角度來理解,兩個(gè)名詞

有共性:

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

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

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

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

五、Monte Carlo

5.1 蒙特卡羅方法引入

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

5.2 概率分布采樣

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

5.3 接受-拒絕采樣

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

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

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

5.5 應(yīng)用場景

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

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

對于蒙特卡洛方法,經(jīng)典的例子就是計(jì)算π\(zhòng)piπ值。在(-1,1)之間隨機(jī)取兩個(gè)數(shù),如果在單位圓內(nèi),則記一次,在圓外則不計(jì)入次數(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由兩個(gè)MC組成,即蒙特卡羅方法(Monte Carlo Simulation,簡稱MC)和馬爾科夫鏈(Markov Chain ,也簡稱MC)。
Monte Carlo (蒙特卡羅)的核心是尋找一個(gè)隨機(jī)的序列。

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

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

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

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

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

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

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

6.1 動因

MCMC通常用于解決高維度的積分和最優(yōu)化問題,這兩種問題也是機(jī)器學(xué)習(xí)、物理、統(tǒng)計(jì)、計(jì)量等領(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)建一個(gè)以Beta分布為其平穩(wěn)分布(stationary distribution)的馬爾科夫鏈的算法,從而使取樣變得簡單,因?yàn)槲覀兛梢詮囊粋€(gè)均勻分布中取樣(這相對容易)。

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

這類算法,一個(gè)典型的代表就是Metropolis-Hastings Algorithm算法。

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

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

一、穩(wěn)定分布

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

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

二、細(xì)致平衡

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

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


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

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

6.4 Metropolis-Hastings (MH) Algorithm算法

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

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



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

(1)隨機(jī)選一個(gè)點(diǎn),然后將船開過去。

(2)用棍子測量湖的深度。

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

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

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

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

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

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

6.6 算法介紹



6.7 MCMC 采樣

6.7.1 背景

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

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

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

6.7.2 隨機(jī)矩陣

在MCMC采樣中先隨機(jī)一個(gè)狀態(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算法在高維時(shí)同樣適用

6.7.4.1 基本概念

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

6.7.4.2 優(yōu)勢

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

6.7.4.3 要點(diǎn)

  • 生成樣本值序列;樣本值產(chǎn)生得越多,這些值的分布就越近似于P(x)P(x)P(x)
  • 迭代產(chǎn)生樣本值:下一個(gè)樣本的分布僅僅取決于當(dāng)前樣本值(馬爾科夫鏈特性)
  • 接受/拒絕概率:接受計(jì)算出的值為下一個(gè)樣本值/拒絕并重復(fù)使用當(dāng)前樣本值;基于P(x)P(x)P(x),接受概率通過比較f(xt)f(x_t)f(xt?)(當(dāng)前值)和f(x′)f(x')f(x)(備選值)得出

6.7.5 Metropolis Algorithm (對稱分布)

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

6.7.5.1 初始化:

6.7.5.2 迭代ttt

6.7.6 缺陷

  • 樣本自相關(guān):相鄰的樣本會相互相關(guān),雖然可以通過每nnn步取樣的方式來減少相關(guān)性,但這樣的后果就是很難讓樣本近似于目標(biāo)分布P(x)P(x)P(x)。
  • 自相關(guān)性可通過增加步調(diào)長度(jumping width,與jumping function g(x∣y)g(x|y)g(xy)的方差有關(guān))來控制,但同時(shí)也增加了拒絕備選樣本的幾率α\alphaα。
  • 過大或過小的jumping size會導(dǎo)致slow-mixing Markov Chain, 即高度自相關(guān)的一組樣本,以至于我們需要得到非常大的樣本量nnn才能得到目標(biāo)分布P(x)P(x)P(x)
    • 初始值的選擇:盡管Markov Chain最后都會收斂到目標(biāo)分布,然而初始值的選擇直接影響到運(yùn)算時(shí)間,尤其是把初始值選在在了“低密度”區(qū)域。因此,選擇初值時(shí),最好加入一個(gè)“burn-in period”(預(yù)燒期,預(yù)選期)。

    6.7.7 優(yōu)勢

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

    6.7.8 衍生算法


    6.7.9 小結(jié)

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

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

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

    6.8 Gibbs 采樣


    6.8.1 二維的流程

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

    6.8.2 多維

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

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

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

    6.8.3 小結(jié)

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

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

    6.9 其他算法

    除了最常見的MH那幾個(gè)算法,后來還有很多新的比較驚艷的算法出現(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)計(jì)之都-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)得到采樣值,只要我們能算出某個(gè)與p(x)成比例的值。這一點(diǎn)很有用,因?yàn)樵陬愃曝惾~斯統(tǒng)計(jì)的許多問題中,最難的部分是計(jì)算歸一化因子,也就是貝葉斯定理中的分母。Metropolis-Hastings算法的步驟如下:

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

    強(qiáng)烈建議閱讀以下博文,作者用一個(gè)簡單的例子實(shí)現(xiàn)了metropolis方法,并將其用于求解后驗(yàn)分布,文中用非常好看的圖展示了采樣的過程,同時(shí)簡單討論了最初選取的步長是如何影響結(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)容還不錯(cuò),歡迎將生活随笔推薦給好友。