Data.Random.Distribution.Ziggurat

A generic "ziggurat algorithm" implementation. Fairly rough right now.

There is a lot of room for improvement in findBin0 especially. It needs a fair amount of cleanup and elimination of redundant calculation, as well as either a justification for using the simple findMinFrom or a proper root-finding algorithm.

It would also be nice to add (preferably by pulling in an external package) support for numerical integration and differentiation, so that tables can be derived from only a PDF (if the end user is willing to take the performance and accuracy hit for the convenience).

Documentation

data Ziggurat v t

A data structure containing all the data that is needed to implement Marsaglia & Tang's "ziggurat" algorithm for sampling certain kinds of random distributions.

The documentation here is probably not sufficient to tell a user exactly how to build one of these from scratch, but it is not really intended to be. There are several helper functions that will build Ziggurats. The pathologically curious may wish to read the runZiggurat source. That is the ultimate specification of the semantics of all these fields.

Constructors

Ziggurat

Fields

zTable_xs :: !(v t)

The X locations of each bin in the distribution. Bin 0 is the infinite one.

In the case of bin 0, the value given is sort of magical - x[0] is defined to be V/f(R). It's not actually the location of any bin, but a value computed to make the algorithm more concise and slightly faster by not needing to specially-handle bin 0 quite as often. If you really need to know why it works, see the runZiggurat source or "the literature" - it's a fairly standard setup.

zTable_y_ratios :: !(v t)

The ratio of each bin's Y value to the next bin's Y value

zTable_ys :: !(v t)

The Y value (zFunc x) of each bin

zGetIU :: !(RVar (Int, t))

An RVar providing a random tuple consisting of:

• a bin index, uniform over [0,c) :: Int (where c is the number of bins in the tables)
• a uniformly distributed fractional value, from -1 to 1 if not mirrored, from 0 to 1 otherwise.

This is provided as a single RVar because it can be implemented more efficiently than naively sampling 2 separate values - a single random word (64 bits) can be efficiently converted to a double (using 52 bits) and a bin number (using up to 12 bits), for example.

zTailDist :: RVar t

The distribution for the final "virtual" bin (the ziggurat algorithm does not handle distributions that wander off to infinity, so another distribution is needed to handle the last "bin" that stretches to infinity)

zUniform :: !(t -> t -> RVar t)

A copy of the uniform RVar generator for the base type, so that Distribution Uniform t is not needed when sampling from a Ziggurat (makes it a bit more self-contained).

zFunc :: !(t -> t)

The (one-sided antitone) PDF, not necessarily normalized

zMirror :: !Bool

A flag indicating whether the distribution should be mirrored about the origin (the ziggurat algorithm it its native form only samples from one-sided distributions. By mirroring, we can extend it to symmetric distributions such as the normal distribution)

Instances

(Num t, Ord t, Vector v t) => Distribution (Ziggurat v) t

mkZigguratRec :: (RealFloat t, Vector v t, Distribution Uniform t) => Bool -> (t -> t) -> (t -> t) -> (t -> t) -> t -> Int -> RVar (Int, t) -> Ziggurat v t

Build a lazy recursive ziggurat. Uses a lazily-constructed ziggurat as its tail distribution (with another as its tail, ad nauseum).

Arguments:

• flag indicating whether to mirror the distribution
• the (one-sided antitone) PDF, not necessarily normalized
• the inverse of the PDF
• the integral of the PDF (definite, from 0)
• the estimated volume under the PDF (from 0 to +infinity)
• the chunk size (number of bins in each layer). 64 seems to perform well in practice.
• an RVar providing the zGetIU random tuple

mkZiggurat :: (RealFloat t, Vector v t, Distribution Uniform t) => Bool -> (t -> t) -> (t -> t) -> (t -> t) -> t -> Int -> RVar (Int, t) -> (t -> RVar t) -> Ziggurat v t

Build the tables to implement the "ziggurat algorithm" devised by Marsaglia & Tang, attempting to automatically compute the R and V values.

Arguments are the same as for mkZigguratRec, with an additional argument for the tail distribution as a function of the selected R value.

mkZiggurat_ :: (RealFloat t, Vector v t, Distribution Uniform t) => Bool -> (t -> t) -> (t -> t) -> Int -> t -> t -> RVar (Int, t) -> RVar t -> Ziggurat v t

Build the tables to implement the "ziggurat algorithm" devised by Marsaglia & Tang, attempting to automatically compute the R and V values.

Arguments:

• flag indicating whether to mirror the distribution
• the (one-sided antitone) PDF, not necessarily normalized
• the inverse of the PDF
• the number of bins
• R, the x value of the first bin
• V, the volume of each bin
• an RVar providing the zGetIU random tuple
• an RVar sampling from the tail (the region where x > R)

findBin0 :: RealFloat b => Int -> (b -> b) -> (b -> b) -> (b -> b) -> b -> (b, b)

I suspect this isn't completely right, but it works well so far. Search the distribution for an appropriate R and V.

Arguments:

• Number of bins
• target function (one-sided antitone PDF, not necessarily normalized)
• function inverse
• function definite integral (from 0 to _)
• estimate of total volume under function (integral from 0 to infinity)

Result: (R,V)

runZiggurat :: (Num a, Ord a, Vector v a) => Ziggurat v a -> RVar a

Sample from the distribution encoded in a Ziggurat data structure.