Menu
Search

Library Reference

overview

Report a bug
If you spot a problem with this page, click here to create a GitHub issue.
Improve this page
Quickly fork, edit online, and submit a pull request for this page. Requires a signed-in GitHub account. This works well for small changes. If you'd like to make larger changes you may want to consider using a local clone.

mir.random.flex

Flex module that allows to sample from arbitrary random distributions.
License:
Boost License 1.0.
Authors:
Sebastian Wilzbach, Ilya Yaroshenko
This module is available in the extended configuration.
The Transformed Density Rejection with Inflection Points (Flex) algorithm can sample from arbitrary distributions given (1) its log-density function f, (2) its first two derivatives and (3) a partitioning into intervals with at most one inflection point.
These can be easily found by plotting f''. Inflection points can be identified by observing at which points f'' is 0 and an inflection interval which is defined by two inflection points can either be:
  • g(x) is entirely concave (f'' is entirely negative)
  • g(x) is entirely convex (f'' is entirely positive)
  • g(x) contains one inflection point (f'' intersects the x-axis once)
It is not important to identify the exact inflection points, but the user input requires:
  • Continuous density function g(x).
  • Continuous differentiability of g(x) except in a finite number of points which need to have a one-sided derivative.
  • Doubled continuous differentiability of g(x) except in a finite number of points which need to be inflection points.
  • At most one inflection point per interval
Internally the Flex algorithm transforms the distribution with a special transformation function and constructs for every interval a linear hat function that majorizes the pdf and a linear squeeze function that is majorized by the pdf from the user-defined, mutually-exclusive partitioning.

Efficiency rho

In further steps the algorithm splits those intervals until a chosen efficiency rho between the ratio of the sum of all hat areas to the sum of all squeeze areas is reached. For example an efficiency of 1.1 means that 10 / 110 of all drawn uniform numbers don't match the target distribution and need be resampled. A higher efficiency constructs more intervals, and thus requires more iterations and a longer setup phase, but increases the speed of sampling.

Unbounded intervals

In each unbounded interval the transformation and thus the density must be concave and strictly monotone.

Transformation function (T_c)

The Flex algorithm uses a family of T_c transformations:
  • `T_0(x) = log(x)
  • `T_c(x) = sign(c) * x^c
Thus c has the following properties
  • Decreasing c may decrease the number of inflection points
  • For unbounded domains, c > -1 is required
  • For unbounded densities, c must be sufficiently small, but should be great than -1. A common choice is -0.5
  • c=0 is the pure log transformation and thus decreases the vulnerability for under- and overflows

References Botts, Carsten, Wolfgang Hörmann, and Josef Leydold. " Transformed density rejection with inflection points." Statistics and Computing 23.2 (2013): 251-260.

auto flex(S, F0, F1, F2)(in F0 f0, in F1 f1, in F2 f2, S c, S[] points, S rho = 1.1, int maxApproxPoints = 1000)
if (isFloatingPoint!S);

auto flex(S, Pdf, F0, F1, F2)(in Pdf pdf, in F0 f0, in F1 f1, in F2 f2, S c, S[] points, S rho = 1.1, int maxApproxPoints = 1000)
if (isFloatingPoint!S);

auto flex(S, F0, F1, F2)(in F0 f0, in F1 f1, in F2 f2, S[] cs, S[] points, S rho = 1.1, int maxApproxPoints = 1000)
if (isFloatingPoint!S);

auto flex(S, Pdf, F0, F1, F2)(in Pdf pdf, in F0 f0, in F1 f1, in F2 f2, S[] cs, S[] points, S rho = 1.1, int maxApproxPoints = 1000)
if (isFloatingPoint!S);

auto flex(S, Pdf)(in Pdf pdf, in FlexInterval!S[] intervals)
if (isFloatingPoint!S);
The Transformed Density Rejection with Inflection Points (Flex) algorithm can sample from arbitrary distributions given its density function f, its first two derivatives and a partitioning into intervals with at most one inflection point. The partitioning needs to be mutually exclusive and sorted.
Parameters:
Pdf pdf probability density function of the distribution
F0 f0 logarithmic pdf
F1 f1 first derivative of logarithmic pdf
F2 f2 second derivative of logarithmic pdf
S c T_c family value
S[] cs T_c family array
S[] points non-overlapping partitioning with at most one inflection point per interval
S rho efficiency of the Flex algorithm
int maxApproxPoints maximal number of points to use for the hat/squeeze approximation
Returns:
Flex Generator.
struct Flex(S, Pdf) if (isFloatingPoint!S);
Data body of the Flex algorithm. Can be used to sample from the distribution.
Examples: 
import mir.math : approxEqual;
import std.meta : AliasSeq;
import mir.random.engine.xorshift : Xorshift;
auto gen = Xorshift(42);
alias S = double;
auto f0 = (S x) => -x^^4 + 5 * x^^2 - 4;
auto f1 = (S x) => 10 * x - 4 * x^^3;
auto f2 = (S x) => 10 - 12 * x^^2;
S[] points = [-3, -1.5, 0, 1.5, 3];

auto tf = flex(f0, f1, f2, 1.5, points, 1.1);
auto value = tf(gen);
const @property const(FlexInterval!S[]) intervals();
Generated partition points
S opCall(RNG)(ref RNG rng)
if (isRandomEngine!RNG);
Sample a value from the distribution.
Parameters:
RNG rng random number generator to use
Returns:
Array of length n with the samples
struct FlexInterval(S) if (isFloatingPoint!S);
Reduced version of Interval. Contains only the necessary information needed in the generation phase.
S lx;
Left position of the interval
S rx;
Right position of the interval
S c;
T_c family of the interval
LinearFun!S hat;
The majorizing linear hat function
LinearFun!S squeeze;
The linear squeeze function which is majorized by log(f(x))
S hatArea;
The total area that is spanned by the hat function in this interval
S squeezeArea;
The total area that is spanned by the squeeze function in this interval
FlexInterval!S[] flexIntervals(S, F0, F1, F2)(in F0 f0, in F1 f1, in F2 f2, in S[] cs, in S[] points, in S rho = 1.1, in int maxApproxPoints = 1000);
Calculate the intervals for the Flex algorithm for a T_c family given its density function, the first two derivatives and a valid start partitioning. The Flex algorithm will try to split the intervals until a chosen efficiency rho is reached.
Parameters:
F0 f0 probability density function of the distribution
F1 f1 first derivative of f0
F1 f1 second derivative of f0
S[] cs T_c family (single value or array)
S[] points non-overlapping partitioning with at most one inflection point per interval
S rho efficiency of the Flex algorithm
int maxApproxPoints maximal number of splitting points before Flex is aborted
Returns:
Array of FlexInterval's.
S flexInverse(bool common = false, S)(in S x, in S c);
Compute inverse transformation of a T_c family given point x. Based on Table 1, column 3 of Botts et al. (2013).
Parameters:
common whether c be 0, -0.5, -1 or 1
S x value to transform
S c T_c family to use for the transformation
Returns:
Flex-inversed value of x