xt = (-1:0.1:6)';
ft = sin( xt*2 );

x = [rand(80,1);rand(20,1)*3+1];
x = sort(x);
fx = sin(x*2);
yx = fx + randn(size(fx))*0.3;

%---
for ii=1:4
    switch ii
        case 1
          z = ([0  5])';
        case 2
          z = ([0 3 5])';
        case 3
          z = (0:1:5)';
        case 4
          z = (0:0.5:5)';
    end
sigma = 1;
Kxx = gk(x,x,sigma);
Kzz = gk(z,z,sigma);
Kxz = gk(x,z,sigma);
Kxt = gk(x,xt,sigma);
Ktt = gk(xt,xt,sigma);
Kzt = gk(z,xt,sigma);

ss = 1;
xt_hat = Kxt' * inv(Kxx+ss*eye(length(Kxx))) * yx;
Sigma_xt = Ktt - Kxt' * inv(Kxx+ss*eye(length(Kxx)))*Kxt;
w = sqrt(diag(Sigma_xt));

Kzzi = inv(Kzz);
Sinv = (1/sigma)*Kzzi*Kxz'*Kxz*Kzzi+Kzzi;
S = inv(Sinv);
m = S*Kzzi*Kxz'*yx;

figure(1)
subplot(2,2,ii)
hold off
plot( xt, ft, 'w-' )
hold on
xt0 = [xt;xt(end:-1:1)];
yt0 = [xt_hat+w;xt_hat(end:-1:1)-w(end:-1:1)];

patch( xt0, yt0, 'b',...
    'EdgeColor','none','FaceAlpha',0.2 )
plot( xt, xt_hat, 'b-')
plot( x, yx, '.')
plot( x, x*0-1.5, 'b+', [-1,6],-1.5+[0,0],'b-')
ylim([-2,2.5])

if ii>1

s = sqrt(diag(S));
for i=1:length(z)
%    plot( z(i)+[0,0], fz(i)+[s,-s], 'r-')
    plot( z(i)+[0,0], m(i)+[s(i),-s(i)], 'r-')
end

plot( z, z*0-1.7, 'r+', [-1,6],-1.7+[0,0],'r-')

ft_hat = Kzt' * Kzzi * m;
Sigma_ft = Ktt - Kzt'*Kzzi*Kzt...
    + Kzt' * Kzzi * S * Kzzi * Kzt;
w = sqrt(diag(Sigma_ft));
xt0 = [xt;xt(end:-1:1)];
ft0 = [ft_hat+w;ft_hat(end:-1:1)-w(end:-1:1)];
patch( xt0, ft0, 'r',...
    'EdgeColor','none','FaceAlpha',0.2 )
plot( xt, ft_hat, 'r-')
end
if ii~=1
    text(2,2,sprintf('M=%d',length(z)-1))
end
set(gca, 'YTick',[-1.7, -1.5, -1:1:2], 'YTickLabel',...
    {'z','x','-1','0','1','2'})
xlim([-1,5])
end
saveas(gcf,'gp1d.png','png')

