New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

mir.random.nonuniform: Add Ziggurat method for Normal & Exponential #261

Closed
wants to merge 6 commits into
base: master
from

Conversation

Projects
None yet
5 participants
@wilzbach
Member

wilzbach commented Jul 19, 2016

Adds the Ziggurat sampling algorithm for Normal & Exponential distribution.

  1. Paper

Marsaglia, George, and Wai Wan Tsang. "The ziggurat method for generating random variables."
Journal of statistical software 5.8 (2000): 1-7.

  1. Thoughts
  1. Distribution plots (dub ./examples/nonuniform_plot.d):

image

image

Ping @joseph-wakeling-sociomantic @9il

@wilzbach wilzbach referenced this pull request Jul 19, 2016

Closed

mir.random.nonuniform #262

0 of 7 tasks complete
}
};
return Ziggurat!(T, fallback, R, true)(pdf, invPdf, 128, rightEnd, T(9.91256303526217e-3));

This comment has been minimized.

@wilzbach

wilzbach Jul 19, 2016

Member

I actually would like to run the initialization in CTFE as it will never change, but exp uses inline assembler which isn't supported in CTFE :/
Has anyone an idea?

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Jul 19, 2016

Well, start by filing an issue against phobos asking for a CTFE'able exp.

BTW where on earth does this magic constant 9.91256...e-3 come from? I would suggest making it a named manifest constant.

This comment has been minimized.

@wilzbach

wilzbach Jul 19, 2016

Member

Well, start by filing an issue against phobos asking for a CTFE'able exp.

Ok thanks - done. I will test whether copying the non-inline version from Phobos works.

BTW where on earth does this magic constant 9.91256...e-3 come from? I would suggest making it a named manifest constant.

It also comes from [Marsaglia00] - there it's called v.
It's the area of every block and thus depends on k and the distribution.

-> I will declare it more explicitly.

/// precalculate scaling to R.max for x_i
T[] xScaled;
/// precalculate pdf value for x_i
T[] fs;

This comment has been minimized.

@wilzbach

wilzbach Jul 19, 2016

Member

"An Improved Ziggurat Method to Generate Normal Random Samples" postulates that saving the fs isn't necessary to which I couldn't follow yet.

DRanU() returns a uniform random number, U (0, 1), and IRanU() returns 32-bit
unsigned random integer. DRanNormalTail is implemented as a separate function:
it gets called only rarely, so that efficiency does not matter. For the same reason, it
is not necessary to avoid a call to exp() when checking for the wedges (this could
be achieved by precomputing the function values f (xi). (page 7)

Imho saving doesn't waste much space and saves time

T function(T x) invPdf;
/// precalculate difference x_i / x_{i+1}
T[] xDiv;

This comment has been minimized.

@wilzbach

wilzbach Jul 19, 2016

Member

As the array size is known at compile-time, should we use T[k]?

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Jul 19, 2016

Don't see why not, but right now k is not provided at compile time AFAICS?

This comment has been minimized.

@wilzbach

wilzbach Jul 20, 2016

Member

Don't see why not, but right now k is not provided at compile time AFAICS?

Well it's not strictly needed at compile-time (in comparison to the other compile-time parameters), but I don't see any use case where one wants to be able to choose the block size at runtime as the rightEnd and averageArea need to be precomputed by hand anyways.

Marsaglia, George, and Wai Wan Tsang. "The ziggurat method for generating random variables."
Journal of statistical software 5.8 (2000): 1-7.
*/
struct Ziggurat(T, string _fallback, R = uint, bool bothSides)

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Jul 19, 2016

Instead of R I would use UIntType (it matches the typical template-parameter name in both phobos std.random and C++11 <random> for the word-type of the uniform RNG).

However, in this case, is there any possibility to avoid needing the word-type to be known at compile time?

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Jul 19, 2016

Also, note that while a dependency on the word size is potentially more flexible, it's worth considering also just templating the ziggurat on the actual RNG type.

This comment has been minimized.

@wilzbach

wilzbach Jul 20, 2016

Member

However, in this case, is there any possibility to avoid needing the word-type to be known at compile time?

Well you already proposed it.

Also, note that while a dependency on the word size is potentially more flexible, it's worth considering also just templating the ziggurat on the actual RNG type.

Hmm while I understand the motivation of providing a simple API for the user, this would create even more template bloat (for every RNG, RNG template in Ziggurat template) and limit the API more than needed. Hence I think making opCall generic depending on the RNG is probably the better way to go.

}
/// samples a value from the discrete distribution using a custom random generator
T opCall(RNG)(ref RNG gen) const

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Jul 19, 2016

Assuming you do need to know the word-type (your template parameter R) at compile time, you might want to validate that typeof(gen.front) matches it.

Note that in principle at least it's possible to handle different R and generator return type: if R's size is a multiple of the generator's word size, then you can populate it by several calls to the generator; conversely, if R is smaller than the generator's word size, you can use a single RNG variate to provide several of the needed values. But this may be adding excessive complexity.

This comment has been minimized.

@wilzbach

wilzbach Jul 20, 2016

Member

Assuming you do need to know the word-type (your template parameter R) at compile time, you might want to validate that typeof(gen.front) matches it.

Nice idea (done).

Note that in principle at least it's possible to handle different R and generator return type: if R's size is a multiple of the generator's word size, then you can populate it by several calls to the generator; conversely, if R is smaller than the generator's word size, you can use a single RNG variate to provide several of the needed values. But this may be adding excessive complexity.

Good point - how common is it to have something different than uint or ulong?

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Jul 20, 2016

Good point - how common is it to have something different than uint or ulong?

Not common, I think, at least not these days.

import mir.internal.math : exp, log;
auto pdf = (T x) => cast(T) exp(-x);
auto invPdf = (T x) => cast(T) -log(x);

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Jul 19, 2016

I'm not sure I like the explicit casts. Surely it's possible to get the same result here without them?

This comment has been minimized.

@wilzbach

wilzbach Jul 20, 2016

Member

(Idk why my comment wasn't saved). The problem is that exp only accepts and returns real.

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Aug 5, 2016

Just as the normal implementation looks to be missing mean + variance, here the exponential implementation is missing its \lambda control parameter.

}
};
return Ziggurat!(T, fallback, R, false)(pdf, invPdf, 256, T(7.697117470131487), T(3.949659822581572e-3));

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Jul 19, 2016

While calculations using pdf and invPdf can only be done at runtime for now (because of the exp implementation issues), is there any reason why the actual lambdas can't be provided as template parameters? It would make for a more logical design, I think (and also be future-proof against the point when you get a CTFE'able exp).

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Jul 19, 2016

I would also suggest defining explicit templates for Normal(T), Exponential(T) etc. to wrap these Ziggurat instantiations.

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Jul 19, 2016

An alternative would be to define actual Normal(T) and Exponential(T) wrapper structs that just use a Ziggurat instance internally; that might also give you a more future-proof API should you ever wish to rework the internals.

That might be better, because it'll probably be easier for users to debug stuff if they see a type called Normal!double instead of Ziggurat!(LOTS, OF, DIFFERENT, STUFF).

This comment has been minimized.

@wilzbach

wilzbach Jul 19, 2016

Member

I put everything as template argument that is absolutely needed at compile time.
Does it make a huge difference? - I thought in the end we can just use enum z = Ziggurat!..

On July 19, 2016 2:15:03 PM GMT+02:00, Joseph Wakeling notifications@github.com wrote:

  • import mir.internal.math : exp, log;
  • auto pdf = (T x) => cast(T) exp(-x);
  • auto invPdf = (T x) => cast(T) -log(x);
  • // values from [Marsaglia00]
  • enum fallback = q{
  •    T fallback(RNG)(ref RNG gen) const
    
  •    {
    
  •        import std.random : uniform;
    
  •        auto u = uniform!("[]", T, T)(0, 1, gen);
    
  •        return 7.69711 - u;
    
  •    }
    
  • };
  • return Ziggurat!(T, fallback, R, false)(pdf, invPdf, 256,
    T(7.697117470131487), T(3.949659822581572e-3));

While calculations using pdf and invPdf can only be done at runtime
for now (because of the exp implementation issues), is there any
reason why the actual lambdas can't be provided as template parameters?
It would make for a more logical design, I think (and also be
future-proof against the point when you get a CTFE'able exp).


You are receiving this because you authored the thread.
Reply to this email directly or view it on GitHub:
https://github.com/libmir/mir/pull/261/files/6a6565f321a80c67b95a7657bf55cc77989fa09c#r71325812

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Jul 19, 2016

It makes less of a difference if the Ziggurat instance is wrapped away inside a type that's explicitly Normal or Exponential or whatever. But if you're set on returning a raw Ziggurat instantiation from your normal and exponential factory functions, it might be useful to have the PDF and inverse PDF clearly there in the template parameters.

This comment has been minimized.

@wilzbach

wilzbach Jul 20, 2016

Member

An alternative would be to define actual Normal(T) and Exponential(T) wrapper structs that just use a Ziggurat instance internally; that might also give you a more future-proof API should you ever wish to rework the internals.

I really like the idea of having a broad Normal(T, UIntType) API, but with my current attempt (see below - for some reasons the PR didn't go through yesterday) it is yet another function call.

wilzbach@4c2f4bd
9

Should I build a Ziggurat struct with mixins?

@codecov-io

This comment has been minimized.

codecov-io commented Jul 20, 2016

Current coverage is 96.74% (diff: 97.77%)

Merging #261 into master will increase coverage by 0.21%

@@             master       #261   diff @@
==========================================
  Files            19         21     +2   
  Lines          3262       3874   +612   
  Methods           0          0          
  Messages          0          0          
  Branches          0          0          
==========================================
+ Hits           3149       3748   +599   
- Misses          113        126    +13   
  Partials          0          0          

Sunburst

Powered by Codecov. Last update d24907b...b7541d9

@wilzbach

This comment has been minimized.

Member

wilzbach commented Jul 20, 2016

(will post more summaries here soon, here's a brief overview)

Design flaws of the Ziggurat algorithm

two main critics:

  • same variable is used to pick the block & value (last 7 or 8 bits). However 2^50 values are needed to detect this
    • SHR3 not uniform (shouldn't affect us)

An Improved Ziggurat Method to Generate Normal Random Samples

  • use double
  • two random variables (one separate 7bit one to pick the block) -> optimization uses 1 + 1/4 random variables
  • no precomputation of f(x_i)
  • idea: use two random integers for higher-precisions values

wilzbach added some commits Jul 20, 2016

@wilzbach

This comment has been minimized.

Member

wilzbach commented Jul 20, 2016

btw an interesting overview paper is Gaussian Random Number Generators by Thomas et. al.
Summary: Wallace-method is the fastest, but doesn't provide good statistical quality. Ziggurat was the second fastest method among the huge benchmark while passing the chi-squared and achieving good scores at the high sigma test (see table 3 and 4).

Authors: Sebastian Wilzbach
*/
module mir.random.nonuniform;

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Aug 5, 2016

Minor, but didn't we agree to create a mir.random.distribution package to contain everything ... ?

This comment has been minimized.

@wilzbach
gen.popFront();
// TODO: this is a bit biased
size_t i = u & kMask;

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Aug 5, 2016

Given your TODO here, I assume this is the first place you're looking to understand why the plotted distributions don't look quite right .... ?

This comment has been minimized.

@wilzbach

wilzbach Aug 5, 2016

Member

Yes resolving whether we can reuse the random bits or need to generate another 1/4 random variable (one of the optimizations described in later papers), shouldn't matter so much :/

T x, y, u;
do
{
u = uniform!("[]", T, T)(0, 1, gen);

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Aug 5, 2016

Marsaglia & Tang talk about UNI (in their paper) being a generator of "uniform (0, 1) variates". That would suggest over the open rather than closed interval, i.e. uniform!"()" rather than uniform!"[]". Does that make a difference to your distribution results?

// TODO: this is a bit biased
size_t i = u & kMask;
//size_t i = uniform!("[)", size_t, size_t)(0, kMask, gen);

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Aug 5, 2016

Minor: note that uniform!T should give uniform distribution across all possible values of an integral type T. But in this case arguably unnecessary.

import mir.internal.math : exp, log, sqrt;
auto pdf = (T x) => cast(T) exp(T(-0.5) * x * x);
auto invPdf = (T x) => cast(T) sqrt(T(-2) * log(x));

This comment has been minimized.

@joseph-wakeling-sociomantic

joseph-wakeling-sociomantic Aug 5, 2016

Aren't these definitions missing the mean and variance parameters of the normal distribution? I recognize that one can generate any normal distribution once one has variates from N(0, 1), but surely that should be baked in to the implementation?

Also: even allowing for N(0, 1), isn't the PDF missing the divisor by \sqrt{2 * \pi} ... ? & presumably this means the inverse PDF is also missing something correspondingly?

@9il

This comment has been minimized.

Member

9il commented Aug 5, 2016

Does Ziggurat method yield better results for normal distribution then other methods in Atmosphere or dstats?

@wilzbach

This comment has been minimized.

Member

wilzbach commented Aug 5, 2016

Does Ziggurat method yield better results for normal distribution then other methods in Atmosphere or dstats?

What's the best way to compare?

@wilzbach

This comment has been minimized.

Member

wilzbach commented Aug 5, 2016

@joseph-wakeling-sociomantic

NormalDist CDF looks good from what I can judge

image

Exp doesn't

image

@9il

This comment has been minimized.

Member

9il commented Aug 5, 2016

Does Ziggurat method yield better results for normal distribution then other methods in Atmosphere or dstats?
What's the best way to compare?

If you don't know if Ziggurat yield better result, then we do not need Ziggurat

@9il

This comment has been minimized.

Member

9il commented Aug 5, 2016

What's the best way to compare?

There are no the best way. The best option is to read a couple of articles. We need to understand need we Ziggurat or not first, before spend time on it

@joseph-wakeling-sociomantic

This comment has been minimized.

joseph-wakeling-sociomantic commented Aug 5, 2016

If you don't know if Ziggurat yield better result, then we do not need Ziggurat

The existing academic literature would suggest that Ziggurat is a very effective method; the paper by Thomas et al. offers a variety of tests of statistical quality that could be used, IIRC.

Question is, given the timelines, is it worth pushing on with Ziggurat, or would it be better to implement more basic implementations of the various distributions, and return to Ziggurat as a longer-term work?

@9il

This comment has been minimized.

Member

9il commented Aug 5, 2016

Question is, given the timelines, is it worth pushing on with Ziggurat, or would it be better to implement more basic implementations of the various distributions, and return to Ziggurat as a longer-term work?

First question is Ziggurat better for basic distributions than basic implementations of them? If there no strong Yes, then we the next after Tinflex is basic implementations.

@wilzbach

This comment has been minimized.

Member

wilzbach commented Aug 5, 2016

The existing academic literature would suggest that Ziggurat is a very effective method; the paper by Thomas et al. offers a variety of tests of statistical quality that could be used, IIRC.

Yes sorry I should have been more precise. @9il what would be needed to convince you that Ziggurat is better than the algorithms in Atmosphere or dstats?
Is a X^2 test (that's what they use to evaluate the goodness of the fit) ok or should I also do the high-sigma test (more complex)?

Question is, given the timelines, is it worth pushing on with Ziggurat, or would it be better to implement more basic implementations of the various distributions, and return to Ziggurat as a longer-term work?

We have time at least time until October. I would prefer to go with the "better" algorithm and tune it. We already have the basics implementations in dstats against which we can benchmark and compare.

@wilzbach

This comment has been minimized.

Member

wilzbach commented Aug 5, 2016

First question is Ziggurat better for basic distributions than basic implementations of them? If there no strong Yes, then we the next after Tinflex is basic implementations.

Quote from my comment above:

btw an interesting overview paper is Gaussian Random Number Generators by Thomas et. al.
Summary: Wallace-method is the fastest, but doesn't provide good statistical quality. Ziggurat was the second fastest method among the huge benchmark while passing the chi-squared and achieving good scores at the high sigma test (see table 3 and 4).

@9il

This comment has been minimized.

Member

9il commented Aug 5, 2016

btw an interesting overview paper is Gaussian Random Number Generators by Thomas et. al.
Summary: Wallace-method is the fastest, but doesn't provide good statistical quality. Ziggurat was the second fastest method among the huge benchmark while passing the chi-squared and achieving good scores at the high sigma test (see table 3 and 4).

Wallace-method is not specialised for Normal if I am not wrong

@wilzbach

This comment has been minimized.

Member

wilzbach commented Aug 5, 2016

Wallace-method is not specialised for Normal if I am not wrong

Afaik Ziggurat is neither. In the literature it's just commonly only used for Normal and exponential distributions, however it's a general method that works for all monotone decreasing distributions (or if their symmetric half is monotone decreasing)

@joseph-wakeling-sociomantic

This comment has been minimized.

joseph-wakeling-sociomantic commented Aug 5, 2016

First question is Ziggurat better for basic distributions than basic implementations of them?

I would say, "Yes, but." The "but" is because Ziggurat is more complicated to implement correctly (as we're learning). So, in terms of the current project, I would say there's a tradeoff between focusing on Ziggurat correct, versus getting a good variety of basic distributions in place with simpler (but more limited) algorithms.

@joseph-wakeling-sociomantic

This comment has been minimized.

joseph-wakeling-sociomantic commented Aug 5, 2016

however it's a general method that works for all monotone decreasing distributions (or if their symmetric half is monotone decreasing

Yes, exactly.

@wilzbach wilzbach referenced this pull request Aug 10, 2016

Merged

Flex2 - first benchmark #287

@WebDrake

This comment has been minimized.

WebDrake commented Aug 10, 2016

@wilzbach:

We have time at least time until October.

The emails I'm getting from GSoC suggest that we're supposed to be finished up by the end of August, with 23 August as your own deadline for finalizing code and 29 August as the deadline for Ilya and me to submit our final evaluation report?
https://developers.google.com/open-source/gsoc/timeline

@wilzbach

This comment has been minimized.

Member

wilzbach commented Aug 10, 2016

The emails I'm getting from GSoC suggest that we're supposed to be finished up by the end of August, with 23 August as your own deadline for finalizing code and 29 August as the deadline for Ilya and me to submit our final evaluation report?
https://developers.google.com/open-source/gsoc/timeline

Yep, but that doesn't stop me to continue to work (I know that I wasted quite a lot of time)
My submission will only include the (Tin)flex algorithm, but I am still onto the mission to write (the building blocks for) a new & fast std.random for D. Hence I suggested to do it properly and as the benchmark in #286 suggests it's worth it.

@WebDrake

This comment has been minimized.

WebDrake commented Aug 10, 2016

Yep, but that doesn't stop me to continue to work (I know that I wasted quite a lot of time)

It's great that you want to keep working, but I was concerned about the expectations raised in the description of your GSoC project and what the people responsible might expect to see (which is why earlier I raised the possibility of doing some basic implementations of a variety of non-uniform distributions).

@9il you're the primary mentor here, so what are your thoughts?

@9il

This comment has been minimized.

Member

9il commented Aug 10, 2016

@9il you're the primary mentor here, so what are your thoughts?

We already have 2 general purpose Discrete RNG realizations and Tinflex will be ready soon.

Tinflex is hard numeric project without obvious workforce requirements. R version contains a lot of numeric bugs, many of them are fixed in the this project. We can stamp / copy-past boost rng, and this is not problem. But copy-pasting is not related to the proper RNG numbers. First, we need proper shell over std.random and fixed uniform generators. @wilzbach expected that he would add them during GSoC, in the same I didn't expect it.

If @wilzbach implemented a variety of non-uniform distributions but not Tinflex, then I would not be able to consider this GSoC project as finished.

@WebDrake

This comment has been minimized.

WebDrake commented Aug 10, 2016

If @wilzbach implemented a variety of non-uniform distributions but not Tinflex, then I would not be able to consider this GSoC project as finished.

I agree with that; I'm asking whether we should implement a variety of non-uniform distributions instead of (short term) focusing on Ziggurat. Tinflex obviously takes primacy as it is the most significant part of the promised work.

@9il

This comment has been minimized.

Member

9il commented Aug 10, 2016

I agree with that; I'm asking whether we should implement a variety of non-uniform distributions instead of (short term) focusing on Ziggurat. Tinflex obviously takes primacy as it is the most significant part of the promised work.

Yes, I prefer to add more a variety of non-uniform distributions instead of Ziggurat

@9il

This comment has been minimized.

Member

9il commented Nov 27, 2016

Please reopen for mir-random

@9il 9il closed this Nov 27, 2016

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment