Get Complete Project Material File(s) Now! »

## Importance Sampling

The limitation seen in the crude Monte Carlo estimator comes from the fact that the standard deviation of a Bernoulli random variable with parameter p becomes much larger than its expectation when the parameter p is small. In this context, the idea behind Importance Sampling is to consider an other random variable with the same mean but a lower variance. Reader is referred to [Glynn and Iglehart, 1989, Juneja and Shahabuddin, 2006] and references therein for more information about Importance Sampling.

Let Y be a real-valued random variable with distribution µY and finite expectation ˜ Y Y E [Y ] < ∞; and Y be a real-valued random variable with distribution µ ˜ such that µ is ˜ absolutely continuous with respect to µY : ˜ ∀B ∈ B(R), µY (B) = 0 ⇒ µY (B) = 0.

### Adaptive splitting

The optimal settings for the splitting method found in Section 1.3.1 are not usable in practice. Indeed they require to know in advance the target probability p and the cdf of g(X) to define the sequence (Fi)mi=1. Moreover it is not possible to generate iid. conditional samples from scratch. While it could be possible to use an oracle sequence (qi)mi=1 to define suboptimal subsets, or to get such a sequence from a first pilot run [Botev and Kroese, 2008, 2012], an other option is to define the splitting on-the-fly while the algorithm is running. These strategies are referred to as adaptive.

Subset Simulation The basic idea is to select the subsets iteratively with a heuristic approach based on a current population, i.e. a set of iid. samples. The optimality result for the conditional probability p∗0 gives a possible choice for this heuristic method: get qi as the crude Monte Carlo estimator of the p0 quantile, with p0 a given arbitrary value (often set to 0.1). A direct application of this principle leads to the Subset Simulation method (see Algorithm 2).

On lines 4 and 9 of Algorithm 2, qm+1 is the crude Monte Carlo estimator of a quantile of order 1 − p0, i.e. the d1 − p0e ordered statistic of the N0−samples [see for example Arnold et al., 1992]. The sequential approach of this algorithm does not only serve the definition of the subsets but also to sample from the conditional distributions (line 7). Indeed its link with Sequential Monte Carlo methods [Chopin, 2002, Del Moral et al., 2006, Le Gland, 2007, Vergé et al., 2013, Smith et al., 2013, Beskos et al., 2016] and Interacting.

#### Conditional sampling

In practice the adaptive splitting uses Markov chains to simulate according to the condi-tional distributions. At least this will increase the variance of the estimator. Very recently this issue has been tackled by Bréhier et al. [2015a] who show that the estimator remains unbiased in the dynamic case X = (Xt)t using the Markov property of the trajectory and successive appropriate filtrations; and by Cérou and Guyader [2016] who proved the same Central Limit Theorem as their previous result in the ideal case [see Cérou et al., 2012].

In this section we present commonly used tools to perform these simulations for both the static X ∈ X ⊂ Rd and the dynamic case (Xt)t ∈ X ⊂ (Rd)R. While a lot of ongoing research is devoted to the improvement of such methods and more advanced algorithms are now available, this thesis focuses on finding an optimal (in some sense specified throughout the manuscript) estimator given these conditional simulations are possible. In a sense, we focus on optimising the estimators while these researches strive to make them more robust. These two approaches are complementary and this distinction justifies the fact that the skilled reader may find that this section lacks depth. Metropolis-Hastings method When X ⊂ Rd for some d ≥ 1, one looks indeed for a simulator of the truncated distributions µX (· | g(X) > y) for all y in the support of µY . Assuming that µX has a density π with respect to the Lebesgue measure, it means that one wants to simulate X ∼ 1g(x)>yπ(x)/ P [g(X) > y]. A general idea for simulating from a given distribution known up to a given constant is to use the convergence property of a Markov chain to its unique invariant measure. This strategy is referred to as the Metropolis-Hastings algorithm in the literature, from the pioneering work of Metropolis and Ulam [1949] further extended by Hastings [1970]. This presentation follows [Delmas and Jourdain, 2006].

**Eﬃciency of the estimators**

Two diﬀerent notions are usually used to quantify the quality of an estimator in the rare event context: 1) its robustness, i.e. how it behaves when the probability becomes smaller; and 2) its reliability, i.e. which confidence intervals can be built and how reliable they are. Indeed, confidence intervals often require to estimate not only the sought probability but also a variance, and this may result in very bad confidence intervals. These notions will be presented here but the reader interested in more details and examples is referred to [Rubino et al., 2009, Section 4] and [L’Ecuyer et al., 2010] for even more advanced concepts.

Remember that the main problem of this thesis is the estimation of the quantity (1.1), it is: p = P [g(X) > q] for a given q. Assume that we have defined an estimator pb of p for any q, one considers the relative moment of order k of the estimator pb: mk(q) = E h k i. (1.27) pk p b.

**Table of contents :**

Contents

Context

**I Introduction **

**1 Monte Carlo methods for rare events **

1.1 Crude Monte Carlo method

1.1.1 Theoretical definition

1.1.2 Limitation

1.1.3 Practical implementation

1.2 Importance Sampling

1.3 Splitting

1.3.1 Ideal splitting

1.3.2 Adaptive splitting

1.3.3 Conditional sampling

1.4 Nested sampling

1.5 Efficiency of the estimators

**2 Rare event simulation and surrogate models **

2.1 Usual surrogate models

2.1.1 First/Second order reliability method

2.1.2 Support-Vector Machine

2.1.3 Polynomial-Chaos expansion

2.1.4 Kriging

2.2 Design of Experiments

2.2.1 First Design of Experiments

2.2.2 Model-oriented designs

2.2.3 Stepwise Uncertainty Reduction

2.3 Metamodels and estimators

2.3.1 Crude Monte Carlo estimator

2.3.2 Importance sampling-based procedures

2.3.3 Subset Simulation

**II Contribution to rare event simulation **

**3 Point process for rare event simulation **

3.1 Introduction

3.2 The increasing random walk

3.3 Probability estimator

3.3.1 Minimal variance unbiased estimators

3.3.2 Efficiency of the estimator

3.3.3 Confidence intervals

3.4 Quantile estimator

3.4.1 Description of the estimator

3.4.2 Statistical analysis of the estimator

3.5 Discontinuous random variables

3.5.1 The increasing random walk for discontinuous random variables

3.5.2 Probability estimators

3.6 Numerical examples

3.6.1 Discretised random path

3.6.2 Discrete random variable

3.7 Conclusion

**4 Nested sampling and rare event simulation **

4.1 Introduction

4.2 Ideal estimator

4.2.1 Extreme event simulation

4.2.2 Definition of the moment estimator

4.2.3 Comparison with classical Monte Carlo

4.3 Randomised unbiased estimator

4.3.1 Definition

4.3.2 Convergence rate

4.3.3 Optimal randomisation

4.3.4 Geometric randomisation

4.3.5 Parallel implementation

4.4 Application to heavy-tailed random variables

4.4.1 Exact resolution for a Pareto distribution

4.4.2 Comparison of the estimators

4.5 Examples

4.5.1 Practical implementation

4.5.2 Variance increase

4.5.3 Adaptive stopping criteria

4.5.4 Nested sampling with fixed computational budget

4.6 Conclusion

**5 Rare event simulation with random processes **

5.1 Rare event simulation

5.1.1 Augmented problem

5.1.2 Link with other metamodel based algorithms

5.1.3 Uncertainty reduction

5.1.4 Integrated SUR criteria

5.2 Algorithms

5.2.1 SUR criteria estimation

5.2.2 Integrated SUR criteria estimation

5.2.3 Bayesian Moving Particles

5.3 Numerical results

5.3.1 SUR criteria

5.3.2 Statistical error

5.3.3 Industrial problem

5.4 Conclusion

Conclusion and perspectives

**III Appendix **

**A Parallel computation of the estimators **

A.1 Introduction

A.2 Parallel algorithms

A.2.1 Sampling conditional distributions

A.2.2 Batch simulation of random walks

A.2.3 Fixed threshold

A.2.4 Fixed number of terms

A.2.5 Last Particle Algorithm

A.3 Wall-clock time

A.3.1 Fixed threshold algorithm

A.3.2 Fixed number of terms algorithm

A.4 Numerical benchmark of the parallelisation

A.4.1 Presentation of the examples

A.4.2 Estimation of failure probability

A.4.3 Estimation of quantile

A.5 Conclusion on parallel implementation

**Bibliography **