Daniel@0: %DEMGP Demonstrate simple regression using a Gaussian Process. Daniel@0: % Daniel@0: % Description Daniel@0: % The problem consists of one input variable X and one target variable Daniel@0: % T. The values in X are chosen in two separated clusters and the Daniel@0: % target data is generated by computing SIN(2*PI*X) and adding Gaussian Daniel@0: % noise. Two Gaussian Processes, each with different covariance Daniel@0: % functions are trained by optimising the hyperparameters using the Daniel@0: % scaled conjugate gradient algorithm. The final predictions are Daniel@0: % plotted together with 2 standard deviation error bars. Daniel@0: % Daniel@0: % See also Daniel@0: % GP, GPERR, GPFWD, GPGRAD, GPINIT, SCG Daniel@0: % Daniel@0: Daniel@0: % Copyright (c) Ian T Nabney (1996-2001) Daniel@0: Daniel@0: Daniel@0: % Find out if flops is available (i.e. pre-version 6 Matlab) Daniel@0: v = version; Daniel@0: if (str2num(strtok(v, '.')) >= 6) Daniel@0: flops_works = logical(0); Daniel@0: else Daniel@0: flops_works = logical(1); Daniel@0: end Daniel@0: Daniel@0: randn('state', 42); Daniel@0: x = [0.1 0.15 0.2 0.25 0.65 0.7 0.75 0.8 0.85 0.9]'; Daniel@0: ndata = length(x); Daniel@0: t = sin(2*pi*x) + 0.05*randn(ndata, 1); Daniel@0: Daniel@0: xtest = linspace(0, 1, 50)'; Daniel@0: Daniel@0: clc Daniel@0: disp('This demonstration illustrates the use of a Gaussian Process') Daniel@0: disp('model for regression problems. The data is generated from a noisy') Daniel@0: disp('sine function.') Daniel@0: disp(' ') Daniel@0: disp('Press any key to continue.') Daniel@0: pause Daniel@0: Daniel@0: flops(0); Daniel@0: % Initialise the parameters. Daniel@0: net = gp(1, 'sqexp'); Daniel@0: prior.pr_mean = 0; Daniel@0: prior.pr_var = 1; Daniel@0: net = gpinit(net, x, t, prior); Daniel@0: Daniel@0: clc Daniel@0: disp('The first GP uses the squared exponential covariance function.') Daniel@0: disp('The hyperparameters are initialised by sampling from a Gaussian with a') Daniel@0: disp(['mean of ', num2str(prior.pr_mean), ' and variance ', ... Daniel@0: num2str(prior.pr_var), '.']) Daniel@0: disp('After initializing the network, we train it using the scaled conjugate') Daniel@0: disp('gradients algorithm for 20 cycles.') Daniel@0: disp(' ') Daniel@0: disp('Press any key to continue') Daniel@0: pause Daniel@0: Daniel@0: % Now train to find the hyperparameters. Daniel@0: options = foptions; Daniel@0: options(1) = 1; % Display training error values Daniel@0: options(14) = 20; Daniel@0: flops(0) Daniel@0: [net, options] = netopt(net, options, x, t, 'scg'); Daniel@0: if flops_works Daniel@0: sflops = flops; Daniel@0: end Daniel@0: Daniel@0: disp('The second GP uses the rational quadratic covariance function.') Daniel@0: disp('The hyperparameters are initialised by sampling from a Gaussian with a') Daniel@0: disp(['mean of ', num2str(prior.pr_mean), ' and variance ', num2str(prior.pr_var)]) Daniel@0: disp('After initializing the network, we train it using the scaled conjugate') Daniel@0: disp('gradients algorithm for 20 cycles.') Daniel@0: disp(' ') Daniel@0: disp('Press any key to continue') Daniel@0: pause Daniel@0: flops(0) Daniel@0: net2 = gp(1, 'ratquad'); Daniel@0: net2 = gpinit(net2, x, t, prior); Daniel@0: flops(0) Daniel@0: [net2, options] = netopt(net2, options, x, t, 'scg'); Daniel@0: if flops_works Daniel@0: rflops = flops; Daniel@0: end Daniel@0: Daniel@0: disp(' ') Daniel@0: disp('Press any key to continue') Daniel@0: disp(' ') Daniel@0: pause Daniel@0: clc Daniel@0: Daniel@0: fprintf(1, 'For squared exponential covariance function,'); Daniel@0: if flops_works Daniel@0: fprintf(1, 'flops = %d', sflops); Daniel@0: end Daniel@0: fprintf(1, '\nfinal hyperparameters:\n') Daniel@0: format_string = strcat(' bias:\t\t\t%10.6f\n noise:\t\t%10.6f\n', ... Daniel@0: ' inverse lengthscale:\t%10.6f\n vertical scale:\t%10.6f\n'); Daniel@0: fprintf(1, format_string, ... Daniel@0: exp(net.bias), exp(net.noise), exp(net.inweights(1)), exp(net.fpar(1))); Daniel@0: fprintf(1, '\n\nFor rational quadratic covariance function,'); Daniel@0: if flops_works Daniel@0: fprintf(1, 'flops = %d', rflops); Daniel@0: end Daniel@0: fprintf(1, '\nfinal hyperparameters:\n') Daniel@0: format_string = [format_string ' cov decay order:\t%10.6f\n']; Daniel@0: fprintf(1, format_string, ... Daniel@0: exp(net2.bias), exp(net2.noise), exp(net2.inweights(1)), ... Daniel@0: exp(net2.fpar(1)), exp(net2.fpar(2))); Daniel@0: disp(' ') Daniel@0: disp('Press any key to continue') Daniel@0: pause Daniel@0: Daniel@0: disp(' ') Daniel@0: disp('Now we plot the data, underlying function, model outputs and two') Daniel@0: disp('standard deviation error bars on a single graph to compare the results.') Daniel@0: disp(' ') Daniel@0: disp('Press any key to continue.') Daniel@0: pause Daniel@0: cn = gpcovar(net, x); Daniel@0: cninv = inv(cn); Daniel@0: [ytest, sigsq] = gpfwd(net, xtest, cninv); Daniel@0: sig = sqrt(sigsq); Daniel@0: Daniel@0: fh1 = figure; Daniel@0: hold on Daniel@0: plot(x, t, 'ok'); Daniel@0: xlabel('Input') Daniel@0: ylabel('Target') Daniel@0: fplot('sin(2*pi*x)', [0 1], '--m'); Daniel@0: plot(xtest, ytest, '-k'); Daniel@0: plot(xtest, ytest+(2*sig), '-b', xtest, ytest-(2*sig), '-b'); Daniel@0: axis([0 1 -1.5 1.5]); Daniel@0: title('Squared exponential covariance function') Daniel@0: legend('data', 'function', 'GP', 'error bars'); Daniel@0: hold off Daniel@0: Daniel@0: cninv2 = inv(gpcovar(net2, x)); Daniel@0: [ytest2, sigsq2] = gpfwd(net2, xtest, cninv2); Daniel@0: sig2 = sqrt(sigsq2); Daniel@0: fh2 = figure; Daniel@0: hold on Daniel@0: plot(x, t, 'ok'); Daniel@0: xlabel('Input') Daniel@0: ylabel('Target') Daniel@0: fplot('sin(2*pi*x)', [0 1], '--m'); Daniel@0: plot(xtest, ytest2, '-k'); Daniel@0: plot(xtest, ytest2+(2*sig2), '-b', xtest, ytest2-(2*sig2), '-b'); Daniel@0: axis([0 1 -1.5 1.5]); Daniel@0: title('Rational quadratic covariance function') Daniel@0: legend('data', 'function', 'GP', 'error bars'); Daniel@0: hold off Daniel@0: Daniel@0: disp(' ') Daniel@0: disp('Press any key to end.') Daniel@0: pause Daniel@0: close(fh1); Daniel@0: close(fh2); Daniel@0: clear all;