%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plotting samples from a Gaussian process %%
%% Venkatraman Narayanan, 16-831, Fall 2012 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Plotting parameters
eps = 0.01;  % Input sampling resolution 
num_samples = 10;   % Number of samples from the GP
kernel_type = 'indicator'; % gaussian, laplace, linear, indicator
kernel_scale = 0.5; % used only for gaussian and laplace


%% Choice of kernels for the GP
switch kernel_type
    case 'gaussian'
        % Gaussian/RBF kernel
        k = @(x1,x2)(exp(-(norm(x1-x2,2))^2/(2*kernel_scale^2)));
        
    case 'laplace'
        % Laplace kernel
        k = @(x1,x2)(exp(-norm(x1-x2,1)/kernel_scale));
    case 'linear'
        % Linear kernel
        k = @(x1,x2)(x1'*x2);
    case 'indicator'
        % Indicator kernel
        k = @(x1,x2)double(sum(x1~=x2)==0);
end

%% Uniformly sample inputs X, construct K, draw samples from the GP, plot
% Sample points on the real axis
x = 0:eps:1;
N = numel(x);

% Generate giant kernel matrix
S = zeros(N);
for i = 1:N
    for j = i:N
        S(i,j) = k(x(i),x(j));
        if i ~= j
            S(j,i) = S(i,j);
        end
    end
end

% Mean for the GP
mu = zeros(1,N);

% Sample from the GP
y = mvnrnd(mu,S,num_samples);

plot(x,y');axis tight;
hold on; plot(x,zeros(size(x)),'r--','linewidth',2);
title(sprintf('Samples from a zero-mean GP with %s kernel',kernel_type))  


