Daniel@0: function samples=demmet1(plot_wait)
Daniel@0: %DEMMET1 Demonstrate Markov Chain Monte Carlo sampling on a Gaussian.
Daniel@0: %
Daniel@0: %	Description
Daniel@0: %	The problem consists of generating data from a Gaussian in two
Daniel@0: %	dimensions using a Markov Chain Monte Carlo algorithm. The points are
Daniel@0: %	plotted one after another to show the path taken by the chain.
Daniel@0: %
Daniel@0: %	DEMMET1(PLOTWAIT) allows the user to set the time (in a whole number
Daniel@0: %	of seconds) between the plotting of points.  This is passed to PAUSE
Daniel@0: %
Daniel@0: %	See also
Daniel@0: %	DEMHMC1, METROP, GMM, DEMPOT
Daniel@0: %
Daniel@0: 
Daniel@0: %	Copyright (c) Ian T Nabney (1996-2001)
Daniel@0: 
Daniel@0: if nargin == 0 | plot_wait < 0
Daniel@0:   plot_wait = 0; % No wait if not specified or incorrect
Daniel@0: end
Daniel@0: dim = 2;            	% Data dimension
Daniel@0: ncentres = 1;		% Number of centres in mixture model
Daniel@0: 
Daniel@0: seed = 42;              % Seed for random weight initialization.
Daniel@0: randn('state', seed);
Daniel@0: rand('state', seed);
Daniel@0: 
Daniel@0: clc
Daniel@0: disp('This demonstration illustrates the use of the Markov chain Monte Carlo')
Daniel@0: disp('algorithm to sample from a Gaussian distribution.')
Daniel@0: disp('The mean is at [0 0].')
Daniel@0: disp(' ')
Daniel@0: disp('First we set up the parameters of the mixture model we are sampling')
Daniel@0: disp('from.')
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue.')
Daniel@0: pause
Daniel@0: 
Daniel@0: % Set up mixture model to sample from
Daniel@0: mix = gmm(dim, ncentres, 'spherical');
Daniel@0: mix.centres(1, :) = [0 0];
Daniel@0: x = [0 4];  % Start vector
Daniel@0: 
Daniel@0: % Set up vector of options for hybrid Monte Carlo.
Daniel@0: 
Daniel@0: nsamples = 150;		% Number of retained samples.
Daniel@0: 
Daniel@0: options = foptions;     % Default options vector.
Daniel@0: options(1) = 0;		% Switch off diagnostics.
Daniel@0: options(14) = nsamples;	% Number of Monte Carlo samples returned. 
Daniel@0: options(18) = 0.1;
Daniel@0: 
Daniel@0: clc
Daniel@0: disp('Next we take 150 samples from the distribution.')
Daniel@0: disp('Sampling starts at the point [0 4].')
Daniel@0: disp('The new state is accepted if the threshold value is greater than')
Daniel@0: disp('a random number between 0 and 1.')
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue.')
Daniel@0: pause
Daniel@0: 
Daniel@0: [samples, energies] = metrop('dempot', x, options, '', mix);
Daniel@0: 
Daniel@0: clc
Daniel@0: disp('The plot shows the samples generated by the MCMC function in order')
Daniel@0: disp('as an animation to show the path taken by the Markov chain.')
Daniel@0: disp('The different colours are used to show that the first few samples')
Daniel@0: disp('should be discarded as they lie too far from the mean.')
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to continue.')
Daniel@0: pause
Daniel@0: probs = exp(-energies);
Daniel@0: fh1 = figure;
Daniel@0: g1end = floor(nsamples/4);
Daniel@0: 
Daniel@0: for n = 1:nsamples
Daniel@0:   
Daniel@0:   if n < g1end
Daniel@0:     Marker = 'k.';
Daniel@0:     p1 = plot(samples(n,1), samples(n,2), Marker, ...
Daniel@0:       'EraseMode', 'none', 'MarkerSize', 12);
Daniel@0:     if n == 1
Daniel@0:       axis([-3 5 -2 5])
Daniel@0:     end
Daniel@0:   else
Daniel@0:     Marker = 'r.';
Daniel@0:     p2 = plot(samples(n,1), samples(n,2), Marker, ...
Daniel@0:       'EraseMode', 'none', 'MarkerSize', 12);
Daniel@0:   end
Daniel@0:   hold on
Daniel@0:   drawnow;  % Force drawing immediately
Daniel@0:   pause(plot_wait);
Daniel@0: end
Daniel@0: lstrings = char(['Samples 1-' int2str(g1end)], ...
Daniel@0:   ['Samples ' int2str(g1end+1) '-' int2str(nsamples)]);
Daniel@0: legend([p1 p2], lstrings, 1);
Daniel@0: 
Daniel@0: disp(' ')
Daniel@0: disp('Press any key to exit.')
Daniel@0: pause
Daniel@0: close(fh1);
Daniel@0: clear all;
Daniel@0: