clc
clf
%k = @(x,y) 1*x'*y;
k= @(x,y) exp(-((x-y)'*(x-y))/(2.*0.1*0.1));

% sample points
Xtest= (0:.001:1);
n = length(Xtest);

% construct covariance matrix
C = zeros(n,n);
for i=1:n
    for j=1:n
        C(i,j)= k(Xtest(i),Xtest(j));
    end
end

% sample from GP at these points
u = randn(n,1);

[A, S, B]= svd(C);
z = A*sqrt(S)*u;  % z = A S^0.5

figure(2); hold on;
plot(Xtest,z,'.-')
axis([0, 1, -3, 3])

for i=2:5
    u = randn(n,1);
    [A, S, B]= svd(C);
    z = A*sqrt(S)*u;  % z = A S^0.5 u
    %   plot(Xtest,z,'.-')
end

% noiseless training data

Xtrain= [0.2 ; 0.49];  % location
ytrain= [0.1 ; -0.4];  %  sample value

% apply the kernel function to the training points
for i=1:length(Xtrain)
    for j=1:length(Xtrain)
        K(i,j)= k(Xtest(i),Xtest(j));
    end
end

[A, S, B]= svd(K);
L= A*sqrt(S);

% compute mean at the test points
for i=1:length(Xtrain)
    for j=1:length(Xtest)
        K_s(i,j)= k(Xtrain(i), Xtest(j));
    end
end

% Lk = L*K_s;
% mu = Lk'*L*ytrain;

%%%%%%%%%

for i=1:length(Xtrain)
    for j=1:length(Xtrain)
        K11(i,j)= k(Xtrain(i),Xtrain(j));
    end
end

for i=1:length(Xtrain)
    for j=1:length(Xtest)
        K12(i,j)= k(Xtrain(i), Xtest(j));
    end
end

K21= K12';

for i=1:length(Xtest)
    for j=1:length(Xtest)
        K22(i,j)= k(Xtest(i), Xtest(j));
    end
end

K = [K11 K12; K21 K22];

m= K21*inv(K11)*ytrain;

C = zeros(n,n);
for i=1:n
    for j=1:n
        C(i,j)= k(Xtest(i),Xtest(j));
    end
end

C= C - K21*inv(K11)*K12;

u = randn(n,1);
[A, S, B]= svd(C);
z = m + A*sqrt(S)*u;  % z = A S^0.5 u
plot(Xtest,z,'.-')

plot(Xtrain(1),ytrain(1),'s','MarkerFaceColor','r', 'MarkerSize', 10);
plot(Xtrain(2),ytrain(2),'s','MarkerFaceColor','r', 'MarkerSize', 10);
