Generating samples of a multivariate cauchy distribution

Question

The question is very simple: Do you know an efficient algorithm to generate samples of the multivariate cauchy distribution

• http://en.wikipedia.org/wiki/Cauchy_distribution#Multivariate_Cauchy_distribution

The pdf is given by

$f(x_i;\mu_i,\Sigma,k) = \frac{\Gamma\left(\frac{1+k}{2}\right)}{\pi^{(1+k)/2} \left[(x_i-\mu_i)\Sigma_{ij}(x_j-\mu_j) + 1\right]^{(1+k)/2}}$

where the indices $i,j=1,...,k$ and $\Sigma_{ij}$ is a $k\times k$ positive definite matrix.

Using that a multivariate cauchy distribution is really a multivariate Student's distribution with one degree of freedom, I would be happy also knowing how to generate samples of a Student's distribution.

Answer 1

As others have pointed out, a multivariate Cauchy is a multivariate Student's $t$ with 1 degree of freedom.

One straightforward approach would be to generate a $Z \sim \chi^2(1)$, and then $X|Z \sim N_k(\mu,\frac{1}{Z}\Sigma)$.

To generate $\chi^2(1)$, it is probably simplest to generate a single standard normal and square it. To generate a multivariate normal, the usual method is to generate $k$ i.i.d. standard normals, multiply by a Cholesky decomposition of the covariance matrix, and add the mean.

Answer 2

Not really an answer. One neat (Edit: not so neat in this case) way of sampling out of the Cauchy (and, actually, student as well) multivariate distribution exploits its (their) ellipticity. In order to sample an elliptical random vector it is enough to proceed along the following steps:

• generate a point $x$ uniformly on the unit sphere of $\mathbb{R}^{n}$; with $n$ being the dimmension of your r.v.
• independently on $x$, generate a real number $\alpha$ according to some univariate law (in your case, it is Cauchy with $\gamma=1$ Edit: as fedja pointed out, this is not so simple -- actually, one has to sample from absolute value of $n$-variate Cauchy, which seems to be the main obstacle here.)
• compute $\alpha x$ and transform the resulting point according to the matrix $\Sigma$ (via a matrix $A$ such that $AA^{T}=\Sigma$)
• advance the resulting point by your desired mean $\mu$