|
clProbDist
An OpenCL library for probability distributions
|
|
clProbDist
An OpenCL library for probability distributions
|
We introduce clProbDist, an OpenCL library for probability distributions. It provides facilities for evaluating probability density or mass functions, distribution functions and their inverses, and reliability functions, currently for five distributions: normal, lognormal, exponential, gamma and Poisson. Applications include nonuniform variate generation, computing confidence intervals, computing \(p\)-values for statistical tests, computing densities to evaluate the likelihood function of a sample, etc.
This library, in its current state, defines a framework for working with probability distributions in OpenCL, and to show some working examples. Dozens of additional distributions need to be contributed to clProbDist to make it a mature library.
clProbDist allows users to generate nonuniform variates (see Example 2: non-uniform variate generation for illustration) through inversion only. Other (possibly faster) methods, including acceptance-rejection based approaches, need to be implemented somewhere else. The inversion method relies on evaluating the inverse cumulative distribution function (CDF) of a given distribution. For certain distributions, it may be profitable to compute lookup arrays in advance and use them to evaluate the inverse CDF repeatedly with the same parameters, as when generating many random variates from the same distribution. clProbDist provides facilities to do so (see Example 3: non-uniform variate generation revisited for an example).
The API is almost the same for every distribution. It is described in detail in clProbDist_template.h. Each distribution has its own header file (which is normally the lowercase name of the distribution, with a .h or .clh extension on the host or device, respectively), and each function contains the name of the distribution in its prefix (listed in Implemented distributions). In the generic probability distribution API reference, the distribution-specific part of the prefix is not shown.
For example, functions from the Normal distribution API have names that begin with clprobdistNormal and are declared in the normal.h header file. In particular, the probability density function clprobdistDensity(), the cumulative distribution function clprobdistCDF() and its inverse clprobdistInverseCDF() from the generic API, expand to clprobdistNormalDensity(), clprobdistNormalCDF() and clprobdistNormalInverseCDF(), respectively, for the normal distribution. More concretely, to store in p the probability density of x in a normal distribution of mean \(\mu\) = mu and standard deviation \(\sigma\) = sigma, one would invoke
Note that x, mu and sigma in the above are all floating-point values of type cl_double. The last argument can be used for reporting errors, but errors are ignored here by passing NULL.
The following code defines a function (on the host) that computes a normal confidence interval of level level on the mean of a set of n observations of type cl_double stored in the array observations. This assumes, of course, that the mean is normally distributed. The lower and upper bounds are returned in the lower and upper variables, respectively:
normalConfidenceInterval() simply computes the sample average and variance of the observations. The second part illustrates how to evaluate the inverse normal distribution function by invoking clprobdistNormalInverseCDF() with the normal distribution parameters \(\mu\)=mu and \(\sigma\)=sigma (the mean and variance) as its first arguments, for two different quantiles. The last argument is used for reporting errors, which are ignored here by setting it to NULL.
A more efficient implementation would compute the half width of the interval using a single call to clprobdistStdNormalInverseCDF() from the Standard normal distribution API. The standard normal distribution always has zero mean and a unit variance, so functions from its API do not take any distribution parameters. With this, we could replace the last four instructions of the normalConfidenceInterval defined above with:
The complete code for this example can be found in DocsTutorial/example1.c.
In this example, we want to simulate a stochastic system (not completely specified as our focus here is the usage of clProbDist) that uses a number, say 100, of Poisson variates, all with the same mean lambda, as its input. To generate a Poisson variate with the inversion method, we first generate a pseudorandom number uniformly distributed in (0,1) using the MRG31k3p generator from the clRNG library, and we map the output through clprobdistPoissonInverseCDF(). (The user is referred to the clRNG documentation and clRNG tutorial document for instructions on random streams usage.) For simplicity, we assume that the Poisson variates are used sequentially, one at a time, so we do not need to store them all at once in advance. Again, we run everything on the host. The code to average 1024 realizations of the output of this stochastic system would look like this:
Here, the output of simulateOneRun() is simply the sum of the 100 Poisson variates. In practice, it would be more elaborate, but the purpose of this example is to show how to generate nonuniform variates. The call to clprobdistPoissonInverseCDF(lambda, u, NULL) evaluates the inverse Poisson distribution function with mean lambda at quantile u, ignoring errors.
The complete code for this example can be found in DocsTutorial/example2.c.
The above code generates 1024 × 100 Poisson variates with the same mean, and can be very inefficient when lambda is large. A faster approach makes use of a lookup array to avoid repeating costly computations. Such facility is provided by a selection of functions of the Poisson distribution API, illustrated in Example 3: non-uniform variate generation revisited below. For this purpose, we introduce Distribution objects.
The API of each distribution comes in two variants. One uses distribution parameters as the first arguments of the functions; the other replaces the distribution parameters with a distribution object, and appends WithObject to the function names. [note]
The type of a distribution object is specific to the distribution it represents, e.g., the Poisson distribution has objects of type clprobdistPoisson. A distribution object stores the parameters of the distribution it represents, and, for some distributions, it may also store additional data.
In particular, for discrete distributions such as the Poisson one, the distribution object also stores lookup arrays to speed up the evaluation of the probability, of the distribution function and its inverse, and of the reliability function. Building the lookup arrays cost some time, so it is advisable to use a distribution object especially when these functions need to be evaluated many times. Before before using a distribution object with these functions, it must be created using clprobdistCreate(); it can be destroyed with clprobdistDestroy().
For other distributions, it should make very little difference in terms of performance which variant of the distribution API is used (with distribution objects or distribution parameters).
In the generic probability distribution API reference, the generic type name for a distribution object is clprobdistObject.
Here, we modify Example 2: non-uniform variate generation to use a distribution object in order to accelerate the multiple evaluations of the inverse distribution function. First, simulateOneRun() needs to receive a distribution object instead of the mean of the Poisson distribution, so its signature becomes
Next, we replace the call to clprobdistPoissonInverseCDF() to its variant that takes a distribution object:
Finally, we create the distribution object with:
and pass it to simulateOneRun() instead of the value of the mean:
The first argument to clprobdistPoissonCreate() is the value of distribution parameter \(\lambda\). For distributions with multiple parameters, their values would be listed as the first arguments to clprobdistCreate() (expanded to the proper name for the distribution). For instance, to create a distribution object for the normal distribution with mean mu and standard deviation sigma, we would call:
The last two arguments to clprobdistPoissonCreate(), or to clprobdistCreate() in general, can be used to return the size of the created distribution object and the error status. Here, we ignore them by setting them to NULL.
The complete code for this example can be found in DocsTutorial/example3.c.
Using clProbDist from device code is very similar to using it from host code. The device API variant that uses distribution parameters is the same as the host API. For example, to store in p the probability density of x in a normal distribution of mean \(\mu\) = mu and standard deviation \(\sigma\) = sigma, one would include the header file normal.clh and invoke:
Here, x, mu and sigma are all floating-point values of type double.
The device API has no functions to create or destroy distribution objects. Besides that, it is almost the same as the host API. Here is an example kernel that receives a normal distribution object and uses it to invoke clprobdistNormalDensity():
Again, the last argument to clprobdistNormalDensity(), that can be used for error reporting but is ignored here by setting it to NULL. It is possible to store the distribution objects into other memory types than global memory, but this requires additional instructions in the host and kernel codes. This is discussed in Device memory types.
Here, we adapt Example 3: non-uniform variate generation revisited, which runs on the host, to run on the device. First, we translate the computation code into device code, by replacing the extensions of the header names with .clh and by replacing cl_double with double, and by writing a kernel to receive data from the host, call simulateOneRun() and send data back to the host:
The __global OpenCL C keyword in front of const clprobdistPoisson* indicates that the Poisson distribution object is stored in the device's global memory.
The host code is responsible for creating the distribution object and sending it into the device's memory. The latter requires the size of the object to be known, so we allow clprobdistPoissonCreate() to return it in its second argument (dist_buf_size):
After this instruction, dist_buf_size contains the size in bytes of the distribution object pointed to by dist. The size of the array is determined at runtime; it may vary depending on the value of the distribution parameter \(\lambda\), set to 50.0 here. An OpenCL buffer must be created to copy the distribution object into the device's memory:
In the above code snippet, the context variable is an OpenCL context. Finally, the buffer must be set as the first argument of the kernel (stored in the kernel variable):
The complete code for this example can be found in DocsTutorial/example4.c and DocsTutorial/example4_kernel.cl.
For all features of the library to work properly, the CLPROBDIST_ROOT environment variable must be set to point to the installation path of the clProbDist package, that is, the directory under which lies the include/clProbDist subdirectory. Means of setting an environment variable depend on the operating system used.
The API of each distribution assumes that the distribution objects are stored in a specific type of memory. It defaults to global memory, but can be customized by the user by changing the value of the preprocessor symbol CLPROBDIST_<DIST>_OBJ_MEM, where <DIST> is the uppercase name of the distribution, before including the device header file of the distribution, to one of the following values:
For example, to store exponential distribution objects in constant memory, the device code should simply begin with:
For other memory types than global and constant memory, the device API provides a clprobdistCopyOverFromGlobal() function, where the clprobdist prefix has to be expanded with the distribution name as explained in Distributions, prefixes and headers. We give examples below of how to use this function.
To store exponential distribution objects in private memory instead, in addition to replacing CLPROBDIST_MEM_TYPE_CONSTANT with CLPROBDIST_MEM_TYPE_PRIVATE, each work item needs to make a private copy of the global distribution object and use the address of that copy with API calls, e.g.,
To store distribution objects in local memory, the user must first reserve local memory from the host code. For example, with a Poisson distribution object, the host code would include something like the following:
In the kernel, the distribution object must be copied from global memory to local memory, which is shared by all work items from the same work group. Perhaps the simplest way to do so is to have the first work item of each work group do the copy, e.g.,
This is substantially more complicated than using global memory. It may be faster or slower depending on the distribution and its parameters, the OpenCL hardware and the vendor implementation of OpenCL.
The following table lists the distributions that are currently implemented in clProbDist with the type of their support (continuous or discrete), the prefix for functions of the distribution API, and the name of the corresponding header file on the host and on the device.
| Distribution | Support Type | Prefix | Host Header File | Device Header File |
|---|---|---|---|---|
| Normal | continuous | clprobdistNormal | normal.h | normal.clh |
| Standard Normal | continuous | clprobdistStdNormal | normal.h | normal.clh |
| Lognormal | continuous | clprobdistLognormal | lognormal.h | lognormal.clh |
| Exponential | continuous | clprobdistExponential | exponential.h | exponential.clh |
| Gamma | continuous | clprobdistGamma | gamma.h | gamma.clh |
| Poisson | discrete | clprobdistPoisson | poisson.h | poisson.clh |
The API of continuous distributions provide a clprobdistDensity() function to compute the probability density; for discrete distributions, it is replaced by a clprobdistProb() function that computes the probability mass.
The common API for all distributions, in a generic form, is given in clProbDist_template.h. Some parts of the API vary across distributions, for instance, the type of return value of the inverse distribution function: cl_double for continuous distributions and cl_int for discrete distributions. This is in fact the type of the random variable described by the distribution, and its programmatic type is denoted by vartype in the generic API. Also, the number of distribution parameters vary, as well as their individual types and names. The token DIST_PARAMS substitutes for the list of distribution parameters in the generic API. The table below explicits the expansion of vartype and DIST_PARAMS for implemented distributions. The meaning of the distribution parameters is given below for each individual distribution.
In addition to the functions listed with generic API names in clProbDist_template.h, each parameterized distribution exposes functions to retrieve the values of the parameters a distribution object was created with.
With mean \(\mu\) and variance \(\sigma^2\), where \(\sigma > 0\), its density function is
\[ f(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left[ \frac{-(x-\mu)^2}{2\sigma^2} \right] \qquad (-\infty < x < \infty). \]
For the normal distribution, DIST_PARAMS expands to:
where mu \(=\mu\) and sigma \(=\sigma\). The API of the normal distribution adds the functions clprobdistNormalGetMu() and clprobdistNormalGetSigma() to retrieve the value of these parameters.
The CDF is evaluated using Uses the Chebyshev approximation proposed in [6], which gives 16 decimals of precision. The inverse CDF is computed using different rational Chebyshev approximations [2], which give 16 decimal digits of precision for \(2.2 \times 10^{-308} < u < 1\).
The complete API of the normal distribution is described in normal.h.
This special case of the normal distribution with \(\mu = 0\) and \(\sigma = 1\) has distribution function
\[ F(x) = \Phi(x) = \frac{1}{\sqrt{2 \pi}} \int_{-\infty}^x e^{-t^2/2} \,dt \qquad (-\infty < x < \infty). \]
For the standard normal distribution, there are no free parameters, so DIST_PARAMS expands to nothing.
The complete API of the standard normal distribution is also described in normal.h.
With scale parameter \(\mu\) and shape parameter \(\sigma^2\), where \(\sigma > 0\), its density function is
\[ f(x) = \frac{1}{\sqrt{2 \pi} \sigma x} \exp\left[ \frac{-(\ln x - \mu)^2}{2\sigma^2} \right] \qquad (0 < x < \infty). \]
Its distribution function is
\[ F(x) = \Phi((\ln x - \mu) / \sigma) \qquad (0 < x < \infty), \]
where \(\Phi\) is the distribution function for the Standard normal distribution. This inverse distribution function is given by:
\[ F^{-1}(u) = \exp\left[ \mu + \sigma\Phi^{-1}(u) \right] \qquad (0 \leq u < 1). \]
If \(\ln Y\) has a normal distribution, then \(Y\) has a lognormal distribution with the same parameters.
The CDF and inverse CDF are evaluated using the implementations for the normal distribution.
For the lognormal distribution, DIST_PARAMS expands to:
where mu \(=\mu\) and sigma \(=\sigma\). The API of the normal distribution adds the functions clprobdistLognormalGetMu() and clprobdistLognormalGetSigma() to retrieve the value of these parameters.
The complete API of the normal distribution is described in lognormal.h.
With mean \(1/\lambda\), where \(\lambda > 0\), its density, distribution and inverse distribution functions are, respectively,
\begin{align} f(x) &= \lambda e^{-\lambda x} & (x \geq 0) \\ F(x) &= 1 - e^{-\lambda x} & (x \geq 0) \\ F^{-1}(x) &= -\lambda^{-1} \ln(1 - u) & (0 \leq u < 1). \end{align}
The inverse CDF is computed exactly with the above formula.
For the exponential distribution, DIST_PARAMS expands to:
where lambda \(=\lambda\). The API of the exponential distribution adds the function clprobdistExponentialGetLambda() retrieve the value of this parameter.
The complete API of the exponential distribution is described in exponential.h.
With shape parameter \(\alpha > 0\) and scale parameter \(\lambda > 0\), its density function is
\[ f(x) = \frac{\lambda^\alpha x^{\alpha - 1} e^{\lambda x}}{\Gamma(\alpha)} \qquad (x > 0), \]
where
\[ \Gamma(\alpha) = \int_0^\infty x^{\alpha - 1} e^{-x} \,dx \]
is the usual gamma function. In particular, \(\Gamma(n) = (n - 1)!\) when \(n\) is a positive integer.
Several functions of the API for the gamma distribution are evaluated using an approximation with a precision of roughly \(d\) decimal digits, where \(d\) can be specified by the user. For consistency, \(d\) must always be specified, but is ignored by clprobdistGammaDensity(), clprobdistGammaMean(), clprobdistGammaStdDeviation() and clprobdistGammaVariance().
The CDF is approximated with an improved version of the algorithm in [1] . The functions clprobdistGammaCDF() and clprobdistGammaCDFWithObject() try to return \(d\) decimals digits of precision. For \(\alpha\) not too large (e.g., \(\alpha \leq 1000\)), \(d\) gives a good idea of the precision attained. The inverse CDF \(F^{-1}(u)\) is approximated by solving \(F(x) - u = 0\) for \(x\) with the bisection method when \(u \leq 10^{-8}\) or \(\alpha \leq 1.5\), or with the Brent-Dekker method (see [3] and [4]) otherwise.
For the gamma distribution, DIST_PARAMS expands to:
where alpha \(=\alpha\), lambda \(=\lambda\) and decprec \(=d\). The API of the gamma distribution adds the functions clprobdistGammaGetAlpha() and clprobdistGammaGetLambda() to retrieve the value of the first two parameters.
The complete API of the gamma distribution is described in gamma.h.
With mean \(\lambda \geq 0\), its mass and distribution functions are, respectively,
\begin{align} p(x) &= \frac{e^{-\lambda} \lambda^x}{x!} & (x = 0, 1, \ldots) \\ F(x) &= e^{-\lambda} \sum_{j=0}^x \frac{\lambda^j}{j!} & (x = 0, 1, \ldots) \end{align}
When \(\lambda > 200\), to compute \(F(x)\) and \(\bar F(x)\), clprobdistPoissonCDF() and clprobdistPoissonComplCDF() exploit the relationship \(F(x) = 1 - G_{x+1}(\lambda)\), where \(G_{x+1}\) is the gamma distribution function with parameters \((\alpha, \lambda) = (x + 1, 1)\). Otherwise, the probabilities are summed term by term. To evaluate the inverse CDF \(F^{-1}(u)\), clprobdistPoissonInverseCDF() performs a sequential search for the smallest value of \(x\) such that \(F(x) \geq u\), by incrementing \(x\) starting from \(x = 0\). When \(\lambda\) is very large, the probabilities \(p(x)\) are empirically negligible when \(x\) is far from \(\lambda\). So, when \(\lambda \geq 700\), the implementation of clprobdistPoissonInverseCDF() starts with a binary search for the lower boundary \(x_{\mathrm{min}} \in \{0,\dots,\lambda\}\) of an interval that contains the non-negligible probability values. Then, the sequential search begins at \(x = x_{\mathrm{min}}\).
If one has to compute \(p(x)\) or \(F(x)\) for several values of \(x\) with the same \(\lambda\), where \(\lambda\) is not too large, then it is more efficient to instantiate an object and use the API variant that uses distribution objects of type clprobdistPoisson, e.g., clprobdistPoissonProbWithObject() and clprobdistPoissonCDFWithObject(), for which return values are computed in advance and stored in lookup arrays. Probabilities smaller than \(10^{-16}\) are not stored in the Poisson distribution object, but are computed directly each time they are needed (which should be very seldom). The inverse CDF is evaluated by clprobdistPoissonInverseCDFWithObject() with a binary search.
For the Poisson distribution, DIST_PARAMS expands to:
where lambda \(=\lambda\). The API of the exponential distribution adds the function clprobdistPoissonGetLambda() retrieve the value of this parameter.
The complete API of the Poisson distribution is described in poisson.h. The device API expects Poisson distribution objects to be stored into global memory by default. See Device memory types for information on using other types of memory.
WARNING: The complementary distribution function, which can be evaluated with clprobdistPoissonComplCDF() and clprobdistPoissonComplCDFWithObject(), is defined as \(\bar F(j) = \mathbb P[X \geq j]\) (for integers \(j\)), so that for the Poisson distribution clProbDist, \(F(j) + \bar F(j) \neq 1\) since both include the term \(\mathbb P[X = j]\).
1.8.10