In eddelbuettel/rcppziggurat: 'Rcpp' Integration of Different "Ziggurat" Normal RNG Implementations

library(RcppZiggurat)
## this RData file contains the pre-run results 
##
##   [1] "chires"   "norres"   "rspeed"   "stdres"   "zigspeed"
##
## from the two speed comparions (ie rspeed and zigspeed) as well the three 
## simulation results (ie norres, stdres, and chires) 
##
## code for the two speed tests is shown below, the simulations results were run as
##
##   #!/bin/bash
##   Rscript -e 'library(RcppZiggurat); options("mc.cores"=6L); resS <- RcppZiggurat:::standardTest(5e9); saveRDS(resS, "standardTest.rds")' &
##   Rscript -e 'library(RcppZiggurat); options("mc.cores"=6L); resC <- RcppZiggurat:::chisqTest(1e10); saveRDS(resC, "chisqTest.rds")' &
##   Rscript -e 'library(RcppZiggurat); options("mc.cores"=6L); resN <- RcppZiggurat:::normalTest(5e9); saveRDS(resN, "normalTest.rds")' &
##
load("RcppZiggurat.RData")

Introduction

Generating random number for use in simulation is a classic topic in scientific computing and about as old as the field itself. Most algorithms concentrate on the uniform distribution. Its values can be used to generate randomly distributed values from other distributions simply by using the relevant inverse function.

Regarding terminology, all computational algorithms for generation of \textsl{random} numbers are by definition deterministic. Here, we follow standard convention and refer to such numbers as \textsl{pseudo-random} as they can always be recreated given the seed value for a given sequence. We are not concerned in this note with \textsl{quasi-random} numbers (also called low-discrepnancy sequences). We consider this topic to be subset of pseudo-random numbers subject to distributional constraints. Without lack of generality, we will henceforth drop the qualifier \textsl{pseudo} when refering to random numbers.

Due to its importance for many modeling tasks, the Normal distribution has also been studied extensively in order to find suitable direct algorithms which generate stream of normally distributed (pseudo) random numbers. A well-known examples for such algorithms includes the method by Box and Muller \citep{Box+Muller:1958}. A useful recent survey of the field is provided by \cite{Thomas+Luk+Leong+Villasenor:2007}.

In an important paper, \cite{Marsaglia+Tsang:2000} introduced the Ziggurat method. Since its initial publication, this algorithm has become reasonably popular\footnote{Usage of a code search engine such as \url{code.ohloh.net} or \url{codesearch.debian.net} provides a good approximation to the popularity of the Ziggurat algorithm as its name is also a reasonably unique search term within the field of computing.} as it provides a very useful combination of both excellent statistical properties and execution speed. \cite{Thomas+Luk+Leong+Villasenor:2007} conclude their survey by saying that \begin{quote} the Ziggurat method, the second in speed, is about 33\% slower than the [fastest] method but does not suffer from correlation problems. Thus, when maintaining extremely high statistical quality is the first priority, and subject to that constraint, speed is also desired, the Ziggurat method will often be the most appropriate choice. \end{quote}

This paper reexamines the Ziggurat method, provides a new and portable C++ implementation, and compares it to several other Open Source implementations of the underlying algorithm by applying three different statistical tests.

Ziggurat

This sections briefly discusses the key papers related to Ziggurat.

Marsaglia and Tsang

\cite{Marsaglia+Tsang:2000} introduced the Ziggurat method. It provides a fast algorithm for generating both normally and exponentially distributed random numbers. The original paper also contains a corresponding implementation in the C language.

The listing in Figure \ref{fig:ZiggMT} shows this initial implementation. We have removed the code for generating exponentially distributed random numbers as well a comment header from the listing to keep the display more compact. The full listing is also included in the \pkg{RcppZiggurat} package for reference.

\begin{figure*} \begin{center} \begin{small} \begin{Shaded} \begin{Highlighting}[] \PreprocessorTok{#include }\ImportTok{} \AttributeTok{static} \DataTypeTok{unsigned} \DataTypeTok{long} \NormalTok{jz,jsr=}\DecValTok{123456789}\NormalTok{;}

\PreprocessorTok{#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)} \PreprocessorTok{#define UNI (.5 + (signed) SHR3*.2328306e-9)} \PreprocessorTok{#define IUNI SHR3}

\AttributeTok{static} \DataTypeTok{long} \NormalTok{hz;} \AttributeTok{static} \DataTypeTok{unsigned} \DataTypeTok{long} \NormalTok{iz, kn[}\DecValTok{128}\NormalTok{], ke[}\DecValTok{256}\NormalTok{];} \AttributeTok{static} \DataTypeTok{float} \NormalTok{wn[}\DecValTok{128}\NormalTok{],fn[}\DecValTok{128}\NormalTok{], we[}\DecValTok{256}\NormalTok{],fe[}\DecValTok{256}\NormalTok{];}

\PreprocessorTok{#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix())}

\CommentTok{/ nfix() generates variates from the residue when rejection in RNOR occurs. /} \DataTypeTok{float} \NormalTok{nfix(}\DataTypeTok{void}\NormalTok{)} \NormalTok{{} \AttributeTok{const} \DataTypeTok{float} \NormalTok{r = }\FloatTok{3.}\ErrorTok{442620f}\NormalTok{; }\CommentTok{/ The start of the right tail /} \AttributeTok{static} \DataTypeTok{float} \NormalTok{x, y;} \ControlFlowTok{for}\NormalTok{(;;)} \NormalTok{{ x=hzwn[iz]; }\CommentTok{/ iz==0, handles the base strip /} \ControlFlowTok{if}\NormalTok{(iz==}\DecValTok{0}\NormalTok{)} \NormalTok{{ }\ControlFlowTok{do}\NormalTok{{ x=-log(UNI)}\FloatTok{0.2904764}\NormalTok{; y=-log(UNI);} }\CommentTok{/ .2904764 is 1/r /} \ControlFlowTok{while}\NormalTok{(y+y}\DecValTok{0}\NormalTok{)? r+x : -r-x;} \NormalTok{}} \CommentTok{/ iz>0, handle the wedges of other strips /} \ControlFlowTok{if}\NormalTok{( fn[iz]+UNI(fn[iz}\DecValTok{-1}\NormalTok{]-fn[iz]) < exp(-.}\DecValTok{5}\NormalTok{x*x) ) }\ControlFlowTok{return} \NormalTok{x;}

 \CommentTok{/* initiate, try to exit for(;;) for loop*/}
  \NormalTok{hz=SHR3;}
  \NormalTok{iz=hz&}\DecValTok{127}\NormalTok{;}
  \ControlFlowTok{if}\NormalTok{(fabs(hz)<kn[iz]) }\ControlFlowTok{return} \NormalTok{(hz*wn[iz]);}

\NormalTok{}} \NormalTok{}}

\CommentTok{/--------This procedure sets the seed and creates the tables------/} \DataTypeTok{void} \NormalTok{zigset(}\DataTypeTok{unsigned} \DataTypeTok{long} \NormalTok{jsrseed)} \NormalTok{{ }\AttributeTok{const} \DataTypeTok{double} \NormalTok{m1 = }\FloatTok{2147483648.0}\NormalTok{, m2 = }\DecValTok{4294967296}\NormalTok{.;} \DataTypeTok{double} \NormalTok{dn=}\FloatTok{3.442619855899}\NormalTok{,tn=dn,vn=}\FloatTok{9.91256303526217e-3}\NormalTok{, q;} \DataTypeTok{double} \NormalTok{de=}\FloatTok{7.697117470131487}\NormalTok{, te=de, ve=}\FloatTok{3.949659822581572e-3}\NormalTok{;} \DataTypeTok{int} \NormalTok{i;} \NormalTok{jsr^=jsrseed;}

\CommentTok{/ Set up tables for RNOR /} \NormalTok{q=vn/exp(-.}\DecValTok{5}\NormalTok{dndn);} \NormalTok{kn[}\DecValTok{0}\NormalTok{]=(dn/q)*m1;} \NormalTok{kn[}\DecValTok{1}\NormalTok{]=}\DecValTok{0}\NormalTok{;}

\NormalTok{wn[}\DecValTok{0}\NormalTok{]=q/m1;} \NormalTok{wn[}\DecValTok{127}\NormalTok{]=dn/m1;}

\NormalTok{fn[}\DecValTok{0}\NormalTok{]=}\DecValTok{1}\NormalTok{.;} \NormalTok{fn[}\DecValTok{127}\NormalTok{]=exp(-.}\DecValTok{5}\NormalTok{dndn);}

\ControlFlowTok{for}\NormalTok{(i=}\DecValTok{126}\NormalTok{;i>=}\DecValTok{1}\NormalTok{;i--)}
\NormalTok{\{dn=sqrt(}\DecValTok{-2}\NormalTok{.*log(vn/dn+exp(-.}\DecValTok{5}\NormalTok{*dn*dn)));}
 \NormalTok{kn[i}\DecValTok{+1}\NormalTok{]=(dn/tn)*m1;}
 \NormalTok{tn=dn;}
 \NormalTok{fn[i]=exp(-.}\DecValTok{5}\NormalTok{*dn*dn);}
 \NormalTok{wn[i]=dn/m1;}
\NormalTok{\}}

\NormalTok{}} \end{Highlighting} \end{Shaded} \end{small} \caption{Ziggurat code by \cite{Marsaglia+Tsang:2000}. \label{fig:ZiggMT}} \end{center} \end{figure*}

As can be seen from Figure \ref{fig:ZiggMT}, the Ziggurat algorithm is implemented in bare-bones C code using a number of macros, and two helper functions. The helper functions allow setting a seed, and deal with parameter updates needed in about 2.5\% of cases. The actual core component---the function to draw a random number distributed according to thstandard normal distribution---is provided by the macro \texttt{RNOR}. Needless to say, using C macros is no longer cosidered \textsl{de rigeur}. Possible side effects include inadvertent changes in globally visible variables, as well as possible bugs from the macro evaluation.

A more important concern is that this implementation uses \texttt{unsigned long} types, and explicit bit mapping operations. The code was originally developed for 32-bit operating systems where \texttt{int} and \texttt{long} are typically four bytes (or 32 bits) wide. Hence the code does not produce correct results on a 64-bit operating system as (signed and unsigned) \texttt{long} types are typically eight bytes (or 64 bits) wide (whereas \texttt{int} is still 32 bits).

Our modified version introduced below overcomes both issues.

Leong, Zhang et al

\cite{Leong+Zhang+Lee+Luk+Villasenor:2005} show in a comment that the Ziggurat method as introduced by \cite{Marsaglia+Tsang:2000} suffers from another weakness due to the \texttt{SHR3} generator (by Marsaglia). The authors show via a $\chi^2$-test that the generator has a short period of about $2^{32}-1$, or the four byte limit. Replacing it with the \texttt{KISS} generator (also by Marsaglia) improves the performance beyond this limit.

```{Rcpp, eval=FALSE, fip.cap="ABC"}

define MWC ((znew<<16)+wnew )

define SHR3 (jz=jsr, jsr^=(jsr<<13), \

    jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)

define CONG (jcong=69069*jcong+1234567)

define KISS ((MWC^CONG)+SHR3)

define RNOR (hz=KISS, iz=hz&127, \

    (fabs(hz)<kn[iz]) ? hz*wn[iz] : nfix())

<!-- \caption{Improved Ziggurat code by Leong et al}\label{fig:ZiggLZLLV} -->

Following \cite{Leong+Zhang+Lee+Luk+Villasenor:2005}, Ziggurat code should
use an improved uniform generator. Choices are either KISS as suggested
initially, or another trusted (and fast) uniform generator such as the
Mersenne Twister \citep{Matsumoto+Nishimura:1998}. Other Open Source
implementations (such as the ones discussed below) frequently use the
Mersenne Twister as the souce of uniformly distributed random numbers.

## Voss's implementation in GNU GSL

\cite{Voss:2011} provided another Ziggurat implementation for use in the GNU
Scientific Library or GSL \citep{GSL}. It uses the Mersenne Twister generator
by \cite{Matsumoto+Nishimura:1998} which also avoids the issue identified by
\cite{Leong+Zhang+Lee+Luk+Villasenor:2005} in which the originally-used
uniform generator had too short a cycle.

\cite{Voss:2011} notes two differences between his implementation and the
original work by \cite{Marsaglia+Tsang:2000}. First, he uses only 128 instead
of 256 steps which reduces the memory requirements and computational
cost at a possible (though presumably minor) loss of precision.
Second, he uses an exponential distribution with tail rejection
for the base strip, which is motivated by a simpler implementations.  Both of
these aspects could have implications for the statistical properties of the
generator. Voss also appears to be unaware of the work by
\cite{Leong+Zhang+Lee+Luk+Villasenor:2005}, yet sidesteps the issue raised by
these author by relying on the Mersenne Twister generator.

Here, the implementation from the current GSL sources and file
\code{randist/gausszig.c} is used, and adapted to the class strucuture
detailed in section \ref{sec:cpp}.


## Gretl

The Gretl \citep{Cottrell+Lucchetti:Gretl} econometrics program contains
another Open Source implementation of the Ziggurat algorithm. The code
credits the implementation by \cite{Voss:2011} described above.
\cite{Yalta+Schreiber:2012} review the Gretl implementation and performance
of Ziggurat and find it to be satisfactory.

We use the implementation from the file \code{src/lib/random.c} from the
current Gretl sources and adapt to the class strucuture detailed in
section \ref{sec:cpp}.


## QuantLib

The QuantLib library for quantitative finance
\citep{Ametrano+Ballabio+et+al:QuantLib} contains another open source
implementation of Ziggurat. It is provided in the files
\code{ql/experimental/math/zigguratrng.hpp} and
\code{ql/experimental/math/zigguratrng.cpp}. As part of the experimental
section, it is made available for further study and use, but not yet part of
the default build.  As before, we integrate this source into the class
structure used here.


# Speed

## R Generators

Before comparing the speed of the different Ziggurat implementations, it is
also illustrative to compare the different R generators.
Figure \ref{fig:RNormalRNGs} provides a comparison.

<!-- this chunk generated the data plotted in the next chunk-->
```r
library(microbenchmark)
res <- microbenchmark({RNGkind(,"Kinderman-Ramage"); rnorm(1e6)},
                      {RNGkind(,"Ahrens-Dieter"); rnorm(1e6)},
                      {RNGkind(,"Box-Muller"); rnorm(1e6)},
                      {RNGkind(,"Inversion"); rnorm(1e6)},
                      times=100)
levels(res$expr) <- c("KR", "AH", "BM", "Inv")
#saveRDS(res, file="~/git/rcppziggurat/vignettes/Rspeed.rds")

library(lattice)
rdf <- as.data.frame(rspeed)
rdf[,1] <- ordered(rdf[,1], levels=c("AH","KR","Inv","BM"), labels=c("AH","KR","Inv","BM"))
#pdf("plot1.pdf", 8, 8)
bwplot(time/1e6 ~ expr, rdf,
       ylab="Time in msec",
       main="Time for 100 times 1e6 normal draws",
       lattice.options=list(fontsize=list(text=8, points=8)),
       panel=function(...,box.ratio) { panel.violin(..., col="lightgray",
                                                    varwidth=FALSE, box.ratio=box.ratio)} )
                                        #dev.off()

\begin{figure}[!tp] \begin{minipage}{0.475\textwidth} \begin{center} \includegraphics[width=0.95\textwidth]{figures/timingR} \caption{R Normal RNG Generator Performance \label{fig:RNormalRNGs}} \begin{minipage}{0.8\textwidth} \scriptsize Note: Figure shows timings from the \pkg{microbenchmark} package using 100 replications of 1,000,000 draws per generator. Code for the figure is included in the \pkg{RcppZiggurat} package, and the source code for this document. \end{minipage} \end{center} \end{minipage} \begin{minipage}{0.475\textwidth} \begin{center} \includegraphics[width=0.95\textwidth]{figures/timingZiggurat} \caption{Ziggurat and R Normal RNG Generator Performance \label{fig:ZiggRNormalRNGs}} \begin{minipage}{0.8\textwidth} \scriptsize Note: Figure shows timings from the \pkg{microbenchmark} package using 100 replications of 1,000,000 draws per generator. Code for the figure is included in the \pkg{RcppZiggurat} package, and the source code for this document. \end{minipage} \end{center} \end{minipage} \end{figure}

We see that the Box-Muller generator is the slowest by some margin. However, both Kinderman-Ramage and Ahrens-Dieter are faster than the Inversion method chosen as the default in R. So even before considering Ziggurat generators, R users could reap a speed benefit simply by calling \code{RNGkind(,"Ahrens-Dieter")} or \code{RNGkind(,"Kinderman-Ramage")}.

Ziggurat Generators

library(microbenchmark)
library(lattice)
library(RcppZiggurat)
N <- 1e6
res <- microbenchmark(rnorm(N),          # Marsgalia and Tsang, JSS, 2000
                      zrnorm(N),         # based on updated Burkardt implementation
                      zrnormGSL(N),      # GSL's ziggurat by Voss
                      zrnormQL(N),       # QuantLib variant
                      zrnormGl(N),       # Gretl
                      times=100)
levels(res$expr) <- c("RInv", "Zigg", "ZiggGSL", "ZiggQL", "ZiggGretl")
#saveRDS(res, file="~/git/rcppziggurat/vignettes/Zigspeed.rds")

library(lattice)
#res <- readRDS("~/git/rcppziggurat/vignettes/Zigspeed.rds")
zdf <- as.data.frame(zigspeed)
zdf[,1] <- ordered(zdf[,1],
                   levels=c("Zigg", "ZiggGSL", "ZiggQL", "ZiggGretl", "RInv"),
                   labels=c("Zigg", "ZiggGSL", "ZiggQL", "ZiggGretl", "RInv"))
print(bwplot(time/1e6 ~ expr, zdf, ylab="Time in msec",
             main="Time for 100 times 1e6 normal draws",
             panel=function(...,box.ratio) {
                 panel.violin(..., col="lightgray", varwidth=FALSE, box.ratio=box.ratio)
}))

All Ziggurat generators are significantly faster the the default generator in R which uses a inversion method. Among the Ziggurat generators, we notice that approaches using an external uniform number generator (GNU GSL, GNU Gretl, QuantLib) are all slower than our compact and self-contained implementation which is seen as the fastest method.
<!--- %% A comparison to the %% original implementation and the comment follow-up (not shown here but easily %% programmable using the function in the \pkg{RcppZiggurat} package) shows them %% to be of comparable speed.

%% However, as the next section details, the initial implementation should not %% be used. And as discussed below, the bare-bones implementation of the initial %% proposals lack some safeguards suggested by modern software practice such as %% encapsulation within a class, a namespace or both. -->

Accuracy

Standard Test for Uniform RNG draws

Test for random number generators are often focussed on the case of uniform generators which are the most common type of generators. As detailed for example in \cite{Brown+Eddelbuettel+Bauer:DieHarder}, a test proceeds as follow:

\begin{enumerate} \item Take $n$ draws from a $U(0,1)$ distribution (as any given $U(a,b)$ can always be scaled to $U(0,1)$), and then compute the sum of the $n$ values. \item Repeat this $m$ times to create a set of $m$ sums of uniform RNG draws. \item With $n$ large enough, the collection of $m$ results will converge towards normally distributed random variable with a mean of $n/2$ and a standard deviation of $\sqrt{n/12}$ (which is the Irwin-Hall distribution of the sum of uniformly distributed values). \item Given this asymptotic result, one can construct a probability value $p_i$ for each of the $m$ values using the inverse of the Normal distribution using the known mean and standard deviation from the Irwin-Hall distribution. \item We now have $m$ uniformly distributed values. A standard test such as Kolmogorov-Smirnov or Wilcoxon can be used to test against departures from the uniform distribution. \end{enumerate}

Here, we can apply this test for first converting the $N(0,1)$ distributed values produced by the given Ziggurat implementation to $U(0,1)$ distributed values by using the inverse of the normal distribution. We are then ready to simulate and test. Figure \ref{fig:StandardTest} below displays Q-Q plots for the empirical distribution against the uniform, and displays the $p$-values of a Kolmogorov-Smirnov as well as a Wilcoxon test.

RcppZiggurat:::plotAll(stdres)

\begin{figure}[!tp] \begin{minipage}{0.475\textwidth} \begin{center} \includegraphics[width=0.95\textwidth]{figures/standardResults} \caption{Standard Test applied to Ziggurat generators \label{fig:StandardTest}} \begin{minipage}{0.8\textwidth} \scriptsize Note: Code for the figure is included in the \pkg{RcppZiggurat} package, and the source code for this document. \end{minipage} \end{center} \end{minipage} \begin{minipage}{0.475\textwidth} \begin{center} \includegraphics[width=0.95\textwidth]{figures/normalResults} \caption{Normal Test applied to Ziggurat generators \label{fig:NormalTest}} \begin{minipage}{0.8\textwidth} \scriptsize Note: Code for the figure is included in the \pkg{RcppZiggurat} package, and the source code for this document. \end{minipage} \end{center} \end{minipage} \end{figure}

We see that five of the six generators pass the test. In the case of the original \cite{Marsaglia+Tsang:2000} generator, we can see the departure from the expected diagonal clearly once we draw more than $4.2 \times 10^9$ numbers (which is the limit of the representation of an unsigned four-byte nunber). However, only the Kolmogorov-Smirnov test can formally reject; the Wilcoxon test appear to lack sufficient power in this setting. The QuantLib generator is seen as suspicious which $p$-value just below a conventional rejection level.

Normal Test for Normal RNG draws

We can propose a simpler variant of the test outlined in the previous section. As the random numbers we are drawing are following a $N(0,1)$ distribution, the sum of their values follows a $N(0, n)$ distribution. This allows us to skip one inversion step:

\begin{enumerate} \item Take $n$ draws from a $N(0,1)$ distribution and then compute the sum of the $n$ values. \item Repeats this $m$ times to create a set of $m$ sums of (standard) normals RNG draws. \item The collection of the $m$ sums of $n$ normals converges towards a mean of $0$ and a standard deviation of $\sqrt{n}$. \item Given this known result, one can construct a probability value $p_i$ for each of the $m$ values using the inverse of the Normal distribution using the known mean and standard deviation. \item We again have $m$ uniformly distributed values. A standard test such as Kolmogorov-Smirnov or Wilcoxon can be used to test against departures from the uniform distribution. \end{enumerate}

RcppZiggurat:::plotAll(norres)

Results, shown in Figure \ref{fig:NormalTest}, are qualitatively similar to the result discussed above. Kolmogorov-Smirnov rejects for the \cite{Marsaglia+Tsang:2000} generator. However, we note that the Wilcoxon test now has a lower $p$-value---we would now reject at conventional test levels. The QuantLib implementation is now rejected by the Wilcoxon test but not the Kolmogorov-Smirnov.

$\chi^2$ test

Another test variant is the $\chi^2$ test which was also used by \cite{Leong+Zhang+Lee+Luk+Villasenor:2005}. The basic idea is as follow: \begin{enumerate} \item The real line is divided into $B$ bins, equally spaced between (symmetric) values distant enough from zero so that no $N(0,1)$ draw should exceed them. \item Here, we follow \cite{Leong+Zhang+Lee+Luk+Villasenor:2005} and use a range from -7 to 7 with a total of 200 bins. \item A large number of $N(0,1)$ random variates is drawn, and for each of these numbers a counter in the bin corresponding to the draw is increased. \item After the $N$ draws, the empirical distribution is compared to the theoretical (provided by the corresponding value of the Normal density function) using a standard $chi^2$ test. \end{enumerate}

As can be seen in Figure \ref{fig:ChiSqTest}, the original proposal by \cite{Marsaglia+Tsang:2000} fails as was shown by \cite{Leong+Zhang+Lee+Luk+Villasenor:2005}. All other tests pass again.

RcppZiggurat:::plotChiSq(chires)

\begin{figure}[!tp] \begin{center} \includegraphics[width=0.9\linewidth]{figures/chisquareResults} \caption{$\chi^2$ Test applied to Ziggurat generators \label{fig:ChiSqTest}} \medskip \begin{minipage}{0.85\linewidth} \scriptsize Note: Code for the figure is included in the \pkg{RcppZiggurat} package, and the source code for this document. \end{minipage} \end{center} \end{figure}

C++ Implementation

\label{sec:cpp}

Preceding work by \cite{Marsaglia+Tsang:2000} and \cite{Leong+Zhang+Lee+Luk+Villasenor:2005} also contained implementations in the C language. These versions were implemented in just a few lines, and used idioms common to C programmers such as macros and global variables.

C++ programming style permits encapsulation in order to avoid possible collisions and side-effects. Moreover, by using a modest amount of object-oriented programming we can use a class structure with a common base class to express commonalities between the implementations. The following code segment shows the virtual base class used here.

```{Rcpp, eval=FALSE}

ifndef RcppZiggurat__Zigg_h

define RcppZiggurat__Zigg_h

include

include // or cstdint (C++11)

namespace Ziggurat {

class Zigg { public: virtual ~Zigg() {}; virtual void setSeed(const uint32_t s) = 0; // no getSeed() as GSL has none virtual double norm() = 0; }; }

endif

<!-- FIXME  \caption{Ziggurat Base Class in package \pkg{RcppZiggurat}. \label{fig:ZigguratBase}}-->
<!-- FIXME  \end{center}
\end{figure}-->

As shown in the preceding code segment, we provide two user-accessible
functions to obtain a normal random deviate, and to set the seed,
respectively. The actual implementation uses portable types such as
\code{uint32\_t}, an unsigned 32-bit integer provided by the C header file
\code{stdint.h}, which provides correct and identical results on both 32-bit
and 64-bit operating systems.

Each actual implementation can then encapsulate its state variable as a
private variable inaccessible to other functions. Such a small core to each
class also makes it feasible to provide a Ziggurat generator in each thread
in a parallel execution framework.

Our Ziggurat implementation has no external dependencies and can therefore be
used in other projects. The testing framework used for this note has a single
dependency on the GNU GSL as the generator by \cite{Voss:2011} is used via
its GSL implementations.  The generator and testing framework in the
corresponding R package have a build-dependency on R, and are of course
accessed by R. But the generator discussed here could equally well be used in
standalone programs or with other scripting languages.

\section{R Integration}

In the \pkg{RcppZiggurat} package, the Rcpp Attributes
\citep{CRAN:Rcpp:Attributes} feature of the \pkg{Rcpp} C++/R integration
package \citep{CRAN:Rcpp,Eddelbuettel:2013:Rcpp,TAS:Rcpp} are used to access instances of
the corresponding generator class.

<!-- FIXME\begin{figure}[!tp]
  \begin{center}-->
```cpp

#include <Rcpp.h>

#include <Ziggurat.h>

static Ziggurat::Ziggurat::Ziggurat zigg;

// [[Rcpp::export]]
Rcpp::NumericVector zrnorm(int n) {
    Rcpp::NumericVector x(n);
    for (int i=0; i<n; i++) {
        x[i] = zigg.norm();
    }
    return x;
}

// [[Rcpp::export]]
void zsetseed(unsigned long int s) {
    zigg.setSeed(s);
    return;
}

In this particular reference implementation, we have chosen a namespace \code{Ziggurat} for the entire project. Within this namespace, we opted to provided an extra namespace layer for each generator as some of these generators still use global variables---which are therefore shielded in their own namespace. For example, for the \cite{Marsaglia+Tsang:2000} generator, we use \code{Ziggurat::ZigurratMT}. Next is the name of the actual class---which in the case of the reference implementation shown above is once again \code{Ziggurat} leading to the triple use of the term. Actual deployment of the Ziggurat generator (without comparison to other implementations and concerns about interaction between variables, particularly for the older implementations having global variables) can of course be used with a single namespace.

The remainder shows how two functions \code{znorm()} and \code{zsetseed()} are provided via the attribute \code{[[Rcpp::Export]]} as described in the Rcpp Attributes vignette \citep{CRAN:Rcpp:Attributes} of the \pkg{Rcpp} package \citep{CRAN:Rcpp,Eddelbuettel:2013:Rcpp,TAS:Rcpp}.

Conclusion

This note describes the \pkg{RcppZiggurat} package and its new implementation of the Ziggurat generator for normally distributed random numbers. The package is implemented in a way which is portable so that it can be used on 32-bit and 64-bit operating systems, filling a gap left by the original implementations \citep{Marsaglia+Tsang:2000,Leong+Zhang+Lee+Luk+Villasenor:2005}.

By embedding the code in a simple C++ class structure, we can ease testing and comparison of several variants of the algorithm. Our note reconfirmed the findings by \cite{Leong+Zhang+Lee+Luk+Villasenor:2005} of a short cycle due to to the use of an inferior uniform generator. Replacing the generator leads to better performance.

We suggest a new test for random number generators producing deviates distributed according to the standard normal distribution by adapting and simplifying an existing test framework for uniform deviates. Both tests confirm the (previously documented) failure of the original Ziggurat proposal but do not find a problem with any of the new implementations (apart from the still-experimental QuantLib generator).

A key motivation for this work has been a desire to improve the speed of creating standard-normally distributed random numbers in R. We find Ziggurat to be faster than the existing implementations, and hope that this generator will be of use to those generating large numbers of draws.

eddelbuettel/rcppziggurat documentation built on April 23, 2024, 7:13 p.m.