Two random variables $X$ and $Y$ are distributed according to \begin{align} f_{XY}(x,y)= \begin{cases} x+y & 0\leq x \leq 1, 0\leq y \leq 1 \\ 0 & otherwise \end{cases} \end{align}

I want to generate rvs of the given function using python. Is there any way to do it ?

  • $\begingroup$ What does "rvs" stand for? $\endgroup$ Apr 4 at 1:59
  • $\begingroup$ "Random values" $\endgroup$ Apr 5 at 3:52

Rejection sampling is straightforward to implement for this case.

import numpy as np

def rng_xpy(n, rng=None, chunk_size=1024):
    rng = np.random.default_rng(rng)
    rvs = []
    n_drawn = 0
    while n_drawn < n:
        # Draw numbers from U(0, 1) for x, y, and z
        # We draw numbers in chunks for efficiency
        x, y, z = rng.random((3, chunk_size))
        # Scale z to the maximum value of x+y
        z *= 2
        # We only want to use (x, y) values when z is less than x+y
        good_mask = z <= (x + y)
        xy = np.column_stack([x, y])[good_mask]
        n_drawn += len(xy)
    # Because we draw in chunks, we probably drew too many.
    # Assemble the final array and trim it to the required number.
    xy = np.concatenate(rvs, axis=0)[:n]
    return xy
  • $\begingroup$ Thank your for your answer but I dont really understand this part:# Scale z to the maximum value of x+y. can you please elaborate $\endgroup$
    – leo
    Apr 4 at 16:21
  • $\begingroup$ As described in the article, the proposal distribution that we sample (x,y,z) from needs to always be above the target distribution function. We are using a very simple proposal distribution that is just uniformly sampling from the bounding box defined by f(x, y). The maximum value of f(x,y) is at f(1, 1)=1+1=2, so we need to sample z from the U(0,2). Sampling from U(0,1) using the random() method then multiplying by 2 is just a convenience. There are lots of other ways to write that code. $\endgroup$ Apr 5 at 3:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.