https://www.youtube.com/watch?v=16oPvgOd3UI

function GP_1d

kernel=5;

switch kernel

case 1; k = @(x,y) 1*x'*y; % linear

case 2; k = @(x,y) 1*min(x,y); % brownian motion

case 3; k = @(x,y) exp(-100*(x-y)'*(x-y)); % squared

case 4; k = @(x,y) exp(-1*sqrt(x-y)'*(x-y)); % Ornistin

case 5; k = @(x,y) exp(-1*sin(5*pi*(x-y))^2); % periodic

end

% choose points at which to sample

x = (0:.005:1);

n = length(x);

% covariance matrix

C = zeros(n,n);

for i=1:n

for j=1:n

C(i,j) = k (x(i), x(j)) ;

end

end

% sample from gaussian process at this points

u = randn(n,1);

[A, S , B ] = svd(C);

z = A *sqrt(S)*u;

% plot

figure(2); hold on;

plot(x, z, '.-');

axis([0, 1, -2, 2]);

end

============ IN 2D ========

function GP_2d

kernel=3;

switch kernel

case 1; k = @(x,y) 1*x'*y; % linear

case 2; k = @(x,y) 1*min(x,y); % brownian motion

case 3; k = @(x,y) exp(-100*(x-y)'*(x-y)); % squared

case 4; k = @(x,y) exp(-1*sqrt(x-y)'*(x-y)); % Ornistin

case 5; k = @(x,y) exp(-1*sin(5*pi*(x-y))^2 ); % periodic

end

% choose points at which to sample

points = (0:0.05:1)';

[U,V] = meshgrid(points,points);

x = [U(:) V(:)]';

n = size(x,2);

% covariance matrix

C = zeros(n,n);

for i=1:n

for j=1:n

C(i,j) = k (x(:,i), x(:,j)) ;

end

end

% sample from gaussian process at this points

u = randn(n,1);

[A, S , B ] = svd(C);

z = A*sqrt(S)*u;

% plot

figure(2); clf;

Z = reshape(z,sqrt(n) , sqrt(n));

surf(U,V,Z);

end