comparison toolboxes/FullBNT-1.0.7/netlab3.3/demmet1.m @ 0:e9a9cd732c1e tip

first hg version after svn

+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;
+