Generating samples of a multivariate cauchy distribution - MathOverflow

Generating samples of a multivariate cauchy distribution

2011-05-28T19:18:05Z

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 by Peter Sarkoci for Generating samples of a multivariate cauchy distribution

2011-05-28T21:31:53Z

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$

Answer by simon for Generating samples of a multivariate cauchy distribution

2011-05-30T15:40:09Z

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.