1 引言

?蒙特卡洛積分通過生成符合特定分布的隨機變量來近似計算積分值。例如:

?上面的式子可以理解成求函數 f(x)f(x) 的期望因此根據大數定理我們生成符合 p(x)p(x) 分布的 N 個隨機變量然后用下式來近似計算函數的期望等效于上式積分值。

?但是如過分布函數 p(x)過于復雜我們是很難生成符合條件的隨機變量的計算機只能生成符合諸如均勻分布正態分布等簡單分布的隨機變量對于復雜分布我們需要用一些特殊的技巧才能生成其分布這里就用到馬爾可夫鏈了。

2 有限維空間馬爾可夫鏈

? Example: Predicting the weather with a finite state-space Markov chain

? 假如加州大學伯克利分校的天氣有三種晴天霧天和雨天這用來模擬狀態空間的三個離散的值。這里的天氣狀態十分穩定所以伯克利的天氣預報員可以輕松地根據當前的天氣來預測下周的天氣按照下面轉移概率矩陣進行計算

2.1 直接采用MC公式進行求解

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33

clear?all close?all clc P?=?[0.8??0.15?0.05;???%?SUNNY ???0.4??0.5??0.1;????%?FOGGY ???0.1??0.3??0.6];????%?RAINY nWeeks?=?25; ??? %?初始化狀態為RAINY概率為1 X(1,:)?=?[0?0?1]; ??? %?RUN?MARKOV?CHAIN for?iB?=?2:nWeeks ????X(iB,:)?=?X(iB-1,:)*P;?%?TRANSITION end ??? %?DISPLAY figure;?hold?on h(1)?=?plot(1:nWeeks,X(:,1),'r','Linewidth',1);???%?SUNNY概率 h(2)?=?plot(1:nWeeks,X(:,2),'k','Linewidth',1);???%?FOGGY概率 h(3)?=?plot(1:nWeeks,X(:,3),'b','Linewidth',1);???%?RAINY概率 h(4)?=?plot([15?15],[0?1],'g--','Linewidth',1);??%暫態階段熱機階段 hold?off legend(h,?{'Sunny','Foggy','Rainy','Burn?In'}); xlabel('Week') ylabel('p(Weather)') xlim([1,nWeeks]); ylim([0?1]); ??? %?PREDICTIONS fprintf('\np(weather)?in?1?week?-->'),?disp(X(2,:)) fprintf('\np(weather)?in?24?weeks?-->'),?disp(X(24,:)) fprintf('\np(weather)?in?25?weeks?-->'),?disp(X(25,:))

輸出結果：

p(weather) in 1 week --> ? ?0.1000 ? ?0.3000 ? ?0.6000

p(weather) in 24 weeks --> ? ?0.5965 ? ?0.2632 ? ?0.1404

p(weather) in 25 weeks --> ? ?0.5965 ? ?0.2632 ? ?0.1404

2.1 采用MC抽樣進行求解

? 上面的代碼展示了經過多次迭代后馬爾可夫鏈的分布趨于一個穩態分布。我們將利用馬爾可夫鏈的這一性質進行抽樣，從而得到符合這個穩態分布的樣例。

? 之所以討論MCMC方法，就是為了得到符合某種分布的樣例，通常，我們不容易對這種分布進行直接抽樣（知道了分布函數也不行），我們計算機只能對均勻分布、正態分布等一系列簡單分布進行抽樣，對復雜分布是無法抽樣的。所以，我們的思路是：構造這樣的馬爾可夫鏈，使其穩態分布符合我們要抽樣的函數分布，然后沿著馬爾可夫鏈進行抽樣，從而得到符合穩態分布的樣例，也就是間接地得到了符合目標分布函數的樣例。下面我們繼續利用上例的轉移矩陣，實地抽樣，看看得到的樣例分布是否符合其穩態分布。

方法一：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

clear?all close?all clc P?=?[0.8?0.15?0.05;???%?SUNNY ???0.4?0.5??0.1;???%?FOGGY ???0.1?0.3??0.6];??%?RAINY nWeeks?=?25000;?%?為了精確，我們進行25000次抽樣 X?=?zeros(nWeeks,3);?%?x?用來存儲該周天氣 p?=?zeros(nWeeks,3);?%?p?用來存儲到當前周為止，每種天氣出現的頻率 X(1,:)?=?[0,0,1];?%?第一周天氣為雨天 for?i=2:nWeeks ????px?=?X(i-1,:)*P;?%?根據前一周的天氣來計算這一周三種天氣出現的概率 ????%?我們把三種天氣的累計頻率計算出來，比如px=[0.3,0.5,0.2] ????%?那么?cump?=?[0.3,0.8,1]，然后生成0到1之間符合均勻分布的隨機數， ????%?這個隨機數小于0.3的概率是0.3，大于0.3小于0.8的概率是0.5，大于0.8 ????%?小于1的概率是0.2，這樣就可以模擬隨機抽樣的過程了 ????cump?=?cumsum(px); ????%?生成一個在0到1之間均勻分布的隨機數 ????rnd?=?rand; ????%?看這個隨機數落到哪個區間：如果小于第一個數，那么就是晴天，如果在第一個數 ????%?到第二數個之間，就是霧天，如果在第二個數到第三個數之間，就是雨天 ????if(rnd<=cump(1)) ????????X(i,:)?=?[1,0,0]; ????elseif(rnd>cump(1)&&rnd<=cump(2)) ????????X(i,:)?=?[0,1,0]; ????else ????????X(i,:)?=?[0,0,1]; ????end end cumx?=?cumsum(X,1);?%?每一行的元素代表到這一周為止，出現的晴天、霧天、雨天的次數 %?每一行各種天氣的頻數除以當前行對應周數 p?=?cumx./repmat(transpose(1:nWeeks),1,3); hold?on; h(1)?=?plot(1:nWeeks,p(:,1),'r','Linewidth',1); h(2)?=?plot(1:nWeeks,p(:,2),'k','Linewidth',1); h(3)?=?plot(1:nWeeks,p(:,3),'b','Linewidth',1); hold?off legend(h,?{'Sunny','Foggy','Rainy'}); xlabel('Week') ylabel('p(Weather)') xlim([1,nWeeks]); ylim([0?1]); %?PREDICTIONS fprintf('\np(weather)?in?1?week?-->'),?disp(p(2,:)) fprintf('\np(weather)?in?24999?weeks?-->'),?disp(p(nWeeks-1,:)) fprintf('\np(weather)?in?25000?weeks?-->'),?disp(p(nWeeks,:))

輸出結果：

p(weather) in 1 week --> ? ? ? ? 0 ? ?0.5000 ? ?0.5000

p(weather) in 24999 weeks --> ? ?0.6002 ? ?0.2569 ? ?0.1429

p(weather) in 25000 weeks --> ? ?0.6002 ? ?0.2569 ? ?0.1429

方法二：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

clear?all close?all clc P?=?[.8?.15?.05;??%?SUNNY ?????.4?.5??.1;???%?FOGGY ?????.1?.3??.6];??%?RAINY nWeeks?=?25000;?%?為了精確，我們進行25000次抽樣 X?=?ones(1,nWeeks);?%?x?用來存儲該周天氣 p?=?zeros(nWeeks,3);?%?p?用來存儲到當前周為止，每種天氣出現的頻率 X(1)?=?3;?%?第一周天氣為雨天 C=cumsum(P,2); for?i=2:nWeeks???? ????%?生成一個在0到1之間均勻分布的隨機數 ????rnd?=?rand; ????%?看這個隨機數落到哪個區間：如果小于第一個數，那么就是晴天，如果在第一個數 ????%?到第二數個之間，就是霧天，如果在第二個數到第三個數之間，就是雨天 ????if(rnd<=C(X(i-1),1)) ????????X(i)?=?1; ????elseif(rnd>C(X(i-1),1)&&rnd<=C(X(i-1),2)) ????????X(i)?=?2; ????else ????????X(i)?=?3; ????end end X=ind2vec(X); X=X'; cumx?=?cumsum(X,1);?%?每一行的元素代表到這一周為止，出現的晴天、霧天、雨天的次數 %?每一行各種天氣的頻數除以當前行對應周數 p?=?cumx./repmat(transpose(1:nWeeks),1,3); hold?on; h(1)?=?plot(1:nWeeks,p(:,1),'r','Linewidth',1); h(2)?=?plot(1:nWeeks,p(:,2),'k','Linewidth',1); h(3)?=?plot(1:nWeeks,p(:,3),'b','Linewidth',1); hold?off legend(h,?{'Sunny','Foggy','Rainy'}); xlabel('Week') ylabel('p(Weather)') xlim([1,nWeeks]); ylim([0?1]); %?PREDICTIONS fprintf('\np(weather)?in?1?week?-->'),?disp(p(2,:)) fprintf('\np(weather)?in?24999?weeks?-->'),?disp(p(nWeeks-1,:)) fprintf('\np(weather)?in?25000?weeks?-->'),?disp(p(nWeeks,:))

輸出結果：

p(weather) in 1 week --> ? (1,3) ? ? ? ?1

p(weather) in 24999 weeks --> ? (1,1) ? ? ? 0.6071

? ?(1,2) ? ? ? 0.2536

? ?(1,3) ? ? ? 0.1393

p(weather) in 25000 weeks --> ? (1,1) ? ? ? 0.6071

? ?(1,2) ? ? ? 0.2536

? ?(1,3) ? ? ? 0.1393

3 連續狀態空間馬爾可夫鏈

? ?我們可以通過連續狀態空間的馬爾可夫鏈的穩態分布，來對一個連續概率密度函數取樣（通常這個概率密度函數都不容易直接取樣，原因是函數太過復雜）：我們運行充分次數的馬爾可夫鏈從而讓它達到穩態分布，然后保留馬爾可夫鏈到達穩態分布后抽取的樣例。

? ?下面的例子里，我們定義一個連續狀態空間馬爾可夫鏈。轉移算子是正態分布，其方差為1，均值為當前狀態與0點的距離的一半，初始狀態的分布是均值為0方差為1的標準正態分布。

? ?為了保證鏈條已經移動得離初始狀態足夠遠，達到了穩態分布，我們將扔掉馬爾可夫鏈最初的50個抽樣（也就是 burn in 狀態）。我們同時運行多個鏈，從而更加致密地得到穩態分布抽樣。這里我們選擇同時運行5個馬爾可夫鏈。

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49

%?EXAMPLE?OF?CONTINUOUS?STATE-SPACE?MARKOV?CHAIN clear?all close?all clc?? %?INITIALIZE randn('seed',12345) nBurnin?=?50;?%?#?BURNIN nChains?=?5;??%?#?MARKOV?CHAINS ??? %?DEFINE?TRANSITION?OPERATOR %?定義內聯函數,用來產生符合正態分布的隨機變量normrnd(均值,方差,行,列) P?=?inline('normrnd(.5*x,1,1,nChains)','x','nChains'); nTransitions?=?1000; %?將所有狀態初始化為0，抽樣時對這個矩陣進行更新 x?=?zeros(nTransitions,nChains); %初始狀態，由于是簡單的正態分布，所以可以用函數得到 x(1,:)?=?randn(1,nChains); ??? %?RUN?THE?CHAINS for?iT?=?2:nTransitions ????%后一個狀態由前一個狀態決定，這里源代碼等號右邊是?P(x(iT-1),nChains)不準確， ????%因為5個馬爾可夫鏈上，每個抽樣都是由上一個抽樣決定的，而不是由第一個鏈的上一個狀態決定 ????x(iT,:)?=?P(x(iT-1,:),nChains); end ??? %?DISPLAY?BURNIN figure subplot(221);?plot(x(1:100,:));?hold?on; minn?=?min(x(:)); maxx?=?max(x(:)); l?=?line([nBurnin?nBurnin],[minn?maxx],'color','k','Linewidth',2); ylim([minn?maxx]) legend(l,'~Burn-in','Location','SouthEast') title('First?100?Samples');?hold?off ??? %?DISPLAY?ENTIRE?MARKOV?CHAIN subplot(223);?plot(x);hold?on; l?=?line([nBurnin?nBurnin],[minn?maxx],'color','k','Linewidth',2); legend(l,'~Burn-in','Location','SouthEast') title('Entire?Chain'); ??? %?DISPLAY?SAMPLES?FROM?STATIONARY?DISTRIBUTION samples?=?x(nBurnin+1:end,:); subplot(122); %?這里?hist?將樣例的區間分成100份，返回bins為各個區間的中點，counts是落到各個區間的數目 [counts,bins]?=?hist(samples(:),100);?colormap?hot b?=?bar(bins,counts); legend(b,sprintf('Markov?Chain\nSamples')); title(['\mu=',num2str(mean(samples(:))),'?\sigma=',num2str(var(samples(:)))])

輸出結果：

?左上方的輸出，我們可以看到1000次轉移中的前100次，由5個馬爾可夫鏈同時運行得到。burn in 的切割點用黑線顯示。左下方可以看到馬爾可夫鏈的整體轉移過程。右邊的狀態分布直方圖，表示 burn in 之后的抽樣的分布情況，我們可以發現大部分抽樣值都集中在-2到2之間，這個是一個正態分布，均值為0，方差為1.3。

4?結語

? 上面的例子中，我們能夠根據 burn in 期之后的抽樣推斷出馬爾可夫鏈的穩態分布。然而，為了能夠利用馬爾可夫鏈來對目標分布函數進行抽樣，我們需要設計轉移算子，使得經過多次迭代后，鏈最終達到與目標函數吻合的分布。這就是諸多 MCMC 方法，例如 Metropolis sampler，Metropolis-Hastings sampler，以及 Gibbs sampler 所要實現的目標。

參考：

https://theclevermachine.wordpress.com/2012/09/24/a-brief-introduction-to-markov-chains/

https://theclevermachine.wordpress.com/2012/10/05/mcmc-the-metropolis-sampler/

https://theclevermachine.wordpress.com/2012/10/20/mcmc-the-metropolis-hastings-sampler/

http://www.cnblogs.com/yinxiangnan-charles/p/5018876.html

? ? ?本文轉自stock0991 51CTO博客，原文鏈接：http://blog.51cto.com/qing0991/1793526，如需轉載請自行聯系原作者

總結

以上是生活随笔為你收集整理的Markovs Chains采样的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。