I have a pair of continuous random variables $\{X,\,Y\}$ that follow a 2-dim density $f_{XY}$ which is uniform over a region (denoted $\mathcal{R}$) and zero elsewhere. The region $\mathcal{R}$ is the "hill top" of a sinusoidal wave such that

$$\begin{align} 0 &< y < g(x)~ &,&& g(x) &\equiv A\cos x - \lfloor A \rfloor \\ 0 &< x < x_0~ &,&& x_0 &\equiv \cos^{-1}\left( \frac{\lfloor A \rfloor}A\right) \quad \text{where}~ g(x_0) = 0 \end{align}$$

where the amplitude $A > 0$ can be any positive value except integers, and $\lfloor A \rfloor$ is the usual flooring, giving the largest integer (including zero) that is less than $A$.

$X$ and $Y$ are correlated simply because $\mathcal{R}$ is not a rectangle.

Here's My Question:

How to sample uniformly from $\mathcal{R}$ described above?

Currently I'm doing acception/rejection, but as usual it is not very efficient: generate $\{x,y\}$uniformly from the rectangle with width $x_0$ and height $g(0)$, and reject those out side of $\mathcal{R}$ where $y > g(x)$. This rejects about $1/3\ldots$ not too bad, but I wonder if it can be improved.

It's not clear to me if there's a transformation that can decouple $X$ and $Y$.

Another approached I have considered is to discretize the region (for given desired sample size $n$) into regular lattice points, then it becomes equivalent to a discrete sampling from $\{1,2,\ldots, n\}$. However, this feels like cheap shot, appropriate just for a certain kind of computations.

It's been many years since I last used Gibbs sampling so I'm still studying how to implement it.

  • 1
    $\begingroup$ It seems quite easy to find a convex polytope say $\mathcal{P}$ which contains the region and is close to it, e.g. using an enveloppe of the curve $y = g(x)$. Then you could use a software that efficiently samples in $\mathcal{P}$ to get proposals and reject those that fall outside of $\mathcal{R}$. $\endgroup$ – Yves Dec 8 '17 at 15:36
  • 2
    $\begingroup$ Sure, you can improve on the rejection rate--at some cost in extra computation. The extra effort seems unlikely to be worth the tiny additional performance improvement that's possible! $\endgroup$ – whuber Dec 8 '17 at 15:50

Generating a random variable $(X,Y)$ uniformly over that set $\mathcal{R}$ is identical to simulate a random variable $X$ with density the upper bound of $\mathcal{R}$, $g(x)$, and then a Uniform variable $Y\sim\mathcal{U}(0,g(x)$. (In our book we call this equivalence the fundamental theorem of simulation.)

To simulate from $$g(x)\propto(A\cos\{x\}-\lfloor A\rfloor)\mathbb{I}_{(0,x_0(A))}(x)$$one can first attempt to find the associated cdf: for $0<x<x_0(A)$ \begin{align*}G(x)&\propto\int_0^x g(x)\text{d}x=\int_0^x (A\cos\{x\}-\lfloor A\rfloor)\text{d}x\\ &=A\sin\{x\}-\lfloor A\rfloor x\end{align*} but this is not easily invertible, at least analytically. (Numerically, $G$ is an increasing function, hence can be inverted by standard methods.) Note that $$G(x_0)\propto A\sin\{x_0\}-\lfloor A\rfloor x_0$$provides the normalisation constant since $G(x_0)=1$.

For instance, here is an R code that inverts $G$:

while (dif>1e-8){
  if (G(x)>u){ x=(x-dif)*(x>dif)}else{
    x=(x+dif);if (x>x0) x=x0}

By drawing first coordinates by cdf inversions, and second coordinates as $\mathcal{U}(0,g(x))$, one gets uniform points under the curve g: enter image description here

Another approach is to use the upper bound $$\cos x \le 1-\frac{x^2}{2}+\frac{x^4}{4!}$$ and figure out a way to simulate from this polynomial density, a secondary accept-reject step being a solution.

The density $g$ being concave, it can also be sandwiched by piece-wise linear functions going through points $(x,g(x))$, suggesting a Gilks and Wild (1992) ARS algorithm. Without the need for a log-transform.

enter image description here

| cite | improve this answer | |
  • $\begingroup$ This makes a lot of sense. I've had some vague ideas about this, and now I can quote your book! Thanks $\endgroup$ – Lee David Chung Lin Dec 9 '17 at 16:19

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.