Fast multivariate Cauchy sampling?

I've been using the code below to efficiently sample from multivariate Gaussian with diagonal covariance matrix. Now I want to try the same for Cauchy samples, however replacing NormalDistribution with CauchyDistribution fails with CompiledFunction::cfse: Compiled expression {{1.54755,-2.73921},{1.28524,-1.49352}} should be a machine-size real number.

I suspect CauchyDistribution is not supported by compile. Is there a workaround?

Compilation speeds up multivariate Gaussian sampling by an order of magnitude. (related post)

gaussianSampler[mu_, diag_] :=
With[{d = Length[diag]}, Assert[d > 0];
Compile[{{n, _Integer}},
Module[{vals},
vals =
mu + Sqrt[diag]*# & /@
RandomVariate[NormalDistribution[], {Max[n, 1], d}]]]];

sampler = gaussianSampler[{0, 0}, {1, 2}];
sampler[5]

• Erm. Why don't you just use RandomVariate[ MultinormalDistribution[\[Mu], DiagonalMatrix[diag]], {n, d}]? It seems to be at least twice as fast as sampler. Oct 25, 2023 at 22:03
• ConstantArray[mu, n] + RandomVariate[NormalDistribution[], {n, d}] . DiagonalMatrix[diag] would have the same effect. Oct 25, 2023 at 22:11
• So, maybe ConstantArray[mu, n] + RandomVariate[CauchyDistribution[0., 1.], {n, d}].DiagonalMatrix[diag] does what you want? Oct 25, 2023 at 22:12
• Also Maybe helpful here to bypass problems with CauchyDistribution not being compileable:  If $u$ is a uniform r.v. on [0,1], you can get a r.v. $x$ distributed as the standard Cauchy distribution by the transform: ${\displaystyle x=\tan \left(\pi (u-{\frac {1}{2}})\right)}$
– ydd
Oct 25, 2023 at 22:12
• @HenrikSchumacher that solution is slow for large dimensionality. Just forming the diagonal matrix takes longer than sampling Oct 25, 2023 at 22:42

Two approaches with timings...

Set parameters.

n = 100000; (* Sample size *)
d = 1000; (* Dimension *)
mu0 = RandomVariate[NormalDistribution[0, 1], d];
sigma0 = RandomVariate[ChiSquareDistribution[4], d];


Based on @ydd 's comment (assuming all values of sigma0 are positive):

cauchySampler1[mu_, sigma_] := Compile[{{n, _Integer}}, Module[{vals},
vals = mu + sigma*Tan[π (# - 1/2)] & /@ RandomReal[{0, 1}, {Max[n, 1], Length[sigma]}]]]
sampler1 = cauchySampler1[mu0, sigma0];
AbsoluteTiming[data = sampler1[n];]
(* {1.16369, Null} *)


Direct approach:

cauchySampler2[mu_, sigma_, n_] :=
mu + sigma # & /@ RandomVariate[CauchyDistribution[0, 1], {n, Length[sigma]}];
AbsoluteTiming[data = cauchySampler2[mu0, sigma0, n];]
(* {1.99362, Null} *)


Improvements in speed might likely need knowledge as to what you want to do with these samples especially making use of the lack of any dependence among samples.