|
Daniel@0
|
1 % Verify that my code gives the same results as the 1D example at
|
|
Daniel@0
|
2 % http://www.media.mit.edu/physics/publications/books/nmm/files/cwm.m
|
|
Daniel@0
|
3
|
|
Daniel@0
|
4 seed = 0;
|
|
Daniel@0
|
5 rand('state', seed);
|
|
Daniel@0
|
6 randn('state', seed);
|
|
Daniel@0
|
7 x = (-10:10)';
|
|
Daniel@0
|
8 y = double(x > 0);
|
|
Daniel@0
|
9 npts = length(x);
|
|
Daniel@0
|
10 plot(x,y,'+')
|
|
Daniel@0
|
11
|
|
Daniel@0
|
12 nclusters = 4;
|
|
Daniel@0
|
13 nplot = 100;
|
|
Daniel@0
|
14 xplot = 24*(1:nplot)'/nplot - 12;
|
|
Daniel@0
|
15
|
|
Daniel@0
|
16 mux = 20*rand(1,nclusters) - 10;
|
|
Daniel@0
|
17 muy = zeros(1,nclusters);
|
|
Daniel@0
|
18 varx = ones(1,nclusters);
|
|
Daniel@0
|
19 vary = ones(1,nclusters);
|
|
Daniel@0
|
20 pc = 1/nclusters * ones(1,nclusters);
|
|
Daniel@0
|
21
|
|
Daniel@0
|
22
|
|
Daniel@0
|
23 I = repmat(eye(1,1), [1 1 nclusters]);
|
|
Daniel@0
|
24 O = repmat(zeros(1,1), [1 1 nclusters]);
|
|
Daniel@0
|
25 X = x(:)';
|
|
Daniel@0
|
26 Y = y(:)';
|
|
Daniel@0
|
27
|
|
Daniel@0
|
28 % Do 1 iteration of EM
|
|
Daniel@0
|
29
|
|
Daniel@0
|
30 %cwr = cwr_em(X, Y, nclusters, 'muX', mux, 'muY', muy, 'SigmaX', I, 'cov_typeX', 'spherical', 'SigmaY', I, 'cov_typeY', 'spherical', 'priorC', pc, 'weightsY', O, 'init_params', 0, 'clamp_weights', 1, 'max_iter', 1, 'cov_priorX', zeros(1,1,nclusters), 'cov_priorY', zeros(1,1,nclusters));
|
|
Daniel@0
|
31
|
|
Daniel@0
|
32 cwr = cwr_em(X, Y, nclusters, 'muX', mux, 'muY', muy, 'SigmaX', I, 'cov_typeX', 'spherical', 'SigmaY', I, 'cov_typeY', 'spherical', 'priorC', pc, 'weightsY', O, 'create_init_params', 0, 'clamp_weights', 1, 'max_iter', 1);
|
|
Daniel@0
|
33
|
|
Daniel@0
|
34
|
|
Daniel@0
|
35 % Check this matches Gershenfeld's code
|
|
Daniel@0
|
36
|
|
Daniel@0
|
37 % E step
|
|
Daniel@0
|
38 % px(t,c) = prob(x(t) | c)
|
|
Daniel@0
|
39 px = exp(-(kron(x,ones(1,nclusters)) ...
|
|
Daniel@0
|
40 - kron(ones(npts,1),mux)).^2 ...
|
|
Daniel@0
|
41 ./ (2*kron(ones(npts,1),varx))) ...
|
|
Daniel@0
|
42 ./ sqrt(2*pi*kron(ones(npts,1),varx));
|
|
Daniel@0
|
43 py = exp(-(kron(y,ones(1,nclusters)) ...
|
|
Daniel@0
|
44 - kron(ones(npts,1),muy)).^2 ...
|
|
Daniel@0
|
45 ./ (2*kron(ones(npts,1),vary))) ...
|
|
Daniel@0
|
46 ./ sqrt(2*pi*kron(ones(npts,1),vary));
|
|
Daniel@0
|
47 p = px .* py .* kron(ones(npts,1),pc);
|
|
Daniel@0
|
48 pp = p ./ kron(sum(p,2),ones(1,nclusters));
|
|
Daniel@0
|
49
|
|
Daniel@0
|
50 % M step
|
|
Daniel@0
|
51 eps = 0.01;
|
|
Daniel@0
|
52 pc2 = sum(pp)/npts;
|
|
Daniel@0
|
53
|
|
Daniel@0
|
54 mux2 = sum(kron(x,ones(1,nclusters)) .* pp) ...
|
|
Daniel@0
|
55 ./ (npts*pc2);
|
|
Daniel@0
|
56 varx2 = eps + sum((kron(x,ones(1,nclusters)) ...
|
|
Daniel@0
|
57 - kron(ones(npts,1),mux2)).^2 .* pp) ...
|
|
Daniel@0
|
58 ./ (npts*pc2);
|
|
Daniel@0
|
59 muy2 = sum(kron(y,ones(1,nclusters)) .* pp) ...
|
|
Daniel@0
|
60 ./ (npts*pc2);
|
|
Daniel@0
|
61 vary2 = eps + sum((kron(y,ones(1,nclusters)) ...
|
|
Daniel@0
|
62 - kron(ones(npts,1),muy2)).^2 .* pp) ...
|
|
Daniel@0
|
63 ./ (npts*pc2);
|
|
Daniel@0
|
64
|
|
Daniel@0
|
65
|
|
Daniel@0
|
66 denom = (npts*pc2);
|
|
Daniel@0
|
67 % denom(c) = N*pc(c) = w(c) = sum_t pp(c,t)
|
|
Daniel@0
|
68 % since pc(c) = sum_t pp(c,t) / N
|
|
Daniel@0
|
69
|
|
Daniel@0
|
70 cwr_mux = cwr.muX;
|
|
Daniel@0
|
71 assert(approxeq(mux2, cwr_mux))
|
|
Daniel@0
|
72 cwr_SigmaX = squeeze(cwr.SigmaX)';
|
|
Daniel@0
|
73 assert(approxeq(varx2, cwr_SigmaX))
|
|
Daniel@0
|
74
|
|
Daniel@0
|
75 cwr_muy = cwr.muY;
|
|
Daniel@0
|
76 assert(approxeq(muy2, cwr_muy))
|
|
Daniel@0
|
77 cwr_SigmaY = squeeze(cwr.SigmaY)';
|
|
Daniel@0
|
78 assert(approxeq(vary2, cwr_SigmaY))
|
|
Daniel@0
|
79
|
|
Daniel@0
|
80
|
