# Generating random samples from a density function

How does a computer algorithm set up to take as input an arbirary bivariate probability density function, generate pairs of numbers from that distribution? I have found a routine called simcontour that is part of LearnBayes in R that performs that operation.

Denote the joint density you want to sample from by $p_{XY}(x,y)$. Then the marginal densities are easily obtainable:

$$p_{X}(x) = \int_{D_{Y}} p_{XY}(x,y) dy, \ \ \ \ \ \ \ \ \ p_{Y}(y) = \int_{D_{X}} p_{XY}(x,y) dx$$

where $D_X, D_Y$ denote the support of $X, Y$, respectively. From this you can calculate each of the conditional density functions:

$$p_{X|Y=y}(x) = \frac{ p_{XY}(x,y) }{ p_{Y}(y) }, \ \ \ \ \ \ \ \ p_{Y|X=x}(y) = \frac{ p_{XY}(x,y) }{ p_{X}(x) }$$

From here you can execute Gibbs' sampling to sample from $p_{XY}$. That is, start with some arbitrary starting value in the $(x_{0}, y_{0})$. Then, at step $i$,

(1) Sample from $p_{X|Y=y_{i}}(x)$ to obtain $x_{i+1}$

(2) Sample from $p_{Y|X=x_{i+1}}(y)$ to obtain $y_{i+1}$

(3) Repeat a number of times equal to the desired sample size

Each $(x_{i}, y_{i})$ is an approximate draw from the joint distribution of $X$ and $Y$. The conditional densities can be sampled from in steps (1) and (2) using, for example, rejection sampling assuming they do not come from one of the "standard" distributions, in which case you can use a pre-existing function to generate from them.

You can use rejection sampling. Let $f$ be the probability density we want to sample from (it does not matter if it is 1- or 2-dimensional). If we can find a probability distribution $g$, from which we can easily sample from, e.g., a uniform distribution (a good choice, if the support of $f$ is bounded), a Gaussian distribution, or a Cauchy distribution (a last-resort choice, inefficient but useful if $f$ has fat tails), such that there exists a constant $c$ with $f \leq c g$, then we can sample from $f$ as follows. Take a random sample $x$ from $g$, take a random number $\lambda$ uniformly in $[0,1]$; if $\lambda c g(x) < f(x)$, keep the sample, otherwise reject it and try again. Repeat until you have the desired number of samples.

The reason why it works is easily seen on a picture: $(x,\lambda c g(x))$ is a point sampled uniformly in the area under the curve $x \mapsto c g(x)$, and the condition $\lambda c g(x) < f(x)$ asks that it be under the curve $x\mapsto f(x)$.