No BSD License  

Highlights from

image thumbnail




06 Nov 2002 (Updated )

Pattern analysis toolbox.

%DEMMDN1 Demonstrate fitting a multi-valued function using a Mixture Density Network.
%	Description
%	The problem consists of one input variable X and one target variable
%	T with data generated by sampling T at equal intervals and then
%	generating target data by computing T + 0.3*SIN(2*PI*T) and adding
%	Gaussian noise. A Mixture Density Network with 3 centres in the
%	mixture model is trained by minimizing a negative log likelihood
%	error function using the scaled conjugate gradient optimizer.
%	The conditional means, mixing coefficients and variances are plotted
%	as a function of X, and a contour plot of the full conditional
%	density is also generated.
%	See also

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

% Generate the matrix of inputs x and targets t.
seedn = 42;
seed = 42;
randn('state', seedn);
rand('state', seed);
ndata = 300;			% Number of data points.
noise = 0.2;			% Range of noise distribution.
t = [0:1/(ndata - 1):1]';
x = t + 0.3*sin(2*pi*t) + noise*rand(ndata, 1) - noise/2;
axis_limits = [-0.2 1.2 -0.2 1.2];

disp('This demonstration illustrates the use of a Mixture Density Network')
disp('to model multi-valued functions.  The data is generated from the')
disp('mapping x = t + 0.3 sin(2 pi t) + e, where e is a noise term.')
disp('We begin by plotting the data.')
disp(' ')
disp('Press any key to continue')
% Plot the data
fh1 = figure;
p1 = plot(x, t, 'ob');
hold on
disp('Note that for x in the range 0.35 to 0.65, there are three possible')
disp('branches of the function.')
disp(' ')
disp('Press any key to continue')

% Set up network parameters.
nin = 1;			% Number of inputs.
nhidden = 5;			% Number of hidden units.
ncentres = 3;			% Number of mixture components.
dim_target = 1;			% Dimension of target space
mdntype = '0';			% Currently unused: reserved for future use
alpha = 100;			% Inverse variance for weight initialisation
				% Make variance small for good starting point

% Create and initialize network weight vector.
net = mdn(nin, nhidden, ncentres, dim_target, mdntype);
init_options = zeros(1, 18);
init_options(1) = -1;	% Suppress all messages
init_options(14) = 10;  % 10 iterations of K means in gmminit
net = mdninit(net, alpha, t, init_options);

% Set up vector of options for the optimiser.
options = foptions;
options(1) = 1;			% This provides display of error values.
options(14) = 200;		% Number of training cycles. 

disp('We initialise the neural network model, which is an MLP with a')
disp('Gaussian mixture model with three components and spherical variance')
disp('as the error function.  This enables us to model the complete')
disp('conditional density function.')
disp(' ')
disp('Next we train the model for 200 epochs using a scaled conjugate gradient')
disp('optimizer.  The error function is the negative log likelihood of the')
disp('training data.')
disp(' ')
disp('Press any key to continue.')

% Train using scaled conjugate gradients.
[net, options] = netopt(net, options, x, t, 'scg');

disp(' ')
disp('Press any key to continue.')

disp('We can also train a conventional MLP with sum of squares error function.')
disp('This will approximate the conditional mean, which is not always a')
disp('good representation of the data.  Note that the error function is the')
disp('sum of squares error on the training data, which accounts for the')
disp('different values from training the MDN.')
disp(' ')
disp('We train the network with the quasi-Newton optimizer for 80 epochs.')
disp(' ')
disp('Press any key to continue.')
mlp_nhidden = 8;
net2 = mlp(nin, mlp_nhidden, dim_target, 'linear');
options(14) = 80; 
[net2, options] = netopt(net2, options, x, t, 'quasinew');
disp(' ')
disp('Press any key to continue.')

disp('Now we plot the underlying function, the MDN prediction,')
disp('represented by the mode of the conditional distribution, and the')
disp('prediction of the conventional MLP.')
disp(' ')
disp('Press any key to continue.')

% Plot the original function, and the trained network function.
plotvals = [0:0.01:1]';
mixes = mdn2gmm(mdnfwd(net, plotvals));
yplot = t+0.3*sin(2*pi*t);
p2 = plot(yplot, t, '--y');

% Use the mode to represent the function
y = zeros(1, length(plotvals));
priors = zeros(length(plotvals), ncentres);
c = zeros(length(plotvals), 3);
widths = zeros(length(plotvals), ncentres);
for i = 1:length(plotvals)
  [m, j] = max(mixes(i).priors);
  y(i) = mixes(i).centres(j,:);
  c(i,:) = mixes(i).centres';
p3 = plot(plotvals, y, '*r');
p4 = plot(plotvals, mlpfwd(net2, plotvals), 'g');
set(p4, 'LineWidth', 2);
legend([p1 p2 p3 p4], 'data', 'function', 'MDN mode', 'MLP mean', 4);
hold off

disp('We can also plot how the mixture model parameters depend on x.')
disp('First we plot the mixture centres, then the priors and finally')
disp('the variances.')
disp(' ')
disp('Press any key to continue.')
fh2 = figure;
subplot(3, 1, 1)
plot(plotvals, c)
hold on
title('Mixture centres')
legend('centre 1', 'centre 2', 'centre 3')
hold off

priors = reshape([mixes.priors], mixes(1).ncentres, size(mixes, 2))';
%%fh3 = figure;
subplot(3, 1, 2)
plot(plotvals, priors)
hold on
title('Mixture priors')
legend('centre 1', 'centre 2', 'centre 3')
hold off

variances = reshape([mixes.covars], mixes(1).ncentres, size(mixes, 2))';
%%fh4 = figure;
subplot(3, 1, 3)
plot(plotvals, variances)
hold on
title('Mixture variances')
legend('centre 1', 'centre 2', 'centre 3')
hold off

disp('The last figure is a contour plot of the conditional probability')
disp('density generated by the Mixture Density Network.  Note how it')
disp('is well matched to the regions of high data density.')
disp(' ')
disp('Press any key to continue.')
% Contour plot for MDN.
i = 0:0.01:1.0;
j = 0:0.01:1.0;

[I, J] = meshgrid(i,j);
I = I(:);
J = J(:);
li = length(i);
lj = length(j);
Z = zeros(li, lj);
for k = 1:li;
  Z(:,k) = gmmprob(mixes(k), j');
fh5 = figure;
% Set up levels by hand to make a good figure
v = [2 2.5 3 3.5 5:3:18];
contour(i, j, Z, v)
hold on
title('Contour plot of conditional density')
hold off

disp(' ')
disp('Press any key to exit.')
%%clear all;

Contact us