Ziggurat algorithm for normal deviates

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

Ziggurat algorithm for normal deviates

Kakhkhor Abdijalilov
Did anyone implement Ziggurat algorithm already?

I converted one from GSL into C++. In my tests I paired Mersenne
Twister with Acklam's approximation, boost's inverse CDF and Ziggurat
to generate standard normal distribution. Ziggurat method was 3 times
faster than Acklam's approximation and 7 times faster than using
boost's inverse CDF. GSL is licensed under GPL v3, but I adopted only
"matematical formula", the implementation follows the original
Ziggurat paper. My implementation is hardwired to use Mersenne Twister
directly and doesn't rely on GSL. GSL version uses smaller look up
table and exponential distribution for the tails. I might create my
own lookup table, but  it would be exactly the same. The math formula
is from the published paper, not from GSL. Given that the only thing I
took from GSL is precomputed look up table, is my implementation
covered by GPL? What do you say?

------------------------------------------------------------------------------
ThinkGeek and WIRED's GeekDad team up for the Ultimate
GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the
lucky parental unit.  See the prize list and enter to win:
http://p.sf.net/sfu/thinkgeek-promo
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev
Reply | Threaded
Open this post in threaded view
|

Re: Ziggurat algorithm for normal deviates

Ferdinando M. Ametrano-2
On Tue, Jun 15, 2010 at 5:07 PM, Kakhkhor Abdijalilov
<[hidden email]> wrote:
> In my tests I paired Mersenne
> Twister with Acklam's approximation, boost's inverse CDF and Ziggurat
> to generate standard normal distribution.

for those who wonder Acklam's approximation is the (best) inverse CDF
implemented in QuantLib

> Ziggurat method was 3 times
> faster than Acklam's approximation and 7 times faster than using
> boost's inverse CDF.

ok, so Ziggurat is interesting, even if it is only for pseudo-random
and not for quasi-random

> GSL version uses smaller look up
> table and exponential distribution for the tails.

smaller than what ? what is this lookup table ?

> I might create my
> own lookup table, but  it would be exactly the same.

of course, if it is the inverse cdf lookup table it's no problem to
generate it again (even with Acklam improved by using one iteration of
Halley's rational method in order to obtain full machine precision)

> GSL is licensed under GPL v3, but I adopted only
> "matematical formula", the implementation follows the original
> Ziggurat paper. My implementation is hardwired to use Mersenne Twister
> directly and doesn't rely on GSL.
> [...] The math formula
> is from the published paper, not from GSL. Given that the only thing I
> took from GSL is precomputed look up table, is my implementation
> covered by GPL?

I would recompute the lookup table and then it's not covered by GPL

ciao -- Nando

------------------------------------------------------------------------------
ThinkGeek and WIRED's GeekDad team up for the Ultimate
GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the
lucky parental unit.  See the prize list and enter to win:
http://p.sf.net/sfu/thinkgeek-promo
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev
Reply | Threaded
Open this post in threaded view
|

Re: Ziggurat algorithm for normal deviates

Kakhkhor Abdijalilov
OK, it is done. I submitted my implementation to SourceForge patch
system (Ziggurat Algorithm by renorm). The underlying RNG is
MersenneTwisterUniformRng (MT19937). In my tests Ziggurat was about
2-3 times faster than MT19937+Acklam and about 7 times faster than
MT19937+boot's quantile function. I tested Ziggurat with Diehard
battery and everything looked OK.

I also did some changes to MersenneTwisterUniformRng and
InverseCumulativeNormal (Acklam's method). On my platform
MersenneTwisterUniformRng  is slightly faster than boost's mt19937 and
considerable faster that the old MersenneTwisterUniformRng, especially
in debug mode. Apparently, using STL vector instead of plain array
makes a difference in this case.

InverseCumulativeNormal features static std_ic (standard normal
inverse cumulative) method. It computes inverse CDF with mean=0 and
sigma=0. Using std_ic instead of operator() saves us 2 floating point
operation. The performance difference is quite noticeable.

Regards,
Kakhkhor Abdijalilov.

------------------------------------------------------------------------------
ThinkGeek and WIRED's GeekDad team up for the Ultimate
GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the
lucky parental unit.  See the prize list and enter to win:
http://p.sf.net/sfu/thinkgeek-promo
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev
Reply | Threaded
Open this post in threaded view
|

Re: Ziggurat algorithm for normal deviates

Luigi Ballabio
On Wed, 2010-06-16 at 18:22 -0500, Kakhkhor Abdijalilov wrote:
> OK, it is done. I submitted my implementation to SourceForge patch
> system

Ok, thanks.  I'll have a look.

Luigi


--

A debugged program is one for which you have not yet found the
conditions that make it fail.
-- Jerry Ogdin



------------------------------------------------------------------------------
ThinkGeek and WIRED's GeekDad team up for the Ultimate
GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the
lucky parental unit.  See the prize list and enter to win:
http://p.sf.net/sfu/thinkgeek-promo
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev
Reply | Threaded
Open this post in threaded view
|

Re: Ziggurat algorithm for normal deviates

Kakhkhor Abdijalilov
I implemented it from the scratch. No copyright problems, whatsoever.
It follows the recipe from Marsaglia and Tsang's paper, except the
sampling from the tail. My implementation uses the conditional
distribution to sample from the tail of normal density.

Regards,
Kakhkhor Abdijalilov.

------------------------------------------------------------------------------
Beautiful is writing same markup. Internet Explorer 9 supports
standards for HTML5, CSS3, SVG 1.1,  ECMAScript5, and DOM L2 & L3.
Spend less time writing and  rewriting code and more time creating great
experiences on the web. Be a part of the beta today.
http://p.sf.net/sfu/beautyoftheweb
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev