馬爾科夫鏈蒙特卡洛采樣（MCMC）入門 之一

 

目錄

馬爾科夫鏈蒙特卡洛采樣（MCMC）入門 之一

馬爾科夫鏈蒙特卡洛采樣（MCMC）入門 之二

 

1、從隨機變量分布中采樣

研究人員提出的概率模型對於分析方法來說往往過於復雜。越來越多的研究人員依賴數學計算的方法處理復雜的概率模型，研究者通過使用計算的方法，擺脫一些分析技術所需要的不切實際的假設。（如，正態和獨立）

大多數近似方法的關鍵是在於從分布中采樣的能力，我們需要通過采樣來預測特定的模型在某些情況下的行為，並為潛在的變量（參數）找到合適的值以及將模型應用到實驗數據中，大多數采樣方法都是將復雜的分布中抽樣的問題轉化到簡單子問題的采樣分布中。

本章，我們解釋兩種采樣方法：逆變換方法（the inverse transformation method）和拒絕采樣（rejection sampling）。這些方法主要適用於單變量的情況，用於處理輸出單變量的問題。在下一章，我們討論馬爾科夫鏈蒙特卡洛方法（Markov chain Monte Carlo），該方法可以有效的用於多元變量分布采樣。

 

1.1 從標准分布中采樣

有一些分布被經常用到，這些分布被MATLAB作為標准分布實現。在MATLAB統計工具箱（Matlab Statistics Toolbox supports）實現了一系列概率分布。使用MATLAB工具箱可以很方便的計算這些分布的概率密度、累積密度、並從這些分布中取樣隨機值。表1.1列舉了一些MATLAB工具箱中的標准分布。在MATLAB文檔中列舉了更多的分布，這些分布可以用MATLAB模擬。利用在線資源，通常很容易能找到對其他常見分布的支持。

 

  

 

 為了說明如何使用這些函數，Listing 1.1展示了正態分布N(μ,σ)可視化的MATLAB代碼，其中μ=100,σ=15。

舉例

舉個例子，可以想象一下用該正態分布表示觀察到的人群的IQ系數變化。該代碼顯示了了如何展示概率密度和累積密度。它還展示了如何從該分布中抽取隨機值以及如何使用hist函數可視化這些隨機樣本。代碼的輸出結果如 圖1.1 所示。

 %% Explore the Normal distribution N( mu , sigma )
 mu = 100; % the mean
 sigma = 15; % the standard deviation
 xmin = 70; % minimum x value for pdf and cdf plot
 xmax = 130; % maximum x value for pdf and cdf plot
 n = 100; % number of points on pdf and cdf plot
 k = 10000; % number of random draws for histogram

 % create a set of values ranging from xmin to xmax
 x = linspace( xmin , xmax , n );
 p = normpdf( x , mu , sigma ); % calculate the pdf
 c = normcdf( x , mu , sigma ); % calculate the cdf

 figure( 1 ); clf; % create a new figure and clear the contents

 subplot( 1,3,1 );
 plot( x , p , 'k?' );
 xlabel( 'x' ); ylabel( 'pdf' );
 title( 'Probability Density Function' );

 subplot( 1,3,2 );
 plot( x , c , 'k?' );
 xlabel( 'x' ); ylabel( 'cdf' );
 title( 'Cumulative Density Function' );

 % draw k random numbers from a N( mu , sigma ) distribution
 y = normrnd( mu , sigma , k , 1 );

 subplot( 1,3,3 );
 hist( y , 20 );
 xlabel( 'x' ); ylabel( 'frequency' );
 title( 'Histogram of random values' );

 

Listing 1.1: Matlab code to visualize Normal distribution.

 

1.2 從非標准分布中采樣

我們希望MATLAB工具也支持從非標准分布中采樣，這種情況在建模過程中經常出現，因為研究人員可以提出一種新的噪聲過程或已存在分布的組合方式。復雜采樣問題的計算方法通常依賴於我們已經知道如何有效地進行采樣的分布。這些從簡單分布中采樣的隨機值可以被轉換成目標分布需要的值。事實上，這一節我們討論的一些技術是MATLAB的內部分布，如正態分布和指數分布。

1.2.1 用離散變量進行逆變換采樣（Inverse transform sampling）

逆變換采樣（也被成為逆變換方法）即給定累積分布函數的逆，可從任意概率分布中生成隨機數。這個方法是對均勻分布的隨機數字進行采樣（在0到1之間）然后使用逆累積分布函數轉換這些值。該過程的簡單之處就在於，潛在的采樣僅僅依賴對統一的參數進行偏移和變換。該過程可以用於采樣很多不同種類的分布，事實上，MATLAB實現很多隨機變量生成方法也是基於該方法的。

在離散分布中，我們知道每個輸出結果的概率。這種情況下，逆變換方法就需要一個簡單的查找表。

給定一個非標准的離散分布的例子，我們使用一些實驗數據來研究人類如何能產生一致的隨機數（如Treisman and Faulkner, 1987）。在這些實驗中，被測試者會產生大量的隨機數字（0,…,9）。研究人員根據每個隨機數字的相對頻率進行制表。你可能會懷疑實驗對象不會總是產生均勻分布。表1.2 展示了一些典型的數據，其中可以看出一些比較低的和高的數字容易被忽視，而一些特殊數字（如數字4）占過高的比例。由於某種原因，數字0和9從來沒有被產生。在任何情況下，這些數字都是相當典型的，而且證明了人類不能很好地產生均勻分布的隨機數字。

 

 

 假設我們想要模擬這個過程，並根據表1.2 中的概率編寫一個算法采樣數字。因此，程序應該用概率0.2生成數字4，根據概率0.175生成數字5等。例如，Listing1.2中的代碼使用MATLAB內置的函數randsample來實現這個過程。代碼的輸出結果如圖1.2所示。

  % Simulate the distribution observed in the
 % human random digit generation task

 % probabilities for each digit
 theta = [0.000; ... % digit 0
 0.100; ... % digit 1
 0.090; ... % digit 2
 0.095; ... % digit 3
 0.200; ... % digit 4
 0.175; ... % digit 5
 0.190; ... % digit 6
 0.050; ... % digit 7
 0.100; ... % digit 8
 0.000 ] ... % digit 9

 % fix the random number generator
 seed = 1; rand( 'state' , seed );

 % let's say we draw K random values
 K = 10000;
 digitset = 0:9;
 Y = randsample(digitset,K,true,theta);
 % create a new figure
 figure( 1 ); clf;

 % Show the histogram of the simulated draws
 counts = hist( Y , digitset );
 bar( digitset , counts , 'k' );
 xlim( [ 0.5 9.5 ] );
 xlabel( 'Digit' );
 ylabel( 'Frequency' );
 title( 'Distribution of simulated draws of human digit generator' );

Listing 1.2: Matlab code to simulate sampling of random digits.

 

 

 Figure 1.2.1: Illustration of the inverse transform procedure for generating discrete random variables. Note that we plot the cumulative probabilities for each outcome. If we sample a uniform random number of U = 0.8, then this yields a random value of X = 6

如果不采用內置的函數如 randsample 和 mnrnd，而是通過逆變換方法來實現底層的采樣算法對我們更有幫助。我們首先需要計算累計概率分布，換句話說，我們需要知道我們觀察到的結果等於或小於某一特定值的概率。如果F(X)表示累計函數，我們需要計算F(X=x)=p(X≤x)。對於離散分布，計算這個值可以通過簡單的求和。我們的例子的累計概率在表1.2的最后一列中給出。在逆變換算法中，該想法是采樣隨機偏差（0和1之間的隨機數）並將隨機數與表中的累計概率比較。隨機偏差的第一個結果小於（或等於）相關的累計概率與抽樣結果相對應。圖1.2.1展示了一個例子，其中隨機偏差U=0.8，導致采樣結果X=6。這個重復采樣隨機偏差的過程，並與累積分布相比較，就會形成離散變量的逆變換方法的基礎。注意我們應用了一個逆函數，因為做的是逆表的查找。

1.2.2 連續變量的逆變換采樣

逆變換樣方法也可以用於連續分布。一般地，該方法目的是獲得均勻的隨機偏差並且將逆函數應用在隨機偏差的累積分布中。下面，令F(X)是目標變量X的累積密度函數（cumulative density function，CDF），F⁻¹(X)是該函數的逆。我們希望獲得隨機變量X的值，可以通過如下步驟獲取：

1. 獲得均勻分布 U∼Uniform(0,1)

 設

 重復上述采樣過程

可通過一個簡單的例子解釋上述方法。假設我們要從指數分布（exponential distribution）中采樣隨機數。當λ>0時，累積密度函數是F(x∣λ)=1−exp(−x/λ)。用一些簡單的代數方法，就可以求出這個函數的逆 。從而引出了下面的從Exponental(λ)分布中采樣隨機數的步驟：

1. 獲得均勻分布

2. 
   

    設

    重復上述采樣過程

    

 

1.2.3 拒絕采樣

在很多情況下，逆變換采樣方法是不適用的，因為很難計算其累積分布和它的逆。在這種情況下，有其他的可供選擇的采樣方法，如拒絕采樣，以及其他使用馬爾科夫鏈蒙特卡洛的方法，將在下一章進行討論。拒絕采樣的主要優勢在於它不需要“burn-in”過程。和其他方法相反，所有的樣本在采樣過程中獲得，並且可直接作為目標分布的樣本。

拒絕采樣的一般概念可用圖1.5進行解釋。假設我們想在以(0,0)為圓心以1為半徑的圓形內均勻畫點。乍一看，似乎很難均勻地在這個圓內直接進行采樣。但是我們可以應用拒絕采樣的方法進行采樣，首先確定圓外圍的正方形，在正方形中采樣(x,y)的值，剔除滿足式子的所有的取值。注意在這個過程中我們使用了一個簡單的建議分布（q），如均勻分布，作為從更復雜的分布中采樣的基礎。

拒絕采樣允許我們從難以采樣的分布中生成樣本，在這些難以采樣的分布中我們可以計算任何特定樣本的概率。換句話說，假定我們有一個分布p(θ)，並且難以直接從該分布中采樣，但是我們可以計算其特定值的概率密度p(θ)。

 

 

第一件要做的就是建議分布（proposal distribution）。建議分布就是指一個簡單的分布，記為q(θ)，該分布是可以直接進行采樣的。主要思路是，計算同時滿足建議分布和目標分布采樣的概率，然后拒絕相對於建議分布中那些不太可能出現在目標分布下的樣本。

 

 

圖1.6說明了這個過程。我們首先需要找到一個常量c使之對所有的可能的樣本θ滿足cq(θ)≥p(θ)。建議函數q(θ)乘以常量c被稱為比較分布，並且使其總是大於目標分布。要確定常數cc並不容易，我們先假定可以通過微積分確定該常數。

我們首先從均勻分布[0,cq(θ)]中獲取一個數u，換句話說，這是直線段從0到cq(θ)的某個點以θ為建議的比較分布。如果u>p(θ)，我們拒絕這個建議分布采樣得到的值，否則，接受之。如果接受了某個建議，則采樣值θ就是從目標分布p(θ)中獲得的。上述的采樣過程步驟如下：

 

步驟

1. 選擇一個容易采樣的分布q(θ)

2. 確定常量c，使對所有θ的有cq(θ)≥p(θ)

3. 從建議分布q(θ)中采樣一個建議樣本θ

4. 從區間[0,cq(θ)]采樣一個數u

5. 如果u>p(θ)則拒絕，否則接受

6. 重復步驟3,4,5，直到達到要求的樣本數量；每個接受的樣本都是從p(θ)中獲得的

這種算法有效的關鍵就是需要有盡可能多的樣本被接受，這取決於建議分布q(θ)的選擇。建議分布不同於目標分布，它會導致很多拒絕接受的樣本，這將減慢采樣的速度。

 

目錄

馬爾科夫鏈蒙特卡洛采樣（MCMC）入門 之一

馬爾科夫鏈蒙特卡洛采樣（MCMC）入門 之二

 

 

 代碼鏈接： http://pan.baidu.com/s/1qXIWzJu

MCMC matlab tutorial 電子版本：http://pan.baidu.com/s/1i48vyr7

Reference

Mark Steyver. Computational Statistics with Matlab. 2011

 

 

From:

https://mp.weixin.qq.com/s?__biz=MzU2OTA0NzE2NA==&mid=2247483828&idx=1&sn=24dc135cf0d45a1dd276c0247c3e99d4&chksm=fc85e0a7cbf269b1b368465431e8050f1591f8d5611a923de92d999aed47abd8ce5cb0367142&scene=21#wechat_redirect