|
wolffd@0
|
1 function samples=demmet1(plot_wait)
|
|
wolffd@0
|
2 %DEMMET1 Demonstrate Markov Chain Monte Carlo sampling on a Gaussian.
|
|
wolffd@0
|
3 %
|
|
wolffd@0
|
4 % Description
|
|
wolffd@0
|
5 % The problem consists of generating data from a Gaussian in two
|
|
wolffd@0
|
6 % dimensions using a Markov Chain Monte Carlo algorithm. The points are
|
|
wolffd@0
|
7 % plotted one after another to show the path taken by the chain.
|
|
wolffd@0
|
8 %
|
|
wolffd@0
|
9 % DEMMET1(PLOTWAIT) allows the user to set the time (in a whole number
|
|
wolffd@0
|
10 % of seconds) between the plotting of points. This is passed to PAUSE
|
|
wolffd@0
|
11 %
|
|
wolffd@0
|
12 % See also
|
|
wolffd@0
|
13 % DEMHMC1, METROP, GMM, DEMPOT
|
|
wolffd@0
|
14 %
|
|
wolffd@0
|
15
|
|
wolffd@0
|
16 % Copyright (c) Ian T Nabney (1996-2001)
|
|
wolffd@0
|
17
|
|
wolffd@0
|
18 if nargin == 0 | plot_wait < 0
|
|
wolffd@0
|
19 plot_wait = 0; % No wait if not specified or incorrect
|
|
wolffd@0
|
20 end
|
|
wolffd@0
|
21 dim = 2; % Data dimension
|
|
wolffd@0
|
22 ncentres = 1; % Number of centres in mixture model
|
|
wolffd@0
|
23
|
|
wolffd@0
|
24 seed = 42; % Seed for random weight initialization.
|
|
wolffd@0
|
25 randn('state', seed);
|
|
wolffd@0
|
26 rand('state', seed);
|
|
wolffd@0
|
27
|
|
wolffd@0
|
28 clc
|
|
wolffd@0
|
29 disp('This demonstration illustrates the use of the Markov chain Monte Carlo')
|
|
wolffd@0
|
30 disp('algorithm to sample from a Gaussian distribution.')
|
|
wolffd@0
|
31 disp('The mean is at [0 0].')
|
|
wolffd@0
|
32 disp(' ')
|
|
wolffd@0
|
33 disp('First we set up the parameters of the mixture model we are sampling')
|
|
wolffd@0
|
34 disp('from.')
|
|
wolffd@0
|
35 disp(' ')
|
|
wolffd@0
|
36 disp('Press any key to continue.')
|
|
wolffd@0
|
37 pause
|
|
wolffd@0
|
38
|
|
wolffd@0
|
39 % Set up mixture model to sample from
|
|
wolffd@0
|
40 mix = gmm(dim, ncentres, 'spherical');
|
|
wolffd@0
|
41 mix.centres(1, :) = [0 0];
|
|
wolffd@0
|
42 x = [0 4]; % Start vector
|
|
wolffd@0
|
43
|
|
wolffd@0
|
44 % Set up vector of options for hybrid Monte Carlo.
|
|
wolffd@0
|
45
|
|
wolffd@0
|
46 nsamples = 150; % Number of retained samples.
|
|
wolffd@0
|
47
|
|
wolffd@0
|
48 options = foptions; % Default options vector.
|
|
wolffd@0
|
49 options(1) = 0; % Switch off diagnostics.
|
|
wolffd@0
|
50 options(14) = nsamples; % Number of Monte Carlo samples returned.
|
|
wolffd@0
|
51 options(18) = 0.1;
|
|
wolffd@0
|
52
|
|
wolffd@0
|
53 clc
|
|
wolffd@0
|
54 disp('Next we take 150 samples from the distribution.')
|
|
wolffd@0
|
55 disp('Sampling starts at the point [0 4].')
|
|
wolffd@0
|
56 disp('The new state is accepted if the threshold value is greater than')
|
|
wolffd@0
|
57 disp('a random number between 0 and 1.')
|
|
wolffd@0
|
58 disp(' ')
|
|
wolffd@0
|
59 disp('Press any key to continue.')
|
|
wolffd@0
|
60 pause
|
|
wolffd@0
|
61
|
|
wolffd@0
|
62 [samples, energies] = metrop('dempot', x, options, '', mix);
|
|
wolffd@0
|
63
|
|
wolffd@0
|
64 clc
|
|
wolffd@0
|
65 disp('The plot shows the samples generated by the MCMC function in order')
|
|
wolffd@0
|
66 disp('as an animation to show the path taken by the Markov chain.')
|
|
wolffd@0
|
67 disp('The different colours are used to show that the first few samples')
|
|
wolffd@0
|
68 disp('should be discarded as they lie too far from the mean.')
|
|
wolffd@0
|
69 disp(' ')
|
|
wolffd@0
|
70 disp('Press any key to continue.')
|
|
wolffd@0
|
71 pause
|
|
wolffd@0
|
72 probs = exp(-energies);
|
|
wolffd@0
|
73 fh1 = figure;
|
|
wolffd@0
|
74 g1end = floor(nsamples/4);
|
|
wolffd@0
|
75
|
|
wolffd@0
|
76 for n = 1:nsamples
|
|
wolffd@0
|
77
|
|
wolffd@0
|
78 if n < g1end
|
|
wolffd@0
|
79 Marker = 'k.';
|
|
wolffd@0
|
80 p1 = plot(samples(n,1), samples(n,2), Marker, ...
|
|
wolffd@0
|
81 'EraseMode', 'none', 'MarkerSize', 12);
|
|
wolffd@0
|
82 if n == 1
|
|
wolffd@0
|
83 axis([-3 5 -2 5])
|
|
wolffd@0
|
84 end
|
|
wolffd@0
|
85 else
|
|
wolffd@0
|
86 Marker = 'r.';
|
|
wolffd@0
|
87 p2 = plot(samples(n,1), samples(n,2), Marker, ...
|
|
wolffd@0
|
88 'EraseMode', 'none', 'MarkerSize', 12);
|
|
wolffd@0
|
89 end
|
|
wolffd@0
|
90 hold on
|
|
wolffd@0
|
91 drawnow; % Force drawing immediately
|
|
wolffd@0
|
92 pause(plot_wait);
|
|
wolffd@0
|
93 end
|
|
wolffd@0
|
94 lstrings = char(['Samples 1-' int2str(g1end)], ...
|
|
wolffd@0
|
95 ['Samples ' int2str(g1end+1) '-' int2str(nsamples)]);
|
|
wolffd@0
|
96 legend([p1 p2], lstrings, 1);
|
|
wolffd@0
|
97
|
|
wolffd@0
|
98 disp(' ')
|
|
wolffd@0
|
99 disp('Press any key to exit.')
|
|
wolffd@0
|
100 pause
|
|
wolffd@0
|
101 close(fh1);
|
|
wolffd@0
|
102 clear all;
|
|
wolffd@0
|
103
|
