Sample a random non-decreasing sequence

Question

Input: Two integers n and k given in any form that is convenient for your code

Output A random non-decreasing sequence of k integers, each in the range 1 to n. The sample should be chosen uniformly from all non-decreasing sequences of k integers with integers in the range 1 to n.

The output can be in any reasonable format you find convenient.

You can use whatever pseudo-random generator your favorite library/language provides.

We can assume that the integers n,k > 0.

Example

Say n,k = 2. The non-decreasing sequences are

1,1
1,2
2,2

Each sequence should have probability 1/3 of being outputted.

Restriction

Your code should run in no more than a few seconds for k = 20 and n = 100.

What doesn't work

If you just sample each integer randomly from the range 1 to n and then sort the list you won't get a uniform distribution.

Answer 1

Actually, 14 12 bytes

This answer is based on Emigna's 05AB1E answer and the answers on this Math.SE question. Golfing suggestions welcome! Try it online!

;;ra+DR╚HS♀-

Ungolfing

      Implicit input n, then k.
;;    Duplicate k twice.
r     Push range [0...k] for later.
a     Invert the stack. Stack: n, k, k, [0...k]
+DR   Push the range [1..n+k-1].
╚     Shuffle the range. Stack: shuffled_range, k, [0...k]
H     Push the first k elements of shuffled_range. Call this increasing.
S     Sort increasing so the elements are actually increasing.
♀-    Subtract each element of [0...k] from each element of increasing.
      This gives us our non-decreasing sequence.
      Implicit return.

Answer 2

Python, 89 bytes

from random import*
lambda n,k:[x-i for i,x in enumerate(sorted(sample(range(1,n+k),k)))]

Generating an increasing sequence rather than a non-decreasing one would be straightforward: this is just a random subset of k numbers between 1 and n, sorted.

But, we can convert an increasing sequence to a non-decreasing one by shrinking each gap between consecutive numbers by 1. So, a gap of 1 becomes a gap of 0, making equal numbers. To do so, decrease the i'th largest value by i

r[0], r[1], ..., r[n-1]  =>  r[0]-0, r[1]-1, ..., r[n-1]-(n-1)

For the result to be from 1 to n, the input must be from 1 to n+k-1. This gives a bijection between non-decreasing sequences of numbers between 1 and n, to increasing sequences between 1 and n+k-1. The same bijection is used in the stars and bars argument for counting such sequences.

The code uses the python function random.sample, which takes k samples without replacement from the input list. Sorting it gives the increasing sequence.

Answer 3

05AB1E, 13 bytes

+<L.r¹£{¹L<-Ä

Try it online!

Explanation

+<L            # range [1 ... n+k-1]
   .r          # scramble order
     ¹£        # take k numbers
       {       # sort
        ¹L<-   # subtract from their 0-based index
            Ä  # absolute value

Answer 4

Python, 87 bytes

from random import*
f=lambda n,k:k>random()*(n+k-1)and f(n,k-1)+[n]or k*[7]and f(n-1,k)

The probability that the maximum possible value n is included equals k/(n+k-1). To include it, put it at the end of the list, and decrement the needed numbers remaining k. To exclude it, decrement the upper bound n. Then, recurse until no more values are needed (k==0).

Python's random doesn't seem to have a built-in for a Bernoulli variable: 1 with some probability, and 0 otherwise. So, this checks if a random value between 0 and 1 generated by random falls below k/(n+k-1). Python 2 would the ratio as float division, so we instead multiply by the denominator: k>random()*(n+k-1).

Answer 5

JavaScript (Firefox 30+), 74 bytes

(n,k,i=0,j=k)=>[for(_ of Array(q=k+n-1))if(Math.random(++i)<k/q--)i-j+k--]

Explanation

xnor's excellent Python answer contains a very good summary of how/why the technique used here works. The first step is to create the range [1, 2, ..., n + k - 1]:

(n,k,i=0)=>[for(_ of Array(q=k+n-1))++i]

Next we need to take k random items from this range. To do this, we need to select each item with probability s / q, where s is the number of items still needed and q is the number of items left in the range. Since we're using an array comprehension, this is fairly easy:

(n,k,i=0)=>[for(_ of Array(q=k+n-1))if(Math.random(++i)<k/q--)k--&&i]

This gives us a uniformly distributed increasing sequence of numbers. This can be fixed by subtracting the number of items j that we have previously found:

(n,k,i=0,j=0)=>[for(_ of Array(q=k+n-1))if(Math.random(++i)<k/q--)k--&&i-j++]

Finally, by storing k in j, we can incorporate k-- into the expression and save a few bytes:

(n,k,i=0,j=k)=>[for(_ of Array(q=k+n-1))if(Math.random(++i)<k/q--)i-j+k--]

Answer 6

TI-BASIC, 54 Bytes

Prompt N,K
K→dim(L1
While K
If rand<K/(N+K-1
Then
N→L1(K
K-1→K
Else
N-1→N
End
End
Disp L1

Follow's xnor's logic, with a small caveat. We could theoretically shave a byte by doing something like this:

K>rand(N+K-1

But rand( is reserved to create a list of random numbers, so we wouldn't be able to do the desired implicit multiplication to save a byte.

This should run decently fast on an 84+ per the restriction but I'll check to make sure when I can.

Answer 7

PHP, 77 75 73 bytes

foreach(array_rand(range(2,$argv[1]+$k=$argv[2]),$k)as$v)echo$v+1-$i++,_;

Run like this:

php -r 'foreach(array_rand(range(2,$argv[1]+$k=$argv[2]),$k)as$v)echo$v+1-$i++,_;' -- 10 5 2>/dev/null;echo
> 1_4_6_9_9_

Explanation

foreach(                    # Iterate over...
  array_rand(               #   a (sorted) random number of items from...
    range(                  #     an array with items...
      2,                    #       from 2
      $argv[1]+$k=$argv[2]  #       to n + k (set arg 2 to $k)
    ),
    $k                      #     Take k number of items (their keys)
  )
  as $v
)
  echo $v +1 - $i++,"_";    # Print the value subtracted by the index.
                            # Need to add 1, because keys are 0-indexed.

Tweaks

• Saved 2 bytes by removing end() call and set $argv[2] to $k instead to shorten access to arguments
• Saved 2 bytes by removing the index from the foreach, since it's simply an incrementing number. Just increment $i each iteration instead