|
wolffd@0
|
1 %DEMHMC1 Demonstrate Hybrid Monte Carlo sampling on mixture of two Gaussians.
|
|
wolffd@0
|
2 %
|
|
wolffd@0
|
3 % Description
|
|
wolffd@0
|
4 % The problem consists of generating data from a mixture of two
|
|
wolffd@0
|
5 % Gaussians in two dimensions using a hybrid Monte Carlo algorithm with
|
|
wolffd@0
|
6 % persistence. A mixture model is then fitted to the sample to compare
|
|
wolffd@0
|
7 % it with the true underlying generator.
|
|
wolffd@0
|
8 %
|
|
wolffd@0
|
9 % See also
|
|
wolffd@0
|
10 % DEMHMC3, HMC, DEMPOT, DEMGPOT
|
|
wolffd@0
|
11 %
|
|
wolffd@0
|
12
|
|
wolffd@0
|
13 % Copyright (c) Ian T Nabney (1996-2001)
|
|
wolffd@0
|
14
|
|
wolffd@0
|
15
|
|
wolffd@0
|
16 dim = 2; % Data dimension
|
|
wolffd@0
|
17 ncentres = 2; % Number of centres in mixture model
|
|
wolffd@0
|
18
|
|
wolffd@0
|
19 seed = 42; % Seed for random weight initialization.
|
|
wolffd@0
|
20 randn('state', seed);
|
|
wolffd@0
|
21 rand('state', seed);
|
|
wolffd@0
|
22
|
|
wolffd@0
|
23 clc
|
|
wolffd@0
|
24 disp('This demonstration illustrates the use of the hybrid Monte Carlo')
|
|
wolffd@0
|
25 disp('algorithm to sample from a mixture of two Gaussians.')
|
|
wolffd@0
|
26 disp('The means of the two components are [0 0] and [2 2].')
|
|
wolffd@0
|
27 disp(' ')
|
|
wolffd@0
|
28 disp('First we set up the parameters of the mixture model we are sampling')
|
|
wolffd@0
|
29 disp('from.')
|
|
wolffd@0
|
30 disp(' ')
|
|
wolffd@0
|
31 disp('Press any key to continue.')
|
|
wolffd@0
|
32 pause
|
|
wolffd@0
|
33
|
|
wolffd@0
|
34 % Set up mixture model to sample from
|
|
wolffd@0
|
35 mix = gmm(dim, ncentres, 'spherical');
|
|
wolffd@0
|
36 mix.centres(1, :) = [0 0];
|
|
wolffd@0
|
37 mix.centres(2, :) = [2 2];
|
|
wolffd@0
|
38 x = [0 1]; % Start vector
|
|
wolffd@0
|
39
|
|
wolffd@0
|
40 % Set up vector of options for hybrid Monte Carlo.
|
|
wolffd@0
|
41
|
|
wolffd@0
|
42 nsamples = 160; % Number of retained samples.
|
|
wolffd@0
|
43
|
|
wolffd@0
|
44 options = foptions; % Default options vector.
|
|
wolffd@0
|
45 options(1) = 1; % Switch on diagnostics.
|
|
wolffd@0
|
46 options(5) = 1; % Use persistence
|
|
wolffd@0
|
47 options(7) = 50; % Number of steps in trajectory.
|
|
wolffd@0
|
48 options(14) = nsamples; % Number of Monte Carlo samples returned.
|
|
wolffd@0
|
49 options(15) = 30; % Number of samples omitted at start of chain.
|
|
wolffd@0
|
50 options(18) = 0.02;
|
|
wolffd@0
|
51
|
|
wolffd@0
|
52 clc
|
|
wolffd@0
|
53 disp(['Next we take ', num2str(nsamples),' samples from the distribution.'...
|
|
wolffd@0
|
54 , 'The first ', num2str(options(15))])
|
|
wolffd@0
|
55 disp('samples at the start of the chain are omitted. As persistence')
|
|
wolffd@0
|
56 disp('is used, the momentum has a small random component added at each step.')
|
|
wolffd@0
|
57 disp([num2str(options(7)), ...
|
|
wolffd@0
|
58 ' iterations are used at each step and the step size is ',...
|
|
wolffd@0
|
59 num2str(options(18))])
|
|
wolffd@0
|
60 disp('Sampling starts at the point [0 1].')
|
|
wolffd@0
|
61 disp('The new state is accepted if the threshold value is greater than')
|
|
wolffd@0
|
62 disp('a random number between 0 and 1.')
|
|
wolffd@0
|
63 disp(' ')
|
|
wolffd@0
|
64 disp('Negative step numbers indicate samples discarded from the start of the')
|
|
wolffd@0
|
65 disp('chain.')
|
|
wolffd@0
|
66 disp(' ')
|
|
wolffd@0
|
67 disp('Press any key to continue.')
|
|
wolffd@0
|
68 pause
|
|
wolffd@0
|
69
|
|
wolffd@0
|
70 [samples, energies] = hmc('dempot', x, options, 'demgpot', mix);
|
|
wolffd@0
|
71
|
|
wolffd@0
|
72 disp(' ')
|
|
wolffd@0
|
73 disp('Press any key to continue.')
|
|
wolffd@0
|
74 pause
|
|
wolffd@0
|
75 clc
|
|
wolffd@0
|
76 disp('The plot shows the samples generated by the HMC function.')
|
|
wolffd@0
|
77 disp('The different colours are used to show how the samples move from')
|
|
wolffd@0
|
78 disp('one component to the other over time.')
|
|
wolffd@0
|
79 disp(' ')
|
|
wolffd@0
|
80 disp('Press any key to continue.')
|
|
wolffd@0
|
81 pause
|
|
wolffd@0
|
82 probs = exp(-energies);
|
|
wolffd@0
|
83 fh1 = figure;
|
|
wolffd@0
|
84 % Plot data in 4 groups
|
|
wolffd@0
|
85 ngroups = 4;
|
|
wolffd@0
|
86 g1end = floor(nsamples/ngroups);
|
|
wolffd@0
|
87 g2end = floor(2*nsamples/ngroups);
|
|
wolffd@0
|
88 g3end = floor(3*nsamples/ngroups);
|
|
wolffd@0
|
89 p1 = plot(samples(1:g1end,1), samples(1:g1end,2), 'k.', 'MarkerSize', 12);
|
|
wolffd@0
|
90 hold on
|
|
wolffd@0
|
91 lstrings = char(['Samples 1-' int2str(g1end)], ...
|
|
wolffd@0
|
92 ['Samples ' int2str(g1end+1) '-' int2str(g2end)], ...
|
|
wolffd@0
|
93 ['Samples ' int2str(g2end+1) '-' int2str(g3end)], ...
|
|
wolffd@0
|
94 ['Samples ' int2str(g3end+1) '-' int2str(nsamples)]);
|
|
wolffd@0
|
95 p2 = plot(samples(g1end+1:g2end,1), samples(g1end+1:g2end,2), ...
|
|
wolffd@0
|
96 'r.', 'MarkerSize', 12);
|
|
wolffd@0
|
97 p3 = plot(samples(g2end+1:g3end,1), samples(g2end+1:g3end,2), ...
|
|
wolffd@0
|
98 'g.', 'MarkerSize', 12);
|
|
wolffd@0
|
99 p4 = plot(samples(g3end+1:nsamples,1), samples(g3end+1:nsamples,2), ...
|
|
wolffd@0
|
100 'b.', 'MarkerSize', 12);
|
|
wolffd@0
|
101 legend([p1 p2 p3 p4], lstrings, 2);
|
|
wolffd@0
|
102
|
|
wolffd@0
|
103 clc
|
|
wolffd@0
|
104 disp('We now fit a Gaussian mixture model to the sampled data.')
|
|
wolffd@0
|
105 disp('The model has spherical covariance structure and the correct')
|
|
wolffd@0
|
106 disp('number of components.')
|
|
wolffd@0
|
107 disp(' ')
|
|
wolffd@0
|
108 disp('Press any key to continue.')
|
|
wolffd@0
|
109 pause
|
|
wolffd@0
|
110 % Fit a mixture model to the sample
|
|
wolffd@0
|
111 newmix = gmm(dim, ncentres, 'spherical');
|
|
wolffd@0
|
112 options = foptions;
|
|
wolffd@0
|
113 options(1) = -1; % Switch off all diagnostics
|
|
wolffd@0
|
114 options(14) = 5; % Just use 5 iterations of k-means in initialisation
|
|
wolffd@0
|
115 % Initialise the model parameters from the samples
|
|
wolffd@0
|
116 newmix = gmminit(newmix, samples, options);
|
|
wolffd@0
|
117
|
|
wolffd@0
|
118 % Set up vector of options for EM trainer
|
|
wolffd@0
|
119 options = zeros(1, 18);
|
|
wolffd@0
|
120 options(1) = 1; % Prints out error values.
|
|
wolffd@0
|
121 options(14) = 15; % Max. Number of iterations.
|
|
wolffd@0
|
122
|
|
wolffd@0
|
123 disp('We now train the model using the EM algorithm for 15 iterations')
|
|
wolffd@0
|
124 disp(' ')
|
|
wolffd@0
|
125 disp('Press any key to continue')
|
|
wolffd@0
|
126 pause
|
|
wolffd@0
|
127 [newmix, options, errlog] = gmmem(newmix, samples, options);
|
|
wolffd@0
|
128
|
|
wolffd@0
|
129 % Print out model
|
|
wolffd@0
|
130 disp(' ')
|
|
wolffd@0
|
131 disp('The trained model has parameters ')
|
|
wolffd@0
|
132 disp(' Priors Centres Variances')
|
|
wolffd@0
|
133 disp([newmix.priors' newmix.centres newmix.covars'])
|
|
wolffd@0
|
134 disp('Note the close correspondence between these parameters and those')
|
|
wolffd@0
|
135 disp('of the distribution used to generate the data')
|
|
wolffd@0
|
136 disp(' ')
|
|
wolffd@0
|
137 disp(' Priors Centres Variances')
|
|
wolffd@0
|
138 disp([mix.priors' mix.centres mix.covars'])
|
|
wolffd@0
|
139 disp(' ')
|
|
wolffd@0
|
140 disp('Press any key to exit')
|
|
wolffd@0
|
141 pause
|
|
wolffd@0
|
142
|
|
wolffd@0
|
143 close(fh1);
|
|
wolffd@0
|
144 clear all;
|
|
wolffd@0
|
145
|
