Report a bug
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.
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`