Markov Chain Monte Carlo (MCMC) and Gibbs Sampling example(MATLAB)

Definition 1. Markov Chain

 (x0, x1, xk, .. )의 Random Variable이 있을 때, 모든 k에 대해서 P(x_k | x_k-1, x_k2, ... x_0) = P(x_k | x_k-1)을 만족한다. 즉 k번 째 Random Variable(x_k)의 probability는 k-1번째 Random Variable(x_k-1)에만 의존한다. 

 이 때 P_k(x_k | x_k-1)은 transition probability distribution(TPD) at time step k라고 불린다. Markov Chain은 TPD가 모든 time step k에 대해서 같을 때 stationary transition probability를 가진다. 이러한 stationary transition probability를 가지는 Markov Chain은 P(x|y)라고 표기할 수 있다. 

Definition 2

 stationary transition probability distribution P(x|y)를 가지는 Markov Chain의 경우, π(x) = ∫ P(x|y)π(y)dy 를 만족하는 π(x)가 있다면, π(x)는 Markov Chain의 stationary distribution이다. 

Gibbs Sampling

 Gibbs sampling은 Metropolis-Hastings sampling의 특수한 경우이다. 이것은 high dimensional distribution에서 sampling하는데 매우 유용하다. 예를 들어서 x가 m-dimensional하다고 할 때, 각 time step에서, Gibbs sampling은 x의 element들은 순차적으로 얻는다. Proposal Distribution이 x의 모든 element들에 있어서 full conditional하다면, acceptance probability는 항상 1일 것이고, 우리는 accept/reject step을 건너뛰어도 된다. 이것은 매우 큰 강점인데, Metro-Hastings sampling에 비해 update할 때 high rejection rate에 대해서 고려하지 않아도 되기 때문이다. 아래 표는 Gibbs sampling에 대해서 설명하고 있다. 

f(x)는 x의 PDF이고, x = (x_1, x_2, x_3, ... x_m)' 는 m-dimensional 하다. 

n개의 sampling을 얻어내고자 할 때 

---

• for i = 1 : n    % 각 sample에 대해서 
• for j = 1 : m    % x의 각 dimension에 대해서 
• draw x_i{ j } from f(x_i | x_i{1}, x_i{2}, ... ,  x_i{j-1},  x_i-1{j+1}, ... ,  x_i -1 {m} )    
• %  x_i{j} 를 제외한 나머지로부터! 이 때 j 번째 dimension이후로는 i-1번째 element를 이용한다. (아직 안구해져서)
• end for
• end for

---

(x_1, x_2, ... , x_n ) 은 f(x)로부터 얻은 sample 들이다. 

참고 자료

mcmc-gibbs-intro.pdf

매트랩 예제 

결과 

소스

clc;

clear all;

close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% x_1 | x_2 ~ Normal( -1 + (x_2-1)/4 , 15/4 ) ;

% x_2 | x_1 ~ Normal( 1 + (x_1+1)/4, 15/4 ) ;

% 일 때 x_1과 x_2의 joiint distribution을 구하라.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = 10;        % the number of samples

m = 2;          % the dimension of random variabe

xSample = zeros(2, n);

% Set Initial value

xSample(1, 1) = -4;

xSample(2, 1) = -3;

for epoch = 1:4

    n = n*10;

    % Run Gibbs Sampling 

    for i = 2:n

        for j = 1:m

            if j == 1

                xSample(1,i) = -1 + (xSample(2, i-1)-1)/4 + (15/4)*randn(1);

            else

                xSample(2,i) = 1 + (xSample(1, i)+1)/4 + (15/4)*randn(1);

            end

        end

    end

    % Surface Plot 3-d histogram

    % 1. Data to plot

    % x = 10*(rand(10000, 1)-0.5);

    % y = 10*(rand(10000, 1)-0.5);

    x = xSample(1, :)';

    y = xSample(2, :)';

    [sizex1 sizex2] = size(x);

    [sizey1 sizey2] = size(y);

    % 2. Range and resolution of histogram

    

    %figure(epoch);

    subplot(2, 2, epoch);

    LL = -20;     % Lower Limit

    UL = 20;      % Upper Limit

    dt = 0.3;   % Delta t

    SHIFTv = 0;

    SHIFTw = 0;

    v = LL - SHIFTv : dt : UL - SHIFTv ;

    w = LL - SHIFTw : dt : UL - SHIFTw ;

    [X,Y] = meshgrid(v,w);

    % 3. Plot specifics

    pdfz1 = hist3([x y], {v w});

    pdfz1 = pdfz1'/sizex1; %!!!!!!!! Transpose pdfz!!!!

    surf(X , Y, pdfz1)

    xlabel('X');

    ylabel('Y');

    zlabel('Joint Probability Mass');

    str = sprintf(' %d gibbs sampling ', n);

    title(str);

    drawnow();

end