3
$\begingroup$

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]
Improve this question
$\endgroup$
8
  • $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    Commented Oct 25, 2023 at 22:03
  • $\begingroup$ ConstantArray[mu, n] + RandomVariate[NormalDistribution[], {n, d}] . DiagonalMatrix[diag] would have the same effect. $\endgroup$
    – Henrik Schumacher
    Commented Oct 25, 2023 at 22:11
  • $\begingroup$ So, maybe ConstantArray[mu, n] + RandomVariate[CauchyDistribution[0., 1.], {n, d}].DiagonalMatrix[diag] does what you want? $\endgroup$
    – Henrik Schumacher
    Commented Oct 25, 2023 at 22:12
  • 4
    $\begingroup$ 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)}$ $\endgroup$
    – ydd
    Commented Oct 25, 2023 at 22:12
  • $\begingroup$ @HenrikSchumacher that solution is slow for large dimensionality. Just forming the diagonal matrix takes longer than sampling $\endgroup$
    – Yaroslav Bulatov
    Commented Oct 25, 2023 at 22:42

1 Answer 1

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

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

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