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 들이다.
참고 자료
매트랩 예제
결과
소스
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
'Enginius > Machine Learning' 카테고리의 다른 글
Random Forest Description (0) | 2012.06.07 |
---|---|
ANN로 headPoseEstimate하기 (0) | 2012.06.05 |
Complex Zernike moments (0) | 2012.05.21 |
푸리에 변환과 라플라스 변환 (18) | 2012.05.01 |
Hierarchical Temporal Memory Prototype ver1.0.0 (4) | 2012.04.26 |