What is the equivalent for cdfs of MCMC for pdfs?

Question

In conjunction with a Cross Validated question on simulating from a specific copula, that is, a multivariate cdf $C(u_1,\ldots,u_k)$ defined on $[0,1]^k$, I started wondering about the larger picture, namely how, when given such a function, can one figure a generic algorithm to simulate from the corresponding probability distribution?

Obviously, one solution is to differentiate $C$ $k$ times to produce the corresponding pdf $\kappa(u_1,\ldots,u_k)$ and then call a generic MCMC algorithm like Metropolis-Hastings to produce a sample from $C$ (or $\kappa$).

Aside: Another solution is to stick to Archimedian copulas, using the Laplace-Stieljes transform for simulation, but this is not always possible in practice. As I found when trying to solve the aforesaid question.

My question is about avoiding this differentiating step in a generic way, if at all possible.

Answer 1

This is an attempt which I didn't completely work through, but too long for the comments section. It might be useful to put it here as another basic alternative for very low $k$. It does not require explicit differentiation + MCMC (but it does perform numerical differentiation, without MCMC).

Algorithm

For small $\varepsilon > 0$:

1. Draw $u_1 \sim C_1 \equiv C(U_1 = u_1,U_2 = 1,\ldots, U_k = 1)$. This can be easily done by drawing $\eta \sim \text{Uniform}[0,1]$ and computing $C_1^{-1}(\eta)$ (which, if anything, can be easily done numerically). This is a draw from the marginal pdf $u_1 \sim \kappa(u_1)$.
2. For $j = 2\ldots k$
   • Define $$D_j^{(\varepsilon)}(u_j|u_1,\ldots,u_{j-1}) \equiv \Pr\left( u_1 - \frac{\varepsilon}{2} \le U_1 \le u_1 + \frac{\varepsilon}{2} \land \dots \land u_{j-1} - \frac{\varepsilon}{2} \le U_{j-1} \le u_{j-1} + \frac{\varepsilon}{2} \land U_{j} \le u_j \land U_{j+1} \le 1 \dots \land U_{k} \le 1\right),$$ which can be computed as a difference of $C$ evaluated at various points (which in the naive way needs $O(2^{j-1})$ evaluations of $C$ for every evaluation of $D_j^{(\varepsilon)}$). $D_j^{(\varepsilon)}$ is the $\varepsilon$-approximate marginal conditional of $u_j$ given $u_1, \ldots, u_{j-1}$.
   • Draw $u_j \sim D_j^{(\varepsilon)}(u_j|u_1,\ldots,u_{j-1})$ as per point 1, which again should be easy to do with numerical inversion.

Discussion

This algorithm should generate i.i.d. samples from an $\varepsilon$-approximation of $C(u_1,\ldots,u_k)$, where $\varepsilon$ merely depends on numerical precision. There are practical technicalities to refine the approximation and make it numerically stable.

The obvious problem is that computational complexity scales as $O(2^{k})$, so, to put it generously, this is not very general in terms of $k$ (but the example you linked had $k = 3$, so perhaps this method is not completely useless -- I am not familiar with the typical scenario in which you would have access to the cdf). On the other hand, for very low-dimensional distributions it could work, and the cost is compensated by the fact that, unlike the other generic solution of "differentiating + MCMC", there is no need to compute derivatives, samples are i.i.d. and there is no tuning (aside the choice of $\varepsilon$, which should just be something slightly above machine precision). And perhaps there are ways to make this better than the naive approach.

As I mentioned, this is off the top of my head so there might be other issues.