# Identifying subsets with minimal standard deviation

Consider the following:

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

DataAll=StandardDeviation[data];
DataMost=StandardDeviation[Most@data];
DataRest=StandardDeviation[Rest@data];

DataOtherCombinations=
Table[
Sqrt[
(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.)

-
rather than "order must not be changed", I think you mean the numbers must be sequential in the original list... –  rm -rf Apr 16 '12 at 22:43

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]]]];


Usage

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}}}
*)

-
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}}
*)

-
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.

-

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).

-