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