# Generating Rvs for a given PDf in python

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 ?

• What does "rvs" stand for? – Wolfgang Bangerth Apr 4 at 1:59
• "Random values" – Robert Kern Apr 5 at 3:52

## 1 Answer

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]
rvs.append(xy)
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

• 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 – leo Apr 4 at 16:21
• 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. – Robert Kern Apr 5 at 3:52