Hard integrals arise frequently in reinforcement learning, for example when computing expectations in policy evaluation and policy iteration. They are often analytically intractable and typically estimated with Monte Carlo methods, whose sampling contributes to high variance in policy values and gradients. In this work, we propose to replace Monte Carlo samples with low-discrepancy point sets. We combine policy gradient methods with Randomized Quasi-Monte Carlo, yielding variance-reduced formulations of policy gradient and actor-critic algorithms. These formulations are effective for policy evaluation and policy improvement, as they outperform state-of-the-art algorithms on standardized continuous control benchmarks. Our empirical analyses validate the intuition that replacing Monte Carlo with Quasi-Monte Carlo yields significantly more accurate gradient estimates.

Estimating the unknown density from which a given independent sample originates is more difficult than estimating the mean in the sense that, for the best popular nonparametric density estimators, the mean integrated square error converges more slowly than at the canonical rate of [Formula: see text]. When the sample is generated from a simulation model and we have control over how this is done, we can do better. We examine an approach in which conditional Monte Carlo yields, under certain conditions, a random conditional density that is an unbiased estimator of the true density at any point. By averaging independent replications, we obtain a density estimator that converges at a faster rate than the usual ones. Moreover, combining this new type of estimator with randomized quasi–Monte Carlo to generate the samples typically brings a larger improvement on the error and convergence rate than for the usual estimators because the new estimator is smoother as a function of the underlying uniform random numbers. Summary of Contribution: Stochastic simulation is commonly used to estimate the mathematical expectation of some output random variable X together with a confidence interval for this expectation. But the simulations usually provide information to do much more, such as estimating the entire distribution (or density) of X. Histograms are routinely provided by standard simulation software, but they are very primitive density estimators. Kernel density estimators perform better, but they are trickier to use, have bias, and their mean square error converges more slowly than the canonical rate of O(1/n) with n independent samples. In this paper, we explain how to construct unbiased density estimators that converge at the canonical rate and even much faster when combined with randomized quasi–Monte Carlo. The key idea is to use conditional Monte Carlo to hide appropriate information and obtain a computable (random) conditional density, which acts (under certain conditions) as an unbiased density estimator. Moreover, this sample density is typically smoother than the classic density estimators as a function of the underlying uniform random numbers, so it can get along much better with randomized quasi–Monte Carlo methods. This offers an opportunity to further improve the O(1/n) rate. We observe rates near O(1/n2) on some examples, and we give conditions under which this type of rate provably holds. The proposed approach is simple, easy to implement, and extremely effective, so it provides a significant addition to the stochastic simulation toolbox.

We propose and analyze a generalized splitting method to sample approximately from a distribution conditional on the occurrence of a rare event. This has important applications in a variety of contexts in operations research, engineering, and computational statistics. The method uses independent trials starting from a single particle. We exploit this independence to obtain asymptotic and non-asymptotic bounds on the total variation error of the sampler. Our main finding is that the approximation error depends crucially on the relative variability of the number of points produced by the splitting algorithm in one run, and that this relative variability can be readily estimated via simulation. We illustrate the relevance of the proposed method on an application in which one needs to sample (approximately) from an intractable posterior density in Bayesian inference.

We study the lattice structure of random number generators of the MIXMAX family, a class of matrix linear congruential generators that produces a vector of random numbers at each step. The design of these generators was inspired by Kolmogorov K-systems over the unit torus in the real space, for which the transition function is measure preserving and produces a chaotic behavior. In actual implementations, however, the state space is a finite set of rational vectors, and the MIXMAX has a lattice structure just like linear congruential and multiple recursive generators. Its matrix entries were also selected in a special way to allow a fast implementation, and this has an impact on the lattice structure. We study this lattice structure for vectors of successive and nonsuccessive output values in various dimensions. We show in particular that for coordinates at specific lags not too far apart, in three dimensions, or if we construct points of k+2 or more successive values from the beginning of an output vector of size k, all the nonzero points lie in only two hyperplanes. This is reminiscent of the behavior of lagged-Fibonacci and add-with-carry/subtract-with-borrow generators. And even if we skip the output coordinates involved in this bad structure, other highly structured projections often remain, depending on the choice of parameters. We show that empirical statistical tests can easily detect this structure.

We introduce a new software tool and library named Lattice Builder, written in C++, that implements a variety of construction algorithms for good rank-1 lattice rules. It supports exhaustive and random searches, as well as component-by-component (CBC) and random CBC constructions, for any number of points, and for various measures of (non)uniformity of the points. The measures currently implemented are all shift-invariant and represent the worst-case integration error for certain classes of integrands. They include, for example, the weighted Pα square discrepancy, the Rα criterion, and figures of merit based on the spectral test, with projection-dependent weights. Each of these measures can be computed as a finite sum. For the Pα and Rα criteria, efficient specializations of the CBC algorithm are provided for projection-dependent, order-dependent, and product weights. For numbers of points that are integer powers of a prime base, the construction of embedded rank-1 lattice rules is supported through any of these algorithms, and through a fast CBC algorithm, with a variety of possibilities for the normalization of the merit values of individual embedded levels and for their combination into a single merit value. The library is extensible, thanks to the decomposition of the algorithms into decoupled components, which makes it easy to implement new types of weights, new search domains, new figures of merit, and so on.

We examine some properties of the points produced by certain classes of long-period linear multiple recursive random number generators proposed by L.-Y. Deng and his co-authors in several papers. These generators have their parameters selected in special ways to make the implementation faster. We show that as a result, the points produced by these generators have a poor lattice structure, and a poor initialization of the state can have long-lasting impact, because of the limited diffusion capacity of the recurrence.

We propose a novel simulation-based method that exploits a generalized splitting (GS) algorithm to estimate the reliability of a graph (or network), defined here as the probability that a given set of nodes are connected, when each link of the graph fails with a given (small) probability. For large graphs, in general, computing the exact reliability is an intractable problem and estimating it by standard Monte Carlo methods poses serious difficulties, because the unreliability (one minus the reliability) is often a rare-event probability. We show that the proposed GS algorithm can accurately estimate extremely small unreliabilities and we exhibit large examples where it performs much better than existing approaches. It is also flexible enough to dispense with the frequently made assumption of independent edge failures.

We propose a new Monte Carlo method, based on dynamic importance sampling, to estimate the probability that a given set of nodes is connected in a graph (or network) where each link is failed with a given probability. The method generates the link states one by one, using a sampling strategy that approximates an ideal zero-variance importance sampling scheme. The approximation is based on minimal cuts in subgraphs. In an asymptotic rare-event regime where failure probability becomes very small, we prove that the relative error of our estimator remains bounded, and even converges to 0 under additional conditions, when the unreliability of individual links converges to 0. The empirical performance of the new sampling scheme is illustrated by examples.

A popular approach for modeling dependence in a finite-dimensional random vector X with given univariate marginals is via a normal copula that fits the rank or linear correlations for the bivariate marginals of X. In this approach, known as the NORTA method, the normal distribution function is applied to each coordinate of a vector Z of correlated standard normals to produce a vector U of correlated uniform random variables over (0,1); then X is obtained by applying the inverse of the target marginal distribution function for each coordinate of U. The fitting requires finding the appropriate correlation ρ between any two given coordinates of Z that would yield the target rank or linear correlation r between the corresponding coordinates of X. This root-finding problem is easy to solve when the marginals are continuous but not when they are discrete. In this paper, we provide a detailed analysis of this root-finding problem for the case of discrete marginals. We prove key properties of r and of its derivative as a function of ρ. It turns out that the derivative is easier to evaluate than the function itself. Based on that, we propose and compare alternative methods for finding or approximating the appropriate ρ. The case of discrete distributions with unbounded support is covered as well. In our numerical experiments, a derivative-supported method is faster and more accurate than a state-of-the-art, nonderivative-based method. We also characterize the asymptotic convergence rate of the function r (as a function of ρ) to the continuous-marginals limiting function, when the discrete marginals converge to continuous distributions.

The fastest long-period random number generators currently available are based on linear recurrences modulo 2. So far, software that provides multiple disjoint streams and substreams has not been available for these generators because of the lack of efficient jump-ahead facilities. In principle, it suffices to multiply the state (a k-bit vector) by an appropriate k × k binary matrix to find the new state far ahead in the sequence. However, when k is large (e.g., for a generator such as the popular Mersenne twister, for which k = 19,937), this matrix-vector multiplication is slow, and a large amount of memory is required to store the k × k matrix. In this paper, we provide a faster algorithm to jump ahead by a large number of steps in a linear recurrence modulo 2. The method uses much less than the k2 bits of memory required by the matrix method. It is based on polynomial calculus modulo the characteristic polynomial of the recurrence, and uses a sliding window algorithm for the multiplication.

We introduce TestU01 , a software library implemented in the ANSI C language, and offering a collection of utilities for the empirical statistical testing of uniform random number generators (RNGs). It provides general implementations of the classical statistical tests for RNGs, as well as several others tests proposed in the literature, and some original ones. Predefined tests suites for sequences of uniform random numbers over the interval (0, 1) and for bit sequences are available. Tools are also offered to perform systematic studies of the interaction between a specific test and the structure of the point sets produced by a given family of RNGs. That is, for a given kind of test and a given class of RNGs, to determine how large should be the sample size of the test, as a function of the generator's period length, before the generator starts to fail the test systematically. Finally, the library provides various types of generators implemented in generic form, as well as many specific generators proposed in the literature or found in widely used software. The tests can be applied to instances of the generators predefined in the library, or to user-defined generators, or to streams of random numbers produced by any kind of device or stored in files. Besides introducing TestU01 , the article provides a survey and a classification of statistical tests for RNGs. It also applies batteries of tests to a long list of widely used RNGs.

We propose a fast algorithm for computing the inverse symmetrical beta distribution. Four series (two around x = 0 and two around x = 1/2) are used to approximate the distribution function, and its inverse is found via Newton's method. This algorithm can be used to generate beta random variates by inversion and is much faster than currently available general inversion methods for the beta distribution. It turns out to be very useful for generating gamma processes efficiently via bridge sampling.

Fast uniform random number generators with extremely long periods have been defined and implemented based on linear recurrences modulo 2. The twisted GFSR and the Mersenne twister are famous recent examples. Besides the period length, the statistical quality of these generators is usually assessed via their equidistribution properties. The huge-period generators proposed so far are not quite optimal in this respect. In this article, we propose new generators of that form with better equidistribution and “bit-mixing” properties for equivalent period length and speed. The state of our new generators evolves in a more chaotic way than for the Mersenne twister. We illustrate how this can reduce the impact of persistent dependencies among successive output values, which can be observed in certain parts of the period of gigantic generators such as the Mersenne twister.

Linear congruential random-number generators with Mersenne prime modulus and multipliers of the form a = ±2 q ± r have been proposed recently. Their main advantage is the availability of a simple and fast implementation algorithm for such multipliers. This note generalizes this algorithm, points out statistical weaknesses of these multipliers when used in a straightforward manner, and suggests in what context they could be used safely.

We consider an ATM switch model, in which each of N sources is a Markov modulated rate process. We look at some approximations that have been proposed for p, the probability that in steady-state conditions, the buffer content exceeds x. These approximations are compared with the simulation results for p and the cell-loss ratio.

Usually, the t-dimensional spectral test for linear congruential generators examines the lattice structure of all the points formed by taking t successive values in the sequence. In this article, we consider the case where the t values taken are not successive, but separated by lags that are chosen a priori. For certain classes of linear congruential and multiple recursive generators, and for certain choices of the lags, we give lower bounds on the distance between hyperplanes. In some cases, those lower bounds are quite large, even in dimensions as small as t = 3. We give illustrations with specific classes of generators that have been proposed in the literature, and discuss the possible implications.

We discuss the implementation of theoretical tests to assess the structural properties of simple or combined linear congruential and multiple recursive random number generators. In particular, we describe a package implementing the so-called spectral and lattice tests for such generators. Our programs analyze the lattices generated by vectors of successive or nonsuccessive values produced by the generator, analyze the behavior of generators in high dimensions, and deal with moduli of practically unlimited sizes. We give numerical illustrations. We also explain how to build lattice bases in several different cases, e.g., for vectors of far-apart nonsuccessive values, or for sublattices generated by the set of periodic states or by a subcycle of a generator, and, for all these cases, how to increase the dimension of a (perhaps partially reduced) basis.

This is a commentary and a complement to the McGeoch’s feature article. Certain topics of current interest related to modeling, statistical output analysis, efficiency improvement, and random number generation, for stochastic discrete-event simulation is presented.

A portable set of software tools is described for uniform random variates generation. It provides for multiple generators running simultaneously, and each generator has its sequence of numbers partitioned into many long (disjoint) substreams. Simple procedure calls allow the user to make any generator “jump” ahead to the beginning of its next substream, back to the beginning of its current substream, or back to the beginning of its first substream.… Implementation issues are discussed.…A Pascal implementation for 32-bit computers is described. — From the Authors' Abstract

A numerical approach for computing optimal dynamic checkpointing strategies for general rollback and recovery systems is presented. The system is modeled as a Markov renewal decision process. General failure distributions, random checkpointing durations, and reprocessing-dependent recovery times are allowed. The aim is to find a dynamic decision rule to maximize the average system availability over an infinite time horizon. A computational approach to approximate such a rule is proposed. This approach is based on value-iteration stochastic dynamic programming with spline or finite-element approximation of the value and policy functions. Numerical illustrations are provided.>