Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/netlab3.3/demhmc3.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 %DEMHMC3 Demonstrate Bayesian regression with Hybrid Monte Carlo sampling.
wolffd@0 2 %
wolffd@0 3 % Description
wolffd@0 4 % The problem consists of one input variable X and one target variable
wolffd@0 5 % T with data generated by sampling X at equal intervals and then
wolffd@0 6 % generating target data by computing SIN(2*PI*X) and adding Gaussian
wolffd@0 7 % noise. The model is a 2-layer network with linear outputs, and the
wolffd@0 8 % hybrid Monte Carlo algorithm (with persistence) is used to sample
wolffd@0 9 % from the posterior distribution of the weights. The graph shows the
wolffd@0 10 % underlying function, 300 samples from the function given by the
wolffd@0 11 % posterior distribution of the weights, and the average prediction
wolffd@0 12 % (weighted by the posterior probabilities).
wolffd@0 13 %
wolffd@0 14 % See also
wolffd@0 15 % DEMHMC2, HMC, MLP, MLPERR, MLPGRAD
wolffd@0 16 %
wolffd@0 17
wolffd@0 18 % Copyright (c) Ian T Nabney (1996-2001)
wolffd@0 19
wolffd@0 20
wolffd@0 21 % Generate the matrix of inputs x and targets t.
wolffd@0 22 ndata = 20; % Number of data points.
wolffd@0 23 noise = 0.1; % Standard deviation of noise distribution.
wolffd@0 24 nin = 1; % Number of inputs.
wolffd@0 25 nout = 1; % Number of outputs.
wolffd@0 26
wolffd@0 27 seed = 42; % Seed for random number generators.
wolffd@0 28 randn('state', seed);
wolffd@0 29 rand('state', seed);
wolffd@0 30
wolffd@0 31 x = 0.25 + 0.1*randn(ndata, nin);
wolffd@0 32 t = sin(2*pi*x) + noise*randn(size(x));
wolffd@0 33
wolffd@0 34 clc
wolffd@0 35 disp('This demonstration illustrates the use of the hybrid Monte Carlo')
wolffd@0 36 disp('algorithm to sample from the posterior weight distribution of a')
wolffd@0 37 disp('multi-layer perceptron.')
wolffd@0 38 disp(' ')
wolffd@0 39 disp('A regression problem is used, with the one-dimensional data drawn')
wolffd@0 40 disp('from a noisy sine function. The x values are sampled from a normal')
wolffd@0 41 disp('distribution with mean 0.25 and variance 0.01.')
wolffd@0 42 disp(' ')
wolffd@0 43 disp('First we initialise the network.')
wolffd@0 44 disp(' ')
wolffd@0 45 disp('Press any key to continue.')
wolffd@0 46 pause
wolffd@0 47
wolffd@0 48 % Set up network parameters.
wolffd@0 49 nhidden = 5; % Number of hidden units.
wolffd@0 50 alpha = 0.001; % Coefficient of weight-decay prior.
wolffd@0 51 beta = 100.0; % Coefficient of data error.
wolffd@0 52
wolffd@0 53 % Create and initialize network model.
wolffd@0 54
wolffd@0 55 % Initialise weights reasonably close to 0
wolffd@0 56 net = mlp(nin, nhidden, nout, 'linear', alpha, beta);
wolffd@0 57 net = mlpinit(net, 10);
wolffd@0 58
wolffd@0 59 clc
wolffd@0 60 disp('Next we take 100 samples from the posterior distribution. The first')
wolffd@0 61 disp('300 samples at the start of the chain are omitted. As persistence')
wolffd@0 62 disp('is used, the momentum has a small random component added at each step.')
wolffd@0 63 disp('10 iterations are used at each step (compared with 100 in demhmc2).')
wolffd@0 64 disp('The step size is 0.005 (compared with 0.002).')
wolffd@0 65 disp('The new state is accepted if the threshold')
wolffd@0 66 disp('value is greater than a random number between 0 and 1.')
wolffd@0 67 disp(' ')
wolffd@0 68 disp('Negative step numbers indicate samples discarded from the start of the')
wolffd@0 69 disp('chain.')
wolffd@0 70 disp(' ')
wolffd@0 71 disp('Press any key to continue.')
wolffd@0 72 pause
wolffd@0 73
wolffd@0 74 % Set up vector of options for hybrid Monte Carlo.
wolffd@0 75 nsamples = 100; % Number of retained samples.
wolffd@0 76
wolffd@0 77 options = foptions; % Default options vector.
wolffd@0 78 options(1) = 1; % Switch on diagnostics.
wolffd@0 79 options(5) = 1; % Use persistence
wolffd@0 80 options(7) = 10; % Number of steps in trajectory.
wolffd@0 81 options(14) = nsamples; % Number of Monte Carlo samples returned.
wolffd@0 82 options(15) = 300; % Number of samples omitted at start of chain.
wolffd@0 83 options(17) = 0.95; % Alpha value in persistence
wolffd@0 84 options(18) = 0.005; % Step size.
wolffd@0 85
wolffd@0 86 w = mlppak(net);
wolffd@0 87 % Initialise HMC
wolffd@0 88 hmc('state', 42);
wolffd@0 89 [samples, energies] = hmc('neterr', w, options, 'netgrad', net, x, t);
wolffd@0 90
wolffd@0 91 clc
wolffd@0 92 disp('The plot shows the underlying noise free function, the 100 samples')
wolffd@0 93 disp('produced from the MLP, and their average as a Monte Carlo estimate')
wolffd@0 94 disp('of the true posterior average.')
wolffd@0 95 disp(' ')
wolffd@0 96 disp('Press any key to continue.')
wolffd@0 97 pause
wolffd@0 98
wolffd@0 99 nplot = 300;
wolffd@0 100 plotvals = [0 : 1/(nplot - 1) : 1]';
wolffd@0 101 pred = zeros(size(plotvals));
wolffd@0 102 fh1 = figure;
wolffd@0 103 hold on
wolffd@0 104 for k = 1:nsamples
wolffd@0 105 w2 = samples(k,:);
wolffd@0 106 net2 = mlpunpak(net, w2);
wolffd@0 107 y = mlpfwd(net2, plotvals);
wolffd@0 108 % Sum predictions
wolffd@0 109 pred = pred + y;
wolffd@0 110 h4 = plot(plotvals, y, '-r', 'LineWidth', 1);
wolffd@0 111 end
wolffd@0 112 pred = pred./nsamples;
wolffd@0 113 % Plot data
wolffd@0 114 h1 = plot(x, t, 'ob', 'LineWidth', 2, 'MarkerFaceColor', 'blue');
wolffd@0 115 axis([0 1 -3 3])
wolffd@0 116
wolffd@0 117 % Plot function
wolffd@0 118 [fx, fy] = fplot('sin(2*pi*x)', [0 1], '--g');
wolffd@0 119 h2 = plot(fx, fy, '--g', 'LineWidth', 2);
wolffd@0 120 set(gca, 'box', 'on');
wolffd@0 121
wolffd@0 122 % Plot averaged prediction
wolffd@0 123 h3 = plot(plotvals, pred, '-c', 'LineWidth', 2);
wolffd@0 124
wolffd@0 125 lstrings = char('Data', 'Function', 'Prediction', 'Samples');
wolffd@0 126 legend([h1 h2 h3 h4], lstrings, 3);
wolffd@0 127 hold off
wolffd@0 128
wolffd@0 129 disp('Note how the predictions become much further from the true function')
wolffd@0 130 disp('away from the region of high data density.')
wolffd@0 131 disp(' ')
wolffd@0 132 disp('Press any key to exit.')
wolffd@0 133 pause
wolffd@0 134 close(fh1);
wolffd@0 135 clear all;