function cwm(y,nclusters,niter)

% Cluster-weighted modeling algorithm 
% (includes Expectation-Maximization)
% Inputs: 
% 	y = vector of points to fit with a constant local function
%	nclusters = number of clusters
%	niter = number of iterations

% Calculate number of input points
npts = length(y);

% Generate the x axis
x = [1:npts]';
y = y';

% Calculate number of input points
npts = length(x);

% Choose a small epsilon to prevent the variances from shrinking to 0
eps = 0.01;

% Number of points for graphics
nplot = 100;

% Axis for ploting
xplot = npts*(1:nplot)'/nplot;

% Generate random variances
varx = rand(1,nclusters) + eps;
vary = rand(1,nclusters) + eps;

% Generate random distribution of clusters on the X axis
mux = rand(1,nclusters)*npts;

% Generate random coefficients (local models are constant functions) 
fxB = zeros(1,nclusters);

% Give the same weight to each cluster
weights = ones(1,nclusters) / nclusters; 

% Simplify with
NCLUSTERS = ones(1,nclusters);
NPTS = ones(npts,1);

for i=1:niter,

	% Calculate the domaine of influence in the input space (1 dimension)
	X = kron(x,NCLUSTERS);
	MUX = kron(NPTS,mux);
	VARX = kron(NPTS,varx);
	
	inputs = exp(-(X-MUX).^2 ./ (2*VARX)) ./ sqrt(2*pi*VARX);
	
	% Calculate the dependance in the output space
	Y = kron(y,NCLUSTERS);
	FXB = kron(NPTS,fxB);
	VARY = kron(NPTS,vary);
	
	outputs = exp(-(Y-FXB).^2 ./ (2*VARY)) ./ sqrt(2*pi*VARY);
	
	% Calculate posteriors
	WEIGHTS = kron(NPTS,weights);
	
	P = inputs .* outputs .* WEIGHTS;	
	posts = P ./ kron(sum(P,2),NCLUSTERS);

	% Calculate new weights
	weights = sum(posts) / npts;

	% Calculate new mux
	X = kron(x,NCLUSTERS);
	
	mux = sum(X .* posts) ./ (npts*weights);
	
	% Calculate new variance x
	X = kron(x,NCLUSTERS);
	MUX = kron(NPTS,mux);
	
	varx = eps + sum((X - MUX).^2 .* posts) ./ (npts*weights);
	
	% Calculate new coefficients fxB
	Y = kron(y,NCLUSTERS);
	
	fxB = sum(Y .* posts) ./ (npts*weights);

	% Calculate new variance y
	Y = kron(y,NCLUSTERS);
	FXB = kron(NPTS,fxB);
	
	vary = eps + sum((Y - FXB).^2 .* posts) ./ (npts*weights);

end

% Plot
NPLOT = ones(nplot,1);
XPLOT = kron(xplot,NCLUSTERS);
MUX = kron(NPLOT,mux);
VARX = kron(NPLOT,varx);
WEIGHTS = kron(NPLOT,weights);
FXB = kron(NPLOT,fxB);

inputs = exp(-(XPLOT-MUX).^2 ./ (2*VARX)) ./ sqrt(2*pi*VARX) .* WEIGHTS;
Fx = sum(FXB .* inputs, 2) ./ sum(inputs, 2);

plot(xplot,Fx);
hold on
plot(x,y,'k+');
hold off

	
