diff options
| author | ziejd2 | 2017-09-28 15:04:40 -0500 |
|---|---|---|
| committer | ziejd2 | 2017-09-28 15:04:40 -0500 |
| commit | 8070dc963753142bb86c4ed698d91fd623ed28e7 (patch) | |
| tree | d0f6dd8fc46a49b819aa55c1a90faa14d8448883 /sourcecodes/bnt-master/netlab3.3/demgp.m | |
| parent | 7cc31810d53176e805532b2789955f4eedbce6bb (diff) | |
| download | BNW-8070dc963753142bb86c4ed698d91fd623ed28e7.tar.gz | |
BNW using Octave instead of Matlab.
This version of BNW should perform the same as the original version. The only difference is that it uses Octave instead of Matlab when running BayesNet Toolbox during parameter learning. I am calling this BNW_1.02. It can be accessed at: compbio.uthsc.edu/BNW_1.02
Diffstat (limited to 'sourcecodes/bnt-master/netlab3.3/demgp.m')
| -rw-r--r-- | sourcecodes/bnt-master/netlab3.3/demgp.m | 160 |
1 files changed, 160 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/netlab3.3/demgp.m b/sourcecodes/bnt-master/netlab3.3/demgp.m new file mode 100644 index 00000000..65c587a3 --- /dev/null +++ b/sourcecodes/bnt-master/netlab3.3/demgp.m @@ -0,0 +1,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; \ No newline at end of file |
