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

Question

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?

Answer 1

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[123];
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} *)

Answer 2

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[
                {PadRight[ConstantArray[1, len], 2 len]},
                {PadLeft[ConstantArray[1, len], 2 len]},
                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[123];
list=RandomInteger[100,50];

result = splitTotal[list]; //AbsoluteTiming

Mean /@ result

{0.007443, Null}

{1121/25, 1121/25}

Answer 3

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[1234];
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    *)

Answer 4

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][[1]];
      pos2 = Position[d2, x_ /; (x > m2), 1, 1][[1]];
      ,
      pos1 = Position[d1, x_ /; (x > m1), 1, 1][[1]];
      pos2 = Position[d2, x_ /; (x < m2), 1, 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[1];
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}