The problem is to split a list of m*n
numbers into m
subsets of size n
whose means are as equal as possible. Here is an iterative solution. First, partition the sorted list into n
sublists of size m
, then reverse the even-numbered sublists, and transpose the array to get m
sublists of size n
. This will give a reasonably good initial approximation.
Then iterate over pairs of subsets. For each pair, find the split of the 2n
values that will make their means as equal as possible. Continue until the current split is best for m(m-1)/2
pairs in a row. This will not necessarily find the absolute best subsets, but it will find a very good set. Restarting the iterations from the initial approximation many times, but with the order of the numbers in each subset randomized each time, is a convenient way to increase the chance of finding the absolute best set. (There is no way to guarantee that a set is best, short of checking all possible sets.)
To facilitate repeated randomized starts, I have rewritten my earlier code as a function. The last argument is the number of repetitions. It may be omitted, in which case only one start is used. The returned list is {subsets, subset means, standard deviation of subset means}.
EDIT - This revision addresses again the problem (noted by incognito007) that the results are susceptible to floating-point roundoff error. One possible effect is that the While
can loop forever; this actually happened when the previous code was used to partition list2
with m = 3
, n = 6
. The new code handles that case, but I can't guarantee that such looping can never occur when the data are Real
. On the other hand, if all the data are exact then the arithmetic will be exact and the problem can not occur. One solution, exemplified below, is to scale and round the data, then undo the scaling at the end. Besides being safe, this will usually be faster, too.
funk[data_List?(VectorQ[#,NumericQ]&), m_Integer?Positive, n_Integer?Positive,
r:(_Integer?Positive):1] /; Length@data == m*n :=
Module[{b,c,d,e,f,g,h,i,j,k,L, v = Infinity, w},
d = Transpose@MapAt[Reverse,Partition[Sort@data,m],List/@Range[2,n,2]];
c = Prepend[#,1]&/@Permutations@Join[ConstantArray[1,{n-1}],ConstantArray[-1,{n}]];
i = Flatten[Position[#, 1]]&/@c;
j = Flatten[Position[#,-1]]&/@c;
If[And@@ExactNumberQ/@data, w := c.f,
w := Total[f*c,Method->"CompensatedSummation"]; d = N@d; c = N@Transpose@c];
c = Developer`ToPackedArray@c;
Do[g = h = 1; L = 0; e = RandomSample/@d;
While[L < m(m-1)/2,
If[++h > m, If[++g >= m, g = 1]; h = g+1];
f = Developer`ToPackedArray@Flatten@e[[{g,h}]];
k = Ordering[Abs@w,1][[1]];
If[k > 1, e[[g]] = f[[i[[k]]]]; e[[h]] = f[[j[[k]]]]; L = 0]; L++];
k = Variance@Total[e,{2}];
If[k < v, v = k; b = e], {r}];
{b, Mean/@b, Sqrt@v/n}]
list1 = {7.49, 7.56, 7.98, 8.09, 8.16, 8.21, 8.73, 8.64, 8.68,
8.46, 8.57, 8.29, 9.38, 9.43, 8.95, 9.04, 8.9, 9.07};
list2 = {21.08, 21.18, 21.35, 21.79, 21.92, 22.15, 22.38, 22.48, 22.48,
22.51, 22.64, 22.68, 22.75, 22.8, 23.01, 23.28, 23.5, 23.54};
list3 = {36.33, 37.1, 37.19, 37.31, 37.34, 37.61, 37.88, 38.32, 38.42,
38.9, 39.06, 39.12, 39.14, 39.31, 39.39, 39.41, 39.41, 39.43};
First some results with m = 6
, n = 3
, the parameters in the original request.
Timing@funk[list1,6,3,1000]
{3.97, {{{8.9, 8.57, 8.16}, {8.64, 8.68, 8.29}, {8.73, 7.49, 9.38},
{7.56, 8.95, 9.07}, {9.43, 8.21, 7.98}, {8.09, 8.46, 9.04}},
{8.54333, 8.53667, 8.53333, 8.52667, 8.54, 8.53}, 0.0062361}}
Timing[.01*funk[Round[100*list1],6,3,1000]]
{3.22, {{{8.9, 8.57, 8.16}, {8.64, 8.68, 8.29}, {8.73, 7.49, 9.38},
{7.56, 8.95, 9.07}, {9.43, 8.21, 7.98}, {8.09, 8.46, 9.04}},
{8.54333, 8.53667, 8.53333, 8.52667, 8.54, 8.53}, 0.0062361}}
Timing@funk[list2,6,3,1000]
{4.79, {{{22.75, 23.28, 21.18}, {21.08, 22.68, 23.5}, {21.92, 22.51, 22.8},
{23.54, 22.38, 21.35}, {22.64, 22.48, 22.15}, {23.01, 21.79, 22.48}},
{22.4033, 22.42, 22.41, 22.4233, 22.4233, 22.4267}, 0.0091084}}
Timing[.01*funk[Round[100*list2],6,3,1000]]
{3.61, {{{22.75, 23.28, 21.18}, {21.08, 22.68, 23.5}, {21.92, 22.51, 22.8},
{23.54, 22.38, 21.35}, {22.64, 22.48, 22.15}, {23.01, 21.79, 22.48}},
{22.4033, 22.42, 22.41, 22.4233, 22.4233, 22.4267}, 0.0091084}}
Timing@funk[list3,6,3,1000]
{1.83, {{{36.33, 39.31, 39.41}, {39.12, 37.1, 39.06}, {39.14, 37.19, 38.9},
{39.39, 38.42, 37.31}, {38.32, 37.34, 39.41}, {39.43, 37.61, 37.88}},
{38.35, 38.4267, 38.41, 38.3733, 38.3567, 38.3067}, 0.0433803}}
Timing[.01*funk[Round[100*list3],6,3,1000]]
{1.49, {{{39.12, 39.06, 37.1}, {37.19, 39.14, 38.9}, {36.33, 39.31, 39.41},
{39.39, 38.42, 37.31}, {38.32, 37.34, 39.41}, {39.43, 37.61, 37.88}},
{38.4267, 38.41, 38.35, 38.3733, 38.3567, 38.3067}, 0.0433803}}
Now the other way around: m = 3
, n = 6
.
Timing@funk[list1,3,6,1000]
{9.85, {{{7.49, 9.04, 8.09, 9.38, 8.64, 8.57},
{8.68, 8.46, 8.29, 8.73, 7.98, 9.07},
{9.43, 8.21, 8.95, 8.16, 8.9, 7.56}},
{8.535, 8.535, 8.535}, 0.}}
Timing[.01*funk[Round[100*list1],3,6,1000]]
{2.45, {{{8.21, 8.29, 8.73, 9.43, 7.98, 8.57},
{8.9, 8.09, 8.95, 7.56, 9.07, 8.64},
{9.04, 8.68, 7.49, 8.16, 9.38, 8.46}},
{8.535, 8.535, 8.535}, 0}}
Timing@funk[list2,3,6,1000] (* this used to hang *)
{8.02, {{{22.68, 23.54, 21.35, 21.18, 23.28, 22.48},
{21.92, 22.48, 22.51, 21.08, 23.5, 23.01},
{22.38, 22.75, 21.79, 22.15, 22.64, 22.8}},
{22.4183,22.4167, 22.4183}, 0.00096225}}
Timing[.01*funk[Round[100*list2],3,6,1000]]
{2.03, {{{22.15, 22.51, 23.54, 22.48, 21.35, 22.48},
{22.38, 22.64, 21.79, 23.01, 21.18, 23.5},
{22.75, 21.92, 22.68, 21.08, 22.8, 23.28}},
{22.4183,22.4167, 22.4183}, 0.00096225}}
Timing@funk[list3,3,6,1000]
{9.31, {{{39.43, 37.88, 39.14, 36.33, 38.32, 39.12},
{37.34, 39.06, 39.31, 37.19, 38.9, 38.42},
{37.1, 37.31, 39.39, 39.41, 39.41, 37.61}},
{38.37, 38.37, 38.3717}, 0.00096225}}
Timing[.01*funk[Round[100*list3],3,6,1000]]
{2.52, {{{39.14, 37.88, 39.12, 36.33, 38.32, 39.43},
{39.41, 39.31, 37.61, 37.19, 37.31, 39.39},
{37.1, 37.34, 38.9, 38.42, 39.06, 39.41}},
{38.37, 38.37, 38.3717}, 0.00096225}}
aux = {#, Variance[Mean[Transpose[#]]]} & /@ (Partition[#, 2, 2] & /@ Permutations[lst, {6}])
andSortBy[aux, #[[2]] &][[1]]
. As hinted by @PlatoManiac, you need a better way for large lists. $\endgroup$Transpose[FindClusters[list, 3]]
For other data it can produce ragged uneven subdivision, so there should be an even-out-and-sort steps included. $\endgroup$