How to sample uniformly from 2-dim region with sinusoidal shape

Question

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.

Answer 1

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$:

A=pi
tA=trunc(A)
x0=acos(tA/A)
G=function(x){(A*sin(x)-tA*x)/(A*sin(x0)-tA*x0)}
inG=function(u){
dif=1/2;x=u*x0
while (dif>1e-8){
  if (G(x)>u){ x=(x-dif)*(x>dif)}else{
    x=(x+dif);if (x>x0) x=x0}
  dif=dif/2}
  return(x)}

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

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.