clear all;
close all;
clc;
% f sim GP(m,k) here f is a function this is a Gaussian Process
% the function f is distributrd as a GP woth mean m and covariance 
% function k

n = 51;
xi = linspace(-5,5,n)';
mu = 1/4*xi.^2;
Sigma = exp(-1/2*pdist2(xi,xi).^2);

% f sim N(mu,Sigma) here f is a vector this is a Gaussian distribution

% --FROM RASMUSSEN
xs = (-5:0.2:5)'; ns = size(xs,1); keps = 1e-9;
m = inline('0.25*x.^2');
K = inline('exp(-0.5 * ( repmat(p'',size(q)) - repmat(q,size(p'')) ).^2)');
fs = m(xs) + chol(K(xs,xs)+keps*eye(ns))' * randn(ns,1);
plot(xs,fs,'--')
% --
hold on;
fs1 = m(xs) + chol(K(xs,xs)+keps*eye(ns))' * randn(ns,1);
plot(xs,fs1,'--')
fs2 = m(xs) + chol(K(xs,xs)+keps*eye(ns))' * randn(ns,1);
plot(xs,fs2,'--')
for mm=1:10
    fsn = m(xs) + chol(K(xs,xs)+keps*eye(ns))' * randn(ns,1);
    plot(xs,fsn,'--')
end
fs3 = m(xs); 
plot(xs,fs3,'-')
lb = m(xs) - 2*sqrt(diag(Sigma)) ;
ub = m(xs) + 2*sqrt(diag(Sigma)) ;
% plot(xs,lb,'-')
% plot(xs,ub,'-')
color = 'b';
edgecolor ='k';
transparency = 0.2;
fill([xs; flipud(xs)], [lb; flipud(ub)], color, 'EdgeColor',edgecolor,'FaceAlpha',transparency,'EdgeAlpha',transparency);