wolffd@0: %DEMRBF1 Demonstrate simple regression using a radial basis function network.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	The problem consists of one input variable X and one target variable
wolffd@0: %	T with data generated by sampling X at equal intervals and then
wolffd@0: %	generating target data by computing SIN(2*PI*X) and adding Gaussian
wolffd@0: %	noise. This data is the same as that used in demmlp1.
wolffd@0: %
wolffd@0: %	Three different RBF networks (with different activation functions)
wolffd@0: %	are trained in two stages. First, a Gaussian mixture model is trained
wolffd@0: %	using the EM algorithm, and the centres of this model are used to set
wolffd@0: %	the centres of the RBF.  Second, the output weights (and biases) are
wolffd@0: %	determined using the pseudo-inverse of the design matrix.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	DEMMLP1, RBF, RBFFWD, GMM, GMMEM
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: 
wolffd@0: % Generate the matrix of inputs x and targets t.
wolffd@0: randn('state', 42);
wolffd@0: rand('state', 42);
wolffd@0: ndata = 20;			% Number of data points.
wolffd@0: noise = 0.2;			% Standard deviation of noise distribution.
wolffd@0: x = (linspace(0, 1, ndata))';
wolffd@0: t = sin(2*pi*x) + noise*randn(ndata, 1);
wolffd@0: mu = mean(x);
wolffd@0: sigma = std(x);
wolffd@0: tr_in = (x - mu)./(sigma);
wolffd@0: 
wolffd@0: clc
wolffd@0: disp('This demonstration illustrates the use of a Radial Basis Function')
wolffd@0: disp('network 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: % Set up network parameters.
wolffd@0: nin = 1;			% Number of inputs.
wolffd@0: nhidden = 7;			% Number of hidden units.
wolffd@0: nout = 1;			% Number of outputs.
wolffd@0: 
wolffd@0: clc
wolffd@0: disp('We assess the effect of three different activation functions.')
wolffd@0: disp('First we create a network with Gaussian activations.')
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: % Create and initialize network weight and parameter vectors.
wolffd@0: net = rbf(nin, nhidden, nout, 'gaussian');
wolffd@0: 
wolffd@0: disp('A two-stage training algorithm is used: it uses a small number of')
wolffd@0: disp('iterations of EM to position the centres, and then the pseudo-inverse')
wolffd@0: disp('of the design matrix to find the second layer weights.')
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: disp('Error values from EM training.')
wolffd@0: % Use fast training method
wolffd@0: options = foptions;
wolffd@0: options(1) = 1;		% Display EM training
wolffd@0: options(14) = 10;	% number of iterations of EM
wolffd@0: net = rbftrain(net, options, tr_in, t);
wolffd@0: 
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: clc
wolffd@0: disp('The second RBF network has thin plate spline activations.')
wolffd@0: disp('The same centres are used again, so we just need to calculate')
wolffd@0: disp('the second layer weights.')
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: % Create a second RBF with thin plate spline functions
wolffd@0: net2 = rbf(nin, nhidden, nout, 'tps');
wolffd@0: 
wolffd@0: % Re-use previous centres rather than calling rbftrain again
wolffd@0: net2.c = net.c;
wolffd@0: [y, act2] = rbffwd(net2, tr_in);
wolffd@0: 
wolffd@0: % Solve for new output weights and biases from RBF activations
wolffd@0: temp = pinv([act2 ones(ndata, 1)]) * t;
wolffd@0: net2.w2 = temp(1:nhidden, :);
wolffd@0: net2.b2 = temp(nhidden+1, :);
wolffd@0: 
wolffd@0: disp('The third RBF network has r^4 log r activations.')
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: % Create a third RBF with r^4 log r functions
wolffd@0: net3 = rbf(nin, nhidden, nout, 'r4logr');
wolffd@0: 
wolffd@0: % Overwrite weight vector with parameters from first RBF
wolffd@0: net3.c = net.c;
wolffd@0: [y, act3] = rbffwd(net3, tr_in);
wolffd@0: temp = pinv([act3 ones(ndata, 1)]) * t;
wolffd@0: net3.w2 = temp(1:nhidden, :);
wolffd@0: net3.b2 = temp(nhidden+1, :);
wolffd@0: 
wolffd@0: disp('Now we plot the data, underlying function, and network outputs')
wolffd@0: disp('on a single graph to compare the results.')
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue.')
wolffd@0: pause
wolffd@0: % Plot the data, the original function, and the trained network functions.
wolffd@0: plotvals = [x(1):0.01:x(end)]';
wolffd@0: inputvals = (plotvals-mu)./sigma;
wolffd@0: y = rbffwd(net, inputvals);
wolffd@0: y2 = rbffwd(net2, inputvals);
wolffd@0: y3 = rbffwd(net3, inputvals);
wolffd@0: fh1 = figure;
wolffd@0: 
wolffd@0: plot(x, t, 'ob')
wolffd@0: hold on
wolffd@0: xlabel('Input')
wolffd@0: ylabel('Target')
wolffd@0: axis([x(1) x(end) -1.5 1.5])
wolffd@0: [fx, fy] = fplot('sin(2*pi*x)', [x(1) x(end)]);
wolffd@0: plot(fx, fy, '-r', 'LineWidth', 2)
wolffd@0: plot(plotvals, y, '--g', 'LineWidth', 2)
wolffd@0: plot(plotvals, y2, 'k--', 'LineWidth', 2)
wolffd@0: plot(plotvals, y3, '-.c', 'LineWidth', 2)
wolffd@0: legend('data', 'function', 'Gaussian RBF', 'Thin plate spline RBF', ...
wolffd@0:   'r^4 log r RBF');
wolffd@0: hold off
wolffd@0: 
wolffd@0: disp('RBF training errors are');
wolffd@0: disp(['Gaussian ', num2str(rbferr(net, tr_in, t)), ' TPS ',  ...
wolffd@0: num2str(rbferr(net2, tr_in, t)), ' R4logr ', num2str(rbferr(net3, tr_in, t))]);
wolffd@0: 
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to end.')
wolffd@0: pause
wolffd@0: close(fh1);
wolffd@0: clear all;