1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
%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;
|
