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

  • $\begingroup$ What is this for? That is, what is the endgame here - what is then done with these lists? This feels like an x-y question. In any case, what is the maximum cardinality of these lists you want to do this for? $\endgroup$
    – ciao
    Mar 9, 2020 at 17:18

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


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}



{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] 
  • $\begingroup$ Many thanks! This is really simple and neat! Regarding the case where we have a custom probability distribution (such as distcustom example in the post), is there a way we could extend the above module to work with those too? (afaik RandomVariate works best with built-in distributions). If not, would an approach like mine (creating a large list of integers sampled from the dist, and then re-sampling that list) be a possible alternative? Thanks again! $\endgroup$
    – user52181
    Mar 9, 2020 at 20:13
  • 1
    $\begingroup$ @user929304 - see update to answer. $\endgroup$
    – ciao
    Mar 9, 2020 at 20:40
  • $\begingroup$ Thx for the prompt edit! EmpiricalDistribution function, so glad I ask about i! :D After some testing, if I may ask a follow-up question: it is working correctly and also rather speedily, however, as soon as the difference between the length of the lists becomes a bit large, it seems to never be able to find a solution. E.g., st[PoissonDistribution[3], 200, 250] It might be that none in fact exists for those constraints, but would the search in our module benefit from refreshing rl1 as well each round? Since right now, if I'm not mistaken, we initialise rl1 and only change the 2nd. $\endgroup$
    – user52181
    Mar 9, 2020 at 21:30
  • $\begingroup$ @user929304: For your example, it's easy to show the probability of two lists being generated with the same sum is ~2.6 x 10^-6. Given a first list, depending on its sum, the probability of a random second list matching is almost always much less probable, for example the expected sum for list 1 is 600, the probability of the second longer list matching is ~1.6 x 10^-9. Yes, you could choose to generate both and then compare, repeating until a match, but you then severely bias the results. If you're not after pairs that track the distribution, that's a different problem/question to the OP. $\endgroup$
    – ciao
    Mar 9, 2020 at 22:08
  • $\begingroup$ Excellent point, that's clarified the problem for me! Indeed it is crucial to be sampling from the desired distribution. What if, we opted for the 2nd approach? That is, first generating a list of length $n_a+n_b$ sampled from the distribution, then looking for partitions of it into lists a and b that satisfy the length and sum criteria. Do you reckon there's any fundamental difference in this approach? $\endgroup$
    – user52181
    Mar 10, 2020 at 12:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.