# Splitting data into two same-size and (approx.) same-mean groups

I need to split a dataset into two same-size and (approx.) same-mean groups. I've looked for a function to address it, but couldn't find it (something similar to SAS PROC OPTEX). Has anyone dealt with something similar using Mathematica? If so, which steps have you followed?

At first I thought your problem might be a duplicate of this one, but it's a bit more subtle because you don't know the means in advance. Here's a compact way you can go about the problem using WFR function "MaximizeOverPermutations" - although here I get it to do a minimization by negating the score. I search for permutations of the original list that, when cut into two equal length halves, gives the lowest absolute mean difference of these halves:

SeedRandom;
list = RandomInteger[100, 50];

(* cuts a list into two halves *)
cut2[list_] := TakeDrop[list, Floor[Length[list]/2]];

(* the score is the absolute difference of means *)
score[list_] := Abs[(Subtract @@ (Mean /@ cut2[list]))]

(* search for a permutation of the list that maximizes the -ve score (minimizes) *)
result = ResourceFunction["MaximizeOverPermutations"][
-score[list[[#]]] &, Length[list],
Method -> {"MonteCarlo", "Iterations" -> 100000}
];

(* the result is a permutation of Range[n], so permute the list first
and then get two halves with cut2 once more *)
{list1, list2} = cut2[ list[[ result[[1, 1]] ]] ];
N@Mean[list1] (* 44.84 *)
N@Mean[list2] (* 44.84 *)


In general, if you want $$n$$ clusters with means as close as possible, you could minimize the maximum absolute difference of means of sub-lists created with "NearEqualPartition":

NearEqualPartition = ResourceFunction["NearEqualPartition"];
MaximizeOverPermutations = ResourceFunction["MaximizeOverPermutations"];

list = {16, 4, 17, 10, 15, 4, 4, 6, 7, 14, 9, 17, 27, 6, 1, 9, 0, 12, 20, 8, 0, 3, 4, 0, 3, 4};
clusters = 4;

score[list_] :=
Max[Abs@*Subtract @@@ Subsets[
Mean /@ NearEqualPartition[list, clusters], {2}]]

result = MaximizeOverPermutations[
-score[list[[#]]] &, Length[list],
Method -> {"MonteCarlo", "Iterations" -> 10000}
];

lists = NearEqualPartition[list[[result[[1, 1]]]], clusters]
(* {{6, 27, 4, 1, 0, 14, 7}, {15, 12, 4, 3, 10, 9, 6}, {9, 4, 0, 17, 4, 17}, {3, 4, 16, 0, 20, 8}} *)

N[Mean /@ lists]
(* {8.42857, 8.42857, 8.5, 8.5} *)

• A quick test with list = Range; gives {4, 10, 3, 5, 6} and {7, 1, 9, 8, 2} with means 5.6 and 5.4. Maybe knock a zero off the 100000 iterations if you want it to be faster. May 19, 2021 at 13:23

You could use LinearProgramming for this. The idea is to create a cost matrix consisting of two sets of your data, one negative and one positive, and the desired minimization vector has half the ones in each part. A simple example to illustrate. Consider the following set of data:

data = {15, 82, 75, 25}


The LinearProgramming call would be:

Quiet @ LinearProgramming[
{15,82,75,25,-15,-82,-75,-25}, (* two copies, one negative and one positive *)
{
{1,1,1,1,0,0,0,0}, (* total number of positive *)
{0,0,0,0,1,1,1,1}, (* total number of negative *)
{1,0,0,0,1,0,0,0}, (* how often 15 is selected *)
{0,1,0,0,0,1,0,0}, (* how often 82 is selected *)
{0,0,1,0,0,0,1,0}, (* how often 75 is selected *)
{0,0,0,1,0,0,0,1}, (* how often 25 is selected *)
{15,82,75,25,-15,-82,-75,-25} (* the total *)
},
{
{2,0}, (* select two positive *)
{2,0}, (* select two negative *)
{1,0}, (* select 15 exactly once *)
{1,0}, (* select 82 exactly once *)
{1,0}, (* select 75 exactly once *)
{1,0}, (* select 25 exactly once *)
{0,1} (* total must be nonnegative *)
},
{
{0,1}, (* each element of the vector is between 0 and 1 *)
{0,1},
{0,1},
{0,1},
{0,1},
{0,1},
{0,1},
{0,1}
},
Integers
]


{0, 0, 1, 1, 1, 1, 0, 0}

So, the positive set is {75, 25} and the negative set is {15, 82}.

The following function creates the desired LinearProgramming call and returns the result:

splitTotal[data_List]:=Module[{r, len=Length@data},
r=Quiet[
LinearProgramming[
Join[data, -data],
Join[
Join[IdentityMatrix[len], IdentityMatrix[len], 2],
{Join[data, -data]}
],
Join[
{{Ceiling[len/2], 0}},
{{Floor[len/2], 0}},
Table[{1, 0}, {len}],
{{0, 1}}
],
Table[{0, 1},{2 len}],
Integers
]
];
{
Pick[data, Take[r, len], 1],
Pick[data, Take[r, len], 0]
}
]


Using @flinty's example:

SeedRandom;
list=RandomInteger[100,50];

result = splitTotal[list]; //AbsoluteTiming

Mean /@ result


{0.007443, Null}

{1121/25, 1121/25}

Here's a slow way to do an exhaustive search. It does not work for lists that are larger than about 30 elements, but it is guaranteed to give the best bipartition and can be used to check the correctness of faster methods.

Use a random list of length 20:

SeedRandom;
list = RandomVariate[NormalDistribution[], 20]
(*    {-0.508336, -0.0706019, -1.5939, 1.53659, 2.67802,
-1.31683, -1.09456, 0.142948, 0.104166, -1.67311,
0.0357514, 0.170186, 0.934682, -0.826864, -0.090963,
1.03341, -0.363529, 0.188591, -0.136142, -0.0839498}    *)


The mean value of the entire list is

μ = Mean[list]
(*    -0.046722    *)


Enumerate all subsets of the list that contain half of the list's elements, and pick the subset whose mean is closest to $$\mu$$:

half = First@MinimalBy[Subsets[list, {Length[list]/2}],
Abs[Mean[#] - μ] &]
(*     {-0.508336, -0.0706019, -1.09456, 0.104166, 0.0357514,
0.170186, 0.934682, -0.090963, 0.188591, -0.136142}      *)


The elements in the other half of list:

otherhalf = Complement[list, half]
(*    {-1.67311, -1.5939, -1.31683, -0.826864, -0.363529,
-0.0839498, 0.142948, 1.03341, 1.53659, 2.67802}      *)


check the means:

Mean[half]
(*    -0.0467231    *)

Mean[otherhalf]
(*    -0.0467209    *)


Here is a proposal how you could proceed.

Calculate the means of both data sets. If the mean m1 of data set d1 is smaller than the mean m2 of data set d2, then exchange an element of d1 that is smaller than m1 against an element of d2 that is larger than m2. If m1>m2 then do the opposite until both means differ no more than a given tolerance.

sameMean[dist10_, dist20_, tol_] :=
Module[{nmax = 1000, dist1 = dist10, dist2 = dist20, swap, d1, d2,
m1, m2, pos1, pos2},
swap[d10_, d20_, m1_, m2_] := (
d1 = d10;
d2 = d20;
If[m1 < m2,
pos1 = Position[d1, x_ /; (x < m1), 1, 1][];
pos2 = Position[d2, x_ /; (x > m2), 1, 1][];
,
pos1 = Position[d1, x_ /; (x > m1), 1, 1][];
pos2 = Position[d2, x_ /; (x < m2), 1, 1][];
];
{d1[[pos1]], d2[[pos2]]} = {d2[[pos2]], d1[[pos1]]};
{d1, d2}
);

While[
(m1 = Mean[dist1];
m2 = Mean[dist2];
Abs[m1 - m2] > tol)
,
{dist1, dist2} = swap[dist1, dist2, m1, m2];
];
{dist1, dist2}
];


Note, for safety, I put a maximum number of exchanges of 1000 in. To test this, we first need 2 data sets with different means:

SeedRandom;
dist1 = RandomVariate[NormalDistribution[0, 1], 100];
dist2 = RandomVariate[NormalDistribution[1, 1], 100];


With this data we get:

Mean /@ {dist1, dist2}
{dist1, dist2} = sameMean[dist1, dist2, 0.1];
Mean /@ {dist1, dist2}


Means start: {-0.00966945, 1.20263}

Means end: {0.564653, 0.628304}

• What you suggest is similar to a zero-temperature version of @flinty's Monte Carlo search. As such I expect it to get stuck in local optima and not reach the global one. May 19, 2021 at 14:28
• To prevent this, you could choose the position to swap by chance or even better, choose the value that is most distant from the current mean. May 19, 2021 at 16:23
• Yes, you're on track to re-invent the Metropolis–Hastings algorithm :-) May 19, 2021 at 17:59
• Do you think I should claim the name for me? :-) May 19, 2021 at 18:10