function samples=demmet1(plot_wait)
%DEMMET1 Demonstrate Markov Chain Monte Carlo sampling on a Gaussian.
%
%	Description
%	The problem consists of generating data from a Gaussian in two
%	dimensions using a Markov Chain Monte Carlo algorithm. The points are
%	plotted one after another to show the path taken by the chain.
%
%	DEMMET1(PLOTWAIT) allows the user to set the time (in a whole number
%	of seconds) between the plotting of points.  This is passed to PAUSE
%
%	See also
%	DEMHMC1, METROP, GMM, DEMPOT
%

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

if nargin == 0 | plot_wait < 0
  plot_wait = 0; % No wait if not specified or incorrect
end
dim = 2;            	% Data dimension
ncentres = 1;		% Number of centres in mixture model

seed = 42;              % Seed for random weight initialization.
randn('state', seed);
rand('state', seed);

clc
disp('This demonstration illustrates the use of the Markov chain Monte Carlo')
disp('algorithm to sample from a Gaussian distribution.')
disp('The mean is at [0 0].')
disp(' ')
disp('First we set up the parameters of the mixture model we are sampling')
disp('from.')
disp(' ')
disp('Press any key to continue.')
pause

% Set up mixture model to sample from
mix = gmm(dim, ncentres, 'spherical');
mix.centres(1, :) = [0 0];
x = [0 4];  % Start vector

% Set up vector of options for hybrid Monte Carlo.

nsamples = 150;		% Number of retained samples.

options = foptions;     % Default options vector.
options(1) = 0;		% Switch off diagnostics.
options(14) = nsamples;	% Number of Monte Carlo samples returned. 
options(18) = 0.1;

clc
disp('Next we take 150 samples from the distribution.')
disp('Sampling starts at the point [0 4].')
disp('The new state is accepted if the threshold value is greater than')
disp('a random number between 0 and 1.')
disp(' ')
disp('Press any key to continue.')
pause

[samples, energies] = metrop('dempot', x, options, '', mix);

clc
disp('The plot shows the samples generated by the MCMC function in order')
disp('as an animation to show the path taken by the Markov chain.')
disp('The different colours are used to show that the first few samples')
disp('should be discarded as they lie too far from the mean.')
disp(' ')
disp('Press any key to continue.')
pause
probs = exp(-energies);
fh1 = figure;
g1end = floor(nsamples/4);

for n = 1:nsamples
  
  if n < g1end
    Marker = 'k.';
    p1 = plot(samples(n,1), samples(n,2), Marker, ...
      'EraseMode', 'none', 'MarkerSize', 12);
    if n == 1
      axis([-3 5 -2 5])
    end
  else
    Marker = 'r.';
    p2 = plot(samples(n,1), samples(n,2), Marker, ...
      'EraseMode', 'none', 'MarkerSize', 12);
  end
  hold on
  drawnow;  % Force drawing immediately
  pause(plot_wait);
end
lstrings = char(['Samples 1-' int2str(g1end)], ...
  ['Samples ' int2str(g1end+1) '-' int2str(nsamples)]);
legend([p1 p2], lstrings, 1);

disp(' ')
disp('Press any key to exit.')
pause
close(fh1);
clear all;