view toolboxes/FullBNT-1.0.7/netlab3.3/demgp.m @ 0:e9a9cd732c1e tip

first hg version after svn

%DEMGP	Demonstrate simple regression using a Gaussian Process.
%
%	Description
%	The problem consists of one input variable X and one target variable
%	T. The values in X are chosen in two separated clusters and the
%	target data is generated by computing SIN(2*PI*X) and adding Gaussian
%	noise. Two Gaussian Processes, each with different covariance
%	functions are trained by optimising the hyperparameters  using the
%	scaled conjugate gradient algorithm.  The final predictions are
%	plotted together with 2 standard deviation error bars.
%
%	See also
%	GP, GPERR, GPFWD, GPGRAD, GPINIT, SCG
%

%	Copyright (c) Ian T Nabney (1996-2001)


% Find out if flops is available (i.e. pre-version 6 Matlab)
v = version;
if (str2num(strtok(v, '.')) >= 6)
    flops_works = logical(0);
else
    flops_works = logical(1);
end

randn('state', 42);
x = [0.1 0.15 0.2 0.25  0.65 0.7 0.75 0.8 0.85 0.9]';
ndata = length(x);
t = sin(2*pi*x) + 0.05*randn(ndata, 1);

xtest = linspace(0, 1, 50)';

clc
disp('This demonstration illustrates the use of a Gaussian Process')
disp('model for regression problems.  The data is generated from a noisy')
disp('sine function.')
disp(' ')
disp('Press any key to continue.')
pause

flops(0);
% Initialise the parameters.
net = gp(1, 'sqexp');
prior.pr_mean = 0;
prior.pr_var = 1;
net = gpinit(net, x, t, prior);

clc
disp('The first GP uses the squared exponential covariance function.')
disp('The hyperparameters are initialised by sampling from a Gaussian with a')
disp(['mean of ', num2str(prior.pr_mean), ' and variance ', ...
    num2str(prior.pr_var), '.'])
disp('After initializing the network, we train it using the scaled conjugate')
disp('gradients algorithm for 20 cycles.')
disp(' ')
disp('Press any key to continue')
pause

% Now train to find the hyperparameters.
options = foptions;
options(1) = 1;    % Display training error values
options(14) = 20;
flops(0)
[net, options] = netopt(net, options, x, t, 'scg');
if flops_works
    sflops = flops;
end

disp('The second GP uses the rational quadratic covariance function.')
disp('The hyperparameters are initialised by sampling from a Gaussian with a')
disp(['mean of ', num2str(prior.pr_mean), ' and variance ', num2str(prior.pr_var)])
disp('After initializing the network, we train it using the scaled conjugate')
disp('gradients algorithm for 20 cycles.')
disp(' ')
disp('Press any key to continue')
pause
flops(0)
net2 = gp(1, 'ratquad');
net2 = gpinit(net2, x, t, prior);
flops(0)
[net2, options] = netopt(net2, options, x, t, 'scg');
if flops_works
    rflops = flops;
end

disp(' ')
disp('Press any key to continue')
disp(' ')
pause
clc

fprintf(1, 'For squared exponential covariance function,');
if flops_works    
    fprintf(1, 'flops = %d', sflops);
end
fprintf(1, '\nfinal hyperparameters:\n')
format_string = strcat('  bias:\t\t\t%10.6f\n  noise:\t\t%10.6f\n', ...
  '  inverse lengthscale:\t%10.6f\n  vertical scale:\t%10.6f\n');
fprintf(1, format_string, ...
    exp(net.bias), exp(net.noise), exp(net.inweights(1)), exp(net.fpar(1)));
fprintf(1, '\n\nFor rational quadratic covariance function,');
if flops_works
    fprintf(1, 'flops = %d', rflops);
end
fprintf(1, '\nfinal hyperparameters:\n')
format_string = [format_string '  cov decay order:\t%10.6f\n'];
fprintf(1, format_string, ...
      exp(net2.bias), exp(net2.noise), exp(net2.inweights(1)), ...
      exp(net2.fpar(1)), exp(net2.fpar(2)));
disp(' ')
disp('Press any key to continue')
pause

disp(' ')
disp('Now we plot the data, underlying function, model outputs and two')
disp('standard deviation error bars on a single graph to compare the results.')
disp(' ')
disp('Press any key to continue.')
pause
cn = gpcovar(net, x); 
cninv = inv(cn);
[ytest, sigsq] = gpfwd(net, xtest, cninv);
sig = sqrt(sigsq);

fh1 = figure;
hold on
plot(x, t, 'ok');
xlabel('Input')
ylabel('Target')
fplot('sin(2*pi*x)', [0 1], '--m');
plot(xtest, ytest, '-k');
plot(xtest, ytest+(2*sig), '-b', xtest, ytest-(2*sig), '-b');
axis([0 1 -1.5 1.5]);
title('Squared exponential covariance function')
legend('data', 'function', 'GP', 'error bars');
hold off

cninv2 = inv(gpcovar(net2, x));
[ytest2, sigsq2] = gpfwd(net2, xtest, cninv2);
sig2 = sqrt(sigsq2);
fh2 = figure;
hold on
plot(x, t, 'ok');
xlabel('Input')
ylabel('Target')
fplot('sin(2*pi*x)', [0 1], '--m');
plot(xtest, ytest2, '-k');
plot(xtest, ytest2+(2*sig2), '-b', xtest, ytest2-(2*sig2), '-b');
axis([0 1 -1.5 1.5]);
title('Rational quadratic covariance function')
legend('data', 'function', 'GP', 'error bars');
hold off

disp(' ')
disp('Press any key to end.')
pause
close(fh1);
close(fh2);
clear all;