Gaussian Processes demo

gp.m

% All based on Rasmussen and Williams.

N = 1000;
lo=0; hi=5;
x = linspace(lo, hi, N);
%xn = [2, 2.5, 3]';
%yn = [-1.9, -2, -1.9]';
xn = ([-4, -3, -1, 0, 2]'+5)/2;
yn =  [-2,  0,  1, 2, -1]';

% alphas = [100]; aa = length(alphas);
% ells = [1]; ll = length(ells);
alphas = 10.^(-2:2); aa = length(alphas);
ells = 10.^(-2:2); ll = length(ells);

for i=1:3; figure(i); clf; end
for j=1:ll
    for i=1:aa
        alpha = alphas(i); ell = ells(j);
        cov = @(r) alpha*exp(-1/(2*ell^2)*r.^2);

        [P,Q] = ndgrid(x, x);
        KStarStar = cov(P-Q) + 1e-6*eye(N);
        KStarStarChol = chol(KStarStar,'lower');

        % Plot the unconditioned draws.
        figure(1)
        subplot(aa, ll, ll*(i-1)+j)
        %subaxis(aa, ll, j, i, 'Spacing', 0.03, 'Padding', 0, 'Margin', 0.03)
        for k=1:10
            plot(x,KStarStarChol*randn(N,1));
            hold on
        end
        if j==1; ylabel(sprintf('$\\alpha$=%g', alpha), 'interpreter', 'latex'); end
        if i==aa; xlabel(sprintf('$\\ell$=%g', ell), 'interpreter', 'latex'); end
        set(gca,'xtick',[])
        xlim([lo, hi])
        
        % Plot mixtures
        figure(2)
        subplot(aa, ll, ll*(i-1)+j)
        %subaxis(aa, ll, j, i, 'Spacing', 0.03, 'Padding', 0, 'Margin', 0.03)
        s1 = KStarStarChol*randn(N,1);
        s2 = KStarStarChol*randn(N,1);
        for theta=linspace(-pi, pi,10)
            plot(x, s1*cos(theta) + s2*sin(theta))
            hold on
        end
        plot(x, s1, 'r')
        plot(x, s2, 'r')
        if j==1; ylabel(sprintf('$\\alpha$=%g', alpha), 'interpreter', 'latex'); end
        if i==aa; xlabel(sprintf('$\\ell$=%g', ell), 'interpreter', 'latex'); end
        set(gca,'xtick',[])
        xlim([lo, hi])
        
        % Compute predictive distributions
        [P,Q] = ndgrid(xn, xn);
        K = cov(P-Q) + 1e-6*eye(length(xn));
        Kchol = chol(K,'lower');
        % vecAlpha = (K + sigma^2 I)^-1 * yn
        vecAlpha = Kchol' \ (Kchol \ yn);

        [P,Q] = ndgrid(xn, x);
        Kstar = cov(P-Q);
        predictiveMean = Kstar'*vecAlpha;
        
        vecV = Kchol \ Kstar;
        predictiveVar = KStarStar - vecV'*vecV;
        varChol = chol(predictiveVar+1e-6*eye(N), 'lower');
        logLikelihood = -(.5*yn'*vecAlpha + sum(log(diag(Kchol))) + length(xn)/2*log(2*pi));
        
        % Plot samples from predictive distributions
        figure(3)
        %subaxis(aa, ll, j, i, 'Spacing', 0.03, 'Padding', 0, 'Margin', 0.03)
        subplot(aa, ll, ll*(i-1)+j)
        v_ = zeros(N,1);
        for q=1:N
            [P,Q] = ndgrid(xn, x(q));
            vecVV = Kchol\cov(P-Q);
            v_(q) = 1e-6 - vecVV'*vecVV;
        end
        v_ = 1e-6+diag(predictiveVar);
        fill([x'; flipud(x')],[predictiveMean+2*sqrt(v_); flipud(predictiveMean-2*sqrt(v_))], [.5 .5 .5])
        hold on
        for sample=1:10
            plot(x, predictiveMean+varChol*randn(N,1));
        end
        plot(xn, yn, 'xr')
        plot(x, predictiveMean, 'g')
        title(sprintf('LL=%g', logLikelihood))
        if j==1; ylabel(sprintf('$\\alpha$=%g', alpha), 'interpreter', 'latex'); end
        if i==aa; xlabel(sprintf('$\\ell$=%g', ell), 'interpreter', 'latex'); end
        set(gca,'xtick',[])
        xlim([lo, hi])
    end
    for fig=1:3; figure(fig); set(gca, 'xtickMode', 'auto'); end
end
figure(1); save_tight_pdf('GPdraws.pdf', 8, 6)
figure(2); save_tight_pdf('GPblends.pdf', 8, 6)
figure(3); save_tight_pdf('GPconditional.pdf', 8, 6)