吉布斯取樣是生成馬爾科夫鏈的一種方法，生成的馬爾科夫鏈可以用來做蒙特卡洛模擬，從而求得一個較複雜的多元分佈。

吉布斯取樣的具體做法：假設有一個k維的隨機向量，現想要構造一條有n個樣本的k維向量（n樣本馬爾科夫序列），那麼（隨機）初始化一個k維向量，然後固定這個向量其中的k-1個元素，抽取剩下的那個元素（生成給定後驗的隨機數），這樣迴圈k次，就把整個向量更新了一遍，也就是生成了一個新的樣本，把這個整體重複n次就得到了一條馬爾科夫鏈。

在統計學和統計物理學中，gibbs抽樣是馬爾可夫鏈蒙特卡爾理論（MCMC）中用來獲取一系列近似等於指定多維概率分佈（比如2個或者多個隨即變數的聯合概率分佈）觀察樣本的演算法。 MCMC是用於構建Markov chain隨機概率分佈的抽樣的一類演算法。MCMC有很多演算法，其中比較流行的是Metropolis-Hastings Algorithm，Gibbs Sampling是Metropolis-Hastings Algorithm的一種特殊情況。 Markov chain 是一組事件的集合，在這個集合中，事件是一個接一個發生的，並且下一個事件的發生，只由當前發生的事件決定。用數學符號表示就是： 這裡的    不一定是一個數字，它有可能是一個向量，或者一個矩陣，例如我們比較感興趣的問題裡    =（g, u, b）這裡g表示基因的效應，u表示環境效應，b表示固定效應，假設我們研究的一個群體，g，u，b的聯合分佈用π（a）表示。事實上，我們研究QTL，就是要找到π（a），但是有時候π（a）並不是那麼好找的，特別是我們要估計的a的引數的個數多於研究的個體數的時候。用一般的least square往往效果不是那麼好。 解決方案： 用一種叫Markov chain Monte Carlo(MCMC)的方法產生Markov chain，產生的Markov chain    具有如下性質：當t 很大時，比如10000，那麼    ~ π（a），這樣的話如果我們產生一個markov chain：   ，那麼我們取後面9000個樣本的平均： a_hat = (g_hat,u_hat,b_hat) = ∑    / 9000 (i=1001,1002, … 10000)； 這裡g_hat, u_hat, b_hat 就是基因效應，環境效應，以及固定效應的估計值； MCMC演算法的關鍵是兩個函式： 1）q（    ,    ），這個函式決定怎麼基於    得到    ； 2）α（    ,    ），這個函式決定得到的    是否保留； 目的是使得    的分佈收斂於π（a）[1] 。

過程

編輯 一般來說我們通常不知道π（a），但我們可以得到p（g | u , b）,p(u | g , b), p ( b | g, u )即三個變數的posterior distribution。 Step1: 給g, u, b 賦初始值：（g0,u0,b0）; Step2: 利用p (g | u0, b0) 產生g1; Step3: 利用p (u | g1, b0) 產生u1; Step4: 利用p (b | g1, u1) 產生b1; Step5: 重複step2~step5 這樣我們就可以得到一個markov chain    。 這裡的q（      ）= p（g | u , b）* p(u | g , b)* p ( b | g, u ) 注意：Gibbs取樣的目的是獲得一個樣本，不是計算概率，但可以通過其他方法來統計概率

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

        在之前的部落格中，我們對比了在一個多元概率分佈p(x)中，分別採用分組（block-wise）和分成分（component-wise）實現梅特羅波利斯哈斯廷斯演算法。對於多元變數問題中的MCMC演算法，分成分更新權重比分組更新更有效，因為通過使每一個成分/維度獨立於其他成分/維度，我們將更可能去接受一個提議取樣【注，這個proposed sample應該就是前面部落格裡面提到的轉移提議分佈】。然而，提議取樣仍然可能被拒絕，導致有些多餘的計算，因為他們被拒絕了，計算了但是一直未使用。吉布斯取樣是另外一種比較受歡迎的MCMC取樣技術，提供了避免這種多餘計算的方法。就想分成分實現Metropolis Hastings演算法，吉布斯仍然使用分成分更新。然而，不像Metropolis Hastings取樣，所有的提議取樣將被接受，因此不會有多餘的計算。

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

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

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

1.設定t=0

2.生成初始狀態

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

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

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

       這個例子與前面一樣，從2維的正態分佈使用分組和分成分的Metropolis-Hastings演算法取樣。這裡我們展示使用同樣的目標分佈如何實現吉布斯取樣。重複提示一下，目標函式p(x)是一種規範化形式，表示如下：

①均值是

②方差是

     為了使用吉布斯取樣從這個分佈中取樣，我們需要有變數/維度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個符合分佈的點

[html] view plain copy  print?

1. %https://victorfang.wordpress.com/2014/04/29/mcmc-the-gibbs-sampler-simple-example-w-matlab-code/  
2. %seed 用來控制 rand 和 randn   
3. % 如果沒有設定seed，每次執行rand或randn產生的隨機數都是不一樣的  
4. % 用了seed，比如設定rand('seed',0);，那麼每次執行rand產生的隨機數是一樣的，這樣對除錯程式很有幫助  
5. rand('seed' ,12345);  
6. nSamples = 5000;  
7. mu = [0 0]; % TARGET MEAN目標均值  
8. rho(1) = 0.8; % rho_21目標方差  
9. rho(2) = 0.8; % rho_12目標方差  
10. % INITIALIZE THE GIBBS SAMPLER  
11. propSigma = 1; % PROPOSAL VARIANCE  
12. minn = [-3 -3];  
13. maxx = [3 3];  
14. % INITIALIZE SAMPLES  
15. x = zeros(nSamples,2);  
16. x(1,1) = unifrnd(minn(1), maxx(1));%unifrnd生成連續均勻分佈的隨機數  
17. x(1,2) = unifrnd(minn(2), maxx(2));  
18. dims = 1:2; % INDEX INTO EACH DIMENSION  
19. % RUN GIBBS SAMPLER  
20. t = 1;  
21. while t <nSamples%總共取樣出5000個取樣點  
22.     t = t + 1;  
23.     T = [t-1,t];  
24.     for iD = 1:2 % LOOP OVER DIMENSIONS總共兩維，註釋先討論第一維  
25.         % UPDATE SAMPLES  
26.         nIx = dims~=iD; % *NOT* THE CURRENT DIMENSION找到另外一維nIx=[0 1]logical型別  
27.         % CONDITIONAL MEAN  
28.         muCond = mu(iD) + rho(iD)*(x(T(iD),nIx)-mu(nIx));%計算均值=表示式μ(1)+ρ(1)*(x(n,2)-μ(2))其中x(n,2)代表樣本第n個數據的第二維  
29.         % CONDITIONAL VARIANCE  
30.         varCond = sqrt(1-rho(iD)^2);%計算方差  
31.         % DRAW FROM CONDITIONAL  
32.         x(t,iD) = normrnd(muCond,varCond);%正態分佈隨機函式，計算得到當前第t個數據的第1維資料value  
33.     end  
34. end  
35. % DISPLAY SAMPLING DYNAMICS  
36. figure;  
37. h1 = scatter(x(:,1),x(:,2),'r.');%scatter描繪散點圖，x為橫座標，y為縱座標  
38. % CONDITIONAL STEPS/SAMPLES  
39. hold on;%畫出前五十個取樣點  
40. for t = 1:50  
41.     plot([x(t,1),x(t+1,1)],[x(t,2),x(t,2)],'k-');  
42.     plot([x(t+1,1),x(t+1,1)],[x(t,2),x(t+1,2)],'k-');  
43.     h2 = plot(x(t+1,1),x(t+1,2),'ko');  
44. end  
45. h3 = scatter(x(1,1),x(1,2),'go','Linewidth',3);  
46. legend([h1,h2,h3],{'Samples','1st 50 Samples','x(t=0)'},'Location','Northwest')  
47. hold off;  
48. xlabel('x_1');  
49. ylabel('x_2');  
50. axis square