Sample a random non-decreasing sequence

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.

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