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:

Authors:

Sebastian Wilzbach, Ilya Yaroshenko
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)

- 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

### 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

- 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 approximationReturns: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 std.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 abortedReturns: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 1S `x`

value to transform S `c`

T_c family to use for the transformation Returns:Flex-inversed value of`x`

Copyright © 1999-2017 by the D Language Foundation | Page generated by
Ddoc on Thu Mar 30 05:51:42 2017