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
% m(x)=1/4*x^2
% k(x,x')=exp(-1/2*(x-x')^2)

n = 51;
xi = linspace(-5,5,n)';
mu = 1/4*xi.^2;
Sigma = exp(-1/2*pdist2(xi,xi).^2);
% Sigma = exp(-1/2*squareform(pdist(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);

% L=chol(K(xs,xs)+keps*eye(ns))';
% R=randn(ns,1);

% xx=1;
% m(xx)
% K(xx,xx)
% 
% xxx = [-3+m(xx):.1:3+m(xx)];
% yyy = normpdf(xxx,m(xx),sqrt(K(xx,xx)));
% plot(xx+yyy,xxx,'b');
% plot(xx+[yyy(1) yyy(end)],[xxx(1) xxx(end)],'b')
% grid on