Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Recently, Oleksandr kindly showed a list of Mathematica commands that can be compiled. RandomVariate was part of that list. However, whether this can be compiled depends upon the distribution that is being sampled.

Needs["CompiledFunctionTools`"]

cf1 = Compile[{{m, _Real}, {s, _Real}},
  Module[{v1, v2, v3, v4, v5, v6},
   v1 = RandomVariate[NormalDistribution[m, s]];
   v2 = RandomVariate[UniformDistribution[{m, s}]];
   v3 = RandomVariate[GammaDistribution[m, s]];
   v4 = RandomVariate[PoissonDistribution[m]];
   v5 = RandomVariate[ChiSquareDistribution[m]];
   v6 = RandomVariate[ExponentialDistribution[m]];
   {v1, v2, v3, v4, v5, v6}
   ]
  ]

Using CompilePrint shows that RandomVariate can be compiled for the Normal Distribution or the Uniform Distribution and not with some others.

CompilePrint[cf1]

2 arguments
        4 Integer registers
        8 Real registers
        1 Tensor register
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        R0 = A1
        R1 = A2
        Result = T(R1)0

1   R2 = RandomNormal[ R0, R1]]
2   R3 = RandomReal[ R0, R1]]
3   I0 = MainEvaluate[ Function[{m, s}, 
             RandomVariate[GammaDistribution[m, s]]][ R0, R1]]
4   I1 = MainEvaluate[ Function[{m, s}, 
             RandomVariate[PoissonDistribution[m]]][ R0, R1]]
5   I2 = MainEvaluate[ Function[{m, s}, 
             RandomVariate[ChiSquareDistribution[m]]][ R0, R1]]
6   I3 = MainEvaluate[ Function[{m, s}, 
             RandomVariate[ExponentialDistribution[m]]][ R0, R1]]
7   R4 = I0
8   R5 = I1
9   R6 = I2
10  R7 = I3
11  T(R1)0 ={ R2, R3, R4, R5, R6, R7 }
12  Return

Does anyone have a list of all the distributions that can be compiled (including, PDF, CDF and RandomVariate functionality)?

share|improve this question
    
Well, you can check each of them, or semi-automate checking. When one is not compiled, and you try to use it, you'll easily find out anyway: see my answer to that post asking about compilable functions. –  Szabolcs Feb 1 '12 at 19:33
    
@szabolcs There are too many combinations to try out, unfortunately. –  asim Feb 2 '12 at 1:43
3  
Beta distribution can be compiled too.. –  asim Feb 10 '12 at 14:09
    
Since the results of RandomVariate are ALWAYS going to be Reals (for a continuous distribution), and Integers (for most inbuilt discrete distributions), I can't see why there would be any advantage to compiling them, that could not have automatically been built into the function. Does anyone have some timing tests where manual compilation yields substantive benefits to something like: RandomVariate[ dist, {10^6}] ? –  wolfies Aug 5 at 15:37

3 Answers 3

up vote 15 down vote accepted

To my knowledge UniformDistribution and NormalDistribution are the only distributions that are directly compilable for RandomVariate.

Consider that sampling from a UniformDistribution is what RandomReal was originally designed to do. This code is likely written deep down in C and so compiles without any special effort. In order to hook up RandomVariate for uniforms Compile just needs to recognize that this is really just a call to RandomReal.

Now, sampling from a NormalDistribution is so common that it was considered worth the time investment to make it compilable. Notice that the call to RandomVariate actually produces a call to RandomNormal which was almost certainly written for this purpose.

As for other distributions, special code would need to be written for each one in a similar fashion to RandomNormal for them to be "supported" by Compile. Since there are well over 100 of these, it would be a huge undertaking. An argument could be made for doing this for a few distributions but who is to decide which ones are most important?

There is a sunny side. Most distributions have their own dedicated and highly optimized methods for random number generation. Often Compile is used under the hood when machine precision numbers are requested.

Because of this, even if they were directly compilable you probably wouldn't see much of a speed boost since the code is already optimized.

Fortunately Compile can happily handle arrays of numbers. I typically just rely on the optimized code used by RandomVariate to generate the numbers and subsequently pass them in as an argument to the compiled function.

Incidentally, everything I just said about RandomVariate is also true of distribution functions like PDF, CDF, etc. Obviously these are just pure functions (in the univariate case) and unless they are built with some exotic components they should compile assuming you evaluate them before putting them into your compiled function.

share|improve this answer
    
Also of note is that Erf[] is in the list of compilable functions in Olek's answer. So yes, I think the normal distribution is pretty well-covered here. –  J. M. Feb 2 '12 at 6:29
    
I just found out that RandomVariate[BetaDistribution[a,b]] can also be compiled. By the way, it is not always that simple to generate the random numbers before passing them to a compiled function, especially, if the number of draws are endogenously determined within the function, as part of some while loop. –  asim Feb 10 '12 at 14:05
    
@asim this is definitely true. I was unaware the BetaDistribution was compilable, thanks! –  Andy Ross Feb 10 '12 at 16:09
    
@asim I have had situations like this. Typically I generate far more numbers than I think I will need and pass them in. It is faster to generate a bunch in one go than to repeatedly call RandomVariate. The downside is having to leave Compile and generate more if I didn't generate enough numbers in the first go. –  Andy Ross Feb 10 '12 at 23:22

Another option to the answer posted by @Andy Ross cropped up in a recent question of mine about corrupting an image with Poisson noise. In my own answer, I made use of LibraryLink to utilise the distributions built into C++, which you can find here: http://www.cplusplus.com/reference/random/

This was especially useful in my case because Poisson noise in an image otherwise relies on a call to RandomVariate for each pixel (although the built-in Mathematica function for doing this via ImageEffect works differently, but is still slower than my answer).

Repeated calls to RandomVariate rapidly become time-consuming:

(* Generate random list *)
listofnums1d = RandomReal[{0, 100}, {10^6}];

(* Calling RandomVariate each time takes 146.78 seconds *)
RandomVariate[PoissonDistribution[#]] & /@ listofnums1d; // AbsoluteTiming

(* Loading my C++ function via LibraryLink (code below) takes 0.77 seconds *) 
mypoissondistribution[listofnums1d]; // AbsoluteTiming

C++ has a whole host of distributions built into <random>, including Gamma, Exponential, ChiSquare, Cauchy and Weibull.

If repeated calls to RandomVariate are a problem, this could be a way to go.


C++ and LibraryLink code

#include "WolframLibrary.h"
#include "stdio.h"
#include "stdlib.h"
#include <random>
#include <chrono>
EXTERN_C DLLEXPORT mint WolframLibrary_getVersion(){return WolframLibraryVersion;}
EXTERN_C DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {return 0;}
EXTERN_C DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {}

EXTERN_C DLLEXPORT int poissondistribution(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){

    int err; // error code

    MTensor m1; // input tensor
    MTensor m2; // output tensor

    mint const* dims; // dimensions of the tensor

    double* data1; // actual data of the input tensor
    double* data2; // data for the output tensor

    mint i; // bean counter

    m1 = MArgument_getMTensor(Args[0]);
    dims = libData->MTensor_getDimensions(m1);
    err = libData->MTensor_new(MType_Real, 1, dims,&m2);
    data1 = libData->MTensor_getRealData(m1);
    data2 = libData->MTensor_getRealData(m2);

    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
    std::default_random_engine generator (seed); // RNG

    for(i=0;i<dims[0];i++) {
        if(data1[i]==0.)
        {
            data2[i] = 0.; 
            // If value is 0, return 0. 
            // You can avoid this by not passing a value of 0, of course.
        }
        else
        {
            std::poisson_distribution<int> poissondistribution(data1[i]); // Poisson distribution
            data2[i] = poissondistribution(generator); // If value is >0, sample from the Poisson distribution
        }
    }

    MArgument_setMTensor(Res, m2);
    return LIBRARY_NO_ERROR;
}

Put the C++ code into a file (e.g. in the notebook directory) and you can compile it like this.

Needs["CCompilerDriver`"]

distributionlib = 
  CreateLibrary[{FileNameJoin[{NotebookDirectory[], "filename"}]},
  "libraryname", "Debug" -> False, 
  "CompileOptions" -> "-O2 -std=c++11"];

Once the library has been compiled and you stop making changes to the C++ code, you can save time by taking it from the default location Mathematica places it ($LibraryPath, found using FindLibrary[]), and move it into your notebook directory and call it. Remember that the file extension is .dll on Windows, .so on Linux, and .dylib on Mac.

This works because LibraryFunctionLoad[] uses FindLibrary internally.

mypoissondistribution = 
  LibraryFunctionLoad[
   FileNameJoin[{NotebookDirectory[], "filename"}], 
   "functionname", {{Real, 1}}, {Real, 1}];

Note that this code takes and returns a 1D list of values. The code in my linked question about images takes a 2D list.

I believe that you can then include mypoissondistribution in a compiled function without a call to MainEvaluate. Perhaps someone can clarify? I don't know if you need to set "InlineExternalDefinitions" and/or "InlineCompiledFunctions" to True for this to work.

share|improve this answer
1  
This is very neat!! +1 –  brama Aug 6 at 10:28

Here is a quick (and dirty?) implementation of the Laplace distribution. Relationship to the uniform distribution as given on Wikipedia:

randomLaplace = With[{$MachineEpsilon = $MachineEpsilon},
  Compile[{{µ, _Real, 0}, {b, _Real, 0}, {n, _Integer, 0}},
   With[{u = RandomReal[{-1/2, 1/2} (1 - $MachineEpsilon), n]},
    µ - b Sign[u] Log[1 - 2 Abs[u]]
   ]
  ]
 ];

Histogram of compiled Laplace distribution results

It looks fine compared to the built-in version:

Histogram of built-in Laplace distribution results

share|improve this answer
    
+1, nice idea and implementation! –  blochwave Aug 5 at 11:39

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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