PRML 그림 3.7의 재현 코드
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
Reference
이 문제에 관하여(PRML 그림 3.7의 재현 코드), 우리는 이곳에서 더 많은 자료를 발견하고 링크를 클릭하여 보았다 https://qiita.com/naoya_t/items/c31d850e516646f712d3텍스트를 자유롭게 공유하거나 복사할 수 있습니다.하지만 이 문서의 URL은 참조 URL로 남겨 두십시오.
우수한 개발자 콘텐츠 발견에 전념 (Collection and Share based on the CC Protocol.)