PRML図6.9の再現

gistfile1.matlab

1;
clear

%global eta = [1 1]; %% 左図
global eta = [1 0.01]; %% 右図

global theta0 = 0.5;

N = 30;  %% N=30は割と時間がかかる
xl = linspace(-1, 1, N);

xs = zeros(2, N*N);
for x1 = 1:N
  for x2 = 1:N
    j = (x2-1)*N + x1;
    xs(:,j) = [x1; x2];
  endfor
endfor

function [retval] = ard_kernel(x1, x2)
  global eta;
  global theta0;
  sum = eta(1)*(x1(1)-x2(1))^2 + eta(2)*(x1(2)-x2(2))^2;
  retval = theta0 * exp(-1/2*sum);
endfunction

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

k = @ard_kernel;

cov = outer(xs, xs, k);

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

cov += eye(N*N)*beta;

graph = 1;
if graph == 1
  clf;
  hold on;
  axis([-1 1 -1 1 -3 3]);
  colormap(gray)

  R = big_rand(zeros(N*N,1), cov);
  [xx, yy] = meshgrid(xl, xl);
  surf(xx, yy, reshape(R, [N N]))
  view(-20, 30);
  hold off;
endif