title: MeanShift Clustering on dataset
body: <p>I have a numeric dataset and I want to cluster data with a non-parametric algorithm.  Basically, I would like to cluster without specifying the number of clusters for the input.  I am using this code that I accessed through the MathWorks File Exchange network which implements the Mean Shift algorithm.  However, I don't Know how to adapt my data to this code as my dataset has dimensions 516 x 19.</p>
  
  <pre><code>function [clustCent,data2cluster,cluster2dataCell] =MeanShiftCluster(dataPts,bandWidth,plotFlag)
  %UNTITLED2 Summary of this function goes here
  %   Detailed explanation goes here
  
  %perform MeanShift Clustering of data using a flat kernel
  %
  % ---INPUT---
  % dataPts           - input data, (numDim x numPts)
  % bandWidth         - is bandwidth parameter (scalar)
  % plotFlag          - display output if 2 or 3 D    (logical)
  % ---OUTPUT---
  % clustCent         - is locations of cluster centers (numDim x numClust)
  % data2cluster      - for every data point which cluster it belongs to (numPts)
  % cluster2dataCell  - for every cluster which points are in it (numClust)
  % 
  % Bryan Feldman 02/24/06
  % MeanShift first appears in
  % K. Funkunaga and L.D. Hosteler, "The Estimation of the Gradient of a
  % Density Function, with Applications in Pattern Recognition"
  
  
  %*** Check input ****
  if nargin &lt; 2
      error('no bandwidth specified')
  end
  
  if nargin &lt; 3
      plotFlag = true;
      plotFlag = false;
  end
  
  %**** Initialize stuff ***
  %[numPts,numDim] = size(dataPts);
  [numDim,numPts] = size(dataPts);
  numClust        = 0;
  bandSq          = bandWidth^2;
  initPtInds      = 1:numPts
  maxPos          = max(dataPts,[],2);                          %biggest size in each dimension
  minPos          = min(dataPts,[],2);                          %smallest size in each dimension
  boundBox        = maxPos-minPos;                        %bounding box size
  sizeSpace       = norm(boundBox);                       %indicator of size of data space
  stopThresh      = 1e-3*bandWidth;                       %when mean has converged
  clustCent       = [];                                   %center of clust
  beenVisitedFlag = zeros(1,numPts,'uint8');              %track if a points been seen already
  numInitPts      = numPts                               %number of points to posibaly use as initilization points
  clusterVotes    = zeros(1,numPts,'uint16');             %used to resolve conflicts on cluster membership
  
  
  while numInitPts
  
      tempInd         = ceil( (numInitPts-1e-6)*rand)        %pick a random seed point
      stInd           = initPtInds(tempInd)                  %use this point as start of mean
      myMean          = dataPts(:,stInd);                           % intilize mean to this points location
      myMembers       = [];                                   % points that will get added to this cluster                          
      thisClusterVotes    = zeros(1,numPts,'uint16');         %used to resolve conflicts on cluster membership
  
      while 1     %loop untill convergence
  
          sqDistToAll = sum((repmat(myMean,1,numPts) - dataPts).^2);    %dist squared from mean to all points still active
          inInds      = find(sqDistToAll &lt; bandSq);               %points within bandWidth
          thisClusterVotes(inInds) = thisClusterVotes(inInds)+1;  %add a vote for all the in points belonging to this cluster
  
  
          myOldMean   = myMean;                                   %save the old mean
          myMean      = mean(dataPts(:,inInds),2);                %compute the new mean
          myMembers   = [myMembers inInds];                       %add any point within bandWidth to the cluster
          beenVisitedFlag(myMembers) = 1;                         %mark that these points have been visited
  
          %*** plot stuff ****
          if plotFlag
              figure(12345),clf,hold on
              if numDim == 2
                  plot(dataPts(1,:),dataPts(2,:),'.')
                  plot(dataPts(1,myMembers),dataPts(2,myMembers),'ys')
                  plot(myMean(1),myMean(2),'go')
                  plot(myOldMean(1),myOldMean(2),'rd')
                  pause
              end
          end
  
          %**** if mean doesnt move much stop this cluster ***
          if norm(myMean-myOldMean) &lt; stopThresh
  
              %check for merge posibilities
              mergeWith = 0;
              for cN = 1:numClust
                  distToOther = norm(myMean-clustCent(:,cN));     %distance from posible new clust max to old clust max
                  if distToOther &lt; bandWidth/2                    %if its within bandwidth/2 merge new and old
                      mergeWith = cN;
                      break;
                  end
              end
  
  
              if mergeWith &gt; 0    % something to merge
                  clustCent(:,mergeWith)       = 0.5*(myMean+clustCent(:,mergeWith));             %record the max as the mean of the two merged (I know biased twoards new ones)
                  %clustMembsCell{mergeWith}    = unique([clustMembsCell{mergeWith} myMembers]);   %record which points inside 
                  clusterVotes(mergeWith,:)    = clusterVotes(mergeWith,:) + thisClusterVotes;    %add these votes to the merged cluster
              else    %its a new cluster
                  numClust                    = numClust+1                   %increment clusters
                  clustCent(:,numClust)       = myMean;                       %record the mean  
                  %clustMembsCell{numClust}    = myMembers;                    %store my members
                  clusterVotes(numClust,:)    = thisClusterVotes;
              end
  
              break;
          end
  
      end
  
  
      initPtInds      = find(beenVisitedFlag == 0);           %we can initialize with any of the points not yet visited
      numInitPts      = length(initPtInds);                   %number of active points in set
  
  end
  
  [val,data2cluster] = max(clusterVotes,[],1);                %a point belongs to the cluster with the most votes
  
  %*** If they want the cluster2data cell find it for them
  if nargout &gt; 2
      cluster2dataCell = cell(numClust,1);
      for cN = 1:numClust
          myMembers = find(data2cluster == cN);
          cluster2dataCell{cN} = myMembers;
      end
  end
  </code></pre>
  
  <p>This is the test code I am using to try and get the Mean Shift program to work:</p>
  
  <pre><code>clear
  profile on
  
  nPtsPerClust = 250;
  nClust  = 3;
  totalNumPts = nPtsPerClust*nClust;
  m(:,1) = [1 1];
  m(:,2) = [-1 -1];
  m(:,3) = [1 -1];
  var = .6;
  bandwidth = .75;
  clustMed = [];
  %clustCent;
  
  
  x = var*randn(2,nPtsPerClust*nClust);
  %*** build the point set
  for i = 1:nClust
      x(:,1+(i-1)*nPtsPerClust:(i)*nPtsPerClust)       = x(:,1+(i-1)*nPtsPerClust:(i)*nPtsPerClust) + repmat(m(:,i),1,nPtsPerClust);   
  end
  
  tic
  [clustCent,point2cluster,clustMembsCell] = MeanShiftCluster(x,bandwidth);
  toc
  
  numClust = length(clustMembsCell)
  
  
  figure(10),clf,hold on
  cVec = 'bgrcmykbgrcmykbgrcmykbgrcmyk';%, cVec = [cVec cVec];
  for k = 1:min(numClust,length(cVec))
      myMembers = clustMembsCell{k};
      myClustCen = clustCent(:,k);
      plot(x(1,myMembers),x(2,myMembers),[cVec(k) '.'])
      plot(myClustCen(1),myClustCen(2),'o','MarkerEdgeColor','k','MarkerFaceColor',cVec(k), 'MarkerSize',10)
  
  end
  title(['no shifting, numClust:' int2str(numClust)])
  </code></pre>
  
  <p>The test script generates random data <code>X</code>.  In my case. I want to use the matrix <code>D</code> of size 516 x 19 but I am not sure how to adapt my data to this function.  The function is returning results that are not agreeing with my understanding of the algorithm.</p>
  
  <p>Does anyone know how to do this?</p>