Metropolis-Hastingsアルゴリズム ref: http://qiita.com/items/b9aa110b395d0a5b0b0a

mh1.matlab

1;
%
% 図11.9 Metropolisアルゴリズムを用いて2次元ガウス分布からサンプリングする例。
%

function plot_gauss2_contour(mu, Sigma)
  cont = zeros(2,101);
  D = chol(Sigma, "lower");
  for j = 1:100
    theta = pi*2/100 * (j-1);
    cont(:,j) = mu + D*[cos(theta); sin(theta)];
  endfor
  cont(:,101) = cont(:,1);
  plot(cont(1,:), cont(2,:), 'k')
endfunction

function pdf = multival_gaussian_pdf(mu, Sigma)
  D = numel(mu);
  Z = ((2*pi)^(D/2)*sqrt(det(Sigma))); % 正規化定数1/Z
  inv_Sigma = inv(Sigma);

  pdf = @(x) 1/Z*exp(-1/2*(x-mu)'*inv_Sigma*(x-mu));
endfunction


mu = [1.5; 1.5];
Sigma = [1.125 0.875; 0.875 1.125]; % λ1=2, v1=[1; 1], λ2=1/4, v2=[1; -1]
p = multival_gaussian_pdf(mu, Sigma);

clf
hold on
axis([0,3,0,3]);
title("Fig 11.9: Metropolis");

plot_gauss2_contour(mu, Sigma);

T = 150;
z = mu;
for j = 1:T
  z_proposal = z + 0.2*randn(2,1);
  A = min(1, p(z_proposal)/p(z)); % (11.33) 受理確率A(z*, z)
  u = rand(1);
  if A > u
    plot([z(1,1) z_proposal(1,1)], [z(2,1) z_proposal(2,1)], 'g')
    z = z_proposal;
  else
    plot([z(1,1) z_proposal(1,1)], [z(2,1) z_proposal(2,1)], 'r')
  endif
endfor

hold off

mh2.matlab

1;
%%
%% 図11.9 Metropolisアルゴリズムを用いて2次元ガウス分布からサンプリングする例。
%%

function p = multival_gaussian(x, mu, Sigma)
  D = rows(x);
  dx = x - mu;
  Z = ((2*pi)^(D/2)*sqrt(det(Sigma))); % 正規化定数1/Z
  p = 1/Z*exp(-1/2*dx'*inv(Sigma)*dx);
endfunction

function pdf = multival_gaussian_pdf(mu, Sigma)
  D = numel(mu);
  Z = ((2*pi)^(D/2)*sqrt(det(Sigma))); % 正規化定数1/Z
  inv_Sigma = inv(Sigma);

  pdf = @(x) 1/Z*exp(-1/2*(x-mu)'*inv_Sigma*(x-mu));
endfunction

function plot_gauss2_contour(mu, Sigma, space)
%  [x,y] = meshgrid(space);
%  gauss = multival_gaussian_pdf(mu, Sigma);
%  z = arrayfun(@(x,y)gauss([x;y]), x, y);
%  contour(space,space,z);

  cont = zeros(2,101);
  D = chol(Sigma, "lower");
  for j = 1:100
    theta = pi*2/100 * (j-1);
    cont(:,j) = mu + D*[cos(theta); sin(theta)];
  endfor
  cont(:,101) = cont(:,1);
  plot(cont(1,:), cont(2,:), 'k')
endfunction

global s1 = 10;
global s2 = 0.5;
global mu = [1.5; 1.5]*s1/sqrt(2);
% Sigma = [1.125 0.875; 0.875 1.125]; % λ1=2, v1=[1; 1], λ2=1/4, v2=[1; -1]
H = [1 1; 1 -1]/sqrt(2);
global Sigma = H * [s1*s1 0; 0 s2*s2] * H;
global p = multival_gaussian_pdf(mu, Sigma);

function plot_m(rho, caption, A, line, step=1)
  global s1;
  global s2;
  global mu;
  global Sigma;

  printf("%s: ", caption);

  clf
  hold on
  axis([0,3,0,3]*s1/sqrt(2));
  title(caption);

  plot_gauss2_contour(mu, Sigma);

  T = 150 * step;
  acc = 0;
  z = mu;

  for j = 1:T
    z_proposal = z + rho*randn(2,1);
    p_accept = A(z_proposal, z);
    u = rand(1);
    if p_accept > u
      if line
        plot([z(1,1) z_proposal(1,1)], [z(2,1) z_proposal(2,1)], 'g')
      else
        if mod(j, step) == 0
          plot([z_proposal(1,1)], [z_proposal(2,1)], 'x')
        endif
      endif
      z = z_proposal;
      acc += 1;
    else
      if line
        plot([z(1,1) z_proposal(1,1)], [z(2,1) z_proposal(2,1)], 'r')
      endif
    endif
  endfor

  hold off

  printf("acceptance rate = %g¥n", acc/T);
endfunction


function metropolis(rho, caption, skip=0)
  global p;
  A = @(z_,z) min(1, p(z_)/p(z)); % (11.33)
  plot_m(rho, caption, A, true, 1)
endfunction

function metropolis_hastings(rho, caption, skip=0)
  global p;
  q = @(z_,z) multival_gaussian(z_, z, rho^2*eye(2));
  A = @(z_,z) min(1, (p(z_)*q(z,z_)) / (p(z)*q(z_,z))); % (11.44)
  if skip == 0
    plot_m(rho, caption, A, true, 1)
  else
    plot_m(rho, caption, A, false, 1+skip)
  endif
endfunction


metropolis(0.2, "Metropolis")
pause

metropolis_hastings(0.2, "Metropolis-Hastings, with rho << sigma_2")
pause

metropolis_hastings(s1, "Metropolis-Hastings, with rho=sigma_{max}")
pause

metropolis_hastings(s2, "Metropolis-Hastings, with rho=sigma_2")
pause

metropolis_hastings(s2, "Metropolis-Hastings, with rho=sigma_2 (skip=9)", 9)