Results in 2008
In 2008 we have worked on different topics.
Splitting technique for rare event simulation
Splitting can be interpreted as decomposing the rare event as a succession of non rare ones and looking successively if these non-rare events are reached. From a given level/event, if the next event is not reached, so is the rare event, but if it is, we create an artificial drift toward the rare event by splitting (or cloning) the current state and trying to proceed further.
We have worked and are working in several directions. First, we have designed new versions of the procedure by for instance playing russian roulette when the trajectiry is going levels back, in order to reduce the computational time. Providing central limit theorems for the different versions is also of primary importance and a topic we are tackling. Similarly, some of our works are also trying to mix splitting and importance sampling.
Part of the work is also joint with F. Legland, from ASPI project-team.
Zero-Variance simulation
Variance reduction techniques can lower the estimation error in Monte Carlo simulation, sometimes by a large factor, but rarely change the convergence rate of the estimation error. This error usually decreases as the inverse square root of the computational effort, as dictated by the central limit theorem. In theory, there exist simulation estimators with zero variance, i.e., that always provide the exact value. The catch is that these estimators are usually much too difficult (or virtually impossible) to implement. However, there are situations, especially in the context of rare-event simulation, where the zero-variance simulation can be approximated well enough to provide huge efficiency gains. Adaptive versions can even yield a faster convergence rate, including exponential convergence in some cases.
We have provided an overview of these methods and discussed their practicality.
Zero-Variance importance sampling for highly reliable systems
We have designed efficient importance sampling (IS) schemes to estimate the reliability of a class of Markov chain models that includes the highly reliable Markovian systems (HRMS) often used to represent the evolution of multicomponent systems in reliability settings. For these models, there is in fact a zero-variance IS scheme that can be written exactly in terms of a value function that gives the expected cost-to-go (the exact reliability, in our case) from any state of the chain.
This IS scheme is impractical to implement exactly, but it can be approximated by approximating this value function. In our implementation, we start with a simple crude approximation of the value function, we use it in a first-order IS scheme to obtain a better approximation at a few selected states, then we interpolate in between and use this interpolation in our final (second-order) IS scheme.
Our approach outperforms the popular IS heuristics previously proposed for this class of problems. In an asymptotic analysis, we also show that with our approximation, the IS estimator has bounded relative error (BRE) under very mild conditions, and vanishing relative error (VRE).
Rare event simulation for static models
We have considered the problem of evaluating the probability of a rare event in the case of a static model, that is, a model where time does not play an explicit role. More specifically, we address the problem of network reliability estimation using Monte Carlo, where the network is modeled by a graph with nodes and/or edges subject to independent failures. When the reliability of the components (the probability that they are operational) is close to one, the path redundancy in the graph makes that the probability of the existence of at least a path between two selected nodes (for instance) is extremely high, and thus, that the unreliability of the whole network (the probability that the two nodes are not connected) is very small.
We have designed efficient Monte Carlo methods for the estimation of the unreliability of such a network, by combining three approaches. First, as in other works, we use easy-to-get knowledge about the network, namely its path structure, to follow a conditional approach allowing to bound the target. In particular, we show how to derive methods having bounded work-normalized relative variance in the homogeneous-lines case. Second, we explore the Randomized Quasi-Monte Carlo (RQMC) technique in this context, in order to accelerate the estimations. Third, we show that the so-called zero-variance importance sampling approximation, well established in the field of Markovian models, for instance, can be extended to this static model, and we explore the properties of the obtained procedure. The idea here is to introduce artificially a (discrete) time variable given by an order on the set of network components. A first component is chosen at time 1, then a second one at time 2, and so on.
Robutsness properties of estimators
The asymptotic robustness of estimators as a function of a rarity parameter, in the context of rare-event simulation, is often qualified by properties such as bounded relative error (BRE) and logarithmic efficiency (LE), also called asymptotic optimality. However, these properties do not suffice to ensure that moments of order higher than one are well estimated. For example, they do not guarantee that the variance of the empirical variance remains under control as a function of the rarity parameter.
We study generalizations of the BRE and LE properties that take care of this limitation. They are named bounded relative moment of order k (BRM-k) and logarithmic efficiency of order k (LE-k), where k>=1 is an arbitrary real number. We also introduce and examine a stronger notion called vanishing relative centered moment of order k, and exhibit examples where it holds. These properties are of interest for various estimators, including the empirical mean and the empirical variance. We develop (sufficient) Lyapunov-type conditions for these properties in a setting where state-dependent importance sampling (IS) is used to estimate first-passage time probabilities. We show how these conditions can guide us in the design of good IS schemes, that enjoy convenient asymptotic robustness properties, in the context of random walks with light-tailed and heavy-tailed increments. As another illustration, we study the hierarchy between these robustness properties (and a few others) for a model of highly-reliable Markovian system (HRMS) where the goal is to estimate the failure probability of the system. In this setting, for a popular class of IS schemes, we show that BRM-k and LE-k are equivalent and that these properties become strictly stronger when k increases. We also obtain a necessary and sufficient condition for BRM-k in terms of quantities that can be readily verified from the parameters of the model.
Array-RQMC methods
Straightforward simulation, which corresponds to using the Monte Carlo (MC) methods to estimate multivariate integrals, often gives very noisy answers. Quasi-Monte Carlo (QMC) and randomized quasi-Monte Carlo (RQMC) methods can beat MC for estimating the integral of a smooth function in a small or moderate number of dimensions. But when the function depends on the sample path of a Markov chain over a large number of steps, which is often the case in simulation applications, the dimension of the integrand is very large and common wisdom says that QMC and RQMC methods are ineffective.
Our work proposes a new RQMC algorithm specially designed to improve the efficiency of Markov chains simulations over several steps. We provide examples showing that their method can dramatically reduce the variance compared with standard MC. Applications include queueing systems, option pricing in finance, reliability and risk assessment models, and many more.
The key component of the algorithm is to properly sort the simulated chains to benefit from the good distribution of QMC sequences. We are now working on defining ``good'' sortings in multidimensional settings.
Another issue we curretnly tackling is the coverage of the confidence intervals in RQMC methods in general. Indeed, the common prectice is to take a few randomizations to get a confidence interval and consider as many QMC points as possible to get a faster convergence. On the other hand, there is no study or assessment on resulting the coverage. We have made interesting numerical observations and are working on thoretical justifications.
Resuts in 2009
Completing the book on rare event simulation
The beginnong of 2009 focused on the completion of the
book on rare event simulation.
The (joint) contribution from the associated team is composed of the two key chapters describing the methodological aspects:
- importance sampling,
- splitting.
Acceleration technques for reliability analysis of static models
Many dependability analyses are performed using static models, that is, models where time is not an explicit variable. In these models, the system and its components are considered at a fixed point in time, and the word “static” means that the past or future behavior is not relevant for the analysis. Examples of such models are reliability diagrams, or fault trees. The main difficulty when evaluating the dependability of these systems is the combinatorial explosion associated with exact solution techniques. For large and complex models, one may turn to Monte Carlo methods, but these methods have to be modified or adapted in the presence of rare important events, which are commonplace in reliability and dependability systems. This chapter examines a recently proposed method designed to deal with the problem of estimating reliability metrics for highly-dependable systems where the failure of the whole system is a rare event.We focus on the robustness properties of estimators. We also propose improvements to the original technique, including its combination with randomized quasi-Monte Carlo, for which we prove that the variance converges at a faster rate (asymptotically) than for standard Monte Carlo.
Limiting distributions of randomized quasi-Monte Carlo methods
Randomized quasi-Monte Carlo (RQMC) methods estimate the expectation of a random variable by the average of n dependent realizations of it. In general, due to the strong dependence, the estimation error may not obey a central limit theorem. Analysis of RQMC methods have so far focused mostly on the convergence rates of asymptotic worst-case error bounds and variance bounds, when n tends to infinity, but little is known about the limiting distribution of the error. Here we examine this limiting distribution for the special case of a randomly-shifted lattice rule, when the integrand is smooth. We start with simple one-dimensional functions, where we show that the limiting distribution is uniform over a bounded interval if the integrand is non-periodic, and has a square root form over a bounded interval if the integrand is periodic. In higher dimensions, for linear functions, the distribution function of the properly standardized error converges to a spline of degree equal to the dimension.
Quasi-Monte Carlo methods and discrete choice models
Advanced discrete choice models, such as parametric/non-parametric mixed logit and hybrid choice models, have been heavily promoted in research since the 2000 Nobel lecture given by McFadden. The complexity of model formulations has nevertheless been limited by the associated estimation difficulties. Even if important progress has been attained these last years, researchers and practitioners still lack of suitable tools to estimate such complex models and often report high levels of failure. In this paper, we try to take advantage of various approaches proposed in the literature to alleviate the estimation cost and we propose to combine them to construct more efficient estimators.
We consider the adaptive trust-region algorithm designed by Bastin et al., including simulation bias correction and a variant based on retrospective trust-region, that simplifies the implementation of the method. The number of draws used at each iteration is indeed then based on the simulated log-likelihood reduction instead of a model of it. This approach has been limited so far to the use of Monte Carlo draws only. But randomized quasi-Monte Carlo sometimes achieves better accuracy with less draws. The adaptive strategy however relies on error estimation, an approach made possible by using randomized quasi-Monte Carlo sequences. The choice of the quasi-Monte Carlo sequence and the randomization technique may have a large impact on the estimation performance. While Halton sequences are popular in transportation research, they suffer from correlation amongst coordinates that only adequate scrambling can remove, while shuffling techniques have been shown to be not so efficient. Scrambled Halton sequences however suffer from the use of a different basis in each dimension, making difficult to adequately choose a good number of draws to uniformly cover the integration space. This also makes them unsuited to the adaptive strategy. We consider here digital sequences using basis two, so that a number of draws equal to a power of two ensures a good coverage of the integration space. It is then easy to construct an adaptive strategy using only adequate quasi-Monte Carlo samples. A particular but popular kind of digital sequence in base two is Sobol sequences that nevertheless require careful choice of primitive polynomials. Sobol sequences usually achieve best results in numerous applications, but they have only received limited attention in discrete choice modeling.
We also investigate memory consumption and computational time when (quasi-)Monte Carlo draws are stored, which has become a classical way to proceed in discrete choice modeling. The generation of long draws sequences can quickly require a large amount of memory, sometimes making the estimation at a desired accuracy impossible. In this paper we suggest to save memory by generating quasi-Monte Carlo draws on the fly, and check that this way to proceed ensures similar estimation times. It is also known that independent draws across individuals are needed in order to be able to express the limit distribution of the simulated likelihood function. Randomization is a way to transform deterministic samples into random and independent samples, but (usually) with a lower variance than what is obtained with standard Monte-Carlo draws. We consider here the possibility of generating the same quasi-Monte Carlo sequence for each individual, but applying a different randomization on the fly. By doing so, we guarantee independence and we drastically limit memory consumption, even when we decide to store quasi-Monte Carlo draws before the randomization. By analogy to scrambling across observations, we can express this feature as randomization across individuals. Scrambling does not always guarantee statistical independence between the individual draw sets, especially when the only target of the scrambling is to break correlations. Interpreting scrambling as a randomization technique depends on how the scrambling is defined. We finally compare the results obtained with this methodology to the case where a single long sequence of draws divided amongst individuals is generated.
Numerical experiments, on both simulated and real datasets, will allow us to (1) demonstrate the advantages of a method based on these various improvements, (2) estimate complex models based on a large number of observations and (3) improve the accuracy of the estimations when compared to the results obtained with standard techniques.
Results in 2010
Randomized quasi-Monte Carlo methods
We have pursued our work on the limiting distributions of randomized quasi-Monte Carlo methods (see the 2009 results). We have generalized the results about the random shift and more deeply analyzed and illustrated the convergence properties in dimesion larger than one. The results are submitted (under revision) at the Electronic Journal of Statistics.
Zero-Variance Importance Sampling and dependability models
We are also continuing our activity around dependability models estimated by Zero-Variance Importance Sampling approximation. We have first modified our key work for static models, under minor revision at IEEE Transactions on Reliability
We have also proposed to combine two efficient rare event Monte Carlo simulation techniques for the estimation of the connectivity probability of a given set of nodes in a graph when links can fail: approximate zero-variance importance sampling and conditional Monte Carlo taking into account the expected number of samples over which a prespecified set of disjoint minpaths linking the set of nodes fails. Those two methods had been applied separately, but we show here how their combination can be implemented, we derive asymptotic robustness properties of the resulting estimator when reliabilities of individual links go arbitrarily close to one, and we illustrate numerically the gain that can be obtained on large graphs.
New activity around robustness properties of adaptive variance reduction techniques
We are also starting a new activity about the robustness of rare event simulation estimators. Robustness properties for rare event simulation have been studied in the litterature when parameters (that is, importance sampling (IS) parameters, or importance function and levels for splitting) are fixed. But there exist many adaptive techniques where parameters are learned and as a consequence random, since depending on sample values. Cross-Entropy technique is a typical example. In most cases, it is then not possible to assert that a robustness property is verified. The goals are here
- to illustrate that getting a robustness property may be random for adaptive methods;
- to define probabilistic robustness properties;
- to provide conditions under which those probabilitic properties are verified.
Results in 2011
We continue our activities around the robustness of adaptive estimators, partly in collaboration with Ad Ridder (Amsterdam). Other obtained and published results are described now.
Static network reliability estimation and generalized splitting
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 is failed 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.
Static network reliability estimation and importance sampling
We speed up the Monte Carlo simulation of static graph reliability models by adding graph reductions to zero-variance importance sampling (ZVIS) approximation techniques. ZVIS approximation samples the status of links sequentially, and at each step we check if series-parallel reductions can be performed. We present two variants of the algorithm and describe their respective advantages. We show that the method satisfies robustness properties as the reliability of links increases. We illustrate theoretically on small examples and numerically on large ones the gains that can be obtained, both in terms of variance and computational time.
An importance sampling method based on a one-step look-ahead density from a Markov chain
We propose a new importance sampling method that constructs an importance sampling density which approximates the zero-variance sampling density nonparametrically as follows. In a first stage, it generates a sample (possibly approximately) from the zero-variance density using, for example, Markov chain Monte Carlo methodology. In a second stage, the method constructs a kernel density estimator of the zero-variance density based on the sample in the first stage. The most important aspect of the method is that, unlike other kernel estimation methods, the kernel of the estimator is defined as the one-step transition density of a Markov chain whose stationary distribution is the zero-variance one. We give examples where this one-step transition density is available analytically and provide numerical illustrations in which the method performs very well.
Randomized Quasi-Monte Carlo and mixed logit models
We examine the effectiveness of randomized quasi-Monte Carlo (RQMC) techniques to estimate the integrals that express the discrete choice probabilities in a mixed logit model, for which no closed form formula is available. These models are used extensively in travel behavior research. We consider popular RQMC constructions such as randomized Sobol’, Faure, and Halton points, but our main emphasis is on randomly-shifted lattice rules, for which we study how to select the parameters as a function of the considered class of integrands. We compare the effectiveness of all these methods and of standard Monte Carlo (MC) to reduce both the variance and the bias when estimating the log-likelihood function at a given parameter value. In our numerical experiments, randomized lattice rules (with carefully selected parameters) and digital nets are the best performers and they reduce the bias as much as the variance. With panel data, in our examples, the performance of all RQMC methods degrades rapidly when we simultaneously increase the dimension and the number of observations per individual.
Program for 2012
Joint organization of the Analysis track at the Winter Simulation Conference 2012, for the first time it leaves North America (Berlin, Germany).
For rare event simulation we would like to focus on dependent events, a key aspect in the domain. Some dependences have already been modeled and analyzed in queuing systems, thanks to Markov modulated processes or batch arrivals, but not in the dependability context, for which the existing methods have to be adapted or redesigned. Failure propagations are barely considered in highly reliable Markovian systems analysis, and never in static reliability systems, whereas a realistic assumption. Another issue we want to address is the design of adaptive methods, where variance reduction parameters are learned to better approach the ideal zero-variance estimator. Dependences are here not on the model itself, but rather internal to the method because the parameters are random and depend on sample values. Several adaptive methods exist, but we would like to design new ones improving the zero-variance importance sampling approximation we have developped and already the most efficient estimator in the contexts we have considered. Robustness to rarity has to be defined for adaptive methods, in a probabilistic way.
Modeling and simulation of complex revenue management models in telecommunications. The models in this industry are becoming more and more complex. Companies have to simultaneously consider heterogeneous networks (WiFi, WiMAX, fixed, 3G wireless), all potential applications, and the corresponding imbricated demand of customers, and propose offers, potentially grouped and trying to minimize churn. Pricing can be dynamic, flat or usage-based, and also differentiated. Capacity investment has also to be studied at a larger time scale, but can also include leasing from competitors (typically what Mobile Virtual Network Operators are doing). Each decision in a category (capacity investiment, list of offers, price, contracts) has an impact on demand and cannot be handled separately. Due to the complexity of the model, simulation is the relevant tool for the analysis. We expect to draw upon the results obtained in the call centers domain to design a model and propose a simulation based on MC and RQMC methods. Queuing systems will be used to represent demand, but modulated by the contract offers. Churn will also depend on the delivered QoS. Decisions of customers can be modeled by discrete choice models, making also use of our previous joint work on this issue to determine by RQMC the distribution of users' preferences.
