一.概述

????????通過概率統計模擬來進行數值計算的方法統稱為蒙特卡羅(Monte Carlo)方法，而MCMC方法稱為馬爾科夫鏈蒙特卡洛(Markov Chain Monte Carlo)方法。顯然，MCMC法為MC法的一種特例。MCMC法是利用馬爾可夫鏈的細致平衡條件進行采樣，再通過所采樣的樣本進行數值計算的一種方法。

二.實例

? ? ? ? 下面將介紹三個利用MC方法的實例，分別為

? ? ? （1）對圓周率進行數值計算

? ? ? （2）求解定積分的數值解

? ? ? （3）由重要采樣法求解定積分的數值解

1.對圓周率進行數值計算

????????設某二維均勻分布U的正概率密度區域是邊長為2的正方形，如下圖所示

該正方形的內切圓已在圖中標紅，假設向該正方形內均勻的撒一把豆子，理想情況下，豆子落入紅圓內的概率為

其中m為總豆子數，n為落入紅圓內的豆子數。顯然，本次撒豆子的行為服從二維均勻分布U，因此，又可以得到

從而有

????????不難發現，通過撒豆子(生成均勻分布隨機數)，便可計算出的數值，這個方法便是MC方法的一個簡單例子，有大數定律可知，撒的豆子越多，所求得的值便越精確，下面將利用MATLAB編程對上述過程進形模擬，程序如下

%在正方形內生成均勻分布隨機數，即撒豆子 N=1e4; x=unifrnd(-1,1,1,N); y=unifrnd(-1,1,1,N); z=x(x.^2+y.^2<=1);%找出落入內切圓內以及圓上的豆子 p=4*length(z)/length(x);%計算pi值 err=(p-pi)/pi*100;%計算相對誤差 %繪圖 hold on for i=1:length(x)if x(i)^2+y(i)^2<=1plot(x(i),y(i),'.r')elseplot(x(i),y(i),'.b')end end title('樣本點分布圖'); axis square; hold off fprintf('pi值為%1.4f，相對誤差為%1.6f\n',p,err);%顯示結果

【注】 為了計算速度，該程序只生成了1e4個樣本點，樣本點越多結果越精確。

運行結果如下

?

?2.求解定積分的數值解

? ? ? ? 由實例1不難發現MC方法的一個基本思想，即利用樣本量足夠大的古典概型與幾何概型建立等式，進而求得等式中未知參數的數值解，因幾何概型涉及到面積之比，并且定積分的數值解也是一個面積，故利用上述思想亦能對定積分進行求解。以如下定積分為例

顯然，被積函數在[0,]上的最大值為1，因此撒豆子的區域為

如下圖所示，設m為豆子總數，n為落入紅色曲線(被積函數)以下豆子的數量

則待求積分值為

下面將利用MATLAB編程對上述過程進形模擬，程序如下

%在矩形目標區域內撒豆子 x=unifrnd(0,pi,1,1e6); y=rand(1,1e6); %定義被積函數并找出落入被積函數曲線以下的豆子 f=@(x)exp(-x.^2); lg=y<=f(x); xx=x(lg);xxf=x(~lg); yy=y(lg);yyf=y(~lg); inti=pi*length(xx)/length(x); %求解積分值 zjb=integral(f,0,pi); %求解真值 err=(inti-zjb)/zjb*100; %計算相對誤差 %繪制樣本點分布圖 hold on scatter(xx,yy,'.r') scatter(xxf,yyf,'.b') title('樣本點分布圖'); grid on;box on; hold off fprintf('積分值為%1.4f，相對誤差為%1.6f\n',inti,err); %顯示結果

?運行結果如下

?

3.由重要采樣法求解定積分的數值解

? ? ? ? 實例1,2介紹了MC方法的一個基本思想，接下來將通過該實例介紹另一種思想——通過重要采樣法求解定積分的數值解

? ? ? ? ?重要采樣(Importance sampling)法概述：

?因此，可由如下過程計算函數f(x)的積分值

下面將通過MATLAB利用上述方法計算定積分?

MATLAB程序如下

%定義q(x） syms x; q1=3*x.^2; q=@(x)double(subs(q1,x)); %定義q(x）的原函數Q(x) Q=@(x)double(subs(finverse(int(q1)),x)); %利用Q(x)對q(x)進行采樣，采樣結果為z，采樣數量為2e4 r=rand(1,2e4); z=Q(r); f=@(x)sin(x);%定義被積函數 F=@(x)f(x)./q(x);%定義隨機變量的函數 jg=mean(F(z));%求取均值，即積分值 zjg=integral(f,0,1);%求取真值 err=(jg-zjg)/zjg*100;%求取相對誤差 fprintf('積分值為%1.4f，相對誤差為%1.6f\n',jg,err); %顯示結果

?運行結果如下

三.采樣方法（這里僅是對下述采樣方法的概括，文獻[1]中有對各個采樣方法的詳細說明）

? ? ? ? 通過上述簡單實例不難發現，MC法需要對某個目標分布進行抽樣，上述實例中的目標分布均為常見的一維分布，采樣方法較為簡單，如實例3使用的是概率分布采樣（利用均勻分布U(0,1)對CDF的反函數進行采樣），但在目標分布較為復雜的情況下，可利用馬爾可夫鏈的細致平衡條件進行采樣，基于此種采樣的MC法，稱為MCMC方法。

? ? ? ? 下面將介紹三種采樣方法

? ? ? ? （1）拒絕采樣

? ? ? ? （2）M-H采樣

? ? ? ? （3）Gibbs采樣

?????????其中（2）和（3）為基于馬爾可夫鏈的采樣法。

1.拒絕采樣(reject sampling)

????????通過建議分布(proposed distribution)q(x),q(x)的CDF應方便求得,構造接受α=f(x)/(mq(x)),其中f(x)為目標分布的概率密度函數，m為一個常數，目的是保證f(x)≤mq(x)，由于M-H法是該方法的改良版，故不對該采樣方法進行編程操作，該方法的偽代碼[1]如下

2.M-H采樣(用于對一維分布進行采樣)

????????該方法的思路與拒絕采樣類似，均構造了接受率，但M-H采樣是利用馬爾可夫鏈細致平衡條件(i)P(i,j)=(j)p(j,i)對接受率進行構造，具體過程如下

對于細致平衡條件(i)P(i,j)=(j)P(j,i)，引入隨機矩陣Q，使得

P(j,i)=Q(i,j)(i,j)? ?

從而有

(i)Q(i,j)(i,j)=(j)Q(j,i)(j,i)

則

?(i,j)=?(j)Q(j,i)

該方法將對隨時間前進的隨機矩陣Q采樣，所采樣本被保留下來的概率為，在M-H方法中，為了使采樣效率相對較高(避免接受率長時間過小)，采用如下接受率公式

下面將通過MATLAB對M-H采樣進行編程,采樣分布的核密度為p(x),隨機矩陣為q(x),接受率為.

偽代碼如下

MATLAB代碼如下

x(1)=0;%樣本初始值 N=2e5;%采樣次數 p=@(x)0.3*exp(-0.2*x.^2)+0.7*exp(-0.2*(x-10).^2); %p(x)的核密度，采樣過程中目標分布PDF中的常數會被消去，因此用核密度即可。 %進行M-H采樣 for i=1:N-1xs=normrnd(x(i),10);%生成候選樣本pxs=p(xs);%計算該候選樣本處的核密度值px=p(x(i));%計算上一個已采樣本處的核密度值qxs=normpdf(xs,x(i),10);%計算在上一個已采樣本條件下，該候選樣本處隨機矩陣q的概率密度值qx=normpdf(x(i),xs,10);%計算在該候選樣本條件下，上一個已采樣本處隨機矩陣q的概率密度值alpha=min(pxs*qx/(px*qxs),1);%計算接受率%按接受率判斷是否接受該候選樣本u=rand;if u<=alphax(i+1)=xs;elsex(i+1)=x(i);end end %繪制各類對比圖像，對采樣結果進行分析 xx=-10:0.01:20; fp=p(xx); A=sum(fp)*0.01;%A=integral(p,-inf,inf),p是核概率密度函數，下面fp/A的目的是使積分為1，求出pdf。 hold on ecdf(x); c=linspace(-10,20,1000); for i=1:length(c)F(i)=integral(@(x)p(x)/A,-inf,c(i)); end plot(c,F,'r') legend('ecdf','cdf') title('cdf對比') [ni,bin]=histcounts(x); binx=bin(2)-bin(1); hold off figure hold on bar((bin(1)+bin(2))/2:binx:(bin(end-1)+bin(end))/2,ni/N/binx) plot(xx,fp/A,'r') hold off legend('采樣','目標分布') title('pdf對比1') figure hold on ksdensity(x) plot(xx,fp/A,'r') title('pdf對比2') legend('采樣','目標分布') hold off

?運行結果如下

3.Gibbs采樣（用于對2維及2維以上分布進行采樣）

該方法的思路與M-H類似，但接受率恒為1，具體數學推導可查閱文獻[1]，其偽代碼如下

MATLAB程序如下

%% 該程序將利用Gibbs采樣二維正態分布 %樣本初始化 N=3000; %采樣數量 mu=[0,0]; %目標分布均值 rho=[0.8,0.8]; %目標分布相關系數p21與p12 Sigma=[1,0.8;0.8,1];%目標分布方差與相關系數矩陣 %確定采樣區域 low=[-3,-3];%采樣區域下限 up=[3,3];%采樣區域上限， x=zeros(N,2);%樣本初始化（分配空間）,設第一列為x1，第二列為x2 %設置樣本初始值 x(1,1)=unifrnd(low(1),up(1)); x(1,2)=unifrnd(low(2),up(2)); dims=1:2; %樣本維度 %Gibbs程序 for t=2:NT=[t-1,t];%為了給下面的均值以及方差公式代數而生成矩陣位置索引for i=1:2 %采樣維度dim=dims~=i;%為了給下面的均值以及方差公式代數而生成的邏輯數組mud=mu(i)+rho(i)*(x(T(i),dim)-mu(dim));%計算均值vard=sqrt(1-rho(i)^2);%計算方差x(t,i)=normrnd(mud,vard);%產生正態分布隨機數end end %繪制Gibbs采樣的對比圖像，分析采樣效果 h1=scatter(x(:,1),x(:,2),'.r'); hold on; for t = 1:30h2=plot([x(t,1),x(t+1,1)],[x(t,2),x(t,2)],'k-');plot([x(t+1,1),x(t+1,1)],[x(t,2),x(t+1,2)],'k-');plot(x(t+1,1),x(t+1,2),'ko'); end h3=scatter(x(1,1),x(1,2),'go','Linewidth',3); legend([h1,h2,h3],'采樣點','前30個樣本點的采樣路徑','初始采樣點','Location','Northwest') hold off title('采樣情況圖') xlabel('x1'); ylabel('x2'); figure xd=mvnrnd(mu,Sigma,5000); subplot(1,3,1);scatter(xd(:,1),xd(:,2),'.r'); xlabel('x1');ylabel('x2');title('目標分布樣本'); subplot(1,3,2);scatter(x(:,1),x(:,2),'.b'); xlabel('x1');ylabel('x2');title('Gibbs采樣樣本') subplot(1,3,3);hold on scatter(xd(:,1),xd(:,2),'.r');title('目標分布樣本') scatter(x(:,1),x(:,2),'.b');title('Gibbs采樣樣本') xlabel('x1');ylabel('x2');title('兩種樣本對比'); hold off xx1=low(1):0.01:up(1); xx2=low(1):0.01:up(1); [X,Y]=meshgrid(xx1,xx2); Z=zeros(length(X),length(X)); for i=1:length(X)for j=1:length(X)Z(i,j)=mvnpdf([X(i,j),Y(i,j)],mu,Sigma);end end figure subplot(1,2,1); scatter3(x(:,1),x(:,2),mvnpdf([x(:,1),x(:,2)],mu,Sigma),'.r'); xlabel('x1');ylabel('x2');zlabel('z');title('采樣樣本PDF'); subplot(1,2,2); mesh(X,Y,Z); xlabel('x1');ylabel('x2');zlabel('z');title('目標分布PDF'); figure mesh(X,Y,Z); hold on; scatter3(x(:,1),x(:,2),mvnpdf([x(:,1),x(:,2)],mu,Sigma),'.r'); xlabel('x1');ylabel('x2');zlabel('z');title('兩種樣本對比'); hold off rectangle('Position', [0,0,pi,1]) z=linspace(0,pi); hold on plot(z,f(z),'r'); grid on

運行結果如下

?

?

四.參考文獻

?[1]An Introduction to MCMC for Machine Learning

總結

以上是生活随笔為你收集整理的MCMC算法原理及其实例的全部內容，希望文章能夠幫你解決所遇到的問題。

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