An Introduction to Rare Event Simulation and Importance Sampling

Gino Biondini*,1    * Department of Mathematics, State University of New York at Buffalo, Buffalo, New York, USA
1 Corresponding author: email address: biondini@buffalo.edu

1 Introduction: Monte Carlo Methods, Rare Event Simulation, and Variance Reduction Techniques

Since its introduction almost 70 years ago (Metropolis and Ulam, 1949) (see Metropolis, 1987 for a historical review), the Monte Carlo (MC) method has been extensively used in engineering and scientific computing. In their most general interpretation, MC methods are a way to compute integrals. They comprise a collection of techniques for generating random samples on a computer as well as their application to solve a variety of problems. In essence, they involve drawing random or pseudo-random samples from a specific distribution and using them to estimate one or more quantities of interest. Such methods are especially advantageous over numerical quadrature methods when the dimensionality of the problem is large. As a result, and thanks to their flexibility, such methods have found a wide range of applications (e.g., see Fishman, 1996; Fishman, 2006; Kroese et al., 2011; Landau and Binder, 2000).

A common challenge in MC simulation is that of rare event simulation, also referred to as the problem of rare events, where very small probabilities need to be accurately estimated—for example, in reliability analysis, or performance analysis of telecommunication systems. In a nutshell, the problem is that if one needs to quantify the probability of one or more events that occur very rarely, an exceedingly large number of samples are needed even to just produce the desired events, and an even larger number of samples are required to obtain accurate estimates. Other applications that call for rare event simulation are queueing systems (to avoid excessively long waiting times), nuclear physics (avoiding catastrophic accident), security systems (false alarms in radar), material science (technical defects), mathematical science, and insurance.

One approach to overcome the problem of rare events is the use of variance reduction techniques (VRTs) (e.g., see the monographs: Bucklew, 2004; Fishman, 1996; Kroese et al., 2011 for general reviews). The general idea behind all of these techniques is to modify the selection of the random samples in such a way that the desired events occur more frequently than they would normally, while simultaneously taking these changes into account in order to obtain unbiased estimates.

Perhaps the most famous VRT is the importance sampling (IS) (Fishman, 1996; Kroese et al., 2011; Srinivasan, 2002). The main idea behind IS is to select an appropriate biasing distribution (i.e., a change of probability measure) from which to draw the MC samples so that most of the distribution mass falls on the regions of interest. This ensures that many of the MC samples will produce the rare events sought. At the same time, the contribution from each sample is weighted according to the likelihood ratio, which ensures that unbiased estimates are obtained.

Of course, for IS to be effective, a good biasing distribution must be chosen. This requires knowledge of which system configurations are likely to produce the rare events of interest. Even though such knowledge is not always available, in many cases it is enough to leverage what is known about the system’s behavior in order to guide the choice of biasing distribution, and indeed IS has been used with success in a variety of applications (Biondini et al., 2004; Li et al., 2007; Moore et al., 2008; Smith et al., 1997). (Note that, often, exact knowledge of the most likely failure configurations may not be needed, and an approximate knowledge may be sufficient, since the statistical nature of the MC sampling allows one to take into account the contributions of nearby points in sample space.)

Many other VRTs have also been used with success in various applications, such as multicanonical MC methods (Yevick, 2002), Markov chain Monte Carlo (MCMC) methods (Secondini and Forestieri, 2005), and Gibbs sampling. See Fishman (1996), Landau and Binder (2000), and MacKay (2003) for a general overview of these methods. The common thread among those VRTs is that they are adaptive. In essence, such methods attempt to find the important regions of sample space numerically. These methods can be applied to problems for which no good choice of biasing distribution is known. When IS is available, however, it is generally advantageous over other methods, because: (i) IS allows one to compute precise error estimates, if desired; (ii) adaptive methods typically require tweaking certain parameters, on which IS is less dependent; and (iii) IS is usually faster than adaptive methods, since, in adaptive methods, a certain portion of numerical simulations needs to be used to look for the most important regions in state space. Indeed, the speed advantage of IS was verified directly in a few cases by a detailed comparison between different methods (Biondini and Kath, 2005; Lima et al., 2005). We should also mention that it is not always necessary to choose between IS and adaptive VRTs. Indeed, yet another technique which has proven to be especially useful in recent years is the cross-entropy method (de Boer et al., 2005; Rubinstein and Kroese, 2004). While it is a useful VRT on its own right, in some cases, the cross-entropy method can also be combined with IS to afford the user the advantages of both IS and those of adaptive techniques.

The remainder of this chapter aims to put the above discussion on a more precise mathematical setting.

2 MC Methods and the Problem of Rare Events

2.1 MC Estimators

We start with a simple one-dimensional (1D) example. Let X be a random variable (RV) with probability density function (pdf) px(x) (Papoulis, 1991). If one defines Y = y(X), where (x)=∫−∞xpx(x)dx, it is easy to show that Y is uniform in [0,1]. [To see this, note that py(y) dy = px(x) dx, with dy = (dy/dx) dx. But dy/dx = px(x), so py(y) = 1.]

Now suppose that we wish to calculate the probability Q that X falls in a range R of interest, namely =P[X∈R], where ⊂R. We can write Q as

=∫IR(x)px(x)dx.

The function IR(x) is the so-called indicator function (or characteristic function) of the set R: namely, IR(x) = 1 if x ∈ R and IR(x) = 0 otherwise. (Hereafter we will drop the subscript R on I whenever that will not cause ambiguity. Also, integrals without limits are always intended as complete—i.e., over all of sample space—unless specifically noted otherwise.) In particular, we are interested in situations in which it is difficult to compute the above integral analytically.

Making the substitution xy, we can express Q as =∫I(x(y))dy. It is therefore natural to try to estimate Q using a frequency count. That is, we draw N independent identically distributed (i.i.d.) random samples 1,…,YN from a uniform distribution and we write the estimate ^N=F/N, where F is the number of samples which fall in the region of interest. More specifically, the above MC estimator is ^N=(1/N)∑n=1NI(x(Yn)). Equivalently, we can forget about Y and write the above estimator as

^N=1N∑n=1NI(Xn),

where the i.i.d. random samples 1,…,XN are drawn according to the distribution...