mistral (version 2.2.2)

MP: Moving Particles

Description

This function runs the Moving Particles algorithm for estimating extreme probability and quantile.

Usage

MP(
  dimension,
  lsf,
  N = 100,
  N.batch = foreach::getDoParWorkers(),
  p,
  q,
  lower.tail = TRUE,
  Niter_1fold,
  alpha = 0.05,
  compute_confidence = FALSE,
  verbose = 0,
  chi2 = FALSE,
  breaks = N.batch/5,
  ...
)

Value

An object of class list containing the outputs described below:

p

the estimated probability or the reference for the quantile estimate.

q

the estimated quantile or the reference for the probability estimate.

cv

the coefficient of variation of the probability estimator.

ecdf

the empirical cdf.

L

the states of the random walk.

L_max

the farthest state reached by the random process. Validity range for the ecdf is then (-Inf, L_max] or [L_max, Inf).

times

the times of the random process.

Ncall

the total number of calls to the lsf.

X

the N particles in their final state.

y

the value of the lsf(X).

moves

a vector containing the number of moves for each batch.

p_int

a 95% confidence intervalle on the probability estimate.

cov

the coefficient of variation of the estimator

q_int

a 95% confidence intervall on the quantile estimate.

chi2

the output of the chisq.test function.

Arguments

dimension

the dimension of the input space.

lsf

the function defining the RV of interest Y = lsf(X).

N

the total number of particles,

N.batch

the number of parallel batches for the algorithm. Each batch will then have N/N.batch particles. Typically this could be detectCores() or some other machine-derived parameters. Note that N/N.batch has to be an integer.

p

a given probability to estimate the corresponding quantile (as in qxxxx functions).

q

a given quantile to estimate the corresponding probability (as in pxxxx functions).

lower.tail

as for pxxxx functions, TRUE for estimating P(lsf(X) < q), FALSE for P(lsf(X) > q).

Niter_1fold

a function = fun(N) giving the deterministic number of iterations for the first pass.

alpha

when using default Niter_1fold function, this is the risk not to have simulated enough samples to produce a quantile estimator.

compute_confidence

if TRUE, the algorithm runs a little bit longer to produces a 95% interval on the quantile estimator.

verbose

to control level of print (either 0, or 1, or 2).

chi2

for a chi2 test on the number of events.

breaks

for the final histogram is chi2 == TRUE.

...

further arguments past to IRW.

Author

Clement WALTER

Details

MP is a wrap up of IRW for probability and quantile estimation. By construction, the several calls to IRW are parallel (foreach) and so is the algorithm. Especially, with N.batch=1, this is the Last Particle Algorithm, which is a specific version of SubsetSimulation with p_0 = 1-1/N. However, note that this algorithm not only gives a quantile or a probability estimate but also an estimate of the whole cdf until the given threshold q.

The probability estimator only requires to generate several random walks as it is the estimation of the parameter of a Poisson random variable. The quantile estimator is a little bit more complicated and requires a 2-passes algorithm. It is thus not exactly fully parallel as cluster/cores have to communicate after the first pass. During the first pass, particles are moved a given number of times, during the second pass particles are moved until the farthest event reach during the first pass. Hence, the random process is completely simulated until this given state.

For an easy user experiment, all the parameters are defined by default with the optimised values as described in the reference paper (see References below) and a typical use will only specify N and N.batch.

References

• A. Guyader, N. Hengartner and E. Matzner-Lober:
  Simulation and estimation of extreme quantiles and extreme probabilities
  Applied Mathematics and Optimization, 64(2), 171-196.
  

• C. Walter:
  Moving Particles: a parallel optimal Multilevel Splitting method with application in quantiles estimation and meta-model based algorithms
  Structural Safety, 55, 10-25.
  

• E. Simonnet:
  Combinatorial analysis of the adaptive last particle method
  Statistics and Computing, 1-20.
  

See Also

SubsetSimulation MonteCarlo IRW

Examples

if (FALSE) {
# Estimate some probability and quantile with the parabolic lsf
p.est <- MP(2, kiureghian, N = 100, q = 0) # estimate P(lsf(X) < 0)
p.est <- MP(2, kiureghian, N = 100, q = 7.8, lower.tail = FALSE) # estimate P(lsf(X) > 7.8)

q.est <- MP(2, kiureghian, N = 100, p = 1e-3) # estimate q such that P(lsf(X) < q) = 1e-3
q.est <- MP(2, kiureghian, N = 100, p = 1e-3, lower.tail = FALSE) # estimate q such
# that P(lsf(X) > q) = 1e-3


# plot the empirical cdf
plot(xplot <- seq(-3, p.est$L_max, l = 100), sapply(xplot, p.est$ecdf_MP))

# check validity range
p.est$ecdf_MP(p.est$L_max - 1)
# this example will fail because the quantile is greater than the limit
tryCatch({
   p.est$ecdf_MP(p.est$L_max + 0.1)},
   error = function(cond) message(cond))

# Run in parallel
library(doParallel)
registerDoParallel()
p.est <- MP(2, kiureghian, N = 100, q = 0, N.batch = getDoParWorkers())
}