parameters = rand(2,4);
means = rand(4,2)*10;

nb_points = 200;
centroids = rand(2,4)*10;

pointsX = normrnd(parameters(1,1),parameters(2,1),1,nb_points/4) .+means(1,1);
pointsY = normrnd(parameters(1,1),parameters(2,1),1,nb_points/4) .+ means(1,2);

for i=2:4
 pointsX = [pointsX , normrnd(parameters(1,i),parameters(2,i),1,nb_points/4) .+means(i,1)];
 pointsY = [pointsY , normrnd(parameters(1,i),parameters(2,i),1,nb_points/4) .+means(i,2)];
endfor

scatter(pointsX,pointsY,'b','filled');
hold on;


y=zeros(nb_points,1);



for cvg = 1:500
  centromeans = zeros(4,2);
  for i=1:nb_points
    min = sqrt( (centroids(1,1)-pointsX(i))^2 + (centroids(2,1)-pointsY(i))^2 );
    y(i) = 1;
    for j=2:4
      argmin = sqrt( (centroids(1,j)-pointsX(i))^2 + (centroids(2,j)-pointsY(i))^2 );
      if argmin < min
        min = argmin;
        y(i) = j;
      endif
    endfor

    centromeans(y(i),1) = centromeans(y(i),1) + pointsX(i);
    centromeans(y(i),2) = centromeans(y(i),2) + pointsY(i);

  endfor

  for i=1:4
    if(sum(y == i) > 0)
      centroids(1,i) = centromeans(i,1)/sum(y == i);
      centroids(2,i) = centromeans(i,2)/sum(y == i);
    endif
  endfor

endfor

scatter(centroids(1,:),centroids(2,:),'r','filled');
hold on;