// genera k centroides C para particionar X (matriz de vectores columna)
function C = kmeans(X,k)
  // inicializa centroides
  Cn = rand(size(X,'r'),k)  // k centroides en el rango [0,1]
  Cn = Cn.*((max(X,'c')-min(X,'c'))*ones(1,k))  // se expanden
  Cn = Cn + min(X,'c')*ones(1,k)  // se posicionan en el espacio
  C = Cn
  // iteracion
  while %t
  disp(C)
    D = distancias(C,X)  // matriz de ||C|| x ||X||
    [d,i] = min(D,'r')  // encuentra valor minimo e indice de la columna
    for j=1:k  // para cada centroide j
      r = find(i==j)  // casos que caen
      Cn(:,j) = sum(X(:,r),'c')/size(r,'c')  // se actualiza centroide
    end
    if and(C==Cn) then break; end  // condicion de termino
    C = Cn
  end

function D = distancias(C,X)
  k = size(C,'c')
  n = size(X,'c')
  D = zeros(k,n)
  for i=1:k
    D(i,:) = sum((X-C(:,i)*ones(1,n)).^2,'r')
  end

function kmtest()
  X = rand(2,300)-rand(2,300)  // desviaciones aleatorias
  // primer cumulo
  X(1,1:100) = X(1,1:100) + rand()*10
  X(2,1:100) = X(2,1:100) + rand()*10
  // segundo cumulo
  X(1,101:200) = X(1,101:200) + rand()*10
  X(2,101:200) = X(2,101:200) + rand()*10
  // tercer cumulo
  X(1,201:300) = X(1,201:300) + rand()*10
  X(2,201:300) = X(2,201:300) + rand()*10
  plot2d(X(1,:),X(2,:),style=.1)
  C = kmeans(X,3)
  plot2d(C(1,:),C(2,:),style=-5)
  [X,C] = return(X,C)