wolffd@0: % Compare my code with
wolffd@0: % http://www.media.mit.edu/physics/publications/books/nmm/files/index.html
wolffd@0: %
wolffd@0: % cwm.m
wolffd@0: % (c) Neil Gershenfeld  9/1/97
wolffd@0: % 1D Cluster-Weighted Modeling example
wolffd@0: %
wolffd@0: clear all
wolffd@0: figure;
wolffd@0: seed = 0;
wolffd@0: rand('state', seed);
wolffd@0: randn('state', seed);
wolffd@0: x = (-10:10)';
wolffd@0: y = (x > 0);
wolffd@0: npts = length(x);
wolffd@0: plot(x,y,'+')
wolffd@0: xlabel('x')
wolffd@0: ylabel('y')
wolffd@0: nclusters = 4;
wolffd@0: nplot = 100;
wolffd@0: xplot = 24*(1:nplot)'/nplot - 12;
wolffd@0: 
wolffd@0: mux = 20*rand(1,nclusters) - 10;
wolffd@0: muy = zeros(1,nclusters);
wolffd@0: varx = ones(1,nclusters);
wolffd@0: vary = ones(1,nclusters);
wolffd@0: pc = 1/nclusters * ones(1,nclusters);
wolffd@0: niterations = 5;
wolffd@0: eps = 0.01;
wolffd@0:   
wolffd@0: 
wolffd@0: I = repmat(eye(1,1), [1 1 nclusters]);
wolffd@0: O = repmat(zeros(1,1), [1 1 nclusters]);
wolffd@0: X = x(:)';
wolffd@0: Y = y(:)';
wolffd@0: 
wolffd@0: cwr = cwr_em(X, Y, nclusters, 'muX', mux, 'muY', muy,  'SigmaX', I, ...
wolffd@0: 	     'cov_typeX', 'spherical', 'SigmaY', I, 'cov_typeY', 'spherical', ...
wolffd@0: 	     'priorC', pc, 'weightsY', O,  'create_init_params', 0, ...
wolffd@0: 	     'clamp_weights', 1, 'max_iter', niterations, ...
wolffd@0: 	     'cov_priorX', eps*ones(1,1,nclusters), ...
wolffd@0: 	     'cov_priorY', eps*ones(1,1,nclusters));
wolffd@0: 
wolffd@0: 
wolffd@0: % Gershenfeld's EM code
wolffd@0: for step = 1:niterations
wolffd@0:     pplot = exp(-(kron(xplot,ones(1,nclusters)) ...
wolffd@0: 		  - kron(ones(nplot,1),mux)).^2 ...
wolffd@0: 		./ (2*kron(ones(nplot,1),varx))) ...
wolffd@0: 	    ./ sqrt(2*pi*kron(ones(nplot,1),varx)) ...
wolffd@0: 	    .* kron(ones(nplot,1),pc);
wolffd@0:     plot(xplot,pplot,'k');
wolffd@0:     pause(0);
wolffd@0:     px = exp(-(kron(x,ones(1,nclusters)) ...
wolffd@0: 	       - kron(ones(npts,1),mux)).^2 ...
wolffd@0: 	     ./ (2*kron(ones(npts,1),varx))) ...
wolffd@0: 	 ./ sqrt(2*pi*kron(ones(npts,1),varx));
wolffd@0:     py = exp(-(kron(y,ones(1,nclusters)) ...
wolffd@0: 	       - kron(ones(npts,1),muy)).^2 ...
wolffd@0: 	     ./ (2*kron(ones(npts,1),vary))) ...
wolffd@0: 	 ./ sqrt(2*pi*kron(ones(npts,1),vary));
wolffd@0:     p = px .* py .* kron(ones(npts,1),pc);
wolffd@0:     pp = p ./ kron(sum(p,2),ones(1,nclusters));
wolffd@0:     pc = sum(pp)/npts;
wolffd@0:     yfit = sum(kron(ones(npts,1),muy) .* p,2) ...
wolffd@0: 	   ./ sum(p,2);
wolffd@0:     mux = sum(kron(x,ones(1,nclusters)) .* pp) ...
wolffd@0: 	  ./ (npts*pc);
wolffd@0:     varx = eps + sum((kron(x,ones(1,nclusters)) ...
wolffd@0: 		      - kron(ones(npts,1),mux)).^2 .* pp) ...
wolffd@0: 	   ./ (npts*pc);
wolffd@0:     muy = sum(kron(y,ones(1,nclusters)) .* pp) ...
wolffd@0: 	  ./ (npts*pc);
wolffd@0:     vary = eps + sum((kron(y,ones(1,nclusters)) ...
wolffd@0: 		      - kron(ones(npts,1),muy)).^2 .* pp) ...
wolffd@0: 	   ./ (npts*pc);
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % Check equal
wolffd@0: cwr_pc = cwr.priorC';
wolffd@0: assert(approxeq(cwr_pc, pc))
wolffd@0: cwr_mux = cwr.muX;
wolffd@0: assert(approxeq(mux, cwr_mux))
wolffd@0: cwr_SigmaX = squeeze(cwr.SigmaX)';
wolffd@0: assert(approxeq(varx, cwr_SigmaX))
wolffd@0: cwr_muy = cwr.muY;
wolffd@0: assert(approxeq(muy, cwr_muy))
wolffd@0: cwr_SigmaY = squeeze(cwr.SigmaY)';
wolffd@0: assert(approxeq(vary, cwr_SigmaY))
wolffd@0: 
wolffd@0: 
wolffd@0: % Prediction
wolffd@0: 
wolffd@0: X = xplot(:)';
wolffd@0: [cwr_mu, Sigma, post] = cwr_predict(cwr, X);
wolffd@0: cwr_ystd = squeeze(Sigma)';
wolffd@0: 
wolffd@0: % pplot(t,c)
wolffd@0: pplot = exp(-(kron(xplot,ones(1,nclusters)) ...
wolffd@0:    - kron(ones(nplot,1),mux)).^2 ...
wolffd@0:    ./ (2*kron(ones(nplot,1),varx))) ...
wolffd@0:    ./ sqrt(2*pi*kron(ones(nplot,1),varx)) ...
wolffd@0:    .* kron(ones(nplot,1),pc);
wolffd@0: yplot = sum(kron(ones(nplot,1),muy) .* pplot,2) ...
wolffd@0:    ./ sum(pplot,2);
wolffd@0: ystdplot = sum(kron(ones(nplot,1),(muy.^2+vary)) .* pplot,2) ...
wolffd@0:    ./ sum(pplot,2) - yplot.^2;
wolffd@0: 
wolffd@0: 
wolffd@0: % Check equal
wolffd@0: assert(approxeq(yplot(:)', cwr_mu(:)'))
wolffd@0: assert(approxeq(ystdplot, cwr_ystd))
wolffd@0: assert(approxeq(pplot ./ repmat(sum(pplot,2), 1, nclusters),post') )
wolffd@0: 
wolffd@0: plot(xplot,yplot,'k');
wolffd@0: hold on
wolffd@0: plot(xplot,yplot+ystdplot,'k--');
wolffd@0: plot(xplot,yplot-ystdplot,'k--');
wolffd@0: plot(x,y,'k+');
wolffd@0: axis([-12 12 -1 1.1]);
wolffd@0: plot(xplot,.8*pplot/max(max(pplot))-1,'k')
wolffd@0: hold off
wolffd@0: