I am trying to use Gibbs sampling or Metropolis-Hastings to draw samples from the joint distribution$$f(x,y)\propto\exp(-|x|-|y|-a \cdot |x-y|)$$ For this I need the conditional distributions of $x$ and $y$ ($a$ is known and constant). How do I obtain these? Is there a trick / shortcut?

  • $\begingroup$ No trick involved, just apply standard probability: find $f(x|y)$ and $f(y|x)$ from the joint and then simulate from these densities. $\endgroup$
    – Xi'an
    Feb 4, 2019 at 9:56
  • 1
    $\begingroup$ @Xi'an The full conditionals don't seem to be from an obvious density though or am I missing something? $\endgroup$ Feb 4, 2019 at 10:21
  • $\begingroup$ @Xi'an By rejection sampling I presume? $\endgroup$ Feb 4, 2019 at 12:00
  • $\begingroup$ @Xi'an Could you write down these inverse CDF? $\endgroup$ Feb 4, 2019 at 12:08

3 Answers 3


Disclaimer: although there is nothing to complain about Ben's answer (!), except maybe that the normalising constant of the conditional is not of direct use, here is what I wrote while being off-line, so I may as well post it!

The full conditional of $X$ given $Y$ has a density that is proportional to \begin{align} f(x|y) &\propto \exp\{ -|x|-a|y-x|\}\\ &=\begin{cases} \exp\{ -x-a(y-x)\} &\text{ if } 0\le x\le y\\ \exp\{ x-a(y-x)\} &\text{ if } x\le \min(y,0)\\ \exp\{ -x+a(y-x)\} &\text{ if } x\ge \max(0,y)\\ \exp\{ x+a(y-x)\} &\text{ if } y\le x\le 0\\ \end{cases} \\ &=\begin{cases} \exp\{ -(1-a)x-ay\} &\text{ if } 0\le x\le y\\ \exp\{ (1-a)x-ay\} &\text{ if } x\le \min(y,0)\\ \exp\{ -(a+1)x+ay\} &\text{ if } x\ge \max(0,y)\\ \exp\{ (1-a)x+ay\} &\text{ if } y\le x\le 0\\ \end{cases} \\ \end{align} with one of the two conditions $0\le x\le y$ or $y\le x\le 0$ being obviously empty. Each of the (lhs) terms can be normalised into a (rhs) density: \begin{align}\exp\{ -(1-a)x\} \mathbb{I}_{0\le x\le y} &\propto \dfrac{(a-1)\exp\{ -(1-a)x\}}{\exp\{ (a-1)y\}-1}\mathbb{I}_{0\le x\le y}\\ \exp\{ (1+a)x\} \mathbb{I}_{x\le \min(y,0)} &\propto \dfrac{(1+a)\exp\{ (1+a)x\}}{\exp\{ (1+a)\min(y,0)\}} \mathbb{I}_{x\le \min(y,0)}\\ \exp\{ -(1+a)x\} \mathbb{I}_{x\ge \max(0,y)}&\propto \dfrac{(1+a)\exp\{ (1+a)x\}}{\exp\{ (1+a)\max(0,y)\}}\mathbb{I}_{x\ge \max(0,y)}\\ \exp\{ (1-a)x\} \mathbb{I}_{y\le x\le 0} &\propto \dfrac{(1-a)\exp\{ (1-a)x\}}{1-\exp\{ (1-a)y\}}\mathbb{I}_{y\le x\le 0}\\ \end{align} which means that it is directly possible to simulate from each of these densities, by inverse cdf, picking one of the four (three if accounting fo) with probabilities $$\dfrac{\exp\{-ay\}[\exp\{ -(1-a)y\}-1]}{a-1}\mathbb{I}_{y\ge 0}, \dfrac{\exp\{-ay\}\exp\{ (1+a)\min(y,0)\}}{1+a}, \dfrac{\exp\{ay\}\exp\{ (1+a)\max(0,y)\}}{1+a}, \dfrac{\exp\{ay\}[1-\exp\{ (1-a)y\}]}{1-a}\mathbb{I}_{y\le 0}$$ (Obviously, a basic accept-reject algorithm based on a double exponential $\pm\mathcal{E}(1)$ proposal avoids this painful derivation since, when simulating $X\sim\pm\mathcal{E}(1)$ the accept-reject ratio $$\dfrac{\exp\{ -|x|-a|y-x|\}}{\exp\{-|x|\}}$$ is $\exp\{-a|y-x|\}\le 1$ and hence can be compared with a Uniformn draw.) Here is a Gibbs outcome for $a=10$ and $10^3$ Gibbs steps (using accept-reject sampling):

enter image description here

  • 1
    $\begingroup$ +1. After your comment, I realized the integrals were available in closed form, but lost my patience somewhere in the second step. Good on you for sticking with it. Also, I suppose a Laplace proposal centered at x(or y) may be more efficient. $\endgroup$ Feb 4, 2019 at 16:15
  • 1
    $\begingroup$ Hmm, but the support of Exp(1) is a subset of the support of this distribution. So That shouldn't work, right? In other words, you won't ever be able to propose/accept negative values. Am I missing something here? $\endgroup$ Feb 4, 2019 at 16:31
  • 1
    $\begingroup$ My bad: a double Exp(1)... aka Laplace distribution. $\endgroup$
    – Xi'an
    Feb 4, 2019 at 16:54

The conditional density kernels are:

$$\begin{equation} \begin{aligned} f(x|y) &\propto \exp(-|x|-a \cdot |x-y|), \\[6pt] f(y|x) &\propto \exp(-|y|-a \cdot |x-y|). \\[6pt] \end{aligned} \end{equation}$$

The difficulty here is to derive the actual densities that go with these kernels, which takes a bit of algebra. Integration over the full range of values gives the appropriate constant-of-integration:

$$C(a,y) \equiv \int \limits_\mathbb{R} \exp(-|x|-a \cdot |x-y|) \ dx.$$

To evaluate this constant you need to break your integral up into pieces so that you can remove the absolute value signs. I will not do all of this, but I will show you part of it, which will show how it can be done. (For simplicity I am going to assume that $0<a<1$.) In the case where $y>0$ we have:

$$\begin{equation} \begin{aligned} C(a,y) &= \int \limits_{-\infty}^{0} \exp(x-a (y-x)) \ dx \\ &\quad + \int \limits_{0}^{y} \exp(-x-a (y-x)) \ dx \\ &\quad + \int \limits_{y}^{\infty} \exp(-x-a (x-y)) \ dx \\[6pt] &= \frac{1}{1+a} \Bigg[ \exp((1+a)x-ay)) \Bigg]_{x \rightarrow -\infty}^{x=0} \\ &\quad - \frac{1}{1-a} \Bigg[ \exp((a-1)x-ay)) \Bigg]_{x=0}^{x=y} \\ &\quad - \frac{1}{1+a} \Bigg[ \exp(-(1+a)x + ay)) \Bigg]_{x=y}^{x \rightarrow \infty} \\[6pt] &= \frac{1}{1+a} \cdot \exp(-ay) \\ &\quad - \frac{1}{1-a} \cdot \exp(-y) + \frac{1}{1-a} \cdot \exp(-ay) \\ &\quad + \frac{1}{1+a} \cdot \exp(-y) \\[6pt] &= \Big( \frac{1}{1+a} + \frac{1}{1-a} \Big) \exp(-ay) + \Big( \frac{1}{1+a} - \frac{1}{1-a} \Big) \exp(-y) \\[6pt] &= \frac{2}{1-a^2} \cdot \exp(-ay) - \frac{2a}{1-a^2} \cdot \exp(-y). \\[6pt] \end{aligned} \end{equation}$$

So for $y>0$ you have:

$$f(x|y) = \frac{(1-a^2) \exp(-|x|-a \cdot |x-y|)}{2 \exp(-ay) - 2a \exp(-y)}.$$

You will need to derive the corresponding case where $y \leqslant 0$ and then you will have the full density. Both the conditional densities in this problem are the same, so you will then be able to implement the Gibbs sampler by generating from this density.


An alternative to the painful simulation from the exact full conditional distributions is to ressort to slice sampling, that is, to express the density in $(x,y)$ as the marginal of a density in $(x,y,u_1,u_2,u_3)$ as follows: \begin{align*}f(x,y)&\propto\exp(-|x|-|y|-a \cdot |x-y|)\\&=\int_0^\infty\mathbb{I}_{u_1\le\exp\{-|x|\}}\text{d}u_1\int_0^\infty\mathbb{I}_{u_2\le\exp\{-|y|\}}\text{d}u_2\int_0^\infty\mathbb{I}_{u_3\le\exp\{-a|x-y|\}}\text{d}u_3\end{align*} The full conditionals are then all uniforms \begin{align*} x|y,u_1,u_2,u_3&\sim f(x|y,u_1,u_2,u_3)\propto\mathbb{I}_{u_1\le\exp\{-|x|\}}\mathbb{I}_{u_3\le\exp\{-a|x-y|\}}\\ y|x,u_1,u_2,u_3&\sim f(y|x,u_1,u_2,u_3)\propto\mathbb{I}_{u_2\le\exp\{-|y|\}}\mathbb{I}_{u_3\le\exp\{-a|x-y|\}}\\ u_1|x&\sim\mathcal U(0,\exp\{-|x|\})\\ u_3|x,y&\sim\mathcal U(0,\exp\{-a|x-y|\})\\ u_2|y&\sim\mathcal U(0,\exp\{-|x|\})\\ \end{align*} hence straightforward to simulate:

for (t in 2:1e3){
 while (a*abs(xz[t]-yz[t-1])>-log(u3)) xz[t]=runif(1,log(u1),-log(u1))
 while (a*abs(xz[t]-yz[t])>-log(u3)) yz[t]=runif(1,log(u2),-log(u2))

$\qquad\qquad$enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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