MCMC(馬爾可夫鏈蒙特卡洛方法)：the Gibbs Sampler(吉布斯取樣)

        基於兩個標準，吉布斯取樣使用某些類別的問題。給定一個目標分佈p(x)，其中，第一個標準是以其他所有變數聯合起來的聯合分佈為條件的每一個變數的條件分佈有解析(數學)表示式。在形式上，如果目標分佈p(x)是D維的，我們必須有D個獨立的表示式：

        這些表示式的每一個都定義了在知道其他j(j≠i)維的數值的情況下第i維的概率。具有每一個變數的條件分佈代表我們不需要像Metropolis Hastings演算法需要提議分佈或者接受/拒絕標準。因此，當其他變數固定的時候，我們可以簡單的從每一個條件中去取樣。第二個標準就是我們必須能夠從每一個條件分佈中去取樣。如果我們想要去實現一個演算法，這個附加條件是非常明顯的。

   吉布斯取樣的工作方法和分成分Metropolis Hastings演算法很像，除了取締借鑑每一個維度的提議分佈，然後對於接受或者拒絕提議取樣，我們採用簡單地依據變數對應的條件分佈去選取此維度的值。我們會接受所有選取的值。類似分成分Metropolis Hastings演算法，我們依次通過每一個變數，在其它變數固定的時候對它取樣。吉布斯取樣的步驟大致如下：

2.生成初始狀態

3.重複直到t=M
{
對於每一個維度i=1...D
從中得到 

}
       為了對吉布斯取樣有更好的理解，我們下面來實現一下吉布斯取樣，去解決與前面提到過的同樣的多元變數取樣問題。

例子：從二元正態分佈中取樣Example: Sampling from a bivariate a Normal distribution

①均值是

②方差是

     為了使用吉布斯取樣從這個分佈中取樣，我們需要有變數/維度x1和x2的條件分佈：

      是第二個維度的前一個狀態，是從中得到的第一個維度的狀態。有差異的原因是更新x1和x2用的是(t-1)和t時刻的狀態，在上一節中的演算法大綱第三步可以看出來。第t次迭代，我們首先以變數x2的最近狀態即第(t-1)次迭代結果為條件，為x1取樣一種新狀態。然後再以第t次迭代得到的x1的最新狀態為條件取樣得到變數x2。
經過一些數學推導(這裡先跳過，下面會有詳細的過程)，我們發現目標正態分佈的兩個條件分佈是：

        每一個都是單變數的正態分佈，其中均值依賴條件變數的最近狀態的值，方差依賴兩個變數之間的目標方差。
       使用上述描述的變數x1和x2的條件概率，我們下面採用matlab實現吉布斯取樣，輸出的取樣如下：

        觀察上表，發現每一次迭代中吉布斯取樣中的馬爾可夫鏈是如何第一次沿著x1方向前進一步，然後沿著x2的方向前進。這展示了吉布斯取樣以分成分方法一次單獨為每一個變數取樣。

總結Wrapping Up

        吉布斯取樣是為複雜多元概率分佈取樣的一個受歡迎的MCMC方法。然而，吉布斯取樣不能用於一般的抽樣問題。對於許多目標分佈，很難或者不可能去獲取到所有需要的條件分佈的近似表達。在其它情況下，對於所有條件的解析表示式或許存在，但是它或許很難從任意的或者全部的條件分佈去取樣(在這種情況下，使用單變數( univariate sampling methods)取樣比如拒絕抽樣(rejection sampling)和Metropolis型別的MCMC技術去逼近每一個條件的樣本是比較普遍的)。吉布斯取樣是非常受歡迎的貝葉斯方法，模型經常以這種方式設計：所有模型變數的條件表示式非常容易獲取，並且採用一種能夠被高效取樣的眾所周知的形式。
吉布斯取樣，就想很多MCMC方法，有“慢混合(slow mixing)”的問題。慢混合的發生是在潛在的馬爾可夫鏈需要很長時間去充分探索出x的值，從而給出一個更好的p(x)的表徵(characterization)。慢混合是因為一些因素包括馬爾可夫鏈的“隨機走動(random walk)”特性，並且馬爾可夫鏈有“卡住”的趨勢，因為僅僅取樣了x的一個單獨區域，這個區域在p(x)下有很高的概率。這種反應對於多模式(multiple modes)或者重尾(heavy tails)中的分佈進行取樣效果不好，比如混合蒙特卡洛已被髮展成一個包含附加動力學(incorporate additional dynamics)的能提高馬爾可夫鏈路徑效率的方法。將來會討論混合蒙特卡洛方法

matlab程式碼

      實現的應該是給定了一個均值和方差，以及初始的一個樣本點，然後取樣出5000個符合分佈的點

%https://victorfang.wordpress.com/2014/04/29/mcmc-the-gibbs-sampler-simple-example-w-matlab-code/
%seed 用來控制 rand 和 randn 
% 如果沒有設定seed，每次執行rand或randn產生的隨機數都是不一樣的
% 用了seed，比如設定rand('seed',0);，那麼每次執行rand產生的隨機數是一樣的，這樣對除錯程式很有幫助
rand('seed' ,12345);

nSamples = 5000;
 
mu = [0 0]; % TARGET MEAN目標均值
rho(1) = 0.8; % rho_21目標方差
rho(2) = 0.8; % rho_12目標方差
 
% INITIALIZE THE GIBBS SAMPLER
propSigma = 1; % PROPOSAL VARIANCE
minn = [-3 -3];
maxx = [3 3];
 
% INITIALIZE SAMPLES
x = zeros(nSamples,2);
x(1,1) = unifrnd(minn(1), maxx(1));%unifrnd生成連續均勻分佈的隨機數
x(1,2) = unifrnd(minn(2), maxx(2));
 
dims = 1:2; % INDEX INTO EACH DIMENSION
 
% RUN GIBBS SAMPLER
t = 1;
while t < nSamples%總共取樣出5000個取樣點
    t = t + 1;
    T = [t-1,t];
    for iD = 1:2 % LOOP OVER DIMENSIONS總共兩維，註釋先討論第一維
        % UPDATE SAMPLES
        nIx = dims~=iD; % *NOT* THE CURRENT DIMENSION找到另外一維nIx=[0 1]logical型別
        % CONDITIONAL MEAN
        muCond = mu(iD) + rho(iD)*(x(T(iD),nIx)-mu(nIx));%計算均值=表示式μ(1)+ρ(1)*(x(n,2)-μ(2))其中x(n,2)代表樣本第n個數據的第二維
        % CONDITIONAL VARIANCE
        varCond = sqrt(1-rho(iD)^2);%計算方差
        % DRAW FROM CONDITIONAL
        x(t,iD) = normrnd(muCond,varCond);%正態分佈隨機函式，計算得到當前第t個數據的第1維資料value
    end
end
 
% DISPLAY SAMPLING DYNAMICS
figure;
h1 = scatter(x(:,1),x(:,2),'r.');%scatter描繪散點圖，x為橫座標，y為縱座標
 
% CONDITIONAL STEPS/SAMPLES
hold on;%畫出前五十個取樣點
for t = 1:50
    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-');
    h2 = 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],{'Samples','1st 50 Samples','x(t=0)'},'Location','Northwest')
hold off;
xlabel('x_1');
ylabel('x_2');
axis square