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

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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]

Follow up

After OP clarifications it is evident that in such a method every transition between the neighbouring tree vertices is equiprobable. To deviate from this assumption the whole tree needs to be traversed keeping account of probabilities. Let us see what can be done here ... (continuation follows)