% This is a little script to demonstrate how the variance of the
% posterior diminishes at points near the samples. This is based on the
% example in: Carl E. Rasmussen. "Gaussian Processes in Machine
% Learning." It was written by Peter Carbonetto and Nando de Freitas.

clear;
N = 10;        % The number of training points.
s = 0.04;      % Noise variance.
h = 2;         % Kernel parameter (sigma).

% Randomly generate training points on [-5,5].
X = -5 + 10*rand(N,1);
x = (-5:0.1:5)';  % The test points.
n = size(x,1);    % The number of test pts.

% Construct the mean and covariance functions.
%m = inline('0.25*x.^2', 'x');
m = inline('sin(0.9*x)','x');
K = inline(['exp((-1/(h^2))*(repmat(transpose(p),size(q)) - ' ...
	    'repmat(q,size(transpose(p)))).^2)'], 'p', 'q','h');

% Generate random training labels.
F = m(X) + chol(K(X,X,h)+s*eye(N))'*randn(N,1);
M = m(X);

% Demonstrante how to sample functions from the prior:
figure(1)
clf;
L = 20;           % Number of samples
f = zeros(n,L);
for i=1:L,
  f(:,i) = m(x) + sqrtm(K(x,x,h)+s*eye(n))'*randn(n,1);
end;
plot(x,f);
xlabel('time','fontsize',15);
ylabel('functions from the prior');

% COMPUTE MEAN AND VARIANCE;
S = eye(N)*s + K(X,X,h);
y = m(x) + K(X,x,h)*inv(S)*(F - M); 
y = y';
for i = 1:n
  xi   = x(i);  
  c(i) = s + K(xi,xi,h) - K(X,xi,h)*inv(S)*K(xi,X,h);
end

% Plot the mean and 95% confidence intervals.
figure(2)
clf
hold on
plot(x,y-2*c,'g-')
plot(x,y+2*c,'g-')
plot(x,y,'r')
plot(X,F,'bx')
legend('confidence intervals',' ','mean prediction','data','location','southwest');
hold off