PRML図3.7の再現コード

（誰得OctaveコードをQiitaに貼り付けていくプロジェクト）

1;

function g = gaussian(x, mu,sigma)
  g = 1/sqrt(2*pi*sigma^2) * exp(-1/(2*sigma^2)*(x-mu)^2);
endfunction

function g = gaussian2(x, mu,Sigma)
  g = 1/(2*pi*sqrt(det(Sigma))) * exp(-1/2*(x-mu)'*inv(Sigma)*(x-mu));
endfunction

function plotfun2(fun, param)
  x = linspace(-1, 1, 3);
  plot(x, arrayfun(fun, x), param);
endfunction

function plotfun3(fun)
  [xx, yy] = meshgrid(linspace(-1, 1, 51));
  mesh(xx, yy, arrayfun(fun, xx, yy));
  view(-15, 80);
endfunction

function savegif(id)
  filename = sprintf("fig37_%03d.gif", id);
  print(filename, "-dgif");
endfunction

N = 20;
a = [-0.3 0.5];

xs = rand(N, 1)*2-1; % U(x|-1,1)
ts = arrayfun(@(x)[1 x]*a', xs) + 0.2*randn(N, 1);

alpha = 2.0;
beta  = 25; % = (1/0.2)^2

%% 初期値
mu       = [0; 0];
Sigma    = eye(2);
invSigma = inv(Sigma);

for i = 1:N
  printf("w空間の事前分布¥n");
  clf;
  hold on;
  title("prior/posterior");
  xlabel("w0");
  ylabel("w1");
  plot3(a(1), a(2), 0, 'x');
  plotfun3(@(w0,w1) gaussian2([w0;w1],mu,Sigma));
  hold off;
  % savegif((i-1)*3+1);
  pause

  printf("現在の事前分布から w をランダムに6個選んで y(x,w) をプロット¥n");
  clf;
  hold on;
  axis([-1, 1, -1, 1])
  title("6 samples from w");
  xlabel("x");
  ylabel("y");
  for j = 1:6
    %% 二次元ガウス分布 N(_mu,_Sigma) に従う乱数を求める
    L = chol(Sigma);
    w = mu + L*rand(2, 1);
    plotfun2(@(x) [1 x]*w, 'r');
  endfor
  plot(xs(1:i), ts(1:i), 'o');
  drawnow()
  % savegif((i-1)*3+2);
  hold off;
  pause

  xi = xs(i);
  ti = ts(i);
  phi = [1 xi]; % 計画行列 Φ は単純に [1 x] で

  printf("尤度関数 p(t|x, w) を wの関数としてプロット¥n");
  clf;
  hold on;
  title("likelihood");
  xlabel("w0");
  ylabel("w1");
  plotfun3(@(w0,w1) gaussian(ti, phi*[w0;w1], 1/beta));
  % savegif((i-1)*3+3);
  hold off;
  pause

  printf("「事後分布 = 事前分布 × 尤度」¥n");
  % 式 (3.50) (3.51) で分布の μ と Σ を更新
  next_invSigma = invSigma + beta*phi'*phi;
  next_Sigma    = inv(next_invSigma);
  next_mu       = next_Sigma * (invSigma*mu + beta*phi'*ti);

  Sigma    = next_Sigma;
  invSigma = next_invSigma;
  mu       = next_mu;
endfor

アニメーションGIF化

$ convert -delay 100 fig37_*.gif -loop 0 fig37_anim.gif