Posted at

PRML図3.7の再現コード

More than 5 years have passed since last update.

(誰得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