Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Consider the following:

data = {10, 11, 15, 12, 9, 8, 7, 6, 8};


                 (StandardDeviation@data[[1 ;; -i]])^2
                +(StandardDeviation@data[[-i + 1 ;; -1]])^2
    {i, 3, 8}]

ListPlot[Join[{DataAll,DataMost,DataRest},DataOtherCombinations], Joined -> True]

The plot shows, that the combination {{10, 11, 15, 12},{9, 8, 7, 6, 8}} (i=6 in DataOtherCombinations) has the lowest total standard deviation. (I define the total standard deviation as Sqrt[StandardDeviation[a]^2+StandardDeviation[b]^2]).

I was hoping to find the same result with FindClusters but this is not the case:

In[1455]:= FindClusters[data,DistanceFunction->EuclideanDistance]
Out[1455]= {{10, 11, 15, 12, 9, 8, 7, 6, 8}}

Maybe there is a way of customising FindClusters or even better: an approach based on the minimisation of the total standard deviation exists already.

By the way: DataOtherCombinations does not represent all possible combinations. I'm looking for all subsets of data with two restrictions: the order of numbers must not be changed (i.e. the subset must be a sequence of data) and the minimum length of one sequence must be greater than or equal to 3. Hence, another possible combination would be {{10, 11, 15},{12, 9, 8},{7, 6, 8}}. (I tried it with Subsets but this function seems more to return all permutations of all subsets.)

share|improve this question
rather than "order must not be changed", I think you mean the numbers must be sequential in the original list... – R. M. Apr 16 '12 at 22:43
up vote 3 down vote accepted

A general solution:

getMinStdDev[data_?ListQ] := Module[{ip, fp, allowedParts},
  ip = IntegerPartitions[Length@data, Length@data, Range[3, Length@data]];
  fp = Flatten[Permutations /@ ip, 1];
  allowedParts = (data[[#[[1]] ;; #[[2]]]] & /@ 
                  Thread[List[(Accumulate@Join[{1}, #])[[;; -2]],Accumulate@#]]) & /@ fp;
  allowedParts[[Ordering[Total /@ Map[StandardDeviation[#]^2 &, allowedParts, {2}], 1]]]];


getMinStdDev@{10, 11, 10, 15, 14, 16, 9, 8, 7, 6, 8, 7}
getMinStdDev@{10, 11, 15, 12, 9, 8, 7, 6, 8} (*Yours*)
getMinStdDev@Join[RandomInteger[{1,5},5],RandomInteger[{5,10},5], RandomInteger[{10,15},5]]
{{{10, 11, 10}, {15, 14, 16}, {9, 8, 7, 6, 8, 7}}}
{{{10, 11, 15, 12}, {9, 8, 7, 6, 8}}}
{{{3, 3, 5, 3, 2}, {9, 9, 5, 8, 8, 11}, {13, 15, 15, 13}}}
share|improve this answer
Your approach does work, but now I would like to add one restriction: in addition to the total StandardDeviation of allowedPart, the StandardDeviation of each element of allowedPart must be smaller than the StandardDeviation of the data_. I suppose that one needs to implement this additional restriction with the &&-operator(?) – John Apr 17 '12 at 14:52

To reproduce your result example, just specify the number of desired clusters:

data = {10, 11, 15, 12, 9, 8, 7, 6, 8}; 
FindClusters[data, 2]
-> {{10, 11, 15, 12}, {9, 8, 7, 6, 8}}
share|improve this answer
That won't help, as your approach is limited to two clusters for each dataset. I have datasets for which more than two clusters must be identified due to the length of the dataset and my restrictions. For instance: {12,13,15,14,5,4,3,4,5,23,24,22,25}->{{12,13,15,14},{5,4,3,4,5},{23,24,22,25}‌​}. – John Apr 16 '12 at 23:24

You can try to find all the possible lengths of clusters

clusters = Map[FindClusters[data, #] &, Range[2, Length[data] - 1]]

Then remove the sets which have clusters with only 1 or 2 elements

selectclusters = Select[clusters, Min[Length /@ #] >= 3 &]

Then its simply a problem of selecting the set with the minimal stdev

stdevs = Sqrt /@ Total /@ Map[StandardDeviation, selectclusters, {2}]^2
Extract[selectclusters, Position[stdevs, Min@stdevs]]

This works for both the examples you posted. Hope this helps.

edit: There are certain datasets (ex. {19, 10, 7, 1, 3, 9, 10, 19, 5, 9}) Where no clusters with all sets having greater than 3 elements are found so this isn't completely foolproof.

share|improve this answer

First off, it's cleaner to minimize the variances (squares of the standard deviations). Take the square root at the very end if you must.

Second, be aware that Mathematica's Variance (and StandardDeviation) are the so-called "unbiased" calculations used for statistical estimation. They might not be wanted here, because they penalize small clusters (pretty heavily), driving the solution to contain one large cluster in many cases. This suggests implementing a solution in which the objective function can be easily modified.

A recursive solution naturally suggests itself for clustering a list of $n$ values: for each initial interval of length $j\ge 3$, find the optimal clustering for the remaining $n-j$ values (provided there are $3$ or more of them). But don't forget to consider a single cluster! Without any optimization that leads directly to this:

cluster::usage = 
    "cluster[x,k,t] returns an interval partition of x in which \
     each interval has at least k >= 2 elements and \
     ('unbiased') variance not exceeding t (if possible).";
cluster[x_List, k_Integer, t_] := Module[{n = Length[x], options},
   options = Prepend[cluster[Take[x, # - n], k, t], Take[x, #]] & /@ Range[k, n - k];
   options = Append[options, {x}];
   options[[First[Ordering[value[#, t] & /@ options, 1]]]]
cluster[x_List, k_Integer, t_] /; Length[x] < 2 k := {x};

This coding style values clarity because that assists debugging, modification, and integration with other code. It also values flexibility. E.g., hard-coding $k=3$ would speed it up about 25% but it is convenient to allow other values of $k$ (if only for testing).

Oh yes, concerning flexibility: the preceding could drive almost any interval-based clustering algorithm whose objective function is local (that is, a sum over the clusters). To complete the present solution, we need to be specific about value, the objective function:

value::usage = "value[x,t] returns the total variance of a list of intervals \
                (or Infinity if any interval has variance exceeding t).";
value[x_List, t_] := Sum[variance[y, t], {y, x}];

variance::usage = "variance[x,t] measures how 'clustered' a set of numbers x is.";
variance[x_List, t_] := With[{v = Variance[x]}, If[v > t, Infinity, v]];

It's fast and easy to experiment with different measures of "clusteredness"; for instance, try replacing variance with a calculation of the sum of squares of residuals, Variance[x](Length[x]-1).

For an example of usage:

solution = cluster[x = RandomReal[{0, 1}, 30], 3, Variance[x]]

This generates 30 random values and requests clustering into a groups of 3 or more with the variance in each group not to exceed the variance of the entire dataset.

Several aspects of this solution seem worthy of comment:

  1. If no solution is possible, it returns a single cluster. You need to check that (by applying value to it).

  2. The running time grows exponentially. It's going to get bad with lists of length 40 or more.

  3. Running time could be improved with some pruning of the search tree (by tracking the value of the best solution so far).

  4. RAM usage is optimized by this recursive approach (compared to generating all possible partitions of x and then checking them all).

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.