PRML図6.4図6.5の再現

gistfile1.matlab

1;
clear

fig_id = 6;

thetas_6_5 = [ 1 4 0 0;
               9 4 0 0;
               1 64 0 0;
               1 0.25 0 0;
               1 4 10 0;
               1 4 0 5 ];

global thetas = thetas_6_5(fig_id,:);
global sigma = 1.2;
global theta = 1;

axes = [ -1 1 -3 3;
         -1 1 -9 9;
         -1 1 -3 3;
         -1 1 -3 3;
         -1 1 -9 9;
         -1 1 -4 4 ];

xmin = axes(fig_id, 1);
xmax = axes(fig_id, 2);

N = 100;
xs = linspace(xmin, xmax, N);

function [retval] = gauss_kernel(x1, x2)
  global sigma;
  retval = exp( -(x1-x2)^2 / (2*sigma^2) );
endfunction

function [retval] = ornstein_uhlenbeck(x1, x2)
  global theta;
  retval = exp(-theta*abs(x1-x2));
endfunction

function [retval] = my_kernel(x1, x2)
  global thetas;
  retval = thetas(1)*exp(-thetas(2)/2*(x1-x2)^2) + thetas(3) + thetas(4)*x1'*x2;
endfunction

function [R] = big_rand(Mu, Cov)
  N = length(Mu);
  L = chol(Cov, 'lower');
  z = randn(N,1);
  R = Mu + L * z;
endfunction

k = @my_kernel;
%k = @ornstein_uhlenbeck;

cov = outer(xs, xs, k);

bm = min(eig(cov));
if bm < 0
  beta_1 = -bm
else
  beta_1 = 1.0e-14;
endif

cov += eye(N)*beta_1;

graph = 1;
if graph == 1
  clf;
  hold on;
  axis(axes(fig_id, :));

  color = ['r';'g';'b';'k';'m'];
  for i = 1:5
    R = big_rand(zeros(N,1), cov);

    plot(xs, R, color(i))
  endfor

  hold off;
endif