# Efficiently generating samples from an urn with maximum per element constraint?

I need to sample from a distribution that is a hybrid of uniform and hypergeometric in the sense that all elements are sampled uniformly until an element reaches some specified maximum observations, upon which that element will no longer be sampled.

Think of it as an urn with some $N$ distinctly colored balls, and I sample with replacement $M$ times, but any time a ball has been observed $O$ times in a sample, it is not replaced.

A quick-n-dirty implementation is:

genoutcomes[src_, max_, len_] :=
Module[{oc = {}, cnt = ConstantArray[max, Length@src], rc, r = Range@Length@src},
Do[oc = {oc, src[[rc = RandomChoice[Unitize@cnt -> r]]]};
cnt[[rc]]--;, len];
Flatten@oc];


For example, if the urn has 3 differently colored balls, I want to sample 6, with a maximum observation limit for any ball of 2,

genoutcomes[{red, green, blue}, 2, 6]


produces such a sample.

I'd like to get a more efficient mechanism for doing this (pure Mathematica - I can write it in C or Lisp or whatever, but I'm not after that), ideally for multiple samples in bulk vs calling a routine once per sample.

The above is as I said a quick-n-dirty sketch, I'm pondering but it's late and Morpheus is calling, figured it's an interesting puzzle to throw out here....

Edit: From interactions with posters, here's a quick sanity check for testing.

Generating a large batch of results for 3 elements, 5 max, 9 sample length (something like res = Table[genoutcomes[Range@3, 5, 9], 1000000] for my routine, substitute yours), the result of

Count[res, #] & /@ {{x_, x_, x_, x_, x_, ___}, {___, x_, x_, x_, x_, x_}}


should have a ratio of ~5:1.

• You made some comments below demonstrating that proposed solutions are wrong. Could you please post the code that you use for the verification and demonstrate that your own solution passes the test (CC @bobbym) – yarchik Jun 17 '16 at 10:25
• @yarchik - no need for code, simply generate a large sample of results, using e.g. Table and check the distribution. – ciao Jun 17 '16 at 10:43
• It can be I misunderstand your question, for instance why the sum of probabilities should be not unity, but 1/27 + 2/81 + 4/243=19/243. If you would like to get help, a more precise formulation would be appreciated. – yarchik Jun 17 '16 at 18:40
• @yarchik That should have probabilities 1/27, 2/81, 4/243, 8/729, 193/17496 respectively for a state with a 4-run (4 consecutive of the same element) starting at positions 1,2,3,4, and 5 respectively. In other words, {red,red,red,red,blue,blue,green,green} is much more probable (>5x) than {blue,blue,green,green,red,red,red,red}. Your routine does not, for obvious reasons, produce this distribution. Were you to remove the DeleteDuplicates, it gets closer with larger samples, but then why bother, better to just generate them in the first place. Make no mistake, I appreciate your efforts. – ciao Jun 17 '16 at 20:35
• @Mr.Wizard Certainly not. Project ciao. – ciao Jun 18 '16 at 18:58

This approach is quite different from my previous answer. Therefore I decided to separate them for clarity.

Let us write first a function that for a given sample returns the probability to obtain it via the random process described in the original post.

prob[ncol_, max_, len_, conf_] :=
Module[{cnt = ConstantArray[max, ncol], ucnt, p = 1},
Do[
ucnt = Unitize@cnt;
p = p ucnt[[conf[[i]]]]/Total[ucnt];
cnt[[conf[[i]]]]--;, {i, len}];
p];


Taking example given above we generate 10 samples:

a = Table[genoutcomes[Range@3, 5, 9], {10}]
(*{{2, 1, 3, 3, 1, 3, 2, 2, 3}, {2, 2, 1, 3, 3, 3, 1, 2, 1}, {3, 3, 1, 1, 2, 2, 1, 3, 2}, {3, 3, 3, 3, 3, 1, 1, 1, 1}, {2, 1, 3, 3, 3, 3, 2, 1, 3}, {2, 3, 3, 3, 3, 3, 1, 2, 1}, {3, 3, 3, 1, 3, 3, 2, 1, 2}, {1, 2, 2, 3, 1, 2, 2, 3, 3}, {3, 2, 3, 1, 3, 3, 3, 2, 2}, {3, 2, 2, 2, 1, 3, 3, 2, 1}}*)


and compute the corresponding probabilities

prob[3, 5, 9, #] & /@ a
(*{1/19683, 1/19683, 1/19683, 1/3888, 1/19683, 1/5832, 1/5832, 1/19683, 1/8748, 1/19683}*)


So far so good. For small max and len one can now generate all possible configurations and assign the probabilities, and the problem is solved. But what if the number of configurations is huge? I guess for their brute force generation the function in the OP is almost unbeatable. But we rarely need just configurations. Most probably the ultimate goal would be to compute or compare probabilities of certain classes of configurations. Here, I can provide an analytical result.

For illustration, let us consider the example above. For max=5 and len=9 we want to analytically compute the probabilities of following outcomes:

{{x_, x_, x_, x_, x_, ___}, {___, x_, x_, x_, x_, x_}}


This can be accomplished with the following code:

p1 = prob[3, 5, 9, #] & /@ (Join[{1, 1, 1, 1, 1}, #] & /@
Tuples[Range[3], 4]) // Total
(*1/243*)
p2 = prob[3, 5, 9, #] & /@ (Join[#, {1, 1, 1, 1, 1}] & /@
Tuples[Range[3], 4]) // Total
(*16/19683*)


The ratio of probabilities is as requested

p1/p2
(*81/16*)


i.e., close to 5.

• @Mr.Wizard Finally I found time :) – yarchik Jul 11 '16 at 19:39

I did

genoutcomes[src_, numberforeach_, draws_] :=
RandomSample[Table[src, {numberforeach}] // Flatten, draws]

genoutcomes[{red, green, blue}, 10, 30]


yields

{red, blue, green, red, red, red, green, blue, green, blue, red, red, \ red, green, green, red, green, red, green, blue, blue, green, blue, \ green, blue, red, blue, green, blue, blue}

Hope I did not miss the point of the exercise.

• This does not generate a proper distribution. For example, genoutcomesLSC[{1,2,3},4,6] should generate a 4-run for leading/middle/end of sample with probabilities 1/27, 2/81, 4/243 respectively. This routine makes such runs equiprobable. Don't fret, I've been pondering this, as have some seriously sharp posters here, this is one of those simple-but-hard problems. – ciao Jun 17 '16 at 6:46

Three observations:

• Your random process is non-Markovian, i.e. depends on the previous steps

• The number of states is finite, in your example there are only 90 of them.

• In view of (i) and (ii) you have a random walk on a tree.

You can build the tree explicitly for small systems and study it. But you can also use a brute force approach based on your existing method.

1) generate all states

 config = Table[genoutcomes[{red, green, blue}, 2, 6], {1000}] //
DeleteDuplicates


2) Verify that the length is as expected, (notice a rather large sampling ---1000):

Length[config]
(*90*)


3) Generate any big number (let say 4) of samples

RandomChoice[config,4]