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.

  • 1
    $\begingroup$ The link "the Laplace-Stieljes transform" seems to be broken now. $\endgroup$
    – jochen
    Nov 13, 2015 at 18:42

1 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).


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.


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.


Your Answer

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

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