|
wolffd@0
|
1 %DEMGMM5 Demonstrate density modelling with a PPCA mixture model.
|
|
wolffd@0
|
2 %
|
|
wolffd@0
|
3 % Description
|
|
wolffd@0
|
4 % The problem consists of modelling data generated by a mixture of
|
|
wolffd@0
|
5 % three Gaussians in 2 dimensions with a mixture model using full
|
|
wolffd@0
|
6 % covariance matrices. The priors are 0.3, 0.5 and 0.2; the centres
|
|
wolffd@0
|
7 % are (2, 3.5), (0, 0) and (0,2); the variances are (0.16, 0.64) axis
|
|
wolffd@0
|
8 % aligned, (0.25, 1) rotated by 30 degrees and the identity matrix. The
|
|
wolffd@0
|
9 % first figure contains a scatter plot of the data.
|
|
wolffd@0
|
10 %
|
|
wolffd@0
|
11 % A mixture model with three one-dimensional PPCA components is trained
|
|
wolffd@0
|
12 % using EM. The parameter vector is printed before training and after
|
|
wolffd@0
|
13 % training. The parameter vector consists of priors (the column), and
|
|
wolffd@0
|
14 % centres (given as (x, y) pairs as the next two columns).
|
|
wolffd@0
|
15 %
|
|
wolffd@0
|
16 % The second figure is a 3 dimensional view of the density function,
|
|
wolffd@0
|
17 % while the third shows the axes of the 1-standard deviation ellipses
|
|
wolffd@0
|
18 % for the three components of the mixture model together with the one
|
|
wolffd@0
|
19 % standard deviation along the principal component of each mixture
|
|
wolffd@0
|
20 % model component.
|
|
wolffd@0
|
21 %
|
|
wolffd@0
|
22 % See also
|
|
wolffd@0
|
23 % GMM, GMMINIT, GMMEM, GMMPROB, PPCA
|
|
wolffd@0
|
24 %
|
|
wolffd@0
|
25
|
|
wolffd@0
|
26 % Copyright (c) Ian T Nabney (1996-2001)
|
|
wolffd@0
|
27
|
|
wolffd@0
|
28
|
|
wolffd@0
|
29 ndata = 500;
|
|
wolffd@0
|
30 data = randn(ndata, 2);
|
|
wolffd@0
|
31 prior = [0.3 0.5 0.2];
|
|
wolffd@0
|
32 % Mixture model swaps clusters 1 and 3
|
|
wolffd@0
|
33 datap = [0.2 0.5 0.3];
|
|
wolffd@0
|
34 datac = [0 2; 0 0; 2 3.5];
|
|
wolffd@0
|
35 datacov = repmat(eye(2), [1 1 3]);
|
|
wolffd@0
|
36 data1 = data(1:prior(1)*ndata,:);
|
|
wolffd@0
|
37 data2 = data(prior(1)*ndata+1:(prior(2)+prior(1))*ndata, :);
|
|
wolffd@0
|
38 data3 = data((prior(1)+prior(2))*ndata +1:ndata, :);
|
|
wolffd@0
|
39
|
|
wolffd@0
|
40 % First cluster has axis aligned variance and centre (2, 3.5)
|
|
wolffd@0
|
41 data1(:, 1) = data1(:, 1)*0.1 + 2.0;
|
|
wolffd@0
|
42 data1(:, 2) = data1(:, 2)*0.8 + 3.5;
|
|
wolffd@0
|
43 datacov(:, :, 3) = [0.1*0.1 0; 0 0.8*0.8];
|
|
wolffd@0
|
44
|
|
wolffd@0
|
45 % Second cluster has variance axes rotated by 30 degrees and centre (0, 0)
|
|
wolffd@0
|
46 rotn = [cos(pi/6) -sin(pi/6); sin(pi/6) cos(pi/6)];
|
|
wolffd@0
|
47 data2(:,1) = data2(:, 1)*0.2;
|
|
wolffd@0
|
48 data2 = data2*rotn;
|
|
wolffd@0
|
49 datacov(:, :, 2) = rotn' * [0.04 0; 0 1] * rotn;
|
|
wolffd@0
|
50
|
|
wolffd@0
|
51 % Third cluster is at (0,2)
|
|
wolffd@0
|
52 data3(:, 2) = data3(:, 2)*0.1;
|
|
wolffd@0
|
53 data3 = data3 + repmat([0 2], prior(3)*ndata, 1);
|
|
wolffd@0
|
54
|
|
wolffd@0
|
55 % Put the dataset together again
|
|
wolffd@0
|
56 data = [data1; data2; data3];
|
|
wolffd@0
|
57
|
|
wolffd@0
|
58 ndata = 100; % Number of data points.
|
|
wolffd@0
|
59 noise = 0.2; % Standard deviation of noise distribution.
|
|
wolffd@0
|
60 x = [0:1/(2*(ndata - 1)):0.5]';
|
|
wolffd@0
|
61 randn('state', 1);
|
|
wolffd@0
|
62 rand('state', 1);
|
|
wolffd@0
|
63 t = sin(2*pi*x) + noise*randn(ndata, 1);
|
|
wolffd@0
|
64
|
|
wolffd@0
|
65 % Fit three one-dimensional PPCA models
|
|
wolffd@0
|
66 ncentres = 3;
|
|
wolffd@0
|
67 ppca_dim = 1;
|
|
wolffd@0
|
68
|
|
wolffd@0
|
69 clc
|
|
wolffd@0
|
70 disp('This demonstration illustrates the use of a Gaussian mixture model')
|
|
wolffd@0
|
71 disp('with a probabilistic PCA covariance structure to approximate the')
|
|
wolffd@0
|
72 disp('unconditional probability density of data in a two-dimensional space.')
|
|
wolffd@0
|
73 disp('We begin by generating the data from a mixture of three Gaussians and')
|
|
wolffd@0
|
74 disp('plotting it.')
|
|
wolffd@0
|
75 disp(' ')
|
|
wolffd@0
|
76 disp('The first cluster has axis aligned variance and centre (0, 2).')
|
|
wolffd@0
|
77 disp('The variance parallel to the x-axis is significantly greater')
|
|
wolffd@0
|
78 disp('than that parallel to the y-axis.')
|
|
wolffd@0
|
79 disp('The second cluster has variance axes rotated by 30 degrees')
|
|
wolffd@0
|
80 disp('and centre (0, 0). The third cluster has significant variance')
|
|
wolffd@0
|
81 disp('parallel to the y-axis and centre (2, 3.5).')
|
|
wolffd@0
|
82 disp(' ')
|
|
wolffd@0
|
83 disp('Press any key to continue.')
|
|
wolffd@0
|
84 pause
|
|
wolffd@0
|
85
|
|
wolffd@0
|
86 fh1 = figure;
|
|
wolffd@0
|
87 plot(data(:, 1), data(:, 2), 'o')
|
|
wolffd@0
|
88 set(gca, 'Box', 'on')
|
|
wolffd@0
|
89 axis equal
|
|
wolffd@0
|
90 hold on
|
|
wolffd@0
|
91
|
|
wolffd@0
|
92 mix = gmm(2, ncentres, 'ppca', ppca_dim);
|
|
wolffd@0
|
93 options = foptions;
|
|
wolffd@0
|
94 options(14) = 10;
|
|
wolffd@0
|
95 options(1) = -1; % Switch off all warnings
|
|
wolffd@0
|
96
|
|
wolffd@0
|
97 % Just use 10 iterations of k-means in initialisation
|
|
wolffd@0
|
98 % Initialise the model parameters from the data
|
|
wolffd@0
|
99 mix = gmminit(mix, data, options);
|
|
wolffd@0
|
100 disp('The mixture model has three components with 1-dimensional')
|
|
wolffd@0
|
101 disp('PPCA subspaces. The model parameters after initialisation using')
|
|
wolffd@0
|
102 disp('the k-means algorithm are as follows')
|
|
wolffd@0
|
103 disp(' Priors Centres')
|
|
wolffd@0
|
104 disp([mix.priors' mix.centres])
|
|
wolffd@0
|
105 disp(' ')
|
|
wolffd@0
|
106 disp('Press any key to continue')
|
|
wolffd@0
|
107 pause
|
|
wolffd@0
|
108
|
|
wolffd@0
|
109 options(1) = 1; % Prints out error values.
|
|
wolffd@0
|
110 options(14) = 30; % Number of iterations.
|
|
wolffd@0
|
111
|
|
wolffd@0
|
112 disp('We now train the model using the EM algorithm for up to 30 iterations.')
|
|
wolffd@0
|
113 disp(' ')
|
|
wolffd@0
|
114 disp('Press any key to continue.')
|
|
wolffd@0
|
115 pause
|
|
wolffd@0
|
116
|
|
wolffd@0
|
117 [mix, options, errlog] = gmmem(mix, data, options);
|
|
wolffd@0
|
118 disp('The trained model has priors and centres:')
|
|
wolffd@0
|
119 disp(' Priors Centres')
|
|
wolffd@0
|
120 disp([mix.priors' mix.centres])
|
|
wolffd@0
|
121
|
|
wolffd@0
|
122 % Now plot the result
|
|
wolffd@0
|
123 for i = 1:ncentres
|
|
wolffd@0
|
124 % Plot the PC vectors
|
|
wolffd@0
|
125 v = mix.U(:,:,i);
|
|
wolffd@0
|
126 start=mix.centres(i,:)-sqrt(mix.lambda(i))*(v');
|
|
wolffd@0
|
127 endpt=mix.centres(i,:)+sqrt(mix.lambda(i))*(v');
|
|
wolffd@0
|
128 linex = [start(1) endpt(1)];
|
|
wolffd@0
|
129 liney = [start(2) endpt(2)];
|
|
wolffd@0
|
130 line(linex, liney, 'Color', 'k', 'LineWidth', 3)
|
|
wolffd@0
|
131 % Plot ellipses of one standard deviation
|
|
wolffd@0
|
132 theta = 0:0.02:2*pi;
|
|
wolffd@0
|
133 x = sqrt(mix.lambda(i))*cos(theta);
|
|
wolffd@0
|
134 y = sqrt(mix.covars(i))*sin(theta);
|
|
wolffd@0
|
135 % Rotate ellipse axes
|
|
wolffd@0
|
136 rot_matrix = [v(1) -v(2); v(2) v(1)];
|
|
wolffd@0
|
137 ellipse = (rot_matrix*([x; y]))';
|
|
wolffd@0
|
138 % Adjust centre
|
|
wolffd@0
|
139 ellipse = ellipse + ones(length(theta), 1)*mix.centres(i,:);
|
|
wolffd@0
|
140 plot(ellipse(:,1), ellipse(:,2), 'r-')
|
|
wolffd@0
|
141 end
|
|
wolffd@0
|
142
|
|
wolffd@0
|
143 disp(' ')
|
|
wolffd@0
|
144 disp('Press any key to exit')
|
|
wolffd@0
|
145 pause
|
|
wolffd@0
|
146 close (fh1);
|
|
wolffd@0
|
147 clear all; |
