# Sample a random non-decreasing sequence

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.

• Outputting the number of non-decreasing sequences for n,k might make an interesting challenge on its own, if it hasn't already been done... Commented Nov 17, 2016 at 16:35
• @ETHproductions Not really, it's just a binomial (related to this) Commented Nov 17, 2016 at 16:49
• @Sp3000 Ah, OK. It was a fun challenge for me to figure out how to calculate it efficiently though. Commented Nov 17, 2016 at 16:53
• Your requirement that each sequence have equal probability of being output is not possible to satisfy with most garden variety PRNGs which may have just 32 or 48 bit states. According to Wolfram, there are 535 quintillion 20 element subsequences of 1, ..., 100 (didn't check how many of them are nondecreasing). 2^64 is just 18 quintillion. Commented Nov 18, 2016 at 16:03

## 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.

• This is impressive. Could you add an explanation of the method please?
– user9206
Commented Nov 17, 2016 at 17:05
• Yup, busy now, will explain later.
– xnor
Commented Nov 17, 2016 at 17:36
• I counted 90 bytes... (and you can also import* to save 1 byte)
– Rod
Commented Nov 17, 2016 at 17:39
• @Rod Thanks, I forgot about that.
– xnor
Commented Nov 17, 2016 at 18:20
• 87 bytes Commented Jan 31, 2021 at 11:12

# 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


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

• Would numpy help here?
– user9206
Commented Nov 17, 2016 at 19:46
• @Lembik Good thought, but it looks like you'd have to import numpy.random, which is too long.
– xnor
Commented Nov 17, 2016 at 19:49

## 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--]


# TI-BASIC, 53 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
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.

# 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.


# PHP, 7775 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
• First JavaScript and now PHP. All the best scientific programming languages :) Thank you.
– user9206
Commented Nov 18, 2016 at 11:46
• @Lembik, you're welcome. Mind you, it uses a basic PRNG. Do not use for cryptographic stuff. :) Commented Nov 18, 2016 at 13:28