Sampling two integer lists of given sizes from a distribution that have equal sums

Question

(This question emerged from discussions in this post.)

Context and code sample:

I am trying to figure out if there is a way to generate two $a$ and $b,$ comprised of $n_a$ and $n_b$ integers which are sampled from a given discrete distribution dist, such that the sums of the two lists are equal (i.e. Total@a == Total@b).

For the sake of discussion, the distribution can either be a custom one, such as:

distcustom = {0.17525, 0.329672, 0.2882, 0.14761, 0.04771, 0.0101, 0.001362, \ 0.0001141, 0.0000001}

which can be histogram of the integers, here 1 has a probability 0.17525, 2 has a probability 0.329672 etc.

Or a more common one, such as a Poisson distribution:

dist = Table[N@PDF[PoissonDistribution[3], i], {i, 9}]

{0.149361, 0.224042, 0.224042, 0.168031, 0.100819, 0.0504094, \ 0.021604, 0.00810151, 0.0027005}

Sampling integers from given distribution:

To sample lists of desired number of integers from such dists without imposing the sum condition, in other words creating sequences of integers from a distribution, we can do:

Say we want a list of 10 elements:

For a built-in distribution such as PoissonDistribution we can use:

sequence = {};
sequence = RandomVariate[PoissonDistribution[3], 10]

from an input custom distribution/histogram such as distcustom we can use:

repeat[m_, n_Integer?Positive] := Sequence @@ ConstantArray[m, n]
sequence = {};
(*we generate a long list of integers sampled from distcustom, then later we randomsample from it*)
For[i = 1, i <= Length[distcustom], i++,
  tmpval = IntegerPart[Round[10000*distcustom[[i]]]];
  If[tmpval == 0, tmpval = 1];
  sequence = Join[sequence, {repeat[i, tmpval]}];
  ];

---

Two approaches come to mind:

• Generating the a and b lists separately: We could e.g. generate a long list of integers sampledls sampled according to dist, this could for example be a list of a million integers for higher accuracy. Then to create say a, we keep trying to extract n_a element sublists from sampledls and similarly for list b, until we find two sublists that satisfy Total@a == Total@b.

• Partitioning a larger list into the smaller lists of a and b: We generate a list of $n_a+n_b$ integers sampled from dist, let’s denote by ab, then we try various partitions of ab into two lists of a and b having $n_a$ and $n_b$ elements respectively, until we find a partition which satisfies Total@a == Total@b.

Intuitively, the second one seems to be a more efficient approach (as we sample once from the dist, then the computation boils down to creating partitions/pairs of given sizes).

1. Do any of these approaches seem sound? would an approach similar to the second one be indeed the more efficient one to opt for?

2. Is there possibly a simpler way to go about solving this problem by better exploiting the built-in functionalities of Mathematica? To re-iterate the problem again:

Generating two lists $a$ and $b,$ comprised of given $n_a$ and $n_b$ integers which are sampled from a given discrete distribution dist, such that the sums of the two lists are equal (i.e. Total@a ==Total@b).

Answer 1

st[dist_, l1_, l2_] := Module[{rl1 = RandomVariate[dist, l1], rl2, t},
   t = Tr@rl1;
   While[Tr[rl2 = RandomVariate[dist, l2]] != t];
   {rl1, rl2}];

Usage:

Create two lists of length 10 and 15 respectively from a Poisson distribution, such that sums of lists are same.

{list1, list2} = st[PoissonDistribution[3], 10, 15]; 

{list1, list2}

Tr/@%

{{4,7,3,3,1,3,1,5,4,2},{4,0,3,2,4,5,2,1,1,0,5,1,1,1,3}}

{33, 33}

Per question in comments from OP re: their distcustom, such constructions can be accommodated, among other ways, as (here assuming probabilities over 1, 2, ..., j for length j probability list):

st[EmpiricalDistribution[distcustom -> Range@Length@distcustom], 100, 100]